using Oetqf, SpecialFunctions, Optim

Generate the mesh for the transform fault, which is suited for using Okaka, (1992) equation:

mf = gen_mesh(Val(:RectOkada), 80e3, 8e3, 10e3, 2e3, 90.0);

Use Gmsh to generate the mantle mesh, which is suited for using Barbot et al., (2017) equation, with no refinement in x or y direction while cell sizes are 1.5 times progressively larger along z axes:

gen_gmsh_mesh(Val(:BEMHex8Mesh), -40e3, -2.5e3, -8e3, 80e3, 5e3, -22e3, 4, 3, 3;
    output = joinpath(@__DIR__, "mantle.vtk"),
    rfzh = cumprod(ones(3) * 1.5), rfy = 1.0, rfyType = "Bump"
)
ma = gen_mesh(Val(:BEMHex8Mesh), joinpath(@__DIR__, "mantle.vtk"));

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 20%] Meshing curve 3 (Line)
Info    : [ 30%] Meshing curve 4 (Line)
Info    : [ 40%] Meshing curve 6 (Extruded)
Info    : [ 50%] Meshing curve 7 (Extruded)
Info    : [ 50%] Meshing curve 8 (Extruded)
Info    : [ 60%] Meshing curve 9 (Extruded)
Info    : [ 70%] Meshing curve 11 (Extruded)
Info    : [ 80%] Meshing curve 12 (Extruded)
Info    : [ 90%] Meshing curve 16 (Extruded)
Info    : [100%] Meshing curve 20 (Extruded)
Info    : Done meshing 1D (Wall 0.000227002s, CPU 0.000223s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Transfinite)
Info    : [ 20%] Meshing surface 13 (Extruded)
Info    : [ 40%] Meshing surface 17 (Extruded)
Info    : [ 50%] Meshing surface 21 (Extruded)
Info    : [ 70%] Meshing surface 25 (Extruded)
Info    : [ 90%] Meshing surface 26 (Extruded)
Info    : Done meshing 2D (Wall 0.000319003s, CPU 0.00031s)
Info    : Meshing 3D...
Info    : Meshing volume 1 (Extruded)
Info    : Done meshing 3D (Wall 0.000125501s, CPU 0.000122s)
Info    : Optimizing mesh...
Info    : Done optimizing mesh (Wall 3.9e-06s, CPU 4e-06s)
Info    : 80 nodes 150 elements
Info    : Writing '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'...
Info    : Done writing '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'
Info    : Reading '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'...
Info    : Reading 80 points
Info    : Reading 150 cells
Info    : Done reading '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'

Compute the stress Green's function between the two meshes:

λ = μ = 3e10
gffile = joinpath(@__DIR__, "gf.h5")
isfile(gffile) && rm(gffile)
@time gf₁₁ = stress_greens_function(mf, λ, μ; buffer_ratio = 1)
h5write(gffile, "gf₁₁", gf₁₁) # fault -> fault
@time gf₁₂ = stress_greens_function(mf, ma, λ, μ; buffer_ratio = 1, qtype = "Gauss1")
h5write(gffile, "gf₁₂", gf₁₂) # fault -> mantle
@time gf₂₁ = stress_greens_function(ma, mf, λ, μ)
h5write(gffile, "gf₂₁", gf₂₁) # mantle -> fault
@time gf₂₂ = stress_greens_function(ma, λ, μ; qtype = "Gauss1")
h5write(gffile, "gf₂₂", gf₂₂) # mantle -> mantle

  2.148996 seconds (3.05 M allocations: 155.368 MiB, 3.69% gc time, 99.91% compilation time)
  0.490442 seconds (507.99 k allocations: 26.954 MiB, 81.86% compilation time)
  6.090463 seconds (9.78 M allocations: 426.988 MiB, 3.70% gc time, 71.95% compilation time)
Maximum real part of eigval is: 1111120100.9503
 12.751689 seconds (1.33 M allocations: 71.680 MiB, 0.31% gc time, 8.29% compilation time)

