Approximations

approximations

This module provides classes that emulate physical approximations of fluid dynamics systems by exposing methods to calculate specific terms in the corresponding mathematical equations. Users instantiate the appropriate class by providing relevant parameters and pass the instance to other objects, such as solvers. Under the hood, G-ADOPT queries variables and methods from the approximation.

BaseApproximation

Bases: ABC

Base class to provide expressions for the coupled Stokes and Energy system.

The basic assumption is that we are solving (to be extended when needed)

div(dev_stress) + grad p + buoyancy(T, p) * khat = 0
div(rho_continuity * u) = 0
rhocp DT/Dt + linearized_energy_sink(u) * T
  = div(kappa * grad(Tbar + T)) + energy_source(u)

where the following terms are provided by Approximation methods:

• linearized_energy_sink(u) = 0 (BA), Di * rhobar * alphabar * g * w (EBA), or Di * rhobar * alphabar * w (TALA/ALA)
• kappa() is diffusivity or conductivity depending on rhocp()
• Tbar is 0 or reference temperature profile (ALA)
• dev_stress depends on the compressible property (False or True):
  • if compressible then dev_stress = mu * [sym(grad(u) - 2/3 div(u)]
  • if not compressible then dev_stress = mu * sym(grad(u)) and rho_continuity is assumed to be 1

compressible: bool abstractmethod property

Defines compressibility.

Returns:

Type	Description
bool	A boolean signalling if the governing equations are in compressible form.

Tbar: Function abstractmethod property

Defines the reference temperature profile.

Returns:

Type	Description
Function	A Firedrake function for the reference temperature profile.

stress(u) abstractmethod

Defines the deviatoric stress.

Returns:

Type	Description
Expr	A UFL expression for the deviatoric stress.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def stress(self, u: Function) -> ufl.core.expr.Expr:
    """Defines the deviatoric stress.

    Returns:
      A UFL expression for the deviatoric stress.

    """
    pass

buoyancy(p, T) abstractmethod

Defines the buoyancy force.

Returns:

Type	Description
Expr	A UFL expression for the buoyancy term (momentum source in gravity direction).

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def buoyancy(self, p: Function, T: Function) -> ufl.core.expr.Expr:
    """Defines the buoyancy force.

    Returns:
      A UFL expression for the buoyancy term (momentum source in gravity direction).

    """
    pass

rho_continuity() abstractmethod

Defines density.

Returns:

Type	Description
Expr	A UFL expression for density in the mass continuity equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def rho_continuity(self) -> ufl.core.expr.Expr:
    """Defines density.

    Returns:
      A UFL expression for density in the mass continuity equation.

    """
    pass

rhocp() abstractmethod

Defines the volumetric heat capacity.

Returns:

Type	Description
Expr	A UFL expression for the volumetric heat capacity in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def rhocp(self) -> ufl.core.expr.Expr:
    """Defines the volumetric heat capacity.

    Returns:
      A UFL expression for the volumetric heat capacity in the energy equation.

    """
    pass

kappa() abstractmethod

Defines thermal diffusivity.

Returns:

Type	Description
Expr	A UFL expression for thermal diffusivity.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def kappa(self) -> ufl.core.expr.Expr:
    """Defines thermal diffusivity.

    Returns:
      A UFL expression for thermal diffusivity.

    """
    pass

linearized_energy_sink(u) abstractmethod

Defines temperature-related sink terms.

Returns:

Type	Description
Expr	A UFL expression for temperature-related sink terms in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def linearized_energy_sink(self, u) -> ufl.core.expr.Expr:
    """Defines temperature-related sink terms.

    Returns:
      A UFL expression for temperature-related sink terms in the energy equation.

    """
    pass

energy_source(u) abstractmethod

Defines additional terms.

Returns:

Type	Description
Expr	A UFL expression for additional independent terms in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def energy_source(self, u) -> ufl.core.expr.Expr:
    """Defines additional terms.

    Returns:
      A UFL expression for additional independent terms in the energy equation.

    """
    pass

free_surface_terms(p, T, eta, theta_fs) abstractmethod

Defines free surface normal stress in the momentum equation and prefactor multiplying the free surface equation.

The normal stress depends on the density contrast across the free surface. Depending on the variable_rho_fs argument, this contrast either involves the interior reference density or the full interior density that accounts for changes in temperature and composition within the domain. For dimensional simulations, the user should specify delta_rho_fs, defined as the difference between the reference interior density and the exterior density. For non-dimensional simulations, the user should specify RaFS, the equivalent of the Rayleigh number for the density contrast across the free surface.

The prefactor multiplies the free surface equation to ensure the top-right and bottom-left corners of the block matrix remain symmetric. This is similar to rescaling eta -> eta_tilde in Kramer et al. (2012, see block matrix shown in Eq. 23).

Returns:

Type	Description
tuple[Expr]	A UFL expression for the free surface normal stress and a UFL expression for
tuple[Expr]	the free surface equation prefactor.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def free_surface_terms(self, p, T, eta, theta_fs) -> tuple[ufl.core.expr.Expr]:
    """Defines free surface normal stress in the momentum equation and prefactor
    multiplying the free surface equation.

    The normal stress depends on the density contrast across the free surface.
    Depending on the `variable_rho_fs` argument, this contrast either involves the
    interior reference density or the full interior density that accounts for
    changes in temperature and composition within the domain. For dimensional
    simulations, the user should specify `delta_rho_fs`, defined as the difference
    between the reference interior density and the exterior density. For
    non-dimensional simulations, the user should specify `RaFS`, the equivalent of
    the Rayleigh number for the density contrast across the free surface.

    The prefactor multiplies the free surface equation to ensure the top-right and
    bottom-left corners of the block matrix remain symmetric. This is similar to
    rescaling eta -> eta_tilde in Kramer et al. (2012, see block matrix shown in
    Eq. 23).

    Returns:
      A UFL expression for the free surface normal stress and a UFL expression for
      the free surface equation prefactor.

    """
    pass

BoussinesqApproximation(Ra, *, mu=1, rho=1, alpha=1, T0=0, g=1, RaB=0, delta_rho=1, kappa=1, H=0)

Bases: BaseApproximation

Expressions for the Boussinesq approximation.

Density variations are considered small and only affect the buoyancy term. Reference parameters are typically constant. Viscous dissipation is neglected (Di << 1).

Parameters:

Name	Type	Description	Default
Ra	Function | Number	Rayleigh number	required
mu	Function | Number	dynamic viscosity	1
rho	Function | Number	reference density	1
alpha	Function | Number	coefficient of thermal expansion	1
T0	Function | Number	reference temperature	0
g	Function | Number	gravitational acceleration	1
RaB	Function | Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers	0
delta_rho	Function | Number	compositional density difference from the reference density	1
kappa	Function | Number	thermal diffusivity	1
H	Function | Number	internal heating rate	0

Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    Ra: Function | Number,
    *,
    mu: Function | Number = 1,
    rho: Function | Number = 1,
    alpha: Function | Number = 1,
    T0: Function | Number = 0,
    g: Function | Number = 1,
    RaB: Function | Number = 0,
    delta_rho: Function | Number = 1,
    kappa: Function | Number = 1,
    H: Function | Number = 0,
):
    self.Ra = ensure_constant(Ra)
    self.mu = ensure_constant(mu)
    self.rho = ensure_constant(rho)
    self.alpha = ensure_constant(alpha)
    self.T0 = T0
    self.g = ensure_constant(g)
    self.kappa_ref = ensure_constant(kappa)
    self.RaB = RaB
    self.delta_rho = ensure_constant(delta_rho)
    self.H = ensure_constant(H)

