Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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 net
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, ln10, pi4
24 : use utils_lib, only: is_bad, mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: set_net
30 : public :: do_net
31 : public :: do1_net
32 : public :: do_micro_change_net
33 : public :: get_screening_mode
34 : public :: default_set_rate_factors
35 : public :: default_set_op_mono_factors
36 :
37 : contains
38 :
39 66 : subroutine do_net(s, nzlo, nzhi, ierr)
40 : use star_utils, only: start_time, update_time
41 : use net_def, only: net_other_net_derivs
42 : use rates_def, only: rates_other_screening
43 : use alloc
44 : type (star_info), pointer :: s
45 : integer, intent(in) :: nzlo, nzhi
46 : integer, intent(out) :: ierr
47 :
48 : logical, parameter :: use_omp = .true.
49 : integer :: k, op_err
50 : integer(i8) :: time0
51 66 : real(dp) :: total
52 : logical, parameter :: only_dlnT = .false.
53 : logical :: okay, check_op_split_burn
54 :
55 : include 'formats'
56 :
57 66 : ierr = 0
58 :
59 66 : if (s% eps_nuc_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) then
60 0 : do k = nzlo, nzhi
61 0 : s% eps_nuc(k) = 0d0
62 0 : s% d_epsnuc_dlnd(k) = 0d0
63 0 : s% d_epsnuc_dlnT(k) = 0d0
64 0 : s% d_epsnuc_dx(:,k) = 0d0
65 0 : s% eps_nuc_categories(:,k) = 0d0
66 0 : s% dxdt_nuc(:,k) = 0d0
67 0 : s% d_dxdt_nuc_dRho(:,k) = 0d0
68 0 : s% d_dxdt_nuc_dT(:,k) = 0d0
69 0 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
70 0 : s% eps_nuc_neu_total(k) = 0d0
71 0 : if (s% op_split_burn) then
72 0 : s% burn_avg_epsnuc(k) = 0d0
73 0 : s% burn_num_iters(k) = 0
74 : end if
75 : end do
76 : return
77 : end if
78 :
79 66 : rates_other_screening => null()
80 66 : if(s% use_other_screening) then
81 0 : rates_other_screening => s% other_screening
82 : end if
83 :
84 66 : net_other_net_derivs => null()
85 66 : if(s% use_other_net_derivs) then
86 0 : if(index(s% net_name,'approx')>0) then
87 0 : write(*,*) 'use_other_net_derivs does not work with approx nets'
88 0 : ierr = -1
89 0 : return
90 : end if
91 0 : net_other_net_derivs => s% other_net_derivs
92 : end if
93 :
94 66 : check_op_split_burn = s% op_split_burn
95 :
96 66 : if (nzlo == nzhi) then
97 : call do1_net(s, nzlo, s% species, &
98 : s% num_reactions, &
99 0 : check_op_split_burn, ierr)
100 0 : return
101 : end if
102 :
103 66 : if (s% doing_timing) call start_time(s, time0, total)
104 : if (use_omp) then
105 66 : okay = .true.
106 66 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
107 : do k = nzlo, nzhi
108 : if (.not. okay) cycle
109 : op_err = 0
110 : call do1_net(s, k, s% species, &
111 : s% num_reactions, &
112 : check_op_split_burn, op_err)
113 : if (op_err /= 0) okay = .false.
114 : end do
115 : !$OMP END PARALLEL DO
116 66 : if (.not. okay) ierr = -1
117 : else
118 : do k = nzlo, nzhi
119 : call do1_net(s, k, s% species, &
120 : s% num_reactions, &
121 : check_op_split_burn, ierr)
122 : if (ierr /= 0) exit
123 : end do
124 : end if
125 66 : if (s% doing_timing) call update_time(s, time0, total, s% time_nonburn_net)
126 :
127 66 : end subroutine do_net
128 :
129 :
130 79994 : subroutine do1_net(s, k, species, &
131 : num_reactions, check_op_split_burn, ierr)
132 66 : use rates_def, only: std_reaction_Qs, std_reaction_neuQs, i_rate
133 : use net_def, only: Net_Info, net_test_partials, &
134 : net_test_partials_val, net_test_partials_dval_dx, net_test_partials_i, &
135 : net_test_partials_iother, get_net_ptr
136 : use net_lib, only: net_get
137 : use star_utils, only: lookup_nameofvar
138 : use chem_def, only: category_name, i_ni56_co56, i_co56_fe56, &
139 : num_categories, iphoto, category_name
140 : use eos_def, only : i_eta
141 : use utils_lib,only: realloc_double, realloc_double3
142 : type (star_info), pointer :: s
143 : integer, intent(in) :: k, species, num_reactions
144 : logical, intent(in) :: check_op_split_burn
145 : integer, intent(out) :: ierr
146 :
147 : integer :: i, j, kk, screening_mode, sz, i_var, i_var_sink
148 79994 : real(dp) :: log10_rho, log10_T, T, alfa, eps_nuc_factor, &
149 79994 : d_eps_nuc_dRho, d_eps_nuc_dT, tau_gamma, eps_cat_sum
150 79994 : type (Net_Info) :: n
151 79994 : real(dp), pointer :: reaction_neuQs(:)
152 : logical :: clipped_T
153 :
154 : logical, parameter :: dbg = .false.
155 :
156 : include 'formats'
157 :
158 79994 : ierr = 0
159 :
160 : if (check_op_split_burn .and. &
161 0 : s% doing_struct_burn_mix .and. &
162 : s% T_start(k) >= s% op_split_burn_min_T) then ! leave this to do_burn
163 : return
164 : end if
165 :
166 79994 : n% star_id = s% id
167 79994 : n% zone = k
168 :
169 79994 : s% eps_nuc(k) = 0d0
170 79994 : s% d_epsnuc_dlnd(k) = 0d0
171 79994 : s% d_epsnuc_dlnT(k) = 0d0
172 719946 : s% d_epsnuc_dx(:,k) = 0d0
173 1999850 : s% eps_nuc_categories(:,k) = 0d0
174 719946 : s% dxdt_nuc(:,k) = 0d0
175 719946 : s% d_dxdt_nuc_dRho(:,k) = 0d0
176 719946 : s% d_dxdt_nuc_dT(:,k) = 0d0
177 5839562 : s% d_dxdt_nuc_dx(:,:,k) = 0d0
178 79994 : s% eps_nuc_neu_total(k) = 0d0
179 :
180 79994 : if ((s% eps_nuc_factor == 0d0 .and. s% dxdt_nuc_factor == 0d0) .or. &
181 : s% abar(k) > s% max_abar_for_burning) then
182 : return
183 : end if
184 :
185 79994 : log10_rho = s% lnd(k)/ln10
186 79994 : log10_T = s% lnT(k)/ln10
187 79994 : T = s% T(k)
188 :
189 79994 : clipped_T = (s% max_logT_for_net > 0 .and. log10_T > s% max_logT_for_net)
190 : if (clipped_T) then
191 0 : log10_T = s% max_logT_for_net
192 0 : T = exp10(log10_T)
193 : end if
194 :
195 79994 : screening_mode = get_screening_mode(s,ierr)
196 79994 : if (ierr /= 0) then
197 0 : write(*,*) 'unknown string for screening_mode: ' // trim(s% screening_mode)
198 0 : call mesa_error(__FILE__,__LINE__,'do1_net')
199 0 : return
200 : end if
201 :
202 79994 : if (s% reaction_neuQs_factor /= 1d0) then
203 0 : sz = size(std_reaction_neuQs,dim=1)
204 0 : allocate(reaction_neuQs(sz))
205 0 : do j=1,sz
206 0 : reaction_neuQs(j) = std_reaction_neuQs(j)*s% reaction_neuQs_factor
207 : end do
208 : else
209 79994 : reaction_neuQs => std_reaction_neuQs
210 : end if
211 :
212 79994 : if (s% solver_test_net_partials) then
213 : net_test_partials = (k == s% solver_test_partials_k .and. &
214 : s% solver_call_number == s% solver_test_partials_call_number .and. &
215 0 : s% solver_iter == s% solver_test_partials_iter_number)
216 : ! if the test is for a partial wrt an abundance, do this
217 : ! in inlist set solver_test_partials_var_name and solver_test_partials_sink_name
218 : ! set solver_test_partials_equ_name = ''
219 0 : i_var = lookup_nameofvar(s, s% solver_test_partials_var_name)
220 0 : i_var_sink = lookup_nameofvar(s, s% solver_test_partials_sink_name)
221 0 : s% solver_test_partials_var = i_var ! index in vars
222 0 : if (i_var > s% nvar_hydro) then ! index in xa for sink
223 0 : s% solver_test_partials_dx_sink = i_var_sink - s% nvar_hydro
224 : else
225 0 : s% solver_test_partials_dx_sink = 0
226 : end if
227 0 : net_test_partials_i = i_var - s% nvar_hydro ! index in xa for var
228 0 : net_test_partials_iother = i_var_sink - s% nvar_hydro ! index in xa for var
229 : end if
230 :
231 79994 : if (s% use_other_net_get) then
232 : call s% other_net_get( &
233 : s% id, k, &
234 : s% net_handle, .false., n, species, num_reactions, s% xa(1:species,k), &
235 : T, log10_T, s% rho(k), log10_Rho, &
236 : s% abar(k), s% zbar(k), s% z2bar(k), s% ye(k), &
237 : s% eta(k), s% d_eos_dlnT(i_eta,k), s% d_eos_dlnd(i_eta,k), &
238 : s% rate_factors, s% weak_rate_factor, &
239 : std_reaction_Qs, reaction_neuQs, &
240 : s% eps_nuc(k), d_eps_nuc_dRho, d_eps_nuc_dT, s% d_epsnuc_dx(:,k), &
241 : s% dxdt_nuc(:,k), s% d_dxdt_nuc_dRho(:,k), s% d_dxdt_nuc_dT(:,k), s% d_dxdt_nuc_dx(:,:,k), &
242 : screening_mode, s% eps_nuc_categories(:,k), &
243 0 : s% eps_nuc_neu_total(k), ierr)
244 : else
245 : call net_get( &
246 : s% net_handle, .false., n, species, num_reactions, s% xa(1:species,k), &
247 : T, log10_T, s% rho(k), log10_Rho, &
248 : s% abar(k), s% zbar(k), s% z2bar(k), s% ye(k), &
249 : s% eta(k), s% d_eos_dlnT(i_eta,k), s% d_eos_dlnd(i_eta,k), &
250 : s% rate_factors, s% weak_rate_factor, &
251 : std_reaction_Qs, reaction_neuQs, &
252 : s% eps_nuc(k), d_eps_nuc_dRho, d_eps_nuc_dT, s% d_epsnuc_dx(:,k), &
253 : s% dxdt_nuc(:,k), s% d_dxdt_nuc_dRho(:,k), s% d_dxdt_nuc_dT(:,k), s% d_dxdt_nuc_dx(:,:,k), &
254 : screening_mode, s% eps_nuc_categories(:,k), &
255 79994 : s% eps_nuc_neu_total(k), ierr)
256 : end if
257 :
258 :
259 79994 : if (is_bad(s% eps_nuc(k))) then
260 0 : ierr = -1
261 0 : if (s% report_ierr) write(*,*) 'net_get returned bad eps_nuc', ierr
262 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_net')
263 0 : return
264 : end if
265 :
266 79994 : call save_rates(s, n, k, ierr)
267 79994 : if(ierr/=0) return
268 :
269 79994 : if (-k == s% nz) then
270 0 : write(*,1) 'logT', log10_T
271 0 : call mesa_error(__FILE__,__LINE__,'net')
272 : end if
273 : !if (k == 864 .and. log10_T >= 7.522497408d0 .and. log10_T <= 7.5224974089d0) then
274 79994 : if (-k == s% nz) then
275 0 : do j=1,num_categories
276 0 : write(*,2) trim(category_name(j)), j, s% eps_nuc_categories(j,k)
277 : end do
278 0 : write(*,'(A)')
279 0 : write(*,1) 'logRho', log10_Rho
280 0 : write(*,1) 'logT', log10_T
281 0 : write(*,1) 'eps_nuc', s% eps_nuc(k)
282 0 : write(*,1) 'sum(eps_nuc_categories)', sum(s% eps_nuc_categories(:,k))
283 0 : write(*,1) 'sum(eps_nuc_categories)/eps_nuc', &
284 0 : sum(s% eps_nuc_categories(:,k))/s% eps_nuc(k)
285 0 : write(*,2) trim(s% net_name), s% species
286 0 : call mesa_error(__FILE__,__LINE__,'after net_get in star')
287 : end if
288 :
289 79994 : if (s% solver_test_net_partials .and. net_test_partials) then
290 0 : s% solver_test_partials_val = net_test_partials_val
291 0 : s% solver_test_partials_dval_dx = net_test_partials_dval_dx
292 : end if
293 :
294 79994 : if (ierr == 0) then
295 :
296 79994 : if (clipped_T) then
297 0 : d_eps_nuc_dT = 0
298 0 : s% d_dxdt_nuc_dT(1:species,k) = 0
299 : end if
300 :
301 79994 : if (s% nonlocal_NiCo_kap_gamma > 0d0 .and. &
302 : .not. s% nonlocal_NiCo_decay_heat) then
303 : tau_gamma = 0
304 0 : do kk = 1, k
305 0 : tau_gamma = tau_gamma + s% dm(kk)/(pi4*s% rmid(kk)*s% rmid(kk))
306 : end do
307 0 : tau_gamma = tau_gamma*s% nonlocal_NiCo_kap_gamma
308 0 : s% eps_nuc(k) = s% eps_nuc(k)*(1d0 - exp(-tau_gamma))
309 : end if
310 :
311 79994 : if (abs(s% eps_nuc(k)) > s% max_abs_eps_nuc) then
312 0 : s% eps_nuc(k) = sign(s% max_abs_eps_nuc, s% eps_nuc(k))
313 0 : d_eps_nuc_dRho = 0d0
314 0 : d_eps_nuc_dT = 0d0
315 0 : s% d_epsnuc_dx(:,k) = 0d0
316 : end if
317 :
318 1999850 : eps_cat_sum = sum(s% eps_nuc_categories(:,k))
319 79994 : if (abs(eps_cat_sum) < 1d-10) then
320 : alfa = 1d0
321 : else
322 44923 : alfa = s% eps_nuc(k)/eps_cat_sum
323 : end if
324 : if (.false. .and. abs(1d0 - alfa) > 1d-10) then
325 : write(*,3) 'do1_net: sum(categories) /= eps_nuc', k, s% model_number, &
326 : log10(abs(alfa)), s% eps_nuc(k), eps_cat_sum, s% eps_nuc_categories(iphoto,k)
327 : else
328 1999850 : do i=1,num_categories
329 1999850 : s% eps_nuc_categories(i,k) = s% eps_nuc_categories(i,k)*alfa
330 : end do
331 : end if
332 :
333 : end if
334 :
335 79994 : if (s% reaction_neuQs_factor /= 1d0) deallocate(reaction_neuQs)
336 :
337 79994 : if (ierr /= 0) then
338 0 : if (s% report_ierr) then
339 0 : write(*,'(A)')
340 0 : write(*,*) 'do1_net: net_get failure for cell ', k
341 : !return
342 0 : call show_stuff(s,k)
343 : end if
344 0 : if (is_bad_num(s% eps_nuc(k))) then
345 0 : if (s% stop_for_bad_nums) then
346 0 : write(*,2) 'eps_nuc', k, s% eps_nuc(k)
347 0 : call mesa_error(__FILE__,__LINE__,'do1_net')
348 : end if
349 : end if
350 0 : return
351 : end if
352 :
353 79994 : if (is_bad_num(s% eps_nuc(k))) then
354 0 : if (s% stop_for_bad_nums) then
355 0 : write(*,2) 'eps_nuc', k, s% eps_nuc(k)
356 0 : call mesa_error(__FILE__,__LINE__,'do1_net')
357 : end if
358 0 : ierr = -1
359 0 : return
360 : end if
361 :
362 79994 : s% d_epsnuc_dlnd(k) = d_eps_nuc_dRho*s% rho(k)
363 79994 : s% d_epsnuc_dlnT(k) = d_eps_nuc_dT*s% T(k)
364 :
365 79994 : eps_nuc_factor = s% eps_nuc_factor
366 79994 : if (eps_nuc_factor /= 1d0) then
367 0 : s% eps_nuc(k) = s% eps_nuc(k)*eps_nuc_factor
368 0 : s% d_epsnuc_dlnd(k) = s% d_epsnuc_dlnd(k)*eps_nuc_factor
369 0 : s% d_epsnuc_dlnT(k) = s% d_epsnuc_dlnT(k)*eps_nuc_factor
370 0 : s% d_epsnuc_dx(:,k) = s% d_epsnuc_dx(:,k)*eps_nuc_factor
371 0 : s% eps_nuc_categories(:,k) = s% eps_nuc_categories(:,k)*eps_nuc_factor
372 : end if
373 :
374 79994 : if (s% dxdt_nuc_factor /= 1d0) then
375 0 : s% dxdt_nuc(:,k) = s% dxdt_nuc(:,k)*s% dxdt_nuc_factor
376 0 : s% d_dxdt_nuc_dRho(:,k) = s% d_dxdt_nuc_dRho(:,k)*s% dxdt_nuc_factor
377 0 : s% d_dxdt_nuc_dT(:,k) = s% d_dxdt_nuc_dT(:,k)*s% dxdt_nuc_factor
378 0 : s% d_dxdt_nuc_dx(:,:,k) = s% d_dxdt_nuc_dx(:,:,k)*s% dxdt_nuc_factor
379 : end if
380 :
381 79994 : if (is_bad_num(s% eps_nuc(k))) then
382 0 : write(*,*) 'k', k
383 0 : write(*,1) 's% eps_nuc(k)', s% eps_nuc(k)
384 0 : ierr = -1
385 0 : call show_stuff(s,k)
386 0 : write(*,*) '(is_bad_num(s% eps_nuc(k)))'
387 0 : write(*,*) 'failed in do1_net'
388 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_net')
389 0 : return
390 : end if
391 :
392 79994 : if (k == -1) then
393 0 : write(*,'(A)')
394 0 : call show_stuff(s,k)
395 : end if
396 :
397 79994 : if (s% model_number == -1) then
398 0 : write(*,5) 'eps_nuc', k, s% solver_iter, s% model_number, s% solver_adjust_iter, &
399 0 : s% eps_nuc(k)
400 : end if
401 :
402 : if (.false.) then
403 : write(*,'(A)')
404 : call show_stuff(s,k)
405 : write(*,'(A)')
406 : write(*,'(A)')
407 : write(*,1) 's% eps_nuc(k)', s% eps_nuc(k)
408 : write(*,1) 's% d_epsnuc_dlnd(k)', s% d_epsnuc_dlnd(k)
409 : write(*,1) 's% d_epsnuc_dlnT(k)', s% d_epsnuc_dlnT(k)
410 : write(*,'(A)')
411 : write(*,*) 'do1_net'
412 : stop
413 : !ierr = -1
414 : end if
415 :
416 : if (.false.) call show_stuff(s,k)
417 :
418 159988 : end subroutine do1_net
419 :
420 :
421 :
422 0 : subroutine show_stuff(s,k)
423 79994 : use chem_def
424 : use eos_def, only : i_eta
425 : use rates_def
426 : use net_lib, only: get_reaction_id_table_ptr
427 : use num_lib, only: qsort
428 : use eos_def, only: i_eta
429 : type (star_info), pointer :: s
430 : integer, intent(in) :: k
431 :
432 0 : integer, pointer :: reaction_id(:) ! maps net reaction number to reaction id
433 : integer :: i, j, ierr, species, num_reactions
434 : real(dp) :: log10_Rho, log10_T
435 0 : real(dp), pointer :: v(:)
436 0 : integer, pointer :: index(:)
437 0 : real(dp), pointer, dimension(:) :: rate_screened, rate_raw
438 :
439 : include 'formats'
440 :
441 : logical, parameter :: do_sort = .true.
442 :
443 : ierr = 0
444 0 : species = s% species
445 0 : num_reactions = s% num_reactions
446 0 : log10_T = s% lnT(k)/ln10
447 0 : log10_Rho = s% lnd(k)/ln10
448 :
449 0 : call get_reaction_id_table_ptr(s% net_handle, reaction_id, ierr)
450 0 : if (ierr /= 0) return
451 :
452 0 : do i=1,num_reactions
453 0 : if (s% rate_factors(i) /= 1d0) then
454 0 : write(*,2) 'rate factor ' // trim(reaction_Name(reaction_id(i))), &
455 0 : s% rate_factors(i)
456 : end if
457 : end do
458 :
459 0 : i = max(species, num_reactions)
460 0 : allocate(v(i), index(i))
461 0 : write(*,'(A)')
462 : if (.true.) then
463 0 : write(*, *)
464 : if (do_sort) then
465 0 : do j=1,num_reactions
466 0 : v(j) = abs(rate_raw(j))
467 : end do
468 0 : call qsort(index, num_reactions, v)
469 : else
470 : do j=1,num_reactions
471 : index(j) = j
472 : end do
473 : end if
474 :
475 0 : write(*,*) 'reaction rate_raw'
476 0 : do i=1,num_reactions
477 0 : j = index(num_reactions+1-i)
478 0 : write(*,2) trim(reaction_Name(reaction_id(j))), k, rate_raw(j)
479 : end do
480 : end if
481 :
482 : if (.false.) then
483 : write(*,'(A)')
484 : write(*,*) 'screened rates'
485 : do j=1,num_reactions
486 : write(*,3) 'screened rate ' // trim(reaction_Name(reaction_id(j))), &
487 : j, k, rate_screened(j)
488 : end do
489 : end if
490 :
491 : if (.true.) then
492 0 : write(*,'(A)')
493 0 : do j=1,species
494 0 : write(*,2) 'dxdt ' // trim(chem_isos% name(s% chem_id(j))), k, s% dxdt_nuc(j, k)
495 : end do
496 0 : write(*,'(A)')
497 0 : do j=1,species
498 0 : write(*,2) 'd_epsnuc_dx ' // trim(chem_isos% name(s% chem_id(j))), k, s% d_epsnuc_dx(j, k)
499 : end do
500 : end if
501 0 : write(*,'(A)')
502 :
503 : if (.false.) then
504 : write(*,'(A)')
505 : do j=1,species
506 : write(*,2) 'dt*dxdt ' // trim(chem_isos% name(s% chem_id(j))), k, &
507 : s% dt * s% dxdt_nuc(j, k)
508 : end do
509 : end if
510 :
511 : if (.true.) then
512 : if (do_sort) then
513 0 : do j=1,species
514 0 : v(j) = s% xa(j,k)
515 : end do
516 0 : call qsort(index, species, v)
517 : else
518 : do j=1,num_reactions
519 : index(j) = j
520 : end do
521 : end if
522 0 : write(*,'(A)')
523 0 : do i=1,species
524 0 : j = index(species+1-i)
525 : if (.true. .or. s% xa(j,k) > 1d-9) &
526 : write(*,1) 'xin(net_iso(i' // &
527 0 : trim(chem_isos% name(s% chem_id(j))) // '))= ', s% xa(j,k)
528 : end do
529 : end if
530 :
531 : if (.false.) then
532 : do i=1,species
533 : write(*,'(a,i3,a,d26.16)') 'values_for_Xinit(', i, ')= ', s% xa(i,k)
534 : end do
535 : end if
536 :
537 0 : write(*,2) 'k', k
538 0 : write(*,'(A)')
539 0 : write(*,*) 'net_name ', trim(s% net_name)
540 0 : write(*,*) 'species', species
541 0 : write(*,'(A)')
542 0 : write(*,1) 'logT =', log10_T
543 0 : write(*,1) 'T =', s% T(k)
544 0 : write(*,1) 'logRho =', log10_Rho
545 0 : write(*,1) 'rho =', s% rho(k)
546 0 : write(*,'(A)')
547 0 : write(*,1) 'eta =', s% eta(k)
548 0 : write(*,1) 'd_eta_lnT =', s% d_eos_dlnT(i_eta,k)
549 0 : write(*,1) 'd_eta_lnd =', s% d_eos_dlnd(i_eta,k)
550 0 : write(*,'(A)')
551 0 : write(*,1) 'abar =', s% abar(k)
552 0 : write(*,1) 'zbar =', s% zbar(k)
553 0 : write(*,1) 'z2bar =', s% z2bar(k)
554 0 : write(*,1) 'ye =', s% ye(k)
555 0 : write(*,*) 'screening_mode = ' // trim(s% screening_mode)
556 0 : write(*,'(A)')
557 :
558 0 : end subroutine show_stuff
559 :
560 :
561 :
562 80005 : integer function get_screening_mode(s,ierr)
563 0 : use rates_lib, only: screening_option
564 : type (star_info), pointer :: s
565 : integer, intent(out) :: ierr
566 : include 'formats'
567 80005 : ierr = 0
568 80005 : if (s% screening_mode_value >= 0) then
569 80005 : get_screening_mode = s% screening_mode_value
570 : return
571 : end if
572 44 : get_screening_mode = screening_option(s% screening_mode, ierr)
573 44 : if (ierr /= 0) return
574 44 : s% screening_mode_value = get_screening_mode
575 44 : end function get_screening_mode
576 :
577 :
578 0 : subroutine do_micro_change_net(s, new_net_name, ierr)
579 : use net_def
580 : type (star_info), pointer :: s
581 : character (len=*), intent(in) :: new_net_name
582 : integer, intent(out) :: ierr
583 0 : ierr = 0
584 0 : s% net_name = new_net_name
585 0 : call set_net(s, new_net_name, ierr)
586 0 : end subroutine do_micro_change_net
587 :
588 :
589 1 : subroutine set_net(s, new_net_name, ierr)
590 : use net_lib
591 : use utils_lib, only: realloc_double
592 : use alloc, only: update_nvar_allocs, set_chem_names
593 : use chem_def, only: ih1, ihe4
594 : use rates_def
595 : type (star_info), pointer :: s
596 : character (len=*), intent(in) :: new_net_name
597 : integer, intent(out) :: ierr
598 :
599 : integer :: old_num_reactions, old_nvar_chem, old_species
600 : integer, parameter :: num_lowT_rates = 10
601 :
602 : include 'formats'
603 :
604 1 : old_num_reactions = s% num_reactions
605 :
606 1 : if (s% net_handle /= 0) call free_net_handle(s% net_handle)
607 :
608 1 : s% net_handle = alloc_net_handle(ierr)
609 1 : if (ierr /= 0) then
610 0 : write(*,*) 'set_net failed in alloc_net_handle'
611 0 : return
612 : end if
613 1 : call net_ptr(s% net_handle, s% net_rq, ierr)
614 1 : if (ierr /= 0) then
615 0 : write(*,*) 'set_net failed in net_ptr'
616 0 : return
617 : end if
618 :
619 1 : call net_tables(s, ierr)
620 1 : if (ierr /= 0) then
621 0 : write(*,*) 'set_net failed in net_tables'
622 0 : return
623 : end if
624 :
625 1 : old_species = s% species
626 1 : s% species = net_num_isos(s% net_handle, ierr)
627 1 : if (ierr /= 0) then
628 0 : write(*,*) 'set_net failed in net_num_isos'
629 0 : return
630 : end if
631 :
632 1 : s% num_reactions = net_num_reactions(s% net_handle, ierr)
633 1 : if (ierr /= 0) then
634 0 : write(*,*) 'set_net failed in net_num_reactions'
635 0 : return
636 : end if
637 :
638 1 : old_nvar_chem = s% nvar_chem
639 1 : s% nvar_chem = s% species
640 1 : call update_nvar_allocs(s, s% nvar_hydro, old_nvar_chem, ierr)
641 1 : if (ierr /= 0) then
642 0 : write(*,*) 'set_net failed in update_nvar_allocs'
643 0 : return
644 : end if
645 :
646 1 : call get_chem_id_table_ptr(s% net_handle, s% chem_id, ierr)
647 1 : if (ierr /= 0) then
648 0 : write(*,*) 'set_net failed in get_chem_id_table_ptr'
649 0 : return
650 : end if
651 :
652 1 : call get_net_iso_table_ptr(s% net_handle, s% net_iso, ierr)
653 1 : if (ierr /= 0) then
654 0 : write(*,*) 'set_net failed in get_net_iso_table_ptr'
655 0 : return
656 : end if
657 :
658 1 : if (s% net_iso(ih1) == 0 .or. s% net_iso(ihe4) == 0) then
659 0 : write(*,*) 'mesa/star requires both h1 and he4 in net isotopes'
660 0 : write(*,*) 'but they are not included in ' // trim(new_net_name)
661 0 : ierr = -1
662 0 : return
663 : end if
664 :
665 1 : if (associated(s% xa_removed)) deallocate(s% xa_removed)
666 1 : allocate(s% xa_removed(s% species))
667 :
668 1 : call set_chem_names(s)
669 :
670 1 : call s% set_rate_factors(s% id, ierr)
671 1 : if (ierr /= 0) then
672 0 : write(*,*) 'set_net failed in s% set_rate_factors'
673 0 : return
674 : end if
675 :
676 1 : if (associated(s% op_mono_factors)) deallocate(s% op_mono_factors)
677 1 : allocate(s% op_mono_factors(s% species))
678 :
679 1 : call s% set_op_mono_factors(s% id, ierr)
680 1 : if (ierr /= 0) then
681 0 : write(*,*) 'set_net failed in s% set_op_mono_factors'
682 0 : return
683 : end if
684 :
685 1 : s% net_rq% use_3a_fl87 = s% job% use_3a_fl87
686 :
687 1 : s% need_to_setvars = .true.
688 :
689 1 : s% net_rq% fill_arrays_with_nans = s% fill_arrays_with_NaNs
690 :
691 6 : end subroutine set_net
692 :
693 :
694 1 : subroutine net_tables(s, ierr)
695 1 : use net_lib ! setup net
696 : use rates_lib
697 : use rates_def, only: rates_reaction_id_max, rates_other_rate_get
698 : type (star_info), pointer :: s
699 : integer, intent(out) :: ierr
700 : ierr = 0
701 :
702 1 : call net_start_def(s% net_handle, ierr)
703 1 : if (ierr /= 0) then
704 0 : if (s% report_ierr) write(*,*) 'failed in net_start_def'
705 0 : return
706 : end if
707 :
708 1 : if (len_trim(s% net_name) == 0) then
709 0 : write(*,*) 'missing net_name -- please set it and try again'
710 0 : ierr = -1
711 0 : return
712 : end if
713 :
714 1 : call read_net_file(s% net_name, s% net_handle, ierr)
715 1 : if (ierr /= 0) then
716 0 : if (s% report_ierr) write(*,*) 'failed in read_net_file ' // trim(s% net_name)
717 0 : return
718 : end if
719 :
720 1 : call net_finish_def(s% net_handle, ierr)
721 1 : if (ierr /= 0) then
722 0 : if (s% report_ierr) write(*,*) 'failed in net_finish_def'
723 0 : return
724 : end if
725 :
726 1 : if (associated(s% rate_factors)) deallocate(s% rate_factors)
727 1 : allocate(s% rate_factors(rates_reaction_id_max))
728 :
729 1 : call s% set_rate_factors(s% id, ierr)
730 1 : if (ierr /= 0) then
731 0 : if (s% report_ierr) write(*,*) 'failed in set_rate_factors'
732 0 : return
733 : end if
734 :
735 1 : call read_rates_from_files(s% job% reaction_for_special_factor, s% job% filename_of_special_rate, ierr)
736 1 : if (ierr /= 0) then
737 0 : if (s% report_ierr) write(*,*) 'failed in read_rates_from_files'
738 0 : return
739 : end if
740 :
741 1 : call net_set_logTcut(s% net_handle, s% net_logTcut_lo, s% net_logTcut_lim, ierr)
742 1 : if (ierr /= 0) then
743 0 : if (s% report_ierr) write(*,*) 'failed in net_set_logTcut'
744 0 : return
745 : end if
746 :
747 : call net_set_fe56ec_fake_factor( &
748 1 : s% net_handle, s% fe56ec_fake_factor, s% min_T_for_fe56ec_fake_factor, ierr)
749 1 : if (ierr /= 0) then
750 0 : if (s% report_ierr) write(*,*) 'failed in net_set_fe56ec_fake_factor'
751 0 : return
752 : end if
753 :
754 1 : rates_other_rate_get => null()
755 1 : if(s% use_other_rate_get) then
756 0 : rates_other_rate_get => s% other_rate_get
757 : end if
758 :
759 : call net_setup_tables( &
760 1 : s% net_handle, rates_cache_suffix_for_star, ierr)
761 1 : if (ierr /= 0) then
762 0 : if (s% report_ierr) write(*,*) 'failed in net_setup_tables'
763 0 : return
764 : end if
765 :
766 7 : end subroutine net_tables
767 :
768 0 : subroutine default_set_rate_factors(id, ierr)
769 : integer, intent(in) :: id
770 : integer, intent(out) :: ierr
771 : type (star_info), pointer :: s
772 : ierr = 0
773 0 : call get_star_ptr(id, s, ierr)
774 0 : if (ierr /= 0) return
775 0 : s% rate_factors(:) = 1
776 1 : end subroutine default_set_rate_factors
777 :
778 :
779 1 : subroutine default_set_op_mono_factors(id, ierr)
780 : integer, intent(in) :: id
781 : integer, intent(out) :: ierr
782 : type (star_info), pointer :: s
783 : ierr = 0
784 1 : call get_star_ptr(id, s, ierr)
785 1 : if (ierr /= 0) return
786 9 : s% op_mono_factors(:) = 1
787 : end subroutine default_set_op_mono_factors
788 :
789 :
790 79994 : subroutine save_rates(s, n, k, ierr)
791 : use net_def, only: net_Info, Net_General_Info, get_net_ptr
792 : type(star_info),pointer :: s
793 : type(net_info) :: n
794 : integer,intent(in) :: k
795 : integer, intent(inout) :: ierr
796 : type(Net_General_Info), pointer :: g=> null()
797 : integer :: i
798 :
799 : ierr = 0
800 :
801 79994 : call get_net_ptr(s% net_handle, g, ierr)
802 79994 : if(ierr/=0) return
803 :
804 2559808 : do i=1, g% num_reactions
805 2479814 : s% raw_rate(i,k) = n% raw_rate(i)
806 2479814 : s% screened_rate(i,k) = n% screened_rate(i)
807 2479814 : s% eps_nuc_rate(i,k) = n% eps_nuc_rate(i)
808 2559808 : s% eps_neu_rate(i,k) = n% eps_neu_rate(i)
809 : end do
810 :
811 : end subroutine save_rates
812 :
813 : end module net
|
