diff --git a/net/private/net_approx21.f90 b/net/private/net_approx21.f90 index ef2fbb239..621b5559e 100644 --- a/net/private/net_approx21.f90 +++ b/net/private/net_approx21.f90 @@ -1683,7 +1683,17 @@ subroutine approx21_eps_info( & xx = xx + a1*a2 end do eps_total_q = -m3(avo,clight,clight) * xx + ! XXX error here??? eps_total is really negative? eps_total = eps_total_q +! if (n% zone >= 1816 .and. n% zone <= 1819) then +! xx = 0.0_qp +! do i=1,species(plus_co56) +! a1 = dydt(i) +! a2 = mion(i) +! xx = xx + a1*a2 +! write(*,*) n% zone, eps_total, i, a1, a2, xx +! end do +! end if fe56ec_fake_factor = eval_fe56ec_fake_factor( & n% g% fe56ec_fake_factor, n% g% min_T_for_fe56ec_fake_factor, n% temp) diff --git a/net/private/net_eval.f90 b/net/private/net_eval.f90 index d0c020a66..cd30fe010 100644 --- a/net/private/net_eval.f90 +++ b/net/private/net_eval.f90 @@ -384,6 +384,8 @@ subroutine eval_net_approx21_procs(n,just_dxdt, ierr) if (ierr /= 0) return n% eps_nuc = n% eps_total - n% eps_neu_total + ! XXX test: + !n% eps_nuc = n% eps_total do i=1, num_isos n% dxdt(i) = chem_isos% Z_plus_N(g% chem_id(i)) * n% dydt1(i) diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/.gitattributes b/star/dev_cases_test_solver/test_20M_near_cc_approx21/.gitattributes new file mode 100644 index 000000000..b896ae1ca --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/.gitattributes @@ -0,0 +1,2 @@ +plots/* filter=lfs diff=lfs merge=lfs -text +photos/* filter=lfs diff=lfs merge=lfs -text diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/README.rst b/star/dev_cases_test_solver/test_20M_near_cc_approx21/README.rst new file mode 100644 index 000000000..f06a59348 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/README.rst @@ -0,0 +1,45 @@ +.. _20M_pre_ms_to_core_collapse: + +*************************** +20M_pre_ms_to_core_collapse +*************************** + +This test suite evolves a solar metalicity 20 |MSun| model from the pre-ms to core collapse. +For bit for bit convergence, we recomended to run by using the ./run_all script instead of restarting from models, +see https://github.com/MESAHub/mesa/issues/610. + +This test_suite has been tested up to 80 solar masses, up to solar metallicity, with mass loss, and produces reasonable HR-tracks. +Note that for higher masses at solar metallicity, some combination of Pextra_factor, mass-loss, and/or superadiabatic convection reduction (e.g. mlt++) +might be necessary to stabilize the surface and avoid numerical issues. See the 80Msun_zams_to_cc test_suite as an example. + +For production science we recommend adopting tighter mesh and timestep controls, such as those suggested in the comments of inlist_common. + +Physical checks +=============== + +None + +Inlists +======= + +This test case has seven parts. + +* Part 1 (``inlist_make_late_pre_zams``) creates a 20 |Msun|, Z=1.42*10^-2 metallicity, pre-main sequence model and evolves it for 100 years. + +* Part 2 (``inlist_to_zams``) evolves the model to the zero age main sequence. + +* Part 3 (``inlist_to_end_core_he_burn``) takes the model to core helium depletion. + +* Part 4 (``inlist_remove_envelope``) removes the remianing hydrogen envelope. (optional) + +* Part 5 (``inlist_to_end_core_c_burn``) takes the model to core carbon depletion. + +* Part 6 (``inlist_to_lgTmax``) evolves the model until the core temperature reaches log T =9.60 (approximately silicon-shell burning) + +* Part 7 (``inlist_to_cc``) evolves until core collapse. + + + + +Last-Updated: 18Dec2023 by EbF + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/ck b/star/dev_cases_test_solver/test_20M_near_cc_approx21/ck new file mode 100755 index 000000000..ac08f15c6 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/clean b/star/dev_cases_test_solver/test_20M_near_cc_approx21/clean new file mode 100755 index 000000000..95545a5c1 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/clean @@ -0,0 +1,4 @@ +#!/bin/bash + +cd make +make clean diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/history_columns.list b/star/dev_cases_test_solver/test_20M_near_cc_approx21/history_columns.list new file mode 100644 index 000000000..7adf357a6 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/history_columns.list @@ -0,0 +1,1071 @@ +! 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/astero_Teff_sun)^0.5) + !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 80 + mixing_regions 40 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_common b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_common new file mode 100644 index 000000000..c99582cf9 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_common @@ -0,0 +1,403 @@ +! the contents of this are set to match the needs of pre_ms_to_cc_12_20 +! so local inlist_common files show differences from what is used for that case. + +&star_job + show_log_description_at_start = .false. + + change_net = .true. + new_net_name = 'approx21_cr60_plus_co56.net' + dump_missing_metals_into_heaviest = .false. + + ! Example ccsn network for properly capturing the energy generation rate during all stages + !change_net_net = .true. + !new_net_name = 'mesa_206.net' + !dump_missing_metals_into_heaviest = .false. + + + ! if using rotation + change_D_omega_flag = .true. + new_D_omega_flag = .true. + + ! turn on hydrodynamics + change_v_flag = .true. + new_v_flag = .true. + + ! Can be adopted to avoid issues with surface for high Z or during advanced burning + !relax_to_this_tau_factor = 1.5d6 + !dlogtau_factor = 0.1d0 + !relax_initial_tau_factor = .true. + + ! If using a big network, comment out and used a local rate_tables directory + 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' + + show_retry_counts_when_terminate = .true. + show_timestep_limit_counts_when_terminate = .true. + + ! pgstar is recommended for diagnosing issues and understanding the evolution + pgstar_flag = .true. + save_pgstar_files_when_terminate = .true. + +/ ! end of star_job namelist + +&eos + use_Skye = .true. + use_PC = .false. +/ ! end of eos namelist + +&kap + ! OPAL asplund 2009 opacities, Zbase set in inlist_mass_Z_wind_rotation + kap_file_prefix = 'a09' ! 'gs98' 'a09' 'OP_a09' 'OP_gs98' + kap_CO_prefix = 'a09_co' ! 'gs98_co' 'a09_co' + kap_lowT_prefix = 'lowT_fa05_a09p' + use_Type2_opacities = .true. + +/ ! end of kap namelist + +&controls + +! wind + ! Dutch scaling factor set in inlist_mass_Z_wind_rotation + cool_wind_full_on_T = 0.8d4 + hot_wind_full_on_T = 1.2d4 + cool_wind_RGB_scheme = 'Dutch' + cool_wind_AGB_scheme = 'Dutch' + hot_wind_scheme = 'Dutch' + Dutch_wind_lowT_scheme = 'de Jager' + + max_T_center_for_any_mass_loss = 1.1d9 + +! atmosphere + Pextra_factor = 2 ! easy to lower with mass loss!ideally 1.5 is most physical + ! if you run into issues, you can increase Pextra + ! extra pressure helps stabilize the atmosphere during core He burning + + + + atm_option = 'T_tau' + atm_T_tau_relation = 'Eddington' + atm_T_tau_opacity = 'fixed' ! next best is 'iterated' + +! rotation + am_nu_visc_factor = 0 + am_D_mix_factor = 0.0333333333333333d00 + D_DSI_factor = 0 + D_SH_factor = 1 + D_SSI_factor = 1 + D_ES_factor = 1 + D_GSF_factor = 1 + D_ST_factor = 1 + +! mlt + mixing_length_alpha = 1.5 + MLT_option = 'TDC' + + ! If you would like to be bold, try radiative damping, although it + ! is a departure from the mlt limit of 'TDC' + alpha_TDC_DAMPR = 0d0 ! 0d0 is default for mlt limit + + + use_other_alpha_mlt = .true. ! implemented in run_star_extras + x_ctrl(21) = 3.0 ! alpha_H + x_ctrl(22) = 1.5 ! alpha_other + x_ctrl(23) = 0.5 ! ! use alpha_H if cell X >= H_limit; else use alpha_other + x_ctrl(24) = 9d0 ! ! use other_alpha_mlt only if star_mass >= this limit. + + + use_Ledoux_criterion = .true. + semiconvection_option = 'Langer_85 mixing; gradT = gradr' + thermohaline_option = 'Kippenhahn' + + alpha_semiconvection = 1d-2 + thermohaline_coeff = 0 + + mlt_make_surface_no_mixing = .true. + + + +! superadiabatic convection routines, it's a choice: + + ! superadiabatic reduction, implicit, new + use_superad_reduction = .false. + superad_reduction_Gamma_limit = 0.3d0 ! default is 0.5d0 + superad_reduction_Gamma_limit_scale = 5d0 + superad_reduction_Gamma_inv_scale = 5d0 + superad_reduction_diff_grads_limit = 1d-2 ! default is 1d-3 + superad_reduction_limit = -1d0 + + + ! MLT ++, explicit, well tested + okay_to_reduce_gradT_excess = .false. + gradT_excess_f1 = 1d-4 + gradT_excess_f2 = 1d-2 + gradT_excess_lambda1 = -1d0 ! full on + +! mixing + D_omega_mixing_rate = 1d0 + D_omega_mixing_across_convection_boundary = .false. + + + + ! we use step overshooting in H core + overshoot_scheme(1) = 'step' + overshoot_zone_type(1) = 'burn_H' + overshoot_zone_loc(1) = 'core' + overshoot_bdy_loc(1) = 'top' + overshoot_f(1) = 0.345 ! for M>10 + !overshoot_f(1) = 0.21 ! For M<10 + overshoot_f0(1) = 0.01 + + ! exponential in the H core + overshoot_scheme(2) = 'exponential' + overshoot_zone_type(2) = 'burn_He' + overshoot_zone_loc(2) = 'core' + overshoot_bdy_loc(2) = 'top' + overshoot_f(2) = 0.01 + overshoot_f0(2) = 0.005 + + ! we don't want to deal with He/CO core mergers + ! and there is reason to believe there is little + ! inward overshooting in the shell across compositions boundaries + overshoot_scheme(3) = 'none' + overshoot_zone_type(3) = 'burn_He' + overshoot_zone_loc(3) = 'shell' + overshoot_bdy_loc(3) = 'bottom' + + ! a small amount of overshooting on top of any other convective core + ! avoid spurious numerical behavior + ! perfect amount for degenerate flames + overshoot_scheme(4) = 'exponential' + overshoot_zone_type(4) = 'any' + overshoot_zone_loc(4) = 'any' + overshoot_bdy_loc(4) = 'any' + overshoot_f(4) = 0.005d0 + overshoot_f0(4) = 0.001d0 + + +! timesteps + time_delta_coeff = 1.0 + varcontrol_target = 1d-3 + + min_timestep_factor = 0.8d0 + max_timestep_factor = 1.05d0 + timestep_factor_for_retries = 0.75 + + limit_for_rel_error_in_energy_conservation = 1d-7 + hard_limit_for_rel_error_in_energy_conservation = 1d-6 + + never_skip_hard_limits = .true. + min_xa_hard_limit = -1d-5 + min_xa_hard_limit_for_highT = -3d-5 + + delta_lgTeff_limit = 0.01 + delta_lgL_limit = 0.1 + delta_lgL_He_limit = 0.1 + + + ! Recommend decreasing all three Rho, T, Tmax + ! to 1d-3 or lower in production runs + delta_lgRho_cntr_limit = 0.03 + delta_lgRho_cntr_hard_limit = 0.1 + delta_lgRho_limit = 0.1 + + delta_lgT_cntr_limit_only_after_near_zams = .true. + delta_lgT_cntr_limit = 0.002 + delta_lgT_cntr_hard_limit = 0.1 + + delta_lgT_max_limit_only_after_near_zams = .true. + delta_lgT_max_limit = 0.002 + delta_lgT_max_hard_limit = 0.1 + + dX_div_X_limit(2) = -1 ! for he4 + + ! On the changes in total abundance of each isotope + ! one of the most useful timestep controls, period + dX_nuc_drop_limit = 2d-2 ! Recommend decreasing to 1d-3 or lower in a production run + dX_nuc_drop_limit_at_high_T = 2d-2 ! default = -1 = same as dX_nuc_drop_limit + dX_nuc_drop_min_X_limit = 1d-3 ! try decreasing to 1d-4 or 1d-5 in a production run + dX_nuc_drop_max_A_limit = 70 ! try increasing beyond 60 in a big network run + dX_nuc_drop_hard_limit = 1d99 + + + ! delta_XHe_cntr_limit = 0.01d0 + ! delta_XHe_cntr_hard_limit = 0.1d0 + ! delta_lg_XH_cntr_limit = -1 + ! delta_lg_XHe_cntr_limit = -1 + ! delta_lg_XC_cntr_limit = -1 + ! delta_lg_XNe_cntr_limit = -1 + ! delta_lg_XO_cntr_limit = -1 + ! delta_lg_XSi_cntr_limit = -1 + ! delta_XC_cntr_limit = -1 + ! delta_XC_cntr_hard_limit = -1 + ! delta_XNe_cntr_limit = -1 + ! delta_XNe_cntr_hard_limit = -1 + ! delta_XO_cntr_limit = -1 + ! delta_XO_cntr_hard_limit = -1 + ! delta_XSi_cntr_limit = -1 + ! delta_XSi_cntr_hard_limit = -1 + + + + delta_lg_XH_cntr_limit = 0.01d0 + delta_lg_XH_cntr_max = 0.0d0 + delta_lg_XH_cntr_min = -2.0d0 + !delta_lg_XH_cntr_hard_limit = 0.02d0 + + delta_lg_XHe_cntr_limit = 0.01d0 + delta_lg_XHe_cntr_max = 0.0d0 + delta_lg_XHe_cntr_min = -2.0d0 + !delta_lg_XHe_cntr_hard_limit = -1!0.02d0 + +! time step controls below are useful for resolution on fuel depletion + !delta_lg_XC_cntr_limit = 0.01d0 + !delta_lg_XC_cntr_max = 0.0d0 + !delta_lg_XC_cntr_min = -2.0d0 + !delta_lg_XC_cntr_hard_limit = -1!0.2d0 + + !delta_lg_XNe_cntr_limit = 0.01d0 + !delta_lg_XNe_cntr_max = 0.0d0 + !delta_lg_XNe_cntr_min = -2.0d0 + !delta_lg_XNe_cntr_hard_limit = 0.02 !0.1d0!0.02d0 can crash run + !delta_XNe_cntr_limit = 0.01 + + !delta_lg_XO_cntr_limit = 0.01d0 + !delta_lg_XO_cntr_max = 0.0d0 + !delta_lg_XO_cntr_min = -2.0d0 + !delta_lg_XO_cntr_hard_limit = -1 !0.1d0!0.02d0 can crash run + + delta_XSi_cntr_limit = 0.01 + delta_XSi_cntr_hard_limit = -1!0.02 + + !delta_lg_XSi_cntr_limit = 0.01d0 + !delta_lg_XSi_cntr_max = 0.0d0 + !delta_lg_XSi_cntr_min = -2.0d0 + !delta_lg_XSi_cntr_hard_limit = -1 !0.1d0!0.02d0 can crash run + + !dX_limit_species(3) = 'fe56' + !dX_limit(3) = 0.1 + !dX_div_X_limit_min_X(3) = 1d-5 + !dX_div_X_limit(3) = 0.5d0 + + delta_Ye_highT_limit = 1d-3 + + +! mesh + !max_dq= 1d-3 ! or lower + mesh_delta_coeff = 2.5 ! try 1.0 or below in production run + mesh_delta_coeff_for_highT = 1.5 ! try 1.0 or below in production run + logT_max_for_standard_mesh_delta_coeff = 9.0 + logT_min_for_highT_mesh_delta_coeff = 9.5 + !min_dq_for_xa = 1d-4 ! avoid over-resolving composition changes, bad for bit for bit convergence + !remesh_dt_limit = 1728000 ! 20 days. turn off remesh when dt smaller than this + +! solver + + ! damped newton and structure only + scale_max_correction = 0.1d0 + ignore_species_in_max_correction = .true. + + ! might help to turn off gold2 in 8-10 Msun degenerate cores + ! for O-Ne flames. + use_gold2_tolerances = .true. + gold2_tol_max_residual2 = 5d-7 + gold2_tol_max_residual3 = 5d-4 + + + use_gold_tolerances = .true. + gold_tol_max_residual2 = 5d-4 + !gold_tol_max_residual3 = 1d-4 ! Default in controls is 1d-5 + tol_correction_high_T_limit = 1d9 ! Switch to lower tol at high temp for large Mass + solver_iters_timestep_limit = 20 + gold_solver_iters_timestep_limit = 20 + iter_for_resid_tol2 = 10 + + max_abs_rel_run_E_err = 1d-2 + + energy_eqn_option = 'dedt' + max_tries_for_implicit_wind = 0 ! Recommend 10 for a production run + + convergence_ignore_equL_residuals = .true. + make_gradr_sticky_in_solver_iters = .true. + xa_scale = 1d-5 + iter_for_resid_tol2 = 10 + min_timestep_limit = 1d-12 ! (seconds) ! 1d-20 if things are sticky + + warn_rates_for_high_temp = .true. + max_safe_logT_for_rates = 10.5d0 + + when_to_stop_rtol = 1d-3 + when_to_stop_atol = 1d-3 + + sig_min_factor_for_high_Tcenter = 0.01 + Tcenter_min_for_sig_min_factor_full_on = 3.2d9 + Tcenter_max_for_sig_min_factor_full_off = 2.8d9 + + ! Can be helpful to decrease op_split_burn_min_T + ! 4d9 has been well tested with the approx21 network + ! lower to 2.5d9 or worst case 1d9 for large networks to help with + ! numerical stability and speed, see MESA VI (Jermyn 2023) + op_split_burn = .true. + op_split_burn_min_T = 4d9 + burn_steps_limit = 150 + burn_steps_hard_limit = 250 + op_split_burn_eps = 1d-5 + op_split_burn_odescal = 1d-5 + +! output + terminal_show_log_dt = .false. + + !max_model_number = 4000000 ! if you're serious + + photo_interval = 200 !1000 + photo_digits = 8 + profile_interval = 100 + max_num_profile_models = 400000 + history_interval = 1 + write_header_frequency = 10 + terminal_interval = 10 + + !report_solver_progress = .true. ! set true to see info about solver iterations + !report_ierr = .true. ! if true, produce terminal output when have some internal error + + + + + +! GYRE output controls + !write_pulse_data_with_profile = .true. + pulse_data_format = 'GYRE' + + x_logical_ctrl(37) = .false. ! if true, then run GYRE + + x_integer_ctrl(1) = 1 ! output GYRE info at this step interval + x_logical_ctrl(1) = .false. ! save GYRE info whenever save profile + + x_integer_ctrl(2) = 2 ! max number of modes to output per call + x_logical_ctrl(2) = .false. ! output eigenfunction files + + x_integer_ctrl(3) = 0 ! mode l (e.g. 0 for p modes, 1 for g modes) + ! should match gyre.in mode l + x_integer_ctrl(4) = 1 ! order + x_ctrl(1) = 0.158d-05 ! freq ~ this (Hz) + x_ctrl(2) = 0.33d+03 ! growth < this (days) + + +/ +&pgstar + +!pause = .true. + +!pgstar_interval = 1 + +! x-axis limits and properties +Profile_Panels3_xaxis_name = 'mass' +Profile_Panels3_xmin = 0.0 +Profile_Panels3_xmax = 1.6 +Profile_Panels3_xaxis_reversed = .false. + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams new file mode 100644 index 000000000..89ed7ec43 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams @@ -0,0 +1,57 @@ +! inlist_make_late_pre_zams - until radiative core starts to develop + +&star_job + + show_log_description_at_start = .false. + + create_pre_main_sequence_model = .true. + pre_ms_relax_to_start_radiative_core = .true. + + save_model_when_terminate = .true. + save_model_filename = 'late_pre_zams.mod' + required_termination_code_string = 'max_age' + + !pgstar_flag = .true. + +/ ! end of star_job namelist + +&eos +/ ! end of eos namelist + +&kap + kap_file_prefix = 'a09' ! 'gs98' 'a09' 'OP_a09' 'OP_gs98' + kap_CO_prefix = 'a09_co' ! 'gs98_co' 'a09_co' + use_Type2_opacities = .true. +/ ! end of kap namelist + +&controls + + ! limit max_model_number as part of test_suite + max_model_number = 10000 + use_gold2_tolerances = .true. + + max_age = 1d3 + + num_trace_history_values = 4 + trace_history_value_name(1) = 'conv_mx1_bot' + trace_history_value_name(2) = 'conv_mx1_top' + trace_history_value_name(3) = 'rel_E_err' + trace_history_value_name(4) = 'log_rel_run_E_err' + + photo_interval = 200!1000 + photo_digits = 8 + profile_interval = 100 + max_num_profile_models = 400000 + history_interval = 1 + write_header_frequency = 10 + terminal_interval = 10 + + x_integer_ctrl(5) = 1 ! Inlist part number + +/ ! end of controls namelist + +&pgstar + + +/ ! end of pgstar namelist + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams_header new file mode 100644 index 000000000..dd8aee7e7 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_make_late_pre_zams_header @@ -0,0 +1,36 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_make_late_pre_zams' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_make_late_pre_zams' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_make_late_pre_zams' +/ ! end of kap namelist + + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_make_late_pre_zams' +/ ! end of controls namelist + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_make_late_pre_zams' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_mass_Z_wind_rotation b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_mass_Z_wind_rotation new file mode 100644 index 000000000..56468203c --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_mass_Z_wind_rotation @@ -0,0 +1,18 @@ +! inlist_mass_Z_wind_rotation + +&star_job + new_Z = 1.42d-2 + new_omega_div_omega_crit = 0d0 + initial_zfracs = 6 ! set to asplund 2009 abundance fractions +/ ! end of star_job namelist + +&kap + Zbase = 1.42d-2 +/ ! end of kap namelist + +&controls + initial_mass = 20 + initial_z = 1.42d-2 ! ! solar + initial_y = 0.2703d0 ! from asplund et al. 2009, solar. + Dutch_scaling_factor = 0d0 +/ ! end of controls namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_pgstar b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_pgstar new file mode 100644 index 000000000..cd1b476b7 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_pgstar @@ -0,0 +1,684 @@ + +&pgstar + +!pause = .true. + +pgstar_interval = 1 +pgstar_show_age_in_years = .true. +pgstar_show_age_in_seconds = .false. +pgstar_sleep = 0.0 + +! some global grid plot settings at end + +pgstar_show_model_number = .false. +pgstar_show_age = .false. + +!------------------------------------------------------------------------------------ + +Grid1_win_flag = .true. +Grid1_win_width = 12 +Grid1_win_aspect_ratio = 0.666 + +! file output +Grid1_file_flag = .true. +Grid1_file_dir = 'png' +Grid1_file_prefix = 'Grid1_' +Grid1_file_interval = 25 ! output when mod(model_number,Grid1_file_interval)==0 +Grid1_file_width = 27 ! (inches) negative means use same value as for window +Grid1_file_aspect_ratio = -1 ! negative means use same value as for window + +! reset the defaults + +Grid1_plot_name(:) = '' +Grid1_plot_row(:) = 1 ! number from 1 at top +Grid1_plot_rowspan(:) = 1 ! plot spans this number of rows +Grid1_plot_col(:) = 1 ! number from 1 at left +Grid1_plot_colspan(:) = 1 ! plot spans this number of columns +Grid1_plot_pad_left(:) = 0.0 ! fraction of full window width for padding on left +Grid1_plot_pad_right(:) = 0.0 ! fraction of full window width for padding on right +Grid1_plot_pad_top(:) = 0.0 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(:) = 0.0 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(:) = 0.7 ! multiply txt_scale for subplot by this + + +Grid1_title = '' + +Grid1_num_cols = 3 ! divide plotting region into this many equal width cols +Grid1_num_rows = 5 ! divide plotting region into this many equal height rows +Grid1_num_plots = 6 ! <= 10 + + +Grid1_plot_name(1) = 'Text_Summary1' +Grid1_plot_row(1) = 1 ! number from 1 at top +Grid1_plot_rowspan(1) = 1 ! plot spans this number of rows +Grid1_plot_col(1) = 1 ! number from 1 at left +Grid1_plot_colspan(1) = 3 ! plot spans this number of columns + +Grid1_plot_pad_left(1) = -0.03 ! fraction of full window width for padding on left +Grid1_plot_pad_right(1) = 0.0 ! fraction of full window width for padding on right +Grid1_plot_pad_top(1) = -0.06 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(1) = 0.07 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(1) = 1 ! 0.8 ! multiply txt_scale for subplot by this + + +Grid1_plot_name(3) = 'TRho' +Grid1_plot_row(3) = 2 ! number from 1 at top +Grid1_plot_rowspan(3) = 1 ! plot spans this number of rows +Grid1_plot_col(3) = 1 ! number from 1 at left +Grid1_plot_colspan(3) = 1 ! plot spans this number of columns + +Grid1_plot_pad_left(3) = 0.00 ! fraction of full window width for padding on left +Grid1_plot_pad_right(3) = 0.08 ! fraction of full window width for padding on right +Grid1_plot_pad_top(3) = -0.04 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(3) = 0.01 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(3) = 0.7 ! multiply txt_scale for subplot by this + + +Grid1_plot_name(5) = 'Profile_Panels1' +Grid1_plot_row(5) = 3 ! number from 1 at top +Grid1_plot_rowspan(5) = 3 ! plot spans this number of rows +Grid1_plot_col(5) = 1 ! number from 1 at left +Grid1_plot_colspan(5) = 1 ! plot spans this number of columns + +Grid1_plot_pad_left(5) = 0.00 ! fraction of full window width for padding on left +Grid1_plot_pad_right(5) = 0.10 ! fraction of full window width for padding on right +Grid1_plot_pad_top(5) = 0.05 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(5) = 0 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(5) = 0.6 ! multiply txt_scale for subplot by this + + +Grid1_plot_name(2) = 'TRho_Profile' +Grid1_plot_row(2) = 2 ! number from 1 at top +Grid1_plot_rowspan(2) = 1 ! plot spans this number of rows +Grid1_plot_col(2) = 2 ! number from 1 at left +Grid1_plot_colspan(2) = 1 ! plot spans this number of columns + +Grid1_plot_pad_left(2) = -0.01 ! fraction of full window width for padding on left +Grid1_plot_pad_right(2) = 0.06 ! fraction of full window width for padding on right +Grid1_plot_pad_top(2) = -0.04 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(2) = 0.01 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(2) = 0.7 ! multiply txt_scale for subplot by this + + +Grid1_plot_name(4) = 'History_Panels1' +Grid1_plot_row(4) = 3 ! number from 1 at top +Grid1_plot_rowspan(4) = 3 ! plot spans this number of rows +Grid1_plot_col(4) = 2 ! number from 1 at left +Grid1_plot_colspan(4) = 1 ! plot spans this number of columns + +Grid1_plot_pad_left(4) = 0.00 ! fraction of full window width for padding on left +Grid1_plot_pad_right(4) = 0.06 ! fraction of full window width for padding on right +Grid1_plot_pad_top(4) = 0.05 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(4) = 0.0 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(4) = 0.6 ! multiply txt_scale for subplot by this + + +Grid1_plot_name(6) = 'Profile_Panels3' +Grid1_plot_row(6) = 2 ! number from 1 at top +Grid1_plot_rowspan(6) = 4 ! plot spans this number of rows +Grid1_plot_col(6) = 3 ! Number from 1 at left +Grid1_plot_colspan(6) = 1 ! plot spans this number of columns + +Grid1_plot_pad_left(6) = 0.04 ! fraction of full window width for padding on left +Grid1_plot_pad_right(6) = 0.06 ! fraction of full window width for padding on right +Grid1_plot_pad_top(6) = -0.04 ! fraction of full window height for padding at top +Grid1_plot_pad_bot(6) = 0 ! fraction of full window height for padding at bottom +Grid1_txt_scale_factor(6) = 0.6 ! multiply txt_scale for subplot by this + + +!------------------------------------------------------------------------------------ + +Profile_Panels3_win_flag = .false. + +Profile_Panels3_title = '' + +Profile_Panels3_num_panels = 5 + +Profile_Panels3_yaxis_name(1) = 'Abundance' + +Profile_Panels3_yaxis_name(2) = 'Power' + +Profile_Panels3_yaxis_name(3) = 'Mixing' +Mixing_legend_txt_scale_factor = 0.9 + +Profile_Panels3_yaxis_name(4) = 'logRho' +Profile_Panels3_other_yaxis_name(4) = 'v_div_cs' ! 'vel_km_per_s' ! 'entropy' +Profile_Panels3_other_dymin(4) = 0.14 + +Profile_Panels3_yaxis_name(5) = 'logT' +Profile_Panels3_other_yaxis_name(5) = 'burn_num_iters' + +! x-axis limits and properties +Profile_Panels3_xaxis_name = 'mass' +Profile_Panels3_xmin = 0.0 +Profile_Panels3_xmax = -101d0 ! 2.2 +Profile_Panels3_xaxis_reversed = .false. + + +Profile_Panels3_txt_scale = 0.7 + +!Profile_Panels3_xaxis_name = 'zone' +!Profile_Panels3_xmin = 800 +!Profile_Panels3_xmax = 1100 +!Profile_Panels3_xaxis_reversed = .true. + +!Profile_Panels3_show_grid = .true. +Profile_Panels3_show_mix_regions_on_xaxis = .true. + +!------------------------------------------------------------------------------------ + +!Profile_Panels1_win_flag = .true. +!Profile_Panels1_file_flag = .true. + Profile_Panels1_file_dir = 'png' + Profile_Panels1_file_prefix = 'profile_panels1_' + Profile_Panels1_file_interval = 1 + Profile_Panels1_file_width = -1 + Profile_Panels1_file_aspect_ratio = -1 + +Profile_Panels1_title = '' + +Profile_Panels1_txt_scale = 0.7 +Profile_Panels1_num_panels = 4 + +Profile_Panels1_yaxis_name(1) = 'logRho' +Profile_Panels1_dymin(1) = 0.14 +Profile_Panels1_other_yaxis_name(1) = 'logT' +Profile_Panels1_other_dymin(1) = 0.14 + +Profile_Panels1_yaxis_name(2) = 'entropy' +Profile_Panels1_other_yaxis_name(2) = 'pgas_div_p' + +Profile_Panels1_yaxis_name(3) = 'gradT' +Profile_Panels1_other_yaxis_name(3) = 'grada' +Profile_Panels1_same_yaxis_range(3) = .true. +Profile_Panels1_other_dymin(3) = 0.08 + + +Profile_Panels1_yaxis_name(4) = 'log_opacity' +Profile_Panels1_dymin(1) = 0.14 +Profile_Panels1_other_yaxis_name(4) = 'gradr' +Profile_Panels1_other_dymin(1) = 0.14 + +! x-axis limits and properties +Profile_Panels1_xaxis_name = 'logtau' +Profile_Panels1_xmin = -101d0 +Profile_Panels1_xmax = 8.1 +Profile_Panels1_xaxis_reversed = .true. + +!Profile_Panels1_xaxis_name = 'zone' +!Profile_Panels1_xmin = 15 +!Profile_Panels1_xmax = 270 +!Profile_Panels1_xaxis_reversed = .true. + +!Profile_Panels1_show_grid = .true. +Profile_Panels1_show_mix_regions_on_xaxis = .true. + +!------------------------------------------------------------------------------------ + + +!TRho_Profile_win_flag = .true. +TRho_Profile_win_width = 8 +TRho_Profile_win_aspect_ratio = 0.75 ! aspect_ratio = height/width + +! file output +!TRho_Profile_file_flag = .true. +TRho_Profile_file_dir = 'TRho' +TRho_Profile_file_prefix = 'trho_' +TRho_Profile_file_interval = 10 ! output when `mod(model_number,TRho_Profile_file_interval)==0` +TRho_Profile_file_width = -1 ! (inches) negative means use same value as for window +TRho_Profile_file_aspect_ratio = -1 ! negative means use same value as for window + +TRho_Profile_xleft = 0.15 +TRho_Profile_xright = 0.85 +TRho_Profile_ybot = 0.15 +TRho_Profile_ytop = 0.85 +TRho_Profile_txt_scale = 0.7 +TRho_Profile_title = ' ' + +TRho_switch_to_Column_Depth = .false. +TRho_switch_to_mass = .false. + +show_TRho_Profile_legend = .false. + TRho_Profile_legend_coord = 0.07 + TRho_Profile_legend_fjust = 0.0 + TRho_Profile_legend_disp1 = -2.0 + TRho_Profile_legend_del_disp = -1.3 + TRho_Profile_legend_txt_scale = 1.1 + + +show_TRho_Profile_text_info = .false. + TRho_Profile_text_info_xfac = 0.77 ! controls x location + TRho_Profile_text_info_dxfac = 0.02 ! controls x spacing to value from text + TRho_Profile_text_info_yfac = 0.6 ! controls y location of 1st line + TRho_Profile_text_info_dyfac = -0.04 ! controls line spacing + +show_TRho_Profile_mass_locs = .false. +show_TRho_accretion_mesh_borders = .false. +show_TRho_Profile_kap_regions = .false. +show_TRho_Profile_gamma1_4_3rd = .true. +show_TRho_Profile_eos_regions = .false. +show_TRho_Profile_degeneracy_line = .true. +show_TRho_Profile_Pgas_Prad_line = .true. +show_TRho_Profile_burn_lines = .true. +show_TRho_Profile_burn_labels = .true. + +! axis limits +TRho_Profile_xmin = -12.0 +TRho_Profile_xmax = 10.0 +TRho_Profile_ymin = 3.0 +TRho_Profile_ymax = 10.0 + +! these are shown if show_TRho_Profile_mass_locs = .true. +! set all the entries +profile_mass_point_q = -1 +profile_mass_point_color_index = 1 +profile_mass_point_symbol = -6 +profile_mass_point_symbol_scale = 1.7 +profile_mass_point_str = '' +profile_mass_point_str_clr = 1 +profile_mass_point_str_scale = 1.0 + +! set defaults +num_profile_mass_points = 3 ! max is defined in star_def (max_num_profile_mass_points) + +profile_mass_point_q(1) = 0.5 +profile_mass_point_color_index(1) = 1 +profile_mass_point_symbol(1) = -6 +profile_mass_point_str(1) = ' 0.5 M\d\(0844)\u' +profile_mass_point_str_clr(1) = 1 + +profile_mass_point_q(2) = 0.95 +profile_mass_point_color_index(2) = 1 +profile_mass_point_symbol(2) = -6 +profile_mass_point_str(2) = ' 0.95 M\d\(0844)\u' +profile_mass_point_str_clr(3) = 1 + +profile_mass_point_q(3) = 0.999 +profile_mass_point_color_index(3) = 1 +profile_mass_point_symbol(3) = -6 +profile_mass_point_str(3) = ' 0.999 M\d\(0844)\u' +profile_mass_point_str_clr(3) = 1 + +!------------------------------------------------------------------------------------ + + +! Text_Summary windows + +Text_Summary1_win_flag = .false. +Text_Summary1_win_width = 10 +Text_Summary1_win_aspect_ratio = 0.15 + +Text_Summary1_xleft = 0.01 +Text_Summary1_xright = 0.99 +Text_Summary1_ybot = 0.0 +Text_Summary1_ytop = 1.0 +Text_Summary1_txt_scale = 0.95 +Text_Summary1_title = '' + +Text_Summary1_num_rows = 6 ! <= 20 +Text_Summary1_num_cols = 5 ! <= 20 +Text_Summary1_name(:,:) = '' + + +Text_Summary1_name(1,1) = 'model_number' +Text_Summary1_name(1,2) = 'log_dt' +Text_Summary1_name(1,3) = 'Mass' +Text_Summary1_name(1,4) = 'H_cntr' +Text_Summary1_name(1,5) = 'H_rich' + +Text_Summary1_name(2,1) = 'non_fe_core_infall' +Text_Summary1_name(2,2) = 'star_age' +Text_Summary1_name(2,3) = 'lg_Mdot' +Text_Summary1_name(2,4) = 'He_cntr' +Text_Summary1_name(2,5) = 'He_core' + +Text_Summary1_name(3,1) = 'fe_core_infall' +Text_Summary1_name(3,2) = 'gamma1_min' +Text_Summary1_name(3,3) = 'eta_cntr' +Text_Summary1_name(3,4) = 'C_cntr' +Text_Summary1_name(3,5) = 'CO_core' + +Text_Summary1_name(4,1) = 'log_max_T' +Text_Summary1_name(4,2) = 'log_LH' +Text_Summary1_name(4,3) = 'lg_Lnuc' +Text_Summary1_name(4,4) = 'O_cntr' +Text_Summary1_name(4,5) = 'Fe_core' + +Text_Summary1_name(5,1) = 'log_cntr_T' +Text_Summary1_name(5,2) = 'log_LHe' +Text_Summary1_name(5,3) = 'lg_Lneu' +Text_Summary1_name(5,4) = 'Ne_cntr' +Text_Summary1_name(5,5) = 'zones' + +Text_Summary1_name(6,1) = 'log_cntr_Rho' +Text_Summary1_name(6,2) = 'log_LZ' +Text_Summary1_name(6,3) = 'lg_Lphoto' +Text_Summary1_name(6,4) = 'Si_cntr' +Text_Summary1_name(6,5) = 'retries' + + +!------------------------------------------------------------------------------------ + +! Abundance profile plot + +Abundance_win_flag = .false. + +! window properties +Abundance_win_width = 10 +Abundance_win_aspect_ratio = 0.75 + +Abundance_xleft = 0.15 +Abundance_xright = 0.85 +Abundance_ybot = 0.15 +Abundance_ytop = 0.85 +Abundance_txt_scale = 1.1 +Abundance_title = '' + +! isotopes to plot + +Abundance_num_isos_to_show = 20 + +Abundance_which_isos_to_show(1) = 'h1' +Abundance_which_isos_to_show(2) = 'he3' +Abundance_which_isos_to_show(3) = 'he4' +Abundance_which_isos_to_show(4) = 'c12' +Abundance_which_isos_to_show(5) = 'n14' +Abundance_which_isos_to_show(6) = 'o16' +Abundance_which_isos_to_show(7) = 'ne20' +Abundance_which_isos_to_show(8) = 'mg24' +Abundance_which_isos_to_show(9) = 'si28' +Abundance_which_isos_to_show(10) = 's32' +Abundance_which_isos_to_show(11) = 'ar36' +Abundance_which_isos_to_show(12) = 'ca40' +Abundance_which_isos_to_show(13) = 'ti44' +Abundance_which_isos_to_show(14) = 'cr48' +Abundance_which_isos_to_show(15) = 'cr56' +Abundance_which_isos_to_show(16) = 'fe52' +Abundance_which_isos_to_show(17) = 'fe54' +Abundance_which_isos_to_show(18) = 'fe56' +Abundance_which_isos_to_show(19) = 'ni56' +Abundance_which_isos_to_show(20) = 'neut' +!Abundance_which_isos_to_show(22) = 'ne22' + + + +! number and size of isotope labels along curves +num_abundance_line_labels = 4 +Abundance_line_txt_scale_factor = 1.1 + + +! number and size of isotopes on legend +Abundance_legend_max_cnt = 10 +Abundance_legend_txt_scale_factor = 1.3 + +! yaxis limits +Abundance_log_mass_frac_min = -4!-6.4 ! -3.5 +Abundance_log_mass_frac_max = 0.3 + +! file output +Abundance_file_flag = .false. +Abundance_file_dir = 'Abundance' +Abundance_file_prefix = 'abund_' +Abundance_file_width = -1 ! (inches) negative means use same value as for window +Abundance_file_aspect_ratio = -1 ! negative means use same value as for window + + +!------------------------------------------------------------------------------------ + +! power plot + +Power_win_flag = .false. +Power_win_width = 10 +Power_win_aspect_ratio = 0.75 +Power_title = '' + +Power_xleft = 0.15 +Power_xright = 0.85 +Power_ybot = 0.15 +Power_ytop = 0.85 +Power_txt_scale = 1.1 +Power_title = ' ' + +Power_legend_max_cnt = 10 +Power_legend_txt_scale_factor = 1.3 ! relative to other text + +! power yaxis limits -- to override system default selections +Power_ymin = -5.0 ! -101d0 ! only used if /= -101d0 +Power_ymax = 25.0 ! -101d0 ! only used if /= -101d0 + +! file output +Power_file_flag = .false. +Power_file_dir = 'png' +Power_file_prefix = 'power_' +Power_file_interval = 5 ! output when mod(model_number,Power_file_interval)==0 +Power_file_width = -1 ! (inches) negative means use same value as for window +Power_file_aspect_ratio = -1 ! negative means use same value as for window + + +!------------------------------------------------------------------------------------ + +! mixing plot + +Mixing_xmin = 0.0 +Mixing_xmax = 1.6 ! -101d0 +Mixing_legend_txt_scale_factor = 1.4 ! relative to other text + +Mixing_show_rotation_details = .false. + +!Mixing_win_flag = .true. +!Mixing_file_flag = .true. +Mixing_file_dir = 'png' +Mixing_file_prefix = 'mixing_' +Mixing_file_interval = 1 ! output when `mod(model_number,Mixing_file_interval)==0` +Mixing_file_width = -1 ! (inches) negative means use same value as for window +Mixing_file_aspect_ratio = -1 ! negative means use same value as for window + + +!----------------------------------------------------------------------- + +! TRho window + ! history of central temperature vs. density + + TRho_txt_scale = 0.7 + TRho_title = '' + + TRho_logT_min = -101d0 + TRho_logT_max = -101d0 + TRho_logRho_min = -101d0 + TRho_logRho_max = -101d0 + show_TRho_degeneracy_line = .true. + + + +!----------------------------------------------------------------------- + + !# HR window + ! history of `lg_L` vs. `lg_Teff` + + HR_win_flag = .true. + + HR_win_width = 6 + HR_win_aspect_ratio = 0.75 ! aspect_ratio = height/width + + HR_xleft = 0.15 + HR_xright = 0.85 + HR_ybot = 0.15 + HR_ytop = 0.85 + HR_txt_scale = 0.7 !1.0 + HR_title = '' + + ! axis limits -- to override system default selections + HR_logT_min = -101d0 ! only used if /= -101d0 + HR_logT_max = -101d0 ! only used if /= -101d0 + HR_logL_min = -101d0 ! only used if /= -101d0 + HR_logL_max = -101d0 ! only used if /= -101d0 + + HR_logL_margin = 0.1 + HR_logT_margin = 0.1 + HR_dlogT_min = -1 + HR_dlogL_min = -1 + + HR_step_min = -1 ! only plot models with model number >= this + HR_step_max = -1 ! only plot models with model number <= this + + show_HR_classical_instability_strip = .true. + show_HR_Mira_instability_region = .false. + show_HR_WD_instabilities = .false. + + show_HR_target_box = .false. + HR_target_n_sigma = -3 ! -n means show sig 1..n + HR_target_logL = 0 + HR_target_logL_sigma = 0 + HR_target_logT = 0 + HR_target_logT_sigma = 0 + + show_HR_annotation1 = .false. + show_HR_annotation2 = .false. + show_HR_annotation3 = .false. + + HR_fname = '' ! file name for extra HR data + + ! Enables calling a subroutine to add extra information to a plot + ! see `$MESA_DIR/star/other/pgstar_decorator.f90` + HR_use_decorator = .false. + + + ! file output + HR_file_flag = .false. + HR_file_dir = 'png' + HR_file_prefix = 'hr_' + HR_file_interval = 5 ! output when `mod(model_number,HR_file_interval)==0` + HR_file_width = -1 ! (inches) negative means use same value as for window + HR_file_aspect_ratio = -1 ! negative means use same value as for window + +!----------------------------------------------------------------------- + + History_Panels1_title = '' + + History_Panels1_xaxis_name = 'model_number' + History_Panels1_max_width = -1 + + !History_Panels1_xaxis_name = 'star_age' + !History_Panels1_max_width = 10 + + History_Panels1_txt_scale = 0.75 + History_Panels1_xmin = -101d0 + History_Panels1_xmax = -101d0 + History_Panels1_dxmin = -1 + History_Panels1_xaxis_reversed = .false. + History_Panels1_xaxis_log = .false. + History_Panels1_xmargin = 0.0 + + ! :: + + History_Panels1_num_panels = 4 + + ! :: + + History_Panels1_yaxis_name(1) = 'log_L' + History_Panels1_yaxis_reversed(1) = .false. + History_Panels1_ymin(1) = -101d0 + History_Panels1_ymax(1) = -101d0 + History_Panels1_dymin(1) = 0.14 + + ! :: + + History_Panels1_other_yaxis_name(1) = 'log_Teff' + History_Panels1_other_yaxis_reversed(1) = .false. + History_Panels1_other_ymin(1) = -101d0 + History_Panels1_other_ymax(1) = -101d0 + History_Panels1_other_dymin(1) = 0.14 + + ! :: + + History_Panels1_yaxis_name(2) = 'lum_div_Ledd' + History_Panels1_yaxis_reversed(2) = .false. + History_Panels1_ymin(2) = -101d0 + History_Panels1_ymax(2) = -101d0 + History_Panels1_dymin(2) = 0.14 + + ! :: + + History_Panels1_other_yaxis_name(2) = 'log_max_T' ! 'v_surf_km_s' + History_Panels1_other_yaxis_reversed(2) = .false. + History_Panels1_other_ymin(2) = -101d0 + History_Panels1_other_ymax(2) = -101d0 + History_Panels1_other_dymin(2) = 0.14 + + ! :: + + History_Panels1_yaxis_name(3) = 'radius' + History_Panels1_yaxis_reversed(3) = .false. + History_Panels1_ymin(3) = -101d0 + History_Panels1_ymax(3) = -101d0 + History_Panels1_dymin(3) = 0.14 + + ! :: + + History_Panels1_other_yaxis_name(3) = 'log_cntr_Rho' + History_Panels1_other_yaxis_reversed(3) = .false. + History_Panels1_other_ymin(3) = -101d0 + History_Panels1_other_ymax(3) = -101d0 + History_Panels1_other_dymin(3) = 0.14 + ! :: + + History_Panels1_yaxis_name(4) = '' + History_Panels1_yaxis_reversed(4) = .false. + History_Panels1_ymin(4) = -101d0 + History_Panels1_ymax(4) = -101d0 + History_Panels1_dymin(4) = 0.14 + + ! :: + + History_Panels1_other_yaxis_name(4) = 'log_dt' + History_Panels1_other_yaxis_reversed(4) = .false. + History_Panels1_other_ymin(4) = -101d0 + History_Panels1_other_ymax(4) = -101d0 + History_Panels1_other_dymin(4) = 0.14 + + +!----------------------------------------------------------------------- + +! some global grid plot settings + +pgstar_grid_show_title = .true. +pgstar_grid_title_scale = 1.0 +pgstar_grid_title_lw = 3 +pgstar_grid_title_disp = 2.5 ! 1.8 +pgstar_grid_title_coord = 0.5 +pgstar_grid_title_fjust = 0.5 + +pgstar_age_scale = 0.8 +pgstar_age_disp = 3.0 +pgstar_age_coord = 0.0 +pgstar_age_fjust = 0.0 + +pgstar_xaxis_label_scale = 1.3 +pgstar_left_yaxis_label_scale = 1.3 +pgstar_xaxis_label_disp = 2.2 +pgstar_left_yaxis_label_disp = 3.1 +pgstar_right_yaxis_label_disp = 4.1 + +pgstar_model_scale = 0.8 +pgstar_model_disp = 3.0 +pgstar_model_coord = 1.0 +pgstar_model_fjust = 1.0 + +! white_on_black flags -- true means white foreground color on black background +file_white_on_black_flag = .true. +file_device = 'png' ! options 'png' and 'vcps' for png and postscript respectively +file_extension = 'png' ! common names are 'png' and 'ps' + + +!file_white_on_black_flag = .false. +!file_device = 'vcps' ! options 'png' and 'vcps' for png and postscript respectively +!file_extension = 'ps' ! common names are 'png' and 'ps' + + +kipp_win_flag=.true. +kipp_file_flag=.true. +Kipp_mix_interval = 1 +Kipp_show_luminosities = .true. + +/ ! end of pgstar namelist + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope new file mode 100644 index 000000000..6d787b9fe --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope @@ -0,0 +1,79 @@ + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'after_core_he_burn.mod' + + save_model_when_terminate = .true. + save_model_filename = 'removed_envelope.mod' + required_termination_code_string = 'max_model_number' + + set_initial_number_retries = .true. + + ! uncomment the following to remove H env + !relax_initial_mass_to_remove_H_env = .true. + !lg_max_abs_mdot = -2 + + !pgstar_flag = .false. + +/ ! end of star_job namelist + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls + + ! just run for a relatively small number of steps to adjust to removing envelope + max_model_number = 100 + !max_number_retries = 37 + +! wind + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + +! mesh + +! solver + +! output + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_rel_run_E_err' + + !photo_interval = 10 + !profile_interval = 2 + !history_interval = 1 + !terminal_interval = 1 + + x_integer_ctrl(5) = 4 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar + + +Profile_Panels3_xmin = 0 ! 2.5 ! -101d0 +Profile_Panels3_xmax = -101d0 ! +Profile_Panels3_yaxis_name(4) = 'gamma1' + +TRho_Profile_xmin = 2.0 +TRho_Profile_xmax = 7.5 +TRho_Profile_ymin = 8.3 +TRho_Profile_ymax = 9.2 + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope_header new file mode 100644 index 000000000..a71696def --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_remove_envelope_header @@ -0,0 +1,42 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_remove_envelope' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_remove_envelope' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_remove_envelope' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_remove_envelope' +/ ! end of controls namelist + + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_remove_envelope' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc new file mode 100644 index 000000000..52ab1a47d --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc @@ -0,0 +1,202 @@ + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'lgTmax.mod' + + save_model_when_terminate = .true. + save_model_filename = 'final.mod' + required_termination_code_string = 'fe_core_infall_limit' + + +/ ! end of star_job namelist + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls + use_DGESVX_in_bcyclic = .true. + use_equilibration_in_DGESVX = .true. + report_min_rcond_from_DGESXV = .true. + +max_number_retries = 0 + + report_solver_progress = .true. ! set true to see info about solver iterations + !report_ierr = .true. ! if true, produce terminal output when have some internal error + !stop_for_bad_nums = .true. + !trace_evolve = .true. + !fill_arrays_with_NaNs = .true. + + !solver_save_photo_call_number = 0 + ! Saves a photo when solver_call_number = solver_save_photo_call_number - 1 + ! e.g., useful for testing partials to set solver_call_number = solver_test_partials_call_number - 1 + + ! solver_test_partials_call_number = 34 + ! solver_test_partials_k = 1817 + ! solver_test_partials_iter_number = 25 + ! solver_test_partials_dx_0 = 1d-6 + ! solver_test_partials_var_name = 'all' ! 'all' or 'lnd', 'lnT', 'lnR', 'L', 'v', etc. '' means code sets + ! solver_test_partials_equ_name = 'eps_nuc' ! 'all' or 'dlnE_dt', 'dlnd_dt', 'dlnR_dt', 'equL', etc '' means code sets + !solver_test_partials_sink_name = 'si28' ! iso name to use for "sink" to keep sum = 1 + !solver_test_partials_show_dx_var_name = 'h1' + + ! equ name can also be one of these + ! 'lnE', 'lnP', 'grad_ad' to test eos + ! 'eps_nuc' to test net + ! 'non_nuc_neu' to test neu + ! 'gradT', 'mlt_vc' to test mlt + ! 'opacity' to test kap + + !solver_test_partials_write_eos_call_info = .true. + + !solver_test_partials_k_low = -1 + !solver_test_partials_k_high = -1 + + !solver_test_eos_partials = .true. + !solver_test_kap_partials = .true. + !solver_test_net_partials = .true. + !solver_test_atm_partials = .true. + + !report_all_dt_limits = .true. + !report_solver_dt_info = .true. + + !show_mesh_changes = .true. + !mesh_dump_call_number = 5189 + !okay_to_remesh = .false. + + !energy_conservation_dump_model_number = -1 + + ! solver debugging + !solver_check_everything = .true. + + !solver_epsder_struct = 1d-6 + !solver_epsder_chem = 1d-6 + + !report_solver_dt_info = .true. + !report_dX_nuc_drop_dt_limits = .true. + !report_bad_negative_xa = .true. + + + + +! to avoid a host of retries we relax our limits on split_burn +! for speed/stability in test_suite, at reduced time resolution. + op_split_burn = .false. + op_split_burn_min_T = 1d9 + +! prevent development of radial pulses during advanced burning + !drag_coefficient = 1d0 + !min_q_for_drag = 0.8d0 + + set_min_D_mix = .true. + min_D_mix = 1d-2 + + fe_core_infall_limit = 1d7 + + ! limit max_model_number as part of test_suite + max_model_number = 6000 + !max_number_retries = 37 + +ignore_too_large_correction = .true. ! for conv_vel's + + delta_lg_XH_cntr_min = 0.5d0 + delta_lg_XO_cntr_min = 0.5d0 + delta_lg_XC_cntr_min = 0.5d0 + delta_lg_XHe_cntr_min = 0.5d0 + +! wind + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + limit_for_rel_error_in_energy_conservation = 1d-3 + hard_limit_for_rel_error_in_energy_conservation = 1d-2 + + delta_XSi_cntr_limit = 0.0025 + + +! mesh + +! solver + + calculate_Brunt_B = .true. ! needed for tau_conv + max_q_for_conv_timescale = 0.2d0 + max_X_for_conv_timescale = 1d-6 ! must be > 0 + +! GYRE output controls + !write_pulse_data_with_profile = .true. + !x_logical_ctrl(37) = .true. ! if true, then run GYRE + x_integer_ctrl(1) = 20 ! output GYRE info at this step interval + x_logical_ctrl(1) = .false. ! save GYRE info whenever save profile + x_integer_ctrl(2) = 2 ! max number of modes to output per call + x_logical_ctrl(2) = .false. ! output eigenfunction files + x_integer_ctrl(3) = 0 ! mode l (e.g. 0 for p modes, 1 for g modes) + x_integer_ctrl(4) = 1 ! order + x_ctrl(1) = 0.158d-05 ! freq ~ this (Hz) + x_ctrl(2) = 0.33d+03 ! growth < this (days) + +! output + terminal_show_age_units = 'secs' + terminal_show_timestep_units = 'secs' + terminal_show_log_dt = .false. + terminal_show_log_age = .false. + + num_trace_history_values = 5 + trace_history_value_name(1) = 'Fe_core' + trace_history_value_name(2) = 'fe_core_infall' + trace_history_value_name(3) = 'rel_E_err' + trace_history_value_name(4) = 'log_rel_run_E_err' + trace_history_value_name(5) = 'dt_div_max_tau_conv' + + photo_interval = 1 + !profile_interval = 1 + !history_interval = 1 + terminal_interval = 10 + + x_integer_ctrl(5) = 7 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar +!Grid1_file_flag = .true. + +!pause = .true. + +Profile_Panels3_xaxis_name = 'mass' +Profile_Panels3_xmin = 0 ! -101d0 +Profile_Panels3_xmax = 5 ! +Profile_Panels3_xaxis_reversed = .false. + + +Profile_Panels3_yaxis_name(4) = 'gamma1' +Profile_Panels3_dymin(4) = 0.14 +Profile_Panels3_other_yaxis_name(4) = 'vel_km_per_s' +Profile_Panels3_dymin(4) = 0.14 + +TRho_Profile_xmin = -15 !3.0 +TRho_Profile_xmax = 10 +TRho_Profile_ymin = 3.35 !8.5 +TRho_Profile_ymax = 10.1 + +TRho_logT_min = 9.2d0 +TRho_logRho_min = 6d0 +TRho_logRho_max = 10d0 +TRho_logT_max = 10d0 + +Text_Summary1_name(1,2) = 'time_step_sec' +Text_Summary1_name(2,2) = 'star_age_sec' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc_header new file mode 100644 index 000000000..93af411e6 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_cc_header @@ -0,0 +1,42 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_to_cc' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_to_cc' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_to_cc' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_to_cc' +/ ! end of controls namelist + + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_to_cc' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn new file mode 100644 index 000000000..49bf2a878 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn @@ -0,0 +1,78 @@ + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'removed_envelope.mod' + + save_model_when_terminate = .true. + save_model_filename = 'after_core_c_burn.mod' + required_termination_code_string = 'xa_central_lower_limit' + +/ ! end of star_job namelist + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls +! prevent development of radial pulses during advanced burning + !drag_coefficient = 1d0 + !min_q_for_drag = 0.8d0 + + set_min_D_mix = .true. + min_D_mix = 1d-2 + + xa_central_lower_limit_species(1) = 'c12' + xa_central_lower_limit(1) = 1d-3 + + ! limit max_model_number as part of test_suite + max_model_number = 10000 + !max_number_retries = 37 + +! wind + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + +! mesh + +! solver + +! output + terminal_show_age_units = 'years' + terminal_show_timestep_units = 'days' + terminal_show_log_dt = .false. + terminal_show_log_age = .false. + + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_rel_run_E_err' + + x_integer_ctrl(5) = 5 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar + +Profile_Panels3_xmin = 0 ! 2.5 ! -101d0 +Profile_Panels3_xmax = -101d0 ! + +TRho_Profile_xmin = 2.5 +TRho_Profile_xmax = 7.5 +TRho_Profile_ymin = 8.4 +TRho_Profile_ymax = 9.2 + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn_header new file mode 100644 index 000000000..f020a5a14 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_c_burn_header @@ -0,0 +1,44 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_to_end_core_c_burn' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_to_end_core_c_burn' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_to_end_core_c_burn' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_to_end_core_c_burn' +/ ! end of controls namelist + + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_to_end_core_c_burn' + read_extra_pgstar_inlist(3) = .true. + extra_pgstar_inlist_name(3)= 'inlist_to_end_core_c_burn' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn new file mode 100644 index 000000000..aebe741cc --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn @@ -0,0 +1,80 @@ + + +&star_job + + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'zams.mod' + + save_model_when_terminate = .true. + save_model_filename = 'after_core_he_burn.mod' + required_termination_code_string = 'xa_central_lower_limit' + + set_initial_cumulative_energy_error = .true. + new_cumulative_energy_error = 0d0 + + new_rotation_flag = .true. + near_zams_relax_omega_div_omega_crit = .true. + change_rotation_flag = .false. ! rotation off until near zams + num_steps_to_relax_rotation = 50 + +/ ! end of star_job namelist + + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls + + xa_central_lower_limit_species(1) = 'he4' + xa_central_lower_limit(1) = 1d-6 + + ! limit max_model_number as part of test_suite + max_model_number = 3000 + !max_number_retries = 58 + +! wind + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + +! mesh + +! solver + +! output + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_rel_run_E_err' + + x_integer_ctrl(5) = 3 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar + +!pause = .true. +Grid1_plot_name(3) = 'HR' +Profile_Panels3_xmin = -101d0 +Profile_Panels3_xmax = -101d0 +Profile_Panels3_other_yaxis_name(4) = 'entropy' + +TRho_Profile_xmin = -15 +TRho_Profile_xmax = 7.5 +TRho_Profile_ymin = 3.5 +TRho_Profile_ymax = 9.2 + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn_header new file mode 100644 index 000000000..2aca02043 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_end_core_he_burn_header @@ -0,0 +1,41 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_to_end_core_he_burn' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_to_end_core_he_burn' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_to_end_core_he_burn' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_to_end_core_he_burn' +/ ! end of controls namelist + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_to_end_core_he_burn' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax new file mode 100644 index 000000000..5d384c4a5 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax @@ -0,0 +1,86 @@ + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'after_core_c_burn.mod' + + save_model_when_terminate = .true. + save_model_filename = 'lgTmax.mod' + required_termination_code_string = 'log_max_temp_upper_limit' + +/ ! end of star_job namelist + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls +! prevent development of radial pulses during advanced burning + !drag_coefficient = 1d0 + !min_q_for_drag = 0.8d0 + + set_min_D_mix = .true. + min_D_mix = 1d-2 + + log_max_temp_upper_limit = 9.600d0 + + ! limit max_model_number as part of test_suite + max_model_number = 4000 + !max_number_retries = 37 + + limit_for_rel_error_in_energy_conservation = 1d-3 + hard_limit_for_rel_error_in_energy_conservation = 1d-2 + +! wind + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + +! mesh + +! solver + +! output + terminal_show_age_units = 'days' + terminal_show_timestep_units = 'days' + terminal_show_log_dt = .false. + terminal_show_log_age = .false. + + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_rel_run_E_err' + + x_integer_ctrl(5) = 6 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar + +TRho_logT_min = 8.6d0 +TRho_logRho_min = 5d0 +TRho_logRho_max = 10d0 +TRho_logT_max = 10d0 + + +!pause = .true. +Profile_Panels3_xmin = 0 ! -101d0 +Profile_Panels3_xmax = 5 ! +Profile_Panels3_other_yaxis_name(4) = 'vel_km_per_s' + +TRho_Profile_xmin = -12d0 +TRho_Profile_ymin = 3.0d0 + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax_header new file mode 100644 index 000000000..0ab492974 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_lgTmax_header @@ -0,0 +1,42 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_to_lgTmax' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_to_lgTmax' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_to_lgTmax' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_to_lgTmax' +/ ! end of controls namelist + + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_to_lgTmax' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams new file mode 100644 index 000000000..b3bd5e068 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams @@ -0,0 +1,88 @@ + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'late_pre_zams.mod' + + save_model_when_terminate = .true. + save_model_filename = 'zams.mod' + required_termination_code_string = 'xa_central_lower_limit' + + set_initial_cumulative_energy_error = .true. + new_cumulative_energy_error = 0d0 + + new_rotation_flag = .true. + near_zams_relax_omega_div_omega_crit = .true. + change_rotation_flag = .false. ! rotation off until near zams + num_steps_to_relax_rotation = 50 + + pgstar_flag = .false. + +/ ! end of star_job namelist + + +&eos +/ ! end of eos namelist + +&kap +/ ! end of kap namelist + +&controls + + xa_central_lower_limit_species(1) = 'h1' + xa_central_lower_limit(1) = 0.68d0 + + ! limit max_model_number as part of test_suite + max_model_number = 6000 + !max_number_retries = 58 + + !stop_near_zams = .false. + !Lnuc_div_L_zams_limit = 0.4 ! if stop_near_zams + +! wind + ! delay wind until reach zams + Dutch_scaling_factor = 0 + +! atmosphere + +! rotation + +! mlt + +! mixing + +! timesteps + +! mesh + +! solver + solver_iters_timestep_limit = 20 + +! output + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_rel_run_E_err' + + x_integer_ctrl(5) = 2 ! Inlist part number + +/ ! end of controls namelist + + + +&pgstar + +!pause = .true. + +Grid1_plot_name(3) = 'HR' +Profile_Panels3_xmin = -101d0 +Profile_Panels3_xmax = -101d0 +Profile_Panels3_other_yaxis_name(4) = 'entropy' + +!Profile_Panels3_xaxis_name = 'zone' +!Profile_Panels3_xmin = 48 +!Profile_Panels3_xmax = 58 +!Profile_Panels3_xaxis_reversed = .true. + +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams_header b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams_header new file mode 100644 index 000000000..97df1345d --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/inlist_to_zams_header @@ -0,0 +1,41 @@ + +&star_job + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_common' + read_extra_star_job_inlist(2) = .true. + extra_star_job_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_star_job_inlist(3) = .true. + extra_star_job_inlist_name(3) = 'inlist_to_zams' +/ ! end of star_job namelist + +&eos + read_extra_eos_inlist(1) = .true. + extra_eos_inlist_name(1) = 'inlist_common' + read_extra_eos_inlist(3) = .true. + extra_eos_inlist_name(3) = 'inlist_to_zams' +/ ! end of eos namelist + +&kap + read_extra_kap_inlist(1) = .true. + extra_kap_inlist_name(1) = 'inlist_common' + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_mass_Z_wind_rotation' + read_extra_kap_inlist(3) = .true. + extra_kap_inlist_name(3) = 'inlist_to_zams' +/ ! end of kap namelist + +&controls + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_common' + read_extra_controls_inlist(2) = .true. + extra_controls_inlist_name(2)= 'inlist_mass_Z_wind_rotation' + read_extra_controls_inlist(3) = .true. + extra_controls_inlist_name(3)= 'inlist_to_zams' +/ ! end of controls namelist + +&pgstar + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + read_extra_pgstar_inlist(2) = .true. + extra_pgstar_inlist_name(2)= 'inlist_to_zams' +/ ! end of pgstar namelist diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/instructions_for_testing.txt b/star/dev_cases_test_solver/test_20M_near_cc_approx21/instructions_for_testing.txt new file mode 100644 index 000000000..1514d884c --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/instructions_for_testing.txt @@ -0,0 +1,13 @@ +# this is for generating the photo necessary for restarts. a working photo should be available already though. +# 1 -execute "./rn" for one timestep to save photo number x00002000. +# 2 - Make sure "op_split_burn = .false." in inlist_to_cc. Uncomment the first line of the "./test" script and then execute "./test" until model number 2031 is created +# 3 - comment out the first line of "./test" again, and then uncomment the "max_num_retries = 0" control in inlist_to_cc. +# 4 - Comment out first line of "./test", so you are restarting from x00002031 and then execture "./test". + +# After the steps above are complete, make sure "export OMP_NUM_THREADS=1" for testing so the data from jacobian populates the "zone_*.txt"files in order. + + +# once "./test" is run in the last step, it generates five files for each zone across all iterations. +# After this, the plot_eps_nuc_terms.py file can be run to make the plots in the "plots" folder, however this is a work in progress. + +# need more derivatives and information across timesteps. Current output definitions on lines 283+ in star/private/net.f90 in the do1_net subroutine. \ No newline at end of file diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/lgTmax.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/lgTmax.mod new file mode 100644 index 000000000..543e1555d --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/lgTmax.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:808de00ef2b333e36afbbf0a8dda7ebc146e78560f8faf79a8b12ce6288e4acf +size 1428461 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/make/makefile b/star/dev_cases_test_solver/test_20M_near_cc_approx21/make/makefile new file mode 100644 index 000000000..f2017faed --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/make/makefile @@ -0,0 +1,5 @@ + + +STAR = star + +include $(MESA_DIR)/star/work_standard_makefile diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/mk b/star/dev_cases_test_solver/test_20M_near_cc_approx21/mk new file mode 100755 index 000000000..aec7a5195 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/plot_eps_nuc_terms.py b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plot_eps_nuc_terms.py new file mode 100644 index 000000000..edcedb200 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plot_eps_nuc_terms.py @@ -0,0 +1,205 @@ +import numpy as np +import re +import glob + +def read_zone_file(filename): + data = [] + with open(filename, 'r') as file: + content = file.readlines() + + iteration_data = {} + i = 0 + while i < len(content): + line = content[i].strip() + if line.startswith('Iteration:'): + # Start a new iteration data block + if iteration_data: + data.append(iteration_data) + iteration_data = {} + iteration_number = int(line.split(':')[1].strip()) + iteration_data['Iteration'] = iteration_number + i += 1 + elif ':' in line: + key = line.split(':')[0].strip() + if key in ['eps_nuc', 'rho', 'T', 'd_eps_nuc_dRho', 'd_eps_nuc_dT']: + # Scalar variable + value = float(line.split(':')[1].strip()) + iteration_data[key] = value + i += 1 + elif key in ['d_epsnuc_dx', 'dxdt_nuc', 'd_dxdt_nuc_dRho', 'd_dxdt_nuc_dT']: + # Vector variable + i += 1 # Move to the next line where data starts + values = [] + while i < len(content) and not content[i].strip().endswith(':'): + line_values = content[i].strip().split() + values.extend([float(val) for val in line_values]) + i += 1 + iteration_data[key] = np.array(values) + elif key == 'd_dxdt_nuc_dx': + # 2D array variable + i += 1 # Move to the next line where data starts + rows = [] + while i < len(content) and not content[i].strip().startswith(('Iteration:', 'eps_nuc:', 'rho:', 'T:', 'd_eps_nuc_dRho:', 'd_eps_nuc_dT:', 'd_epsnuc_dx:', 'dxdt_nuc:', 'd_dxdt_nuc_dRho:', +'d_dxdt_nuc_dT:')): + line_values = content[i].strip().split() + if line_values: + row = [float(val) for val in line_values] + rows.append(row) + i += 1 + iteration_data[key] = np.array(rows) + else: + # Unknown key + i += 1 + else: + i += 1 + + # Append the last iteration data + if iteration_data: + data.append(iteration_data) + + return data + +def read_all_zones(): + zone_files = sorted(glob.glob('zone_*.txt')) + all_zone_data = {} + for filename in zone_files: + zone_number = re.findall(r'zone_(\d+)\.txt', filename)[0] + zone_data = read_zone_file(filename) + all_zone_data[zone_number] = zone_data + return all_zone_data + +# Example usage +if __name__ == '__main__': + all_data = read_all_zones() + + import matplotlib.pyplot as plt + + plt.figure(figsize=(8, 6)) + for i in [1815,1816,1817,1818,1819]: + # Access data for a specific zone and iteration + zone_number = str(i) # Change as needed + zone_data = all_data[zone_number] + + # For example, get the eps_nuc values over iterations + iterations = [iter_data['Iteration'] for iter_data in zone_data] + eps_nuc_values = [iter_data['eps_nuc'] for iter_data in zone_data] + + # Plotting eps_nuc over iterations for the specified zone + + eps_nuc_values = np.array(eps_nuc_values) + # Initialize an array to store the log-scaled values + eps_nuc_values_log = np.zeros_like(eps_nuc_values) + + # Handle positive terms: directly take the logarithm + positive_mask = eps_nuc_values > 0 + eps_nuc_values_log[positive_mask] = np.log10(eps_nuc_values[positive_mask]) + + # Handle negative terms: take the logarithm of the absolute value and scale by -1 + negative_mask = eps_nuc_values < 0 + eps_nuc_values_log[negative_mask] = -np.log10(np.abs(eps_nuc_values[negative_mask])) + + #name = "zone" +str(i) + plt.plot(iterations,eps_nuc_values_log, marker='o', label = str(i)) + + +# plt.plot(iterations,eps_nuc_values_log, marker='o') +# plt.ylim(-1e21,1e21) + print (iterations) +# print (eps_nuc_values_log) + plt.xlabel('Iteration') + plt.ylabel('eps_nuc') + plt.title(f'eps_nuc over all iterations for core zones') + plt.legend() + plt.grid(True) + plt.savefig('plots/eps_nuc.pdf',dpi = 300) + plt.show() + + + + + + + "plot partials for d_eps_nuc_dT in core zones" + plt.figure(figsize=(8, 6)) + for i in [1815,1816,1817,1818,1819]: + # Access data for a specific zone and iteration + zone_number = str(i) # Change as needed + zone_data = all_data[zone_number] + + # For example, get the eps_nuc values over iterations + iterations = [iter_data['Iteration'] for iter_data in zone_data] + eps_nuc_values = [iter_data['d_eps_nuc_dT'] for iter_data in zone_data] + + # Plotting eps_nuc over iterations for the specified zone + + eps_nuc_values = np.array(eps_nuc_values) + # Initialize an array to store the log-scaled values + eps_nuc_values_log = np.zeros_like(eps_nuc_values) + + # Handle positive terms: directly take the logarithm + positive_mask = eps_nuc_values > 0 + eps_nuc_values_log[positive_mask] = np.log10(eps_nuc_values[positive_mask]) + + # Handle negative terms: take the logarithm of the absolute value and scale by -1 + negative_mask = eps_nuc_values < 0 + eps_nuc_values_log[negative_mask] = -np.log10(np.abs(eps_nuc_values[negative_mask])) + + #name = "zone" +str(i) + plt.plot(iterations,eps_nuc_values_log, marker='o', label = str(i)) + + +# plt.plot(iterations,eps_nuc_values_log, marker='o') +# plt.ylim(-1e21,1e21) + print (iterations) +# print (eps_nuc_values_log) + plt.xlabel('Iteration') + plt.ylabel('d_eps_nuc_dT') + plt.title(f'd_eps_nuc_dT over all iterations for core zones') + plt.legend() + plt.grid(True) + plt.savefig('plots/d_eps_nuc_dT.pdf',dpi = 300) + plt.show() + + + + "plot partials for d_eps_nuc_dRho in core zones" + plt.figure(figsize=(8, 6)) + for i in [1815,1816,1817,1818,1819]: + # Access data for a specific zone and iteration + zone_number = str(i) # Change as needed + zone_data = all_data[zone_number] + + # For example, get the eps_nuc values over iterations + iterations = [iter_data['Iteration'] for iter_data in zone_data] + eps_nuc_values = [iter_data['d_eps_nuc_dRho'] for iter_data in zone_data] + + # Plotting eps_nuc over iterations for the specified zone + + eps_nuc_values = np.array(eps_nuc_values) + # Initialize an array to store the log-scaled values + eps_nuc_values_log = np.zeros_like(eps_nuc_values) + + # Handle positive terms: directly take the logarithm + positive_mask = eps_nuc_values > 0 + eps_nuc_values_log[positive_mask] = np.log10(eps_nuc_values[positive_mask]) + + # Handle negative terms: take the logarithm of the absolute value and scale by -1 + negative_mask = eps_nuc_values < 0 + eps_nuc_values_log[negative_mask] = -np.log10(np.abs(eps_nuc_values[negative_mask])) + + #name = "zone" +str(i) + plt.plot(iterations,eps_nuc_values_log, marker='o', label = str(i)) + + +# plt.plot(iterations,eps_nuc_values_log, marker='o') +# plt.ylim(-1e21,1e21) + print (iterations) +# print (eps_nuc_values_log) + plt.xlabel('Iteration') + plt.ylabel('d_eps_nuc_dRho') + plt.title(f'd_eps_nuc_dRho over all iterations for core zones') + plt.legend() + plt.grid(True) + plt.savefig('plots/d_eps_nuc_dRho.pdf',dpi = 300) + plt.show() + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dRho.pdf b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dRho.pdf new file mode 100644 index 000000000..d7bee0868 Binary files /dev/null and b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dRho.pdf differ diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dT.pdf b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dT.pdf new file mode 100644 index 000000000..0525d2291 Binary files /dev/null and b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/d_eps_nuc_dT.pdf differ diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/eps_nuc.pdf b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/eps_nuc.pdf new file mode 100644 index 000000000..c99e711e9 Binary files /dev/null and b/star/dev_cases_test_solver/test_20M_near_cc_approx21/plots/eps_nuc.pdf differ diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/profile_columns.list b/star/dev_cases_test_solver/test_20M_near_cc_approx21/profile_columns.list new file mode 100644 index 000000000..7f6c0749c --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/re b/star/dev_cases_test_solver/test_20M_near_cc_approx21/re new file mode 100755 index 000000000..c9ef26f96 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/re_nomodfiles b/star/dev_cases_test_solver/test_20M_near_cc_approx21/re_nomodfiles new file mode 100755 index 000000000..565b74208 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/re_nomodfiles @@ -0,0 +1,41 @@ +#!/bin/bash + +#echo $# +#echo $1 +#echo $2 + +photo_directory=photos + +function most_recent_photo { + ls -t "$photo_directory" | head -1 +} + +if [ "$#" -ne 2 ] +then + echo "must pass two arguments, photo string and inlist name" + exit 1 +fi + +if [ $1 = "." ] +then + photo=$(most_recent_photo) +else + photo=$1 +fi + +if [ -z "$photo" ] || ! [ -f "$photo_directory/$photo" ] +then + echo "specified photo does not exist" + exit 1 +fi + +echo "restart from $photo" +if ! cp "$photo_directory/$photo" restart_photo +then + echo "failed to copy photo" + exit 1 +fi + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star $2 +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn new file mode 100755 index 000000000..f8a1b66f4 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn @@ -0,0 +1,56 @@ +#!/bin/bash + +# 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" + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_late_pre_zams.mod late_pre_zams.mod +# else +# do_one inlist_make_late_pre_zams_header late_pre_zams.mod +# cp late_pre_zams.mod standard_late_pre_zams.mod +# fi + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_zams.mod zams.mod +# else +# do_one inlist_to_zams_header zams.mod +# cp zams.mod standard_zams.mod +# fi + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_after_core_he_burn.mod after_core_he_burn.mod +# else +# do_one inlist_to_end_core_he_burn_header after_core_he_burn.mod +# cp after_core_he_burn.mod standard_after_core_he_burn.mod +# fi + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_removed_envelope.mod removed_envelope.mod +# else +# do_one inlist_remove_envelope_header removed_envelope.mod +# cp removed_envelope.mod standard_removed_envelope.mod +# fi + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then +# cp standard_after_core_c_burn.mod after_core_c_burn.mod +# else +# do_one inlist_to_end_core_c_burn_header after_core_c_burn.mod +# cp after_core_c_burn.mod standard_after_core_c_burn.mod +# fi + +# if [ -n "$MESA_SKIP_OPTIONAL" ]; then + cp standard_lgTmax.mod lgTmax.mod +# else +# do_one inlist_to_lgTmax_header lgTmax.mod +# cp lgTmax.mod standard_lgTmax.mod +# fi + +do_one inlist_to_cc_header final.mod + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +echo 'finished' + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn1 b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn1 new file mode 100755 index 000000000..25590040a --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/rn_all b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_all new file mode 100755 index 000000000..1ef23b587 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_all @@ -0,0 +1,31 @@ +#!/bin/bash + +# 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" + +do_one inlist_make_late_pre_zams_header late_pre_zams.mod +cp late_pre_zams.mod standard_late_pre_zams.mod + +do_one inlist_to_zams_header zams.mod +cp zams.mod standard_zams.mod + +do_one inlist_to_end_core_he_burn_header after_core_he_burn.mod +cp after_core_he_burn.mod standard_after_core_he_burn.mod + +do_one inlist_remove_envelope_header removed_envelope.mod +cp removed_envelope.mod standard_removed_envelope.mod + +do_one inlist_to_end_core_c_burn_header after_core_c_burn.mod +cp after_core_c_burn.mod standard_after_core_c_burn.mod + +do_one inlist_to_lgTmax_header lgTmax.mod +cp lgTmax.mod standard_lgTmax.mod + +do_one inlist_to_cc_header final.mod + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +echo 'finished' diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_nomodfiles b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_nomodfiles new file mode 100755 index 000000000..35bd84f4c --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_nomodfiles @@ -0,0 +1,7 @@ +#!/bin/bash +rm -f restart_photo +echo $1 + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star $1 +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_standard b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_standard new file mode 100755 index 000000000..fde5dbfab --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/rn_standard @@ -0,0 +1,56 @@ +#!/bin/bash + +# 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" + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_late_pre_zams.mod late_pre_zams.mod +else + do_one inlist_make_late_pre_zams_header late_pre_zams.mod + cp late_pre_zams.mod standard_late_pre_zams.mod +fi + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_zams.mod zams.mod +else + do_one inlist_to_zams_header zams.mod + cp zams.mod standard_zams.mod +fi + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_after_core_he_burn.mod after_core_he_burn.mod +else + do_one inlist_to_end_core_he_burn_header after_core_he_burn.mod + cp after_core_he_burn.mod standard_after_core_he_burn.mod +fi + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_removed_envelope.mod removed_envelope.mod +else + do_one inlist_remove_envelope_header removed_envelope.mod + cp removed_envelope.mod standard_removed_envelope.mod +fi + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_after_core_c_burn.mod after_core_c_burn.mod +else + do_one inlist_to_end_core_c_burn_header after_core_c_burn.mod + cp after_core_c_burn.mod standard_after_core_c_burn.mod +fi + +if [ -z "$MESA_RUN_OPTIONAL" ]; then + cp standard_lgTmax.mod lgTmax.mod +else + do_one inlist_to_lgTmax_header lgTmax.mod + cp lgTmax.mod standard_lgTmax.mod +fi + +do_one inlist_to_cc_header final.mod + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +echo 'finished' + diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/run_all b/star/dev_cases_test_solver/test_20M_near_cc_approx21/run_all new file mode 100755 index 000000000..e0b9ab8d7 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/run_all @@ -0,0 +1,6 @@ +./rn_nomodfiles inlist_make_late_pre_zams_header +./re_nomodfiles . inlist_to_zams_header +./re_nomodfiles . inlist_to_end_core_he_burn_header +./re_nomodfiles . inlist_to_end_core_c_burn_header +./re_nomodfiles . inlist_to_lgTmax_header +./re_nomodfiles . inlist_to_cc_header diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/src/run.f90 b/star/dev_cases_test_solver/test_20M_near_cc_approx21/src/run.f90 new file mode 100644 index 000000000..815cbeabd --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/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_solver/test_20M_near_cc_approx21/src/run_star_extras.f90 b/star/dev_cases_test_solver/test_20M_near_cc_approx21/src/run_star_extras.f90 new file mode 100644 index 000000000..2cf66b188 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/src/run_star_extras.f90 @@ -0,0 +1,227 @@ +! *********************************************************************** +! +! Copyright (C) 2010 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 + + implicit none + + include "test_suite_extras_def.inc" + +! here are the x controls used below + +!alpha_mlt_routine + !alpha_H = s% x_ctrl(21) + !alpha_other = s% x_ctrl(22) + !H_limit = s% x_ctrl(23) + + 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_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 + s% other_alpha_mlt => alpha_mlt_routine + end subroutine extras_controls + + + subroutine alpha_mlt_routine(id, ierr) + use chem_def, only: ih1 + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + integer :: k, h1 + real(dp) :: alpha_H, alpha_other, H_limit + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + alpha_H = s% x_ctrl(21) + alpha_other = s% x_ctrl(22) + H_limit = s% x_ctrl(23) + h1 = s% net_iso(ih1) + !write(*,1) 'alpha_H', alpha_H + !write(*,1) 'alpha_other', alpha_other + !write(*,1) 'H_limit', H_limit + !write(*,2) 'h1', h1 + !write(*,2) 's% nz', s% nz + if (alpha_H <= 0 .or. alpha_other <= 0 .or. h1 <= 0) return + do k=1,s% nz + if (s% xa(h1,k) >= H_limit) then + s% alpha_mlt(k) = alpha_H + else + s% alpha_mlt(k) = alpha_other + end if + !write(*,2) 'alpha_mlt', k, s% alpha_mlt(k), + end do + !stop + end subroutine alpha_mlt_routine + + + subroutine extras_startup(id, restart, ierr) + integer, intent(in) :: id + logical, intent(in) :: restart + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + call test_suite_startup(s, restart, ierr) + + if (.not. s% x_logical_ctrl(37)) return + + end subroutine extras_startup + + + subroutine extras_after_evolve(id, ierr) + use num_lib, only: find0 + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + real(dp) :: dt, m + integer :: k, nz + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + nz = s% nz + write(*,'(A)') + select case (s% x_integer_ctrl(5)) + case (7) + ! put target info in TestHub output + testhub_extras_names(1) = 'fe_core_mass'; testhub_extras_vals(1) = s% fe_core_mass + + if(s% fe_core_mass < 1d0) then + write(*,1) "Bad fe_core_mass", s%fe_core_mass + else + if(s% fe_core_infall > s% fe_core_infall_limit) then + write(*,'(a)') 'all values are within tolerance' + else + write(*,'(a)') "Bad fe core infall" + end if + end if + end select + call test_suite_after_evolve(s, ierr) + if (.not. s% x_logical_ctrl(37)) return + end subroutine extras_after_evolve + + + ! returns either keep_going, retry, or terminate. + integer function extras_check_model(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + 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 = 0 + 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) + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + 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 + how_many_extra_profile_columns = 1 + 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) + integer, intent(out) :: ierr + type (star_info), pointer :: s + integer :: k + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + names(1) = 'zbar_div_abar' + do k=1,s% nz + vals(k,1) = s% zbar(k)/s% abar(k) + end do + end subroutine data_for_extra_profile_columns + + ! returns either keep_going or terminate. + integer function extras_finish_step(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_finish_step = keep_going + if (.not. s% x_logical_ctrl(37)) return + 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_solver/test_20M_near_cc_approx21/standard_after_core_c_burn.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_after_core_c_burn.mod new file mode 100644 index 000000000..27d843893 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_after_core_c_burn.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0438de9dfa9479211fd27de44844fb5d5e64a2ee279369d95b183034d3912f2f +size 785401 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_after_core_he_burn.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_after_core_he_burn.mod new file mode 100644 index 000000000..b068bd47e --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_after_core_he_burn.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cad819ff0177ba09b07f78a3a5788d4e04627701c186a96c18d3bfbf5925da59 +size 700871 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_late_pre_zams.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_late_pre_zams.mod new file mode 100644 index 000000000..e6246cfc8 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_late_pre_zams.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e2da28aa6d45c618355a983ffab589b5ab5cf12b6b03a4c6de60fac8bb384631 +size 198408 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_lgTmax.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_lgTmax.mod new file mode 100644 index 000000000..543e1555d --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_lgTmax.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:808de00ef2b333e36afbbf0a8dda7ebc146e78560f8faf79a8b12ce6288e4acf +size 1428461 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_removed_envelope.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_removed_envelope.mod new file mode 100644 index 000000000..e2e3b1d99 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_removed_envelope.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:636aa0a7b04e6809ddd105569d27d629451f4a3cefb788b56e9908ccc0aece3c +size 700871 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_zams.mod b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_zams.mod new file mode 100644 index 000000000..9a4ba47f3 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/standard_zams.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:133a5ee96544e548f4094310fe94ef2dbce7246c8bdbe51f9a6a7fb16ce119b0 +size 309031 diff --git a/star/dev_cases_test_solver/test_20M_near_cc_approx21/test b/star/dev_cases_test_solver/test_20M_near_cc_approx21/test new file mode 100755 index 000000000..abe072d53 --- /dev/null +++ b/star/dev_cases_test_solver/test_20M_near_cc_approx21/test @@ -0,0 +1,2 @@ +#./re_nomodfiles x00002100 inlist_to_cc_header +./re_nomodfiles x00002131 inlist_to_cc_header diff --git a/star/private/net.f90 b/star/private/net.f90 index 943f01859..e39f1c933 100644 --- a/star/private/net.f90 +++ b/star/private/net.f90 @@ -154,6 +154,10 @@ subroutine do1_net(s, k, species, & type (Net_Info) :: n real(dp), pointer :: reaction_neuQs(:) logical :: clipped_T + ! Declare variables for writing net data to file. + character(len=20) :: filename + integer :: unit_number, num_isotopes, ii, jj + logical, parameter :: dbg = .false. @@ -257,6 +261,66 @@ subroutine do1_net(s, k, species, & s% dxdt_nuc(:,k), s% d_dxdt_nuc_dRho(:,k), s% d_dxdt_nuc_dT(:,k), s% d_dxdt_nuc_dx(:,:,k), & screening_mode, s% eps_nuc_categories(:,k), & s% eps_nuc_neu_total(k), ierr) + + +! if (k == 1815 .or. k == 1816 .or. k == 1817 .or. k == 1818 .or. k == 1819) then +! write(*,*) 'XXX CELL ', k +! write(*,*) 'eps_nuc', s% eps_nuc(k) +! !write(*,*) 'eps_nuc', s% eps_nuc(k), s% rho(k), T +! !write(*,*) 'eps_nuc', s% eps_nuc(k), s% d_epsnuc_dx(:,k) +! !write(*,*) 'eps_nuc', s% eps_nuc(k), s% abar(k), s% zbar(k), s% z2bar(k), s% ye(k) +! !write(*,*) 'eps_nuc', s% eps_nuc(k), d_eps_nuc_dRho, d_eps_nuc_dT +! !!write(*,*) 'd_eps_nuc_dRho', d_eps_nuc_dRho +! !!write(*,*) 'd_eps_nuc_dT', d_eps_nuc_dT +! !write(*,*) 'd_epsnuc_dx', s% d_epsnuc_dx(:,k) +! !write(*,*) 'dxdt_nuc', s% dxdt_nuc(:,k) +! !write(*,*) 'd_dxdt_nuc_dRho', s% d_dxdt_nuc_dRho(:,k) +! !write(*,*) 'd_dxdt_nuc_dT', s% d_dxdt_nuc_dT(:,k) +! !write(*,*) 'd_dxdt_nuc_dx', s% d_dxdt_nuc_dx(:,:,k) +! end if + + + if (k == 1815 .or. k == 1816 .or. k == 1817 .or. k == 1818 .or. k == 1819) then + ! Construct the filename based on the zone + write(filename, '("zone_",I0,".txt")') k + ! Assign a unique unit number for file I/O + unit_number = 10 + k - 1815 ! Unit numbers from 10 upwards + ! Open the file in append mode + open(unit=unit_number, file=filename, status='unknown', position='append') + ! Write the iteration marker + write(unit_number, '(A,I0)') 'Iteration: ', s%num_solver_iterations +! write(unit_number, '(A,I0)') 'solver_adjust_iter: ', s% solver_adjust_iter + ! Write scalar variables with modified format + write(unit_number, '(A,ES23.15E3)') 'eps_nuc: ', s%eps_nuc(k) + write(unit_number, '(A,ES23.15E3)') 'rho: ', s%rho(k) + write(unit_number, '(A,ES23.15E3)') 'T: ', T + write(unit_number, '(A,ES23.15E3)') 'd_eps_nuc_dRho: ', d_eps_nuc_dRho + write(unit_number, '(A,ES23.15E3)') 'd_eps_nuc_dT: ', d_eps_nuc_dT + ! Determine the number of isotopes + num_isotopes = size(s%d_epsnuc_dx,1) + ! Write vector variables with modified format + write(unit_number, '(A)') 'd_epsnuc_dx: ' + write(unit_number, '(ES23.15E3,1X)') (s%d_epsnuc_dx(ii,k), ii=1,num_isotopes) + write(unit_number, *) + write(unit_number, '(A)') 'dxdt_nuc: ' + write(unit_number, '(ES23.15E3,1X)') (s%dxdt_nuc(ii,k), ii=1,num_isotopes) + write(unit_number, *) + write(unit_number, '(A)') 'd_dxdt_nuc_dRho: ' + write(unit_number, '(ES23.15E3,1X)') (s%d_dxdt_nuc_dRho(ii,k), ii=1,num_isotopes) + write(unit_number, *) + write(unit_number, '(A)') 'd_dxdt_nuc_dT: ' + write(unit_number, '(ES23.15E3,1X)') (s%d_dxdt_nuc_dT(ii,k), ii=1,num_isotopes) + write(unit_number, *) + ! Write the 2D array variable + write(unit_number, '(A)') 'd_dxdt_nuc_dx: ' + do ii = 1, num_isotopes + write(unit_number, '(ES23.15E3,1X)') (s%d_dxdt_nuc_dx(ii,jj,k), jj=1,num_isotopes) + end do + write(unit_number, *) ! New line + ! Close the file + close(unit_number) + end if + end if