Line data Source code
1 : module test_atm_support
2 :
3 : use const_def
4 : use math_lib
5 : use atm_def
6 : use atm_lib
7 : use chem_lib, only: composition_info
8 : use chem_def
9 : use eos_def
10 : use eos_lib
11 : use kap_def
12 : use kap_lib
13 :
14 : implicit none
15 :
16 : logical, parameter :: SKIP_PARTIALS = .TRUE.
17 :
18 : logical :: test_verbosely
19 : integer :: eos_handle, kap_handle
20 : real(dp) :: cgrav, Pextra_factor
21 : ! composition info
22 : integer, parameter :: species = 7
23 : integer, pointer :: chem_id(:), net_iso(:)
24 : real(dp) :: X, Y, Z, XC, XN, XO, xa(species), abar, zbar, z53bar
25 :
26 : integer :: ierr, max_iters, max_steps, iters
27 : real(dp) :: logg, Teff, M, R, L, kap_guess, tau, tau_base, tau_phot, &
28 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, &
29 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, &
30 : atol, rtol, kap, T, err, P, Prad, &
31 : dlnT_dlnkap, dlnP_dlnkap, &
32 : logTeff, log_gsurf, log_Rsurf, log_M, g, &
33 : T_eq, kap_v, gamma, T_int
34 :
35 : contains
36 :
37 6218 : subroutine do_test_atm( &
38 : test_verbosely_in, cgrav_in, eos_handle_in, kap_handle_in)
39 : logical, intent(in) :: test_verbosely_in
40 : real(dp), intent(in) :: cgrav_in
41 : integer, intent(in) :: eos_handle_in, kap_handle_in
42 :
43 2 : test_verbosely = test_verbosely_in
44 2 : cgrav = cgrav_in
45 2 : eos_handle = eos_handle_in
46 2 : kap_handle = kap_handle_in
47 :
48 2 : call test_table(ATM_TABLE_TAU_1M1, 'tau=0.1')
49 2 : call test_table(ATM_TABLE_TAU_1, 'tau=1')
50 2 : call test_table(ATM_TABLE_TAU_10, 'tau=10')
51 2 : call test_table(ATM_TABLE_TAU_100, 'tau=100')
52 :
53 2 : call test_table(ATM_TABLE_WD_TAU_25, 'WD tau=25')
54 2 : call test_table(ATM_TABLE_DB_WD_TAU_25, 'DB WD tau=25')
55 2 : call test_table(ATM_TABLE_PHOTOSPHERE, 'photosphere')
56 :
57 2 : call test_T_tau_varying(ATM_T_TAU_EDDINGTON, 'Eddington', -1._dp)
58 2 : call test_T_tau_varying(ATM_T_TAU_KRISHNA_SWAMY, 'Krishna-Swamy', -1._dp)
59 2 : call test_T_tau_varying(ATM_T_TAU_SOLAR_HOPF, 'solar Hopf', -1._dp)
60 2 : call test_T_tau_varying(ATM_T_TAU_TRAMPEDACH_SOLAR, 'Trampedach solar', -1._dp)
61 2 : call test_T_tau_varying(ATM_T_TAU_EDDINGTON, 'Eddington', 100._dp)
62 :
63 2 : call test_T_tau_uniform('fixed', 100._dp)
64 2 : call test_T_tau_uniform('iterated', 150._dp)
65 :
66 2 : call test_irradiated()
67 :
68 2 : end subroutine do_test_atm
69 :
70 : !****
71 :
72 14 : subroutine test_table(table_id, label)
73 :
74 : integer, intent(in) :: table_id
75 : character(*), intent(in) :: label
76 :
77 : include 'formats'
78 :
79 14 : ierr = 0
80 :
81 14 : if (test_verbosely) then
82 7 : write (*, *)
83 7 : write (*, *) 'test_table: ', label
84 : end if
85 :
86 14 : if (table_id == ATM_TABLE_WD_TAU_25) then
87 2 : logg = 7.1d0
88 2 : Teff = 5000
89 2 : M = 0.8d0*Msun
90 2 : R = sqrt(cgrav*M/exp10(logg))
91 : !write(*,*) 'R/Rsun', R/Rsun
92 2 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
93 : !write(*,*) 'L/Lsun', L/Lsun
94 12 : elseif (table_id == ATM_TABLE_DB_WD_TAU_25) then
95 2 : logg = 8.0d0
96 2 : Teff = 25000
97 2 : M = 0.8d0*Msun
98 2 : R = sqrt(cgrav*M/exp10(logg))
99 : !write(*,*) 'R/Rsun', R/Rsun
100 2 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
101 : !write(*,*) 'L/Lsun', L/Lsun else
102 : else
103 10 : M = Msun
104 10 : R = Rsun
105 10 : L = Lsun
106 10 : Teff = atm_Teff(L, R)
107 : end if
108 :
109 14 : Z = 0.02d0
110 :
111 : call atm_eval_table( &
112 : L, R, M, cgrav, table_id, Z, SKIP_PARTIALS, &
113 : Teff, &
114 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
115 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
116 14 : ierr)
117 14 : if (ierr /= 0) then
118 0 : if (test_verbosely) write (*, *) 'failed in atm_eval_table'
119 0 : call mesa_error(__FILE__, __LINE__)
120 : end if
121 :
122 14 : T = exp(lnT)
123 14 : P = exp(lnP)
124 :
125 14 : Prad = radiation_pressure(T)
126 :
127 14 : if (test_verbosely) write (*, *)
128 14 : if (test_verbosely) write (*, '(99a16)') &
129 7 : 'T', 'log_T', 'log P', 'M/Msun', 'L/Lsun', 'R/Rsun', 'logPgas'
130 14 : if (test_verbosely) write (*, '(i16,99f16.8)') &
131 7 : floor(0.5d0 + T), lnT/ln10, log10(P), M/Msun, L/Lsun, R/Rsun, log10(P - Prad)
132 14 : if (test_verbosely) write (*, *)
133 :
134 14 : end subroutine test_table
135 :
136 : !****
137 :
138 10 : subroutine test_T_tau_varying(T_tau_id, label, tau_base_in)
139 :
140 : integer, intent(in) :: T_tau_id
141 : character(*), intent(in) :: label
142 : real(dp), intent(in) :: tau_base_in
143 :
144 10 : real(dp) :: errtol
145 : integer :: max_steps
146 :
147 : include 'formats'
148 :
149 10 : ierr = 0
150 :
151 10 : if (tau_base_in < 0._dp) then
152 8 : call atm_get_T_tau_base(T_tau_id, tau_base, ierr)
153 8 : if (ierr /= 0) then
154 0 : if (test_verbosely) write (*, *) 'failed in atm_eval_T_tau_varying'
155 0 : return
156 : end if
157 : else
158 2 : tau_base = tau_base_in
159 : end if
160 :
161 10 : if (test_verbosely) then
162 5 : write (*, *)
163 5 : write (*, *) 'test_T_tau_varying: ', label
164 : end if
165 :
166 : ! TEST SOLAR VALUES
167 10 : logTeff = log10(5776d0)
168 10 : log_gsurf = log10(cgrav*Msun/(Rsun*Rsun))
169 10 : log_Rsurf = log10(Rsun)
170 10 : log_M = log10(Msun)
171 :
172 10 : Z = 0.02d0
173 10 : X = 0.70d0
174 10 : XC = 3.2724592105263235d-03
175 10 : XN = 9.5023842105263292d-04
176 10 : XO = 8.8218000000000601d-03
177 10 : call set_composition()
178 :
179 10 : R = exp10(log_Rsurf)
180 10 : M = exp10(log_M)
181 :
182 10 : g = exp10(log_gsurf)
183 10 : Teff = exp10(logTeff)
184 10 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
185 :
186 10 : errtol = 1.E-9_dp
187 10 : max_steps = 500
188 :
189 : call atm_eval_T_tau_varying( &
190 : tau_base, L, R, M, cgrav, &
191 : T_tau_id, eos_proc, kap_proc, &
192 : errtol, max_steps, SKIP_PARTIALS, &
193 : Teff, &
194 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
195 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
196 10 : ierr)
197 10 : if (ierr /= 0) then
198 0 : if (test_verbosely) write (*, *) 'failed in atm_eval_T_tau_varying'
199 0 : call mesa_error(__FILE__, __LINE__)
200 : end if
201 :
202 10 : T = exp(lnT)
203 10 : P = exp(lnP)
204 :
205 10 : if (test_verbosely) write (*, *)
206 10 : if (test_verbosely) write (*, '(99a16)') 'T', 'log T', 'log P', 'M/Msun', 'L/Lsun', 'X', 'Z'
207 10 : if (test_verbosely) write (*, '(i16,99f16.8)') &
208 5 : floor(0.5d0 + T), log10(T), log10(P), M/Msun, L/Lsun, X, Z
209 10 : if (test_verbosely) write (*, *)
210 :
211 10 : deallocate (chem_id, net_iso)
212 :
213 : end subroutine test_T_tau_varying
214 :
215 : !****
216 :
217 4 : subroutine test_T_tau_uniform(T_tau_opacity, tau_base_in)
218 :
219 : character(*), intent(in) :: T_tau_opacity
220 : real(dp), intent(in) :: tau_base_in
221 :
222 4 : real(dp) :: tau_base
223 4 : real(dp) :: errtol
224 :
225 : include 'formats'
226 :
227 4 : if (test_verbosely) then
228 2 : write (*, *)
229 2 : write (*, *) 'test_T_tau_uniform: ', TRIM(T_tau_opacity)
230 : end if
231 :
232 4 : tau_base = tau_base_in
233 4 : if (tau_base <= 0._dp) then
234 0 : tau_base = 2._dp/3._dp
235 : end if
236 :
237 4 : Pextra_factor = 1._dp
238 :
239 2 : select case (T_tau_opacity)
240 :
241 : case ('fixed')
242 :
243 2 : M = 1.9892000000000002D+32
244 2 : R = 6.3556231577545586D+10
245 2 : L = 2.4015399190199118D+32
246 2 : Teff = atm_Teff(L, R)
247 :
248 2 : kap_guess = 5.8850802481174469D-02
249 :
250 2 : max_iters = 0
251 :
252 : case ('iterated')
253 :
254 2 : logTeff = log10(5776._dp)
255 2 : log_gsurf = log10(cgrav*Msun/(Rsun*Rsun))
256 2 : log_Rsurf = log10(Rsun)
257 2 : log_M = log10(Msun)
258 :
259 2 : R = exp10(log_Rsurf)
260 2 : M = exp10(log_M)
261 2 : Teff = exp10(logTeff)
262 2 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
263 :
264 2 : kap_guess = 0.5d0
265 :
266 2 : max_iters = 30
267 :
268 : case default
269 :
270 0 : write (*, *) 'unknown value for T_tau_opacity '//trim(T_tau_opacity)
271 :
272 4 : stop
273 :
274 : end select
275 :
276 4 : Z = 0.02d0
277 4 : X = 0.70d0
278 4 : XC = 3.2724592105263235d-03
279 4 : XN = 9.5023842105263292d-04
280 4 : XO = 8.8218000000000601d-03
281 4 : call set_composition()
282 :
283 4 : errtol = 1.E-6_dp
284 :
285 : call atm_eval_T_tau_uniform( &
286 : tau_base, L, R, M, cgrav, kap_guess, Pextra_factor, &
287 : ATM_T_TAU_EDDINGTON, eos_proc, kap_proc, errtol, max_iters, SKIP_PARTIALS, &
288 : Teff, kap, &
289 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
290 : lnP, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
291 4 : ierr)
292 4 : if (ierr /= 0) then
293 0 : if (test_verbosely) write (*, *) 'failed in atm_eval_T_tau_uniform'
294 0 : call mesa_error(__FILE__, __LINE__)
295 : end if
296 4 : if (test_verbosely) write (*, 1) 'atm_get_grey logP surf', lnP/ln10
297 4 : if (test_verbosely) write (*, 1) 'atm_get_grey logT surf', lnT/ln10
298 4 : if (test_verbosely) write (*, 1) 'atm_get_grey tau_phot', tau_base
299 4 : if (test_verbosely) write (*, *)
300 :
301 4 : deallocate (chem_id, net_iso)
302 :
303 4 : end subroutine test_T_tau_uniform
304 :
305 : !****
306 :
307 6 : subroutine test_irradiated()
308 :
309 2 : real(dp) :: errtol
310 : type(Kap_General_Info), pointer :: rq
311 :
312 : include 'formats'
313 :
314 2 : if (test_verbosely) then
315 1 : write (*, *)
316 1 : write (*, *) 'test_irradiated'
317 : end if
318 :
319 2 : ierr = 0
320 :
321 2 : X = 0.70d0
322 2 : Z = 1d-2
323 2 : XC = 3.2724592105263235d-03
324 2 : XN = 9.5023842105263292d-04
325 2 : XO = 8.8218000000000601d-03
326 :
327 2 : call set_composition()
328 :
329 : ! at these conditions, the appropriate lowT opacity table is Freedman11
330 : ! this is around logR = 3.5, way off the default tables (max logR = 1)
331 2 : call kap_ptr(kap_handle, rq, ierr)
332 2 : if (ierr /= 0) return
333 2 : rq%kap_lowT_option = kap_lowT_Freedman11
334 :
335 : ! must set up tables again after changing options
336 2 : call kap_setup_tables(kap_handle, ierr)
337 :
338 2 : errtol = 1.E-6_dp
339 2 : max_iters = 30
340 :
341 2 : T_eq = 1000.
342 2 : kap_v = 4.E-3_dp
343 2 : kap_guess = 1.5d-2
344 2 : gamma = 0._dp
345 2 : P = 1d6
346 2 : M = 1.5d0*M_jupiter
347 : tau = 10 ! just a guess for use in getting R
348 2 : R = sqrt(cgrav*M*tau/(P*kap_guess)) ! g = P*kap/tau = G*M/R^2
349 2 : T_int = 900
350 2 : L = pi*crad*clight*R*R*T_int*T_int*T_int*T_int
351 :
352 : call atm_eval_irradiated( &
353 : L, R, M, cgrav, T_eq, P, kap_guess, kap_v, gamma, &
354 : eos_proc, kap_proc, errtol, max_iters, SKIP_PARTIALS, &
355 : Teff, kap, tau, &
356 : lnT, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
357 2 : ierr)
358 2 : if (ierr /= 0) then
359 0 : if (test_verbosely) write (*, *) 'bad return from atm_eval_irradiated'
360 0 : call mesa_error(__FILE__, __LINE__)
361 : end if
362 :
363 2 : if (test_verbosely) write (*, *)
364 2 : if (test_verbosely) write (*, 1) 'test_grey_irradiated: kap', kap
365 2 : if (test_verbosely) write (*, 1) 'M/M_jupiter', M/M_jupiter
366 2 : if (test_verbosely) write (*, 1) 'R/R_jupiter', R/R_jupiter
367 2 : if (test_verbosely) write (*, 1) 'R/Rsun', R/Rsun
368 2 : if (test_verbosely) write (*, 1) 'kap_v', kap_v
369 2 : if (test_verbosely) write (*, 1) 'P', P
370 2 : if (test_verbosely) write (*, *)
371 :
372 2 : if (test_verbosely) write (*, 1) 'tau', tau
373 2 : if (test_verbosely) write (*, 1) 'Teff', Teff
374 2 : if (test_verbosely) write (*, 1) 'T_eq', T_eq
375 2 : if (test_verbosely) write (*, 1) 'T_int', T_int
376 2 : if (test_verbosely) write (*, 1) 'T', exp(lnT)
377 2 : if (test_verbosely) write (*, 1) 'logT', lnT/ln10
378 2 : if (test_verbosely) write (*, *)
379 :
380 2 : deallocate (chem_id, net_iso)
381 :
382 : end subroutine test_irradiated
383 :
384 : !****
385 :
386 16 : subroutine set_composition()
387 16 : real(dp) :: xz, z2bar, ye, mass_correction, sumx, &
388 384 : dabar_dx(species), dzbar_dx(species), dmc_dx(species)
389 : integer :: i
390 : type(Kap_General_Info), pointer :: rq
391 16 : call kap_ptr(kap_handle, rq, ierr)
392 16 : rq%Zbase = Z
393 16 : Y = 1 - (X + Z)
394 16 : allocate (chem_id(species), net_iso(num_chem_isos), stat=ierr)
395 16 : if (ierr /= 0) then
396 0 : write (*, *) 'allocate failed'
397 0 : call mesa_error(__FILE__, __LINE__)
398 : end if
399 128 : chem_id(:) = [ih1, ihe4, ic12, in14, io16, ine20, img24]
400 125712 : net_iso(:) = 0
401 128 : forall (i=1:species) net_iso(chem_id(i)) = i
402 128 : xa(:) = [X, Y, xc, xn, xo, 0d0, 0d0]
403 128 : xa(species) = 1 - sum(xa(:))
404 : call composition_info( &
405 : species, chem_id, xa, X, Y, xz, abar, zbar, z2bar, z53bar, ye, &
406 16 : mass_correction, sumx, dabar_dx, dzbar_dx, dmc_dx)
407 16 : end subroutine set_composition
408 :
409 : !****
410 :
411 6216 : subroutine eos_proc( &
412 : lnP, lnT, &
413 6216 : lnRho, res, dres_dlnRho, dres_dlnT, &
414 : ierr)
415 :
416 : use eos_def, only: num_eos_basic_results
417 : use eos_lib, only: eosPT_get, radiation_pressure
418 :
419 : real(dp), intent(in) :: lnP
420 : real(dp), intent(in) :: lnT
421 : real(dp), intent(out) :: lnRho
422 : real(dp), intent(out) :: res(:)
423 : real(dp), intent(out) :: dres_dlnRho(:)
424 : real(dp), intent(out) :: dres_dlnT(:)
425 : integer, intent(out) :: ierr
426 :
427 6216 : real(dp) :: T, P, Prad, Pgas, logPgas, rho
428 6216 : real(dp) :: logRho, dlnRho_dlnPgas, dlnRho_dlnT
429 136752 : real(dp), dimension(num_eos_d_dxa_results, species) :: dres_dxa
430 :
431 6216 : T = exp(lnT)
432 6216 : P = exp(lnP)
433 :
434 6216 : Prad = radiation_pressure(T)
435 6216 : Pgas = max(1E-99_dp, P - Prad)
436 6216 : logPgas = log10(Pgas)
437 :
438 : call eosPT_get( &
439 : eos_handle, &
440 : species, chem_id, net_iso, xa, &
441 : Pgas, logPgas, T, lnT/ln10, &
442 : Rho, logRho, dlnRho_dlnPgas, dlnRho_dlnT, &
443 6216 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
444 :
445 6216 : lnRho = logRho*ln10
446 :
447 6216 : end subroutine eos_proc
448 :
449 : !****
450 :
451 6216 : subroutine kap_proc( &
452 6216 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
453 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
454 : ierr)
455 :
456 6216 : use kap_def, only: num_kap_fracs
457 : use kap_lib, only: kap_get
458 : use eos_def, only: i_lnfree_e, i_eta
459 :
460 : real(dp), intent(in) :: lnRho
461 : real(dp), intent(in) :: lnT
462 : real(dp), intent(in) :: res(:)
463 : real(dp), intent(in) :: dres_dlnRho(:)
464 : real(dp), intent(in) :: dres_dlnT(:)
465 : real(dp), intent(out) :: kap
466 : real(dp), intent(out) :: dlnkap_dlnRho
467 : real(dp), intent(out) :: dlnkap_dlnT
468 : integer, intent(out) :: ierr
469 :
470 80808 : real(dp) :: kap_fracs(num_kap_fracs), dlnkap_dxa(species)
471 :
472 : call kap_get( &
473 : kap_handle, species, chem_id, net_iso, xa, &
474 6216 : lnRho/ln10, lnT/ln10, res(i_lnfree_e), dres_dlnRho(i_lnfree_e), dres_dlnT(i_lnfree_e), &
475 6216 : res(i_eta), dres_dlnRho(i_eta), dres_dlnT(i_eta), &
476 12432 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, dlnkap_dxa, ierr)
477 :
478 6216 : end subroutine kap_proc
479 :
480 : end module test_atm_support
|
