API overview

This is a detailed overview of Aragog's API for the user's reference. If you want to understand the underlying model, please visit the model overview.

Directory structure

aragog
├── cfg
│   ├── abe_liquid.cfg
│   ├── abe_mixed.cfg
│   ├── abe_mixed_init.cfg
│   ├── abe_mixed_lookup.cfg
│   ├── abe_solid.cfg
│   └── __init__.py
├── cli.py
├── core.py
├── data.py
├── __init__.py
├── interfaces.py
├── mesh.py
├── output.py
├── parser.py
├── phase.py
├── solver.py
└── utilities.py

Here, the cfg subdirectory contains config files used as input for Aragog, that can be adjusted or added to.

High-level architecture

Aragog is structured around a small number of components:

Configuration and scaling (parser.py)
Mesh and equation-of-state (EOS) setup (mesh.py)
Phase/property evaluators (phase.py, interfaces.py)
State evaluation (fluxes, heating) (solver.py)
Time integration (solver.py → scipy.integrate.solve_ivp)
Postprocessing and export (output.py)

Module reference

aragog.cli

`cli`

`all()`

Download all lookup table data.

Source code in aragog/cli.py

@click.command()
def all():
    """Download all lookup table data."""
    from .data import DownloadLookupTableData

    DownloadLookupTableData()

`cli()`

Source code in aragog/cli.py

@click.group()
def cli():
    pass

`download()`

Download data

Source code in aragog/cli.py

@click.group()
def download():
    """Download data"""
    pass

`env()`

Show environment variables and locations

Source code in aragog/cli.py

@click.command()
def env():
    """Show environment variables and locations"""
    from .data import FWL_DATA_DIR

    click.echo(f"FWL_DATA location: {FWL_DATA_DIR}")

aragog.data

`data`

`FWL_DATA_DIR = Path(os.environ.get('FWL_DATA', platformdirs.user_data_dir('fwl_data')))` `module-attribute`

`basic_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Wolf_Bower+2018')` `module-attribute`

`full_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_400GPa', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Monteux+600', 'Melting_curves/Monteux-600', 'Melting_curves/Wolf_Bower+2018')` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`project_id = 'phsxf'` `module-attribute`

`DownloadLookupTableData(fname='')`

Download lookup table data

Inputs : - fname (optional) : folder name, i.e. "1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018" if not provided download all the folder list

Source code in aragog/data.py

def DownloadLookupTableData(fname: str = ""):
    """
    Download lookup table data

    Inputs :
        - fname (optional) :    folder name, i.e. "1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018"
                                if not provided download all the folder list
    """

    osf = OSF()
    project = osf.project(project_id)
    storage = project.storage("osfstorage")

    data_dir = GetFWLData() / "interior_lookup_tables"
    data_dir.mkdir(parents=True, exist_ok=True)

    # If no folder specified download all basic list
    if not fname:
        folder_list = basic_list
    elif fname in full_list:
        folder_list = [fname]
    else:
        raise ValueError(f"Unrecognised folder name: {fname}")

    for folder in folder_list:
        folder_dir = data_dir / folder
        max_tries = 2 # Maximum download attempts, could be a function argument

        if not folder_dir.exists():
            logger.info(f"Downloading interior lookup table data to {data_dir}")
            for i in range(max_tries):
                logger.info(f"Attempt {i + 1} of {max_tries}")
                success = False

                try:
                    download_zenodo_folder(folder = folder, data_dir=data_dir)
                    success = True
                except RuntimeError as e:
                    logger.error(f"Zenodo download failed: {e}")
                    folder_dir.rmdir()

                if not success:
                    try:
                        download_OSF_folder(storage=storage, folders=folder, data_dir=data_dir)
                        success = True
                    except RuntimeError as e:
                        logger.error(f"OSF download failed: {e}")

                if success:
                    break

                if i < max_tries - 1:
                    logger.info("Retrying download...")
                    sleep(5) # Wait 5 seconds before retrying
                else:
                    logger.error("Max retries reached. Download failed.")

    return

`GetFWLData()`

Get path to FWL data directory on the disk

Source code in aragog/data.py

def GetFWLData() -> Path:
    """
    Get path to FWL data directory on the disk
    """
    return Path(FWL_DATA_DIR).absolute()

`download_OSF_folder(*, storage, folders, data_dir)`

Download a specific folder in the OSF repository

Inputs : - storage : OSF storage name - folders : folder names to download - data_dir : local repository where data are saved

Source code in aragog/data.py

def download_OSF_folder(*, storage, folders: list[str], data_dir: Path):
    """
    Download a specific folder in the OSF repository

    Inputs :
        - storage : OSF storage name
        - folders : folder names to download
        - data_dir : local repository where data are saved
    """
    for file in storage.files:
        for folder in folders:
            if not file.path[1:].startswith(folder):
                continue
            parts = file.path.split("/")[1:]
            target = Path(data_dir, *parts)
            target.parent.mkdir(parents=True, exist_ok=True)
            logger.info(f"Downloading {file.path}...")
            with open(target, "wb") as f:
                file.write_to(f)
            break

`download_zenodo_folder(folder, data_dir)`

Download a specific Zenodo record into specified folder

Inputs : - folder : str Folder name to download - folder_dir : Path local repository where data are saved

Source code in aragog/data.py

def download_zenodo_folder(folder: str, data_dir: Path):
    """
    Download a specific Zenodo record into specified folder

    Inputs :
        - folder : str
            Folder name to download
        - folder_dir : Path
            local repository where data are saved
    """

    folder_dir = data_dir / folder
    folder_dir.mkdir(parents=True)
    zenodo_id = get_zenodo_record(folder)
    cmd = [
            "zenodo_get", zenodo_id,
            "-o", folder_dir
        ]
    out = os.path.join(GetFWLData(), "zenodo.log")
    logger.debug("    logging to %s"%out)
    with open(out,'w') as hdl:
        sp.run(cmd, check=True, stdout=hdl, stderr=hdl)

`get_zenodo_record(folder)`

Get Zenodo record ID for a given folder.

Inputs : - folder : str Folder name to get the Zenodo record ID for

Returns : - str | None : Zenodo record ID or None if not found

Source code in aragog/data.py

def get_zenodo_record(folder: str) -> str | None:
    """
    Get Zenodo record ID for a given folder.

    Inputs :
        - folder : str
            Folder name to get the Zenodo record ID for

    Returns :
        - str | None : Zenodo record ID or None if not found
    """
    zenodo_map = {
        '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018': '15877374',
        '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_400GPa': '15877424',
        '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa': '17417017',
        'Melting_curves/Monteux+600': '15728091',
        'Melting_curves/Monteux-600': '15728138',
        'Melting_curves/Wolf_Bower+2018': '15728072',
    }
    return zenodo_map.get(folder, None)

aragog.core

`core`

Core classes and functions

`logger = logging.getLogger(name)` `module-attribute`

`BoundaryConditions(_parameters, _mesh)` `dataclass`

Boundary conditions

Args: parameters: Parameters mesh: Mesh

`apply_flux_boundary_conditions(state)`

Applies the boundary conditions to the state.

Args: state: The state to apply the boundary conditions to

Source code in aragog/core.py

def apply_flux_boundary_conditions(self, state: State) -> None:
    """Applies the boundary conditions to the state.

    Args:
        state: The state to apply the boundary conditions to
    """
    self.apply_flux_inner_boundary_condition(state)
    self.apply_flux_outer_boundary_condition(state)
    logger.debug("temperature = %s", state.temperature_basic)
    logger.debug("heat_flux = %s", state.heat_flux)

`apply_flux_inner_boundary_condition(state)`

Applies the flux boundary condition to the state at the inner boundary.

Args: state: The state to apply the boundary conditions to

Equivalent to CORE_BC in C code. 1: Simple core cooling 2: Prescribed heat flux 3: Prescribed temperature

Source code in aragog/core.py

def apply_flux_inner_boundary_condition(self, state: State) -> None:
    """Applies the flux boundary condition to the state at the inner boundary.

    Args:
        state: The state to apply the boundary conditions to

    Equivalent to CORE_BC in C code.
        1: Simple core cooling
        2: Prescribed heat flux
        3: Prescribed temperature
    """
    if self._settings.inner_boundary_condition == 1:
        self.core_cooling(state)
    elif self._settings.inner_boundary_condition == 2:
        state.heat_flux[0, :] = self._settings.inner_boundary_value
    elif self._settings.inner_boundary_condition == 3:
        pass
        # raise NotImplementedError
    else:
        msg: str = (
            f"inner_boundary_condition = {self._settings.inner_boundary_condition} is unknown"
        )
        raise ValueError(msg)

`apply_flux_outer_boundary_condition(state)`

Applies the flux boundary condition to the state at the outer boundary.

Args: state: The state to apply the boundary conditions to

Equivalent to SURFACE_BC in C code. 1: Grey-body atmosphere 2: Zahnle steam atmosphere 3: Couple to atmodeller 4: Prescribed heat flux 5: Prescribed temperature

Source code in aragog/core.py

def apply_flux_outer_boundary_condition(self, state: State) -> None:
    """Applies the flux boundary condition to the state at the outer boundary.

    Args:
        state: The state to apply the boundary conditions to

    Equivalent to SURFACE_BC in C code.
        1: Grey-body atmosphere
        2: Zahnle steam atmosphere
        3: Couple to atmodeller
        4: Prescribed heat flux
        5: Prescribed temperature
    """
    if self._settings.outer_boundary_condition == 1:
        self.grey_body(state)
    elif self._settings.outer_boundary_condition == 2:
        raise NotImplementedError
    elif self._settings.outer_boundary_condition == 3:
        msg: str = "Requires coupling to atmodeller"
        logger.error(msg)
        raise NotImplementedError(msg)
    elif self._settings.outer_boundary_condition == 4:
        state.heat_flux[-1, :] = self._settings.outer_boundary_value
    elif self._settings.outer_boundary_condition == 5:
        pass
    else:
        msg: str = (
            f"outer_boundary_condition = {self._settings.outer_boundary_condition} is unknown"
        )
        raise ValueError(msg)

`apply_temperature_boundary_conditions(temperature, temperature_basic, dTdr)`

Conforms the temperature and dTdr at the basic nodes to temperature boundary conditions.

Args: temperature: Temperature at the staggered nodes temperature_basic: Temperature at the basic nodes dTdr: Temperature gradient at the basic nodes

Source code in aragog/core.py

def apply_temperature_boundary_conditions(
    self, temperature: npt.NDArray, temperature_basic: npt.NDArray, dTdr: npt.NDArray
) -> None:
    """Conforms the temperature and dTdr at the basic nodes to temperature boundary conditions.

    Args:
        temperature: Temperature at the staggered nodes
        temperature_basic: Temperature at the basic nodes
        dTdr: Temperature gradient at the basic nodes
    """
    # Core-mantle boundary
    if self._settings.inner_boundary_condition == 3:
        temperature_basic[0, :] = self._settings.inner_boundary_value
        dTdr[0, :] = (
            2 * (temperature[0, :] - temperature_basic[0, :])
            / self._mesh.basic.delta_mesh[0]
            * self._mesh.dxidr[0]
        )
    # Surface
    if self._settings.outer_boundary_condition == 5:
        temperature_basic[-1, :] = self._settings.outer_boundary_value
        dTdr[-1, :] = (
            2 * (temperature_basic[-1, :] - temperature[-1, :])
            / self._mesh.basic.delta_mesh[-1]
            * self._mesh.dxidr[-1]
        )

`apply_temperature_boundary_conditions_melt(melt_fraction, melt_fraction_basic, dphidr)`

Conforms the melt fraction gradient dphidr at the basic nodes to temperature boundary conditions.

Args: melt_fraction: Melt fraction at the staggered nodes melt_fraction_basic: Melt fraction at the basic nodes dphidr: Melt fraction gradient at the basic nodes

Source code in aragog/core.py

def apply_temperature_boundary_conditions_melt(
    self, melt_fraction: npt.NDArray, melt_fraction_basic: npt.NDArray, dphidr: npt.NDArray
) -> None:
    """Conforms the melt fraction gradient dphidr at the basic nodes
       to temperature boundary conditions.

    Args:
        melt_fraction: Melt fraction at the staggered nodes
        melt_fraction_basic: Melt fraction at the basic nodes
        dphidr: Melt fraction gradient at the basic nodes
    """
    # Core-mantle boundary
    if self._settings.inner_boundary_condition == 3:
        dphidr[0, :] = (
            2 * (melt_fraction[0, :] - melt_fraction_basic[0, :])
            / self._mesh.basic.delta_mesh[0]
            * self._mesh.dxidr[0]
        )
    # Surface
    if self._settings.outer_boundary_condition == 5:
        dphidr[-1, :] = (
            2 * (melt_fraction_basic[-1, :] - melt_fraction[-1, :])
            / self._mesh.basic.delta_mesh[-1]
            * self._mesh.dxidr[-1]
        )

`core_cooling(state)`

Applies a core cooling heat flux according to Eq. (37) of Bower et al., 2018

Args: state: The state to apply the boundary condition to

Source code in aragog/core.py

def core_cooling(self, state: State) -> None:
    """Applies a core cooling heat flux according to Eq. (37) of Bower et al., 2018

    Args:
        state: The state to apply the boundary condition to
    """
    core_capacity: float = (
        4
        / 3
        * np.pi
        * np.power(self._mesh.basic.radii[0], 3)
        * self._mesh.settings.core_density
        * self._settings.core_heat_capacity
    )
    cell_capacity = self._mesh.basic.volume[0] * state.capacitance_staggered()[0, :]
    radius_ratio: float = self._mesh.basic.radii[1] / self._mesh.basic.radii[0]
    alpha = np.power(radius_ratio, 2) / ((cell_capacity / (core_capacity * 1.147)) + 1)

    state.heat_flux[0, :] = alpha * state.heat_flux[1, :]

`grey_body(state)`

Applies a grey body flux at the surface.

Args: state: The state to apply the boundary conditions to

Source code in aragog/core.py

def grey_body(self, state: State) -> None:
    """Applies a grey body flux at the surface.

    Args:
        state: The state to apply the boundary conditions to
    """
    state.heat_flux[-1, :] = (
        self._settings.emissivity
        * self._settings.scalings_.stefan_boltzmann_constant
        * (np.power(state.top_temperature, 4) - self._settings.equilibrium_temperature**4)
    )

`InitialCondition(_parameters, _mesh, _phases)` `dataclass`

Initial condition

Args: parameters: Parameters mesh: Mesh phases: PhaseEvaluatorCollection

`get_adiabat(pressure_basic)`

Gets an adiabatic temperature profile by integrating the adiatiabatic temperature gradient dTdPs from the surface. Uses the set surface temperature.

Args: Pressure field on the basic nodes

Returns: Adiabatic temperature profile for the staggered nodes

Source code in aragog/core.py

def get_adiabat(self, pressure_basic) -> npt.NDArray:
    """Gets an adiabatic temperature profile by integrating
       the adiatiabatic temperature gradient dTdPs from the surface.
       Uses the set surface temperature.

    Args:
        Pressure field on the basic nodes

    Returns:
        Adiabatic temperature profile for the staggered nodes
    """

    def adiabat_ode(P,T):
        self._phases.active.set_pressure(P)
        self._phases.active.set_temperature(T)
        self._phases.active.update()
        return self._phases.active.dTdPs()

    # flip the pressure field top to bottom
    pressure_basic = np.flip(pressure_basic)

    sol = solve_ivp(
         adiabat_ode, (pressure_basic[0], pressure_basic[-1]),
         [self._settings.surface_temperature], t_eval=pressure_basic,
         method='RK45', rtol=1e-6, atol=1e-9)

    # flip back the temperature field from bottom to top
    temperature_basic = np.flip(sol.y[0])

    # Return temperature field at staggered nodes
    return self._mesh.quantity_at_staggered_nodes(temperature_basic)

`get_linear()`

Gets a linear temperature profile

Returns: Linear temperature profile for the staggered nodes Only works for uniform spatial mesh.

Source code in aragog/core.py

def get_linear(self) -> npt.NDArray:
    """Gets a linear temperature profile

    Returns:
        Linear temperature profile for the staggered nodes
        Only works for uniform spatial mesh.
    """
    temperature_basic: npt.NDArray = np.linspace(
        self._settings.basal_temperature,
        self._settings.surface_temperature,
        self._mesh.basic.number_of_nodes,
    )
    return self._mesh.quantity_at_staggered_nodes(temperature_basic)

`Mesh(parameters)`

A staggered mesh.

The basic mesh is used for the flux calculations and the staggered mesh is used for the volume calculations.

Args: parameters: Parameters

Source code in aragog/mesh.py

def __init__(self, parameters: Parameters):

    # STEP 1: Set up the basic mesh
    self.settings: _MeshParameters = parameters.mesh
    basic_coordinates: npt.NDArray = self.get_constant_spacing()
    if self.settings.eos_method == 1:
        self.eos = AdamsWilliamsonEOS(
            self.settings, basic_coordinates
        )
    elif self.settings.eos_method == 2:
        self.eos = UserDefinedEOS(
            self.settings, basic_coordinates
        )
    else:
        msg: str = (f"Unknown method to initialize Equation of State")
        raise ValueError(msg)
    if parameters.mesh.mass_coordinates:
        self._planet_density: float = self.get_planet_density(basic_coordinates)
        basic_mass_coordinates: npt.NDArray = (
            self.get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates))
        logger.debug("Basic mass coordinates = %s", basic_mass_coordinates)
    else:
        basic_mass_coordinates = basic_coordinates
    self.basic: FixedMesh = FixedMesh(
        self.settings,
        basic_coordinates,
        basic_mass_coordinates,
        np.max(basic_coordinates),
        np.min(basic_coordinates)
    )

    # STEP 2: Set up the staggered mesh
    staggered_mass_coordinates: npt.NDArray = (
        self.basic.mass_radii[:-1] + 0.5 * self.basic.delta_mesh)
    if parameters.mesh.mass_coordinates:
        staggered_coordinates: npt.NDArray = (
            self.get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates))
    else:
        staggered_coordinates = staggered_mass_coordinates
    self.staggered: FixedMesh = FixedMesh(
        self.settings,
        staggered_coordinates,
        staggered_mass_coordinates,
        self.basic.outer_boundary,
        self.basic.inner_boundary,
    )
    self.eos.set_staggered_pressure(self.staggered.radii)

    # STEP 3: Set up the transform matrices
    if parameters.mesh.mass_coordinates:
        self._dxidr: npt.NDArray = self.get_dxidr_basic()
    else:
        self._dxidr: npt.NDArray = np.ones_like(self.basic.radii)
    self._d_dr_transform: npt.NDArray = self._get_d_dr_transform_matrix()
    self._quantity_transform: npt.NDArray = self._get_quantity_transform_matrix()

`dxidr` `property`

dxi/dr at basic nodes

`d_dr_at_basic_nodes(staggered_quantity)`

Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

Args: staggered_quantity: A quantity defined at the staggered nodes.

Returns: d/dr at the basic nodes

Source code in aragog/mesh.py

def d_dr_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes.

    Returns:
        d/dr at the basic nodes
    """
    d_dr_at_basic_nodes: npt.NDArray = self._d_dr_transform.dot(staggered_quantity)
    logger.debug("d_dr_at_basic_nodes = %s", d_dr_at_basic_nodes)

    return d_dr_at_basic_nodes

`get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates)`

Computes the basic mass coordinates from basic spatial coordinates.

Args: Basic spatial coordinates

Returns: Basic mass coordinates

Source code in aragog/mesh.py

def get_basic_mass_coordinates_from_spatial_coordinates(self, basic_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the basic mass coordinates from basic spatial coordinates.

    Args:
        Basic spatial coordinates

    Returns:
        Basic mass coordinates
    """

    # Set the mass coordinates at the inner boundary from the core mass
    basic_mass_coordinates = np.zeros_like(basic_coordinates)
    basic_mass_coordinates[:,:] = (
        self.settings.core_density / self._planet_density
        * np.power(basic_coordinates[0,:], 3.0)
    )

    # Get mass coordinates by adding individual cell contributions to the mantle mass
    basic_volumes = (np.power(basic_coordinates[1:,:],3.0)
        - np.power(basic_coordinates[:-1,:],3.0))
    for i in range(1, self.settings.number_of_nodes):
        basic_mass_coordinates[i:,:] += (
            self.staggered_effective_density[i-1,:] * basic_volumes[i-1,:] / self._planet_density
        )

    return np.power(basic_mass_coordinates, 1.0/3.0)

`get_constant_spacing()`

Constant radius spacing across the mantle

Returns: Radii with constant spacing as a column vector

Source code in aragog/mesh.py

def get_constant_spacing(self) -> npt.NDArray:
    """Constant radius spacing across the mantle

    Returns:
        Radii with constant spacing as a column vector
    """
    radii: npt.NDArray = np.linspace(
        self.settings.inner_radius, self.settings.outer_radius, self.settings.number_of_nodes
    )
    radii = np.atleast_2d(radii).T

    return radii

`get_dxidr_basic()`

Computes dxidr at basic nodes.

Source code in aragog/mesh.py

def get_dxidr_basic(self) -> npt.NDArray:
    """Computes dxidr at basic nodes."""

    dxidr = (
        self.eos.basic_density
        / self._planet_density
        * np.power(self.basic.radii,2.0)
        / np.power(self.basic.mass_radii,2.0)
    )

    return dxidr

`get_planet_density(basic_coordinates)`

Computes the planet density.

Args: Basic spatial coordinates

Returns: Planet effective density

Source code in aragog/mesh.py

def get_planet_density(self, basic_coordinates: npt.NDArray) -> float:
    """Computes the planet density.

    Args:
        Basic spatial coordinates

    Returns:
        Planet effective density
    """
    core_mass = self.settings.core_density *  np.power(basic_coordinates[0,0], 3.0)
    basic_volumes = (np.power(basic_coordinates[1:,0],3.0)
        - np.power(basic_coordinates[:-1,0],3.0))
    mantle_mass = np.sum(
        self.staggered_effective_density[:,0] * basic_volumes
    )
    planet_density = (
        (core_mass + mantle_mass) / (np.power(basic_coordinates[-1,0],3.0)))
    return planet_density.item()

`get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates)`

Computes the staggered spatial coordinates from staggered mass coordinates.

Args: Staggered mass coordinates

Returns: Staggered spatial coordinates

Source code in aragog/mesh.py

