List of all functions
Here an overview of all functions:
GeoModBox.DataFields — Type
DataFields()Structure to initialize some data arrays.
The default values are:
Q = zeros(1,1)
T = zeros(1,1)
T0 = zeros(1,1)
T_ex = zeros(1,1)
T_exo = zeros(1,1)
ηc = zeros(1,1)
η_ex = zeros(1,1)
ηv = zeros(1,1)
ρ = zeros(1,1)
ρ_ex = zeros(1,1)
cp = zeros(1,1)
vx = zeros(1,1)
vy = zeros(1,1)
Pt = zeros(1,1)
vxc = zeros(1,1)
vyc = zeros(1,1)
vc = zeros(1,1)
wt = zeros(1,1)
wte = zeros(1,1)
wtv = zeros(1,1)
ΔTtop = zeros(1)
ΔTbot = zeros(1)
Tmax = 0.0
Tmin = 0.0
Tmean = 0.0Examples
julia> D = DataFields()
DataFields([0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [0.0], [0.0], 0.0, 0.0, 0.0)GeoModBox.Geometry — Type
Geometry()Structure to initialize the geometry of the model domain.
The default values are:
xmin = 0.0
xmax = 1.0
ymin = -1.0
ymax = 0.0Examples
julia> M = Geometry()
Geometry(0.0, 1.0, -1.0, 0.0)GeoModBox.GridSpacing — Type
GridSpacing()Structure to initialize the horizontal and vertical grid spacing.
The default values are:
x = 0.0
y = 0.0Examples
julia> Δ = GridSpacing()
GridSpacing(0.0, 0.0)GeoModBox.Physics — Type
Physics()Structure to initialize some physical constants in SI units.
The default values are:
g = 9.81 # Gravitational acceleration [ m/s² ]
ρ₀ = 3300.0 # Reference density [ kg/m³ ]
k = 4.125 # Thermal conductivity [ W/m/K ]
cp = 1250.0 # Specific heat capacity [ J/kg/K ]
α = 2.0e-5 # Thermal expansion coefficient [ 1/K ]
Q₀ = 0.0 # Heat production rate [ W/m³ ]
η₀ = 3.947725485e23 # Reference viscosity [ Pa s ]
κ = k/ρ₀/cp # Thermal Diffusivity [ m²/s ]
ΔT = 2500 # Temperature difference [ K ]
Ttop = 273.15 # Temperature at the top [ K ]
Tbot = Ttop + ΔT # Temperature at the bottom [ K ]
Ra = -9999 # Rayleigh numberThe Rayleigh number is set to a negative value, which results into the calculation of the basal Rayleigh number based on the given physical constants.
Examples
julia> P = Physics()
Physics(9.81, 3300.0, 4.125, 1250.0, 2.0e-5, 0.0, 3.947725485e23, 1.0e-6, 2500.0, 273.15, 2773.15, -9999.0)GeoModBox.TimeParameter — Type
TimeParameter()Structure to initialize the time parameters.
The default values are:
const year = 365.25*3600*24 # Seconds per year
tmax = 1000.0 # [ Ma ]
Δfacc = 0.9 # Courant time factor
Δfacd = 0.9 # Diffusion time factor
Δ = 0.0 # Absolute time step
Δc = 0.0 # Courant time step
Δd = 0.0 # Diffusion time stability criterion
itmax = 8000 # Maximum iterationsExamples
julia> T = TimeParameter
TimeParameter(3.15576e7, 1000.0, 0.9, 0.9, 0.0, 0.0, 0.0, 8000)GeoModBox.Dani_Solution_vec! — Method
Dani_Solution_vec!(type,D,M,x,y,rad,mus_i,NC,NV)GeoModBox.Scaling.Constants — Type
Constants()Structure to initialize the scaling constants.
The default values are:
hsc = 0.0
tsc = 0.0
vsc = 0.0
τsc = 0.0
Tsc = 0.0
Qsc = 0.0Examples
julia> S = Constants()
Constants(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)GeoModBox.Scaling.ScaleParameters! — Method
ScaleParameters!(S,M,Δ,T,P,D)Structure to scale the parmeters using the scaling constants.
The scaling laws are:
# Geometry ---
xmin /= S.hsc # Minimum x-coordinate
xmax /= S.hsc # Maximum x-coordinate
ymin /= S.hsc # Minimum y-coordinate
ymax /= S.hsc # Minimum y-coordinate
x /= S.hsc # x-coordinate
y /= S.hsc # y-coordinate
# Time ---
tmax /= S.tsc # Maximum time
Δc /= S.tsc # Courant time step
Δd /= S.tsc # Diffusion time step
Δ /= S.tsc # Total time step
# Temperature etc. ---
Ttop = (P.Ttop - 273.15)/S.Tsc # Temperature at the top
Tbot = (P.Tbot - 273.15)/S.Tsc # Temperature at the bottom
@. D.Q /= S.Qsc
# Viscosity ---
@. D.η /= P.η₀Examples
julia> M = Geometry(0,1000,-1000,0)
Geometry(0.0, 1000.0, -1000.0, 0.0)
julia> P = Physics()
Physics(9.81, 3300.0, 4.125, 1250.0, 2.0e-5, 0.0, 3.947725485e23, 1.0e-6, 2500.0, 273.15, 2773.15, -9999.0)
julia> S = ScalingConstants!(M,P)
Constants(1000.0, 1.0e12, 9.999999999999999e-10, 3.947725485e14, 2500.0, 0.0103125)
julia> NC = (x=100,y=100,)
(x = 100, y = 100)
julia> NV = (
x = NC.x + 1,
y = NC.y + 1,
)
(x = 101, y = 101)
julia> Δ = GridSpacing(
x = (M.xmax - M.xmin)/NC.x,
y = (M.ymax - M.ymin)/NC.y,
)
GridSpacing(10.0, 10.0)
julia> D = DataFields(
Q = zeros(Float64,(NC...)),
T = zeros(Float64,(NC...)),
T_ex = zeros(Float64,(NC.x+2,NC.y+2)),
ρ = ones(Float64,(NC...)),
cp = ones(Float64,(NC...)),
vx = zeros(Float64,(NV.x,NV.y+1)),
vy = zeros(Float64,(NV.x+1,NV.y)),
Pt = zeros(Float64,(NC...)),
vxc = zeros(Float64,(NC...)),
vyc = zeros(Float64,(NC...)),
vc = zeros(Float64,(NC...)),
ΔTtop = zeros(Float64,NC.x),
ΔTbot = zeros(Float64,NC.x),
Tmax = 0.0,
Tmin = 0.0,
Tmean = 0.0,
)
DataFields([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0;;], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0;;], [0.0;;], [0.0;;], [0.0;;], [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [0.0;;], [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0;;], [0.0;;], [0.0;;], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 0.0, 0.0, 0.0)
julia> T = TimeParameter(
tmax = 1000000.0, # [ Ma ]
Δfacc = 1.0, # Courant time factor
Δfacd = 1.0, # Diffusion time factor
itmax = 8000, # Maximum iterations
)
TimeParameter(3.15576e7, 1.0e6, 1.0, 1.0, 0.0, 0.0, 0.0, 8000)
julia> T.tmax = T.tmax*1e6*T.year # [ s ]
3.15576e19
julia> T.Δc = T.Δfacc * minimum((Δ.x,Δ.y)) /
(sqrt(maximum(abs.(D.vx))^2 + maximum(abs.(D.vy))^2))
Inf
julia> T.Δd = T.Δfacd * (1.0 / (2.0 * P.κ *(1.0/Δ.x^2 + 1/Δ.y^2)))
2.5e7
julia> T.Δ = minimum([T.Δd,T.Δc])
2.5e7
ScaleParameters!(S,M,Δ,T,P,D)GeoModBox.Scaling.ScalingConstants! — Method
ScalingConstants!(M,P)Structure to initialize and calcuate the scaling constants in SI units.
The scaling constants are calculated by:
hsc = (ymax-ymin) # Length scale [ m ]
tsc = (ymax-ymin)^2 / κ # Time scale [ s ]
vsc = κ / (ymax-ymin) # Velocity scale [ m/s ]
τsc = (η₀ * κ)/(ymax-ymin) # Stress scale [ Pa ]
Tsc = ΔT # Temperature scale [ K ]
Qsc = (ΔT*κ*ρ₀*cp)/(ymax-ymin)^2 # Heat source scale [ w/m³ ]Examples
julia> M = Geometry(0,1000,-1000,0)
Geometry(0.0, 1000.0, -1000.0, 0.0)
julia> P = Physics()
Physics(9.81, 3300.0, 4.125, 1250.0, 2.0e-5, 0.0, 3.947725485e23, 1.0e-6, 2500.0, 273.15, 2773.15, -9999.0)
julia> ScalingConstants!(M,P)
Constants(1000.0, 1.0e12, 9.999999999999999e-10, 3.947725485e14, 2500.0, 0.0103125)GeoModBox.HeatEquation.OneD.AssembleMatrix1D — Method
AssembleMatrix1D(ρ, cp, k, Δx, Δt, nc, BC, K;C=0 )Function to build the coefficient matrix K for the unknown centroid temperature field in the 2D heat diffusion equation assuming variable thermal parameter and radiogenig heating only.
The coefficient matrix is build using a five-point finite difference stencil, resulting in a five, non-zero diagonal matrix to solve the system of equations.
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
k : Thermal conductivity [ W/m/K ]
Δx : Tuple or structure containing the horizontal grid resolution
Δt : Time step [ s ]
nc : Tuple or structure containing the number of centroids in the horizontal direction
BC : Tuple for the boundary condition
K : Coefficient matrixOptional input values: C : Constant defining the residual for a certain finite difference discretization method, i.e.: C = 0 -> implicit, backward Euler discretization (default) C = 0.5 -> Crank-Nicolson discretization C = 1 -> explicit, forward Euler discretization
GeoModBox.HeatEquation.OneD.AssembleMatrix1Dc! — Method
AssembleMatrix1Dc!( κ, Δx, Δt, nc, BC, K; C )Setup the coefficient matrix for the linear system of equations.
GeoModBox.HeatEquation.OneD.BackwardEuler1Dc! — Method
BackwardEuler1Dc!( implicit, κ, Δx, Δt, nc, BC , K)Solves the onedimensional heat diffusion equation assuming no internal heating and constant thermal parameters using an implicit, backward euler finite difference scheme.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
implicit : Tuple, containing the current temperature array T0 and
the new temperature array T
κ : Diffusivity [ m²/s ]
Δt : Time step [ s ]
nc : Number of central nodes
Δx : Grid spacing [ m ]
BC : Tuple for the boundary condition
K : Coefficient matrix for linear system of equationsGeoModBox.HeatEquation.OneD.CNA1Dc! — Method
CNA1Dc!( cna, κ, Δx, Δt, nc, BC, K1, K2 )Solves the onedimensional heat diffusion equation assuming no internal heating and constant thermal parameters using Crank-Nicolson finite difference scheme.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
cna : Tuple, containing the current temperature array T0 and
the new temperature array T
κ : Diffusivity [ m²/s ]
Δt : Time step [ s ]
nc : Number of central nodes
Δx : Grid spacing [ m ]
BC : Tuple for the boundary condition
K1 : Coefficient matrix for the unknow variables
K2 : Coefficient matrix for the know variablesGeoModBox.HeatEquation.OneD.ComputeResiduals1D! — Method
ComputeResiduals1D!(R, T, T_ex, T0, T_ex0, Q, ∂T, q, ρ, Cp, k, BC, Δx, Δt;C=0)Function to calculate the residual of the one dimensional heat diffusion equation assuming variable thermal parameters and radiogenic heating only. The residual is required for defection correction to solve the linear system of equations.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions. The caluclated residual is used in an iteration loop to calculate the correction term of the initial temperature guess. The coefficient matrix is build via AssembleMatrix1D().
R : 1D array for the residual
T : 1D centroid temperature field of the next time step
T_ex : 1D extended centroid temperature field including the ghost nodes
T_0 : 1D centroid temperature field of the current time step
T_ex0 : 1D extended centroid temperature field of the current time step
Q : Volumetric heat production rate [ W/m³ ]
∂T : Tuple containing the first space derivatives of the temperature
in the horizontal direction
q : Tuple or structure containing the 2D horizontal heat flux
ρ : Density [ kg/m³ ]
Cp : Specific heat capacity [ J/kg/K ]
k : Thermal conductivity [ W/m/K ]
BC : Tuple for the boundary condition
Δx : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
Δt : Time step [ s ]Optional input values: C : Constant defining the residual for a certain finite difference discretization method, i.e.: C = 0 -> implicit, backward Euler discretization (default) C = 0.5 -> Crank-Nicolson discretization C = 1 -> explicit, forward Euler discretization
GeoModBox.HeatEquation.OneD.ComputeResiduals1Dc! — Method
ComputeResiduals1Dc!( cna, κ, Δx, Δt, nc, BC, K1, K2 )Computes the residual of the onedimensional heat diffusion equation assuming no internal heating and constant thermal parameters.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
dc : Tuple, containing the current temperature array T,
the temperature array with ghost nodes T_ex,
the partial derivatives ∂2T∂x2, and the
residual R
κ : Diffusivity [ m²/s ]
Δx : Grid spacing [ m ]
Δt : Time step [ s ]
BC : Tuple for the boundary conditionGeoModBox.HeatEquation.OneD.ForwardEuler1D! — Method
ForwardEuler1D!( explicit, κ, Δx, Δt, nc, BC)Solves the onedimensional heat diffusion equation assuming internal heating and variable thermal parameters using an explicit, forward euler finite difference scheme.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
T : Tuple, containing the regular temperature array T and
array containing the ghost nodes T_ex
Py : Tuple, containing the thermal parameters ρ, k, cp, and H [ W/kg ]
Δt : Time step [ s ]
Δy : Grid spacing [ m ]
nc : Number of central nodes
BC : Tuple for the boundary conditionGeoModBox.HeatEquation.OneD.ForwardEuler1Dc! — Method
ForwardEuler1Dc!( explicit, κ, Δx, Δt, nc, BC)Solves the onedimensional heat diffusion equation assuming no internal heating and constant thermal parameters using an explicit, forward euler finite difference scheme.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
explicit : Tuple, containing the regular temperature array T and
array containing the ghost nodes T_ex
κ : Diffusivity [ m²/s ]
Δt : Time step [ s ]
nc : Number of central nodes
Δx : Grid spacing [ m ]
BC : Tuple for the boundary conditionGeoModBox.HeatEquation.TwoD.ADI2Dc! — Method
ADI2Dc!(T, κ, Δx, Δy, Δt, ρ, cp, NC, BC)Solves the two dimensional heat diffusion equation assuming constant thermal parameters using the alternating-direction implicit finite difference scheme for a linear problem, i.e. a left-matrix division.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Boundary conditions are directly implemented within the system of equations.
T : Tuple, containing the regular temperature array T,
the extended array containing the ghost nodes T_ex, and
the volumentric heat production rate Q
κ : Diffusivity [ m²/s ]
Δx : Horizontal grid spacing [ m ]
Δy : Vertical grid spacing [ m ]
Δt : Time step [ s ]
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
NC : Tuple containing the numer of centroids in x- and y-direction
BC : Tuple for the boundary conditionGeoModBox.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)Function to build the coefficient matrix K for the unknown centroid temperature field in the 2D heat diffusion equation assuming variable thermal parameter and radiogenig heating only.
The coefficient matrix is build using a five-point finite difference stencil, resulting in a five, non-zero diagonal matrix to solve the system of equations.
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
k : Thermal conductivity [ W/m/K ]
BC : Tuple for the boundary condition
Num : Tuple or structure containing the global numbering of the centroids
nc : Tuple or structure containing the number of centroids in the horizontal and vertical direction
Δ : Tuple or structure containint the horizontal and vertical grid resolution
Δt : Time step [ s ]Optional input values: C : Constant defining the residual for a certain finite difference discretization method, i.e.: C = 0 -> implicit, backward Euler discretization (default) C = 0.5 -> Crank-Nicolson discretization C = 1 -> explicit, forward Euler discretization
GeoModBox.HeatEquation.TwoD.AssembleMatrix2Dc — Method
AssembleMatrix2Dc( κ, BC, Num, nc, Δ, Δt )Function to build the coefficient matrix K for the unknown centroid temperature field in the 2D heat diffusion equation assuming constant thermal parameter and radiogenig heating only.
The coefficient matrix is build using a five-point finite difference stencil, resulting in a five, non-zero diagonal matrix to solve the system of equations.
κ : Diffusivity [ m²/s ]
BC : Tuple for the boundary condition
Num : Tuple or structure containing the global numbering of the centroids
nc : Tuple or structure containing the number of centroids in the horizontal and vertical direction
Δ : Tuple or structure containint the horizontal and vertical grid resolution
Δt : Time step [ s ]GeoModBox.HeatEquation.TwoD.BackwardEuler2Dc! — Method
BackwardEuler2Dc!(D, κ, Δx, Δy, Δt, ρ, cp, NC, BC, rhs, K, Num)Solves the two dimensional heat diffusion equation assuming constant thermal parameters using an implicit, backward Euler finite difference scheme for a linear problem, i.e. a left-matrix division.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Boundary conditions are directly implemented within the system of equations.
D : Tuple, containing the regular temperature array T and
the extended array containing the ghost nodes T_ex
κ : Diffusivity [ m²/s ]
Δx : Horizontal grid spacing [ m ]
Δy : Vertical grid spacing [ m ]
Δt : Time step [ s ]
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
NC : Tuple containing the numer of centroids in x- and y-direction
BC : Tuple for the boundary condition
rhs : Known right-hand side vector
K : Coefficient matrix in sparse format
Num : Tuple or structure containing the global centroid numberingGeoModBox.HeatEquation.TwoD.CNA2Dc! — Method
CNA2Dc!(D, κ, Δx, Δy, Δt, ρ, cp, NC, BC, rhs, K1, K2, Num)Solves the two dimensional heat diffusion equation assuming constant thermal parameters using the Crank-Nicolson finite difference scheme for a linear problem, i.e. a left-matrix division.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Boundary conditions are directly implemented within the system of equations.
D : Tuple, containing the regular temperature array T and
the extended array containing the ghost nodes T_ex
κ : Diffusivity [ m²/s ]
Δx : Horizontal grid spacing [ m ]
Δy : Vertical grid spacing [ m ]
Δt : Time step [ s ]
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
NC : Tuple containing the numer of centroids in x- and y-direction
BC : Tuple for the boundary condition
rhs : Known right-hand side vector
K1 : Coefficient matrix in sparse format for the unknown temperature
K1 : Coefficient matrix in sparse format for the known temperature
Num : Tuple or structure containing the global centroid numberingGeoModBox.HeatEquation.TwoD.ComputeResiduals2D! — Method
ComputeResiduals2D!(R, T, T_ex, T0, T_ex0, Q, ∂T, q, ρ, Cp, k, BC, Δ, Δt;C=0)Function to calculate the residual of the two dimensional heat diffusion equation assuming variable thermal parameters and radiogenic heating only. The residual is required for defection correction to solve the linear system of equations.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions. The caluclated residual is used in an iteration loop to calculate the correction term of the initial temperature guess. The coefficient matrix is build via AssembleMatrix2D().
R : 2D array for the residual
T : 2D centroid temperature field of the next time step
T_ex : 2D extended centroid temperature field including the ghost nodes
T_0 : 2D centroid temperature field of the current time step
T_ex0 : 2D extended centroid temperature field of the current time step
Q : Volumetric heat production rate [ W/m³ ]
∂T : Tuple containing the first space derivatives of the temperature
in the horizontal and vertical direction
q : Tuple or structure containing the 2D horizontal and vertical heat flux
ρ : Density [ kg/m³ ]
Cp : Specific heat capacity [ J/kg/K ]
k : Thermal conductivity [ W/m/K ]
BC : Tuple for the boundary condition
Δ : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
Δt : Time step [ s ]Optional input values: C : Constant defining the residual for a certain finite difference discretization method, i.e.: C = 0 -> implicit, backward Euler discretization (default) C = 0.5 -> Crank-Nicolson discretization C = 1 -> explicit, forward Euler discretization
GeoModBox.HeatEquation.TwoD.ComputeResiduals2Dc! — Method
ComputeResiduals2Dc!(R, T, T_ex, T0, T_ex0, ∂2T, Q, κ, BC, Δ, Δt;C=0)Function to calculate the residual of the two dimensional heat diffusion equation assuming constant thermal parameters and radiogenic heating only. The residual is required for defection correction to solve the linear system of equations.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions. The caluclated residual is used in an iteration loop to calculate the correction term of the initial temperature guess. The coefficient matrix is build via AssembleMatrix2Dc().
R : 2D array for the residual
T : 2D centroid temperature field of the next time step
T_ex : 2D extended centroid temperature field including the ghost nodes
T_0 : 2D centroid temperature field of the current time step
T_ex0 : 2D extended centroid temperature field of the current time step
∂2T : Tuple containing the second space derivatives of the temperature
in the horizontal and vertical direction
Q : Volumetric heat production rate [ W/m³ ]
κ : Diffusivity [ m²/s ]
BC : Tuple for the boundary condition
Δ : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
Δt : Time step [ s ]Optional input values: C : Constant defining the residual for a certain finite difference discretization method, i.e.: C = 0 -> implicit, backward Euler discretization (default) C = 0.5 -> Crank-Nicolson discretization C = 1 -> explicit, forward Euler discretization
GeoModBox.HeatEquation.TwoD.ForwardEuler2Dc! — Method
ForwardEuler2Dc!(D, κ, Δx, Δy, Δt, ρ, cp, NC, BC)Solves the two dimensional heat diffusion equation assuming constant thermal parameters using an explicit, forward Euler finite difference scheme.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Using central temperature nodes requires external ghost nodes, which are used to define the boundary conditions.
D : Tuple, containing the regular temperature array T and
the extended array containing the ghost nodes T_ex
κ : Diffusivity [ m²/s ]
Δx : Horizontal grid spacing [ m ]
Δy : Vertical grid spacing [ m ]
Δt : Time step [ s ]
ρ : Density [ kg/m³ ]
cp : Specific heat capacity [ J/kg/K ]
NC : Tuple containing the numer of centroids in x- and y-direction
BC : Tuple for the boundary conditionGeoModBox.HeatEquation.TwoD.Poisson2D! — Method
Poisson2D!( T, Q, kx, ky, Δx, Δy, NC, BC, K, rhs, Num )Function to solve the two dimensional poisson equation of the steady state heat diffusion equation assuming variable thermal parameters and radiogenic heat sources only.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Boundary conditions are directly implemented within the system of equations.
T : Tuple or structure containing the centroid temperature field.
Q : Tuple or structure containing the centroid volumentric heat prodcution rate [ W/m³ ].
kx : Thermal conductivity in the horizontal direction [ W/m/K ].
ky : Thermal conductivity in the vertical direction [ W/m/K ].
Δx : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
Δy : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
NC : Tuple or structure containing the number of horizontal and vertical
centroids.
BC : Tuple or structure containint the boundary conditions.
K : Coefficient matrix in sparse format
rhs : Known right-hand side vector
Num : Tuple or structure containing the global centroid numberingGeoModBox.HeatEquation.TwoD.Poisson2Dc! — Method
Poisson2Dc!(D,NC,P,BC,Δ,K,rhs,Num)Function to solve the two dimensional poisson equation of the steady state heat diffusion equation assuming constant thermal parameters and radiogenic heat sources only.
The temperature is defined on central nodes and the heat flux on the vertices. Boundary conditions are currently limited to Dirichlet and Neumann. Boundary conditions are directly implemented within the system of equations.
D : Tuple or structure containing the centroid temperature field
and volumetic heat production rate.
NC : Tuple or structure containing the number of horizontal and vertical
centroids.
P : Tuple or structure containing the thermal conductivity [ W/m/K ].
BC : Tuple or structure containint the boundary conditions.
Δ : Tuple or strucutre containint the horizontal and vertical grid resolution [ m ]
K : Coefficient matrix in sparse format
rhs : Known right-hand side vector
Num : Tuple or structure containing the global centroid numberingGeoModBox.AdvectionEquation.OneD.RK4O1D! — Method
RK4O1D!( x, Δt, vx, xmin, xmax )Function to advect a tracer in 1-D for one time step using fourth order Runge-Kutta.
The final position of the tracer is defined by:
x .+ = (x1 + 2.0 * (x2 + x3) + x4) / 6.0where
x1 = Δt * vx
x2 = Δt * vx
x3 = Δt * vx
x4 = Δt * vxExamples
julia> xminx,xmax = 0,100
(0, 100)
julia> x = collect(xmin:0.1:xmax)
1001-element Vector{Float64}:
0.0
0.1
0.2
0.3
0.4
⋮
99.7
99.8
99.9
100.0
julia> x = RK4O1D!(x,0.1,1,xmin,xmax)GeoModBox.AdvectionEquation.OneD.lax1D! — Method
lax1D!( A, vx, Δt, Δx )Function to advect a property using the Lax-Friedrich method in one dimension.
A : 1-D array of the advected property
vx : Horizontal velocity
Δt : Time step
Δx : Grid spacingGeoModBox.AdvectionEquation.OneD.semilag1D! — Method
semilag1D!( A, xc, vx, Δt, Δx )Function to advect a property using the semi-lagrangian method in one dimension.
A : 1-D array of the advected property
xc : x-coordinates of each centroid
vx : Horizontal velocity
Δt : Time step
Δx : Horizontal velocityGeoModBox.AdvectionEquation.OneD.slf1D! — Method
slf1D!( A, Aold2, vx, Δt, Δx )Function to advect a property using the staggered leaped frog method in one dimension.
A : 1-D array of the advected property
Aold2 : 1-D array of the advected property at the previous time step
vx : Horizontal velocity
Δt : Time step
Δx : Grid spacingGeoModBox.AdvectionEquation.OneD.upwind1D! — Method
upwind1D!( A, vx, Δt, Δx )Function to advect a tracer using the upwind method in one dimension.
A : 1-D Array of the advected property
vx : Horizontal velocity
Δt : Time step
Δx : Grid spacingGeoModBox.AdvectionEquation.TwoD.semilagc2D! — Method
semilagc2D!()Function to advect a property using the semi-lagrange method in two dimensions.
P : 2-D array of the advecte property on the centroids
P_ex : 2-D array including the ghost nodes
vxc : 2-D array of the horizontal velocity on the centroids
vyc : 2-D array of the vertical velocity on the centroids
vxo : 2-D array of the horizontal velocity of the previous time step on the centroids
vyo : 2-D array of the vertical velocity of the previous time step on the centroids
x : Structure or Tuple containing the 2-D horizontal coordinates of the centroids
y : Structure or Tuple containing the 2-D vertical coordinates of the centroids
Δt : Time stepThe function calculates a mean velocity of the previous and current velocity field, if needed, and advects the property using a mid-point iteration scheme.
GeoModBox.AdvectionEquation.TwoD.slfc2D! — Method
slfc2D!Function to advect a property using the staggered-leaped frog method in two dimensions.
P : 2-D array of the advecte property on the centroids
P_ex : 2-D array including the ghost nodes
P_exo : 2-D array of the previous time step including the ghost nodes
vxc : 2-D array of the horizontal velocity on the centroids
vyc : 2-D array of the vertical velocity on the centroids
NC : Structure or Tuple containing the number of centroids
Δt : Time step
Δx : Horizontal grid spacing
Δy : Vertical grid spacingGeoModBox.AdvectionEquation.TwoD.upwindc2D! — Method
upwindc2D!(P,P_ex,vxc,vyc,NC,Δt,Δx,Δy)Function to advect a property using the upwind method in two dimensions.
P : 2-D array of the advected property on the centroids
P_ex : 2-D array including the ghost nodes
vxc : 2-D array of the horizontal velocity on the centroids
vyc : 2-D array of the vertical velocity on the centroids
NC : Structure or Tuple containing the number of centroids
Δt : Time step
Δx : Horizontal grid spacing
Δy : Vertical grid spacingGeoModBox.InitialCondition.IniPhase! — Method
IniPhase!(type,D,M,x,y,NC;phase=0)Function to setup an initial phase field for a two dimensional problem. The phase ID is defined on a the centroids of a regular finite difference grid. The function is outdated. Better to use the tracers to define the phase distributions and to inerpolate the phase ID from the tracers to the grid.
type : Parameter defining the phase type
D : Structure or tuple containing the field arrays
M : Structure or tuple containing the geometry
x : Structure or tuple containing the x-coordinates
y : Structure or tuple containing the y-coordinates
NC : Structure or tuple containing the centroids parameterCertain default values can be modified as well:
phase : Vector containing the phase ID numbers (e.g. phase=[0,1])Currently, only one initial phase configuration is available:
1) A rectangular shaped anomaly (block)Example:
IniPhase!(:block,D,M,x,y,NC;phase=[0 1])GeoModBox.InitialCondition.IniTemperature! — Method
IniTemperature!(type,M,NC,D,x,y;Tb=1000,Ta=1200,Ampl=200,σ=0.05)Function to setup an initial temperature condition for a two dimensional problem. The temperature is defined on the centroids of a regular finite difference grid.
type : Parameter defining the anomaly type
M : Structure or tuple containing the geometry
NC : Structure or tuple containing the centroids parameter
D : Structure or tuple containing the field arrays
x : Structure or tuple containing the x-coordinates
y : Structure or tuple containing the y-coordinatesCertain default values can be modified as well:
Tb : Scalar value for the background (or top) temperature
Ta : Scalar value for the maximum (bottom or anomaly) temperature
σ : Width of the Gaussian temperature anomalyPossible anomaly types are:
1) A circular anomaly with a constant background (circle)
2) A Gaussian anomaly (gaussian)
3) A rectangular shaped anomaly with a constant background (block)
4) A linear increasing temperature with depth (linear)
5) A linear increasing temperature with depth including an elliptical anomaly (lineara)Example:
IniTemperature!(:circle,M,NC,D,x,y;Tb=1200,Ta=0)GeoModBox.InitialCondition.IniVelocity! — Method
IniVelocity!(type,D,NC,Δ,M,x,y)Function to setup an initial velocity field for a two dimensional problem. The velocity is defined on a staggered finite difference grid, that is inbetween the vertices of the regular grid.
type : Parameter defining the velocity type
D : Structure or tuple containing the field arrays
BC : Structure or tuple containing the velocity boundary conditions
NV : Structure or tuple containing the vertices parameter
Δ : Structure or tuple containing the grid resolution
M : Structure or tuple containing the geometry
x : Structure or tuple containing the x-coordinates
y : tructure or tuple containing the y-coordinatesCertain default values can be modified as well:
ε : Background strain rate for pure shear or simple shearThe following velocity configurations are currently supported:
1) A rigid-body rotation (RigidBody)
2) A shear cell (ShearCell)
3) Simple Shear (SimpleShear)
4) Pure Shear (PureShear)For Pure Shear and Simple Shear the boundary velocity values are updated accordingly.
Example:
IniVelocity!(:PureShear,D,VBC,NV,Δ,M,x,y;ε=1e-15)GeoModBox.Tracers.OneD.Itp1D_Centers2Markers! — Method
Itp1D_Centers2Markers!( Tm, xm, Tc, xc, Δx, xmin )GeoModBox.Tracers.OneD.Itp1D_Markers2Centers! — Method
Itp1D_Markers2Centers!( Tc, xc, Tm, xm, dx, xmin )GeoModBox.Tracers.TwoD.Markers — Type
Ma = Markers(xmi,ymi,phase)Function to allocate the marker arrays to advect the phase ID.
xmi : Initial horizontal marker coordinates
ymi : Initial vertical marker coordinates
phase : Initial marker phase IDExample: Ma = TMarkers( vec(xmi), vec(ymi), zeros(Int64, nmark) )
GeoModBox.Tracers.TwoD.TMarkers — Type
Ma = TMarkers(xmi,ymi,T,phase)Function to allocate the marker arrays to advect the absolute temperature.
xmi : Initial horizontal marker coordinates
ymi : Initial vertical marker coordinates
T : Initial marker temperature
phase : Initial marker phase IDExample: Ma = TMarkers( vec(xmi), vec(ymi), zeros(Float64, nmark), zeros(Int64, nmark) )
GeoModBox.Tracers.TwoD.AdvectTracer2D — Method
AdvectTracer2D(Ma,nmark,D,x,y,dt,Δ,NC,rkw,rkv,style)Function to advect tracer in 2D using forth order Runge-Kutta.
Ma : Tuple or structure containing the marker information
nmark : Total number of tracers
D : Tuple or structure containing the velocities
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
dt : Time step
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroids
rkw : Runge-Kutta coefficients for velocity averaging
rkv : Runge-Kutta coefficients for time steppingCertain default values can be modified as well:
style : Defines the calculation of the velocity, default 1, (1,2,3)GeoModBox.Tracers.TwoD.CountMPC — Method
CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,NV,it)Function to count the marker per cell.
Ma : Tuple or structure containing the marker properties
nmark : Total number of markers
MPC : Tuple or structure containing marker per node per thread information
M : Tuple or structure containing the model geometry
x : Tuple or sturcture containing the x-coordinates
y : Tuple or sturcture containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroids
NV : Tuple or structure containing the number of verticesExample:
CountMPC(Ma,nmark,MPC,M,x,y,Δ,NC,NV,)GeoModBox.Tracers.TwoD.FromCtoM — Method
FromCtoM(Prop, Ma, x, y, Δ, NC)Function to interpolate the centroid property onto the tracers.
Prop : Property array on the centroids
Ma : Tuple or structure containing marker information
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroidsExample: for k = 1:nmark Ma.T[k] = FromCtoM(D.T_ex, k, Ma, x, y, Δ, NC) end
GeoModBox.Tracers.TwoD.IniTracer2D — Method
IniTracer2D(Aparam,nmx,nmy,Δ,M,NC,noise,ini,phase;λ,δA,ellA,ellB,α)Function to setup the initial tracer distribution.
Aparam : Defines if temperature (thermal) or phase (phase) is advected
nmx : Number of horizontal tracers per cell
nmy : Number of vertical tracers per cell
Δ : Structure or tuple containing the grid resolution
M : Structure or tuple containing the geometry
NC : Structure or tuple containing centroids parameter
noise : add noise; 1 - yes, 0 - no
ini : Initial phase distribution (e.g., block)
phase : Vector with phase IDs, (e.g. [0,1])Certain default values can be modified as well:
λ : Wavelength [m] for a cosine perturbation, e.g. for the RTI
δA : Amplitude [m] of the perturbation
ellA : Major half axis [m] of the elliptical inclusion
ellB : Minor half axis [m] of the elliptical inclusion
α : Rotation angle [°] of the elliptical inclusionExample: 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)Function to interpolate the marker property to the centroids using a weighted bilinear interpolation scheme.
Ma : Tuple or structure containing the marker information
nmark : Total number of markers
PC_th : Property on the centroids per thread
PC : Property on the centroids
weight_th : Weight on the centroids per thread
weight : Weight on the centroids
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
param : Switch to define the advected property: :thermal or :phase
param2 : Array containing the value for each phase, e.g. [ρ₁ ρ₂]Certain default values can be modified as well:
avgm : Averaging scheme: arithmetic (arith), geometric (geom), or harmonic (harm)Example:
Markers2Cells(Ma,nmark,MAVG.PC_th,D.ρ_ex,MAVG.wte_th,D.wte,x,y,Δ,Aparam,ρ)
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)Function to interpolate the marker property to the vertices using a weighted bilinear interpolation scheme.
Ma : Tuple or structure containing the marker information
nmark : Total number of markers
PG_th : Property on the vertices per thread
PG : Property on the vertices
weight_th : Weight on the centroids per thread
weight : Weight on the centroids
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
param : Switch to define the advected property: :thermal or :phase
param2 : Array containing the value for each phase, e.g. [ρ₁ ρ₂]Certain default values can be modified as well:
avgm : Averaging scheme: arithmetic (arith), geometric (geom), or harmonic (harm)Example:
Markers2Vertices(Ma,nmark,MAVG.PV_th,D.ηv,MAVG.wtv_th,D.wtv,x,y,Δ,Aparam,η)GeoModBox.Tracers.TwoD.VxFromVxNodes — Method
VxFromVxNodes(Vx, k, Ma, x, y, Δ, NC, new)Function to calculate horizontal marker velocities from the vx staggered nodes. Modified from https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
Vx : Horizontal velocity field
k : Marker number
Ma : Tuple or structure containing marker information
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroids
new : Switch to activate updated velocity calculationExample:
vxm = VxFromVxNodes(D.vx, k, Ma, x, y, Δ, NC, 0)GeoModBox.Tracers.TwoD.VxVyFromPrNodes — Method
VxVyFromPrNodes(Vxp ,Vyp, k, Ma, x, y, Δ, NC )Function to calculate horizontal and vertical marker velocities from the centroids. Modified from https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
Vxp : Horizontal velocity field at the centroids
Vyp : Vertical velocity field at the centroids
k : Marker number
Ma : Tuple or structure containing marker information
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroidsExample:
vxp, vyp = VxVyFromPrNodes(D.vxc ,D.vyc, k, Ma, x, y, Δ, NC)GeoModBox.Tracers.TwoD.VyFromVyNodes — Method
VyFromVyNodes(Vy, k, Ma, x, y, Δ, NC, new)Function to calculate vertical marker velocities from the vy staggered nodes. Modified from https://github.com/tduretz/M2DptJulia/blob/master/Markers2D/MainTarasv6Hackathon.jl
Vy : Vertical velocity field
k : Marker number
Ma : Tuple or structure containing marker information
x : Tuple or structure containing the x-coordinates
y : Tuple or structure containing the y-coordinates
Δ : Tuple or structure containing the grid resolution
NC : Tuple or structure containing the number of centroids
new : Switch to activate updated velocity calculationExample:
vym = VyFromVyNodes(D.vy, k, Ma, x, y, Δ, NC, 0)GeoModBox.MomentumEquation.OneD.AssembleStokesMatrix1D — Method
AssembleStokesMatrix1D(nc, η, Δy, BC, K)GeoModBox.MomentumEquation.OneD.ComputeStokesResiduals1D! — Method
ComputeStokesResiduals1D!( D, ∂P∂x, Δy, BC)GeoModBox.MomentumEquation.OneD.Stokes_1D_direct — Method
Stokes_1D_direct(vₓ,η,Δy,nc,BC,K,rhs)GeoModBox.MomentumEquation.TwoD.Assembly — Method
Assembly(NC, NV, Δ, ηc, ηv, BC, Num)GeoModBox.MomentumEquation.TwoD.Assemblyc — Method
Assemblyc(NC, NV, Δ, η, BC, Num)GeoModBox.MomentumEquation.TwoD.Residuals2D! — Method
Residuals2D!(D,BC,ε,τ,divV,Δ,ηc,ηv,g,Fm,FPt)GeoModBox.MomentumEquation.TwoD.Residuals2Dc! — Method
Residuals2Dc!(D,BC,ε,τ,divV,Δ,η,g,Fm,FPt)GeoModBox.MomentumEquation.TwoD.updaterhs — Method
updaterhs(NC, NV, Δ, ηc, ηv, ρ, g, BC, Num)GeoModBox.MomentumEquation.TwoD.updaterhsc — Method
updaterhsc(NC, NV, Δ, η, ρ, g, BC, Num)