List of all functions
Here an overview of all functions:
GeoModBox.DataFields — Type
DataFields()Structure containing the primary model fields and auxiliary arrays used by the numerical solvers.
DataFields provides a common container for thermal, mechanical, material, and interpolation-related quantities. By default, all array fields are initialized as 1 × 1 arrays and all scalar diagnostics are set to zero. For an actual model, the required fields should be allocated with dimensions consistent with the chosen finite-difference grid before the corresponding solver is called.
Fields
Thermal fields
- `Q` : Volumetric heat-production field.
- `Hs` : Shear-heating field.
- `T` : Temperature at the cell centers.
- `T0` : Temperature at the previous or reference time level.
- `T_ex` : Extended temperature field including ghost cells.
- `T_ex0` : Extended temperature field at the previous or reference time level.
- `Told_ex` : Previously stored extended temperature field.Viscosity fields
- `ηc` : Viscosity at the cell centers.
- `η_ex`: Extended cell-centered viscosity field including ghost cells.
- `ηv` : Viscosity at the grid vertices.Density and phase fields
- `ρ` : Density at the cell centers.
- `ρ_ex`: Extended cell-centered density field.
- `p` : Phase field at the cell centers.
- `p_ex`: Extended cell-centered phase field.
- `pv` : Phase field at the grid vertices.Heat-capacity field
- `cp` : Specific heat capacity at the cell centers.Velocity and pressure fields
- `vx` : Horizontal velocity on the staggered grid.
- `vy` : Vertical velocity on the staggered grid.
- `Pt` : Pressure at the cell centers.
- `vxc` : Horizontal velocity interpolated to the cell centers.
- `vyc` : Vertical velocity interpolated to the cell centers.
- `vxco`: Previously stored cell-centered horizontal velocity.
- `vyco`: Previously stored cell-centered vertical velocity.
- `vc` : Cell-centered velocity magnitude.Interpolation weights
- `wt` : Generic cell-centered interpolation weights.
- `wte` : Interpolation weights on the extended cell-centered grid.
- `wtv` : Interpolation weights at the grid vertices.Thermal diagnostics
- `ΔTtop` : Temperature difference or temperature gradient evaluated along
the upper boundary.
- `ΔTbot` : Temperature difference or temperature gradient evaluated along
the lower boundary.
- `Tmax` : Maximum temperature.
- `Tmin` : Minimum temperature.
- `Tmean` : Mean temperature.Default values
All matrix fields are initialized as
zeros(1, 1)all vector fields as
zeros(1)and all scalar diagnostics as
0.0Notes
DataFields is designed as a flexible container. Individual examples and solvers may use only a subset of the available fields. The fields required by a specific model should therefore be allocated explicitly with dimensions matching the corresponding staggered-grid locations.
Example
NC = (x = 100, y = 50)
NV = (x = NC.x + 1, y = NC.y + 1)
D = DataFields(
T = zeros(Float64, NC.x, NC.y),
T_ex = zeros(Float64, NC.x + 2, NC.y + 2),
vx = zeros(Float64, NV.x, NV.y + 1),
vy = zeros(Float64, NV.x + 1, NV.y),
Pt = zeros(Float64, NC.x, NC.y),
ηc = ones(Float64, NC.x, NC.y),
ηv = ones(Float64, NV.x, NV.y),
)GeoModBox.Geometry — Type
Geometry()Structure containing the spatial limits of the two-dimensional model domain.
The geometry is defined by the minimum and maximum coordinates in the horizontal and vertical directions. The fields can be specified during construction or modified afterward because Geometry is a mutable structure.
Fields
xmin: Minimum horizontal coordinate.xmax: Maximum horizontal coordinate.ymin: Minimum vertical coordinate.ymax: Maximum vertical coordinate.
Default values
xmin = 0.0
xmax = 1.0
ymin = -1.0
ymax = 0.0Example
M = Geometry()which returns
Geometry(0.0, 1.0, -1.0, 0.0)A custom model geometry can be initialized using keyword arguments:
M = Geometry(
xmin = -50.0e3,
xmax = 50.0e3,
ymin = 0.0,
ymax = 100.0e3,
)GeoModBox.GridSpacing — Type
GridSpacing()Structure containing the grid spacing of a two-dimensional Cartesian finite-difference grid.
The structure stores the horizontal and vertical grid spacing used throughout the numerical discretization. The values are typically computed from the model geometry and the number of grid cells in each coordinate direction.
Fields
- `x`: Horizontal grid spacing, Δx.
- `y`: Vertical grid spacing, Δy.Default values
x = 0.0
y = 0.0Notes
For a regular Cartesian grid, the grid spacing is commonly computed as
\[\Delta x = \frac{x_{\max}-x_{\min}}{N_x},\]
\[\Delta y = \frac{y_{\max}-y_{\min}}{N_y},\]
where Nₓ and Nᵧ denote the number of cells in the horizontal and vertical directions, respectively.
Example
Δ = GridSpacing()
NC = (x = 100, y = 50)
Δ = GridSpacing(
x = (M.xmax - M.xmin) / NC.x,
y = (M.ymax - M.ymin) / NC.y,
)GeoModBox.Physics — Type
Physics()Structure containing the physical parameters required for thermal and thermo-mechanical simulations.
The structure stores the material properties and physical constants used throughout GeoModBox.jl. All quantities are specified in SI units and can be modified during construction using keyword arguments.
Fields
- `g` : Gravitational acceleration [m/s²].
- `ρ₀` : Reference density [kg/m³].
- `k` : Thermal conductivity [W/m/K].
- `cp` : Specific heat capacity [J/kg/K].
- `α` : Thermal expansion coefficient [1/K].
- `Q₀` : Volumetric heat-production rate [W/m³].
- `η₀` : Reference dynamic viscosity [Pa·s].
- `κ` : Thermal diffusivity [m²/s].
- `ΔT` : Reference temperature difference [K].
- `Ttop`: Upper boundary temperature [K].
- `Tbot`: Lower boundary temperature [K].
- `Ra` : Rayleigh number.
- `RG` : Universal gas constant [J/mol/K].Default values
g = 9.81
ρ₀ = 3300.0
k = 4.125
cp = 1250.0
α = 2.0e-5
Q₀ = 0.0
η₀ = 3.947725485e23
κ = k / (ρ₀ * cp)
ΔT = 2500.0
Ttop = 273.15
Tbot = Ttop + ΔT
Ra = -9999.0
RG = 8.314Notes
If Ra is negative (default), the basal Rayleigh number is computed from the remaining physical parameters during model initialization. A positive value may be supplied when prescribing a Rayleigh number directly.
Example
P = Physics()
P = Physics(
η₀ = 1.0e21,
ΔT = 1300.0,
Q₀ = 1.0e-6,
)GeoModBox.TimeParameter — Type
TimeParameter()Structure containing the parameters controlling the temporal discretization of a simulation.
The structure stores the simulation end time, the stability factors used to compute the Courant and diffusion time steps, the current time-step size, and the maximum number of time steps. These parameters are shared by the thermal, advection, and thermo-mechanical solvers.
Fields
- `year` : Number of seconds in one year.
- `tmax` : Maximum simulation time.
- `Δfacc` : Safety factor for the Courant (advection) stability criterion.
- `Δfacd` : Safety factor for the diffusion stability criterion.
- `Δ` : Current time step.
- `Δc` : Time step determined from the Courant stability criterion.
- `Δd` : Time step determined from the diffusion stability criterion.
- `itmax` : Maximum number of iterations (time steps).Default values
year = 365.25 * 24 * 3600
tmax = 1000.0 # [Ma]
Δfacc = 0.9
Δfacd = 0.9
Δ = 0.0
Δc = 0.0
Δd = 0.0
itmax = 8000Notes
The stable time step is typically chosen as
\[\Delta t = \min(\Delta t_c,\Delta t_d),\]
where Δt_c is obtained from the Courant criterion and Δt_d from the diffusion stability criterion. The factors Δfacc and Δfacd provide an additional safety margin.
Although tmax is initialized in millions of years (Ma), it is commonly converted to seconds before the simulation begins and subsequently non-dimensionalized if required.
Example
T = TimeParameter(
tmax = 50.0,
Δfacc = 0.8,
Δfacd = 0.9,
itmax = 10000,
)
T.tmax *= 1e6 * T.yearGeoModBox.Dani_Solution_vec! — Method
Dani_Solution_vec!(type, D, M, x, y, rad, mus_i, NC, NV)Compute the analytical pressure and velocity solution for a circular viscous inclusion embedded in an infinite viscous matrix undergoing either pure shear or simple shear deformation.
The analytical solution is evaluated on the staggered finite-difference grid used by GeoModBox.jl. Pressure is computed at the cell centers, whereas the velocity components are evaluated at the corresponding staggered velocity locations. In addition, the analytical boundary velocities required by the numerical Stokes solver are extracted from the analytical solution.
The implementation is based on the complex-variable formulation of the Eshelby inclusion problem presented by Schmid (2002) and was vectorized by Thibault Duretz and Ludovic Räss.
Arguments
- `type` : Deformation mode:
- `:PureShear`
- `:SimpleShear`
- `D` : Data structure containing the analytical pressure, velocity, and
boundary-condition arrays.
- `M` : Model geometry.
- `x` : Horizontal coordinates of the staggered grid.
- `y` : Vertical coordinates of the staggered grid.
- `rad` : Radius of the circular inclusion.
- `mus_i` : Inclusion viscosity normalized by the matrix viscosity.
- `NC` : Number of cell centers in each coordinate direction.
- `NV` : Number of grid vertices in each coordinate direction.Output
The function updates the analytical solution stored in D:
- `D.Pa` : Pressure at the cell centers.
- `D.Vxa` : Horizontal velocity on the staggered `vₓ` grid.
- `D.Vya` : Vertical velocity on the staggered `vᵧ` grid.
- `D.Vx_W`, `D.Vx_E`: Horizontal boundary velocities.
- `D.Vy_S`, `D.Vy_N`: Vertical boundary velocities.
- `D.Vx_S`, `D.Vx_N`: Horizontal velocities on the non-conforming boundary
nodes.
- `D.Vy_W`, `D.Vy_E`: Vertical velocities on the non-conforming boundary
nodes.The arrays contained in D are modified in place.
Notes
The matrix viscosity is normalized to unity (ηₘ = 1), and the inclusion viscosity is specified through the viscosity ratio mus_i.
This analytical solution is primarily intended for benchmarking the Stokes solver by comparing the numerical and analytical pressure and velocity fields.
Example
Dani_Solution_vec!(
:PureShear,
D,
M,
x,
y,
0.2,
1.0e3,
NC,
NV,
)References
Schmid, D. W. (2002). Finite and infinite heterogeneities under pure and simple shear (Doctoral dissertation, ETH Zurich).
The implementation is based on the original MATLAB routine CYLPMATRIX.m by Dani Schmid and the vectorized Julia version by Thibault Duretz and Ludovic Räss.
GeoModBox.Scaling.Constants — Type
Constants()Structure containing the characteristic scaling quantities used for non-dimensionalization and post-processing.
The structure stores the reference scales for length, time, velocity, stress, temperature, and volumetric heat production. By default, all values are initialized to zero and can be assigned by the user depending on the chosen scaling.
Fields
- `hsc`: Characteristic length scale.
- `tsc`: Characteristic time scale.
- `vsc`: Characteristic velocity scale.
- `τsc`: Characteristic stress scale.
- `Tsc`: Characteristic temperature scale.
- `Qsc`: Characteristic volumetric heat-production scale.Default values
hsc = 0.0
tsc = 0.0
vsc = 0.0
τsc = 0.0
Tsc = 0.0
Qsc = 0.0Example
S = Constants()
S.hsc = 100e3 # m S.tsc = 1.0e6 # yr S.vsc = 0.1 # cm/yr S.τsc = 1.0e8 # Pa S.Tsc = 1300.0 # °C S.Qsc = 1.0e-6 # W/m³
GeoModBox.Scaling.ScaleParameters! — Method
ScaleParameters!(S, M, Δ, T, P, D)Scale the model geometry, time parameters, physical parameters, and field variables using the characteristic scales stored in S.
The function converts the dimensional model setup to its non-dimensional form. The characteristic scales are typically obtained from ScalingConstants!() and are based on the model height and the thermal diffusion time. The scaling is performed in place by modifying the corresponding structures.
The following quantities are non-dimensionalized:
Geometry
\[x=\frac{x}{h_{\mathrm{sc}}}, \qquad y=\frac{y}{h_{\mathrm{sc}}}\]
including the model boundaries and the grid spacing.
Time
\[t=\frac{t}{t_{\mathrm{sc}}}\]
including the maximum simulation time and all time-step estimates.
Temperature
The boundary temperatures are converted from Kelvin to non-dimensional temperature using
\[T=\frac{T-273.15}{T_{\mathrm{sc}}}.\]
Volumetric heat production
\[Q=\frac{Q}{Q_{\mathrm{sc}}}.\]
The reference viscosity scaling
\[\eta=\frac{\eta}{\eta_0}\]
can optionally be applied if viscosity is stored in dimensional units.
Arguments
- `S`: Structure containing the characteristic scaling constants.
- `M`: Model geometry.
- `Δ`: Grid spacing.
- `T`: Time-parameter structure.
- `P`: Physical-parameter structure.
- `D`: Structure containing the model fields.Notes
The function modifies all supplied structures in place. After calling ScaleParameters!, the governing equations can be solved in their non-dimensional form.
Typically, the scaling workflow is
S = ScalingConstants!(M, P)
ScaleParameters!(S, M, Δ, T, P, D)Example
M = Geometry(0, 1000, -1000, 0)
P = Physics()
S = ScalingConstants!(M, P)
ScaleParameters!(S, M, Δ, T, P, D)GeoModBox.Scaling.ScalingConstants! — Method
ScalingConstants!(M, P)Compute the characteristic scaling constants used for non-dimensionalization.
The scaling is based on the model height and the thermal diffusion time, which are commonly used as characteristic scales in mantle convection problems. The resulting structure contains the reference scales for length, time, velocity, stress, temperature, and volumetric heat production.
The scaling constants are defined as
\[h_{\mathrm{sc}} = H, t_{\mathrm{sc}} = \frac{H^2}{\kappa}, v_{\mathrm{sc}} = \frac{\kappa}{H}, au_{\mathrm{sc}} = \frac{\eta_0 \kappa}{H^2}, T_{\mathrm{sc}} = \Delta T, Q_{\mathrm{sc}} = \frac{\Delta T \kappa \rho_0 c_p}{H^2},\]
where - H = ymax - ymin is the model height, - κ is the thermal diffusivity, - η₀ is the reference viscosity, - ΔT is the reference temperature difference, - ρ₀ is the reference density, and - cₚ is the specific heat capacity.
Arguments
- M: Model geometry.
- P: Structure containing the physical parameters.Returns
Returns a Constants structure containing the characteristic scaling constants in SI units.
Example
M = Geometry(0, 1000, -1000, 0) P = Physics() S = ScalingConstants!(M, P)
S.vsc
Notes
The returned scaling constants can be used to convert dimensional model parameters to their non-dimensional counterparts, or to convert non-dimensional simulation results back to SI units for analysis and visualization.
GeoModBox.HeatEquation.OneD.AssembleMatrix1D — Method
AssembleMatrix1D(ρ, cp, k, Δx, Δt, nc, BC, K; C=0.0)Assembles the coefficient matrix for the one-dimensional transient heat equation with variable thermal properties.
Temperature is defined at the cell centroids, while thermal conductivity is defined at the cell boundaries. The conductive term is discretized in conservative flux-divergence form using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the matrix coefficients at the western and eastern boundaries.
The temporal discretization is controlled by C. For C < 1, the matrix contains the implicit contribution of the conductive term. For C = 1, the conductive contribution vanishes and the matrix contains only the transient storage term.
Arguments
ρ : Density defined at the cell centroids.
cp : Specific heat capacity defined at the cell centroids.
k : Thermal conductivity defined at the cell boundaries.
Δx : Grid spacing.
Δt : Time step.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
K : Coefficient matrix to be assembled.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.Notes
The matrix corresponds to the current-time contribution of the generalized time-discretization scheme used by ComputeResiduals1D!.
The assembled matrix is tridiagonal because each centroid temperature is coupled only to its western and eastern neighbours. Boundary conditions modify the coefficients in the first and last matrix rows.
The matrix is modified in place and finalized using flush!(K).
GeoModBox.HeatEquation.OneD.AssembleMatrix1Dc! — Method
AssembleMatrix1Dc!(κ, Δx, Δt, nc, BC, K; C=0.0)Assembles the coefficient matrix for the one-dimensional transient heat equation with constant thermal diffusivity.
Temperature is defined at cell centroids, and the spatial diffusion term is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the matrix coefficients at the western and eastern boundaries.
The temporal discretization is controlled by C. For C < 1, the matrix contains the implicit contribution of the diffusion term. For C = 1, the diffusion contribution vanishes and the matrix reduces to the temporal term.
Arguments
κ : Thermal diffusivity.
Δx : Grid spacing.
Δt : Time step.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
K : Coefficient matrix to be assembled in place.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.Notes
The matrix corresponds to the current-time contribution of the generalized time-discretization scheme used by ComputeResiduals1Dc!. Backward Euler and Crank–Nicolson require a diffusion contribution in the matrix, whereas Forward Euler is fully explicit.
The assembled matrix is tridiagonal. Boundary conditions modify the diagonal coefficients of the first and last equations. After assembly, flush!(K) is called to finalize the sparse matrix.
GeoModBox.HeatEquation.OneD.BackwardEuler1Dc! — Method
BackwardEuler1Dc!(
implicit, κ, Δx, Δt, nc, BC, K, rhs;
Q=0.0, ρ=3200.0, cp=1200.0
)Solves the one-dimensional transient heat equation using an implicit Backward Euler finite-difference scheme with constant thermal diffusivity.
Temperature is defined at cell centroids, while the diffusion term is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrix and right-hand-side vector.
Optional internal heating can be included through the volumetric heat-production term Q.
Arguments
implicit : Structure or tuple containing `T`, the temperature array
defined on the centroids. The array is updated in place with
the solution at the new time level.
κ : Thermal diffusivity.
Δx : Grid spacing.
Δt : Time step.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
K : Coefficient matrix for the linear system of equations.
rhs : Right-hand-side vector.Keyword Arguments
Q : Volumetric heat-production rate (default: `0.0`).
ρ : Density (default: `3200.0`).
cp : Specific heat capacity (default: `1200.0`).Notes
Backward Euler is first-order accurate in time and unconditionally stable for the linear diffusion equation. The spatial discretization is second-order accurate.
The coefficient matrix and right-hand-side vector are assembled inside the routine. The temperature at the new time level is obtained by solving the resulting linear system, and implicit.T is updated in place.
GeoModBox.HeatEquation.OneD.CNA1Dc! — Method
CNA1Dc!(
cna, κ, Δx, Δt, nc, BC, K1, K2;
Q=0.0, ρ=3200.0, cp=1200.0
)Solves the one-dimensional transient heat equation using the Crank–Nicolson finite-difference scheme with constant thermal diffusivity.
Temperature is defined at cell centroids, while the diffusion term is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrices and right-hand-side vector.
The Crank–Nicolson method combines the diffusion operators evaluated at the current and new time levels with equal weighting, resulting in a second-order accurate temporal discretization. Optional internal heating can be included through the volumetric heat-production term Q.
Arguments
cna : Structure or tuple containing `T`, the temperature array
defined on the centroids. The array is updated in place with
the solution at the new time level.
κ : Thermal diffusivity.
Δx : Grid spacing.
Δt : Time step.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
K1 : Coefficient matrix associated with the unknown temperature
field at the new time level.
K2 : Coefficient matrix associated with the known temperature
field at the previous time level.Keyword Arguments
Q : Volumetric heat-production rate (default: `0.0`).
ρ : Density (default: `3200.0`).
cp : Specific heat capacity (default: `1200.0`).Notes
The Crank–Nicolson scheme is second-order accurate in both space and time. For linear diffusion problems it is unconditionally stable, although large time steps may produce weak temporal oscillations.
The coefficient matrices K1 and K2 are assembled inside the routine. The right-hand-side vector is constructed from the previous temperature field and the volumetric heat-production term before solving the resulting linear system. The computed temperature field overwrites cna.T.
GeoModBox.HeatEquation.OneD.ComputeResiduals1D! — Method
ComputeResiduals1D!(
R, T, T_ex, T0, T_ex0, Q, ∂T, q, ρ, Cp, k, BC, Δx, Δt;
C=0.0
)Computes the residual of the one-dimensional transient heat equation with variable thermal properties and volumetric heat production.
The conductive heat flux is discretized in conservative form, with temperature defined at the cell centroids and thermal conductivity at the cell boundaries. Ghost nodes are used to impose Dirichlet and Neumann boundary conditions.
The routine evaluates the residual for the generalized θ-method, allowing Backward Euler, Crank–Nicolson, and Forward Euler time discretizations through the weighting parameter C. The residual is intended for defect-correction iterations and is used together with the coefficient matrix assembled by AssembleMatrix1D!.
Arguments
R : Residual vector defined at the cell centroids.
T : Temperature at the current iteration or new time level.
T_ex : Extended current-temperature array including ghost nodes.
T0 : Temperature at the previous time level.
T_ex0 : Extended previous-temperature array including ghost nodes.
Q : Volumetric heat-production rate defined at the cell centroids.
∂T : Structure or tuple containing the temperature gradients
`∂x` and `∂x0`.
q : Structure or tuple containing the conductive heat fluxes
`x` and `x0`.
ρ : Density defined at the cell centroids.
Cp : Specific heat capacity defined at the cell centroids.
k : Thermal conductivity defined at the cell boundaries.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
Δx : Grid spacing.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.Notes
The residual is evaluated in conservative form as
- the transient storage term,
- the divergence of the conductive heat flux, and
- the volumetric heat-production term.
Backward Euler evaluates the conductive heat flux entirely at the new time level, Crank–Nicolson averages the previous and current time levels, and Forward Euler uses only the previous time level.
The routine updates the ghost-node values, temperature gradients, and heat fluxes before evaluating the residual.
GeoModBox.HeatEquation.OneD.ComputeResiduals1Dc! — Method
ComputeResiduals1Dc!(
R, T, T_ex, T0, T_ex0, ∂2T, κ, BC, Δx, Δt;
C=0.0, Q=0.0, ρ=3200.0, cp=1200.0
)Computes the residual of the one-dimensional transient heat equation for constant thermal parameters.
Temperature is defined at cell centroids, while second spatial derivatives are evaluated using central finite differences and ghost nodes. Dirichlet and Neumann boundary conditions are supported at the western and eastern boundaries.
The temporal discretization is controlled by C and can represent Backward Euler, Crank–Nicolson, or Forward Euler time integration. Optional internal heating can be included through the volumetric heat-production term Q.
Arguments
R : Residual vector defined on the centroids.
T : Temperature at the current iteration or new time level.
T_ex : Extended current-temperature array including ghost nodes.
T0 : Temperature at the previous time level.
T_ex0 : Extended previous-temperature array including ghost nodes.
∂2T : Structure or tuple containing the second spatial derivatives
`∂x2` and `∂x20`.
κ : Thermal diffusivity.
BC : Structure or tuple defining the boundary-condition types and
values at the western and eastern boundaries.
Δx : Grid spacing.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.
Q : Volumetric heat-production rate (default: `0.0`).
ρ : Density (default: `3200.0`).
cp : Specific heat capacity (default: `1200.0`).Notes
The residual is evaluated as a weighted combination of the diffusion term at the previous and current time levels. Backward Euler and Crank–Nicolson require the current-temperature contribution, whereas Forward Euler uses only the previous temperature field.
Backward Euler is first-order accurate in time, Crank–Nicolson is second-order accurate in time, and Forward Euler is first-order accurate in time. The spatial discretization is second-order accurate.
GeoModBox.HeatEquation.OneD.ForwardEuler1D! — Method
ForwardEuler1D!(T, Py, Δt, Δy, nc, BC)Solves the one-dimensional transient heat equation using an explicit Forward Euler finite-difference scheme with variable thermal parameters and internal heat production.
Temperature is defined at cell centroids, while thermal conductivity is defined at the cell boundaries. The conductive heat flux is therefore discretized in conservative form. Density, specific heat capacity, and heat production are defined at the cell centroids.
Dirichlet and Neumann boundary conditions are supported at the southern and northern boundaries. Ghost nodes are used to impose the boundary conditions.
Arguments
T : Structure or tuple containing:
`T`, the temperature array defined at the cell centroids, and
`T_ex`, the extended temperature array including ghost nodes.
Both arrays are updated in place.
Py : Structure or tuple containing the thermal parameters:
`k`, thermal conductivity defined at the cell boundaries;
`ρ`, density defined at the cell centroids;
`cp`, specific heat capacity defined at the cell centroids;
`H`, specific heat-production rate defined at the cell
centroids.
Δt : Time step.
Δy : Grid spacing.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types and
values at the southern and northern boundaries.Notes
Scalar values may be supplied for k, ρ, cp, and H. In this case, the values are expanded internally to arrays of the required size.
The heat-production parameter H is expressed per unit mass. The corresponding volumetric heat-production rate is therefore
Q = ρ * H.The method is first-order accurate in time and second-order accurate in space. Because the temperature update is explicit, the timestep must satisfy the diffusive stability condition. For variable thermal parameters, the timestep should be selected using the largest local thermal diffusivity.
After each timestep, the updated centroid temperatures are copied to the interior of T.T_ex. The ghost-node values are updated again when the routine is called for the next timestep.
GeoModBox.HeatEquation.OneD.ForwardEuler1Dc! — Method
ForwardEuler1Dc!(explicit, κ, Δx, Δt, nc, BC;
Q=zeros(nc), ρ=3200.0, cp=1200.0)Solves the one-dimensional transient heat equation using an explicit Forward Euler finite-difference scheme.
Temperature is defined at cell centroids, while the heat flux is evaluated at the cell boundaries. Ghost nodes are used to impose Dirichlet or Neumann boundary conditions. Thermal diffusivity is assumed to be constant.
Optional internal heating can be included through the volumetric heat production term Q.
Arguments
explicit : Structure or tuple containing:
`T`, the temperature array on the centroids, and
`T_ex`, the extended temperature array including ghost nodes.
κ : Thermal diffusivity.
Δx : Grid spacing.
Δt : Time step.
nc : Number of centroid nodes.
BC : Structure or tuple defining the boundary-condition types
and values at the western and eastern boundaries.Keyword Arguments
Q : Volumetric heat production rate defined on the centroids
(default: `zeros(nc)`).
ρ : Density (default: `3200.0`).
cp : Specific heat capacity (default: `1200.0`).Notes
The temperature update is explicit and second-order accurate in space but first-order accurate in time. The timestep must therefore satisfy the diffusive stability condition.
After each timestep, the updated centroid temperatures are copied back to the interior of explicit.T_ex. The ghost-node values are updated again when the routine is called for the next timestep.
GeoModBox.HeatEquation.TwoD.ADI2Dc! — Method
ADI2Dc!(
T, κ, Δx, Δy, Δt, NC, BC;
Q=zeros(NC...), ρ=3300.0, cp=1200.0
)Solves the two-dimensional transient heat equation using a second-order Alternating Direction Implicit (ADI) finite-difference scheme with constant thermal diffusivity.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the linear systems at the western, eastern, southern, and northern boundaries.
Optional volumetric heat production may be included.
Arguments
T : Structure or tuple containing:
`T`, the temperature field defined at the cell centroids, and
`T_ex`, the extended temperature field including ghost nodes.
Both arrays are updated in place.
κ : Thermal diffusivity.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.
Δt : Time step.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.Keyword Arguments
Q : Volumetric heat-production rate (default: `zeros(NC...)`).
ρ : Density (default: `3300.0`).
cp : Specific heat capacity (default: `1200.0`).Notes
The ADI method advances the solution over one time step by performing two successive implicit half-steps:
- An implicit solve in the vertical direction while treating the horizontal diffusion term explicitly.
- An implicit solve in the horizontal direction while treating the vertical diffusion term explicitly.
The symmetric combination of the two directional half-steps yields second-order accuracy in time. The spatial diffusion operators are discretized using second-order central finite differences, making the method second-order accurate in both space and time.
For the linear two-dimensional diffusion equation, the ADI scheme is unconditionally stable. Each half-step requires the solution of a system that is implicit in only one spatial direction, rather than a fully coupled two-dimensional system.
The temperature field is updated in place, and the interior of T.T_ex is synchronized with the new temperature field after the second half-step.
GeoModBox.HeatEquation.TwoD.AnalyticalSolution2D! — Method
AnalyticalSolution!(Te, x, y, t)Calls 2D analytical solution for a 2D diffusion problem:
Te : 2D matrix
x : x coordinate array
y : y coordinate array
t : timeExamples
julia> a = zeros(2,2); x = [0, 1]; y = [0, 1];
julia> AnalyticalSolution!(a, x, y, 0.0)
2×2 Matrix{Float64}:
100.0 6.04202e-67
6.04202e-67 3.6506e-135GeoModBox.HeatEquation.TwoD.AssembleMatrix2D — Method
AssembleMatrix2D(ρ, cp, k, BC, Num, nc, Δ, Δt; C=0.0)Assembles the coefficient matrix for the conservative two-dimensional transient heat equation with spatially variable thermal properties.
Temperature is defined at the cell centroids, while thermal conductivity is defined at the cell faces. The conductive term is discretized in flux form using second-order finite differences, resulting in a conservative discretization of the variable-conductivity diffusion operator.
Dirichlet and Neumann boundary conditions are incorporated directly into the matrix coefficients at the western, eastern, southern, and northern boundaries.
Arguments
ρ : Density defined at the cell centroids.
cp : Specific heat capacity defined at the cell centroids.
k : Structure or tuple containing the face-centered thermal
conductivities `x` and `y`.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
Num : Structure or tuple containing the global numbering of the
centroid nodes.
nc : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
Δ : Structure or tuple containing the horizontal and vertical grid
spacings `x` and `y`.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.Notes
The matrix corresponds to the current-time contribution of the generalized time-discretization scheme used by ComputeResiduals2D!.
The transient storage term is weighted by the local volumetric heat capacity
ρ * cp,while the conductive coefficients are determined from the thermal conductivities at the western, eastern, southern, and northern cell faces.
The spatial operator is represented by a five-point finite-difference stencil. For an interior cell, the corresponding matrix row couples the central temperature to its four neighboring centroid temperatures. This results in a sparse matrix with five non-zero diagonals for the adopted global numbering.
For C = 1, the implicit conductive contribution vanishes, and the matrix contains only the transient storage term. Volumetric heat sources do not enter the coefficient matrix and are instead included in the residual or right-hand-side vector.
The matrix is assembled internally as an ExtendableSparseMatrix and finalized using flush!(K) before being returned.
GeoModBox.HeatEquation.TwoD.AssembleMatrix2Dc — Method
AssembleMatrix2Dc(κ, BC, Num, nc, Δ, Δt; C=0.0)Assembles the coefficient matrix for the two-dimensional transient heat equation with constant thermal diffusivity.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the matrix coefficients at the western, eastern, southern, and northern boundaries.
The temporal discretization is controlled by C. For C < 1, the matrix contains the implicit contribution of the diffusion operator. For C = 1, the diffusion contribution vanishes and the matrix contains only the transient storage term.
Arguments
κ : Thermal diffusivity.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
Num : Structure or tuple containing the global numbering of the
centroid nodes.
nc : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
Δ : Structure or tuple containing the horizontal and vertical grid
spacings `x` and `y`.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.Notes
The matrix corresponds to the current-time contribution of the generalized time-discretization scheme used by ComputeResiduals2Dc!.
The diffusion operator is discretized using the classical five-point finite difference stencil, resulting in a sparse matrix with five non-zero diagonals. Boundary conditions modify the coefficients of the equations adjacent to the domain boundaries.
The matrix is assembled internally as an ExtendableSparseMatrix and finalized using flush!(K) before being returned.
GeoModBox.HeatEquation.TwoD.BackwardEuler2Dc! — Method
BackwardEuler2Dc!(
D, κ, Δx, Δy, Δt, NC, BC, rhs, K, Num;
Q=zeros(NC...), ρ=3300.0, cp=1200.0, Qₛ=zeros(NC...)
)Solves the two-dimensional transient heat equation using an implicit Backward Euler finite-difference scheme with constant thermal diffusivity.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrix and right-hand-side vector at the western, eastern, southern, and northern boundaries.
Optional volumetric heat-production and shear-heating terms may be included. The resulting linear system is solved using a direct left-matrix division.
Arguments
D : Structure or tuple containing:
`T`, the temperature field defined at the cell centroids, and
`T_ex`, the extended temperature field including ghost nodes.
Both arrays are updated in place.
κ : Thermal diffusivity.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.
Δt : Time step.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
rhs : Right-hand-side vector.
K : Sparse coefficient matrix for the linear system.
Num : Structure or tuple containing the global numbering of the
centroid nodes.Keyword Arguments
Q : Volumetric heat-production rate (default: `zeros(NC...)`).
ρ : Density (default: `3300.0`).
cp : Specific heat capacity (default: `1200.0`).
Qₛ : Volumetric shear-heating rate (default: `zeros(NC...)`).Notes
Backward Euler is first-order accurate in time and unconditionally stable for the linear diffusion equation. The spatial discretization is second-order accurate.
The diffusion operator is represented by a five-point finite-difference stencil. Boundary conditions modify both the matrix coefficients and the right-hand-side vector for the equations adjacent to the domain boundaries.
The coefficient matrix and right-hand-side vector are assembled inside the routine. The temperature at the new time level is obtained from
K * Tⁿ⁺¹ = rhsand written to D.T. The interior of D.T_ex is then updated with the new temperature field.
GeoModBox.HeatEquation.TwoD.CNA2Dc! — Method
CNA2Dc!(
D, κ, Δx, Δy, Δt, NC, BC, rhs, K1, K2, Num;
Q=zeros(NC...), ρ=3300.0, cp=1200.0, Qₛ=zeros(NC...)
)Solves the two-dimensional transient heat equation using the Crank–Nicolson finite-difference scheme with constant thermal diffusivity.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrices and right-hand-side vector at the western, eastern, southern, and northern boundaries.
Optional volumetric heat-production and shear-heating terms may be included. The resulting linear system is solved using a direct left-matrix division.
Arguments
D : Structure or tuple containing:
`T`, the temperature field defined at the cell centroids, and
`T_ex`, the extended temperature field including ghost nodes.
Both arrays are updated in place.
κ : Thermal diffusivity.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.
Δt : Time step.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
rhs : Right-hand-side vector.
K1 : Sparse coefficient matrix associated with the unknown
temperature field at the new time level.
K2 : Sparse coefficient matrix associated with the known
temperature field at the previous time level.
Num : Structure or tuple containing the global numbering of the
centroid nodes.Keyword Arguments
Q : Volumetric heat-production rate (default: `zeros(NC...)`).
ρ : Density (default: `3300.0`).
cp : Specific heat capacity (default: `1200.0`).
Qₛ : Volumetric shear-heating rate (default: `zeros(NC...)`).Notes
The Crank–Nicolson scheme is second-order accurate in both space and time. For linear diffusion problems it is unconditionally stable, although large time steps may introduce weak temporal oscillations.
The diffusion operator is represented by a five-point finite-difference stencil. Boundary conditions modify both the coefficient matrices and the right-hand-side vector for the equations adjacent to the domain boundaries.
The matrices K1 and K2 represent the implicit and explicit contributions of the generalized time discretization, respectively. The right-hand-side vector is assembled from the previous temperature field and the volumetric heat-source terms before solving
K1 * Tⁿ⁺¹ = rhsThe computed temperature field is written to D.T, and the interior of D.T_ex is updated accordingly.
GeoModBox.HeatEquation.TwoD.ComputeResiduals2D! — Method
ComputeResiduals2D!(
R, T, T_ex, T0, T_ex0, Q, ∂T, q, ρ, Cp, k, BC, Δ, Δt;
C=0.0, Qₛ=0.0
)Computes the residual of the conservative two-dimensional transient heat equation with variable thermal conductivity and optional volumetric heat sources.
Temperature is defined at the cell centroids, while heat fluxes are evaluated at the cell faces using Fourier's law. The divergence of the heat flux is then computed using second-order finite differences, providing a conservative discretization of the diffusion operator.
Ghost nodes are used to impose Dirichlet and Neumann boundary conditions at the western, eastern, southern, and northern boundaries.
The temporal discretization is controlled by C and can represent Backward Euler, Crank–Nicolson, or Forward Euler time integration. The residual is intended for defect-correction iterations and is used together with the coefficient matrix assembled by AssembleMatrix2D.
Arguments
R : Residual field defined at the cell centroids.
T : Temperature field at the current iteration or new time level.
T_ex : Extended current-temperature field including ghost nodes.
T0 : Temperature field at the previous time level.
T_ex0 : Extended previous-temperature field including ghost nodes.
Q : Volumetric heat-production rate.
∂T : Structure or tuple containing the temperature gradients
`∂x` and `∂y`.
q : Structure or tuple containing the heat fluxes
(`x`, `y`, `x0`, and `y0`).
ρ : Density.
Cp : Specific heat capacity.
k : Thermal conductivity defined on the cell faces.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
Δ : Structure or tuple containing the horizontal and vertical grid
spacings `x` and `y`.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.
Qₛ : Volumetric shear-heating rate (default: `0.0`).Notes
The residual is evaluated from
- the transient heat-storage term,
- the divergence of the conductive heat flux,
- the volumetric heat-production term, and
- the volumetric shear-heating term.
The heat flux is computed from Fourier's law,
q = -k ∇T,using thermal conductivities defined at the cell faces. The divergence of the heat flux is then evaluated at the cell centroids, resulting in a conservative finite-difference discretization that naturally accommodates spatially varying thermal conductivity.
Backward Euler evaluates the conductive fluxes entirely at the new time level, Crank–Nicolson averages the previous and current time levels, and Forward Euler evaluates the conductive fluxes at the previous time level.
GeoModBox.HeatEquation.TwoD.ComputeResiduals2Dc! — Method
ComputeResiduals2Dc!(
R, T, T_ex, T0, T_ex0, ∂2T, κ, BC, Δ, Δt;
C=0.0, Q=0.0, ρ=3300.0, cp=1200.0, Qₛ=0.0
)Computes the residual of the two-dimensional transient heat equation with constant thermal diffusivity and optional volumetric heat sources.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Ghost nodes are used to impose Dirichlet and Neumann boundary conditions at the western, eastern, southern, and northern boundaries.
The temporal discretization is controlled by C and can represent Backward Euler, Crank–Nicolson, or Forward Euler time integration. The residual is intended for defect-correction iterations and is used together with the coefficient matrix assembled by AssembleMatrix2Dc.
Arguments
R : Residual field defined at the cell centroids.
T : Temperature field at the current iteration or new time level.
T_ex : Extended current-temperature field including ghost nodes.
T0 : Temperature field at the previous time level.
T_ex0 : Extended previous-temperature field including ghost nodes.
∂2T : Structure or tuple containing the second spatial derivatives
`∂x2`, `∂y2`, `∂x20`, and `∂y20`.
κ : Thermal diffusivity.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
Δ : Structure or tuple containing the horizontal and vertical grid
spacings `x` and `y`.
Δt : Time step.Keyword Arguments
C : Temporal weighting parameter (default: `0.0`):
`0.0` for Backward Euler,
`0.5` for Crank–Nicolson,
`1.0` for Forward Euler.
Q : Volumetric heat-production rate (default: `0.0`).
ρ : Density (default: `3300.0`).
cp : Specific heat capacity (default: `1200.0`).
Qₛ : Volumetric shear-heating rate (default: `0.0`).Notes
The residual is evaluated as the sum of
- the transient temperature-change term,
- the weighted two-dimensional diffusion term, and
- the volumetric heat-source terms.
Backward Euler evaluates the diffusion operator entirely at the new time level, Crank–Nicolson averages the previous and current time levels, and Forward Euler uses only the previous time level.
The routine updates the ghost-node values and second spatial derivatives before evaluating the residual.
GeoModBox.HeatEquation.TwoD.ForwardEuler2Dc! — Method
ForwardEuler2Dc!(
D, κ, Δx, Δy, Δt, NC, BC;
Q=zeros(NC...), ρ=3300.0, cp=1200.0, Qₛ=zeros(NC...)
)Solves the two-dimensional transient heat equation using an explicit Forward Euler finite-difference scheme with constant thermal diffusivity.
Temperature is defined at the cell centroids, while the diffusion term is discretized using second-order central finite differences. Ghost nodes are used to impose Dirichlet and Neumann boundary conditions at the western, eastern, southern, and northern boundaries.
Optional volumetric and shear-heating source terms may be included.
Arguments
D : Structure or tuple containing:
`T`, the temperature field defined at the cell centroids, and
`T_ex`, the extended temperature field including ghost nodes.
Both arrays are updated in place.
κ : Thermal diffusivity.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.
Δt : Time step.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.Keyword Arguments
Q : Volumetric heat-production rate (default: `zeros(NC...)`).
ρ : Density (default: `3300.0`).
cp : Specific heat capacity (default: `1200.0`).
Qₛ : Volumetric shear-heating rate (default: `zeros(NC...)`).Notes
The method is first-order accurate in time and second-order accurate in space. Since the temperature update is explicit, the timestep must satisfy the diffusive stability condition in two dimensions.
For a uniform grid (Δx = Δy = h), the stability criterion is
Δt ≤ h² / (4κ).For non-uniform grids, the timestep should satisfy
Δt ≤ 1 / (2κ(1/Δx² + 1/Δy²)).After each timestep, the updated centroid temperatures are copied to the interior of D.T_ex. The ghost-node values are updated again when the routine is called for the next timestep.
GeoModBox.HeatEquation.TwoD.Poisson2D! — Method
Poisson2D!(T, Q, kx, ky, Δx, Δy, NC, BC, K, rhs, Num)Solves the conservative two-dimensional steady-state heat equation with spatially variable thermal conductivity and optional volumetric heat production.
Temperature is defined at the cell centroids, while the horizontal and vertical thermal conductivities are defined at the corresponding cell faces. The conductive term is discretized in flux form using second-order finite differences.
Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrix and right-hand-side vector at the western, eastern, southern, and northern boundaries. The resulting linear system is solved using a direct left-matrix division.
Arguments
T : Temperature field defined at the cell centroids. The array is
updated in place with the steady-state solution.
Q : Volumetric heat-production rate defined at the cell centroids.
kx : Thermal conductivity defined at the vertical cell faces.
ky : Thermal conductivity defined at the horizontal cell faces.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
K : Sparse coefficient matrix.
rhs : Right-hand-side vector.
Num : Structure or tuple containing the global numbering of the
centroid nodes.Notes
This routine solves the steady-state heat equation in conservative form, where thermal conductivity may vary spatially.
The horizontal and vertical conductive fluxes are evaluated using the face-centered conductivity fields kx and ky. Their divergence is then discretized at the cell centroids using a five-point finite-difference stencil. This formulation conserves conductive heat flux and is suitable for material interfaces with discontinuous thermal conductivity.
Dirichlet and Neumann boundary conditions modify both the matrix coefficients and the right-hand-side vector for equations adjacent to the domain boundaries.
After assembly, the steady-state temperature field is obtained from
KT=rhs,
and written to T.
GeoModBox.HeatEquation.TwoD.Poisson2Dc! — Method
Poisson2Dc!(D, NC, P, BC, Δ, K, rhs, Num)Solves the two-dimensional steady-state heat equation (Poisson equation) assuming constant thermal conductivity and volumetric heat production.
Temperature is defined at the cell centroids, and the diffusion operator is discretized using second-order central finite differences. Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrix and right-hand-side vector at the western, eastern, southern, and northern boundaries.
The resulting linear system is solved using a direct left-matrix division.
Arguments
D : Structure or tuple containing the centroid temperature field
`T` and the volumetric heat-production rate `Q`.
NC : Structure or tuple containing the number of centroid nodes
in the x- and y-directions.
P : Structure or tuple containing the thermal conductivity `k`.
BC : Structure or tuple defining the boundary-condition types and
values at the western, eastern, southern, and northern
boundaries.
Δ : Structure or tuple containing the horizontal and vertical grid
spacings.
K : Sparse coefficient matrix.
rhs : Right-hand-side vector.
Num : Structure or tuple containing the global numbering of the
centroid nodes.Notes
This routine solves the steady-state heat equation where the thermal conductivity is assumed to be constant throughout the domain.
The diffusion operator is discretized using a classical five-point finite difference stencil, resulting in a sparse matrix with five non-zero diagonals for the adopted global numbering.
Dirichlet and Neumann boundary conditions are incorporated directly into the coefficient matrix and right-hand-side vector. After assembly, the steady-state temperature field is obtained from
KT=rhs,
and written to D.T.
GeoModBox.AdvectionEquation.OneD.RK4O1D! — Method
RK4O1D!(x, Δt, vx, xmin, xmax)Advects Lagrangian tracers in one dimension using the classical fourth-order Runge–Kutta (RK4) method.
The tracer trajectories are integrated explicitly over one time step assuming a constant advection velocity. Under this assumption, all four Runge–Kutta stages evaluate the same velocity, reducing the method to an exact integration of the particle trajectories for uniform flow. Tracers leaving one side of the computational domain are reinserted at the opposite side, resulting in periodic boundary conditions.
Arguments
x : Vector containing the tracer coordinates.
Δt : Time step.
vx : Constant horizontal advection velocity.
xmin : Minimum x-coordinate of the computational domain.
xmax : Maximum x-coordinate of the computational domain.Notes
The routine assumes a spatially and temporally constant velocity field during the current time step. Periodic boundary conditions are enforced by wrapping tracers that leave the computational domain to the opposite boundary.
GeoModBox.AdvectionEquation.OneD.lax1D! — Method
lax1D!(A, vx, Δt, Δx)Advects a scalar property using the first-order Lax–Friedrichs finite-difference scheme on a one-dimensional grid.
The Lax–Friedrichs method replaces the property at each grid point by the average of its neighboring values before applying a central-difference approximation of the advection term. This additional averaging introduces artificial numerical diffusion, which stabilizes the explicit scheme and suppresses non-physical oscillations.
Arguments
A : Vector containing the advected property.
vx : Constant horizontal advection velocity.
Δt : Time step.
Δx : Grid spacing.Notes
The advection update is performed explicitly using a combination of neighbor averaging and central finite differences. The implementation assumes a spatially constant velocity during the current time step.
Compared to the first-order upwind scheme, the Lax–Friedrichs method is symmetric with respect to the flow direction but generally introduces stronger numerical diffusion.
For stability, the time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition.
GeoModBox.AdvectionEquation.OneD.semilag1D! — Method
semilag1D!( A, xc, vx, Δt )Advects a scalar property using a semi-Lagrangian method on a one-dimensional grid.
The departure point of each centroid is calculated by tracing the trajectory backward over one time step using the prescribed horizontal velocity. The advected property is then evaluated at the departure points using linear spline interpolation.
Arguments
A : Vector containing the advected property.
xc : Vector containing the centroid x-coordinates.
vx : Horizontal advection velocity.
Δt : Time step.
Δx : Grid spacing.Notes
The implementation assumes that the velocity remains constant during the current time step. Linear interpolation is used to evaluate the property at the departure points.
Unlike explicit Eulerian advection schemes, the semi-Lagrangian method is not restricted by the conventional advective CFL stability condition. However, the time step still affects the trajectory and interpolation accuracy.
GeoModBox.AdvectionEquation.OneD.slf1D! — Method
slf1D!(A, Aold2, vx, Δt, Δx)Advects a scalar property using the staggered leapfrog (SLF) finite-difference scheme on a one-dimensional grid.
The staggered leapfrog method is an explicit second-order accurate advection scheme in both space and time. The spatial derivative is approximated using a central finite difference, while the leapfrog time integration combines the solution from the current and previous time steps. Compared to first-order schemes, the staggered leapfrog method exhibits significantly lower numerical diffusion but may develop dispersive oscillations near sharp gradients.
Arguments
A : Vector containing the advected property at the current time step.
Aold2 : Vector containing the advected property from the previous time step.
vx : Constant horizontal advection velocity.
Δt : Time step.
Δx : Grid spacing.Notes
The advection update is performed explicitly using a central finite-difference approximation of the spatial derivative and leapfrog time integration. The implementation assumes a spatially constant velocity during the current time step.
For stability, the time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition.
GeoModBox.AdvectionEquation.OneD.upwind1D! — Method
upwind1D!(A, vx, Δt, Δx)Advects a scalar property using a first-order upwind finite-difference scheme on a one-dimensional grid.
The spatial derivative is approximated using a one-sided finite difference, with the upwind direction selected according to the sign of the constant advection velocity. The method is explicit, robust, and monotonic, although it introduces numerical diffusion that increases with the Courant number and the number of time steps.
Arguments
A : Vector containing the advected property.
vx : Constant horizontal advection velocity.
Δt : Time step.
Δx : Grid spacing.Notes
The advection update is performed explicitly using first-order upwind finite differences. The implementation assumes a spatially constant velocity during the current time step.
For stability, the time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition.
GeoModBox.AdvectionEquation.TwoD.slfc2D! — Method
slfc2D!(P, P_ex, P_exo, vxc, vyc, NC, Δt, Δx, Δy)Advects a scalar property using the staggered leapfrog (SLF) finite-difference scheme on a two-dimensional finite-difference grid.
The staggered leapfrog method is an explicit second-order accurate advection scheme in both space and time. Spatial derivatives are approximated using central finite differences, while the leapfrog time integration combines the property fields from the current and previous time steps. Compared to the first-order upwind scheme, the staggered leapfrog method exhibits significantly lower numerical diffusion but may develop dispersive oscillations in regions with strong gradients.
Arguments
P : 2-D array containing the property defined on the centroids.
P_ex : 2-D array containing the property including ghost nodes at the current time step.
P_exo : 2-D array containing the property including ghost nodes from the previous time step.
vxc : Horizontal velocity defined on the centroids.
vyc : Vertical velocity defined on the centroids.
NC : Structure or tuple containing the number of centroids.
Δt : Time step.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.Notes
The advection update is performed explicitly using central finite differences for the spatial derivatives and a leapfrog time integration. After each advection step, the current property field is copied to P_exo, and the newly computed solution is written back to the interior of P_ex. The ghost nodes can subsequently be updated by the appropriate boundary-condition routine before the next time step.
For stability, the time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition.
GeoModBox.AdvectionEquation.TwoD.upwindc2D! — Method
upwindc2D!(P, P_ex, vxc, vyc, NC, Δt, Δx, Δy)Advects a scalar property using a first-order upwind finite-difference scheme on a two-dimensional finite-difference grid.
The spatial derivatives are evaluated using one-sided finite differences selected according to the sign of the local velocity components. This upwind discretization provides a robust and monotonic solution, although it introduces numerical diffusion, particularly when the flow is not aligned with the computational grid.
Arguments
P : 2-D array containing the property defined on the centroids.
P_ex : 2-D array containing the property including ghost nodes.
vxc : Horizontal velocity defined on the centroids.
vyc : Vertical velocity defined on the centroids.
NC : Structure or tuple containing the number of centroids.
Δt : Time step.
Δx : Horizontal grid spacing.
Δy : Vertical grid spacing.Notes
The advection update is performed explicitly using first-order upwind differences in both coordinate directions. After each advection step, the interior values are copied back to the extended array containing the ghost nodes, allowing the boundary conditions to be updated before the next time step.
For stability, the time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition.
GeoModBox.InitialCondition.IniPhase! — Method
IniPhase!(type, D, M, x, y, NC; phase=(0, 1))Initialize a two-dimensional phase field at the cell centers of a regular finite-difference grid.
This function defines phases directly on the Eulerian grid and is retained primarily for compatibility with older examples. For models involving material transport, phases should preferably be initialized on tracers and subsequently interpolated to the grid.
Arguments
- `type`: Symbol selecting the initial phase distribution.
- `D` : Data structure containing the cell-centered phase field `p` and the
extended phase field `p_ex`.
- `M` : Model geometry containing `xmin`, `xmax`, `ymin`, and `ymax`.
- `x` : Horizontal cell-center coordinates.
- `y` : Vertical cell-center coordinates.
- `NC` : Number of cell centers in the horizontal and vertical directions.Keyword arguments
- `phase=(0, 1)`: Two phase values given as
`(background_phase, anomaly_phase)`.Supported phase distributions
- `:block`: Rectangular anomaly located between 40% and 60% of the model
width and between 10% and 30% of the model height below the upper
boundary.The interior phase field is assigned to D.p. The corresponding extended field D.p_ex is then updated using constant extrapolation into the ghost cells.
Example
```julia IniPhase!(:block,D,M,x,y,NC;phase = (0, 1))
GeoModBox.InitialCondition.IniTemperature! — Method
IniTemperature!(type,M,NC,D,x,y;Tb = 600.0,Ta = 1200.0,σ = 0.1)Initialize a two-dimensional temperature field on the extended cell-centered finite-difference grid.
The selected initial condition is assigned to D.T_ex, including its ghost-cell layers. The interior temperature field D.T is subsequently updated from
D.T_ex[2:end-1, 2:end-1]If available, the auxiliary fields D.T0, D.Told_ex, D.T_exD0, and D.T_ex0 are initialized from the resulting temperature field.
Arguments
- `type`: Symbol selecting the initial temperature distribution.
- `M` : Model geometry containing `xmin`, `xmax`, `ymin`, and `ymax`.
- `NC` : Number of cell centers in the horizontal and vertical directions.
- `D` : Data structure containing the temperature fields.
- `x` : Horizontal coordinates, including the extended centroid coordinates
`x.ce`.
- `y` : Vertical coordinates, including the extended centroid coordinates
`y.ce`.Keyword arguments
- `Tb=600.0` : Background temperature for localized anomalies or bottom
temperature for vertically varying profiles.
- `Ta=1200.0` : Anomaly temperature for localized anomalies or top
temperature for vertically varying profiles.
- `σ=0.1` : Nondimensional width of the Gaussian anomaly.Initial-condition types
- `:const` : Spatially uniform temperature equal to `Tb`.
- `:circle` : Circular or elliptical region with temperature `Ta` embedded
in a background temperature `Tb`.
- `:gaussian` : Gaussian temperature anomaly with background temperature
`Tb` and peak temperature `Ta`.
- `:block` : Rectangular region with temperature `Ta` embedded in a
background temperature `Tb`.
- `:linear` : Linear conductive temperature profile between the prescribed
bottom and top temperatures.
- `:lineara` : Linear conductive temperature profile with an additional
elliptical temperature anomaly.
- `:blankenbach`: Conductive temperature profile with the small sinusoidal
perturbation used to initialize the Blankenbach convection benchmark.Example
IniTemperature!(:circle,M,NC,D,x,y;Tb = 1200.0,Ta = 0.0,)GeoModBox.InitialCondition.IniVelocity! — Method
IniVelocity!(type,D,BC,NV,Δ,M,x,y;ε = 1e-15)Initialize a predefined two-dimensional velocity field or update the velocity boundary values on a staggered finite-difference grid.
For the rigid-body and shear-cell configurations, the complete staggered velocity fields D.vx and D.vy are initialized. For simple-shear and pure-shear configurations, the corresponding entries of BC.val are updated and written to the velocity nodes and ghost nodes located along the model boundaries. Interior velocity values are retained in these cases, allowing the function to update deformation boundary conditions after a grid-remeshing step without erasing the current Stokes solution.
Arguments
- `type`: Symbol selecting the velocity configuration.
- `D` : Data structure containing the staggered velocity fields `vx` and
`vy`.
- `BC` : Velocity boundary-condition structure containing `type` and `val`.
- `NV` : Number of grid vertices in the horizontal and vertical directions.
- `Δ` : Grid spacing containing `x` and `y`.
- `M` : Model geometry containing `xmin`, `xmax`, `ymin`, and `ymax`.
- `x` : Horizontal coordinates at cell centers, vertices, and staggered
velocity locations.
- `y` : Vertical coordinates at cell centers, vertices, and staggered
velocity locations.Keyword arguments
- `ε=1e-15` : Background deformation rate in s⁻¹ used by the simple- and
pure-shear velocity configurations.Velocity configurations
- `:RigidBody` : Clockwise rigid-body rotation within a circular region.
The velocity decreases to zero outside the prescribed rotating region.
- `:ShearCell` : Analytic cellular velocity field used for advection tests.
- `:SimpleShear`: Simple-shear deformation with horizontal velocity varying
linearly with vertical position and zero vertical velocity.
- `:PureShear` : Domain-centered incompressible pure shear,\[ v_x = -\dot{\varepsilon}(x-x_c), \qquad v_y = \dot{\varepsilon}(y-y_c). ``` - `:ShearBandPS`: Origin-centered pure-shear field used by the deforming shear-band experiment, ```math v_x = -\dot{\varepsilon}x, \qquad v_y = \dot{\varepsilon}y. ``` For `:SimpleShear`, `:PureShear`, and `:ShearBandPS`, the boundary types in `BC.type` must be configured consistently with the selected deformation field before calling this function. # Returns Returns the modified data and boundary-condition structures as `(D, BC)`. # Example \]
julia D, VBC = IniVelocity!(:PureShear,D,VBC,NV,Δ,M,x,y;ε = 1e-15) ```
GeoModBox.Tracers.OneD.Itp1D_Centers2Markers! — Method
Itp1D_Centers2Markers!(Tm, xm, Tc, xc, Δx)Interpolate a one-dimensional cell-centered field to marker positions using linear interpolation.
For each marker coordinate in xm, the function identifies the two neighboring cell centers in xc and interpolates the corresponding values of Tc. The interpolated marker values are written in place to Tm.
Arguments
Tm: Array receiving the interpolated marker values.xm: Marker coordinates.Tc: Field values defined at the cell centers.xc: Cell-center coordinates.Δx: Uniform spacing between adjacent cell centers.
Notes
The arrays Tm and xm must have the same length, while Tc and xc must contain the same number of entries. At least two cell centers are required.
Marker positions outside the interval spanned by xc are linearly extrapolated using the first or last pair of cell centers.
Example
xc = [0.5, 1.5, 2.5, 3.5]
Tc = [10.0, 20.0, 30.0, 40.0]
xm = [1.0, 2.0, 3.0]
Tm = similar(xm)
Itp1D_Centers2Markers!(Tm, xm, Tc, xc, 1.0)GeoModBox.Tracers.OneD.Itp1D_Markers2Centers! — Method
Itp1D_Markers2Centers!(Tc, xc, Tm, xm, Δx)Interpolate values from one-dimensional markers to cell centers using linear weighted averaging.
Each marker contributes to the two adjacent cell centers according to its relative position between them. Marker contributions and interpolation weights are accumulated separately, after which the value at each cell center is obtained by dividing the weighted sum by the total weight.
Arguments
Tc: Cell-centered field receiving the interpolated marker values.xc: Coordinates of the cell centers.Tm: Values carried by the markers.xm: Marker coordinates.Δx: Uniform spacing between adjacent cell centers.
Notes
The arrays Tm and xm must have the same length, and Tc and xc must contain the same number of entries. At least two cell centers are required.
Markers located outside the interval spanned by xc are assigned to the nearest cell center. Cell centers receiving no marker contributions retain their previous value in Tc.
Example
xc = [0.5, 1.5, 2.5, 3.5]
xm = [1.0, 2.0, 3.0]
Tm = [10.0, 20.0, 30.0]
Tc = zeros(length(xc))GeoModBox.Tracers.TwoD.Markers — Type
Ma = Markers(x, y, phase)Container storing the Lagrangian marker positions and phase IDs.
The structure is primarily used for marker-based material advection. Each marker stores its Cartesian coordinates together with an integer phase ID, allowing the phase distribution to evolve independently of the Eulerian computational grid.
Fields
- `x` : Horizontal marker coordinates.
- `y` : Vertical marker coordinates.
- `phase` : Integer phase identifier carried by each marker.All arrays must have identical length corresponding to the total number of markers.
GeoModBox.Tracers.TwoD.TMarkers — Type
Ma = TMarkers(x, y, T, phase)Container storing the Lagrangian marker properties used by the advection routines.
Each marker stores its Cartesian coordinates, temperature, and phase ID. The structure is primarily used for marker-based temperature advection and for tracking material interfaces during thermo-mechanical simulations.
Fields
- `x` : Horizontal marker coordinates.
- `y` : Vertical marker coordinates.
- `T` : Marker temperature.
- `phase` : Integer phase identifier carried by each marker.All arrays must have identical length corresponding to the total number of markers.
GeoModBox.Tracers.TwoD.AdvectTracer2D — Method
AdvectTracer2D(
Ma,
nmark,
D,
x,
y,
dt,
Δ,
NC,
rkw,
rkv;
style = 1,
)Advect two-dimensional Lagrangian tracers using a fourth-order Runge–Kutta scheme.
For each active tracer, the velocity is interpolated from the Eulerian grid at the intermediate Runge–Kutta positions. The stage velocities are combined using the coefficients in rkw, while rkv defines the intermediate positions used during the four Runge–Kutta stages.
Only tracers with a non-negative phase ID are advected.
Arguments
Ma: Marker structure containing the tracer coordinates and phase IDs.nmark: Total number of tracers.D: Data structure containing the staggered and cell-centered velocity fields.x: Horizontal coordinates of the staggered and cell-centered grids.y: Vertical coordinates of the staggered and cell-centered grids.dt: Time-step size.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.rkw: Runge–Kutta weights used to combine the stage velocities.rkv: Runge–Kutta coefficients used to calculate the intermediate tracer positions.
Keyword arguments
style=1: Velocity-interpolation method:1: Bilinear interpolation from the staggeredvₓandvᵧnodes.2: Combined interpolation using staggered and cell-centered velocities.3: Interpolation from the staggered velocity nodes with the optional higher-order correction enabled.
Returns
Returns the modified marker structure Ma. The tracer coordinates Ma.x and Ma.y are updated in place.
Notes
Markers with
Ma.phase[k] < 0are treated as inactive and are not advected.
The velocity-interpolation routines are evaluated at the intermediate Runge–Kutta positions by temporarily updating the tracer coordinates during each stage.
Example
rkw = (1.0 / 6.0) .* [1.0, 2.0, 2.0, 1.0] rkv = (1.0 / 2.0) .* [1.0, 1.0, 2.0, 2.0]
AdvectTracer2D( Ma, nmark, D, x, y, dt, Δ, NC, rkw, rkv; style = 1, ) ```
GeoModBox.Tracers.TwoD.CountMPC — Method
CountMPC(Ma, nmark, MPC, M, x, y, Δ, NC, NV; verbose=false)Count the number of active markers associated with each cell center and grid vertex.
Markers located outside the current model domain are deactivated by setting their phase ID to -1. The remaining active markers are counted separately for the cell-centered and vertex-centered grids. Thread-local counting arrays are used during the parallel loops and are subsequently summed to obtain the total marker counts stored in MPC.c and MPC.v.
Arguments
Ma: Marker structure containing the marker coordinates and phase IDs.nmark: Total number of markers.MPC: Structure containing the marker-count arrays:MPC.c: Total number of active markers associated with each cell center.MPC.v: Total number of active markers associated with each grid vertex.MPC.th: Thread-local marker counts for the cell-centered grid.MPC.thv: Thread-local marker counts for the vertex-centered grid.
M: Model geometry containing the domain limits.x: Horizontal coordinates of the cell centers and vertices.y: Vertical coordinates of the cell centers and vertices.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.NV: Number of grid vertices in each coordinate direction.
Keyword arguments
verbose=false: Print the number of markers located outside the model domain.
Notes
Only markers with
Ma.phase[k] >= 0are included in the marker counts. Markers outside the domain are assigned a phase ID of -1 and are excluded from subsequent interpolation and advection operations.
The function modifies Ma.phase, MPC.c, MPC.v, MPC.th, and MPC.thv in place.
Example
CountMPC( Ma, nmark, MPC, M, x, y, Δ, NC, NV; verbose = false, )
GeoModBox.Tracers.TwoD.FromCtoM — Method
FromCtoM(Prop, k, Ma, x, y, Δ, NC)Interpolate a cell-centered property to the position of marker k.
The function identifies the four surrounding cell centers enclosing the marker position and evaluates the property using bilinear interpolation. Marker coordinates located close to the model boundaries are restricted to valid interpolation indices to ensure that only existing grid values are used.
This routine can be used to initialize or update marker properties from Eulerian fields, such as temperature, density, viscosity, or any other scalar quantity defined at the cell centers.
Arguments
Prop: Property field defined at the cell centers.k: Index of the marker for which the property is evaluated.Ma: Marker structure containing the marker coordinates.x: Horizontal coordinates of the extended cell centers.y: Vertical coordinates of the extended cell centers.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.
Returns
Returns the interpolated property value at the position of marker k.
Example
for k = 1:nmark
Ma.T[k] = FromCtoM(D.T_ex, k, Ma, x, y, Δ, NC)
endGeoModBox.Tracers.TwoD.IniTracer2D — Method
IniTracer2D(
Aparam,
nmx,
nmy,
Δ,
M,
NC,
noise,
ini,
phase;
xc,
yc,
λ,
δA,
ellA,
ellB,
α,
)Initialize a regular two-dimensional distribution of Lagrangian tracers and assign their initial phase IDs.
The tracers are distributed uniformly within each finite-difference cell. Depending on Aparam, the returned marker structure stores either phase information only or both phase and temperature. An optional random perturbation can be applied to the initial tracer coordinates.
Arguments
- `Aparam` : Defines the marker properties:
- `:phase` initializes markers carrying phase IDs.
- `:thermal` initializes markers carrying phase IDs and temperature.
- `nmx` : Number of tracers per cell in the horizontal direction.
- `nmy` : Number of tracers per cell in the vertical direction.
- `Δ` : Grid spacing containing the horizontal and vertical increments.
- `M` : Model geometry containing the domain limits.
- `NC` : Number of finite-difference cells in the horizontal and vertical
directions.
- `noise` : Controls whether random perturbations are added to the tracer
positions. Use `1` to add noise and `0` to retain the regular distribution.
- `ini` : Symbol selecting the initial phase distribution.
- `phase` : Collection containing the phase IDs. `phase[1]` represents the
background or matrix phase, while `phase[2]` represents the anomaly,
inclusion, or second layer.Keyword arguments
- `xc=(M.xmax-M.xmin)` : Horizontal reference coordinate used by the
Rayleigh–Taylor instability setup.
- `yc=(M.ymax-M.ymin)/2`: Vertical reference coordinate used by the
Rayleigh–Taylor instability setup.
- `λ=1.0e3` : Wavelength of the cosine perturbation used for `:RTI`.
- `δA=5e2/15` : Amplitude of the perturbation used for `:RTI`.
- `ellA=100.0` : Major semiaxis of the elliptical inclusion or radius of the
circular shear-band inclusion.
- `ellB=100.0` : Minor semiaxis of the elliptical inclusion.
- `α=0.0` : Rotation angle of the elliptical inclusion in degrees.Initial phase distributions
- `:block` : Rectangular anomaly embedded in a homogeneous background phase.
- `:RTI` : Two-layer configuration with a cosine perturbation of the tracer
positions for a Rayleigh–Taylor instability.
- `:Inclusion` : Elliptical inclusion located at the center of the model
domain. The inclusion may be rotated by the angle `α`.
- `:ShearBandSetting` : Circular weak inclusion centered at the lower
boundary for the thermo-mechanical shear-localization experiment.Returns
Returns either a Markers or TMarkers structure containing the initialized tracer coordinates, phase IDs, and, where applicable, temperature values.
Example
Ma = IniTracer2D(
:phase,
nmx,
nmy,
Δ,
M,
NC,
1,
:RTI,
[0, 1],
)GeoModBox.Tracers.TwoD.Markers2Cells — Method
Markers2Cells(
Ma,
nmark,
PC_th,
PC,
weight_th,
weight,
x,
y,
Δ,
param,
param2;
avgm = :arith,
)Interpolate a marker property to the extended cell-centered grid using weighted bilinear interpolation.
Each active marker contributes to the four surrounding cell centers according to its relative position within the corresponding grid cell. Task-local property and weight arrays are used to accumulate the marker contributions in parallel. The local arrays are subsequently summed and normalized to obtain the interpolated cell-centered field.
The function can interpolate either marker temperature directly or a phase-dependent material property specified through param2.
Arguments
Ma: Marker structure containing marker coordinates, phase IDs, and, where applicable, marker temperatures.nmark: Total number of markers.PC_th: Collection of local property arrays used by the parallel
interpolation tasks defined on the extended cell-centered grid.
PC: Extended cell-centered array receiving the interpolated property.weight_th: Collection of local property interpolation-weight arrays.weight: Extended cell-centered array receiving the accumulated weights.x: Horizontal coordinates of the extended cell centers.y: Vertical coordinates of the extended cell centers.Δ: Grid spacing in the horizontal and vertical directions.param: Property to interpolate::thermalinterpolatesMa.T.:phaseinterpolates the phase-dependent values contained inparam2.
param2: Collection containing the material-property value assigned to each phase. Phase IDpaccesses the valueparam2[p+1]. This argument is not used whenparam == :thermal.
Keyword arguments
avgm=:arith: Averaging method used for phase-dependent properties::arith: Weighted arithmetic average.:harm: Weighted harmonic average.:geom: Weighted geometric average.
For thermal interpolation, the marker temperatures are combined using the weighted arithmetic average.
Notes
The number of arrays in PC_th and weight_th must be at least nthreads(:default).
Grid points receiving no marker contributions retain their previous values from PC. The function modifies PC, weight, PC_th, and weight_th in place.
Example
Markers2Cells(
Ma,
nmark,
MAVG.PC_th,
D.ρ_ex,
MAVG.wte_th,
D.wte,
x,
y,
Δ,
:phase,
ρ;
avgm = :arith,
)
D.ρ .= D.ρ_ex[2:end-1, 2:end-1]GeoModBox.Tracers.TwoD.Markers2Vertices — Method
Markers2Vertices(
Ma,
nmark,
PG_th,
PG,
weight_th,
weight,
x,
y,
Δ,
param,
param2;
avgm = :arith,
)Interpolate a marker property to the grid vertices using weighted bilinear interpolation.
Each marker contributes to the four surrounding vertices according to its relative position within the corresponding grid cell. Task-local property and weight arrays are used to accumulate the marker contributions in parallel. The local arrays are subsequently summed and normalized to obtain the interpolated vertex field.
The function can interpolate either marker temperature directly or a phase-dependent material property specified through param2.
Arguments
- `Ma`: Marker structure containing marker coordinates, phase IDs, and,
where applicable, marker temperatures.
- `nmark`: Total number of markers.
- `PG_th`: Collection of local property arrays used by the parallel
interpolation tasks defined at the grid vertices.
- `PG`: Vertex-centered array receiving the interpolated property.
- `weight_th`: Collection of local property interpolation-weight arrays
defined at the vertices.
- `weight`: Vertex-centered array receiving the accumulated interpolation
weights.
- `x`: Horizontal coordinates of the grid vertices.
- `y`: Vertical coordinates of the grid vertices.
- `Δ`: Grid spacing in the horizontal and vertical directions.
- `param`: Property to interpolate:
- `:thermal` interpolates `Ma.T`.
- `:phase` interpolates the phase-dependent values contained in `param2`.
- `param2`: Collection containing the material-property value assigned to
each phase. Phase ID `p` accesses the value `param2[p+1]`. This argument
is not used when `param == :thermal`.Keyword arguments
- `avgm=:arith`: Averaging method used for phase-dependent properties:
- `:arith`: Weighted arithmetic average.
- `:harm`: Weighted harmonic average.
- `:geom`: Weighted geometric average.For thermal interpolation, marker temperatures are combined using the weighted arithmetic average.
Notes
The number of arrays in PG_th and weight_th must be at least nthreads(:default).
Vertices receiving no marker contributions retain their previous values from PG. The function modifies PG, weight, PG_th, and weight_th in place.
Example
Markers2Vertices(
Ma,
nmark,
MAVG.PV_th,
D.ηv,
MAVG.wtv_th,
D.wtv,
x,
y,
Δ,
:phase,
η;
avgm = :arith,
)GeoModBox.Tracers.TwoD.VxFromVxNodes — Method
VxFromVxNodes(Vx, k, Ma, x, y, Δ, NC, new)Interpolate the horizontal velocity component from the staggered vₓ nodes to the position of marker k.
The function first determines the surrounding staggered velocity nodes and performs bilinear interpolation in the horizontal and vertical directions. Marker coordinates located close to the model boundaries are restricted to valid interpolation indices.
An optional higher-order velocity correction can be activated through new. The implementation is adapted from the marker-advection routines in M2Dpt_Julia.
Arguments
Vx: Horizontal velocity field defined at the staggeredvₓnodes.k: Index of the marker for which the velocity is evaluated.Ma: Marker structure containing the marker coordinates.x: Horizontal coordinates of the staggered velocity nodes.y: Vertical coordinates of the extended cell centers.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.new: Switch controlling the interpolation formulation:0: Use bilinear interpolation.1: Apply the additional higher-order velocity correction where sufficient neighboring nodes are available.
Returns
Returns the interpolated horizontal marker velocity vxm.
Example
vxm = VxFromVxNodes(
D.vx,
k,
Ma,
x,
y,
Δ,
NC,
0,
)Reference
The interpolation formulation is modified from the marker-advection implementation in M2Dpt_Julia:
https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
GeoModBox.Tracers.TwoD.VxVyFromPrNodes — Method
VxVyFromPrNodes(Vxp, Vyp, k, Ma, x, y, Δ, NC)Interpolate the horizontal and vertical velocity components from the cell-centered grid to the position of marker k.
The function identifies the four surrounding cell-center nodes and applies bilinear interpolation to both velocity components. Marker coordinates close to the model boundaries are restricted to valid interpolation indices.
Arguments
Vxp: Horizontal velocity field defined at the cell centers.Vyp: Vertical velocity field defined at the cell centers.k: Index of the marker for which the velocity is evaluated.Ma: Marker structure containing the marker coordinates.x: Horizontal coordinates of the extended cell centers.y: Vertical coordinates of the extended cell centers.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.
Returns
Returns the interpolated horizontal and vertical marker velocities as
vxm, vymExample
vxm, vym = VxVyFromPrNodes( D.vxc, D.vyc, k, Ma, x, y, Δ, NC, )
Reference
The interpolation formulation is modified from the marker-advection implementation in M2Dpt_Julia:
https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
GeoModBox.Tracers.TwoD.VyFromVyNodes — Method
VyFromVyNodes(Vy, k, Ma, x, y, Δ, NC, new)Interpolate the vertical velocity component from the staggered vᵧ nodes to the position of marker k.
The function first identifies the surrounding staggered velocity nodes and computes the marker velocity using bilinear interpolation. Marker coordinates close to the model boundaries are clamped to valid interpolation indices to ensure that only existing grid nodes are used.
An optional higher-order interpolation correction can be enabled through new. The implementation is adapted from the marker-advection routines in M2Dpt_Julia.
Arguments
Vy: Vertical velocity field defined at the staggeredvᵧnodes.k: Index of the marker for which the velocity is evaluated.Ma: Marker structure containing the marker coordinates.x: Horizontal coordinates of the extended cell centers.y: Vertical coordinates of the staggered velocity nodes.Δ: Grid spacing in the horizontal and vertical directions.NC: Number of cell centers in each coordinate direction.new: Switch controlling the interpolation formulation:0: Use bilinear interpolation.1: Apply the additional higher-order velocity correction where sufficient neighboring nodes are available.
Returns
Returns the interpolated vertical marker velocity vym.
Example
vym = VyFromVyNodes(
D.vy,
k,
Ma,
x,
y,
Δ,
NC,
0,
)Reference
The interpolation formulation is modified from the marker-advection implementation in M2Dpt_Julia:
https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
GeoModBox.MomentumEquation.OneD.AssembleStokesMatrix1D — Method
AssembleStokesMatrix1D(nc, η, Δy, BC, K)Assemble the sparse coefficient matrix for the one-dimensional incompressible Stokes equation
∂/∂y(η ∂vₓ/∂y) = ∂P/∂x.The function constructs the variable-viscosity finite-difference stencil and incorporates the prescribed Dirichlet or Neumann boundary conditions into the matrix coefficients. A system with Neumann conditions at both boundaries is rejected because the resulting matrix is singular.
Arguments
- `nc` : Number of computational cells.
- `η` : Dynamic viscosity defined at the cell faces.
- `Δy` : Grid spacing in the y-direction.
- `BC` : Boundary condition structure specifying the southern and
northern boundary types and values.
- `K` : Extendable sparse matrix used to store the coefficient matrix.Returns
The finalized sparse coefficient matrix.
GeoModBox.MomentumEquation.OneD.ComputeStokesResiduals1D! — Method
ComputeStokesResiduals1D!(D, ∂P∂x, Δy, BC)Compute the residual of the one-dimensional incompressible Stokes equation
-∂P∂x + ∂τxy/∂y = 0,where the shear stress is given by
τxy = η ∂vₓ/∂y.The function applies the prescribed Dirichlet or Neumann boundary conditions using ghost nodes, computes the velocity gradient, shear stress, and stress divergence, and stores the resulting residual in D.R.
Arguments
- `D` : Structure containing the velocity field, viscosity, stress,
derivatives, and residual arrays.
- `∂P∂x`: Prescribed pressure gradient (scalar or array).
- `Δy` : Grid spacing in the y-direction.
- `BC` : Boundary condition structure specifying the boundary types
and values.Returns
Nothing. The residual and intermediate quantities are updated in-place.
GeoModBox.MomentumEquation.OneD.Stokes_1D_direct — Method
Stokes_1D_direct(vₓ, η, Δy, nc, BC, K, rhs)Solve the one-dimensional incompressible Stokes equation using a direct sparse linear solver.
The function assembles the finite-difference coefficient matrix, incorporates the prescribed Dirichlet or Neumann boundary conditions, modifies the right-hand side accordingly, and solves the resulting linear system for the horizontal velocity. A system with Neumann conditions at both boundaries is rejected because the resulting matrix is singular.
Arguments
- `vₓ` : Array storing the computed horizontal velocity.
- `η` : Dynamic viscosity defined at the cell faces.
- `Δy` : Grid spacing in the y-direction.
- `nc` : Number of computational cells.
- `BC` : Boundary condition structure specifying the southern and
northern boundary types and values.
- `K` : Extendable sparse matrix used to assemble the coefficient matrix.
- `rhs` : Right-hand side vector representing the pressure gradient and
any additional body forces.Returns
The computed horizontal velocity vₓ.
GeoModBox.MomentumEquation.TwoD.Assembly — Method
Assembly(NC, NV, Δ, ηc, ηv, BC, Num)Assemble the sparse coefficient matrix for the two-dimensional incompressible Stokes equations with spatially variable viscosity.
The function constructs the finite-difference stencils for the horizontal and vertical momentum equations and for the continuity equation on a staggered grid. Cell-centered viscosities are used for the normal stress components, whereas vertex-centered viscosities are used for the shear stress components. The prescribed boundary conditions are incorporated into the matrix coefficients, and one pressure degree of freedom is fixed to remove the pressure nullspace.
Arguments
- `NC` : Number of computational cells in the x- and y-directions.
- `NV` : Number of velocity nodes in the x- and y-directions.
- `Δ` : Grid spacing in the x- and y-directions.
- `ηc` : Cell-centered dynamic viscosity field.
- `ηv` : Vertex-centered dynamic viscosity field.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `Num` : Structure containing the global equation numbers for `vₓ`,
`vᵧ`, and pressure.Returns
The finalized sparse coefficient matrix.
GeoModBox.MomentumEquation.TwoD.Assemblyc — Method
Assemblyc(NC, NV, Δ, η, BC, Num)Assemble the sparse coefficient matrix for the two-dimensional incompressible Stokes equations with constant viscosity.
The function constructs the finite-difference stencils for the horizontal and vertical momentum equations and for the continuity equation on a staggered grid. The prescribed boundary conditions are incorporated into the matrix coefficients, and one pressure degree of freedom is fixed to remove the pressure nullspace.
Arguments
- `NC` : Number of computational cells in the x- and y-directions.
- `NV` : Number of velocity nodes in the x- and y-directions.
- `Δ` : Grid spacing in the x- and y-directions.
- `η` : Constant dynamic viscosity.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `Num` : Structure containing the global equation numbers for `vₓ`,
`vᵧ`, and pressure.Returns
The finalized sparse coefficient matrix.
GeoModBox.MomentumEquation.TwoD.CheckStokesBoundaryTypes — Method
CheckStokesBoundaryTypes(BC)Validate the boundary condition types used by the two-dimensional Stokes solver.
The function verifies that the boundary conditions prescribed on the western, eastern, southern, and northern boundaries are among the supported types. An ArgumentError is thrown if an unsupported boundary condition is encountered.
Supported boundary condition types are: - :freeslip - :noslip - :const - :ps
Arguments
- `BC`: Boundary condition structure containing the boundary types.Returns
Nothing. Throws an ArgumentError if an invalid boundary condition type is detected.
GeoModBox.MomentumEquation.TwoD.Residuals2D! — Method
Residuals2D!(D, BC, ε, τ, divV, Δ, ηc, ηv, g, Fm, FPt)Compute the momentum and continuity residuals for the two-dimensional incompressible Stokes equations with spatially variable viscosity.
The function first applies the prescribed velocity boundary conditions. It then evaluates the velocity divergence, strain-rate components, viscous stresses, momentum residuals, and continuity residual on the staggered grid. Cell-centered viscosities are used for the normal stress components, whereas vertex-centered viscosities are used for the shear stress component.
The residuals are computed consistently with the variable-viscosity operator assembled by Assembly and are intended for use in the defect-correction solver. All output fields are updated in place.
Arguments
- `D` : Structure containing the velocity, pressure, and density fields.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `ε` : Structure storing the strain-rate components.
- `τ` : Structure storing the viscous stress components.
- `divV`: Cell-centered velocity-divergence field.
- `Δ` : Grid spacing in the x- and y-directions.
- `ηc` : Cell-centered dynamic viscosity field.
- `ηv` : Vertex-centered dynamic viscosity field.
- `g` : Gravitational acceleration.
- `Fm` : Structure storing the horizontal and vertical momentum residuals.
- `FPt` : Array storing the continuity residual.Returns
Nothing. The strain-rate, stress, divergence, and residual fields are updated in place.
GeoModBox.MomentumEquation.TwoD.Residuals2Dc! — Method
Residuals2Dc!(D, BC, ε, τ, divV, Δ, η, g, Fm, FPt)Compute the momentum and continuity residuals for the two-dimensional incompressible Stokes equations with constant viscosity.
The function first applies the prescribed velocity boundary conditions. It then evaluates the velocity divergence, strain-rate components, viscous stresses, momentum residuals, and continuity residual on the staggered grid. The residuals are computed consistently with the constant-viscosity operator assembled by Assemblyc and are intended for use in the defect-correction solver.
All output fields are updated in place.
Arguments
- `D` : Structure containing the velocity, pressure, and density fields.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `ε` : Structure storing the strain-rate components.
- `τ` : Structure storing the viscous stress components.
- `divV`: Cell-centered velocity-divergence field.
- `Δ` : Grid spacing in the x- and y-directions.
- `η` : Constant dynamic viscosity.
- `g` : Gravitational acceleration.
- `Fm` : Structure storing the horizontal and vertical momentum residuals.
- `FPt` : Array storing the continuity residual.Returns
Nothing. The strain-rate, stress, divergence, and residual fields are updated in place.
GeoModBox.MomentumEquation.TwoD.updaterhs — Method
updaterhs(NC, NV, Δ, ηc, ηv, ρ, g, BC, Num)Assemble the right-hand-side vector for the two-dimensional incompressible Stokes equations with spatially variable viscosity.
The function incorporates the contributions from prescribed velocity boundary conditions and the gravitational body force into the linear system. Vertex-centered viscosities are used for the boundary terms associated with the shear stress components. The returned right-hand-side vector is compatible with the coefficient matrix assembled by Assembly.
Arguments
- `NC` : Number of computational cells in the x- and y-directions.
- `NV` : Number of velocity nodes in the x- and y-directions.
- `Δ` : Grid spacing in the x- and y-directions.
- `ηc` : Cell-centered dynamic viscosity field.
- `ηv` : Vertex-centered dynamic viscosity field.
- `ρ` : Cell-centered density field.
- `g` : Gravitational acceleration.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `Num` : Structure containing the global equation numbers for `vₓ`,
`vᵧ`, and pressure.Returns
The right-hand-side vector of the linear Stokes system.
GeoModBox.MomentumEquation.TwoD.updaterhsc — Method
updaterhsc(NC, NV, Δ, η, ρ, g, BC, Num)Assemble the right-hand-side vector for the two-dimensional incompressible Stokes equations with constant viscosity.
The function incorporates the contributions from prescribed velocity boundary conditions and the gravitational body force into the linear system. The returned right-hand-side vector is compatible with the coefficient matrix assembled by Assemblyc.
Arguments
- `NC` : Number of computational cells in the x- and y-directions.
- `NV` : Number of velocity nodes in the x- and y-directions.
- `Δ` : Grid spacing in the x- and y-directions.
- `η` : Constant dynamic viscosity.
- `ρ` : Cell-centered density field.
- `g` : Gravitational acceleration.
- `BC` : Boundary condition structure specifying the boundary types and
values.
- `Num` : Structure containing the global equation numbers for `vₓ`,
`vᵧ`, and pressure.Returns
The right-hand-side vector of the linear Stokes system.