Example 2D
Note
This example is adapted from Wei, 2016 AGU
Tip
It will automatically use parallel scheme if nprocs() ≂̸ 1 when building stiffness tensor. To do so:
1 2 3 | using Distributed addprocs(4); # add # of cores you desire using JuEQ |
Define Parameters¶
First, we load the package and define some basic parameters:
1 2 3 4 5 6 7 8 9 | using JuEQ using Plots ms2mmyr = 365 * 86400 * 1e3 ρ = 2670.0 # kg/m³ cs = 3044.0 # m/s vpl = 100.0 # mm/yr v0 = 3.2e4 # mm/yr f0 = 0.6; |
1 | 0.6 |
Then we come to parameters implicit by above:
1 2 3 4 | μ = 0.3 # Bar·km/mm λ = μ # poisson material α = (λ + μ) / (λ + 2μ) η = μ / 2(cs * 1e-3 * 365 * 86400); # Bar·yr/mm |
1 | 1.562571878306402e-9 |
Create a fault:
1 | fa = fault(StrikeSlipFault, (80., 10.)); |
1 | PlaneFaultDomain{StrikeSlipFault,2,Float64}(90.0, (80.0, 10.0))
|
Generate grids:
1 | gd = discretize(fa; nx=160, nξ=20, buffer_ratio=1); |
1 | JuEQ.BEMGrid_2D{Array{Float64,1},Array{Array{Float64,1},1},Float64,Int64}([-39.75, -39.25, -38.75, -38.25, -37.75, -37.25, -36.75, -36.25, -35.75, -35.25 … 35.25, 35.75, 36.25, 36.75, 37.25, 37.75, 38.25, 38.75, 39.25, 39.75], [-0.25, -0.75, -1.25, -1.75, -2.25, -2.75, -3.25, -3.75, -4.25, -4.75, -5.25, -5.75, -6.25, -6.75, -7.25, -7.75, -8.25, -8.75, -9.25, -9.75], 0.5, 0.5, 160, 20, Array{Float64,1}[[-40.0, -39.5], [-39.5, -39.0], [-39.0, -38.5], [-38.5, -38.0], [-38.0, -37.5], [-37.5, -37.0], [-37.0, -36.5], [-36.5, -36.0], [-36.0, -35.5], [-35.5, -35.0] … [35.0, 35.5], [35.5, 36.0], [36.0, 36.5], [36.5, 37.0], [37.0, 37.5], [37.5, 38.0], [38.0, 38.5], [38.5, 39.0], [39.0, 39.5], [39.5, 40.0]], Array{Float64,1}[[-0.5, 0.0], [-1.0, -0.5], [-1.5, -1.0], [-2.0, -1.5], [-2.5, -2.0], [-3.0, -2.5], [-3.5, -3.0], [-4.0, -3.5], [-4.5, -4.0], [-5.0, -4.5], [-5.5, -5.0], [-6.0, -5.5], [-6.5, -6.0], [-7.0, -6.5], [-7.5, -7.0], [-8.0, -7.5], [-8.5, -8.0], [-9.0, -8.5], [-9.5, -9.0], [-10.0, -9.5]], [-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0], [-0.25, -0.75, -1.25, -1.75, -2.25, -2.75, -3.25, -3.75, -4.25, -4.75, -5.25, -5.75, -6.25, -6.75, -7.25, -7.75, -8.25, -8.75, -9.25, -9.75], 1)
|
Tip
It is recommended (from Yajing Liu's personal communication) to add buffer zones adjacent the horizontal edges to immitate zero dislocation at the ridge region. Basically, it affects how the stiffness tensor are periodically summed. To what extent it alters the results remains further testing.
Under the hood, it shall impose buffer areas on both sides of along-strike, each of which has a length of bufferratio/2*fa[:x]. Thus, the stiffness contributions falling into those buffer zone shall be neglected, which is equivalent to impose zero-slip correspondingly.
Time for us to establish frictional parameters profile:
1 2 3 4 5 6 7 8 9 10 11 | a = 0.015 .* ones(gd.nx, gd.nξ) b = 0.0115 .* ones(gd.nx, gd.nξ) left_patch = @. -25. ≤ gd.x ≤ -5. right_patch = @. 5. ≤ gd.x ≤ 25. vert_patch = @. -6. ≤ gd.z ≤ -1. b[xor.(left_patch, right_patch), vert_patch] .= 0.0185 amb = a - b σmax = 500. σ = [min(σmax, 15. + 180. * z) for z in -gd.z] σ = Matrix(repeat(σ, 1, gd.nx)') L = 12.; |
1 | 12.0 |
Check our profile:
1 2 3 4 5 6 7 8 9 10 11 12 13 | p1 = heatmap(amb', xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2), yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]), yflip=true, color=:isolum, aspect_ratio=2, title="a-b" ); p2 = heatmap(σ', xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2), yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]), yflip=true, color=:isolum, aspect_ratio=2, title="\\sigma" ); plot(p1, p2, layout=(2, 1)) |
Documenter.Documents.RawHTML("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n\n")
Construct Model¶
Construct our material property profile:
1 | mp = properties(fa, gd, [:a=>a, :b=>b, :L=>L, :σ=>σ, :η=>η, :k=>[:λ=>λ, :μ=>μ], :vpl=>vpl, :f0=>f0, :v0=>v0]); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | [ Info: Calculating stiffness tensor ...
[ Info: Fault material properties establised.
JuEQ.PlaneMaterialProperties{2,Float64,Array{Float64,2},Array{Complex{Float64},3}}
dims: Tuple{Int64,Int64}
a: Array{Float64}((160, 20)) [0.015 0.015 … 0.015 0.015; 0.015 0.015 … 0.015 0.015; … ; 0.015 0.015 … 0.015 0.015; 0.015 0.015 … 0.015 0.015]
b: Array{Float64}((160, 20)) [0.0115 0.0115 … 0.0115 0.0115; 0.0115 0.0115 … 0.0115 0.0115; … ; 0.0115 0.0115 … 0.0115 0.0115; 0.0115 0.0115 … 0.0115 0.0115]
L: Array{Float64}((160, 20)) [12.0 12.0 … 12.0 12.0; 12.0 12.0 … 12.0 12.0; … ; 12.0 12.0 … 12.0 12.0; 12.0 12.0 … 12.0 12.0]
k: Array{Complex{Float64}}((160, 20, 20)) Complex{Float64}[0.254648+0.0im -0.152789+0.0im … -0.00055943+0.0im -0.000503351+0.0im; 0.254915-9.92045e-18im -0.152567-7.18284e-19im … -0.000455488-2.85874e-21im -0.00040226-9.21148e-21im; … ; 0.78401-3.89143e-18im -0.0149563-6.43374e-18im … 6.78809e-13-1.71483e-21im 2.84531e-14+1.35689e-23im; 0.784076+2.99306e-19im -0.0149479-1.53641e-18im … 4.33558e-13-1.04805e-22im -1.86728e-13+5.91773e-22im]
Complex{Float64}[-0.152789+0.0im 0.371059+0.0im … -0.000569356+0.0im -0.00051138+0.0im; -0.15257-2.79182e-18im 0.37131-9.40545e-18im … -0.000469-6.77626e-21im -0.000413734-5.39984e-21im; … ; -0.0173791+2.19565e-18im 0.826915+8.36764e-19im … -4.6572e-13-1.34349e-21im 6.73039e-13-7.19082e-22im; -0.0173699-1.04872e-18im 0.826977+2.85134e-18im … -4.86674e-13-4.91837e-23im 4.39738e-13-6.26032e-22im]
Complex{Float64}[-0.0363781+0.0im -0.133387+0.0im … -0.000589806+0.0im -0.000527873+0.0im; -0.0361831-2.1684e-19im -0.133181-7.99599e-19im … -0.000492091-2.75286e-21im -0.00043281-3.59989e-21im; … ; 9.42205e-5-4.75247e-19im -0.0138966-1.67687e-18im … 4.78959e-13+1.37684e-21im -4.54393e-13+4.76101e-22im; 9.44e-5-6.84032e-19im -0.0138888-2.03806e-18im … 8.00646e-14-1.39889e-22im -4.8123e-13-6.10271e-23im]
...
Complex{Float64}[-0.000625436+0.0im -0.000637855+0.0im … -0.127397+0.0im -0.0255343+0.0im; -0.000518554-1.27055e-20im -0.000534698-1.27055e-21im … -0.127238-1.28749e-19im -0.0253991-3.45589e-19im; … ; 9.26806e-14-1.29914e-21im 2.0563e-13+1.61014e-21im … -0.0138824+1.01786e-18im 0.000228985-5.75375e-19im; 4.21404e-14-7.87087e-23im -2.05234e-13+5.16603e-22im … -0.0138746-2.00471e-18im 0.000228943-1.22109e-20im]
Complex{Float64}[-0.00055943+0.0im -0.000569356+0.0im … 0.381902+0.0im -0.12739+0.0im; -0.00045553-4.5528e-21im -0.000469034-6.24687e-21im … 0.382104+1.51788e-18im -0.127231-1.81604e-18im; … ; 1.53612e-13-4.10456e-22im 8.4465e-14+1.09244e-21im … 0.827193+2.16326e-18im -0.0138824-2.40971e-18im; -9.174e-14+3.52978e-23im 6.35933e-14-1.94583e-22im … 0.827255-3.09464e-18im -0.0138746-6.87684e-19im]
Complex{Float64}[-0.000503351+0.0im -0.00051138+0.0im … -0.12739+0.0im 0.381909+0.0im; -0.000402298-1.27055e-21im -0.000413764-7.41154e-21im … -0.127231+7.72494e-19im 0.38211-1.14248e-17im; … ; 1.45108e-13-2.9691e-22im 1.478e-13-5.9501e-22im … -0.0138824+4.56513e-18im 0.827193+8.6732e-18im; -7.00794e-14+6.1893e-22im -8.55222e-14-9.09477e-24im … -0.0138746-4.34537e-19im 0.827255-2.18242e-18im]
σ: Array{Float64}((160, 20)) [60.0 150.0 … 500.0 500.0; 60.0 150.0 … 500.0 500.0; … ; 60.0 150.0 … 500.0 500.0; 60.0 150.0 … 500.0 500.0]
η: Array{Float64}((160, 20)) [1.56257e-9 1.56257e-9 … 1.56257e-9 1.56257e-9; 1.56257e-9 1.56257e-9 … 1.56257e-9 1.56257e-9; … ; 1.56257e-9 1.56257e-9 … 1.56257e-9 1.56257e-9; 1.56257e-9 1.56257e-9 … 1.56257e-9 1.56257e-9]
vpl: Float64 100.0
f0: Float64 0.6
v0: Float64 32000.0
|
Provide the initial condition:
1 2 3 | vinit = vpl .* ones(gd.nx, gd.nξ) θ0 = L ./ vinit ./ 1.1 u0 = cat(vinit, θ0, dims=3); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | 160×20×2 Array{Float64,3}:
[:, :, 1] =
100.0 100.0 100.0 100.0 100.0 … 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 … 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
⋮ ⋱ ⋮
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 … 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
[:, :, 2] =
0.109091 0.109091 0.109091 0.109091 … 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 … 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
⋮ ⋱
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 … 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
0.109091 0.109091 0.109091 0.109091 0.109091 0.109091 0.109091
|
Get our ODEs problem:
1 | prob = EarthquakeCycleProblem(gd, mp, u0, (0., 18.); se=DieterichStateLaw(), fform=CForm()); |
1 2 3 4 5 | ODEProblem with uType Array{Float64,3} and tType Float64. In-place: true
timespan: (0.0, 18.0)
u0: [100.0 100.0 … 100.0 100.0; 100.0 100.0 … 100.0 100.0; … ; 100.0 100.0 … 100.0 100.0; 100.0 100.0 … 100.0 100.0]
[0.109091 0.109091 … 0.109091 0.109091; 0.109091 0.109091 … 0.109091 0.109091; … ; 0.109091 0.109091 … 0.109091 0.109091; 0.109091 0.109091 … 0.109091 0.109091]
|
Solve Model¶
Solve the model:
1 | sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3009-element Array{Float64,1}:
0.0
0.023024248421406265
0.0351255887693807
0.05876569503928833
0.07976766186029392
0.11039745942303938
0.1434255742236179
0.18505221605685096
0.23039961851701168
0.28149227315967174
⋮
17.620910529073324
17.66877212637217
17.71657790765609
17.764324436550062
17.81200818883823
17.859625294155794
17.90717089472236
17.954637490041815
18.0
u: 3009-element Array{Array{Float64,3},1}:
[100.0 100.0 … 100.0 100.0; 100.0 100.0 … 100.0 100.0; … ; 100.0 100.0 … 100.0 100.0; 100.0 100.0 … 100.0 100.0]
[0.109091 0.109091 … 0.109091 0.109091; 0.109091 0.109091 … 0.109091 0.109091; … ; 0.109091 0.109091 … 0.109091 0.109091; 0.109091 0.109091 … 0.109091 0.109091]
[98.9146 98.7561 … 98.6482 98.6672; 98.7364 98.6541 … 98.614 98.6385; … ; 98.7364 98.6541 … 98.614 98.6385; 98.9146 98.7561 … 98.6482 98.6672]
[0.111114 0.111126 … 0.111134 0.111132; 0.111127 0.111133 … 0.111136 0.111134; … ; 0.111127 0.111133 … 0.111136 0.111134; 0.111114 0.111126 … 0.111134 0.111132]
[98.5394 98.2571 … 98.0281 98.0706; 98.2146 98.0464 … 97.9516 98.0067; … ; 98.2146 98.0464 … 97.9516 98.0067; 98.5394 98.2571 … 98.0281 98.0706]
[0.112104 0.112139 … 0.112163 0.112159; 0.112143 0.112162 … 0.112171 0.112166; … ; 0.112143 0.112162 … 0.112171 0.112166; 0.112104 0.112139 … 0.112163 0.112159]
[98.0049 97.5137 … 96.9797 97.0899; 97.3976 97.0433 … 96.7821 96.9264; … ; 97.3976 97.0433 … 96.7821 96.9264; 98.0049 97.5137 … 96.9797 97.0899]
[0.11387 0.113979 … 0.114077 0.114058; 0.113999 0.114068 … 0.114111 0.114086; … ; 0.113999 0.114068 … 0.114111 0.114086; 0.11387 0.113979 … 0.114077 0.114058]
[97.6637 97.0474 … 96.2126 96.4019; 96.8364 96.3324 … 95.874 96.1241; … ; 96.8364 96.3324 … 95.874 96.1241; 97.6637 97.0474 … 96.2126 96.4019]
[0.115256 0.115451 … 0.11566 0.115617; 0.115497 0.115636 … 0.115739 0.115681; … ; 0.115497 0.115636 … 0.115739 0.115681; 0.115256 0.115451 … 0.11566 0.115617]
[97.3008 96.5932 … 95.3364 95.6632; 96.2099 95.5445 … 94.7547 95.1918; … ; 96.2099 95.5445 … 94.7547 95.1918; 97.3008 96.5932 … 95.3364 95.6632]
[0.116988 0.117317 … 0.117761 0.117658; 0.117432 0.117697 … 0.117945 0.117808; … ; 0.117432 0.117697 … 0.117945 0.117808; 0.116988 0.117317 … 0.117761 0.117658]
[97.0249 96.3026 … 94.6681 95.1597; 95.7147 94.9515 … 93.7973 94.4637; … ; 95.7147 94.9515 … 93.7973 94.4637; 97.0249 96.3026 … 94.6681 95.1597]
[0.118521 0.118979 … 0.119739 0.119542; 0.119207 0.119618 … 0.120092 0.119827; … ; 0.119207 0.119618 … 0.120092 0.119827; 0.118521 0.118979 … 0.119739 0.119542]
[96.7876 96.1141 … 94.1599 94.865; 95.275 94.4751 … 92.9191 93.8905; … ; 95.275 94.4751 … 92.9191 93.8905; 96.7876 96.1141 … 94.1599 94.865]
[0.120043 0.120617 … 0.12181 0.121452; 0.121037 0.121613 … 0.122447 0.121959; … ; 0.121037 0.121613 … 0.122447 0.121959; 0.120043 0.120617 … 0.12181 0.121452]
[96.6216 96.0368 … 93.9387 94.8624; 94.958 94.1876 … 92.3239 93.6197; … ; 94.958 94.1876 … 92.3239 93.6197; 96.6216 96.0368 … 93.9387 94.8624]
[0.121286 0.121926 … 0.123552 0.122981; 0.122592 0.123301 … 0.124566 0.123773; … ; 0.122592 0.123301 … 0.124566 0.123773; 0.121286 0.121926 … 0.123552 0.122981]
[96.5105 96.0338 … 93.9958 95.1269; 94.7385 94.0453 … 92.0305 93.6516; … ; 94.7385 94.0453 … 92.0305 93.6516; 96.5105 96.0338 … 93.9958 95.1269]
[0.122295 0.122944 … 0.12493 0.124097; 0.123904 0.124696 … 0.126408 0.125225; … ; 0.123904 0.124696 … 0.126408 0.125225; 0.122295 0.122944 … 0.12493 0.124097]
⋮
[96.8206 96.8583 … 98.6429 99.0521; 95.1492 95.1577 … 97.9583 98.6168; … ; 95.1492 95.1577 … 97.9583 98.6168; 96.8206 96.8583 … 98.6429 99.0521]
[0.124076 0.124025 … 0.121731 0.121202; 0.126331 0.126321 … 0.122628 0.121765; … ; 0.126331 0.126321 … 0.122628 0.121765; 0.124076 0.124025 … 0.121731 0.121202]
[96.8731 96.9101 … 98.6743 99.0733; 95.2289 95.2372 … 98.007 98.6486; … ; 95.2289 95.2372 … 98.007 98.6486; 96.8731 96.9101 … 98.6743 99.0733]
[0.124021 0.123972 … 0.121698 0.12118; 0.126246 0.126236 … 0.122576 0.121731; … ; 0.126246 0.126236 … 0.122576 0.121731; 0.124021 0.123972 … 0.121698 0.12118]
[96.9287 96.9648 … 98.7066 99.0951; 95.3133 95.3211 … 98.0566 98.6812; … ; 95.3133 95.3211 … 98.0566 98.6812; 96.9287 96.9648 … 98.7066 99.0951]
[0.123962 0.123913 … 0.121664 0.121157; 0.126153 0.126143 … 0.122522 0.121696; … ; 0.126153 0.126143 … 0.122522 0.121696; 0.123962 0.123913 … 0.121664 0.121157]
[96.9873 97.0224 … 98.7393 99.1173; 95.4023 95.4095 … 98.107 98.7143; … ; 95.4023 95.4095 … 98.107 98.7143; 96.9873 97.0224 … 98.7393 99.1173]
[0.123899 0.123851 … 0.121627 0.121132; 0.126053 0.126044 … 0.122465 0.121659; … ; 0.126053 0.126044 … 0.122465 0.121659; 0.123899 0.123851 … 0.121627 0.121132]
[97.0489 97.0828 … 98.7724 99.1398; 95.4958 95.5021 … 98.1578 98.7478; … ; 95.4958 95.5021 … 98.1578 98.7478; 97.0489 97.0828 … 98.7724 99.1398]
[0.123831 0.123784 … 0.121589 0.121107; 0.125947 0.125938 … 0.122406 0.121621; … ; 0.125947 0.125938 … 0.122406 0.121621; 0.123831 0.123784 … 0.121589 0.121107]
[97.1134 97.1459 … 98.8058 99.1626; 95.5937 95.5988 … 98.2088 98.7816; … ; 95.5937 95.5988 … 98.2088 98.7816; 97.1134 97.1459 … 98.8058 99.1626]
[0.123759 0.123713 … 0.12155 0.121081; 0.125835 0.125827 … 0.122345 0.121581; … ; 0.125835 0.125827 … 0.122345 0.121581; 0.123759 0.123713 … 0.12155 0.121081]
[97.1809 97.2118 … 98.8392 99.1855; 95.6961 95.6998 … 98.2598 98.8155; … ; 95.6961 95.6998 … 98.2598 98.8155; 97.1809 97.2118 … 98.8392 99.1855]
[0.123683 0.123639 … 0.121511 0.121054; 0.125716 0.125709 … 0.122284 0.121541; … ; 0.125716 0.125709 … 0.122284 0.121541; 0.123683 0.123639 … 0.121511 0.121054]
[97.2514 97.2805 … 98.8728 99.2087; 95.803 95.805 … 98.3106 98.8495; … ; 95.803 95.805 … 98.3106 98.8495; 97.2514 97.2805 … 98.8728 99.2087]
[0.123604 0.123561 … 0.121471 0.121027; 0.125592 0.125586 … 0.122222 0.121501; … ; 0.125592 0.125586 … 0.122222 0.121501; 0.123604 0.123561 … 0.121471 0.121027]
[97.3218 97.3491 … 98.9049 99.231; 95.91 95.9099 … 98.3592 98.8821; … ; 95.91 95.9099 … 98.3592 98.8821; 97.3218 97.3491 … 98.9049 99.231]
[0.123524 0.123484 … 0.121432 0.121001; 0.125468 0.125464 … 0.122162 0.121462; … ; 0.125468 0.125464 … 0.122162 0.121462; 0.123524 0.123484 … 0.121432 0.121001]
|
Sanity Check of Results¶
Take a look at the max velocity:
1 2 | maxv = max_velocity(sol) plot(sol.t, log10.(maxv / ms2mmyr), xlabel="Time (year)", ylabel="Max Velocity (log10 (m/s))", label="") |
Documenter.Documents.RawHTML("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n\n")
View some snapshots to see the rupture (quasi-dynamic) patterns:
1 2 3 4 5 6 7 8 9 | ind = argmax(maxv) myplot = (ind) -> heatmap(log10.(sol.u[ind][:,:,1]./ms2mmyr)', xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2), yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]), yflip=true, color=:isolum, aspect_ratio=2, title="t = $(sol.t[ind])") snaps = [myplot(i) for i in ind-700: 200: ind+500] plot(snaps..., layout=(length(snaps), 1), size=(600, 1800)) |
Documenter.Documents.RawHTML("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n\n")
This page was generated using Literate.jl.