Development manual

Contributing

If you are interested in contributing to the model, please contact the developers using the information on the Home page.

Coding style

Indentation uses 4 spaces, no tabs.
Function names should be lowercase, with words separated by underscores .
Lines should aim to have a length of no more than 92 characters.
All functions should have docstrings, ideally with Arguments and Returns listed.
More comments are always better, even if they seem redundant.
Use type hinting where possible.
Print statements should be made through the logger where possible.

Code reference

Obliqua.solid0d.compute_solid_lovenumbers — Method

compute_solid_lovenumbers(μc, mass_tot, R, n)

Calculate k2 lovenumbers in solid.

Arguments

μc::precc : Complex shear modulus.
mass_tot::prec : Total mass of the planet.
R::prec : Outer radius of the planet.
n::Int=2 : Power of the radial factor (goes with (r/a)^{n}, since r<<a only n=2 contributes significantly).

Returns

k2_T::precc : Complex Tidal k2 Lovenumber.
k2_L::precc : Complex Load k2 Lovenumber.

source

Obliqua.solid0d.mean_cmu — Method

mean_cmu(μc, r)

Calculate the Hill-averaged complex shear modulus in a spherical segment (Hill+1952).

Arguments

μc::Array{precc,1} : Complex shear modulus profile of the planet.
r::Array{prec,1} : Radial positions of layers, from bottom to top of segment.

Returns

mean_μc::precc : Mean complex shear modulus in segment.

Notes

(DOI 10.1088/0370-1298/65/5/307)

source

Obliqua.solid1d.Ynm — Method

Ynm(n, m, theta, phi)

Compute the spherical harmonic Ynm for given n, m, theta, and phi.

Arguments

n::Int : Tidal degree.
m::Int : Tidal order.
theta::AbstractArray : Array of colatitudes in radians.
phi::LinearAlgebra.Adjoint : Array of longitudes in radians.

Returns

Ynm::Array{ComplexF64,2} : 2D array of spherical harmonic values for each combination of theta and phi.

source

Obliqua.solid1d.compute_M — Method

compute_M(ω, r, ρ, g, μ, K, n, ρ_core, μ_core, κ_core; core="liquid")

Compute the M matrix, which is used to propagate the solution across the entire interior. This is used in the compute_y function.

Arguments

ω::prec : Forcing frequency.
r::Array{prec,2} : 2D array of layer boundaries.
ρ::Array{prec,1} : 1D array of layer densities.
g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
μ::Array{precc,1} : 1D array of layer shear moduli.
K::Array{prec,1} : 1D array of layer bulk moduli.
n::Int : Tidal degree.
ρ_core::prec : Density of the core, which is used to compute the starting vector for the numerical integration across the interior.
μ_core::prec : Shear modulus of the core.
κ_core::prec : Bulk modulus of the core.

Keyword Arguments

core::String="liquid" : Type of core, either "liquid" or "solid". This is used to compute the starting vector for the numerical integration across the interior.

Returns

M::Array{precc,2} : 3x3 M matrix, which is used to propagate the solution across the entire interior.
y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.

source

Obliqua.solid1d.compute_strain_ten! — Method

compute_strain_ten!(ϵ, y, n, rr, ρr, gr, μr, Ksr)

Calculate the strain tensor ϵ at a particular radial level.

Arguments

ϵ::SubArray : 3D array to store the strain tensor at a particular radial level, with dimensions corresponding to latitude, longitude, and the 6 independent components of the strain tensor.
y::SubArray : 1D array of the solution vector y at a particular radial level, with 6 components.
n::Int : Tidal degree.
rr::Float64 : Radius at which to compute the strain tensor.
ρr::Float64 : Density at radius rr.
gr::Float64 : Gravity at radius rr.
μr::ComplexF64 : Shear modulus at radius rr.
Ksr::Float64 : Bulk modulus at radius rr.

source

Obliqua.solid1d.compute_y — Method

compute_y(r, g, M, y1_4, n; load=false)

Compute the solution vector y across the entire interior, given the M matrix and the y1_4 solutions across each layer. This is used to compute the strain tensor and heating profile.

Arguments

r::Array{prec,2} : 2D array of layer boundaries.
g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
M::Array{precc,2} : 3x3 M matrix, which is used to propagate the solution across the entire interior.
y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.
n::Int : Tidal degree.

Keyword Arguments

load::Bool=false : If true, compute the solution for a loaded problem.

Returns

y::Array{ComplexF64,3} : 3D array of the solution vector y across the interior.

source

Obliqua.solid1d.convert_params_to_prec — Method

convert_params_to_prec(r, ρ, g, μ, κs)

