GeoParams.jl API

Features

  • [x] Functors to available creep laws
  • [x] Elasticity struct
  • [x] Interoperability with ParallelStencil.jl
  • [x] compute_viscosity() agnostic to input parameters
  • [x] composable custom creep law

The following API is created using the docstrings avaliable within the ParallelStencil.jl package, and is included for convinience here for look-up purpose. Please refer to the official repository and ask the creator of the packages if anything is unclear.

<!–

GeoParams.ConstTempType
ConstTemp(T=1000)

Sets a constant temperature inside the box

Parameters:

  • T : the value
GeoParams.HalfspaceCoolTempType
HalfspaceCoolTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)

Sets a halfspace temperature structure in plate

Parameters:

  • Tsurface: surface temperature [C]
  • Tmantle: mantle temperature [C]
  • Age: Thermal Age of plate [Myrs]
  • Adiabat: Mantle Adiabat [K/km]
GeoParams.LinTempType
LinTemp(Ttop=0, Tbot=1000)

Sets a linear temperature structure from top to bottom

Parameters:

  • Ttop: the value @ the top
  • Tbot: the value @ the bottom
GeoParams.CompTempStructMethod
CompTempStruct(Z, s::AbstractTempStruct)

Returns a temperature vector that matches the coordinate vector and temperature structure that were passed in

Parameters:

  • Z: vector with depth coordinates
  • s: Temperature structure (CostTemp, LinTemp, HalfspaceCoolTemp)
GeoParams.LithPresMethod
LithPres(MatParam, Phases, ρ, T, dz, g)

Iteratively solves a 1D lithostatic pressure profile (compatible with temperature- and pressure-dependent densities)

Parameters:

  • MatParam: a tuple of materials (including the following properties: Phase, Density)
  • Phases: vector with the distribtion of phases
  • ρ: density vector for initial guess(can be zeros)
  • T: temperature vector
  • dz: grid spacing
  • g: gravitational accelaration
GeoParams.StrengthEnvelopeCompMethod
StrengthEnvelopeComp(MatParam::NTuple{N, AbstractMaterialParamsStruct}, Thickness::Vector{U}, TempType::AbstractTempStruct=LinTemp(0C, 800C), ε=1e-15/s, nz::Int64=101) where {N, U}

Calculates a 1D strength envelope. Pressure used for Drucker Prager plasticity is lithostatic. To visualize the results in a GUI, use StrengthEnvelopePlot.

Parameters:

  • MatParam: a tuple of materials (including the following properties: Phase, Density, CreepLaws, Plasticity)
  • Thickness: a vector listing the thicknesses of the respective layers (should carry units)
  • TempType: the type of temperature profile (ConstTemp(), LinTemp(), HalfspaceCoolTemp())
  • ε: background strainrate
  • nz: optional argument controlling the number of points along the profile (default = 101)

Example:

