diff --git a/README.rst b/README.rst index f4546ba..70bfdf5 100644 --- a/README.rst +++ b/README.rst @@ -10,7 +10,7 @@ More details about the MESA implementation can be found in `Section 10.2 of Jerm Installation ------------ -First ensure that you have installed `MESA `_ and have the +First ensure that you have installed `MESA `_ and have the environment variable ``MESA_DIR`` pointing to your MESA installation. Currently supported MESA versions: @@ -85,4 +85,10 @@ eos inlist The ``eos`` inlist is the same as the MESA EOS inlist and can contain anything that inlist specifies. It is the only inlist that allows for nesting of other inlists. -See `MESA's eos options `_ for the full set of supported options. +See `MESA's eos options `_ for the full set of supported options. + + +nuclear inlist +~~~~~~~~~~~~~~ + +The ``nuclear`` inliust contains many of the nuclear related inlist options from MESA's `star_job `_ and `controls `_ inlists. diff --git a/defaults/bbq.defaults b/defaults/bbq.defaults index c82a922..d3760d4 100644 --- a/defaults/bbq.defaults +++ b/defaults/bbq.defaults @@ -3,12 +3,6 @@ ! Nuclear network in MESA format. If path is not fully specified then we look in $MESA_DIR/data/net_data/nets net_name = 'approx21.net' - ! Screening mode - screening_mode = 'chugunov' - - ! Scale factor for weak reactions - weak_rate_factor = 1 - ! Solver tolerances max_steps = 1000000 eps = 1d-8 diff --git a/defaults/nuclear.defaults b/defaults/nuclear.defaults new file mode 100644 index 0000000..3e7dba9 --- /dev/null +++ b/defaults/nuclear.defaults @@ -0,0 +1,384 @@ + ! Users can also provide tabulated rates for any of the reactions. + ! Tabulated rates automatically take priority over any other options for the reaction. + ! e.g., if you provide a rate table for c12ag, those rates will be used + + ! To provide tabulated rates: + ! create a file of (T8, rate) pairs as in ``data/rates_data/rate_tables`` + ! You can give as many pairs as you want with any spacing in T8. + ! The first uncommented line of the file should be a number giving the + ! total number of (T8, rate) pairs in the subsequent lines. + ! The following lines are your specified values of T8 and rate separated + ! by a single space, one pair per line. + ! Add the filename to ``rate_list.txt`` along with the name of the rate you + ! want it to govern, either in ``data/rates_data/rate_tables`` or in a local + ! directory specified with the ``rate_tables_dir`` control. + ! Be aware that if you choose to put the modified ``rate_list.txt`` in + ! ``data/rates_data/rate_tables`` rather than a local directory, + ! your custom tabulated rate will override the rate for that reaction + ! for all future MESA runs. + + ! If the reaction you wish to control does not already have a + ! name that MESA will recognize, you will also need to add it to + ! the file specified by ``net_reaction_filename`` (defaults to reactions.list). + ! The default version of this file is located + ! in ``data/rates_data``. If you place a modified copy of this file + ! in your work directory, it will take precedence. + + + ! num_special_rate_factors + ! ~~~~~~~~~~~~~~~~~~~~~~~~ + ! reaction_for_special_factor + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! special_rate_factor + ! ~~~~~~~~~~~~~~~~~~~ + ! filename_of_special_rate + ! ~~~~~~~~~~~~~~~~~~~~~~~~ + + ! For using other special rate factors. + ! ``num_special_rate_factors`` must be <= ``max_num_special_rate_factors``. + + ! :: + + num_special_rate_factors = 0 + reaction_for_special_factor(:) = '' + special_rate_factor(:) = 1 + + ! If set, we read this filename from the local work directory, then try ``data/rates_data/rate_tables``. + ! This enables using the MESA provide custom rate tables with out changing the default for everyone. + ! Note we still multiple the loaded rate by special_rate_factor so leave that as 1 if you want the rate to be unchanged. + + ! For instance setting: + ! num_special_rate_factors = 1 + ! reaction_for_special_factor(1) = 'r_c12_ag_o16' + ! special_rate_factor(1) = 1 + ! filename_of_special_rate(1) = 'r_c12_ag_o16_kunz.txt' + + ! Recovers the old set_rate_c12ag = 'Kunz' option + + filename_of_special_rate(:) = '' + + + ! use_3a_fl87 + ! ~~~~~~~~~~~ + + ! If true then the triple alpha reaction is taken from Fushiki and Lamb, Apj, 317, 368-388, 1987 + ! This replaces the old set_rate_3a = 'Fl87' option + + use_3a_fl87 = .false. + + ! net_reaction_filename + ! ~~~~~~~~~~~~~~~~~~~~~ + + ! Looks first in current directory, then in ``mesa_data_dir/rates_data``. + + ! :: + + net_reaction_filename = 'reactions.list' + + + ! jina_reaclib_filename + ! ~~~~~~~~~~~~~~~~~~~~~ + + ! Empty string means use current standard version. + ! Which is defined in rates/public/rates_def.f90 as reaclib_filename + ! and is currently 'jina_reaclib_results_20171020_default' + + ! Else give name of file in directory ``mesa/data/rates_data``, + ! e.g., ``jina_reaclib_results_20130213default2`` + ! (which is an 18.8 MB file of rates data). + ! To use previous version, set to ``jina_reaclib_results_v2.2``. + + ! If you change reaclib version, you should clear the cache + ! after making the change in order to ensure that cached + ! rates from the default reaclib version are not being + ! read. (You can use the script ``empty_caches`` in + ! ``$MESA_DIR``.) + + ! In order to avoid this caching issue, one can also specify + ! a local rates cache directory via the control + ! ``rates_cache_dir``. + + ! :: + + jina_reaclib_filename = '' + + + ! jina_reaclib_min_T9 + ! ~~~~~~~~~~~~~~~~~~~ + + ! set jina reaclib rates to zero for T9 <= this. + ! if this control is <= 0, then use the standard default from rates. + ! need <= 3d-3 for pre-ms li7 burning + ! if change this, must remove old cached rates from data/rates_data/cache + + ! :: + + jina_reaclib_min_T9 = -1 + + + ! rate_tables_dir + ! ~~~~~~~~~~~~~~~ + + ! When MESA looks for the files ``rate_list.txt`` and ``weak_rate_list.txt``, + ! it will look in a local directory with this name first. + ! If doesn't find one, it will use the one in ``data/rates_data/rate_tables``. + + ! :: + + rate_tables_dir = 'rate_tables' + + + ! rate_cache_suffix + ! ~~~~~~~~~~~~~~~~~ + + ! If this not empty, then use it when creating names + ! for cache files for reaction rates from ``rate_tables_dir``. + ! If empty, the suffix will be '0'. + + ! :: + + rate_cache_suffix = '0' + + + ! T9_weaklib_full_off + ! ~~~~~~~~~~~~~~~~~~~ + ! T9_weaklib_full_on + ! ~~~~~~~~~~~~~~~~~~ + + ! Weak rates blend weaklib and reaclib according to temperature. + ! These can be used to overwrite the defaults in ``mesa/rates/public/rates_def`` + + ! + ``T9_weaklib_full_off`` : use pure reaclib for T <= this (ignore if <= 0) + ! + ``T9_weaklib_full_on`` : use pure weaklib for T >= this (ignore if <= 0) + + ! :: + + T9_weaklib_full_off = 0.01d0 + T9_weaklib_full_on = 0.02d0 + + + ! weaklib_blend_hi_Z + ! ~~~~~~~~~~~~~~~~~~ + + ! Ignore if <= 0. + ! Blend for intermediate temperatures. + ! For high Z elements, switch to reaclib at temp where no longer fully ionized. + ! As rough approximation for this, we switch at Fe to higher values of T9. + + ! :: + + weaklib_blend_hi_Z = 26 + + + ! T9_weaklib_full_off_hi_Z + ! ~~~~~~~~~~~~~~~~~~~~~~~~ + ! T9_weaklib_full_on_hi_Z + ! ~~~~~~~~~~~~~~~~~~~~~~~ + + ! If input element has Z >= ``weaklib_blend_hi_Z``, then use the following T9 limits: + + ! + ``T9_weaklib_full_off_hi_Z`` : use pure reaclib for T <= this (ignore if <= 0) + ! + ``T9_weaklib_full_on_hi_Z`` : use pure weaklib for T >= this (ignore if <= 0) + + ! :: + + T9_weaklib_full_off_hi_Z = 0.063d0 + T9_weaklib_full_on_hi_Z = 0.073d0 + + ! controls for other weak rate sources + ! ____________________________________ + + + ! use_suzuki_weak_rates + ! ~~~~~~~~~~~~~~~~~~~~~ + + ! If this is true, use the A=17-28 weak reaction rates from + + ! Suzuki, Toki, and Nomoto (2016) + ! Electron-capture and $\beta$-decay rates for sd-shell nuclei in stellar environments relevant to high-density O-Ne-Mg cores + ! http://adsabs.harvard.edu/abs/2016ApJ...817..163S + + ! If you make use of these rates, please cite the above paper. + + ! :: + + use_suzuki_weak_rates = .false. + + + ! use_special_weak_rates + ! ~~~~~~~~~~~~~~~~~~~~~~ + + ! If this is true, calculate special weak rates using the + ! approach described in Section 8 of Paxton et al. (2015). + + ! :: + + use_special_weak_rates = .false. + + + ! special_weak_states_file + ! ~~~~~~~~~~~~~~~~~~~~~~~~ + + ! File specifying which states to include + + ! Provide the low-lying energy levels of a given nucleus. + ! These are needed to calculate the partition function + ! and to indicate which states have allowed transitions. + ! Each isotope should have an entry of the form + + ! :: + + ! + ! + ! ... + ! + + ! where E = energy, J = spin. + + ! :: + + special_weak_states_file = 'special_weak_rates.states' + + + ! special_weak_transitions_file + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! File specifying to include + + ! These are the transitions for electron capture / beta decay + ! reactions that should be used. + + ! Each reaction should have and entry of the form + + ! :: + + ! + ! + ! ... + ! + + ! where si / sf are the n-th parent / daughter state, counting + ! in the order that you specified in the states file. logft is + ! the comparative half-life of that transition. + + ! :: + + special_weak_transitions_file = 'special_weak_rates.transitions' + + + ! ion_coulomb_corrections + ! ~~~~~~~~~~~~~~~~~~~~~~~ + + ! select which expression for the ion chemical potential to use + ! to calculate the energy shift associated with changing ion charge + + ! + 'none': no corrections + ! + 'DGC1973': Dewitt, Graboske, & Cooper, M. S. 1973, ApJ, 181, 439 + ! + 'I1993': Ichimaru, 1993, Reviews of Modern Physics, 65, 255 + ! + 'PCR2009': Potekhin, Chabrier, & Rogers, 2009, Phys. Rev. E, 79, 016411 + + ! :: + + ion_coulomb_corrections = 'none' + + + ! electron_coulomb_corrections + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! select which expression to use to calculate the shift in the + ! electron chemical potential at the location of the nucleus + + ! + 'none': no corrections + ! + 'ThomasFermi': Thomas-Fermi theory + ! + 'Itoh2002': Itoh et al., 2002, ApJ, 579, 380 + + ! :: + + electron_coulomb_corrections = 'none' + + + ! net_logTcut_lo + ! ~~~~~~~~~~~~~~ + + ! strong rates are zero ``logT < logTcut_lo`` + ! use default from net if this is <= 0 + + ! :: + + net_logTcut_lo = -1 + + + ! net_logTcut_lim + ! ~~~~~~~~~~~~~~~ + + ! strong rates cutoff smoothly for ``logT < logTcut_lim`` + ! use default from net if this is <= 0 + + ! :: + + net_logTcut_lim = -1 + + ! max_logT_for_net + ! ~~~~~~~~~~~~~~~~ + + ! :: + + max_logT_for_net = 10.2d0 + + + ! screening_mode + ! ~~~~~~~~~~~~~~ + + ! + empty string means no screening + ! + ``' extended'`` : + ! extends the Graboske (1973) method using results from Alastuey and Jancovici (1978), + ! along with plasma parameters from Itoh et al (1979) for strong screening. + ! + ``'salpeter'`` : + ! weak screening only. following Salpeter (1954), + ! with equations (4-215) and (4-221) of Clayton (1968). + ! + ``'chugunov'`` : + ! based on code from Sam Jones + ! Implements screening from Chugunov et al (2007) for weak and strong screening + ! MESA versions <=11435 used extended as the default value + + ! :: + + screening_mode = 'chugunov' + + ! fe56ec_fake_factor + ! ~~~~~~~~~~~~~~~~~~ + ! min_T_for_fe56ec_fake_factor + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Multiplier on ni56 electron capture rate to take isotopes in hardwired networks + ! to more neutron rich isotopes. + + ! :: + + fe56ec_fake_factor = 1d-7 + min_T_for_fe56ec_fake_factor = 3d9 + + + ! warn_rates_for_high_temp + ! ~~~~~~~~~~~~~~~~~~~~~~~~ + + ! If true then when any zone tries to evaluate a rate above ``max_safe_logT_for_rates`` + ! it generates a warning message. The code will cap the rate at the value + ! for ``max_safe_logT_for_rates`` whether the warning is on or not. + ! 10d0 is a sensible default for the max temperature, as that is where the partition tables + ! and polynomial fits to the rates are valid until. + ! warning messages include the text "rates have been truncated" and "WARNING: evaluating rates". + + ! :: + + warn_rates_for_high_temp = .true. + max_safe_logT_for_rates = 10d0 + + ! weak_rate_factor + ! ~~~~~~~~~~~~~~~~ + + ! all weak rates are multiplied by this factor + + ! :: + + weak_rate_factor = 1 diff --git a/defaults/rates.defaults b/defaults/rates.defaults deleted file mode 100644 index e69de29..0000000 diff --git a/inlist b/inlist index 9158df4..98f862b 100644 --- a/inlist +++ b/inlist @@ -26,4 +26,9 @@ &rates -/ \ No newline at end of file +/ + + +&nuclear + +/ diff --git a/src/ctrls.f90 b/src/ctrls.f90 index 9e88832..d32c10d 100644 --- a/src/ctrls.f90 +++ b/src/ctrls.f90 @@ -2,109 +2,70 @@ module ctrls use const_def use net_def use net_lib + use rates_def implicit none + integer, parameter :: max_num_special_rate_factors = 20 - type rate_t - integer :: eos_handle, net_handle - integer :: species, num_reactions - type (Net_Info) :: netinfo_target - type (Net_Info), pointer :: netinfo =>null() - type (Net_General_Info), pointer :: g =>null() + type state_t + include 'private/state.inc' + end type state_t - real(dp), pointer :: rate_factors(:) =>null() ! (num_reactions) - integer, pointer :: net_reaction_ptr(:) =>null() + type nuclear_t + include 'private/nuclear.inc' - integer, pointer :: reaction_id(:) - integer, dimension(:), pointer :: net_iso =>null(), chem_id =>null() - - real(dp), pointer :: ending_x(:) =>null() ! (species) - integer :: nfcn ! number of function evaluations - integer :: njac ! number of jacobian evaluations - integer :: nstep ! number of computed steps - integer :: naccpt ! number of accepted steps - integer :: nrejct ! number of rejected steps - - real(dp), dimension(:), pointer :: & - rate_raw =>null(), rate_raw_dT =>null(), rate_raw_dRho =>null(), & - rate_screened =>null(), rate_screened_dT =>null(), rate_screened_dRho =>null() - - integer :: screening_opt - - integer :: neut_id, prot_id - - end type rate_t + integer :: screening_opt + end type nuclear_t type bbq_t - character(len=strlen) :: net_name,screening_mode,iso_list_filename - integer :: max_steps - real(dp) :: weak_rate_factor, eps,odescal,stptry - character(len=strlen) :: inlist - logical :: use_input_file,use_random_sampling,use_profile - logical :: just_write_isos,write_iso_list - - type(rate_t) :: state + include 'private/bbq.inc' + + type(state_t) :: state + type(nuclear_t) :: nuclear end type bbq_t type profile_t - character(len=strlen) :: input_filename,output_filename,input_composition_filename - logical :: reflective_boundaries,write_comp_every_loop - integer :: num_loops + include 'private/profile.inc' end type profile_t type random_t - real(dp) :: log_time_min, log_time_max - real(dp) :: log_temp_min, log_temp_max - real(dp) :: log_rho_min, log_rho_max - real(dp) :: log_xa_min, log_xa_max - real(dp) :: neut_prot_limit_frac - character(len=strlen) :: output_starting_filename,output_ending_filename - integer :: num_samples, seed + include 'private/random.inc' end type random_t type sample_t - character(len=strlen) :: input_filename,output_filename - logical :: uniform_composition + include 'private/sample.inc' end type sample_t - + contains - subroutine load_bbq_inputs(inlist, options, ierr) - character(len=*), intent(in) :: inlist + subroutine load_bbq_inputs(inlist_in, options, ierr) + character(len=*), intent(in) :: inlist_in type(bbq_t), intent(inout) :: options integer, intent(out) :: ierr integer :: unit - - character(len=strlen) :: net_name, screening_mode,iso_list_filename - integer :: max_steps - real(dp) :: weak_rate_factor, eps,odescal,stptry - logical :: use_input_file,use_random_sampling,use_profile - logical :: just_write_isos,write_iso_list - integer :: flush_freq + include 'private/bbq.inc' - namelist /bbq/ net_name, screening_mode, weak_rate_factor, & + namelist /bbq/ net_name, & eps, odescal, stptry, max_steps, & use_input_file, use_random_sampling, use_profile,& iso_list_filename, just_write_isos, write_iso_list include '../defaults/bbq.defaults' - open(newunit=unit,file=inlist,status='old',action='read') + open(newunit=unit,file=inlist_in,status='old',action='read') read(unit,nml=bbq) close(unit) - options% inlist = inlist + options% inlist = inlist_in options% net_name = net_name - options% screening_mode = screening_mode options% iso_list_filename = iso_list_filename options% max_steps = max_steps - options% weak_rate_factor = weak_rate_factor options% eps = eps options% odescal = odescal options% stptry = stptry @@ -123,9 +84,7 @@ subroutine load_profile_inputs(inlist, options, ierr) integer, intent(out) :: ierr integer :: unit - character(len=strlen) :: input_filename,output_filename,input_composition_filename - logical :: reflective_boundaries,write_comp_every_loop - integer :: num_loops + include 'private/profile.inc' namelist /profile/ input_filename,output_filename, reflective_boundaries, num_loops,& input_composition_filename,write_comp_every_loop @@ -154,13 +113,7 @@ subroutine load_random_inputs(inlist, options, ierr) integer, intent(out) :: ierr integer :: unit - real(dp) :: log_time_min, log_time_max - real(dp) :: log_temp_min, log_temp_max - real(dp) :: log_rho_min, log_rho_max - real(dp) :: log_xa_min, log_xa_max - real(dp) :: neut_prot_limit_frac - character(len=strlen) :: output_starting_filename,output_ending_filename - integer :: num_samples, seed + include 'private/random.inc' namelist /random/ log_time_min, log_time_max, log_temp_min, log_temp_max, & @@ -200,8 +153,7 @@ subroutine load_sampling_inputs(inlist, options, ierr) integer, intent(out) :: ierr integer :: unit - character(len=strlen) :: input_filename,output_filename - logical :: uniform_composition + include 'private/sample.inc' namelist /sampling/ input_filename, output_filename, uniform_composition @@ -221,4 +173,87 @@ subroutine load_sampling_inputs(inlist, options, ierr) end subroutine load_sampling_inputs + subroutine load_nuclear_inputs(inlist, options, ierr) + character(len=*), intent(in) :: inlist + type(nuclear_t), intent(inout) :: options + integer, intent(out) :: ierr + integer :: unit + + include 'private/nuclear.inc' + + namelist /nuclear/ & + num_special_rate_factors,& + special_rate_factor,& + reaction_for_special_factor,& + filename_of_special_rate,& + weaklib_blend_hi_Z,& + T9_weaklib_full_off,& + T9_weaklib_full_on, & + T9_weaklib_full_off_hi_Z,& + T9_weaklib_full_on_hi_Z,& + use_suzuki_weak_rates,& + use_3a_fl87,& + use_special_weak_rates,& + special_weak_states_file,& + special_weak_transitions_file,& + ion_coulomb_corrections,& + electron_coulomb_corrections,& + net_reaction_filename,& + jina_reaclib_filename,& + jina_reaclib_min_T9,& + rate_tables_dir,& + rate_cache_suffix,& + screening_mode,& + net_logTcut_lo,& + net_logTcut_lim,& + max_logT_for_net,& + weak_rate_factor, & + fe56ec_fake_factor,& + min_T_for_fe56ec_fake_factor,& + warn_rates_for_high_temp,& + max_safe_logT_for_rates + + + include '../defaults/nuclear.defaults' + + open(newunit=unit,file=inlist,status='old',action='read') + read(unit,nml=nuclear) + close(unit) + + + options% num_special_rate_factors = num_special_rate_factors + options% special_rate_factor = special_rate_factor + options% reaction_for_special_factor = reaction_for_special_factor + options% filename_of_special_rate = filename_of_special_rate + options% weaklib_blend_hi_Z = weaklib_blend_hi_Z + options% T9_weaklib_full_off = T9_weaklib_full_off + options% T9_weaklib_full_on = T9_weaklib_full_on + options% T9_weaklib_full_off_hi_Z = T9_weaklib_full_off_hi_Z + options% T9_weaklib_full_on_hi_Z = T9_weaklib_full_on_hi_Z + options% use_suzuki_weak_rates = use_suzuki_weak_rates + options% use_3a_fl87 = use_3a_fl87 + options% use_special_weak_rates = use_special_weak_rates + options% special_weak_states_file = special_weak_states_file + options% special_weak_transitions_file = special_weak_transitions_file + options% ion_coulomb_corrections = ion_coulomb_corrections + options% electron_coulomb_corrections = electron_coulomb_corrections + options% net_reaction_filename = net_reaction_filename + options% jina_reaclib_filename = jina_reaclib_filename + options% jina_reaclib_min_T9 = jina_reaclib_min_T9 + options% rate_tables_dir = rate_tables_dir + options% rate_cache_suffix = rate_cache_suffix + options% screening_mode = screening_mode + options% net_logTcut_lo = net_logTcut_lo + options% net_logTcut_lim = net_logTcut_lim + options% max_logT_for_net = max_logT_for_net + options% weak_rate_factor = weak_rate_factor + options% fe56ec_fake_factor = fe56ec_fake_factor + options% min_T_for_fe56ec_fake_factor = min_T_for_fe56ec_fake_factor + options% warn_rates_for_high_temp = warn_rates_for_high_temp + options% max_safe_logT_for_rates = max_safe_logT_for_rates + + end subroutine load_nuclear_inputs + + + end module ctrls diff --git a/src/lib_bbq.f90 b/src/lib_bbq.f90 index 64e65ed..e50d621 100644 --- a/src/lib_bbq.f90 +++ b/src/lib_bbq.f90 @@ -52,6 +52,9 @@ subroutine read_inlist(bbq_in) call load_bbq_inputs(inlist, bbq_in, ierr) if(ierr/=0) return + call load_nuclear_inputs(bbq_in% inlist, bbq_in% nuclear, ierr) + if(ierr/=0) return + end subroutine read_inlist subroutine setup_eos(bbq_in) @@ -76,10 +79,13 @@ subroutine setup_eos(bbq_in) end subroutine setup_eos subroutine load_libs(bbq_in) - type(bbq_t) :: bbq_in + type(bbq_t),target :: bbq_in + type(nuclear_t), pointer :: nuc => null() integer :: ierr character (len=64) :: my_mesa_dir + nuc => bbq_in% nuclear + my_mesa_dir = '' call const_init(my_mesa_dir,ierr) if (ierr /= 0) then @@ -96,13 +102,17 @@ subroutine load_libs(bbq_in) call mesa_error(__FILE__,__LINE__) end if - call rates_init('reactions.list', '', 'rate_tables', .false., .false., & - '', '', '',ierr) + call rates_init(nuc% net_reaction_filename, nuc% jina_reaclib_filename,& + nuc% rate_tables_dir, & + nuc% use_suzuki_weak_rates, nuc% use_special_weak_rates, & + nuc% special_weak_states_file, nuc% special_weak_transitions_file,& + '',& + ierr) if (ierr /= 0) call mesa_error(__FILE__,__LINE__) - call rates_warning_init(.true., 10d0) + call rates_warning_init(nuc% warn_rates_for_high_temp ,nuc% max_safe_logT_for_rates) - bbq_in% state% screening_opt = screening_option(bbq_in% screening_mode, ierr) + nuc% screening_opt = screening_option(nuc% screening_mode, ierr) if (ierr /= 0) then write(*,*) 'read_inlist failed setting screening mode' call mesa_error(__FILE__,__LINE__) @@ -113,7 +123,8 @@ end subroutine load_libs subroutine net_setup(bbq_in) type(bbq_t),target :: bbq_in integer :: ierr, j - type(rate_t), pointer :: s => null() + type(state_t), pointer :: s => null() + type(nuclear_t), pointer :: nuc => null() call read_inlist(bbq_in) @@ -122,6 +133,7 @@ subroutine net_setup(bbq_in) call setup_eos(bbq_in) s => bbq_in% state + nuc => bbq_in% nuclear allocate(s% net_iso(num_chem_isos), & s% chem_id(num_chem_isos)) @@ -163,12 +175,6 @@ subroutine net_setup(bbq_in) s% num_reactions = s% g% num_reactions allocate(s% reaction_id(s% num_reactions)) - - call net_setup_tables(s% net_handle, cache_suffix, ierr) - if (ierr /= 0) then - write(*,*) 'net_setup_tables failed' - call mesa_error(__FILE__,__LINE__) - end if call get_chem_id_table(s% net_handle, s% species, s% chem_id, ierr) if (ierr /= 0) then @@ -189,12 +195,6 @@ subroutine net_setup(bbq_in) call mesa_error(__FILE__,__LINE__) end if - call net_set_fe56ec_fake_factor(s% net_handle, 1d-7,3d9, ierr) - if (ierr /= 0) then - write(*,*) 'net_set_fe56ec_fake_factor failed' - call mesa_error(__FILE__,__LINE__) - end if - s% netinfo => s% netinfo_target allocate( & @@ -206,13 +206,43 @@ subroutine net_setup(bbq_in) call mesa_error(__FILE__,__LINE__) end if - s% rate_factors(:) = 1 + call set_rate_factors(s, nuc, ierr) + if (ierr /= 0) then + write(*,*) 'failed in set_rate_factors' + call mesa_error(__FILE__,__LINE__) + end if + + call set_global_nuc_opts(s, nuc) + call get_net_reaction_table_ptr(s% net_handle, & s% net_reaction_ptr, ierr) if (ierr /= 0) then write(*,*) 'bad net? get_net_reaction_table_ptr failed' return + end if + + call set_rate_factors(s, nuc, ierr) + if (ierr /= 0) then + write(*,*) 'failed in set_rate_factors' + call mesa_error(__FILE__,__LINE__) + end if + + call set_global_nuc_opts(s, nuc) + + call net_set_fe56ec_fake_factor(s% net_handle,& + nuc% fe56ec_fake_factor,& + nuc% min_T_for_fe56ec_fake_factor,& + ierr) + if (ierr /= 0) then + write(*,*) 'net_set_fe56ec_fake_factor failed' + call mesa_error(__FILE__,__LINE__) + end if + + call net_setup_tables(s% net_handle, nuc% rate_cache_suffix, ierr) + if (ierr /= 0) then + write(*,*) 'net_setup_tables failed' + call mesa_error(__FILE__,__LINE__) end if s% neut_id = -1 @@ -225,10 +255,96 @@ subroutine net_setup(bbq_in) end subroutine net_setup + subroutine set_rate_factors(state, nuc, ierr) + type(state_t) :: state + type(nuclear_t) :: nuc + integer, intent(out) :: ierr + integer :: j, i, ir + logical :: error + + state% rate_factors(:) = 1 + if (nuc% num_special_rate_factors <= 0) return + + do i=1,nuc% num_special_rate_factors + if (len_trim(nuc% reaction_for_special_factor(i)) == 0) cycle + ir = rates_reaction_id(nuc% reaction_for_special_factor(i)) + j = 0 + if (ir > 0) j = state% net_reaction_ptr(ir) + if (j <= 0) then + write(*,*) 'Failed to find reaction_for_special_factor ' // & + trim(nuc% reaction_for_special_factor(i)), & + j, nuc% special_rate_factor(i) + error = .true. + cycle + end if + state% rate_factors(j) = nuc% special_rate_factor(i) + write(*,*) 'set special rate factor for ' // & + trim(nuc% reaction_for_special_factor(i)), & + j, nuc% special_rate_factor(i) + end do + + if(error) call mesa_error(__FILE__,__LINE__) + + + call read_rates_from_files(nuc% reaction_for_special_factor, nuc% filename_of_special_rate, ierr) + if (ierr /= 0) then + write(*,*) 'failed in read_rates_from_files' + return + end if + + end subroutine set_rate_factors + + + subroutine set_global_nuc_opts(state, nuc) + use rates_def + type(state_t) :: state + type(nuclear_t) :: nuc + integer :: ierr + + + if (abs(nuc% T9_weaklib_full_off - T9_weaklib_full_off) > 1d-6) then + write(*,*) 'set T9_weaklib_full_off', nuc% T9_weaklib_full_off + T9_weaklib_full_off = nuc% T9_weaklib_full_off + end if + + if (abs(nuc% T9_weaklib_full_on - T9_weaklib_full_on) > 1d-6) then + write(*,*) 'set T9_weaklib_full_on', nuc% T9_weaklib_full_on + T9_weaklib_full_on = nuc% T9_weaklib_full_on + end if + + if (nuc% weaklib_blend_hi_Z /= weaklib_blend_hi_Z) then + write(*,*) 'set weaklib_blend_hi_Z', nuc% weaklib_blend_hi_Z + weaklib_blend_hi_Z = nuc% weaklib_blend_hi_Z + end if + + if (abs(nuc% T9_weaklib_full_off_hi_Z - T9_weaklib_full_off_hi_Z) > 1d-6) then + write(*,*) 'set T9_weaklib_full_off_hi_Z', nuc% T9_weaklib_full_off_hi_Z + T9_weaklib_full_off_hi_Z = nuc% T9_weaklib_full_off_hi_Z + end if + + if (abs(nuc% T9_weaklib_full_on_hi_Z - T9_weaklib_full_on_hi_Z) > 1d-6) then + write(*,*) 'set T9_weaklib_full_on_hi_Z', nuc% T9_weaklib_full_on_hi_Z + T9_weaklib_full_on_hi_Z = nuc% T9_weaklib_full_on_hi_Z + end if + + ! set up coulomb corrections for the special weak rates + which_mui_coulomb = get_mui_value(nuc% ion_coulomb_corrections) + which_vs_coulomb = get_vs_value(nuc% electron_coulomb_corrections) + + + call net_set_logTcut(state% net_handle, nuc% net_logTcut_lo, nuc% net_logTcut_lim, ierr) + if (ierr /= 0) return + + end subroutine set_global_nuc_opts + + + + subroutine do_burn(in, out, bbq_in, ierr) type(inputs_t),intent(in) :: in type(outputs_t),intent(out) :: out - type(bbq_t) :: bbq_in + type(bbq_t),target :: bbq_in + type(nuclear_t), pointer :: nuc=>null() integer, intent(out) :: ierr real(dp) :: eta, logRho, logT, Rho, T, xsum, eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, avg_eps_nuc @@ -250,6 +366,8 @@ subroutine do_burn(in, out, bbq_in, ierr) log10Ts_f =>null(), log10Rhos_f =>null(), etas_f =>null(), log10Ps_f =>null() + nuc => bbq_in% nuclear + logT = in% logT T = exp10(logT) logRho = in% logRho @@ -299,9 +417,9 @@ subroutine do_burn(in, out, bbq_in, ierr) bbq_in% state% net_handle, bbq_in% state% eos_handle, bbq_in% state% species, & bbq_in% state% num_reactions, 0d0, times(1), in% xa, & num_times, times, log10Ts_f1, log10Rhos_f1, etas_f1, dxdt_source_term, & - bbq_in% state% rate_factors, bbq_in% weak_rate_factor, & + bbq_in% state% rate_factors, nuc% weak_rate_factor, & std_reaction_Qs, std_reaction_neuQs, & - bbq_in% state% screening_opt, & + nuc% screening_opt, & bbq_in% stptry, bbq_in% max_steps, bbq_in% eps, bbq_in% odescal, & .true., .false., burn_dbg, burn_finish_substep, & bbq_in% state% ending_x, out% eps_nuc_categories, & diff --git a/src/private/bbq.inc b/src/private/bbq.inc new file mode 100644 index 0000000..9b3e477 --- /dev/null +++ b/src/private/bbq.inc @@ -0,0 +1,6 @@ +character(len=strlen) :: net_name,iso_list_filename +integer :: max_steps +real(dp) :: eps,odescal,stptry +character(len=strlen) :: inlist +logical :: use_input_file,use_random_sampling,use_profile +logical :: just_write_isos,write_iso_list diff --git a/src/private/nuclear.inc b/src/private/nuclear.inc new file mode 100644 index 0000000..7ba2bcc --- /dev/null +++ b/src/private/nuclear.inc @@ -0,0 +1,31 @@ +integer :: num_special_rate_factors +real(dp) :: special_rate_factor(max_num_special_rate_factors) +character(len=maxlen_reaction_Name) :: & + reaction_for_special_factor(max_num_special_rate_factors) +character(len=strlen) :: filename_of_special_rate(max_num_special_rate_factors) +integer :: weaklib_blend_hi_Z +real(dp) :: & + T9_weaklib_full_off, T9_weaklib_full_on, & + T9_weaklib_full_off_hi_Z, T9_weaklib_full_on_hi_Z + +logical :: use_suzuki_weak_rates, use_3a_fl87 + +logical :: use_special_weak_rates +character (len=1000) :: special_weak_states_file, special_weak_transitions_file +character (len=strlen) :: ion_coulomb_corrections, electron_coulomb_corrections + +character (len=strlen) :: net_reaction_filename, jina_reaclib_filename +real(dp) :: jina_reaclib_min_T9 + +character (len=strlen) :: rate_tables_dir, rate_cache_suffix + +character (len=strlen) :: screening_mode +real(dp) :: net_logTcut_lo, net_logTcut_lim + +real(dp) :: max_logT_for_net + +real(dp) :: weak_rate_factor, & + fe56ec_fake_factor, min_T_for_fe56ec_fake_factor + +logical :: warn_rates_for_high_temp +real(dp) :: max_safe_logT_for_rates diff --git a/src/private/profile.inc b/src/private/profile.inc new file mode 100644 index 0000000..dda3e98 --- /dev/null +++ b/src/private/profile.inc @@ -0,0 +1,3 @@ +character(len=strlen) :: input_filename,output_filename,input_composition_filename +logical :: reflective_boundaries,write_comp_every_loop +integer :: num_loops diff --git a/src/private/random.inc b/src/private/random.inc new file mode 100644 index 0000000..5db1a2e --- /dev/null +++ b/src/private/random.inc @@ -0,0 +1,7 @@ +real(dp) :: log_time_min, log_time_max +real(dp) :: log_temp_min, log_temp_max +real(dp) :: log_rho_min, log_rho_max +real(dp) :: log_xa_min, log_xa_max +real(dp) :: neut_prot_limit_frac +character(len=strlen) :: output_starting_filename,output_ending_filename +integer :: num_samples, seed diff --git a/src/private/sample.inc b/src/private/sample.inc new file mode 100644 index 0000000..cca9803 --- /dev/null +++ b/src/private/sample.inc @@ -0,0 +1,2 @@ +character(len=strlen) :: input_filename,output_filename +logical :: uniform_composition diff --git a/src/private/state.inc b/src/private/state.inc new file mode 100644 index 0000000..a67c9fa --- /dev/null +++ b/src/private/state.inc @@ -0,0 +1,25 @@ +integer :: eos_handle, net_handle +integer :: species, num_reactions + +type (Net_Info) :: netinfo_target +type (Net_Info), pointer :: netinfo =>null() +type (Net_General_Info), pointer :: g =>null() + +real(dp), pointer :: rate_factors(:) =>null() ! (num_reactions) +integer, pointer :: net_reaction_ptr(:) =>null() + +integer, pointer :: reaction_id(:) +integer, dimension(:), pointer :: net_iso =>null(), chem_id =>null() + +real(dp), pointer :: ending_x(:) =>null() ! (species) +integer :: nfcn ! number of function evaluations +integer :: njac ! number of jacobian evaluations +integer :: nstep ! number of computed steps +integer :: naccpt ! number of accepted steps +integer :: nrejct ! number of rejected steps + +real(dp), dimension(:), pointer :: & + rate_raw =>null(), rate_raw_dT =>null(), rate_raw_dRho =>null(), & + rate_screened =>null(), rate_screened_dT =>null(), rate_screened_dRho =>null() + +integer :: neut_id, prot_id diff --git a/tests/big_net/inlist b/tests/big_net/inlist index dd3d8d2..6be0032 100644 --- a/tests/big_net/inlist +++ b/tests/big_net/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'mesa_495.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -40,3 +38,7 @@ &eos / + +&nuclear + +/ diff --git a/tests/big_net/test b/tests/big_net/test index 6af6286..82822b3 100755 --- a/tests/big_net/test +++ b/tests/big_net/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/isos/inlist b/tests/isos/inlist index c2f75b8..734d2e2 100644 --- a/tests/isos/inlist +++ b/tests/isos/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'approx21.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! What mode to run just_write_isos=.true. @@ -28,3 +26,7 @@ &eos / + +&nuclear + +/ diff --git a/tests/isos/test b/tests/isos/test index 6af6286..82822b3 100755 --- a/tests/isos/test +++ b/tests/isos/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/ni56/inlist b/tests/ni56/inlist index bf38a92..29449dd 100644 --- a/tests/ni56/inlist +++ b/tests/ni56/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'ni56.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -38,3 +36,7 @@ &eos / + +&nuclear + +/ diff --git a/tests/ni56/test b/tests/ni56/test index 6af6286..82822b3 100755 --- a/tests/ni56/test +++ b/tests/ni56/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/profile/inlist b/tests/profile/inlist index 3abed52..383f58e 100644 --- a/tests/profile/inlist +++ b/tests/profile/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'approx21.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -42,3 +40,7 @@ &eos / + +&nuclear + +/ diff --git a/tests/profile/test b/tests/profile/test index 6af6286..82822b3 100755 --- a/tests/profile/test +++ b/tests/profile/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/random/inlist b/tests/random/inlist index 526187c..459fa1f 100644 --- a/tests/random/inlist +++ b/tests/random/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'approx21.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -54,17 +52,15 @@ &profile ! For use_profile - input_filename = '' - input_composition_filename = '' - output_filename = '' - reflective_boundaries=.true. - num_loops = 1 +/ + +&eos / -&eos +&nuclear / diff --git a/tests/random/test b/tests/random/test index 6af6286..82822b3 100755 --- a/tests/random/test +++ b/tests/random/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/sample/inlist b/tests/sample/inlist index 88aa3a0..6620e58 100644 --- a/tests/sample/inlist +++ b/tests/sample/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'approx21.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -36,3 +34,8 @@ &eos / + + +&nuclear + +/ diff --git a/tests/sample/test b/tests/sample/test index 6af6286..82822b3 100755 --- a/tests/sample/test +++ b/tests/sample/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck diff --git a/tests/sample_test_suite/inlist b/tests/sample_test_suite/inlist index a4908d7..37a344b 100644 --- a/tests/sample_test_suite/inlist +++ b/tests/sample_test_suite/inlist @@ -2,8 +2,6 @@ ! Physcis options net_name = 'approx21.net' - screening_mode = 'chugunov' - weak_rate_factor = 1 ! Solver tolerances max_steps = 1000000 @@ -33,3 +31,7 @@ &eos / + +&nuclear + +/ diff --git a/tests/sample_test_suite/test b/tests/sample_test_suite/test index 6af6286..82822b3 100755 --- a/tests/sample_test_suite/test +++ b/tests/sample_test_suite/test @@ -1,4 +1,6 @@ #!/bin/bash -../../bbq inlist +if ! ../../bbq inlist;then + exit 1 +fi ./ck