Convert input parameters into the required precision.

Arguments

r::Array{Float64,2} : 2D array of layer boundaries.
ρ::Array{Float64,1} : 1D array of layer densities.
g::Array{Float64,2} : 2D array of gravity values at the layer boundaries.
μ::Array{Float64,1} : 1D array of layer shear moduli.
κs::Array{Float64,1} : 1D array of layer bulk moduli.

Returns

Tuple of converted parameters in the required precision.

source

Obliqua.solid1d.define_spherical_grid — Method

define_spherical_grid(res, n, m)

Create the spherical grid of angular resolution res in degrees. This is used for numerical integrations over solid angle. A new grid can easily be defined by recalling the function with a new res.

Arguments

res::Float64 : Desired angular resolution in degrees.
n::Int : Tidal degree.
m::Int : Tidal order.

Notes

The grid is internal to solid1d, but can be accessed with

solid1d.clats[:] # colatitude grid
solid1d.lons[:]  # longitude grid

source

Obliqua.solid1d.expand_layers — Method

expand_layers(r; nr::Int=80)

Discretize the primary layers given by r into nr discrete secondary layers.

Arguments

r::Array{prec,1} : 1D array of primary layer boundaries.

Keyword Arguments

nr::Int=80 : Number of secondary layers to discretize.

Returns

rs::Array{prec,2} : 2D array of secondary layer boundaries/

source

Obliqua.solid1d.get_B! — Method

get_B!(B, ω, r1, r2, g1, g2, ρ, μ, K, n)

Compute the 6x6 numerical integrator matrix, which integrates dy/dr from r1 to r2 for the solid-body problem. B here represnts the RK4 integrator, given by Eq. S5.5 in Hay et al., (2025).

Arguments

B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.
ω::prec : Forcing frequency.
r1::prec : Starting radius for integration.
r2::prec : Ending radius for integration.
g1::prec : Gravity at radius r1.
g2::prec : Gravity at radius r2.
ρ::prec : Density at radius r.
μ::precc : Shear modulus at radius r.
K::prec : Bulk modulus at radius r.
n::Int : Tidal degree.

Notes

See also get_B

source

Obliqua.solid1d_mush.get_B — Method

get_B(ω, r1, r2, g1, g2, ρ, μ, K, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n)

Compute the 8x8 numerical integrator matrix, which integrates dy/dr from r1 to r2 for the two-phase problem.

Arguments

ω::prec : Forcing frequency.
r1::prec : Starting radius for integration.
r2::prec : Ending radius for integration.
g1::prec : Gravity at radius r1.
g2::prec : Gravity at radius r2.
ρ::prec : Density at radius r.
μ::precc : Shear modulus at radius r.
K::prec : Bulk modulus at radius r.
ρₗ::prec : Liquid density at radius r.
Kl::prec : Liquid bulk modulus at radius r.
Kd::prec : Drained bulk modulus at radius r.
α::prec : Biot coefficient at radius r.
ηₗ::prec : Liquid viscosity at radius r.
ϕ::prec : Porosity at radius r.
k::prec : Permeability at radius r.
n::Int : Tidal degree.

Returns

B::Array{precc,2} : 8x8 numerical integrator matrix for integrating dy/dr from r1 to r2 for the two-phase problem.

Notes

See 'get_B!' for definition.

source

Obliqua.solid1d_mush.get_B — Method

get_B(ω, r1, r2, g1, g2, ρ, μ, K, n)

Compute the 6x6 numerical integrator matrix, which integrates dy/dr from r1 to r2 for the solid-body problem.

Arguments

ω::prec : Forcing frequency.
r1::prec : Starting radius for integration.
r2::prec : Ending radius for integration.
g1::prec : Gravity at radius r1.
g2::prec : Gravity at radius r2.
ρ::prec : Density at radius r.
μ::precc : Shear modulus at radius r.
K::prec : Bulk modulus at radius r.
n::Int : Tidal degree.

Returns

B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.

Notes

See 'get_B!' for definition.

source

Obliqua.solid1d_mush.get_B_product! — Method

get_B_product!(Bprod2, ω, r, ρ, g, μ, K, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n)

Compute the product of the 8x8 B matrices within a primary layer. This is used to propgate the y solution across a single two-phase primary layer. Bprod is denoted by D in Eq. S5.14 in Hay et al., (2025).

Arguments

