Line data Source code
1 : module test_kap_support
2 :
3 : use chem_def
4 : use chem_lib
5 : use eos_def
6 : use eos_lib
7 : use kap_def
8 : use kap_lib
9 : use const_def, only: dp, ln10, arg_not_provided
10 : use math_lib
11 : use utils_lib, only: mesa_error
12 :
13 : implicit none
14 :
15 : logical, parameter :: use_shared_data_dir = .true. ! if false, then test using local version data
16 : !logical, parameter :: use_shared_data_dir = .false.
17 : logical, parameter :: use_cache = .true.
18 : logical, parameter :: show_info = .false.
19 :
20 : character(len=32) :: my_mesa_dir
21 : integer, parameter :: ionmax = 8
22 :
23 : real(dp) :: abar, zbar, z2bar, z53bar, ye, mass_correction, sumx
24 : integer, parameter :: species = 8
25 : integer, parameter :: h1 = 1, he4 = 2, c12 = 3, n14 = 4, o16 = 5, ne20 = 6, mg24 = 7, fe56 = 8
26 : integer, pointer, dimension(:) :: net_iso, chem_id
27 : real(dp) :: xa(species)
28 :
29 : real(dp) :: X, Y, Z, Zbase
30 :
31 : contains
32 :
33 4 : subroutine Do_One(quietly)
34 :
35 : logical, intent(in) :: quietly
36 : integer :: ierr
37 :
38 2 : call setup(quietly)
39 :
40 2 : allocate (net_iso(num_chem_isos), chem_id(species))
41 :
42 15714 : net_iso(:) = 0
43 :
44 2 : chem_id(h1) = ih1; net_iso(ih1) = h1
45 2 : chem_id(he4) = ihe4; net_iso(ihe4) = he4
46 2 : chem_id(c12) = ic12; net_iso(ic12) = c12
47 2 : chem_id(n14) = in14; net_iso(in14) = n14
48 2 : chem_id(o16) = io16; net_iso(io16) = o16
49 2 : chem_id(ne20) = ine20; net_iso(ine20) = ne20
50 2 : chem_id(mg24) = img24; net_iso(img24) = mg24
51 2 : chem_id(fe56) = ife56; net_iso(ife56) = fe56
52 :
53 2 : call test1(quietly, 1, 'fixed metals', ierr)
54 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
55 :
56 2 : call test1(quietly, 2, 'C/O enhanced', ierr)
57 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
58 :
59 : !call test1(quietly, 3, 'op_mono', ierr)
60 : !if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
61 :
62 2 : call test1(quietly, 4, 'AESOPUS', ierr)
63 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
64 :
65 2 : call kap_shutdown
66 :
67 2 : end subroutine Do_One
68 :
69 0 : subroutine test1_op_mono(quietly, test_str)
70 : logical, intent(in) :: quietly
71 : character(len=*), intent(in) :: test_str
72 : real(dp) :: &
73 0 : zbar, Z, xh, XC, XN, XO, XNe, kap1, &
74 0 : xmass(ionmax), &
75 0 : frac_Type2, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
76 0 : logT, logRho, kap, log10kap, dlnkap_dlnRho, dlnkap_dlnT
77 :
78 : logical :: CO_enhanced
79 : logical, parameter :: dbg = .false.
80 : integer :: ierr
81 0 : real(dp) :: chem_factors(ionmax)
82 :
83 : include 'formats'
84 :
85 0 : ierr = 0
86 :
87 0 : call setup_op_mono(ierr)
88 0 : if (ierr /= 0) then
89 0 : write (*, *) 'failed in setup_op_mono'
90 0 : return
91 : end if
92 :
93 0 : lnfree_e = 0; d_lnfree_e_dlnRho = 0; d_lnfree_e_dlnT = 0
94 0 : xc = 0d0
95 0 : xn = 0d0
96 0 : xo = 0d0
97 0 : xne = 0d0
98 :
99 0 : logT = 5.4d0
100 0 : logRho = -5.7d0
101 0 : z = 0.02d0
102 0 : xh = 0.65d0
103 :
104 : !call get_composition_info(Z, xh, abar, zbar, chem_id, xmass)
105 :
106 0 : chem_factors(:) = 1d0 ! scale factors for element opacity
107 :
108 0 : frac_Type2 = 0d0
109 0 : call test_op_mono(0, ierr)
110 0 : if (ierr /= 0) return
111 :
112 0 : log10kap = safe_log10(kap)
113 :
114 0 : if (.not. quietly) then
115 0 : write (*, *) trim(test_str)
116 0 : write (*, '(A)')
117 0 : call show_args
118 0 : call show_results
119 : end if
120 :
121 : ! test element factors with pure Fe56
122 0 : write (*, 1) 'pure fe56; factors all 1.0'
123 : !call get_pure_fe56_composition_info(abar, zbar, chem_id, xmass, fe56)
124 0 : call test_op_mono(fe56, ierr)
125 0 : if (ierr /= 0) return
126 0 : write (*, '(A)')
127 0 : kap1 = kap
128 :
129 0 : chem_factors(fe56) = 1.75d0
130 0 : write (*, 1) 'pure fe56; fe56 factor increased', chem_factors(fe56)
131 0 : call test_op_mono(fe56, ierr)
132 0 : if (ierr /= 0) return
133 0 : write (*, '(A)')
134 0 : write (*, 1) 'new/old', kap/kap1, kap, kap1
135 0 : write (*, '(A)')
136 :
137 : contains
138 :
139 0 : subroutine setup_op_mono(ierr)
140 : integer, intent(out) :: ierr
141 : character(len=256) :: op_mono_data_path, op_mono_data_cache_filename
142 :
143 0 : ierr = 0
144 :
145 : call GET_ENVIRONMENT_VARIABLE( &
146 0 : "MESA_OP_MONO_DATA_PATH", op_mono_data_path, status=ierr, trim_name=.true.)
147 0 : if (ierr /= 0) then
148 0 : write (*, *) 'failed to get environment variable MESA_OP_MONO_DATA_PATH'
149 0 : return
150 : end if
151 : call GET_ENVIRONMENT_VARIABLE( &
152 : "MESA_OP_MONO_DATA_CACHE_FILENAME", op_mono_data_cache_filename, &
153 0 : status=ierr, trim_name=.true.)
154 0 : if (ierr /= 0) then
155 0 : write (*, *) 'failed to get environment variable MESA_OP_MONO_DATA_CACHE_FILENAME'
156 0 : return
157 : end if
158 :
159 0 : write (*, *) 'call load_op_mono_data'
160 : call load_op_mono_data( &
161 0 : op_mono_data_path, op_mono_data_cache_filename, ierr)
162 0 : if (ierr /= 0) then
163 0 : write (*, *) 'failed in load_op_mono_data'
164 0 : write (*, *) 'my_mesa_dir '//trim(my_mesa_dir)
165 0 : write (*, *) 'op_mono_data_path '//trim(op_mono_data_path)
166 0 : write (*, *) 'op_mono_data_cache_filename '//trim(op_mono_data_cache_filename)
167 0 : return
168 : end if
169 :
170 : end subroutine setup_op_mono
171 :
172 0 : subroutine test_op_mono(fe56, ierr)
173 : use const_def, only: Lsun, Rsun, pi
174 : integer, intent(in) :: fe56
175 : integer, intent(out) :: ierr
176 :
177 : real, pointer :: &
178 0 : umesh(:), semesh(:), ff(:, :, :, :), ta(:, :, :, :), rs(:, :, :)
179 : integer :: kk, nel, nptot, ipe, nrad, iz(ionmax), iZ_rad(ionmax)
180 0 : real(dp), dimension(ionmax) :: fap, fac, &
181 0 : lgrad
182 0 : real(dp) :: flux, L, r
183 : logical, parameter :: screening = .true.
184 : !logical, parameter :: screening = .false.
185 :
186 : include 'formats'
187 :
188 0 : ierr = 0
189 :
190 : !write(*,*) 'call get_op_mono_params'
191 0 : call get_op_mono_params(nptot, ipe, nrad)
192 : allocate ( &
193 : umesh(nptot), semesh(nptot), ff(nptot, ipe, 4, 4), &
194 : ta(nptot, nrad, 4, 4), &
195 0 : rs(nptot, 4, 4), stat=ierr)
196 0 : if (ierr /= 0) return
197 :
198 : !write(*,*) 'call get_op_mono_args'
199 : call get_op_mono_args( &
200 : ionmax, xmass, 0d0, chem_id, chem_factors, &
201 0 : nel, iz, fap, fac, ierr)
202 0 : if (ierr /= 0) then
203 0 : write (*, *) 'error in get_op_mono_args, ierr = ', ierr
204 0 : return
205 : end if
206 :
207 0 : L = 11d0*Lsun
208 0 : r = 0.15d0*Rsun
209 0 : flux = L/(4*pi*r*r)
210 0 : if (fe56 > 0) then
211 0 : kk = 1
212 0 : iZ_rad(1) = iZ(fe56)
213 : else
214 0 : kk = ionmax
215 0 : iZ_rad(:) = iZ(:)
216 : end if
217 :
218 : call op_mono_get_radacc( &
219 : ! input
220 : kk, iZ_rad, ionmax, iZ, fap, fac, &
221 : flux, logT, logRho, screening, &
222 : ! output
223 : log10kap, &
224 : lgrad, &
225 : ! work arrays
226 : umesh, semesh, ff, ta, rs, &
227 0 : ierr)
228 :
229 0 : deallocate (umesh, semesh, ff, rs, ta)
230 0 : if (ierr /= 0) then
231 0 : write (*, *) 'error in op_mono_get_radacc, ierr = ', ierr
232 0 : return
233 : end if
234 :
235 0 : kap = exp10(log10kap)
236 :
237 0 : if (fe56 > 0) then
238 0 : write (*, '(A)')
239 0 : write (*, 1) 'grad', exp10(lgrad(kk))
240 0 : write (*, 1) 'lgrad', lgrad(kk)
241 : end if
242 :
243 0 : write (*, '(A)')
244 0 : write (*, 1) 'kap', kap
245 0 : write (*, 1) 'log10kap', log10kap
246 0 : write (*, 1) 'dlnkap_dlnRho', dlnkap_dlnRho
247 0 : write (*, 1) 'dlnkap_dlnT', dlnkap_dlnT
248 0 : write (*, 1) 'sum fap', sum(fap(:))
249 0 : write (*, '(A)')
250 :
251 0 : end subroutine test_op_mono
252 :
253 0 : subroutine show_args
254 : 1 format(a40, 1pe26.16)
255 0 : write (*, *) 'CO_enhanced', CO_enhanced
256 0 : write (*, 1) 'logT', logT
257 0 : write (*, 1) 'logRho', logRho
258 0 : write (*, 1) 'Z', Z
259 0 : write (*, 1) 'Zbase', Zbase
260 0 : write (*, 1) 'zbar', zbar
261 0 : write (*, 1) 'xh', xh
262 0 : write (*, 1) 'xc', xc
263 0 : write (*, 1) 'xn', xn
264 0 : write (*, 1) 'xo', xo
265 0 : write (*, 1) 'xne', xne
266 0 : write (*, 1) 'lnfree_e', lnfree_e
267 0 : write (*, '(A)')
268 0 : end subroutine show_args
269 :
270 0 : subroutine show_results
271 : use utils_lib
272 : 1 format(a40, 1pe26.16)
273 0 : write (*, 1) 'log10kap', log10kap
274 0 : write (*, 1) 'dlnkap_dlnRho', dlnkap_dlnRho
275 0 : write (*, 1) 'dlnkap_dlnT', dlnkap_dlnT
276 0 : write (*, '(A)')
277 0 : write (*, 1) 'kap', kap
278 0 : write (*, 1) 'dkap_dlnd', dlnkap_dlnRho*kap
279 0 : write (*, 1) 'dkap_dlnT', dlnkap_dlnT*kap
280 0 : write (*, '(A)')
281 0 : write (*, 1) 'frac_Type2', frac_Type2
282 0 : write (*, '(A)')
283 0 : if (is_bad(log10kap)) then
284 0 : write (*, *) 'bad log10kap'
285 : end if
286 0 : end subroutine show_results
287 :
288 : end subroutine test1_op_mono
289 :
290 12 : subroutine test1(quietly, which, test_str, ierr)
291 : use kap_def, only: Kap_General_Info, num_kap_fracs, i_frac_Type2
292 : logical, intent(in) :: quietly
293 : integer, intent(in) :: which
294 : character(len=*), intent(in) :: test_str
295 : integer, intent(out) :: ierr
296 : real(dp) :: &
297 6 : XC, XN, XO, XNe, &
298 6 : fC, fN, fO, fNe, dXC, dXO, &
299 12 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
300 12 : eta, d_eta_dlnRho, d_eta_dlnT, &
301 18 : logT, logRho, kap, log10kap, dlnkap_dlnRho, dlnkap_dlnT
302 78 : real(dp) :: kap_fracs(num_kap_fracs), dlnkap_dxa(species)
303 :
304 : ! eos results
305 480 : real(dp), dimension(num_eos_basic_results) :: res, deos_dlnd, deos_dlnT
306 150 : real(dp), dimension(num_eos_d_dxa_results, species) :: deos_dxa
307 :
308 : character(len=64) :: inlist
309 : integer :: eos_handle, kap_handle
310 : type(Kap_General_Info), pointer :: rq
311 :
312 : logical :: CO_enhanced
313 : logical, parameter :: dbg = .false.
314 : include 'formats'
315 :
316 6 : ierr = 0
317 :
318 6 : xa = 0
319 6 : X = 0; Z = 0; xc = 0; xn = 0; xo = 0; xne = 0
320 :
321 0 : select case (which)
322 : case (0) ! special test
323 :
324 : case (1) ! fixed
325 :
326 2 : inlist = 'inlist_test_fixed'
327 :
328 2 : CO_enhanced = .false.
329 2 : logT = 6d0
330 2 : logRho = -6d0
331 2 : X = 0.7d0
332 2 : Z = 0.018d0
333 :
334 2 : xa(h1) = X
335 2 : xa(he4) = 1d0 - X - Z
336 2 : xa(fe56) = Z
337 :
338 : case (2) ! co
339 :
340 2 : inlist = 'inlist_test_co'
341 :
342 2 : CO_enhanced = .true.
343 2 : logT = 6d0
344 2 : logRho = -6d0
345 2 : Zbase = 0.018d0
346 2 : dXC = 0.021d0
347 2 : dXO = 0.019d0
348 2 : X = 0.0d0
349 2 : fC = 0.173312d0
350 2 : fN = 0.053152d0
351 2 : fO = 0.482398d0
352 2 : fNe = 0.098668d0
353 2 : Z = Zbase + dXC + dXO
354 2 : xc = dXC + fC*Zbase
355 2 : xn = fN*Zbase
356 2 : xo = dXO + fO*Zbase
357 2 : xne = fNe*Zbase
358 :
359 2 : xa(h1) = X
360 2 : xa(he4) = 1d0 - X - Z
361 2 : xa(c12) = xc
362 2 : xa(n14) = xn
363 2 : xa(o16) = xo
364 2 : xa(ne20) = xne
365 :
366 : case (3) ! OP
367 :
368 0 : call mesa_error(__FILE__, __LINE__)
369 :
370 : case (4) ! AESOPUS
371 :
372 2 : inlist = 'inlist_aesopus'
373 2 : CO_enhanced = .true.
374 :
375 : ! conditions from RG surface
376 2 : logT = 3.5571504546260235D+000
377 2 : logRho = -8.2496430699014667D+000
378 2 : Z = 2.0074120713487353D-002
379 2 : X = 6.7888662523180188D-001
380 2 : xc = 2.9968458709806432D-003
381 2 : xn = 1.5270900145630591D-003
382 2 : xo = 9.3376263240514384D-003
383 :
384 2 : xa(h1) = X
385 2 : xa(he4) = 1d0 - X - Z
386 2 : xa(c12) = xc
387 2 : xa(n14) = xn
388 6 : xa(o16) = xo
389 :
390 : end select
391 :
392 : call basic_composition_info( &
393 : species, chem_id, xa, X, Y, Z, abar, zbar, z2bar, z53bar, &
394 6 : ye, mass_correction, sumx)
395 :
396 6 : eos_handle = alloc_eos_handle_using_inlist(inlist, ierr)
397 6 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
398 :
399 : call eosDT_get( &
400 : eos_handle, species, chem_id, net_iso, xa, &
401 : exp10(logRho), logRho, exp10(logT), logT, &
402 6 : res, deos_dlnd, deos_dlnT, deos_dxa, ierr)
403 :
404 6 : lnfree_e = res(i_lnfree_e)
405 6 : d_lnfree_e_dlnRho = deos_dlnd(i_lnfree_e)
406 6 : d_lnfree_e_dlnT = deos_dlnT(i_lnfree_e)
407 :
408 6 : eta = res(i_eta)
409 6 : d_eta_dlnRho = deos_dlnd(i_eta)
410 6 : d_eta_dlnT = deos_dlnT(i_eta)
411 :
412 6 : kap_handle = alloc_kap_handle_using_inlist(inlist, ierr)
413 6 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
414 :
415 6 : call kap_ptr(kap_handle, rq, ierr)
416 :
417 : call kap_get( &
418 : kap_handle, species, chem_id, net_iso, xa, logRho, logT, &
419 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
420 : eta, d_eta_dlnRho, d_eta_dlnT, &
421 6 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, dlnkap_dxa, ierr)
422 :
423 6 : log10kap = safe_log10(kap)
424 :
425 12 : if (.not. quietly) then
426 3 : write (*, *) 'test number', which
427 3 : write (*, *) trim(test_str)
428 3 : write (*, '(A)')
429 3 : call show_args
430 3 : call show_results
431 : end if
432 :
433 : contains
434 :
435 3 : subroutine show_args
436 : 1 format(a40, 1pe26.16)
437 3 : write (*, *) 'CO_enhanced', CO_enhanced
438 3 : write (*, 1) 'logT', logT
439 3 : write (*, 1) 'logRho', logRho
440 3 : write (*, 1) 'Z', Z
441 3 : write (*, 1) 'Zbase', rq%Zbase
442 3 : write (*, 1) 'zbar', zbar
443 3 : write (*, 1) 'X', X
444 3 : write (*, 1) 'xc', xc
445 3 : write (*, 1) 'xn', xn
446 3 : write (*, 1) 'xo', xo
447 3 : write (*, 1) 'xne', xne
448 3 : write (*, 1) 'lnfree_e', lnfree_e
449 3 : write (*, '(A)')
450 3 : end subroutine show_args
451 :
452 3 : subroutine show_results
453 : use utils_lib
454 : 1 format(a40, 1pe26.16)
455 3 : write (*, 1) 'log10kap', log10kap
456 3 : write (*, 1) 'dlnkap_dlnRho', dlnkap_dlnRho
457 3 : write (*, 1) 'dlnkap_dlnT', dlnkap_dlnT
458 3 : write (*, '(A)')
459 3 : write (*, 1) 'kap', kap
460 3 : write (*, 1) 'dkap_dlnd', dlnkap_dlnRho*kap
461 3 : write (*, 1) 'dkap_dlnT', dlnkap_dlnT*kap
462 3 : write (*, '(A)')
463 3 : write (*, 1) 'frac_Type2', kap_fracs(i_frac_Type2)
464 3 : write (*, '(A)')
465 3 : if (is_bad(log10kap)) then
466 0 : write (*, *) 'bad log10kap'
467 : end if
468 3 : end subroutine show_results
469 :
470 : end subroutine test1
471 :
472 2 : subroutine setup(quietly)
473 : use chem_lib
474 : use const_lib, only: const_init
475 : logical, intent(in) :: quietly
476 :
477 : integer :: ierr
478 : logical, parameter :: use_cache = .true.
479 :
480 2 : my_mesa_dir = '../..'
481 :
482 2 : call const_init(my_mesa_dir, ierr)
483 2 : if (ierr /= 0) then
484 0 : write (*, *) 'const_init failed'
485 0 : call mesa_error(__FILE__, __LINE__)
486 : end if
487 :
488 2 : call math_init()
489 :
490 2 : call chem_init('isotopes.data', ierr)
491 2 : if (ierr /= 0) then
492 0 : write (*, *) 'chem_init failed'
493 0 : call mesa_error(__FILE__, __LINE__)
494 : end if
495 :
496 2 : call eos_init('', use_cache, ierr)
497 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
498 :
499 2 : call kap_init(use_cache, '', ierr)
500 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
501 :
502 2 : end subroutine setup
503 :
504 : end module test_kap_support
|
