Go to the documentation of this file.
69 real(
r_kind),
intent(in) :: time
70 real(
r_kind),
intent(in) :: temp
71 integer,
intent(in),
optional :: idx
72 real(
r_kind),
intent(in) :: rho
73 real(
r_kind),
intent(in) :: ye
74 real(
r_kind),
intent(in) :: rkm
75 real(
r_kind),
intent(out),
optional :: rat_calc
78 logical :: is_same_step
80 real(
r_kind) :: eta_e,eta_p
82 integer :: reac_origin
96 is_same_step = .false.
103 if (is_same_step .eqv. .false.)
then
118 call chempot(1.d9*temp,rho,ye,eta_e,eta_p)
119 mue = eta_e * temp * 1.d9 *
unit%k_MeV
125 call nuflux(time, rkm)
132 if (
present(rrate_array).eqv..false.)
return
136 src = rrate_array(idx)%source
137 reac_origin = rrate_array(idx)%reac_src
142 if (reac_origin .eq.
rrs_twr)
then
145 else if (reac_origin .eq.
rrs_tabl)
then
148 else if (reac_origin .eq.
rrs_nu)
then
151 else if (reac_origin .eq.
rrs_detb)
then
154 else if (reac_origin .eq.
rrs_fiss)
then
157 else if (reac_origin .eq.
rrs_wext)
then
160 else if (reac_origin .eq.
rrs_aext)
then
163 else if (reac_origin .eq.
rrs_reacl)
then
168 "'."//new_line(
"A")//
"Reaction source is: "//src,&
169 "calculate_reaction_rate", 280011)
176 rrate_array(idx)%cached_p = rat_calc
179 if(
iscreen .eqv. .true.) rat_calc=rat_calc*dexp(
hv(idx))
182 rat_calc = rat_calc * rrate_array(idx)%one_over_n_fac
187 if (trim(adjustl(src)) .eq.
'ec') rat_calc = rat_calc * elnd
190 select case (rrate_array(idx)%group)
192 if (rrate_array(idx)%is_reverse)
then
193 rat_calc = rat_calc * (
pf(rrate_array(idx)%parts(2))*
pf(rrate_array(idx)%parts(3))* &
194 pf(rrate_array(idx)%parts(4))*
pf(rrate_array(idx)%parts(5)))/ &
195 pf(rrate_array(idx)%parts(1))
198 rat_calc = rat_calc * rho
199 if (rrate_array(idx)%is_reverse)
then
200 rat_calc = rat_calc * (
pf(rrate_array(idx)%parts(3))*
pf(rrate_array(idx)%parts(4)) * &
201 pf(rrate_array(idx)%parts(5))*
pf(rrate_array(idx)%parts(6))) / &
202 (
pf(rrate_array(idx)%parts(1))*
pf(rrate_array(idx)%parts(2)))
205 rat_calc = rat_calc * rho**2
206 if (rrate_array(idx)%is_reverse)
then
207 rat_calc = rat_calc * (
pf(rrate_array(idx)%parts(4))*
pf(rrate_array(idx)%parts(5)) &
208 / (
pf(rrate_array(idx)%parts(1))*
pf(rrate_array(idx)%parts(2)) &
209 *
pf(rrate_array(idx)%parts(3))))
212 rat_calc = rat_calc * rho**3
213 if (rrate_array(idx)%is_reverse)
then
214 rat_calc = rat_calc * (
pf(rrate_array(idx)%parts(5))*
pf(rrate_array(idx)%parts(6))) &
215 / (
pf(rrate_array(idx)%parts(1))*
pf(rrate_array(idx)%parts(2)) &
216 *
pf(rrate_array(idx)%parts(3))*
pf(rrate_array(idx)%parts(4)))
222 "calculate_reaction_rate",280009)
233 rrate_array(idx)%cached = rat_calc
288 subroutine abchange (time, itemp, rho, Ye, rkm, Y, dYdt, evolution_mode)
293 real(
r_kind),
intent(in) :: time
294 real(
r_kind),
intent(in) :: itemp
295 real(
r_kind),
intent(in) :: rho
296 real(
r_kind),
intent(in) :: ye
297 real(
r_kind),
intent(in) :: rkm
298 real(
r_kind),
dimension(:),
intent(in) :: y
299 real(
r_kind),
dimension(:),
intent(inout) :: dydt
300 integer,
intent(in) :: evolution_mode
330 if ((evolution_mode.eq.em_nse).and.(
rrate(i)%is_weak.eqv..false.)) cycle outer
349 select case (
rrate(i)%group)
352 if (
rrate(i)%parts(j) .eq. 0)
exit
353 dydt(
rrate(i)%parts(j)) = dydt(
rrate(i)%parts(j)) + rat * &
358 if (
rrate(i)%parts(j) .eq. 0)
exit
359 dydt(
rrate(i)%parts(j)) = dydt(
rrate(i)%parts(j)) + rat * &
360 rrate(i)%ch_amount(j) * y(
rrate(i)%parts(1)) * &
365 if (
rrate(i)%parts(j) .eq. 0)
exit
366 dydt(
rrate(i)%parts(j)) = dydt(
rrate(i)%parts(j)) + rat * &
367 rrate(i)%ch_amount(j) * y(
rrate(i)%parts(1)) * &
372 if (
rrate(i)%parts(j) .eq. 0)
exit
373 dydt(
rrate(i)%parts(j)) = dydt(
rrate(i)%parts(j)) + rat * &
374 rrate(i)%ch_amount(j) * y(
rrate(i)%parts(1)) * &
375 y(
rrate(i)%parts(2)) * y(
rrate(i)%parts(3)) * &
386 if ((rat .ge. infty) .or. (rat .ne. rat))
then
389 new_line(
"A")//
"was infinity or NaN.",&
395 if ((fissflag.ne.0) .and. (evolution_mode.ne.em_nse))
then
413 if ((rat .ge. infty) .or. (rat .ne. rat))
then
416 " was infinity or NaN.",
"abchange",280004)
418 if (y(
fissrate(i)%fissparts(1)) .eq. 0)
then
431 if ((rat .ge. infty) .or. (rat .ne. rat))
then
434 " was infinity or NaN.",
"abchange",280005)
468 real(
r_kind),
intent(in) :: rat
470 real(
r_kind),
dimension(net_size),
intent(in) :: y
473 integer :: zero_count
480 if ((rrate_in%reac_type .eq.
rrt_nf) .or. &
481 (rrate_in%reac_type .eq.
rrt_bf) .or. &
482 (rrate_in%reac_type .eq.
rrt_sf))
then
489 if (y(rrate_in%parts(i)) .le. 0) zero_count = zero_count+1
490 if (zero_count .gt. 1)
exit
492 if (zero_count .gt. 1) res = .true.
522 subroutine jacobi_init (time, itemp, rho, rkm, Y, Y_p, dYdt, rhs, h, evolution_mode)
532 real(
r_kind),
intent(in) :: time
533 real(
r_kind),
intent(in) :: itemp
534 real(
r_kind),
intent(in) :: rho
535 real(
r_kind),
intent(in) :: rkm
536 real(
r_kind),
dimension(:),
intent(in) :: y
537 real(
r_kind),
dimension(:),
intent(in) :: y_p
538 real(
r_kind),
dimension(:),
intent(inout) :: dydt
539 real(
r_kind),
dimension(:),
intent(inout) :: rhs
540 real(
r_kind),
intent(in) :: h
541 integer,
intent(in) :: evolution_mode
549 character*1,
dimension(6) :: descra
555 info_entry(
"jacobi_init")
561 descra = (/
"G",
" ",
" ",
"F",
" ",
" "/)
566 ye= sum(isotope(1:net_size)%p_nr*y(1:net_size))
574 elseif (ye.lt.0)
then
576 'jacobi_init',280007)
580 if(dlog10(rho*ye).ne.dlog10(rho*ye))
then
583 "Density is : "//trim(adjustl(
num_to_str(rho)))//
" [g/ccm]."//new_line(
"A")//&
585 'jacobi_init',280009)
589 if (itemp.le.freeze_rate_temp)
then
590 temp = freeze_rate_temp
592 print *,
"Freezing rates at "//
num_to_str(freeze_rate_temp)//
" GK"
606 vals(
dia(i)) = 1.d0/h
608 elseif(solver==1)
then
625 if ((evolution_mode.eq.em_nse).and.(rr_tmp%is_weak.eqv..false.)) cycle outer
634 select case (rr_tmp%group)
638 if (rr_tmp%parts(j) .eq. 0)
exit
639 dydt(rr_tmp%parts(j)) = dydt(rr_tmp%parts(j)) + rat * &
640 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1))
642 vals(rr_tmp%cscf_ind(1,j)) = vals(rr_tmp%cscf_ind(1,j)) - rat * &
648 if (rr_tmp%parts(j) .eq. 0)
exit
649 dydt(rr_tmp%parts(j)) = dydt(rr_tmp%parts(j)) + rat * &
650 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
653 vals(rr_tmp%cscf_ind(1,j)) = vals(rr_tmp%cscf_ind(1,j)) - rat * &
654 rr_tmp%ch_amount(j) * y(rr_tmp%parts(2))
655 vals(rr_tmp%cscf_ind(2,j)) = vals(rr_tmp%cscf_ind(2,j)) - rat * &
656 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1))
661 if (rr_tmp%parts(j) .eq. 0)
exit
662 dydt(rr_tmp%parts(j)) = dydt(rr_tmp%parts(j)) + rat * &
663 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
664 y(rr_tmp%parts(2)) * y(rr_tmp%parts(3))
666 vals(rr_tmp%cscf_ind(1,j)) = vals(rr_tmp%cscf_ind(1,j)) - rat * &
667 rr_tmp%ch_amount(j) * y(rr_tmp%parts(2)) * &
669 vals(rr_tmp%cscf_ind(2,j)) = vals(rr_tmp%cscf_ind(2,j)) - rat * &
670 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
672 vals(rr_tmp%cscf_ind(3,j)) = vals(rr_tmp%cscf_ind(3,j)) - rat * &
673 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
679 if (rr_tmp%parts(j) .eq. 0)
exit
680 dydt(rr_tmp%parts(j)) = dydt(rr_tmp%parts(j)) + rat * &
681 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
682 y(rr_tmp%parts(2)) * y(rr_tmp%parts(3)) * y(rr_tmp%parts(4))
684 vals(rr_tmp%cscf_ind(1,j)) = vals(rr_tmp%cscf_ind(1,j)) - rat * &
685 rr_tmp%ch_amount(j) * y(rr_tmp%parts(2)) * &
686 y(rr_tmp%parts(3)) * y(rr_tmp%parts(4))
687 vals(rr_tmp%cscf_ind(2,j)) = vals(rr_tmp%cscf_ind(2,j)) - rat * &
688 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
689 y(rr_tmp%parts(3)) * y(rr_tmp%parts(4))
690 vals(rr_tmp%cscf_ind(3,j)) = vals(rr_tmp%cscf_ind(3,j)) - rat * &
691 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
692 y(rr_tmp%parts(2)) * y(rr_tmp%parts(4))
693 vals(rr_tmp%cscf_ind(4,j)) = vals(rr_tmp%cscf_ind(4,j)) - rat * &
694 rr_tmp%ch_amount(j) * y(rr_tmp%parts(1)) * &
695 y(rr_tmp%parts(2)) * y(rr_tmp%parts(3))
701 ,
"jacobi_init",280009)
705 rrate(i)%cached = rat
710 if ((rat .gt. infty) .or. (rat .ne. rat))
then
714 new_line(
"A")//
"was infinity or NaN."//new_line(
"A")//&
715 "Logarithm of screening factor was "//
num_to_str(
hv(i)),&
716 "jacobi_init",280003)
720 new_line(
"A")//
"was infinity or NaN.",&
721 "jacobi_init",280003)
731 if ((fissflag.ne.0) .and. (evolution_mode.ne.em_nse))
then
732 fissloop:
do i=1,
nfiss
743 if (verbose_level .ge. 2)
then
744 if (
fissrate(i)%fissparts(1) .ne. ineu)
then
745 write(*,*)
"Fission rate is overflowing. Fissioning nucleus is "//&
746 trim(adjustl(isotope(
fissrate(i)%fissparts(1))%name))
748 write(*,*)
"Fission rate is overflowing. Fissioning nucleus is "//&
749 trim(adjustl(isotope(
fissrate(i)%fissparts(2))%name))
763 if ((y(
fissrate(i)%fissparts(1)) .le. 0) .and. (solver .eq. 0))
then
785 if ((y(
fissrate(i)%fissparts(1)) * y(
fissrate(i)%fissparts(2)) .le. 0) .and. (solver .eq. 0))
then
815 if (verbose_level .ge. 2)
then
817 if (dydt(i).ne.dydt(i))
then
818 call raise_exception(
"dYdt of "//trim(adjustl(net_names(i)))//
" is NaN.",&
819 "jacobi_init",280010)
823 rhs = dydt + d*(y_p - y)
825 call mkl_dcscmv(
'T',net_size,net_size,1.d0,descra,vals,
rows,
pt_b, &
827 elseif(solver==1)
then
830 dydt = dydt + d*(y_p-y)
832 info_exit(
"jacobi_init")
integer iwformat
defines format of the weak rates (0 = tabulated, 1 = log<ft>)
real(r_kind) function, dimension(net_size), public get_predictor_dydt()
Determines the predicted dydt_next.
Module to calculate electron screening.
real(r_kind), dimension(:), allocatable, public vals
contains non-zero matrix elements (of the Jacobian)
real(r_kind) freeze_rate_temp
Tmperature at which rates get frozen (for reacl. rates this should be 1d-2GK)
character *5, dimension(:), allocatable, public net_names
list of isotopes contained in the network
type(reactionrate_type), dimension(:), allocatable, public rrate
array containing all reaction rates used in the network
subroutine, public calculate_nu_rate(rrate, rat_calc)
Calculates the cross section times the neutrino flux.
Module that deals with reaclib reaction rates.
Takes care of the temperature and entropy updates due to the nuclear energy release.
real(r_kind), parameter, private rate_min_cutoff
Minimum cutoff for the reaction rates.
real(r_kind), private old_dens
character(50) function, public reaction_string(reac)
Return a string to represent a given reaction.
Module to deal with fission reactions.
real(r_kind), dimension(:), allocatable, public pf
partition functions
character(:) function, allocatable, public int_to_str(num)
Converts a given integer to a string.
integer, dimension(:), allocatable, public dia
dia(j) gives the index into vals that contains the (j,j) diagonal element of the matrix
subroutine, public raise_exception(msg, sub, error_code)
Raise a exception with a given error message.
logical function skip_rate(rat, rrate_in, Y)
Function to decide whether a rate can contribute or not.
real(r_kind), private old_temp
integer, public nfiss
Amount of fission rates.
integer, dimension(:), allocatable, public pt_e
pt_e(j)-pt_b(j) gives the index into vals that contains the LAST non-zero element of column j of the ...
logical use_timmes_mue
Use electron chemical potentials from timmes EOS for theoretical weak rates.
subroutine, public abchange(time, itemp, rho, Ye, rkm, Y, dYdt, evolution_mode)
This subroutine calculates the change in abundances (dYdt) due to reaction rates at given temperature...
real(r_kind), private old_time
type(isotope_type), dimension(:), allocatable, public isotope
all nuclides used in the network
real(r_kind), public mue
chemical potential of the electron
Module with helper functions such as rate tabulation, partition functions, and theoretical weak rates...
Contains variables and parameters related to neutrino fluxes.
integer, dimension(:,:), allocatable, public jind
jind(i,j) contains the index of jacobian entry in vals
subroutine, public calculate_detailed_balance(temp, rrate_array, inverse_rate, rate)
Calculate the inverse rate from the forward one via detailed balance.
logical, public weak
switch for weak rates
real(r_kind), private old_ye
Module that deals with inverse reaction rates.
real(r_kind) function, public get_l1()
The current l_1 (in Fortran: l(2))
integer, public net_size
total number of isotopes (network size)
Calculate the effective phase space integral needed for the log(ft) weak rates.
type(unit_type), public unit
constants and unit conversion factors
real(r_kind), parameter, private rate_max_cutoff
Maximum cutoff for the reaction rates.
real(r_kind) function, dimension(net_size), public get_predictor_y()
Determines the predicted y_next.
type(fissionrate_type), dimension(:), allocatable, public fissrate
Array storing fission reactions.
real(r_kind), private old_rad
character(:) function, allocatable, public num_to_str(num)
Converts a given real to a string with format "(1pE10.2)".
subroutine reset_qdot(time)
Reset the qdot variable.
logical, public tabulated
switch for tabulated rates
integer heating_mode
Mode for heating: 0 - no heating, 1 - heating using an entropy equation, 2 - heating from the energy ...
Contains subroutines for sparse matrix assembly and the solver core.
integer, dimension(:), allocatable, public rows
rows(i) is the number of the row in the matrix that contains the i-th value in vals
subroutine, public tabulated_index(temp)
Set tab_index for a given temperature.
subroutine chempot(temp, den, ye, etaele, etapos)
Given a temperature temp [K], density den [g/cm**3], and a composition characterized by y_e,...
real(r_kind), dimension(:), allocatable, public hv
Screening correction.
subroutine calculate_qdot(rrate, Y, h)
logical, private freeze_rate_indicator
Indicator for debug statement.
integer, parameter, public fiss_neglect
Amount of fission fragments not to be neglected.
subroutine, public jacobi_init(time, itemp, rho, rkm, Y, Y_p, dYdt, rhs, h, evolution_mode)
This subroutine calculates the entries in the jacobian and the right-hand side of the numerical integ...
integer function get_nr_reactants(group)
Get number of reactants of the reaction based on the reaclib chapter.
Contains subroutines to assemble the system Jacobian and changes in the abundances.
real(r_kind), dimension(9), public t9_pow
Powers of T, used for the rates.
subroutine, public screen(t9, rho, n, ye, mode)
This function calculates the screening coefficients hv.
This module contains everything for the tabulated rates that can replace reaclib rates.
integer, dimension(:), allocatable, public pt_b
pt_b(j) gives the index into vals that contains the FIRST non-zero element of column j of the matrix
subroutine, public inter_partf(temp, interpol)
Calculates partition function.
Contains types and objects shared between multiple modules.
subroutine, public calculate_tab_rate(rrate, temp, rat_calc)
Calculates the tabulated rate.
integer, public nreac
total number of reactions
subroutine calculate_reacl_rate(rrate, rat_calc)
Calculate a reaclib rate.
This module contains everything for the theoretical weak rates that are added into the rate array.
subroutine, public calc_t9_pow(t9)
A helper to compute powers of temperature used in the reaction rates.
subroutine, private calculate_reaction_rate(time, temp, rho, Ye, rkm, rrate_array, idx, rat_calc)
Subroutine to calculate the reaction rates.
integer nuflag
defines type of neutrino reactions used
logical, public iscreen
Flag whether screening is enabled or not.
gear_module contains adaptive high-order Gear solver
Contains all runtime parameters as well as phys and math constants.
subroutine, public calculate_twr_rate(rrate, temp, rho, Ye, rat_calc, nuloss)
Calculate the theoretical weak rate.
Subroutines needed to initialise the network.