def get_staggered_spatial_coordinates_from_mass_coordinates(self, staggered_mass_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the staggered spatial coordinates from staggered mass coordinates.

    Args:
        Staggered mass coordinates

    Returns:
        Staggered spatial coordinates
    """

    # Initialise the staggered spatial coordinate to the inner boundary
    staggered_coordinates = (np.ones_like(staggered_mass_coordinates)
        * np.power(self.settings.inner_radius,3.0))

    # Add first half cell contribution
    staggered_coordinates += (
        self._planet_density
         * (np.power(staggered_mass_coordinates[0,:],3.0)
        - np.power(self.basic.mass_radii[0,:], 3.0))
        / self.staggered_effective_density[0,:]
    )

    # Get spatial coordinates by adding individual cell contributions to the mantle mass
    shell_effective_density = 0.5*(
        self.staggered_effective_density[1:,:] + self.staggered_effective_density[:-1,:])
    shell_mass_volumes = (np.power(staggered_mass_coordinates[1:,:],3.0)
        - np.power(staggered_mass_coordinates[:-1,:],3.0))
    for i in range(1,self.settings.number_of_nodes-1):
        staggered_coordinates[i:,:] += (
            self._planet_density
            * shell_mass_volumes[i-1,:]
            / shell_effective_density[i-1,:]
        )

    return np.power(staggered_coordinates, 1.0/3.0)

`quantity_at_basic_nodes(staggered_quantity)`

Determines a quantity at the basic nodes that is defined at the staggered nodes.

Uses backward and forward differences at the inner and outer radius, respectively, to obtain the quantity values of the basic nodes at the innermost and outermost nodes. When using temperature boundary conditions, values at outer boundaries will be overwritten. When using flux boundary conditions, values at outer boundaries will be used to provide estimate of individual components of heat fluxes though the total heat flux is imposed.

Args: staggered_quantity: A quantity defined at the staggered nodes

Returns: The quantity at the basic nodes

Source code in aragog/mesh.py

def quantity_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the basic nodes that is defined at the staggered nodes.

    Uses backward and forward differences at the inner and outer radius, respectively, to
    obtain the quantity values of the basic nodes at the innermost and outermost nodes.
    When using temperature boundary conditions, values at outer boundaries will be overwritten.
    When using flux boundary conditions, values at outer boundaries will be used to provide
    estimate of individual components of heat fluxes though the total heat flux is imposed.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes

    Returns:
        The quantity at the basic nodes
    """
    quantity_at_basic_nodes: npt.NDArray = self._quantity_transform.dot(staggered_quantity)
    logger.debug("quantity_at_basic_nodes = %s", quantity_at_basic_nodes)

    return quantity_at_basic_nodes

`quantity_at_staggered_nodes(basic_quantity)`

Determines a quantity at the staggered nodes that is defined at the basic nodes.

Staggered nodes are always located at cell centers, whatever the mesh.

Args: basic_quantity: A quantity defined at the basic nodes

Returns: The quantity at the staggered nodes

Source code in aragog/mesh.py

def quantity_at_staggered_nodes(self, basic_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the staggered nodes that is defined at the basic nodes.

    Staggered nodes are always located at cell centers, whatever the mesh.

    Args:
        basic_quantity: A quantity defined at the basic nodes

    Returns:
        The quantity at the staggered nodes
    """
    quantity_at_staggered_nodes: npt.NDArray = 0.5 * (
        basic_quantity[:-1, ...] + basic_quantity[1:, ...])
    logger.debug("quantity_at_staggered_nodes = %s", quantity_at_staggered_nodes)

    return quantity_at_staggered_nodes

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`PhaseEvaluatorCollection(parameters)` `dataclass`

A collection of phase evaluators

Creates the phase evaluators and selects the active phase based on configuration data.

Args: parameters: Parameters

Attributes: liquid: Liquid evaluator solid: Solid evaluator mixed: Mixed evaluator composite: Composite evaluator active: The active evaluator, which is defined by configuration data

`State(parameters, _evaluator)` `dataclass`

Stores and updates the state at temperature and pressure.

Args: parameters: Parameters evaluator: Evaluator

Attributes: evaluator: Evaluator critical_reynolds_number: Critical Reynolds number gravitational_separation_flux: Gravitational separation flux at the basic nodes heating: Total internal heat production at the staggered nodes (power per unit mass) heating_radio: Radiogenic heat production at the staggered nodes (power per unit mass) heating_tidal: Tidal heat production at the staggered nodes (power per unit mass) heat_flux: Heat flux at the basic nodes (power per unit area) inviscid_regime: True if the flow is inviscid and otherwise False, at the basic nodes inviscid_velocity: Inviscid velocity at the basic nodes is_convective: True if the flow is convecting and otherwise False, at the basic nodes reynolds_number: Reynolds number at the basic nodes super_adiabatic_temperature_gradient: Super adiabatic temperature gradient at the basic nod temperature_basic: Temperature at the basic nodes temperature_staggered: Temperature at the staggered nodes bottom_temperature: Temperature at the bottom basic node top_temperature: Temperature at the top basic node viscous_regime: True if the flow is viscous and otherwise False, at the basic nodes viscous_velocity: Viscous velocity at the basic nodes

`critical_reynolds_number` `property`

Critical Reynolds number from Abe (1993)

`heat_flux` `property`

The total heat flux according to the fluxes specified in the configuration.

`heating` `property`

The power generation according to the heat sources specified in the configuration.

`heating_dilatation` `property`

The heat source through dilation/compression.

`heating_radio` `property`

The radiogenic power generation.

`heating_tidal` `property`

The tidal power generation.

`mass_flux` `property`

The total melt mass flux according to the fluxes specified in the configuration.

`conductive_heat_flux()`

Conductive heat flux:

\[ q_{cond} = -k \frac{\partial T}{\partial r} \]

where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.

Source code in aragog/solver.py

def conductive_heat_flux(self) -> npt.NDArray:
    r"""Conductive heat flux:

    $$
    q_{cond} = -k \frac{\partial T}{\partial r}
    $$

    where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.
    """
    conductive_heat_flux: npt.NDArray = self.phase_basic.thermal_conductivity() * -self.dTdr()

    return conductive_heat_flux

`convective_heat_flux()`

Convective heat flux:

\[ q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r} - \left( \frac{\partial T}{\partial r} \right)_S \right) \]

where $\rho$ is density, $c_p$ is heat capacity at constant pressure, $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and $S$ is entropy.

Source code in aragog/solver.py

def convective_heat_flux(self) -> npt.NDArray:
    r"""Convective heat flux:

    $$
    q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r}
        - \left( \frac{\partial T}{\partial r} \right)_S \right)
    $$ 

    where $\rho$ is density, $c_p$ is heat capacity at constant pressure,
    $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and
    $S$ is entropy.
    """
    convective_heat_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.heat_capacity()
        * self.eddy_diffusivity()
        * -self._super_adiabatic_temperature_gradient
    )

    return convective_heat_flux

`dilatation_heating()`

Dilatation/compression heating (power per unit mass)

Returns: Dilatation/compression heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def dilatation_heating(self) -> npt.NDArray:
    """Dilatation/compression heating (power per unit mass)

    Returns:
        Dilatation/compression heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    mass_flux_staggered: npt.NDArray = self._evaluator.mesh.quantity_at_staggered_nodes(self._mass_flux)

    dilatation_volume_source: npt.NDArray = (
        self.phase_staggered.gravitational_acceleration()
        * self.phase_staggered.delta_specific_volume()
        * mass_flux_staggered
    )

    return dilatation_volume_source

`gravitational_separation_mass_flux()`

Gravitational separation mass flux:

\[ j_{grav} = \rho \phi (1 - \phi) v_{rel} \]

where $\rho$ is density, $\phi$ is melt fraction, and $v_{rel}$ is relative velocity.

Source code in aragog/solver.py

def gravitational_separation_mass_flux(self) -> npt.NDArray:
    r"""Gravitational separation mass flux:

    $$
    j_{grav} = \rho \phi (1 - \phi) v_{rel}
    $$ 

    where $\rho$ is density, $\phi$ is melt fraction, and
    $v_{rel}$ is relative velocity.
    """
    gravitational_separation_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.melt_fraction()
        * (1.0 - self.phase_basic.melt_fraction())
        * self.phase_basic.relative_velocity()
    )
    return gravitational_separation_mass_flux

`mixing_mass_flux()`

Mixing mass flux:

\[ j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r} \]

where $\rho$ is density, $\kappa_h$ is eddy diffusivity, $\phi$ is melt mass fraction, and $r$ is radius.

Source code in aragog/solver.py

def mixing_mass_flux(self) -> npt.NDArray:
    r"""Mixing mass flux:

    $$
    j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r}
    $$

    where $\rho$ is density, $\kappa_h$ is eddy diffusivity,
    $\phi$ is melt mass fraction, and $r$ is radius.
    """
    mixing_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.eddy_diffusivity()
        * -self.dphidr()
    )
    return mixing_mass_flux

`radiogenic_heating(time)`

Radiogenic heating (constant with radius)

Args: time: Time

Returns: Radiogenic heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def radiogenic_heating(self, time: float) -> npt.NDArray:
    """Radiogenic heating (constant with radius)

    Args:
        time: Time

    Returns:
        Radiogenic heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    # Total heat production at a given time (power per unit mass)
    radio_heating_float: float = 0
    for radionuclide in self._evaluator.radionuclides:
        radio_heating_float += radionuclide.get_heating(time)

    # Convert to 1D array (assuming abundances are constant)
    return radio_heating_float * np.ones_like(self.temperature_staggered)

`tidal_heating()`

Tidal heating at each layer of the mantle.

Args: time: Time

Returns: Tidal heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def tidal_heating(self) -> npt.NDArray:
    """Tidal heating at each layer of the mantle.

    Args:
        time: Time

    Returns:
        Tidal heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    length = len(self._settings.tidal_array)

    # Must have correct shape
    if length == 1:
        # scalar
        out = np.ones_like(self.temperature_staggered) * self._settings.tidal_array[0]

    elif length == len(self.temperature_staggered):
        # vector with correct length
        out = np.array([self._settings.tidal_array]).T

    else:
        # vector with invalid length
        raise ValueError(f"Tidal heating array has invalid length {length}")

    return out

`update(temperature, time)`

Updates the state.

The evaluation order matters because we want to minimise the number of evaluations.

Args: temperature: Temperature at the staggered nodes pressure: Pressure at the staggered nodes time: Time

Source code in aragog/solver.py

def update(self, temperature: npt.NDArray, time: FloatOrArray) -> None:
    """Updates the state.

    The evaluation order matters because we want to minimise the number of evaluations.

    Args:
        temperature: Temperature at the staggered nodes
        pressure: Pressure at the staggered nodes
        time: Time
    """
    logger.debug("Updating the state")

    logger.debug("Setting the temperature profile")
    self._temperature_staggered = temperature
    self._temperature_basic = self._evaluator.mesh.quantity_at_basic_nodes(temperature)
    logger.debug("temperature_basic = %s", self.temperature_basic)
    self._dTdr = self._evaluator.mesh.d_dr_at_basic_nodes(temperature)
    logger.debug("dTdr = %s", self.dTdr())
    self._evaluator.boundary_conditions.apply_temperature_boundary_conditions(
        temperature, self._temperature_basic, self.dTdr()
    )

    logger.debug("Setting the melt fraction profile")
    self.phase_staggered.set_temperature(temperature)
    self.phase_staggered.update()
    self.phase_basic.set_temperature(self._temperature_basic)
    self.phase_basic.update()
    phase_to_use = self._evaluator._parameters.phase_mixed.phase
    if phase_to_use == "mixed" or phase_to_use == "composite":
        self._dphidr = self._evaluator.mesh.d_dr_at_basic_nodes(
            self.phase_staggered.melt_fraction()
        )
        logger.debug("dphidr = %s", self.dphidr())
        self._evaluator.boundary_conditions.apply_temperature_boundary_conditions_melt(
            self.phase_staggered.melt_fraction(), self.phase_basic.melt_fraction(), self._dphidr
        )
    else:
        self._dphidr = np.zeros_like(self._dTdr)

    self._super_adiabatic_temperature_gradient = self.dTdr() - self.phase_basic.dTdrs()
    self._is_convective = self._super_adiabatic_temperature_gradient < 0
    velocity_prefactor: npt.NDArray = (
        self.phase_basic.gravitational_acceleration()
        * self.phase_basic.thermal_expansivity()
        * -self._super_adiabatic_temperature_gradient
    )
    # Viscous velocity
    self._viscous_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_cubed
    ) / (18 * self.phase_basic.kinematic_viscosity())
    self._viscous_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    # Inviscid velocity
    self._inviscid_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_squared
    ) / 16
    self._inviscid_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    self._inviscid_velocity[self._is_convective] = np.sqrt(
        self._inviscid_velocity[self._is_convective]
    )
    # Reynolds number
    self._reynolds_number = (
        self._viscous_velocity
        * self._evaluator.mesh.basic.mixing_length
        / self.phase_basic.kinematic_viscosity()
    )
    # Eddy diffusivity
    self._eddy_diffusivity = np.where(
        self.viscous_regime, self._viscous_velocity, self._inviscid_velocity
    )
    self._eddy_diffusivity *= self._evaluator.mesh.basic.mixing_length

    # Heat flux (power per unit area)
    self._heat_flux = np.zeros_like(self.temperature_basic)
    self._mass_flux = np.zeros_like(self.temperature_basic)
    if self._settings.conduction:
        self._heat_flux += self.conductive_heat_flux()
    if self._settings.convection:
        self._heat_flux += self.convective_heat_flux()
    if self._settings.gravitational_separation:
        self._mass_flux += self.gravitational_separation_mass_flux()
    if self._settings.mixing:
        self._mass_flux += self.mixing_mass_flux()
    self._heat_flux += self._mass_flux * self.phase_basic.latent_heat()

    # Heating (power per unit mass)
    self._heating = np.zeros_like(self.temperature_staggered)
    self._heating_radio = np.zeros_like(self.temperature_staggered)
    self._heating_dilatation = np.zeros_like(self.temperature_staggered)
    self._heating_tidal = np.zeros_like(self.temperature_staggered)

    if self._settings.radionuclides:
        self._heating_radio = self.radiogenic_heating(time)
        self._heating += self._heating_radio

    if self._settings.dilatation:
        self._heating_dilatation = self.dilatation_heating()
        self._heating += self._heating_dilatation

    if self._settings.tidal:
        self._heating_tidal = self.tidal_heating()
        self._heating += self._heating_tidal

`_BoundaryConditionsParameters(outer_boundary_condition, outer_boundary_value, inner_boundary_condition, inner_boundary_value, emissivity, equilibrium_temperature, core_heat_capacity)` `dataclass`

Stores parameters in the boundary_conditions section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.equilibrium_temperature /= self.scalings_.temperature
    self.core_heat_capacity /= self.scalings_.heat_capacity
    self._scale_inner_boundary_condition()
    self._scale_outer_boundary_condition()

`_InitialConditionParameters(initial_condition=1, surface_temperature=4000, basal_temperature=4000, init_file='')` `dataclass`

Stores the settings in the initial_condition section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes.

Initial condition method 1: Linear profile 2: User-defined temperature field (from file) 3: Adiabatic profile

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Initial condition method
        1: Linear profile
        2: User-defined temperature field (from file)
        3: Adiabatic profile

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.surface_temperature /= self.scalings_.temperature
    self.basal_temperature /= self.scalings_.temperature

    if self.initial_condition == 2:
        if self.init_file == "":
            msg: str = (f"you must provide an initial temperature file")
            raise ValueError(msg)
        self.init_temperature = np.loadtxt(self.init_file)
        self.init_temperature /= self.scalings_.temperature

aragog.interfaces

`interfaces`

Interfaces

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`MixedPhaseEvaluatorProtocol`

Bases: PhaseEvaluatorProtocol, Protocol

Mixed phase evaluator protocol

`PhaseEvaluatorABC`

Bases: ABC

Phase evaluator ABC

`dTdPs()`

TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2

Source code in aragog/interfaces.py

def dTdPs(self) -> npt.NDArray:
    """TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2"""
    dTdPs: npt.NDArray = (
        self.thermal_expansivity() * self.temperature / (self.density() * self.heat_capacity())
    )

    return dTdPs

`set_pressure(pressure)`

Sets the pressure.

Source code in aragog/interfaces.py

def set_pressure(self, pressure: npt.NDArray) -> None:
    """Sets the pressure."""
    logger.debug("set_pressure = %s", pressure)
    self.pressure = pressure

`set_temperature(temperature)`

Sets the temperature.

Source code in aragog/interfaces.py

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Sets the temperature."""
    logger.debug("set_temperature = %s", temperature)
    self.temperature = temperature

`update()`

Updates quantities to avoid repeat, possibly expensive, calculations.

Source code in aragog/interfaces.py

def update(self) -> None:
    """Updates quantities to avoid repeat, possibly expensive, calculations."""

`PhaseEvaluatorProtocol`

Bases: Protocol

Phase evaluator protocol

`PropertyProtocol`

Bases: Protocol

Property protocol

aragog.mesh

`mesh`

Mesh

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`AdamsWilliamsonEOS(settings, basic_radii)`

Bases: EOS

Adams-Williamson equation of state (EOS).

EOS due to adiabatic self-compression from the definition of the adiabatic bulk modulus:

\[ \left( \frac{{d\rho}}{{dP}} \right)_S = \frac{{\rho}}{{K_S}} \]

where $\rho$ is density, $K_S$ the adiabatic bulk modulus, and $S$ is entropy.

Source code in aragog/mesh.py

def __init__(
    self,
    settings: _MeshParameters,
    basic_radii: npt.NDArray,
):
    self._settings: _MeshParameters = settings
    self._basic_radii: npt.NDArray = basic_radii
    self._outer_boundary = np.max(basic_radii)
    self._inner_boundary = np.min(basic_radii)
    self._surface_density: float = self._settings.surface_density
    self._gravitational_acceleration: float = self._settings.gravitational_acceleration
    self._adiabatic_bulk_modulus: float = self._settings.adiabatic_bulk_modulus
    self._basic_pressure = self.get_pressure_from_radii(basic_radii)
    self._basic_density = self.get_density_from_radii(basic_radii)
    self._staggered_effective_density = self.get_effective_density(basic_radii)

`basic_density` `property`

Density at basic nodes

`basic_pressure` `property`

Pressure at basic nodes

`staggered_effective_density` `property`

Effective density at staggered nodes

`staggered_pressure` `property`

Pressure at staggered nodes

`get_density(pressure)`

Computes density from pressure:

\[ \rho(P) = \rho_s \exp(P/K_S) \]

where $\rho$ is density, $P$ is pressure, $\rho_s$ is surface density, and $K_S$ is adiabatic bulk modulus.

Args: pressure: Pressure

Returns: Density

Source code in aragog/mesh.py

def get_density(self, pressure: FloatOrArray) -> npt.NDArray:
    r"""Computes density from pressure:

    $$
    \rho(P) = \rho_s \exp(P/K_S)
    $$

    where $\rho$ is density, $P$ is pressure, $\rho_s$ is surface density,
    and $K_S$ is adiabatic bulk modulus.

    Args:
        pressure: Pressure

    Returns:
        Density
    """
    density: npt.NDArray = self._surface_density * np.exp(
        pressure / self._adiabatic_bulk_modulus
    )

    return density

`get_density_from_radii(radii)`

Computes density from radii:

$$ \rho(r) = \frac{\rho_s K_S}{K_S + \rho_s g (r-r_s)} $$ where $\rho$ is density, $r$ is radius, $\rho_s$ is surface density, $K_S$ is adiabatic bulk modulus, and $r_s$ is surface radius.

Args: radii: Radii

Returns Density

Source code in aragog/mesh.py

def get_density_from_radii(self, radii: FloatOrArray) -> FloatOrArray:
    r"""Computes density from radii:

    $$
        \rho(r) = \frac{\rho_s K_S}{K_S + \rho_s g (r-r_s)}
    $$
    where $\rho$ is density, $r$ is radius, $\rho_s$ is surface density,
    $K_S$ is adiabatic bulk modulus, and $r_s$ is surface radius.

    Args:
        radii: Radii

    Returns
        Density
    """
    density: FloatOrArray = (self._surface_density * self._adiabatic_bulk_modulus) / (
        self._adiabatic_bulk_modulus
        + self._surface_density
        * self._gravitational_acceleration
        * (radii - self._outer_boundary)
    )

    return density

`get_effective_density(radii)`

Computes effective density using density rho(r) integration over a spherical shell bounded by radii

Args: radii: Radii array

Returns: Effective Density array

Source code in aragog/mesh.py

def get_effective_density(self, radii) -> npt.NDArray:
    r"""
    Computes effective density using density rho(r) integration
    over a spherical shell bounded by radii

    Args:
        radii: Radii array

    Returns:
        Effective Density array
    """

    mass_shell = self.get_mass_within_shell(radii)
    volume_shell = 4 / 3 * np.pi  \
        * (np.power(radii[1:],3.0) - np.power(radii[:-1],3.0))

    return mass_shell/volume_shell

`get_mass_element(radii)`

Computes the mass element:

\[ \frac{\delta m}{\delta r} = 4 \pi r^2 \rho \]

where $\delta m$ is the mass element, $r$ is radius, and $\rho$ is density.

Args: radii: Radii

Returns: The mass element at radii

Source code in aragog/mesh.py

def get_mass_element(self, radii: FloatOrArray) -> npt.NDArray:
    r"""Computes the mass element:

    $$
    \frac{\delta m}{\delta r} = 4 \pi r^2 \rho
    $$

    where $\delta m$ is the mass element, $r$ is radius, and $\rho$ is
    density.

    Args:
        radii: Radii

    Returns:
        The mass element at radii
    """
    mass_element: npt.NDArray = (
        4 * np.pi * np.square(radii) * self.get_density_from_radii(radii)
    )

    return mass_element

`get_mass_within_radii(radii)`

Computes mass within radii:

$$ m(r) = \int 4 \pi r^2 \rho dr $$ where $m$ is mass, $r$ is radius, and $\rho$ is density.

The integral was evaluated using WolframAlpha.

Args: radii: Radii

Returns: Mass within radii

Source code in aragog/mesh.py

def get_mass_within_radii(self, radii: FloatOrArray) -> npt.NDArray:
    r"""Computes mass within radii:

    $$
    m(r) = \int 4 \pi r^2 \rho dr
    $$ 
    where $m$ is mass, $r$ is radius, and $\rho$ is density.

    The integral was evaluated using WolframAlpha.

    Args:
        radii: Radii

    Returns:
        Mass within radii
    """
    a: float = self._surface_density
    b: float = self._adiabatic_bulk_modulus
    c: float = self._gravitational_acceleration
    d: float = self._outer_boundary
    beta: float = b / (a * c) - d

    def mass_integral(radii_: FloatOrArray) -> npt.NDArray:
        """Mass within radii including arbitrary constant of integration.

        Args:
            radii_: Radii

        Returns:
            Mass within radii
        """

        mass: npt.NDArray = (
            4
            * np.pi
            * b
            / c
            * (
                -1.5 * beta * beta
                - beta * radii_
                + 0.5 * radii_ * radii_
                + beta * beta * np.log(abs(beta + radii_))
            )
        )
        # + constant

        return mass

    mass: npt.NDArray = mass_integral(radii) - mass_integral(self._inner_boundary)

    return mass

`get_mass_within_shell(radii)`

Computes the mass within spherical shells bounded by radii.

Args: radii: Radii

Returns: Mass within the bounded spherical shells

Source code in aragog/mesh.py

def get_mass_within_shell(self, radii: npt.NDArray) -> npt.NDArray:
    """Computes the mass within spherical shells bounded by radii.

    Args:
        radii: Radii

    Returns:
        Mass within the bounded spherical shells
    """
    mass: npt.NDArray = self.get_mass_within_radii(radii[1:]) - self.get_mass_within_radii(
        radii[:-1]
    )

    return mass

`get_pressure_from_radii(radii)`

Computes pressure from radii:

\[ P(r) = -K_S \ln \left( 1 + \frac{\rho_s g (r-r_s)}{K_S} \right) \]

where $r$ is radius, $K_S$ is adiabatic bulk modulus, $P$ is pressure, $\rho_s$ is surface density, $g$ is gravitational acceleration, and $r_s$ is surface radius.

Args: radii: Radii

Returns: Pressure

Source code in aragog/mesh.py

def get_pressure_from_radii(self, radii: FloatOrArray) -> npt.NDArray:
    r"""Computes pressure from radii:

    $$
    P(r) = -K_S \ln \left( 1 + \frac{\rho_s g (r-r_s)}{K_S} \right)
    $$

    where $r$ is radius, $K_S$ is adiabatic bulk modulus, $P$ is pressure,
    $\rho_s$ is surface density, $g$ is gravitational acceleration, and
    $r_s$ is surface radius.

    Args:
        radii: Radii

    Returns:
        Pressure
    """
    pressure: npt.NDArray = -self._adiabatic_bulk_modulus * np.log(
        (
            self._adiabatic_bulk_modulus
            + self._surface_density
            * self._gravitational_acceleration
            * (radii - self._outer_boundary)
        )
        / self._adiabatic_bulk_modulus
    )

    return pressure

`get_pressure_gradient(pressure)`

Computes the pressure gradient:

\[ \frac{dP}{dr} = -g \rho \]

where $\rho$ is density, $P$ is pressure, and $g$ is gravitational acceleration.

Args: pressure: Pressure

Returns: Pressure gradient

Source code in aragog/mesh.py

def get_pressure_gradient(self, pressure: FloatOrArray) -> npt.NDArray:
    r"""Computes the pressure gradient:

    $$
    \frac{dP}{dr} = -g \rho
    $$ 

    where $\rho$ is density, $P$ is pressure, and  $g$ is gravitational
    acceleration.

    Args:
        pressure: Pressure

    Returns:
        Pressure gradient
    """
    dPdr: npt.NDArray = -self._gravitational_acceleration * self.get_density(pressure)

    return dPdr

`get_radii_from_pressure(pressure)`

Computes radii from pressure:

\[ P(r) = \int \frac{dP}{dr} dr = \int -g \rho_s \exp(P/K_S) dr \]

And apply the boundary condition $P=0$ at $r=r_s$ to get:

\[ r(P) = \frac{K_s \left( \exp(-P/K_S)-1 \right)}{\rho_s g} + r_s \]

where $r$ is radius, $K_S$ is adiabatic bulk modulus, $P$ is pressure, $\rho_s$ is surface density, $g$ is gravitational acceleration, and $r_s$ is surface radius.

Args: pressure: Pressure

Returns: Radii

Source code in aragog/mesh.py

def get_radii_from_pressure(self, pressure: FloatOrArray) -> npt.NDArray:
    r"""Computes radii from pressure:

    $$
    P(r) = \int \frac{dP}{dr} dr = \int -g \rho_s \exp(P/K_S) dr
    $$ 

    And apply the boundary condition $P=0$ at $r=r_s$ to get:

    $$
    r(P) = \frac{K_s \left( \exp(-P/K_S)-1 \right)}{\rho_s g} + r_s
    $$

    where $r$ is radius, $K_S$ is adiabatic bulk modulus, $P$ is pressure,
    $\rho_s$ is surface density, $g$ is gravitational acceleration, and
    $r_s$ is surface radius.

    Args:
        pressure: Pressure

    Returns:
        Radii
    """
    radii: npt.NDArray = (
        self._adiabatic_bulk_modulus
        * (np.exp(-pressure / self._adiabatic_bulk_modulus) - 1)
        / (self._surface_density * self._gravitational_acceleration)
        + self._outer_boundary
    )

    return radii

`set_staggered_pressure(staggered_radii)`

Set staggered pressure based on staggered radii.

Source code in aragog/mesh.py

def set_staggered_pressure(self, staggered_radii: npt.NDArray,) -> None:
    """Set staggered pressure based on staggered radii."""
    self._staggered_pressure = self.get_pressure_from_radii(staggered_radii)

`EOS`

Bases: ABC

Generic EOS class

`FixedMesh(settings, radii, mass_radii, outer_boundary, inner_boundary)` `dataclass`

A fixed mesh

Args: settings: Mesh parameters radii: Radii of the mesh mass_radii: Mass coordinates of the mesh outer_boundary: Outer boundary for computing depth below the surface inner_boundary: Inner boundary for computing height above the base

Attributes: settings: Mesh parameters radii: Radii of the mesh mass_radii: Mass coordinates of the mesh outer_boundary: Outer boundary for computing depth below the surface inner_boundary: Inner boundary for computing height above the base area: Surface area delta_mesh: Delta radii in mass coordinates depth: Depth below the outer boundary height: Height above the inner boundary mixing_length: Mixing length mixing_length_squared: Mixing length squared mixing_length_cubed: Mixing length cubed number_of_nodes: Number of nodes volume: Volume of the spherical shells defined between neighbouring radii total_volume: Total volume

`area` `cached` `property`

Area

`Mesh(parameters)`

A staggered mesh.

The basic mesh is used for the flux calculations and the staggered mesh is used for the volume calculations.

Args: parameters: Parameters

Source code in aragog/mesh.py

def __init__(self, parameters: Parameters):

    # STEP 1: Set up the basic mesh
    self.settings: _MeshParameters = parameters.mesh
    basic_coordinates: npt.NDArray = self.get_constant_spacing()
    if self.settings.eos_method == 1:
        self.eos = AdamsWilliamsonEOS(
            self.settings, basic_coordinates
        )
    elif self.settings.eos_method == 2:
        self.eos = UserDefinedEOS(
            self.settings, basic_coordinates
        )
    else:
        msg: str = (f"Unknown method to initialize Equation of State")
        raise ValueError(msg)
    if parameters.mesh.mass_coordinates:
        self._planet_density: float = self.get_planet_density(basic_coordinates)
        basic_mass_coordinates: npt.NDArray = (
            self.get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates))
        logger.debug("Basic mass coordinates = %s", basic_mass_coordinates)
    else:
        basic_mass_coordinates = basic_coordinates
    self.basic: FixedMesh = FixedMesh(
        self.settings,
        basic_coordinates,
        basic_mass_coordinates,
        np.max(basic_coordinates),
        np.min(basic_coordinates)
    )

    # STEP 2: Set up the staggered mesh
    staggered_mass_coordinates: npt.NDArray = (
        self.basic.mass_radii[:-1] + 0.5 * self.basic.delta_mesh)
    if parameters.mesh.mass_coordinates:
        staggered_coordinates: npt.NDArray = (
            self.get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates))
    else:
        staggered_coordinates = staggered_mass_coordinates
    self.staggered: FixedMesh = FixedMesh(
        self.settings,
        staggered_coordinates,
        staggered_mass_coordinates,
        self.basic.outer_boundary,
        self.basic.inner_boundary,
    )
    self.eos.set_staggered_pressure(self.staggered.radii)

    # STEP 3: Set up the transform matrices
    if parameters.mesh.mass_coordinates:
        self._dxidr: npt.NDArray = self.get_dxidr_basic()
    else:
        self._dxidr: npt.NDArray = np.ones_like(self.basic.radii)
    self._d_dr_transform: npt.NDArray = self._get_d_dr_transform_matrix()
    self._quantity_transform: npt.NDArray = self._get_quantity_transform_matrix()

`dxidr` `property`

dxi/dr at basic nodes

`d_dr_at_basic_nodes(staggered_quantity)`

Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

Args: staggered_quantity: A quantity defined at the staggered nodes.

Returns: d/dr at the basic nodes

Source code in aragog/mesh.py

def d_dr_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes.

    Returns:
        d/dr at the basic nodes
    """
    d_dr_at_basic_nodes: npt.NDArray = self._d_dr_transform.dot(staggered_quantity)
    logger.debug("d_dr_at_basic_nodes = %s", d_dr_at_basic_nodes)

    return d_dr_at_basic_nodes