julia> using GLMakie
julia> MatParam = (SetMaterialParams(Name="UC", Phase=1, Density=ConstantDensity(ρ=2700kg/m^3), CreepLaws = SetDislocationCreep("Wet Quartzite | Ueda et al. (2008)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)),
                   SetMaterialParams(Name="MC", Phase=2, Density=Density=ConstantDensity(ρ=2900kg/m^3), CreepLaws = SetDislocationCreep("Plagioclase An75 | Ji and Zhao (1993)"), Plasticity = DruckerPrager(ϕ=20.0, C=10MPa)),
                   SetMaterialParams(Name="LC", Phase=3, Density=PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-10/Pa), CreepLaws = SetDislocationCreep("Maryland strong diabase | Mackwell et al. (1998)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)));
julia> Thickness = [15,10,15]*km;

julia> StrengthEnvelopeComp(MatParam, Thickness, LinTemp(), ε=1e-15/s)
GeoParams.StrengthEnvelopePlotMethod
StrengthEnvelopePlot(MatParam, Thickness; TempType, nz)

Requires GLMakie

Creates a GUI that plots a 1D strength envelope. In the GUI, temperature profile and strain rate can be adjusted. The Drucker-Prager plasticity uses lithostatic pressure.

Parameters:

  • MatParam: a tuple of materials (including the following properties: Phase, Density, CreepLaws, Plasticity)
  • Thickness: a vector listing the thicknesses of the respective layers (should carry units)
  • TempType: the type of temperature profile (LinTemp=default, HalfspaceCoolTemp, ConstTemp)
  • nz: optional argument controlling the number of points along the profile (default = 101)

Example:

julia> using GLMakie
julia> MatParam = (SetMaterialParams(Name="UC", Phase=1, Density=ConstantDensity(ρ=2700kg/m^3), CreepLaws = SetDislocationCreep("Wet Quartzite | Ueda et al. (2008)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)),
                   SetMaterialParams(Name="MC", Phase=2, Density=Density=ConstantDensity(ρ=2900kg/m^3), CreepLaws = SetDislocationCreep("Plagioclase An75 | Ji and Zhao (1993)"), Plasticity = DruckerPrager(ϕ=20.0, C=10MPa)),
                   SetMaterialParams(Name="LC", Phase=3, Density=PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-10/Pa), CreepLaws = SetDislocationCreep("Maryland strong diabase | Mackwell et al. (1998)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)));
julia> Thickness = [15,10,15]*km;

julia> StrengthEnvelopePlot(MatParam, Thickness, LinTemp())
GeoParams.compute_p_τij!Method
compute_p_τij!(Txx, Tyy, Txy, Tii, η_vep, P, Exx, Eyy, Exy, P_o, Txx_o, Tyy_o, Txy_o, phase, MatParam, dt)

Computes 2D deviatoric stress components (Txx,Tyy,Txy) and pressure P given deviatoric strainrate components (Exx,Eyy,Exy) and old deviatoric stresses (Txx_o, Tyy_o, Txy_o) and old pressure P_o (only used for viscoelastic cases). Also returned are Tii (second invariant of the deviatoric stress tensor), and η_vep the viscoelastoplastic effective viscosity. Also required as input is MatParam, the material parameters for every phase and phase, an integer array of size(Exx) that indicates the phase of every point.

This function assumes that strainrate points are collocated and that Exx,Eyy,Exy are at the same points.

GeoParams.compute_p_τijMethod
P,τij, τII = compute_p_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple{N2,T}, P_old::T, args, τij_old::NTuple{3,T}, phase::I)

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a collocated grid case with a single phase phase.

GeoParams.compute_p_τijMethod
τij, τII, η_eff = compute_p_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple, args, τij_old::NTuple, phases::NTuple)

This computes pressure p and deviatoric stress components τij, their second invariant τII, and effective viscosity η_eff for given deviatoric strainrates εij, old stresses τij_old, phases (integer) for every point and arguments args. It handles staggered grids of various taste

GeoParams.compute_p_τijMethod
p,τij,τII = compute_p_τij(v, εij::NTuple{n,T}, P_old::T, args,  τij_old::NTuple{N,T})

Computes pressure p and deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a case that all points are collocated and we have a single phase.

GeoParams.compute_p_τijMethod
p, τij, τII = compute_p_τij(v, εij::NTuple{N,Union{T,NTuple{4,T}}}, P_old::T, args, τij_old::NTuple{3,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a staggered grid case with a single phase.

GeoParams.compute_p_τij_stagcenter!Method
compute_τij_stagcenter!(Txx, Tyy, Txy, Tii, η_vep,  P, Exx, Eyy, Exyv, P_o, Txx_o, Tyy_o, Txyv_o, phase_center, phase_vertex, MatParam, dt)

Updates deviatoric stresses on a staggered grid in a centered based manner (averages strainrates/old stresses from vertices -> centers). Take care of the sizes of the input matrixes:

  • Txx,Tyy,Txy,Tii,P: 2D matrixes of size (nx,ny) which describe updated deviatoric stress components (x,y,xy, 2nd invariant) at the center point
  • η_vep: viscoelastoplastic viscosity @ center (nx,ny)
  • Exx,Eyy,P_o: deviatoric strain rate components & old pressure @ center
  • Exy: shear strain rate @ vertices (nx+1,ny+1)
  • Txx_o, Tyy_o: deviatoric stress normal components of last timestep at center (nx,ny)
  • Txy_o: deviatoric stress shear components at vertices (nx+1,ny+1)
  • phase_center: integer array with phase at center points (nx,ny)
  • phase_vertex: integer array with phase at vertex points (nx+1,ny+1)
  • MatParam: material parameter array that includes a CompositeRheology field which specifies the rheology for different phases
  • dt: timestep
GeoParams.compute_τijMethod
τij, τII = compute_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple{N2,T}, args, τij_old::NTuple{3,T}, phase::I)

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a collocated grid case with a single phase phase.

GeoParams.compute_τijMethod
τij, τII, η_eff = compute_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple, args, τij_old::NTuple, phases::NTuple)

This computes deviatoric stress components τij, their second invariant τII, and effective viscosity η_eff for given deviatoric strainrates εij, old stresses τij_old, phases (integer) for every point and arguments args. This handles various staggered grid arrangements; if staggered components are given as NTuple{4,T}, they will be averaged. Note that the phase of all staggered points should be the same.

GeoParams.compute_τijMethod
τij,τII = compute_τij(v, εij::NTuple{n,T}, args, τij_old::NTuple{N,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a case that all points are collocated and we have a single phase.

GeoParams.compute_τijMethod
τij, τII = compute_τij(v, εij::NTuple{N,Union{T,NTuple{4,T}}}, args, τij_old::NTuple{3,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a staggered grid case with a single phase.

GeoParams.compute_τij_stagcenter!Method
compute_τij_stagcenter!(Txx, Tyy, Txy, Tii, η_vep,  Exx, Eyy, Exyv, P, Txx_o, Tyy_o, Txyv_o, phase_center, phase_vertex, MatParam, dt)

Updates deviatoric stresses on a staggered grid in a centered based manner (averages strainrates/old stresses from vertices -> centers). Take care of the sizes of the input matrixes:

  • Txx,Tyy,Txy,Tii: 2D matrixes of size (nx,ny) which describe updated deviatoric stress components (x,y,xy, 2nd invariant) at the center point
  • η_vep: viscoelastoplastic viscosity @ center (nx,ny)
  • Exx,Eyy,P: deviatoric strain rate components & pressure @ center
  • Exy: shear strain rate @ vertices (nx+1,ny+1)
  • Txx_o, Tyy_o: deviatoric stress normal components of last timestep at center (nx,ny)
  • Txy_o: deviatoric stress shear components at vertices (nx+1,ny+1)
  • phase_center: integer array with phase at center points (nx,ny)
  • phase_vertex: integer array with phase at vertex points (nx+1,ny+1)
  • MatParam: material parameter array that includes a CompositeRheology field which specifies the rheology for different phases
  • dt: timestep
GeoParams.second_invariant_staggeredMethod
second_invariant_staggered(Axx::NTuple{4,T}, Ayy::NTuple{4,T}, Azz::NTuple{4,T}, Aij::NTuple{3,T}) where {T}

Computes the second invariant of the 2D tensor A when its diagonal components need to be maped from cell center to cell vertex. Axx, Ayy, and Azz are tuples containinig the diagonal terms of A at the cell centers around the i-th vertex., and Aij is a tuple that contains the off-diagonal components at the i-th vertex.

GeoParams.second_invariant_staggeredMethod
second_invariant_staggered(Axx::NTuple{4,T}, Ayy::NTuple{4,T}, Axy::Number) where {T}

Computes the second invariant of the 2D tensor A when its diagonal components need to be maped from cell center to cell vertex. Axx, and Ayy are tuples containinig the diagonal terms of A at the cell centers around the i-th vertex., and Axy is the xy component at the i-th vertex.

GeoParams.second_invariant_staggeredMethod
second_invariant_staggered(Aii::NTuple{3,T}, Ayz::NTuple{4,T}, Axz::NTuple{4,T}, Axy::NTuple{4,T}) where {T}

Computes the second invariant of the 2D tensor A when its off-diagonal components need to be maped from cell center to cell vertex. Aii is a tuple containinig the diagonal terms of A at the i-th vertex, and Ayz, Axz, and Axy are tuples that contain the off-diagonal components of the tensor at the cell centers around the i-th vertex.

GeoParams.second_invariant_staggeredMethod
second_invariant_staggered(Aii::NTuple{2,T}, Axy::NTuple{4,T}) where {T}

Computes the second invariant of the 2D tensor A when its off-diagonal components need to be maped from cell center to cell vertex. Aii is a tuple containinig the diagonal terms of A at the i-th vertex, and Axy is a tuple that contains A_xy at the cell centers around the i-th vertex.

GeoParams.time_p_τII_0D!Method
time_p_τII_0D!(P_vec::Vector{T}, τ_vec::Vector{T}, v::CompositeRheology, εII_vec::Vector{T}, εvol_vec::Vector{T}, args, t_vec::AbstractVector{T}) where {T}

Computes stress-time evolution for a 0D (homogeneous) setup with given shear and volumetric strainrate vectors (which can vary with time).

GeoParams.time_p_τII_0DMethod
t_vec, P_vec, τ_vec = time_p_τII_0D(v::CompositeRheology, εII::Number, εvol::Number, args; t=(0.,100.), τ0=0., nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, shear strainrate εII and volumetric strainrate εvol, as well as rheology arguments args. The initial stress τ0, the time range t and the number of timesteps nt can be modified

GeoParams.time_τII_0D!Method
time_τII_0D!(τ_vec::Vector{NTuple{N,T}}, τII_vec::Vector{T}, v::Union{CompositeRheology,Tuple, Parallel}, ε_vec::Vector{NTuple{N,T}}, args, t_vec::AbstractVector{T}; verbose=false) where {N,T}

Computes stress-time evolution for a 0D (homogeneous) setup with given strainrate tensor (which can vary with time).

GeoParams.time_τII_0D!Method
time_τII_0D!(τ_vec::Vector{T}, v::CompositeRheology, εII_vec::Vector{T}, args, t_vec::AbstractVector{T}) where {T}

Computes stress-time evolution for a 0D (homogeneous) setup with given strainrate vector (which can vary with time).

GeoParams.time_τII_0DMethod
t_vec, τ_vec = time_τII_0D(v::CompositeRheology, εII::Number, args; t=(0.,100.), τ0=0., nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, strainrate εII and rheology arguments args. The initial stress τ0, the time range t and the number of timesteps nt can be modified

GeoParams.time_τII_0DMethod
t_vec, τ_vec = time_τII_0D(v::CompositeRheology, ε::NTuple{N,T}, args; t=(0.,100.), τ0=NTuple{N,T}, nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, strainrate tensor ε and rheology arguments args. The initial deviatoric stress tensor τ, the time range t and the number of timesteps nt can be modified

–>