Bprod2::Array{precc} : 8x8x(nr-1)x(nlayers-1) array to store the B products across each secondary layer within each primary layer.
ω::prec : Forcing frequency.
r::SubArray{prec} : 1D array of layer boundaries.
ρ::prec : 1D array of layer densities.
g::SubArray{prec} : 1D array of gravity values at the layer boundaries.
μ::precc : 1D array of layer shear moduli.
K::prec : 1D array of layer bulk moduli.
ρₗ::prec : 1D array of liquid densities at layer boundaries.
Kl::prec : 1D array of liquid bulk moduli at layer boundaries.
Kd::prec : 1D array of drained bulk moduli at layer boundaries.
α::prec : 1D array of Biot coefficients at layer boundaries.
ηₗ::prec : 1D array of liquid viscosities at layer boundaries.
ϕ::prec : 1D array of porosities at layer boundaries.
k::prec : 1D array of permeabilities at layer boundaries.
n::Int : Tidal degree.

source

Obliqua.solid1d_mush.get_g — Method

get_g(r, ρ, m_core)

Compute the radial gravity structure associated with a density profile r at intervals given by r.

Arguments

r::Array{prec,2} : 2D array of layer boundaries.
ρ::Array{prec,1} : 1D array of layer densities. The length of ρ must be equal to the number of columns in r.
m_core::prec : Mass of the planetary core.

Returns

g::Array{prec,2} : 2D array of gravity values at the layer boundaries. The dimensions of g are the same as r.

Notes

r must be be a 2D array, with index 1 representing the top radius of secondary layers, and index 2 representing the top radius of primary layers.

source

Obliqua.solid1d_mush.get_heating_profile — Method

get_heating_profile(y, r, ρ, g, μ, Ks, ω, ρl, Kl, Kd, α, ηl, ϕ, k, n; lay=nothing)

Get the radial volumetric heating for two-phase tides and eccentricity forcing, assuming synchronous rotation. Heating rate is computed with numerical integration using the solution y returned by compute_y, using Eq. 2.39a/b/c integrated over solid angle. The heating profile for a specific layer is specified with lay, otherwise all layers will be caclulated.

Arguments

y::Array{ComplexF64} : 4D array of the solution vector y across the interior, returned by compute_y.
r::Matrix : 2D array of layer boundaries.
ρ::AbstractVector : 1D array of layer densities.
g::Matrix : 2D array of gravity values at the layer boundaries.
μ::AbstractVector : 1D array of layer shear moduli.
Ks::AbstractVector : 1D array of layer bulk moduli for shear dissipation.
ω::prec : Tidal frequency in radians per second.
ρl::AbstractVector : 1D array of liquid densities at layer boundaries.
Kl::AbstractVector : 1D array of liquid bulk moduli at layer boundaries.
Kd::AbstractVector : 1D array of drained bulk moduli at layer boundaries.
α::AbstractVector : 1D array of Biot coefficients at layer boundaries.
ηl::AbstractVector : 1D array of liquid viscosities at layer boundaries.
ϕ::AbstractVector : 1D array of porosities at layer boundaries.
k::AbstractVector : 1D array of permeabilities at layer boundaries.
n::Int : Tidal degree.

Keyword Arguments

lay::Int=nothing : If specified, compute the heating profile for only the layer corresponding to this index. Otherwise, compute for all layers.

Returns

Eμ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to shear, in W.
Eμ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to shear, with dimensions corresponding to the secondary layer and primary layer indices.
Eκ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to compaction, in W.
Eκ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to compaction, with dimensions corresponding to the secondary layer and primary layer indices.
El_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to Darcy flow, in W.
El_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to Darcy flow, with dimensions corresponding to the secondary layer and primary layer indices.

source

Obliqua.solid1d_relax.compute_y — Method

compute_y(r, ρ, g, μ, K, ω, n, R, ρ_core, μ_core, κ_core, M_tot; core="liquid")

Compute the solution y to the solid-body problem using a relaxation method. This function performs the forward-backward relaxation scheme described in the main text of N. Kobayashi (2006), where we first solve the radial system of ODEs to obtain the solution at the surface, and then perform back-substitution to compute the solution at all interior points.

Arguments

r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Shear modulus of the core, used for core boundary conditions.
κ_core::prec : Bulk modulus of the core, used for core boundary conditions.
M_tot::prec : Total mass of the planet, used for non-dimensionalization.

Keyword Arguments

`core::String="liquid" : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.

Returns

y_t::Matrix{ComplexF64} : 6xN matrix of the solution at all radial grid points, where N is the number of radial layers. Each column corresponds to a radial grid point, and each row corresponds to a state variable (displacements, stresses, potential).
y_l::Matrix{ComplexF64} : 6x1 matrix of the solution at the surface for the load problem.

source

Obliqua.solid1d_relax.core_boundary — Method

core_boundary(R, ids, r, ρ, g, μ, K, ω, ρ_core, μ_core, κ_core, core, n; G0=1)

Perform the forward-backward relaxation step at the core boundary. This function implements the recursion described in N. Kobayashi (2007) for the initial step of the relaxation scheme, where we apply the core boundary condition and get the first solution for the first layer above the core.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
`μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Shear modulus of the core, used for core boundary conditions.
κ_core::prec : Bulk modulus of the core, used for core boundary conditions.
core::String : Type of core boundary condition to apply.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.

