diff --git a/input/escape.cfg b/input/escape.cfg index fe0150dc..1de713ba 100644 --- a/input/escape.cfg +++ b/input/escape.cfg @@ -29,17 +29,17 @@ P_top = 1.0e-5 # bar, Pressure at TOA ##### PROTEUS settings iter_max = 9000 -log_level = INFO +log_level = DEBUG # Output subdirectory name (relative to output folder) -dir_output = earth_demo +dir_output = escape # Choose times time_star = 100.0e6 # yr, time since star formation time_target = 4.567e+9 # yr, target time for MO evolution # SOCRATES spectral file to use (relative to FWL_DATA folder) -spectral_file = spectral_files/Frostflow/256/Frostflow.sf +spectral_file = spectral_files/Frostflow/48/Frostflow.sf # Stellar heating toggle, 0: disabled | 1: enabled stellar_heating = 1 @@ -80,8 +80,7 @@ F_crit = 0.2 # Critical flux, below which the model w # Atmospheric escape escape_model = 1 # Escape model to be used, 0: None | 1: ZEPHYRUS | 2: Dummy escape_stop = 3e-4 # Terminate when atm mass drops below this fraction of its initial mass -escape_dummy_rate = 0.0 # Bulk escape rate for dummy escape model [kg s-1] -escape_el_rate = 0.15 # Efficiency factor for energy-limited escape (varies typically 0.1 < epsilon < 0.6) [dimensionless] +efficiency_factor = 0.50 # Efficiency factor for energy-limited escape (varies typically 0.1 < epsilon < 0.6) [dimensionless] escape_el_tidal_correction = 0 # Tidal correction factor, 0: None | 1: yes [dimensionless] @@ -140,7 +139,7 @@ F_atm = 8.0E4 fO2_shift_IW = 4 # Enable solving for initial partial pressures (0: off | 1: on) -solvevol_use_params = 0 +solvevol_use_params = 1 # Parameters used to solve for initial partial pressures (when solvepp_use_params = 1) Phi_global = 1.0 # Mantle melt fraction initial guess @@ -152,16 +151,16 @@ sulfur_ppmw = 200.0 # Sulfur inventory in ppmw relative to m # Prescribed injected partial pressures [bar] # Summed with solvepp results when solvepp_enabled = 1 H2O_included = 1 -H2O_initial_bar = 30.0 +H2O_initial_bar = 0.0 CO2_included = 1 -CO2_initial_bar = 30.0 +CO2_initial_bar = 0.0 N2_included = 1 -N2_initial_bar = 30.0 +N2_initial_bar = 0.0 S2_included = 1 -S2_initial_bar = 30.0 +S2_initial_bar = 0.0 SO2_included = 0 SO2_initial_bar = 0.0 diff --git a/src/proteus/utils/escape.py b/src/proteus/utils/escape.py index 7b88d8d0..f83c7915 100644 --- a/src/proteus/utils/escape.py +++ b/src/proteus/utils/escape.py @@ -11,6 +11,10 @@ log = logging.getLogger("PROTEUS") +# Define global variables +Fxuv_star_SI_full = None +age_star_mors = None + def RunDummyEsc(hf_row:dict, dt:float, phi_bulk:float): """Run dummy escape model. @@ -92,19 +96,23 @@ def RunZEPHYRUS(hf_row, dt, M_star,Omega_star,tidal_contribution, semi_major_axi out : dict Dictionary of updated total elemental mass inventories [kg] """ + global Fxuv_star_SI_full + global age_star_mors log.info("Running EL escape (ZEPHYRUS) ...") ## Step 1 : Load stellar evolution track + compute EL escape log.info("Step 1 : Load stellar evolution track + compute EL escape ") - # Get the age of the star at time t to compute XUV flux at that time - age_star = hf_row["age_star"] # [years] - star = mors.Star(Mstar=M_star, Age=age_star/1e6, Omega=Omega_star) # Load the stellar evolution track from MORS - age_star_mors = star.Tracks['Age'] # [Myrs] + # Get the age of the star at time t to compute XUV flux at that time + age_star = hf_row["age_star"] # [years] + + if (Fxuv_star_SI_full is None): + star = mors.Star(Mstar=M_star, Age=age_star/1e6, Omega=Omega_star) + age_star_mors = star.Tracks['Age'] + Fxuv_star_SI_full = ((star.Tracks['Lx'] + star.Tracks['Leuv']) / (4 * np.pi * (semi_major_axis * 1e2)**2)) * ergcm2stoWm2 # Interpolating the XUV flux at the age of the star - Fxuv_star_SI_full = ((star.Tracks['Lx'] + star.Tracks['Leuv']) / (4 * np.pi * (semi_major_axis * 1e2)**2)) * ergcm2stoWm2 Fxuv_star_SI = np.interp(age_star, age_star_mors * 1e6, Fxuv_star_SI_full) # Interpolate to get Fxuv at age_star log.info(f"Interpolated Fxuv_star_SI at age_star = {age_star} years is {Fxuv_star_SI}") diff --git a/start_proteus.py b/start_proteus.py index 27daccdd..0bf9a467 100755 --- a/start_proteus.py +++ b/start_proteus.py @@ -476,7 +476,7 @@ def main(): OPTIONS['mean_distance']*AU, OPTIONS["eccentricity"], hf_row["M_planet"], - OPTIONS["escape_el_rate"], + OPTIONS["efficiency_factor"], hf_row["R_planet"], hf_row["R_planet"]) @@ -497,10 +497,6 @@ def main(): # update total elemental inventory solvevol_target[e] = esc_m - # print info to user - log.debug(" escape %s: m=%.2e kg, dm=%.2e (%.3f%%)"% - (e, esc_m, esc_dm, 100*esc_dm/esc_m)) - # do not allow negative masses solvevol_target[e] = max(0.0, solvevol_target[e]) diff --git a/tools/GridPROTEUS.py b/tools/GridPROTEUS.py old mode 100755 new mode 100644