ExtendedBoussinesqApproximation(Ra, Di, *, H=None, **kwargs)

Bases: BoussinesqApproximation

Expressions for the extended Boussinesq approximation.

Extends the Boussinesq approximation by including viscous dissipation and work against gravity (both scaled with Di).

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
mu		dynamic viscosity	required
H	Optional[Number]	volumetric heat production	None

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	coefficient of thermal expansion
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	thermal diffusivity

Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

Source code in g-adopt/gadopt/approximations.py

def __init__(self, Ra: Number, Di: Number, *, H: Optional[Number] = None, **kwargs):
    super().__init__(Ra, **kwargs)
    self.Di = Di
    self.H = H

TruncatedAnelasticLiquidApproximation(Ra, Di, *, Tbar=0, cp=1, **kwargs)

Bases: ExtendedBoussinesqApproximation

Truncated Anelastic Liquid Approximation

Compressible approximation. Excludes linear dependence of density on pressure.

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
Tbar	Function | Number	reference temperature. In the diffusion term we use Tbar + T (i.e. T is the pertubartion)	0
cp	Function | Number	reference specific heat at constant pressure	1

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	reference thermal expansion coefficient
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	diffusivity
mu	Number	viscosity used in viscous dissipation
H	Number	volumetric heat production

Note

Other keyword arguments may be depth-dependent, but default to 1 if not supplied.

