Merge pull request #272 from FormingWorlds/hn/radius

Set interior structure by: total planet mass, radius

input/aragog.toml

@@ -114,12 +114,11 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-    mass        = 1.0       # M_earth
+    set_by      = 'radius_int' # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 1.0       # M_earth
+    radius_int  = 0.9       # R_earth
     corefrac    = 0.55      # non-dim., radius fraction
 
-    # No module for specifically for internal structure
-    module  = "none"
-
 # Atmosphere - physics table
 [atmos_clim]
     prevent_warming = true      # do not allow the planet to heat up

input/default.toml

@@ -114,13 +114,10 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-
-    # Enclosed mass at the atmosphere-interior interface
-    mass        = 1.0       # M_earth
-    corefrac    = 0.55      # non-dim., radius fraction
-
-    # No module for specifically for internal structure
-    module  = "none"
+    set_by      = 'mass_tot'    # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 1.0           # Total planet mass [M_earth]
+    radius_int  = 0.9           # Radius at mantle-atmosphere boundary [R_earth]
+    corefrac    = 0.55          # non-dim., radius fraction
 
 # Atmosphere - physics table
 [atmos_clim]

input/dummy.toml

@@ -114,12 +114,11 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-    mass        = 1.0       # M_earth
+    set_by      = 'radius_int' # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 1.0       # M_earth
+    radius_int  = 0.9       # R_earth
     corefrac    = 0.55      # non-dim., radius fraction
 
-    # No module for specifically for internal structure
-    module  = "none"
-
 # Atmosphere - physics table
 [atmos_clim]
     prevent_warming = true      # do not allow the planet to heat up

input/hd63433d.toml

@@ -114,12 +114,11 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-    mass        = 1.0       # M_earth
+    set_by      = 'radius_int' # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 1.0       # M_earth
+    radius_int  = 1.073       # R_earth
     corefrac    = 0.55      # non-dim., radius fraction
 
-    # No module for specifically for internal structure
-    module  = "none"
-
 # Atmosphere - physics table
 [atmos_clim]
     prevent_warming = true      # do not allow the planet to heat up

input/k218b.toml

@@ -114,12 +114,11 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-    mass        = 8.63       # M_earth
+    set_by      = 'mass_tot' # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 8.63       # M_earth
+    radius_int  = 1.0       # R_earth
     corefrac    = 0.75      # non-dim., radius fraction
 
-    # No module for specifically for internal structure
-    module  = "none"
-
 # Atmosphere - physics table
 [atmos_clim]
     prevent_warming = true      # do not allow the planet to heat up

input/trappist1c.toml

@@ -114,12 +114,11 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Planetary structure - physics table
 [struct]
-    mass        = 1.308       # M_earth
+    set_by      = 'mass_tot' # Variable to set interior structure: 'radius_int' or 'mass_tot'
+    mass_tot    = 1.308       # M_earth
+    radius_int  = 1.0       # R_earth
     corefrac    = 0.55      # non-dim., radius fraction
 
-    # No module for specifically for internal structure
-    module  = "none"
-
 # Atmosphere - physics table
 [atmos_clim]
     prevent_warming = true      # do not allow the planet to heat up

src/proteus/atmos_clim/agni.py

@@ -363,7 +363,6 @@ def run_agni(atmos, loops_total:int, dirs:dict, config:Config, hf_row:dict):
         easy_start  = False
         dx_max      = config.interior.spider.tsurf_atol+5.0
         ls_increase = 1.02
-        reset_vmrs  = True
         perturb_all = True
         max_steps   = 100
 
@@ -397,7 +396,7 @@ def run_agni(atmos, loops_total:int, dirs:dict, config:Config, hf_row:dict):
                             conv_atol=config.atmos_clim.agni.solution_atol,
                             conv_rtol=config.atmos_clim.agni.solution_rtol,
 
-                            method=1, ls_increase=ls_increase, reset_vmrs=reset_vmrs,
+                            method=1, ls_increase=ls_increase, rainout=True,
                             dx_max=dx_max, ls_method=linesearch, easy_start=easy_start,
                             perturb_all=perturb_all,
 

src/proteus/config/_struct.py

@@ -3,21 +3,24 @@
 from attrs import define, field
 from attrs.validators import gt, in_, lt
 
-from ._converters import none_if_none
-
 
 @define
 class Struct:
     """Planetary structure (mass, radius).
 
-    mass: float
-        Dry mass of the planet's interior [M_earth]
+    set_by: str
+        What defines the interior structure? Choices: "mass_tot", "radius_int".
+    mass_tot: float
+        Total mass of the planet [M_earth]
+    radius_int: float
+        Radius of the atmosphere-mantle boundary [R_earth]
     corefrac: float
         Fraction of the planet's interior radius corresponding to the core.
-    module: str | None
-        Select internal structure module to use. Not used currently.
     """
-    mass: float = field(validator=(gt(0),lt(20)))
-    corefrac: float = field(validator=(gt(0), lt(1)))
 
