From 73051852fbf14cb909fda11c8566730d7b729da1 Mon Sep 17 00:00:00 2001 From: rdc49 Date: Mon, 26 Jan 2026 10:22:19 +0000 Subject: [PATCH 1/4] commit before pulling from main branch --- input/chili/intercomp/_base.grid.toml | 6 +- input/chili/intercomp/_base.toml | 9 +- input/chili/protocol/earth.toml | 82 ++++- input/chili/protocol/tr1b.toml | 82 ++++- input/chili/readme.md | 1 - input/demos/escape.toml | 3 +- input/ensembles/SNS_0.5%_sanity_check.toml | 34 +++ input/ensembles/SNS_0.5%_sweep.toml | 40 +++ input/ensembles/solubility_test.toml | 22 ++ input/planets/GDSP_fiducial.toml | 281 ++++++++++++++++++ input/planets/earth.toml | 4 +- input/planets/fiducial_sub_Neptune.toml | 11 +- .../fiducial_sub_Neptune_Tang_comparison.toml | 281 ++++++++++++++++++ input/planets/gj9827d.toml | 3 +- input/planets/hd63433d.toml | 5 +- input/planets/l9859b.toml | 4 +- input/planets/l9859c.toml | 4 +- input/planets/l9859d.toml | 3 +- input/planets/toi6255.toml | 3 +- input/planets/trappist1c.toml | 3 +- 20 files changed, 840 insertions(+), 41 deletions(-) create mode 100644 input/ensembles/SNS_0.5%_sanity_check.toml create mode 100644 input/ensembles/SNS_0.5%_sweep.toml create mode 100644 input/ensembles/solubility_test.toml create mode 100644 input/planets/GDSP_fiducial.toml create mode 100644 input/planets/fiducial_sub_Neptune_Tang_comparison.toml diff --git a/input/chili/intercomp/_base.grid.toml b/input/chili/intercomp/_base.grid.toml index 2483e9835..a81279d9d 100644 --- a/input/chili/intercomp/_base.grid.toml +++ b/input/chili/intercomp/_base.grid.toml @@ -1,15 +1,15 @@ output = "chili_base_grid/" symlink = "" ref_config = "input/chili/intercomp/_base.toml" -use_slurm = false +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 15adab058..1f7919a4d 100644 --- a/input/chili/intercomp/_base.toml +++ b/input/chili/intercomp/_base.toml @@ -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 @@ -117,8 +116,8 @@ spectral_bands = "48" p_top = 1e-04 p_obs = 0.02 surf_material = "surface_albedos/Hammond24/lunarmarebasalt.dat" -num_levels = 45 -chemistry = "eq" +num_levels = 50 +chemistry = "none" solve_energy = true solution_atol = 0.08 solution_rtol = 0.1 @@ -133,7 +132,7 @@ ls_default = 1 max_steps = 70 perturb_all = true mlt_criterion = "s" -fastchem_floor = 500.0 +fastchem_floor = 273.0 ini_profile = "isothermal" [atmos_clim.dummy] diff --git a/input/chili/protocol/earth.toml b/input/chili/protocol/earth.toml index b8d476b1e..ad8e18b80 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" @@ -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..880c5ef9c 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" @@ -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/escape.toml b/input/demos/escape.toml index 7307b54a7..b4a966af6 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 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/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/planets/GDSP_fiducial.toml b/input/planets/GDSP_fiducial.toml new file mode 100644 index 000000000..d79ab7b39 --- /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 = 1.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 = 4 # 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 = 3546.0 # Hydrogen inventory in ppmw relative to mantle mass + + CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system + # C_ppmw = 10.0 # Carbon inventory in ppmw relative to mantle mass + + NH_ratio = 0.0 # 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 82167268b..6ac238ba7 100644 --- a/input/planets/earth.toml +++ b/input/planets/earth.toml @@ -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 @@ -85,7 +84,6 @@ version = "2.0" [interior.dummy] ini_tmagma = 2000.0 - [interior.aragog] ini_tmagma = 3000.0 diff --git a/input/planets/fiducial_sub_Neptune.toml b/input/planets/fiducial_sub_Neptune.toml index 979fddc67..d935cb3f1 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 @@ -192,7 +191,7 @@ version = "2.0" rheo_phi_loc = 0.4 # Centre of rheological transition rheo_phi_wid = 0.15 # Width of rheological transition - module = "aragog" # Which interior module to use + module = "spider" # Which interior module to use [interior.spider] num_levels = 110 # Number of SPIDER grid levels @@ -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/fiducial_sub_Neptune_Tang_comparison.toml b/input/planets/fiducial_sub_Neptune_Tang_comparison.toml new file mode 100644 index 000000000..ad1248c31 --- /dev/null +++ b/input/planets/fiducial_sub_Neptune_Tang_comparison.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 = "fiducial_sub_Neptune" + 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 = true + phi_crit = 0.01 # non-dim., model will terminate when global melt fraction < phi_crit + + # radiative equilibrium + [params.stop.radeqm] + enabled = true + atol = 0.2 # 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 = 10 # 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 = 4.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 = false # 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 = 60.0 # Hydrogen inventory in units of equivalent Earth oceans + H_ppmw = 14180.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/gj9827d.toml b/input/planets/gj9827d.toml index 139b265c7..4c475ca0d 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] diff --git a/input/planets/hd63433d.toml b/input/planets/hd63433d.toml index 576e04635..90a3a915b 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] diff --git a/input/planets/l9859b.toml b/input/planets/l9859b.toml index c445a7d96..024c0c3e0 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] diff --git a/input/planets/l9859c.toml b/input/planets/l9859c.toml index 0fb6bcf2a..db69e5ecc 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 diff --git a/input/planets/l9859d.toml b/input/planets/l9859d.toml index d07d652c5..41b09cc84 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] 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..856d84581 100644 --- a/input/planets/trappist1c.toml +++ b/input/planets/trappist1c.toml @@ -23,8 +23,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] From dd817c2280f92c198b07523ebba52ca30b5225c0 Mon Sep 17 00:00:00 2001 From: rdc49 Date: Mon, 26 Jan 2026 17:42:29 +0000 Subject: [PATCH 2/4] added git ignore to some .toml files --- input/planets/GDSP_fiducial.toml | 281 ------------------ .../fiducial_sub_Neptune_Tang_comparison.toml | 281 ------------------ 2 files changed, 562 deletions(-) delete mode 100644 input/planets/GDSP_fiducial.toml delete mode 100644 input/planets/fiducial_sub_Neptune_Tang_comparison.toml diff --git a/input/planets/GDSP_fiducial.toml b/input/planets/GDSP_fiducial.toml deleted file mode 100644 index d79ab7b39..000000000 --- a/input/planets/GDSP_fiducial.toml +++ /dev/null @@ -1,281 +0,0 @@ -# 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 = 1.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 = 4 # 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 = 3546.0 # Hydrogen inventory in ppmw relative to mantle mass - - CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system - # C_ppmw = 10.0 # Carbon inventory in ppmw relative to mantle mass - - NH_ratio = 0.0 # 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/fiducial_sub_Neptune_Tang_comparison.toml b/input/planets/fiducial_sub_Neptune_Tang_comparison.toml deleted file mode 100644 index ad1248c31..000000000 --- a/input/planets/fiducial_sub_Neptune_Tang_comparison.toml +++ /dev/null @@ -1,281 +0,0 @@ -# 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 = "fiducial_sub_Neptune" - 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 = true - phi_crit = 0.01 # non-dim., model will terminate when global melt fraction < phi_crit - - # radiative equilibrium - [params.stop.radeqm] - enabled = true - atol = 0.2 # 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 = 10 # 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 = 4.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 = false # 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 = 60.0 # Hydrogen inventory in units of equivalent Earth oceans - H_ppmw = 14180.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" From 5f202b4cd3c214e3fd31049602dc3a3d4d222364 Mon Sep 17 00:00:00 2001 From: rdc49 Date: Tue, 14 Apr 2026 14:01:11 +0100 Subject: [PATCH 3/4] working with dummy atm and no radio --- input/all_options.toml | 105 +--- input/chili/intercomp/_base.grid.toml | 2 +- input/chili/intercomp/_base.toml | 42 +- input/chili/protocol/earth.toml | 2 +- input/chili/protocol/tr1b.toml | 2 +- input/demos/aragog.toml | 2 +- input/demos/escape.toml | 2 +- input/demos/sat_physical.toml | 2 +- input/ensembles/GDSP_initial_grid_2.toml | 36 ++ input/ensembles/boundary_interior_test.toml | 31 + input/ensembles/example.infer.toml | 28 +- input/ensembles/ocean_formation_grid.toml | 31 + input/ensembles/stellar_spectra_test_bb.toml | 29 + .../ensembles/stellar_spectra_test_mors.toml | 29 + input/minimal.toml | 3 +- input/planets/GDSP_fiducial.toml | 281 +++++++++ input/planets/earth.toml | 30 +- input/planets/fiducial_sub_Neptune.toml | 2 +- input/planets/gj9827d.toml | 2 +- input/planets/hd63433d.toml | 2 +- input/planets/l9859b.toml | 2 +- input/planets/l9859c.toml | 2 +- input/planets/l9859d.toml | 2 +- input/planets/trappist1c.toml | 19 +- src/proteus/atmos_clim/dummy.py | 5 +- src/proteus/atmos_clim/wrapper.py | 9 +- src/proteus/config/_interior.py | 110 +++- src/proteus/interior/boundary.py | 568 ++++++++++++++++++ src/proteus/interior/wrapper.py | 43 +- src/proteus/proteus.py | 14 +- 30 files changed, 1258 insertions(+), 179 deletions(-) create mode 100644 input/ensembles/GDSP_initial_grid_2.toml create mode 100644 input/ensembles/boundary_interior_test.toml create mode 100644 input/ensembles/ocean_formation_grid.toml create mode 100644 input/ensembles/stellar_spectra_test_bb.toml create mode 100644 input/ensembles/stellar_spectra_test_mors.toml create mode 100644 input/planets/GDSP_fiducial.toml create mode 100644 src/proteus/interior/boundary.py diff --git a/input/all_options.toml b/input/all_options.toml index 4511c4544..4ee41fe9e 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] @@ -112,22 +110,11 @@ version = "2.0" module = "mors" [star.mors] - rot_pcntle = 50.0 # rotation percentile - rot_period = "none" # rotation period [days] - tracks = "spada" # evolution tracks: spada | baraffe - age_now = 4.567 # current age of star [Gyr] - spectrum_source = "solar" # Spectrum source: 'solar' for solar spectra; 'muscles' for MUSCLES spectra; 'phoenix' for synthetic PHOENIX spectrum; see https://proteus-framework.org/PROTEUS/data.html#stellar-spectra - star_name = "sun" # star name, relevant for when spectrum_source = 'solar' (use e.g. 'sun' or 'Sun0.6Ga') or when spectrum_source = 'muscles' (use e.g. 'trappist-1' or 'gj1214'). Not relevent when spectrum_source = 'phoenix'. - star_path = "none" # optional override star path to custom stellar spectrum - - # PHOENIX parameters, only relevant if spectrum_source = "phoenix". Defaults to solar (0.0). - phoenix_FeH = 0.0 # metallicity [Fe/H] - phoenix_alpha = 0.0 # alpha enhancement [α/M] - - # if None, calculated by mors - 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] + rot_pcntle = 50.0 # rotation percentile + rot_period = "none" # rotation period [days] + tracks = "spada" # evolution tracks: spada | baraffe + age_now = 4.567 # current age of star [Gyr] + spec = "stellar_spectra/Named/sun.txt" # path to stellar spectrum file [star.dummy] radius = 1.0 # Constant stellar radius [R_sun] @@ -137,11 +124,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 = "lovepy" # module used to calculate tidal heating @@ -183,7 +170,7 @@ version = "2.0" relative_tolerance = 1e-5 # relative tolerance for solve_ivp absolute_tolerance = 1e-6 # absolute tolerance for solve_ivp target_surface_pressure = 101325 # target surface pressure - pressure_tolerance = 1e9 # tolerance surface pressure + pressure_tolerance = 1e11 # tolerance surface pressure max_iterations_pressure = 200 # max. iterations for the innermost loop pressure_adjustment_factor = 1.1 # factor for adjusting the pressure in the innermost loop @@ -196,15 +183,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 @@ -216,29 +202,19 @@ version = "2.0" solution_atol = 0.5 # solver absolute tolerance solution_rtol = 0.15 # solver relative tolerance overlap_method = "ee" # gas overlap method - surf_roughness = 1e-3 # characteristic surface roughness [m] - surf_windspeed = 2.0 # characteristic surface wind speed [m/s]. - rainout = true # include volatile condensation/evaporation aloft - latent_heat = false # include latent heat release when `rainout=true`? - oceans = true # form liquid oceans at planet surface? - convection = true # include convective heat transport, with MLT - conduction = true # include conductive heat transport, with Fourier's law - sens_heat = true # include sensible heat flux near surface, with TKE scheme + condensation = false # volatile condensation + latent_heat = false # include latent heat release by phase change? 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 @@ -258,7 +234,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] @@ -268,26 +244,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] @@ -297,26 +253,19 @@ version = "2.0" rheo_phi_loc = 0.6 # Centre of rheological transition rheo_phi_wid = 0.2 # Width of rheological transition melting_dir = "Monteux-600" # Name of folder constaining melting curves - lookup_dir = "1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa" # Name of folder with EOS tables, etc. - module = "aragog" # Which interior module to use [interior.spider] - num_levels = 100 # Number of SPIDER grid levels - mixing_length = 2 # Mixing length parameterization - tolerance = 1.0e-8 # absolute solver tolerance - tolerance_rel = 1.0e-8 # relative solver tolerance - solver_type = "bdf" # SUNDIALS solver method - tsurf_atol = 20.0 # tsurf_poststep_change - tsurf_rtol = 0.02 # 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] - conduction = true # enable conductive heat transfer - convection = true # enable convective heat transfer - gravitational_separation = true # enable gravitational separation - mixing = true # enable mixing - matprop_smooth_width = 1e-2 # melt-fraction window width over which to smooth material properties + num_levels = 100 # Number of SPIDER grid levels + mixing_length = 2 # Mixing length parameterization + tolerance = 1.0e-8 # absolute solver tolerance + tolerance_rel = 1.0e-8 # relative solver tolerance + solver_type = "bdf" # SUNDIALS solver method + tsurf_atol = 20.0 # tsurf_poststep_change + tsurf_rtol = 0.02 # 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] logging = "ERROR" # Aragog log verbosity @@ -442,7 +391,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 a81279d9d..1b4a553c2 100644 --- a/input/chili/intercomp/_base.grid.toml +++ b/input/chili/intercomp/_base.grid.toml @@ -1,7 +1,7 @@ output = "chili_base_grid/" symlink = "" ref_config = "input/chili/intercomp/_base.toml" -use_slurm = true +use_slurm = false 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) diff --git a/input/chili/intercomp/_base.toml b/input/chili/intercomp/_base.toml index 1f7919a4d..de11d8066 100644 --- a/input/chili/intercomp/_base.toml +++ b/input/chili/intercomp/_base.toml @@ -5,24 +5,23 @@ version = "2.0" [params.out] path = "chili_base" logging = "INFO" -plot_fmt = "pdf" -write_mod = 2 -plot_mod = 30 +plot_fmt = "png" +write_mod = 1 +plot_mod = 20 archive_mod = "none" remove_sf = true [params.dt] -starspec = 1e9 +starspec = 10000000.0 starinst = 10.0 method = "adaptive" -minimum = 1000.0 +minimum = 10000.0 minimum_rel = 1e-05 maximum = 10000000.0 initial = 30.0 - [params.dt.adaptive] -atol = 0.02 -rtol = 0.09 +atol = 0.03 +rtol = 0.10 [params.stop] strict = false @@ -98,7 +97,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 @@ -113,27 +112,24 @@ 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 = 50 +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" -fastchem_floor = 273.0 -ini_profile = "isothermal" +mlt_criterion = "l" +fastchem_floor = 100.0 [atmos_clim.dummy] gamma = 0.01 @@ -158,9 +154,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 ad8e18b80..5641b6c2b 100644 --- a/input/chili/protocol/earth.toml +++ b/input/chili/protocol/earth.toml @@ -137,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 diff --git a/input/chili/protocol/tr1b.toml b/input/chili/protocol/tr1b.toml index 880c5ef9c..e58551456 100644 --- a/input/chili/protocol/tr1b.toml +++ b/input/chili/protocol/tr1b.toml @@ -137,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 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 b4a966af6..27178e1d5 100644 --- a/input/demos/escape.toml +++ b/input/demos/escape.toml @@ -95,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/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/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 3c37215f2..f0d128be9 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 b8605b578..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 @@ -60,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 @@ -97,12 +85,6 @@ version = "2.0" [interior.aragog] ini_tmagma = 3000.0 - H_radio = 8e-12 - mantle_cp = 2000.0 - - [interior.spider] - ini_entropy = 2900.0 - num_levels = 60 [outgas] fO2_shift_IW = 4 @@ -114,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 d935cb3f1..f88288257 100644 --- a/input/planets/fiducial_sub_Neptune.toml +++ b/input/planets/fiducial_sub_Neptune.toml @@ -153,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] diff --git a/input/planets/gj9827d.toml b/input/planets/gj9827d.toml index 4c475ca0d..1f372a218 100644 --- a/input/planets/gj9827d.toml +++ b/input/planets/gj9827d.toml @@ -101,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 90a3a915b..28c7109d4 100644 --- a/input/planets/hd63433d.toml +++ b/input/planets/hd63433d.toml @@ -60,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 024c0c3e0..9377a0a8b 100644 --- a/input/planets/l9859b.toml +++ b/input/planets/l9859b.toml @@ -90,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 db69e5ecc..c55b82cd4 100644 --- a/input/planets/l9859c.toml +++ b/input/planets/l9859c.toml @@ -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 41b09cc84..9de5ad7f8 100644 --- a/input/planets/l9859d.toml +++ b/input/planets/l9859d.toml @@ -94,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/trappist1c.toml b/input/planets/trappist1c.toml index 856d84581..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 # ---------------------------------------------------- @@ -51,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] @@ -75,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 0d660894d..9a4296896 100644 --- a/src/proteus/atmos_clim/dummy.py +++ b/src/proteus/atmos_clim/dummy.py @@ -55,7 +55,10 @@ 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 + 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 57d31a4d5..7c4adee87 100644 --- a/src/proteus/atmos_clim/wrapper.py +++ b/src/proteus/atmos_clim/wrapper.py @@ -74,7 +74,11 @@ def run_atmosphere(atmos_o:Atmos_t, config:Config, dirs:dict, loop_counter:dict, hf_row["albedo_pl"] = 0.0 # Handle new surface temperature - if config.atmos_clim.surf_state == 'mixed_layer': + 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': @@ -131,7 +135,8 @@ def run_atmosphere(atmos_o:Atmos_t, config:Config, dirs:dict, loop_counter:dict, 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 574147e58..ce5928899 100644 --- a/src/proteus/config/_interior.py +++ b/src/proteus/config/_interior.py @@ -17,6 +17,20 @@ def valid_spider(instance, attribute, value): if not (spider.conduction or spider.convection or spider.mixing or spider.gravitational_separation): 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") @@ -162,6 +176,93 @@ def valid_interiordummy(instance, attribute, value): if tliq <= tsol: 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. + activation_energy: float + Activation energy for solid viscosity [J/mol]. + solidus_parameterisation: str + Solidus parameterisation to use. + critical_rayleigh_number: float + Critical Rayleigh number for onset of convection [-]. + dynamic_viscosity: float + Reference solid viscosity prefactor [Pa s]. + 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]. + viscosity_activation_temp: float + VFT characteristic temperature for melt viscosity [K]. + viscosity_prefactor: float + VFT prefactor for melt viscosity [Pa s]. + use_aggregate_viscosity: bool + Use aggregate viscosity formulation instead of separate phase viscosities. + 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))) + + # Material constants (defaults match boundary.py) + solidus_parameterisation: str = field(default="Pierru_2022", validator=in_(("dry_peridotite", "ET08", "dry_delta", + "wet_Katz","Monteux2020","Pierru_2022"))) + + activation_energy: float = field(default=3.5e5, validator=gt(0)) # J/mol + critical_rayleigh_number: float = field(default=1.1e3, validator=gt(0)) # - + dynamic_viscosity: float = field(default=3.8e9, validator=gt(0)) # Pa s + 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 + viscosity_activation_temp: float = field(default=4600.0, validator=gt(0)) # K + viscosity_prefactor: float = field(default=2.4e-4, validator=gt(0)) # Pa s + + # Aggregate viscosity parameters + use_aggregate_viscosity: bool = field(default=True) + 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. @@ -228,16 +329,17 @@ class Interior: Set of lookup data files to use in the model (e.g. equations of state). """ - module: str = field(validator=in_(('spider', 'aragog', 'dummy'))) + module: str = field(validator=in_(('spider', 'aragog', 'dummy', 'boundary'))) melting_dir: str = field(default="Monteux-600",validator=valid_path) lookup_dir: str = field(default="1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa", validator=valid_path) 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..6d09e6d63 --- /dev/null +++ b/src/proteus/interior/boundary.py @@ -0,0 +1,568 @@ +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.solidus_parameterisation = config.interior.boundary.solidus_parameterisation + 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 + self.boundary_thickness = config.atmos_clim.surface_d + + # Material constants + self.activation_energy = config.interior.boundary.activation_energy # J/mol + self.critical_rayleigh_number = config.interior.boundary.critical_rayleigh_number # dimensionless + self.dynamic_viscosity = config.interior.boundary.dynamic_viscosity # Pa s + 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 + self.viscosity_activation_temp = config.interior.boundary.viscosity_activation_temp # K + self.viscosity_prefactor = config.interior.boundary.viscosity_prefactor # Pa s + + #atm data + self.gamma = config.atmos_clim.dummy.gamma + self.zenith_angle = config.orbit.zenith_angle + self.albedo_pl = float(hf_row["albedo_pl"]) + self.inst_sf = config.orbit.s0_factor + self.albedo_s = config.atmos_clim.surf_greyalbedo + self.finst = hf_row["F_ins"] + + # Aggregate viscosity parameters + self.use_aggregate_viscosity = config.interior.boundary.use_aggregate_viscosity # boolean + 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 + if self.use_aggregate_viscosity: + eta = self.viscosity_aggregate(phi) + else: + # Original behavior: use separate viscosities + if phi < self.critical_melt_fraction: + eta = self.viscosity_solid(T_p) + else: + eta = self.viscosity_magma_ocean(T_p, 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 76c0b3a48..5ccd49059 100644 --- a/src/proteus/interior/wrapper.py +++ b/src/proteus/interior/wrapper.py @@ -8,6 +8,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.utils.constants import M_earth, R_earth, const_G, element_list from proteus.utils.helper import UpdateStatusfile @@ -38,6 +39,8 @@ def get_nlevb(config:Config): return int(config.interior.spider.num_levels) case "aragog": return int(config.interior.aragog.num_levels) + case "boundary": + return 2 case "dummy": return 2 raise ValueError(f"Invalid interior module selected '{config.interior.module}'") @@ -94,13 +97,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: + + 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") @@ -118,7 +130,8 @@ def determine_interior_radius_with_zalmoxis(dirs:dict, config:Config, hf_all:pd. int_o = Interior_t(get_nlevb(config)) int_o.ic = 1 zalmoxis_solver(config, outdir, hf_row) - 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(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict, outdir:str): ''' @@ -133,6 +146,8 @@ def solve_structure(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict, # 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 + 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 @@ -155,7 +170,7 @@ def solve_structure(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict, def run_interior(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict, - interior_o:Interior_t, verbose:bool=True): + interior_o:Interior_t, atmos_o: Atmos_t=None, verbose:bool=True): """Run interior mantle evolution model. Parameters @@ -170,6 +185,8 @@ def run_interior(dirs:dict, config:Config, 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. """ @@ -199,6 +216,14 @@ def run_interior(dirs:dict, config:Config, 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 39854490e..3b2e3e31b 100644 --- a/src/proteus/proteus.py +++ b/src/proteus/proteus.py @@ -14,6 +14,7 @@ import proteus.utils.archive as archive from proteus.config import read_config_object from proteus.utils.constants import ( + M_earth, vap_list, vol_list, ) @@ -359,12 +360,23 @@ 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() run_interior(self.directories, self.config, - self.hf_all, self.hf_row, self.interior_o) + self.hf_all, self.hf_row, self.interior_o, self.atmos_o) # Advance current time in main loop according to interior step + log.info("Current time: %.2f years"%self.hf_row["Time"]) + log.info("Advancing time by %.2f years"%self.interior_o.dt) self.hf_row["Time"] += self.interior_o.dt # in years self.hf_row["age_star"] += self.interior_o.dt # in years From e986c34fc730b2d7c5b8e480c41dde842f952a3c Mon Sep 17 00:00:00 2001 From: rdc49 Date: Tue, 14 Apr 2026 15:08:32 +0100 Subject: [PATCH 4/4] bug fixes --- input/all_options.toml | 19 +++++++++++++ src/proteus/config/_interior.py | 23 --------------- src/proteus/interior/boundary.py | 48 +++++++++----------------------- src/proteus/interior/wrapper.py | 2 +- 4 files changed, 33 insertions(+), 59 deletions(-) diff --git a/input/all_options.toml b/input/all_options.toml index 9395f699e..d98b803ca 100644 --- a/input/all_options.toml +++ b/input/all_options.toml @@ -336,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)] diff --git a/src/proteus/config/_interior.py b/src/proteus/config/_interior.py index 4d1a2a37d..925f5b4f5 100644 --- a/src/proteus/config/_interior.py +++ b/src/proteus/config/_interior.py @@ -220,14 +220,8 @@ class InteriorBoundary: Mantle liquidus temperature [K]. critical_melt_fraction: float Critical melt fraction for rheological transition. - activation_energy: float - Activation energy for solid viscosity [J/mol]. - solidus_parameterisation: str - Solidus parameterisation to use. critical_rayleigh_number: float Critical Rayleigh number for onset of convection [-]. - dynamic_viscosity: float - Reference solid viscosity prefactor [Pa s]. heat_fusion_silicate: float Latent heat of fusion for silicates [J/kg]. nusselt_exponent: float @@ -240,12 +234,6 @@ class InteriorBoundary: Thermal diffusivity [m^2/s]. thermal_expansivity: float Thermal expansivity [1/K]. - viscosity_activation_temp: float - VFT characteristic temperature for melt viscosity [K]. - viscosity_prefactor: float - VFT prefactor for melt viscosity [Pa s]. - use_aggregate_viscosity: bool - Use aggregate viscosity formulation instead of separate phase viscosities. transition_width: float Width of viscosity transition in melt fraction space [-]. eta_solid_const: float @@ -265,24 +253,15 @@ class InteriorBoundary: critical_melt_fraction: float = field(default=0.4, validator=(gt(0), lt(1))) - # Material constants (defaults match boundary.py) - solidus_parameterisation: str = field(default="Pierru_2022", validator=in_(("dry_peridotite", "ET08", "dry_delta", - "wet_Katz","Monteux2020","Pierru_2022"))) - - activation_energy: float = field(default=3.5e5, validator=gt(0)) # J/mol critical_rayleigh_number: float = field(default=1.1e3, validator=gt(0)) # - - dynamic_viscosity: float = field(default=3.8e9, validator=gt(0)) # Pa s 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 - viscosity_activation_temp: float = field(default=4600.0, validator=gt(0)) # K - viscosity_prefactor: float = field(default=2.4e-4, validator=gt(0)) # Pa s # Aggregate viscosity parameters - use_aggregate_viscosity: bool = field(default=True) 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 @@ -363,8 +342,6 @@ class Interior: module: str = field(validator=in_(('spider', 'aragog', 'dummy', 'boundary'))) melting_dir: str = field(default="Monteux-600",validator=valid_path) - lookup_dir: str = field(default="1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa", - 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) diff --git a/src/proteus/interior/boundary.py b/src/proteus/interior/boundary.py index 6d09e6d63..e5d5f5614 100644 --- a/src/proteus/interior/boundary.py +++ b/src/proteus/interior/boundary.py @@ -90,45 +90,30 @@ def __init__(self, config: Config, dirs: dict, hf_row: dict, hf_all: 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.solidus_parameterisation = config.interior.boundary.solidus_parameterisation 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 - self.boundary_thickness = config.atmos_clim.surface_d # Material constants - self.activation_energy = config.interior.boundary.activation_energy # J/mol self.critical_rayleigh_number = config.interior.boundary.critical_rayleigh_number # dimensionless - self.dynamic_viscosity = config.interior.boundary.dynamic_viscosity # Pa s - 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 - self.viscosity_activation_temp = config.interior.boundary.viscosity_activation_temp # K - self.viscosity_prefactor = config.interior.boundary.viscosity_prefactor # Pa s - - #atm data - self.gamma = config.atmos_clim.dummy.gamma - self.zenith_angle = config.orbit.zenith_angle - self.albedo_pl = float(hf_row["albedo_pl"]) - self.inst_sf = config.orbit.s0_factor - self.albedo_s = config.atmos_clim.surf_greyalbedo - self.finst = hf_row["F_ins"] + 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.use_aggregate_viscosity = config.interior.boundary.use_aggregate_viscosity # boolean 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 + 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 + 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, @@ -191,14 +176,7 @@ def rayleigh_number(self, T_p: float, T_surf: float, phi: float) -> float: Rayleigh number [dimensionless] """ # Determine viscosity based on configuration - if self.use_aggregate_viscosity: - eta = self.viscosity_aggregate(phi) - else: - # Original behavior: use separate viscosities - if phi < self.critical_melt_fraction: - eta = self.viscosity_solid(T_p) - else: - eta = self.viscosity_magma_ocean(T_p, phi) + eta = self.viscosity_aggregate(phi) # Calculate Rayleigh number Ra = ((self.bulk_density * self.surface_gravity * self.thermal_expansivity * diff --git a/src/proteus/interior/wrapper.py b/src/proteus/interior/wrapper.py index 08aacf6db..3a4c8b3dd 100644 --- a/src/proteus/interior/wrapper.py +++ b/src/proteus/interior/wrapper.py @@ -146,7 +146,7 @@ def _resid(x): update_gravity(hf_row) else: - #fill change this once interior module works with zalmoxis + # *** 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)