Returns

C1l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the C1 matrix for the next iteration.
D2l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the D2 matrix for the next iteration.

source

Obliqua.solid1d_relax.doublefactorial — Method

doublefactorial(n)

Compute the double factorial of an integer n, defined as n!! = n * (n-2) * (n-4) * ... until 1 or 0.

Arguments

n::Integer : The integer for which to compute the double factorial. Must be non-negative.

Returns

result::Integer : The double factorial of n.

source

Obliqua.solid1d_relax.get_core_bc! — Method

get_core_bc!(ω, r, ρ, g, μ, K, type, n; G0=1)

Get the core boundary condition matrix B for the solid-body problem. The core boundary conditions are derived from the requirement that the radial stress at the core-mantle boundary must balance the tidal potential, and that the tangential stresses must vanish.

Arguments

ω::prec : Forcing frequency.
r::prec : Radial position of the core-mantle boundary.
ρ::prec : Average core density.
g::prec : Gravity at the core-mantle boundary.
μ::prec : Average core shear modulus.
K::prec : Average core bulk modulus.
type::String : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.

Returns

B::Array{precc,2} : 3x6 matrix representing the linear relationship between the state variables at the core and the boundary conditions.

source

Obliqua.solid1d_relax.get_g — Method

get_g(r, ρ, m_core)

Compute the radial gravity structure associated with a density profile r at intervals given by r.

Arguments

r::Array{prec,1} : 1D array of layer boundaries.
ρ::Array{prec,1} : 1D array of layer densities. The length of ρ must be equal to the number of columns in r.
m_core::prec : Mass of the core.

Returns

g::Array{prec,1} : 1D array of gravity values at the layer boundaries. The dimensions of g are the same as r.
M_enc::prec : Total mass enclosed within the outermost layer boundary.

source

Obliqua.solid1d_relax.get_surface_bc! — Method

get_surface_bc!(R, g, n, U, U_prime, tau, P; G0=1)

