Aragog constant

docs/installation.md

@@ -107,6 +107,13 @@ You will need to setup Python (>=3.10) on your system. This can be done via brew
         pip install -e CALLIOPE/.
         ```
 
+    - ARAGOG interior model
+
+        ```console
+        git clone git@github.com:ExPlanetology/aragog.git
+        pip install -e aragog/.
+        ```
+
     - ZEPHYRUS escape model
 
         ```console

examples/default/init_coupler.toml

@@ -186,7 +186,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

examples/hd63433d/init_coupler.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

input/default.toml

@@ -186,7 +186,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

input/dummy.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

input/escape.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

input/hd63433d.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

input/jgr_grid.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 280       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]

src/proteus/config/_interior.py

@@ -66,7 +66,10 @@ class Aragog:
 
     Attributes
     ----------
-    some_parameter: str
-        Not used.
+    num_levels: int
+        Number of Aragog grid levels (basic mesh).
+    tolerance: float
+        Solver tolerance.
     """
-    some_parameter: str
+    num_levels: int  = field(validator=ge(40))
+    tolerance: float

src/proteus/interior/aragog.py

@@ -0,0 +1,203 @@
+# Aragog interior module
+from __future__ import annotations
+
+import logging
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pandas as pd
+from aragog import Output, Solver
+from aragog.parser import (
+    Parameters,
+    _BoundaryConditionsParameters,
+    _EnergyParameters,
+    _InitialConditionParameters,
+    _MeshParameters,
+    _PhaseMixedParameters,
+    _PhaseParameters,
+    _Radionuclide,
+    _ScalingsParameters,
+    _SolverParameters,
+)
+
+from proteus.interior.timestep import next_step
+from proteus.utils.constants import R_earth, secs_per_year
+
+if TYPE_CHECKING:
+    from proteus.config import Config
+
+aragog_solver = None
+log = logging.getLogger("fwl."+__name__)
+
+# Run the Aragog interior module
+def RunAragog(config:Config, dirs:dict, IC_INTERIOR:int, hf_row:dict, hf_all:pd.DataFrame):
+
+    global aragog_solver
+
+    log.info("Running Aragog...")
+
+    # Compute time step
+    if IC_INTERIOR==1:
+        dt = 0.0
+    else:
+        step_sf = 1.0 # dt scale factor
+        dt = next_step(config, dirs, hf_row, hf_all, step_sf)
+
+    # Setup Aragog parameters from options at first iteration
+    if (aragog_solver is None):
+        SetupAragogSolver(config, hf_row)
+    # Else only update varying parameters
+    else:
+        UpdateAragogSolver(dt, hf_row)
+
+    # Run Aragog solver
+    aragog_solver.initialize()
+    aragog_solver.solve()
+
+    # Get Aragog output
+    output = GetAragogOutput(hf_row)
+    sim_time = aragog_solver.parameters.solver.end_time
+
+    return sim_time, output
+
+def SetupAragogSolver(config:Config, hf_row:dict):
+
+    global aragog_solver
+
+    scalings = _ScalingsParameters(
+            radius = R_earth, # scaling radius [m]
+            temperature = 4000, # scaling temperature [K]
+            density = 4000, # scaling density
+            time = secs_per_year, # scaling time [sec]
+            )
+
+    solver = _SolverParameters(
+            start_time = 0,
+            end_time = 0,
+            atol = config.interior.aragog.tolerance,
+            rtol = config.interior.aragog.tolerance,
+            )
+
+    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 = 3, # 3 = prescribed temperature
+            inner_boundary_value = 4000, # core temperature [K]
+            emissivity = 1, # only used in gray body BC, outer_boundary_condition = 1
+            equilibrium_temperature = 273, # only used in gray body BC, outer_boundary_condition = 1
+            core_radius = config.struct.corefrac * hf_row["R_int"], # not used now
+            core_density = 10738.332568062382, # not used now
+            core_heat_capacity = 880, # not used now
+            )
+
+    mesh = _MeshParameters(
+            outer_radius = hf_row["R_int"], # planet radius [m]
+            inner_radius = config.struct.corefrac * hf_row["R_int"], # core radius [m]
+            number_of_nodes = config.interior.aragog.num_levels, # basic nodes
+            mixing_length_profile = "constant",
+            surface_density = 4090, # AdamsWilliamsonEOS parameter [kg/m3]
+            gravitational_acceleration = hf_row["gravity"], # [m/s-2]
+            adiabatic_bulk_modulus = 260E9, # AdamsWilliamsonEOS parameter [Pa]
+            )
+
+    energy = _EnergyParameters(
+            conduction = True,
+            convection = True,
+            gravitational_separation = False,
+            mixing = False,
+            radionuclides = False,
+            tidal = False,
+            )
+
+    initial_condition = _InitialConditionParameters(
+            surface_temperature = 4000, # initial top temperature (K)
+            basal_temperature = 4000, # initial bottom temperature (K)
+            )
+
+    phase_liquid = _PhaseParameters(
+            density = 4000,
+            viscosity = 1E2,
+            heat_capacity = 1000,
+            melt_fraction = 1,
+            thermal_conductivity = 4,
+            thermal_expansivity = 1.0E-5,
+            )
+
+    phase_solid = _PhaseParameters(
+            density = 4200,
+            viscosity = 1E21,
+            heat_capacity = 1000,
+            melt_fraction = 0,
+            thermal_conductivity = 4,
+            thermal_expansivity = 1.0E-5,
+            )
+
+    phase_mixed = _PhaseMixedParameters(
+            latent_heat_of_fusion = 4e6,
+            rheological_transition_melt_fraction = 0.4,
+            rheological_transition_width = 0.15,
+            solidus = "aragog/data/test/solidus_1d_lookup.dat",
+            liquidus = "aragog/data/test/liquidus_1d_lookup.dat",
+            phase = "mixed",
+            phase_transition_width = 0.1,
+            grain_size = config.interior.grain_size,
+            )
+
+    radionuclides = _Radionuclide(
+            name = "U235",
+            t0_years = 4.55E9,
+            abundance = 0.0072045,
+            concentration = 0.031,
+            heat_production = 5.68402E-4,
+            half_life_years = 704E6,
+            )
+
+    param = Parameters(
+            boundary_conditions = boundary_conditions,
+            energy = energy,
+            initial_condition = initial_condition,
+            mesh = mesh,
+            phase_solid = phase_solid,
+            phase_liquid = phase_liquid,
+            phase_mixed = phase_mixed,
+            radionuclides = radionuclides,
+            scalings = scalings,
+            solver = solver,
+            )
+
+    aragog_solver = Solver(param)
+
+def UpdateAragogSolver(dt:float, hf_row:dict):
+
+    # Set solver time
+    # hf_row["Time"] is in yr so do not need to scale as long as scaling time is secs_per_year
+    aragog_solver.parameters.solver.start_time = hf_row["Time"]
+    aragog_solver.parameters.solver.end_time = hf_row["Time"] + dt
+
+    # Update initial condition (temperature field from previous run)
+    Tfield: np.array = aragog_solver.temperature_staggered[:, -1]
+    Tfield = Tfield / aragog_solver.parameters.scalings.temperature
+    aragog_solver.parameters.initial_condition.from_field = True
+    aragog_solver.parameters.initial_condition.init_temperature = Tfield
+
+    # Update boundary conditions
+    aragog_solver.parameters.boundary_conditions.outer_boundary_value = hf_row["F_atm"]
+
+    return
+
+def GetAragogOutput(hf_row:dict):
+
+    aragog_output: Output = Output(aragog_solver)
+    output = {}
+
+    output["M_mantle"] = aragog_output.mantle_mass
+    output["T_magma"] = aragog_output.solution_top_temperature
+    output["Phi_global"] = float(aragog_output.melt_fraction_global[-1])
+    output["RF_depth"] = aragog_output.rheological_front
+    output["F_int"] = aragog_output.convective_heat_flux_basic[-1,-1] # Need to be revised for consistency
+
+    output["M_mantle_liquid"] = output["M_mantle"] * output["Phi_global"]
+    output["M_mantle_solid"] = output["M_mantle"] - output["M_mantle_liquid"]
+    output["M_core"] = hf_row["M_core"]
+
+    return output

src/proteus/interior/wrapper.py

@@ -31,8 +31,11 @@ def run_interior(dirs:dict, config:Config, loop_counter:dict, IC_INTERIOR:int, h
         sim_time, output = ReadSPIDER(dirs, config, hf_row["R_int"])
 
     elif config.interior.module == 'aragog':
-        # Not supported
-        raise NotImplementedError("Aragog interface not yet implemented")
+        # Import
+        from proteus.interior.aragog import RunAragog
+
+        # Run Aragog
+        sim_time, output = RunAragog(config, dirs, IC_INTERIOR, hf_row, hf_all)
 
     elif config.interior.module == 'dummy':
         # Import
@@ -53,6 +56,8 @@ def run_interior(dirs:dict, config:Config, loop_counter:dict, IC_INTERIOR:int, h
 
     log.debug("    T_magma    = %.1f K"%hf_row["T_magma"])
     log.debug("    Phi_global = %.3f  "%hf_row["Phi_global"])
+    log.debug("    F_int      = %.2e" %hf_row["F_int"])
+    log.debug("    RF_depth   = %.2e" %hf_row["RF_depth"])
 
     # Time step size
     dt = float(sim_time) - hf_row["Time"]

src/proteus/utils/coupler.py

@@ -255,11 +255,11 @@ def UpdatePlots( hf_all:pd.DataFrame, dirs:dict, config:Config, end=False, num_s
 
     # Check model configuration
     dummy_atm = config.atmos_clim.module == 'dummy'
-    dummy_int = config.interior.module == 'dummy'
+    spider = config.interior.module == 'spider'
     escape    = config.escape.module is not None
 
     # Get all output times
-    if dummy_int:
+    if not spider:
         output_times = []
     else:
         from proteus.interior.spider import get_all_output_times
@@ -279,7 +279,7 @@ def UpdatePlots( hf_all:pd.DataFrame, dirs:dict, config:Config, end=False, num_s
         nc_times = [int(f.split("/")[-1].split("_atm")[0]) for f in ncs]
 
         # Check intersection of atmosphere and interior data
-        if dummy_int:
+        if not spider:
             output_times = nc_times
         else:
             output_times = sorted(list(set(output_times) & set(nc_times)))
@@ -294,7 +294,7 @@ def UpdatePlots( hf_all:pd.DataFrame, dirs:dict, config:Config, end=False, num_s
         log.debug("Snapshots to plot:" + str(plot_times))
 
     # Interior profiles
-    if not dummy_int:
+    if spider:
         jsons = read_jsons(output_dir, plot_times)
         plot_interior(output_dir, plot_times, jsons, config.params.out.plot_fmt)
 
@@ -306,7 +306,7 @@ def UpdatePlots( hf_all:pd.DataFrame, dirs:dict, config:Config, end=False, num_s
         plot_atmosphere(output_dir, plot_times, ncdfs, config.params.out.plot_fmt)
 
         # Atmosphere and interior, stacked
-        if not dummy_int:
+        if spider:
             plot_stacked(output_dir, plot_times, jsons, ncdfs, config.params.out.plot_fmt)
 
         # Flux profiles
@@ -324,7 +324,7 @@ def UpdatePlots( hf_all:pd.DataFrame, dirs:dict, config:Config, end=False, num_s
         plot_sflux(output_dir,          plot_format=config.params.out.plot_fmt)
         plot_sflux_cross(output_dir,    plot_format=config.params.out.plot_fmt)
 
-        if not dummy_int:
+        if spider:
             plot_interior_cmesh(output_dir, plot_format=config.params.out.plot_fmt)
 
         if not dummy_atm:

tests/data/integration/dummy/dummy.toml

@@ -184,7 +184,8 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
-        some_parameter = "some_value"
+        num_levels      = 220       # Number of Aragog grid levels
+        tolerance       = 1.0e-10   # solver tolerance
 
 # Outgassing - physics table
 [outgas]