Set up the rate-and-state friction parameters in the fault:

cs = 3044.14 # m/s
vpl = 140e-3 / 365 / 86400 # 140 mm/yr
v0 = 1e-6
f0 = 0.6
μ = 3e10
η = μ / 2cs # radiation damping
ν = λ / 2(λ + μ)
avw = 0.015
abvw = 0.0047
Dc = 8e-3
σmax = 5e7
a = ones(mf.nx, mf.nξ) .* avw
b = ones(mf.nx, mf.nξ) .* (avw - abvw)
L = ones(mf.nx, mf.nξ) .* Dc
σ = [min(σmax, 1.5e6 + 18.0e3 * z) for z in -mf.z] # Pa
σ = repeat(σ, 1, mf.nx)' |> Matrix # Pa
left_patch = @. -25.e3 ≤ mf.x ≤ -5.e3
right_patch = @. 5.e3 ≤ mf.x ≤ 25.e3
vert_patch = @. -6.e3 ≤ mf.z ≤ -1e3
b[xor.(left_patch, right_patch), vert_patch] .= avw + abvw # assign velocity weakening
pf = RateStateQuasiDynamicProperty(a, b, L, σ, η, vpl, f0, v0)
save_property(joinpath(@__DIR__, "para-fault.bson"), pf);

Set up rheology parameters in the mantle assuming power-law viscosity with lab-derived results:

A_wet_dis = 3e1
Q_wet_dis = 480e3
V_wet_dis = 11e-6
m_wet_dis = 0
r_wet_dis = 1.2
n_wet_dis = 3.5
grain_size = 10000.0 # μm
COH = 1000 # ppm / HSi
𝙍 = 8.314 # gas contant
crust_depth = 7e3
κ = 8e-7
𝚃(z) = 1673 * erf(z / sqrt(4κ * 1e6 * 365 * 86400)) # 1 Myr OTF
𝙿(z) = 2800 * 9.8 * crust_depth + 3300 * 9.8 * (z - crust_depth)
prefactor_dis(z) = A_wet_dis / (1e6)^n_wet_dis * COH^r_wet_dis * grain_size^m_wet_dis * exp(-(Q_wet_dis + 𝙿(z) * V_wet_dis) / 𝙍 / 𝚃(z))
rel_dϵ = [0.0, -1e-12, 0.0, 0.0, 0.0, 0.0]
amplifier = 1e0
γ_dis = prefactor_dis.(-ma.cz) .* amplifier
pa = PowerLawViscosityProperty(γ_dis, ones(length(ma.cz)) * (n_wet_dis - 1), rel_dϵ) # notice to save `n-1` instead of `n` where `n` refers the stress power
save_property(joinpath(@__DIR__, "para-mantle" * ".bson"), pa);

Set up initial conditions on the fault with an offset between left and right half fault:

vinit = pf.vpl .* ones(size(pf.a))
θinit = pf.L ./ vinit
θinit[1: size(θinit, 1) >> 1, :] ./= 1.1
θinit[size(θinit, 1) >> 1 + 1: end, :] ./= 2.5
δinit = zeros(size(pf.a));

Set up initial conditions in the mantle

ϵinit = zeros(length(pa.γ), 6)
P = map(z -> 2800 * 9.8 * crust_depth + 3300 * 9.8 * (z - crust_depth), -ma.cz) # change the depth of crust
σinit = repeat(P, 1, 6)
σinit[:,3] .= 0.0 # xz
σinit[:,5] .= 0.0 # yz
target(i) = x -> (pa.γ[i] * (sqrt(2) * x) ^ (pa.n[i]) * x - abs(pa.dϵ₀[2])) ^ 2
σxyinit = -map(i -> Optim.minimizer(optimize(target(i), 1e1, 1e14)), 1: length(pa.γ))
reldϵ = map(i -> pa.γ[i] * (sqrt(2) * abs(σxyinit[i])) ^ (pa.n[i]) * σxyinit[i], 1: length(pa.γ))
@assert all(isapprox.(reldϵ, pa.dϵ₀[2]; rtol=1e-3))
σinit[:,2] .= σxyinit;

