Hello,

I am working on a project to model dehydration of some rocks. I am planning on adapting this model https://github.com/rmskarbek/PorosityWaves to different conditions. However, before moving on to that I wanted to see if I could replicate the results of the original MATLAB code in Julia using DiffEq.

The original MATLAB model uses **ode15s** as the solver and takes around 1-2 hours solving for a timespan of 120,000 years.

I rewrote the model in Julia to the best of my ability, but I am unable to solve the system in a reasonable amount of time for the original timespan of 0-120,000 years. I have been able to solve it over the timespan of 0-12,000 years fairly rapidly using **FBDF**. But the moment I increase the tspan to 120,000 years the solver progress slows to a crawl and I haven’t had the patience to wait for it to complete. I have even tried running it over night without it finishing, but ideally it should be faster than the MATLAB implementation.

Some reports using @btime to show how the evaluation time behaves with increasing tspan:

tspan = (0, 120000) → Never Finishes

tspan = (0, 12000) → 25.844 s (484132372 allocations: 23.47 GiB)

tspan = (0, 1200) → 9.151 s (108983605 allocations: 9.39 GiB)

tspan = (0, 120) → 6.795 s (59317085 allocations: 7.22 GiB)

tspan = (0, 12) → 5.266 s (38895461 allocations: 5.89 GiB)

Also, at first I tried using **QBDF** and **QNDF** as solvers, but they both suggest it will take ~2e12 hours to complete on the tspan of 0-12,000 years. I am not sure why this is the case.

Here is the code:

