diff --git a/input/all_options.toml b/input/all_options.toml index 1f412c9e3..d98b803ca 100644 --- a/input/all_options.toml +++ b/input/all_options.toml @@ -24,7 +24,6 @@ # [observe] synthetic observations # ---------------------------------------------------- - version = "2.0" # Parameters @@ -100,7 +99,6 @@ version = "2.0" spin_enabled = true offset_spin = 0 # correction to calculated Breakup period [s] - # ---------------------------------------------------- # Star [star] @@ -128,7 +126,6 @@ version = "2.0" phoenix_radius = "none" # Stellar radius [R_sun] used for PHOENIX spectrum scaling phoenix_log_g = "none" # Stellar surface gravity [dex] phoenix_Teff = "none" # Stellar effective temperature [K] - [star.dummy] radius = 1.0 # Constant stellar radius [R_sun] calculate_radius = false # Calculate star radius using scaling from Teff? @@ -137,11 +134,11 @@ version = "2.0" # Orbital system [orbit] instellation_method = 'sma' # whether to define orbit using semi major axis ('sma') or instellation flux ('inst') - instellationflux = 1.0 # instellation flux received from the planet in [Earth units] - semimajoraxis = 1.0 # initial semi-major axis of planet's orbit [AU] - eccentricity = 0.0 # initial eccentricity of planet's orbit [dimensionless] - zenith_angle = 48.19 # characteristic zenith angle [degrees] - s0_factor = 0.375 # instellation scale factor [dimensionless] + instellationflux = 1.0 # instellation flux received from the planet in Earth units. + semimajoraxis = 1.0 # AU + eccentricity = 0.0 # dimensionless + zenith_angle = 48.19 # degrees + s0_factor = 0.375 # dimensionless evolve = false # whether to evolve the SMaxis and eccentricity module = "none" # module used to calculate tidal heating @@ -173,8 +170,7 @@ version = "2.0" update_dtmagma_frac = 0.03 # trigger on relative T_magma change (3%) update_dphi_abs = 0.05 # trigger on absolute Phi_global change - [struct.zalmoxis] - # Per-layer EOS: ":". See Zalmoxis docs for valid options. + [struct.zalmoxis]# Per-layer EOS: ":". See Zalmoxis docs for valid options. # Tabulated: Seager2007:iron, Seager2007:MgSiO3, WolfBower2018:MgSiO3, RTPress100TPa:MgSiO3, Seager2007:H2O # Analytic: Analytic:iron, Analytic:MgSiO3, Analytic:MgFeSiO3, Analytic:H2O, Analytic:graphite, Analytic:SiC core_eos = "Seager2007:iron" # core layer EOS @@ -212,15 +208,14 @@ version = "2.0" cloud_alpha = 0.0 # condensate retention fraction (1 -> fully retained) surf_state = "fixed" # surface scheme: "mixed_layer" | "fixed" | "skin" surf_greyalbedo = 0.1 # surface grey albedo - albedo_pl = 0.0 # Enforced input value of bond albedo; float (0 to 1) or str (path to CSV) - rayleigh = true # enable rayleigh scattering + albedo_pl = 0.0 # Bond albedo (scattering) + rayleigh = false # enable rayleigh scattering tmp_minimum = 0.5 # temperature floor on solver tmp_maximum = 5000.0 # temperature ceiling on solver module = "agni" # Which atmosphere module to use [atmos_clim.agni] - verbosity = 1 # output verbosity for agni (0:none, 1:info, 2:debug) p_top = 1.0e-5 # bar, top of atmosphere grid pressure p_obs = 0.02 # bar, level probed in transmission spectral_group = "Honeyside" # which gas opacities to include @@ -245,19 +240,15 @@ version = "2.0" sens_heat = true # include sensible heat flux near surface, with TKE scheme real_gas = false # use real-gas equations of state psurf_thresh = 0.1 # bar, surface pressure where we switch to 'transparent' mode - dx_max_ini = 300.0 # initial maximum temperature step [kelvin] allowed by solver dx_max = 35.0 # maximum temperature step [kelvin] allowed by solver max_steps = 70 # max steps allowed by solver during each iteration perturb_all = true # updated entire jacobian each step? - mlt_criterion = "s" # MLT convection stability criterion; (l)edoux or (s)chwarzschild + mlt_criterion = "l" # MLT convection stability criterion; (l)edoux or (s)chwarzschild fastchem_floor = 150.0 # Minimum temperature allowed to be sent to FC fastchem_maxiter_chem = 60000 # Maximum FC iterations (chemistry) fastchem_maxiter_solv = 20000 # Maximum FC iterations (internal solver) fastchem_xtol_chem = 1e-4 # FC solver tolerance (chemistry) fastchem_xtol_elem = 1e-4 # FC solver tolerance (elemental) - ini_profile = 'isothermal' # Initial guess for temperature profile shape - ls_default = 2 # Default linesearch method (0:none, 1:gs, 2:bt) - fdo = 2 # finite-difference order (options: 2, 4) [atmos_clim.janus] p_top = 1.0e-6 # bar, top of atmosphere grid pressure @@ -277,7 +268,7 @@ version = "2.0" [escape] module = "zephyrus" # Which escape module to use - reservoir = "outgas" # Reservoir that sets escaping composition when fractionation disabled + reservoir = "outgas" # Reservoir that sets escaping gas composition [escape.zephyrus] Pxuv = 5e-5 # Pressure at which XUV radiation become opaque in the planetary atmosphere [bar] @@ -287,26 +278,6 @@ version = "2.0" [escape.dummy] rate = 0.0 # Bulk unfractionated escape rate [kg s-1] - [escape.boreas] - fractionate = true # Include fractionation in outflow? - efficiency = 0.1 # Escape efficiency factor - sigma_H = 1.89e-18 # H absorption cross-section in XUV [cm2] - sigma_O = 2.00e-18 # O absorption ^ - sigma_C = 2.50e-18 # C absorption ^ - sigma_N = 3.00e-18 # N absorption ^ - sigma_S = 6.00e-18 # S absorption ^ - kappa_H2 = 0.01 # H2 opacity in IR, grey [cm2 g-1] - kappa_H2O = 1.0 # H2O opacity ^ - kappa_O2 = 1.0 # O2 opacity ^ - kappa_CO2 = 1.0 # CO2 opacity ^ - kappa_CO = 1.0 # CO opacity ^ - kappa_CH4 = 1.0 # CH4 opacity ^ - kappa_N2 = 1.0 # N2 opacity ^ - kappa_NH3 = 1.0 # NH3 opacity ^ - kappa_H2S = 1.0 # H2S opacity ^ - kappa_SO2 = 1.0 # SO2 opacity ^ - kappa_S2 = 1.0 # S2 opacity ^ - # Interior - physics table [interior] grain_size = 0.1 # crystal settling grain size [m] @@ -318,7 +289,6 @@ version = "2.0" melting_dir = "Monteux-600" # Melting curve set used by all interior modules (Zalmoxis, Aragog, SPIDER) eos_dir = "WolfBower2018_MgSiO3" # Dynamic EOS for SPIDER + Aragog (in FWL_DATA/interior_lookup_tables/EOS/dynamic/) - module = "spider" # Which interior module to use [interior.spider] @@ -366,6 +336,25 @@ version = "2.0" mantle_cp = 1792.0 # Mantle heat capacity [J K-1 kg-1] H_radio = 0.0 # Radiogenic heating rate [W/kg] + [interior.boundary] + rtol = 0.0001 # ODE solver relative tolerance. + atol = 1e-08 # ODE solver absolute tolerance. + T_surf_0 = 3900.0 # Initial surface temperature [K] + T_p_0 = 4000.0 # Initial mantle potential temperature [K] + T_solidus = 1420.0 # Mantle solidus temperature [K] + T_liquidus = 2020.0 # Mantle liquidus temperature [K] + critical_melt_fraction = 0.4 # Critical melt fraction for rheological transition. + critical_rayleigh_number = 1100.0 # Critical Rayleigh number for onset of convection + heat_fusion_silicate = 400000.0 # Latent heat of fusion for silicates [J/kg] + nusselt_exponent = 0.33 # Nusselt-Rayleigh scaling exponent + silicate_heat_capacity = 1200.0 # Silicate heat capacity [J/kg/K] + thermal_conductivity = 4.2 # Thermal conductivity [W/m/K] + thermal_diffusivity = 1e-06 # Thermal diffusivity [m^2/s] + thermal_expansivity = 2e-05 # Thermal expansivity [1/K] + transition_width = 0.15 # Width of viscosity transition in melt fraction space + eta_solid_const = 1e+21 # Constant solid viscosity for aggregate formulation [Pa s] + eta_melt_const = 100.0 # Constant melt viscosity for aggregate formulation [Pa s] + # Outgassing - physics table [outgas] fO2_shift_IW = 4 # atmosphere/interior boundary oxidation state [log10(ΔIW)] @@ -461,7 +450,7 @@ version = "2.0" ini_mix = "profile" # Initial mixing ratios (profile, outgas) fix_surf = false # Fixed surface mixing ratios network = "SNCHO" # Class of chemical network to use (CHO, NCHO, SNCHO) - save_frames = false # Plot frames during iterations + save_frames = true # Plot frames during iterations yconv_cri = 0.05 # Convergence criterion, value of mixing ratios slope_cri = 0.0001 # Convergence criterion, rate of change of mixing ratios diff --git a/input/chili/intercomp/_base.grid.toml b/input/chili/intercomp/_base.grid.toml index dc7f06383..a81279d9d 100644 --- a/input/chili/intercomp/_base.grid.toml +++ b/input/chili/intercomp/_base.grid.toml @@ -3,13 +3,13 @@ symlink = "" ref_config = "input/chili/intercomp/_base.toml" use_slurm = true -max_jobs = 9 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_jobs = 16 # maximum number of concurrent tasks (e.g. 500 on Habrok) max_days = 1 # maximum number of days to run (e.g. 1) max_mem = 3 # maximum memory per CPU in GB (e.g. 3) ["delivery.elements.H_kg"] method = "direct" - values = [1.60e20, 7.80e20, 16e20] + values = [1.60e20, 4.70e20, 7.80e20, 16e20] ["delivery.elements.C_kg"] method = "direct" diff --git a/input/chili/intercomp/_base.toml b/input/chili/intercomp/_base.toml index ef195be23..0fcacc2c5 100644 --- a/input/chili/intercomp/_base.toml +++ b/input/chili/intercomp/_base.toml @@ -12,7 +12,7 @@ archive_mod = "none" remove_sf = true [params.dt] -starspec = 1e9 +starspec = 10000000.0 starinst = 10.0 method = "adaptive" minimum = 100.0 @@ -64,11 +64,10 @@ age_ini = 0.05 [star.mors] age_now = 4.567 +spec = "stellar_spectra/Named/sun.txt" rot_pcntle = 50.0 rot_period = "none" tracks = "spada" -star_name = "sun" -spectrum_source = "solar" [star.dummy] Teff = 5772.0 @@ -99,7 +98,7 @@ radius_int = "none" [atmos_clim] module = "agni" -surf_state = "skin" +surf_state = "fixed" prevent_warming = true surface_d = 0.01 surface_k = 2.0 @@ -114,22 +113,20 @@ tmp_maximum = 5000.0 [atmos_clim.agni] spectral_group = "Dayspring" spectral_bands = "48" -p_top = 1e-04 +p_top = 1e-05 p_obs = 0.02 -surf_material = "surface_albedos/Hammond24/lunarmarebasalt.dat" -num_levels = 45 -chemistry = "eq" +surf_material = "greybody" +num_levels = 42 +chemistry = "none" solve_energy = true -solution_atol = 0.08 +solution_atol = 0.1 solution_rtol = 0.1 overlap_method = "ee" -rainout = false +condensation = false latent_heat = false -real_gas = true +real_gas = false psurf_thresh = 0.1 -dx_max = 15.0 -dx_max_ini = 200.0 -ls_default = 1 +dx_max = 10.0 max_steps = 70 perturb_all = true mlt_criterion = "s" @@ -159,9 +156,9 @@ rheo_phi_loc = 0.5 rheo_phi_wid = 0.2 [interior.spider] -ini_entropy = 3900.0 +ini_entropy = 3000.0 ini_dsdr = -4.698e-06 -num_levels = 150 +num_levels = 140 mixing_length = 2 tolerance = 1e-08 tolerance_rel = 1e-08 diff --git a/input/chili/protocol/earth.toml b/input/chili/protocol/earth.toml index b8d476b1e..5641b6c2b 100644 --- a/input/chili/protocol/earth.toml +++ b/input/chili/protocol/earth.toml @@ -69,11 +69,15 @@ age_ini = 0.05 [star.mors] age_now = 4.567 +spec = "stellar_spectra/Named/sun.txt" rot_pcntle = 50.0 rot_period = "none" tracks = "spada" -star_name = "sun" -spectrum_source = "solar" + +[star.dummy] +Teff = 5772.0 +radius = 1.0 +calculate_radius = false [orbit] module = "none" @@ -89,6 +93,15 @@ semimajoraxis_sat = 300000000.0 instellation_method = "sma" instellationflux = 1.0 +[orbit.dummy] +H_tide = 1e-07 +Phi_tide = "<0.3" +Imk2 = 0.0 + +[orbit.lovepy] +visc_thresh = 1000000000.0 +ncalc = 1000 + [struct] corefrac = 0.55 module = "self" @@ -124,7 +137,7 @@ solve_energy = true solution_atol = 0.1 solution_rtol = 0.1 overlap_method = "ee" -rainout = false +condensation = false latent_heat = false real_gas = false psurf_thresh = 0.1 @@ -134,6 +147,20 @@ perturb_all = true mlt_criterion = "l" fastchem_floor = 100.0 +[atmos_clim.janus] +spectral_group = "none" +spectral_bands = "none" +p_top = 1e-05 +p_obs = 0.002 +F_atm_bc = 0 +num_levels = 90 +tropopause = "none" +overlap_method = "ee" + +[atmos_clim.dummy] +gamma = 0.01 +height_factor = 3.0 + [atmos_chem] module = "vulcan" when = "manually" @@ -163,6 +190,9 @@ Pxuv = 5e-05 efficiency = 0.1 tidal = true +[escape.dummy] +rate = 0.0 + [interior] module = "spider" melting_dir = "Monteux-600" @@ -184,6 +214,36 @@ solver_type = "bdf" tsurf_atol = 10.0 tsurf_rtol = 0.02 +[interior.aragog] +logging = "ERROR" +ini_tmagma = "none" +basal_temperature = 7000 +init_file = "none" +num_levels = 100 +initial_condition = 1 +tolerance = 1e-10 +inner_boundary_condition = 1 +inner_boundary_value = 4000 +conduction = true +convection = true +gravitational_separation = false +mixing = false +dilatation = false +mass_coordinates = false +tsurf_poststep_change = 30 +event_triggering = true +bulk_modulus = 260000000000.0 + +[interior.dummy] +ini_tmagma = 3300.0 +tmagma_atol = 30.0 +tmagma_rtol = 0.05 +mantle_tliq = 2700.0 +mantle_tsol = 1700.0 +mantle_rho = 4550.0 +mantle_cp = 1792.0 +H_radio = 0.0 + [outgas] fO2_shift_IW = 4.0 module = "calliope" @@ -205,6 +265,8 @@ rtol = 0.0001 xtol = 1e-06 solubility = true +[outgas.atmodeller] +some_parameter = "some_value" [delivery] module = "none" @@ -230,5 +292,21 @@ SH_ratio = 0.0 S_kg = 0.0 S_ppmw = 0.0 +[delivery.volatiles] +H2O = 0 +CO2 = 0 +N2 = 0 +S2 = 0 +SO2 = 0 +H2S = 0 +NH3 = 0 +H2 = 0 +CH4 = 0 +CO = 0 + [observe] synthesis = "none" + +[observe.platon] +downsample = 8 +clip_vmr = 1e-08 diff --git a/input/chili/protocol/tr1b.toml b/input/chili/protocol/tr1b.toml index ebd393e88..e58551456 100644 --- a/input/chili/protocol/tr1b.toml +++ b/input/chili/protocol/tr1b.toml @@ -69,11 +69,15 @@ age_ini = 0.05 [star.mors] age_now = 7.6 +spec = "stellar_spectra/Named/trappist-1.txt" rot_pcntle = 50.0 rot_period = "none" tracks = "spada" -star_name = "trappist-1" -spectrum_source = "muscles" + +[star.dummy] +Teff = 5772.0 +radius = 1.0 +calculate_radius = false [orbit] module = "none" @@ -89,6 +93,15 @@ semimajoraxis_sat = 300000000.0 instellation_method = "sma" instellationflux = 1.0 +[orbit.dummy] +H_tide = 1e-07 +Phi_tide = "<0.3" +Imk2 = 0.0 + +[orbit.lovepy] +visc_thresh = 1000000000.0 +ncalc = 1000 + [struct] corefrac = 0.55 module = "self" @@ -124,7 +137,7 @@ solve_energy = true solution_atol = 0.1 solution_rtol = 0.1 overlap_method = "ee" -rainout = false +condensation = false latent_heat = false real_gas = false psurf_thresh = 0.1 @@ -134,6 +147,19 @@ perturb_all = true mlt_criterion = "l" fastchem_floor = 100.0 +[atmos_clim.janus] +spectral_group = "none" +spectral_bands = "none" +p_top = 1e-05 +p_obs = 0.002 +F_atm_bc = 0 +num_levels = 90 +tropopause = "none" +overlap_method = "ee" + +[atmos_clim.dummy] +gamma = 0.01 +height_factor = 3.0 [atmos_chem] module = "vulcan" @@ -164,6 +190,9 @@ Pxuv = 5e-05 efficiency = 0.1 tidal = true +[escape.dummy] +rate = 0.0 + [interior] module = "spider" melting_dir = "Monteux-600" @@ -185,6 +214,36 @@ solver_type = "bdf" tsurf_atol = 10.0 tsurf_rtol = 0.02 +[interior.aragog] +logging = "ERROR" +ini_tmagma = "none" +basal_temperature = 7000 +init_file = "none" +num_levels = 100 +initial_condition = 1 +tolerance = 1e-10 +inner_boundary_condition = 1 +inner_boundary_value = 4000 +conduction = true +convection = true +gravitational_separation = false +mixing = false +dilatation = false +mass_coordinates = false +tsurf_poststep_change = 30 +event_triggering = true +bulk_modulus = 260000000000.0 + +[interior.dummy] +ini_tmagma = 3300.0 +tmagma_atol = 30.0 +tmagma_rtol = 0.05 +mantle_tliq = 2700.0 +mantle_tsol = 1700.0 +mantle_rho = 4550.0 +mantle_cp = 1792.0 +H_radio = 0.0 + [outgas] fO2_shift_IW = 4.0 module = "calliope" @@ -206,6 +265,9 @@ rtol = 0.0001 xtol = 1e-06 solubility = true +[outgas.atmodeller] +some_parameter = "some_value" + [delivery] module = "none" initial = "elements" @@ -230,5 +292,21 @@ SH_ratio = 0.0 S_kg = 0.0 S_ppmw = 0.0 +[delivery.volatiles] +H2O = 0 +CO2 = 0 +N2 = 0 +S2 = 0 +SO2 = 0 +H2S = 0 +NH3 = 0 +H2 = 0 +CH4 = 0 +CO = 0 + [observe] synthesis = "none" + +[observe.platon] +downsample = 8 +clip_vmr = 1e-08 diff --git a/input/chili/readme.md b/input/chili/readme.md index fe642fc96..8466c4452 100644 --- a/input/chili/readme.md +++ b/input/chili/readme.md @@ -1,6 +1,5 @@ Configuration files for the CHILI model intercomparison project. Website: https://nexss.info/cuisines/chili-mip/ -GitHub for CHILI: https://github.com/projectcuisines/chili The `protocol/` subfolder contains the two configurations used in the protocol paper. diff --git a/input/demos/aragog.toml b/input/demos/aragog.toml index 8a58fcef7..cd8955c1a 100644 --- a/input/demos/aragog.toml +++ b/input/demos/aragog.toml @@ -71,7 +71,7 @@ version = "2.0" solution_atol = 1e-3 # solver absolute tolerance solution_rtol = 2e-2 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = true # use real-gas equations of state [atmos_clim.dummy] diff --git a/input/demos/escape.toml b/input/demos/escape.toml index 7307b54a7..27178e1d5 100644 --- a/input/demos/escape.toml +++ b/input/demos/escape.toml @@ -57,8 +57,7 @@ version = "2.0" rot_pcntle = 50.0 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 4.567 # Gyr, current age of star used for scaling - star_name = "sun" - spectrum_source = "solar" + spec = "stellar_spectra/Named/sun.txt" # stellar spectrum [orbit] semimajoraxis = 0.5 # AU @@ -96,7 +95,7 @@ version = "2.0" solution_atol = 2.0 solution_rtol = 0.13 overlap_method = "ee" - rainout = false + condensation = false real_gas = true [escape] diff --git a/input/demos/sat_physical.toml b/input/demos/sat_physical.toml index 7b92031f8..6b20ff818 100644 --- a/input/demos/sat_physical.toml +++ b/input/demos/sat_physical.toml @@ -125,7 +125,7 @@ version = "2.0" solution_atol = 1e-3 # solver absolute tolerance solution_rtol = 2e-2 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = false # use real-gas equations of state [escape] diff --git a/input/ensembles/GDSP_initial_grid_2.toml b/input/ensembles/GDSP_initial_grid_2.toml new file mode 100644 index 000000000..f05486fca --- /dev/null +++ b/input/ensembles/GDSP_initial_grid_2.toml @@ -0,0 +1,36 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "GDPS_initial_grid_2/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/GDSP_fiducial.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 25 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +["delivery.elements.H_ppmw"] + method = "logspace" + start = 709 + stop = 7090 + count = 5 + +["orbit.instellationflux"] + method = "logspace" + start = 1 + stop = 100 + count = 10 diff --git a/input/ensembles/SNS_0.5%_sanity_check.toml b/input/ensembles/SNS_0.5%_sanity_check.toml new file mode 100644 index 000000000..2dbb34640 --- /dev/null +++ b/input/ensembles/SNS_0.5%_sanity_check.toml @@ -0,0 +1,34 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "SNS_0.5%_sanity_check/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/fiducial_sub_Neptune.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 80 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +["delivery.elements.H_ppmw"] + method = "direct" + values = [3545] + +["orbit.instellationflux"] + method = "logspace" + start = 1 + stop = 50 + count = 30 diff --git a/input/ensembles/SNS_0.5%_sweep.toml b/input/ensembles/SNS_0.5%_sweep.toml new file mode 100644 index 000000000..afda50709 --- /dev/null +++ b/input/ensembles/SNS_0.5%_sweep.toml @@ -0,0 +1,40 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "SNS_0.5%_sweep_test/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/fiducial_sub_Neptune.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 80 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +["delivery.elements.H_ppmw"] + method = "direct" + values = [2128] + +["star.dummy.Teff"] + method = "linspace" + start = 4900 + stop = 6300 + count = 8 + +["orbit.instellationflux"] + method = "linspace" + start = 11 + stop = 12 + count = 10 diff --git a/input/ensembles/boundary_interior_test.toml b/input/ensembles/boundary_interior_test.toml new file mode 100644 index 000000000..a5863f030 --- /dev/null +++ b/input/ensembles/boundary_interior_test.toml @@ -0,0 +1,31 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "boundary_interior_test/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/nogit_minimal_boundary_dummy_atm.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 10 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +# Planet mass set directly +["atmos_clim.dummy.gamma"] + method = "linspace" + start = 0.1 + stop = 1.0 + count = 10 diff --git a/input/ensembles/example.infer.toml b/input/ensembles/example.infer.toml index 9ae566e40..623f0260c 100644 --- a/input/ensembles/example.infer.toml +++ b/input/ensembles/example.infer.toml @@ -1,34 +1,36 @@ # Config file for inference scheme # Set seed for reproducibility -seed = 1 +seed = 2 # Path to output folder where inference will be saved (relative to PROTEUS output folder) -output = "infer_toy/BO/" +output = "infer_demo/" # Path to base (reference) config file relative to PROTEUS root folder ref_config = "input/demos/dummy.toml" # Method for initialising the inference scheme (one of these must be 'none') -init_samps = 10 # Number of random samples if starting from scratch. +init_samps = 2 # Number of random samples if starting from scratch. init_grid = 'none' # grid_demo/' # Path pre-computed grid (relative to PROTEUS output folder) # Parameters for Bayesian optimisation -n_workers = 6 # Number of parallel workers -kernel = "MAT3/2" # Kernel type for GP, "RBF" | "MAT1/2" | "MAT3/2" | "MAT5/2" -acqf = "LogEI" # Acquisition function, "UCB" | "LogEI" | "E-LogEI" -n_steps = 50 # Total number of evaluations (i.e. BO steps) +n_workers = 7 # Number of parallel workers +kernel = "MAT" # Kernel type for GP, "RBF" | "MAT" +acqf = "LogEI" # Acquisition function, "UCB" | "LogEI" +n_steps = 30 # Total number of evaluations (i.e. BO steps) +n_restarts = 10 # GP optimization restarts +n_samples = 1000 # Raw samples for acquisition optimization # Parameters to optimize (with bounds) [parameters] "struct.mass_tot" = [0.7, 3.0] "struct.corefrac" = [0.3, 0.9] -# "delivery.elements.H_ppmw" = [1e3, 2e4] -# "outgas.fO2_shift_IW" = [-3.0, 5.0] -# "escape.dummy.rate" = [1e5, 1e7] +"delivery.elements.H_ppmw" = [1e3, 2e4] +"outgas.fO2_shift_IW" = [-3.0, 5.0] +"escape.dummy.rate" = [1e5, 1e7] # Target observables to match by optimisation [observables] -"R_obs" = 7965615.42040000017732382 -# "H2O_vmr" = 0.69457002234999998 -"P_surf" = 16333.19793399999980466 +"R_obs" = 6.4781639626e+6 +"H2O_vmr" = 7.5032517428e-1 +"P_surf" = 3.8198606763e+3 diff --git a/input/ensembles/ocean_formation_grid.toml b/input/ensembles/ocean_formation_grid.toml new file mode 100644 index 000000000..4bca045cb --- /dev/null +++ b/input/ensembles/ocean_formation_grid.toml @@ -0,0 +1,31 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "GDPS_initial_grid/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/GDSP_fiducial.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 20 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +# Hydrogen inventory set by arange +["delivery.elements.H_oceans"] + method = "arange" + start = 1 + stop = 20 + step = 1 diff --git a/input/ensembles/solubility_test.toml b/input/ensembles/solubility_test.toml new file mode 100644 index 000000000..251a8978e --- /dev/null +++ b/input/ensembles/solubility_test.toml @@ -0,0 +1,22 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "SNS_solubility_test/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/fiducial_sub_Neptune.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 80 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +["outgas.calliope.solubility"] + method = "direct" + values = [true,false] diff --git a/input/ensembles/stellar_spectra_test_bb.toml b/input/ensembles/stellar_spectra_test_bb.toml new file mode 100644 index 000000000..14a2bd03a --- /dev/null +++ b/input/ensembles/stellar_spectra_test_bb.toml @@ -0,0 +1,29 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "stellar_spectra_test_bb/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/nogit_fiducial_sub_Neptune.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 10 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +# Planet mass set directly +["star.dummy.Teff"] + method = "direct" + values = [3467,5212,5772] diff --git a/input/ensembles/stellar_spectra_test_mors.toml b/input/ensembles/stellar_spectra_test_mors.toml new file mode 100644 index 000000000..3cb66dc17 --- /dev/null +++ b/input/ensembles/stellar_spectra_test_mors.toml @@ -0,0 +1,29 @@ +# Config file for running a grid of forward models + +# Path to output folder where grid will be saved (relative to PROTEUS output folder) +output = "stellar_spectra_test_mors/" + +# Make `output` a symbolic link to this absolute location. To disable: set to empty string. +symlink = "" + +# Path to base (reference) config file relative to PROTEUS root folder +ref_config = "input/planets/nogit_fiducial_sub_Neptune_MORS.toml" + +# Use SLURM? +use_slurm = false + +# Execution limits +max_jobs = 10 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) +max_mem = 3 # maximum memory per CPU in GB (e.g. 3) + +# Now define grid axes... +# Each axis must be a new section (table) in this file. +# Each table corresponds to the name of the parameter to be varied. +# Each table name must be written in double quotes. +# See examples below + +# Planet mass set directly +["star.mors.spec"] + method = "direct" + values = ["stellar_spectra/Named/gj849.txt","stellar_spectra/Named/hd97658.txt","stellar_spectra/Named/sun.txt"] diff --git a/input/minimal.toml b/input/minimal.toml index b3137bd1e..973142e33 100644 --- a/input/minimal.toml +++ b/input/minimal.toml @@ -19,8 +19,7 @@ version = "2.0" [star.mors] rot_pcntle = 50.0 # Rotation percentile age_now = 4.567 # Current age of star used for scaling [Gyr] - star_name = "sun" - spectrum_source = "solar" + spec = "stellar_spectra/Named/sun.txt" # Stellar spectrum [orbit] semimajoraxis = 1.0 # [AU] diff --git a/input/planets/GDSP_fiducial.toml b/input/planets/GDSP_fiducial.toml new file mode 100644 index 000000000..3499a423b --- /dev/null +++ b/input/planets/GDSP_fiducial.toml @@ -0,0 +1,281 @@ +# PROTEUS configuration file (version 2.0) + +# Root tables should be physical, with the exception of "params" +# Software related options should go within the appropriate physical table + +# The general structure is: +# [root] metadata +# [params] parameters for code execution, output files, time-stepping, convergence +# [star] stellar parameters, model selection +# [orbit] planetary orbital parameters +# [struct] planetary structure (mass, radius) +# [atmos] atmosphere parameters, model selection +# [escape] escape parameters, model selection +# [interior] magma ocean model selection and parameters +# [outgas] outgassing parameters (fO2) and included volatiles +# [delivery] initial volatile inventory, and delivery model selection + +# ---------------------------------------------------- +# Metadata +version = "2.0" + +# ---------------------------------------------------- +# Parameters +[params] + # output files + [params.out] + path = "GDPS_fiducial" + logging = "DEBUG" + plot_mod = 1 # Plotting frequency, 0: wait until completion | n: every n iterations + plot_fmt = "png" # Plotting image file format, "png" or "pdf" recommended + write_mod = 1 # Write CSV frequency, 0: wait until completion | n: every n iterations + + # time-stepping + [params.dt] + minimum = 3e2 # yr, minimum time-step + maximum = 1e7 # yr, maximum time-step + initial = 1e3 # yr, inital step size + starspec = 3e8 # yr, interval to re-calculate the stellar spectrum + starinst = 1e2 # yr, interval to re-calculate the instellation + method = "adaptive" # proportional | adaptive | maximum + + [params.dt.proportional] + propconst = 52.0 # Proportionality constant + + [params.dt.adaptive] + atol = 0.03 # Step size atol + rtol = 0.15 # Step size rtol + + # Termination criteria + # Set enabled=true/false in each section to enable/disable that termination criterion + [params.stop] + + # Require criteria to be satisfied twice before model will exit? + strict = false + + # required number of iterations + [params.stop.iters] + enabled = true + minimum = 5 + maximum = 9000 + + # required time constraints + [params.stop.time] + enabled = true + minimum = 1.0e3 # yr, model will certainly run to t > minimum + maximum = 2.0e+10 # yr, model will terminate when t > maximum + + # solidification + [params.stop.solid] + enabled = false + phi_crit = 0.01 # non-dim., model will terminate when global melt fraction < phi_crit + + # radiative equilibrium + [params.stop.radeqm] + enabled = true + atol = 0.1 # absolute tolerance [W m-2] + rtol = 1e-3 # relative tolerance + + [params.stop.escape] + enabled = false + mass_frac = 3e-4 # Stop when atm_mass < this frac of initial mass + + +# ---------------------------------------------------- +# Star +[star] + + # Physical parameters + mass = 0.465 # M_sun + age_ini = 0.100 # Gyr, model initialisation/start age + + module = "dummy" + [star.mors] + rot_pctle = 50.0 # rotation percentile + tracks = "spada" # evolution tracks: spada | baraffe + age_now = 2.4 # Gyr, current age of star used for scaling + spec = "stellar_spectra/Named/gj849.txt" # stellar spectrum + + [star.dummy] + calculate_radius = true # calculate radius using empirical mass-luminosity and mass-radius relations + Teff = 4500.0 # K + +# Orbital system +[orbit] + instellation_method = 'inst' # whether to define orbit using semi major axis ('sma') or instellation flux ('inst') + instellationflux = 1 # instellation flux received from the planet in Earth units. + semimajoraxis = 0.15 # AU + eccentricity = 0.0 # dimensionless + zenith_angle = 54.74 # degrees + s0_factor = 0.25 # dimensionless + + module = "none" + + [orbit.dummy] + H_tide = 1e-11 # Fixed tidal power density [W kg-1] + Phi_tide = "<0.3" # Tidal heating applied when inequality locally satisfied + + [orbit.lovepy] + visc_thresh = 1e9 # Minimum viscosity required for heating [Pa s] + +# Planetary structure - physics table +[struct] + set_by = 'mass_tot' # Variable to set interior structure: 'radius_int' or 'mass_tot' + mass_tot = 5.0 # M_earth + corefrac = 0.55 # non-dim., radius fraction + core_density = 10738.33 # Core density [kg m-3] + core_heatcap = 880.0 # Core specific heat capacity [J K-1 kg-1] + +# Atmosphere - physics table +[atmos_clim] + prevent_warming = true # do not allow the planet to heat up + surface_d = 0.01 # m, conductive skin thickness + surface_k = 2.0 # W m-1 K-1, conductive skin thermal conductivity + cloud_enabled = false # enable water cloud radiative effects + cloud_alpha = 0.0 # condensate retention fraction (1 -> fully retained) + surf_state = "skin" # surface scheme: "mixed_layer" | "fixed" | "skin" + surf_greyalbedo = 0.0 # surface grey albedo + albedo_pl = 0.0 # Enforced Bond albedo (do not use with `rayleigh = true`) + rayleigh = true # Enable rayleigh scattering + tmp_minimum = 0.5 # temperature floor on solver + tmp_maximum = 5000.0 # temperature ceiling on solver + + module = "agni" # Which atmosphere module to use + + [atmos_clim.agni] + p_top = 1.0e-5 # bar, top of atmosphere grid pressure + spectral_group = "Dayspring" # which gas opacities to include + spectral_bands = "48" # how many spectral bands? + num_levels = 30 # Number of atmospheric grid levels + chemistry = "none" # "none" | "eq" + surf_material = "greybody" # surface material file for scattering + solve_energy = true # solve for energy-conserving atmosphere profile + solution_atol = 1e-2 # solver absolute tolerance + solution_rtol = 1e-1 # solver relative tolerance + overlap_method = "ee" # gas overlap method + condensation = true # volatile condensation + real_gas = true # use real-gas equations of state + + [atmos_clim.janus] + p_top = 1.0e-5 # bar, top of atmosphere grid pressure + p_obs = 1.0e-3 # bar, observed pressure level + spectral_group = "Dayspring" # which gas opacities to include + spectral_bands = "48" # how many spectral bands? + F_atm_bc = 0 # measure outgoing flux at: (0) TOA | (1) Surface + num_levels = 90 # Number of atmospheric grid levels + tropopause = "none" # none | skin | dynamic + overlap_method = "ee" # gas overlap method + + [atmos_clim.dummy] + gamma = 0.3 # atmosphere opacity between 0 and 1 + +# Volatile escape - physics table +[escape] + + module = "none" # Which escape module to use + + [escape.zephyrus] + Pxuv = 1e-2 # Pressure at which XUV radiation become opaque in the planetary atmosphere [bar] + efficiency = 1.0 # Escape efficiency factor + tidal = false # Tidal contribution enabled + + [escape.dummy] + rate = 0.0 # Bulk unfractionated escape rate [kg s-1] + +# Interior - physics table +[interior] + grain_size = 0.1 # crystal settling grain size [m] + F_initial = 1e5 # Initial heat flux guess [W m-2] + radiogenic_heat = true # enable radiogenic heat production + tidal_heat = false # enable tidal heat production + rheo_phi_loc = 0.4 # Centre of rheological transition + rheo_phi_wid = 0.15 # Width of rheological transition + + module = "spider" # Which interior module to use + + [interior.spider] + num_levels = 110 # Number of SPIDER grid levels + mixing_length = 2 # Mixing length parameterization + tolerance = 1.0e-7 # solver tolerance + tsurf_atol = 20.0 # tsurf_poststep_change + tsurf_rtol = 0.01 # tsurf_poststep_change_frac + ini_entropy = 3000.0 # Surface entropy conditions [J K-1 kg-1] + ini_dsdr = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1] + + [interior.aragog] + num_levels = 110 # Number of Aragog grid levels + tolerance = 1.0e-7 # solver tolerance + ini_tmagma = 3000.0 # Initial magma surface temperature [K] + + [interior.dummy] + ini_tmagma = 3500.0 # Initial magma surface temperature [K] + +# Outgassing - physics table +[outgas] + fO2_shift_IW = -3 # log10(ΔIW), atmosphere/interior boundary oxidation state + + module = "calliope" # Which outgassing module to use + + [outgas.calliope] + include_H2O = true # Include H2O compound + include_CO2 = true # Include CO2 compound + include_N2 = true # Include N2 compound + include_S2 = true # Include S2 compound + include_SO2 = true # Include SO2 compound + include_H2 = true # Include H2 compound + include_CH4 = true # Include CH4 compound + include_CO = true # Include CO compound + T_floor = 700.0 # Temperature floor applied to outgassing calculation [K]. + + [outgas.atmodeller] + some_parameter = "some_value" + +# Volatile delivery - physics table +[delivery] + + # Radionuclide parameters + radio_tref = 4.55 # Reference age for concentrations [Gyr] + radio_K = 310.0 # ppmw of potassium (all isotopes) + radio_U = 0.031 # ppmw of uranium (all isotopes) + radio_Th = 0.124 # ppmw of thorium (all isotopes) + + # Which initial inventory to use? + initial = 'elements' # "elements" | "volatiles" + + # No module for accretion as of yet + module = "none" + + # Set initial volatile inventory by planetary element abundances + [delivery.elements] + use_metallicity = false # whether or not to specify the elemental abundances in terms of solar metallicity + metallicity = 1000 # metallicity relative to solar metallicity + + # H_oceans = 1 # Hydrogen inventory in units of equivalent Earth oceans + H_ppmw = 709.0 # Hydrogen inventory in ppmw relative to mantle mass + + CH_ratio = 0.32 # C/H mass ratio in mantle/atmosphere system + # C_ppmw = 10.0 # Carbon inventory in ppmw relative to mantle mass + + NH_ratio = 0.09 # N/H mass ratio in mantle/atmosphere system + # N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass + + SH_ratio = 0.0 # S/H mass ratio in mantle/atmosphere system + # S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mantle mass + + + # Set initial volatile inventory by partial pressures in atmosphere + [delivery.volatiles] + H2O = 0.0 # partial pressure of H2O + CO2 = 0.0 # partial pressure of CO2 + N2 = 0.0 # etc + S2 = 0.0 + SO2 = 0.0 + H2 = 0.0 + CH4 = 0.0 + CO = 0.0 + +[observe] + synthesis = "none" + +[atmos_chem] + module = "none" diff --git a/input/planets/earth.toml b/input/planets/earth.toml index 25e2d474f..79661172e 100644 --- a/input/planets/earth.toml +++ b/input/planets/earth.toml @@ -12,8 +12,8 @@ version = "2.0" [params.dt] starspec = 1e8 starinst = 1e2 - minimum = 1e5 - initial = 1e5 + minimum = 1e6 + initial = 2e6 [params.dt.adaptive] atol = 0.03 @@ -46,8 +46,7 @@ version = "2.0" [star.mors] rot_pcntle = 50.0 age_now = 4.57 # Gyr - star_name = "sun" - spectrum_source = "solar" + spec = "stellar_spectra/Named/sun.txt" [orbit] semimajoraxis = 1.0 @@ -61,32 +60,20 @@ version = "2.0" corefrac = 0.55 [atmos_clim] - surf_state = "skin" + surf_state = "fixed" rayleigh = true module = "agni" - albedo_pl = 0.0 # test guess, contribution from clouds + albedo_pl = 0.2 # test guess, contribution from clouds [atmos_clim.agni] p_top = 1e-4 - num_levels = 40 - dx_max = 100.0 - sens_heat = false - cloud_enabled = false - rainout = false - oceans = true + condensation = true spectral_group = "Dayspring" spectral_bands = "48" - solution_atol = 2.0 - ls_default = 2 [escape] module = "zephyrus" - [escape.zephyrus] - Pxuv = 1e-4 # Pressure at which XUV radiation become opaque in the planetary atmosphere [bar] - efficiency = 0.15 # Escape efficiency factor - tidal = true # Tidal contribution enabled - [interior] radiogenic_heat = false tidal_heat = false @@ -94,13 +81,10 @@ version = "2.0" module = "dummy" [interior.dummy] - ini_tmagma = 3000.0 - H_radio = 8e-12 - mantle_cp = 2000.0 + ini_tmagma = 2000.0 - [interior.spider] - ini_entropy = 2900.0 - num_levels = 60 + [interior.aragog] + ini_tmagma = 3000.0 [outgas] fO2_shift_IW = 4 @@ -112,7 +96,7 @@ version = "2.0" module = "none" [delivery.elements] - H_oceans = 2.0 # Hydrogen inventory in units of equivalent Earth oceans + H_oceans = 1.0 # Hydrogen inventory in units of equivalent Earth oceans # H_ppmw = 109.0 # Hydrogen inventory in ppmw relative to mantle mass CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system diff --git a/input/planets/fiducial_sub_Neptune.toml b/input/planets/fiducial_sub_Neptune.toml index a56921284..f88288257 100644 --- a/input/planets/fiducial_sub_Neptune.toml +++ b/input/planets/fiducial_sub_Neptune.toml @@ -94,8 +94,7 @@ version = "2.0" rot_pctle = 50.0 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 2.4 # Gyr, current age of star used for scaling - star_name = "gj849" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/gj849.txt" # stellar spectrum [star.dummy] calculate_radius = true # calculate radius using empirical mass-luminosity and mass-radius relations @@ -104,7 +103,7 @@ version = "2.0" # Orbital system [orbit] instellation_method = 'inst' # whether to define orbit using semi major axis ('sma') or instellation flux ('inst') - instellationflux = 100 # instellation flux received from the planet in Earth units. + instellationflux = 1 # instellation flux received from the planet in Earth units. semimajoraxis = 0.15 # AU eccentricity = 0.0 # dimensionless zenith_angle = 54.74 # degrees @@ -154,7 +153,7 @@ version = "2.0" solution_atol = 1e-2 # solver absolute tolerance solution_rtol = 1e-1 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = true # volatile condensation + condensation = true # volatile condensation real_gas = true # use real-gas equations of state [atmos_clim.janus] @@ -204,7 +203,7 @@ version = "2.0" ini_dsdr = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1] [interior.aragog] - num_levels = 100 # Number of Aragog grid levels + num_levels = 110 # Number of Aragog grid levels tolerance = 1.0e-7 # solver tolerance ini_tmagma = 3000.0 # Initial magma surface temperature [K] @@ -252,7 +251,7 @@ version = "2.0" metallicity = 1000 # metallicity relative to solar metallicity # H_oceans = 60.0 # Hydrogen inventory in units of equivalent Earth oceans - H_ppmw = 709.0 # Hydrogen inventory in ppmw relative to mantle mass + H_ppmw = 3546.0 # Hydrogen inventory in ppmw relative to mantle mass CH_ratio = 0.32 # C/H mass ratio in mantle/atmosphere system # C_ppmw = 10.0 # Carbon inventory in ppmw relative to mantle mass diff --git a/input/planets/gj9827d.toml b/input/planets/gj9827d.toml index 139b265c7..1f372a218 100644 --- a/input/planets/gj9827d.toml +++ b/input/planets/gj9827d.toml @@ -64,8 +64,7 @@ version = "2.0" rot_pcntle = 50.0 # rotation rate tracks = "spada" # evolution tracks: spada | baraffe age_now = 5.465 # Gyr, current age of star used for scaling - star_name = "hd85512" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/hd85512.txt" # stellar spectrum # Orbital system [orbit] @@ -102,7 +101,7 @@ version = "2.0" solution_atol = 1e-1 # solver absolute tolerance solution_rtol = 8e-2 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = false # use real-gas equations of state [atmos_clim.dummy] diff --git a/input/planets/hd63433d.toml b/input/planets/hd63433d.toml index 576e04635..28c7109d4 100644 --- a/input/planets/hd63433d.toml +++ b/input/planets/hd63433d.toml @@ -23,10 +23,7 @@ version = "2.0" rot_pcntle = 50.0 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 0.414 # Gyr, current age of star used for scaling - star_name = "sun" - spectrum_source = "solar" - - + spec = "stellar_spectra/Named/sun.txt" # stellar spectrum # Orbital system [orbit] @@ -63,7 +60,7 @@ version = "2.0" solution_atol = 1e-2 # solver absolute tolerance solution_rtol = 5e-2 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = true # volatile condensation + condensation = true # volatile condensation real_gas = true # use real-gas equations of state # Volatile escape - physics table diff --git a/input/planets/l9859b.toml b/input/planets/l9859b.toml index c445a7d96..9377a0a8b 100644 --- a/input/planets/l9859b.toml +++ b/input/planets/l9859b.toml @@ -52,9 +52,7 @@ version = "2.0" rot_period = 80.9 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 4.94 # Gyr, current age of star used for scaling - star_name = "l-98-59" - spectrum_source = "muscles" - + spec = "stellar_spectra/Named/l-98-59.txt" # stellar spectrum [orbit] @@ -92,7 +90,7 @@ version = "2.0" solution_atol = 0.5 # solver absolute tolerance solution_rtol = 0.12 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = true # use real-gas equations of state [escape] diff --git a/input/planets/l9859c.toml b/input/planets/l9859c.toml index 0fb6bcf2a..c55b82cd4 100644 --- a/input/planets/l9859c.toml +++ b/input/planets/l9859c.toml @@ -52,8 +52,8 @@ version = "2.0" rot_period = 80.9 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 4.94 # Gyr, current age of star used for scaling - star_name = "l-98-59" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/l-98-59.txt" # stellar spectrum + [orbit] semimajoraxis = 0.0304 # AU @@ -89,7 +89,7 @@ version = "2.0" solution_atol = 0.5 # solver absolute tolerance solution_rtol = 0.12 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = true # use real-gas equations of state [escape] diff --git a/input/planets/l9859d.toml b/input/planets/l9859d.toml index d07d652c5..9de5ad7f8 100644 --- a/input/planets/l9859d.toml +++ b/input/planets/l9859d.toml @@ -57,8 +57,7 @@ version = "2.0" rot_period = 77.5 # rotation period tracks = "spada" # evolution tracks: spada | baraffe age_now = 2.5 # Gyr, current age of star used for scaling - star_name = "l-98-59" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/l-98-59.txt" # stellar spectrum [orbit] @@ -95,7 +94,7 @@ version = "2.0" solution_atol = 2.0 # solver absolute tolerance solution_rtol = 0.14 # solver relative tolerance overlap_method = "ee" # gas overlap method - rainout = false # volatile condensation + condensation = false # volatile condensation real_gas = false # use real-gas equations of state [escape] diff --git a/input/planets/toi6255.toml b/input/planets/toi6255.toml index 966150e3f..15298bf4e 100644 --- a/input/planets/toi6255.toml +++ b/input/planets/toi6255.toml @@ -25,8 +25,7 @@ version = "2.0" rot_period = "none" # rotation period [days] tracks = "spada" # evolution tracks: spada | baraffe age_now = 6 # Gyr, current age of star used for scaling - star_name = "gj176" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/gj176.txt" # stellar spectrum # Orbital system [orbit] diff --git a/input/planets/trappist1c.toml b/input/planets/trappist1c.toml index a67fdb7ed..f096d6b38 100644 --- a/input/planets/trappist1c.toml +++ b/input/planets/trappist1c.toml @@ -7,8 +7,7 @@ version = "2.0" [params] # output files [params.out] - logging = "DEBUG" - path = "trappist1c" + path = "hd63433d" plot_mod = 5 # Plotting frequency, 0: wait until completion | n: every n iterations # ---------------------------------------------------- @@ -23,8 +22,7 @@ version = "2.0" rot_pcntle = 50.0 # rotation percentile tracks = "spada" # evolution tracks: spada | baraffe age_now = 7.6 # Gyr, current age of star used for scaling - star_name = "trappist-1" - spectrum_source = "muscles" + spec = "stellar_spectra/Named/trappist-1.txt" # stellar spectrum # Orbital system [orbit] @@ -52,17 +50,17 @@ version = "2.0" [atmos_clim.agni] p_top = 1.0e-5 # bar, top of atmosphere grid pressure - spectral_group = "Dayspring" # which gas opacities to include - spectral_bands = "48" # how many spectral bands? - num_levels = 40 # Number of atmospheric grid levels - chemistry = "eq" # "none" | "eq" + spectral_group = "Honeyside" # which gas opacities to include + spectral_bands = "256" # how many spectral bands? + num_levels = 60 # Number of atmospheric grid levels + chemistry = "none" # "none" | "eq" surf_material = "greybody" # surface material file for scattering solve_energy = true # solve for energy-conserving atmosphere profile solution_atol = 1e-2 # solver absolute tolerance solution_rtol = 5e-2 # solver relative tolerance - overlap_method = "ee" # gas overlap method - rainout = true # volatile condensation - real_gas = false # use real-gas equations of state + overlap_method = "ee" # gas overlap method + condensation = true # volatile condensation + real_gas = true # use real-gas equations of state # Volatile escape - physics table [escape] @@ -76,7 +74,7 @@ version = "2.0" module = "spider" # Which interior module to use [interior.spider] - num_levels = 200 # Number of SPIDER grid levels + num_levels = 220 # Number of SPIDER grid levels tolerance = 1.0e-10 # solver tolerance tsurf_atol = 20.0 # tsurf_poststep_change tsurf_rtol = 0.01 # tsurf_poststep_change_frac diff --git a/src/proteus/atmos_clim/dummy.py b/src/proteus/atmos_clim/dummy.py index c6f00749f..d4a475771 100644 --- a/src/proteus/atmos_clim/dummy.py +++ b/src/proteus/atmos_clim/dummy.py @@ -56,8 +56,11 @@ def _calc_fluxes(x): # fixed T_Surf if config.atmos_clim.surf_state == 'fixed': - log.info('Calculating fluxes with dummy atmosphere') - T_surf_atm = T_magma + log.info("Calculating fluxes with dummy atmosphere") + if config.interior.module == 'boundary': + T_surf_atm = hf_row["T_surf"] + else: + T_surf_atm = T_magma fluxes = _calc_fluxes(T_surf_atm) # conductive lid diff --git a/src/proteus/atmos_clim/wrapper.py b/src/proteus/atmos_clim/wrapper.py index bdff334a7..2567d523a 100644 --- a/src/proteus/atmos_clim/wrapper.py +++ b/src/proteus/atmos_clim/wrapper.py @@ -87,8 +87,12 @@ def run_atmosphere( hf_row['albedo_pl'] = 0.0 # Handle new surface temperature - if config.atmos_clim.surf_state == 'mixed_layer': - hf_row['T_surf'] = ShallowMixedOceanLayer(hf_all.iloc[-1].to_dict(), hf_row) + if config.interior.module=='boundary': + #surface temperature is already calculated by boundary interior module + pass + + elif config.atmos_clim.surf_state == 'mixed_layer': + hf_row["T_surf"] = ShallowMixedOceanLayer(hf_all.iloc[-1].to_dict(), hf_row) elif config.atmos_clim.surf_state == 'fixed': hf_row['T_surf'] = hf_row['T_magma'] @@ -147,7 +151,9 @@ def run_atmosphere( if no_atm: activate_julia(dirs, config.atmos_clim.agni.verbosity) # surface temperature guess - hf_row['T_surf'] = hf_row['T_magma'] + + if config.interior.module != 'boundary': + hf_row["T_surf"] = hf_row["T_magma"] else: # Remove old spectral file if it exists safe_rm(spfile_path) diff --git a/src/proteus/config/_interior.py b/src/proteus/config/_interior.py index 61e02548d..925f5b4f5 100644 --- a/src/proteus/config/_interior.py +++ b/src/proteus/config/_interior.py @@ -23,6 +23,20 @@ def valid_spider(instance, attribute, value): raise ValueError('Must enable at least one energy transport term in SPIDER') +def valid_interiorboundary(instance, attribute, value): + if instance.module != "boundary": + return + + tsol = instance.boundary.T_solidus + tliq = instance.boundary.T_liquidus + if tliq <= tsol: + raise ValueError(f"Boundary liquidus ({tliq}K) must be greater than solidus ({tsol}K)") + + t_surf_0 = instance.boundary.T_surf_0 + t_p_0 = instance.boundary.T_p_0 + if t_surf_0 > t_p_0: + raise ValueError(f"Initial surface temperature ({t_surf_0}K) must be less than or equal to initial potential temperature ({t_p_0}K)") + def valid_path(instance, attribute, value): if not isinstance(value, str) or not value.strip(): raise ValueError(f"'{attribute.name}' must be a non-empty string") @@ -186,6 +200,72 @@ def valid_interiordummy(instance, attribute, value): raise ValueError(f'Dummy liquidus ({tliq}K) must be greater than solidus ({tsol}K)') +@define +class InteriorBoundary: + """Parameters for Boundary interior module. Default values taken from Schaefer et al. 2016 (https://iopscience.iop.org/article/10.3847/0004-637X/829/2/63/pdf). + + Attributes + ---------- + rtol: float + ODE solver relative tolerance. + atol: float + ODE solver absolute tolerance. + T_surf_0: float + Initial surface temperature [K] for boundary solver. + T_p_0: float + Initial mantle potential temperature [K] for boundary solver. + T_solidus: float + Mantle solidus temperature [K]. + T_liquidus: float + Mantle liquidus temperature [K]. + critical_melt_fraction: float + Critical melt fraction for rheological transition. + critical_rayleigh_number: float + Critical Rayleigh number for onset of convection [-]. + heat_fusion_silicate: float + Latent heat of fusion for silicates [J/kg]. + nusselt_exponent: float + Nusselt-Rayleigh scaling exponent [-]. + silicate_heat_capacity: float + Silicate heat capacity [J/kg/K]. + thermal_conductivity: float + Thermal conductivity [W/m/K]. + thermal_diffusivity: float + Thermal diffusivity [m^2/s]. + thermal_expansivity: float + Thermal expansivity [1/K]. + transition_width: float + Width of viscosity transition in melt fraction space [-]. + eta_solid_const: float + Constant solid viscosity for aggregate formulation [Pa s]. + eta_melt_const: float + Constant melt viscosity for aggregate formulation [Pa s]. + """ + + rtol: float = field(default=1e-6, validator=gt(0)) + atol: float = field(default=1e-9, validator=gt(0)) + + T_surf_0: float = field(default=3999.0, validator=ge(0)) + T_p_0: float = field(default=4000.0, validator=ge(0)) + + T_solidus: float = field(default=1420.0, validator=ge(0)) + T_liquidus: float = field(default=2020.0, validator=gt(0)) + + critical_melt_fraction: float = field(default=0.4, validator=(gt(0), lt(1))) + + critical_rayleigh_number: float = field(default=1.1e3, validator=gt(0)) # - + heat_fusion_silicate: float = field(default=4.0e5, validator=gt(0)) # J/kg + nusselt_exponent: float = field(default=0.33, validator=gt(0)) # - + silicate_heat_capacity: float = field(default=1.2e3, validator=gt(0)) # J/kg/K + thermal_conductivity: float = field(default=4.2, validator=gt(0)) # W/m/K + thermal_diffusivity: float = field(default=1e-6, validator=gt(0)) # m^2/s + thermal_expansivity: float = field(default=2e-5, validator=gt(0)) # 1/K + + # Aggregate viscosity parameters + transition_width: float = field(default=0.15, validator=(gt(0), lt(1))) # - + eta_solid_const: float = field(default=1e21, validator=gt(0)) # Pa s + eta_melt_const: float = field(default=1e2, validator=gt(0)) # Pa s + @define class InteriorDummy: """Parameters for Dummy interior module. @@ -260,15 +340,16 @@ class Interior: Zalmoxis derives its EOS paths from struct.zalmoxis config instead. """ - module: str = field(validator=in_(('spider', 'aragog', 'dummy'))) - melting_dir: str = field(default='Monteux-600', validator=valid_path) + module: str = field(validator=in_(('spider', 'aragog', 'dummy', 'boundary'))) + melting_dir: str = field(default="Monteux-600",validator=valid_path) eos_dir: str = field(default='WolfBower2018_MgSiO3', validator=valid_path) - radiogenic_heat: bool = field(default=True) - tidal_heat: bool = field(default=True) + radiogenic_heat: bool = field(default=True) + tidal_heat: bool = field(default=True) - spider: Spider = field(factory=Spider, validator=valid_spider) - aragog: Aragog = field(factory=Aragog, validator=valid_aragog) - dummy: InteriorDummy = field(factory=InteriorDummy, validator=valid_interiordummy) + spider: Spider = field(factory=Spider, validator=valid_spider) + aragog: Aragog = field(factory=Aragog, validator=valid_aragog) + dummy: InteriorDummy = field(factory=InteriorDummy, validator=valid_interiordummy) + boundary: InteriorBoundary = field(factory=InteriorBoundary, validator=valid_interiorboundary) grain_size: float = field(default=0.1, validator=gt(0)) F_initial: float = field(default=1e3, validator=gt(0)) diff --git a/src/proteus/interior/boundary.py b/src/proteus/interior/boundary.py new file mode 100644 index 000000000..e5d5f5614 --- /dev/null +++ b/src/proteus/interior/boundary.py @@ -0,0 +1,546 @@ +from __future__ import annotations + +import csv +import logging +from typing import TYPE_CHECKING + +import numpy as np +import pandas as pd +from scipy.integrate import solve_ivp + +from proteus.atmos_clim.common import Atmos_t +from proteus.interior.common import Interior_t +from proteus.interior.timestep import next_step +from proteus.utils.constants import ( + M_earth, + const_G, + radnuc_data, + secs_per_year, +) + +if TYPE_CHECKING: + from proteus.config import Config + + +class LineLimitFilter(logging.Filter): + def __init__(self, log_file: str, max_lines: int) -> None: + super().__init__() + self.log_file = log_file + self.max_lines = max_lines + self.line_count = self._count_lines() + + def _count_lines(self) -> int: + try: + with open(self.log_file, 'r', encoding='utf-8') as handle: + return sum(1 for _ in handle) + except FileNotFoundError: + return 0 + + def filter(self, record: logging.LogRecord) -> bool: + if self.line_count >= self.max_lines: + return False + message = record.getMessage() + lines_to_add = message.count("\n") + 1 + if self.line_count + lines_to_add > self.max_lines: + return False + self.line_count += lines_to_add + return True + +class BoundaryRunner(): + + def __init__(self, config: Config, dirs: dict, hf_row: dict, hf_all: + pd.DataFrame, interior_o: Interior_t, atmos_o: Atmos_t): + + self.curr_time = hf_row["Time"] * secs_per_year + self.dt = self.compute_time_step(config, dirs, hf_row, hf_all, interior_o) * secs_per_year + self.iteration = 1 if hf_all is None else len(hf_all) + + self.planet_radius = hf_row["R_int"] + self.planet_mass = config.struct.mass_tot * M_earth + self.core_mass = hf_row["M_core"] + self.m_atm = hf_row["M_atm"] + self.f_atm = hf_row["F_atm"] + + cp_layer = getattr(getattr(atmos_o, "_atm", None), "layer_cp", None) + if cp_layer is not None and len(cp_layer) > 2 and not pd.isna(cp_layer[2]): + self.atmosphere_heat_capacity = cp_layer[2] # J/kg/K + else: + self.atmosphere_heat_capacity = 1.7e4 # J/kg/K for H2 at 2000K + + self.core_radius = config.struct.corefrac * self.planet_radius + + self.mantle_radius = self.planet_radius - self.core_radius + self.mantle_volume = (4/3) * np.pi * (self.planet_radius**3 - self.core_radius**3) + self.mantle_mass = (self.planet_mass - self.core_mass) + self.bulk_density = self.mantle_mass / self.mantle_volume + + self.surface_gravity = const_G * self.planet_mass / self.planet_radius**2 + self.scale_length_mantle = self.mantle_radius / 3 + + self.rtol = config.interior.boundary.rtol + self.atol = config.interior.boundary.atol + + if interior_o.ic == 2: + self.T_p_0 = hf_row.get("T_magma", 0.0) + self.T_surf_0 = hf_row.get("T_surf", 0.0) + else: + self.T_p_0 = config.interior.boundary.T_p_0 + self.T_surf_0 = config.interior.boundary.T_surf_0 + + if self.T_surf_0 > self.T_p_0: + self.T_surf_0 = self.T_p_0 - 1.0 # Ensure initial surface temperature does not exceed potential temperature + + self.T_solidus = config.interior.boundary.T_solidus + self.T_liquidus = config.interior.boundary.T_liquidus + self.critical_melt_fraction = config.interior.boundary.critical_melt_fraction + + # Material constants + self.critical_rayleigh_number = config.interior.boundary.critical_rayleigh_number # dimensionless + self.heat_fusion_silicate = config.interior.boundary.heat_fusion_silicate # J/kg + self.nusselt_exponent = config.interior.boundary.nusselt_exponent # dimensionless + self.silicate_heat_capacity = config.interior.boundary.silicate_heat_capacity # J/kg/K + self.thermal_conductivity = config.interior.boundary.thermal_conductivity # W/m/K + self.thermal_diffusivity = config.interior.boundary.thermal_diffusivity # m^2/s + self.thermal_expansivity = config.interior.boundary.thermal_expansivity # 1/K + + # Aggregate viscosity parameters + self.transition_width = config.interior.boundary.transition_width # dimensionless + self.eta_solid_const = config.interior.boundary.eta_solid_const # Pa s + self.eta_melt_const = config.interior.boundary.eta_melt_const # Pa s + + # Radioactive heating parameters + self.use_radiogenic_heating = config.interior.radiogenic_heat + self.radio_tref = config.delivery.radio_tref + self.U_abun = config.delivery.radio_U * 1e-6 # Convert ppm to kg/kg + self.Th_abun = config.delivery.radio_Th * 1e-6 # Convert ppm to kg/kg + self.K_abun = config.delivery.radio_K * 1e-6 # Convert ppm to kg/kg + + @staticmethod + def compute_time_step(config: Config, dirs: dict, hf_row: dict, + hf_all: pd.DataFrame, interior_o: Interior_t) -> float: + if interior_o.ic == 1: + return 0.0 + else: + step_sf = 1.0 # dt scale factor + return next_step(config, dirs, hf_row, hf_all, step_sf) + + def viscosity_aggregate(self, phi: float) -> float: + """ + Calculate the aggregate viscosity using a smooth transition function. + + This function blends between solid and magma ocean viscosities using a + hyperbolic tangent transition function centered at the critical melt fraction. + + Parameters + ---------- + phi : float + Melt fraction [0-1] + + Returns + ------- + float + Aggregate dynamic viscosity [Pa s] + """ + # Use constant viscosities from config + eta_solid = self.eta_solid_const + eta_magma = self.eta_melt_const + + # Calculate transition parameter + y = (phi - self.critical_melt_fraction) / self.transition_width + + # Calculate transition function (0 to 1) + z = 0.5 * (1 + np.tanh(y)) + + # Calculate aggregate viscosity using logarithmic interpolation + log_eta = z * np.log10(eta_magma) + (1 - z) * np.log10(eta_solid) + eta = 10**log_eta + + return eta + + def rayleigh_number(self, T_p: float, T_surf: float, phi: float) -> float: + """ + Calculate the Rayleigh number for mantle convection. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + T_surf : float + Surface temperature of the planet [K] + phi : float + Melt fraction [0-1] + + Returns + ------- + float + Rayleigh number [dimensionless] + """ + # Determine viscosity based on configuration + eta = self.viscosity_aggregate(phi) + + # Calculate Rayleigh number + Ra = ((self.bulk_density * self.surface_gravity * self.thermal_expansivity * + np.abs(T_p - T_surf) * self.scale_length_mantle**3) / + (eta * self.thermal_diffusivity)) + + return Ra + + def q_m(self, T_p: float, T_surf: float, phi: float) -> float: + """ + Calculate the convective heat flux from the mantle. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + T_surf : float + Surface temperature of the planet [K] + phi : float + Melt fraction [0-1] + + Returns + ------- + float + Convective heat flux [W/m^2] + """ + Ra = self.rayleigh_number(T_p, T_surf, phi) + + # Nusselt number scaling relation + Nu = (Ra / self.critical_rayleigh_number)**self.nusselt_exponent + + # Calculate convective heat flux + q_m_val = Nu * self.thermal_conductivity * np.abs(T_p - T_surf) / self.scale_length_mantle + + return q_m_val + + def radioactive_heating(self, t: float) -> float: + """ + Calculate the volumetric radioactive heating rate as a function of time. + + Parameters + ---------- + t : float + Time [s] + + Returns + ------- + float + Radioactive heating rate per unit volume [W/m^3] + """ + # Get radioactive constants from radnuc_data + secs_per_year = 3.154e7 # s + + # Radioactive isotope properties from radnuc_data + # Convert half-lives to decay constants: lambda = ln(2) / half-life + K_decay_constant = np.log(2) / (radnuc_data['k40']['halflife'] * secs_per_year) # 1/s + K_heat_production = radnuc_data['k40']['heatprod'] # W/kg + + Th_decay_constant = np.log(2) / (radnuc_data['th232']['halflife'] * secs_per_year) # 1/s + Th_heat_production = radnuc_data['th232']['heatprod'] # W/kg + + Ur_decay_constant = np.log(2) / (radnuc_data['u238']['halflife'] * secs_per_year) # 1/s + Ur_heat_production = radnuc_data['u238']['heatprod'] # W/kg + + # Calculate heating contributions from each isotope + K_term = self.K_abun * K_heat_production * np.exp( + K_decay_constant * (self.radio_tref * secs_per_year - t) + ) + Ur_term = self.U_abun * Ur_heat_production * np.exp( + Ur_decay_constant * (self.radio_tref * secs_per_year - t) + ) + Th_term = self.Th_abun * Th_heat_production * np.exp( + Th_decay_constant * (self.radio_tref * secs_per_year - t) + ) + H_total = K_term + Ur_term + Th_term + + if self.use_radiogenic_heating: + return H_total + else: + return 0.0 + + def melt_fraction(self, T_p: float) -> float: + """ + Calculate the melt fraction from potential temperature. + + Parameters + ---------- + T_p : float + Potential temperature [K] + + Returns + ------- + float + Melt fraction [0-1] + """ + phi = (T_p - self.T_solidus) / (self.T_liquidus - self.T_solidus) + return np.clip(phi, 0.0, 1.0) + + def r_s(self, T_p) -> float: + """ + Calculate the solidification radius directly from the potential temperature. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + + Returns + ------- + float + Solidification radius [m] + """ + phi = self.melt_fraction(T_p) + + if phi >= 1.0: + return self.core_radius # Fully molten, solidification radius at core-mantle boundary + elif phi <= 0.0: + return self.planet_radius # Fully solid, solidification radius at surface + else: + # Linear interpolation between core and surface based on melt fraction + return (self.planet_radius**3 - phi * (self.planet_radius**3 - self.core_radius**3))**(1/3) + + def drs_dTp(self, T_p: float) -> float: + """ + Calculate the derivative of solidification radius with respect to potential temperature. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + + Returns + ------- + float + Derivative dr_s/dT_p [m/K] + """ + + # Compute r_s + r_s = self.r_s(T_p) + + # Compute dr_s/dT_p using the chain rule as derived in drs_dt + volume_diff = self.planet_radius**3 - self.core_radius**3 + T_range = self.T_liquidus - self.T_solidus + + dr_s_dT_p = -(volume_diff) / (3 * T_range * r_s**2) + + return dr_s_dT_p + + def dT_pdt(self, T_p: float, T_surf: float, t: float) -> float: + """ + Calculate the rate of change of potential temperature of the mantle. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + T_surf : float + Surface temperature of the planet [K] + t : float + Time [s] + + Returns + ------- + float + Rate of change of potential temperature [K/s] + """ + # Calculate melt fraction + phi = self.melt_fraction(T_p) + + # Calculate convective heat flux + q_m_val = self.q_m(T_p, T_surf, phi) + + # Calculate radioactive heating + Q_val = self.radioactive_heating(t) + + # Energy balance numerator: heat loss - radiogenic heating + numerator = (-4 * np.pi * self.planet_radius**2 * q_m_val + + (4/3) * np.pi * self.bulk_density * Q_val * + (self.planet_radius**3 - self.core_radius**3)) + + r_s_val = self.r_s(T_p) # Update r_s_val based on current T_p + + # Calculate dr_s/dT_p for the latent heat term + dr_s_dT_p = self.drs_dTp(T_p) + + # Energy balance denominator: sensible heat + latent heat + if r_s_val < self.planet_radius: + denominator = ((4/3) * np.pi * self.bulk_density * self.silicate_heat_capacity * + (self.planet_radius**3 - r_s_val**3) - + 4 * np.pi * r_s_val**2 * self.bulk_density * self.heat_fusion_silicate * dr_s_dT_p) + else: + denominator = self.silicate_heat_capacity * self.mantle_mass + + dT_pdt_val = numerator / denominator + + return dT_pdt_val + + def dT_surfdt(self, T_p: float, T_surf: float) -> float: + """ + Calculate the rate of change of surface temperature of the planet. + + Parameters + ---------- + T_p : float + Potential temperature of the mantle [K] + T_surf : float + Surface temperature of the planet [K] + + Returns + ------- + dT_surfdt_val : float + Rate of change of surface temperature [K/s] + """ + phi = self.melt_fraction(T_p) + q_m_val = self.q_m(T_p, T_surf, phi) + + delta = self.thermal_conductivity * (T_p - T_surf) / q_m_val + + numerator = 4 * np.pi * self.planet_radius**2 * (q_m_val - self.f_atm) + denominator = self.atmosphere_heat_capacity * self.m_atm + \ + (4/3) * np.pi * self.silicate_heat_capacity * self.bulk_density * (self.planet_radius**3 - (self.planet_radius-delta)**3) + + dT_surfdt_val = numerator / denominator + + return dT_surfdt_val + + def thermal_rhs(self, t: float, y: list) -> list: + """ + Right-hand side function for the coupled thermal evolution ODEs. + like scipy.integrate.solve_ivp. + + Parameters + ---------- + t : float + Time [s] + y : list or array-like + State vector containing [T_p, T_surf, r_s] where: + - T_p is the mantle potential temperature [K] + - T_surf is the surface temperature [K] + + Returns + ------- + list + Time derivatives [dT_p/dt, dT_surf/dt, drs/dt] where: + - dT_p/dt is the rate of change of potential temperature [K/s] + - dT_surf/dt is the rate of change of surface temperature [K/s] + """ + T_p, T_surf = y + + if T_surf>T_p: + T_surf = T_p-1.0 # Ensure surface temperature does not exceed potential temperature + + dTp = self.dT_pdt(T_p, T_surf, t) + dTs = self.dT_surfdt(T_p, T_surf) + + return [dTp, dTs] + + def run_solver(self, hf_row: dict, interior_o: Interior_t, dirs: dict) -> tuple: + """ + Run the thermal evolution solver for a single timestep. + + Parameters + ---------- + interior_o : Interior_t + Interior model object containing structural and thermodynamic state + dirs : dict + Dictionary of directory paths for output files + + Returns + ------- + tuple + A tuple containing: + - sim_time : float + Final simulation time after integration [s] + - output : dict + Dictionary containing thermal evolution results + """ + # Set up CSV logging for step diagnostics. + output_dir = dirs.get('output', '.') + csv_log_file = f"{output_dir}/boundary_solver_debug.csv" + self._run_solver_call_count = self.iteration + + csv_columns = [ + "call_index", + "step_index", + "n_steps", + "time_s", + "dt_s", + "T_p_K", + "T_surf_K", + "dT_p_dt_K_per_s", + "dT_surf_dt_K_per_s", + "convective_heat_loss_term_W", + "radiogenic_heating_term_W", + "F_atmosphere_W_per_m2", + ] + csv_needs_header = not pd.io.common.file_exists(csv_log_file) + + y0 = [self.T_p_0, self.T_surf_0] + t_span = (self.curr_time, self.curr_time + self.dt) + sol = solve_ivp(self.thermal_rhs, t_span, y0, method='BDF', rtol=self.rtol, atol=self.atol, dense_output=True) + + with open(csv_log_file, mode='a', encoding='utf-8', newline='') as handle: + writer = csv.DictWriter(handle, fieldnames=csv_columns) + if csv_needs_header: + writer.writeheader() + + for i, t in enumerate(sol.t): + timestep = t - sol.t[i-1] if i > 0 else t - self.curr_time + T_p = sol.y[0, i] + T_surf = sol.y[1, i] + dT_pdt_val = self.dT_pdt(T_p, T_surf, t) + dT_surfdt_val = self.dT_surfdt(T_p, T_surf) + q_m_val = self.q_m(T_p, T_surf, self.melt_fraction(T_p)) + q_m_term = 4 * np.pi * self.planet_radius**2 * q_m_val + q_rad = self.radioactive_heating(t) + q_rad_term = (4 / 3) * np.pi * self.bulk_density * q_rad * (self.planet_radius**3 - self.core_radius**3) + + writer.writerow( + { + "call_index": self._run_solver_call_count, + "step_index": i + 1, + "n_steps": len(sol.t), + "time_s": t, + "dt_s": timestep, + "T_p_K": T_p, + "T_surf_K": T_surf, + "dT_p_dt_K_per_s": dT_pdt_val, + "dT_surf_dt_K_per_s": dT_surfdt_val, + "convective_heat_loss_term_W": q_m_term, + "radiogenic_heating_term_W": q_rad_term, + "F_atmosphere_W_per_m2": self.f_atm, + } + ) + + # Extract results + T_p_final = sol.y[0, -1] + T_surf_final = sol.y[1, -1] + r_s_final = self.r_s(T_p_final) + sim_time = (self.curr_time + self.dt)/secs_per_year # convert back to years + q_m_final = self.q_m(T_p_final, T_surf_final, self.melt_fraction(T_p_final)) + phi_final = self.melt_fraction(T_p_final) + f_radio_final = self.radioactive_heating(self.curr_time + self.dt) + + m_liquid = (4/3) * np.pi * self.bulk_density * (self.planet_radius**3 - r_s_final**3) + m_solid = self.mantle_mass - m_liquid + + if T_surf_final > T_p_final: + T_surf_final = T_p_final - 1.0 # Ensure surface temperature does not exceed potential temperature + + output = { + "T_magma": T_p_final, + "T_pot": T_p_final, + "T_surf": T_surf_final, + "F_int": q_m_final, + "Phi_global": phi_final, + "Phi_global_vol": phi_final, + "F_radio": f_radio_final/(4*np.pi*self.planet_radius**2), + "RF_depth": r_s_final/self.planet_radius, + "M_mantle_liquid": m_liquid, + "M_mantle_solid": m_solid, + "F_tidal": 0.0, + "M_mantle": self.mantle_mass, + "M_core": self.core_mass + } + + return sim_time, output diff --git a/src/proteus/interior/wrapper.py b/src/proteus/interior/wrapper.py index 932ff8767..3a4c8b3dd 100644 --- a/src/proteus/interior/wrapper.py +++ b/src/proteus/interior/wrapper.py @@ -11,6 +11,7 @@ import pandas as pd import scipy.optimize as optimise +from proteus.atmos_clim.common import Atmos_t from proteus.interior.common import Interior_t from proteus.outgas.wrapper import calc_target_elemental_inventories from proteus.utils.constants import M_earth, R_earth, const_G, element_list @@ -67,7 +68,9 @@ def get_nlevb(config: Config): return int(config.interior.spider.num_levels) case 'aragog': return int(config.interior.aragog.num_levels) - case 'dummy': + case "boundary": + return 2 + case "dummy": return 2 raise ValueError(f"Invalid interior module selected '{config.interior.module}'") @@ -132,20 +135,22 @@ def _resid(x): case _: rtol = 1e-7 - # Find the radius - r = optimise.root_scalar( - _resid, - method='secant', - xtol=1e3, - rtol=rtol, - maxiter=10, - x0=hf_row['R_int'], - x1=hf_row['R_int'] * 1.02, - ) - hf_row['R_int'] = float(r.root) - calculate_core_mass(hf_row, config) - run_interior(dirs, config, hf_all, hf_row, int_o) - update_gravity(hf_row) + if config.interior.module != "boundary": + + # Find the radius + r = optimise.root_scalar(_resid, method='secant', xtol=1e3, rtol=rtol, maxiter=10, + x0=hf_row["R_int"], x1=hf_row["R_int"]*1.02) + hf_row["R_int"] = float(r.root) + calculate_core_mass(hf_row, config) + run_interior(dirs, config, hf_all, hf_row, int_o) + update_gravity(hf_row) + + else: + # *** will change this once interior module works with zalmoxis *** + hf_row["R_int"] = R_earth + hf_row["M_tot"] = M_target + calculate_core_mass(hf_row, config) + hf_row["gravity"] = 9.81 # Result log.info('Found solution for interior structure') @@ -220,7 +225,8 @@ def determine_interior_radius_with_zalmoxis( # by the finally block above), not the overridden 'adiabatic'. This is # correct: the Zalmoxis solver already used the adiabatic mode to compute # the structure, and run_interior (SPIDER/ARAGOG) manages its own T(r). - run_interior(dirs, config, hf_all, hf_row, int_o) + if config.interior.module!="boundary": + run_interior(dirs, config, hf_all, hf_row, int_o) def solve_structure( @@ -238,7 +244,9 @@ def solve_structure( # We might need here to setup a determine_interior_mass function as mass calculation depends on gravity if config.struct.set_by == 'radius_int': # radius defines interior structure - hf_row['R_int'] = config.struct.radius_int * R_earth + if config.interior.module=="boundary": + raise ValueError("Must set structure by 'mass_tot' if boundary interior module is used") + hf_row["R_int"] = config.struct.radius_int * R_earth calculate_core_mass(hf_row, config) # initial guess for mass, which will be updated by the interior model hf_row['M_int'] = 1.2 * M_earth @@ -270,13 +278,13 @@ def solve_structure( else: log.error('Invalid constraint on interior structure: %s' % config.struct.set_by) - def run_interior( dirs: dict, config: Config, hf_all: pd.DataFrame, hf_row: dict, interior_o: Interior_t, + atmos_o: Atmos_t=None, verbose: bool = True, ): """Run interior mantle evolution model. @@ -293,6 +301,8 @@ def run_interior( Dictionary of current runtime variables interior_o : Interior_t Interior struct. + atmos_o : Atmos_t + Atmosphere struct (only required for boundary module). verbose : bool Verbose printing enabled. """ @@ -322,6 +332,14 @@ def run_interior( # Run Aragog sim_time, output = AragogRunnerInstance.run_solver(hf_row, interior_o, dirs) + elif config.interior.module == 'boundary': + from proteus.interior.boundary import BoundaryRunner + BoundaryRunnerInstance = BoundaryRunner(config, dirs, hf_row, hf_all, + interior_o, atmos_o) + # Run the boundary interior module + sim_time, output = BoundaryRunnerInstance.run_solver(hf_row, interior_o, + dirs) + elif config.interior.module == 'dummy': # Import from proteus.interior.dummy import run_dummy_int diff --git a/src/proteus/proteus.py b/src/proteus/proteus.py index 8a75e930b..5fd40e54e 100644 --- a/src/proteus/proteus.py +++ b/src/proteus/proteus.py @@ -13,7 +13,11 @@ import proteus.utils.archive as archive from proteus.config import read_config_object -from proteus.utils.constants import vap_list, vol_list +from proteus.utils.constants import ( + M_earth, + vap_list, + vol_list, +) from proteus.utils.helper import ( CleanDir, PrintHalfSeparator, @@ -405,6 +409,15 @@ def start(self, *, resume: bool = False, offline: bool = False): ) ############### INTERIOR + if self.config.interior.module == "boundary" and self.loops["total"] == 0: + + self.hf_row["M_mantle"] = self.config.struct.mass_tot * M_earth - self.hf_row["M_core"] + self.hf_row["Phi_global"] = 1.0 + self.hf_row["T_magma"] = self.config.interior.boundary.T_p_0 + + calc_target_elemental_inventories(self.directories, self.config, self.hf_row) + run_outgassing(self.directories, self.config, self.hf_row) + PrintHalfSeparator() # Evolve interior