`get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates)`

Computes the basic mass coordinates from basic spatial coordinates.

Args: Basic spatial coordinates

Returns: Basic mass coordinates

Source code in aragog/mesh.py

def get_basic_mass_coordinates_from_spatial_coordinates(self, basic_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the basic mass coordinates from basic spatial coordinates.

    Args:
        Basic spatial coordinates

    Returns:
        Basic mass coordinates
    """

    # Set the mass coordinates at the inner boundary from the core mass
    basic_mass_coordinates = np.zeros_like(basic_coordinates)
    basic_mass_coordinates[:,:] = (
        self.settings.core_density / self._planet_density
        * np.power(basic_coordinates[0,:], 3.0)
    )

    # Get mass coordinates by adding individual cell contributions to the mantle mass
    basic_volumes = (np.power(basic_coordinates[1:,:],3.0)
        - np.power(basic_coordinates[:-1,:],3.0))
    for i in range(1, self.settings.number_of_nodes):
        basic_mass_coordinates[i:,:] += (
            self.staggered_effective_density[i-1,:] * basic_volumes[i-1,:] / self._planet_density
        )

    return np.power(basic_mass_coordinates, 1.0/3.0)

`get_constant_spacing()`

Constant radius spacing across the mantle

Returns: Radii with constant spacing as a column vector

Source code in aragog/mesh.py

def get_constant_spacing(self) -> npt.NDArray:
    """Constant radius spacing across the mantle

    Returns:
        Radii with constant spacing as a column vector
    """
    radii: npt.NDArray = np.linspace(
        self.settings.inner_radius, self.settings.outer_radius, self.settings.number_of_nodes
    )
    radii = np.atleast_2d(radii).T

    return radii

`get_dxidr_basic()`

Computes dxidr at basic nodes.

Source code in aragog/mesh.py

def get_dxidr_basic(self) -> npt.NDArray:
    """Computes dxidr at basic nodes."""

    dxidr = (
        self.eos.basic_density
        / self._planet_density
        * np.power(self.basic.radii,2.0)
        / np.power(self.basic.mass_radii,2.0)
    )

    return dxidr

`get_planet_density(basic_coordinates)`

Computes the planet density.

Args: Basic spatial coordinates

Returns: Planet effective density

Source code in aragog/mesh.py

def get_planet_density(self, basic_coordinates: npt.NDArray) -> float:
    """Computes the planet density.

    Args:
        Basic spatial coordinates

    Returns:
        Planet effective density
    """
    core_mass = self.settings.core_density *  np.power(basic_coordinates[0,0], 3.0)
    basic_volumes = (np.power(basic_coordinates[1:,0],3.0)
        - np.power(basic_coordinates[:-1,0],3.0))
    mantle_mass = np.sum(
        self.staggered_effective_density[:,0] * basic_volumes
    )
    planet_density = (
        (core_mass + mantle_mass) / (np.power(basic_coordinates[-1,0],3.0)))
    return planet_density.item()

`get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates)`

Computes the staggered spatial coordinates from staggered mass coordinates.

Args: Staggered mass coordinates

Returns: Staggered spatial coordinates

Source code in aragog/mesh.py

def get_staggered_spatial_coordinates_from_mass_coordinates(self, staggered_mass_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the staggered spatial coordinates from staggered mass coordinates.

    Args:
        Staggered mass coordinates

    Returns:
        Staggered spatial coordinates
    """

    # Initialise the staggered spatial coordinate to the inner boundary
    staggered_coordinates = (np.ones_like(staggered_mass_coordinates)
        * np.power(self.settings.inner_radius,3.0))

    # Add first half cell contribution
    staggered_coordinates += (
        self._planet_density
         * (np.power(staggered_mass_coordinates[0,:],3.0)
        - np.power(self.basic.mass_radii[0,:], 3.0))
        / self.staggered_effective_density[0,:]
    )

    # Get spatial coordinates by adding individual cell contributions to the mantle mass
    shell_effective_density = 0.5*(
        self.staggered_effective_density[1:,:] + self.staggered_effective_density[:-1,:])
    shell_mass_volumes = (np.power(staggered_mass_coordinates[1:,:],3.0)
        - np.power(staggered_mass_coordinates[:-1,:],3.0))
    for i in range(1,self.settings.number_of_nodes-1):
        staggered_coordinates[i:,:] += (
            self._planet_density
            * shell_mass_volumes[i-1,:]
            / shell_effective_density[i-1,:]
        )

    return np.power(staggered_coordinates, 1.0/3.0)

`quantity_at_basic_nodes(staggered_quantity)`

Determines a quantity at the basic nodes that is defined at the staggered nodes.

Uses backward and forward differences at the inner and outer radius, respectively, to obtain the quantity values of the basic nodes at the innermost and outermost nodes. When using temperature boundary conditions, values at outer boundaries will be overwritten. When using flux boundary conditions, values at outer boundaries will be used to provide estimate of individual components of heat fluxes though the total heat flux is imposed.

Args: staggered_quantity: A quantity defined at the staggered nodes

Returns: The quantity at the basic nodes

Source code in aragog/mesh.py

def quantity_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the basic nodes that is defined at the staggered nodes.

    Uses backward and forward differences at the inner and outer radius, respectively, to
    obtain the quantity values of the basic nodes at the innermost and outermost nodes.
    When using temperature boundary conditions, values at outer boundaries will be overwritten.
    When using flux boundary conditions, values at outer boundaries will be used to provide
    estimate of individual components of heat fluxes though the total heat flux is imposed.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes

    Returns:
        The quantity at the basic nodes
    """
    quantity_at_basic_nodes: npt.NDArray = self._quantity_transform.dot(staggered_quantity)
    logger.debug("quantity_at_basic_nodes = %s", quantity_at_basic_nodes)

    return quantity_at_basic_nodes

`quantity_at_staggered_nodes(basic_quantity)`

Determines a quantity at the staggered nodes that is defined at the basic nodes.

Staggered nodes are always located at cell centers, whatever the mesh.

Args: basic_quantity: A quantity defined at the basic nodes

Returns: The quantity at the staggered nodes

Source code in aragog/mesh.py

def quantity_at_staggered_nodes(self, basic_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the staggered nodes that is defined at the basic nodes.

    Staggered nodes are always located at cell centers, whatever the mesh.

    Args:
        basic_quantity: A quantity defined at the basic nodes

    Returns:
        The quantity at the staggered nodes
    """
    quantity_at_staggered_nodes: npt.NDArray = 0.5 * (
        basic_quantity[:-1, ...] + basic_quantity[1:, ...])
    logger.debug("quantity_at_staggered_nodes = %s", quantity_at_staggered_nodes)

    return quantity_at_staggered_nodes

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`UserDefinedEOS(settings, basic_radii)`

Bases: EOS

User defined equation of state (EOS).

Pressure field and effective density field on staggered nodes provided by the user.

Source code in aragog/mesh.py

def __init__(
    self,
    settings: _MeshParameters,
    basic_radii: npt.NDArray,
):
    self._interp_pressure = PchipInterpolator(settings.eos_radius, settings.eos_pressure)
    self._interp_density = PchipInterpolator(settings.eos_radius, settings.eos_density)
    self._basic_pressure = self._interp_pressure(basic_radii).reshape(-1,1)
    basic_effective_density = self._interp_density(basic_radii).reshape(-1,1)
    self._staggered_effective_density = 0.5*(
        basic_effective_density[:-1, :] + basic_effective_density[1:, :])
    # Assumes density and effective density are the same at basic nodes
    self._basic_density = basic_effective_density

`basic_density` `property`

Effective density at basic nodes

`basic_pressure` `property`

Pressure at basic nodes

`staggered_effective_density` `property`

Effective density at staggered nodes

`staggered_pressure` `property`

Pressure at staggered nodes

`set_staggered_pressure(staggered_radii)`

Set staggered pressure based on staggered radii.

Source code in aragog/mesh.py

def set_staggered_pressure(self, staggered_radii: npt.NDArray,) -> None:
    """Set staggered pressure based on staggered radii."""
    self._staggered_pressure = self._interp_pressure(staggered_radii).reshape(-1,1)

`_MeshParameters(outer_radius, inner_radius, number_of_nodes, mixing_length_profile, core_density, eos_method=1, surface_density=4000, gravitational_acceleration=9.81, adiabatic_bulk_modulus=260000000000.0, mass_coordinates=False, eos_file='')` `dataclass`

Stores parameters in the mesh section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.outer_radius /= self.scalings_.radius
    self.inner_radius /= self.scalings_.radius
    self.core_density /= self.scalings_.density
    self.surface_density /= self.scalings_.density
    self.gravitational_acceleration /= self.scalings_.gravitational_acceleration
    self.adiabatic_bulk_modulus /= self.scalings_.pressure

    if self.eos_method == 2:
        if self.eos_file == "":
            msg: str = (f"you must provide a file for setting up equation of state")
            raise ValueError(msg)
        arr = np.loadtxt(self.eos_file)
        self.eos_radius = arr[:,0] / self.scalings_.radius
        self.eos_pressure = arr[:,1] / self.scalings_.pressure
        self.eos_density = arr[:,2] / self.scalings_.density
        self.eos_gravity = arr[:,3] / self.scalings_.gravitational_acceleration
        # Check that provided eos radius roughly match with Aragog mesh
        if ((self.eos_radius[0] < self.inner_radius) or
            (self.eos_radius[-1] > self.outer_radius) or
            (self.eos_radius[-1]-self.eos_radius[0]) < 0.75*(self.outer_radius-self.inner_radius)):
            msg: str = (f"Radius array in EOS file: Values out of range.")
            raise ValueError(msg)

`is_monotonic_increasing(some_array)`

Returns True if an array is monotonically increasing, otherwise returns False.

Source code in aragog/utilities.py

def is_monotonic_increasing(some_array: npt.NDArray) -> bool:
    """Returns True if an array is monotonically increasing, otherwise returns False."""
    return np.all(np.diff(some_array) > 0)

aragog.output

`output`

Output

`FloatOrArray = float | npt.NDArray` `module-attribute`

`version = '26.01.06'` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`Evaluator(_parameters)` `dataclass`

Contains classes that evaluate quantities necessary to compute the interior evolution.

Args: _parameters: Parameters

Attributes: boundary_conditions: Boundary conditions initial_condition: Initial condition mesh: Mesh phases: Evaluators for all phases radionuclides: Radionuclides

`Output(solver)`

Stores inputs and outputs of the models.

Source code in aragog/output.py

def __init__(self, solver: Solver):
    self.solver: Solver = solver
    self.parameters: Parameters = solver.parameters
    self.solution: OptimizeResult = self.solver.solution
    self.evaluator: Evaluator = self.solver.evaluator
    self.state: State = self.solver.state

`conductive_heat_flux_basic` `property`

Conductive heat flux

`convective_heat_flux_basic` `property`

Convective heat flux

`core_mass` `property`

Core mass computed with constant density

`dTdr` `property`

dTdr

`dTdrs` `property`

dTdrs

`density_basic` `property`

Density

`gravitational_separation_heat_flux_basic` `property`

Gravitational separation heat flux

`heat_capacity_basic` `property`

Heat capacity

`heating` `property`

Internal heat generation at staggered nodes

`heating_dilatation` `property`

Internal heat generation from dilatation/compression at staggered nodes

`heating_radio` `property`

Internal heat generation from radioactive decay at staggered nodes

`heating_tidal` `property`

Internal heat generation from tidal heat dissipation at staggered nodes

`liquidus_K_staggered` `property`

Liquidus

`log10_viscosity_basic` `property`

Viscosity of the basic mesh

`log10_viscosity_staggered` `property`

Viscosity of the staggered mesh

`mantle_mass` `property`

Mantle mass computed from the AdamsWilliamsonEOS

`mass_radii_km_basic` `property`

Mass radii of the basic mesh in km

`mass_radii_km_staggered` `property`

Mass radii of the staggered mesh in km

`mass_staggered` `property`

Mass of each layer on staggered mesh

`melt_fraction_basic` `property`

Melt fraction on the basic mesh

`melt_fraction_global` `property`

Volume-averaged melt fraction

`melt_fraction_staggered` `property`

Melt fraction on the staggered mesh

`mixing_heat_flux_basic` `property`

Convective mixing heat flux

`pressure_GPa_basic` `property`

Pressure of the basic mesh in GPa

`pressure_GPa_staggered` `property`

Pressure of the staggered mesh in GPa

`radii_km_basic` `property`

Radii of the basic mesh in km

`radii_km_staggered` `property`

Radii of the staggered mesh in km

`rheological_front` `property`

Rheological front at the last solve iteration given user defined threshold. It is defined as a dimensionless distance with respect to the outer radius.

`shape_basic` `property`

Shape of the basic data

`shape_staggered` `property`

Shape of the staggered data

`solidus_K_staggered` `property`

Solidus

`solution_top_temperature` `property`

Solution (last iteration) temperature at the top of the domain (planet surface)

`super_adiabatic_temperature_gradient_basic` `property`

Super adiabatic temperature gradient

`temperature_K_basic` `property`

Temperature of the basic mesh in K

`temperature_K_staggered` `property`

Temperature of the staggered mesh in K

`thermal_expansivity_basic` `property`

Thermal expansivity

`times` `property`

Times in years

`total_heat_flux_basic` `property`

Conductive heat flux

`plot(num_lines=11, figsize=(25, 10))`

Plots the solution with labelled lines according to time.

Args: num_lines: Number of lines to plot. Defaults to 11. figsize: Size of the figure. Defaults to (25, 10).

Source code in aragog/output.py

def plot(self, num_lines: int = 11, figsize: tuple = (25, 10)) -> None:
    """Plots the solution with labelled lines according to time.

    Args:
        num_lines: Number of lines to plot. Defaults to 11.
        figsize: Size of the figure. Defaults to (25, 10).
    """
    assert self.solution is not None

    self.state.update(self.solution.y, self.solution.t)

    _, axs = plt.subplots(1, 10, sharey=True, figsize=figsize)

    # Ensure there are at least 2 lines to plot (first and last).
    num_lines = max(2, num_lines)

    # Calculate the time step based on the total number of lines.
    time_step: float = self.time_range / (num_lines - 1)

    # plot temperature
    try:
        axs[0].scatter(self.liquidus_K_staggered, self.pressure_GPa_staggered)
        axs[0].scatter(self.solidus_K_staggered, self.pressure_GPa_staggered)
    except AttributeError:
        pass

    # Plot the first line.
    def plot_times(ax, x: npt.NDArray, y: npt.NDArray) -> None:
        # If `x` is a float, create an array of the same length as `y` filled with `x`
        if np.isscalar(x):
            x = np.full((len(y), len(self.times)), x, dtype=np.float64)
        elif x.ndim == 1:
            # If `x` is 1D, reshape it or repeat it across the time dimension
            x = np.tile(x[:, np.newaxis], (1, len(self.times)))

        label_first: str = f"{self.times[0]:.2f}"
        ax.plot(x[:, 0], y, label=label_first)

        # Loop through the selected lines and plot each with a label.
        for i in range(1, num_lines - 1):
            desired_time: float = self.times[0] + i * time_step
            # Find the closest available time step.
            closest_time_index: np.intp = np.argmin(np.abs(self.times - desired_time))
            time: float = self.times[closest_time_index]
            label: str = f"{time:.2f}"  # Create a label based on the time.
            ax.plot(
                x[:, closest_time_index],
                y,
                label=label,
            )

        # Plot the last line.
        times_end: float = self.times[-1]
        label_last: str = f"{times_end:.2f}"
        ax.plot(x[:, -1], y, label=label_last)

    plot_times(axs[0], self.temperature_K_basic, self.pressure_GPa_basic)
    axs[0].set_ylabel("Pressure (GPa)")
    axs[0].set_xlabel("Temperature (K)")
    axs[0].set_title("Temperature")

    plot_times(axs[1], self.melt_fraction_staggered, self.pressure_GPa_staggered)
    axs[1].set_xlabel("Melt fraction")
    axs[1].set_title("Melt fraction")

    plot_times(axs[2], self.log10_viscosity_basic, self.pressure_GPa_basic)
    axs[2].set_xlabel("Log10(viscosity)")
    axs[2].set_title("Log10(viscosity)")

    plot_times(axs[3], self.convective_heat_flux_basic, self.pressure_GPa_basic)
    axs[3].set_xlabel("Convective heat flux")
    axs[3].set_title("Convective heat flux")

    plot_times(
        axs[4],
        self.super_adiabatic_temperature_gradient_basic,
        self.pressure_GPa_basic,
    )
    axs[4].set_xlabel("Super adiabatic temperature gradient")
    axs[4].set_title("Super adiabatic temperature gradient")

    plot_times(
        axs[5],
        self.dTdr,
        self.pressure_GPa_basic,
    )
    axs[5].set_xlabel("dTdr")
    axs[5].set_title("dTdr")

    plot_times(
        axs[6],
        self.dTdrs,
        self.pressure_GPa_basic,
    )
    axs[6].set_xlabel("dTdrs")
    axs[6].set_title("dTdrs")

    plot_times(
        axs[7],
        self.density_basic,
        self.pressure_GPa_basic,
    )
    axs[7].set_xlabel("Density")
    axs[7].set_title("Density")

    plot_times(
        axs[8],
        self.heat_capacity_basic,
        self.pressure_GPa_basic,
    )
    axs[8].set_xlabel("Heat capacity")
    axs[8].set_title("Heat capacity")

    plot_times(
        axs[9],
        self.thermal_expansivity_basic,
        self.pressure_GPa_basic,
    )
    axs[9].set_xlabel("Thermal expansivity")
    axs[9].set_title("Thermal expansivity")

    # Shrink current axis by 20% to allow space for the legend.
    # box = ax.get_position()
    # ax.set_position((box.x0, box.y0, box.width * 0.8, box.height))

    # axs[0].grid(True)

    legend = axs[2].legend()  # (loc="center left", bbox_to_anchor=(1, 0.5))
    legend.set_title("Time (yr)")
    plt.gca().invert_yaxis()

    plt.show()

`write_at_time(file_path, tidx=-1, compress=False)`

Write the state of the model at a particular time to a NetCDF4 file on the disk.

Args: file_path: Path to the output file tidx: Index on the time axis at which to access the data compress: Whether to compress the data

Source code in aragog/output.py

def write_at_time(self, file_path: str, tidx: int = -1, compress: bool = False) -> None:
    """Write the state of the model at a particular time to a NetCDF4 file on the disk.

    Args:
        file_path: Path to the output file
        tidx: Index on the time axis at which to access the data
        compress: Whether to compress the data
    """

    logger.debug("Writing i=%d NetCDF file to %s", tidx, file_path)

    # Update the state
    assert self.solution is not None
    self.state.update(self.solution.y, self.solution.t)

    # Open the dataset
    ds: nc.Dataset = nc.Dataset(file_path, mode="w")

    # Metadata
    ds.description = "Aragog output data"
    ds.argog_version = __version__

    # Function to save scalar quantities
    def _add_scalar_variable(key: str, value: float, units: str):
        ds.createVariable(key, np.float64)
        ds[key][0] = value.item() if hasattr(value, 'item') else value
        ds[key].units = units

    # Save scalar quantities
    _add_scalar_variable("time", self.times[tidx], "yr")
    _add_scalar_variable("phi_global", self.melt_fraction_global, "")
    _add_scalar_variable("mantle_mass", self.mantle_mass, "kg")
    _add_scalar_variable("rheo_front", self.rheological_front, "")

    # Create dimensions for mesh quantities
    ds.createDimension("basic", self.shape_basic[0])
    ds.createDimension("staggered", self.shape_staggered[0])

    # Function to save mesh quantities
    def _add_mesh_variable(key: str, some_property: Any, units: str):
        if key[-2:] == "_b":
            mesh = "basic"
        elif key[-2:] == "_s":
            mesh = "staggered"
        else:
            raise KeyError(f"NetCDF variable name must end in _b or _s: {key}")
        ds.createVariable(
            key,
            np.float64,
            (mesh,),
            compression="zlib" if compress else None,
            shuffle = True,
            complevel = 4
        )
        ds[key][:] = some_property[:, tidx]
        ds[key].units = units

    # Save mesh quantities
    _add_mesh_variable("radius_b", self.radii_km_basic, "km")
    _add_mesh_variable("pres_b", self.pressure_GPa_basic, "GPa")
    _add_mesh_variable("temp_b", self.temperature_K_basic, "K")
    _add_mesh_variable("phi_b", self.melt_fraction_basic, "")
    _add_mesh_variable("Fcond_b", self.conductive_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fconv_b", self.convective_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fgrav_b", self.gravitational_separation_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fmix_b", self.mixing_heat_flux_basic, "W m-2")
    _add_mesh_variable("Ftotal_b", self.total_heat_flux_basic, "W m-2")
    _add_mesh_variable("log10visc_b", self.log10_viscosity_basic, "Pa s")
    _add_mesh_variable("log10visc_s", self.log10_viscosity_staggered, "Pa s")
    _add_mesh_variable("density_b", self.density_basic, "kg m-3")
    _add_mesh_variable("heatcap_b", self.heat_capacity_basic, "J kg-1 K-1")
    _add_mesh_variable("mass_s", self.mass_staggered, "kg")
    _add_mesh_variable("Hradio_s", self.heating_radio, "W kg-1")
    _add_mesh_variable("Hvol_s", self.heating_dilatation, "W kg-1")
    _add_mesh_variable("Htidal_s", self.heating_tidal, "W kg-1")
    _add_mesh_variable("Htotal_s", self.heating, "W kg-1")

    # Close the dataset
    ds.close()

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`Solver(param)`

Solves the interior dynamics

Args: filename: Filename of a file with configuration settings root: Root path to the flename

Attributes: filename: Filename of a file with configuration settings root: Root path to the filename. Defaults to empty parameters: Parameters evaluator: Evaluator state: State

Source code in aragog/solver.py

def __init__(self, param: Parameters):
    logger.info("Creating an Aragog model")
    self.parameters: Parameters = param
    self.evaluator: Evaluator
    self.state: State
    self._solution: OptimizeResult

`solution` `property`

The solution.

`temperature_basic` `property`

Temperature of the basic mesh in K

`temperature_staggered` `property`

Temperature of the staggered mesh in K

`dTdt(time, temperature)`

dT/dt at the staggered nodes

Args: time: Time temperature: Temperature at the staggered nodes

Returns: dT/dt at the staggered nodes

Source code in aragog/solver.py

def dTdt(
    self,
    time: npt.NDArray | float,
    temperature: npt.NDArray,
) -> npt.NDArray:
    """dT/dt at the staggered nodes

    Args:
        time: Time
        temperature: Temperature at the staggered nodes

    Returns:
        dT/dt at the staggered nodes
    """
    logger.debug("temperature passed into dTdt = %s", temperature)
    # logger.debug("temperature.shape = %s", temperature.shape)
    self.state.update(temperature, time)
    heat_flux: npt.NDArray = self.state.heat_flux
    # logger.debug("heat_flux = %s", heat_flux)
    self.evaluator.boundary_conditions.apply_flux_boundary_conditions(self.state)
    # logger.debug("heat_flux = %s", heat_flux)
    # logger.debug("mesh.basic.area.shape = %s", self.data.mesh.basic.area.shape)

    energy_flux: npt.NDArray = heat_flux * self.evaluator.mesh.basic.area
    # logger.debug("energy_flux size = %s", energy_flux.shape)

    delta_energy_flux: npt.NDArray = np.diff(energy_flux, axis=0)
    # logger.debug("delta_energy_flux size = %s", delta_energy_flux.shape)
    # logger.debug("capacitance = %s", self.state.phase_staggered.capacitance.shape)
    # FIXME: Update capacitance for mixed phase (enthalpy of fusion contribution)
    capacitance: npt.NDArray = (
        self.state.capacitance_staggered() * self.evaluator.mesh.basic.volume
    )

    # Heating rate (dT/dt) from flux divergence (power per unit area)
    dTdt: npt.NDArray = -delta_energy_flux / capacitance
    logger.debug("dTdt (fluxes only) = %s", dTdt)

    # Additional heating rate (dT/dt) from internal heating (power per unit mass)
    dTdt += self.state.heating * (
        self.state.phase_staggered.density() / self.state.capacitance_staggered()
    )

    logger.debug("dTdt (with internal heating) = %s", dTdt)

    return dTdt

`from_file(filename, root=Path())` `classmethod`

Parses a configuration file

Args: filename: Filename root: Root of the filename

Returns: Parameters

Source code in aragog/solver.py

@classmethod
def from_file(cls, filename: str | Path, root: str | Path = Path()) -> Self:
    """Parses a configuration file

    Args:
        filename: Filename
        root: Root of the filename

    Returns:
        Parameters
    """
    configuration_file: Path = Path(root) / Path(filename)
    logger.info("Parsing configuration file = %s", configuration_file)
    parameters: Parameters = Parameters.from_file(configuration_file)

    return cls(parameters)

`initialize()`

Initializes the model.

Source code in aragog/solver.py

def initialize(self) -> None:
    """Initializes the model."""
    logger.info("Initializing %s", self.__class__.__name__)
    self.evaluator = Evaluator(self.parameters)
    self.state = State(self.parameters, self.evaluator)

`make_tsurf_event()`

Creates a temperature event function for use with an ODE solver to monitor changes in the surface temperature.The event triggers when the change exceeds the threshold, allowing the solver to stop integration.

Returns: The event has the attributes: - terminal = True: Integration stops when the event is triggered. - direction = -1: Only triggers when the function is decreasing through zero.

Source code in aragog/solver.py

def make_tsurf_event(self):
    """
    Creates a temperature event function for use with an ODE solver to monitor changes 
    in the surface temperature.The event triggers when the change exceeds the 
    threshold, allowing the solver to stop integration.

    Returns:
        The event has the attributes:
            - terminal = True: Integration stops when the event is triggered.
            - direction = -1: Only triggers when the function is decreasing through zero.
    """
    tsurf_initial = [None]

    def tsurf_event(time: float, temperature: npt.NDArray) -> float:
        """
        Event function to detect when surface temperature changes beyond a specified threshold.

        Args:
            time (float): Current time.
            temperature (np.ndarray): Current temperature profile.

        Returns:
            float: The difference between the threshold and the actual change in surface 
                temperature. When this value crosses zero from above, the event is triggered.
        """
        tsurf_current = temperature[-1] * self.parameters.scalings.temperature
        tsurf_threshold = self.parameters.solver.tsurf_poststep_change * self.parameters.scalings.temperature

        if tsurf_initial[0] is None:
            tsurf_initial[0] = tsurf_current
            return 1.0  

        delta = abs(tsurf_current - tsurf_initial[0])

        return tsurf_threshold - delta  

    tsurf_event.terminal = self.parameters.solver.event_triggering
    tsurf_event.direction = -1

    return tsurf_event

`reset()`

This function initializes the model, while keeping the previous state of the PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data when running Aragog multiple times.

Source code in aragog/solver.py

def reset(self) -> None:
    """This function initializes the model, while keeping the previous state of the
    PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data
    when running Aragog multiple times.
    """
    logger.info("Resetting %s", self.__class__.__name__)
    # Update the Evaluator object except the phase properties
    self.evaluator.mesh = Mesh(self.parameters)
    self.evaluator.boundary_conditions = BoundaryConditions(self.parameters, self.evaluator.mesh)
    self.evaluator.initial_condition = InitialCondition(self.parameters, self.evaluator.mesh, self.evaluator.phases)
    # Reinstantiate the solver state
    self.state = State(self.parameters, self.evaluator)

`State(parameters, _evaluator)` `dataclass`

Stores and updates the state at temperature and pressure.

Args: parameters: Parameters evaluator: Evaluator

Attributes: evaluator: Evaluator critical_reynolds_number: Critical Reynolds number gravitational_separation_flux: Gravitational separation flux at the basic nodes heating: Total internal heat production at the staggered nodes (power per unit mass) heating_radio: Radiogenic heat production at the staggered nodes (power per unit mass) heating_tidal: Tidal heat production at the staggered nodes (power per unit mass) heat_flux: Heat flux at the basic nodes (power per unit area) inviscid_regime: True if the flow is inviscid and otherwise False, at the basic nodes inviscid_velocity: Inviscid velocity at the basic nodes is_convective: True if the flow is convecting and otherwise False, at the basic nodes reynolds_number: Reynolds number at the basic nodes super_adiabatic_temperature_gradient: Super adiabatic temperature gradient at the basic nod temperature_basic: Temperature at the basic nodes temperature_staggered: Temperature at the staggered nodes bottom_temperature: Temperature at the bottom basic node top_temperature: Temperature at the top basic node viscous_regime: True if the flow is viscous and otherwise False, at the basic nodes viscous_velocity: Viscous velocity at the basic nodes

`critical_reynolds_number` `property`

Critical Reynolds number from Abe (1993)

`heat_flux` `property`

The total heat flux according to the fluxes specified in the configuration.

`heating` `property`

The power generation according to the heat sources specified in the configuration.

`heating_dilatation` `property`

The heat source through dilation/compression.

`heating_radio` `property`

The radiogenic power generation.

`heating_tidal` `property`

The tidal power generation.

`mass_flux` `property`

The total melt mass flux according to the fluxes specified in the configuration.

`conductive_heat_flux()`

Conductive heat flux:

\[ q_{cond} = -k \frac{\partial T}{\partial r} \]

where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.

Source code in aragog/solver.py

def conductive_heat_flux(self) -> npt.NDArray:
    r"""Conductive heat flux:

    $$
    q_{cond} = -k \frac{\partial T}{\partial r}
    $$

    where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.
    """
    conductive_heat_flux: npt.NDArray = self.phase_basic.thermal_conductivity() * -self.dTdr()

    return conductive_heat_flux

`convective_heat_flux()`

Convective heat flux:

\[ q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r} - \left( \frac{\partial T}{\partial r} \right)_S \right) \]