Get the surface boundary condition vector b and matrix BN for the solid-body problem. The surface boundary conditions are determined by setting, respectively (U, U', tau, P) to (1,0,0,0) for tidal Love number and (0,1,0,0) for load Love number in system.

https://hal.science/hal-03421553/document

Arguments

R::prec : Mantle radius, used for surface boundary conditions.
g::prec : Gravity at the surface, used for surface boundary conditions.
n::Int : Tidal degree.
U::Int : Tidal potential at the surface.
U_prime::Int : Radial derivative of the tidal potential at the surface.
tau::Int : Tangential tidal stress at the surface.
P::Int : Surface mass load at the surface.

Keyword Arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.

Returns

BN::Array{precc,2} : 3x6 matrix representing the linear relationship between the state variables at the surface and the boundary conditions.
b::Vector{precc} : Vector of length 6 representing the inhomogeneous part of the surface boundary conditions.

source

Obliqua.solid1d_relax.propagate_solid — Method

propagate_solid(R, B, Cn_l, Dnp_l, ids, r, ρ, g, μ, K, ω, n; G0=1)

Perform the forward-backward relaxation step for the solid propagation segments. This function implements the recursion described in N. Kobayashi (2007) for the segments of the radial grid that correspond to the solid layers, where we propagate the solution up to the surface using the 6x6 system of equations.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
B::Vector{Matrix{precc}} : Vector of 8x1 matrices representing the inhomogeneous terms of the ODE system at each radial layer.
Cn_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Cn matrix from the previous step.
Dnp_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Dnp matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.

Returns

Cn_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Cn matrix for the next iteration.
Dnp_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Dnp matrix for the next iteration.

source

Obliqua.solid1d_relax.resample_profiles — Method

resample_profiles(radius, rho, visc, shear, bulk, m_core, dr_min, dr_max)

Resample the input profiles onto a new grid with ncalc points. The new grid is generated using a stretched and refined scheme, which allows for better resolution in regions of interest (e.g., near layer boundaries).

Arguments

radius::Vector{prec} : Original radius profile (layer boundaries).
rho::Vector{prec} : Original density profile (defined at layer centers).
visc::Vector{prec} : Original viscosity profile (defined at layer centers).
shear::Vector{precc} : Original shear modulus profile (defined at layer centers).
bulk::Vector{prec} : Original bulk modulus profile (defined at layer centers).
m_core::prec : Mass of the core, used for gravity calculations.
Δr_min::Int64 : Minimum grid spacing for the new grid.
Δr_max::Int64 : Maximum grid spacing for the new grid.

Returns

Tuple of resampled profiles on the new grid:

r_new_b::Vector{prec} : New radius profile at layer boundaries.
ρ_new::Vector{prec} : New density profile at layer centers.
η_new::Vector{prec} : New viscosity profile at layer centers.
μ_new::Vector{precc} : New shear modulus profile at layer centers.
κ_new::Vector{prec} : New bulk modulus profile at layer centers.
g_new::Vector{prec} : New gravity profile at layer centers.
M_tot::Float64 : Total mass enclosed within the outermost layer boundary.

source

Obliqua.solid1d_relax.solve_radial_system — Method

solve_radial_system(r, ρ, g, μ, K, ω, n, R_planet, ρ_core, μ_core, κ_core, M_tot; core="liquid")

Solve the radial system of ODEs for the solid-body problem using a relaxation method. This function implements the forward-backward relaxation scheme described in the main text of N. Kobayashi (2006).

Arguments

r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Shear modulus of the core, used for core boundary conditions.
κ_core::prec : Bulk modulus of the core, used for core boundary conditions.
M_tot::prec : Total mass of the planet, used for gravity calculations.

Keyword Arguments

`core::String="liquid" : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.

Returns

y_t::Vector{precc} : Vector of length 6 representing the tidal solution at the surface (radius = R_planet). This includes the displacements, stresses, and potential at the surface.
y_l::Vector{precc} : Vector of length 6 representing the load solution at the surface (radius = R_planet). This includes the displacements, stresses, and potential at the surface.
R::Vector{Matrix{precc}} : Vector of 6x6 matrices representing the coefficients of the ODE system at each radial layer.
S::Vector{Matrix{precc}} : Vector of 6x6 matrices representing the normalization.

source

Obliqua.solid1d_relax.surface_boundary — Method

surface_boundary(R, CNm_l, DN_l, ids, r, ρ, g, μ, K, ω, n; G0=1)

Perform the forward-backward relaxation step at the surface boundary. This function implements the recursion described in N. Kobayashi (2007) for the final step of the relaxation scheme, where we apply the surface boundary condition and solve for the final solution at the surface, using the 6x6 system of equations.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
CNm_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the CNm matrix from the previous step.
DN_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the DN matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.

Returns

y_t::Matrix{precc} : 6x1 matrix representing the solution at the surface for the tidal problem.
y_l::Matrix{precc} : 6x1 matrix representing the solution at the surface for the load problem.

source

Obliqua.solid1d_mush_relax.compute_y — Method

compute_y(r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n, ρ_core, μ_core, κ_core, M_tot; core="liquid")

Arguments

r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.
n::Int : Tidal degree.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Complex shear modulus of the core, used for core boundary conditions.
κ_core::prec : Complex bulk modulus of the core, used for core boundary conditions.
M_tot::prec : Total mass of the planet.

Keyword Arguments

`core::String="liquid" : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.

Returns

y_t::Matrix{ComplexF64} : 8xN matrix of the solution at all radial grid points, where N is the number of radial layers. Each column corresponds to a radial grid point, and each row corresponds to a state variable (displacements, stresses, potential).
y_l::Matrix{ComplexF64} : 8x1 matrix of the solution at the surface for the load problem.

source

Obliqua.solid1d_mush_relax.core_boundary — Method

core_boundary(R, ids, r, ρ, g, μ, K, ω, ρ_core, μ_core, κ_core, core, n; G0=1, Y=[1,2,3,4,5,6])

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Complex shear modulus of the core, used for core boundary conditions.
κ_core::prec : Complex bulk modulus of the core, used for core boundary conditions.
core::String : Type of core boundary condition to apply.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6] : Ordering of the solution vector components.

Returns

C1l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the C1 matrix for the next iteration.
D2l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the D2 matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.core_boundary_mush — Method

core_boundary_mush(R, ids, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, ρ_core, μ_core, κ_core, core, n; G0=1, Y=[1,2,3,4,5,6,7,8])

Perform the forward-backward relaxation step at the core boundary for the two-phase problem. This function implements the recursion described in N. Kobayashi (2007) for the initial step of the relaxation scheme, where we apply the core boundary condition and get the first solution for the first layer above the core, but now using the full 8x8 system that includes the porous layer components.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Complex shear modulus of the core, used for core boundary conditions.
κ_core::prec : Complex bulk modulus of the core, used for core boundary conditions.
core::String : Type of core boundary condition to apply.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6,7,8]: Ordering of the solution vector components.

