Select better solver method

Note that when moving to ComponentArrays.jl the testcase Hammel Loam NLAYER=103 showed oscillations with the previous solver method of Tsit5(). This lead to failure of the test.

Selecting an alternative solver method for stiff problems such as ImplicitEuler(autodiff=false), Rodas4P(autodiff=false), or Rosenbrock23(autodiff=false); seem to be much more efficient and result in better accuracy. Below code was used to to identify the issue.

# Define simulation model by reading in system definition and input data
model = SPAC("test/test-assets/Hammel-2001/input-files/",
                "Hammel_loam-NLayer-103-RESET=FALSE";
                simulate_isotopes = false);

# Prepare simulation by discretizing spatial domain
soil_discretization = discretize_soil(model.continuousIC.soil)
zspan = (maximum(model.soil_horizons.Upper_m),
            minimum(model.soil_horizons.Lower_m))
@assert zspan[1] == soil_discretization.Upper_m[1] "Soil discretization is not compatible with provided soil horizons"
@assert zspan[2] == soil_discretization.Lower_m[end] "Soil discretization is not compatible with provided soil horizons"
simulation = LWFBrook90.discretize(model;tspan=(0,2.5));
# Solve ODE:
# LWFBrook90.simulate!(simulation)
# solve_LWFB90(simulation.ODEProblem)
function norm_to_use(u::Real,            t) norm(u) end
function norm_to_use(u::ComponentVector, t)
    DiffEqBase.ODE_DEFAULT_NORM(reduce(vcat, [u.GWAT.mm,
                                                u.INTS.mm,
                                                u.INTR.mm,
                                                u.SNOW.mm,
                                                u.CC.mm,
                                                u.SNOWLQ.mm,
                                                u.SWATI.mm,
                                                u.RWU.mm,
                                                u.XYLEM.mm,
                                                u.TRANI.mm,
                                                u.aux.θ, u.aux.ψ, u.aux.K,
                                                u.accum
            ]),
        t)
end
function unstable_check_function(dt,u::ComponentVector,p,t)
    any(isnan, u.SWATI.mm) # Do it for SWATI (only amt)
end
# source: https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Fixed-Stepsize-Usage
# Force the same time steps as with sol_working:
# If tstops is set without dt, then the algorithm will step directly to each value in tstops
tspan = simulation.ODEProblem.tspan
@time sol_current = solve(simulation.ODEProblem,
            # ImplicitEuler(autodiff=false);  # for stiff problems (~6.7s)
            # Rodas4P(autodiff=false);        # for stiff problems (~3.2s)
            Rosenbrock23(autodiff=false);   # for stiff problems   (~2.3s)
            # Tsit5(); # Tsit5 recommended for non-stiff problems
            progress = true,
            adaptive = true, internalnorm = norm_to_use # fix adaptivity norm for NAs
            # save_everystep = true,
            # tstops = sol_working.t, #[0.223, 0.4],
            # tstops = tspan[1]:0.00001:tspan[2], adaptive = false
    );
plot(sol_current.t, reduce(hcat, [sol_current.u[t].SWATI.mm[idx] for t in eachindex(sol_current)])', color = "orange", line = :dot)

fabern / LWFBrook90.jl

Select better solver method #44