where $\rho$ is density, $c_p$ is heat capacity at constant pressure, $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and $S$ is entropy.

Source code in aragog/solver.py

def convective_heat_flux(self) -> npt.NDArray:
    r"""Convective heat flux:

    $$
    q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r}
        - \left( \frac{\partial T}{\partial r} \right)_S \right)
    $$ 

    where $\rho$ is density, $c_p$ is heat capacity at constant pressure,
    $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and
    $S$ is entropy.
    """
    convective_heat_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.heat_capacity()
        * self.eddy_diffusivity()
        * -self._super_adiabatic_temperature_gradient
    )

    return convective_heat_flux

`dilatation_heating()`

Dilatation/compression heating (power per unit mass)

Returns: Dilatation/compression heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def dilatation_heating(self) -> npt.NDArray:
    """Dilatation/compression heating (power per unit mass)

    Returns:
        Dilatation/compression heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    mass_flux_staggered: npt.NDArray = self._evaluator.mesh.quantity_at_staggered_nodes(self._mass_flux)

    dilatation_volume_source: npt.NDArray = (
        self.phase_staggered.gravitational_acceleration()
        * self.phase_staggered.delta_specific_volume()
        * mass_flux_staggered
    )

    return dilatation_volume_source

`gravitational_separation_mass_flux()`

Gravitational separation mass flux:

\[ j_{grav} = \rho \phi (1 - \phi) v_{rel} \]

where $\rho$ is density, $\phi$ is melt fraction, and $v_{rel}$ is relative velocity.

Source code in aragog/solver.py

def gravitational_separation_mass_flux(self) -> npt.NDArray:
    r"""Gravitational separation mass flux:

    $$
    j_{grav} = \rho \phi (1 - \phi) v_{rel}
    $$ 

    where $\rho$ is density, $\phi$ is melt fraction, and
    $v_{rel}$ is relative velocity.
    """
    gravitational_separation_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.melt_fraction()
        * (1.0 - self.phase_basic.melt_fraction())
        * self.phase_basic.relative_velocity()
    )
    return gravitational_separation_mass_flux

`mixing_mass_flux()`

Mixing mass flux:

\[ j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r} \]

where $\rho$ is density, $\kappa_h$ is eddy diffusivity, $\phi$ is melt mass fraction, and $r$ is radius.

Source code in aragog/solver.py

def mixing_mass_flux(self) -> npt.NDArray:
    r"""Mixing mass flux:

    $$
    j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r}
    $$

    where $\rho$ is density, $\kappa_h$ is eddy diffusivity,
    $\phi$ is melt mass fraction, and $r$ is radius.
    """
    mixing_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.eddy_diffusivity()
        * -self.dphidr()
    )
    return mixing_mass_flux

`radiogenic_heating(time)`

Radiogenic heating (constant with radius)

Args: time: Time

Returns: Radiogenic heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def radiogenic_heating(self, time: float) -> npt.NDArray:
    """Radiogenic heating (constant with radius)

    Args:
        time: Time

    Returns:
        Radiogenic heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    # Total heat production at a given time (power per unit mass)
    radio_heating_float: float = 0
    for radionuclide in self._evaluator.radionuclides:
        radio_heating_float += radionuclide.get_heating(time)

    # Convert to 1D array (assuming abundances are constant)
    return radio_heating_float * np.ones_like(self.temperature_staggered)

`tidal_heating()`

Tidal heating at each layer of the mantle.

Args: time: Time

