Hello everyone,
I’m currently working on solving an ODE problem using Julia, but I’ve run into some issues. Before diving into possible modifications to my model, I’d like to first rule out any potential problems arising from the code settings or the solver.
Below is a snippet of the key parts of my code. If anyone needs more information to pinpoint the problem, please let me know, and I’ll gladly provide it.
# assemble the problem
uinit = ArrayPartition(vinit, θinit, δinit,𝓅init) # in this order
prob = assemble(gf, pf, dila, uinit, (0.0, 200.0 * 365 * 86400)) # `ODEProblem`
# setup saving scheme
handler(u::ArrayPartition, t, integrator) = (u.x[1], u.x[2], integrator(integrator.t, Val{1}).x[4], u.x[3], u.x[4],integrator(integrator.t, Val{1}).x[1],integrator(integrator.t, Val{1}).x[2]) # velocity, state variable, slip, dilatancy respectively
output = joinpath(@__DIR__, "output-dila.h5")
FFTW.set_num_threads(Threads.nthreads()) # full threads for FFT
@time sol = wsolve(
prob,
VCABM5(), # ODE solver
output, # path to the solution file
100, # the buffer size of the solution to be written
handler, # the function to retrieve the desired solution components from the integrator
["v", "θ", "d𝓅", "δ", "𝓅","dv","dθ"], # name corresponding to each components in `handler`
"t"; # name for the time step
reltol=1e-6, abstol=1e-6, dtmax=0.2*365*86400, dt=1e-2, maxiters=1e9, # ODE solver options
stride=100, # downsampling rate in saving the output
force=true, # whether overwrite existing solution file
)
The model stops at 16 year. The issue I have is as below.
Warning: dt(9.5367431640625e-7) <= dtmin(9.5367431640625e-7) at t=4.921133612180672e8. Aborting. There is either an error in your model specification or the true solution is unstable.
Then I lower the dtmin to see what will happen.
@time sol = wsolve(
prob,
VCABM5(), # ODE solver
output, # path to the solution file
100, # the buffer size of the solution to be written
handler, # the function to retrieve the desired solution components from the integrator
["v", "θ", "d𝓅", "δ", "𝓅","dv","dθ"], # name corresponding to each components in `handler`
"t"; # name for the time step
reltol=1e-5, abstol=1e-6, dtmax=0.2*365*86400, dtmin=1e-8, dt=1e-2, maxiters=1e9, # ODE solver options
stride=100, # downsampling rate in saving the output
force=true, # whether overwrite existing solution file
)
I got result like below. It does extend my computation to 30 years but it stoped at the end and stuck without any information showing.
I also changed my solver to Vern() and the result is even worse. And I got error message that shows instability is detected. The result is below. Time step seems normal.
Any insights or suggestions are greatly appreciated!
The whole code is below:
using Oetqf, JLD2
# setup the fault mesh
# Length: 1000 km; Depth: 20 km; grid space in both direction: 500 m; fault angle: 90 degree
mf = gen_mesh(Val(:RectOkada), 80e3, 10e3, 80e3/512, 10e3/64, 90.0)
save(joinpath(@__DIR__, "mesh.jld2"), Dict("fault" => mf))
# setup the Green's function parameters
λ = μ = 3e10 # Pa
# compute the Green's function
@time gf = stress_greens_function(mf, λ, μ; buffer_ratio=1)
# save the Green's function kernel
gffile = joinpath(@__DIR__, "gf-1km.h5")
isfile(gffile) && rm(gffile)
h5write(gffile, "gf", gf)
using Oetqf
using Printf, PyPlot, JLD2, ProgressMeter
# using MKL
# load previous fault mesh
mf = load(joinpath(@__DIR__, "mesh.jld2"), "fault")
# setup basic fault parameters
begin
a = 0.015
ab = 0.0035 # a - b
L = 4e-3 # critical distance (m)
σmax = 45e6 # confining pressure (Pa)
f₀ = 0.6 # reference friction coefficient
V₀ = 1e-6 # reference velocity (m/s)
μ = 3e10 # shear modulus (Pa)
Cs = 3044.14 # shear wave velocity (m/s)
Vpl = 140e-3 / 365 / 86400 # plate rate (m/s)
η = μ / 2Cs # radiation damping (Pa s/m)
hˢ = 2μ * (a + ab) * L / π / ab ^ 2 / σmax # nucleation size
@printf("Nucleation size: %.2f m\n", hˢ)
Λ₀ = μ * L / (a + ab) / σmax #cohesive zone
@info "cohesive zone size: $(Λ₀), resolved by $(Λ₀ / mf.Δx) cells."
φ = 0.05 #default is 0.05 fault porsity
tₚ = 1e9 * φ # s 1e9 * φ = diffustion time scale
ϵ = 2e-4 #backgroud dilatancy efficient
ϵd = 1e-6#rupture patch dilatancy efficient
β = φ * (5e-4 + 1e-2) * 1e-6#φ * (βf + βφ) = gouge bulk compressibility /Pa
𝓅init = 2e5
end
# resolution assertion
# @assert hˢ ≥ 8mf.Δx && hˢ ≥ 8mf.Δξ
# setup the concrete fault parameters
begin
a′ = ones(mf.nx, mf.nξ) .* a
b′ = ones(mf.nx, mf.nξ) .* (a + ab)
L′ = ones(mf.nx, mf.nξ) .* L
#σ = [min(σmax, 1.5e6 + 18.0e3 * z) for z in -mf.z]
#σ′ = Matrix(repeat(σ, 1, mf.nx)')
σ′ = ones(mf.nx, mf.nξ) .* σmax
#dilatancy parameters
ϵ′ = ones(mf.nx, mf.nξ) .* ϵ
tₚ′ = ones(mf.nx, mf.nξ) .* tₚ
β′ = ones(mf.nx, mf.nξ) .* β
p₀′ = zeros(mf.nx, mf.nξ)
# rupture patch
x = 10e3
dila_patch = @. -x ≤ mf.x ≤ x
vertical_patch = @. (-1e3 - 5e3) ≤ mf.z ≤ -1e3
ϵ′[xor.(dila_patch), vertical_patch] .= ϵd
pf = RateStateQuasiDynamicProperty(a′, b′, L′, σ′, η, Vpl, f₀, V₀)
save_property(joinpath(@__DIR__, "step-2-fault-parameters.bson"), pf)
dila = DilatancyProperty(tₚ′, ϵ′, β′, p₀′)
save_property(joinpath(@__DIR__,"step-2-fault-dila-parameters.bson"),dila)
begin # visually check the parameter setting
fig, ax = plt.subplots(1, 1, figsize=(10, 4))
img = ax.imshow((pf.a - pf.b)', extent=[-500, 500, -20, 0], cmap=plt.get_cmap("bwr_r"), vmin=-ab*1.2, vmax=ab*1.2)
ax.set_aspect(10)
ax.set_xlabel("Along-strike (km)")
ax.set_ylabel("Depth (km)")
cax = fig.add_axes([0.4, 0.8, 0.3, 0.01])
cbar = fig.colorbar(img, cax=cax, orientation="horizontal")
cbar.ax.set_title(L"a-b")
fig.savefig(joinpath(@__DIR__, "step-2-para-ab.pdf"), bbox_inches="tight")
end
end
# setup initial condition
vinit = Vpl .* ones(mf.nx, mf.nξ)
θinit = pf.L ./ vinit
# offset two halves
θinit[1: size(θinit, 1) >> 1, :] ./= 1.1
θinit[size(θinit, 1) >> 1 + 1: end, :] ./= 1.2
δinit = zeros(mf.nx, mf.nξ)
𝓅init = p₀′
#𝓅init = ones(mf.nx, mf.nξ) .* 𝓅init
# load the pre-computed Green's function
gf = h5read(joinpath(@__DIR__, "gf-1km.h5"), "gf")
# assemble the problem
uinit = ArrayPartition(vinit, θinit, δinit,𝓅init) # in this order
prob = assemble(gf, pf, dila, uinit, (0.0, 200.0 * 365 * 86400)) # `ODEProblem`
# setup saving scheme
handler(u::ArrayPartition, t, integrator) = (u.x[1], u.x[2], integrator(integrator.t, Val{1}).x[4], u.x[3], u.x[4],integrator(integrator.t, Val{1}).x[1],integrator(integrator.t, Val{1}).x[2]) # velocity, state variable, slip, dilatancy respectively
output = joinpath(@__DIR__, "output-dila.h5")
FFTW.set_num_threads(Threads.nthreads()) # full threads for FFT
@time sol = wsolve(
prob,
VCABM5(), # ODE solver
output, # path to the solution file
100, # the buffer size of the solution to be written
handler, # the function to retrieve the desired solution components from the integrator
["v", "θ", "d𝓅", "δ", "𝓅","dv","dθ"], # name corresponding to each components in `handler`
"t"; # name for the time step
reltol=1e-5, abstol=1e-6, dtmax=0.2*365*86400, dtmin=1e-8, dt=1e-2, maxiters=1e9, # ODE solver options
stride=100, # downsampling rate in saving the output
force=true, # whether overwrite existing solution file
)
The code package can be found here: GitHub - RiveHe/Oetqf.jl at upstream