Returns

C1l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the C1 matrix for the next iteration.
D2l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the D2 matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.get_core_bc! — Method

get_core_bc!(ω, r, ρ, g, μ, K, type, n; G0=1, Y=[1,2,3,4,5,6])

Arguments

ω::Float64 : Forcing frequency.
r::prec : Radial position of the core-mantle boundary.
ρ::prec : Average core density.
g::prec : Gravity at the core-mantle boundary.
μ::prec : Average core shear modulus.
K::prec : Average core bulk modulus.
type::String : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.
Y::Vector{Int}=[1,2,3,4,5,6] : Indices for the state variables (default is for standard case).

Returns

B::Array{precc,2} : 3x6 matrix representing the linear relationship between the state variables at the core and the boundary conditions.

source

Obliqua.solid1d_mush_relax.get_g — Method

get_g(r, ρ, m_core)

Compute the radial gravity structure associated with a density profile r at intervals given by r.

Arguments

r::Array{prec,1} : 1D array of layer boundaries.
ρ::Array{prec,1} : 1D array of layer densities. The length of ρ must be equal to the number of columns in r.
m_core::prec : Mass of the core.

Returns

g::Array{prec,1} : 1D array of gravity values at the layer boundaries. The dimensions of g are the same as r.
M_enc::prec : Total mass enclosed within the outermost layer boundary.

source

Obliqua.solid1d_mush_relax.get_surface_bc! — Method

get_surface_bc!(R, g, n, U, U_prime, tau, P; G0=1, Y=[1,2,3,4,5,6,7,8])

https://hal.science/hal-03421553/document

Arguments

R::prec : Planetary radius, used for surface boundary conditions.
g::prec : Gravity at the surface, used for surface boundary conditions.
n::Int : Tidal degree.
U::Int : Tidal potential at the surface.
U_prime::Int : Radial derivative of the tidal potential at the surface.
tau::Int : Tangential tidal stress at the surface.
P::Int : Surface mass load at the surface.

Keyword Arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.
Y::Vector{Int}=[1,2,3,4,5,6,7,8] : Indices for the state variables (default is for standard case).

Returns

B::Array{precc,2} : 3x6 matrix representing the linear relationship between the state variables at the surface and the boundary conditions.
b::Vector{precc} : Vector of length 6 representing the inhomogeneous part of the surface boundary conditions.

source

Obliqua.solid1d_mush_relax.interface_mush_solid — Method

interface_mush_solid(R, Cn_l, Dnp_l, ids, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n; G0=1, Y=[1,2,3,4,5,6,7,8])

Perform the forward-backward relaxation step at the interface between the mushy layer and the solid layer. This function implements the recursion described in N. Kobayashi (2007) for the transition from the 8x8 system to the 6x6 system.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
Cn_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Cn matrix from the previous step.
Dnp_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Dnp matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.

keyword arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.
Y::Vector{Int}=1:8 : Vector of column indices corresponding to the 6x6 system variables in the 8x8 system. This allows for

Returns

Cn_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Cn matrix for the next iteration.
Dnp_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Dnp matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.interface_solid_mush — Method

interface_solid_mush(R, Cn_l, Dnp_l, ids, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n; G0=1, Y=[1,2,3,4,5,6,7,8])

Perform the forward-backward relaxation step at the interface between the solid layer and the mushy layer. This function implements the recursion described in N. Kobayashi (2007) for the transition from the 6x6 system to the 8x8 system.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
Cn_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Cn matrix from the previous step.
Dnp_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Dnp matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.
n::Int : Tidal degree.

keyword arguments

G0::prec=1 : Gravitational constant scale for non-dimensionalization.
Y::Vector{Int}=1:8 : Vector of column indices corresponding to the 6x6 system variables in the 8x8 system. This allows for

Returns

Cn_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Cn matrix for the next iteration.
Dnp_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Dnp matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.propagate_mush — Method

propagate_mush(R, Cn_l, Dnp_l, ids, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n; G0=1, Y=[1,2,3,4,5,6,7,8])

Perform the forward-backward relaxation step for the mushy layer propagation segment. This function implements the recursion described in N. Kobayashi (2007) for the segment of the radial grid that corresponds to the mushy layer, where we propagate the solution up to the surface using the full 8x8 system of equations that includes the porous layer components.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
Cn_l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the Cn matrix from the previous step.
Dnp_l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the Dnp matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6,7,8] : Indices for the state variables (default is for standard case).