Returns: Tidal heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def tidal_heating(self) -> npt.NDArray:
    """Tidal heating at each layer of the mantle.

    Args:
        time: Time

    Returns:
        Tidal heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    length = len(self._settings.tidal_array)

    # Must have correct shape
    if length == 1:
        # scalar
        out = np.ones_like(self.temperature_staggered) * self._settings.tidal_array[0]

    elif length == len(self.temperature_staggered):
        # vector with correct length
        out = np.array([self._settings.tidal_array]).T

    else:
        # vector with invalid length
        raise ValueError(f"Tidal heating array has invalid length {length}")

    return out

`update(temperature, time)`

Updates the state.

The evaluation order matters because we want to minimise the number of evaluations.

Args: temperature: Temperature at the staggered nodes pressure: Pressure at the staggered nodes time: Time

Source code in aragog/solver.py

def update(self, temperature: npt.NDArray, time: FloatOrArray) -> None:
    """Updates the state.

    The evaluation order matters because we want to minimise the number of evaluations.

    Args:
        temperature: Temperature at the staggered nodes
        pressure: Pressure at the staggered nodes
        time: Time
    """
    logger.debug("Updating the state")

    logger.debug("Setting the temperature profile")
    self._temperature_staggered = temperature
    self._temperature_basic = self._evaluator.mesh.quantity_at_basic_nodes(temperature)
    logger.debug("temperature_basic = %s", self.temperature_basic)
    self._dTdr = self._evaluator.mesh.d_dr_at_basic_nodes(temperature)
    logger.debug("dTdr = %s", self.dTdr())
    self._evaluator.boundary_conditions.apply_temperature_boundary_conditions(
        temperature, self._temperature_basic, self.dTdr()
    )

    logger.debug("Setting the melt fraction profile")
    self.phase_staggered.set_temperature(temperature)
    self.phase_staggered.update()
    self.phase_basic.set_temperature(self._temperature_basic)
    self.phase_basic.update()
    phase_to_use = self._evaluator._parameters.phase_mixed.phase
    if phase_to_use == "mixed" or phase_to_use == "composite":
        self._dphidr = self._evaluator.mesh.d_dr_at_basic_nodes(
            self.phase_staggered.melt_fraction()
        )
        logger.debug("dphidr = %s", self.dphidr())
        self._evaluator.boundary_conditions.apply_temperature_boundary_conditions_melt(
            self.phase_staggered.melt_fraction(), self.phase_basic.melt_fraction(), self._dphidr
        )
    else:
        self._dphidr = np.zeros_like(self._dTdr)

    self._super_adiabatic_temperature_gradient = self.dTdr() - self.phase_basic.dTdrs()
    self._is_convective = self._super_adiabatic_temperature_gradient < 0
    velocity_prefactor: npt.NDArray = (
        self.phase_basic.gravitational_acceleration()
        * self.phase_basic.thermal_expansivity()
        * -self._super_adiabatic_temperature_gradient
    )
    # Viscous velocity
    self._viscous_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_cubed
    ) / (18 * self.phase_basic.kinematic_viscosity())
    self._viscous_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    # Inviscid velocity
    self._inviscid_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_squared
    ) / 16
    self._inviscid_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    self._inviscid_velocity[self._is_convective] = np.sqrt(
        self._inviscid_velocity[self._is_convective]
    )
    # Reynolds number
    self._reynolds_number = (
        self._viscous_velocity
        * self._evaluator.mesh.basic.mixing_length
        / self.phase_basic.kinematic_viscosity()
    )
    # Eddy diffusivity
    self._eddy_diffusivity = np.where(
        self.viscous_regime, self._viscous_velocity, self._inviscid_velocity
    )
    self._eddy_diffusivity *= self._evaluator.mesh.basic.mixing_length

    # Heat flux (power per unit area)
    self._heat_flux = np.zeros_like(self.temperature_basic)
    self._mass_flux = np.zeros_like(self.temperature_basic)
    if self._settings.conduction:
        self._heat_flux += self.conductive_heat_flux()
    if self._settings.convection:
        self._heat_flux += self.convective_heat_flux()
    if self._settings.gravitational_separation:
        self._mass_flux += self.gravitational_separation_mass_flux()
    if self._settings.mixing:
        self._mass_flux += self.mixing_mass_flux()
    self._heat_flux += self._mass_flux * self.phase_basic.latent_heat()

    # Heating (power per unit mass)
    self._heating = np.zeros_like(self.temperature_staggered)
    self._heating_radio = np.zeros_like(self.temperature_staggered)
    self._heating_dilatation = np.zeros_like(self.temperature_staggered)
    self._heating_tidal = np.zeros_like(self.temperature_staggered)

    if self._settings.radionuclides:
        self._heating_radio = self.radiogenic_heating(time)
        self._heating += self._heating_radio

    if self._settings.dilatation:
        self._heating_dilatation = self.dilatation_heating()
        self._heating += self._heating_dilatation

    if self._settings.tidal:
        self._heating_tidal = self.tidal_heating()
        self._heating += self._heating_tidal

aragog.parser

`parser`

Parses the configuration file and scales and stores the parameters.

`logger = logging.getLogger(name)` `module-attribute`

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`_BoundaryConditionsParameters(outer_boundary_condition, outer_boundary_value, inner_boundary_condition, inner_boundary_value, emissivity, equilibrium_temperature, core_heat_capacity)` `dataclass`

Stores parameters in the boundary_conditions section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.equilibrium_temperature /= self.scalings_.temperature
    self.core_heat_capacity /= self.scalings_.heat_capacity
    self._scale_inner_boundary_condition()
    self._scale_outer_boundary_condition()

`_EnergyParameters(conduction, convection, gravitational_separation, mixing, radionuclides, dilatation, tidal, tidal_array=(lambda: np.array([0.0], dtype=float))())` `dataclass`

Stores parameters in the energy section

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.tidal_array /= self.scalings_.power_per_mass

`_InitialConditionParameters(initial_condition=1, surface_temperature=4000, basal_temperature=4000, init_file='')` `dataclass`

Stores the settings in the initial_condition section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes.

Initial condition method 1: Linear profile 2: User-defined temperature field (from file) 3: Adiabatic profile

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Initial condition method
        1: Linear profile
        2: User-defined temperature field (from file)
        3: Adiabatic profile

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.surface_temperature /= self.scalings_.temperature
    self.basal_temperature /= self.scalings_.temperature

    if self.initial_condition == 2:
        if self.init_file == "":
            msg: str = (f"you must provide an initial temperature file")
            raise ValueError(msg)
        self.init_temperature = np.loadtxt(self.init_file)
        self.init_temperature /= self.scalings_.temperature

`_MeshParameters(outer_radius, inner_radius, number_of_nodes, mixing_length_profile, core_density, eos_method=1, surface_density=4000, gravitational_acceleration=9.81, adiabatic_bulk_modulus=260000000000.0, mass_coordinates=False, eos_file='')` `dataclass`

Stores parameters in the mesh section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.outer_radius /= self.scalings_.radius
    self.inner_radius /= self.scalings_.radius
    self.core_density /= self.scalings_.density
    self.surface_density /= self.scalings_.density
    self.gravitational_acceleration /= self.scalings_.gravitational_acceleration
    self.adiabatic_bulk_modulus /= self.scalings_.pressure

    if self.eos_method == 2:
        if self.eos_file == "":
            msg: str = (f"you must provide a file for setting up equation of state")
            raise ValueError(msg)
        arr = np.loadtxt(self.eos_file)
        self.eos_radius = arr[:,0] / self.scalings_.radius
        self.eos_pressure = arr[:,1] / self.scalings_.pressure
        self.eos_density = arr[:,2] / self.scalings_.density
        self.eos_gravity = arr[:,3] / self.scalings_.gravitational_acceleration
        # Check that provided eos radius roughly match with Aragog mesh
        if ((self.eos_radius[0] < self.inner_radius) or
            (self.eos_radius[-1] > self.outer_radius) or
            (self.eos_radius[-1]-self.eos_radius[0]) < 0.75*(self.outer_radius-self.inner_radius)):
            msg: str = (f"Radius array in EOS file: Values out of range.")
            raise ValueError(msg)

`_PhaseMixedParameters(latent_heat_of_fusion, rheological_transition_melt_fraction, rheological_transition_width, solidus, liquidus, phase, phase_transition_width, grain_size)` `dataclass`

Stores settings in the phase_mixed section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.latent_heat_of_fusion /= self.scalings_.latent_heat_per_mass
    self.grain_size /= self.scalings_.radius

`_PhaseParameters(density, heat_capacity, melt_fraction, thermal_conductivity, thermal_expansivity, viscosity)` `dataclass`

Stores settings in a phase section in the configuration data.

This is used to store settings from phase_liquid and phase_solid.

`scale_attributes(scalings)`

Scales the attributes if they are numbers.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes if they are numbers.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    cls_fields: tuple[Field, ...] = fields(self.__class__)
    for field_ in cls_fields:
        value: Any = getattr(self, field_.name)
        try:
            scaling: float = getattr(self.scalings_, field_.name)
            scaled_value = value / scaling
            setattr(self, field_.name, scaled_value)
            logger.info(
                "%s is a number (value = %s, scaling = %s, scaled_value = %s)",
                field_.name,
                value,
                scaling,
                scaled_value,
            )
        except AttributeError:
            logger.info("No scaling found for %s", field_.name)
        except TypeError:
            logger.info(
                "%s is a string (path to a filename) so the data will be scaled later",
                field_.name,
            )

`_Radionuclide(name, t0_years, abundance, concentration, heat_production, half_life_years)` `dataclass`

Stores the settings in a radionuclide section in the configuration data.

`get_heating(time)`

Radiogenic heating

Args: time: Time

Returns: Radiogenic heat production (power per unit mass) as a float if time is a float, otherwise a numpy row array where each entry in the row is associated with a single time in the time array.

Source code in aragog/parser.py

def get_heating(self, time: npt.NDArray | float) -> npt.NDArray | float:
    """Radiogenic heating

    Args:
        time: Time

    Returns:
        Radiogenic heat production (power per unit mass) as a float if time is a float,
            otherwise a numpy row array where each entry in the row is associated
            with a single time in the time array.
    """
    arg: npt.NDArray | float = np.log(2) * (self.t0_years - time) / self.half_life_years
    heating: npt.NDArray | float = (
        self.heat_production * self.abundance * self.concentration * np.exp(arg)
    )

    return heating

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.t0_years /= self.scalings_.time_years
    self.concentration *= 1e-6  # to mass fraction
    self.heat_production /= self.scalings_.power_per_mass
    self.half_life_years /= self.scalings_.time_years

`_ScalingsParameters(radius=1, temperature=1, density=1, time=1)` `dataclass`

Stores parameters in the scalings section in the configuration data. All units are SI.

Args: radius: Radius in metres. Defaults to 1. temperature: Temperature in Kelvin. Defaults to 1. density: Density in kg/m^3. Defaults to 1. time: Time in seconds. Defaults to 1.

Attributes: radius, m temperature, K density, kg/m^3 time, s area, m^2 kinetic_energy_per_volume, J/m^3 gravitational_acceleration, m/s^2 heat_capacity, J/kg/K heat_flux, W/m^2 latent_heat_per_mass, J/kg power_per_mass, W/kg power_per_volume, W/m^3 pressure, Pa temperature_gradient, K/m thermal_expansivity, 1/K thermal_conductivity, W/m/K velocity, m/s viscosity, Pa s time_years, years stefan_boltzmann_constant (non-dimensional)

`_SolverParameters(start_time, end_time, atol, rtol, tsurf_poststep_change=30.0, event_triggering=False)` `dataclass`

Stores settings in the solver section in the configuration data.

`_get_dataclass_from_section_name()`

Maps the section names in the configuration data to the dataclasses that stores the data.

Source code in aragog/parser.py

def _get_dataclass_from_section_name() -> dict[str, Any]:
    """Maps the section names in the configuration data to the dataclasses that stores the data."""
    mapping: dict[str, Any] = {
        "scalings": _ScalingsParameters,
        "solver": _SolverParameters,
        "boundary_conditions": _BoundaryConditionsParameters,
        "mesh": _MeshParameters,
        "energy": _EnergyParameters,
        "initial_condition": _InitialConditionParameters,
        "phase_liquid": _PhaseParameters,
        "phase_solid": _PhaseParameters,
        "phase_mixed": _PhaseMixedParameters,
        # radionuclides are dealt with separately
    }

    return mapping

aragog.phase

`phase`

A phase defines the equation of state (EOS) and transport properties.

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`CompositePhaseEvaluator(parameters)`

Bases: PhaseEvaluatorABC

Evaluates the EOS and transport properties of a composite phase.

This combines the single phase evaluators for the liquid and solid regions with the mixed phase evaluator for the mixed phase region. This ensure that the phase properties are computed correctly for all temperatures and pressures.

Args: parameters: Parameters

Source code in aragog/phase.py

def __init__(self, parameters: Parameters):
    gravitational_acceleration: PropertyProtocol = setup_gravitational_acceleration(parameters)
    self._liquid: PhaseEvaluatorProtocol = SinglePhaseEvaluator(
        parameters.phase_liquid, gravitational_acceleration
    )
    self._solid: PhaseEvaluatorProtocol = SinglePhaseEvaluator(
        parameters.phase_solid, gravitational_acceleration
    )
    self._mixed: MixedPhaseEvaluator = MixedPhaseEvaluator(parameters)

`delta_specific_volume()`

Difference of specific volume between melt and solid

Source code in aragog/phase.py

@override
def delta_specific_volume(self) -> npt.NDArray:
    """Difference of specific volume between melt and solid"""
    return self._delta_specific_volume

`heat_capacity()`

Heat capacity

Source code in aragog/phase.py

@override
def heat_capacity(self) -> npt.NDArray:
    """Heat capacity"""
    return self._heat_capacity

`latent_heat()`

Latent heat of fusion

Source code in aragog/phase.py

@override
def latent_heat(self) -> FloatOrArray:
    """Latent heat of fusion"""
    return self._mixed.latent_heat()

`melt_fraction()`

Melt fraction

Source code in aragog/phase.py

@override
def melt_fraction(self) -> npt.NDArray:
    """Melt fraction"""
    return self._mixed.melt_fraction()

`relative_velocity()`

Relative velocity between melt and solid

Source code in aragog/phase.py

@override
def relative_velocity(self) -> npt.NDArray:
    """Relative velocity between melt and solid"""
    return self._relative_velocity

`set_pressure(pressure)`

Sets pressure and updates quantities that only depend on pressure

Source code in aragog/phase.py

@override
def set_pressure(self, pressure: npt.NDArray) -> None:
    """Sets pressure and updates quantities that only depend on pressure"""
    super().set_pressure(pressure)
    self._solid.set_pressure(pressure)
    self._liquid.set_pressure(pressure)
    self._mixed.set_pressure(pressure)

`thermal_conductivity()`

Thermal conductivity

Source code in aragog/phase.py

@override
def thermal_conductivity(self) -> npt.NDArray:
    """Thermal conductivity"""
    return self._thermal_conductivity

`thermal_expansivity()`

Thermal expansivity

Source code in aragog/phase.py

@override
def thermal_expansivity(self) -> npt.NDArray:
    """Thermal expansivity"""
    return self._thermal_expansivity

`viscosity()`

Viscosity

Source code in aragog/phase.py

def viscosity(self) -> npt.NDArray:
    """Viscosity"""
    return self._viscosity

`ConstantProperty(name, *, value)` `dataclass`

Bases: PropertyProtocol

A property with a constant value

Args: name: Name of the property value: The constant value

Attributes: name: Name of the property value: The constant value ndim: Number of dimensions, which is equal to zero for a constant property

`LookupProperty1D(name, *, value)` `dataclass`

Bases: PropertyProtocol

A property from a 1-D lookup

Args: name: Name of the property value: A 2-D array, with x values in the first column and y values in the second column.

Attributes: name: Name of the property value: A 2-D array ndim: Number of dimensions, which is equal to one for a 1-D lookup

`gradient(pressure)`

Computes the gradient

Source code in aragog/phase.py

def gradient(self, pressure: npt.NDArray) -> npt.NDArray:
    """Computes the gradient"""
    return np.interp(pressure, self.value[:, 0], self._gradient)

`LookupProperty2D(name, *, value)` `dataclass`

Bases: PropertyProtocol

A property from a 2-D lookup

Args: name: Name of the property value: The 2-D array

Attributes: name: Name of the property value: The 2-D array ndim: Number of dimensions, which is equal to two for a 2-D lookup

`prepare_data_for_spline(data)`

Ensure your data is on a regular grid for RectBivariateSpline

Source code in aragog/phase.py

def prepare_data_for_spline(self, data):
    """Ensure your data is on a regular grid for RectBivariateSpline"""
    # Extract x, y, and z values
    x_values = np.unique(data[:, 0])  # Unique pressure values
    y_values = np.unique(data[:, 1])  # Unique temperature values

    # Create a grid for z values
    z_values = np.full((x_values.size, y_values.size), np.nan)

    # Find the indices of the x and y values in the unique arrays
    x_indices = np.searchsorted(x_values, data[:, 0])
    y_indices = np.searchsorted(y_values, data[:, 1])

    # Fill the z_values grid
    z_values[x_indices, y_indices] = data[:, 2]

    return x_values, y_values, z_values

`MixedPhaseEvaluator(parameters)`

Bases: PhaseEvaluatorABC

Evaluates the EOS and transport properties of a mixed phase.

This only computes quantities within the mixed phase region between the solidus and the liquidus. Computing quantities outside of this region will give incorrect results.

Args: parameters: Parameters

Attributes: settings: Mixed phase parameters

Source code in aragog/phase.py

def __init__(self, parameters: Parameters):
    gravitational_acceleration: PropertyProtocol = setup_gravitational_acceleration(parameters)
    self.settings: _PhaseMixedParameters = parameters.phase_mixed
    self._liquid: PhaseEvaluatorProtocol = SinglePhaseEvaluator(
        parameters.phase_liquid, gravitational_acceleration
    )
    self._solid: PhaseEvaluatorProtocol = SinglePhaseEvaluator(
        parameters.phase_solid, gravitational_acceleration
    )
    self._solidus: LookupProperty1D = self._get_melting_curve_lookup(
        "solidus", self.settings.solidus
    )
    self._liquidus: LookupProperty1D = self._get_melting_curve_lookup(
        "liquidus", self.settings.liquidus
    )
    self._grain_size: float = self.settings.grain_size
    self._latent_heat: float = self.settings.latent_heat_of_fusion

`dTdPs()`

TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2

Source code in aragog/interfaces.py

def dTdPs(self) -> npt.NDArray:
    """TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2"""
    dTdPs: npt.NDArray = (
        self.thermal_expansivity() * self.temperature / (self.density() * self.heat_capacity())
    )

    return dTdPs

`delta_specific_volume()`

Difference of specific volume between melt and solid

Source code in aragog/phase.py

@override
def delta_specific_volume(self) -> npt.NDArray:
    """Difference of specific volume between melt and solid"""
    return self._delta_specific_volume

`heat_capacity()`

Heat capacity of the mixed phase :cite:p:{Equation 4,}SOLO07

Source code in aragog/phase.py

@override
def heat_capacity(self) -> FloatOrArray:
    """Heat capacity of the mixed phase :cite:p:`{Equation 4,}SOLO07`"""
    return self._heat_capacity

`latent_heat()`

Latent heat of fusion

Source code in aragog/phase.py

@override
def latent_heat(self) -> float:
    """Latent heat of fusion"""
    return self._latent_heat

`liquidus()`

Liquidus

Source code in aragog/phase.py

def liquidus(self) -> npt.NDArray:
    """Liquidus"""
    liquidus: npt.NDArray = self._liquidus.eval(self.pressure)
    logger.debug("liquidus = %s", liquidus)

    return liquidus

`liquidus_gradient()`

Liquidus gradient

Source code in aragog/phase.py

def liquidus_gradient(self) -> npt.NDArray:
    """Liquidus gradient"""
    return self._liquidus.gradient(self.pressure)

`melt_fraction()`

Melt fraction of the mixed phase

The melt fraction is always between zero and one.

Source code in aragog/phase.py

@override
def melt_fraction(self) -> npt.NDArray:
    """Melt fraction of the mixed phase

    The melt fraction is always between zero and one.
    """
    return self._melt_fraction

`melt_fraction_no_clip()`

Melt fraction without clipping

Source code in aragog/phase.py

def melt_fraction_no_clip(self) -> npt.NDArray:
    """Melt fraction without clipping"""
    return self._melt_fraction_no_clip

`porosity()`

Porosity of the mixed phase, that is the volume fraction occupied by the melt

Source code in aragog/phase.py

def porosity(self) -> npt.NDArray:
    """Porosity of the mixed phase, that is the volume fraction occupied by the melt"""
    return self._porosity

`relative_velocity()`

Relative velocity between melt and solid

Source code in aragog/phase.py

@override
def relative_velocity(self) -> npt.NDArray:
    """Relative velocity between melt and solid"""
    return self._relative_velocity

`set_pressure(pressure)`

Sets pressure and updates quantities that only depend on pressure

Source code in aragog/phase.py

@override
def set_pressure(self, pressure: npt.NDArray) -> None:
    """Sets pressure and updates quantities that only depend on pressure"""
    super().set_pressure(pressure)
    # Sets the temperature of the solid and liquid phases to the appropriate melting curve in
    # order to evaluate mixed properties
    self._solid.set_temperature(self.solidus())
    self._solid.set_pressure(pressure)
    self._liquid.set_temperature(self.liquidus())
    self._liquid.set_pressure(pressure)
    self._delta_density = self._solid.density() - self._liquid.density()
    self._delta_specific_volume = 1.0/self._liquid.density() - 1.0/self._solid.density()
    self._delta_fusion = self.liquidus() - self.solidus()
    # Heat capacity of the mixed phase :cite:p:`{Equation 4,}SOLO07`
    self._heat_capacity = self.latent_heat() / self.delta_fusion()

`set_temperature(temperature)`

Sets the temperature.

Source code in aragog/interfaces.py

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Sets the temperature."""
    logger.debug("set_temperature = %s", temperature)
    self.temperature = temperature

`solidus()`

Solidus

Source code in aragog/phase.py

def solidus(self) -> npt.NDArray:
    """Solidus"""
    solidus: npt.NDArray = self._solidus.eval(self.pressure)
    logger.debug("solidus = %s", solidus)

    return solidus

`solidus_gradient()`

Solidus gradient

Source code in aragog/phase.py

def solidus_gradient(self) -> npt.NDArray:
    """Solidus gradient"""
    return self._solidus.gradient(self.pressure)

`MixedPhaseEvaluatorProtocol`

Bases: PhaseEvaluatorProtocol, Protocol

Mixed phase evaluator protocol

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`PhaseEvaluatorABC`

Bases: ABC

Phase evaluator ABC

`dTdPs()`

TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2

Source code in aragog/interfaces.py

def dTdPs(self) -> npt.NDArray:
    """TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2"""
    dTdPs: npt.NDArray = (
        self.thermal_expansivity() * self.temperature / (self.density() * self.heat_capacity())
    )

    return dTdPs

`set_pressure(pressure)`

Sets the pressure.

Source code in aragog/interfaces.py

def set_pressure(self, pressure: npt.NDArray) -> None:
    """Sets the pressure."""
    logger.debug("set_pressure = %s", pressure)
    self.pressure = pressure

`set_temperature(temperature)`

Sets the temperature.

Source code in aragog/interfaces.py

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Sets the temperature."""
    logger.debug("set_temperature = %s", temperature)
    self.temperature = temperature

`update()`

Updates quantities to avoid repeat, possibly expensive, calculations.

Source code in aragog/interfaces.py

def update(self) -> None:
    """Updates quantities to avoid repeat, possibly expensive, calculations."""

`PhaseEvaluatorCollection(parameters)` `dataclass`

A collection of phase evaluators

Creates the phase evaluators and selects the active phase based on configuration data.

Args: parameters: Parameters

Attributes: liquid: Liquid evaluator solid: Solid evaluator mixed: Mixed evaluator composite: Composite evaluator active: The active evaluator, which is defined by configuration data

`PhaseEvaluatorProtocol`

Bases: Protocol

Phase evaluator protocol

`PropertyProtocol`

Bases: Protocol

Property protocol

`SinglePhaseEvaluator(settings, gravitational_acceleration)`

Bases: PhaseEvaluatorABC

Contains the objects to evaluate the EOS and transport properties of a phase.

Args: settings: Phase parameters gravitational_acceleration: PropertyProtocol

Source code in aragog/phase.py

def __init__(self, settings: _PhaseParameters, gravitational_acceleration: PropertyProtocol):
    self._settings: _PhaseParameters = settings
    cls_fields: tuple[Field, ...] = fields(self._settings)
    for field_ in cls_fields:
        name: str = field_.name
        private_name: str = f"_{name}"
        value = getattr(self._settings, field_.name)

        if is_number(value):
            # Numbers have already been scaled by the parser
            setattr(self, private_name, ConstantProperty(name=name, value=value))

        elif is_file(value):
            with open(value, encoding="utf-8") as infile:
                logger.debug("%s is a file = %s", name, value)
                header = infile.readline()
                col_names = header[1:].split()
            value_array: npt.NDArray = np.loadtxt(value, ndmin=2)
            logger.debug("before scaling, value_array = %s", value_array)
            # Scale lookup data
            for nn, col_name in enumerate(col_names):
                logger.info("Scaling %s from %s", col_name, value)
                value_array[:, nn] /= getattr(self._settings.scalings_, col_name)
            logger.debug("after scaling, value_array = %s", value_array)
            ndim = value_array.shape[1]
            logger.debug("ndim = %d", ndim)
            if ndim == 2:
                setattr(self, private_name, LookupProperty1D(name=name, value=value_array))
            elif ndim == 3:
                setattr(self, private_name, LookupProperty2D(name=name, value=value_array))
            else:
                raise ValueError(f"Lookup data must have 2 or 3 dimensions, not {ndim}")
        else:
            logger.info("Cannot interpret value (%s): not a number or a file", value)

    self._gravitational_acceleration = gravitational_acceleration

`dTdPs()`

TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2

Source code in aragog/interfaces.py

def dTdPs(self) -> npt.NDArray:
    """TODO: Update reference to sphinx: Solomatov (2007), Treatise on Geophysics, Eq. 3.2"""
    dTdPs: npt.NDArray = (
        self.thermal_expansivity() * self.temperature / (self.density() * self.heat_capacity())
    )

    return dTdPs

`set_pressure(pressure)`

Sets the pressure.

Source code in aragog/interfaces.py

def set_pressure(self, pressure: npt.NDArray) -> None:
    """Sets the pressure."""
    logger.debug("set_pressure = %s", pressure)
    self.pressure = pressure

`set_temperature(temperature)`

Sets the temperature.

Source code in aragog/interfaces.py

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Sets the temperature."""
    logger.debug("set_temperature = %s", temperature)
    self.temperature = temperature

`update()`

Updates quantities to avoid repeat, possibly expensive, calculations.

Source code in aragog/interfaces.py

def update(self) -> None:
    """Updates quantities to avoid repeat, possibly expensive, calculations."""

`_PhaseMixedParameters(latent_heat_of_fusion, rheological_transition_melt_fraction, rheological_transition_width, solidus, liquidus, phase, phase_transition_width, grain_size)` `dataclass`

Stores settings in the phase_mixed section in the configuration data.

`scale_attributes(scalings)`

Scales the attributes

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.latent_heat_of_fusion /= self.scalings_.latent_heat_per_mass
    self.grain_size /= self.scalings_.radius

`_PhaseParameters(density, heat_capacity, melt_fraction, thermal_conductivity, thermal_expansivity, viscosity)` `dataclass`

Stores settings in a phase section in the configuration data.

This is used to store settings from phase_liquid and phase_solid.

`scale_attributes(scalings)`

Scales the attributes if they are numbers.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes if they are numbers.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    cls_fields: tuple[Field, ...] = fields(self.__class__)
    for field_ in cls_fields:
        value: Any = getattr(self, field_.name)
        try:
            scaling: float = getattr(self.scalings_, field_.name)
            scaled_value = value / scaling
            setattr(self, field_.name, scaled_value)
            logger.info(
                "%s is a number (value = %s, scaling = %s, scaled_value = %s)",
                field_.name,
                value,
                scaling,
                scaled_value,
            )
        except AttributeError:
            logger.info("No scaling found for %s", field_.name)
        except TypeError:
            logger.info(
                "%s is a string (path to a filename) so the data will be scaled later",
                field_.name,
            )

`combine_properties(weight, property1, property2)`

Linear weighting of two quantities.

Args: weight: The weight to apply to property1 property1: The value of the first property property2: The value of the second property

Returns: The combined (weighted) property

Source code in aragog/utilities.py

def combine_properties(weight: Any, property1: Any, property2: Any) -> Any:
    """Linear weighting of two quantities.

    Args:
        weight: The weight to apply to property1
        property1: The value of the first property
        property2: The value of the second property

    Returns:
        The combined (weighted) property
    """
    out = weight * property1 + (1.0 - weight) * property2

    return out

`is_file(value)`

Checks if value is a file.

Args: value: Object to be checked

Returns: True if the value is a file, otherwise False

Source code in aragog/utilities.py

def is_file(value: Any) -> bool:
    """Checks if value is a file.

    Args:
        value: Object to be checked

    Returns:
        True if the value is a file, otherwise False
    """
    if isinstance(value, (str, Path)):
        return Path(value).is_file()

    return False

`is_number(value)`

Checks if value is a number.

Args: value: Object to be checked

Returns: True if the value is a number, otherwise False

Source code in aragog/utilities.py

def is_number(value: Any) -> bool:
    """Checks if value is a number.

    Args:
        value: Object to be checked

    Returns:
        True if the value is a number, otherwise False
    """
    try:
        float(value)
        return True

    except (TypeError, ValueError):
        return False

`setup_gravitational_acceleration(parameters)`

Sets up the gravitational acceleration property.

Args: parameters: Parameters

Returns: gravitational_acceleration: PropertyProtocol

Source code in aragog/phase.py

def setup_gravitational_acceleration(parameters: Parameters):
    """Sets up the gravitational acceleration property.

    Args:
        parameters: Parameters

    Returns:
        gravitational_acceleration: PropertyProtocol
    """
    if parameters.mesh.eos_method == 1:
        # AdamsWilliamson EOS method = constant gravitational acceleration
        gravitational_acceleration = ConstantProperty(
            "gravitational_acceleration", value=parameters.mesh.gravitational_acceleration
        )
    elif parameters.mesh.eos_method == 2:
        # User defined EOS method = space dependent gravitational acceleration (1D lookup)
        # Re-interpolate EOS-data to match standard lookup pressure grid
        interp_gravity = PchipInterpolator(
            np.flip(parameters.mesh.eos_pressure),
            np.flip(parameters.mesh.eos_gravity))
        num_points: int = 138
        max_pressure: float = 1.37e11 / parameters.scalings.pressure
        lookup_data: npt.NDArray = np.zeros((num_points, 2), dtype=float)
        lookup_data[:,0] = np.linspace(0.0, max_pressure, num=num_points)
        lookup_data[:,1] = interp_gravity(lookup_data[:,0])
        gravitational_acceleration = LookupProperty1D(
            "gravitational_acceleration",
            value=lookup_data,
        )
    else:
        raise ValueError(f"EOS method = {parameters.mesh.eos_method} is not a valid selection")

    return gravitational_acceleration

`tanh_weight(value, threshold, width)`

Computes the tanh weight for viscosity profile and smoothing.

Args: value: Value threshold: Threshold value width: Width of smoothing

Returns: weight

Source code in aragog/utilities.py

def tanh_weight(value: FloatOrArray, threshold: float, width: float) -> npt.NDArray:
    """Computes the tanh weight for viscosity profile and smoothing.

    Args:
        value: Value
        threshold: Threshold value
        width: Width of smoothing

    Returns:
        weight
    """
    arg: FloatOrArray = (value - threshold) / width
    weight: npt.NDArray = 0.5 * (1.0 + np.tanh(arg))

    return weight

aragog.solver

`solver`

Solver

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`BoundaryConditions(_parameters, _mesh)` `dataclass`

Boundary conditions

Args: parameters: Parameters mesh: Mesh

`apply_flux_boundary_conditions(state)`

Applies the boundary conditions to the state.

Args: state: The state to apply the boundary conditions to

Source code in aragog/core.py

