From 8aa9d9d7f9e62c42bc4849ce8d5f87a249f60c11 Mon Sep 17 00:00:00 2001 From: EbF Date: Wed, 27 Dec 2023 12:01:20 -0500 Subject: [PATCH 1/3] full test, fix checksum: reduce photo interval --- star/test_suite/make_pre_ccsn_13bvn/inlist_to_cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/star/test_suite/make_pre_ccsn_13bvn/inlist_to_cc b/star/test_suite/make_pre_ccsn_13bvn/inlist_to_cc index 2676ca09b..27b0006af 100644 --- a/star/test_suite/make_pre_ccsn_13bvn/inlist_to_cc +++ b/star/test_suite/make_pre_ccsn_13bvn/inlist_to_cc @@ -73,7 +73,7 @@ ! FOR DEBUGGING - !photo_interval = 10 + photo_interval = 50 !profile_interval = 1 !history_interval = 1 terminal_interval = 1 From 0938e547946e757cd38770993c97474162770127 Mon Sep 17 00:00:00 2001 From: EbF Date: Wed, 27 Dec 2023 12:54:43 -0500 Subject: [PATCH 2/3] [ci skip] add v_drag and test through ppisn pulses --- .../dev_TDC_through_ppisn/README.rst | 21 + .../dev_TDC_through_ppisn/README_DEV | 7 + .../dev_TDC_through_ppisn/ck | 7 + .../dev_TDC_through_ppisn/clean | 4 + .../dev_TDC_through_ppisn/he_dep.mod | 3 + .../history_columns.list | 1073 ++++++++++++ .../inlist_after_first_pulse | 17 + .../dev_TDC_through_ppisn/inlist_extra | 14 + .../dev_TDC_through_ppisn/inlist_hydro_off | 26 + .../dev_TDC_through_ppisn/inlist_hydro_on | 36 + .../dev_TDC_through_ppisn/inlist_pgstar | 210 +++ .../dev_TDC_through_ppisn/inlist_ppisn | 395 +++++ .../dev_TDC_through_ppisn/inlist_pulses | 62 + .../inlist_pulses_header | 41 + .../dev_TDC_through_ppisn/inlist_to_he_dep | 47 + .../inlist_to_he_dep_header | 41 + .../dev_TDC_through_ppisn/make/makefile | 5 + .../dev_TDC_through_ppisn/mk | 13 + .../profile_columns.list | 962 +++++++++++ .../dev_TDC_through_ppisn/re | 33 + .../dev_TDC_through_ppisn/rn | 25 + .../dev_TDC_through_ppisn/rn1 | 7 + .../dev_TDC_through_ppisn/src/run.f90 | 15 + .../src/run_star_extras.f90 | 1498 +++++++++++++++++ .../dev_TDC_through_ppisn/standard_he_dep.mod | 3 + .../dev_TDC_through_ppisn/zams.mod | 3 + star/private/ctrls_io.f90 | 10 + star/private/hydro_riemann.f90 | 22 + star_data/private/star_controls.inc | 6 +- 29 files changed, 4605 insertions(+), 1 deletion(-) create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/README.rst create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/README_DEV create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/ck create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/clean create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/he_dep.mod create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/history_columns.list create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_after_first_pulse create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_extra create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_off create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_on create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pgstar create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_ppisn create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses_header create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep_header create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/make/makefile create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/mk create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/profile_columns.list create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/re create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn create mode 100755 star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn1 create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run.f90 create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run_star_extras.f90 create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/standard_he_dep.mod create mode 100644 star/dev_cases_test_TDC/dev_TDC_through_ppisn/zams.mod diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README.rst b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README.rst new file mode 100644 index 000000000..a25931921 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README.rst @@ -0,0 +1,21 @@ +.. _ppisn: + +***** +ppisn +***** + +This test case evolves a very massive helium star from the He-ZAMS +up to the ocurrence of a pulsational pair-instability event (see |Marchant2019|). + +.. |Marchant2019| replace:: `Marchant et al. 2019 `__ + +Initialization of the model +=========================== +The initial mass of the helium star is set in ``inlist_extra`` + +.. literalinclude:: ../../../star/test_suite/ppisn/inlist_extra + +In this case we use a :math:`72 M_\odot` + +Last-Updated: 2019-11-12 (mesa r12413) by Pablo Marchant + diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README_DEV b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README_DEV new file mode 100644 index 000000000..64734bd30 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/README_DEV @@ -0,0 +1,7 @@ +Test cases with names beginning in "DEV_" are under development. +This means that they are not ready for consideration by general users as a "how to" guide for research to be published. + +Play but don't publish! + +The situation is similar to that for options in controls_dev.defaults and star_job_dev.defaults. + diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/ck b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/ck new file mode 100755 index 000000000..ac08f15c6 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/ck @@ -0,0 +1,7 @@ +#!/bin/bash + +# this provides the definition of check_one +# check_one +source "${MESA_DIR}/star/test_suite/test_suite_helpers" + +check_one diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/clean b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/clean new file mode 100755 index 000000000..95545a5c1 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/clean @@ -0,0 +1,4 @@ +#!/bin/bash + +cd make +make clean diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/he_dep.mod b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/he_dep.mod new file mode 100644 index 000000000..0058b23ab --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/he_dep.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d44651e90a73db4a753274968e3e2649ed1cecb47fbdd72a0d915a3514ef6329 +size 1243628 diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/history_columns.list b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/history_columns.list new file mode 100644 index 000000000..03687af8c --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/history_columns.list @@ -0,0 +1,1073 @@ +! history_columns.list -- determines the contents of star history logs +! you can use a non-standard version by setting history_columns_file in your inlist + +! units are cgs unless otherwise noted. + +! reorder the following names as desired to reorder columns. +! comment out the name to omit a column (fewer columns => less IO => faster running). +! remove '!' to restore a column. + +! if you have a situation where you want a non-standard set of columns, +! make a copy of this file, edit as desired, and give the new filename in your inlist +! as history_columns_file. if you are just adding columns, you can 'include' this file, +! and just list the additions in your file. note: to include the standard default +! version, use include '' -- the 0 length string means include the default file. + +! blank lines and comments can be used freely. +! if a column name appears more than once in the list, only the first occurrence is used. + +! if you need to have something added to the list of options, let me know.... + + +! the first few lines of the log file contain a few items: + + ! version_number -- for the version of mesa being used + ! burn_min1 -- 1st limit for reported burning, in erg/g/s + ! burn_min2 -- 2nd limit for reported burning, in erg/g/s + + +!# other files + +! note: you can include another list by doing +! include 'filename' +! include '' means include the default standard list file + +! the following lines of the log file contain info about 1 model per row + +!---------------------------------------------------------------------------------------------- + +!# general info about the model + + model_number ! counting from the start of the run + num_zones ! number of zones in the model + + !## age + + star_age ! elapsed simulated time in years since the start of the run + !star_age_sec ! elapsed simulated time in seconds since the start of the run + !star_age_min ! elapsed simulated time in minutes since the start of the run + !star_age_hr ! elapsed simulated time in hours since the start of the run + !star_age_day ! elapsed simulated time in days since the start of the run + !day ! elapsed simulated time in days since the start of the run + + !log_star_age + !log_star_age_sec + + !## timestep + + !time_step ! timestep in years since previous model + !time_step_sec ! timestep in seconds since previous model + !time_step_days + log_dt ! log10 time_step in years + !log_dt_sec ! log10 time_step in seconds + !log_dt_days ! log10 time_step in days + + !## mass + + star_mass ! in Msun units + !log_star_mass + + !star_gravitational_mass ! star_mass is baryonic mass + !star_mass_grav_div_mass + + !delta_mass ! star_mass - initial_mass in Msun units + log_xmstar ! log10 mass exterior to M_center (grams) + + !## mass change + + star_mdot ! d(star_mass)/dt (in msolar per year) + log_abs_mdot ! log10(abs(star_mdot)) (in msolar per year) + + !## imposed surface conditions + !Tsurf_factor + !tau_factor + !tau_surface + + !## imposed center conditions + !m_center + !m_center_gm + !r_center + !r_center_cm + !r_center_km + !L_center + !log_L_center + !log_L_center_ergs_s + !v_center + !v_center_kms + + !logt_max + +!---------------------------------------------------------------------------------------------- + +!# mixing and convection + + !max_conv_vel_div_csound + !max_gradT_div_grada + !max_gradT_sub_grada + !min_log_mlt_Gamma + + + !## mixing regions + + mass_conv_core ! (Msun) mass coord of top of convective core. 0 if core is not convective + + ! mx1 refers to the largest (by mass) convective region. + ! mx2 is the 2nd largest. + + ! conv_mx1_top and conv_mx1_bot are the region where mixing_type == convective_mixing. + ! mx1_top and mx1_bot are the extent of all kinds of mixing, convective and other. + + ! values are m/Mstar + conv_mx1_top + conv_mx1_bot + conv_mx2_top + conv_mx2_bot + mx1_top + mx1_bot + mx2_top + mx2_bot + + ! radius -- values are radii in Rsun units + !conv_mx1_top_r + !conv_mx1_bot_r + !conv_mx2_top_r + !conv_mx2_bot_r + !mx1_top_r + !mx1_bot_r + !mx2_top_r + !mx2_bot_r + + ! you might want to get a more complete list of mixing regions by using the following + + !mixing_regions ! note: this includes regions where the mixing type is no_mixing. + + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives the mixing type as defined in const/public/const_def.f90. + + ! the second column for a region gives the m/mstar location of the top of the region + ! entries for extra columns after the last region in the star will have an invalid mixing_type value of -1. + ! mstar is the total mass of the star, so these locations range from 0 to 1 + ! all regions are include starting from the center, so the bottom of one region + ! is the top of the previous one. since we start at the center, the bottom of the 1st region is 0. + + ! the columns in the log file will have names like 'mix_type_1' and 'mix_qtop_1' + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + + + !mix_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'mix_relr_type_1' and 'mix_relr_top_1' + + + !## conditions at base of largest convection zone (by mass) + !cz_bot_mass ! mass coordinate of base (Msun) + !cz_mass ! mass coordinate of base (Msun) -- same as cz_bot_mass + !cz_log_xmass ! mass exterior to base (g) + !cz_log_xmsun ! mass exterior to base (Msun) + !cz_xm ! mass exterior to base (Msun) + !cz_logT + !cz_logRho + !cz_logP + !cz_bot_radius ! Rsun + !cz_log_column_depth + !cz_log_radial_depth + !cz_luminosity ! Lsun + !cz_opacity + !cz_log_tau + !cz_eta + !cz_log_eps_nuc ! log10(ergs/g/s) + !cz_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_csound + !cz_scale_height + !cz_grav + + !cz_omega + !cz_omega_div_omega_crit + + !cz_zone + + ! mass fractions at base of largest convection zone (by mass) + !cz_log_xa h1 + !cz_log_xa he4 + + !## conditions at top of largest convection zone (by mass) + !cz_top_mass ! mass coordinate of top (Msun) + !cz_top_log_xmass ! mass exterior to top (g) + !cz_top_log_xmsun ! mass exterior to top (Msun) + !cz_top_xm ! mass exterior to top (Msun) + !cz_top_logT + !cz_top_logRho + !cz_top_logP + !cz_top_radius ! Rsun + !cz_top_log_column_depth + !cz_top_log_radial_depth + !cz_top_luminosity ! Lsun + !cz_top_opacity + !cz_top_log_tau + !cz_top_eta + !cz_top_log_eps_nuc ! log10(ergs/g/s) + !cz_top_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_top_csound + !cz_top_scale_height + !cz_top_grav + + !cz_top_omega + !cz_top_omega_div_omega_crit + + !cz_top_zone + !cz_top_zone_logdq + + ! mass fractions at top of largest convection zone (by mass) + !cz_top_log_xa h1 + !cz_top_log_xa he4 + +!---------------------------------------------------------------------------------------------- + +!# nuclear reactions + + !## integrated quantities + + !power_h_burn ! total thermal power from PP and CNO, excluding neutrinos (in Lsun units) + !power_he_burn ! total thermal power from triple-alpha, excluding neutrinos (in Lsun units) + !power_photo + !power_z_burn + !log_power_nuc_burn ! total thermal power from all burning, excluding photodisintegrations + log_LH ! log10 power_h_burn + log_LHe ! log10 power_he_burn + log_LZ ! log10 total burning power including LC, but excluding LH and LHe and photodisintegrations + log_Lnuc ! log(LH + LHe + LZ) + !log_Lnuc_ergs_s + !log_Lnuc_sub_log_L + !lnuc_photo + + !extra_L ! integral of extra_heat in Lsun units + !log_extra_L ! log10 extra_L + + !## neutrino losses + log_Lneu ! log10 power emitted in neutrinos, nuclear and thermal (in Lsun units) + !log_Lneu_nuc ! log10 power emitted in neutrinos, nuclear sources only (in Lsun units) + !log_Lneu_nonnuc ! log10 power emitted in neutrinos, thermal sources only (in Lsun units) + + !mass_loc_of_max_eps_nuc ! (in Msun units) + !mass_ext_to_max_eps_nuc ! (in Msun units) + !eps_grav_integral ! (in Lsun units) + !log_abs_Lgrav ! log10 abs(eps_grav_integral) (in Lsun units) + + !## information about reactions (by category) + + ! log10 total luminosity for reaction categories (Lsun units) + + pp + cno + tri_alpha + !c_alpha + !n_alpha + !o_alpha + !ne_alpha + !na_alpha + !mg_alpha + !si_alpha + !s_alpha + !ar_alpha + !ca_alpha + !ti_alpha + !fe_co_ni + !c12_c12 + !c12_o16 + !o16_o16 + !photo + !pnhe4 + !other + + !## information about individual reactions + + ! adds columns for all of the reactions that are in the current net + ! Note that if using op_split_burn=.true. then zones which have been split will report 0 for thier rates + !add_raw_rates ! raw reaction rates, reactions/second + !add_screened_rates ! screened reaction rates reactions/second + !add_eps_nuc_rates ! Nuclear energy (minus neutrino losses) released erg/s + !add_eps_neu_rates ! Neutrino losses erg/s + + ! individual reactions (as many as desired) + ! use list_net_reactions = .true. in star_job to list all reactions in the current net + ! reactions/second + !raw_rate r_h1_h1_ec_h2 + !raw_rate r_h1_h1_wk_h2 + + + + !## nuclear reactions at center + + ! center log10 burn erg/g/s for reaction categories + + !c_log_eps_burn cno + !c_log_eps_burn tri_alfa + + ! center d_eps_nuc_dlnd for reaction categories + + !c_d_eps_dlnd cno + !c_d_eps_dlnd tri_alfa + + ! center d_eps_nuc_dlnT for reaction categories + + !c_d_eps_dlnT cno + !c_d_eps_dlnT tri_alfa + + !## regions of strong nuclear burning + + ! 2 zones where eps_nuc > burn_min1 erg/g/s + ! for each zone have 4 numbers: start1, start2, end2, end1 + ! start1 is mass of inner edge where first goes > burn_min1 (or -20 if none such) + ! start2 is mass of inner edge where first zone reaches burn_min2 erg/g/sec (or -20 if none such) + ! end2 is mass of outer edge where first zone drops back below burn_min2 erg/g/s + ! end1 is mass of outer edge where first zone ends (i.e. eps_nuc < burn_min1) + ! similar for the second zone + + epsnuc_M_1 ! start1 for 1st zone + epsnuc_M_2 ! start2 + epsnuc_M_3 ! end2 + epsnuc_M_4 ! end1 + + epsnuc_M_5 ! start1 for 2nd zone + epsnuc_M_6 ! start2 + epsnuc_M_7 ! end2 + epsnuc_M_8 ! end1 + + + ! you might want to get a more complete list of burning regions by using the following + + !burning_regions + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives int(sign(val)*log10(max(1,abs(val)))) + ! where val = ergs/gm/sec nuclear energy minus all neutrino losses. + ! the second column for a region gives the q location of the top of the region + ! entries for extra columns after the last region in the star will have a value of -9999 + ! all regions are included starting from the center, so the bottom of one region + ! is the top of the previous one. + ! since we start at the center, the bottom of the 1st region is q=0 and top of last is q=1. + + ! the columns in the log file will have names like 'burn_type_1' and 'burn_qtop_1' + + !burn_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'burn_relr_type_1' and 'burn_relr_top_1' + + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + +!---------------------------------------------------------------------------------------------- + +!# information about core and envelope + + !## helium core + he_core_mass + he_core_radius + !he_core_lgT + !he_core_lgRho + !he_core_L + !he_core_v + !he_core_omega + !he_core_omega_div_omega_crit + !he_core_k + + !## CO core + co_core_mass + !CO_core + co_core_radius + !co_core_lgT + !co_core_lgRho + !co_core_L + !co_core_v + !co_core_omega + !co_core_omega_div_omega_crit + !co_core_k + + !## ONe core + one_core_mass + !one_core_radius + !one_core_lgT + !one_core_lgRho + !one_core_L + !one_core_v + !one_core_omega + !one_core_omega_div_omega_crit + !one_core_k + + !## iron core + fe_core_mass + !fe_core_radius + !fe_core_lgT + !fe_core_lgRho + !fe_core_L + !fe_core_v + !fe_core_omega + !fe_core_omega_div_omega_crit + !fe_core_k + + !## neuton rich core + neutron_rich_core_mass + !neutron_rich_core_radius + !neutron_rich_core_lgT + !neutron_rich_core_lgRho + !neutron_rich_core_L + !neutron_rich_core_v + !neutron_rich_core_omega + !neutron_rich_core_omega_div_omega_crit + !neutron_rich_core_k + + !## envelope + + !envelope_mass ! = star_mass - he_core_mass + !envelope_fraction_left ! = envelope_mass / (initial_mass - he_core_mass) + + !h_rich_layer_mass ! = star_mass - he_core_mass + !he_rich_layer_mass ! = he_core_mass - c_core_mass + !co_rich_layer_mass + +!---------------------------------------------------------------------------------------------- + +!# timescales + + !dynamic_timescale ! dynamic timescale (seconds) -- estimated by 2*pi*sqrt(r^3/(G*m)) + !kh_timescale ! kelvin-helmholtz timescale (years) + !mdot_timescale ! star_mass/abs(star_mdot) (years) + !kh_div_mdot_timescales ! kh_timescale/mdot_timescale + !nuc_timescale ! nuclear timescale (years) -- proportional to mass divided by luminosity + + !dt_cell_collapse ! min time for any cell to collapse at current velocities + !dt_div_dt_cell_collapse + + !dt_div_max_tau_conv ! dt/ maximum conv timescale + !dt_div_min_tau_conv ! dt/ minimum conv timescale + + + !min_dr_div_cs ! min over all cells of dr/csound (seconds) + !min_dr_div_cs_k ! location of min + !log_min_dr_div_cs ! log10 min dr_div_csound (seconds) + !min_dr_div_cs_yr ! min over all cells of dr/csound (years) + !log_min_dr_div_cs_yr ! log10 min dr_div_csound (years) + !dt_div_min_dr_div_cs + !log_dt_div_min_dr_div_cs + + !min_t_eddy ! minimum value of scale_height/conv_velocity + +!---------------------------------------------------------------------------------------------- + +!# conditions at or near the surface of the model + + !## conditions at the photosphere + !effective_T + !Teff + log_Teff ! log10 effective temperature + ! Teff is calculated using Stefan-Boltzmann relation L = 4 pi R^2 sigma Teff^4, + ! where L and R are evaluated at the photosphere (tau_factor < 1) + ! or surface of the model (tau_factor >= 1) when photosphere is not inside the model. + + !photosphere_black_body_T + !photosphere_cell_T ! temperature at model location closest to the photosphere, not necessarily Teff + !photosphere_cell_log_T + !photosphere_cell_density + !photosphere_cell_log_density + !photosphere_cell_opacity + !photosphere_cell_log_opacity + !photosphere_L ! Lsun units + !photosphere_log_L ! Lsun units + !photosphere_r ! Rsun units + !photosphere_log_r ! Rsun units + !photosphere_m ! Msun units + !photosphere_v_km_s + !photosphere_cell_k + !photosphere_column_density + !photosphere_csound + !photosphere_log_column_density + !photosphere_opacity + !photosphere_v_div_cs + !photosphere_xm + !photosphere_cell_free_e + !photosphere_cell_log_free_e + !photosphere_logg + !photosphere_T + + !## conditions at or near the surface of the model (outer edge of outer cell) + + !luminosity ! luminosity in Lsun units + !luminosity_ergs_s ! luminosity in cgs units + log_L ! log10 luminosity in Lsun units + !log_L_ergs_s ! log10 luminosity in cgs units + !radius ! Rsun + log_R ! log10 radius in Rsun units + !radius_cm + !log_R_cm + + log_g ! log10 gravity + gravity + !log_Ledd + !log_L_div_Ledd ! log10(L/Leddington) + !lum_div_Ledd + !log_surf_optical_depth + !surface_optical_depth + + !log_surf_cell_opacity ! old name was log_surf_opacity + !log_surf_cell_P ! old name was log_surf_P + !log_surf_cell_pressure ! old name was log_surf_pressure + !log_surf_cell_density ! old name was log_surf_density + !log_surf_cell_temperature ! old name was log_surf_temperature + !surface_cell_temperature ! old name was surface_temperature + !log_surf_cell_z ! old name was log_surf_z + !surface_cell_entropy ! in units of kerg per baryon + ! old name was surface_entropy + + v_surf ! (cm/s) + !v_surf_km_s ! (km/s) + v_div_csound_surf ! velocity divided by sound speed at outermost grid point + !v_div_csound_max ! max value of velocity divided by sound speed at face + !v_div_vesc + !v_phot_km_s + !v_surf_div_escape_v + + !v_surf_div_v_kh ! v_surf/(photosphere_r/kh_timescale) + + !surf_avg_j_rot + !surf_avg_omega + !surf_avg_omega_crit + !surf_avg_omega_div_omega_crit + !surf_avg_v_rot ! km/sec rotational velocity at equator + !surf_avg_v_crit ! critical rotational velocity at equator + !surf_avg_v_div_v_crit + !surf_avg_Lrad_div_Ledd + !surf_avg_logT + !surf_avg_logRho + !surf_avg_opacity + + ! Gravity Darkening, reports the surface averaged L/Lsun and Teff (K) caused by + ! gravity darkening in rotating stars. Based on the model of Espinosa Lara & Rieutord (2011) + ! 'polar' refers to the line of sight being directed along the rotation axis of the star + ! 'equatorial' refers to the line of sight coincident with the stellar equator + !grav_dark_L_polar !Lsun + !grav_dark_Teff_polar !K + !grav_dark_L_equatorial !Lsun + !grav_dark_Teff_equatorial !K + + !surf_escape_v ! cm/s + + !v_wind_Km_per_s ! Km/s + ! = 1d-5*s% opacity(1)*max(0d0,-s% mstar_dot)/ & + ! (4*pi*s% photosphere_r*Rsun*s% tau_base) + ! Lars says: + ! wind_mdot = 4*pi*R^2*rho*v_wind + ! tau = integral(opacity*rho*dr) from R to infinity + ! so tau = opacity*wind_mdot/(4*pi*R*v_wind) at photosphere + ! or v_wind = opacity*wind_mdot/(4*pi*R*tau) at photosphere + + !rotational_mdot_boost ! factor for increase in mass loss mdot due to rotation + !log_rotational_mdot_boost ! log factor for increase in mass loss mdot due to rotation + !surf_r_equatorial_div_r_polar + !surf_r_equatorial_div_r + !surf_r_polar_div_r + +!---------------------------------------------------------------------------------------------- + +!# conditions near center + + log_center_T ! temperature + log_center_Rho ! density + log_center_P ! pressure + + ! shorter names for above + log_cntr_P + log_cntr_Rho + log_cntr_T + + !center_T ! temperature + !center_Rho ! density + !center_P ! pressure + + !center_degeneracy ! the electron chemical potential in units of k*T + !center_gamma ! plasma interaction parameter + center_mu + center_ye + center_abar + !center_zbar + + !center_eps_grav + + !center_non_nuc_neu + !center_eps_nuc + !d_center_eps_nuc_dlnT + !d_center_eps_nuc_dlnd + !log_center_eps_nuc + + !center_entropy ! in units of kerg per baryon + !max_entropy ! in units of kerg per baryon + !fe_core_infall + !non_fe_core_infall + !non_fe_core_rebound + !max_infall_speed + + !compactness_parameter ! (m/Msun)/(R(m)/1000km) for m = 2.5 Msun + !compactness + !m4 ! Mass co-ordinate where entropy=4 + ! mu4 is sensitive to the choice of how much dm/dr you average over, thus we average dm and dr over M(entropy=4) and M(entropy=4)+0.3Msun + !mu4 ! dM(Msun)/dr(1000km) where entropy=4 + + + !center_omega + !center_omega_div_omega_crit + +!---------------------------------------------------------------------------------------------- + +!# abundances + + !species ! size of net + + !## mass fractions near center + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_center_abundances + !add_log_center_abundances + + ! individual central mass fractions (as many as desired) + center h1 + center he4 + center c12 + center o16 + + ! individual log10 central mass fractions (as many as desired) + !log_center h1 + !log_center he4 + ! etc. + + + !## mass fractions near surface + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_surface_abundances + !add_log_surface_abundances + + ! individual surface mass fractions (as many as desired) + !surface h1 + !surface he4 + surface c12 + surface o16 + ! etc. + + ! individual log10 surface mass fractions (as many as desired) + + !log_surface h1 + !log_surface he4 + + + !## mass fractions for entire star + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_average_abundances + !add_log_average_abundances + + ! individual average mass fractions (as many as desired) + !average h1 + !average he4 + ! etc. + + ! individual log10 average mass fractions (as many as desired) + !log_average h1 + !log_average he4 + ! etc. + + + !## mass totals for entire star (in Msun units) + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_total_mass + !add_log_total_mass + + ! individual mass totals for entire star (as many as desired) + total_mass h1 + total_mass he4 + ! etc. + + ! individial log10 mass totals for entire star (in Msun units) + !log_total_mass h1 + !log_total_mass he4 + ! etc. + +!---------------------------------------------------------------------------------------------- + +!# info at specific locations + + !## info at location of max temperature + !max_T + !log_max_T + + +!---------------------------------------------------------------------------------------------- + +!# information about shocks + + !## info about outermost outward moving shock + ! excluding locations with q > max_q_for_outer_mach1_location + ! returns values at location of max velocity + !shock_mass ! baryonic (Msun) + !shock_mass_gm ! baryonic (grams) + !shock_q + !shock_radius ! (Rsun) + !shock_radius_cm ! (cm) + !shock_velocity + !shock_csound + !shock_v_div_cs + !shock_lgT + !shock_lgRho + !shock_lgP + !shock_gamma1 + !shock_entropy + !shock_tau + !shock_k + !shock_pre_lgRho + +!---------------------------------------------------------------------------------------------- + +!# asteroseismology + + !delta_nu ! large frequency separation for p-modes (microHz) + ! 1e6/(seconds for sound to cross diameter of star) + !delta_Pg ! g-mode period spacing for l=1 (seconds) + ! sqrt(2) pi^2/(integral of brunt_N/r dr) + !log_delta_Pg + !nu_max ! estimate from scaling relation (microHz) + ! nu_max = nu_max_sun * M/Msun / ((R/Rsun)^2 (Teff/Teff_sun)^0.5) + ! with nu_max_sun = 3100 microHz, Teff_sun = 5777 + !nu_max_3_4th_div_delta_nu ! nu_max^0.75/delta_nu + !acoustic_cutoff ! 0.5*g*sqrt(gamma1*rho/P) at surface + !acoustic_radius ! integral of dr/csound (seconds) + !ng_for_nu_max ! = 1 / (nu_max*delta_Pg) + ! period for g-mode with frequency nu_max = nu_max_ng*delta_Pg + !gs_per_delta_nu ! delta_nu / (nu_max**2*delta_Pg) + ! number of g-modes per delta_nu at nu_max + + !int_k_r_dr_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=1 + !int_k_r_dr_2pt0_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=1 + !int_k_r_dr_0pt5_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=1 + !int_k_r_dr_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=2 + !int_k_r_dr_2pt0_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=2 + !int_k_r_dr_0pt5_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=2 + !int_k_r_dr_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=3 + !int_k_r_dr_2pt0_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=3 + !int_k_r_dr_0pt5_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=3 + +!---------------------------------------------------------------------------------------------- + +!# energy information + + !total_energy ! at end of step + !log_total_energy ! log(abs(total_energy)) + !total_energy_after_adjust_mass ! after mass adjustments + + ! shorter versions of above + !tot_E + !log_tot_E + + + !total_gravitational_energy + !log_total_gravitational_energy ! log(abs(total_gravitational_energy)) + !total_gravitational_energy_after_adjust_mass + + ! shorter versions of above + !tot_PE + !log_tot_PE + + !total_internal_energy + !log_total_internal_energy + !total_internal_energy_after_adjust_mass + + ! shorter versions of above + !tot_IE + !log_tot_IE + + !total_radial_kinetic_energy + !log_total_radial_kinetic_energy + !total_radial_kinetic_energy_after_adjust_mass + + ! shorter versions of above (does not include rot KE) + !tot_KE + !log_tot_KE + + !total_turbulent_energy + !log_total_turbulent_energy + !total_turbulent_energy_after_adjust_mass + !tot_Et + !log_tot_Et + + !total_energy_foe + + !tot_IE_div_IE_plus_KE + !total_IE_div_IE_plus_KE + + !total_entropy + !total_eps_grav + + !total_energy_sources_and_sinks ! for this step + !total_nuclear_heating + !total_non_nuc_neu_cooling + !total_irradiation_heating + !total_extra_heating ! extra heat integrated over the model times dt (erg) + !total_WD_sedimentation_heating + + !rel_run_E_err + + !rel_E_err + !abs_rel_E_err + !log_rel_E_err + + !tot_e_equ_err + !tot_e_err + + + !error_in_energy_conservation ! for this step + ! = total_energy - (total_energy_start + total_energy_sources_and_sinks) + !cumulative_energy_error ! = sum over all steps of abs(error_in_energy_conservation) + !rel_cumulative_energy_error ! = cumulative_energy_error/total_energy + !log_rel_cumulative_energy_error ! = log10 of rel_cumulative_energy_error + !log_rel_run_E_err ! shorter name for rel_cumulative_energy_error + + !rel_error_in_energy_conservation ! = error_in_energy_conservation/total_energy + !log_rel_error_in_energy_conservation + + !virial_thm_P_avg + !virial_thm_rel_err + !work_inward_at_center + !work_outward_at_surface + + +!---------------------------------------------------------------------------------------------- + + !# rotation + + !total_angular_momentum + log_total_angular_momentum + !i_rot_total ! moment of inertia + + !total_rotational_kinetic_energy + !log_total_rotational_kinetic_energy + !total_rotational_kinetic_energy_after_adjust_mass + +!---------------------------------------------------------------------------------------------- + +!# velocities + + !avg_abs_v_div_cs + !log_avg_abs_v_div_cs + !max_abs_v_div_cs + !log_max_abs_v_div_cs + + !avg_abs_v + !log_avg_abs_v + !max_abs_v + !log_max_abs_v + + !u_surf + !u_surf_km_s + !u_div_csound_surf + !u_div_csound_max + + !infall_div_cs + +!---------------------------------------------------------------------------------------------- + +!# misc + + !e_thermal ! sum over all zones of Cp*T*dm + + !## eos + !logQ_max ! logQ = logRho - 2*logT + 12 + !logQ_min + !gamma1_min + + !## core mixing + !mass_semiconv_core + + !## H-He boundary + + !diffusion_time_H_He_bdy + !temperature_H_He_bdy + + + !## optical depth and opacity + + !one_div_yphot + !log_one_div_yphot + + !log_min_opacity + !min_opacity + + !log_tau_center + + !log_max_tau_conv + !max_tau_conv + !log_min_tau_conv + !min_tau_conv + + !tau_qhse_yrs + + !## other + + !Lsurf_m + !dlnR_dlnM + !h1_czb_mass ! location (in Msun units) of base of 1st convection zone above he core + !kh_mdot_limit + !log_cntr_dr_cm + !min_Pgas_div_P + !surf_c12_minus_o16 ! this is useful for seeing effects of dredge up on AGB + !surf_num_c12_div_num_o16 + + !phase_of_evolution ! Integer mapping to the type of evolution see star_data/public/star_data_def.inc for definitions + + !## MLT++ + !gradT_excess_alpha + !gradT_excess_min_beta + !gradT_excess_max_lambda + + !max_L_rad_div_Ledd + !max_L_rad_div_Ledd_div_phi_Joss + + + !## RTI + !rti_regions + + !## Ni & Co + !total_ni_co_56 + + + !## internal structure constants + + ! this is evaluated assuming a spherical star and does not account for rotation + !apsidal_constant_k2 + + +!---------------------------------------------------------------------------------------------- + +!# accretion + + !k_below_const_q + !q_below_const_q + !logxq_below_const_q + + !k_const_mass + !q_const_mass + !logxq_const_mass + + !k_below_just_added + !q_below_just_added + !logxq_below_just_added + + !k_for_test_CpT_absMdot_div_L + !q_for_test_CpT_absMdot_div_L + !logxq_for_test_CpT_absMdot_div_L + +!---------------------------------------------------------------------------------------------- + +!# Color output + + ! Outputs the bolometric correction (bc) for the star in filter band ``filter'' (case sensitive) + !bc filter + + ! Outputs the absolute magnitude for the star in filter band ``filter'' (case sensitive) + !abs_mag filter + + ! Adds all the bc's to the output + !add_bc + + ! Adds all the absolute magnitudes to the output + !add_abs_mag + + ! Outputs luminosity in filter band ``filter'' (erg s^-1) (case sensitive) + ! lum_band filter + + ! Adds all the filter band luminosities to the output (erg s^-1) + ! add_lum_band + + ! Outputs log luminosity in filter band ``filter'' (log erg s^-1) (case sensitive) + ! log_lum_band filter + + ! Adds all the filter band luminosities to the output (log erg s^-1) + ! add_log_lum_band + +!---------------------------------------------------------------------------------------------- + +!# RSP + + !rsp_DeltaMag ! absolute magnitude difference between minimum and maximum light (mag) + !rsp_DeltaR ! R_max - R_min difference in the max and min radius (Rsun) + !rsp_GREKM ! fractional growth of the kinetic energy per pulsation period ("nonlinear growth rate") - see equation 5 in MESA5 + !rsp_num_periods ! Count of the number of pulsation cycles completed + !rsp_period_in_days ! Running period, ie., period between two consecutive values of R_max (days) + !rsp_phase ! Running pulsation phase for a cycle + +!---------------------------------------------------------------------------------------------- +!# debugging + + !## retries + num_retries ! total during the run + + !## solver iterations + + num_iters ! same as num_solver_iterations + !num_solver_iterations ! iterations at this step + !total_num_solver_iterations ! total iterations during the run + !avg_num_solver_iters + + !rotation_solver_steps + + !diffusion_solver_steps + !diffusion_solver_iters + + !avg_setvars_per_step + !avg_skipped_setvars_per_step + !avg_solver_setvars_per_step + + !burn_solver_maxsteps + + !total_num_solver_calls_converged + !total_num_solver_calls_failed + !total_num_solver_calls_made + !total_num_solver_relax_calls_converged + !total_num_solver_relax_calls_failed + !total_num_solver_relax_calls_made + !total_num_solver_relax_iterations + + !total_step_attempts + !total_step_redos + !total_step_retries + !total_steps_finished + !total_steps_taken + + !TDC_num_cells + + !## Relaxation steps + !total_relax_step_attempts + !total_relax_step_redos + !total_relax_step_retries + !total_relax_steps_finished + !total_relax_steps_taken + + !## conservation during mesh adjust + !log_mesh_adjust_IE_conservation + !log_mesh_adjust_KE_conservation + !log_mesh_adjust_PE_conservation + + !## amr + !num_hydro_merges + !num_hydro_splits + + !## timing + !elapsed_time ! time since start of run (seconds) + + !## Extras + burning_regions 40 + mixing_regions 40 + mix_relr_regions 40 diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_after_first_pulse b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_after_first_pulse new file mode 100644 index 000000000..c3bee1f55 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_after_first_pulse @@ -0,0 +1,17 @@ + +&star_job + +/ !end of star_job namelist + + +&controls + !inlist is currently unused, things can be added here + !to be loaded after the first pulse begins + max_years_for_timestep = 0.5d0 + !max_q_for_convection_with_hydro_on = 0.9999d0 + +/ ! end of controls namelist + +&pgstar + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_extra b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_extra new file mode 100644 index 000000000..23b059254 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_extra @@ -0,0 +1,14 @@ +&kap +Zbase = 0.00142d0 +/ ! end of kap namelist + +&star_job +new_omega_div_omega_crit = 0.10d0 +/ ! end of controls namelist + +&controls +initial_mass = 72.50d0 +initial_Y = 0.99858d0 +initial_Z = 0.00142d0 +initial_he3 = 0d0 +/ ! end of controls namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_off b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_off new file mode 100644 index 000000000..d887bf748 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_off @@ -0,0 +1,26 @@ + +&star_job + +/ !end of star_job namelist + + +&controls + use_other_wind = .true. + + use_dPrad_dm_form_of_T_gradient_eqn = .true. + use_momentum_outer_BC = .false. + use_split_merge_amr = .false. + + max_surface_cell_dq = 1d-12 + + max_mdot_redo_cnt = 200 + + restore_mesh_on_retry = .false. + + MLT_option = 'TDC' + +/ ! end of controls namelist + +&pgstar + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_on b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_on new file mode 100644 index 000000000..ac232a095 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_hydro_on @@ -0,0 +1,36 @@ + +&star_job + +/ !end of star_job namelist + + +&controls + use_other_wind = .false. + + !use_ODE_var_eqn_pairing = .true. + use_dPrad_dm_form_of_T_gradient_eqn = .true. + use_momentum_outer_BC = .true. + + use_split_merge_amr = .true. + + max_surface_cell_dq = 1d0 + split_merge_amr_dq_min = 1d-16 + + max_mdot_redo_cnt = 0 + + restore_mesh_on_retry = .true. + num_steps_to_hold_mesh_after_retry = 5 + + MLT_option = 'TDC' + + + !v_drag is set in run_star_extras + q_for_v_drag_full_off = 0.98d0 + q_for_v_drag_full_on = 0.99d0 + max_surface_cell_dq = 1d0 + +/ ! end of controls namelist + +&pgstar + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pgstar b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pgstar new file mode 100644 index 000000000..e2bdeb174 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pgstar @@ -0,0 +1,210 @@ + +&pgstar + pgstar_interval = 1 + !pause = .true. + + pgstar_age_disp = 2.5 + pgstar_model_disp = 2.5 + + !### scale for axis labels + pgstar_xaxis_label_scale = 1.3 + pgstar_left_yaxis_label_scale = 1.3 + pgstar_right_yaxis_label_scale = 1.3 + + Grid2_win_flag = .true. + !Grid2_win_flag = .false. + + Grid2_win_width = 13 !15 + Grid2_win_aspect_ratio = 0.65 ! aspect_ratio = height/width + + Grid2_plot_name(4) = 'Mixing' + + Grid2_num_cols = 7 ! divide plotting region into this many equal width cols + Grid2_num_rows = 8 ! divide plotting region into this many equal height rows + + Grid2_num_plots = 6 ! <= 10 + + Grid2_plot_name(1) = 'TRho_Profile' + Grid2_plot_row(1) = 1 ! number from 1 at top + Grid2_plot_rowspan(1) = 3 ! plot spans this number of rows + Grid2_plot_col(1) = 1 ! number from 1 at left + Grid2_plot_colspan(1) = 2 ! plot spans this number of columns + Grid2_plot_pad_left(1) = -0.05 ! fraction of full window width for padding on left + Grid2_plot_pad_right(1) = 0.01 ! fraction of full window width for padding on right + Grid2_plot_pad_top(1) = 0.00 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(1) = 0.05 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(1) = 0.65 ! multiply txt_scale for subplot by this + + + Grid2_plot_name(5) = 'Kipp' + Grid2_plot_row(5) = 4 ! number from 1 at top + Grid2_plot_rowspan(5) = 3 ! plot spans this number of rows + Grid2_plot_col(5) = 1 ! number from 1 at left + Grid2_plot_colspan(5) = 2 ! plot spans this number of columns + Grid2_plot_pad_left(5) = -0.05 ! fraction of full window width for padding on left + Grid2_plot_pad_right(5) = 0.01 ! fraction of full window width for padding on right + Grid2_plot_pad_top(5) = 0.03 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(5) = 0.0 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(5) = 0.65 ! multiply txt_scale for subplot by this + Kipp_title = '' + Kipp_show_mass_boundaries = .false. + + Grid2_plot_name(6) = 'History_Panels1' + Grid2_plot_row(6) = 6 ! number from 1 at top + Grid2_plot_rowspan(6) = 3 ! plot spans this number of rows + Grid2_plot_col(6) = 6 ! number from 1 at left + Grid2_plot_colspan(6) = 2 ! plot spans this number of columns + !Grid2_plot_pad_left(6) = 0.00 ! fraction of full window width for padding on left + !Grid2_plot_pad_right(6) = 0.05 ! fraction of full window width for padding on right + !Grid2_plot_pad_top(6) = 0.03 ! fraction of full window height for padding at top + !Grid2_plot_pad_bot(6) = 0.0 ! fraction of full window height for padding at bottom + !Grid2_txt_scale_factor(6) = 0.65 ! multiply txt_scale for subplot by this + + Grid2_plot_pad_left(6) = 0.05 ! fraction of full window width for padding on left + Grid2_plot_pad_right(6) = 0.03 ! fraction of full window width for padding on right + Grid2_plot_pad_top(6) = 0.0 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(6) = 0.0 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(6) = 0.65 ! multiply txt_scale for subplot by this + + History_Panels1_title = '' + History_Panels1_num_panels = 2 + + History_Panels1_xaxis_name='model_number' + History_Panels1_max_width = -1 ! only used if > 0. causes xmin to move with xmax. + + History_Panels1_yaxis_name(1) = 'log_center_T' + History_Panels1_yaxis_reversed(1) = .false. + History_Panels1_ymin(1) = -101d0 ! only used if /= -101d0 + History_Panels1_ymax(1) = -101d0 ! only used if /= -101d0 + History_Panels1_dymin(1) = 1 + + History_Panels1_other_yaxis_name(1) = 'gamma_integral' + History_Panels1_other_yaxis_reversed(1) = .false. + History_Panels1_other_ymin(1) = -0.06d0 ! only used if /= -101d0 + History_Panels1_other_ymax(1) = 0.1d0 ! only used if /= -101d0 + History_Panels1_other_dymin(1) = 0.05d0 + + History_Panels1_yaxis_name(2) = 'log_Lnuc' + History_Panels1_yaxis_reversed(2) = .false. + History_Panels1_ymin(2) = 0d0 ! only used if /= -101d0 + History_Panels1_ymax(2) = 20d0 ! only used if /= -101d0 + History_Panels1_dymin(2) = 1 + + History_Panels1_other_yaxis_name(2) = 'log_Lneu' + History_Panels1_other_yaxis_reversed(2) = .false. + History_Panels1_other_ymin(2) = 0d0 ! only used if /= -101d0 + History_Panels1_other_ymax(2) = 20d0 ! only used if /= -101d0 + History_Panels1_other_dymin(2) = 1 + + Grid2_plot_name(2) = 'Text_Summary1' + Grid2_plot_row(2) = 7 ! number from 1 at top + Grid2_plot_rowspan(2) = 2 ! plot spans this number of rows + Grid2_plot_col(2) = 1 ! number from 1 at left + Grid2_plot_colspan(2) = 4 ! plot spans this number of columns + Grid2_plot_pad_left(2) = -0.08 ! fraction of full window width for padding on left + Grid2_plot_pad_right(2) = -0.10 ! fraction of full window width for padding on right + Grid2_plot_pad_top(2) = 0.08 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(2) = -0.04 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(2) = 0.19 ! multiply txt_scale for subplot by this + + Text_Summary1_name(2,1) = 'star_age' + Text_Summary1_name(3,1) = 'time_step_sec' + Text_Summary1_name(8,1) = 'yr_since_coll' + Text_Summary1_name(3,2) = 'co_core_mass' + Text_Summary1_name(4,2) = 'M_below_vesc' + Text_Summary1_name(5,2) = 'gamma_integral' + Text_Summary1_name(6,2) = 'log_R_vesc_098' + Text_Summary1_name(7,2) = 'log_R_vesc_095' + Text_Summary1_name(8,2) = 'log_R_vesc_090' + + Text_Summary1_name(8,4) = 'center_si28' + + Grid2_plot_name(3) = 'Profile_Panels3' + Profile_Panels3_title = 'Abundance-Power-Mixing' + Profile_Panels3_num_panels = 3 + Profile_Panels3_yaxis_name(1) = 'Abundance' + Profile_Panels3_yaxis_name(2) = 'Power' + Profile_Panels3_yaxis_name(3) = 'Mixing' + + Profile_Panels3_xaxis_name = 'mass' + Profile_Panels3_xaxis_reversed = .false. + Profile_Panels3_xmin = -101d0 ! only used if /= -101d0 + Profile_Panels3_xmax = -101d0 ! 10 ! only used if /= -101d0 + + Grid2_plot_row(3) = 1 ! number from 1 at top + Grid2_plot_rowspan(3) = 6 ! plot spans this number of rows + Grid2_plot_col(3) = 3 ! plot spans this number of columns + Grid2_plot_colspan(3) = 3 ! plot spans this number of columns + + Grid2_plot_pad_left(3) = 0.09 ! fraction of full window width for padding on left + Grid2_plot_pad_right(3) = 0.07 ! fraction of full window width for padding on right + Grid2_plot_pad_top(3) = 0.0 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(3) = 0.0 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(3) = 0.65 ! multiply txt_scale for subplot by this + + Grid2_plot_name(4) = 'Profile_Panels1' + Grid2_plot_row(4) = 1 ! number from 1 at top + Grid2_plot_rowspan(4) = 5 ! plot spans this number of rows + Grid2_plot_col(4) = 6 ! number from 1 at left + Grid2_plot_colspan(4) = 2 ! plot spans this number of columns + Grid2_plot_pad_left(4) = 0.05 ! fraction of full window width for padding on left + Grid2_plot_pad_right(4) = 0.03 ! fraction of full window width for padding on right + Grid2_plot_pad_top(4) = 0.0 ! fraction of full window height for padding at top + Grid2_plot_pad_bot(4) = 0.07 ! fraction of full window height for padding at bottom + Grid2_txt_scale_factor(4) = 0.65 ! multiply txt_scale for subplot by this + + + Abundance_line_txt_scale_factor = 1.1 ! relative to other text + Abundance_legend_txt_scale_factor = 1.1 ! relative to other text + Abundance_legend_max_cnt = 0 + Abundance_log_mass_frac_min = -3.5 ! only used if < 0 + + Profile_Panels1_title = '' + + Profile_Panels1_xaxis_name = 'zone' + Profile_Panels1_xaxis_reversed = .true. + Profile_Panels1_xmin = -101d0 ! only used if /= -101d0 + Profile_Panels1_xmax = -101d0 ! only used if /= -101d0 + Profile_Panels1_num_panels = 3 + + Profile_Panels1_yaxis_name(1) = 'vel_km_per_s' + Profile_Panels1_other_yaxis_name(1) = 'gamma1' + Profile_Panels1_dymin(1) = 50 + Profile_Panels1_other_ymax(1) = 1.41 + Profile_Panels1_other_ymin(1) = 1.28 + + Profile_Panels1_yaxis_name(2) = 'omega' + Profile_Panels1_other_yaxis_name(2) = 'fp_rot' ! 'vel_km_per_s' + Profile_Panels1_other_ymax(2) = -101d0 ! -101d0 + Profile_Panels1_other_ymin(2) = -101d0 ! -101d0 + !Profile_Panels1_other_ymargin(2) = 0.1 ! 10 + !Profile_Panels1_other_dymin(2) = 0.5 + + + Profile_Panels1_yaxis_name(3) = 'logRho' + Profile_Panels1_ymax(3) = -101d0 + Profile_Panels1_ymin(3) = -2 + Profile_Panels1_dymin(3) = 0.1 + + Profile_Panels1_other_yaxis_name(3) = 'log_rel_E_err'!'radius' !'spin_parameter' ! 'mass' + Profile_Panels1_other_ymax(3) = -101d0!2!5d0 + Profile_Panels1_other_ymin(3) = -101d0!-2d0 + !Profile_Panels1_other_dymin(3) = 0.15 + + show_TRho_Profile_kap_regions = .false. + show_TRho_Profile_gamma1_4_3rd = .true. + TRho_Profile_xmin = -30 + !TRho_Profile_xmax = 10 + TRho_Profile_ymin = 2.5 + !TRho_Profile_ymax = 10 + + Grid2_file_flag = .true. + Grid2_file_dir = 'png' + Grid2_file_prefix = 'grid_' + Grid2_file_interval = 10 ! 1 ! output when mod(model_number,Grid2_file_interval)==0 + Grid2_file_width = -1 ! negative means use same value as for window + Grid2_file_aspect_ratio = -1 ! negative means use same value as for window + + !show_TRho_Profile_eos_regions = .true. + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_ppisn b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_ppisn new file mode 100644 index 000000000..4823f3f15 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_ppisn @@ -0,0 +1,395 @@ + +&kap + + ! no H, so do blend higher up + kap_blend_logT_upper_bdy = 4.1d0 + kap_blend_logT_lower_bdy = 4.0d0 + + ! We use 10% solar, with the definition of Asplund+ 2009 + ! Zbase defined in inlist_extra + use_Type2_opacities = .true. + + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_extra' +/ + +&eos + + +/ ! end of eos namelist + +&star_job + + timescale_for_relax_composition = 1d5 + + !change_initial_net = .true. + !although not used, run_star_extras reads new_net_name during pulses + new_net_name = "approx21_plus_co56.net" + adv_net = "approx21_plus_co56.net" + + initial_zfracs = 6 + + ! this are options for the relaxation of the model in-between pulses + timescale_for_relax_entropy = 1d-15 + max_dt_for_relax_entropy = 5d-17 + num_timescales_for_relax_entropy = 30 + max_steps_to_relax_entropy = 10000 + relax_omega_max_yrs_dt = 1d-5 + + set_initial_cumulative_energy_error = .true. + new_cumulative_energy_error = 0d0 + + change_rotation_flag = .true. + new_rotation_flag = .true. + change_w_div_wc_flag = .true. + new_w_div_wc_flag = .true. + + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_extra' + + pgstar_flag = .true. + save_pgstar_files_when_terminate = .true. + + num_special_rate_factors = 2 + reaction_for_special_factor(1) = 'r_c12_ag_o16' + special_rate_factor(1) = 1 + filename_of_special_rate(1) = 'c12ag_deboer_sigma_0p0_2000_Tgrid.dat' + + reaction_for_special_factor(2) = 'r_he4_he4_he4_to_c12' + special_rate_factor(2) = 1 + filename_of_special_rate(2) = 'r_he4_he4_he4_to_c12_cf88.txt' + +/ !end of star_job namelist + + +&controls + +! atmosphere + ! Pextra_factor = 1 + ! ! need this extra pressure to stabilize the atmosphere during core He burning + + ! atm_option = 'fixed_Tsurf' !'T_tau' + ! !atm_T_tau_relation = 'Eddington' + ! !atm_T_tau_opacity = 'fixed' + ! atm_fixed_Tsurf = 160000 + + + + okay_to_reduce_gradT_excess = .false. + gradT_excess_f1 = 1d-4 + gradT_excess_f2 = 1d-2 + gradT_excess_lambda1 = -1d0 ! full on + + use_superad_reduction = .false. + superad_reduction_Gamma_limit = 0.3d0 + superad_reduction_Gamma_limit_scale = 5d0 + superad_reduction_Gamma_inv_scale = 5d0 + superad_reduction_diff_grads_limit = 1d-2 + superad_reduction_limit = -1d0 + + + fill_arrays_with_nans = .false. + report_solver_progress = .true. + report_ierr = .false. ! if true, produce terminal output when have some internal error + + max_resid_jump_limit = 1d99 + + okay_to_remove_mixing_singleton = .false. + + num_trace_history_values = 2 + trace_history_value_name(1) = 'log_rel_run_E_err' + trace_history_value_name(2) = 'rel_E_err' + warn_when_large_rel_run_E_err = 1d-2 + energy_eqn_option = 'dedt' + + + limit_for_rel_error_in_energy_conservation = 1d1 + hard_limit_for_rel_error_in_energy_conservation = 1d2 + + ! adjustments to the newton solver + solver_max_tries_before_reject = 30 + max_tries_for_retry = 30 + max_tries_after_5_retries = 40 + max_tries_after_10_retries = 40 + max_tries_after_20_retries = 40 + corr_coeff_limit = 1d-2 + + + !relax_use_gold_tolerances = .false. + use_gold_tolerances = .true. + use_gold2_tolerances = .false. + maxT_for_gold_tolerances = 3.7d9 ! this is to ease core-collapse (~logT=9.7) + gold_iter_for_resid_tol2 = 10 + gold_iter_for_resid_tol3 = 10 + gold_tol_residual_norm3 = 1d-6 + gold_tol_max_residual3 = 1d-3 + gold_solver_iters_timestep_limit = 20 + solver_iters_timestep_limit = 50 + ignore_too_large_correction = .true. + scale_max_correction = 0.1d0 + !corr_coeff_limit = 0.2d0 + ignore_min_corr_coeff_for_scale_max_correction = .true. + ignore_species_in_max_correction = .true. + + ! Can be helpful to decrease op_split_burn_min_T + ! to 3d9 or even 1d9 for large networks to avoid crashes + op_split_burn = .false. + op_split_burn_min_T = 3.0d9 + burn_steps_limit = 150 + burn_steps_hard_limit = 250 + op_split_burn_eps = 1d-5 + op_split_burn_odescal = 1d-5 + + + + mlt_make_surface_no_mixing = .true. + convergence_ignore_equL_residuals = .true. + make_gradr_sticky_in_solver_iters = .true. + xa_scale = 1d-5 + iter_for_resid_tol2 = 10 + + + + ! during pulses very small cells near the surface can sometimes exceed + ! the speed of light. This has no impact in the behaviour of the bulk + ! of the star so we don't want to have a retry if that happens + retry_for_v_above_clight = .false. + + ! limit max_model_number as part of test_suite + max_model_number = 99999999!10000 ! RECOMMENDED -1 + + relax_max_number_retries = 99999999!1000 + max_number_retries = 99999999 !500 ! RECOMMENDED 5000 + + ! in principle this is the only thing that needs changing + ! it is set in inlist_extra + !initial_mass = 72 + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_extra' + + ! our wind implementation follows Brott+ 2011 + use_other_wind = .true. + + ! when using hydro, we reduce the rotational corrections near the surface + use_other_eval_fp_ft = .true. + + ! During hydro we use es scaterring opacity on the very outermost cell + use_other_kap = .true. + + ! convection controls + use_ledoux_criterion = .true. + mixing_length_alpha = 2d0 + alpha_semiconvection = 1d0 + thermohaline_coeff = 0d0 + num_cells_for_smooth_gradL_composition_term = 0 + max_conv_vel_div_csound = 1d0 + + max_abs_du_div_cs_for_convection = 0.03d0 + max_v_div_cs_for_convection = 1d1 + max_v_for_convection = 1d4 + + ! overshoot controls + ! we only include a bit of exponential overshooting to smooth things out + + overshoot_scheme(1) = 'exponential' + overshoot_zone_type(1) = 'any' + overshoot_zone_loc(1) = 'any' + overshoot_bdy_loc(1) = 'any' + overshoot_f(1) = 0.01 + overshoot_f0(1) = 0.005 + + ! termination conditions + ! make this a bit higher to avoid interrupting pair-instability collapse + fe_core_infall_limit = 8d8 + + ! timestep options + varcontrol_target = 5d-4 + max_timestep_factor = 1.2d0 + min_timestep_factor = 0.8d0 + dX_nuc_drop_limit = 5d-2 + dX_nuc_drop_limit_at_high_T = 1d-2 ! for center logT > 9.45 + delta_Ye_highT_limit = 1d-3 + dX_nuc_drop_max_A_limit = 52 + dX_nuc_drop_min_X_limit = 1d-4 + dX_nuc_drop_hard_limit = 1d99 + delta_lgTeff_limit = 1d0 + delta_lgL_limit = -1!1d99 + delta_lgL_hard_limit = -1!1d200 + delta_lgR_limit = 0.025d0 + delta_lgR_hard_limit = 0.05d0 + delta_lgL_He_limit = -1d0 + lgL_nuc_burn_min = 4d0 + retry_hold = 0 + + ! limit for changes in central abundances, RECOMMENDED 0.001d0 for all + delta_XH_cntr_limit = 0.01d0 + delta_XHe_cntr_limit = 0.01d0 + delta_XC_cntr_limit = 0.01d0 + delta_XO_cntr_limit = 0.01d0 + + ! extra controls for timestep + delta_lg_star_mass_limit = 2d-3 ! RECOMMENDED 2d-3 + delta_lgRho_cntr_limit = 0.005d0 ! RECOMMENDED 0.0025d0 + delta_lgRho_cntr_hard_limit = 0.01d0 ! RECOMMENDED 0.005d0 + dt_div_min_dr_div_cs_limit = 5d0 ! this value is adjusted in run_star_extras ! RECOMMENDED 3d0 + min_timestep_limit = 1d-20 ! (seconds) + relax_hard_limits_after_retry = .false. + + ! mesh controls + okay_to_remesh = .true. + max_dq = 1d-3 ! RECOMMENDED 1d-3 + mesh_delta_coeff = 0.8d0 ! RECOMMENDED 0.8d0 + !! this one is turned on in run_star_extras + !use_split_merge_amr = .true. + split_merge_amr_log_zoning = .true. + split_merge_amr_nz_baseline = 6000 ! RECOMMENDED 6000 + split_merge_amr_MaxLong = 1.25d0 + split_merge_amr_MaxShort = 2.5d0 + split_merge_amr_max_iters = 1000 + split_merge_amr_okay_to_split_nz = .false. + merge_amr_ignore_surface_cells = .false. + merge_amr_max_abs_du_div_cs = 0.05d0 + merge_amr_du_div_cs_limit_only_for_compression = .true. + split_merge_amr_avoid_repeated_remesh = .true. + + + + ! rotational mixing coeffs + am_nu_ST_factor = 1.0 + D_visc_factor = 0.0 + am_nu_SH_factor = 0.0 + D_ST_factor = 0.0 + D_SH_factor = 0.0 + D_GSF_factor = 1.0 + D_ES_factor = 1.0 + D_SSI_factor = 1.0 + D_DSI_factor = 1.0 + am_D_mix_factor = 0.0333333d0 + am_gradmu_factor = 0.1d0 + num_cells_for_smooth_gradL_composition_term = 2 + + ! use implicit wind close to critical + surf_avg_tau_min = 0 + surf_avg_tau = 10 + !max_mdot_redo_cnt = 200 ! this is set in inlist_hydro_on and inlist_hydro_off + min_years_dt_for_redo_mdot = 0 + surf_omega_div_omega_crit_limit = 0.98d0 + surf_omega_div_omega_crit_tol = 0.02d0 + rotational_mdot_boost_fac = 1d10 + rotational_mdot_kh_fac = 1d10 + mdot_revise_factor = 1.1 + implicit_mdot_boost = 0.05 + ! this needs to be relaxed just to avoid a crash when hydro Riemann is turned on + angular_momentum_error_retry = 1d99 + angular_momentum_error_warn = 1d-10 + + ! Fixing the position of the Lagrangian region of the mesh helps + ! convergence near the Eddington limit + max_logT_for_k_below_const_q = 100 + max_q_for_k_below_const_q = 0.99 + min_q_for_k_below_const_q = 0.99 + max_logT_for_k_const_mass = 100 + max_q_for_k_const_mass = 0.98 + min_q_for_k_const_mass = 0.98 + + + history_interval = 1 + max_num_profile_models = 10000 ! 100 ! RECOMMENDED 10000 + + warn_rates_for_high_temp = .false. + + ! Here begins the inferno of additional parameters to control turning hydro on/off and + ! relaxing the star to a lower mass + ! 'Abandon hope all ye who enter here' + ! also given are the variable names to which each is assigned in run_star_extras + + ! in order for Riemann hydro to be turned on, the following two conditions need to be met + x_ctrl(1) = 0.01d0 ! the integrated gamma1 must be smaller than this (min_gamma_sub_43_for_hydro) + + ! parameters controlling swtich between hydrostatic and Riemann hydro + ! after the pulse starts, we estimate the dynamical timescale of the star + ! (excluding some of the outer layers), and model the dynamical evolution + ! for at least a fixed number of these timescales before turning off hydro + x_ctrl(2) = 50d0 ! velocity in km/s at any zone for which we consider the system to be undergoing a pulse + ! time will start counting from this point before switching off hydro + ! (max_v_for_pulse) + ! NOTE: we only consider the velocities for q < x_ctrl(7) + x_ctrl(3) = 0.9d0 ! dynamical timescale will be estimated for the inside of the + ! star, up to a mass fraction equal to this value of the CO core mass. + ! If there is no CO core yet (unlikely but just to be sure), + ! the He core mass is used. The timescale is calculated as 1/sqrt(G*), + ! where =M/(4 pi R^3/3) up to the specified q + ! (q_for_dyn_ts) + x_ctrl(4) = 50d0 ! number of dynamical timescales to evolve before turning off hydro + ! (num_dyn_ts_for_relax) + x_ctrl(5) = 0.99d0 ! after the set number of timescales have passed, before turning off hydro + ! two more conditions must be met. This is to allow more time for additional + ! shocks to develop on the outermost layers, which can add additional heat there. + ! These conditions must be met at all points inside a given mass fraction of + ! what would remain of the star, and that fraction is given by this option. + ! For example, for x_ctrl(7) = 0.95, if a star with 55 Msun would eject 5 Msun + ! leaving a 50 Msun star, then the additional conditions need to be satisfied + ! in the inner 47.5 Msun of the star. + ! (q_for_relax_check) + ! The additional conditions are as follows + x_ctrl(6) = 20d0 ! absolute velocities must be below this in km/s (max_v_for_relax) + x_ctrl(7) = 0.5d0 ! the ratio of the (absolute) velocity to the sound speed must be below this (max_machn_for_relax) + x_ctrl(8) = 11d0 ! do not relax if log10 Lneu/Lsun is bigger than this (max_Lneu_for_relax) + x_ctrl(9) = 10d0 ! do not relax if log10 Lnuc/Lsun is bigger than this (max_Lnuc_for_relax) + + x_integer_ctrl(1) = 1 !10 ! conditions for relax must be satisfied for this number of consecutive steps + ! before the outer layers are removed + ! (num_steps_before_relax) + + x_ctrl(10) = 1d4 !5d2 ! With Riemann hydro, limit the timestep to this (in seconds) before the onset of a pulse. + ! prevents surface cells from going wild + ! (max_dt_before_pulse) + ! RECOMMENDED 1d4 + + x_ctrl(11) = 9d0 ! turn off wind mass loss if log10(Lneu/Lsun) exceeds this (max_Lneu_for_mass_loss) + + x_ctrl(12) = 0.025d0 ! 0.05d0 ! log(L_nuc) timestep limit is set to this, hard limit to double this value. + ! This is turned off during relax and near core collapse + ! (delta_lgLnuc_limit) + ! RECOMMENDED 0.025d0 + + x_ctrl(13) = 10d0 ! Ignore the L_nuc limit of x_ctrl(18) if L_photo exceeds this value, otherwise step to step + ! variations can lead to tiny timesteps + ! When we reach this value we switch the timestep limit for lgLnuc to lgLphot instead + ! (max_Lphoto_for_lgLnuc_limit) + + x_ctrl(14) = 25d0 ! If L_photo is above this limit, then the timestep limit given by x_ctrl(18) is not used for L_photo either + + x_ctrl(15) = 9d0 ! After first pulse, turn on v_flag if riemann hydro is off and logT_center falls below this. + ! (logT_for_v_flag) + x_ctrl(16) = 9d0 ! After first pulse, turn on v_flag if riemann hydro is off and log10(Lneu/Lsun) exceeds this. + ! (logLneu_for_v_flag) + + x_logical_ctrl(1) = .false.! If true, terminate after evolving 100d after first relax + ! meant for use in the test_suite + ! (stop_100d_after_pulse) + ! RECOMMENDED .false. + + x_ctrl(18) = 1d4 ! After a pulse begins, limit the timestep to this (in seconds). Ignored after outer layers expand + ! beyond 1d4 Rsun when x_logical_ctrl(1) is true + ! (max_dt_during_pulse) + ! RECOMMENDED 1d4 + + + x_ctrl(19) = 2d4 ! if surface velocity exceeds this, switch to fixed_vsurf BC with the current vsurf + ! (vsurf_for_fixed_bc) + +! remove extended layers +x_logical_ctrl(2) = .true. + + photo_digits = 8 + photo_interval = 50 + profile_interval = 50 + history_interval = 10 + terminal_interval = 1!10 + +/ ! end of controls namelist + +&pgstar + +/ diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses new file mode 100644 index 000000000..7b24a30ed --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses @@ -0,0 +1,62 @@ +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_ppisn' +/ + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_ppisn' +/ + +&star_job + create_pre_main_sequence_model = .false. + + save_model_when_terminate = .true. + save_model_filename = 'final.mod' + required_termination_code_string = 'Successful test: evolved 100 days past first relax' + + load_saved_model = .true. + load_model_filename = 'he_dep.mod' + + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_ppisn' + + change_initial_v_flag = .true. + new_v_flag = .true. + + +/ !end of star_job namelist + + +&controls + + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_ppisn' + + x_logical_ctrl(2) = .true. + + ! See Farag et al 2022 + set_min_d_mix=.true. + min_d_mix=1d-2 + + MLT_option = 'TDC' + + + !v_drag is set in run_star_extras + q_for_v_drag_full_off = 0.98d0 + q_for_v_drag_full_on = 0.99d0 + max_surface_cell_dq = 1d0 + + + !drag_coefficient = 0.8d0 + !min_q_for_drag = 0.98d0 + +!Pextra_factor = 2 +/ ! end of controls namelist + +&pgstar + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses_header b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses_header new file mode 100644 index 000000000..67f9191a3 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_pulses_header @@ -0,0 +1,41 @@ + +&star_job + + + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_pulses' + +/ ! end of star_job namelist + + +&eos + + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_pulses' + +/ ! end of eos namelist + + +&kap + + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_pulses' + +/ ! end of kap namelist + + + +&controls + + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_pulses' + +/ ! end of controls namelist + + +&pgstar + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pulses' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep new file mode 100644 index 000000000..9f939f3c3 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep @@ -0,0 +1,47 @@ +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_ppisn' +/ + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_ppisn' +/ + +&star_job + create_pre_main_sequence_model = .true. + + save_model_when_terminate = .true. + save_model_filename = 'he_dep.mod' + required_termination_code_string = 'xa_central_lower_limit' + + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_ppisn' + +/ !end of star_job namelist + + +&controls + + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_ppisn' + + ! stop when the center mass fraction of h1 drops below this limit + xa_central_lower_limit_species(1) = 'he4' + xa_central_lower_limit(1) = 1d-3 + + ! relax this or you get tiny timesteps during pre He-ZAMS + limit_for_rel_error_in_energy_conservation = -1d0 + hard_limit_for_rel_error_in_energy_conservation = -1d0 + max_abs_rel_run_E_err = 1d99 + + x_logical_ctrl(2) = .false. + +/ ! end of controls namelist + +&pgstar + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep_header b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep_header new file mode 100644 index 000000000..fbbfc3430 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/inlist_to_he_dep_header @@ -0,0 +1,41 @@ + +&star_job + + + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_to_he_dep' + +/ ! end of star_job namelist + + +&eos + + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_to_he_dep' + +/ ! end of eos namelist + + +&kap + + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_to_he_dep' + +/ ! end of kap namelist + + + +&controls + + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_to_he_dep' + +/ ! end of controls namelist + + +&pgstar + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_to_he_dep' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/make/makefile b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/make/makefile new file mode 100644 index 000000000..f2017faed --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/make/makefile @@ -0,0 +1,5 @@ + + +STAR = star + +include $(MESA_DIR)/star/work_standard_makefile diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/mk b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/mk new file mode 100755 index 000000000..aec7a5195 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/mk @@ -0,0 +1,13 @@ +#!/bin/bash + +function check_okay { + if [ $? -ne 0 ] + then + echo + echo "FAILED" + echo + exit 1 + fi +} + +cd make; make; check_okay diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/profile_columns.list b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/profile_columns.list new file mode 100644 index 000000000..4846a1b3b --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/profile_columns.list @@ -0,0 +1,962 @@ +! profile_columns.list -- determines the contents of star model profiles +! you can use a non-standard version by setting profile_columns_file in your inlist + +! units are cgs unless otherwise noted. + +! reorder the following names as desired to reorder columns. +! comment out the name to omit a column (fewer columns => less IO => faster running). +! remove '!' to restore a column. + +! if you have a situation where you want a non-standard set of columns, +! make a copy of this file, edit as desired, and give the new filename in your inlist +! as profile_columns_file. if you are just adding columns, you can 'include' this file, +! and just list the additions in your file. note: to include the standard default +! version, use include '' -- the 0 length string means include the default file. + +! if you need to have something added to the list of options, let me know.... + +! the first few lines of the profile contain general info about the model. +! for completeness, those items are described at the end of this file. + + +! note: you can include another list by doing +! include 'filename' +! include '' means include the default standard list file + + +! the following lines of the profile contain info for 1 zone per row, surface to center. + +! minimal set of enabled columns: + + zone ! numbers start with 1 at the surface + mass ! m/Msun. mass coordinate of outer boundary of cell. + logR ! log10(radius/Rsun) at outer boundary of zone + logT ! log10(temperature) at center of zone + logRho ! log10(density) at center of zone + logP ! log10(pressure) at center of zone + x_mass_fraction_H + y_mass_fraction_He + z_mass_fraction_metals + + +! everything below this line is deactivated + + +!# Structure + !logM ! log10(m/Msun) + !log_mass + !dm ! cell mass (grams) + !dm_bar ! boundary mass (grams) average of adjacent dm's + !logdq ! log10(dq) + !log_dq + !dq_ratio ! dq(k-1)/dq(k) + q ! fraction of star mass interior to outer boundary of this zone + !log_q ! log10(q) + !xq + + !grav ! gravitational acceleration (cm sec^2) + !log_g ! log10 gravitational acceleration (cm sec^2) + !g_div_r ! grav/radius (sec^2) + !r_div_g ! radius/grav (sec^-2) + !cgrav_factor ! = cgrav(k)/standard_cgrav + !vel_km_per_s ! velocity at outer boundary of zone (km/s) -- 0 if no velocity variable + + radius ! radius at outer boundary of zone (in Rsun units) + !radius_cm ! radius at outer boundary of zone (in centimeters) + !radius_km ! radius at outer boundary of zone (in kilometers) + !logR_cm ! log10 radius at outer boundary of zone (in centimeters) + !rmid ! radius at center by mass of zone (in Rsun units) + !r_div_R ! fraction of total radius + + velocity ! velocity at outer boundary of zone (cm/s) -- 0 if no velocity variable + !v_div_r ! velocity divided by radius + !v_times_t_div_r + !rho_times_r3 ! at face + !log_rho_times_r3 ! at face + !scale_height ! in Rsun units + pressure_scale_height ! in Rsun units + + !m_div_r ! gm/cm + !dmbar_m_div_r + !log_dmbar_m_div_r + !mass_grams ! mass coordinate of outer boundary of cell in grams + !mmid ! mass at midpoint of cell (average of mass coords of the cell boundaries) Msun units. + + !m_grav ! total enclosed gravitational mass. Msun units. + !m_grav_div_m_baryonic ! mass_gravitational/mass at cell boundary + !mass_correction_factor ! dm_gravitational/dm (dm is baryonic mass of cell) + + !xm ! mass exterior to point (Msun units) + dq ! mass of zone as a fraction of total star mass + !logxq ! log10(1-q) + !logxm ! log10(xm) + + !xr ! radial distance from point to surface (Rsun) + !xr_cm ! radial distance from point to surface (cm) + !xr_div_R ! radial distance from point to surface in units of star radius + !log_xr ! log10 radial distance from point to surface (Rsun) + !log_xr_cm ! log10 radial distance from point to surface (cm) + !log_xr_div_R ! log10 radial distance from point to surface in units of star radius + + !dr ! r(outer edge) - r(inner edge); radial extent of cell in cm. + !log_dr ! log10 cell width (cm) + !dv ! v(inner edge) - v(outer edge); rate at which delta_r is shrinking (cm/sec). + + !dt_dv_div_dr ! dt*dv/dr; need to have this << 1 for every cell + !dr_div_R ! cell width divided by star R + !log_dr_div_R ! log10 cell width divided by star R + !dr_div_rmid ! cell width divided by rmid + !log_dr_div_rmid ! log(dr_div_rmid) + + !dr_div_cs ! cell sound crossing time (sec) + !log_dr_div_cs ! log10 cell sound crossing time (sec) + !dr_div_cs_yr ! cell sound crossing time (years) + !log_dr_div_cs_yr ! log10 cell sound crossing time (years) + + !acoustic_radius ! sound time from center to outer cell boundary (sec) + !log_acoustic_radius ! log10(acoustic_radius) (sec) + !acoustic_depth ! sound time from surface to outer cell boundary (sec) + !log_acoustic_depth ! log10(acoustic_depth) (sec) + !acoustic_r_div_R_phot + + !cell_collapse_time ! only set if doing explicit hydro + ! time (seconds) for cell inner edge to catch cell outer edge at current velocities + ! 0 if distance between inner and outer is increasing + !log_cell_collapse_time ! log of cell_collapse_time + + !compression_gradient + + + +!# Thermodynamics + !temperature ! temperature at center of zone + !logT_face ! log10(temperature) at outer boundary of zone + !logT_bb ! log10(black body temperature) at outer boundary of zone + !logT_face_div_logT_bb + + energy ! internal energy (ergs/g) + !logE ! log10(specific internal energy) at center of zone + !rho ! density + !density ! rho + + entropy ! specific entropy divided by (avo*kerg) + !logS ! log10(specific entropy) + !logS_per_baryon ! log10(specific entropy per baryon / kerg) + + !pressure ! total pressure at center of zone (pgas + prad) + !prad ! radiation pressure at center of zone + !pgas ! gas pressure at center of zone (electrons and ions) + logPgas ! log10(pgas) + pgas_div_ptotal ! pgas/pressure + + eta ! electron degeneracy parameter (eta >> 1 for significant degeneracy) + mu ! mean molecular weight per gas particle (ions + free electrons) + + grada ! dlnT_dlnP at constant S + !dE_dRho ! at constant T + !cv ! specific heat at constant volume + !cp ! specific heat at constant total pressure + + !log_CpT + gamma1 ! dlnP_dlnRho at constant S + !gamma3 ! gamma3 - 1 = dlnT_dlnRho at constant S + !gam ! plasma interaction parameter (> 160 or so means starting crystallization) + free_e ! free_e is mean number of free electrons per nucleon + !logfree_e ! log10(free_e), free_e is mean number of free electrons per nucleon + !chiRho ! dlnP_dlnRho at constant T + !chiT ! dlnP_dlnT at constant Rho + + csound ! sound speed + !log_csound + !csound_face ! sound speed (was previously called csound_at_face) + !cs_at_cell_bdy ! sound speed at cell boundary (csound is at cell center) + !v_div_cs ! velocity divided by sound speed + v_div_csound ! velocity divided by sound speed + !div_v + + !thermal_time_to_surface ! in seconds + !log_thermal_time_to_surface + !t_rad + !log_t_rad + !log_t_sound + !log_t_thermal + + !eos_phase + !eos_frac_OPAL_SCVH + !eos_frac_HELM + !eos_frac_Skye + !eos_frac_PC + !eos_frac_FreeEOS + !eos_frac_CMS + !eos_frac_ideal + + !pgas_div_p + !prad_div_pgas + !prad_div_pgas_div_L_div_Ledd + !pressure_scale_height_cm + + !eps_grav_composition_term + eps_grav_plus_eps_mdot + + !chiRho_for_partials + !chiT_for_partials + !rel_diff_chiRho_for_partials + !rel_diff_chiT_for_partials + + !latent_ddlnRho + !latent_ddlnT + + !log_P_face + !log_Ptrb + !log_cp_T_div_t_sound + + !QQ + + +!# Mass accretion + eps_grav ! -T*ds/dt (negative for expansion) + !log_abs_eps_grav_dm_div_L + !log_abs_v ! log10(abs(velocity)) (cm/s) + !log_mdot_cs + !log_mdot_v + eps_mdot + !env_eps_grav + !xm_div_delta_m + !log_xm_div_delta_m + + +!# Nuclear energy generation + !signed_log_eps_grav ! sign(eps_grav)*log10(max(1,abs(eps_grav))) + !signed_log_eps_nuc + net_nuclear_energy ! erg/gm/s from nuclear reactions minus all neutrino losses + ! The value plotted is net_nuclear_energy = sign(val)*log10(max(1,abs(val))) + ! where val = net nuclear energy minus all neutrino losses. + net_energy ! net_energy + eps_grav. + ! The value plotted is net_energy = sign(val)*log10(max(1,abs(val))) + ! where val = net nuclear energy plus eps_grav minus all neutrino losses. + eps_nuc_plus_nuc_neu + !eps_nuc_minus_non_nuc_neu + !eps_nuc_start + + eps_nuc ! ergs/g/sec from nuclear reactions (including losses to reaction neutrinos) + !log_abs_eps_nuc + !d_lnepsnuc_dlnd + !d_epsnuc_dlnd + !deps_dlnd_face + ! (was previously called deps_dlnd_at_face) + !d_lnepsnuc_dlnT + !d_epsnuc_dlnT + !deps_dlnT_face + ! (was previously called deps_dlnT_at_face) + !eps_nuc_neu_total ! erg/gm/sec as neutrinos from nuclear reactions + + non_nuc_neu ! non-nuclear-reaction neutrino losses + !nonnucneu_plas ! plasmon neutrinos (for collective reactions like gamma_plasmon => nu_e + nubar_e) + !nonnucneu_brem ! bremsstrahlung (for reactions like e- + (z,a) => e- + (z,a) + nu + nubar) + !nonnucneu_phot ! photon neutrinos (for reactions like e- + gamma => e- + nu_e + nubar_e) + !nonnucneu_pair ! pair production (for reactions like e+ + e- => nu_e + nubar_e) + !nonnucneu_reco ! recombination neutrinos (for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e) + + ! ergs/g/sec for reaction categories + add_reaction_categories ! this adds all the reaction categories + ! NOTE: you can list specific categories by giving their names (from chem_def) + pp + cno + tri_alpha + !c_alpha + !n_alpha + !o_alpha + !ne_alpha + !na_alpha + !mg_alpha + !si_alpha + !s_alpha + !ar_alpha + !ca_alpha + !ti_alpha + !fe_co_ni + !c12_c12 + !c12_o16 + !o16_o16 + !photo + !pnhe4 + !other + + ! adds columns for all of the reactions that are in the current net + ! Note that if using op_split_burn=.true. then zones which have been split will report 0 for thier rates + !add_raw_rates ! raw reaction rates, reactions/second + !add_screened_rates ! screened reaction rates reactions/second + !add_eps_nuc_rates ! Nuclear energy (minus neutrino losses) released erg/s + !add_eps_neu_rates ! Neutrino losses erg/s + + ! individual reactions (as many as desired) + ! use list_net_reactions = .true. in star_job to list all reactions in the current net + ! reactions/second + !raw_rate r_h1_h1_ec_h2 + !raw_rate r_h1_h1_wk_h2 + + !burn_num_iters ! Number of split_burn iterations taken + !burn_avg_epsnuc + !log_burn_avg_epsnuc + +!# Composition + !x_mass_fraction_H + !y_mass_fraction_He + !z_mass_fraction_metals + abar ! average atomic weight (g/mole) + !zbar ! average charge + !z2bar ! average charge^2 + ye ! average charge per baryon = proton fraction + + x ! hydrogen mass fraction + !log_x + y ! helium mass fraction + !log_y + z ! metallicity + !log_z ! metallicity + + add_abundances ! this adds all of the isos that are in the current net + ! NOTE: you can list specific isotopes by giving their names (from chem_def) + !h1 + !he3 + !he4 + !c12 + !n14 + !o16 + + !add_log_abundances ! this adds log10 of all of the isos that are in the current net + ! NOTE: you can list specific isotopes by giving their names (from chem_def) + !log h1 + !log he3 + !log he4 + !log c12 + !log n14 + !log o16 + + ! log concentration of species + ! concentration = number density / number density of electrons + ! Ci = (Xi/Ai) / sum(Zi*Xi/Ai) [see Thoul et al, ApJ 421:828-842, 1994] + !log_concentration h1 + !log_concentration he4 + + + ! typical charge for given species + ! (used by diffusion) + !typical_charge he4 + !typical_charge c12 + !typical_charge fe52 + + ! ionization state for given species + ! (same as typical charge, except that it's unsmoothed) + !ionization he4 + !ionization c12 + !ionization fe52 + + !cno_div_z ! abundance of c12, n14, and o16 as a fraction of total z + + + + +!# Opacity + !opacity ! opacity measured at center of zone + log_opacity ! log10(opacity) + !dkap_dlnrho_face ! partial derivative of opacity wrt. ln rho (at T=const) at outer edge of cell + ! (was previously called dkap_dlnrho_at_face) + !dkap_dlnT_face ! partial derivative of opacity wrt. ln T (at rho=const) at outer edge of cell + ! (was previously called dkap_dlnT_at_face) + !kap_frac_lowT ! fraction of opacity from lowT tables + !kap_frac_highT ! fraction of opacity from highT tables + !kap_frac_Type2 ! fraction of opacity from Type2 tables + !kap_frac_Compton ! fraction of opacity from Compton_Opacity + !kap_frac_op_mono ! fraction of opacity from OP mono + + !log_kap + !log_kap_times_factor + + !log_c_div_tau + !xtau + !xlogtau + !logtau_sub_xlogtau + +!# Luminosity + luminosity ! luminosity at outer boundary of zone (in Lsun units) + !logL ! log10(max(1d-2,L/Lsun)) + !log_Lrad + log_Ledd ! log10(Leddington/Lsun) -- local Ledd, 4 pi clight G m / kap + !log_L_div_Ledd ! log10(max(1d-12,L/Leddington)) + log_Lrad_div_Ledd + !log_Lrad_div_L + !signed_log_power ! sign(L)*log10(max(1,abs(L))) + + !lum_adv + !lum_conv + !lum_conv_MLT + !lum_div_Ledd + !lum_erg_s + !lum_plus_lum_adv + !lum_rad + + !log_L_div_CpTMdot + !log_abs_lum_erg_s + + !L + !Lc + !Lc_div_L + !Lr + !Lr_div_L + !Lt + !Lt_div_L + +!# Energetics + !total_energy ! specific total energy of cell (ergs/g). internal+potential+kinetic+rotation. + !cell_specific_IE + !cell_specific_KE + !cell_IE_div_IE_plus_KE + !cell_KE_div_IE_plus_KE + + !cell_ie_div_star_ie + !cell_internal_energy_fraction + !cell_internal_energy_fraction_start + !cell_specific_PE + + !log_cell_ie_div_star_ie + !log_cell_specific_IE + + ergs_eps_grav_plus_eps_mdot + ergs_error + !ergs_error_integral + ergs_mdot + !ergs_rel_error_integral + !dm_eps_grav + + !dE + + !etrb + !log_etrb + !extra_grav + log_rel_E_err + + !total_energy_sign + +!# Convection + !mlt_mixing_length ! mixing length for mlt (cm) + mlt_mixing_type ! value returned by mlt + !mlt_Pturb + !alpha_mlt + + !conv_vel ! convection velocity (cm/sec) + log_conv_vel ! log10 convection velocity (cm/sec) + + !conv_L_div_L + !log_conv_L_div_L + !lum_conv_div_lum_rad + !lum_rad_div_L_Edd + !lum_conv_div_lum_Edd + !lum_conv_div_L + !lum_rad_div_L + !lum_rad_div_L_Edd_sub_fourPrad_div_PchiT ! density increases outward if this is > 0 + ! see Joss, Salpeter, and Ostriker, "Critical Luminosity", ApJ 181:429-438, 1973. + + gradT ! mlt value for required temperature gradient dlnT/dlnP + + gradr ! dlnT/dlnP required for purely radiative transport + !grad_temperature ! smoothed dlnT/dlnP at cell boundary + !grad_density ! smoothed dlnRho/dlnP at cell boundary + + gradL ! gradient for Ledoux criterion for convection + !sch_stable ! 1 if grada > gradr, 0 otherwise + !ledoux_stable ! 1 if gradL > gradr, 0 otherwise + + !grada_sub_gradT + gradT_sub_grada ! gradT-grada at cell boundary + gradT_div_grada ! gradT/grada at cell boundary + + !gradr_sub_gradT + !gradT_sub_gradr ! gradT-gradr at cell boundary + !gradT_div_gradr ! gradT/gradr at cell boundary + + !log_gradT_div_gradr ! log10 gradT/gradr at cell boundary + log_mlt_Gamma ! convective efficiency + conv_vel_div_csound ! convection velocity divided by sound speed + !conv_vel_div_L_vel ! L_vel is velocity needed to carry L by convection; L = 4*pi*r^2*rho*vel**3 + log_mlt_D_mix ! log10 diffusion coefficient for mixing from mlt (cm^2/sec) + + !gradr_div_grada ! gradr/grada_face; > 1 => Schwarzschild unstable for convection + !gradr_sub_grada ! gradr - grada_face; > 0 => Schwarzschild unstable for convection + + !gradL_sub_gradr + !gradP_div_rho + !gradT_excess_effect + !gradT_rel_err + !gradT_sub_a + !grada_face + !grada_sub_gradr + !diff_grads + !log_diff_grads + + !mlt_D + !mlt_Gamma + !mlt_Y_face + !mlt_Zeta + !mlt_gradT + !mlt_log_abs_Y + !mlt_vc + !log_mlt_vc + !dvc_dt_TDC_div_g + + !superad_reduction_factor + !conv_vel_div_mlt_vc + + !log_Lconv + !log_Lconv_div_L + +!# Mixing + mixing_type ! mixing types are defined in mesa/const/public/const_def + log_D_mix ! log10 diffusion coefficient for mixing in units of cm^2/second (Eulerian) + !log_D_mix_non_rotation + !log_D_mix_rotation + + log_D_conv ! D_mix for regions where mix_type = convective_mixing + !log_D_leftover ! D_mix for regions where mix_type = leftover_convective_mixing + log_D_semi ! D_mix for regions where mix_type = semiconvective_mixing + log_D_ovr ! D_mix for regions where mix_type = overshoot_mixing + log_D_thrm ! D_mix for regions where mix_type = thermohaline_mixing + !log_D_minimum ! D_mix for regions where mix_type = minimum_mixing + !log_D_rayleigh_taylor ! D_mix for regions where mix_type = rayleigh_taylor_mixing + !log_D_anon ! D_mix for regions where mix_type = anonymous_mixing + !log_D_omega + + !log_sig_mix ! sig(k) is mixing flow across face k in (gm sec^1) + ! sig(k) = D_mix*(4*pi*r(k)**2*rho_face)**2/dmavg + + !dominant_isoA_for_thermohaline + !dominant_isoZ_for_thermohaline + !gradL_composition_term + + !mix_type + + + +!# Optical Depth + tau ! optical depth + !log_column_depth ! log10 column depth, exterior mass / area (g cm^-2) + !log_radial_depth ! log10 radial distance to surface (cm) + !logtau ! log10(optical depth) at cell face + !tau_eff ! tau that gives the local P == P_atm if this location at surface + ! tau_eff = kap*(P/g - Pextra_factor*(L/M)/(6*pi*clight*cgrav)) + !tau_eff_div_tau + + + +!# Rotation + omega ! angular velocity = j_rot/i_rot + !log_omega + log_j_rot + !log_J_div_M53 ! J is j*1e-15 integrated from center; M53 is m^(5/3) + log_J_inside ! J_inside is j_rot integrated from center + !shear ! -dlnomega/dlnR + !log_abs_shear ! log10(abs(dlnomega/dlnR)) + !richardson_number + i_rot ! specific moment of inertia at cell boundary + !j_rot ! specific angular momentum at cell boundary + !v_rot ! rotation velocity at cell boundary (km/sec) + !w_div_w_crit_roche !ratio of rotational velocity to keplerian at the equator + !without the contribution from the Eddington factor + fp_rot ! rotation factor for pressure + ft_rot ! rotation factor for temperature + !ft_rot_div_fp_rot ! gradr factor + + !log_am_nu_non_rot ! log10(am_nu_non_rot) + !log_am_nu_rot ! log10(am_nu_rot) + log_am_nu ! log10(am_nu_non_rot + am_nu_rot) + + !r_polar ! (Rsun) + !log_r_polar ! log10 (Rsun) + !r_equatorial ! (Rsun) + !log_r_equatorial ! log10 (Rsun) + !r_e_div_r_p ! equatorial/r_polar + !omega_crit ! breakup angular velocity = sqrt(G M / equatorial^3) + !omega_div_omega_crit + + !am_log_nu_omega ! for diffusion of omega + !am_log_nu_j ! for diffusion of angular momentum + + !am_log_nu_rot ! diffusion of angular momentum driven by rotation + !am_log_nu_non_rot ! diffusion driven by other sources, e.g. convection + + !am_log_sig_omega ! for diffusion of omega + !am_log_sig_j ! for diffusion of angular momentum + !am_log_sig ! == am_log_sig_omega + + !am_log_D_visc ! diffusion coeff for kinematic viscosity + !am_log_D_DSI ! diffusion coeff for dynamical shear instability + !am_log_D_SH ! diffusion coeff for Solberg-Hoiland instability + !am_log_D_SSI ! diffusion coeff for secular shear instability + !am_log_D_ES ! diffusion coeff for Eddington-Sweet circulation + !am_log_D_GSF ! diffusion coeff for Goldreich-Schubert-Fricke instability + !am_log_D_ST ! Spruit dynamo mixing diffusivity + !am_log_nu_ST ! Spruit dynamo effective viscosity + + !dynamo_log_B_r ! (Gauss) + !dynamo_log_B_phi ! (Gauss) + + !am_domega_dlnR + !log_abs_dlnR_domega + + !w_div_w_crit_roche2 + + +!# Diffusion + ! electric field from element diffusion calculation + !e_field + !log_e_field + + ! gravitational field from element diffusion calculation + !g_field_element_diffusion + !log_g_field_element_diffusion + + !eE_div_mg_element_diffusion + !log_eE_div_mg_element_diffusion + + ! element diffusion velocity for species + !edv h1 + !edv he4 + !edv o16 + + ! Energy generated by Ne22 sedimentation. + !eps_WD_sedimentation + !log_eps_WD_sedimentation + + !eps_diffusion + !log_eps_diffusion + + !diffusion_D h1 ! self diffusion coeff + !diffusion_dX h1 ! change in h1 mass fraction from diffusion + !diffusion_dX he4 ! change in he4 mass fraction from diffusion + !diffusion_dX n20 ! change in n20 mass fraction from diffusion + + !v_rad h1 ! velocity from radiative levitation + !v_rad he4 ! velocity from radiative levitation + !v_rad ne20 ! velocity from radiative levitation + + !log_g_rad h1 ! log10 acceleration from radiative levitation + !log_g_rad he4 ! log10 acceleration from radiative levitation + !log_g_rad ne20 ! log10 acceleration from radiative levitation + +!# Phase Separation + !eps_phase_separation + +!# Oscillations + brunt_N2 ! brunt-vaisala frequency squared + brunt_N2_structure_term + brunt_N2_composition_term + !log_brunt_N2_structure_term + !log_brunt_N2_composition_term + !brunt_A ! = N^2*r/g + !brunt_A_div_x2 ! x = r(k)/r(1) + !brunt_N2_dimensionless ! N2 in units of 3GM/R^3 + !brunt_N_dimensionless ! N in units of sqrt(3GM/R^3) + !brunt_frequency ! cycles per day + !brunt_N ! sqrt(abs(brunt_N2)) + !log_brunt_N ! log10(brunt_N) + !log_brunt_N2 ! log10(brunt_N2) + !log_brunt_N2_dimensionless ! log10(brunt_N2_dimensionless) + + !brunt_B ! smoothed numerical difference + !brunt_nonB ! = grada - gradT + !log_brunt_B ! smoothed numerical difference + !log_brunt_nonB ! = grada - gradT + + !sign_brunt_N2 ! sign of brunt_N2 (+1 for Ledoux stable; -1 for Ledoux unstable) + !brunt_nu ! brunt_frequency in microHz + !log_brunt_nu ! brunt_frequency in microHz + + !lamb_S ! lamb frequency for l=1: S = sqrt(2)*csound/r (rad/s) + !lamb_S2 ! squared lamb frequency for l=1: S2 = 2*(csound/r)^2 (rad^2/s^2) + + !lamb_Sl1 ! lamb frequency for l=1; = sqrt(2)*csound/r (microHz) + !lamb_Sl2 ! lamb frequency for l=2; = sqrt(6)*csound/r (microHz) + !lamb_Sl3 ! lamb frequency for l=3; = sqrt(12)*csound/r (microHz) + !lamb_Sl10 ! lamb frequency for l=10; = sqrt(110)*csound/r (microHz) + + !log_lamb_Sl1 ! log10(lamb_Sl1) + !log_lamb_Sl2 ! log10(lamb_Sl2) + !log_lamb_Sl3 ! log10(lamb_Sl3) + !log_lamb_Sl10 ! log10(lamb_Sl10) + + !brunt_N_div_r_integral ! integral from center of N*dr/r + !k_r_integral ! integral from center of k_r*dr + !brunt_N2_sub_omega2 + !sl2_sub_omega2 + + +!# RSP + + !rsp_Chi ! dlnP_dlnRho + !rsp_Et ! Specific turbulent energy + !rsp_logEt ! Log specific turbulent energy + !rsp_erad ! Specific internal (radiative) energy + !rsp_log_erad ! Log specific internal (radiative) energy + !rsp_Hp_face ! Pressure scale height at cell face + !rsp_Lc ! Convective luminosity + !rsp_Lc_div_L ! Convective luminosity div total luminosity + !rsp_Lr ! Radiative luminosity + !rsp_Lr_div_L ! Radiative luminosity div total luminosity + !rsp_Lt ! Turbulent luminosity + !rsp_Lt_div_L ! Turbulent luminosity div total luminosity + !rsp_Pt ! Turbulent pressure, p_t, see Table 1 in MESA5 + !rsp_Uq ! Viscous momentum transfer rate, U_q, see Table 1 in MESA5 + !rsp_Eq ! Viscous energy transfer rate, epsilon_q, see Table 1 in MESA5 + !rsp_Pvsc ! Artificial viscosity, p_av, see Table 1 in MESA5 + !rsp_gradT ! Temperature gradient + !rsp_Y_face ! Superadiabatic gradient at cell face, Y_sag, see Table 1 in MESA5 + !rsp_damp ! Turbulent dissipation, D, see Table 1 in MESA5 + !rsp_dampR ! Radiative cooling, D_r, see Table 1 in MESA5 + !rsp_sink ! Sum of turbulent dissipation and radiative cooling terms + !rsp_src ! Source function, S, see Table 1 in MESA5 + !rsp_src_snk ! Convective coupling, C, see Table 1 in MESA5 + !rsp_heat_exchange_timescale ! 1d0/(clight * opacity * density) + !rsp_log_heat_exchange_timescale + !rsp_log_dt_div_heat_exchange_timescale ! Ratio of time step to heat exchange timescale + !w + !log_w + + !COUPL + !DAMP + !DAMPR + !SOURCE + !Chi + !Eq + !Hp_face + !PII_face + !Ptrb + !Pvsc + !Uq + !Y_face + +!# RTI + + !RTI_du_diffusion_kick + !alpha_RTI + !boost_for_eta_RTI + !dedt_RTI + !dudt_RTI + !eta_RTI + !log_alpha_RTI + !log_boost_for_eta_RTI + !log_eta_RTI + !log_etamid_RTI + !log_lambda_RTI_div_Hrho + !log_sig_RTI + !log_sigmid_RTI + !log_source_RTI + !log_source_minus_alpha_RTI + !log_source_plus_alpha_RTI + !source_minus_alpha_RTI + !source_plus_alpha_RTI + !lambda_RTI + +!# Hydrodynamics + + + !v + !v_div_v_escape + !v_div_vesc + !v_kms + !log_v_escape + + !u + !u_face + + !P_face + + +!# Extras + !extra_heat + !extra_L ! extra_heat integrated from center (Lsun) + !log_extra_L ! log10 integrated from center (Lsun) + !log_irradiation_heat + + !extra_jdot ! set in other_torque routine + !extra_omegadot ! set in other_torque routine + + !extra_opacity_factor ! set in other_opacity_factor routine + + ! diffusion factor profile for species, set in other_diffusion_factor routine + !extra_diffusion_factor h1 + !extra_diffusion_factor he4 + !extra_diffusion_factor o16 + + + +!# Miscellaneous + + !dlog_h1_dlogP ! (log(h1(k)) - log(h1(k-1)))/(log(P(k)) - log(P(k-1))) + !dlog_he3_dlogP + !dlog_he4_dlogP + !dlog_c12_dlogP + !dlog_c13_dlogP + !dlog_n14_dlogP + !dlog_o16_dlogP + !dlog_ne20_dlogP + !dlog_mg24_dlogP + !dlog_si28_dlogP + + !dlog_pp_dlogP + !dlog_cno_dlogP + !dlog_3alf_dlogP + + !dlog_burn_c_dlogP + !dlog_burn_n_dlogP + !dlog_burn_o_dlogP + + !dlog_burn_ne_dlogP + !dlog_burn_na_dlogP + !dlog_burn_mg_dlogP + + !dlog_cc_dlogP + !dlog_co_dlogP + !dlog_oo_dlogP + + !dlog_burn_si_dlogP + !dlog_burn_s_dlogP + !dlog_burn_ar_dlogP + !dlog_burn_ca_dlogP + !dlog_burn_ti_dlogP + !dlog_burn_cr_dlogP + !dlog_burn_fe_dlogP + + !dlog_pnhe4_dlogP + !dlog_photo_dlogP + !dlog_other_dlogP + + !logR_kap ! logR = logRho - 3*logT + 18 ; used in kap tables + !logW ! logW = logPgas - 4*logT + !logQ ! logQ = logRho - 2*logT + 12 + !logV ! logV = logRho - 0.7*logE + 20 + + !log_CpT_absMdot_div_L ! log10(s% Cp(k)*s% T(k)*abs(s% mstar_dot)/s% L(k)) + + !delta_r ! r - r_start, change during step + !delta_L ! L - L_start, change during step + !delta_cell_vol ! cell_vol - cell_vol_start, change during step + !delta_entropy ! entropy - entropy_start, change during step (does not include effects of diffusion) + !delta_T ! T - T_start, change during step + !delta_rho ! rho - rho_start, change during step + !delta_eps_nuc ! eps_nuc - eps_nuc_start, change during step + !delta_mu ! mu - mu_start, change during step + + !zFe ! mass fraction of "Fe" = Fe+Co+Ni + !log_zFe + !dPdr_dRhodr_info + !log_sig_raw_mix + + !d_u_div_rmid + !d_u_div_rmid_start + !d_v_div_r_dm + !d_v_div_r_dr + + !dlnP_dlnR + !dlnRho_dlnR + !dlnRho_dr + !dlnX_dr + !dlnY_dr + !dlogR + !dPdr_div_grav + !dPdr_info + !dRhodr_info + !dRstar_div_dr + !dr_ratio + !dm_eps_grav + !dr_ratio + !dt_cs_div_dr + !dt_div_tau_conv + !dt_times_conv_vel_div_mixing_length + !log_dt_cs_div_dr + !log_dt_div_tau_conv + !log_dt_times_conv_vel_div_mixing_length + !log_du_kick_div_du + !du + !dvdt_dPdm + !dvdt_grav + + !tau_conv + !tau_cool + !tau_epsnuc + !tau_qhse + + !max_abs_xa_corr + + !tdc_num_iters + + !k + + +! the first few lines of the profile contain general info about the model. +! for completeness, those items are described here. + + ! initial mass and Z + ! initial_mass + ! initial_z + ! general properties of the current state + ! model_number + ! num_zones + ! star_age + ! time_step + ! properties at the photosphere + ! Teff + ! photosphere_L + ! photosphere_r + ! properties at the outermost zone of the model + ! log_surface_L + ! log_surface_radius + ! log_surface_temp + ! properties near the center of the model + ! log_center_temp + ! log_center_density + ! log_center_P + ! center_eta + ! abundances near the center + ! center_h1 + ! center_he3 + ! center_he4 + ! center_c12 + ! center_n14 + ! center_o16 + ! center_ne20 + ! information about total mass + ! star_mass + ! star_mdot + ! star_mass_h1 + ! star_mass_he3 + ! star_mass_he4 + ! star_mass_c12 + ! star_mass_n14 + ! star_mass_o16 + ! star_mass_ne20 + ! locations of abundance transitions + ! he_core_mass + ! c_core_mass + ! o_core_mass + ! si_core_mass + ! fe_core_mass + ! location of optical depths 10 and 100 + ! tau10_mass + ! tau10_radius + ! tau100_mass + ! tau100_radius + ! time scales + ! dynamic_time + ! kh_timescale + ! nuc_timescale + ! various kinds of total power + ! power_nuc_burn + ! power_h_burn + ! power_he_burn + ! power_neu + ! a few control parameter values + ! h1_boundary_limit + ! he4_boundary_limit + ! c12_boundary_limit + ! burn_min1 + ! burn_min2 diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/re b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/re new file mode 100755 index 000000000..c9ef26f96 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/re @@ -0,0 +1,33 @@ +#!/bin/bash + +shopt -u expand_aliases + +photo_directory=photos + +function most_recent_photo { + ls -tp "$photo_directory" | grep -v / | head -1 +} + +if [ $# -eq 0 ] +then + photo=$(most_recent_photo) +else + photo=$1 +fi + +if [ -z "$photo" ] || ! [ -f "$photo_directory/$photo" ] +then + echo "specified photo ($photo) does not exist" + exit 1 +fi + +echo "restart from $photo" +if ! cp "$photo_directory/$photo" restart_photo +then + echo "failed to copy photo ($photo)" + exit 1 +fi + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn new file mode 100755 index 000000000..14c5878ce --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn @@ -0,0 +1,25 @@ +#!/bin/bash + +# uncomment the following line to skip the optional inlists +# MESA_SKIP_OPTIONAL=t + +# this provides the definition of do_one (run one part of test) +# do_one [inlist] [output model] [LOGS directory] +source "${MESA_DIR}/star/test_suite/test_suite_helpers" + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +# check if can skip building starting model +#if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_he_dep.mod he_dep.mod +#else +# do_one inlist_to_he_dep_header he_dep.mod +# cp he_dep.mod standard_he_dep.mod +#fi + +do_one inlist_pulses_header final.mod + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +echo 'finished' + diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn1 b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn1 new file mode 100755 index 000000000..25590040a --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/rn1 @@ -0,0 +1,7 @@ +#!/bin/bash + +rm -f restart_photo + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run.f90 b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run.f90 new file mode 100644 index 000000000..112660f7c --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run.f90 @@ -0,0 +1,15 @@ +program run + use run_star_support, only: do_read_star_job + use run_star, only: do_run_star + + implicit none + + integer :: ierr + + ierr = 0 + call do_read_star_job('inlist', ierr) + if (ierr /= 0) stop 1 + + call do_run_star + +end program diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run_star_extras.f90 b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run_star_extras.f90 new file mode 100644 index 000000000..2b9d1a9ad --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/src/run_star_extras.f90 @@ -0,0 +1,1498 @@ +! *********************************************************************** +! +! Copyright (C) 2010-2019 The MESA Team +! +! this file is part of mesa. +! +! mesa is free software; you can redistribute it and/or modify +! it under the terms of the gnu general library public license as published +! by the free software foundation; either version 2 of the license, or +! (at your option) any later version. +! +! mesa is distributed in the hope that it will be useful, +! but without any warranty; without even the implied warranty of +! merchantability or fitness for a particular purpose. see the +! gnu library general public license for more details. +! +! you should have received a copy of the gnu library general public license +! along with this software; if not, write to the free software +! foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa +! +! *********************************************************************** + + module run_star_extras + + use star_lib + use star_def + use const_def + use math_lib + use auto_diff + use chem_def + use utils_lib + use rates_def, only: i_rate + + use interp_1d_def, only: pm_work_size + use interp_1d_lib, only: interp_pm, interp_values, interp_value + + implicit none + + include "test_suite_extras_def.inc" + + logical :: dbg = .false. + + logical :: pgstar_flag + integer :: pgstar_interval, pgstar_file_interval + + !!!!!!!!!!!!!!!!!!!!!!!!! + ! here are definitions for the indices of all s% xtra, s% lxtra and s% ixtra variables used in the run + !!!!!!!!!!!!!!!!!!!!!!!!! + + ! lxtra(lx_hydro_on) is true while hydro is on + integer, parameter :: lx_hydro_on = 1 + ! lxtra(lx_hydro_has_been_on) is true after hydro is turned on for the first time + integer, parameter :: lx_hydro_has_been_on = 2 + ! lxtra(lx_have_saved_gamma_prof) is true, after a profile is saved when gamma_integral=0 + ! One profile is saved during each hydro phase, this is turned + ! back to false after a relax + integer, parameter :: lx_have_saved_gamma_prof = 3 + ! lx_have_reached_gamma_limit is same as before, but it's kept as .true. after the first time + ! gamma_integral=0. Used to count time from first pulse. + integer, parameter :: lx_have_reached_gamma_limit = 4 + ! lxtra(lx_he_zams) indicates if star has reached he zams + integer, parameter :: lx_he_zams = 5 + + ! xtra(x_time_start_pulse) contains the time at which a pulse starts + integer, parameter :: x_time_start_pulse = 1 + ! xtra(x_dyn_time) contains the estimate for the dynamical time when that happens + integer, parameter :: x_dyn_time = 2 + ! xtra(x_gamma_int_bound) contains gamma_integral in bound regions + integer, parameter :: x_gamma_int_bound = 3 + ! xtra(x_time_since_first_gamma_zero) contains time in seconds since first time gamma_integral=0 + integer, parameter :: x_time_since_first_gamma_zero = 4 + ! xtra(x_star_age_at_relax) Stores s% star_age at the point the last relax was made + integer, parameter :: x_star_age_at_relax = 8 + + ! ixtra(ix_steps_met_relax_cond) counts the steps during which the conditions to relax a model are met + integer, parameter :: ix_steps_met_relax_cond = 1 + ! ixtra(ix_num_relaxations) counts the number of times the star has relaxed + integer, parameter :: ix_num_relaxations = 2 + ! ixtra(ix_steps_since_relax) counts the number of steps since last relax + integer, parameter :: ix_steps_since_relax = 3 + ! ixtra(ix_steps_since_hydro_on) counts the number of steps since hydro was turned on + integer, parameter :: ix_steps_since_hydro_on = 4 + + + ! xtra(x_vsurf_at_remove_surface) Store velocity the first time remove_surface is used + ! after that point we fix the surface velocity to this + integer, parameter :: x_vsurf_at_remove_surface = 7 + + ! lxtra(lx_using_bb_bcs) indicates if black body BCs are being used + integer, parameter :: lx_using_bb_bcs = 7 + + ! xtra(x_tau_factor) stores tau_factor which is set by remove_surface + integer, parameter :: x_tau_factor = 6 + + ! xtra(x_Tsurf_for_atm) Surface T for fixed Tsurf atm option + ! this is used when layers at large radii are removed + integer, parameter :: x_Tsurf_for_atm = 5 + + !!!!!!!!!!!!!!!!!!!!!!!!! + ! These variables are loaded up from x_ctrl, x_integer_ctrl and x_logical_ctrl + ! values specified on inlist_ppisn + !!!!!!!!!!!!!!!!!!!!!!!!! + + real(dp) :: min_gamma_sub_43_for_hydro + real(dp) :: max_v_for_pulse, q_for_dyn_ts, num_dyn_ts_for_relax + real(dp) :: q_for_relax_check, max_v_for_relax, max_machn_for_relax, & + max_Lneu_for_relax, max_Lnuc_for_relax + integer :: num_steps_before_relax + logical :: remove_extended_layers !, in_inlist_pulses + logical :: in_inlist_pulses + real(dp) :: max_dt_before_pulse + real(dp) :: max_Lneu_for_mass_loss + real(dp) :: delta_lgLnuc_limit, max_Lphoto_for_lgLnuc_limit, max_Lphoto_for_lgLnuc_limit2 + real(dp) :: delta_lgRho_cntr_hard_limit, dt_div_min_dr_div_cs_limit + real(dp) :: logT_for_v_flag, logLneu_for_v_flag + logical :: stop_100d_after_pulse + + real(dp) :: max_dt_during_pulse + real(dp) :: vsurf_for_fixed_bc + + contains + + include "test_suite_extras.inc" + + + subroutine extras_controls(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + s% extras_startup => extras_startup + s% extras_start_step => extras_start_step + s% extras_check_model => extras_check_model + s% extras_finish_step => extras_finish_step + s% extras_after_evolve => extras_after_evolve + s% how_many_extra_history_columns => how_many_extra_history_columns + s% data_for_extra_history_columns => data_for_extra_history_columns + s% how_many_extra_profile_columns => how_many_extra_profile_columns + s% data_for_extra_profile_columns => data_for_extra_profile_columns + + ! we turn pgstar on and off at some points, so we store the original setting + pgstar_flag = s% job% pgstar_flag + pgstar_interval = s% pg% pgstar_interval + pgstar_file_interval = s% pg% Grid2_file_interval + + ! this is optional + s% other_wind => brott_wind + s% other_adjust_mdot => my_adjust_mdot + s% other_before_struct_burn_mix => my_before_struct_burn_mix + s% other_eval_fp_ft => my_other_eval_fp_ft + s% other_kap_get => my_other_kap_get + + ! store user provided options from the inlist + min_gamma_sub_43_for_hydro = s% x_ctrl(1) + max_v_for_pulse = s% x_ctrl(2) + q_for_dyn_ts = s% x_ctrl(3) + num_dyn_ts_for_relax = s% x_ctrl(4) + q_for_relax_check = s% x_ctrl(5) + max_v_for_relax = s% x_ctrl(6) + max_machn_for_relax = s% x_ctrl(7) + max_Lneu_for_relax = s% x_ctrl(8) + max_Lnuc_for_relax = s% x_ctrl(9) + num_steps_before_relax = s% x_integer_ctrl(1) + in_inlist_pulses = s% x_logical_ctrl(2) + max_dt_before_pulse = s% x_ctrl(10) + max_Lneu_for_mass_loss = s% x_ctrl(11) + delta_lgLnuc_limit = s% x_ctrl(12) + max_Lphoto_for_lgLnuc_limit = s% x_ctrl(13) + max_Lphoto_for_lgLnuc_limit2 = s% x_ctrl(14) + logT_for_v_flag = s% x_ctrl(15) + logLneu_for_v_flag = s% x_ctrl(16) + stop_100d_after_pulse = s% x_logical_ctrl(1) + remove_extended_layers = s% x_logical_ctrl(2) + max_dt_during_pulse = s% x_ctrl(18) + vsurf_for_fixed_bc = s% x_ctrl(19) + + ! we store the value given in inlist_ppisn and deactivate it at + ! high T + delta_lgRho_cntr_hard_limit = s% delta_lgRho_cntr_hard_limit + ! we also store dt_div_min_dr_div_cs_limit, we keep it at a + ! high value until the onset of a pulse to prevent unnecesarily + ! small timesteps before a pulsation + dt_div_min_dr_div_cs_limit = s% dt_div_min_dr_div_cs_limit + + end subroutine extras_controls + + subroutine brott_wind(id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z, w, ierr) + use star_def + integer, intent(in) :: id + real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z ! surface values (cgs) + ! NOTE: surface is outermost cell. not necessarily at photosphere. + ! NOTE: don't assume that vars are set at this point. + ! so if you want values other than those given as args, + ! you should use values from s% xh(:,:) and s% xa(:,:) only. + ! rather than things like s% Teff or s% lnT(:) which have not been set yet. + real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0) + integer, intent(out) :: ierr + + integer :: h1, he4 + real(dp) :: Xs, Ys, Z_div_Z_solar, Teff_jump, alfa, L1, M1, R1, T1, & + vink_wind, nieu_wind, hamann_wind + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + L1 = Lsurf + M1 = Msurf + R1 = Rsurf + T1 = Tsurf + + h1 = s% net_iso(ih1) + he4 = s% net_iso(ihe4) + Xs = s% xa(h1,1) + Ys = s% xa(he4,1) + ! Z=0.0142 is Z from Asplund et al. 2009 + Z_div_Z_solar = s% kap_rq% Zbase/0.0142d0 + ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature + Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10(Z_div_Z_solar))) + + vink_wind = 0d0 + nieu_wind = 0d0 + hamann_wind = 0d0 + w = 0 + + call eval_Vink_wind(vink_wind) + call eval_Nieuwenhuijzen_wind(nieu_wind) + call eval_Hamann_wind(hamann_wind) + + ! use 1/10 hamann + hamann_wind = hamann_wind/10d0 + + if (T1 < Teff_jump) then + ! low T wind + w = max(vink_wind, nieu_wind) + else + ! high T wind + alfa = 0d0 + if (Xs > 0.7d0) then + alfa = 1d0 + else if (Xs > 0.4d0 .and. Xs < 0.7d0) then + alfa = (Xs - 0.4d0)/0.3d0 + end if + w = alfa * vink_wind + (1d0-alfa) * hamann_wind + end if + + ierr = 0 + + contains + + subroutine eval_Vink_wind(w) + real(dp), intent(inout) :: w + real(dp) :: alfa, w1, w2, logMdot, dT, vinf_div_vesc + + ! alfa = 1 for hot side, = 0 for cool side + if (T1 > 27500d0) then + alfa = 1 + else if (T1 < 22500d0) then + alfa = 0 + else + dT = 100d0 + if (T1 > Teff_jump + dT) then + alfa = 1 + else if (T1 < Teff_jump - dT) then + alfa = 0 + else + alfa = (T1 - (Teff_jump - dT)) / (2*dT) + end if + end if + + if (alfa > 0) then ! eval hot side wind (eqn 24) + vinf_div_vesc = 2.6d0 ! this is the hot side galactic value + vinf_div_vesc = vinf_div_vesc*pow(Z_div_Z_solar,0.13d0) ! corrected for Z + logMdot = & + - 6.697d0 & + + 2.194d0*log10(L1/Lsun/1d5) & + - 1.313d0*log10(M1/Msun/30) & + - 1.226d0*log10(vinf_div_vesc/2d0) & + + 0.933d0*log10(T1/4d4) & + - 10.92d0*pow2(log10(T1/4d4)) & + + 0.85d0*log10(Z_div_Z_solar) + w1 = exp10(logMdot) + else + w1 = 0 + end if + + if (alfa < 1) then ! eval cool side wind (eqn 25) + vinf_div_vesc = 1.3d0 ! this is the cool side galactic value + vinf_div_vesc = vinf_div_vesc*pow(Z_div_Z_solar,0.13d0) ! corrected for Z + logMdot = & + - 6.688d0 & + + 2.210d0*log10(L1/Lsun/1d5) & + - 1.339d0*log10(M1/Msun/30) & + - 1.601d0*log10(vinf_div_vesc/2d0) & + + 1.07d0*log10(T1/2d4) & + + 0.85d0*log10(Z_div_Z_solar) + w2 = exp10(logMdot) + else + w2 = 0 + end if + + w = alfa*w1 + (1 - alfa)*w2 + + end subroutine eval_Vink_wind + + subroutine eval_Nieuwenhuijzen_wind(w) + ! Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134 (eqn 2) + real(dp), intent(out) :: w + real(dp) :: log10w + include 'formats' + log10w = -14.02d0 & + +1.24d0*log10(L1/Lsun) & + +0.16d0*log10(M1/Msun) & + +0.81d0*log10(R1/Rsun) & + +0.85d0*log10(Z_div_Z_solar) + w = exp10(log10w) + end subroutine eval_Nieuwenhuijzen_wind + + subroutine eval_Hamann_wind(w) + ! Hamann, W.-R.; Koesterke, L.; Wessolowski, U. 1995, A&A, 299, 151 + real(dp), intent(out) :: w + real(dp) :: log10w + include 'formats' + log10w = -11.95d0 & + +1.5d0*log10(L1/Lsun) & + -2.85d0*Xs & + + 0.85d0*log10(Z_div_Z_solar) + w = exp10(log10w) + end subroutine eval_Hamann_wind + + end subroutine brott_wind + + subroutine my_adjust_mdot(id, ierr) + use star_def + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + real(dp) :: Lrad_div_Ledd + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + if (s% generations > 2) then + write(*,*) "check mdots", s% mstar_dot, s% mstar_dot_old + if (abs(s% mstar_dot) > 1.05d0*abs(s% mstar_dot_old)) then + s% mstar_dot = 1.05d0*s% mstar_dot_old + else if (abs(s% mstar_dot) < 0.95d0*abs(s% mstar_dot_old)) then + s% mstar_dot = 0.95d0*s% mstar_dot_old + end if + end if + end subroutine my_adjust_mdot + + subroutine my_other_eval_fp_ft( & + id, nz, xm, r, rho, aw, ft, fp, r_polar, r_equatorial, report_ierr, ierr) + use num_lib + integer, intent(in) :: id + integer, intent(in) :: nz + real(dp), intent(in) :: aw(:), r(:), rho(:), xm(:) ! (nz) + type(auto_diff_real_star_order1), intent(out) :: ft(:), fp(:) ! (nz) + real(dp), intent(inout) :: r_polar(:), r_equatorial(:) ! (nz) + logical, intent(in) :: report_ierr + integer, intent(out) :: ierr + + type (star_info), pointer :: s + integer :: j + real(dp) :: alpha + + type(auto_diff_real_1var_order1) :: A_omega,fp_numerator, ft_numerator, w, w2, w3, w4, w5, w6, lg_one_sub_w4, & + fp_temp, ft_temp + + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + +!$OMP PARALLEL DO PRIVATE(j, A_omega, fp_numerator, ft_numerator, fp_temp, ft_temp, w, w2, w3, w4, w5, w6, lg_one_sub_w4) SCHEDULE(dynamic,2) + do j=1, s% nz + !Compute fp, ft, re and rp using fits to the Roche geometry of a single star. + !by this point in the code, w_div_w_crit_roche is set + w = s% w_div_w_crit_roche(j) + w%d1val1 = 1d0 + + w2 = pow2(w) + w4 = pow4(w) + w6 = pow6(w) + lg_one_sub_w4 = log(1d0-w4) + A_omega = (1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4) + fp_numerator = (1d0-two_thirds*w2-0.06837d0*w4-0.2495d0*w6) + ft_numerator = (1d0+0.2185d0*w4-0.1109d0*w6) + !fits for fp, ft + fp_temp = fp_numerator/A_omega + ft_temp = ft_numerator/A_omega + !re and rp can be derived analytically from w_div_wcrit + r_equatorial(j) = r(j)*(1d0+w2% val/6d0-0.0002507d0*w4% val+0.06075d0*w6% val) + r_polar(j) = r_equatorial(j)/(1d0+0.5d0*w2% val) + ! Be sure they are consistent with r_Phi + r_equatorial(j) = max(r_equatorial(j),r(j)) + r_polar(j) = min(r_polar(j),r(j)) + + fp(j) = 0d0 + ft(j) = 0d0 + fp(j)% val = fp_temp% val + ft(j)% val = ft_temp% val + if (s% w_div_wc_flag) then + fp(j)% d1Array(i_w_div_wc_00) = fp_temp% d1val1 + ft(j)% d1Array(i_w_div_wc_00) = ft_temp% d1val1 + end if + end do +!$OMP END PARALLEL DO + + if (s% u_flag) then + !make fp and ft 1 in the outer 0.001 mass fraction of the star. softly turn to zero from the outer 0.002 + do j=1, s% nz + if (s% q(j) > 0.999) then + fp(j) = 1d0 + ft(j) = 1d0 + else if (s% q(j) > 0.998) then + alpha = (1d0-(s% q(j)-0.998)/(0.001)) + fp(j) = fp(j)*alpha + 1d0*(1-alpha) + ft(j) = ft(j)*alpha + 1d0*(1-alpha) + end if + end do + end if + + end subroutine my_other_eval_fp_ft + + + subroutine my_other_kap_get( & + id, k, handle, species, chem_id, net_iso, xa, & + log10_rho, log10_T, & + lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, & + eta, d_eta_dlnRho, d_eta_dlnT, & + kap_fracs, kap, dln_kap_dlnRho, dln_kap_dlnT, dln_kap_dxa, ierr) + + use kap_def, only: num_kap_fracs + use kap_lib + + ! INPUT + integer, intent(in) :: id ! star id if available; 0 otherwise + integer, intent(in) :: k ! cell number or 0 if not for a particular cell + integer, intent(in) :: handle ! from alloc_kap_handle + integer, intent(in) :: species + integer, pointer :: chem_id(:) ! maps species to chem id + ! index from 1 to species + ! value is between 1 and num_chem_isos + integer, pointer :: net_iso(:) ! maps chem id to species number + ! index from 1 to num_chem_isos (defined in chem_def) + ! value is 0 if the iso is not in the current net + ! else is value between 1 and number of species in current net + real(dp), intent(in) :: xa(:) ! mass fractions + real(dp), intent(in) :: log10_rho ! density + real(dp), intent(in) :: log10_T ! temperature + real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT + ! free_e := total combined number per nucleon of free electrons and positrons + real(dp), intent(in) :: eta, d_eta_dlnRho, d_eta_dlnT + ! eta := electron degeneracy parameter + + ! OUTPUT + real(dp), intent(out) :: kap_fracs(num_kap_fracs) + real(dp), intent(out) :: kap ! opacity + real(dp), intent(out) :: dln_kap_dlnRho ! partial derivative at constant T + real(dp), intent(out) :: dln_kap_dlnT ! partial derivative at constant Rho + real(dp), intent(out) :: dln_kap_dxa(:) ! partial derivative w.r.t. to species + integer, intent(out) :: ierr ! 0 means AOK. + + type (star_info), pointer :: s + real(dp) :: velocity + real(dp) :: radius, logR + real(dp) :: logT_alt, inv_diff + real(dp) :: log_kap, alpha + + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + kap = 0; dln_kap_dlnRho = 0; dln_kap_dlnT = 0; dln_kap_dxa = 0 + velocity = 0 + radius = 0 + + !if (k==1 .and. s% u_flag .and. .not. is_nan(s% lnR_start(1))) then !very surface cell can go mad, things are more stable if we fix opacity + ! if (s% xh_start(s% i_u,1)>sqrt(2*s% cgrav(1)*s% m(1)/exp(s% lnR_start(1)))) then + if (k==1 .and. s% u_flag) then !very surface cell can go mad, things are more stable if we fix opacity + ! this is to support restarts, as xh_start and r_start are + ! not properly set when model loads + if (s% solver_iter > 0) then + velocity = s% xh_start(s% i_u,1) + radius = s% r_start(1) + else + velocity = s% xh(s% i_u,1) + radius = s% r(1) + end if + if (velocity>sqrt(2*s% cgrav(1)*s% m(1)/radius)) then + kap = 0.2d0*(1 + s% X(1)) + dln_kap_dlnRho = 0d0 + dln_kap_dlnT = 0d0 + return + else + call kap_get( & + s% kap_handle, species, chem_id, net_iso, xa, & + log10_rho, log10_T, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, & + eta, d_eta_dlnRho, d_eta_dlnT, & + kap_fracs, kap, dln_kap_dlnRho, dln_kap_dlnT, dln_kap_dxa, ierr) + end if + else + call kap_get( & + s% kap_handle, species, chem_id, net_iso, xa, & + log10_rho, log10_T, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, & + eta, d_eta_dlnRho, d_eta_dlnT, & + kap_fracs, kap, dln_kap_dlnRho, dln_kap_dlnT, dln_kap_dxa, ierr) + end if + + + end subroutine my_other_kap_get + + subroutine extras_startup(id, restart, ierr) + integer, intent(in) :: id + logical, intent(in) :: restart + integer, intent(out) :: ierr + type (star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + call test_suite_startup(s, restart, ierr) + if (.not. restart) then + s% lxtra(lx_hydro_on) = .false. + s% lxtra(lx_hydro_has_been_on) = .false. + s% lxtra(lx_have_saved_gamma_prof) = .false. + s% lxtra(lx_have_reached_gamma_limit) = .false. + s% lxtra(lx_he_zams) = .false. + s% xtra(x_time_start_pulse) = -1d0 + s% xtra(x_dyn_time) = -1d0 + s% xtra(x_gamma_int_bound) = -1d0 + s% xtra(x_time_since_first_gamma_zero) = 0d0 + s% ixtra(ix_steps_met_relax_cond) = 0 + s% ixtra(ix_num_relaxations) = 0 + s% ixtra(ix_steps_since_relax) = 0 + s% ixtra(ix_steps_since_hydro_on) = 0 + s% xtra(x_star_age_at_relax) = -1d0 + + s% xtra(x_vsurf_at_remove_surface) = 0d0 + s% lxtra(lx_using_bb_bcs) = .false. + s% xtra(x_tau_factor) = 0d0 + s% xtra(x_Tsurf_for_atm) = 0d0 + + + ! to avoid showing pgstar stuff during initial model creation + s% pg% pgstar_interval = 100000000 + s% pg% Grid2_file_interval = 100000000 + else ! it is a restart + if (s% lxtra(lx_hydro_on)) then + call star_read_controls(id, 'inlist_hydro_on', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + call star_set_u_flag(id, .true., ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + else if (s% lxtra(lx_hydro_has_been_on)) then + call star_read_controls(id, 'inlist_hydro_off', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + end if + if (s% lxtra(lx_hydro_has_been_on)) then + call star_read_controls(id, 'inlist_after_first_pulse', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + end if + end if + end subroutine extras_startup + + + subroutine extras_after_evolve(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + real(dp) :: dt + character (len=strlen) :: test + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + call test_suite_after_evolve(s, ierr) + end subroutine extras_after_evolve + + + ! returns either keep_going, retry, or terminate. + integer function extras_check_model(id) + integer, intent(in) :: id + integer :: ierr, k + real(dp) :: max_v + type (star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_check_model = keep_going + + end function extras_check_model + + + integer function how_many_extra_history_columns(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_history_columns = 22 + end function how_many_extra_history_columns + + + subroutine data_for_extra_history_columns(id, n, names, vals, ierr) + integer, intent(in) :: id, n + character (len=maxlen_history_column_name) :: names(n) + real(dp) :: vals(n), v_esc + integer, intent(out) :: ierr + type (star_info), pointer :: s + integer :: k, k0 + integer, parameter :: h_relax_count = 1, & + h_log_R_098 = 2, & + h_log_R_095 = 3, & + h_log_R_090 = 4, & + h_helium_core_mass = 5, & + h_co_core_mass = 6, & + h_log_R_100 = 7, & + h_log_R_vesc = 8, & + h_log_R_vesc_098 = 9, & + h_log_R_vesc_095 = 10, & + h_log_R_vesc_090 = 11, & + h_M_below_vesc = 12, & + h_u_flag = 13, & + h_U_ejecta = 14, & + h_T_ejecta = 15, & + h_Omega_ejecta = 16, & + h_gamma_integral = 17, & + h_max_velocity = 18, & + h_min_velocity = 19, & + h_yr_since_coll = 20, & + h_total_J = 21, & + h_total_J_bel_vesc = 22 + + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + names(h_relax_count) = "relax_count" + vals(h_relax_count) = s% ixtra(ix_num_relaxations) + + names(h_log_R_098) = "log_R_098" ! Radius of 98% of the mass of the star + names(h_log_R_095) = "log_R_095" ! of 95% + names(h_log_R_090) = "log_R_090" ! and 90% + names(h_helium_core_mass) = "helium_core_mass" + names(h_co_core_mass) = "co_core_mass" + names(h_log_R_100) = "log_R_100" ! Radius of outermost layer + names(h_log_R_vesc) = "log_R_vesc" ! Radius of layers just below the escape velocity + names(h_log_R_vesc_098) = "log_R_vesc_098" ! Radius of 98% of the mass just below the esc. vel. + names(h_log_R_vesc_095) = "log_R_vesc_095" ! of 95% + names(h_log_R_vesc_090) = "log_R_vesc_090" ! and 90% + names(h_M_below_vesc) = "M_below_vesc" ! Mass of layers below vesc + names(h_u_flag) = "u_flag" ! 1 if hydro is on + names(h_U_ejecta) = "U_ejecta" ! internal energy of layers above vesc + names(h_T_ejecta) = "T_ejecta" ! kinetic energy of layers above vesc + names(h_Omega_ejecta) = "Omega_ejecta" ! gravitational energy of layers above vesc + names(h_gamma_integral) = "gamma_integral" ! integral of gamma1-4/3 in bound layers + names(h_max_velocity) = "max_velocity" + names(h_min_velocity) = "min_velocity" + names(h_yr_since_coll) = "yr_since_coll" + names(h_total_J) = "total_J" + names(h_total_J_bel_vesc) = "total_J_bel_vesc" + + vals(h_gamma_integral) = s% xtra(x_gamma_int_bound) + + if (s% u_flag) then + vals(h_max_velocity) = maxval(s% u(1:s% nz)) + else + vals(h_max_velocity) = 0 + end if + + if (s% u_flag) then + vals(h_min_velocity) = minval(s% u(1:s% nz)) + else + vals(h_min_velocity) = 0 + end if + + vals(h_yr_since_coll) = s% xtra(x_time_since_first_gamma_zero)/secyer + + vals(h_log_R_098) = -100 + vals(h_log_R_095) = -100 + vals(h_log_R_090) = -100 + vals(h_helium_core_mass) = -100 + vals(h_co_core_mass) = 0 + vals(h_log_R_100) = log10(s% r(1)/Rsun) + vals(h_log_R_vesc_098) = -100 + vals(h_log_R_vesc_095) = -100 + vals(h_log_R_vesc_090) = -100 + vals(h_U_ejecta) = 0d0 + vals(h_T_ejecta) = 0d0 + vals(h_Omega_ejecta) = 0d0 + do k=1, s% nz + if (s% q(k) < 0.98d0 .and. vals(h_log_R_098) == -100) then + vals(h_log_R_098) = log10(s% r(k)/Rsun) + else if (s% q(k) < 0.95d0 .and. vals(h_log_R_095) == -100) then + vals(h_log_R_095) = log10(s% r(k)/Rsun) + else if (s% q(k) < 0.9d0 .and. vals(h_log_R_090) == -100) then + vals(h_log_R_090) = log10(s% r(k)/Rsun) + end if + if (vals(h_helium_core_mass) < 0d0 .and. s% X(k) < 0.01d0) then + vals(h_helium_core_mass) = s% m(k)/Msun + end if + if (vals(h_co_core_mass) <= 0d0 .and. s% Y(k) < 0.01d0) then + vals(h_co_core_mass) = s% m(k)/Msun + end if + end do + + if (s% u_flag) then + do k = s% nz, 1, -1 + v_esc = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) + if (s% u(k) > v_esc) then + exit + end if + end do + if (k==0) then + k=1 + end if + else + k=1 + end if + + if (.not. s% rotation_flag) then + vals(h_total_J) = 0d0 + vals(h_total_J_bel_vesc) = 0d0 + else + vals(h_total_J) = dot_product(s% dm_bar(1:s% nz), s% j_rot(1:s% nz)) + vals(h_total_J_bel_vesc) = dot_product(s% dm_bar(k:s% nz), s% j_rot(k:s% nz)) + end if + + vals(h_log_R_vesc) = log10(s% r(k)/Rsun) + vals(h_M_below_vesc) = s% m(k)/Msun + ! get internal radii + do k0=k, s% nz + if (s% q(k0)/s% q(k) < 0.98d0 .and. vals(h_log_R_vesc_098) == -100) then + vals(h_log_R_vesc_098) = log10(s% r(k0)/Rsun) + else if (s% q(k0)/s% q(k) < 0.95d0 .and. vals(h_log_R_vesc_095) == -100) then + vals(h_log_R_vesc_095) = log10(s% r(k0)/Rsun) + else if (s% q(k0)/s% q(k) < 0.9d0 .and. vals(h_log_R_vesc_090) == -100) then + vals(h_log_R_vesc_090) = log10(s% r(k0)/Rsun) + end if + end do + ! get energies + if (k>1) then + do k0 = 1, k + vals(h_U_ejecta) = vals(h_U_ejecta) + s% dm(k0)*s% energy(k0) + if (s% u_flag) then + vals(h_T_ejecta) = vals(h_T_ejecta) + 0.5d0*s% dm(k0)*s% u(k0)*s% u(k0) + end if + vals(h_Omega_ejecta) = vals(h_Omega_ejecta) - s% dm_bar(k0)*s% cgrav(k0)*s% m(k0)/s% r(k0) + end do + end if + + if (s% u_flag) then + vals(h_u_flag) = 1d0 + else + vals(h_u_flag) = 0d0 + end if + + end subroutine data_for_extra_history_columns + + + integer function how_many_extra_profile_columns(id) + use star_def, only: star_info + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + if (s% u_flag) then + how_many_extra_profile_columns = 8 + else + how_many_extra_profile_columns = 1 + end if + end function how_many_extra_profile_columns + + + subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr) + use star_def, only: star_info, maxlen_profile_column_name + use const_def, only: dp + integer, intent(in) :: id, n, nz + character (len=maxlen_profile_column_name) :: names(n) + real(dp) :: vals(nz,n), alpha, J_inside + integer, intent(out) :: ierr + type (star_info), pointer :: s + integer :: k + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + J_inside = 0d0 + if (s% u_flag) then + names(1) = "vesc" + names(2) = "v_div_vesc" + names(3) = "specific_grav_e" + names(4) = "specific_kin_e" + names(5) = "specific_thermal_e" + names(6) = "total_specific_e" + names(7) = "mlt_vc" + names(8) = "spin_parameter" + + do k = s% nz, 1, -1 + vals(k,1) = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) + vals(k,2) = s% u(k)/vals(k,1) + vals(k,3) = -s% cgrav(k)*s% m(k)/s% r(k) + vals(k,4) = 0.5d0*s% u(k)*s% u(k) + vals(k,5) = two_thirds*avo*kerg*s% T(k)/(2*s% mu(k)*s% rho(k)) & + + crad*pow4(s% T(k))/s% rho(k) + vals(k,6) = vals(k,3) + vals(k,4) + vals(k,5) + vals(k,7) = s% mlt_vc(k) + if (s% rotation_flag) then + J_inside = J_inside + s% j_rot(k)*s% dm_bar(k) + vals(k,8) = J_inside*clight/(pow2(s% m(k))*standard_cgrav) + else + vals(k,8) = 0d0 + end if + end do + else + names(1) = "spin_parameter" + + if (s% rotation_flag) then + do k = s% nz, 1, -1 + J_inside = J_inside + s% j_rot(k)*s% dm_bar(k) + vals(k,1) = J_inside*clight/(pow2(s% m(k))*standard_cgrav) + end do + else + vals(:,1) = 0d0 + end if + end if + + + end subroutine data_for_extra_profile_columns + + + integer function extras_start_step(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + integer :: k, k0, k1, num_pts, species, model_number, num_trace_history_values + real(dp) :: v_esc, time, gamma1_integral, integral_norm, tdyn, & + max_center_cell_dq, avg_v_div_vesc, energy_removed_layers, dt_next, dt, & + max_years_for_timestep, omega_crit, & + denergy + real(dp) :: core_mass, rmax, alfa, log10_r, lburn_div_lsurf + logical :: just_did_relax + character (len=200) :: fname + include 'formats' + extras_start_step = terminate + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + !this is used to ensure we read the right inlist options + s% use_other_before_struct_burn_mix = .true. + + ! be sure power info is stored + call star_set_power_info(s) + + if (.not. in_inlist_pulses .and. .not. s% lxtra(lx_he_zams)) then + lburn_div_lsurf = abs(s% L_nuc_burn_total*Lsun/s% L(1)) + if (lburn_div_lsurf > 0d0) then + if((abs(log10(lburn_div_lsurf))) < 0.01 .and. & + (s% star_age > 1d3 .or. s% center_he4 < 0.98d0)) then + s% use_other_before_struct_burn_mix = .false. + call star_relax_uniform_omega(id, 1, s% job% new_omega_div_omega_crit,& + s% job% num_steps_to_relax_rotation, 1d0, ierr) + s% use_other_before_struct_burn_mix = .true. + if (ierr/=0) then + write(*,*) "error when relaxing omega at he ZAMS" + stop + end if + s% lxtra(lx_he_zams) = .true. + end if + end if + end if + + ! Ignore energy checks before first time hydro is turned on + ! otherwise need small steps during core helium burning and it + ! slows down the test_suite + if (.not. s% lxtra(lx_hydro_has_been_on)) then + s% cumulative_energy_error = 0d0 + s% cumulative_extra_heating = 0d0 + if(mod(s%model_number, s%terminal_interval) == 0) then + write(*,*) & + "Setting energy conservation error to zero until hydro is turned on for the first time" + end if + end if + + if (pgstar_flag) then + s% pg% pgstar_interval = pgstar_interval + s% pg% Grid2_file_interval = pgstar_file_interval + end if + just_did_relax = .false. + if (s% u_flag) then ! get point where v v_esc) then + exit + end if + end do + if (k>0 .and. k < s% nz) k = k+1 + else + k = 0 + end if + + ! s% xtra(x_gamma_int_bound) stores the gamma integral + ! during pulses integrate gamma only inside bound regions + ! otherwise compute it for the whole star + s% xtra(x_gamma_int_bound) = -1d0 + + ! can be adjusted below if nearing breakout + s% profile_interval = 100 + + if (s% u_flag .and. k > 0 .and. s% xtra(x_time_start_pulse) > 0d0) then + + ! check energy and average escape velocity of outer layers + avg_v_div_vesc = 0d0 + energy_removed_layers = 0d0 + do k0 = 1, k + avg_v_div_vesc = avg_v_div_vesc + s% dm(k0)*s% u(k0)/sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0))) + energy_removed_layers = energy_removed_layers + & + 0.5d0*s% dm(k0)*s% u(k0)*s% u(k0) - s% dm(k0)*s% cgrav(k0)*s% m(k0)/s% r(k0) & + +s% energy(k0)*s% dm(k0) + end do + ! adjust location of boundary to remove by considering also + ! material below the escape velocity that has a positive net + ! total specific energy. + if (energy_removed_layers > 0d0) then ! possible to eject material + if(mod(s%model_number, s%terminal_interval) == 0) then + write(*,*) "k, q, energy_removed_layers before adjustment is", k, s% q(k), energy_removed_layers + end if + do k0 = k+1, s% nz + denergy = & + 0.5d0*s% u(k0)*s% u(k0) - s% cgrav(k0)*s% m(k0)/s% r(k0) & + +s% energy(k0) + if (denergy < 0d0) then + k = k0-1 + exit + else + energy_removed_layers = energy_removed_layers+denergy*s% dm(k0) + avg_v_div_vesc = avg_v_div_vesc + & + s% dm(k0)*s% u(k0)/sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0))) + end if + end do + if(mod(s%model_number, s%terminal_interval) == 0) then + write(*,*) "k, q, energy_removed_layers after adjustment is", k, s% q(k), energy_removed_layers + end if + end if + + avg_v_div_vesc = avg_v_div_vesc/(s% m(1) - s% m(k)) + + ! compute gamma integral in what will remain of the star + integral_norm = 0d0 + gamma1_integral = 0d0 + do k0=k,s% nz + integral_norm = integral_norm + (s% Pgas(k0)+s% Prad(k0))*s% dm(k0)/s% rho(k0) + gamma1_integral = gamma1_integral + & + (s% gamma1(k0)-4.d0/3.d0)*(s% Pgas(k0)+s% Prad(k0))*s% dm(k0)/s% rho(k0) + end do + gamma1_integral = gamma1_integral/max(1d-99,integral_norm) + s% xtra(x_gamma_int_bound) = gamma1_integral + + do k1 = 1, s% nz + if (s% q(k1) < 0.9d0) then + !increase profile density near breakout, check that at q=0.9 velocities are + ! above 500 km/s, and that at the surface they're below that. + ! for the first pulse, increase profile density from the onset of instability + ! to breakout + if ((s% u(k1)>5d7 .and. s% u(1)<5d7) & + .or. (s% ixtra(ix_num_relaxations) == 0 .and. gamma1_integral < 0d0 .and. s% u(1)<5d7)) then + s% profile_interval = 10 + else + s% profile_interval = 100 + end if + exit + end if + end do + + ! To relax the star after a pulse, we check for a series of conditions that + ! must apply in layers below q=q_for_relax_check + do k0 = k, s% nz + if (s% q(k0) < q_for_relax_check*s% q(k)) then + exit + end if + end do + + if(mod(s%model_number, s%terminal_interval) == 0) then + write(*,*) 'Layers above q=', s% q(k), 'will be removed' + write(*,*) 'checking for conditions inside q=', q_for_relax_check, 'of material that will remain' + write(*,*) 'check time left', & + s% xtra(x_time_start_pulse) + s% xtra(x_dyn_time)*num_dyn_ts_for_relax - s% time + write(*,*) 'max vel inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)))/1e5 + write(*,*) 'max c/cs inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) + write(*,*) 'average v/vesc outside cutoff', avg_v_div_vesc + write(*,*) 'Kinetic plus potential energy outside cutoff', energy_removed_layers + write(*,*) 'mass inside cutoff', k, s% m(k)/Msun + write(*,*) 'relax counter', s% ixtra(ix_steps_met_relax_cond) + write(*,*) 'log max v [km/s]=', safe_log10(maxval(s% u(1:s% nz))/1e5) + end if + ! If enough time has passed after a pulse, then turn off Riemann hydro + ! also, wait for the inner q_for_relax_check of what would remain of the star to be below a threshold + ! in v/cs, and below a threshold in v + ! Verify also that the conditions to turn on Riemann hydro are not still satisfied + ! For details on all these options check inlist_project + + ! ignore if s% q(k) < 1d-3, in that case it's very likely a PISN + if (s% q(k) > 1d-3) then + if (s% lxtra(lx_hydro_on) .and. s% xtra(x_time_start_pulse) > 0d0 & + .and. s% time > s% xtra(x_time_start_pulse) + s% xtra(x_dyn_time)*num_dyn_ts_for_relax & + .and. maxval(abs(s% u(k0:s% nz)))/1d5 < max_v_for_relax & + .and. maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) < max_machn_for_relax & + .and. gamma1_integral > min_gamma_sub_43_for_hydro & + .and. log10(s% power_neutrinos) < max_Lneu_for_relax & + .and. log10(s% power_nuc_burn) < max_Lnuc_for_relax & + .and. energy_removed_layers > 0d0) then + s% ixtra(ix_steps_met_relax_cond) = s% ixtra(ix_steps_met_relax_cond) + 1 + else + s% ixtra(ix_steps_met_relax_cond) = 0 + end if + else + ! escape velocity reached within a tiny fraction of the + ! core. Before marking as PISN verify if any cell above + ! this is below the escape velocity + do k0 = k, 1, -1 + v_esc = sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0))) + if (s% u(k0) < v_esc) then + write(*,*) "Likely PISN", s% q(k0), s% u(k0)/v_esc + exit + end if + end do + ! If all cells were above the escape velocity, we + ! mark this as a PISN + if (k0 == 0) then + extras_start_step = terminate + write(fname, fmt="(a19)") 'LOGS/pisn_prof.data' + call star_write_profile_info(id, fname, ierr) + write(*,*) "Entire star is expanding above the escape velocity, PISN!" + return + end if + end if + if(s% ixtra(ix_steps_met_relax_cond) >= num_steps_before_relax .and. s% q(k) > 1d-3) then + write(*,*) "Relaxing model to lower mass!" + s% ixtra(ix_num_relaxations) = s% ixtra(ix_num_relaxations) + 1 + s% xtra(x_star_age_at_relax) = s% star_age + + !save a profile just before relaxation + write(fname, fmt="(a18,i0.3,a5)") 'LOGS/prerelax_prof', s% ixtra(ix_num_relaxations), '.data' + call star_write_profile_info(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + s% lxtra(lx_hydro_on) = .false. + s% lxtra(lx_have_saved_gamma_prof) = .false. + + s% ixtra(ix_steps_since_relax) = 0 + + write(*,*) "removing cells", k, s% m(k)/Msun + + max_center_cell_dq = s% max_center_cell_dq + s% max_center_cell_dq = s% dq(s% nz) + dt = s% dt + dt_next = s% dt_next + + call star_read_controls(id, 'inlist_hydro_off', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + call star_set_u_flag(id, .false., ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + call my_before_struct_burn_mix(s% id, s% dt, extras_start_step) + + ! to avoid showing pgstar stuff during relax + s% pg% pgstar_interval = 100000000 + s% pg% Grid2_file_interval = 100000000 + s% job% pgstar_flag = .false. + + max_years_for_timestep = s% max_years_for_timestep + s% max_years_for_timestep = 1d0 + + s% delta_lgL_nuc_limit = -1d0 + s% delta_lgL_nuc_hard_limit = -1d0 + s% use_other_before_struct_burn_mix = .false. + s% timestep_hold = 0 + + call star_relax_to_star_cut(s% id, k, .true., .true., .true., ierr) + if (ierr /= 0) then + write(*,*) "error when removing mass through star_relax_to_star_cut", ierr + stop + end if + + s% max_years_for_timestep = max_years_for_timestep + s% photo_interval = 100 + s% dt_next = min(1d2, dt_next) + s% dt = min(1d2, dt) + + s% ixtra(ix_steps_met_relax_cond) = 0 + + s% max_center_cell_dq = max_center_cell_dq + if (pgstar_flag) then + s% pg% pgstar_interval = pgstar_interval + s% pg% Grid2_file_interval = pgstar_file_interval + s% job% pgstar_flag = .true. + end if + just_did_relax = .true. + + !save a profile and a model just after relaxation + write(fname, fmt="(a19,i0.3,a5)") 'LOGS/postrelax_prof', s% ixtra(ix_num_relaxations), '.data' + call star_write_profile_info(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + write(fname, fmt="(a20,i0.3,a4)") 'LOGS/postrelax_model', s% ixtra(ix_num_relaxations), '.mod' + call star_write_model(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + s% use_other_before_struct_burn_mix = .true. + end if + end if + + ! turn on Riemman hydro if close to instability + integral_norm = 0d0 + gamma1_integral = 0d0 + do k=1,s% nz + integral_norm = integral_norm + (s% Pgas(k)+s% Prad(k))*s% dm(k)/s% rho(k) + gamma1_integral = gamma1_integral + & + (s% gamma1(k)-4.d0/3.d0)*(s% Pgas(k)+s% Prad(k))*s% dm(k)/s% rho(k) + end do + gamma1_integral = gamma1_integral/max(1d-99,integral_norm) + if (s% xtra(x_gamma_int_bound) == -1d0) then + ! if xtra(x_gamma_int_bound) is different from -1d0 it means it was computed for the bound material, + ! use that instead + s% xtra(x_gamma_int_bound) = gamma1_integral + end if + ! Save profile if gamma1_integral becomes negative + if (gamma1_integral < 0d0 .and. .not. s% lxtra(lx_have_saved_gamma_prof)) then + s% lxtra(lx_have_saved_gamma_prof) = .true. + write(fname, fmt="(a16,i0.3,a5)") 'LOGS/gamma1_prof', s% ixtra(ix_num_relaxations)+1, '.data' + call star_write_profile_info(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + end if + if(mod(s%model_number, s%terminal_interval) == 0) then + write(*,*) "check gamma integral", gamma1_integral + end if + + if (just_did_relax) then + extras_start_step = keep_going + return + end if + + if (s% ixtra(ix_steps_since_relax) > 10 .and. .not. just_did_relax .and. .not. s% lxtra(lx_hydro_on) & + .and. gamma1_integral < min_gamma_sub_43_for_hydro) then + write(*,*) "Turning on Riemann hydro!" + + call star_read_controls(id, 'inlist_hydro_on', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + call star_read_controls(id, 'inlist_after_first_pulse', ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + s% dt_next = min(1d2,s% dt_next) + s% dt = min(1d2,s% dt) + + !ensure implicit wind is not used + s% was_in_implicit_wind_limit = .false. + + write(fname, fmt="(a18,i0.3,a5)") 'LOGS/prehydro_prof', s% ixtra(ix_num_relaxations)+1, '.data' + call star_write_profile_info(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + write(fname, fmt="(a19,i0.3,a4)") 'LOGS/prehydro_model', s% ixtra(ix_num_relaxations)+1, '.mod' + call star_write_model(id, fname, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + call star_set_v_flag(id, .false., ierr)!TODO store v + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + + !s% job% pgstar_flag = .false. + call star_set_u_flag(id, .true., ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + s% lxtra(lx_hydro_on) = .true. + s% lxtra(lx_hydro_has_been_on) = .true. + s% xtra(x_time_start_pulse) = -1 + s% ixtra(ix_steps_since_hydro_on) = 0 + + !force initial velocities to zero to prevent issues in outer layers + s% xh(s% i_u,:) = 0d0 + s% xh_old(s% i_u,:) = 0d0 + s% generations = 1 + end if + + ! check time to relax model + ! the star is considered to be in a pulse if the velocity is high enough + ! only check velocity in q_for_relax_check of the mass of the star, this avoids + ! considering as a pulse cases where the surface goes a bit nuts + do k0 = 1, s% nz + if (s% q(k0) < q_for_relax_check) then + exit + end if + end do + if (s% lxtra(lx_hydro_on) .and. s% xtra(x_time_start_pulse) < 0d0) then + if ((maxval(abs(s% xh(s% i_u, k0:s% nz)))/1d5 > max_v_for_pulse .or. gamma1_integral < 1d-3)) then + core_mass = s% m(1) + !if(s% o_core_mass > 0d0) then + ! core_mass = s% o_core_mass*Msun + !else if(s% c_core_mass > 0d0) then + if(s% co_core_mass > 0d0) then + core_mass = s% co_core_mass*Msun + else if (s% he_core_mass > 0d0) then + core_mass = s% he_core_mass*Msun + end if + do k=1, s% nz + if (s% m(k)/core_mass < q_for_dyn_ts) then + exit + end if + end do + tdyn = 1d0/sqrt(s% m(k)/(4d0/3d0*pi*pow3(s% r(k)))*standard_cgrav) + s% xtra(x_time_start_pulse) = s% time + s% xtra(x_dyn_time) = tdyn + write(*,*) "reached high velocities", maxval(abs(s% xh(s% i_u, k0:s% nz)))/1e5, & + "stopping in at least", num_dyn_ts_for_relax*tdyn, "seconds" + end if + end if + + !After a relax, wait for ten days before turning on v_flag + !this avoids the surface post-relax from going crazy + if((logT_for_v_flag < log10(s% T(s% nz)) .or. logLneu_for_v_flag < safe_log10(s% power_neutrinos)) & + .and. .not. s% u_flag .and. .not. s% v_flag) then + write(*,*) "log10 central T has lowered below logT_for_v_flag, turn on v_flag" + call star_set_v_flag(id, .true., ierr) + s% dt_next = min(1d2, s% dt_next) + s% dt = min(1d2, s% dt) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + end if + + !Always call this at the end to ensure we are using the correct + !inlists + call my_before_struct_burn_mix(s% id, s% dt, extras_start_step) + + extras_start_step = keep_going + end function extras_start_step + + subroutine my_before_struct_burn_mix(id, dt, res) + use const_def, only: dp + use star_def + integer, intent(in) :: id + real(dp), intent(in) :: dt + integer, intent(out) :: res ! keep_going, redo, retry, terminate + real(dp) :: power_photo, v_esc + integer :: ierr, k + type (star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + !read proper options when hydro is on/off + !do this to ensure proper behaviour of retries + if(s% u_flag) then + call star_read_controls(id, 'inlist_hydro_on', ierr) + if (s% xtra(x_time_start_pulse) > 0d0) then + s% v_drag = s% u(1) + s% v_drag_factor = 1d0 + write(*,*) 'using drag' + if (.not. s% lxtra(lx_using_bb_bcs)) then + s% max_timestep = max_dt_during_pulse + else + s% max_timestep = 1d99 + end if + s% max_q_for_convection_with_hydro_on = 1d0!0.999d0 + do k = s% nz, 1, -1 + v_esc = sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k))) + if (s% u(k) > 4*v_esc) then ! 4 in farmer inlist + exit + end if + end do + if (k > 1) then + s% max_q_for_convection_with_hydro_on = s% q(k) + end if + else + s% max_timestep = max_dt_before_pulse + !!this prevents surface from going mad right when hydro is turned on + !right after a relax + s% v_drag = 0d0 + s% v_drag_factor = 1d0 + !s% max_timestep = 1d99 + s% max_q_for_convection_with_hydro_on = 0.999d0 + end if + + do k=1, s% nz + s% xh(s% i_u,k) = min(s% xh(s% i_u,k), 1d5*vsurf_for_fixed_bc) + s% u(k) = s% xh(s% i_u,k) + end do + + ! use fixed_vsurf if surface v remains too high + if (s% xh(s% i_u,1) >= 1d5*vsurf_for_fixed_bc) then + s% use_fixed_vsurf_outer_BC = .true. + s% use_momentum_outer_BC = .false. + s% fixed_vsurf = vsurf_for_fixed_bc*1d5 + else + if (s% lxtra(lx_using_bb_bcs)) then + s% use_fixed_vsurf_outer_BC = .true. + s% use_momentum_outer_BC = .false. + s% fixed_vsurf = s% xtra(x_vsurf_at_remove_surface) + else + s% use_fixed_vsurf_outer_BC = .false. + s% use_momentum_outer_BC = .true. + end if + end if + + if (maxval(abs(s% xh(s% i_u,:s% nz)))<2d7) then + s% dt_div_min_dr_div_cs_limit = 2d1 + else + s% dt_div_min_dr_div_cs_limit = dt_div_min_dr_div_cs_limit + end if + + else + s% max_timestep = 1d99 + call star_read_controls(id, 'inlist_hydro_off', ierr) + end if + !ensure we are using the proper BCs and options every time before struct_burn_mix + if (s% lxtra(lx_using_bb_bcs)) then + !s% use_atm_PT_at_center_of_surface_cell = .true. + s% use_momentum_outer_BC = .true. ! offset_P_to_center_cell = .false. + !s% offset_P_to_center_cell = .false. + s% atm_option = "fixed_Tsurf" + s% atm_fixed_Tsurf = s% xtra(x_Tsurf_for_atm) + s% tau_factor = s% xtra(x_tau_factor) + s% force_tau_factor = s% xtra(x_tau_factor) + s% delta_lgL_limit = -1d0 + s% delta_lgL_limit_L_min = 1d99 + s% delta_lgL_limit_L_min = 1d99 + else + !s% use_atm_PT_at_center_of_surface_cell = .false. + s% use_momentum_outer_BC = .false. ! offset_P_to_center_cell = .true. + !s% offset_P_to_center_cell = .true. + s% atm_option = 'T_tau' + s% atm_T_tau_relation = 'Eddington' + s% atm_T_tau_opacity = 'fixed' + s% tau_factor = 1d0 + !s% Pextra_factor = 1d0 + s% force_tau_factor = 1d0 + s% delta_lgL_limit = 0.25d0 + s% delta_lgL_limit_L_min = 1d99!-100 + s% delta_lgL_limit_L_min = 1d99!-100 + end if + + + + + + + + !ignore L_nuc limit if L_phot is too high or if we just did a relax + !(ixtra(ix_steps_since_relax) is set to zero right after a relax) + + !when L_phot exceeds max_Lphoto_for_lgLnuc_limit, the timestep limit is applied + !to L_phot instead + + ! be sure power info is stored + call star_set_power_info(s) + + power_photo = dot_product(s% dm(1:s% nz), s% eps_nuc_categories(iphoto,1:s% nz))/Lsun + if (safe_log10(abs(power_photo)) > max_Lphoto_for_lgLnuc_limit2) then + s% delta_lgL_nuc_limit = -1d0 + s% delta_lgL_nuc_hard_limit = -1d0 + s% delta_lgL_power_photo_limit = -1d0 + s% delta_lgL_power_photo_hard_limit = -1d0 + else + if (s% ixtra(ix_steps_since_relax) == 0 & + .or. safe_log10(abs(power_photo)) > max_Lphoto_for_lgLnuc_limit) then + s% delta_lgL_nuc_limit = -1d0 + s% delta_lgL_nuc_hard_limit = -1d0 + else + s% delta_lgL_nuc_limit = delta_lgLnuc_limit + s% delta_lgL_nuc_hard_limit = 2d0*delta_lgLnuc_limit + end if + if (safe_log10(abs(power_photo)) > max_Lphoto_for_lgLnuc_limit) then + s% delta_lgL_power_photo_limit = delta_lgLnuc_limit + s% delta_lgL_power_photo_hard_limit = 2d0*delta_lgLnuc_limit + else + s% delta_lgL_power_photo_limit = -1d0 + s% delta_lgL_power_photo_hard_limit = -1d0 + end if + end if + if (s% use_fixed_vsurf_outer_BC) then + write(*,*) "Using fixed_vsurf", s% fixed_vsurf/1e5 + end if + !ignore winds if neutrino luminosity is too high or for a few steps after + !a relax + if(s% ixtra(ix_steps_since_relax) < 50 & + .or. safe_log10(s% power_neutrinos) > max_Lneu_for_mass_loss & + .or. s% u_flag) then + s% use_other_wind = .false. + s% was_in_implicit_wind_limit = .false. + else + s% use_other_wind = .true. + end if + + if (maxval(s% T(1:s% nz)) > 9.8) then + s% delta_lgRho_cntr_hard_limit = -1d0 + else + s% delta_lgRho_cntr_hard_limit = delta_lgRho_cntr_hard_limit + end if + + ! reading inlists can turn this flag off for some reason + s% use_other_before_struct_burn_mix = .true. + + res = keep_going + end subroutine my_before_struct_burn_mix + + subroutine null_binary_controls(id, binary_id, ierr) + integer, intent(in) :: id, binary_id + integer, intent(out) :: ierr + ierr = 0 + end subroutine null_binary_controls + + ! returns either keep_going or terminate. + integer function extras_finish_step(id) + use run_star_support + integer, intent(in) :: id + integer :: ierr,k + real(dp) :: max_vel_inside + type (star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + extras_finish_step = keep_going + + !count time since first collapse + if (s% lxtra(lx_have_reached_gamma_limit)) then + s% xtra(x_time_since_first_gamma_zero) = & + s% xtra(x_time_since_first_gamma_zero) + s% dt + end if + + + ! if (s% conv_vel_flag) then + ! do k=1, s% nz + ! if (s% conv_vel(k) < 1d-5 .and. s% mlt_vc(k) <= 0d0) then + ! s% xh(s% i_ln_cvpv0, k) = 0d0 + ! s% conv_vel(k) = 0d0 + ! end if + ! end do + ! end if + + if (remove_extended_layers .and. s% u_flag & + .and. s% log_surface_radius > 4.0d0) then + do k = 1, s% nz + if (log10(s% r(k)/Rsun) < 4.0d0) then + exit + end if + end do + s% atm_option = "fixed_Tsurf" + s% atm_fixed_Tsurf = s% T(k) + s% xtra(x_Tsurf_for_atm) = s% T(k) + s% xtra(x_vsurf_at_remove_surface) = max(s% u(k),2*sqrt(2d0*standard_cgrav*s% m(k)/(pow(10d0,4.0d0)*Rsun))) + s% lxtra(lx_using_bb_bcs) = .true. + + write(*,*) "Removing surface layers", k, s% m(k)/Msun + call star_remove_surface_at_cell_k(s% id, k, ierr) + if (dbg) write(*,*) "check ierr", ierr + if (ierr /= 0) return + s% xtra(x_tau_factor) = s% tau_factor + end if + + s% ixtra(ix_steps_since_relax) = s% ixtra(ix_steps_since_relax) + 1 + s% ixtra(ix_steps_since_hydro_on) = s% ixtra(ix_steps_since_hydro_on) + 1 + + if (s% ixtra(ix_num_relaxations) == 1 .and. stop_100d_after_pulse & + .and. s% star_age - s% xtra(x_star_age_at_relax) > 100d0/dayyer) then + !for the test_suite, terminate at the onset of the second pulse + extras_finish_step = terminate + s% termination_code = t_xtra1 + termination_code_str(t_xtra1) = "Successful test: evolved 100 days past first relax" + return + end if + + if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step + + end function extras_finish_step + + end module run_star_extras + diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/standard_he_dep.mod b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/standard_he_dep.mod new file mode 100644 index 000000000..0058b23ab --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/standard_he_dep.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d44651e90a73db4a753274968e3e2649ed1cecb47fbdd72a0d915a3514ef6329 +size 1243628 diff --git a/star/dev_cases_test_TDC/dev_TDC_through_ppisn/zams.mod b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/zams.mod new file mode 100644 index 000000000..2d558da66 --- /dev/null +++ b/star/dev_cases_test_TDC/dev_TDC_through_ppisn/zams.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0c89bfcb53e4c027cfb33bec543defcc88d508e71c3792cf2b1b4add6a9278d7 +size 381596 diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 2673028f9..543b16c8d 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -352,6 +352,7 @@ module ctrls_io steps_before_use_TDC, use_P_d_1_div_rho_form_of_work_when_time_centering_velocity, compare_TDC_to_MLT, & velocity_logT_lower_bound, max_dt_yrs_for_velocity_logT_lower_bound, velocity_q_upper_bound, & drag_coefficient, min_q_for_drag, & + v_drag_factor, v_drag, q_for_v_drag_full_off, q_for_v_drag_full_on, & retry_for_v_above_clight, & ! hydro solver @@ -1863,6 +1864,11 @@ subroutine store_controls(s, ierr) s% drag_coefficient = drag_coefficient s% min_q_for_drag = min_q_for_drag + s% v_drag_factor = v_drag_factor + s% v_drag = v_drag + s% q_for_v_drag_full_off = q_for_v_drag_full_off + s% q_for_v_drag_full_on = q_for_v_drag_full_on + s% velocity_logT_lower_bound = velocity_logT_lower_bound s% max_dt_yrs_for_velocity_logT_lower_bound = max_dt_yrs_for_velocity_logT_lower_bound @@ -3534,6 +3540,10 @@ subroutine set_controls_for_writing(s, ierr) drag_coefficient = s% drag_coefficient min_q_for_drag = s% min_q_for_drag + v_drag_factor = s% v_drag_factor + v_drag = s% v_drag + q_for_v_drag_full_off = s% q_for_v_drag_full_off + q_for_v_drag_full_on = s% q_for_v_drag_full_on velocity_logT_lower_bound = s% velocity_logT_lower_bound max_dt_yrs_for_velocity_logT_lower_bound = s% max_dt_yrs_for_velocity_logT_lower_bound diff --git a/star/private/hydro_riemann.f90 b/star/private/hydro_riemann.f90 index 22f3a5a1d..502f6fde1 100644 --- a/star/private/hydro_riemann.f90 +++ b/star/private/hydro_riemann.f90 @@ -96,6 +96,7 @@ subroutine do1_dudt_eqn( & type(accurate_auto_diff_real_star_order1) :: sum_ad real(dp) :: dt, dm, ie_plus_ke, scal, residual logical :: dbg, do_diffusion, test_partials + real(dp) :: v_drag, drag_factor, drag_fraction include 'formats' dbg = .false. @@ -135,6 +136,27 @@ subroutine do1_dudt_eqn( & dudt_expected_ad = sum_ad dudt_expected_ad = dudt_expected_ad/dm + ! implement drag + drag_factor = s% v_drag_factor + v_drag = s% v_drag + if (s% q(k) < s% q_for_v_drag_full_off) then + drag_fraction = 0d0 + else if (s% q(k) > s% q_for_v_drag_full_on) then + drag_fraction = 1d0 + else + drag_fraction = (s% q(k) - s% q_for_v_drag_full_off)& + /(s% q_for_v_drag_full_on - s% q_for_v_drag_full_off) + end if + drag_factor = drag_factor*drag_fraction + + if (drag_factor > 0d0) then + if (s% u(k) > v_drag) then + dudt_expected_ad = dudt_expected_ad - drag_factor*pow2(s% u(k) - v_drag)/s% r(k) + else if (s% u(k) < -v_drag) then + dudt_expected_ad = dudt_expected_ad + drag_factor*pow2(s% u(k) + v_drag)/s% r(k) + end if + end if + ! make residual units be relative difference in energy ie_plus_ke = s% energy_start(k) + 0.5d0*s% u_start(k)*s% u_start(k) scal = dt*max(abs(s% u_start(k)),s% csound_start(k))/ie_plus_ke diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index 971b51eeb..f02dada81 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -832,7 +832,11 @@ real(dp) :: & drag_coefficient, & - min_q_for_drag + min_q_for_drag, & + v_drag_factor, & + v_drag, & + q_for_v_drag_full_off, & + q_for_v_drag_full_on real(dp) :: & RTI_A, RTI_B, RTI_C, RTI_D, & From 7190be18b2d68dd3aa51d8fd9bc66ceaacfb14fa Mon Sep 17 00:00:00 2001 From: Eb Date: Wed, 31 Jan 2024 11:59:58 -0700 Subject: [PATCH 3/3] fix rti bugs and add documentation for drag --- star/defaults/controls.defaults | 31 +++++++++++++++++++++++++++---- star/private/mix_info.f90 | 2 +- star/private/remove_shells.f90 | 12 ++++++++++-- 3 files changed, 38 insertions(+), 7 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index ed0d2b4c2..7c4237d62 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2983,7 +2983,7 @@ conv_premix_dump_snapshots = .false. - ! Rayleigh Taylor Instability + ! Rayleigh-Taylor Instability ! ____________________________ ! derived from Paul Duffell's code RT1D. @@ -8196,8 +8196,11 @@ use_dPrad_dm_form_of_T_gradient_eqn = .false. use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false. - - + + + ! Hydrodynamic drag + ! ================= + ! drag_coefficient ! ~~~~~~~~~~~~~~~~ ! min_q_for_drag @@ -8219,12 +8222,32 @@ drag_coefficient = 0d0 min_q_for_drag = 0d0 + + ! v_drag_factor + ! ~~~~~~~~~~~~~ + ! v_drag + ! ~~~~~~ + ! q_for_v_drag_full_off + ! ~~~~~~~~~~~~~~~~~~~~~ + ! q_for_v_drag_full_on + ! ~~~~~~~~~~~~~~~~~~~~ - ! for hydro comparison tests (e.g., Sedov) + ! Implemented only for u_flag right now. Adds a pseudo drag term of the form + ! -v_drag_factor*(v-v_drag)^2/r, can be used damp velocities in outer layers + ! of a star. Effect is full on for q>q_for_v_drag_full_on and full off for + ! q < q_for_v_drag_full_off. + + ! :: + + v_drag_factor = 0d0 + v_drag = 0d0 + q_for_v_drag_full_off = 0.95d0 + q_for_v_drag_full_on = 0.96d0 ! Rayleigh-Taylor Instability ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! for hydro comparison tests (e.g., Sedov) ! RTI_A ! ~~~~~ ! RTI_B diff --git a/star/private/mix_info.f90 b/star/private/mix_info.f90 index b27f174fb..2d4247de0 100644 --- a/star/private/mix_info.f90 +++ b/star/private/mix_info.f90 @@ -2016,7 +2016,7 @@ subroutine set_RTI_mixing_info(s, ierr) r = s% r(k) s% eta_RTI(k) = C*alpha_face*cs*r - if (is_bad(s% eta_RTI(k))) then + if (is_bad(s% eta_RTI(k)) .and. s% q(k) <= s% alpha_RTI_src_max_q) then ierr = -1 return if (s% stop_for_bad_nums) then diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index 746595a66..69fcc3a40 100644 --- a/star/private/remove_shells.f90 +++ b/star/private/remove_shells.f90 @@ -1138,7 +1138,7 @@ subroutine do_relax_to_star_cut( & use interp_1d_def, only: pm_work_size use interp_1d_lib, only: interp_pm, interp_values, interp_value use adjust_xyz, only: change_net - use set_flags, only: set_v_flag, set_u_flag, set_rotation_flag + use set_flags, only: set_v_flag, set_u_flag, set_RTI_flag, set_rotation_flag use rotation_mix_info, only: set_rotation_mixing_info use hydro_rotation, only: set_i_rot, set_rotation_info use relax, only: do_relax_composition, do_relax_angular_momentum, do_relax_entropy @@ -1150,7 +1150,7 @@ subroutine do_relax_to_star_cut( & ! determines if we turn off non_nuc_neu and eps_nuc for entropy relax integer, intent(out) :: ierr - logical :: v_flag, u_flag, rotation_flag + logical :: v_flag, u_flag, RTI_flag, rotation_flag type (star_info), pointer :: s character (len=net_name_len) :: net_name integer :: model_number, num_trace_history_values, photo_interval @@ -1240,6 +1240,7 @@ subroutine do_relax_to_star_cut( & if (dbg) write(*,*) "set_v_flag ierr", ierr v_flag = .true. end if + u_flag = .false. if (s% u_flag) then call set_u_flag(id, .false., ierr) @@ -1247,6 +1248,13 @@ subroutine do_relax_to_star_cut( & u_flag = .true. end if + RTI_flag = .false. + if (s% RTI_flag) then + call set_RTI_flag(id, .false., ierr) + if (dbg) write(*,*) "set_RTI_flag ierr", ierr + RTI_flag = .true. + end if + if (s% rotation_flag) then call set_rotation_flag(id, .false., ierr) if (dbg) write(*,*) "set_rotation_flag ierr", ierr