```
using DifferentialEquations
using ProgressLogging, SparseArrays, Statistics, Plots, Polynomials, IfElse, BenchmarkTools
# DeltG of Rxn.
function DeltaG(T, P_t)
## For Gibbs Calculation
a0 = 65945.7007021266
a1 = -415.847466841169
a2 = 0.654522687965514
dT = 58
a0 = a0 - a1*dT - 2*a2*dT + a2*dT^2
a1 = a1-2*a2*dT
g0 = -643.932375104237
return @. -g0*((1/(2*a2))*(-a1 + sqrt(a1^2-4*a2*(a0-(P_t/1e5))))-(T-273))
end
## Initial Conditions
function initial(phi_i, N, depth, depth_0, rho_s, rho_f)
# Plate Interface
plateInt = [0.2709, -9.8967e3]
distance_0 = (depth_0 - plateInt[2])/plateInt[1]
# Geometry
betaT = atan(plateInt[1])
v_plate = 0.037 / (365*24*3600)
v_x = v_plate * cos(betaT) # horizontal velocity
stressFit = [8693.11664571161, -365056544.708814]
dsigmadt = stressFit[1] * cos(betaT) * v_plate
#Temperature
tempDepth = [200/30000, (400*50000 - 600*20000)/30000]
dTdz = 500 / 25000 # geotherm
dTempdt = tempDepth[1] * plateInt[1] * v_x
const_vars = tempDepth[1] * (plateInt[1] * distance_0 + plateInt[2]) + tempDepth[2] + 273
Temp = dTdz .* depth
Temp = Temp .+ const_vars
# Stress overburden - polynomial coefficients are reverse in julia compared to matlab polyval
stressPoly = Polynomial([stressFit[2], stressFit[1]])
overburden = stressPoly(distance_0)
# porosity
porosity = phi_i .* ones(N)
# Stress in column
lithostatic = zeros(N)
bulkDensity = @. rho_s - (rho_s - rho_f)*porosity
lithostatic[1] = 9.81 * (depth[1])* bulkDensity[1]
for i in range(2,N, step = 1)
lithostatic[i] = 9.81 * (depth[i] - depth[i-1]) * bulkDensity[i] + lithostatic[i - 1]
end
hydrostatic = @. 9.81 * rho_f * (depth + depth_0)
# Total pore pressure
pressure_0 = (overburden * ones(N)) .+ lithostatic
# Excess Pore Pressure above hydrostatic (flow pressure)
pressure_i = pressure_0 .- hydrostatic
# Calc overpressure gradient
dp_lithdz = 9.81*(bulkDensity .- rho_f)
dp_lithdz = mean(dp_lithdz)
return pressure_i, pressure_0, hydrostatic, Temp, dTempdt, dsigmadt, dp_lithdz
end
## Diff Mats
function ddz(z0, zL, n, v)
## This is from MatLab MatMol group (2009)
## Second-order directional differencing for MOL
dz = (zL-z0)/(n-1)
r2fdz = 1 / (2*dz)
# Finite difference for positive v (indicates flow direction)
if v > 0
D = spdiagm(-2 => ones(n-2), -1 => -4*ones(n-1), 0 => 3 * ones(n))
## Boundaries
D[1, 1:3] = [-3, 4, -1]
D[2, 1:3] = [-1, 0, 1]
end
if v < 0
D = spdiagm(0 => -3*ones(n), 1 => 4*ones(n-1), 2 => -ones(n-2))
## Boundaries
D[n-1, (n-2):n] = [-1, 0, +1]
D[n, (n-2):n] = [1, -4, 3]
end
D = r2fdz*D
return D
end
function d2dz2(z0, zL, n)
## From Matlab MatMol group (2009)
## Four point second order approx. of the second derivative
## See the MatMol 2009 docs to see the derivation
dz = (zL-z0)/(n-1)
rdzs = 1 / (2*dz)
# interior
D = spdiagm(-1 => ones(n-1), 0 => -2*ones(n), 1 => ones(n-1))
# boundaries
D[1, 1:4] = [2, -5, 4, -1]
D[n, (n-3):n] = [-1, 4, -5, 2]
D = rdzs*D
return D
end
## Sparsity Pattern
function sparsity_pattern(N)
iD = 1:N
iSup = 1:N-1
iSupSup = 1:N-2
# Pressure section
P_p = sparse(iD, iD, ones(N)) + sparse(iSup, iSup .+ 1, ones(N-1), N, N) +
sparse(iSup.+1, iSup, ones(N-1), N, N) +
sparse(iSupSup, iSupSup .+ 2, ones(N-2), N, N)
P_p[N, N-2] = 1
P_p[N, N-3] = 1
P_T = sparse(iD, iD, ones(N))
P_phi = sparse(iD, iD, ones(N)) + sparse(iSup, iSup .+ 1, ones(N-1), N, N) +
sparse(iSupSup, iSupSup .+ 2, ones(N-2), N, N)
P_phi[N-1, N-2] = 1
P_phi[N, N-1] = 1
P_phi[N, N-2] = 1
P_m = P_T
Jac_P = sparse_hcat(P_p, P_T, P_phi, P_m)
# Temperature equation
T_p = sparse(iD, iD, [ones(N-1); 0]) + sparse(iSup, iSup .+ 1, ones(N-1), N, N) +
sparse(iSup.+1, iSup, [ones(N-2); 0], N, N)
T_T = T_p
T_phi = P_phi
T_m = sparse(iD, iD, [ones(N-1); 0])
Jac_T = sparse_hcat(T_p, T_T, T_phi, T_m)
# Porosity section is equal to the pressure section
Jac_m = sparse_hcat(P_T, P_T, spzeros(N,N), P_T)
# Combine all the sections
S = sparse_vcat(Jac_P, Jac_T, Jac_P, Jac_m)
S = convert(SparseMatrixCSC{Bool, Int64}, S)
return S
end
# Model
function dudt!(du, u, p, t)
# parameters
(m_i, M_H2O, A, c0, hydro, n_r, E_a, R, T0, litho, doubleSnake,
ρ_f, η, k_i, c_m, μ, n, dsigmadt, P_init_b, dz, dp_lithdz, alpha,
K, G, K_f, C_s, K_T, L, dTempdt, ϕ_0, D1, DD2) = p
# variables
p_ex = @view u[:, 1]
T = @view u[:, 2]
ϕ = @view u[:, 3]
m_d = @view u[:, 4]
dpdz = D1 * p_ex
dpdzz = DD2 * p_ex
dϕdz = D1 * ϕ
dTdzz = DD2 * T
### Dehydration
#dM_d/dt
@. du[:,4] .= (1/m_i) * M_H2O * A * c0 * (m_i - m_d) *
(IfElse.ifelse(DeltaG(T, p_ex + hydro) < 0, abs(DeltaG(T, p_ex + hydro))^n_r, 0)) *
exp((-E_a/R)*(1/T - 1/T0))
# last 1000 nodes do not react
du[end-1000:end, 4] .= 0
### Pressure
p_eff = @. litho + dsigmadt*t - (p_ex + hydro) # Effective Pressure
dPHIdt = @. (doubleSnake/ρ_f) * du[:,4] - (3/4)*ϕ* p_eff / η
#dP/dt
@. du[begin:end-1000, 1] = (k_i * c_m/μ) * ((ϕ[begin:end-1000] / ϕ_0)^n) * ((n/ϕ[begin:end-1000]) * dϕdz[begin:end-1000] * dpdz[begin:end-1000] + dpdzz[begin:end-1000]) +
c_m*(du[begin:end-1000,4]/ρ_f - dPHIdt[begin:end-1000]) + dsigmadt
@. du[end-1000+1:end, 1] = (k_i * c_m/μ) * ((ϕ[end-1000+1:end] / ϕ_0)^n) * ((n/ϕ[end-1000+1:end]) * dϕdz[end-1000+1:end] * dpdz[end-1000+1:end] + dpdzz[end-1000+1:end]) + dsigmadt
# Upper absorbing boundary
du[1,1] = @. (k_i*c_m/μ) * ((ϕ[1]/ϕ_0)^n) * (n/ϕ[1]) * dϕdz[1] * dpdz[1] + (c_m/η)*(3/4)*ϕ[1]*p_eff[1] + dsigmadt
# bottom boundary
pN = @. P_init_b + dz * dp_lithdz + dsigmadt*t
dpdzzN = (p_ex[end-1] - 2*p_ex[end] + pN)/ dz^2
du[end, 1] = @. (k_i*c_m/μ) * ((ϕ[end]/ ϕ_0)^n) * ((n/ϕ[end]) * dϕdz[end] * dpdz[end] + dpdzzN) +
c_m*(du[end,4]/ ρ_f - dPHIdt[end]) + dsigmadt
## Porosity
#dϕ/dt
@.du[1:end-1000,3] = -(alpha/(K+4*G/3))*dsigmadt +
(1/c_m - ϕ[1:end-1000]/K_f)*du[1:end-1000,1] +
dPHIdt[1:end-1000]
@. du[end-1000+1:end,3] = -(alpha/(K+4*G/3))*dsigmadt +
(1/c_m - ϕ[end-1000+1:end]/K_f)*du[end-1000+1:end,1]
## Temperature
# conduction only
@. du[:,2] .= (1/C_s) * (K_T * dTdzz - L*du[:, 4])
# Upper boundary
du[1,2] = -(L/C_s) * du[1,4]
# Lower boundary constant
du[end,2] = dTempdt
#done
nothing
end
### Parameters
m_i = 320 # Initial bound water
C_s = 2.5e6 # volumetric heat capacity
k_i = 1e-20 # background permeability
K_f = 2.9e9 # fluid bulk modulus
K_T = 2.25 # thermal conductivity
n = 3 # k-ϕ powerlaw exponent
v_p = 7000 # p-wave velocity
v_s = v_p / 2.35 # s-wave velocity
α = 0.2 # Biot parameter
η = 1e15 # bulk viscosity
μ = 2e-4 # fluid viscosity
ν = 0.3894 # Poisson's ratio
ϕ_i = 0.02 # initial porosity
ϕ_0 = 0.015 # reference porosity
ρ_f = 1e3 # fluid density
ρ_s = 3e3 # solid density
## Computed Parameters
ρ_b = (1-ϕ_i)*ρ_s + ϕ_i*ρ_f # Bulk Density
G = ρ_b * v_s^2 # Shear modulus
K_u = (2*(1+ν)*G)/(3*(1-2*ν)) # Undrained bulk modulus
# Drained bulk modulus
K = (1/(2*ϕ_i))*((K_f+K_u)*ϕ_i - α*(1+ϕ_i)*K_f +
(4*(α-1)*(ϕ_i-α)*ϕ_i*K_f*K_u +
(ϕ_i*(K_f+K_u)-α*(1+ϕ_i)*K_f)^2)^(1/2))
# Elastic Parameter
c_m = ((K_u - K)*(K + 4*G/3))/(α^2 * (K_u + 4*G/3))
## Dehydration Parameters
A = 201 # Specific surface area of RLM
c0 = 1.29e-20 # Ref. reaction constant
E_a = 2e4 # Activation energy
L = 500e3 # latent heat of Rxn.
M_H2O = 0.018 # H2O molar mass
n_r = 3.64 # Rxn. order
R = 1.987 # gas constant
T0 = 633 # Reference Temperature
doubleSnake = 0.1 # Plastic compaction parameter
### Discritization
N = 10000 # dehydrating nodes
NE = Int(0.1 * N) # Non dehydrating nodes
dz = 1000/N # distance between nodes in meters
N = 2*N + NE # Number of nodes in the column
depth = collect(dz:dz:dz*N)
### Set up problem
# Initial conditions
P_i, litho, hydro, T_i, dTdt, dsigmadt, dP_lithdz = initial(ϕ_i, N, depth, 28.345e3, ρ_s, ρ_f)
Phi_i = ϕ_i * ones(N)
M_di = zeros(N)
u0 = cat(P_i, T_i, Phi_i, M_di; dims = 2)
v = -1 # upwards flow
D1 = ddz(depth[1], depth[N], N, v)
DD2 = d2dz2(depth[1], depth[N], N)
# Parameters
p = (m_i, M_H2O, A, c0, hydro, n_r, E_a, R,
T0, litho, doubleSnake, ρ_f, η, k_i, c_m,
μ, n, dsigmadt, P_i[end], dz, dP_lithdz, α, K,
G, K_f, C_s, K_T, L, dTdt, ϕ_0, D1, DD2)
# Time -- need to be able to run for 120000 years
tspan = (0.0, 31536000 * 120000)
# Build sparsity pattern
jac_sparsity = sparsity_pattern(N)
f = ODEFunction(dudt! ; jac_prototype = float.(jac_sparsity))
# Make problem
dehydrate_prob = ODEProblem(f, u0, tspan, p)
sol = solve(dehydrate_prob, FBDF(), saveat = tspan[2]/2000, save_everystep = false, abstol = 1e-12, reltol = 1e-10, progress = true, maxiters = Int(1e6))
```

I would appreciate any help on how I can make the larger time spans possible to solve for. I figure if I cannot make the original model solve over these timespans, then the model I will develop after this will be likely worse.

Thanks!