def apply_flux_boundary_conditions(self, state: State) -> None:
    """Applies the boundary conditions to the state.

    Args:
        state: The state to apply the boundary conditions to
    """
    self.apply_flux_inner_boundary_condition(state)
    self.apply_flux_outer_boundary_condition(state)
    logger.debug("temperature = %s", state.temperature_basic)
    logger.debug("heat_flux = %s", state.heat_flux)

`apply_flux_inner_boundary_condition(state)`

Applies the flux boundary condition to the state at the inner boundary.

Args: state: The state to apply the boundary conditions to

Equivalent to CORE_BC in C code. 1: Simple core cooling 2: Prescribed heat flux 3: Prescribed temperature

Source code in aragog/core.py

def apply_flux_inner_boundary_condition(self, state: State) -> None:
    """Applies the flux boundary condition to the state at the inner boundary.

    Args:
        state: The state to apply the boundary conditions to

    Equivalent to CORE_BC in C code.
        1: Simple core cooling
        2: Prescribed heat flux
        3: Prescribed temperature
    """
    if self._settings.inner_boundary_condition == 1:
        self.core_cooling(state)
    elif self._settings.inner_boundary_condition == 2:
        state.heat_flux[0, :] = self._settings.inner_boundary_value
    elif self._settings.inner_boundary_condition == 3:
        pass
        # raise NotImplementedError
    else:
        msg: str = (
            f"inner_boundary_condition = {self._settings.inner_boundary_condition} is unknown"
        )
        raise ValueError(msg)

`apply_flux_outer_boundary_condition(state)`

Applies the flux boundary condition to the state at the outer boundary.

Args: state: The state to apply the boundary conditions to

Equivalent to SURFACE_BC in C code. 1: Grey-body atmosphere 2: Zahnle steam atmosphere 3: Couple to atmodeller 4: Prescribed heat flux 5: Prescribed temperature

Source code in aragog/core.py

def apply_flux_outer_boundary_condition(self, state: State) -> None:
    """Applies the flux boundary condition to the state at the outer boundary.

    Args:
        state: The state to apply the boundary conditions to

    Equivalent to SURFACE_BC in C code.
        1: Grey-body atmosphere
        2: Zahnle steam atmosphere
        3: Couple to atmodeller
        4: Prescribed heat flux
        5: Prescribed temperature
    """
    if self._settings.outer_boundary_condition == 1:
        self.grey_body(state)
    elif self._settings.outer_boundary_condition == 2:
        raise NotImplementedError
    elif self._settings.outer_boundary_condition == 3:
        msg: str = "Requires coupling to atmodeller"
        logger.error(msg)
        raise NotImplementedError(msg)
    elif self._settings.outer_boundary_condition == 4:
        state.heat_flux[-1, :] = self._settings.outer_boundary_value
    elif self._settings.outer_boundary_condition == 5:
        pass
    else:
        msg: str = (
            f"outer_boundary_condition = {self._settings.outer_boundary_condition} is unknown"
        )
        raise ValueError(msg)

`apply_temperature_boundary_conditions(temperature, temperature_basic, dTdr)`

Conforms the temperature and dTdr at the basic nodes to temperature boundary conditions.

Args: temperature: Temperature at the staggered nodes temperature_basic: Temperature at the basic nodes dTdr: Temperature gradient at the basic nodes

Source code in aragog/core.py

def apply_temperature_boundary_conditions(
    self, temperature: npt.NDArray, temperature_basic: npt.NDArray, dTdr: npt.NDArray
) -> None:
    """Conforms the temperature and dTdr at the basic nodes to temperature boundary conditions.

    Args:
        temperature: Temperature at the staggered nodes
        temperature_basic: Temperature at the basic nodes
        dTdr: Temperature gradient at the basic nodes
    """
    # Core-mantle boundary
    if self._settings.inner_boundary_condition == 3:
        temperature_basic[0, :] = self._settings.inner_boundary_value
        dTdr[0, :] = (
            2 * (temperature[0, :] - temperature_basic[0, :])
            / self._mesh.basic.delta_mesh[0]
            * self._mesh.dxidr[0]
        )
    # Surface
    if self._settings.outer_boundary_condition == 5:
        temperature_basic[-1, :] = self._settings.outer_boundary_value
        dTdr[-1, :] = (
            2 * (temperature_basic[-1, :] - temperature[-1, :])
            / self._mesh.basic.delta_mesh[-1]
            * self._mesh.dxidr[-1]
        )

`apply_temperature_boundary_conditions_melt(melt_fraction, melt_fraction_basic, dphidr)`

Conforms the melt fraction gradient dphidr at the basic nodes to temperature boundary conditions.

Args: melt_fraction: Melt fraction at the staggered nodes melt_fraction_basic: Melt fraction at the basic nodes dphidr: Melt fraction gradient at the basic nodes

Source code in aragog/core.py

def apply_temperature_boundary_conditions_melt(
    self, melt_fraction: npt.NDArray, melt_fraction_basic: npt.NDArray, dphidr: npt.NDArray
) -> None:
    """Conforms the melt fraction gradient dphidr at the basic nodes
       to temperature boundary conditions.

    Args:
        melt_fraction: Melt fraction at the staggered nodes
        melt_fraction_basic: Melt fraction at the basic nodes
        dphidr: Melt fraction gradient at the basic nodes
    """
    # Core-mantle boundary
    if self._settings.inner_boundary_condition == 3:
        dphidr[0, :] = (
            2 * (melt_fraction[0, :] - melt_fraction_basic[0, :])
            / self._mesh.basic.delta_mesh[0]
            * self._mesh.dxidr[0]
        )
    # Surface
    if self._settings.outer_boundary_condition == 5:
        dphidr[-1, :] = (
            2 * (melt_fraction_basic[-1, :] - melt_fraction[-1, :])
            / self._mesh.basic.delta_mesh[-1]
            * self._mesh.dxidr[-1]
        )

`core_cooling(state)`

Applies a core cooling heat flux according to Eq. (37) of Bower et al., 2018

Args: state: The state to apply the boundary condition to

Source code in aragog/core.py

def core_cooling(self, state: State) -> None:
    """Applies a core cooling heat flux according to Eq. (37) of Bower et al., 2018

    Args:
        state: The state to apply the boundary condition to
    """
    core_capacity: float = (
        4
        / 3
        * np.pi
        * np.power(self._mesh.basic.radii[0], 3)
        * self._mesh.settings.core_density
        * self._settings.core_heat_capacity
    )
    cell_capacity = self._mesh.basic.volume[0] * state.capacitance_staggered()[0, :]
    radius_ratio: float = self._mesh.basic.radii[1] / self._mesh.basic.radii[0]
    alpha = np.power(radius_ratio, 2) / ((cell_capacity / (core_capacity * 1.147)) + 1)

    state.heat_flux[0, :] = alpha * state.heat_flux[1, :]

`grey_body(state)`

Applies a grey body flux at the surface.

Args: state: The state to apply the boundary conditions to

Source code in aragog/core.py

def grey_body(self, state: State) -> None:
    """Applies a grey body flux at the surface.

    Args:
        state: The state to apply the boundary conditions to
    """
    state.heat_flux[-1, :] = (
        self._settings.emissivity
        * self._settings.scalings_.stefan_boltzmann_constant
        * (np.power(state.top_temperature, 4) - self._settings.equilibrium_temperature**4)
    )

`Evaluator(_parameters)` `dataclass`

Contains classes that evaluate quantities necessary to compute the interior evolution.

Args: _parameters: Parameters

Attributes: boundary_conditions: Boundary conditions initial_condition: Initial condition mesh: Mesh phases: Evaluators for all phases radionuclides: Radionuclides

`InitialCondition(_parameters, _mesh, _phases)` `dataclass`

Initial condition

Args: parameters: Parameters mesh: Mesh phases: PhaseEvaluatorCollection

`get_adiabat(pressure_basic)`

Gets an adiabatic temperature profile by integrating the adiatiabatic temperature gradient dTdPs from the surface. Uses the set surface temperature.

Args: Pressure field on the basic nodes

Returns: Adiabatic temperature profile for the staggered nodes

Source code in aragog/core.py

def get_adiabat(self, pressure_basic) -> npt.NDArray:
    """Gets an adiabatic temperature profile by integrating
       the adiatiabatic temperature gradient dTdPs from the surface.
       Uses the set surface temperature.

    Args:
        Pressure field on the basic nodes

    Returns:
        Adiabatic temperature profile for the staggered nodes
    """

    def adiabat_ode(P,T):
        self._phases.active.set_pressure(P)
        self._phases.active.set_temperature(T)
        self._phases.active.update()
        return self._phases.active.dTdPs()

    # flip the pressure field top to bottom
    pressure_basic = np.flip(pressure_basic)

    sol = solve_ivp(
         adiabat_ode, (pressure_basic[0], pressure_basic[-1]),
         [self._settings.surface_temperature], t_eval=pressure_basic,
         method='RK45', rtol=1e-6, atol=1e-9)

    # flip back the temperature field from bottom to top
    temperature_basic = np.flip(sol.y[0])

    # Return temperature field at staggered nodes
    return self._mesh.quantity_at_staggered_nodes(temperature_basic)

`get_linear()`

Gets a linear temperature profile

Returns: Linear temperature profile for the staggered nodes Only works for uniform spatial mesh.

Source code in aragog/core.py

def get_linear(self) -> npt.NDArray:
    """Gets a linear temperature profile

    Returns:
        Linear temperature profile for the staggered nodes
        Only works for uniform spatial mesh.
    """
    temperature_basic: npt.NDArray = np.linspace(
        self._settings.basal_temperature,
        self._settings.surface_temperature,
        self._mesh.basic.number_of_nodes,
    )
    return self._mesh.quantity_at_staggered_nodes(temperature_basic)

`Mesh(parameters)`

A staggered mesh.

The basic mesh is used for the flux calculations and the staggered mesh is used for the volume calculations.

Args: parameters: Parameters

Source code in aragog/mesh.py

def __init__(self, parameters: Parameters):

    # STEP 1: Set up the basic mesh
    self.settings: _MeshParameters = parameters.mesh
    basic_coordinates: npt.NDArray = self.get_constant_spacing()
    if self.settings.eos_method == 1:
        self.eos = AdamsWilliamsonEOS(
            self.settings, basic_coordinates
        )
    elif self.settings.eos_method == 2:
        self.eos = UserDefinedEOS(
            self.settings, basic_coordinates
        )
    else:
        msg: str = (f"Unknown method to initialize Equation of State")
        raise ValueError(msg)
    if parameters.mesh.mass_coordinates:
        self._planet_density: float = self.get_planet_density(basic_coordinates)
        basic_mass_coordinates: npt.NDArray = (
            self.get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates))
        logger.debug("Basic mass coordinates = %s", basic_mass_coordinates)
    else:
        basic_mass_coordinates = basic_coordinates
    self.basic: FixedMesh = FixedMesh(
        self.settings,
        basic_coordinates,
        basic_mass_coordinates,
        np.max(basic_coordinates),
        np.min(basic_coordinates)
    )

    # STEP 2: Set up the staggered mesh
    staggered_mass_coordinates: npt.NDArray = (
        self.basic.mass_radii[:-1] + 0.5 * self.basic.delta_mesh)
    if parameters.mesh.mass_coordinates:
        staggered_coordinates: npt.NDArray = (
            self.get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates))
    else:
        staggered_coordinates = staggered_mass_coordinates
    self.staggered: FixedMesh = FixedMesh(
        self.settings,
        staggered_coordinates,
        staggered_mass_coordinates,
        self.basic.outer_boundary,
        self.basic.inner_boundary,
    )
    self.eos.set_staggered_pressure(self.staggered.radii)

    # STEP 3: Set up the transform matrices
    if parameters.mesh.mass_coordinates:
        self._dxidr: npt.NDArray = self.get_dxidr_basic()
    else:
        self._dxidr: npt.NDArray = np.ones_like(self.basic.radii)
    self._d_dr_transform: npt.NDArray = self._get_d_dr_transform_matrix()
    self._quantity_transform: npt.NDArray = self._get_quantity_transform_matrix()

`dxidr` `property`

dxi/dr at basic nodes

`d_dr_at_basic_nodes(staggered_quantity)`

Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

Args: staggered_quantity: A quantity defined at the staggered nodes.

Returns: d/dr at the basic nodes

Source code in aragog/mesh.py

def d_dr_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines d/dr at the basic nodes of a quantity defined at the staggered nodes.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes.

    Returns:
        d/dr at the basic nodes
    """
    d_dr_at_basic_nodes: npt.NDArray = self._d_dr_transform.dot(staggered_quantity)
    logger.debug("d_dr_at_basic_nodes = %s", d_dr_at_basic_nodes)

    return d_dr_at_basic_nodes

`get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates)`

Computes the basic mass coordinates from basic spatial coordinates.

Args: Basic spatial coordinates

Returns: Basic mass coordinates

Source code in aragog/mesh.py

def get_basic_mass_coordinates_from_spatial_coordinates(self, basic_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the basic mass coordinates from basic spatial coordinates.

    Args:
        Basic spatial coordinates

    Returns:
        Basic mass coordinates
    """

    # Set the mass coordinates at the inner boundary from the core mass
    basic_mass_coordinates = np.zeros_like(basic_coordinates)
    basic_mass_coordinates[:,:] = (
        self.settings.core_density / self._planet_density
        * np.power(basic_coordinates[0,:], 3.0)
    )

    # Get mass coordinates by adding individual cell contributions to the mantle mass
    basic_volumes = (np.power(basic_coordinates[1:,:],3.0)
        - np.power(basic_coordinates[:-1,:],3.0))
    for i in range(1, self.settings.number_of_nodes):
        basic_mass_coordinates[i:,:] += (
            self.staggered_effective_density[i-1,:] * basic_volumes[i-1,:] / self._planet_density
        )

    return np.power(basic_mass_coordinates, 1.0/3.0)

`get_constant_spacing()`

Constant radius spacing across the mantle

Returns: Radii with constant spacing as a column vector

Source code in aragog/mesh.py

def get_constant_spacing(self) -> npt.NDArray:
    """Constant radius spacing across the mantle

    Returns:
        Radii with constant spacing as a column vector
    """
    radii: npt.NDArray = np.linspace(
        self.settings.inner_radius, self.settings.outer_radius, self.settings.number_of_nodes
    )
    radii = np.atleast_2d(radii).T

    return radii

`get_dxidr_basic()`

Computes dxidr at basic nodes.

Source code in aragog/mesh.py

def get_dxidr_basic(self) -> npt.NDArray:
    """Computes dxidr at basic nodes."""

    dxidr = (
        self.eos.basic_density
        / self._planet_density
        * np.power(self.basic.radii,2.0)
        / np.power(self.basic.mass_radii,2.0)
    )

    return dxidr

`get_planet_density(basic_coordinates)`

Computes the planet density.

Args: Basic spatial coordinates

Returns: Planet effective density

Source code in aragog/mesh.py

def get_planet_density(self, basic_coordinates: npt.NDArray) -> float:
    """Computes the planet density.

    Args:
        Basic spatial coordinates

    Returns:
        Planet effective density
    """
    core_mass = self.settings.core_density *  np.power(basic_coordinates[0,0], 3.0)
    basic_volumes = (np.power(basic_coordinates[1:,0],3.0)
        - np.power(basic_coordinates[:-1,0],3.0))
    mantle_mass = np.sum(
        self.staggered_effective_density[:,0] * basic_volumes
    )
    planet_density = (
        (core_mass + mantle_mass) / (np.power(basic_coordinates[-1,0],3.0)))
    return planet_density.item()

`get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates)`

Computes the staggered spatial coordinates from staggered mass coordinates.

Args: Staggered mass coordinates

Returns: Staggered spatial coordinates

Source code in aragog/mesh.py

def get_staggered_spatial_coordinates_from_mass_coordinates(self, staggered_mass_coordinates: npt.NDArray) -> npt.NDArray:
    """Computes the staggered spatial coordinates from staggered mass coordinates.

    Args:
        Staggered mass coordinates

    Returns:
        Staggered spatial coordinates
    """

    # Initialise the staggered spatial coordinate to the inner boundary
    staggered_coordinates = (np.ones_like(staggered_mass_coordinates)
        * np.power(self.settings.inner_radius,3.0))

    # Add first half cell contribution
    staggered_coordinates += (
        self._planet_density
         * (np.power(staggered_mass_coordinates[0,:],3.0)
        - np.power(self.basic.mass_radii[0,:], 3.0))
        / self.staggered_effective_density[0,:]
    )

    # Get spatial coordinates by adding individual cell contributions to the mantle mass
    shell_effective_density = 0.5*(
        self.staggered_effective_density[1:,:] + self.staggered_effective_density[:-1,:])
    shell_mass_volumes = (np.power(staggered_mass_coordinates[1:,:],3.0)
        - np.power(staggered_mass_coordinates[:-1,:],3.0))
    for i in range(1,self.settings.number_of_nodes-1):
        staggered_coordinates[i:,:] += (
            self._planet_density
            * shell_mass_volumes[i-1,:]
            / shell_effective_density[i-1,:]
        )

    return np.power(staggered_coordinates, 1.0/3.0)

`quantity_at_basic_nodes(staggered_quantity)`

Determines a quantity at the basic nodes that is defined at the staggered nodes.

Uses backward and forward differences at the inner and outer radius, respectively, to obtain the quantity values of the basic nodes at the innermost and outermost nodes. When using temperature boundary conditions, values at outer boundaries will be overwritten. When using flux boundary conditions, values at outer boundaries will be used to provide estimate of individual components of heat fluxes though the total heat flux is imposed.

Args: staggered_quantity: A quantity defined at the staggered nodes

Returns: The quantity at the basic nodes

Source code in aragog/mesh.py