Returns

Cn_l::Matrix{precc} : Updated 4x8 matrix representing the "stored" lower half of the Cn matrix for the next iteration.
Dnp_l::Matrix{precc} : Updated 4x8 matrix representing the "stored" lower half of the Dnp matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.propagate_solid — Method

propagate_solid(R, Cn_l, Dnp_l, ids, r, ρ, g, μ, K, ω, n; G0=1, Y=[1,2,3,4,5,6])

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
Cn_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Cn matrix from the previous step.
Dnp_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the Dnp matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6] : Ordering of the solution vector components.

Returns

Cn_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Cn matrix for the next iteration.
Dnp_l::Matrix{precc} : Updated 3x6 matrix representing the "stored" lower half of the Dnp matrix for the next iteration.

source

Obliqua.solid1d_mush_relax.resample_profiles — Method

resample_profiles(radius, rho, visc, shear, bulk_s, bulk_l, bulk_d, alpha, visc_l, phi, k, m_core, dr_min, dr_max)

Arguments

radius::Vector{prec} : Original radius profile (layer boundaries).
rho::Vector{prec} : Original density profile (defined at layer centers).
visc::Vector{prec} : Original viscosity profile (defined at layer centers).
shear::Vector{precc} : Original shear modulus profile (defined at layer centers).
bulk_s::Vector{prec} : Original solid bulk modulus profile (defined at layer centers).
bulk_l::Vector{prec} : Original liquid bulk modulus profile (defined at layer centers).
bulk_d::Vector{prec} : Original deep bulk modulus profile (defined at layer centers).
alpha::Vector{prec} : Original alpha profile (defined at layer centers).
visc_l::Vector{prec} : Original liquid viscosity profile (defined at layer centers).
phi::Vector{prec} : Original phi profile (defined at layer centers).
m_core::prec : Mass of the core, used for gravity calculations.
Δr_min::Int64 : Minimum grid spacing for the new grid.
Δr_max::Int64 : Maximum grid spacing for the new grid.

Returns

Tuple of resampled profiles on the new grid:

r_new_b::Vector{prec} : New radius profile at layer boundaries.
ρ_new::Vector{prec} : New density profile at layer centers.
η_new::Vector{prec} : New viscosity profile at layer centers.
μ_new::Vector{precc} : New shear modulus profile at layer centers.
κ_new::Vector{prec} : New bulk modulus profile at layer centers.
φ_new::Vector{prec} : New phi profile at layer centers.
g_new::Vector{prec} : New gravity profile at layer centers.

source

Obliqua.solid1d_mush_relax.solve_radial_system — Method

solve_radial_system(r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, n, ρ_core, μ_core, κ_core, M_tot; core="liquid")

Arguments

r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{precc} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
ρₗ::Vector{prec} : Vector of liquid densities at the layer centers.
Kl::Vector{prec} : Vector of liquid bulk moduli at the layer centers.
Kd::Vector{prec} : Vector of drained bulk moduli at the layer centers.
α::Vector{prec} : Vector of Biot coefficients at the layer centers.
ηₗ::Vector{prec} : Vector of liquid viscosities at the layer centers.
ϕ::Vector{prec} : Vector of porosities at the layer centers.
k::Vector{prec} : Vector of permeabilities at the layer centers.
n::Int : Tidal degree.
ρ_core::prec : Density of the core, used for core boundary conditions.
μ_core::prec : Shear modulus of the core, used for core boundary conditions.
κ_core::prec : Bulk modulus of the core, used for core boundary conditions.
M_tot::prec : Total mass of the body, used for non-dimensionalization.

Keyword Arguments