Source code in g-adopt/gadopt/approximations.py

def __init__(self,
             Ra: Number,
             Di: Number,
             *,
             Tbar: Function | Number = 0,
             cp: Function | Number = 1,
             **kwargs):
    super().__init__(Ra, Di, **kwargs)
    self.Tbar = Tbar
    self.cp = cp

AnelasticLiquidApproximation(Ra, Di, *, chi=1, gamma0=1, cp0=1, cv0=1, **kwargs)

Bases: TruncatedAnelasticLiquidApproximation

Anelastic Liquid Approximation

Compressible approximation. Includes linear dependence of density on pressure.

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
chi	Function | Number	reference isothermal compressibility	1
gamma0	Function | Number	Gruneisen number (in pressure-dependent buoyancy term)	1
cp0	Function | Number	specific heat at constant pressure, reference for entire Mantle (in pressure-dependent buoyancy term)	1
cv0	Function | Number	specific heat at constant volume, reference for entire Mantle (in pressure-dependent buoyancy term)	1

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	reference thermal expansion coefficient
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	diffusivity
mu	Number	viscosity used in viscous dissipation
H	Number	volumetric heat production
Tbar	Number	reference temperature. In the diffusion term we use Tbar + T (i.e. T is the pertubartion)
cp	Number	reference specific heat at constant pressure

Source code in g-adopt/gadopt/approximations.py

def __init__(self,
             Ra: Number,
             Di: Number,
             *,
             chi: Function | Number = 1,
             gamma0: Function | Number = 1,
             cp0: Function | Number = 1,
             cv0: Function | Number = 1,
             **kwargs):
    super().__init__(Ra, Di, **kwargs)
    # Dynamic pressure contribution towards buoyancy
    self.chi = chi
    self.gamma0, self.cp0, self.cv0 = gamma0, cp0, cv0

SmallDisplacementViscoelasticApproximation(density, shear_modulus, viscosity, *, g=1)

Expressions for the small displacement viscoelastic approximation.

By assuming a small displacement, we can linearise the problem, assuming a perturbation away from a reference state. We assume a Maxwell viscoelastic rheology, i.e. stress is the same but viscous and elastic strains combine linearly. We follow the approach by Zhong et al. 2003 redefining the problem in terms of incremental displacement, i.e. velocity * dt where dt is the timestep. This produces a mixed stokes system for incremental displacement and pressure which can be solved in the same way as mantle convection (where unknowns are velocity and pressure), with a modfied viscosity and stress term accounting for the deviatoric stress at the previous timestep.

Zhong, Shijie, Archie Paulson, and John Wahr. "Three-dimensional finite-element modelling of Earth’s viscoelastic deformation: effects of lateral variations lithospheric thickness." Geophysical Journal International 155.2 (2003): 679-695.

N.b. that the implentation currently assumes all terms are dimensional.

Parameters:

Name	Type	Description	Default
density	Function | Number	background density	required
shear_modulus	Function | Number	shear modulus	required
viscosity	Function | Number	viscosity	required
g	Function | Number	gravitational acceleration	1

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    density: Function | Number,
    shear_modulus: Function | Number,
    viscosity: Function | Number,
    *,
    g: Function | Number = 1,
):

    self.density = ensure_constant(density)
    self.shear_modulus = ensure_constant(shear_modulus)
    self.viscosity = ensure_constant(viscosity)
    self.g = ensure_constant(g)

    self.maxwell_time = viscosity / shear_modulus