def quantity_at_basic_nodes(self, staggered_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the basic nodes that is defined at the staggered nodes.

    Uses backward and forward differences at the inner and outer radius, respectively, to
    obtain the quantity values of the basic nodes at the innermost and outermost nodes.
    When using temperature boundary conditions, values at outer boundaries will be overwritten.
    When using flux boundary conditions, values at outer boundaries will be used to provide
    estimate of individual components of heat fluxes though the total heat flux is imposed.

    Args:
        staggered_quantity: A quantity defined at the staggered nodes

    Returns:
        The quantity at the basic nodes
    """
    quantity_at_basic_nodes: npt.NDArray = self._quantity_transform.dot(staggered_quantity)
    logger.debug("quantity_at_basic_nodes = %s", quantity_at_basic_nodes)

    return quantity_at_basic_nodes

`quantity_at_staggered_nodes(basic_quantity)`

Determines a quantity at the staggered nodes that is defined at the basic nodes.

Staggered nodes are always located at cell centers, whatever the mesh.

Args: basic_quantity: A quantity defined at the basic nodes

Returns: The quantity at the staggered nodes

Source code in aragog/mesh.py

def quantity_at_staggered_nodes(self, basic_quantity: npt.NDArray) -> npt.NDArray:
    """Determines a quantity at the staggered nodes that is defined at the basic nodes.

    Staggered nodes are always located at cell centers, whatever the mesh.

    Args:
        basic_quantity: A quantity defined at the basic nodes

    Returns:
        The quantity at the staggered nodes
    """
    quantity_at_staggered_nodes: npt.NDArray = 0.5 * (
        basic_quantity[:-1, ...] + basic_quantity[1:, ...])
    logger.debug("quantity_at_staggered_nodes = %s", quantity_at_staggered_nodes)

    return quantity_at_staggered_nodes

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

Assembles all the parameters.

The parameters in each section are scaled here to ensure that all the parameters are scaled (non-dimensionalised) consistently with each other.

`from_file(*filenames)` `classmethod`

Parses the parameters in a configuration file(s)

Args: *filenames: Filenames of the configuration data

Source code in aragog/parser.py

@classmethod
def from_file(cls, *filenames) -> Self:
    """Parses the parameters in a configuration file(s)

    Args:
        *filenames: Filenames of the configuration data
    """
    parser: ConfigParser = ConfigParser()
    parser.read(*filenames)

    init_dict: dict[str, Any] = {}
    for section_name, dataclass_ in _get_dataclass_from_section_name().items():
        init_dict[section_name] = parser.parse_section(
            using_dataclass=dataclass_, section_name=section_name
        )
    radionuclides: list[_Radionuclide] = []
    for radionuclide_section in cls.radionuclide_sections(parser):
        radionuclide = parser.parse_section(
            using_dataclass=_Radionuclide, section_name=radionuclide_section
        )
        radionuclides.append(radionuclide)

    init_dict["radionuclides"] = radionuclides

    return cls(**init_dict)  # Unpacking gives required arguments so pylint: disable=E1125

`radionuclide_sections(parser)` `staticmethod`

Section names relating to radionuclides

Sections relating to radionuclides must have the prefix radionuclide_

Source code in aragog/parser.py

@staticmethod
def radionuclide_sections(parser: ConfigParser) -> list[str]:
    """Section names relating to radionuclides

    Sections relating to radionuclides must have the prefix radionuclide_
    """
    return [
        parser[section].name
        for section in parser.sections()
        if section.startswith("radionuclide_")
    ]

`PhaseEvaluatorCollection(parameters)` `dataclass`

A collection of phase evaluators

Creates the phase evaluators and selects the active phase based on configuration data.

Args: parameters: Parameters

Attributes: liquid: Liquid evaluator solid: Solid evaluator mixed: Mixed evaluator composite: Composite evaluator active: The active evaluator, which is defined by configuration data

`PhaseEvaluatorProtocol`

Bases: Protocol

Phase evaluator protocol

`Solver(param)`

Solves the interior dynamics

Args: filename: Filename of a file with configuration settings root: Root path to the flename

Attributes: filename: Filename of a file with configuration settings root: Root path to the filename. Defaults to empty parameters: Parameters evaluator: Evaluator state: State

Source code in aragog/solver.py

def __init__(self, param: Parameters):
    logger.info("Creating an Aragog model")
    self.parameters: Parameters = param
    self.evaluator: Evaluator
    self.state: State
    self._solution: OptimizeResult

`solution` `property`

The solution.

`temperature_basic` `property`

Temperature of the basic mesh in K

`temperature_staggered` `property`

Temperature of the staggered mesh in K

`dTdt(time, temperature)`

dT/dt at the staggered nodes

Args: time: Time temperature: Temperature at the staggered nodes

Returns: dT/dt at the staggered nodes

Source code in aragog/solver.py

def dTdt(
    self,
    time: npt.NDArray | float,
    temperature: npt.NDArray,
) -> npt.NDArray:
    """dT/dt at the staggered nodes

    Args:
        time: Time
        temperature: Temperature at the staggered nodes

    Returns:
        dT/dt at the staggered nodes
    """
    logger.debug("temperature passed into dTdt = %s", temperature)
    # logger.debug("temperature.shape = %s", temperature.shape)
    self.state.update(temperature, time)
    heat_flux: npt.NDArray = self.state.heat_flux
    # logger.debug("heat_flux = %s", heat_flux)
    self.evaluator.boundary_conditions.apply_flux_boundary_conditions(self.state)
    # logger.debug("heat_flux = %s", heat_flux)
    # logger.debug("mesh.basic.area.shape = %s", self.data.mesh.basic.area.shape)

    energy_flux: npt.NDArray = heat_flux * self.evaluator.mesh.basic.area
    # logger.debug("energy_flux size = %s", energy_flux.shape)

    delta_energy_flux: npt.NDArray = np.diff(energy_flux, axis=0)
    # logger.debug("delta_energy_flux size = %s", delta_energy_flux.shape)
    # logger.debug("capacitance = %s", self.state.phase_staggered.capacitance.shape)
    # FIXME: Update capacitance for mixed phase (enthalpy of fusion contribution)
    capacitance: npt.NDArray = (
        self.state.capacitance_staggered() * self.evaluator.mesh.basic.volume
    )

    # Heating rate (dT/dt) from flux divergence (power per unit area)
    dTdt: npt.NDArray = -delta_energy_flux / capacitance
    logger.debug("dTdt (fluxes only) = %s", dTdt)

    # Additional heating rate (dT/dt) from internal heating (power per unit mass)
    dTdt += self.state.heating * (
        self.state.phase_staggered.density() / self.state.capacitance_staggered()
    )

    logger.debug("dTdt (with internal heating) = %s", dTdt)

    return dTdt

`from_file(filename, root=Path())` `classmethod`

Parses a configuration file

Args: filename: Filename root: Root of the filename

Returns: Parameters

Source code in aragog/solver.py

@classmethod
def from_file(cls, filename: str | Path, root: str | Path = Path()) -> Self:
    """Parses a configuration file

    Args:
        filename: Filename
        root: Root of the filename

    Returns:
        Parameters
    """
    configuration_file: Path = Path(root) / Path(filename)
    logger.info("Parsing configuration file = %s", configuration_file)
    parameters: Parameters = Parameters.from_file(configuration_file)

    return cls(parameters)

`initialize()`

Initializes the model.

Source code in aragog/solver.py

def initialize(self) -> None:
    """Initializes the model."""
    logger.info("Initializing %s", self.__class__.__name__)
    self.evaluator = Evaluator(self.parameters)
    self.state = State(self.parameters, self.evaluator)

`make_tsurf_event()`

Creates a temperature event function for use with an ODE solver to monitor changes in the surface temperature.The event triggers when the change exceeds the threshold, allowing the solver to stop integration.

Returns: The event has the attributes: - terminal = True: Integration stops when the event is triggered. - direction = -1: Only triggers when the function is decreasing through zero.

Source code in aragog/solver.py

def make_tsurf_event(self):
    """
    Creates a temperature event function for use with an ODE solver to monitor changes 
    in the surface temperature.The event triggers when the change exceeds the 
    threshold, allowing the solver to stop integration.

    Returns:
        The event has the attributes:
            - terminal = True: Integration stops when the event is triggered.
            - direction = -1: Only triggers when the function is decreasing through zero.
    """
    tsurf_initial = [None]

    def tsurf_event(time: float, temperature: npt.NDArray) -> float:
        """
        Event function to detect when surface temperature changes beyond a specified threshold.

        Args:
            time (float): Current time.
            temperature (np.ndarray): Current temperature profile.

        Returns:
            float: The difference between the threshold and the actual change in surface 
                temperature. When this value crosses zero from above, the event is triggered.
        """
        tsurf_current = temperature[-1] * self.parameters.scalings.temperature
        tsurf_threshold = self.parameters.solver.tsurf_poststep_change * self.parameters.scalings.temperature

        if tsurf_initial[0] is None:
            tsurf_initial[0] = tsurf_current
            return 1.0  

        delta = abs(tsurf_current - tsurf_initial[0])

        return tsurf_threshold - delta  

    tsurf_event.terminal = self.parameters.solver.event_triggering
    tsurf_event.direction = -1

    return tsurf_event

`reset()`

This function initializes the model, while keeping the previous state of the PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data when running Aragog multiple times.

Source code in aragog/solver.py

def reset(self) -> None:
    """This function initializes the model, while keeping the previous state of the
    PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data
    when running Aragog multiple times.
    """
    logger.info("Resetting %s", self.__class__.__name__)
    # Update the Evaluator object except the phase properties
    self.evaluator.mesh = Mesh(self.parameters)
    self.evaluator.boundary_conditions = BoundaryConditions(self.parameters, self.evaluator.mesh)
    self.evaluator.initial_condition = InitialCondition(self.parameters, self.evaluator.mesh, self.evaluator.phases)
    # Reinstantiate the solver state
    self.state = State(self.parameters, self.evaluator)

`State(parameters, _evaluator)` `dataclass`

Stores and updates the state at temperature and pressure.

Args: parameters: Parameters evaluator: Evaluator

Attributes: evaluator: Evaluator critical_reynolds_number: Critical Reynolds number gravitational_separation_flux: Gravitational separation flux at the basic nodes heating: Total internal heat production at the staggered nodes (power per unit mass) heating_radio: Radiogenic heat production at the staggered nodes (power per unit mass) heating_tidal: Tidal heat production at the staggered nodes (power per unit mass) heat_flux: Heat flux at the basic nodes (power per unit area) inviscid_regime: True if the flow is inviscid and otherwise False, at the basic nodes inviscid_velocity: Inviscid velocity at the basic nodes is_convective: True if the flow is convecting and otherwise False, at the basic nodes reynolds_number: Reynolds number at the basic nodes super_adiabatic_temperature_gradient: Super adiabatic temperature gradient at the basic nod temperature_basic: Temperature at the basic nodes temperature_staggered: Temperature at the staggered nodes bottom_temperature: Temperature at the bottom basic node top_temperature: Temperature at the top basic node viscous_regime: True if the flow is viscous and otherwise False, at the basic nodes viscous_velocity: Viscous velocity at the basic nodes

`critical_reynolds_number` `property`

Critical Reynolds number from Abe (1993)

`heat_flux` `property`

The total heat flux according to the fluxes specified in the configuration.

`heating` `property`

The power generation according to the heat sources specified in the configuration.

`heating_dilatation` `property`

The heat source through dilation/compression.

`heating_radio` `property`

The radiogenic power generation.

`heating_tidal` `property`

The tidal power generation.

`mass_flux` `property`

The total melt mass flux according to the fluxes specified in the configuration.

`conductive_heat_flux()`

Conductive heat flux:

\[ q_{cond} = -k \frac{\partial T}{\partial r} \]

where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.

Source code in aragog/solver.py

def conductive_heat_flux(self) -> npt.NDArray:
    r"""Conductive heat flux:

    $$
    q_{cond} = -k \frac{\partial T}{\partial r}
    $$

    where $k$ is thermal conductivity, $T$ is temperature, and $r$ is radius.
    """
    conductive_heat_flux: npt.NDArray = self.phase_basic.thermal_conductivity() * -self.dTdr()

    return conductive_heat_flux

`convective_heat_flux()`

Convective heat flux:

\[ q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r} - \left( \frac{\partial T}{\partial r} \right)_S \right) \]

where $\rho$ is density, $c_p$ is heat capacity at constant pressure, $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and $S$ is entropy.

Source code in aragog/solver.py

def convective_heat_flux(self) -> npt.NDArray:
    r"""Convective heat flux:

    $$
    q_{conv} = -\rho c_p \kappa_h \left( \frac{\partial T}{\partial r}
        - \left( \frac{\partial T}{\partial r} \right)_S \right)
    $$ 

    where $\rho$ is density, $c_p$ is heat capacity at constant pressure,
    $\kappa_h$ is eddy diffusivity, $T$ is temperature, $r$ is radius, and
    $S$ is entropy.
    """
    convective_heat_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.heat_capacity()
        * self.eddy_diffusivity()
        * -self._super_adiabatic_temperature_gradient
    )

    return convective_heat_flux

`dilatation_heating()`

Dilatation/compression heating (power per unit mass)

Returns: Dilatation/compression heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def dilatation_heating(self) -> npt.NDArray:
    """Dilatation/compression heating (power per unit mass)

    Returns:
        Dilatation/compression heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    mass_flux_staggered: npt.NDArray = self._evaluator.mesh.quantity_at_staggered_nodes(self._mass_flux)

    dilatation_volume_source: npt.NDArray = (
        self.phase_staggered.gravitational_acceleration()
        * self.phase_staggered.delta_specific_volume()
        * mass_flux_staggered
    )

    return dilatation_volume_source

`gravitational_separation_mass_flux()`

Gravitational separation mass flux:

\[ j_{grav} = \rho \phi (1 - \phi) v_{rel} \]

where $\rho$ is density, $\phi$ is melt fraction, and $v_{rel}$ is relative velocity.

Source code in aragog/solver.py

def gravitational_separation_mass_flux(self) -> npt.NDArray:
    r"""Gravitational separation mass flux:

    $$
    j_{grav} = \rho \phi (1 - \phi) v_{rel}
    $$ 

    where $\rho$ is density, $\phi$ is melt fraction, and
    $v_{rel}$ is relative velocity.
    """
    gravitational_separation_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.phase_basic.melt_fraction()
        * (1.0 - self.phase_basic.melt_fraction())
        * self.phase_basic.relative_velocity()
    )
    return gravitational_separation_mass_flux

`mixing_mass_flux()`

Mixing mass flux:

\[ j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r} \]

where $\rho$ is density, $\kappa_h$ is eddy diffusivity, $\phi$ is melt mass fraction, and $r$ is radius.

Source code in aragog/solver.py

def mixing_mass_flux(self) -> npt.NDArray:
    r"""Mixing mass flux:

    $$
    j_{cm} = -\rho \kappa_h \frac{\partial \phi}{\partial r}
    $$

    where $\rho$ is density, $\kappa_h$ is eddy diffusivity,
    $\phi$ is melt mass fraction, and $r$ is radius.
    """
    mixing_mass_flux: npt.NDArray = (
        self.phase_basic.density()
        * self.eddy_diffusivity()
        * -self.dphidr()
    )
    return mixing_mass_flux

`radiogenic_heating(time)`

Radiogenic heating (constant with radius)

Args: time: Time

