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 3109 : 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 1 : test_verbosely = test_verbosely_in
44 1 : cgrav = cgrav_in
45 1 : eos_handle = eos_handle_in
46 1 : kap_handle = kap_handle_in
47 :
48 1 : call test_table(ATM_TABLE_TAU_1M1, 'tau=0.1')
49 1 : call test_table(ATM_TABLE_TAU_1, 'tau=1')
50 1 : call test_table(ATM_TABLE_TAU_10, 'tau=10')
51 1 : call test_table(ATM_TABLE_TAU_100, 'tau=100')
52 :
53 1 : call test_table(ATM_TABLE_WD_TAU_25, 'WD tau=25')
54 1 : call test_table(ATM_TABLE_DB_WD_TAU_25, 'DB WD tau=25')
55 1 : call test_table(ATM_TABLE_PHOTOSPHERE, 'photosphere')
56 :
57 1 : call test_T_tau_varying(ATM_T_TAU_EDDINGTON, 'Eddington', -1._dp)
58 1 : call test_T_tau_varying(ATM_T_TAU_KRISHNA_SWAMY, 'Krishna-Swamy', -1._dp)
59 1 : call test_T_tau_varying(ATM_T_TAU_SOLAR_HOPF, 'solar Hopf', -1._dp)
60 1 : call test_T_tau_varying(ATM_T_TAU_TRAMPEDACH_SOLAR, 'Trampedach solar', -1._dp)
61 1 : call test_T_tau_varying(ATM_T_TAU_EDDINGTON, 'Eddington', 100._dp)
62 :
63 1 : call test_T_tau_uniform('fixed', 100._dp)
64 1 : call test_T_tau_uniform('iterated', 150._dp)
65 :
66 1 : call test_irradiated()
67 :
68 1 : end subroutine do_test_atm
69 :
70 : !****
71 :
72 7 : subroutine test_table(table_id, label)
73 :
74 : integer, intent(in) :: table_id
75 : character(*), intent(in) :: label
76 :
77 : include 'formats'
78 :
79 7 : ierr = 0
80 :
81 7 : if (test_verbosely) then
82 7 : write (*, *)
83 7 : write (*, *) 'test_table: ', label
84 : end if
85 :
86 7 : if (table_id == ATM_TABLE_WD_TAU_25) then
87 1 : logg = 7.1d0
88 1 : Teff = 5000
89 1 : M = 0.8d0*Msun
90 1 : R = sqrt(cgrav*M/exp10(logg))
91 : !write(*,*) 'R/Rsun', R/Rsun
92 1 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
93 : !write(*,*) 'L/Lsun', L/Lsun
94 6 : elseif (table_id == ATM_TABLE_DB_WD_TAU_25) then
95 1 : logg = 8.0d0
96 1 : Teff = 25000
97 1 : M = 0.8d0*Msun
98 1 : R = sqrt(cgrav*M/exp10(logg))
99 : !write(*,*) 'R/Rsun', R/Rsun
100 1 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
101 : !write(*,*) 'L/Lsun', L/Lsun else
102 : else
103 5 : M = Msun
104 5 : R = Rsun
105 5 : L = Lsun
106 5 : Teff = atm_Teff(L, R)
107 : end if
108 :
109 7 : 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 7 : ierr)
117 7 : 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 7 : T = exp(lnT)
123 7 : P = exp(lnP)
124 :
125 7 : Prad = radiation_pressure(T)
126 :
127 7 : if (test_verbosely) write (*, *)
128 7 : if (test_verbosely) write (*, '(99a16)') &
129 7 : 'T', 'log_T', 'log P', 'M/Msun', 'L/Lsun', 'R/Rsun', 'logPgas'
130 7 : 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 7 : if (test_verbosely) write (*, *)
133 :
134 7 : end subroutine test_table
135 :
136 : !****
137 :
138 5 : 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 : real(dp) :: errtol
145 : integer :: max_steps
146 :
147 : include 'formats'
148 :
149 5 : ierr = 0
150 :
151 5 : if (tau_base_in < 0._dp) then
152 4 : call atm_get_T_tau_base(T_tau_id, tau_base, ierr)
153 4 : 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 1 : tau_base = tau_base_in
159 : end if
160 :
161 5 : if (test_verbosely) then
162 5 : write (*, *)
163 5 : write (*, *) 'test_T_tau_varying: ', label
164 : end if
165 :
166 : ! TEST SOLAR VALUES
167 5 : logTeff = log10(5776d0)
168 5 : log_gsurf = log10(cgrav*Msun/(Rsun*Rsun))
169 5 : log_Rsurf = log10(Rsun)
170 5 : log_M = log10(Msun)
171 :
172 5 : Z = 0.02d0
173 5 : X = 0.70d0
174 5 : XC = 3.2724592105263235d-03
175 5 : XN = 9.5023842105263292d-04
176 5 : XO = 8.8218000000000601d-03
177 5 : call set_composition()
178 :
179 5 : R = exp10(log_Rsurf)
180 5 : M = exp10(log_M)
181 :
182 5 : g = exp10(log_gsurf)
183 5 : Teff = exp10(logTeff)
184 5 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
185 :
186 5 : errtol = 1.E-9_dp
187 5 : 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 5 : ierr)
197 5 : 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 5 : T = exp(lnT)
203 5 : P = exp(lnP)
204 :
205 5 : if (test_verbosely) write (*, *)
206 5 : if (test_verbosely) write (*, '(99a16)') 'T', 'log T', 'log P', 'M/Msun', 'L/Lsun', 'X', 'Z'
207 5 : if (test_verbosely) write (*, '(i16,99f16.8)') &
208 5 : floor(0.5d0 + T), log10(T), log10(P), M/Msun, L/Lsun, X, Z
209 5 : if (test_verbosely) write (*, *)
210 :
211 5 : deallocate (chem_id, net_iso)
212 :
213 : end subroutine test_T_tau_varying
214 :
215 : !****
216 :
217 2 : 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 : real(dp) :: tau_base
223 : real(dp) :: errtol
224 :
225 : include 'formats'
226 :
227 2 : if (test_verbosely) then
228 2 : write (*, *)
229 2 : write (*, *) 'test_T_tau_uniform: ', TRIM(T_tau_opacity)
230 : end if
231 :
232 2 : tau_base = tau_base_in
233 2 : if (tau_base <= 0._dp) then
234 0 : tau_base = 2._dp/3._dp
235 : end if
236 :
237 2 : Pextra_factor = 1._dp
238 :
239 1 : select case (T_tau_opacity)
240 :
241 : case ('fixed')
242 :
243 1 : M = 1.9892000000000002D+32
244 1 : R = 6.3556231577545586D+10
245 1 : L = 2.4015399190199118D+32
246 1 : Teff = atm_Teff(L, R)
247 :
248 1 : kap_guess = 5.8850802481174469D-02
249 :
250 1 : max_iters = 0
251 :
252 : case ('iterated')
253 :
254 1 : logTeff = log10(5776._dp)
255 1 : log_gsurf = log10(cgrav*Msun/(Rsun*Rsun))
256 1 : log_Rsurf = log10(Rsun)
257 1 : log_M = log10(Msun)
258 :
259 1 : R = exp10(log_Rsurf)
260 1 : M = exp10(log_M)
261 1 : Teff = exp10(logTeff)
262 1 : L = pi*crad*clight*R*R*Teff*Teff*Teff*Teff
263 :
264 1 : kap_guess = 0.5d0
265 :
266 1 : max_iters = 30
267 :
268 : case default
269 :
270 0 : write (*, *) 'unknown value for T_tau_opacity '//trim(T_tau_opacity)
271 :
272 2 : stop
273 :
274 : end select
275 :
276 2 : Z = 0.02d0
277 2 : X = 0.70d0
278 2 : XC = 3.2724592105263235d-03
279 2 : XN = 9.5023842105263292d-04
280 2 : XO = 8.8218000000000601d-03
281 2 : call set_composition()
282 :
283 2 : 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 2 : ierr)
292 2 : 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 2 : if (test_verbosely) write (*, 1) 'atm_get_grey logP surf', lnP/ln10
297 2 : if (test_verbosely) write (*, 1) 'atm_get_grey logT surf', lnT/ln10
298 2 : if (test_verbosely) write (*, 1) 'atm_get_grey tau_phot', tau_base
299 2 : if (test_verbosely) write (*, *)
300 :
301 2 : deallocate (chem_id, net_iso)
302 :
303 2 : end subroutine test_T_tau_uniform
304 :
305 : !****
306 :
307 3 : subroutine test_irradiated()
308 :
309 : real(dp) :: errtol
310 : type(Kap_General_Info), pointer :: rq
311 :
312 : include 'formats'
313 :
314 1 : if (test_verbosely) then
315 1 : write (*, *)
316 1 : write (*, *) 'test_irradiated'
317 : end if
318 :
319 1 : ierr = 0
320 :
321 1 : X = 0.70d0
322 1 : Z = 1d-2
323 1 : XC = 3.2724592105263235d-03
324 1 : XN = 9.5023842105263292d-04
325 1 : XO = 8.8218000000000601d-03
326 :
327 1 : 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 1 : call kap_ptr(kap_handle, rq, ierr)
332 1 : if (ierr /= 0) return
333 1 : rq%kap_lowT_option = kap_lowT_Freedman11
334 :
335 : ! must set up tables again after changing options
336 1 : call kap_setup_tables(kap_handle, ierr)
337 :
338 1 : errtol = 1.E-6_dp
339 1 : max_iters = 30
340 :
341 1 : T_eq = 1000.
342 1 : kap_v = 4.E-3_dp
343 1 : kap_guess = 1.5d-2
344 1 : gamma = 0._dp
345 1 : P = 1d6
346 1 : M = 1.5d0*M_jupiter
347 : tau = 10 ! just a guess for use in getting R
348 1 : R = sqrt(cgrav*M*tau/(P*kap_guess)) ! g = P*kap/tau = G*M/R^2
349 1 : T_int = 900
350 1 : 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 1 : ierr)
358 1 : 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 1 : if (test_verbosely) write (*, *)
364 1 : if (test_verbosely) write (*, 1) 'test_grey_irradiated: kap', kap
365 1 : if (test_verbosely) write (*, 1) 'M/M_jupiter', M/M_jupiter
366 1 : if (test_verbosely) write (*, 1) 'R/R_jupiter', R/R_jupiter
367 1 : if (test_verbosely) write (*, 1) 'R/Rsun', R/Rsun
368 1 : if (test_verbosely) write (*, 1) 'kap_v', kap_v
369 1 : if (test_verbosely) write (*, 1) 'P', P
370 1 : if (test_verbosely) write (*, *)
371 :
372 1 : if (test_verbosely) write (*, 1) 'tau', tau
373 1 : if (test_verbosely) write (*, 1) 'Teff', Teff
374 1 : if (test_verbosely) write (*, 1) 'T_eq', T_eq
375 1 : if (test_verbosely) write (*, 1) 'T_int', T_int
376 1 : if (test_verbosely) write (*, 1) 'T', exp(lnT)
377 1 : if (test_verbosely) write (*, 1) 'logT', lnT/ln10
378 1 : if (test_verbosely) write (*, *)
379 :
380 1 : deallocate (chem_id, net_iso)
381 :
382 : end subroutine test_irradiated
383 :
384 : !****
385 :
386 8 : subroutine set_composition()
387 : real(dp) :: xz, z2bar, ye, mass_correction, sumx, &
388 : dabar_dx(species), dzbar_dx(species), dmc_dx(species)
389 : integer :: i
390 : type(Kap_General_Info), pointer :: rq
391 8 : call kap_ptr(kap_handle, rq, ierr)
392 8 : rq%Zbase = Z
393 8 : Y = 1 - (X + Z)
394 8 : allocate (chem_id(species), net_iso(num_chem_isos), stat=ierr)
395 8 : if (ierr /= 0) then
396 0 : write (*, *) 'allocate failed'
397 0 : call mesa_error(__FILE__, __LINE__)
398 : end if
399 64 : chem_id(:) = [ih1, ihe4, ic12, in14, io16, ine20, img24]
400 62856 : net_iso(:) = 0
401 64 : forall (i=1:species) net_iso(chem_id(i)) = i
402 64 : xa(:) = [X, Y, xc, xn, xo, 0d0, 0d0]
403 64 : xa(species) = 1 - sum(xa(:))
404 : call composition_info( &
405 : species, chem_id, xa, X, Y, xz, abar, zbar, z2bar, z53bar, ye, &
406 8 : mass_correction, sumx, dabar_dx, dzbar_dx, dmc_dx)
407 8 : end subroutine set_composition
408 :
409 : !****
410 :
411 3108 : subroutine eos_proc( &
412 : lnP, lnT, &
413 3108 : 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 : real(dp) :: T, P, Prad, Pgas, logPgas, rho
428 : real(dp) :: logRho, dlnRho_dlnPgas, dlnRho_dlnT
429 : real(dp), dimension(num_eos_d_dxa_results, species) :: dres_dxa
430 :
431 3108 : T = exp(lnT)
432 3108 : P = exp(lnP)
433 :
434 3108 : Prad = radiation_pressure(T)
435 3108 : Pgas = max(1E-99_dp, P - Prad)
436 3108 : 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 3108 : res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)
444 :
445 3108 : lnRho = logRho*ln10
446 :
447 3108 : end subroutine eos_proc
448 :
449 : !****
450 :
451 3108 : subroutine kap_proc( &
452 3108 : lnRho, lnT, res, dres_dlnRho, dres_dlnT, &
453 : kap, dlnkap_dlnRho, dlnkap_dlnT, &
454 : ierr)
455 :
456 3108 : 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 : 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 3108 : lnRho/ln10, lnT/ln10, res(i_lnfree_e), dres_dlnRho(i_lnfree_e), dres_dlnT(i_lnfree_e), &
475 3108 : res(i_eta), dres_dlnRho(i_eta), dres_dlnT(i_eta), &
476 6216 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, dlnkap_dxa, ierr)
477 :
478 3108 : end subroutine kap_proc
479 :
480 : end module test_atm_support
|
