diff --git a/star/Makefile b/star/Makefile index c248e3d6b..33fb44b3c 100644 --- a/star/Makefile +++ b/star/Makefile @@ -98,6 +98,7 @@ SRCS := \ private/init_model.f90 \ private/ionization_potentials.f90 \ private/kap_support.f90 \ + private/mlt_tdc_face_support.f90 \ private/magnetic_diffusion.f90 \ private/mass_utils.f90 \ private/mesh_adjust.f90 \ diff --git a/star/defaults/controls_dev.defaults b/star/defaults/controls_dev.defaults index 59666ad67..e578abdfe 100644 --- a/star/defaults/controls_dev.defaults +++ b/star/defaults/controls_dev.defaults @@ -200,6 +200,25 @@ use_TDC_enthalpy_flux_limiter = .false. + ! use_face_values_eos_and_kap_mlt_tdc + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! MESA interpolates primitive solver variables such as T, rho, L, and + ! composition, together with constructed quantities such as opacity, + ! pressure, and other EOS results, onto cell faces for use in MLT and + ! TDC. If use_face_values_eos_and_kap_mlt_tdc = .true., this option + ! instead reconstructs only the primitive face variables and then makes + ! additional EOS and opacity calls to obtain the thermodynamic quantities + ! corresponding to that face state. Despite the cost of the extra EOS and + ! opacity work, this can improve accuracy on coarse meshes, especially in + ! unresolved opacity bumps in pulsating envelopes, and is a potential + ! candidate for becoming the default mode, pending further testing. + + ! :: + + use_face_values_eos_and_kap_mlt_tdc = .false. + + ! include_mlt_Pturb_in_thermodynamic_gradients ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -218,6 +237,7 @@ ! If doing velocity time_centering, include time centering in mlt. Depending on the time ! centering flags adopted, include_mlt_in_velocity_time_centering includes time centering ! for P, L, and r in geff, gradr, grada, and mlt. + ! Currently not supported together with ``use_face_values_eos_and_kap_mlt_tdc``. ! :: diff --git a/star/defaults/profile_columns.list b/star/defaults/profile_columns.list index 2bf7338c2..188d2e45e 100644 --- a/star/defaults/profile_columns.list +++ b/star/defaults/profile_columns.list @@ -740,6 +740,17 @@ !Uq !Y_face + !mlt_tdc_T_face + !mlt_tdc_rho_face + !mlt_tdc_P_face + !mlt_tdc_Cp_face + !mlt_tdc_ChiRho_face + !mlt_tdc_ChiT_face + !mlt_tdc_grada_face + !mlt_tdc_opacity_face + !mlt_tdc_scale_height_face + !mlt_tdc_gradr_face + !# RTI !RTI_du_diffusion_kick diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/Makefile b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/Makefile new file mode 100644 index 000000000..c779b9d6d --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/Makefile @@ -0,0 +1 @@ +include $(MESA_DIR)/star/work/Makefile \ No newline at end of file diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/ck b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/ck new file mode 100755 index 000000000..78ca63a8c --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/ck @@ -0,0 +1,7 @@ +#!/usr/bin/env 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_TDC_Pulsation/dev_TDC_Cepheid_9M/clean b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/clean new file mode 100755 index 000000000..9dff1275b --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/clean @@ -0,0 +1,4 @@ +#!/usr/bin/env bash + +cd make +make clean diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/gyre.in b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/gyre.in new file mode 100644 index 000000000..10381b567 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/gyre.in @@ -0,0 +1,34 @@ +&model + add_center = .FALSE. +/ + +&mode + l = 0 +/ + +&osc + inner_bound = 'ZERO_R' ! for envelopes + nonadiabatic = .TRUE. +/ + +&rot +/ + +&num + diff_scheme = 'MAGNUS_GL2'!'COLLOC_GL2' +/ + +&scan + grid_type = 'LINEAR' + freq_min = 0.5!1d-4 + freq_max = 5.0!10.0 + freq_min_units = 'ACOUSTIC_DELTA' + freq_max_units = 'ACOUSTIC_DELTA' + n_freq = 50 +/ + +&grid + w_osc = 10 + w_exp = 2 + w_ctr = 10 +/ diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/history_columns.list b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/history_columns.list new file mode 100644 index 000000000..f461f8f35 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/history_columns.list @@ -0,0 +1,1073 @@ +! history_columns.list -- determines the contents of star history logs +! you can use a non-standard version by setting history_columns_file in your inlist + +! units are cgs unless otherwise noted. + +! reorder the following names as desired to reorder columns. +! comment out the name to omit a column (fewer columns => less IO => faster running). +! remove '!' to restore a column. + +! if you have a situation where you want a non-standard set of columns, +! make a copy of this file, edit as desired, and give the new filename in your inlist +! as history_columns_file. if you are just adding columns, you can 'include' this file, +! and just list the additions in your file. note: to include the standard default +! version, use include '' -- the 0 length string means include the default file. + +! blank lines and comments can be used freely. +! if a column name appears more than once in the list, only the first occurrence is used. + +! if you need to have something added to the list of options, let me know.... + + +! the first few lines of the log file contain a few items: + + ! version_number -- for the version of mesa being used + ! burn_min1 -- 1st limit for reported burning, in erg/g/s + ! burn_min2 -- 2nd limit for reported burning, in erg/g/s + + +!# other files + +! note: you can include another list by doing +! include 'filename' +! include '' means include the default standard list file + +! the following lines of the log file contain info about 1 model per row + +!---------------------------------------------------------------------------------------------- + +!# general info about the model + + model_number ! counting from the start of the run + num_zones ! number of zones in the model + + !## age + + star_age ! elapsed simulated time in years since the start of the run + !star_age_sec ! elapsed simulated time in seconds since the start of the run + !star_age_min ! elapsed simulated time in minutes since the start of the run + !star_age_hr ! elapsed simulated time in hours since the start of the run + !star_age_day ! elapsed simulated time in days since the start of the run + day ! elapsed simulated time in days since the start of the run + + log_star_age + log_star_age_sec + + !## timestep + + time_step ! timestep in years since previous model + time_step_sec ! timestep in seconds since previous model + time_step_days + log_dt ! log10 time_step in years + log_dt_sec ! log10 time_step in seconds + !log_dt_days ! log10 time_step in days + + !## mass + + star_mass ! in Msun units + !log_star_mass + + !star_gravitational_mass ! star_mass is baryonic mass + !star_mass_grav_div_mass + + !delta_mass ! star_mass - initial_mass in Msun units + log_xmstar ! log10 mass exterior to M_center (grams) + + !## mass change + + star_mdot ! d(star_mass)/dt (in msolar per year) + log_abs_mdot ! log10(abs(star_mdot)) (in msolar per year) + + !## imposed surface conditions + !Tsurf_factor + !tau_factor + !tau_surface + + !## imposed center conditions + !m_center + !m_center_gm + !r_center + !r_center_cm + !r_center_km + !L_center + !log_L_center + !log_L_center_ergs_s + !v_center + !v_center_kms + + !logt_max + +!---------------------------------------------------------------------------------------------- + +!# mixing and convection + + !max_conv_vel_div_csound + !max_gradT_div_grada + !max_gradT_sub_grada + !min_log_mlt_Gamma + + + !## mixing regions + + mass_conv_core ! (Msun) mass coord of top of convective core. 0 if core is not convective + + ! mx1 refers to the largest (by mass) convective region. + ! mx2 is the 2nd largest. + + ! conv_mx1_top and conv_mx1_bot are the region where mixing_type == convective_mixing. + ! mx1_top and mx1_bot are the extent of all kinds of mixing, convective and other. + + ! values are m/Mstar + conv_mx1_top + conv_mx1_bot + conv_mx2_top + conv_mx2_bot + mx1_top + mx1_bot + mx2_top + mx2_bot + + ! radius -- values are radii in Rsun units + !conv_mx1_top_r + !conv_mx1_bot_r + !conv_mx2_top_r + !conv_mx2_bot_r + !mx1_top_r + !mx1_bot_r + !mx2_top_r + !mx2_bot_r + + ! you might want to get a more complete list of mixing regions by using the following + + !mixing_regions ! note: this includes regions where the mixing type is no_mixing. + + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives the mixing type as defined in const/public/const_def.f90. + + ! the second column for a region gives the m/mstar location of the top of the region + ! entries for extra columns after the last region in the star will have an invalid mixing_type value of -1. + ! mstar is the total mass of the star, so these locations range from 0 to 1 + ! all regions are include starting from the center, so the bottom of one region + ! is the top of the previous one. since we start at the center, the bottom of the 1st region is 0. + + ! the columns in the log file will have names like 'mix_type_1' and 'mix_qtop_1' + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + + + !mix_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'mix_relr_type_1' and 'mix_relr_top_1' + + + !## conditions at base of largest convection zone (by mass) + !cz_bot_mass ! mass coordinate of base (Msun) + !cz_mass ! mass coordinate of base (Msun) -- same as cz_bot_mass + !cz_log_xmass ! mass exterior to base (g) + !cz_log_xmsun ! mass exterior to base (Msun) + !cz_xm ! mass exterior to base (Msun) + !cz_logT + !cz_logRho + !cz_logP + !cz_bot_radius ! Rsun + !cz_log_column_depth + !cz_log_radial_depth + !cz_luminosity ! Lsun + !cz_opacity + !cz_log_tau + !cz_eta + !cz_log_eps_nuc ! log10(ergs/g/s) + !cz_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_csound + !cz_scale_height + !cz_grav + + !cz_omega + !cz_omega_div_omega_crit + + !cz_zone + + ! mass fractions at base of largest convection zone (by mass) + !cz_log_xa h1 + !cz_log_xa he4 + + !## conditions at top of largest convection zone (by mass) + !cz_top_mass ! mass coordinate of top (Msun) + !cz_top_log_xmass ! mass exterior to top (g) + !cz_top_log_xmsun ! mass exterior to top (Msun) + !cz_top_xm ! mass exterior to top (Msun) + !cz_top_logT + !cz_top_logRho + !cz_top_logP + !cz_top_radius ! Rsun + !cz_top_log_column_depth + !cz_top_log_radial_depth + !cz_top_luminosity ! Lsun + !cz_top_opacity + !cz_top_log_tau + !cz_top_eta + !cz_top_log_eps_nuc ! log10(ergs/g/s) + !cz_top_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_top_csound + !cz_top_scale_height + !cz_top_grav + + !cz_top_omega + !cz_top_omega_div_omega_crit + + !cz_top_zone + !cz_top_zone_logdq + + ! mass fractions at top of largest convection zone (by mass) + !cz_top_log_xa h1 + !cz_top_log_xa he4 + +!---------------------------------------------------------------------------------------------- + +!# nuclear reactions + + !## integrated quantities + + !power_h_burn ! total thermal power from PP and CNO, excluding neutrinos (in Lsun units) + !power_he_burn ! total thermal power from triple-alpha, excluding neutrinos (in Lsun units) + !power_photo + !power_z_burn + !log_power_nuc_burn ! total thermal power from all burning, excluding photodisintegrations + log_LH ! log10 power_h_burn + log_LHe ! log10 power_he_burn + log_LZ ! log10 total burning power including LC, but excluding LH and LHe and photodisintegrations + log_Lnuc ! log(LH + LHe + LZ) + !log_Lnuc_ergs_s + !log_Lnuc_sub_log_L + !lnuc_photo + + !extra_L ! integral of extra_heat in Lsun units + !log_extra_L ! log10 extra_L + + !## neutrino losses + log_Lneu ! log10 power emitted in neutrinos, nuclear and thermal (in Lsun units) + !log_Lneu_nuc ! log10 power emitted in neutrinos, nuclear sources only (in Lsun units) + !log_Lneu_nonnuc ! log10 power emitted in neutrinos, thermal sources only (in Lsun units) + + !mass_loc_of_max_eps_nuc ! (in Msun units) + !mass_ext_to_max_eps_nuc ! (in Msun units) + !eps_grav_integral ! (in Lsun units) + !log_abs_Lgrav ! log10 abs(eps_grav_integral) (in Lsun units) + + !## information about reactions (by category) + + ! log10 total luminosity for reaction categories (Lsun units) + + pp + cno + tri_alpha + !c_alpha + !n_alpha + !o_alpha + !ne_alpha + !na_alpha + !mg_alpha + !si_alpha + !s_alpha + !ar_alpha + !ca_alpha + !ti_alpha + !fe_co_ni + !c12_c12 + !c12_o16 + !o16_o16 + !photo + !pnhe4 + !other + + !## information about individual reactions + + ! adds columns for all of the reactions that are in the current net + ! Note that if using op_split_burn=.true. then zones which have been split will report 0 for thier rates + !add_raw_rates ! raw reaction rates, reactions/second + !add_screened_rates ! screened reaction rates reactions/second + !add_eps_nuc_rates ! Nuclear energy (minus neutrino losses) released erg/s + !add_eps_neu_rates ! Neutrino losses erg/s + + ! individual reactions (as many as desired) + ! use list_net_reactions = .true. in star_job to list all reactions in the current net + ! reactions/second + !raw_rate r_h1_h1_ec_h2 + !raw_rate r_h1_h1_wk_h2 + + + + !## nuclear reactions at center + + ! center log10 burn erg/g/s for reaction categories + + !c_log_eps_burn cno + !c_log_eps_burn tri_alfa + + ! center d_eps_nuc_dlnd for reaction categories + + !c_d_eps_dlnd cno + !c_d_eps_dlnd tri_alfa + + ! center d_eps_nuc_dlnT for reaction categories + + !c_d_eps_dlnT cno + !c_d_eps_dlnT tri_alfa + + !## regions of strong nuclear burning + + ! 2 zones where eps_nuc > burn_min1 erg/g/s + ! for each zone have 4 numbers: start1, start2, end2, end1 + ! start1 is mass of inner edge where first goes > burn_min1 (or -20 if none such) + ! start2 is mass of inner edge where first zone reaches burn_min2 erg/g/sec (or -20 if none such) + ! end2 is mass of outer edge where first zone drops back below burn_min2 erg/g/s + ! end1 is mass of outer edge where first zone ends (i.e. eps_nuc < burn_min1) + ! similar for the second zone + + epsnuc_M_1 ! start1 for 1st zone + epsnuc_M_2 ! start2 + epsnuc_M_3 ! end2 + epsnuc_M_4 ! end1 + + epsnuc_M_5 ! start1 for 2nd zone + epsnuc_M_6 ! start2 + epsnuc_M_7 ! end2 + epsnuc_M_8 ! end1 + + + ! you might want to get a more complete list of burning regions by using the following + + !burning_regions + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives int(sign(val)*log10(max(1,abs(val)))) + ! where val = ergs/gm/sec nuclear energy minus all neutrino losses. + ! the second column for a region gives the q location of the top of the region + ! entries for extra columns after the last region in the star will have a value of -9999 + ! all regions are included starting from the center, so the bottom of one region + ! is the top of the previous one. + ! since we start at the center, the bottom of the 1st region is q=0 and top of last is q=1. + + ! the columns in the log file will have names like 'burn_type_1' and 'burn_qtop_1' + + !burn_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'burn_relr_type_1' and 'burn_relr_top_1' + + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + +!---------------------------------------------------------------------------------------------- + +!# information about core and envelope + + !## helium core + he_core_mass + he_core_radius + !he_core_lgT + !he_core_lgRho + !he_core_L + !he_core_v + !he_core_omega + !he_core_omega_div_omega_crit + !he_core_k + + !## CO core + co_core_mass + !CO_core + co_core_radius + !co_core_lgT + !co_core_lgRho + !co_core_L + !co_core_v + !co_core_omega + !co_core_omega_div_omega_crit + !co_core_k + + !## ONe core + one_core_mass + !one_core_radius + !one_core_lgT + !one_core_lgRho + !one_core_L + !one_core_v + !one_core_omega + !one_core_omega_div_omega_crit + !one_core_k + + !## iron core + fe_core_mass + !fe_core_radius + !fe_core_lgT + !fe_core_lgRho + !fe_core_L + !fe_core_v + !fe_core_omega + !fe_core_omega_div_omega_crit + !fe_core_k + + !## neutron rich core + neutron_rich_core_mass + !neutron_rich_core_radius + !neutron_rich_core_lgT + !neutron_rich_core_lgRho + !neutron_rich_core_L + !neutron_rich_core_v + !neutron_rich_core_omega + !neutron_rich_core_omega_div_omega_crit + !neutron_rich_core_k + + !## envelope + + !envelope_mass ! = star_mass - he_core_mass + !envelope_fraction_left ! = envelope_mass / (initial_mass - he_core_mass) + + !h_rich_layer_mass ! = star_mass - he_core_mass + !he_rich_layer_mass ! = he_core_mass - c_core_mass + !co_rich_layer_mass + +!---------------------------------------------------------------------------------------------- + +!# timescales + + !dynamic_timescale ! dynamic timescale (seconds) -- estimated by 2*pi*sqrt(r^3/(G*m)) + !kh_timescale ! kelvin-helmholtz timescale (years) + !mdot_timescale ! star_mass/abs(star_mdot) (years) + !kh_div_mdot_timescales ! kh_timescale/mdot_timescale + !nuc_timescale ! nuclear timescale (years) -- proportional to mass divided by luminosity + + !dt_cell_collapse ! min time for any cell to collapse at current velocities + !dt_div_dt_cell_collapse + + !dt_div_max_tau_conv ! dt/ maximum conv timescale + !dt_div_min_tau_conv ! dt/ minimum conv timescale + + + !min_dr_div_cs ! min over all cells of dr/csound (seconds) + !min_dr_div_cs_k ! location of min + !log_min_dr_div_cs ! log10 min dr_div_csound (seconds) + !min_dr_div_cs_yr ! min over all cells of dr/csound (years) + !log_min_dr_div_cs_yr ! log10 min dr_div_csound (years) + !dt_div_min_dr_div_cs + !log_dt_div_min_dr_div_cs + + !min_t_eddy ! minimum value of scale_height/conv_velocity + +!---------------------------------------------------------------------------------------------- + +!# conditions at or near the surface of the model + + !## conditions at the photosphere + !effective_T + !Teff + log_Teff ! log10 effective temperature + ! Teff is calculated using Stefan-Boltzmann relation L = 4 pi R^2 sigma Teff^4, + ! where L and R are evaluated at the photosphere (tau_factor < 1) + ! or surface of the model (tau_factor >= 1) when photosphere is not inside the model. + + !photosphere_black_body_T + !photosphere_cell_T ! temperature at model location closest to the photosphere, not necessarily Teff + !photosphere_cell_log_T + !photosphere_cell_density + !photosphere_cell_log_density + !photosphere_cell_opacity + !photosphere_cell_log_opacity + !photosphere_L ! Lsun units + !photosphere_log_L ! Lsun units + !photosphere_r ! Rsun units + !photosphere_log_r ! Rsun units + !photosphere_m ! Msun units + !photosphere_v_km_s + !photosphere_cell_k + !photosphere_column_density + !photosphere_csound + !photosphere_log_column_density + !photosphere_opacity + !photosphere_v_div_cs + !photosphere_xm + !photosphere_cell_free_e + !photosphere_cell_log_free_e + !photosphere_logg + !photosphere_T + + !## conditions at or near the surface of the model (outer edge of outer cell) + + luminosity ! luminosity in Lsun units + !luminosity_ergs_s ! luminosity in cgs units + log_L ! log10 luminosity in Lsun units + !log_L_ergs_s ! log10 luminosity in cgs units + radius ! Rsun + log_R ! log10 radius in Rsun units + !radius_cm + !log_R_cm + + log_g ! log10 gravity + gravity + log_Ledd + log_L_div_Ledd ! log10(L/Leddington) + lum_div_Ledd + !log_surf_optical_depth + !surface_optical_depth + + !log_surf_cell_opacity ! old name was log_surf_opacity + !log_surf_cell_P ! old name was log_surf_P + !log_surf_cell_pressure ! old name was log_surf_pressure + !log_surf_cell_density ! old name was log_surf_density + !log_surf_cell_temperature ! old name was log_surf_temperature + !surface_cell_temperature ! old name was surface_temperature + !log_surf_cell_z ! old name was log_surf_z + !surface_cell_entropy ! in units of kerg per baryon + ! old name was surface_entropy + + v_surf ! (cm/s) + v_surf_km_s ! (km/s) + v_div_csound_surf ! velocity divided by sound speed at outermost grid point + !v_div_csound_max ! max value of velocity divided by sound speed at face + !v_div_vesc + !v_phot_km_s + !v_surf_div_escape_v + + !v_surf_div_v_kh ! v_surf/(photosphere_r/kh_timescale) + + !surf_avg_j_rot + !surf_avg_omega + !surf_avg_omega_crit + !surf_avg_omega_div_omega_crit + !surf_avg_v_rot ! km/sec rotational velocity at equator + !surf_avg_v_crit ! critical rotational velocity at equator + !surf_avg_v_div_v_crit + !surf_avg_Lrad_div_Ledd + !surf_avg_logT + !surf_avg_logRho + !surf_avg_opacity + + ! Gravity Darkening, reports the surface averaged L/Lsun and Teff (K) caused by + ! gravity darkening in rotating stars. Based on the model of Espinosa Lara & Rieutord (2011) + ! 'polar' refers to the line of sight being directed along the rotation axis of the star + ! 'equatorial' refers to the line of sight coincident with the stellar equator + !grav_dark_L_polar !Lsun + !grav_dark_Teff_polar !K + !grav_dark_L_equatorial !Lsun + !grav_dark_Teff_equatorial !K + + !surf_escape_v ! cm/s + + !v_wind_Km_per_s ! Km/s + ! = 1d-5*s% opacity(1)*max(0d0,-s% mstar_dot)/ & + ! (4*pi*s% photosphere_r*Rsun*s% tau_base) + ! Lars says: + ! wind_mdot = 4*pi*R^2*rho*v_wind + ! tau = integral(opacity*rho*dr) from R to infinity + ! so tau = opacity*wind_mdot/(4*pi*R*v_wind) at photosphere + ! or v_wind = opacity*wind_mdot/(4*pi*R*tau) at photosphere + + !rotational_mdot_boost ! factor for increase in mass loss mdot due to rotation + !log_rotational_mdot_boost ! log factor for increase in mass loss mdot due to rotation + !surf_r_equatorial_div_r_polar + !surf_r_equatorial_div_r + !surf_r_polar_div_r + +!---------------------------------------------------------------------------------------------- + +!# conditions near center + + log_center_T ! temperature + log_center_Rho ! density + log_center_P ! pressure + + ! shorter names for above + log_cntr_P + log_cntr_Rho + log_cntr_T + + !center_T ! temperature + !center_Rho ! density + !center_P ! pressure + + !center_degeneracy ! the electron chemical potential in units of k*T + !center_gamma ! plasma interaction parameter + center_mu + center_ye + center_abar + !center_zbar + + !center_eps_grav + + !center_non_nuc_neu + !center_eps_nuc + !d_center_eps_nuc_dlnT + !d_center_eps_nuc_dlnd + !log_center_eps_nuc + + !center_entropy ! in units of kerg per baryon + !max_entropy ! in units of kerg per baryon + !fe_core_infall + !non_fe_core_infall + !non_fe_core_rebound + !max_infall_speed + + !compactness_parameter ! (m/Msun)/(R(m)/1000km) for m = 2.5 Msun + !compactness + !m4 ! Mass co-ordinate where entropy=4 + ! mu4 is sensitive to the choice of how much dm/dr you average over, thus we average dm and dr over M(entropy=4) and M(entropy=4)+0.3Msun + !mu4 ! dM(Msun)/dr(1000km) where entropy=4 + + + !center_omega + !center_omega_div_omega_crit + +!---------------------------------------------------------------------------------------------- + +!# abundances + + !species ! size of net + + !## mass fractions near center + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_center_abundances + !add_log_center_abundances + + ! individual central mass fractions (as many as desired) + center h1 + center he4 + center c12 + center o16 + + ! individual log10 central mass fractions (as many as desired) + !log_center h1 + !log_center he4 + ! etc. + + + !## mass fractions near surface + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_surface_abundances + !add_log_surface_abundances + + ! individual surface mass fractions (as many as desired) + !surface h1 + !surface he4 + surface c12 + surface o16 + ! etc. + + ! individual log10 surface mass fractions (as many as desired) + + !log_surface h1 + !log_surface he4 + + + !## mass fractions for entire star + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_average_abundances + !add_log_average_abundances + + ! individual average mass fractions (as many as desired) + !average h1 + !average he4 + ! etc. + + ! individual log10 average mass fractions (as many as desired) + !log_average h1 + !log_average he4 + ! etc. + + + !## mass totals for entire star (in Msun units) + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_total_mass + !add_log_total_mass + + ! individual mass totals for entire star (as many as desired) + total_mass h1 + total_mass he4 + ! etc. + + ! individial log10 mass totals for entire star (in Msun units) + !log_total_mass h1 + !log_total_mass he4 + ! etc. + +!---------------------------------------------------------------------------------------------- + +!# info at specific locations + + !## info at location of max temperature + !max_T + !log_max_T + + +!---------------------------------------------------------------------------------------------- + +!# information about shocks + + !## info about outermost outward moving shock + ! excluding locations with q > max_q_for_outer_mach1_location + ! returns values at location of max velocity + !shock_mass ! baryonic (Msun) + !shock_mass_gm ! baryonic (grams) + !shock_q + !shock_radius ! (Rsun) + !shock_radius_cm ! (cm) + !shock_velocity + !shock_csound + !shock_v_div_cs + !shock_lgT + !shock_lgRho + !shock_lgP + !shock_gamma1 + !shock_entropy + !shock_tau + !shock_k + !shock_pre_lgRho + +!---------------------------------------------------------------------------------------------- + +!# asteroseismology + + delta_nu ! large frequency separation for p-modes (microHz) + ! 1e6/(seconds for sound to cross diameter of star) + delta_Pg ! g-mode period spacing for l=1 (seconds) + ! sqrt(2) pi^2/(integral of brunt_N/r dr) + !log_delta_Pg + nu_max ! estimate from scaling relation (microHz) + ! nu_max = nu_max_sun * M/Msun / ((R/Rsun)^2 (Teff/Teff_sun)^0.5) + ! with nu_max_sun = 3100 microHz, Teff_sun = 5777 + !nu_max_3_4th_div_delta_nu ! nu_max^0.75/delta_nu + acoustic_cutoff ! 0.5*g*sqrt(gamma1*rho/P) at surface + acoustic_radius ! integral of dr/csound (seconds) + !ng_for_nu_max ! = 1 / (nu_max*delta_Pg) + ! period for g-mode with frequency nu_max = nu_max_ng*delta_Pg + !gs_per_delta_nu ! delta_nu / (nu_max**2*delta_Pg) + ! number of g-modes per delta_nu at nu_max + + !int_k_r_dr_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=1 + !int_k_r_dr_2pt0_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=1 + !int_k_r_dr_0pt5_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=1 + !int_k_r_dr_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=2 + !int_k_r_dr_2pt0_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=2 + !int_k_r_dr_0pt5_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=2 + !int_k_r_dr_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=3 + !int_k_r_dr_2pt0_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=3 + !int_k_r_dr_0pt5_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=3 + +!---------------------------------------------------------------------------------------------- + +!# energy information + + !total_energy ! at end of step + !log_total_energy ! log(abs(total_energy)) + !total_energy_after_adjust_mass ! after mass adjustments + + ! shorter versions of above + !tot_E + !log_tot_E + + + !total_gravitational_energy + !log_total_gravitational_energy ! log(abs(total_gravitational_energy)) + !total_gravitational_energy_after_adjust_mass + + ! shorter versions of above + !tot_PE + !log_tot_PE + + !total_internal_energy + !log_total_internal_energy + !total_internal_energy_after_adjust_mass + + ! shorter versions of above + !tot_IE + !log_tot_IE + + !total_radial_kinetic_energy + !log_total_radial_kinetic_energy + !total_radial_kinetic_energy_after_adjust_mass + + ! shorter versions of above (does not include rot KE) + !tot_KE + !log_tot_KE + + !total_turbulent_energy + !log_total_turbulent_energy + !total_turbulent_energy_after_adjust_mass + !tot_Et + !log_tot_Et + + !total_energy_foe + + !tot_IE_div_IE_plus_KE + !total_IE_div_IE_plus_KE + + !total_entropy + !total_eps_grav + + !total_energy_sources_and_sinks ! for this step + !total_nuclear_heating + !total_non_nuc_neu_cooling + !total_irradiation_heating + !total_extra_heating ! extra heat integrated over the model times dt (erg) + !total_WD_sedimentation_heating + + !rel_run_E_err + + !rel_E_err + !abs_rel_E_err + !log_rel_E_err + + !tot_e_equ_err + !tot_e_err + + + !error_in_energy_conservation ! for this step + ! = total_energy - (total_energy_start + total_energy_sources_and_sinks) + !cumulative_energy_error ! = sum over all steps of abs(error_in_energy_conservation) + !rel_cumulative_energy_error ! = cumulative_energy_error/total_energy + log_rel_cumulative_energy_error ! = log10 of rel_cumulative_energy_error + log_rel_run_E_err ! shorter name for rel_cumulative_energy_error + + !rel_error_in_energy_conservation ! = error_in_energy_conservation/total_energy + log_rel_error_in_energy_conservation + + !virial_thm_P_avg + !virial_thm_rel_err + !work_inward_at_center + !work_outward_at_surface + + +!---------------------------------------------------------------------------------------------- + + !# rotation + + !total_angular_momentum + log_total_angular_momentum + !i_rot_total ! moment of inertia + + !total_rotational_kinetic_energy + !log_total_rotational_kinetic_energy + !total_rotational_kinetic_energy_after_adjust_mass + +!---------------------------------------------------------------------------------------------- + +!# velocities + + !avg_abs_v_div_cs + !log_avg_abs_v_div_cs + !max_abs_v_div_cs + !log_max_abs_v_div_cs + + !avg_abs_v + !log_avg_abs_v + !max_abs_v + !log_max_abs_v + + !u_surf + !u_surf_km_s + !u_div_csound_surf + !u_div_csound_max + + !infall_div_cs + +!---------------------------------------------------------------------------------------------- + +!# misc + + !e_thermal ! sum over all zones of Cp*T*dm + + !## eos + !logQ_max ! logQ = logRho - 2*logT + 12 + !logQ_min + !gamma1_min + + !## core mixing + !mass_semiconv_core + + !## H-He boundary + + !diffusion_time_H_He_bdy + !temperature_H_He_bdy + + + !## optical depth and opacity + + !one_div_yphot + !log_one_div_yphot + + !log_min_opacity + !min_opacity + + !log_tau_center + + !log_max_tau_conv + !max_tau_conv + !log_min_tau_conv + !min_tau_conv + + !tau_qhse_yrs + + !## other + + !Lsurf_m + !dlnR_dlnM + !h1_czb_mass ! location (in Msun units) of base of 1st convection zone above he core + !kh_mdot_limit + !log_cntr_dr_cm + !min_Pgas_div_P + !surf_c12_minus_o16 ! this is useful for seeing effects of dredge up on AGB + !surf_num_c12_div_num_o16 + + !phase_of_evolution ! Integer mapping to the type of evolution see star_data/public/star_data_def.inc for definitions + + !## MLT++ + !gradT_excess_alpha + !gradT_excess_min_beta + !gradT_excess_max_lambda + + !max_L_rad_div_Ledd + !max_L_rad_div_Ledd_div_phi_Joss + + + !## RTI + !rti_regions + + !## Ni & Co + !total_ni_co_56 + + + !## internal structure constants + + ! this is evaluated assuming a spherical star and does not account for rotation + !apsidal_constant_k2 + + +!---------------------------------------------------------------------------------------------- + +!# accretion + + !k_below_const_q + !q_below_const_q + !logxq_below_const_q + + !k_const_mass + !q_const_mass + !logxq_const_mass + + !k_below_just_added + !q_below_just_added + !logxq_below_just_added + + !k_for_test_CpT_absMdot_div_L + !q_for_test_CpT_absMdot_div_L + !logxq_for_test_CpT_absMdot_div_L + +!---------------------------------------------------------------------------------------------- + +!# Color output + + ! Outputs the bolometric correction (bc) for the star in filter band ``filter'' (case sensitive) + !bc filter + + ! Outputs the absolute magnitude for the star in filter band ``filter'' (case sensitive) + !abs_mag filter + + ! Adds all the bc's to the output + add_bc + + ! Adds all the absolute magnitudes to the output + add_abs_mag + + ! Outputs luminosity in filter band ``filter'' (erg s^-1) (case sensitive) + !lum_band filter + + ! Adds all the filter band luminosities to the output (erg s^-1) + add_lum_band + + ! Outputs log luminosity in filter band ``filter'' (log erg s^-1) (case sensitive) + !log_lum_band filter + + ! Adds all the filter band luminosities to the output (log erg s^-1) + add_log_lum_band + +!---------------------------------------------------------------------------------------------- + +!# RSP + + !rsp_DeltaMag ! absolute magnitude difference between minimum and maximum light (mag) + !rsp_DeltaR ! R_max - R_min difference in the max and min radius (Rsun) + !rsp_GREKM ! fractional growth of the kinetic energy per pulsation period ("nonlinear growth rate") - see equation 5 in MESA5 + !rsp_num_periods ! Count of the number of pulsation cycles completed + !rsp_period_in_days ! Running period, ie., period between two consecutive values of R_max (days) + !rsp_phase ! Running pulsation phase for a cycle + +!---------------------------------------------------------------------------------------------- +!# debugging + + !## retries + num_retries ! total during the run + + !## solver iterations + + num_iters ! same as num_solver_iterations + !num_solver_iterations ! iterations at this step + !total_num_solver_iterations ! total iterations during the run + !avg_num_solver_iters + + !rotation_solver_steps + + !diffusion_solver_steps + !diffusion_solver_iters + + !avg_setvars_per_step + !avg_skipped_setvars_per_step + !avg_solver_setvars_per_step + + !burn_solver_maxsteps + + !total_num_solver_calls_converged + !total_num_solver_calls_failed + !total_num_solver_calls_made + !total_num_solver_relax_calls_converged + !total_num_solver_relax_calls_failed + !total_num_solver_relax_calls_made + !total_num_solver_relax_iterations + + !total_step_attempts + !total_step_redos + !total_step_retries + !total_steps_finished + !total_steps_taken + + !TDC_num_cells + + !## Relaxation steps + !total_relax_step_attempts + !total_relax_step_redos + !total_relax_step_retries + !total_relax_steps_finished + !total_relax_steps_taken + + !## conservation during mesh adjust + !log_mesh_adjust_IE_conservation + !log_mesh_adjust_KE_conservation + !log_mesh_adjust_PE_conservation + + !## amr + !num_hydro_merges + !num_hydro_splits + + !## timing + !elapsed_time ! time since start of run (seconds) + + !## Extras + burning_regions 40 + mixing_regions 40 + mix_relr_regions 40 diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pgstar b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pgstar new file mode 100644 index 000000000..4a101208a --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pgstar @@ -0,0 +1,751 @@ + +&pgstar + +!pause = .true. + +pgstar_interval = 200 ! making this too small slows the model down. +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 = 10 ! 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) = 'HR' +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) = 'conv_vel_div_csound'!'conv_vel'!'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) = 'v_div_vesc'!'lum_div_Ledd'!'log_csound'!'logT' +Profile_Panels3_other_yaxis_name(5) = 'vel_km_per_s'!'burn_num_iters' + Profile_Panels3_yaxis_log(5) = .false. + +! x-axis limits and properties +Profile_Panels3_xaxis_name = 'logtau'!'zone'!'radius' +Profile_Panels3_xmin = -101d0 +Profile_Panels3_xmax = -101d0 !1700! -101d0 ! 2.2 +Profile_Panels3_xaxis_reversed = .true. + + +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) = 'Eq' +!Profile_Panels1_dymin(1) = 0.14 +Profile_Panels1_other_yaxis_name(1) = 'Uq'!'lum_conv' +!Profile_Panels1_other_dymin(1) = 0.14 + +Profile_Panels1_yaxis_name(2) = 'logT'!'logdq'!'radius'!'entropy' +Profile_Panels1_other_yaxis_name(2) = 'mlt_Y_face'! 'Lc_div_L'!'log_dt_cs_div_dr'!'pgas_div_p' +Profile_Panels1_other_yaxis_log(2) = .false. + +!Profile_Panels1_ymax(2) = 210 +!Profile_Panels1_ymin(2) = 160 + +Profile_Panels1_yaxis_name(3) = 'Lc_div_L'!'gradT' +Profile_Panels1_other_yaxis_name(3) = 'conv_vel'!'lum_div_Ledd'!'grada' +Profile_Panels1_same_yaxis_range(3) = .false. +Profile_Panels1_other_dymin(3) = 0.08 +!Profile_Panels1_yaxis_log(3) = .false. +Profile_Panels1_ymax(3) =-101d0 + +Profile_Panels1_yaxis_name(4) = 'log_opacity' +Profile_Panels1_ymin(4) = 0 +Profile_Panels1_dymin(4) = 0.14 +Profile_Panels1_other_yaxis_name(4) = 'gradr' +Profile_Panels1_other_dymin(4) = 0.14 + +! x-axis limits and properties +Profile_Panels1_xaxis_name = 'logtau'!'zone' +Profile_Panels1_xmin = -101d0!-101d0 +Profile_Panels1_xmax = 5.5d0!-101d0!9d0 !-101d0 !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 = .true. +show_TRho_accretion_mesh_borders = .false. +show_TRho_Profile_kap_regions = .false. +show_TRho_Profile_gamma1_4_3rd = .true. +show_TRho_Profile_eos_regions = .true. +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 = -15.0 +TRho_Profile_xmax = 0!5d0!10.0 +TRho_Profile_ymin = 2.5 +TRho_Profile_ymax = 7.5!8.5d0!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) = 'KE_growth_avg' +Text_Summary1_name(3,2) = 'growth' +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) = 'period'!'log_max_T' +Text_Summary1_name(4,2) = 'log_LH' +Text_Summary1_name(4,3) = 'lg_Lnuc_tot' +Text_Summary1_name(4,4) = 'O_cntr' +Text_Summary1_name(4,5) = 'radius' + +Text_Summary1_name(5,1) = 'num_periods' +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) = 'retries' +Text_Summary1_name(6,5) = 'log_cntr_T' + + +!------------------------------------------------------------------------------------ + +! 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 = 12 + 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 + + + + History_Panels1_xaxis_name = 'yr_since_coll' + + + ! axis limits -- to override system default selections + HR_logT_min = -101d0!3.3 !-101d0 ! only used if /= -101d0 + HR_logT_max = -101d0!3.9!-101d0 ! only used if /= -101d0 + HR_logL_min = -101d0!4.4!-101d0 ! only used if /= -101d0 + HR_logL_max = -101d0!5.5!-101d0 ! only used if /= -101d0 + + ! 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 = 'hr_png' + HR_file_prefix = 'hr_' + HR_file_interval = 10 ! output when `mod(model_number,HR_file_interval)==0` + HR_file_width = 27 ! (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 = 'day'!'model_number' + History_Panels1_max_width = 50! 10000 + + !History_Panels1_xaxis_name = 'star_age' + !History_Panels1_max_width = 10 + + History_Panels1_txt_scale = 0.75 + History_Panels1_xmin = -101d0!1300!200!500 + 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) = 'growth'!'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) = 'luminosity'!'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) = 'v_surf_km_s'!'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) = 'KE_growth_avg'!'log_dt' + 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) = 'time_step_sec' + 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_white_on_black_flag = .false. +!file_device = 'vcps' ! options 'png' and 'vcps' for png and postscript respectively + + +kipp_win_flag=.false. +kipp_file_flag=.false. +Kipp_mix_interval = 1 +Kipp_show_luminosities = .true. + + + +! history tracks for pulsations + + +! history tracks + History_Track1_file_flag = .false. + History_Track2_file_flag = .false. + + +History_Track1_win_flag = .true. +History_Track1_file_interval = 50 +History_Track1_win_width = 12 +History_Track1_win_aspect_ratio = 0.75 + +History_Track1_xname = 'v_surf_km_s' +History_Track1_yname = 'log_L' +History_Track1_xaxis_label = 'Vsurf' +History_Track1_yaxis_label = 'log L/L\d\(2281)' +History_Track1_reverse_xaxis = .false. + + +!History_Track1_xmin = -50d0 +!History_Track1_xmax = 50d0 +!History_Track1_ymin = 3.50d0 +!History_Track1_ymax = 3.98d0 + + +History_Track2_win_flag = .true. +History_Track2_file_interval = 50 + +History_Track2_win_width = 12 +History_Track2_win_aspect_ratio = 0.75 + +History_Track2_xname = 'radius' !'v_surf_km_s' +History_Track2_yname = 'log_L' +History_Track2_xaxis_label = 'Radius' +History_Track2_yaxis_label = 'log L/L\d\(2281)' +History_Track2_reverse_xaxis = .false. + +!History_Track2_xmin = 72d0 +!History_Track2_xmax = 96d0 +!History_Track2_ymin = 3.50d0 +!History_Track2_ymax = 3.98d0 + +/ ! end of pgstar namelist + diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses new file mode 100644 index 000000000..28bd13e75 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses @@ -0,0 +1,245 @@ +&kap +cubic_interpolation_in_X = .false. +cubic_interpolation_in_Z = .false. + +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 = .false. ! if envelope model only +Zbase = 0.004d0 + +/ + + +&eos +/ + +&star_job + + save_model_when_terminate = .true. + save_model_filename = 'final.mod' + !required_termination_code_string = 'Successful test: evolve through 5 pulsation cycles.' + + load_saved_model = .true. + load_model_filename = 'rsp_cepheid_9M_cycle40.mod' + +! turn off RSP and turn on v_flag in mesa star. + + change_initial_RSP_flag = .true. + new_RSP_flag = .false. + + change_v_flag = .true. + new_v_flag = .true. + + relax_to_this_tau_factor = 1d-3 + dlogtau_factor = 0.1d0 + relax_tau_factor = .false. + + set_initial_dt = .true. + years_for_initial_dt = -1 + seconds_for_initial_dt = 1d2 + + set_initial_cumulative_energy_error = .true. + new_cumulative_energy_error = 0d0 + + set_initial_age = .true. + initial_age = 0 + + set_initial_model_number = .true. + initial_model_number = 0 + + pgstar_flag = .true. + +! on first pass, uncomment below if you are loading from and evolution model. + !remove_initial_center_by_temperature = 2d6 + !remove_center_set_zero_v_center = .true. ! set center v = 0 + +/ !end of star_job namelist + + +&controls + tol_max_correction = 1d99 + + ! let the rse know we are in inlist_pulses, useful for multi_inlist setups. + x_logical_ctrl(22) = .true. ! flag for in inlist_pulses + + + terminal_show_age_units = 'days' + terminal_show_timestep_units = 'secs' + terminal_show_log_dt = .false. + terminal_show_log_age = .false. + + +! GYRE output controls, make square with other things below. + !write_pulse_data_with_profile = .true. + pulse_data_format = 'GYRE' + + x_logical_ctrl(37) = .true. ! if true, then run GYRE + + x_integer_ctrl(1) = 1000 ! output GYRE info at this step interval + x_logical_ctrl(1) = .false. ! save GYRE info whenever save profile + + x_integer_ctrl(2) = 3 ! 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) + + + ! limit max model number as part of test suite. + max_model_number = 2000 !3000 + + ! controls for analyzing pulsations using GYRE in MESA + x_integer_ctrl(1) = 1000 ! gyre interval to check + + x_logical_ctrl(7) = .true. ! doing pulses + x_integer_ctrl(7) = -1 ! 3 ! which period to check (<= 0 means don't check any) + x_ctrl(7) = 12d0 ! expected period (in days) + x_ctrl(8) = -1 ! min_deltaR_for_periods (Rsun) + x_ctrl(9) = 1d0 ! KE_growth_avg_abs_frac_new ! for averaging growth rates across steps. + x_ctrl(10) = 0.1 ! min_period_div_target + + + ! new rsp style meshing, thanks to Bill P. + x_logical_ctrl(23) = .false. ! .true. = Remesh for TDC + TDC_hydro_use_mass_interp_face_values = .false. + TDC_hydro_nz = 150 + TDC_hydro_nz_outer = 40 + TDC_hydro_T_anchor = 11d3 + TDC_hydro_dq_1_factor = 2d0 + remesh_for_TDC_pulsations_log_core_zoning = .false. ! .false means do rsp style core + + okay_to_remesh = .false. ! freeze mesh after initial remesh. + x_logical_ctrl(24) = .true. ! if true turn off remesh at the following model number + x_ctrl(12) = 200 ! model number to turn off remesh ( only if if okay_to_remesh = .true.) + + P_function_weight = 0 + min_surface_cell_dq = 1d-8 + max_surface_cell_dq = 1d0 + max_num_merge_surface_cells = 5 + +! GYRE set starting velocities, kick! +! kick when true and not restarting. + x_logical_ctrl(5) = .false. ! to turn on gyre kick + x_ctrl(11) = 200! kick model at this model number + + x_ctrl(4) = 0d0 ! fraction_1st_overtone (order 2) + x_ctrl(5) = 0d0 ! fraction_2nd_overtone (order 3) + x_ctrl(6) = 1d0 ! initial vsurf (kms) + + +! turn off burning and mixing, if evolving envelope model + eps_nuc_factor = 0 + non_nuc_neu_factor = 0 + dxdt_nuc_factor = 0 + mix_factor = 0 + +! timesteps for saturation + + ! TDC Pulsation timestepping. + x_ctrl(13) = 100! model number to drop timestep + x_ctrl(17) = 2d3 ! dt before pulse, seconds + x_ctrl(18) = 4.5d3 ! After a pulse begins, limit the timestep to this (in seconds). + + ! have used these values to do run to saturation, but may not be necessary ~ Bill P. + dt_div_min_dr_div_cs_limit = 2d0! set less than 1d0 in high L/M stars. + dt_div_min_dr_div_cs_hard_limit = 100d0! i don't like hard limits ~ EbF + min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit = 0.8d0 + min_q_for_dt_div_min_dr_div_cs_limit = 0d0 + max_q_for_dt_div_min_dr_div_cs_limit = 1d0 + min_k_for_dt_div_min_dr_div_cs_limit = 1 + min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit = 0.001d0 + + ! artificial viscosity if necessary, v_flag only + use_Pvsc_art_visc = .true. + Pvsc_cq = 4.0d0 + Pvsc_zsh = 0.1d0 + + ! velocity time centering for v_flag only. + steps_before_use_velocity_time_centering = 300 ! no v centering when commented + use_P_d_1_div_rho_form_of_work_when_time_centering_velocity = .true. + + include_P_in_velocity_time_centering = .true. ! set to false for u_flag + P_theta_for_velocity_time_centering = 0.5d0 + + include_L_in_velocity_time_centering = .true. + L_theta_for_velocity_time_centering = 0.5d0 + + include_mlt_in_velocity_time_centering = .false. + + set_rho_to_dm_div_dV = .false. + +! OUTER BC for TDC Pulsations + use_RSP_L_eqn_outer_BC = .true. ! uses RSP2_Lsurf_factor + RSP2_Lsurf_factor = 0.5d0 + use_zero_Pgas_outer_BC = .true. + use_compression_outer_BC = .false. + use_momentum_outer_BC = .false. ! adopt this if not using rsp and zero P bc. + +! Convection model + steps_before_use_TDC = 0 + MLT_option = 'TDC' + include_mlt_corr_to_TDC = .false. ! true K. 1986 model, no mlt limiting + use_face_values_eos_and_kap_mlt_tdc = .true. ! build smarter interpolant for face quantities + TDC_include_eturb_in_energy_equation = .true. + use_rsp_form_of_scale_height = .true. + mixing_length_alpha = 1.5d0 + TDC_alpha_M = 0.25d0 + use_ledoux_criterion = .true. + alt_scale_height_flag = .false. ! ignore eggleton, not compatible with rsp_form_of_Hp + TDC_num_innermost_cells_forced_nonturbulent = 2 ! for envelope models only. + +! controls for shock capturing, relaxed for pulsations + ! main purpose is to force radiative in shock face + max_abs_du_div_cs_for_convection = 10d0 + max_v_div_cs_for_convection = 1d2 + max_v_for_convection = 1d4 + +! output + log_directory = 'LOGS_pulsation' + + photo_digits = 8 + photo_interval = 1000 + profile_interval = 10000 + history_interval = 1 + terminal_interval = 10 + max_num_profile_models = 1000000 ! 100 ! RECOMMENDED 10000 + +! terminal output + num_trace_history_values = 2 + trace_history_value_name(1) = 'log_rel_run_E_err' + trace_history_value_name(2) = 'rel_E_err' + warn_when_large_rel_run_E_err = 1d-2 + +! solver, relax solver hard for summer school + use_gold2_tolerances = .true. + !gold_iter_for_resid_tol2 = 3 ! default is 5 + !gold_tol_residual_norm3 = 1d-4 + !gold_tol_max_residual3= 1d-2 + !make_gradr_sticky_in_solver_iters = .true. + +! Debugging + report_solver_progress = .false. + report_ierr = .false. ! if true, produce terminal output when have some internal error + +/ ! end of controls namelist + +&pgstar + +Grid1_file_dir = 'png_pulsation' + +pgstar_interval = 100 + + +! axis limits -- to override system default selections +!HR_logT_min = 3.6 !-101d0 ! only used if /= -101d0 +!HR_logT_max = 3.75!-101d0 ! only used if /= -101d0 +!HR_logL_min = 3.85!-101d0 ! only used if /= -101d0 +!HR_logL_max = 3.6!-101d0 ! only used if /= -101d0 + +HR_file_interval = 50 ! output when `mod(model_number,HR_file_interval)==0` + +/ ! end of pgstar namelist diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses_header b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses_header new file mode 100644 index 000000000..be75026dd --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/inlist_pulses_header @@ -0,0 +1,50 @@ + +&star_job + + read_extra_star_job_inlist(1) = .true. + extra_star_job_inlist_name(1) = 'inlist_pulses' + +/ ! end of star_job namelist + + +&eos + + read_extra_eos_inlist(2) = .true. + extra_eos_inlist_name(2) = 'inlist_pulses' + +/ ! end of eos namelist + + +&kap + + read_extra_kap_inlist(2) = .true. + extra_kap_inlist_name(2) = 'inlist_pulses' + +/ ! end of kap namelist + + +&colors + + read_extra_colors_inlist(1) = .false. + extra_colors_inlist_name(1) = 'inlist_pulses' + +/ ! end of colors namelist + + +&controls + + read_extra_controls_inlist(1) = .true. + extra_controls_inlist_name(1)= 'inlist_pulses' + +/ ! end of controls namelist + + +&pgstar + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pulses' + + read_extra_pgstar_inlist(1) = .true. + extra_pgstar_inlist_name(1)= 'inlist_pgstar' + +/ ! end of pgstar namelist diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/profile_columns.list b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/profile_columns.list new file mode 100644 index 000000000..7087d449c --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/profile_columns.list @@ -0,0 +1,966 @@ +! 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 + dwork_dm + + !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 + Frad_div_cUrad ! Frad/(C*Urad), must be < 1 to not violate c. + !flux_limit_lambda + !flux_limit_R + !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_TDC_Pulsation/dev_TDC_Cepheid_9M/re b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/re new file mode 100755 index 000000000..7093783a3 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/re @@ -0,0 +1,31 @@ +#!/usr/bin/env 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 + +make run diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/re_nomodfiles b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/re_nomodfiles new file mode 100755 index 000000000..d447f6665 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/re_nomodfiles @@ -0,0 +1,41 @@ +#!/usr/bin/env 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_TDC_Pulsation/dev_TDC_Cepheid_9M/rn b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn new file mode 100755 index 000000000..910c56f82 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn @@ -0,0 +1,14 @@ +#!/usr/bin/env 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_pulses_header final.mod + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" + +echo 'finished' + diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn1 b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn1 new file mode 100755 index 000000000..7eb455c6e --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn1 @@ -0,0 +1,7 @@ +#!/usr/bin/env bash + +rm -f restart_photo + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./build/bin/star +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn_nomodfiles b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn_nomodfiles new file mode 100755 index 000000000..f2217502a --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rn_nomodfiles @@ -0,0 +1,7 @@ +#!/usr/bin/env bash +rm -f restart_photo +echo $1 + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./build/bin/star $1 +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rsp_cepheid_9M_cycle40.mod b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rsp_cepheid_9M_cycle40.mod new file mode 100644 index 000000000..c5ca7a952 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/rsp_cepheid_9M_cycle40.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1504617dac77bedac0a5f6e99edb7c60f52fe209f3633ea74811f2724a315206 +size 80845 diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run.f90 b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run.f90 new file mode 100644 index 000000000..76d423f1a --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/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 run diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras.f90 b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras.f90 new file mode 100644 index 000000000..1d91f08e5 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras.f90 @@ -0,0 +1,321 @@ +! *********************************************************************** +! +! Copyright (C) 2010-2025 The MESA Team +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License +! as published by the Free Software Foundation, +! either version 3 of the License, or (at your option) any later version. +! +! This program 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with this program. If not, see . +! +! *********************************************************************** + +module run_star_extras + + use star_lib + use star_def + use const_def + use math_lib + use auto_diff + use chem_def + use utils_lib + use gyre_mesa_m + + use interp_1d_def, only: pm_work_size + use interp_1d_lib, only: interp_pm, interp_values, interp_value + + implicit none + + include 'run_star_extras_TDC_pulsation_defs.inc' + + logical :: dbg = .false. + + !!!!!!!!!!!!!!!!!!!!!!!!! + ! These variables are loaded up from x_ctrl, x_integer_ctrl and x_logical_ctrl + ! values specified on inlist_common, inlist_pulses + !!!!!!!!!!!!!!!!!!!!!!!!! + + logical :: in_inlist_pulses, remesh_for_envelope_model, turn_off_remesh + integer :: kick_model_number, timestep_drop_model_number, turn_off_remesh_model_number + integer :: initial_model_number + real(dp) :: max_dt_before_pulse, max_dt_during_pulse + +contains + + include 'run_star_extras_TDC_pulsation.inc' + + subroutine extras_controls(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type(star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + s%extras_startup => extras_startup + s%extras_start_step => extras_start_step + s%extras_check_model => extras_check_model + s%extras_finish_step => extras_finish_step + s%extras_after_evolve => extras_after_evolve + s%how_many_extra_history_columns => how_many_extra_history_columns + s%data_for_extra_history_columns => data_for_extra_history_columns + s%how_many_extra_profile_columns => how_many_extra_profile_columns + s%data_for_extra_profile_columns => data_for_extra_profile_columns + + ! pulsation info + s%other_photo_write => photo_write + s%other_photo_read => photo_read + + ! store user provided options from the inlist + + in_inlist_pulses = s%x_logical_ctrl(22) + max_dt_before_pulse = s%x_ctrl(17) + max_dt_during_pulse = s%x_ctrl(18) + remesh_for_envelope_model = s%x_logical_ctrl(23) + turn_off_remesh = s%x_logical_ctrl(24) + kick_model_number = s%x_ctrl(11) + timestep_drop_model_number = s%x_ctrl(13) + turn_off_remesh_model_number = s%x_ctrl(12) + end subroutine extras_controls + + + subroutine extras_startup(id, restart, ierr) + integer, intent(in) :: id + logical, intent(in) :: restart + integer, intent(out) :: ierr + type(star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + call TDC_pulsation_extras_startup(id, restart, ierr) + + ! Initialize GYRE + + call init('gyre.in') + + ! Set constants + + call set_constant('G_GRAVITY', standard_cgrav) + call set_constant('C_LIGHT', clight) + call set_constant('A_RADIATION', crad) + + call set_constant('M_SUN', Msun) + call set_constant('R_SUN', Rsun) + call set_constant('L_SUN', Lsun) + + call set_constant('GYRE_DIR', TRIM(mesa_dir)//'/build/gyre/src') + + if (.not. restart .and. in_inlist_pulses) then + initial_model_number = s% model_number + end if + !initial_model_number = 0 ! since we are setting model # to 0 in inlist_pulses + + ! for rsp style mesh + if (.not. restart .and. in_inlist_pulses .and. remesh_for_envelope_model) then + call remesh_for_TDC_pulsation(id, ierr) + end if + end subroutine extras_startup + + subroutine extras_after_evolve(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type(star_info), pointer :: s + real(dp) :: dt + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + if (.not. s%x_logical_ctrl(37)) return + call final() + end subroutine extras_after_evolve + + ! returns either keep_going, retry, or terminate. + integer function extras_check_model(id) + integer, intent(in) :: id + integer :: ierr, k + real(dp) :: max_v + type(star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_check_model = keep_going + + end function extras_check_model + + integer function how_many_extra_history_columns(id) + integer, intent(in) :: id + integer :: ierr + type(star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_history_columns = TDC_pulsation_how_many_extra_history_columns(id) + end function how_many_extra_history_columns + + subroutine data_for_extra_history_columns(id, n, names, vals, ierr) + integer, intent(in) :: id, n + character(len=maxlen_history_column_name) :: names(n) + real(dp) :: vals(n), v_esc + integer, intent(out) :: ierr + type(star_info), pointer :: s + integer :: k, k0 + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + call TDC_pulsation_data_for_extra_history_columns(id, n, names, vals, ierr) + 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 = TDC_pulsation_how_many_extra_profile_columns(id) + + 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 + call TDC_pulsation_data_for_extra_profile_columns(id, n, nz, names, vals, ierr) + + end subroutine data_for_extra_profile_columns + + integer function extras_start_step(id) + integer, intent(in) :: id + integer :: ierr + type(star_info), pointer :: s + include 'formats' + extras_start_step = terminate + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + ! we want to ignore T gradient equation for a few steps after remesh + if (s%model_number < initial_model_number + 10 .and. in_inlist_pulses) then + s%convergence_ignore_equL_residuals = .true. + else if (in_inlist_pulses) then + s%convergence_ignore_equL_residuals = .false. + end if + + if (s%model_number == kick_model_number .and. in_inlist_pulses & + .and. s%x_logical_ctrl(5)) then + + ! if v= 0, turn on v so we can kick + if (.not. s%v_flag .and. .not. s%u_flag) then + call star_set_v_flag(id, .true., ierr) + write (*,*) 'turning on v_flag hydro for kick' + end if + + call gyre_in_mesa_extras_set_velocities(s, .false., ierr) + write (*, *) 'kick' + write (*, *) 'kick' + write (*, *) 'kick' + write (*, *) 'kick' + write (*, *) 'kick' + + end if + + call my_before_struct_burn_mix(s%id, s%dt, extras_start_step) + + extras_start_step = keep_going + end function extras_start_step + + subroutine my_before_struct_burn_mix(id, dt, res) + use const_def, only: dp + use star_def + integer, intent(in) :: id + real(dp), intent(in) :: dt + integer, intent(out) :: res ! keep_going, redo, retry, terminate + real(dp) :: power_photo, v_esc + integer :: ierr, k + type(star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + if (in_inlist_pulses) then + if (s%model_number > timestep_drop_model_number) then + s%max_timestep = max_dt_during_pulse + else + s%max_timestep = max_dt_before_pulse + end if + + ! time step control on pulsations + if (period > 0d0 .and. period/s%max_timestep < 600 .and. & + s%model_number > timestep_drop_model_number) then + s%max_timestep = period/600d0 + end if + + if (s%model_number > turn_off_remesh_model_number .and. turn_off_remesh) then + s%okay_to_remesh = .false. + end if + end if + + res = keep_going + end subroutine my_before_struct_burn_mix + + subroutine null_binary_controls(id, binary_id, ierr) + integer, intent(in) :: id, binary_id + integer, intent(out) :: ierr + ierr = 0 + end subroutine null_binary_controls + + ! returns either keep_going or terminate. + integer function extras_finish_step(id) + use run_star_support + use math_lib + integer, intent(in) :: id + integer :: ierr, k + real(dp) :: max_vel_inside, vesc_for_cell, vesc_surf !check_avg_v_div_vesc + type(star_info), pointer :: s + include 'formats' + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + extras_finish_step = keep_going + extras_finish_step = TDC_pulsation_extras_finish_step(id) + +! if (.not. s% x_logical_ctrl(37)) return +! extras_finish_step = gyre_in_mesa_extras_finish_step(id) + + if (extras_finish_step == terminate) s%termination_code = t_extras_finish_step + + end function extras_finish_step + + + subroutine photo_write(id, iounit) + integer, intent(in) :: id, iounit + call TDC_pulsation_photo_write(id, iounit) + end subroutine photo_write + + subroutine photo_read(id, iounit, ierr) + integer, intent(in) :: id, iounit + integer, intent(out) :: ierr + call TDC_pulsation_photo_read(id, iounit, ierr) + end subroutine photo_read + +end module run_star_extras diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation.inc b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation.inc new file mode 100644 index 000000000..e1d2ca8ca --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation.inc @@ -0,0 +1,383 @@ + + + subroutine TDC_pulsation_photo_write(id, iounit) + integer, intent(in) :: id, iounit + write(iounit) num_periods, run_num_steps_end_prev, & + run_num_iters_end_prev, run_num_retries_end_prev, & + period, KE_growth, KE_growth_avg, prev_KE_max, & + delta_R_growth, delta_R_growth_avg, prev_delta_R, & + period_max_v_div_vesc, period_max_v_div_cs, period_delta_R, & + period_delta_Teff, period_delta_logTeff, & + period_delta_logL, period_delta_Mag, & + time_started, v_div_cs_max, v_div_vesc_max, & + KE_min, KE_max, R_min, R_max, L_min, L_max, T_min, T_max, & + best_period, best_model_number, best_order, best_cycles_to_double + end subroutine TDC_pulsation_photo_write + + + subroutine TDC_pulsation_photo_read(id, iounit, ierr) + integer, intent(in) :: id, iounit + integer, intent(out) :: ierr + ierr = 0 + read(iounit, iostat=ierr) num_periods, run_num_steps_end_prev, & + run_num_iters_end_prev, run_num_retries_end_prev, & + period, KE_growth, KE_growth_avg, prev_KE_max, & + delta_R_growth, delta_R_growth_avg, prev_delta_R, & + period_max_v_div_vesc, period_max_v_div_cs, period_delta_R, & + period_delta_Teff, period_delta_logTeff, & + period_delta_logL, period_delta_Mag, & + time_started, v_div_cs_max, v_div_vesc_max, & + KE_min, KE_max, R_min, R_max, L_min, L_max, T_min, T_max, & + best_period, best_model_number, best_order, best_cycles_to_double + end subroutine TDC_pulsation_photo_read + + + subroutine TDC_pulsation_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 + if (.not. restart) then + num_periods = 0 + run_num_steps_end_prev = 0 + run_num_iters_end_prev = 0 + run_num_retries_end_prev = 0 + period = 0 + KE_growth = 0 + KE_growth_avg = 0 + prev_KE_max = 0 + delta_R_growth = 0 + delta_R_growth_avg = 0 + prev_delta_R = 0 + period_max_v_div_cs = 0 + period_max_v_div_vesc = 0 + period_delta_R = 0 + period_delta_Teff = 0 + period_delta_logTeff = 0 + period_delta_logL = 0 + period_delta_Mag = 0 + time_started = 0 + v_div_cs_max = 0 + v_div_vesc_max = 0 + KE_min = 0 + KE_max = 0 + R_min = 0 + R_max = 0 + L_min = 0 + L_max = 0 + T_min = 0 + T_max = 0 + best_period = 0 + best_model_number = 0 + best_order = 0 + best_cycles_to_double = 0 + end if + end subroutine TDC_pulsation_extras_startup + + + integer function TDC_pulsation_extras_finish_step(id) + use chem_def + integer, intent(in) :: id + type (star_info), pointer :: s + integer :: ierr, gyre_interval, test_period + real(dp) :: target_period + logical :: doing_pulses + include 'formats' + + TDC_pulsation_extras_finish_step = terminate + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + TDC_pulsation_extras_finish_step = keep_going + + gyre_interval = s% x_integer_ctrl(1) + if (gyre_interval > 0) then + if (MOD(s% model_number, gyre_interval) == 0) & + call get_gyre_info_for_this_step + if (TDC_pulsation_extras_finish_step == terminate) & + s% termination_code = t_extras_finish_step + if (TDC_pulsation_extras_finish_step /= keep_going) return + end if + + doing_pulses = s% x_logical_ctrl(7) + if (.not. doing_pulses) return + target_period = s% x_ctrl(7) + if (target_period <= 0d0) return + if (.not. get_period_info()) return + + test_period = s% x_integer_ctrl(7) + if (num_periods < test_period .or. test_period <= 0) return + + ! have finished test run + call report_test_results + TDC_pulsation_extras_finish_step = terminate + + contains + + subroutine get_gyre_info_for_this_step + integer :: i + TDC_pulsation_extras_finish_step = gyre_in_mesa_extras_finish_step(id) + if (s% ixtra3_array(1) > 0) then ! unpack the GYRE results + do i=1,s% ixtra3_array(1) + if (s% xtra1_array(i) == 0d0 .or. & + (s% ixtra1_array(i) /= s% x_integer_ctrl(4) .and. s% x_integer_ctrl(4) > 0)) cycle + if (s% xtra3_array(i) > 0d0 .and. & + (best_cycles_to_double == 0d0 .or. s% xtra3_array(i) < best_cycles_to_double)) then + !best_growth = s% xtra1_array(i) + best_period = 1d0/s% xtra2_array(i) ! xtra2_array = freq (s^-1) + best_period = best_period/(24*3600) ! change to days + best_cycles_to_double = s% xtra3_array(i) + best_order = s% ixtra1_array(i) + best_model_number = s% model_number + end if + end do + end if + end subroutine get_gyre_info_for_this_step + + logical function get_period_info() + real(dp) :: v_surf, v_surf_start, KE, min_period, time_ended, & + delta_R, min_deltaR_for_periods, growth_avg_frac_new, & + min_period_div_target, cs, vesc + include 'formats' + get_period_info = .false. + + if (s% r(1) < R_min) R_min = s% r(1) + if (s% r(1) > R_max) R_max = s% r(1) + if (s% L(1) < L_min) L_min = s% L(1) + if (s% L(1) > L_max) L_max = s% L(1) + if (s% Teff < T_min) T_min = s% Teff + if (s% Teff > T_max) T_max = s% Teff + KE = s% total_radial_kinetic_energy_end + if (KE > KE_max) KE_max = KE + if (KE < KE_min) KE_min = KE + + if (s% v_flag) then + v_surf = s% v(1) + v_surf_start = s% v_start(1) + else if (s% u_flag) then + v_surf = s% u(1) ! s% u_face_val(1) + v_surf_start = s% u_start(1) ! s% u_face_start(1) + else ! + v_surf = 0d0 + v_surf_start = 0d0 + !stop 'TDC_pulsation_extras_finish_step: both v_flag and u_flag are false' + end if + cs = s% csound(1) + if (v_surf > v_div_cs_max*cs) v_div_cs_max = v_surf/cs + vesc = sqrt(2*s% cgrav(1)*s% m(1)/(s% r(1))) + if (v_surf > v_div_vesc_max*vesc) v_div_vesc_max = v_surf/vesc + + ! period is completed when v_surf goes from positive to negative during step + if (v_surf > 0d0 .or. v_surf_start < 0d0) return + + if (time_started == 0) then ! start of 1st cycle + time_started = s% time + run_num_steps_end_prev = s% model_number + run_num_iters_end_prev = s% total_num_solver_iterations + run_num_retries_end_prev = s% num_retries + prev_KE_max = 0d0 + call init_min_max_info + write(*,*) 'first maximum radius, period calculations starting at model, day', & + s% model_number, s% time/(24*3600) + return + end if + + delta_R = R_max - R_min + min_deltaR_for_periods = s% x_ctrl(8)*Rsun + if (min_deltaR_for_periods > 0d0) then + if (delta_R < min_deltaR_for_periods) return ! filter out glitches + end if + + time_ended = s% time + if (abs(v_surf - v_surf_start) > 1d-10) & ! tweak the end time to match when v_surf == 0 + time_ended = s% time - v_surf*s% dt/(v_surf - v_surf_start) + min_period_div_target = s% x_ctrl(10) + min_period = target_period*(24*3600)*min_period_div_target + if (min_period > 0d0 .and. & + time_ended - time_started < min_period) return ! filter out glitches + + period = time_ended - time_started + num_periods = num_periods + 1 + + if (num_periods > 1) then + growth_avg_frac_new = s% x_ctrl(9) + KE_growth = (KE_max - prev_KE_max)/prev_KE_max + KE_growth_avg = growth_avg_frac_new*KE_growth + & + (1d0 - growth_avg_frac_new)*KE_growth_avg + delta_R_growth = (delta_R - prev_delta_R)/prev_delta_R + delta_R_growth_avg = growth_avg_frac_new*delta_R_growth + & + (1d0 - growth_avg_frac_new)*delta_R_growth_avg + end if + + period_delta_Teff = T_max - T_min + period_delta_logTeff = log10(T_max/T_min) + period_delta_R = R_max - R_min + period_delta_logL = log10(L_max/L_min) + period_delta_Mag = 2.5d0*period_delta_logL + period_max_v_div_cs = v_div_cs_max + period_max_v_div_vesc = v_div_vesc_max + prev_KE_max = KE_max + prev_delta_R = period_delta_R + ! 1 2 3 4 5 6 7 8 9 + write(*,'(i4,a14,i6,a13,f8.3,a13,f9.3,a9,f9.4,a15,f10.4,a13,f9.4,a13,f10.4,a11,f9.4,a13,f9.4)') & + num_periods, & + 'steps/cycle', s% model_number - run_num_steps_end_prev, & ! 1 a14,i6 + 'iters/step', & + dble(s% total_num_solver_iterations - run_num_iters_end_prev)/ & + dble(s% model_number - run_num_steps_end_prev), & ! 2 a13,f8.3 + 'period (d)', period/(24*3600), & ! 3 a13,f9.3 + 'growth', delta_R_growth_avg, & ! 4 a9,f9.4 + 'delta R/Rsun', period_delta_R/Rsun, & ! 5 a15,f10.4 + 'delta logL', period_delta_logL, & ! 6 a13,f9.4 + 'delta Teff', period_delta_Teff, & ! 7 a13,f10.4 + 'max v/cs', period_max_v_div_cs, & ! 8 a11,f9.4 + 'max v/vesc', period_max_v_div_vesc ! 9 a13,f9.4 + + time_started = time_ended + run_num_steps_end_prev = s% model_number + run_num_iters_end_prev = s% total_num_solver_iterations + run_num_retries_end_prev = s% num_retries + call init_min_max_info + get_period_info = .true. + + end function get_period_info + + subroutine init_min_max_info + v_div_cs_max = 0d0 + v_div_vesc_max = 0d0 + KE_min = 1d99 + KE_max = -1d99 + R_min = 1d99 + R_max = -1d99 + L_min = 1d99 + L_max = -1d99 + T_min = 1d99 + T_max = -1d99 + end subroutine init_min_max_info + + subroutine report_test_results + real(dp) :: rel_run_E_err + write(*,*) + write(*,*) + write(*,*) + rel_run_E_err = s% cumulative_energy_error/s% total_energy + write(*,*) 'rel_run_E_err', rel_run_E_err + if (s% total_energy /= 0d0 .and. abs(rel_run_E_err) > 1d-5) then + write(*,*) '*** BAD rel_run_E_error ***', & + s% cumulative_energy_error/s% total_energy + else if (abs(period/(24*3600) - target_period) > 1d-2) then + write(*,*) '*** BAD period ***', period/(24*3600) - target_period, & + period/(24*3600), target_period + else + write(*,*) 'good match for period', & + period/(24*3600), target_period + end if + write(*,*) + write(*,*) + write(*,*) + end subroutine report_test_results + + end function TDC_pulsation_extras_finish_step + + + include 'gyre_in_mesa_extras_finish_step.inc' + + + subroutine TDC_pulsation_extras_after_evolve(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + ierr = 0 + call final() + end subroutine TDC_pulsation_extras_after_evolve + + + integer function TDC_pulsation_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 + TDC_pulsation_how_many_extra_history_columns = 11 + end function TDC_pulsation_how_many_extra_history_columns + + + subroutine TDC_pulsation_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 + integer :: i + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + i = 1 + names(i) = 'num_periods'; vals(i) = num_periods; i=i+1 + names(i) = 'period'; vals(i) = period/(24*3600); i=i+1 + names(i) = 'growth'; vals(i) = delta_R_growth_avg; i=i+1 + names(i) = 'max_v_div_cs'; vals(i) = period_max_v_div_cs; i=i+1 + names(i) = 'max_v_div_vesc'; vals(i) = period_max_v_div_vesc; i=i+1 + names(i) = 'delta_R'; vals(i) = period_delta_R/Rsun; i=i+1 + names(i) = 'delta_Teff'; vals(i) = period_delta_Teff; i=i+1 + names(i) = 'delta_logTeff'; vals(i) = period_delta_logTeff; i=i+1 + names(i) = 'delta_logL'; vals(i) = period_delta_logL; i=i+1 + names(i) = 'delta_Mag'; vals(i) = period_delta_Mag; i=i+1 + names(i) = 'KE_growth_avg'; vals(i) = KE_growth_avg; i = i+1 + + end subroutine TDC_pulsation_data_for_extra_history_columns + + + integer function TDC_pulsation_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 + TDC_pulsation_how_many_extra_profile_columns = 0 ! 6 + end function TDC_pulsation_how_many_extra_profile_columns + + + subroutine TDC_pulsation_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 + return + + + call star_ptr(id, s, ierr) + if (ierr /= 0) return + + names(1) = 'xtra1' + names(2) = 'xtra2' + names(3) = 'xtra3' + names(4) = 'xtra4' + names(5) = 'xtra5' + names(6) = 'xtra6' + + do k=1,nz + vals(k,1) = s% xtra1_array(k) + vals(k,2) = s% xtra2_array(k) + vals(k,3) = s% xtra3_array(k) + vals(k,4) = s% xtra4_array(k) + vals(k,5) = s% xtra5_array(k) + vals(k,6) = s% xtra6_array(k) + end do + + end subroutine TDC_pulsation_data_for_extra_profile_columns + + + include 'gyre_in_mesa_extras_set_velocities.inc' + + diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation_defs.inc b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation_defs.inc new file mode 100644 index 000000000..b7b6e47b6 --- /dev/null +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_9M/src/run_star_extras_TDC_pulsation_defs.inc @@ -0,0 +1,16 @@ + + ! GYRE "best" info + real(dp) :: best_period, best_cycles_to_double + integer :: best_model_number, best_order + + ! summary info at time of recently completely period + integer :: num_periods, run_num_steps_end_prev, & + run_num_iters_end_prev, run_num_retries_end_prev + real(dp) :: period, KE_growth, KE_growth_avg, prev_KE_max, & + delta_R_growth, delta_R_growth_avg, prev_delta_R, & + period_max_v_div_vesc, period_max_v_div_cs, period_delta_R, & + period_delta_Teff, period_delta_logTeff, & + period_delta_logL, period_delta_Mag + ! info for period in progress + real(dp) :: time_started, v_div_cs_max, v_div_vesc_max, & + KE_min, KE_max, R_min, R_max, L_min, L_max, T_min, T_max diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 675345b21..b94e58e80 100644 --- a/star/private/alloc.f90 +++ b/star/private/alloc.f90 @@ -1086,6 +1086,28 @@ subroutine star_info_arrays(s, c_in, action_in, ierr) if (failed('mlt_D_ad')) exit call do1_ad(s% mlt_Gamma_ad, c% mlt_Gamma_ad) if (failed('mlt_Gamma_ad')) exit + call do1_logical(s% have_mlt_tdc_face_state, c% have_mlt_tdc_face_state) + if (failed('have_mlt_tdc_face_state')) exit + call do1_ad(s% mlt_tdc_T_face_ad, c% mlt_tdc_T_face_ad) + if (failed('mlt_tdc_T_face_ad')) exit + call do1_ad(s% mlt_tdc_rho_face_ad, c% mlt_tdc_rho_face_ad) + if (failed('mlt_tdc_rho_face_ad')) exit + call do1_ad(s% mlt_tdc_P_face_ad, c% mlt_tdc_P_face_ad) + if (failed('mlt_tdc_P_face_ad')) exit + call do1_ad(s% mlt_tdc_Cp_face_ad, c% mlt_tdc_Cp_face_ad) + if (failed('mlt_tdc_Cp_face_ad')) exit + call do1_ad(s% mlt_tdc_ChiRho_face_ad, c% mlt_tdc_ChiRho_face_ad) + if (failed('mlt_tdc_ChiRho_face_ad')) exit + call do1_ad(s% mlt_tdc_ChiT_face_ad, c% mlt_tdc_ChiT_face_ad) + if (failed('mlt_tdc_ChiT_face_ad')) exit + call do1_ad(s% mlt_tdc_grada_face_ad, c% mlt_tdc_grada_face_ad) + if (failed('mlt_tdc_grada_face_ad')) exit + call do1_ad(s% mlt_tdc_opacity_face_ad, c% mlt_tdc_opacity_face_ad) + if (failed('mlt_tdc_opacity_face_ad')) exit + call do1_ad(s% mlt_tdc_scale_height_face_ad, c% mlt_tdc_scale_height_face_ad) + if (failed('mlt_tdc_scale_height_face_ad')) exit + call do1_ad(s% mlt_tdc_gradr_face_ad, c% mlt_tdc_gradr_face_ad) + if (failed('mlt_tdc_gradr_face_ad')) exit call do1_ad(s% PII_ad, c% PII_ad) if (failed('PII_ad')) exit @@ -1198,6 +1220,8 @@ subroutine star_info_arrays(s, c_in, action_in, ierr) if (failed('Peos_start')) exit call do1(s% Peos_face_start, c% Peos_face_start) if (failed('Peos_face_start')) exit + call do1(s% mlt_tdc_P_face_start, c% mlt_tdc_P_face_start) + if (failed('mlt_tdc_P_face_start')) exit call do1(s% lnT_start, c% lnT_start) if (failed('lnT_start')) exit call do1(s% energy_start, c% energy_start) diff --git a/star/private/conv_premix.f90 b/star/private/conv_premix.f90 index 25bc99e6c..b93a5cfe1 100644 --- a/star/private/conv_premix.f90 +++ b/star/private/conv_premix.f90 @@ -1339,6 +1339,13 @@ subroutine update_model_ (s, update_mode, kc_t, kc_b) kf_t = kc_t kf_b = kc_b + 1 + if (s% use_face_values_eos_and_kap_mlt_tdc) then + ! update_model_ changed the local composition and refreshed EOS and kap + ! on kc_t:kc_b, so the cached face thermo bundle must be rebuilt for the + ! interior faces before set_mlt_vars uses it. + s% have_mlt_tdc_face_state(kf_t+1:kf_b-1) = .false. + end if + call set_mlt_vars(s, kf_t+1, kf_b-1, ierr) if (ierr /= 0) then write(*,*) 'conv_premix: failed in call to set_mlt_vars during update_model_' diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index f9cae2aaa..21d1d2bfb 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -111,7 +111,9 @@ module ctrls_io TDC_use_density_form_for_eddy_viscosity, & TDC_num_innermost_cells_forced_nonturbulent, TDC_num_outermost_cells_forced_nonturbulent, & include_mlt_Pturb_in_thermodynamic_gradients, & - include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter, TDC_include_eturb_in_energy_equation, & + include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter, & + use_face_values_eos_and_kap_mlt_tdc, & + TDC_include_eturb_in_energy_equation, & use_rsp_form_of_scale_height, include_mlt_in_velocity_time_centering, & TDC_hydro_use_mass_interp_face_values, TDC_hydro_nz, TDC_hydro_nz_outer, TDC_hydro_T_anchor, TDC_hydro_dq_1_factor, & @@ -2097,6 +2099,7 @@ subroutine store_controls(s, ierr) s% include_mlt_Pturb_in_thermodynamic_gradients = include_mlt_Pturb_in_thermodynamic_gradients s% include_mlt_corr_to_TDC = include_mlt_corr_to_TDC s% use_TDC_enthalpy_flux_limiter = use_TDC_enthalpy_flux_limiter + s% use_face_values_eos_and_kap_mlt_tdc = use_face_values_eos_and_kap_mlt_tdc s% TDC_include_eturb_in_energy_equation = TDC_include_eturb_in_energy_equation s% use_rsp_form_of_scale_height = use_rsp_form_of_scale_height s% include_mlt_in_velocity_time_centering = include_mlt_in_velocity_time_centering @@ -3808,6 +3811,7 @@ subroutine set_controls_for_writing(s, ierr) include_mlt_Pturb_in_thermodynamic_gradients = s% include_mlt_Pturb_in_thermodynamic_gradients include_mlt_corr_to_TDC = s% include_mlt_corr_to_TDC use_TDC_enthalpy_flux_limiter = s% use_TDC_enthalpy_flux_limiter + use_face_values_eos_and_kap_mlt_tdc = s% use_face_values_eos_and_kap_mlt_tdc TDC_include_eturb_in_energy_equation = s% TDC_include_eturb_in_energy_equation use_rsp_form_of_scale_height = s% use_rsp_form_of_scale_height include_mlt_in_velocity_time_centering = s% include_mlt_in_velocity_time_centering diff --git a/star/private/hydro_temperature.f90 b/star/private/hydro_temperature.f90 index c6e84ef2b..c6f15a31b 100644 --- a/star/private/hydro_temperature.f90 +++ b/star/private/hydro_temperature.f90 @@ -21,6 +21,7 @@ module hydro_temperature use star_private_def use const_def, only: dp, ln10, pi4, crad, clight, convective_mixing + use mlt_tdc_face_support, only: get_face_eos_kap_ad use utils_lib, only: mesa_error, is_bad use auto_diff use auto_diff_support @@ -51,6 +52,7 @@ subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr) type(auto_diff_real_star_order1) :: L_ad, r_00, area, area2, Lrad_ad, & kap_00, kap_m1, kap_face, d_P_rad_expected_ad, T_m1, T4_m1, T_00, T4_00, & P_rad_m1, P_rad_00, d_P_rad_actual_ad, resid + type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face type(auto_diff_real_star_order1) :: flxR, flxLambda integer :: i_equL @@ -88,9 +90,19 @@ subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr) Lrad_ad = L_ad end if - kap_00 = wrap_kap_00(s,k) - kap_m1 = wrap_kap_m1(s,k) - kap_face = alfa*kap_00 + beta*kap_m1 + if (s% use_face_values_eos_and_kap_mlt_tdc) then + if (s% have_mlt_tdc_face_state(k)) then + kap_face = s% mlt_tdc_opacity_face_ad(k) + else + call get_face_eos_kap_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, kap_face, ierr) + if (ierr /= 0) return + end if + else + kap_00 = wrap_kap_00(s,k) + kap_m1 = wrap_kap_m1(s,k) + kap_face = alfa*kap_00 + beta*kap_m1 + end if if (kap_face%val < s% min_kap_for_dPrad_dm_eqn) & kap_face = s% min_kap_for_dPrad_dm_eqn @@ -396,8 +408,9 @@ subroutine eval_dlnPdm_qhse(s, k, & ! calculate the expected dlnPdm for HSE type(auto_diff_real_star_order1), intent(out) :: dlnPdm_qhse, Ppoint integer, intent(out) :: ierr - real(dp) :: alfa - type(auto_diff_real_star_order1) :: grav, area, P00, Pm1, inv_R2, mlt_Ptrb00, mlt_Ptrbm1 + real(dp) :: alfa, P_theta + type(auto_diff_real_star_order1) :: grav, area, P00, Pm1, inv_R2, mlt_Ptrb00, mlt_Ptrbm1, mlt_Ptrb_face + type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face include 'formats' ierr = 0 @@ -408,37 +421,67 @@ subroutine eval_dlnPdm_qhse(s, k, & ! calculate the expected dlnPdm for HSE ! for rotation, multiply gravity by factor fp. MESA 2, eqn 22. call expected_HSE_grav_term(s, k, grav, area, ierr) ! note that expected_HSE_grav_term is negative - ! mlt_pturb in thermodynamic gradients does not currently support time centering because it is timelagged. - ! replace mlt_vc check with s% mlt_vc_old(k) >0 check. - if ((s% have_mlt_vc .and. s% okay_to_set_mlt_vc) .and. s% include_mlt_Pturb_in_thermodynamic_gradients & - .and. s% mlt_Pturb_factor > 0d0) then - if (k ==1) then - mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0 - mlt_Ptrbm1 = 0d0 + if (s% using_velocity_time_centering .and. & + s% include_P_in_velocity_time_centering .and. & + s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then + P_theta = s% P_theta_for_velocity_time_centering + else + P_theta = 1d0 + end if + + if (s% use_face_values_eos_and_kap_mlt_tdc) then + if (s% have_mlt_tdc_face_state(k)) then + rho_face = s% mlt_tdc_rho_face_ad(k) + Ppoint = s% mlt_tdc_P_face_ad(k) else - mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0 - mlt_Ptrbm1 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_m1(s,k)/3d0 + call get_face_eos_kap_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face, ierr) + if (ierr /= 0) return + Ppoint = P_face + end if + if (P_theta /= 1d0) then + Ppoint = P_theta*Ppoint + (1d0 - P_theta)*s% mlt_tdc_P_face_start(k) + end if + if ((s% have_mlt_vc .and. s% okay_to_set_mlt_vc) .and. s% include_mlt_Pturb_in_thermodynamic_gradients & + .and. s% mlt_Pturb_factor > 0d0) then + ! Keep the lagged convective velocity, but form the pressure term from the same + ! face density used by the reconstructed face thermodynamic quantities. + mlt_Ptrb_face = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*rho_face/3d0 + Ppoint = Ppoint + mlt_Ptrb_face + end if + else + ! mlt_pturb in thermodynamic gradients does not currently support time centering because it is timelagged. + ! replace mlt_vc check with s% mlt_vc_old(k) >0 check. + if ((s% have_mlt_vc .and. s% okay_to_set_mlt_vc) .and. s% include_mlt_Pturb_in_thermodynamic_gradients & + .and. s% mlt_Pturb_factor > 0d0) then + if (k ==1) then + mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0 + mlt_Ptrbm1 = 0d0 + else + mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0 + mlt_Ptrbm1 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_m1(s,k)/3d0 + end if + else ! no mlt_pturb + mlt_Ptrb00 = 0d0 + mlt_Ptrbm1 = 0d0 end if - else ! no mlt_pturb - mlt_Ptrb00 = 0d0 - mlt_Ptrbm1 = 0d0 - end if - P00 = wrap_Peos_00(s,k) + P00 = wrap_Peos_00(s,k) - ! mlt Pturb doesn't support time centering yet. - if (s% using_velocity_time_centering) P00 = 0.5d0*(P00 + s% Peos_start(k)) + ! mlt Pturb doesn't support time centering yet. + if (P_theta /= 1d0) P00 = P_theta*P00 + (1d0 - P_theta)*s% Peos_start(k) - if (k == 1) then - Pm1 = 0d0 - Ppoint = P00 + mlt_Ptrb00 - else - Pm1 = wrap_Peos_m1(s,k) - if (s% using_velocity_time_centering) Pm1 = 0.5d0*(Pm1 + s% Peos_start(k-1)) ! pm1 wasn't time centered until now - Pm1 = Pm1 + mlt_Ptrbm1 ! include mlt Ptrb in k-1 - P00 = P00 + mlt_Ptrb00 ! include mlt Ptrb in k - alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) - Ppoint = alfa*P00 + (1d0-alfa)*Pm1 + if (k == 1) then + Pm1 = 0d0 + Ppoint = P00 + mlt_Ptrb00 + else + Pm1 = wrap_Peos_m1(s,k) + if (P_theta /= 1d0) Pm1 = P_theta*Pm1 + (1d0 - P_theta)*s% Peos_start(k-1) + Pm1 = Pm1 + mlt_Ptrbm1 ! include mlt Ptrb in k-1 + P00 = P00 + mlt_Ptrb00 ! include mlt Ptrb in k + alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) + Ppoint = alfa*P00 + (1d0-alfa)*Pm1 + end if end if dlnPdm_qhse = grav/(area*Ppoint) ! note that expected_HSE_grav_term is negative diff --git a/star/private/hydro_vars.f90 b/star/private/hydro_vars.f90 index 23a6d4843..270e070e3 100644 --- a/star/private/hydro_vars.f90 +++ b/star/private/hydro_vars.f90 @@ -564,6 +564,8 @@ subroutine set_hydro_vars( & if (.not. skip_mlt .and. .not. s% RSP_flag) then + s% have_mlt_tdc_face_state(1:s%nz) = .false. + if (.not. skip_mixing_info) then if (s% make_gradr_sticky_in_solver_iters) then s% fixed_gradr_for_rest_of_solver_iters(nzlo:nzhi) = .false. diff --git a/star/private/mlt_tdc_face_support.f90 b/star/private/mlt_tdc_face_support.f90 new file mode 100644 index 000000000..5fb41409f --- /dev/null +++ b/star/private/mlt_tdc_face_support.f90 @@ -0,0 +1,513 @@ +! *********************************************************************** +! +! Copyright (C) 2010-2025 The MESA Team +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License +! as published by the Free Software Foundation, +! either version 3 of the License, or (at your option) any later version. +! +! This program 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with this program. If not, see . +! +! *********************************************************************** + +module mlt_tdc_face_support + + use star_private_def + use const_def, only: dp, ln10, pi4, clight, crad + use auto_diff + use kap_support, only: get_kap + + implicit none + + private + public :: get_mlt_face_state_ad + public :: get_face_eos_kap_ad + public :: get_face_scale_height_ad + +contains + + ! Returns the MLT/TDC face thermodynamic state as + ! auto_diff_real_star_order1 quantities, either from recomputed face + ! EOS/opacity data or from the stored face quantities. + subroutine get_mlt_face_state_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, & + opacity_face, scale_height_face, gradr_face, ierr) + use star_utils, only: get_T_face, get_Peos_face, get_kap_face, get_rho_face, & + get_ChiRho_face, get_ChiT_face, get_Cp_face, get_grada_face, get_scale_height_face, get_gradr_face + + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(out) :: & + T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, & + opacity_face, scale_height_face, gradr_face + integer, intent(out) :: ierr + + ierr = 0 + if (s%use_face_values_eos_and_kap_mlt_tdc) then + call ensure_mlt_tdc_face_state_ad(s, k, ierr) + if (ierr /= 0) return + T_face = s%mlt_tdc_T_face_ad(k) + rho_face = s%mlt_tdc_rho_face_ad(k) + P_face = s%mlt_tdc_P_face_ad(k) + Cp_face = s%mlt_tdc_Cp_face_ad(k) + ChiRho_face = s%mlt_tdc_ChiRho_face_ad(k) + ChiT_face = s%mlt_tdc_ChiT_face_ad(k) + grada_face = s%mlt_tdc_grada_face_ad(k) + opacity_face = s%mlt_tdc_opacity_face_ad(k) + scale_height_face = s%mlt_tdc_scale_height_face_ad(k) + gradr_face = s%mlt_tdc_gradr_face_ad(k) + else + T_face = get_T_face(s, k) + P_face = get_Peos_face(s, k) + opacity_face = get_kap_face(s, k) + rho_face = get_rho_face(s, k) + ChiRho_face = get_ChiRho_face(s, k) + ChiT_face = get_ChiT_face(s, k) + Cp_face = get_Cp_face(s, k) + grada_face = get_grada_face(s, k) + scale_height_face = get_scale_height_face(s, k) + gradr_face = get_gradr_face(s, k) + end if + end subroutine get_mlt_face_state_ad + + + ! Ensures that the recomputed MLT and TDC face thermodynamic quantities have + ! been assembled and cached for face k. + subroutine ensure_mlt_tdc_face_state_ad(s, k, ierr) + type(star_info), pointer :: s + integer, intent(in) :: k + integer, intent(out) :: ierr + + type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, Cp_face + type(auto_diff_real_star_order1) :: ChiRho_face, ChiT_face, grada_face, opacity_face + type(auto_diff_real_star_order1) :: scale_height_face, gradr_face + + ierr = 0 + if (s%have_mlt_tdc_face_state(k)) return + + call build_mlt_tdc_face_state_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, & + opacity_face, scale_height_face, gradr_face, ierr) + if (ierr /= 0) return + + s%mlt_tdc_T_face_ad(k) = T_face + s%mlt_tdc_rho_face_ad(k) = rho_face + s%mlt_tdc_P_face_ad(k) = P_face + s%mlt_tdc_Cp_face_ad(k) = Cp_face + s%mlt_tdc_ChiRho_face_ad(k) = ChiRho_face + s%mlt_tdc_ChiT_face_ad(k) = ChiT_face + s%mlt_tdc_grada_face_ad(k) = grada_face + s%mlt_tdc_opacity_face_ad(k) = opacity_face + s%mlt_tdc_scale_height_face_ad(k) = scale_height_face + s%mlt_tdc_gradr_face_ad(k) = gradr_face + s%have_mlt_tdc_face_state(k) = .true. + end subroutine ensure_mlt_tdc_face_state_ad + + ! Reconstructs the face composition from either the current or the + ! start-of-step composition and renormalizes xa_face. + subroutine get_face_composition(s, k, use_starting_comp, zbar_face, xa_face, ierr) + use star_utils, only: get_face_weights + + type(star_info), pointer :: s + integer, intent(in) :: k + logical, intent(in) :: use_starting_comp + real(dp), intent(out) :: zbar_face + real(dp), intent(out) :: xa_face(:) + integer, intent(out) :: ierr + + real(dp) :: alfa, beta, sum_xa + + ierr = 0 + if (k == 1) then + alfa = 1d0 + beta = 0d0 + else + call get_face_weights(s, k, alfa, beta) + end if + + if (use_starting_comp) then + zbar_face = alfa*s%zbar_start(k) + beta*s%zbar_start(max(1, k-1)) + xa_face(1:s%species) = alfa*s%xa_start(1:s%species, k) + if (k > 1) xa_face(1:s%species) = xa_face(1:s%species) + beta*s%xa_start(1:s%species, k-1) + else + zbar_face = alfa*s%zbar(k) + beta*s%zbar(max(1, k-1)) + xa_face(1:s%species) = alfa*s%xa(1:s%species, k) + if (k > 1) xa_face(1:s%species) = xa_face(1:s%species) + beta*s%xa(1:s%species, k-1) + end if + + sum_xa = sum(xa_face) + if (sum_xa <= 0d0) then + ierr = -1 + if (s%report_ierr) then + !$OMP critical (mlt_tdc_face_report_ierr) + write(*,*) 'get_face_composition: sum_xa <= 0 for k', k + !$OMP end critical (mlt_tdc_face_report_ierr) + end if + return + end if + xa_face = xa_face/sum_xa + end subroutine get_face_composition + + + ! Builds the recomputed face EOS input state by wrapping T and rho to the + ! face, reconstructing the face composition, and evaluating the EOS there. + subroutine get_face_eos_inputs( & + s, k, T_face, rho_face, eos_res, d_dlnd, d_dlnT, ierr) + use eos_support, only: get_eos + use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results + use star_utils, only: get_rho_face, get_T_face + + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(out) :: T_face, rho_face + real(dp), intent(out) :: eos_res(num_eos_basic_results), d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results) + integer, intent(out) :: ierr + + real(dp) :: log10_T, log10_rho, zbar_face + real(dp) :: eos_xa_face(s%species) + real(dp) :: d_dxa(num_eos_d_dxa_results, s%species) + + ierr = 0 + T_face = get_T_face(s, k) + rho_face = get_rho_face(s, k) + if (T_face%val <= 0d0 .or. rho_face%val <= 0d0) then + ierr = -1 + if (s%report_ierr) then + !$OMP critical (mlt_tdc_face_report_ierr) + write(*,*) 'get_face_eos_inputs: bad face T or rho for k', k, T_face%val, rho_face%val + !$OMP end critical (mlt_tdc_face_report_ierr) + end if + return + end if + + call get_face_composition(s, k, .false., zbar_face, eos_xa_face, ierr) + if (ierr /= 0) return + + log10_T = log10(T_face%val) + log10_rho = log10(rho_face%val) + + call get_eos( & + s, k, eos_xa_face, rho_face%val, log10_rho, T_face%val, log10_T, & + eos_res, d_dlnd, d_dlnT, d_dxa, ierr) + if (ierr /= 0) then + if (s%report_ierr) call write_face_eos_call_info(s, k, T_face, rho_face, zbar_face, eos_xa_face) + return + end if + end subroutine get_face_eos_inputs + + + ! Interpolates extra_opacity_factor to the face and applies the existing + ! logT taper used to turn that factor on and off. + subroutine get_face_opacity_factor(s, k, log10_T, opacity_factor_face) + use star_utils, only: get_face_weights + + type(star_info), pointer :: s + integer, intent(in) :: k + real(dp), intent(in) :: log10_T + real(dp), intent(out) :: opacity_factor_face + + real(dp) :: alfa, beta + + if (k == 1) then + alfa = 1d0 + beta = 0d0 + else + call get_face_weights(s, k, alfa, beta) + end if + + opacity_factor_face = alfa*s%extra_opacity_factor(k) + if (k > 1) opacity_factor_face = opacity_factor_face + beta*s%extra_opacity_factor(k-1) + if (s%min_logT_for_opacity_factor_off > 0) then + if (log10_T >= s%max_logT_for_opacity_factor_off .or. & + log10_T <= s%min_logT_for_opacity_factor_off) then + opacity_factor_face = 1d0 + else if (log10_T > s%max_logT_for_opacity_factor_on) then + opacity_factor_face = 1d0 + (opacity_factor_face - 1d0)* & + (log10_T - s%max_logT_for_opacity_factor_off)/ & + (s%max_logT_for_opacity_factor_on - s%max_logT_for_opacity_factor_off) + else if (log10_T < s%min_logT_for_opacity_factor_on) then + opacity_factor_face = 1d0 + (opacity_factor_face - 1d0)* & + (log10_T - s%min_logT_for_opacity_factor_off)/ & + (s%min_logT_for_opacity_factor_on - s%min_logT_for_opacity_factor_off) + end if + end if + end subroutine get_face_opacity_factor + + + ! Returns the cached or newly built face EOS and opacity state as + ! auto_diff_real_star_order1 quantities for the MLT/TDC solve, + ! instead of using the stored face quantities. + subroutine get_face_eos_kap_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face, ierr) + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(out) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face + integer, intent(out) :: ierr + + ierr = 0 + call ensure_mlt_tdc_face_state_ad(s, k, ierr) + if (ierr /= 0) return + + T_face = s%mlt_tdc_T_face_ad(k) + rho_face = s%mlt_tdc_rho_face_ad(k) + P_face = s%mlt_tdc_P_face_ad(k) + Cp_face = s%mlt_tdc_Cp_face_ad(k) + ChiRho_face = s%mlt_tdc_ChiRho_face_ad(k) + ChiT_face = s%mlt_tdc_ChiT_face_ad(k) + grada_face = s%mlt_tdc_grada_face_ad(k) + opacity_face = s%mlt_tdc_opacity_face_ad(k) + end subroutine get_face_eos_kap_ad + + + ! Builds the full set of recomputed face thermodynamic quantities for the + ! MLT and TDC solve from one EOS call and one opacity call at face k. + subroutine build_mlt_tdc_face_state_ad( & + s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, & + opacity_face, scale_height_face, gradr_face, ierr) + use eos_def, only: num_eos_basic_results, i_lnPgas, i_grad_ad, i_gamma1, i_Cp, i_chiRho, i_chiT, i_eta, i_lnfree_e + use kap_def, only: num_kap_fracs + + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(out) :: & + T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, & + opacity_face, scale_height_face, gradr_face + integer, intent(out) :: ierr + + real(dp) :: log10_T, log10_rho, kap_zbar_face, opacity_factor_face + real(dp) :: eos_res(num_eos_basic_results), d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results) + real(dp) :: dlnT_face(auto_diff_star_num_vars), dlnd_face(auto_diff_star_num_vars) + real(dp) :: kap, dlnkap_dlnd, dlnkap_dlnT + real(dp) :: kap_fracs(num_kap_fracs), kap_xa_face(s%species) + type(auto_diff_real_star_order1) :: Pgas_face, gamma1_face, mlt_Pturb_ad, alpha + + ierr = 0 + call get_face_eos_inputs( & + s, k, T_face, rho_face, eos_res, d_dlnd, d_dlnT, ierr) + if (ierr /= 0) return + log10_T = log10(T_face%val) + log10_rho = log10(rho_face%val) + call set_face_log_partials(T_face, rho_face, dlnT_face, dlnd_face) + + call set_face_ad_from_log(eos_res(i_lnPgas), d_dlnd(i_lnPgas), d_dlnT(i_lnPgas), dlnd_face, dlnT_face, Pgas_face) + P_face = Pgas_face + crad*pow4(T_face)/3d0 + call set_face_ad_from_value(eos_res(i_Cp), d_dlnd(i_Cp), d_dlnT(i_Cp), dlnd_face, dlnT_face, Cp_face) + call set_face_ad_from_value(eos_res(i_chiRho), d_dlnd(i_chiRho), d_dlnT(i_chiRho), dlnd_face, dlnT_face, ChiRho_face) + call set_face_ad_from_value(eos_res(i_chiT), d_dlnd(i_chiT), d_dlnT(i_chiT), dlnd_face, dlnT_face, ChiT_face) + call set_face_ad_from_value(eos_res(i_grad_ad), d_dlnd(i_grad_ad), d_dlnT(i_grad_ad), dlnd_face, dlnT_face, grada_face) + call set_face_ad_from_value(eos_res(i_gamma1), d_dlnd(i_gamma1), d_dlnT(i_gamma1), dlnd_face, dlnT_face, gamma1_face) + + call get_face_composition(s, k, s% use_starting_composition_for_kap, kap_zbar_face, kap_xa_face, ierr) + if (ierr /= 0) return + call get_face_opacity_factor(s, k, log10_T, opacity_factor_face) + + call get_kap( & + s, k, kap_zbar_face, kap_xa_face, log10_rho, log10_T, & + eos_res(i_lnfree_e), d_dlnd(i_lnfree_e), d_dlnT(i_lnfree_e), & + eos_res(i_eta), d_dlnd(i_eta), d_dlnT(i_eta), & + kap_fracs, kap, dlnkap_dlnd, dlnkap_dlnT, ierr) + if (ierr /= 0) then + if (s%report_ierr) call write_face_kap_call_info( & + s, k, T_face, rho_face, kap_zbar_face, kap_xa_face, opacity_factor_face) + return + end if + if (is_bad_num(kap) .or. kap <= 0d0) then + ierr = -1 + if (s%report_ierr) then + !$OMP critical (mlt_tdc_face_report_ierr) + write(*,*) 'get_face_eos_kap_ad: bad face opacity for k', k, kap + !$OMP end critical (mlt_tdc_face_report_ierr) + call write_face_kap_call_info(s, k, T_face, rho_face, kap_zbar_face, kap_xa_face, opacity_factor_face) + end if + return + end if + + kap = kap*opacity_factor_face + if (s%opacity_max > 0d0 .and. kap > s%opacity_max) then + kap = s%opacity_max + dlnkap_dlnd = 0d0 + dlnkap_dlnT = 0d0 + end if + if (s%opacity_min > 0d0 .and. kap < s%opacity_min) then + kap = s%opacity_min + dlnkap_dlnd = 0d0 + dlnkap_dlnT = 0d0 + end if + call set_face_ad_from_value(kap, kap*dlnkap_dlnd, kap*dlnkap_dlnT, dlnd_face, dlnT_face, opacity_face) + + if (s% have_mlt_vc .and. s% okay_to_set_mlt_vc .and. s% include_mlt_Pturb_in_thermodynamic_gradients & + .and. s% mlt_Pturb_factor > 0d0 .and. k > 1) then + mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*rho_face/3d0 + alpha = mlt_Pturb_ad/(P_face*gamma1_face) + grada_face = grada_face*(P_face + mlt_Pturb_ad)/(P_face*(1d0 + alpha)) + end if + + call set_scale_height_from_face_state(s, k, P_face, rho_face, scale_height_face) + call set_gradr_from_face_state(s, k, P_face, opacity_face, T_face, gradr_face) + end subroutine build_mlt_tdc_face_state_ad + + + ! Returns the cached or newly built face pressure scale height as an + ! auto_diff_real_star_order1 quantity from the face EOS state. + subroutine get_face_scale_height_ad(s, k, scale_height_face, ierr) + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(out) :: scale_height_face + integer, intent(out) :: ierr + + ierr = 0 + call ensure_mlt_tdc_face_state_ad(s, k, ierr) + if (ierr /= 0) return + scale_height_face = s%mlt_tdc_scale_height_face_ad(k) + end subroutine get_face_scale_height_ad + + + subroutine set_scale_height_from_face_state(s, k, P_face, rho_face, scale_height_face) + use auto_diff_support, only: wrap_r_00 + + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(in) :: P_face, rho_face + type(auto_diff_real_star_order1), intent(out) :: scale_height_face + + real(dp) :: G + type(auto_diff_real_star_order1) :: grav, scale_height2 + + G = s%cgrav(k) + grav = G*s%m_grav(k)/pow2(wrap_r_00(s,k)) + scale_height_face = P_face/(grav*rho_face) + if (s%alt_scale_height_flag) then + scale_height2 = sqrt(P_face/G)/rho_face + if (scale_height2 < scale_height_face) scale_height_face = scale_height2 + end if + end subroutine set_scale_height_from_face_state + + + subroutine set_gradr_from_face_state(s, k, P_face, opacity_face, T_face, gradr_face) + use auto_diff_support, only: wrap_L_00 + + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(in) :: P_face, opacity_face, T_face + type(auto_diff_real_star_order1), intent(out) :: gradr_face + + real(dp) :: L_theta + type(auto_diff_real_star_order1) :: L_face, Pr_face + + if (s%include_mlt_in_velocity_time_centering) then + if (s%using_velocity_time_centering .and. & + s%include_L_in_velocity_time_centering .and. & + s%lnT(k) <= s%max_logT_for_include_P_and_L_in_velocity_time_centering*ln10) then + L_theta = s%L_theta_for_velocity_time_centering + else + L_theta = 1d0 + end if + L_face = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s%L_start(k) + else + L_face = wrap_L_00(s, k) + end if + + Pr_face = crad*pow4(T_face)/3d0 + gradr_face = P_face*opacity_face*L_face/(4d0*pi4*clight*s%m_grav(k)*s%cgrav(k)*Pr_face) + end subroutine set_gradr_from_face_state + + ! Precomputes dlnT_face and dlnd_face for converting scalar d/dlnT and + ! d/dlnd microphysics partials into star-order1 autodiff derivatives. + subroutine set_face_log_partials(T_face, rho_face, dlnT_face, dlnd_face) + type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face + real(dp), intent(out) :: dlnT_face(auto_diff_star_num_vars), dlnd_face(auto_diff_star_num_vars) + + dlnT_face = T_face%d1Array/T_face%val + dlnd_face = rho_face%d1Array/rho_face%val + end subroutine set_face_log_partials + + + ! Converts a scalar value with d/dlnd and d/dlnT partials into an + ! auto_diff_real_star_order1 quantity using the face chain rule. + subroutine set_face_ad_from_value(value, dvalue_dlnd, dvalue_dlnT, dlnd_face, dlnT_face, quantity_ad) + real(dp), intent(in) :: value, dvalue_dlnd, dvalue_dlnT + real(dp), intent(in) :: dlnd_face(auto_diff_star_num_vars), dlnT_face(auto_diff_star_num_vars) + type(auto_diff_real_star_order1), intent(out) :: quantity_ad + + quantity_ad = 0d0 + quantity_ad%val = value + quantity_ad%d1Array = dvalue_dlnd*dlnd_face + dvalue_dlnT*dlnT_face + end subroutine set_face_ad_from_value + + + ! Same as set_face_ad_from_value, but for a quantity returned in + ! logarithmic form by the microphysics routine. + subroutine set_face_ad_from_log(log_value, dlog_dlnd, dlog_dlnT, dlnd_face, dlnT_face, quantity_ad) + real(dp), intent(in) :: log_value, dlog_dlnd, dlog_dlnT + real(dp), intent(in) :: dlnd_face(auto_diff_star_num_vars), dlnT_face(auto_diff_star_num_vars) + type(auto_diff_real_star_order1), intent(out) :: quantity_ad + real(dp) :: value + + value = exp(log_value) + call set_face_ad_from_value(value, value*dlog_dlnd, value*dlog_dlnT, dlnd_face, dlnT_face, quantity_ad) + end subroutine set_face_ad_from_log + + + ! Writes the recomputed face EOS inputs that were passed to get_eos. + subroutine write_face_eos_call_info(s, k, T_face, rho_face, zbar_face, xa_face) + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face + real(dp), intent(in) :: zbar_face + real(dp), intent(in) :: xa_face(:) + + integer :: j + include 'formats' + + !$OMP critical (mlt_tdc_face_eos_call_info) + write(*,'(A)') + write(*,*) 'face EOS input info for k', k + write(*,1) 'T_face', T_face%val + write(*,1) 'rho_face', rho_face%val + write(*,1) 'log10_T_face', log10(T_face%val) + write(*,1) 'log10_rho_face', log10(rho_face%val) + write(*,1) 'zbar_face', zbar_face + write(*,1) 'sum(xa_face)', sum(xa_face) + do j = 1, s%species + write(*,2) 'xa_face ' // trim(s%nameofequ(j+s%nvar_hydro)), j, xa_face(j) + end do + !$OMP end critical (mlt_tdc_face_eos_call_info) + end subroutine write_face_eos_call_info + + + ! Writes the recomputed face opacity inputs that were passed to get_kap. + subroutine write_face_kap_call_info(s, k, T_face, rho_face, zbar_face, xa_face, opacity_factor_face) + type(star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face + real(dp), intent(in) :: zbar_face + real(dp), intent(in) :: xa_face(:) + real(dp), intent(in) :: opacity_factor_face + + integer :: j + include 'formats' + + !$OMP critical (mlt_tdc_face_kap_call_info) + write(*,'(A)') + write(*,*) 'face opacity input info for k', k + write(*,1) 'T_face', T_face%val + write(*,1) 'rho_face', rho_face%val + write(*,1) 'log10_T_face', log10(T_face%val) + write(*,1) 'log10_rho_face', log10(rho_face%val) + write(*,1) 'zbar_face', zbar_face + write(*,1) 'opacity_factor_face', opacity_factor_face + write(*,1) 'sum(xa_face)', sum(xa_face) + do j = 1, s%species + write(*,2) 'xa_face ' // trim(s%nameofequ(j+s%nvar_hydro)), j, xa_face(j) + end do + !$OMP end critical (mlt_tdc_face_kap_call_info) + end subroutine write_face_kap_call_info + +end module mlt_tdc_face_support diff --git a/star/private/phase_separation.f90 b/star/private/phase_separation.f90 index 09b25374a..9d2c8de59 100644 --- a/star/private/phase_separation.f90 +++ b/star/private/phase_separation.f90 @@ -387,6 +387,13 @@ subroutine update_model_ (s, kc_t, kc_b, do_brunt) kf_t = kc_t kf_b = kc_b + 1 + if (s% use_face_values_eos_and_kap_mlt_tdc) then + ! update_model_ changed the local composition and refreshed EOS and + ! kap on kc_t:kc_b, so the cached face thermo bundle must be rebuilt + ! for the interior faces before set_mlt_vars uses it. + s% have_mlt_tdc_face_state(kf_t+1:kf_b-1) = .false. + end if + call set_mlt_vars(s, kf_t+1, kf_b-1, ierr) if (ierr /= 0) then write(*,*) 'phase_separation: failed in call to set_mlt_vars during update_model_' diff --git a/star/private/profile_getval.f90 b/star/private/profile_getval.f90 index 471399f94..383e6a2ad 100644 --- a/star/private/profile_getval.f90 +++ b/star/private/profile_getval.f90 @@ -279,7 +279,7 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) d_dlnR00, d_dlnRp1, d_dv00, d_dvp1 integer :: j, nz, ionization_k, klo, khi, i, ii, ierr real(dp) :: f, lgT, full_on, full_off, am_nu_factor - logical :: rsp_or_w + logical :: rsp_or_w, face_cache_active include 'formats' if (s% rotation_flag) then @@ -304,6 +304,7 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) int_flag = .false. rsp_or_w = s% RSP_flag .or. s% RSP2_flag + face_cache_active = s% use_face_values_eos_and_kap_mlt_tdc if (c > extra_offset) then i = c - extra_offset @@ -1919,6 +1920,26 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) if (rsp_or_w) val = s% Lt(k) case(p_Lt_div_L) if (rsp_or_w) val = s% Lt(k)/s% L(k) + case(p_mlt_tdc_T_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_T_face_ad(k)% val + case(p_mlt_tdc_rho_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_rho_face_ad(k)% val + case(p_mlt_tdc_P_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_P_face_ad(k)% val + case(p_mlt_tdc_Cp_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_Cp_face_ad(k)% val + case(p_mlt_tdc_ChiRho_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_ChiRho_face_ad(k)% val + case(p_mlt_tdc_ChiT_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_ChiT_face_ad(k)% val + case(p_mlt_tdc_grada_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_grada_face_ad(k)% val + case(p_mlt_tdc_opacity_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_opacity_face_ad(k)% val + case(p_mlt_tdc_scale_height_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_scale_height_face_ad(k)% val + case(p_mlt_tdc_gradr_face) + if (face_cache_active .and. s% have_mlt_tdc_face_state(k)) val = s% mlt_tdc_gradr_face_ad(k)% val case(p_rsp_Et) diff --git a/star/private/star_profile_def.f90 b/star/private/star_profile_def.f90 index a51dd25c6..945c60217 100644 --- a/star/private/star_profile_def.f90 +++ b/star/private/star_profile_def.f90 @@ -599,7 +599,18 @@ module star_profile_def integer, parameter :: p_Lt = p_Lc_div_L + 1 integer, parameter :: p_Lt_div_L = p_Lt + 1 - integer, parameter :: p_rsp_log_erad = p_Lt_div_L + 1 + integer, parameter :: p_mlt_tdc_T_face = p_Lt_div_L + 1 + integer, parameter :: p_mlt_tdc_rho_face = p_mlt_tdc_T_face + 1 + integer, parameter :: p_mlt_tdc_P_face = p_mlt_tdc_rho_face + 1 + integer, parameter :: p_mlt_tdc_Cp_face = p_mlt_tdc_P_face + 1 + integer, parameter :: p_mlt_tdc_ChiRho_face = p_mlt_tdc_Cp_face + 1 + integer, parameter :: p_mlt_tdc_ChiT_face = p_mlt_tdc_ChiRho_face + 1 + integer, parameter :: p_mlt_tdc_grada_face = p_mlt_tdc_ChiT_face + 1 + integer, parameter :: p_mlt_tdc_opacity_face = p_mlt_tdc_grada_face + 1 + integer, parameter :: p_mlt_tdc_scale_height_face = p_mlt_tdc_opacity_face + 1 + integer, parameter :: p_mlt_tdc_gradr_face = p_mlt_tdc_scale_height_face + 1 + + integer, parameter :: p_rsp_log_erad = p_mlt_tdc_gradr_face + 1 integer, parameter :: p_rsp_erad = p_rsp_log_erad + 1 integer, parameter :: p_rsp_logEt = p_rsp_erad + 1 integer, parameter :: p_rsp_Et = p_rsp_logEt + 1 @@ -1277,6 +1288,16 @@ subroutine profile_column_names_init(ierr) profile_column_name(p_Lc_div_L) = 'Lc_div_L' profile_column_name(p_Lt) = 'Lt' profile_column_name(p_Lt_div_L) = 'Lt_div_L' + profile_column_name(p_mlt_tdc_T_face) = 'mlt_tdc_T_face' + profile_column_name(p_mlt_tdc_rho_face) = 'mlt_tdc_rho_face' + profile_column_name(p_mlt_tdc_P_face) = 'mlt_tdc_P_face' + profile_column_name(p_mlt_tdc_Cp_face) = 'mlt_tdc_Cp_face' + profile_column_name(p_mlt_tdc_ChiRho_face) = 'mlt_tdc_ChiRho_face' + profile_column_name(p_mlt_tdc_ChiT_face) = 'mlt_tdc_ChiT_face' + profile_column_name(p_mlt_tdc_grada_face) = 'mlt_tdc_grada_face' + profile_column_name(p_mlt_tdc_opacity_face) = 'mlt_tdc_opacity_face' + profile_column_name(p_mlt_tdc_scale_height_face) = 'mlt_tdc_scale_height_face' + profile_column_name(p_mlt_tdc_gradr_face) = 'mlt_tdc_gradr_face' profile_column_name(p_rsp_Et) = 'rsp_Et' profile_column_name(p_rsp_logEt) = 'rsp_logEt' diff --git a/star/private/struct_burn_mix.f90 b/star/private/struct_burn_mix.f90 index 387bda8d2..de0ef0b8f 100644 --- a/star/private/struct_burn_mix.f90 +++ b/star/private/struct_burn_mix.f90 @@ -268,9 +268,12 @@ end function do_rsp_step subroutine save_start_values(s, ierr) use hydro_rsp2, only: set_etrb_start_vars use star_utils, only: eval_total_energy_integrals, set_luminosity_by_category, get_Peos_face_val + use mlt_tdc_face_support, only: get_face_eos_kap_ad type (star_info), pointer :: s integer, intent(out) :: ierr integer :: k, j + type(auto_diff_real_star_order1) :: & + T_face_ad, rho_face_ad, P_face_ad, Cp_face_ad, ChiRho_face_ad, ChiT_face_ad, grada_face_ad, opacity_face_ad include 'formats' ierr = 0 @@ -298,6 +301,18 @@ subroutine save_start_values(s, ierr) s% lnPeos_start(k) = s% lnPeos(k) s% Peos_start(k) = s% Peos(k) s% Peos_face_start(k) = get_Peos_face_val(s,k) + if (s% use_face_values_eos_and_kap_mlt_tdc) then + if (s% have_mlt_tdc_face_state(k)) then + s% mlt_tdc_P_face_start(k) = s% mlt_tdc_P_face_ad(k)%val + else + call get_face_eos_kap_ad( & + s, k, T_face_ad, rho_face_ad, P_face_ad, Cp_face_ad, ChiRho_face_ad, ChiT_face_ad, grada_face_ad, opacity_face_ad, ierr) + if (ierr /= 0) return + s% mlt_tdc_P_face_start(k) = P_face_ad%val + end if + else + s% mlt_tdc_P_face_start(k) = s% Peos_face_start(k) + end if s% lnPgas_start(k) = s% lnPgas(k) s% energy_start(k) = s% energy(k) s% lnR_start(k) = s% lnR(k) @@ -1209,5 +1224,3 @@ end subroutine burn1_zone end module struct_burn_mix - - diff --git a/star/private/tdc_hydro.f90 b/star/private/tdc_hydro.f90 index 3b10b394c..cf714d4bb 100644 --- a/star/private/tdc_hydro.f90 +++ b/star/private/tdc_hydro.f90 @@ -26,6 +26,7 @@ module tdc_hydro use auto_diff_support use accurate_sum_auto_diff_star_order1 use star_utils + use mlt_tdc_face_support, only: get_face_scale_height_ad implicit none @@ -55,10 +56,17 @@ subroutine set_viscosity_vars_TDC(s, ierr) return end if - !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) + !$OMP PARALLEL DO PRIVATE(k,op_err,x) SCHEDULE(dynamic,2) do k = 1, s%nz - ! Hp_face(k) <= 0 means it needs to be set. e.g., after read file - if (s%Hp_face(k) <= 0) then + if (s%use_face_values_eos_and_kap_mlt_tdc) then + x = get_TDC_Hp_face(s, k, op_err) + if (op_err /= 0) then + !$OMP ATOMIC WRITE + ierr = op_err + else if (s%Hp_face(k) <= 0d0) then + s%Hp_face(k) = x%val + end if + else if (s%Hp_face(k) <= 0d0) then ! this scale height for face is already calculated in TDC s%Hp_face(k) = get_scale_height_face_val(s, k) ! because this is called before s% scale_height(k) is updated in mlt_vars. end if @@ -68,18 +76,27 @@ subroutine set_viscosity_vars_TDC(s, ierr) if (s%report_ierr) write (*, 2) 'failed in set_viscosity_vars_TDC loop 1', s%model_number return end if - !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) + !$OMP PARALLEL DO PRIVATE(k,op_err,x) SCHEDULE(dynamic,2) do k = 1, s%nz x = compute_Chi_div_w_face(s, k, op_err) ! Sets Chi_face - if (op_err /= 0) ierr = op_err + if (op_err /= 0) then + !$OMP ATOMIC WRITE + ierr = op_err + end if x = compute_tdc_Eq_div_w_face(s, k, op_err) ! Sets Eq_face - if (op_err /= 0) ierr = op_err + if (op_err /= 0) then + !$OMP ATOMIC WRITE + ierr = op_err + end if if (s% v_flag) then x = compute_tdc_Uq_face(s, k, op_err) else if (s% u_flag) then x = compute_tdc_Uq_dm_cell(s, k, op_err) end if - if (op_err /= 0) ierr = op_err + if (op_err /= 0) then + !$OMP ATOMIC WRITE + ierr = op_err + end if end do !$OMP END PARALLEL DO if (ierr /= 0) then @@ -104,14 +121,33 @@ subroutine get_TDC_alfa_beta_face_weights(s, k, alfa, beta) end subroutine get_TDC_alfa_beta_face_weights - function wrap_Hp_cell(s, k) result(Hp_cell) ! cm , different than rsp2 + function get_TDC_Hp_face(s, k, ierr) result(Hp_face) type(star_info), pointer :: s integer, intent(in) :: k + integer, intent(out) :: ierr + type(auto_diff_real_star_order1) :: Hp_face + + ierr = 0 + if (s%use_face_values_eos_and_kap_mlt_tdc) then + call get_face_scale_height_ad(s, k, Hp_face, ierr) + else + Hp_face = get_scale_height_face(s, k) + end if + end function get_TDC_Hp_face + + + function wrap_Hp_cell(s, k, ierr) result(Hp_cell) ! cm , different than rsp2 + type(star_info), pointer :: s + integer, intent(in) :: k + integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: Hp1, Hp0, Hp_cell - Hp0 = get_scale_height_face(s,k) + ierr = 0 + Hp0 = get_TDC_Hp_face(s, k, ierr) + if (ierr /= 0) return Hp1 = 0d0 if (k+1 < s%nz) then - Hp1 = shift_p1(get_scale_height_face(s,k+1)) + Hp1 = shift_p1(get_TDC_Hp_face(s, k+1, ierr)) + if (ierr /= 0) return end if Hp_cell = 0.5d0*(Hp0 + Hp1) !0.5d0*(wrap_Hp_00(s, k) + wrap_Hp_p1(s, k)) @@ -127,7 +163,8 @@ function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm include 'formats' ierr = 0 - Hp_cell = wrap_Hp_cell(s, k) + Hp_cell = wrap_Hp_cell(s, k, ierr) + if (ierr /= 0) return return ! below is skipped, for now. d_00 = wrap_d_00(s, k) @@ -256,7 +293,7 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then Chi_face = 0d0 else - Hp_face = get_scale_height_face(s,k) !Hp_cell_for_Chi(s, k, ierr) + Hp_face = get_TDC_Hp_face(s, k, ierr) if (ierr /= 0) return if (s%TDC_use_density_form_for_eddy_viscosity) then ! new density derivative form @@ -370,8 +407,8 @@ function compute_tdc_Uq_face(s, k, ierr) result(Uq_face) !(v_flag only) ! cm s^ else r_00 = wrap_opt_time_center_r_00(s, k) - ! which do we adopt? - Chi_00 = compute_Chi_cell(s, k, ierr) ! s% Chi_ad(k) XXX + Chi_00 = compute_Chi_cell(s, k, ierr) + if (ierr /= 0) return if (k > 1) then Chi_m1 = shift_m1(compute_Chi_cell(s, k-1, ierr)) diff --git a/star/private/turb_info.f90 b/star/private/turb_info.f90 index 52913e5e8..2ec8fe6dc 100644 --- a/star/private/turb_info.f90 +++ b/star/private/turb_info.f90 @@ -22,6 +22,7 @@ module turb_info use star_private_def use const_def, only: dp, i8, ln10, pi4, no_mixing, convective_mixing, crystallized, phase_separation_mixing + use mlt_tdc_face_support, only: get_mlt_face_state_ad use num_lib use utils_lib use auto_diff_support @@ -90,7 +91,8 @@ subroutine do1_mlt_2(s, k, & real(dp) :: crystal_pad logical :: no_mix type(auto_diff_real_star_order1) :: & - grada_face_ad, scale_height_ad, gradr_ad, rho_face_ad, & + T_face_ad, P_face_ad, opacity_face_ad, rho_face_ad, chiRho_face_ad, chiT_face_ad, Cp_face_ad, & + grada_face_ad, scale_height_ad, gradr_ad, & gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad include 'formats' @@ -118,9 +120,15 @@ subroutine do1_mlt_2(s, k, & gradL_composition_term = 0d0 end if - grada_face_ad = get_grada_face(s,k) - scale_height_ad = get_scale_height_face(s,k) - gradr_ad = get_gradr_face(s,k) + ! Assemble the full set of face thermodynamic quantities for the + ! MLT/TDC solve. + ! This helper either computes EOS and kap on faces directly or + ! recomputes them from face wrapped primitives, and returns consistent + ! T, P, rho, kap, Cp, chiRho, chiT, grada, gradr, and scale height values. + call get_mlt_face_state_ad( & + s, k, T_face_ad, rho_face_ad, P_face_ad, Cp_face_ad, chiRho_face_ad, chiT_face_ad, & + grada_face_ad, opacity_face_ad, scale_height_ad, gradr_ad, ierr) + if (ierr /= 0) return if (s% rotation_flag .and. s% mlt_use_rotation_correction) then gradr_factor = s% ft_rot(k)/s% fp_rot(k)*s% gradr_factor(k) @@ -223,6 +231,7 @@ subroutine do1_mlt_2(s, k, & end if call do1_mlt_eval(s, k, s% MLT_option, gradL_composition_term, & + T_face_ad, P_face_ad, opacity_face_ad, rho_face_ad, chiRho_face_ad, chiT_face_ad, Cp_face_ad, & gradr_ad, grada_face_ad, scale_height_ad, mixing_length_alpha, & mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, ierr) if (ierr /= 0) then diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index af4758200..66b2982d0 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -110,6 +110,7 @@ end subroutine get_gradT subroutine do1_mlt_eval( & s, k, MLT_option, gradL_composition_term, & + T_in, P_in, opacity_in, rho_in, chiRho_in, chiT_in, Cp_in, & gradr_in, grada, scale_height, mixing_length_alpha, & mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) use chem_def, only: ih1 @@ -118,7 +119,9 @@ subroutine do1_mlt_eval( & type (star_info), pointer :: s integer, intent(in) :: k character (len=*), intent(in) :: MLT_option - type(auto_diff_real_star_order1), intent(in) :: gradr_in, grada, scale_height + type(auto_diff_real_star_order1), intent(in) :: & + T_in, P_in, opacity_in, rho_in, chiRho_in, chiT_in, Cp_in, & + gradr_in, grada, scale_height real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha integer, intent(out) :: mixing_type type(auto_diff_real_star_order1), intent(out) :: & @@ -133,7 +136,7 @@ subroutine do1_mlt_eval( & ierr = 0 - P = get_Peos_face(s,k) ! if u_flag, should this be P_face_ad? (time centered in riemann) + P = P_in ! if u_flag, should this be P_face_ad? (time centered in riemann) if (s% include_mlt_in_velocity_time_centering) then ! could be cleaner with a wrapper for time_centered P and L if (s% using_velocity_time_centering .and. & @@ -152,7 +155,11 @@ subroutine do1_mlt_eval( & L_theta = 1d0 end if L = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s% L_start(k) - P = P_theta*P + (1d0-P_theta)*s% Peos_face_start(k) + if (s% use_face_values_eos_and_kap_mlt_tdc) then + P = P_theta*P + (1d0-P_theta)*s% mlt_tdc_P_face_start(k) + else + P = P_theta*P + (1d0-P_theta)*s% Peos_face_start(k) + end if r = wrap_opt_time_center_r_00(s,k) else L = wrap_L_00(s,k) @@ -161,14 +168,14 @@ subroutine do1_mlt_eval( & gradr = gradr_in cgrav = s% cgrav(k) m = s% m_grav(k) - T = get_T_face(s,k) - opacity = get_kap_face(s,k) - rho = get_rho_face(s,k) + T = T_in + opacity = opacity_in + rho = rho_in rho_start = get_rho_start_face(s,k) dV = 1d0/rho - 1d0/rho_start ! both variables are face wrapped. - chiRho = get_ChiRho_face(s,k) - chiT = get_ChiT_face(s,k) - Cp = get_Cp_face(s,k) + chiRho = chiRho_in + chiT = chiT_in + Cp = Cp_in energy = get_e_face(s,k) iso = s% dominant_iso_for_thermohaline(k) XH1 = s% xa(s% net_iso(ih1),k) diff --git a/star_data/private/star_controls_dev.inc b/star_data/private/star_controls_dev.inc index 4b8d983e1..d778d43bc 100644 --- a/star_data/private/star_controls_dev.inc +++ b/star_data/private/star_controls_dev.inc @@ -4,6 +4,7 @@ logical :: TDC_use_density_form_for_eddy_viscosity logical :: include_mlt_Pturb_in_thermodynamic_gradients logical :: use_TDC_enthalpy_flux_limiter + logical :: use_face_values_eos_and_kap_mlt_tdc logical :: include_mlt_in_velocity_time_centering logical :: use_hydro_merge_limits_in_mesh_plan diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 1991b965d..07904cb3a 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -378,6 +378,12 @@ type(auto_diff_real_star_order1), pointer, dimension(:) :: & gradT_ad, Y_face_ad, gradr_ad, mlt_vc_ad, grada_face_ad, & gradL_ad, scale_height_ad, Lambda_ad, mlt_D_ad, mlt_Gamma_ad + logical, pointer :: have_mlt_tdc_face_state(:) + type(auto_diff_real_star_order1), pointer, dimension(:) :: & + mlt_tdc_T_face_ad, mlt_tdc_rho_face_ad, mlt_tdc_P_face_ad, & + mlt_tdc_Cp_face_ad, mlt_tdc_ChiRho_face_ad, mlt_tdc_ChiT_face_ad, & + mlt_tdc_grada_face_ad, mlt_tdc_opacity_face_ad, & + mlt_tdc_scale_height_face_ad, mlt_tdc_gradr_face_ad logical, pointer :: fixed_gradr_for_rest_of_solver_iters(:) @@ -812,6 +818,7 @@ real(dp), pointer :: lnd_start(:) ! (nz) real(dp), pointer :: Peos_start(:) ! (nz) real(dp), pointer :: Peos_face_start(:) ! (nz) + real(dp), pointer :: mlt_tdc_P_face_start(:) ! (nz) real(dp), pointer :: lnPeos_start(:) ! (nz) real(dp), pointer :: lnPgas_start(:) ! (nz) real(dp), pointer :: lnT_start(:) ! (nz)