Returns: Radiogenic heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def radiogenic_heating(self, time: float) -> npt.NDArray:
    """Radiogenic heating (constant with radius)

    Args:
        time: Time

    Returns:
        Radiogenic heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    # Total heat production at a given time (power per unit mass)
    radio_heating_float: float = 0
    for radionuclide in self._evaluator.radionuclides:
        radio_heating_float += radionuclide.get_heating(time)

    # Convert to 1D array (assuming abundances are constant)
    return radio_heating_float * np.ones_like(self.temperature_staggered)

`tidal_heating()`

Tidal heating at each layer of the mantle.

Args: time: Time

Returns: Tidal heating (power per unit mass) at each layer of the staggered mesh, at a given point in time.

Source code in aragog/solver.py

def tidal_heating(self) -> npt.NDArray:
    """Tidal heating at each layer of the mantle.

    Args:
        time: Time

    Returns:
        Tidal heating (power per unit mass) at each layer of the staggered
            mesh, at a given point in time.
    """

    length = len(self._settings.tidal_array)

    # Must have correct shape
    if length == 1:
        # scalar
        out = np.ones_like(self.temperature_staggered) * self._settings.tidal_array[0]

    elif length == len(self.temperature_staggered):
        # vector with correct length
        out = np.array([self._settings.tidal_array]).T

    else:
        # vector with invalid length
        raise ValueError(f"Tidal heating array has invalid length {length}")

    return out

`update(temperature, time)`

Updates the state.

The evaluation order matters because we want to minimise the number of evaluations.

Args: temperature: Temperature at the staggered nodes pressure: Pressure at the staggered nodes time: Time

Source code in aragog/solver.py

def update(self, temperature: npt.NDArray, time: FloatOrArray) -> None:
    """Updates the state.

    The evaluation order matters because we want to minimise the number of evaluations.

    Args:
        temperature: Temperature at the staggered nodes
        pressure: Pressure at the staggered nodes
        time: Time
    """
    logger.debug("Updating the state")

    logger.debug("Setting the temperature profile")
    self._temperature_staggered = temperature
    self._temperature_basic = self._evaluator.mesh.quantity_at_basic_nodes(temperature)
    logger.debug("temperature_basic = %s", self.temperature_basic)
    self._dTdr = self._evaluator.mesh.d_dr_at_basic_nodes(temperature)
    logger.debug("dTdr = %s", self.dTdr())
    self._evaluator.boundary_conditions.apply_temperature_boundary_conditions(
        temperature, self._temperature_basic, self.dTdr()
    )

    logger.debug("Setting the melt fraction profile")
    self.phase_staggered.set_temperature(temperature)
    self.phase_staggered.update()
    self.phase_basic.set_temperature(self._temperature_basic)
    self.phase_basic.update()
    phase_to_use = self._evaluator._parameters.phase_mixed.phase
    if phase_to_use == "mixed" or phase_to_use == "composite":
        self._dphidr = self._evaluator.mesh.d_dr_at_basic_nodes(
            self.phase_staggered.melt_fraction()
        )
        logger.debug("dphidr = %s", self.dphidr())
        self._evaluator.boundary_conditions.apply_temperature_boundary_conditions_melt(
            self.phase_staggered.melt_fraction(), self.phase_basic.melt_fraction(), self._dphidr
        )
    else:
        self._dphidr = np.zeros_like(self._dTdr)

    self._super_adiabatic_temperature_gradient = self.dTdr() - self.phase_basic.dTdrs()
    self._is_convective = self._super_adiabatic_temperature_gradient < 0
    velocity_prefactor: npt.NDArray = (
        self.phase_basic.gravitational_acceleration()
        * self.phase_basic.thermal_expansivity()
        * -self._super_adiabatic_temperature_gradient
    )
    # Viscous velocity
    self._viscous_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_cubed
    ) / (18 * self.phase_basic.kinematic_viscosity())
    self._viscous_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    # Inviscid velocity
    self._inviscid_velocity = (
        velocity_prefactor * self._evaluator.mesh.basic.mixing_length_squared
    ) / 16
    self._inviscid_velocity[~self.is_convective] = 0  # Must be super-adiabatic
    self._inviscid_velocity[self._is_convective] = np.sqrt(
        self._inviscid_velocity[self._is_convective]
    )
    # Reynolds number
    self._reynolds_number = (
        self._viscous_velocity
        * self._evaluator.mesh.basic.mixing_length
        / self.phase_basic.kinematic_viscosity()
    )
    # Eddy diffusivity
    self._eddy_diffusivity = np.where(
        self.viscous_regime, self._viscous_velocity, self._inviscid_velocity
    )
    self._eddy_diffusivity *= self._evaluator.mesh.basic.mixing_length

    # Heat flux (power per unit area)
    self._heat_flux = np.zeros_like(self.temperature_basic)
    self._mass_flux = np.zeros_like(self.temperature_basic)
    if self._settings.conduction:
        self._heat_flux += self.conductive_heat_flux()
    if self._settings.convection:
        self._heat_flux += self.convective_heat_flux()
    if self._settings.gravitational_separation:
        self._mass_flux += self.gravitational_separation_mass_flux()
    if self._settings.mixing:
        self._mass_flux += self.mixing_mass_flux()
    self._heat_flux += self._mass_flux * self.phase_basic.latent_heat()

    # Heating (power per unit mass)
    self._heating = np.zeros_like(self.temperature_staggered)
    self._heating_radio = np.zeros_like(self.temperature_staggered)
    self._heating_dilatation = np.zeros_like(self.temperature_staggered)
    self._heating_tidal = np.zeros_like(self.temperature_staggered)

    if self._settings.radionuclides:
        self._heating_radio = self.radiogenic_heating(time)
        self._heating += self._heating_radio

    if self._settings.dilatation:
        self._heating_dilatation = self.dilatation_heating()
        self._heating += self._heating_dilatation

    if self._settings.tidal:
        self._heating_tidal = self.tidal_heating()
        self._heating += self._heating_tidal

`_EnergyParameters(conduction, convection, gravitational_separation, mixing, radionuclides, dilatation, tidal, tidal_array=(lambda: np.array([0.0], dtype=float))())` `dataclass`

Stores parameters in the energy section

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.tidal_array /= self.scalings_.power_per_mass

`_Radionuclide(name, t0_years, abundance, concentration, heat_production, half_life_years)` `dataclass`

Stores the settings in a radionuclide section in the configuration data.

`get_heating(time)`

Radiogenic heating

Args: time: Time

Returns: Radiogenic heat production (power per unit mass) as a float if time is a float, otherwise a numpy row array where each entry in the row is associated with a single time in the time array.

Source code in aragog/parser.py

def get_heating(self, time: npt.NDArray | float) -> npt.NDArray | float:
    """Radiogenic heating

    Args:
        time: Time

    Returns:
        Radiogenic heat production (power per unit mass) as a float if time is a float,
            otherwise a numpy row array where each entry in the row is associated
            with a single time in the time array.
    """
    arg: npt.NDArray | float = np.log(2) * (self.t0_years - time) / self.half_life_years
    heating: npt.NDArray | float = (
        self.heat_production * self.abundance * self.concentration * np.exp(arg)
    )

    return heating

`scale_attributes(scalings)`

Scales the attributes.

Args: scalings: scalings

Source code in aragog/parser.py

def scale_attributes(self, scalings: _ScalingsParameters) -> None:
    """Scales the attributes.

    Args:
        scalings: scalings
    """
    self.scalings_ = scalings
    self.t0_years /= self.scalings_.time_years
    self.concentration *= 1e-6  # to mass fraction
    self.heat_production /= self.scalings_.power_per_mass
    self.half_life_years /= self.scalings_.time_years

aragog.utilities

`utilities`

Utilities

`FloatOrArray = float | npt.NDArray` `module-attribute`

`MultiplyT = TypeVar('MultiplyT', float, npt.NDArray, pd.Series, pd.DataFrame)` `module-attribute`

`combine_properties(weight, property1, property2)`

Linear weighting of two quantities.

Args: weight: The weight to apply to property1 property1: The value of the first property property2: The value of the second property

Returns: The combined (weighted) property

Source code in aragog/utilities.py

def combine_properties(weight: Any, property1: Any, property2: Any) -> Any:
    """Linear weighting of two quantities.

    Args:
        weight: The weight to apply to property1
        property1: The value of the first property
        property2: The value of the second property

    Returns:
        The combined (weighted) property
    """
    out = weight * property1 + (1.0 - weight) * property2

    return out

`is_file(value)`

Checks if value is a file.

Args: value: Object to be checked

Returns: True if the value is a file, otherwise False

Source code in aragog/utilities.py

def is_file(value: Any) -> bool:
    """Checks if value is a file.

    Args:
        value: Object to be checked

    Returns:
        True if the value is a file, otherwise False
    """
    if isinstance(value, (str, Path)):
        return Path(value).is_file()

    return False

`is_monotonic_increasing(some_array)`

Returns True if an array is monotonically increasing, otherwise returns False.

Source code in aragog/utilities.py

def is_monotonic_increasing(some_array: npt.NDArray) -> bool:
    """Returns True if an array is monotonically increasing, otherwise returns False."""
    return np.all(np.diff(some_array) > 0)

`is_number(value)`

Checks if value is a number.

Args: value: Object to be checked

Returns: True if the value is a number, otherwise False

Source code in aragog/utilities.py

def is_number(value: Any) -> bool:
    """Checks if value is a number.

    Args:
        value: Object to be checked

    Returns:
        True if the value is a number, otherwise False
    """
    try:
        float(value)
        return True

    except (TypeError, ValueError):
        return False

`profile_decorator(func)`

Decorator to profile a function

Source code in aragog/utilities.py

def profile_decorator(func):
    """Decorator to profile a function"""

    @wraps(func)
    def wrapper(*args, **kwargs):
        with Profile() as profile:
            result = func(*args, **kwargs)
        stats = Stats(profile).strip_dirs().sort_stats(SortKey.TIME)
        stats.print_stats()
        return result

    return wrapper

`tanh_weight(value, threshold, width)`

Computes the tanh weight for viscosity profile and smoothing.

Args: value: Value threshold: Threshold value width: Width of smoothing

Returns: weight

Source code in aragog/utilities.py

def tanh_weight(value: FloatOrArray, threshold: float, width: float) -> npt.NDArray:
    """Computes the tanh weight for viscosity profile and smoothing.

    Args:
        value: Value
        threshold: Threshold value
        width: Width of smoothing

    Returns:
        weight
    """
    arg: FloatOrArray = (value - threshold) / width
    weight: npt.NDArray = 0.5 * (1.0 + np.tanh(arg))

    return weight

aragog (package)

`aragog`

Package level variables and initialises the package logger

`CFG_DATA = importlib.resources.files(f'{package}.cfg')` `module-attribute`

`version = '26.01.06'` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`Output(solver)`

Stores inputs and outputs of the models.

Source code in aragog/output.py

def __init__(self, solver: Solver):
    self.solver: Solver = solver
    self.parameters: Parameters = solver.parameters
    self.solution: OptimizeResult = self.solver.solution
    self.evaluator: Evaluator = self.solver.evaluator
    self.state: State = self.solver.state

`conductive_heat_flux_basic` `property`

Conductive heat flux

`convective_heat_flux_basic` `property`

Convective heat flux

`core_mass` `property`

Core mass computed with constant density

`dTdr` `property`

dTdr

`dTdrs` `property`

dTdrs

`density_basic` `property`

Density

`gravitational_separation_heat_flux_basic` `property`

Gravitational separation heat flux

`heat_capacity_basic` `property`

Heat capacity

`heating` `property`

Internal heat generation at staggered nodes

`heating_dilatation` `property`

Internal heat generation from dilatation/compression at staggered nodes

`heating_radio` `property`

Internal heat generation from radioactive decay at staggered nodes

`heating_tidal` `property`

Internal heat generation from tidal heat dissipation at staggered nodes

`liquidus_K_staggered` `property`

Liquidus

`log10_viscosity_basic` `property`

Viscosity of the basic mesh

`log10_viscosity_staggered` `property`

Viscosity of the staggered mesh

`mantle_mass` `property`

Mantle mass computed from the AdamsWilliamsonEOS

`mass_radii_km_basic` `property`

Mass radii of the basic mesh in km

`mass_radii_km_staggered` `property`

Mass radii of the staggered mesh in km

`mass_staggered` `property`

Mass of each layer on staggered mesh

`melt_fraction_basic` `property`

Melt fraction on the basic mesh

`melt_fraction_global` `property`

Volume-averaged melt fraction

`melt_fraction_staggered` `property`

Melt fraction on the staggered mesh

`mixing_heat_flux_basic` `property`

Convective mixing heat flux

`pressure_GPa_basic` `property`

Pressure of the basic mesh in GPa

`pressure_GPa_staggered` `property`

Pressure of the staggered mesh in GPa

`radii_km_basic` `property`

Radii of the basic mesh in km

`radii_km_staggered` `property`

Radii of the staggered mesh in km

`rheological_front` `property`

Rheological front at the last solve iteration given user defined threshold. It is defined as a dimensionless distance with respect to the outer radius.

`shape_basic` `property`

Shape of the basic data

`shape_staggered` `property`

Shape of the staggered data

`solidus_K_staggered` `property`

Solidus

`solution_top_temperature` `property`

Solution (last iteration) temperature at the top of the domain (planet surface)

`super_adiabatic_temperature_gradient_basic` `property`

Super adiabatic temperature gradient

`temperature_K_basic` `property`

Temperature of the basic mesh in K

`temperature_K_staggered` `property`

Temperature of the staggered mesh in K

`thermal_expansivity_basic` `property`

Thermal expansivity

`times` `property`

Times in years

`total_heat_flux_basic` `property`

Conductive heat flux

`plot(num_lines=11, figsize=(25, 10))`

Plots the solution with labelled lines according to time.

Args: num_lines: Number of lines to plot. Defaults to 11. figsize: Size of the figure. Defaults to (25, 10).

Source code in aragog/output.py

def plot(self, num_lines: int = 11, figsize: tuple = (25, 10)) -> None:
    """Plots the solution with labelled lines according to time.

    Args:
        num_lines: Number of lines to plot. Defaults to 11.
        figsize: Size of the figure. Defaults to (25, 10).
    """
    assert self.solution is not None

    self.state.update(self.solution.y, self.solution.t)

    _, axs = plt.subplots(1, 10, sharey=True, figsize=figsize)

    # Ensure there are at least 2 lines to plot (first and last).
    num_lines = max(2, num_lines)

    # Calculate the time step based on the total number of lines.
    time_step: float = self.time_range / (num_lines - 1)

    # plot temperature
    try:
        axs[0].scatter(self.liquidus_K_staggered, self.pressure_GPa_staggered)
        axs[0].scatter(self.solidus_K_staggered, self.pressure_GPa_staggered)
    except AttributeError:
        pass

    # Plot the first line.
    def plot_times(ax, x: npt.NDArray, y: npt.NDArray) -> None:
        # If `x` is a float, create an array of the same length as `y` filled with `x`
        if np.isscalar(x):
            x = np.full((len(y), len(self.times)), x, dtype=np.float64)
        elif x.ndim == 1:
            # If `x` is 1D, reshape it or repeat it across the time dimension
            x = np.tile(x[:, np.newaxis], (1, len(self.times)))

        label_first: str = f"{self.times[0]:.2f}"
        ax.plot(x[:, 0], y, label=label_first)

        # Loop through the selected lines and plot each with a label.
        for i in range(1, num_lines - 1):
            desired_time: float = self.times[0] + i * time_step
            # Find the closest available time step.
            closest_time_index: np.intp = np.argmin(np.abs(self.times - desired_time))
            time: float = self.times[closest_time_index]
            label: str = f"{time:.2f}"  # Create a label based on the time.
            ax.plot(
                x[:, closest_time_index],
                y,
                label=label,
            )

        # Plot the last line.
        times_end: float = self.times[-1]
        label_last: str = f"{times_end:.2f}"
        ax.plot(x[:, -1], y, label=label_last)

    plot_times(axs[0], self.temperature_K_basic, self.pressure_GPa_basic)
    axs[0].set_ylabel("Pressure (GPa)")
    axs[0].set_xlabel("Temperature (K)")
    axs[0].set_title("Temperature")

    plot_times(axs[1], self.melt_fraction_staggered, self.pressure_GPa_staggered)
    axs[1].set_xlabel("Melt fraction")
    axs[1].set_title("Melt fraction")

    plot_times(axs[2], self.log10_viscosity_basic, self.pressure_GPa_basic)
    axs[2].set_xlabel("Log10(viscosity)")
    axs[2].set_title("Log10(viscosity)")

    plot_times(axs[3], self.convective_heat_flux_basic, self.pressure_GPa_basic)
    axs[3].set_xlabel("Convective heat flux")
    axs[3].set_title("Convective heat flux")

    plot_times(
        axs[4],
        self.super_adiabatic_temperature_gradient_basic,
        self.pressure_GPa_basic,
    )
    axs[4].set_xlabel("Super adiabatic temperature gradient")
    axs[4].set_title("Super adiabatic temperature gradient")

    plot_times(
        axs[5],
        self.dTdr,
        self.pressure_GPa_basic,
    )
    axs[5].set_xlabel("dTdr")
    axs[5].set_title("dTdr")

    plot_times(
        axs[6],
        self.dTdrs,
        self.pressure_GPa_basic,
    )
    axs[6].set_xlabel("dTdrs")
    axs[6].set_title("dTdrs")

    plot_times(
        axs[7],
        self.density_basic,
        self.pressure_GPa_basic,
    )
    axs[7].set_xlabel("Density")
    axs[7].set_title("Density")

    plot_times(
        axs[8],
        self.heat_capacity_basic,
        self.pressure_GPa_basic,
    )
    axs[8].set_xlabel("Heat capacity")
    axs[8].set_title("Heat capacity")

    plot_times(
        axs[9],
        self.thermal_expansivity_basic,
        self.pressure_GPa_basic,
    )
    axs[9].set_xlabel("Thermal expansivity")
    axs[9].set_title("Thermal expansivity")

    # Shrink current axis by 20% to allow space for the legend.
    # box = ax.get_position()
    # ax.set_position((box.x0, box.y0, box.width * 0.8, box.height))

    # axs[0].grid(True)

    legend = axs[2].legend()  # (loc="center left", bbox_to_anchor=(1, 0.5))
    legend.set_title("Time (yr)")
    plt.gca().invert_yaxis()

    plt.show()

`write_at_time(file_path, tidx=-1, compress=False)`

Write the state of the model at a particular time to a NetCDF4 file on the disk.

Args: file_path: Path to the output file tidx: Index on the time axis at which to access the data compress: Whether to compress the data

Source code in aragog/output.py

def write_at_time(self, file_path: str, tidx: int = -1, compress: bool = False) -> None:
    """Write the state of the model at a particular time to a NetCDF4 file on the disk.

    Args:
        file_path: Path to the output file
        tidx: Index on the time axis at which to access the data
        compress: Whether to compress the data
    """

    logger.debug("Writing i=%d NetCDF file to %s", tidx, file_path)

    # Update the state
    assert self.solution is not None
    self.state.update(self.solution.y, self.solution.t)

    # Open the dataset
    ds: nc.Dataset = nc.Dataset(file_path, mode="w")

    # Metadata
    ds.description = "Aragog output data"
    ds.argog_version = __version__

    # Function to save scalar quantities
    def _add_scalar_variable(key: str, value: float, units: str):
        ds.createVariable(key, np.float64)
        ds[key][0] = value.item() if hasattr(value, 'item') else value
        ds[key].units = units

    # Save scalar quantities
    _add_scalar_variable("time", self.times[tidx], "yr")
    _add_scalar_variable("phi_global", self.melt_fraction_global, "")
    _add_scalar_variable("mantle_mass", self.mantle_mass, "kg")
    _add_scalar_variable("rheo_front", self.rheological_front, "")

    # Create dimensions for mesh quantities
    ds.createDimension("basic", self.shape_basic[0])
    ds.createDimension("staggered", self.shape_staggered[0])

    # Function to save mesh quantities
    def _add_mesh_variable(key: str, some_property: Any, units: str):
        if key[-2:] == "_b":
            mesh = "basic"
        elif key[-2:] == "_s":
            mesh = "staggered"
        else:
            raise KeyError(f"NetCDF variable name must end in _b or _s: {key}")
        ds.createVariable(
            key,
            np.float64,
            (mesh,),
            compression="zlib" if compress else None,
            shuffle = True,
            complevel = 4
        )
        ds[key][:] = some_property[:, tidx]
        ds[key].units = units

    # Save mesh quantities
    _add_mesh_variable("radius_b", self.radii_km_basic, "km")
    _add_mesh_variable("pres_b", self.pressure_GPa_basic, "GPa")
    _add_mesh_variable("temp_b", self.temperature_K_basic, "K")
    _add_mesh_variable("phi_b", self.melt_fraction_basic, "")
    _add_mesh_variable("Fcond_b", self.conductive_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fconv_b", self.convective_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fgrav_b", self.gravitational_separation_heat_flux_basic, "W m-2")
    _add_mesh_variable("Fmix_b", self.mixing_heat_flux_basic, "W m-2")
    _add_mesh_variable("Ftotal_b", self.total_heat_flux_basic, "W m-2")
    _add_mesh_variable("log10visc_b", self.log10_viscosity_basic, "Pa s")
    _add_mesh_variable("log10visc_s", self.log10_viscosity_staggered, "Pa s")
    _add_mesh_variable("density_b", self.density_basic, "kg m-3")
    _add_mesh_variable("heatcap_b", self.heat_capacity_basic, "J kg-1 K-1")
    _add_mesh_variable("mass_s", self.mass_staggered, "kg")
    _add_mesh_variable("Hradio_s", self.heating_radio, "W kg-1")
    _add_mesh_variable("Hvol_s", self.heating_dilatation, "W kg-1")
    _add_mesh_variable("Htidal_s", self.heating_tidal, "W kg-1")
    _add_mesh_variable("Htotal_s", self.heating, "W kg-1")

    # Close the dataset
    ds.close()

`Solver(param)`

Solves the interior dynamics

Args: filename: Filename of a file with configuration settings root: Root path to the flename

Attributes: filename: Filename of a file with configuration settings root: Root path to the filename. Defaults to empty parameters: Parameters evaluator: Evaluator state: State

Source code in aragog/solver.py

def __init__(self, param: Parameters):
    logger.info("Creating an Aragog model")
    self.parameters: Parameters = param
    self.evaluator: Evaluator
    self.state: State
    self._solution: OptimizeResult

`solution` `property`

The solution.

`temperature_basic` `property`

Temperature of the basic mesh in K

`temperature_staggered` `property`

Temperature of the staggered mesh in K

`dTdt(time, temperature)`

dT/dt at the staggered nodes

Args: time: Time temperature: Temperature at the staggered nodes

Returns: dT/dt at the staggered nodes

Source code in aragog/solver.py

def dTdt(
    self,
    time: npt.NDArray | float,
    temperature: npt.NDArray,
) -> npt.NDArray:
    """dT/dt at the staggered nodes

    Args:
        time: Time
        temperature: Temperature at the staggered nodes

    Returns:
        dT/dt at the staggered nodes
    """
    logger.debug("temperature passed into dTdt = %s", temperature)
    # logger.debug("temperature.shape = %s", temperature.shape)
    self.state.update(temperature, time)
    heat_flux: npt.NDArray = self.state.heat_flux
    # logger.debug("heat_flux = %s", heat_flux)
    self.evaluator.boundary_conditions.apply_flux_boundary_conditions(self.state)
    # logger.debug("heat_flux = %s", heat_flux)
    # logger.debug("mesh.basic.area.shape = %s", self.data.mesh.basic.area.shape)

    energy_flux: npt.NDArray = heat_flux * self.evaluator.mesh.basic.area
    # logger.debug("energy_flux size = %s", energy_flux.shape)

    delta_energy_flux: npt.NDArray = np.diff(energy_flux, axis=0)
    # logger.debug("delta_energy_flux size = %s", delta_energy_flux.shape)
    # logger.debug("capacitance = %s", self.state.phase_staggered.capacitance.shape)
    # FIXME: Update capacitance for mixed phase (enthalpy of fusion contribution)
    capacitance: npt.NDArray = (
        self.state.capacitance_staggered() * self.evaluator.mesh.basic.volume
    )

    # Heating rate (dT/dt) from flux divergence (power per unit area)
    dTdt: npt.NDArray = -delta_energy_flux / capacitance
    logger.debug("dTdt (fluxes only) = %s", dTdt)

    # Additional heating rate (dT/dt) from internal heating (power per unit mass)
    dTdt += self.state.heating * (
        self.state.phase_staggered.density() / self.state.capacitance_staggered()
    )

    logger.debug("dTdt (with internal heating) = %s", dTdt)

    return dTdt

`from_file(filename, root=Path())` `classmethod`

Parses a configuration file

Args: filename: Filename root: Root of the filename

Returns: Parameters

Source code in aragog/solver.py

@classmethod
def from_file(cls, filename: str | Path, root: str | Path = Path()) -> Self:
    """Parses a configuration file

    Args:
        filename: Filename
        root: Root of the filename

    Returns:
        Parameters
    """
    configuration_file: Path = Path(root) / Path(filename)
    logger.info("Parsing configuration file = %s", configuration_file)
    parameters: Parameters = Parameters.from_file(configuration_file)

    return cls(parameters)

`initialize()`

Initializes the model.

Source code in aragog/solver.py

def initialize(self) -> None:
    """Initializes the model."""
    logger.info("Initializing %s", self.__class__.__name__)
    self.evaluator = Evaluator(self.parameters)
    self.state = State(self.parameters, self.evaluator)

`make_tsurf_event()`

Creates a temperature event function for use with an ODE solver to monitor changes in the surface temperature.The event triggers when the change exceeds the threshold, allowing the solver to stop integration.

Returns: The event has the attributes: - terminal = True: Integration stops when the event is triggered. - direction = -1: Only triggers when the function is decreasing through zero.

Source code in aragog/solver.py

def make_tsurf_event(self):
    """
    Creates a temperature event function for use with an ODE solver to monitor changes 
    in the surface temperature.The event triggers when the change exceeds the 
    threshold, allowing the solver to stop integration.

    Returns:
        The event has the attributes:
            - terminal = True: Integration stops when the event is triggered.
            - direction = -1: Only triggers when the function is decreasing through zero.
    """
    tsurf_initial = [None]

    def tsurf_event(time: float, temperature: npt.NDArray) -> float:
        """
        Event function to detect when surface temperature changes beyond a specified threshold.

        Args:
            time (float): Current time.
            temperature (np.ndarray): Current temperature profile.

        Returns:
            float: The difference between the threshold and the actual change in surface 
                temperature. When this value crosses zero from above, the event is triggered.
        """
        tsurf_current = temperature[-1] * self.parameters.scalings.temperature
        tsurf_threshold = self.parameters.solver.tsurf_poststep_change * self.parameters.scalings.temperature

        if tsurf_initial[0] is None:
            tsurf_initial[0] = tsurf_current
            return 1.0  

        delta = abs(tsurf_current - tsurf_initial[0])

        return tsurf_threshold - delta  

    tsurf_event.terminal = self.parameters.solver.event_triggering
    tsurf_event.direction = -1

    return tsurf_event

`reset()`

This function initializes the model, while keeping the previous state of the PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data when running Aragog multiple times.

Source code in aragog/solver.py

def reset(self) -> None:
    """This function initializes the model, while keeping the previous state of the
    PhaseEvaluatorCollection object. This avoids multiple loads of lookup table data
    when running Aragog multiple times.
    """
    logger.info("Resetting %s", self.__class__.__name__)
    # Update the Evaluator object except the phase properties
    self.evaluator.mesh = Mesh(self.parameters)
    self.evaluator.boundary_conditions = BoundaryConditions(self.parameters, self.evaluator.mesh)
    self.evaluator.initial_condition = InitialCondition(self.parameters, self.evaluator.mesh, self.evaluator.phases)
    # Reinstantiate the solver state
    self.state = State(self.parameters, self.evaluator)

`DownloadLookupTableData(fname='')`

Download lookup table data

Inputs : - fname (optional) : folder name, i.e. "1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018" if not provided download all the folder list

Source code in aragog/data.py

def DownloadLookupTableData(fname: str = ""):
    """
    Download lookup table data

    Inputs :
        - fname (optional) :    folder name, i.e. "1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018"
                                if not provided download all the folder list
    """

    osf = OSF()
    project = osf.project(project_id)
    storage = project.storage("osfstorage")

    data_dir = GetFWLData() / "interior_lookup_tables"
    data_dir.mkdir(parents=True, exist_ok=True)

    # If no folder specified download all basic list
    if not fname:
        folder_list = basic_list
    elif fname in full_list:
        folder_list = [fname]
    else:
        raise ValueError(f"Unrecognised folder name: {fname}")

    for folder in folder_list:
        folder_dir = data_dir / folder
        max_tries = 2 # Maximum download attempts, could be a function argument

        if not folder_dir.exists():
            logger.info(f"Downloading interior lookup table data to {data_dir}")
            for i in range(max_tries):
                logger.info(f"Attempt {i + 1} of {max_tries}")
                success = False

                try:
                    download_zenodo_folder(folder = folder, data_dir=data_dir)
                    success = True
                except RuntimeError as e:
                    logger.error(f"Zenodo download failed: {e}")
                    folder_dir.rmdir()

                if not success:
                    try:
                        download_OSF_folder(storage=storage, folders=folder, data_dir=data_dir)
                        success = True
                    except RuntimeError as e:
                        logger.error(f"OSF download failed: {e}")

                if success:
                    break

                if i < max_tries - 1:
                    logger.info("Retrying download...")
                    sleep(5) # Wait 5 seconds before retrying
                else:
                    logger.error("Max retries reached. Download failed.")

    return

`aragog_file_logger(console_level=logging.INFO, file_level=logging.DEBUG, log_dir=os.getcwd())`

Sets up console logging and file logging according to arguments.

Returns: A logger

Source code in aragog/__init__.py

def aragog_file_logger(
    console_level=logging.INFO, file_level=logging.DEBUG, log_dir=os.getcwd()
) -> logging.Logger:
    """Sets up console logging and file logging according to arguments.

    Returns:
        A logger
    """
    # Console logger
    package_logger: logging.Logger = logging.getLogger(__name__)
    package_logger.setLevel(file_level)
    package_logger.handlers = []
    console_handler: logging.Handler = logging.StreamHandler()
    console_formatter: logging.Formatter = simple_formatter()
    console_handler.setFormatter(console_formatter)
    console_handler.setLevel(console_level)
    logger.addHandler(console_handler)
    # File logger
    log_file = os.path.join(log_dir, f"{__package__}.log")
    file_handler: logging.Handler = logging.FileHandler(log_file)
    file_formatter: logging.Formatter = complex_formatter()
    file_handler.setFormatter(file_formatter)
    file_handler.setLevel(logging.DEBUG)
    package_logger.addHandler(file_handler)

    return package_logger

`complex_formatter()`

Complex formatter for logging

Returns: Formatter for logging

Source code in aragog/__init__.py

def complex_formatter() -> logging.Formatter:
    """Complex formatter for logging

    Returns:
        Formatter for logging
    """
    fmt: str = "[%(asctime)s - %(name)-30s - %(lineno)03d - %(levelname)-9s - %(funcName)s()]"
    fmt += " - %(message)s"
    datefmt: str = "%Y-%m-%d %H:%M:%S"
    formatter: logging.Formatter = logging.Formatter(fmt, datefmt=datefmt)

    return formatter

`debug_logger()`

Sets up debug logging to the console.

Returns: A logger

Source code in aragog/__init__.py

def debug_logger() -> logging.Logger:
    """Sets up debug logging to the console.

    Returns:
        A logger
    """
    package_logger: logging.Logger = logging.getLogger(__name__)
    package_logger.setLevel(logging.DEBUG)
    logger.handlers = []
    console_handler: logging.Handler = logging.StreamHandler()
    console_formatter: logging.Formatter = simple_formatter()
    console_handler.setFormatter(console_formatter)
    logger.addHandler(console_handler)

    return package_logger

`simple_formatter()`

Simple formatter for logging

Returns: Formatter for logging

Source code in aragog/__init__.py

def simple_formatter() -> logging.Formatter:
    """Simple formatter for logging

    Returns:
        Formatter for logging
    """
    fmt: str = "[%(asctime)s - %(name)-30s - %(levelname)-9s] - %(message)s"
    datefmt: str = "%H:%M:%S"
    formatter: logging.Formatter = logging.Formatter(fmt, datefmt=datefmt)

    return formatter

API overview

Directory structure

High-level architecture

Module reference

aragog.cli

cli

all()

cli()

download()

env()

aragog.data

data

FWL_DATA_DIR = Path(os.environ.get('FWL_DATA', platformdirs.user_data_dir('fwl_data'))) module-attribute

basic_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Wolf_Bower+2018') module-attribute

full_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_400GPa', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Monteux+600', 'Melting_curves/Monteux-600', 'Melting_curves/Wolf_Bower+2018') module-attribute

logger = logging.getLogger(__name__) module-attribute

project_id = 'phsxf' module-attribute

DownloadLookupTableData(fname='')

GetFWLData()

download_OSF_folder(*, storage, folders, data_dir)

download_zenodo_folder(folder, data_dir)

get_zenodo_record(folder)

aragog.core

core

logger = logging.getLogger(__name__) module-attribute

BoundaryConditions(_parameters, _mesh) dataclass

apply_flux_boundary_conditions(state)

apply_flux_inner_boundary_condition(state)

apply_flux_outer_boundary_condition(state)

apply_temperature_boundary_conditions(temperature, temperature_basic, dTdr)

apply_temperature_boundary_conditions_melt(melt_fraction, melt_fraction_basic, dphidr)

core_cooling(state)

grey_body(state)

InitialCondition(_parameters, _mesh, _phases) dataclass

get_adiabat(pressure_basic)

get_linear()

Mesh(parameters)

dxidr property

d_dr_at_basic_nodes(staggered_quantity)

get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates)

get_constant_spacing()

get_dxidr_basic()

get_planet_density(basic_coordinates)

get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates)

quantity_at_basic_nodes(staggered_quantity)

quantity_at_staggered_nodes(basic_quantity)

Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver) dataclass

from_file(*filenames) classmethod

radionuclide_sections(parser) staticmethod

PhaseEvaluatorCollection(parameters) dataclass

State(parameters, _evaluator) dataclass

critical_reynolds_number property

heat_flux property

heating property

heating_dilatation property

heating_radio property

heating_tidal property

mass_flux property

conductive_heat_flux()

convective_heat_flux()

dilatation_heating()

gravitational_separation_mass_flux()

mixing_mass_flux()

radiogenic_heating(time)

tidal_heating()

update(temperature, time)

_BoundaryConditionsParameters(outer_boundary_condition, outer_boundary_value, inner_boundary_condition, inner_boundary_value, emissivity, equilibrium_temperature, core_heat_capacity) dataclass

scale_attributes(scalings)

_InitialConditionParameters(initial_condition=1, surface_temperature=4000, basal_temperature=4000, init_file='') dataclass

scale_attributes(scalings)

aragog.interfaces

interfaces

FloatOrArray = float | npt.NDArray module-attribute

logger = logging.getLogger(__name__) module-attribute

MixedPhaseEvaluatorProtocol

PhaseEvaluatorABC

dTdPs()

set_pressure(pressure)

set_temperature(temperature)

update()

`cli`

`all()`

`cli()`

`download()`

`env()`

`data`

`FWL_DATA_DIR = Path(os.environ.get('FWL_DATA', platformdirs.user_data_dir('fwl_data')))` `module-attribute`

`basic_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Wolf_Bower+2018')` `module-attribute`

`full_list = ('1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_400GPa', '1TPa-dK09-elec-free/MgSiO3_Wolf_Bower_2018_1TPa', 'Melting_curves/Monteux+600', 'Melting_curves/Monteux-600', 'Melting_curves/Wolf_Bower+2018')` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`project_id = 'phsxf'` `module-attribute`

`DownloadLookupTableData(fname='')`

`GetFWLData()`

`download_OSF_folder(*, storage, folders, data_dir)`

`download_zenodo_folder(folder, data_dir)`

`get_zenodo_record(folder)`

`core`

`logger = logging.getLogger(name)` `module-attribute`

`BoundaryConditions(_parameters, _mesh)` `dataclass`

`apply_flux_boundary_conditions(state)`

`apply_flux_inner_boundary_condition(state)`

`apply_flux_outer_boundary_condition(state)`

`apply_temperature_boundary_conditions(temperature, temperature_basic, dTdr)`

`apply_temperature_boundary_conditions_melt(melt_fraction, melt_fraction_basic, dphidr)`

`core_cooling(state)`

`grey_body(state)`

`InitialCondition(_parameters, _mesh, _phases)` `dataclass`

`get_adiabat(pressure_basic)`

`get_linear()`

`Mesh(parameters)`

`dxidr` `property`

`d_dr_at_basic_nodes(staggered_quantity)`

`get_basic_mass_coordinates_from_spatial_coordinates(basic_coordinates)`

`get_constant_spacing()`

`get_dxidr_basic()`

`get_planet_density(basic_coordinates)`

`get_staggered_spatial_coordinates_from_mass_coordinates(staggered_mass_coordinates)`

`quantity_at_basic_nodes(staggered_quantity)`

`quantity_at_staggered_nodes(basic_quantity)`

`Parameters(*, boundary_conditions, energy, initial_condition, mesh, phase_solid, phase_liquid, phase_mixed, radionuclides, scalings, solver)` `dataclass`

`from_file(*filenames)` `classmethod`

`radionuclide_sections(parser)` `staticmethod`

`PhaseEvaluatorCollection(parameters)` `dataclass`

`State(parameters, _evaluator)` `dataclass`

`critical_reynolds_number` `property`

`heat_flux` `property`

`heating` `property`

`heating_dilatation` `property`

`heating_radio` `property`

`heating_tidal` `property`

`mass_flux` `property`

`conductive_heat_flux()`

`convective_heat_flux()`

`dilatation_heating()`

`gravitational_separation_mass_flux()`

`mixing_mass_flux()`

`radiogenic_heating(time)`

`tidal_heating()`

`update(temperature, time)`

`_BoundaryConditionsParameters(outer_boundary_condition, outer_boundary_value, inner_boundary_condition, inner_boundary_value, emissivity, equilibrium_temperature, core_heat_capacity)` `dataclass`

`scale_attributes(scalings)`

`_InitialConditionParameters(initial_condition=1, surface_temperature=4000, basal_temperature=4000, init_file='')` `dataclass`

`scale_attributes(scalings)`

`interfaces`

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`MixedPhaseEvaluatorProtocol`

`PhaseEvaluatorABC`

`dTdPs()`

`set_pressure(pressure)`

`set_temperature(temperature)`

`update()`

`PhaseEvaluatorProtocol`

`PropertyProtocol`

`mesh`

`FloatOrArray = float | npt.NDArray` `module-attribute`

`logger = logging.getLogger(name)` `module-attribute`

`AdamsWilliamsonEOS(settings, basic_radii)`

`basic_density` `property`

`basic_pressure` `property`