Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module micro
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, ln10, crad, qe, avo, kerg, one_third, four_thirds_pi
24 : use star_utils, only: foreach_cell
25 : use utils_lib, only: is_bad
26 :
27 : implicit none
28 :
29 : private
30 : public :: set_micro_vars
31 : public :: set_eos_with_mask
32 : public :: do_eos_for_cell
33 : public :: store_eos_for_cell
34 : public :: do_kap_for_cell
35 : public :: shutdown_microphys
36 :
37 : logical, parameter :: dbg = .false.
38 : logical :: initiaze_kap_grid = .true.
39 : real(dp), public, save :: fk_pcg_old(17)
40 :
41 : contains
42 :
43 67 : subroutine set_micro_vars( &
44 : s, nzlo, nzhi, skip_eos, skip_net, skip_neu, skip_kap, ierr)
45 :
46 : use kap_support, only: prepare_kap
47 : use star_utils, only: start_time, update_time
48 : use net, only: do_net
49 : use chem_def, only: icno, ipp, chem_isos
50 : use rates_def, only: i_rate
51 : use kap_lib
52 :
53 :
54 : type (star_info), pointer :: s
55 : integer, intent(in) :: nzlo, nzhi
56 : logical, intent(in) :: skip_eos, skip_net, skip_neu, skip_kap
57 : integer, intent(out) :: ierr
58 :
59 : integer :: k, op_err, i
60 : integer(i8) :: time0
61 66 : real(dp) :: total, alfa, beta
62 : character(len=4) :: e_name
63 1188 : real(dp) :: fk(17), delta
64 :
65 :
66 :
67 : include 'formats'
68 :
69 66 : ierr = 0
70 : if (dbg) then
71 : write(*,'(A)')
72 : write(*,*) 'set_micro_vars'
73 : write(*,*) 'nzlo', nzlo
74 : write(*,*) 'nzhi', nzhi
75 : write(*,*) 'skip_net', skip_net
76 : write(*,*) 'skip_kap', skip_kap
77 : end if
78 :
79 66 : if (.not. skip_eos) then
80 66 : call set_eos(ierr)
81 66 : if (ierr /= 0) return
82 : end if
83 :
84 80060 : do k=nzlo, nzhi
85 80060 : if (s% csound_start(k) < 0) then
86 40733 : s% csound_start(k) = s% csound(k)
87 : end if
88 : end do
89 :
90 80060 : do k=nzlo, nzhi
91 80060 : if (k == 1) then
92 66 : s% rho_face(k) = s% rho(k)
93 66 : if (.not. s% u_flag) s% P_face_ad(k)%val = s% Peos(k)
94 66 : s% csound_face(1) = s% csound(1)
95 : else
96 79928 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
97 79928 : beta = 1 - alfa
98 79928 : s% rho_face(k) = alfa*s% rho(k) + beta*s% rho(k-1)
99 79928 : if (.not. s% u_flag) s% P_face_ad(k)%val = alfa*s% Peos(k) + beta*s% Peos(k-1)
100 79928 : s% csound_face(k) = alfa*s% csound(k) + beta*s% csound(k-1)
101 : end if
102 : end do
103 :
104 66 : if (.not. (skip_kap .and. skip_neu)) then
105 :
106 66 : if (.not. skip_kap .and. s% op_mono_method == 'hu') then
107 66 : call prepare_kap(s, ierr)
108 66 : if (ierr /= 0) return
109 : end if
110 66 : if (.not. skip_kap .and. s% op_mono_method == 'mombarg' .and. &
111 : s% high_logT_op_mono_full_on - s% low_logT_op_mono_full_on > 0 ) then
112 0 : fk = 0
113 0 : do i=1, s% species
114 0 : e_name = chem_isos% name(s% chem_id(i))
115 0 : if (e_name == 'h1') fk(1) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
116 0 : if (e_name == 'he4') fk(2) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
117 0 : if (e_name == 'c12') fk(3) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
118 0 : if (e_name == 'n14') fk(4) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
119 0 : if (e_name == 'o16') fk(5) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
120 0 : if (e_name == 'ne20')fk(6) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
121 0 : if (e_name == 'na23')fk(7) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
122 0 : if (e_name == 'mg24')fk(8) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
123 0 : if (e_name == 'al27')fk(9) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
124 0 : if (e_name == 'si28')fk(10) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
125 0 : if (e_name == 's32') fk(11) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
126 0 : if (e_name == 'ar40')fk(12) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
127 0 : if (e_name == 'ca40')fk(13) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
128 0 : if (e_name == 'cr52')fk(14) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
129 0 : if (e_name == 'mn55')fk(15) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
130 0 : if (e_name == 'fe56')fk(16) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
131 0 : if (e_name == 'ni58')fk(17) = s% xa(i,s% nz)/ chem_isos% W(s% chem_id(i))
132 : end do
133 0 : fk = fk / sum(fk)
134 0 : delta = MAXVAL(ABS(fk - fk_pcg_old)/fk_pcg_old)
135 :
136 0 : if (initiaze_kap_grid) then
137 0 : call call_load_op_master(s% emesh_data_for_op_mono_path, ierr)
138 0 : write(*,*) 'Computing kappa grid for initial mixture.'
139 0 : call call_compute_kappa_grid_mombarg(fk, ierr)
140 : !write(*,*) 'Finished computing grid for initial mixture.'
141 0 : fk_pcg_old = fk
142 0 : initiaze_kap_grid = .false.
143 0 : else if (delta > 1d-4) THEN
144 0 : write(*,*) 'Computing kappa grid for core mixture.'
145 0 : write(*,*) delta
146 0 : call call_compute_kappa_grid_mombarg(fk, ierr)
147 : !write(*,*) 'Finished computing grid for core mixture.'
148 0 : fk_pcg_old = fk
149 : end if
150 : end if
151 :
152 66 : if (s% use_other_opacity_factor) then
153 0 : call s% other_opacity_factor(s% id, ierr)
154 0 : if (ierr /= 0) return
155 : else
156 80060 : s% extra_opacity_factor(1:s% nz) = s% opacity_factor
157 : end if
158 :
159 66 : if (s% doing_timing) call start_time(s, time0, total)
160 :
161 :
162 :
163 66 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
164 : do k = nzlo, nzhi
165 : op_err = 0
166 : call do1_neu_kap(s,k,op_err)
167 : if (op_err /= 0) ierr = op_err
168 : end do
169 : !$OMP END PARALLEL DO
170 66 : if (s% doing_timing) call update_time(s, time0, total, s% time_neu_kap)
171 66 : if (ierr /= 0) then
172 0 : if (s% report_ierr) write(*,*) 'do1_neu_kap returned ierr', ierr
173 0 : return
174 : end if
175 :
176 : end if
177 :
178 66 : if (ierr /= 0) return
179 :
180 132 : if (.not. skip_net) then
181 :
182 : if (dbg) write(*,*) 'micro: call do_net'
183 66 : call do_net(s, nzlo, nzhi, ierr)
184 : if (dbg) write(*,*) 'micro: done do_net'
185 66 : if (ierr /= 0) then
186 0 : if (s% report_ierr) write(*,*) 'do_net returned ierr', ierr
187 0 : return
188 : end if
189 :
190 : end if
191 :
192 : contains
193 :
194 80192 : subroutine set_eos(ierr)
195 : integer, intent(out) :: ierr
196 : include 'formats'
197 66 : ierr = 0
198 : if (dbg) write(*,*) 'call do_eos'
199 66 : if (s% doing_timing) call start_time(s, time0, total)
200 : call do_eos(s,nzlo,nzhi,ierr)
201 66 : if (s% doing_timing) call update_time(s, time0, total, s% time_eos)
202 66 : if (ierr /= 0) then
203 0 : if (s% report_ierr) write(*,*) 'do_eos returned ierr', ierr
204 0 : return
205 : end if
206 66 : end subroutine set_eos
207 :
208 : subroutine debug(str)
209 : use chem_def
210 : character (len=*), intent(in) :: str
211 : integer :: k, j
212 : include 'formats'
213 : k = 1469
214 : do j=1,1 !s% species
215 : if (.true. .or. s% xa(j,k) > 1d-9) &
216 : write(*,1) trim(str) // ' xin(net_iso(i' // &
217 : trim(chem_isos% name(s% chem_id(j))) // '))= ', s% xa(j,k)
218 : end do
219 : end subroutine debug
220 :
221 79994 : subroutine do1_neu_kap(s,k,ierr)
222 : use neu_def,only:Tmin_neu
223 : use neu, only: do_neu_for_cell, do_clear_neu_for_cell
224 : type (star_info), pointer :: s
225 : integer, intent(in) :: k
226 : integer, intent(out) :: ierr
227 79994 : if (s% T(k) >= Tmin_neu .and. s% non_nuc_neu_factor > 0d0) then
228 26684 : call do_neu_for_cell(s,k,ierr)
229 : else
230 53310 : call do_clear_neu_for_cell(s,k,ierr)
231 : end if
232 79994 : if (ierr /= 0) return
233 79994 : if (skip_kap) return
234 79994 : call do_kap_for_cell(s,k,ierr)
235 79994 : if (ierr /= 0) return
236 79994 : end subroutine do1_neu_kap
237 :
238 : end subroutine set_micro_vars
239 :
240 :
241 66 : subroutine do_eos(s,nzlo,nzhi,ierr)
242 :
243 : type (star_info), pointer :: s
244 : integer, intent(in) :: nzlo, nzhi
245 : integer, intent(out) :: ierr
246 : logical, parameter :: use_omp = .true.
247 : include 'formats'
248 :
249 : ierr = 0
250 :
251 : if (dbg) write(*,*) 'do_eos call foreach_cell', nzlo, nzhi
252 :
253 66 : call foreach_cell(s,nzlo,nzhi,use_omp,do_eos_for_cell,ierr)
254 :
255 : end subroutine do_eos
256 :
257 :
258 0 : subroutine set_eos_with_mask (s,nzlo,nzhi,mask,ierr)
259 :
260 : type (star_info), pointer :: s
261 : integer, intent(in) :: nzlo, nzhi
262 : logical, intent(in) :: mask(:)
263 : integer, intent(out) :: ierr
264 : include 'formats'
265 : integer :: k
266 0 : logical :: lerr(s%nz)
267 :
268 0 : ierr = 0
269 :
270 : if (dbg) write(*,*) 'set_eos_with_mask'
271 :
272 0 : !$OMP PARALLEL DO PRIVATE(ierr) SCHEDULE(dynamic,2)
273 : do k = nzlo, nzhi
274 : if (mask(k)) then
275 : call do_eos_for_cell(s, k, ierr)
276 : lerr(k) = ierr /= 0
277 : else
278 : lerr(k) = .FALSE.
279 : end if
280 : end do
281 : !$OMP END PARALLEL DO
282 :
283 0 : if (ANY(lerr(nzlo:nzhi))) then
284 0 : ierr = 1
285 : else
286 0 : ierr = 0
287 : end if
288 :
289 0 : end subroutine set_eos_with_mask
290 :
291 :
292 79994 : subroutine do_eos_for_cell(s,k,ierr)
293 :
294 : use chem_def
295 : use chem_lib
296 : use eos_def
297 : use eos_support
298 : use net_def, only: net_general_info
299 : use star_utils, only: write_eos_call_info, lookup_nameofvar
300 :
301 : type (star_info), pointer :: s
302 : integer, intent(in) :: k
303 : integer, intent(out) :: ierr
304 :
305 6399520 : real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
306 1999850 : real(dp), dimension(num_eos_d_dxa_results,s% species) :: d_dxa
307 79994 : real(dp) :: sumx, logRho, logT, logPgas
308 79994 : integer, pointer :: net_iso(:)
309 : integer :: species
310 : real(dp), parameter :: epsder = 1d-4, Z_limit = 0.5d0
311 : real(dp), parameter :: LOGRHO_TOL = 1d-8, LOGPGAS_TOL = 1d-8
312 :
313 : logical, parameter :: testing = .false.
314 :
315 : include 'formats'
316 :
317 79994 : ierr = 0
318 :
319 79994 : net_iso => s% net_iso
320 79994 : species = s% species
321 : call basic_composition_info( &
322 : species, s% chem_id, s% xa(1:species,k), s% X(k), s% Y(k), s% Z(k), &
323 : s% abar(k), s% zbar(k), s% z2bar(k), s% z53bar(k), &
324 79994 : s% ye(k), s% mass_correction(k), sumx)
325 :
326 79994 : logT = s% lnT(k)/ln10
327 :
328 79994 : if (s%fix_Pgas) then
329 :
330 0 : logPgas = s% lnPgas(k)/ln10
331 :
332 : call solve_eos_given_PgasT( &
333 : s, k, s% xa(:,k), &
334 : logT, logPgas, s% lnd(k)/ln10, LOGRHO_TOL, LOGPGAS_TOL, &
335 0 : logRho, res, d_dlnd, d_dlnT, d_dxa, ierr)
336 0 : if (ierr /= 0) then
337 0 : if (s% report_ierr) then
338 0 : write(*,*) 'do_eos_for_cell: solve_eos_given_PT ierr', ierr
339 : end if
340 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
341 0 : return
342 : end if
343 :
344 0 : s% lnd(k) = logRho*ln10
345 0 : s% rho(k) = exp(s% lnd(k))
346 :
347 : end if
348 :
349 79994 : logRho = s% lnd(k)/ln10
350 :
351 79994 : if (s% solver_test_eos_partials) then
352 : eos_test_partials = (k == s% solver_test_partials_k .and. &
353 : s% solver_call_number == s% solver_test_partials_call_number .and. &
354 0 : s% solver_iter == s% solver_test_partials_iter_number )
355 : end if
356 :
357 : call get_eos( &
358 : s, k, s% xa(:,k), &
359 : s% rho(k), logRho, s% T(k), logT, &
360 : res, s% d_eos_dlnd(:,k), s% d_eos_dlnT(:,k), &
361 79994 : s% d_eos_dxa(:,:,k), ierr)
362 79994 : if (ierr /= 0) then
363 0 : if (s% report_ierr) then
364 0 : write(*, *) s% retry_message
365 0 : write(*,*) 'do_eos_for_cell: get_eos ierr', ierr
366 : end if
367 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
368 0 : return
369 : end if
370 :
371 79994 : if (s% solver_test_eos_partials .and. eos_test_partials) then
372 0 : s% solver_test_partials_val = eos_test_partials_val
373 0 : s% solver_test_partials_dval_dx = eos_test_partials_dval_dx
374 : end if
375 :
376 79994 : s% lnPgas(k) = res(i_lnPgas)
377 79994 : s% Pgas(k) = exp(s% lnPgas(k))
378 :
379 79994 : call store_stuff(ierr)
380 79994 : if (ierr /= 0) return
381 :
382 : if (s% solver_test_partials_write_eos_call_info .and. &
383 : k == s% solver_test_partials_k .and. &
384 79994 : s% solver_call_number == s% solver_test_partials_call_number .and. &
385 159988 : s% solver_iter == s% solver_test_partials_iter_number) then
386 0 : call write_eos_call_info(s,k)
387 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell: write_eos_call_info')
388 : end if
389 :
390 : contains
391 :
392 79994 : subroutine store_stuff(ierr)
393 :
394 : integer, intent(out) :: ierr
395 :
396 : include 'formats'
397 :
398 : call store_eos_for_cell(s, k, &
399 : res, s% d_eos_dlnd(:,k), s% d_eos_dlnT(:,k), &
400 79994 : s% d_eos_dxa(:,:,k), ierr)
401 :
402 79994 : if (ierr /= 0) then
403 0 : if (s% report_ierr) then
404 0 : call write_eos_call_info(s,k)
405 0 : write(*,2) 'store_eos_for_cell failed', k
406 : end if
407 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_eos_for_cell')
408 0 : return
409 : end if
410 :
411 79994 : end subroutine store_stuff
412 :
413 : end subroutine do_eos_for_cell
414 :
415 :
416 159988 : subroutine store_eos_for_cell(s, k, res, d_dlnd, d_dlnT, d_dxa, ierr)
417 :
418 : use eos_def
419 : use chem_def
420 : use star_utils, only: eval_csound, write_eos_call_info
421 :
422 : type (star_info), pointer :: s
423 : integer, intent(in) :: k
424 : real(dp), intent(in), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT
425 : real(dp), intent(in) :: d_dxa(num_eos_d_dxa_results,s% species)
426 : integer, intent(out) :: ierr
427 :
428 : integer :: i, j
429 79994 : real(dp) :: T, rho
430 :
431 : include 'formats'
432 :
433 79994 : ierr = 0
434 :
435 : ! check values
436 2159838 : do i = 1, num_eos_basic_results
437 2159838 : if (is_bad(res(i) + d_dlnd(i) + d_dlnT(i))) then
438 0 : ierr = -1
439 0 : if (s% report_ierr) then
440 0 : !$OMP critical (micro_crit0)
441 0 : write(*,2) trim(eosDT_result_names(i)), k, res(i)
442 0 : write(*,2) 'd_dlnd ' // trim(eosDT_result_names(i)), k, d_dlnd(i)
443 0 : write(*,2) 'd_dlnT ' // trim(eosDT_result_names(i)), k, d_dlnT(i)
444 0 : write(*,'(A)')
445 0 : call write_eos_call_info(s,k)
446 : !$OMP end critical (micro_crit0)
447 : end if
448 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'store_eos_for_cell')
449 0 : return
450 : end if
451 : end do
452 :
453 79994 : T = s% T(k)
454 79994 : rho = s% rho(k)
455 79994 : if (T > 1d15 .or. rho > 1d15) then
456 0 : ierr = -1
457 0 : if (s% report_ierr) then
458 0 : write(*,4) 'bad T or rho for eos', k, s% solver_iter, s% model_number
459 0 : write(*,2) 'T', k, T
460 0 : write(*,2) 'rho', k, rho
461 : end if
462 0 : return
463 : end if
464 79994 : s% Prad(k) = crad * T*T*T*T / 3
465 79994 : s% Peos(k) = s% Prad(k) + s% Pgas(k)
466 79994 : s% lnPeos(k) = log(s% Peos(k))
467 79994 : s% lnS(k) = res(i_lnS)
468 79994 : s% lnE(k) = res(i_lnE)
469 79994 : s% energy(k) = exp(s% lnE(k))
470 79994 : s% entropy(k) = exp(s% lnS(k))
471 79994 : s% grada(k) = res(i_grad_ad)
472 79994 : s% dE_dRho(k) = res(i_dE_drho)
473 79994 : s% Cv(k) = res(i_Cv)
474 79994 : s% Cp(k) = res(i_Cp)
475 79994 : s% chiRho(k) = res(i_chiRho)
476 79994 : s% chiT(k) = res(i_chiT)
477 79994 : s% gamma1(k) = res(i_gamma1)
478 79994 : s% gamma3(k) = res(i_gamma3)
479 79994 : s% eta(k) = res(i_eta)
480 : s% gam(k) = s% z53bar(k)*qe*qe * &
481 79994 : pow(four_thirds_pi*avo*rho*s% zbar(k)/s% abar(k),one_third) / (kerg*T)
482 79994 : s% mu(k) = res(i_mu)
483 79994 : s% lnfree_e(k) = res(i_lnfree_e)
484 :
485 79994 : s% phase(k) = res(i_phase)
486 79994 : s% latent_ddlnT(k) = res(i_latent_ddlnT)
487 79994 : s% latent_ddlnRho(k) = res(i_latent_ddlnRho)
488 :
489 79994 : s% eos_frac_OPAL_SCVH(k) = res(i_frac_OPAL_SCVH)
490 79994 : s% eos_frac_HELM(k) = res(i_frac_HELM)
491 79994 : s% eos_frac_Skye(k) = res(i_frac_Skye)
492 79994 : s% eos_frac_PC(k) = res(i_frac_PC)
493 79994 : s% eos_frac_FreeEOS(k) = res(i_frac_FreeEOS)
494 79994 : s% eos_frac_CMS(k) = res(i_frac_CMS)
495 79994 : s% eos_frac_ideal(k) = res(i_frac_ideal)
496 :
497 79994 : s% chiRho_for_partials(k) = s% Pgas(k)*d_dlnd(i_lnPgas)/s% Peos(k)
498 79994 : s% chiT_for_partials(k) = (s% Pgas(k)*d_dlnT(i_lnPgas) + 4d0*s% Prad(k))/s% Peos(k)
499 79994 : s% dE_drho_for_partials(k) = d_dlnd(i_lnE)*s% energy(k)/s% rho(k)
500 79994 : s% Cv_for_partials(k) = d_dlnT(i_lnE)*s% energy(k)/s% T(k)
501 79994 : s% dS_drho_for_partials(k) = d_dlnd(i_lnS)*s% entropy(k)/s% rho(k)
502 79994 : s% dS_dT_for_partials(k) = d_dlnT(i_lnS)*s% entropy(k)/s% T(k)
503 719946 : do j=1, s% species
504 639952 : s% dlnE_dxa_for_partials(j,k) = d_dxa(i_lnE,j)
505 719946 : s% dlnPeos_dxa_for_partials(j,k) = s% Pgas(k)*d_dxa(i_lnPgas,j)/s% Peos(k)
506 : end do
507 :
508 79994 : s% QQ(k) = s% chiT(k)/(s% rho(k)*s% T(k)*s% chiRho(k)) ! thermal expansion coefficient
509 79994 : s% csound(k) = eval_csound(s,k,ierr)
510 :
511 : ! check some key values
512 : if (s% gamma1(k) <= 0 .or. &
513 : s% grada(k) <= 0 .or. &
514 : s% csound(k) <= 0 .or. is_bad(s% csound(k)) .or. &
515 : s% chiT(k) <= 0 .or. &
516 : s% chiRho(k) <= 0 .or. &
517 79994 : s% Cv(k) <= 0 .or. &
518 : s% Cp(k) <= 0) then
519 0 : ierr = -1
520 : end if
521 :
522 79994 : if (ierr /= 0) then
523 0 : if (s% report_ierr) then
524 0 : !$OMP critical (micro_crit1)
525 0 : write(*,2) 's% Cp(k)', k, s% Cp(k)
526 0 : write(*,2) 's% csound(k)', k, s% csound(k)
527 0 : write(*,2) 's% lnPeos(k)', k, s% lnPeos(k)
528 0 : write(*,2) 's% gam(k)', k, s% gam(k)
529 0 : write(*,2) 's% Peos(k)', k, s% Peos(k)
530 0 : write(*,2) 's% Pgas(k)', k, s% Pgas(k)
531 0 : write(*,2) 's% rho(k)', k, s% rho(k)
532 0 : write(*,2) 's% T(k)', k, s% T(k)
533 0 : write(*,2) 'logRho', k, s% lnd(k)/ln10
534 0 : write(*,2) 'logT', k, s% lnT(k)/ln10
535 0 : write(*,2) 'abar', k, s% abar(k)
536 0 : write(*,2) 'zbar', k, s% zbar(k)
537 0 : write(*,*)
538 0 : call write_eos_call_info(s,k)
539 : !$OMP end critical (micro_crit1)
540 : end if
541 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'store_eos_for_cell')
542 0 : ierr = -1
543 : end if
544 :
545 79994 : end subroutine store_eos_for_cell
546 :
547 :
548 79994 : subroutine do_kap_for_cell(s,k,ierr)
549 :
550 79994 : use net_def,only:net_general_info
551 : use rates_def, only:i_rate
552 : use chem_def
553 : use kap_def
554 : use kap_lib
555 : use kap_support, only: fraction_of_op_mono, get_kap
556 : use eos_def, only : i_lnfree_e, i_eta
557 : use eos_lib, only: Radiation_Pressure
558 : use star_utils, only: get_XYZ
559 :
560 : type (star_info), pointer :: s
561 : integer, intent(in) :: k
562 : integer, intent(out) :: ierr
563 :
564 : real(dp) :: &
565 159988 : log10_rho, log10_T, dlnkap_dlnd, dlnkap_dlnT, zbar, &
566 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
567 79994 : eta, d_eta_dlnRho, d_eta_dlnT, opacity_factor
568 :
569 399970 : real(dp), dimension(num_kap_fracs) :: kap_fracs
570 :
571 79994 : real(dp), pointer :: xa(:)
572 : logical :: test_partials
573 :
574 : include 'formats'
575 :
576 79994 : ierr = 0
577 :
578 79994 : log10_rho = s% lnd(k)/ln10
579 79994 : log10_T = s% lnT(k)/ln10
580 :
581 79994 : lnfree_e = s% lnfree_e(k)
582 79994 : d_lnfree_e_dlnRho = s% d_eos_dlnd(i_lnfree_e,k)
583 79994 : d_lnfree_e_dlnT = s% d_eos_dlnT(i_lnfree_e,k)
584 :
585 79994 : eta = s% eta(k)
586 79994 : d_eta_dlnRho = s% d_eos_dlnd(i_eta,k)
587 79994 : d_eta_dlnT = s% d_eos_dlnT(i_eta,k)
588 :
589 79994 : if (s% use_starting_composition_for_kap) then
590 0 : xa(1:s% species) => s% xa_start(1:s% species,k)
591 0 : zbar = s% zbar_start(k)
592 : else
593 79994 : xa(1:s% species) => s% xa(1:s% species,k)
594 79994 : zbar = s% zbar(k)
595 : end if
596 :
597 : call get_kap( &
598 : s, k, zbar, xa, log10_rho, log10_T, &
599 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
600 : eta, d_eta_dlnRho, d_eta_dlnT, &
601 79994 : kap_fracs, s% opacity(k), dlnkap_dlnd, dlnkap_dlnT, ierr)
602 :
603 :
604 : ! unpack fractions
605 79994 : s% kap_frac_lowT(k) = kap_fracs(i_frac_lowT)
606 79994 : s% kap_frac_highT(k) = kap_fracs(i_frac_highT)
607 79994 : s% kap_frac_Type2(k) = kap_fracs(i_frac_Type2)
608 79994 : s% kap_frac_Compton(k) = kap_fracs(i_frac_Compton)
609 :
610 79994 : if (is_bad_num(s% opacity(k)) .or. ierr /= 0) then
611 0 : if (s% report_ierr) then
612 0 : write(*,*) 'do_kap_for_cell: get_kap ierr', ierr
613 0 : !$omp critical (star_kap_get)
614 0 : call show_stuff()
615 : !$omp end critical (star_kap_get)
616 : end if
617 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_kap_for_cell')
618 0 : ierr = -1
619 0 : return
620 : end if
621 :
622 79994 : opacity_factor = s% extra_opacity_factor(k)
623 79994 : if (s% min_logT_for_opacity_factor_off > 0) then
624 0 : if (log10_T >= s% max_logT_for_opacity_factor_off .or. &
625 : log10_T <= s% min_logT_for_opacity_factor_off) then
626 : opacity_factor = 1
627 0 : else if (log10_T > s% max_logT_for_opacity_factor_on) then
628 : opacity_factor = 1 + (opacity_factor-1)* &
629 : (log10_T - s% max_logT_for_opacity_factor_off)/ &
630 0 : (s% max_logT_for_opacity_factor_on - s% max_logT_for_opacity_factor_off)
631 0 : else if (log10_T < s% min_logT_for_opacity_factor_on) then
632 : opacity_factor = 1 + (opacity_factor-1)* &
633 : (log10_T - s% min_logT_for_opacity_factor_off)/ &
634 0 : (s% min_logT_for_opacity_factor_on - s% min_logT_for_opacity_factor_off)
635 : end if
636 : end if
637 :
638 79994 : s% opacity(k) = s% opacity(k)*opacity_factor
639 79994 : s% d_opacity_dlnd(k) = s% opacity(k)*dlnkap_dlnd
640 79994 : s% d_opacity_dlnT(k) = s% opacity(k)*dlnkap_dlnT
641 :
642 79994 : if (s% opacity(k) > s% opacity_max .and. s% opacity_max > 0) then
643 0 : s% opacity(k) = s% opacity_max
644 0 : s% d_opacity_dlnd(k) = 0
645 0 : s% d_opacity_dlnT(k) = 0
646 : end if
647 :
648 79994 : if (s% opacity(k) < s% opacity_min .and. s% opacity_min > 0) then
649 0 : s% opacity(k) = s% opacity_min
650 0 : s% d_opacity_dlnd(k) = 0
651 0 : s% d_opacity_dlnT(k) = 0
652 : end if
653 :
654 79994 : if (is_bad_num(s% opacity(k))) then
655 0 : if (s% stop_for_bad_nums) then
656 0 : !$omp critical (star_kap_get_bad_num)
657 0 : call show_stuff()
658 0 : call mesa_error(__FILE__,__LINE__,'do_kap_for_cell')
659 : !$omp end critical (star_kap_get_bad_num)
660 : end if
661 0 : ierr = -1
662 0 : return
663 : end if
664 :
665 :
666 : !test_partials = (k == s% solver_test_partials_k)
667 79994 : test_partials = .false.
668 :
669 159988 : if (test_partials) then
670 : s% solver_test_partials_val = s% opacity(k)
671 : s% solver_test_partials_var = s% i_lnd
672 : s% solver_test_partials_dval_dx = s% d_opacity_dlnd(k)
673 : write(*,*) 'do_kap_for_cell', s% solver_test_partials_var
674 : end if
675 :
676 :
677 :
678 : contains
679 :
680 0 : subroutine show_stuff()
681 79994 : use star_utils, only: write_eos_call_info
682 :
683 : include 'formats'
684 :
685 : real(dp) :: xc, xo, xh, xhe, Z
686 : integer :: i, iz
687 :
688 0 : write(*,*) 'eos info'
689 0 : call write_eos_call_info(s,k)
690 :
691 0 : write(*,*) 'kap info'
692 0 : call get_XYZ(s, s% xa(:,k), xh, xhe, Z)
693 0 : xc = 0; xo = 0
694 0 : do i=1, s% species
695 0 : iz = chem_isos% Z(s% chem_id(i))
696 0 : select case(iz)
697 : case (6)
698 0 : xc = xc + s% xa(i,k)
699 : case (8)
700 0 : xo = xo + s% xa(i,k)
701 : end select
702 : end do
703 0 : write(*,2) 'show opacity info', k
704 0 : write(*,*) 'kap_option ' // trim(kap_option_str(s% kap_rq% kap_option))
705 0 : write(*,*) 'kap_CO_option ' // trim(kap_CO_option_str(s% kap_rq% kap_CO_option))
706 0 : write(*,*) 'kap_lowT_option ' // trim(kap_lowT_option_str(s% kap_rq% kap_lowT_option))
707 0 : write(*,1) 'logT = ', log10_T
708 0 : write(*,1) 'logRho = ', log10_rho
709 0 : write(*,1) 'z = ', z
710 0 : write(*,1) 'Zbase = ', s% kap_rq% Zbase
711 0 : write(*,1) 'xh = ', xh
712 0 : write(*,1) 'xc = ', xc
713 0 : write(*,1) 'xo = ', xo
714 0 : write(*,1) 'lnfree_e = ', lnfree_e
715 0 : write(*,1) 'd_lnfree_e_dlnRho = ', d_lnfree_e_dlnRho
716 0 : write(*,1) 'd_lnfree_e_dlnT = ', d_lnfree_e_dlnT
717 0 : write(*,1) 'abar = ', s% abar(k)
718 0 : write(*,1) 'zbar = ', s% zbar(k)
719 0 : write(*,'(A)')
720 0 : write(*,*) 'cubic_interpolation_in_X = ', s% kap_rq% cubic_interpolation_in_X
721 0 : write(*,*) 'cubic_interpolation_in_Z = ', s% kap_rq% cubic_interpolation_in_Z
722 0 : write(*,*) 'include_electron_conduction = ', s% kap_rq% include_electron_conduction
723 0 : write(*,*) 'use_Zbase_for_Type1 = ', s% kap_rq% use_Zbase_for_Type1
724 0 : write(*,*) 'use_Type2_opacities = ', s% kap_rq% use_Type2_opacities
725 0 : write(*,1) 'kap_Type2_full_off_X = ', s% kap_rq% kap_Type2_full_off_X
726 0 : write(*,1) 'kap_Type2_full_on_X = ', s% kap_rq% kap_Type2_full_on_X
727 0 : write(*,1) 'kap_Type2_full_off_dZ = ', s% kap_rq% kap_Type2_full_off_dZ
728 0 : write(*,1) 'kap_Type2_full_on_dZ = ', s% kap_rq% kap_Type2_full_on_dZ
729 0 : write(*,'(A)')
730 0 : write(*,'(A)')
731 0 : write(*,1) 'rho = ', s% rho(k)
732 0 : write(*,1) 'lnrho = ', s% lnd(k)
733 0 : write(*,1) 'T = ', s% T(k)
734 0 : write(*,1) 'lnT = ', s% lnT(k)
735 0 : write(*,1) 'logKap = ', log10(s% opacity(k))
736 0 : write(*,1) 'opacity = ', s% opacity(k)
737 0 : write(*,1) 'dlnkap_dlnd = ', dlnkap_dlnd
738 0 : write(*,1) 'dlnkap_dlnT = ', dlnkap_dlnT
739 0 : write(*,1) 'd_opacity_dlnd = ', s% d_opacity_dlnd(k)
740 0 : write(*,1) 'd_opacity_dlnT = ', s% d_opacity_dlnT(k)
741 0 : write(*,'(A)')
742 0 : write(*,1) 'logQ = ', s% lnd(k)/ln10 - 2*s% lnT(k)/ln10 + 12
743 0 : write(*,'(A)')
744 0 : write(*,1) 'kap_frac_lowT', s% kap_frac_lowT(k)
745 0 : write(*,1) 'kap_frac_highT', s% kap_frac_highT(k)
746 0 : write(*,1) 'kap_frac_Type2', s% kap_frac_Type2(k)
747 0 : write(*,1) 'kap_frac_Compton', s% kap_frac_Compton(k)
748 0 : write(*,'(A)')
749 0 : write(*,1) 'extra_opacity_factor', s% extra_opacity_factor(k)
750 0 : write(*,'(A)')
751 0 : write(*,'(A)')
752 0 : end subroutine show_stuff
753 :
754 : end subroutine do_kap_for_cell
755 :
756 :
757 1 : subroutine shutdown_microphys
758 : use eos_lib, only: eos_shutdown
759 : use kap_lib, only: kap_shutdown
760 : use net_lib
761 1 : call eos_shutdown
762 1 : call kap_shutdown
763 1 : call net_shutdown
764 1 : end subroutine shutdown_microphys
765 :
766 : end module micro
|