Assemble the problem:

uinit = ArrayPartition(vinit, θinit, ϵinit, σinit, δinit)
prob = assemble(gf₁₁, gf₁₂, gf₂₁, gf₂₂, pf, pa, uinit, (0.0, 0.1 * 365 * 86400));

Set up the saving scheme and solve the equation:

handler(u::ArrayPartition, t, integrator) = (u.x[1], u.x[2], integrator(integrator.t, Val{1}).x[3], u.x[3], u.x[4], u.x[5])
output = joinpath(@__DIR__, "output.h5")
@time sol = wsolve(prob, VCABM5(), output, 100, handler, ["v", "θ", "dϵ", "ϵ", "σ", "δ"], "t";
    reltol=1e-6, abstol=1e-8, dtmax=0.2*365*86400, dt=1e-8, maxiters=1e9, stride=100, force=true
)

retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 3.1536e6
u: 2-element Vector{ArrayPartition{Float64, NTuple{5, Matrix{Float64}}}}:
 ([4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; … ; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9], [1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6; 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6; … ; 720822.857142857 720822.857142857 720822.857142857 720822.857142857; 720822.857142857 720822.857142857 720822.857142857 720822.857142857], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [2.993126315789474e8 -1.159793132258934e6 … 0.0 2.993126315789474e8; 4.865442105263158e8 -279723.0841659712 … 0.0 4.865442105263158e8; … ; 4.865442105263158e8 -279723.0841659712 … 0.0 4.865442105263158e8; 7.673915789473685e8 -235923.5595516861 … 0.0 7.673915789473685e8], [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; … ; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0])
 ([4.0767159594235405e-9 4.0726161646218115e-9 4.075442712505631e-9 4.088964986854444e-9; 3.4430994936084665e-9 3.417383507316748e-9 3.4394142405451146e-9 4.062644180138336e-9; … ; 7.586803545971388e-10 7.350000917671253e-10 7.495634426596897e-10 2.050210814376089e-9; 2.0738117751487143e-9 2.038088908694873e-9 2.04703978043994e-9 2.0960945203688157e-9], [1.859041830005468e6 1.8598256458685251e6 1.8593496130393536e6 1.857018956937463e6; 2.0017790267709633e6 2.0064139281263128e6 2.002447943695447e6 1.861281287169844e6; … ; 2.87470435919983e6 2.885326685849754e6 2.878660033158451e6 2.27726833669272e6; 2.270762721125289e6 2.279965835489089e6 2.2775965956855654e6 2.264780686989732e6], [-1.9490353620031998e-10 -3.1502254797285903e-6 … 1.1170093924188181e-10 -5.760585512089933e-12; -3.143400768082191e-10 -3.147648023537708e-6 … 4.658202394508259e-10 -4.579661829808136e-11; … ; -5.615414365772426e-10 -3.130506294386605e-6 … -9.986930109540351e-10 -6.37906389157408e-10; -1.2895019450552975e-10 -3.144725879035746e-6 … -6.507466234411018e-10 -2.713827658666108e-10], [2.993123106670303e8 -1.1588113187960377e6 … 109.26822528313411 2.993125101945185e8; 4.865440819229595e8 -279332.7091137653 … 107.15415362778771 4.865441499927153e8; … ; 4.865439287083998e8 -278353.0326551578 … -211.77387653432862 4.865439457601052e8; 7.673915013824514e8 -235463.2839044049 … -117.59525367783492 7.673914887318368e8], [0.0133529772462193 0.013347861796139792 0.01335086397318031 0.01336568039662874; 0.012402031284312586 0.01237569089985771 0.01239823817095998 0.013339120424513825; … ; 0.005481590408943952 0.005435821491326209 0.005464692151618529 0.008821136598718868; 0.008856835911829546 0.008807659334967277 0.008820414519594021 0.0088895007217099])

This page was generated using Literate.jl.