Merge pull request #291 from FormingWorlds/ls/structure-aragog

Ls/structure aragog

pyproject.toml

@@ -49,7 +49,7 @@ dependencies = [
     "fwl-mors>=24.11.18",
     "fwl-calliope>=24.9.10",
     "fwl-zephyrus>=24.10.15",
-    "aragog==0.1.6a0",
+    "aragog>=0.1.7a0",
     "cmcrameri",
     "juliacall",
     "matplotlib",

src/proteus/config/_struct.py

@@ -16,6 +16,8 @@ class Struct:
         Radius of the atmosphere-mantle boundary [R_earth]
     corefrac: float
         Fraction of the planet's interior radius corresponding to the core.
+    core_density: float
+        Density of the planet's core, optional [kg/m3]
     """
 
     set_by: str = field(validator=in_(('mass_tot','radius_int')))
@@ -24,3 +26,4 @@ class Struct:
     radius_int: float = field(validator=(gt(0),lt(10)))
 
     corefrac: float = field(validator=(gt(0), lt(1)))
+    core_density: float = field(default = 10738.33)

src/proteus/interior/aragog.py

@@ -99,12 +99,12 @@ def SetupAragogSolver(config:Config, hf_row:dict):
     boundary_conditions = _BoundaryConditionsParameters(
             outer_boundary_condition = 4, # 4 = prescribed heat flux
             outer_boundary_value = hf_row["F_atm"], # first guess surface heat flux [W/m2]
-            inner_boundary_condition = 1, # 3 = prescribed temperature
+            inner_boundary_condition = 1, # 1 = core cooling model, 3 = prescribed temperature
             inner_boundary_value = 4000, # core temperature [K], if inner_boundary_condition = 3
             emissivity = 1, # only used in gray body BC, outer_boundary_condition = 1
             equilibrium_temperature = hf_row["T_eqm"], # only used in gray body BC, outer_boundary_condition = 1
-            core_density = 10738.332568062382, # not used now
-            core_heat_capacity = 880, # not used now
+            core_density = config.struct.core_density, # used if inner_boundary_condition = 1
+            core_heat_capacity = 880, # used if inner_boundary_condition = 1
             )
 
     mesh = _MeshParameters(
@@ -253,7 +253,6 @@ def GetAragogOutput(hf_row:dict):
 
     output["M_mantle_liquid"] = output["M_mantle"] * output["Phi_global"]
     output["M_mantle_solid"] = output["M_mantle"] - output["M_mantle_liquid"]
-    output["M_core"] = aragog_output.core_mass
 
     # Calculate surface area
     radii = aragog_output.radii_km_basic * 1e3 # [m]

src/proteus/interior/dummy.py

@@ -8,30 +8,13 @@
 import pandas as pd
 
 from proteus.interior.timestep import next_step
-from proteus.utils.constants import M_earth, R_earth, secs_per_year
+from proteus.utils.constants import secs_per_year
 
 if TYPE_CHECKING:
     from proteus.config import Config
 
 log = logging.getLogger("fwl."+__name__)
 
-def calculate_simple_core_mass(radius:float, corefrac:float):
-    '''
-    Calculate the mass of a planet's core based on the density of Earth's.
-    '''
-
-    earth_fr = 0.55     # earth core radius fraction
-    earth_fm = 0.325    # earth core mass fraction  (https://arxiv.org/pdf/1708.08718.pdf)
-
-    # Estimate core density
-    core_rho = (3.0 * earth_fm * M_earth) / (4.0 * np.pi * ( earth_fr * R_earth )**3.0 )  # core density [kg m-3]
-
-    # Get core mass
-    core_mass = core_rho * 4.0/3.0 * np.pi * (radius * corefrac )**3.0
-
-    return core_mass
-
-
 def calculate_simple_mantle_mass(radius:float, corefrac:float)->float:
     '''
     A very simple interior structure model.
@@ -67,7 +50,6 @@ def RunDummyInt(config:Config, dirs:dict, IC_INTERIOR:int, hf_row:dict, hf_all:p
     output["F_int"] = hf_row["F_atm"]
 
     # Interior structure
-    output["M_core"]   = calculate_simple_core_mass(hf_row["R_int"], config.struct.corefrac)
     output["M_mantle"] = calculate_simple_mantle_mass(hf_row["R_int"], config.struct.corefrac)
 
     # Parameters
@@ -90,7 +72,7 @@ def _calc_phi(tmp:float):
         return  (tmp-tmp_sol)/(tmp_liq-tmp_sol )
 
     # Interior heat capacity [J K-1]
-    cp_int = cp_m*output["M_mantle"] + cp_c*output["M_core"]
+    cp_int = cp_m*output["M_mantle"] + cp_c*hf_row["M_core"]
 
     # Subtract tidal contribution to the total heat flux.
     #    This heat energy is generated only in the mantle, not in the core.

src/proteus/interior/spider.py

@@ -325,7 +325,7 @@ def _try_spider( dirs:dict, config:Config,
 
     # core-mantle boundary condition
     call_sequence.extend(["-CORE_BC",  "1"])
-    call_sequence.extend(["-rho_core", "10738.332568062382"]) # core density
+    call_sequence.extend(["-rho_core", "%.6e"%(config.struct.core_density)]) # core density
     call_sequence.extend(["-cp_core",  "880.0"]) # core heat capacity
 
     # surface boundary condition
@@ -478,7 +478,6 @@ def ReadSPIDER(dirs:dict, config:Config, R_int:float):
     output["M_mantle_liquid"] = float(data_a[0])
     output["M_mantle_solid"]  = float(data_a[1])
     output["M_mantle"]        = float(data_a[2])
-    output["M_core"]          = float(data_a[3])
 
     # Surface properties
     output["T_magma"]         = float(data_a[4])

src/proteus/interior/wrapper.py

@@ -4,6 +4,7 @@
 import logging
 from typing import TYPE_CHECKING
 
+import numpy as np
 import pandas as pd
 import scipy.optimize as optimise
 
@@ -21,6 +22,12 @@ def update_gravity(hf_row:dict):
     '''
     hf_row["gravity"] = const_G * hf_row["M_int"] / (hf_row["R_int"] * hf_row["R_int"])
 
+def calculate_core_mass(hf_row:dict, config:Config):
+    '''
+    Calculate the core mass of the planet.
+    '''
+    hf_row["M_core"] = config.struct.core_density * 4.0/3.0 * np.pi * (hf_row["R_int"] * config.struct.corefrac)**3.0
+
 def determine_interior_radius(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict):
     '''
     Determine the interior radius (R_int) of the planet.
@@ -32,14 +39,10 @@ def determine_interior_radius(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_
 
     log.info("Using %s interior module to solve strcture"%config.interior.module)
 
-    solve_g = True
-    if config.interior.module == 'aragog':
-        log.warning("Cannot solve for interior structure with aragog")
-        solve_g = False
-
     # Initial guess for interior radius and gravity
     IC_INTERIOR = 1
     hf_row["R_int"]   = R_earth
+    calculate_core_mass(hf_row, config)
     hf_row["gravity"] = 9.81
 
     # Target mass
@@ -53,9 +56,9 @@ def _resid(x):
         log.debug("Try R = %.2e m = %.3f R_earth"%(x,x/R_earth))
 
         # Use interior model to get dry mass from radius
+        calculate_core_mass(hf_row, config)
         run_interior(dirs, config, IC_INTERIOR, hf_all, hf_row, verbose=False)
-        if solve_g:
-            update_gravity(hf_row)
+        update_gravity(hf_row)
 
         # Calculate residual
         res = hf_row["M_tot"] - M_target
@@ -76,9 +79,9 @@ def _resid(x):
     r = optimise.root_scalar(_resid, method='secant', xtol=1e3, rtol=rtol, maxiter=12,
                                     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, IC_INTERIOR, hf_all, hf_row)
-    if solve_g:
-        update_gravity(hf_row)
+    update_gravity(hf_row)
 
     # Result
     log.info("Found solution for interior structure")
@@ -95,10 +98,12 @@ def solve_structure(dirs:dict, config:Config, hf_all:pd.DataFrame, hf_row:dict):
     solved as an inverse problem for now.
     '''
 
-    # Trivial case of setting it by the interior radius
+    # Set total mass by radius
+    # We might need here to setup a determine_interior_mass function as mass calculation depends on gravity
     if config.struct.set_by == 'radius_int':
         # radius defines interior structure
         hf_row["R_int"] = config.struct.radius_int * R_earth
+        calculate_core_mass(hf_row, config)
         # initial guess for mass, which will be updated by the interior model
         hf_row["M_int"] = config.struct.mass_tot * M_earth
         update_gravity(hf_row)

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 run_interior, solve_structure, update_gravity
+from proteus.interior.wrapper import run_interior, solve_structure
 from proteus.orbit.wrapper import run_orbit
 from proteus.outgas.wrapper import calc_target_elemental_inventories, run_outgassing
 from proteus.star.wrapper import (
@@ -260,22 +260,15 @@ def start(self, *, resume: bool = False):
                 )
             )
 
-            ############### INTERIOR AND STRUCTURE
+            ############### INTERIOR
 
             PrintHalfSeparator()
 
-            # Solve interior structure
-            if self.loops["total"] < self.loops["init_loops"]:
-                solve_structure(self.directories, self.config, self.hf_all, self.hf_row)
-
             # Run interior model
             self.dt = run_interior(self.directories, self.config,
                                     self.IC_INTERIOR, self.hf_all,  self.hf_row)
 
 
-            # Update surface gravity
-            update_gravity(self.hf_row)
-
             # Advance current time in main loop according to interior step
             self.hf_row["Time"]     += self.dt    # in years
             self.hf_row["age_star"] += self.dt    # in years