Below I have the important parts of the scripts for a model that determines the geomorphic behaviors of a tidal marsh and lagoon system. Seen below there are 3 “behaviors” in a vector continuous callback function and everything works perfectly, the question is, “How can I determine which one of the outs in the condition function occurred for each model run?” so I can say “For this parameter set the marsh drowned, filled, or eroded.”
# Computational Parameters
years = 150
t0 = 1
tf = years*365*24*60*60
tspan = (t0,tf)
tyears = range(t0,years,years)
n = 40
# Model Parameters
β = 10^10
Bpeak = 2.500
Co = range(0,150*(10^-3),n)
χref = 0.15
Dmin = 0.0
g = 9.80171
ka = 2.0
ke = 0.1/(365*24*60*60)
ko = 0.001
λ = 0.0001
νGp = 0.0138
ϕ = 0.377
r = 1.4
R = range(0,15*(10^-3)/(365*24*60*60),n)
ρ = 1000
ρo = 1000
Tω = 12.5*(60*60)
τcr = 0.1
Uref = 8
ws = 0.5*10^-3
x = 10
# Initial Conditions
bf0 = 9000
dm0 = (0.7167*r-0.483)/2
L0 = 10000
trajectories = zeros(n,n,4)
times = zeros(n,n)
# n^2 Model Runs
for i in eachindex(R)
for j in eachindex(Co)
local sol = ode(i,j)
trajectories[i,j,:] = sol[:,length(sol)]
times[i,j] = length(sol)
end
end
This is the script I have that runs a model n^2 times to create a regime diagram of different geomorphological behaviors of the system I am interested in. Below is the callback function and solving function I am using. Everything works well.
function behavior()
function condition(out, u, t, integrator)
out[1] = u[3]-1.4 # Drown
out[2] = 0.99*u[4]-u[1] # Erode
out[3] = u[1]-0.01*u[4] # Fill
end
function affect!(integrator, idx)
if idx == 1
terminate!(integrator)
elseif idx == 2
terminate!(integrator)
elseif idx == 3
terminate!(integrator)
end
end
cb = VectorContinuousCallback(condition,affect!,3)
return cb
end
function ode(i,j)
p = [β Bpeak Co[j] χref Dmin g ka ke ko λ νGp ϕ r R[i] ρ ρo Tω τcr Uref ws x bf0]
prob_df0 = NonlinearProblem(lagoon_depth,[2],p)
df0 = solve(prob_df0, LevenbergMarquardt(), abstol=10^-28, reltol=10^-28, maxiters=10000)
u0 = [bf0 df0 dm0 L0]
prob = ODEProblem(system_equations,u0,tspan,p)
sol = solve(prob,Rosenbrock23(),abstol=10^-6,reltol=10^-6,callback = behavior(),saveat = range(t0,tf,years))
return sol
end