-    module: str | None = field(validator=in_((None,)), converter=none_if_none)
+    set_by: str = field(validator=in_(('mass_tot','radius_int')))
+
+    mass_tot: float = field(validator=(gt(0),lt(20)))
+    radius_int: float = field(validator=(gt(0),lt(10)))
+
+    corefrac: float = field(validator=(gt(0), lt(1)))

src/proteus/interior/wrapper.py

@@ -7,7 +7,7 @@
 import pandas as pd
 import scipy.optimize as optimise
 
-from proteus.utils.constants import M_earth, R_earth, const_G
+from proteus.utils.constants import M_earth, R_earth, const_G, element_list
 from proteus.utils.helper import UpdateStatusfile
 
 if TYPE_CHECKING:
@@ -28,8 +28,6 @@ def determine_interior_radius(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_
     This uses the interior model's hydrostatic integration to estimate the planet's
     interior mass from a given radius. The radius is then adjusted until the interior mass
     achieves the target mass provided by the user in the config file.
-
-    For the dummy interior, the radius is simply specified by the user in the config file.
     '''
 
     log.info("Using %s interior module to solve strcture"%config.interior.module)
@@ -45,21 +43,23 @@ def determine_interior_radius(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_
     hf_row["gravity"] = 9.81
 
     # Target mass
-    M_target = config.struct.mass * M_earth
+    M_target = config.struct.mass_tot * M_earth
 
-    # We need to solve for the state M_int = config.struct.mass
+    # We need to solve for the state hf_row[M_tot] = config.struct.mass_tot
     # This function takes R_int as the input value, and returns the mass residual
     def _resid(x):
         hf_row["R_int"] = x
 
         log.debug("Try R = %.2e m = %.3f R_earth"%(x,x/R_earth))
 
+        # Use interior model to get dry mass from radius
         run_interior(dirs, config, IC_INTERIOR, hf_all, hf_row, verbose=False)
         if solve_g:
             update_gravity(hf_row)
 
-        res = hf_row["M_int"] - M_target
-        log.debug("    yields M = %.5e kg , resid = %.3e kg"%(hf_row["M_int"], res))
+        # Calculate residual
+        res = hf_row["M_tot"] - M_target
+        log.debug("    yields M = %.5e kg , resid = %.3e kg"%(hf_row["M_tot"], res))
 
         return res
 
@@ -77,14 +77,36 @@ def _resid(x):
                                     x0=hf_row["R_int"], x1=hf_row["R_int"]*1.02)
     hf_row["R_int"] = float(r.root)
     run_interior(dirs, config, IC_INTERIOR, hf_all, hf_row)
-    # update_gravity(hf_row)
+    if solve_g:
+        update_gravity(hf_row)
 
     # Result
     log.info("Found solution for interior structure")
-    log.info("M_int: %.1e kg = %.3f M_earth"%(hf_row["M_int"], hf_row["M_int"]/M_earth))
+    log.info("M_tot: %.1e kg = %.3f M_earth"%(hf_row["M_tot"], hf_row["M_tot"]/M_earth))
     log.info("R_int: %.1e m  = %.3f R_earth"%(hf_row["R_int"], hf_row["R_int"]/R_earth))
     log.info(" ")
 
+def solve_structure(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict):
+    '''
+    Solve for the planet structure based on the method set in the configuration file.
+
+    If the structure is set by the radius, then this is trivial because the radius is used
+    as an input to the interior modules anyway. If the structure is set by mass, then it is
+    solved as an inverse problem for now.
+    '''
+
+    # Trivial case of setting it by the interior radius
+    if config.struct.set_by == 'radius_int':
+        hf_row["R_int"] = config.struct.radius_int * R_earth
+
+    # Set by total mass (mantle + core + volatiles)
+    elif config.struct.set_by == 'mass_tot':
+        determine_interior_radius(dirs, config, hf_all, hf_row)
+
+    # Otherwise, error
+    else:
+        log.error("Invalid constraint on interior structure: %s"%config.struct.set_by)
+
 
 def run_interior(dirs:dict, config:Config, IC_INTERIOR:int,
                  hf_all:pd.DataFrame, hf_row:dict, verbose:bool=True):
@@ -124,6 +146,10 @@ def run_interior(dirs:dict, config:Config, IC_INTERIOR:int,
         if k in hf_row.keys():
             hf_row[k] = output[k]
 
+    # Ensure values are >= 0
+    for k in ("M_mantle","M_mantle_liquid","M_mantle_solid","M_core","Phi_global"):
+        hf_row[k] = max(hf_row[k], 0.0)
+
     # Check that the new temperature is remotely reasonable
     if not (0 < hf_row["T_magma"] < 1e6):
         log.error("T_magma is out of range: %g K"%hf_row["T_magma"])
@@ -133,6 +159,15 @@ def run_interior(dirs:dict, config:Config, IC_INTERIOR:int,
     # Update dry interior mass
     hf_row["M_int"] = hf_row["M_mantle"]  + hf_row["M_core"]
 
+    # Update total planet mass
+    M_volatiles = 0.0
+    for e in element_list:
+        if (e == 'O'):
+            # do not include oxygen, because it varies over time in order to set fO2.
+            continue
+        M_volatiles += hf_row[e+"_kg_total"]
+    hf_row["M_tot"] = hf_row["M_int"] + M_volatiles
+
     # Prevent increasing melt fraction
     if config.atmos_clim.prevent_warming and (IC_INTERIOR==2):
         hf_row["Phi_global"] = min(hf_row["Phi_global"], hf_all.iloc[-1]["Phi_global"])

src/proteus/plot/cpl_population.py

@@ -117,7 +117,8 @@ def plot_population_mass_radius(hf_all:pd.DataFrame, output_dir: str, fwl_dir:st
 
 
 def plot_population_time_density(hf_all:pd.DataFrame, output_dir: str, fwl_dir:str,
-                                 plot_format:str, t0: float=100.0, xmin:float=1e5):
+                                 plot_format:str, t0: float=100.0,
+                                 xmin:float=1e5, ymin:float=0.1, ymax:float=10.0):
     """
     Plot planetary evolution on population time-density diagram
     """
@@ -143,6 +144,8 @@ def plot_population_time_density(hf_all:pd.DataFrame, output_dir: str, fwl_dir:s
     # Make sure that we actually plot the simulation data
     if xmin > np.amax(time):
         xmin = time[0]
+    ymin = min(ymin, np.amin(sim_rho[2:])/2)
+    ymax = max(ymin, np.amax(sim_rho[2:])*2)
 
     # Create plot
     scale = 1.0
@@ -171,7 +174,7 @@ def plot_population_time_density(hf_all:pd.DataFrame, output_dir: str, fwl_dir:s
     # Save figure
     ax.set_ylabel(r"Bulk density [g cm$^{-3}$]")
     ax.set_yscale("log")
-    ax.set_ylim(bottom=0.1, top=10)
+    ax.set_ylim(bottom=ymin, top=ymax)
 
     ax.set_xlabel(r"Time [yr]")
     ax.set_xscale("log")

src/proteus/proteus.py

@@ -10,7 +10,7 @@
 from proteus.atmos_clim import RunAtmosphere
 from proteus.config import read_config_object
 from proteus.escape.wrapper import RunEscape
-from proteus.interior.wrapper import determine_interior_radius, run_interior, update_gravity
+from proteus.interior.wrapper import run_interior, solve_structure, update_gravity
 from proteus.orbit.wrapper import update_period, update_separation
 from proteus.outgas.wrapper import calc_target_elemental_inventories, run_outgassing
 from proteus.star.wrapper import (
@@ -21,7 +21,6 @@
     write_spectrum,
 )
 from proteus.utils.constants import (
-    M_earth,
     gas_list,
     vap_list,
     vol_list,
@@ -180,11 +179,8 @@ def start(self, *, resume: bool = False):
             self.hf_row["F_int"] = self.hf_row["F_atm"]
             self.hf_row["T_eqm"] = 2000.0
 
-            # Planet interior mass
-            self.hf_row["M_int"] = self.config.struct.mass * M_earth
-
-            # Calculate R_int from M_int
-            determine_interior_radius(self.directories, self.config, self.hf_all, self.hf_row)
+            # Solve interior structure
+            solve_structure(self.directories, self.config, self.hf_all, self.hf_row)
 
             # Store partial pressures and list of included volatiles
             inc_gases = []
@@ -263,10 +259,9 @@ def start(self, *, resume: bool = False):
 
             PrintHalfSeparator()
 
-            # Calculate R_int from M_int
+            # Solve interior structure
             if self.hf_row["Time"] < 100:
-                determine_interior_radius(self.directories, self.config,
-                                            self.hf_all, self.hf_row)
+                solve_structure(self.directories, self.config, self.hf_all, self.hf_row)
 
             # Run interior model
             self.dt = run_interior(self.directories, self.config,

src/proteus/utils/coupler.py

@@ -155,9 +155,6 @@ def print_module_configuration(dirs:dict, config:Config, config_path:str):
     if config.interior.module == "spider":
         log.info("  - PETSc         version " + _get_petsc_version())
 
-    # Structure module
-    log.info("Structure module  %s" % config.struct.module)
-
     # Atmosphere module
     write = "Atmos_clim module %s" % config.atmos_clim.module
     match config.atmos_clim.module:
@@ -253,6 +250,9 @@ def GetHelpfileKeys():
             # Dry interior radius (calculated) and mass (from config)
             "R_int", "M_int", # [m], [kg]
 
+            # Total planet mass (mantle + core + volatiles)
+            "M_tot", # [kg]
+
             # Temperatures
             "T_surf", "T_magma", "T_eqm", "T_skin", # all [K]