`core::String="liquid" : Type of core boundary condition to apply. Options are "liquid" for a fluid core, "solid" for a solid core, and "inertial" for a core with inertial response.

Returns

y_t::Vector{precc} : Vector of length 8 representing the tidal solution at the top of the mantle. This includes the displacements, stresses, and potential at the surface.
y_l::Vector{precc} : Vector of length 8 representing the load solution at the top of the mantle. This includes the displacements, stresses, and potential at the surface.
R::Vector{Matrix{precc}} : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
Y_inv::Vector{Int} : Vector of inverse ordering indices for the 8x8 case.
transitions::Vector{Int} : Indices of the interface layers (the ones closer to the core and surface).

source

Obliqua.solid1d_mush_relax.surface_boundary — Method

surface_boundary(R, CNm_l, DN_l, ids, r, ρ, g, μ, K, ω, n; G0=1, Y=[1,2,3,4,5,6])

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
CNm_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the CNm matrix from the previous step.
DN_l::Matrix{precc} : 3x6 matrix representing the "stored" lower half of the DN matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6,7,8]: Ordering of the solution vector components.

Returns

y_t::Matrix{precc} : 6x1 matrix representing the solution at the surface for the tidal problem.
y_l::Matrix{precc} : 6x1 matrix representing the solution at the surface for the load problem.

source

Obliqua.solid1d_mush_relax.surface_boundary_mush — Method

surface_boundary_mush(R, CNm_l, DN_l, ids, r, ρ, g, μ, K, ω, n; G0=1, Y=[1,2,3,4,5,6,7,8])

Perform the forward-backward relaxation step at the surface boundary for the two-phase problem. This function implements the recursion described in N. Kobayashi (2007) for the final step of the relaxation scheme, where we apply the surface boundary condition and solve for the final solution at the surface, using the full 8x8 system of equations that includes the porous layer components.

Arguments

R::Vector : Vector of 8x8 matrices representing the coefficients of the ODE system at each radial layer.
CNm_l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the CNm matrix from the previous step.
DN_l::Matrix{precc} : 4x8 matrix representing the "stored" lower half of the DN matrix from the previous step.
ids::Tuple{Int, Int} : Tuple containing the start and end indices of the current segment in the radial grid.
r::Vector{prec} : Vector of radial grid points (layer centers).
ρ::Vector{prec} : Vector of densities at the layer centers.
g::Vector{prec} : Vector of gravity values at the layer centers.
μ::Vector{prec} : Vector of shear moduli at the layer centers.
K::Vector{prec} : Vector of bulk moduli at the layer centers.
ω::prec : Angular frequency of the tidal forcing.
n::Int : Tidal degree.

Keyword Arguments

G0::prec=1 : Gravitational constant used for non-dimensional scaling.
Y::Vector{Int}=[1,2,3,4,5,6,7,8]: Ordering of the solution vector components.

Returns

y_t::Matrix{precc} : 8x1 matrix representing the solution at the surface for the tidal problem.
y_l::Matrix{precc} : 8x1 matrix representing the solution at the surface for the load problem.

source

Obliqua.fluid0d.compute_fluid_lovenumbers — Method

compute_fluid_lovenumbers(omega, R, H_magma, g, ρ_ratio, n, σ_R)

Calculate k2 lovenumbers in fluid.

Arguments

omega::prec : Forcing frequency range.
R::prec : Outer radius of fluid segment in mantle.
H_magma::prec : Height of fluid segment in mantle.
g::prec : Surface gravity at top of fluid segment in mantle.
ρ_ratio::prec : Density contrast between current (fluid) and lower (non-fluid) layer.
n::Int=2 : Power of the radial factor (goes with (r/a)^{n}, since r<<a only n=2 contributes significantly).
sigma_R::Float64=1e-3 : Rayleigh drag coefficient.

Returns

k2_T::precc : Complex Tidal k2 Lovenumber.
k2_L::precc : Complex Load k2 Lovenumber.

source

Obliqua.fluid0d.mean_rho — Method

mean_rho(ρ, r)

Calculate the mean density in segment.

Arguments

ρ::Array{prec,1} : Density profile of the planet.
r::Array{prec,1} : Radial positions of layers, from bottom to top of segment.

Returns

mean_density::prec : Mean density in segment.

source

Obliqua.Hansen.get_hansen — Method

get_hansen(ecc)

Wrapper used in your original code: computes Hansen coefficients for given eccentricity. Relies on global n, m, kmin, kmax.

source

Obliqua.Hansen.get_k_range — Method

get_k_range(e, n, m)

Compute k-range for Hansen coefficients given eccentricity e and tidal degree n and order m.

Arguments

e::Float64 : Eccentricity of the orbit.
n::Int : Tidal degree l.
m::Int : Tidal order m.

Returns

k_min::Int : Minimum k-index for Hansen coefficients.
k_max::Int : Maximum k-index for Hansen coefficients.

source

Obliqua.Hansen.hansen_fft — Method

hansen_fft(n, m, e, kmin, kmax; N=nothing)

Compute Hansen coefficients X_k^{n,m}(e) using FFT on mean anomaly.

Returns a tuple (k, Xkm), where:

k is a vector of integer k-indices
Xkm are the corresponding Hansen coefficients (real part)

source

Obliqua.Hansen.kepler_newton — Method

kepler_newton(M, e)

Solve Kepler’s equation E - e*sin(E) = M for the eccentric anomaly E, using Newton iteration. M may be scalar or array.

source

Obliqua.Hansen.nextpow2_int — Method

nextpow2_int(x)

Return the integer p such that 2^p ≥ x.

source