diff --git a/docs/make.jl b/docs/make.jl index 574afa472..e6aee0a27 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -90,7 +90,8 @@ makedocs(modules = [PositiveIntegrators, Base.get_extension(PositiveIntegrators, "Experimental order of convergence" => "convergence.md", "NPZD model" => "npzd_model_benchmark.md", "Robertson problem" => "robertson_benchmark.md", - "Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md" + "Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md", + "Spatially heterogeneous diffusion equation" => "diffusion_benchmark.md" ], "Troubleshooting, FAQ" => "faq.md", "API reference" => "api_reference.md", diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md index 3bb67575d..9b8fb3089 100644 --- a/docs/src/api_reference.md +++ b/docs/src/api_reference.md @@ -16,12 +16,15 @@ PDSProblem ```@docs prob_pds_bertolazzi prob_pds_brusselator +prob_pds_diffusion +prob_ode_diffusion prob_pds_linmod prob_pds_linmod_inplace prob_pds_minmapk prob_pds_nonlinmod prob_pds_npzd prob_pds_robertson +prob_pds_saceirqd prob_pds_sir prob_pds_stratreac prob_ode_stratreac_scaled diff --git a/docs/src/diffusion_benchmark.md b/docs/src/diffusion_benchmark.md new file mode 100644 index 000000000..210dd345e --- /dev/null +++ b/docs/src/diffusion_benchmark.md @@ -0,0 +1,250 @@ +# [Benchmark: Solution of the Diffusion problem](@id benchmark-diffusion) + +We consider the test problem [`prob_pds_diffusion`](@ref) of the spatially heterogeneous diffusion equation to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl), especially on larger-scale problems. + +```@example DIFFU +using OrdinaryDiffEqFIRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqSDIRK +using PositiveIntegrators + +# select spatially heterogeneous diffusion problem +prob = prob_pds_diffusion +nothing # hide +``` + +To keep the following code as clear as possible, we define a helper function `diffusion_plot` that we use for plotting. + +```@example DIFFU +using Plots + +function diffusion_plot(sol, ref_sol=nothing, title_str="") + N = length(sol.u[1]) + x = range(0.0, 1.0; length=N) + colors = palette(:default)[1:3]' + + lbls = ["t=$(sol.t[1])" "t=$(sol.t[2])" "t=$(sol.t[3])"] + + if ref_sol !== nothing + p1 = plot(x, ref_sol.u[1], linestyle=:dash, linewidth=2, color=colors[2]) + p2 = plot(x, ref_sol.u[2], linestyle=:dash, linewidth=2, color=colors[2]) + p3 = plot(x, ref_sol.u[3], linestyle=:dash, linewidth=2, color=colors[2]) + else + p1= () + p2 = () + p3 = () + end + + plot!(p1, x, sol.u[1], linewidth=0.5, color=colors[1], title=lbls[1]) + plot!(p2, x, sol.u[2], linewidth=0.5, color=colors[1], title=lbls[2]) + plot!(p3, x, sol.u[3], linewidth=0.5, color=colors[1], title=lbls[3]) + + p = plot(p1, p2, p3, layout=(1,3), size=(1200,350), plot_title=title_str, xlabel="x", ylabel="u", legend=false) + + return p +end + + +nothing # hide +``` + +## Work-Precision diagrams + +In the following we show several work-precision diagrams, which compare different methods with respect to computing time and the respective error. + +Since the diffusion problem is stiff, we need to use a suited implicit scheme to compute a reference solution, see the [solver guide](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/#Stiff-Problems). + +```@example DIFFU +# select solver to compute reference solution +alg_ref = RadauIIA5() +nothing # hide +``` + +We use the functions [`work_precision_adaptive`](@ref) and [`work_precision_adaptive!`](@ref) to compute the data for the diagrams. +Furthermore, the following absolute and relative tolerances are used. + +```@example DIFFU +# set absolute and relative tolerances +abstols = 1.0 ./ 10.0 .^ (2:1:10) +reltols = abstols .* 10.0 +nothing # hide +``` + +### Relative maximum error at the final time + +In this section the chosen error is the relative maximum error at the final time ``t = 60.0``. + +```@example DIFFU +# select relative maximum error at the end of the problem's time span. +compute_error = rel_max_error_tend +nothing # hide +``` + +We start with a comparison of different adaptive MPRK schemes. +```@example DIFFU +# choose methods to compare +algs = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0); + MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0); + MPDeC(2); MPDeC(3); MPDeC(4); MPDeC(5); MPDeC(6); MPDeC(7); MPDeC(8); MPDeC(9); MPDeC(10)] +labels = ["MPRK22(0.5)"; "MPPRK22(2/3)"; "MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; + "MPRK43I(1.0, 0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"; + "MPDeC(2)"; "MPDeC(3)"; "MPDeC(4)"; "MPDeC(5)"; "MPDeC(6)"; "MPDeC(7)"; "MPDeC(8)"; "MPDeC(9)"; "MPDeC(10)"] + +# compute work-precision data +wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, labels; title = "Diffusion benchmark", legend = :outerright, + color = permutedims([ + repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)..., repeat([5], 9)... ]), + xlims = (4*10^-11, 5*10^-2), xticks = 10.0 .^ (-10:1:-2), + ylims = (2*10^-3, 2*10^1), yticks = 10.0 .^ (-3:1:1), minorticks = 10) +``` + +For comparisons with other schemes we choose `MPRK22(1.0)`, `SSPMPRK22(0.5, 1.0)`, `MPRK43I(0.5, 0.75)` and `MPDeC(10)`. + +```@example DIFFU +# compute reference solution for plotting +saveat = (0.05, 1.0, 60.0) +ref_sol = solve(prob, RadauIIA5(); abstol = 1e-14, reltol = 1e-13, saveat=saveat); + +# compute solutions with loose tolerances +abstol = 1e-2 +reltol = 1e-1 +sol_MPRK22_1 = solve(prob, MPRK22(1.0); abstol, reltol, saveat=saveat) +sol_SSPMPRK22 = solve(prob, SSPMPRK22(0.5, 1.0); abstol, reltol, saveat=saveat) +sol_MPRK43 = solve(prob, MPRK43I(0.5, 0.75); abstol, reltol, saveat=saveat) +sol_MPDeC10 = solve(prob, MPDeC(10); abstol, reltol, saveat=saveat) + + +p1 = diffusion_plot(sol_MPRK22_1, ref_sol, "MPRK22(1.0)"); +p2 = diffusion_plot(sol_SSPMPRK22, ref_sol, "SSPMPRK(0.5, 1.0)"); +p3 = diffusion_plot(sol_MPRK43, ref_sol, "MPRK43I(0.5, 0.75)"); +p4 = diffusion_plot(sol_MPDeC10, ref_sol, "MPDeC(10)"); +plot(p1, p2, p3, p4, layout=(4,1), size=(1200, 1400)) +``` + +Next, we compare these four schemes with a selection of second- and third-order stiff solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Due to the special structure of the problem, which involves tridiagonal matrices, the classical solvers become inefficient when applied to the PDS formulation. Therefore, for these solvers we instead solve the problem in its original form directly. + +```@example DIFFU +# set the problem for the classical solvers +prob_classic = prob_ode_diffusion + + +# select reference MPRK methods +algs1 = [MPRK22(1.0); SSPMPRK22(0.5, 1.0); MPRK43I(0.5, 0.75); MPDeC(10)] +labels1 = ["MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; "MPRK43I(0.5,0.75)"; "MPDeC(10)"] + +# select methods from OrdinaryDiffEq +algs2 = [TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23()] +labels2 = ["TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23"] + +# compute work-precision data +wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) +# add work-precision data +work_precision_adaptive!(wp, prob_classic, algs2, labels2, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, [labels1; labels2]; title = "Diffusion benchmark", legend = :outerright, + color = permutedims([repeat([1], 2)..., 3, 5, repeat([6], 3)..., repeat([7], 4)...]), + xlims = (10^-11, 10^-1), xticks = 10.0 .^ (-11:1:-1), + ylims = (10^-3, 10^3), yticks = 10.0 .^ (-3:1:3), minorticks = 10) +``` + +In addition, we compare the selected MPRK schemes to some [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) of higher order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). + +```@example DIFFU +algs3 = [Rodas5P(); Rodas4P(); RadauIIA5()] +labels3 = ["Rodas5P"; "Rodas4P"; "RadauIIA5"] + +# compute work-precision data +wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) +# add work-precision data with isoutofdomain = isnegative +work_precision_adaptive!(wp, prob_classic, algs3, labels3, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, [labels1; labels3]; title = "Diffusion benchmark", legend = :bottomleft, + color = permutedims([repeat([1],2)..., 3, 5, repeat([4], 2)..., 6]), + xlims = (10^-12, 10^-2), xticks = 10.0 .^ (-12:1:-2), + ylims = (10^-3, 10^2), yticks = 10.0 .^ (-3:1:2), minorticks = 10) +``` + + +### Relative maximum error over all time steps + +In this section we do not compare the relative maximum errors at the final time ``t = 60.0}``, but the relative maximum errors over all time steps. + +```@example DIFFU +# select relative maximum error at the end of the problem's time span. +compute_error = rel_max_error_overall +nothing # hide +``` + +First, we compare different MPRK schemes. + +```@example DIFFU +# compute work-precision data +wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, labels; title = "Diffusion benchmark", legend = :outerright, + color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)..., repeat([5], 9)... ]), + xlims = (10^-9, 10^0), xticks = 10.0 .^ (-9:1:2), + ylims = (10^-3, 2*10^1), yticks = 10.0 .^ (-3:1:1), minorticks = 10) +``` + +We choose the second-order scheme `MPRK22(1.0)` and the third-order scheme `MPRK43I(0.5, 0.75)` for comparison with solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). + +```@example DIFFU +# select reference MPRK methods +algs1 = [MPRK22(1.0); MPRK43I(0.5, 0.75)] +labels1 = ["MPRK22(1.0)"; "MPRK43I(0.5,0.75)"] + +# compute work-precision data +wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) +# add work-precision data with isoutofdomain = isnegative +work_precision_adaptive!(wp, prob_classic, algs2, labels2, abstols, reltols, alg_ref; adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, [labels1; labels2]; title = "Diffusion benchmark", legend = :bottomleft, + color = permutedims([1, 3, repeat([5], 3)..., repeat([6], 4)...]), + xlims = (10^-9, 10^0), xticks = 10.0 .^ (-9:1:3), + ylims = (10^-3, 10^3), yticks = 10.0 .^ (-3:1:3), minorticks = 10) +``` + +Finally, we compare `MPRK43I(0.5, 0.75)` and `MPRK22(1.0)` to [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) of higher order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). + +```@example DIFFU +# compute work-precision data +wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; + adaptive_ref = true, compute_error) +# add work-precision data with isoutofdomain = isnegative +work_precision_adaptive!(wp, prob_classic, algs3, labels3, abstols, reltols, alg_ref;adaptive_ref = true, compute_error) + +# plot work-precision diagram +plot(wp, [labels1; labels3]; title = "Diffusion benchmark", legend = :topright, + color = permutedims([1, 3, repeat([4], 2)..., 5]), + xlims = (10^-10, 10^0), xticks = 10.0 .^ (-10:1:0), + ylims = (10^-3, 2*10^1), yticks = 10.0 .^ (-3:1:1), minorticks = 10) +``` + +## Package versions + +These results were obtained using the following versions. +```@example DIFFU +using InteractiveUtils +versioninfo() +println() + +using Pkg +Pkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", + "OrdinaryDiffEqFIRK", "OrdinaryDiffEqRosenbrock", + "OrdinaryDiffEqSDIRK"], + mode=PKGMODE_MANIFEST) +nothing # hide +``` diff --git a/src/PDSProblemLibrary.jl b/src/PDSProblemLibrary.jl index a706bf5f1..768f7e384 100644 --- a/src/PDSProblemLibrary.jl +++ b/src/PDSProblemLibrary.jl @@ -597,3 +597,260 @@ There are two independent linear invariants, e.g. ``u_1+u_4+u_6=1.75`` and ``u_2 prob_pds_minmapk = PDSProblem(P_minmapk, D_minmapk, u0, tspan; std_rhs = f_minmapk, linear_invariants = @SMatrix[1.0 0.0 0.0 1.0 0.0 1.0; 0.0 1.0 1.0 1.0 1.0 0.0]) + +# SACEIRQD Covid-19 model +function f_saceirqd(u, p, t) + Npop = 6.046e7 + alpha = 0.0194 + beta = 7.567 + mu = 2.278e-6 + eta = 9.180e-7 + sigma = 1.4633e-3 + tau = 1.109e-4 + xi = 0.263 + gamma = 0.021 + delta = 0.077 + lambda = 6.2800e-04 + Kd = 0.0013 + + # infection-like term + inf_term = (beta * u[5] + sigma * u[2]) / Npop + eta + + return @SVector [-alpha * u[1] - inf_term * u[1]; + xi * u[4] - tau * u[2]; + alpha * u[1] - mu * u[3]; + inf_term * u[1] + mu * u[3] - (gamma + xi) * u[4]; + tau * u[2] + gamma * u[4] - delta * u[5]; + lambda * u[7]; + delta * u[5] - (lambda + Kd) * u[7]; + Kd * u[7]] +end + +function P_saceirqd(u, p, t) + Npop = 6.046e7 + alpha = 0.0194 + beta = 7.567 + mu = 2.278e-6 + eta = 9.180e-7 + sigma = 1.4633e-3 + tau = 1.109e-4 + xi = 0.263 + gamma = 0.021 + delta = 0.077 + lambda = 6.2800e-04 + Kd = 0.0013 + + # P[i,j] is flux from compartment j -> i + P41 = (u[1] / Npop) * (beta * u[5] + sigma * u[2]) + eta * u[1] + + return @SMatrix [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 xi*u[4] 0.0 0.0 0.0 0.0; + alpha*u[1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0; + P41 0.0 mu*u[3] 0.0 0.0 0.0 0.0 0.0; + 0.0 tau*u[2] 0.0 gamma*u[4] 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0 lambda*u[7] 0.0; + 0.0 0.0 0.0 0.0 delta*u[5] 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0 Kd*u[7] 0.0] +end + +# initial value (from SirdTest with sost = 1e-10) +u0_saceirqd = @SVector [6.046e7 - (4e-10 + 3.0); 1e-10; 1e-10; 1.0; 1.0; 1e-10; 1.0; + 1e-10] + +tspan_saceirqd = (0.0, 180.0) + +""" + prob_pds_saceirqd + +Positive and conservative autonomous nonlinear PDS +```math +\\begin{aligned} +u_1' &= -\\alpha u_1 - \\left(\\frac{\\beta}{N_{\\mathrm{pop}}}u_5 + \\frac{\\sigma}{N_{\\mathrm{pop}}}u_2 + \\eta\\right)u_1,\\\\ +u_2' &= \\xi u_4 - \\tau u_2,\\\\ +u_3' &= \\alpha u_1 - \\mu u_3,\\\\ +u_4' &= \\left(\\frac{\\beta}{N_{\\mathrm{pop}}}u_5 + \\frac{\\sigma}{N_{\\mathrm{pop}}}u_2 + \\eta\\right)u_1 + + \\mu u_3 - (\\gamma + \\xi)u_4,\\\\ +u_5' &= \\tau u_2 + \\gamma u_4 - \\delta u_5,\\\\ +u_6' &= \\lambda u_7,\\\\ +u_7' &= \\delta u_5 - (\\lambda + K_d)u_7,\\\\ +u_8' &= K_d u_7, +\\end{aligned} +``` + +with constants + +```math +\\begin{aligned} +N_{\\mathrm{pop}} &= 6.046\\cdot 10^{7}, \\\\ +\\alpha &= 0.0194, \\\\ +\\beta &= 7.567, \\\\ +\\mu &= 2.278\\cdot 10^{-6}, \\\\ +\\eta &= 9.180\\cdot 10^{-7}, \\\\ +\\sigma &= 1.4633\\cdot 10^{-3}, \\\\ +\\tau &= 1.109\\cdot 10^{-4}, \\\\ +\\xi &= 0.263, \\\\ +\\gamma &= 0.021, \\\\ +\\delta &= 0.077, \\\\ +\\lambda &= 6.2800\\cdot 10^{-4}, \\\\ +K_d &= 0.0013. +\\end{aligned} +``` + +The initial value is ``\\mathbf{u}_0 = (6.046\\cdot 10^{7}-(4\\cdot 10^{-10}+3),\\,10^{-10},\\,10^{-10},\\,1,\\,1,\\,10^{-10},\\,1,\\,10^{-10})^T`` and the time domain ``(0.0, 180.0)``. +There is one independent linear invariant, namely total population ``u_1+u_2+u_3+u_4+u_5+u_6+u_7+u_8 = N_{\\mathrm{pop}}``. + +## References + +- D. Sen and D. Sen. + "Use of a modified SIRD model to analyze COVID-19 data." + Industrial & Engineering Chemistry Research 60(11) (2021): 4251–4260. + [DOI: 10.1021/acs.iecr.0c04754](https://doi.org/10.1021/acs.iecr.0c04754) +- Giuseppe Izzo, Eleonora Messina, Mario Pezzella, and Antonia Vecchio. + "Modified Patankar Linear Multistep Methods for Production-Destruction Systems." + Journal of Scientific Computing 102 (2025): 87. + [DOI: 10.1007/s10915-025-02804-5](https://doi.org/10.1007/s10915-025-02804-5) +""" +prob_pds_saceirqd = ConservativePDSProblem(P_saceirqd, u0_saceirqd, + tspan_saceirqd, std_rhs = f_saceirqd, + linear_invariants = @SMatrix[1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0]) + +# diffusion problem +function f_diffusion!(du, u, p, t) + K = p.K + invdx2 = p.invdx2 + N = length(u) + + @inbounds begin + du[1] = (K[2] * u[2] - K[1] * u[1]) * invdx2 + + for i in 2:(N - 1) + du[i] = (K[i - 1] * u[i - 1] + K[i + 1] * u[i + 1] - 2 * K[i] * u[i]) * invdx2 + end + + du[N] = (K[N - 1] * u[N - 1] - K[N] * u[N]) * invdx2 + end + return nothing +end + +function P_diffusion!(P::Tridiagonal, u, p, t) + K = p.K + invdx2 = p.invdx2 + N = length(u) + + fill!(P.dl, zero(eltype(P))) + fill!(P.d, zero(eltype(P))) + fill!(P.du, zero(eltype(P))) + + @inbounds for i in 1:(N - 1) + P.du[i] = K[i + 1] * u[i + 1] * invdx2 + P.dl[i] = K[i] * u[i] * invdx2 + end + + return nothing +end + +N_diffusion = 2000 +L_diffusion = 1.0 +dx_diffusion = L_diffusion / N_diffusion +invdx2_diffusion = 1.0 / (dx_diffusion^2) +x_diffusion = collect(range(dx_diffusion / 2, L_diffusion - dx_diffusion / 2, + length = N_diffusion)) + +D0 = 1e-2 +kfun = x -> 1e-5 + + (x - 2 * L_diffusion / 3) .^ 2 .* D0 .* + atan(0.5 * (2 * x - L_diffusion * 1.5 * 2)) ./ + (0.5 * (2 * x - L_diffusion * 1.5 * 2)) +K_diffusion = kfun.(x_diffusion) + +f0 = x -> 2 * (1 - sin(pi * (x * pi / 2 - 0.25))^2) +u0_diffusion = f0.(x_diffusion) + +tspan_diffusion = (0.0, 60.0) + +p_diffusion = (K = K_diffusion, invdx2 = invdx2_diffusion) + +p_prototype_diffusion = Tridiagonal(zeros(eltype(u0_diffusion), N_diffusion - 1), + zeros(eltype(u0_diffusion), N_diffusion), + zeros(eltype(u0_diffusion), N_diffusion - 1)) + +#TODO Docs must be revised. What is K? What is f? +""" + prob_pds_diffusion + +Positive and conservative autonomous nonlinear production–destruction system +obtained from a finite-volume discretization of a one-dimensional diffusion equation +with spatially varying diffusion coefficient. + +```math +\\begin{aligned} +u_i' &= \\sum_{j=1}^{N} \\bigl( P_{ij}(u) - P_{ji}(u) \\bigr), \\qquad i = 1,\\dots,N,\\\\ +P_{i,i+1}(u) &= \\frac{1}{\\Delta x^2} K_{i+1} u_{i+1},\\qquad +P_{i+1,i}(u) = \\frac{1}{\\Delta x^2} K_i u_i, +\\end{aligned} +``` + +with ``P_{i,j}(u)=0`` otherwise. + +The grid consists of N = 2000 cells with width ``\\Delta x = 10^{-2}`` +and centers ``x_i = (i-\\tfrac12)\\Delta x`` (``L = 1``). +The initial value is ``\\mathbf{u}_0 = (u_1^0,\\dots,u_N^0)^T`` with +``u_i^0 = f(x_i)``, and the time domain ``(0.0, 60.0)``. + +There is one independent linear invariant, namely +``\\sum_{i=1}^{N} u_i = \\text{const}.`` + +## References + +- Giuseppe Izzo, Eleonora Messina, Mario Pezzella, and Antonia Vecchio. + "Modified Patankar Linear Multistep Methods for Production-Destruction Systems." + Journal of Scientific Computing* 102 (2025): 87. + [DOI: 10.1007/s10915-025-02804-5](https://doi.org/10.1007/s10915-025-02804-5) +""" +prob_pds_diffusion = ConservativePDSProblem(P_diffusion!, + u0_diffusion, + tspan_diffusion, + p_diffusion; + p_prototype = p_prototype_diffusion, + std_rhs = ODEFunction(f_diffusion!; + jac_prototype = p_prototype_diffusion), + linear_invariants = ones(1, N_diffusion)) + +#TODO Docs must be revised. What is K? What is f? f_diffusion! +""" + prob_ode_diffusion + +Positive and conservative autonomous nonlinear system of ordinary differential equations +obtained from a finite-volume discretization of a one-dimensional diffusion equation +with spatially varying diffusion coefficient. + +```math +\\begin{aligned} +u_i' &= \\sum_{j=1}^{N} \\bigl( P_{ij}(u) - P_{ji}(u) \\bigr), \\qquad i = 1,\\dots,N,\\\\ +P_{i,i+1}(u) &= \\frac{1}{\\Delta x^2} K_{i+1} u_{i+1},\\qquad +P_{i+1,i}(u) = \\frac{1}{\\Delta x^2} K_i u_i, +\\end{aligned} +``` + +with ``P_{i,j}(u)=0`` otherwise. + +The grid consists of N = 2000 cells with width ``\\Delta x = 10^{-2}`` +and centers ``x_i = (i-\\tfrac12)\\Delta x`` (``L = 1``). +The initial value is ``\\mathbf{u}_0 = (u_1^0,\\dots,u_N^0)^T`` with +``u_i^0 = f(x_i)``, and the time domain ``(0.0, 60.0)``. + +There is one independent linear invariant, namely +``\\sum_{i=1}^{N} u_i = \\text{const}.`` + +## References + +- Giuseppe Izzo, Eleonora Messina, Mario Pezzella, and Antonia Vecchio. + "Modified Patankar Linear Multistep Methods for Production-Destruction Systems." + Journal of Scientific Computing* 102 (2025): 87. + [DOI: 10.1007/s10915-025-02804-5](https://doi.org/10.1007/s10915-025-02804-5) +""" +prob_ode_diffusion = ODEProblem(ODEFunction(f_diffusion!; + jac_prototype = p_prototype_diffusion), + u0_diffusion, + tspan_diffusion, + p_diffusion) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index eba70eb39..a521a6ea0 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -6,14 +6,15 @@ using Statistics: median using SparseArrays: SparseArrays, AbstractSparseMatrix, SparseMatrixCSC, issparse, nonzeros, nzrange, rowvals, spdiagm -using StaticArrays: SVector, SMatrix, StaticArray, StaticMatrix, @SVector, @SMatrix, MMatrix +using StaticArrays: SVector, SMatrix, StaticArray, StaticMatrix, @SVector, @SMatrix, + MMatrix using FastBroadcast: @.. using MuladdMacro: @muladd using Reexport: @reexport -@reexport using SciMLBase: ODEProblem, init, solve +@reexport using SciMLBase: ODEProblem, ODEFunction, init, solve using SciMLBase: AbstractODEFunction, NullParameters, FullSpecialize, isinplace @@ -30,7 +31,7 @@ import SciMLBase: interp_summary using OrdinaryDiffEqCore: OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, False, - _vec + _vec, @cache import OrdinaryDiffEqCore: alg_order, isfsal, calculate_residuals, calculate_residuals!, alg_cache, get_tmp_cache, @@ -47,11 +48,13 @@ export ConservativePDSFunction, ConservativePDSProblem export MPE, MPRK22, MPRK43I, MPRK43II export SSPMPRK22, SSPMPRK43 export MPDeC +export MPLM22, MPLM33, MPLM43, MPLM54, MPLM75, MPLM106 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, prob_pds_bertolazzi, prob_pds_npzd, prob_pds_stratreac, prob_pds_minmapk export prob_ode_stratreac_scaled, linear_invariants_stratreac_scaled +export prob_pds_saceirqd, prob_pds_diffusion, prob_ode_diffusion export isnegative, isnonnegative export work_precision_adaptive, work_precision_adaptive!, work_precision_fixed, @@ -77,6 +80,9 @@ include("sspmprk.jl") # MPDeC methods include("mpdec.jl") +# MPLM methods +include("mplm.jl") + # interpolation for dense output include("interpolation.jl") diff --git a/src/interpolation.jl b/src/interpolation.jl index d3be83e31..18b6216c3 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -47,7 +47,8 @@ const MPRKCaches = Union{MPEConstantCache, MPECache, MPRK43ConstantCache, MPRK43Cache, SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK43ConstantCache, SSPMPRK43Cache, - MPDeCConstantCache, MPDeCCache, MPDeCConservativeCache} + MPDeCConstantCache, MPDeCCache, MPDeCConservativeCache, + MPLM22oopCache, MPLM33oopCache, MPLM43oopCache} function interp_summary(::Type{cacheType}, dense::Bool) where {cacheType <: MPRKCaches} diff --git a/src/mplm.jl b/src/mplm.jl new file mode 100644 index 000000000..b4257bf45 --- /dev/null +++ b/src/mplm.jl @@ -0,0 +1,2265 @@ +######################################################################################## +### Structs and caches ################################################################# +######################################################################################## +abstract type MPLMMutableCache <: OrdinaryDiffEqMutableCache end +get_fsalfirstlast(cache::MPLMMutableCache, rate_prototype) = (nothing, nothing) + +#### MPLM22 ############################################################################ +struct MPLM22{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM22) = 2 +isfsal(::MPLM22) = false +alg_extrapolates(alg::MPLM22) = false # TODO: Should probably be false + +#TODO: Check if OrdinaryDiffEqConstantCache is correct supertype +@cache mutable struct MPLM22oopCache{uType, T} <: OrdinaryDiffEqConstantCache + uprevprev::uType + step::Int + small_constant::T +end + +@cache mutable struct MPLM22Cache{uType, dType, T, PType, F} <: MPLMMutableCache + uprevprev::uType + step::Int + small_constant::T + tmp::uType + P::PType + P2::PType + D::dType + D2::dType + σ::uType + linsolve::F +end + +get_tmp_cache(integrator, ::MPLM22, cache::MPLMMutableCache) = (cache.σ,) + +function MPLM22(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM22(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM22, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + MPLM22oopCache(u, 1, alg.small_constant_function(uEltypeNoUnits)) +end + +# In-place +function alg_cache(alg::MPLM22, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + uprevprev = zero(u) + uprevprev[1] = 111.0 + uprevprev[2] = 222.0 + + step = 1 + small_constant = alg.small_constant_function(uEltypeNoUnits) + tmp = zero(u) + P = p_prototype(u, f) + # We use P2 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + P2 = p_prototype(u, f) + σ = zero(u) + + if f isa ConservativePDSFunction + # The right hand side of the linear system is always uprev. But using + # tmp instead of uprev for the rhs we allow `alias_b=true`. uprev must + # not be altered, since it is needed to compute the adaptive time step + # size. + linprob = LinearProblem(P2, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPLM22Cache(uprevprev, step, small_constant, tmp, P, P2, nothing, nothing, σ, + linsolve) + elseif f isa PDSFunction + linprob = LinearProblem(P2, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPLM22Cache(uprevprev, step, small_constant, tmp, P, P2, + similar(u), # D + similar(u), # D2 + σ, + linsolve) + else + throw(ArgumentError("MPLM22 can only be applied to production-destruction systems")) + end +end + +#### MPLM33 ############################################################################ +struct MPLM33{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM33) = 3 +alg_extrapolates(alg::MPLM33) = true # TODO: Should probably be false + +@cache mutable struct MPLM33oopCache{uType, PType, dType, T, T2} <: + OrdinaryDiffEqConstantCache + uprevprev::uType + uprev3::uType + P2::PType + P3::PType + d2::dType + d3::dType + αβ::NTuple{6, T} + step::Int + small_constant::T2 +end + +function MPLM33(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM33(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM33, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + α1 = zero(uEltypeNoUnits) + α2 = zero(uEltypeNoUnits) + α3 = one(uEltypeNoUnits) + β1 = 9 / 4 * one(uEltypeNoUnits) + β2 = zero(uEltypeNoUnits) + β3 = 3 / 4 * one(uEltypeNoUnits) + αβ = (α1, α2, α3, β1, β2, β3) + + # TODO: This is currently necessary to get the correct type of P (d is of type rateType) + P, d = evaluate_pds(f, u, p, t) + # TODO: integrator_stats_nf = 1 + + MPLM33oopCache(u, u, P, P, d, d, αβ, 1, alg.small_constant_function(uEltypeNoUnits)) +end + +#### MPLM43 ############################################################################ +struct MPLM43{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM43) = 3 +alg_extrapolates(alg::MPLM43) = true # TODO: Should probably be false + +@cache mutable struct MPLM43oopCache{uType, PType, dType, T, T2} <: + OrdinaryDiffEqConstantCache + uprevprev::uType + uprev3::uType + uprev4::uType + P2::PType + P3::PType + P4::PType + d2::dType + d3::dType + d4::dType + αβ::NTuple{8, T} + step::Int + small_constant::T2 +end + +function MPLM43(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM43(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM43, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + # TODO: This is currently necessary to get the correct type of P (d is of type rateType) + P, d = evaluate_pds(f, u, p, t) + # TODO: integrator_stats_nf = 1 + + α1 = 1 / 4 * one(uEltypeNoUnits) + α2 = zero(uEltypeNoUnits) + α3 = 3 / 4 * one(uEltypeNoUnits) + α4 = zero(uEltypeNoUnits) + β1 = 35 / 18 * one(uEltypeNoUnits) + β2 = 1 / 3 * one(uEltypeNoUnits) + β3 = zero(uEltypeNoUnits) + β4 = 2 / 9 * one(uEltypeNoUnits) + αβ = (α1, α2, α3, α4, β1, β2, β3, β4) + MPLM43oopCache(u, u, u, P, P, P, d, d, d, αβ, 1, + alg.small_constant_function(uEltypeNoUnits)) +end + +#### MPLM54 ############################################################################ +struct MPLM54{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM54) = 4 +alg_extrapolates(alg::MPLM54) = true # TODO: Should probably be false + +@cache mutable struct MPLM54oopCache{uType, PType, dType, T, T2} <: + OrdinaryDiffEqConstantCache + uprevprev::uType + uprev3::uType + uprev4::uType + uprev5::uType + P2::PType + P3::PType + P4::PType + P5::PType + d2::dType + d3::dType + d4::dType + d5::dType + αβ::NTuple{10, T} + step::Int + small_constant::T2 +end + +function MPLM54(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM54(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM54, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + # TODO: This is currently necessary to get the correct type of P (d is of type rateType) + P, d = evaluate_pds(f, u, p, t) + # TODO: integrator_stats_nf = 1 + + α1 = zero(uEltypeNoUnits) + α2 = zero(uEltypeNoUnits) + α3 = zero(uEltypeNoUnits) + α4 = zero(uEltypeNoUnits) + α5 = one(uEltypeNoUnits) + β1 = 225 / 96 * one(uEltypeNoUnits) + β2 = zero(uEltypeNoUnits) + β3 = 50 / 96 * one(uEltypeNoUnits) + β4 = 200 / 96 * one(uEltypeNoUnits) + β5 = 5 / 96 * one(uEltypeNoUnits) + αβ = (α1, α2, α3, α4, α5, β1, β2, β3, β4, β5) + MPLM54oopCache(u, u, u, u, P, P, P, P, d, d, d, d, αβ, 1, + alg.small_constant_function(uEltypeNoUnits)) +end + +#### MPLM75 ############################################################################ +struct MPLM75{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM75) = 5 +alg_extrapolates(alg::MPLM75) = true # TODO: Should probably be false + +@cache mutable struct MPLM75oopCache{uType, PType, dType, T, T2} <: + OrdinaryDiffEqConstantCache + uprevprev::uType + uprev3::uType + uprev4::uType + uprev5::uType + uprev6::uType + uprev7::uType + P2::PType + P3::PType + P4::PType + P5::PType + P6::PType + P7::PType + d2::dType + d3::dType + d4::dType + d5::dType + d6::dType + d7::dType + αβ::NTuple{14, T} + step::Int + small_constant::T2 +end + +function MPLM75(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM75(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM75, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + # TODO: This is currently necessary to get the correct type of P (d is of type rateType) + P, d = evaluate_pds(f, u, p, t) + # TODO: integrator_stats_nf = 1 + + α1 = zero(uEltypeNoUnits) + α2 = zero(uEltypeNoUnits) + α3 = zero(uEltypeNoUnits) + α4 = zero(uEltypeNoUnits) + α5 = zero(uEltypeNoUnits) + α6 = zero(uEltypeNoUnits) + α7 = one(uEltypeNoUnits) + β1 = 12 / 5 * one(uEltypeNoUnits) + β2 = zero(uEltypeNoUnits) + β3 = 197 / 720 * one(uEltypeNoUnits) + β4 = 701 / 360 * one(uEltypeNoUnits) + β5 = 43 / 30 * one(uEltypeNoUnits) + β6 = 107 / 360 * one(uEltypeNoUnits) + β7 = 467 / 720 * one(uEltypeNoUnits) + αβ = (α1, α2, α3, α4, α5, α6, α7, β1, β2, β3, β4, β5, β6, β7) + MPLM75oopCache(u, u, u, u, u, u, P, P, P, P, P, P, + d, d, d, d, d, d, αβ, 1, alg.small_constant_function(uEltypeNoUnits)) +end +#### MPLM106 ########################################################################### +struct MPLM106{F, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +alg_order(alg::MPLM106) = 6 +alg_extrapolates(alg::MPLM106) = true # TODO: Should probably be false + +@cache mutable struct MPLM106oopCache{uType, PType, dType, T, T2} <: + OrdinaryDiffEqConstantCache + uprevprev::uType + uprev3::uType + uprev4::uType + uprev5::uType + uprev6::uType + uprev7::uType + uprev8::uType + uprev9::uType + uprev10::uType + P2::PType + P3::PType + P4::PType + P5::PType + P6::PType + P7::PType + P8::PType + P9::PType + P10::PType + d2::dType + d3::dType + d4::dType + d5::dType + d6::dType + d7::dType + d8::dType + d9::dType + d10::dType + αβ::NTuple{20, T} + step::Int + small_constant::T2 +end + +function MPLM106(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPLM106(linsolve, small_constant_function) +end + +function alg_cache(alg::MPLM106, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, + verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + # TODO: This is currently necessary to get the correct type of P (d is of type rateType) + P, d = evaluate_pds(f, u, p, t) + # TODO: integrator_stats_nf = 1 + + α1 = zero(uEltypeNoUnits) + α2 = zero(uEltypeNoUnits) + α3 = zero(uEltypeNoUnits) + α4 = zero(uEltypeNoUnits) + α5 = zero(uEltypeNoUnits) + α6 = zero(uEltypeNoUnits) + α7 = zero(uEltypeNoUnits) + α8 = zero(uEltypeNoUnits) + α9 = zero(uEltypeNoUnits) + α10 = one(uEltypeNoUnits) + + β1 = 11125 / 4536 * one(uEltypeNoUnits) + β2 = zero(uEltypeNoUnits) + β3 = zero(uEltypeNoUnits) + β4 = 50 / 27 * one(uEltypeNoUnits) + β5 = 85 / 36 * one(uEltypeNoUnits) + β6 = zero(uEltypeNoUnits) + β7 = zero(uEltypeNoUnits) + β8 = 125 / 63 * one(uEltypeNoUnits) + β9 = 25 / 24 * one(uEltypeNoUnits) + β10 = 25 / 81 * one(uEltypeNoUnits) + αβ = (α1, α2, α3, α4, α5, α6, α7, α8, α9, α10, + β1, β2, β3, β4, β5, β6, β7, β8, β9, β10) + MPLM106oopCache(u, u, u, u, u, u, u, u, u, + P, P, P, P, P, P, P, P, P, + d, d, d, d, d, d, d, d, d, + αβ, 1, alg.small_constant_function(uEltypeNoUnits)) +end + +function initialize!(integrator, + cache::Union{MPLM22oopCache, MPLM33oopCache, MPLM43oopCache, + MPLM54oopCache, MPLM75oopCache, MPLM106oopCache, + MPLMMutableCache}) +end + +######################################################################################## +### perform_step! ###################################################################### +######################################################################################## +#### MPLM22 ############################################################################ +@muladd function start_MPLM22_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + dt = dt / 4 + + u = uprev + @inbounds for _ in 1:4 + u = perform_step_MPE_oop(P, d, dt, u, small_constant, linsolve) + + P, d = evaluate_pds(f, u, p, t) + + t += dt + end + + # 4 function evals and 4 solves + nf = 4 + ns = 4 + + return u, t, nf, ns +end + +@muladd function start_MPLM22!(u, P, d, t, dt, uprev, σ, f, p, small_constant, linsolve) + dt = dt / 4 + + u .= uprev + @inbounds for _ in 1:4 + perform_step_MPE!(u, P, d, dt, u, σ, small_constant, linsolve) + + evaluate_pds!(P, d, f, u, p, t) + + t += dt + end + + # 4 function evals and 4 solves + nf = 4 + ns = 4 + + return nf, ns +end + +#TODO Use αβ in MPLM22 +@muladd function perform_step_MPLM22_oop(P, d, dt, uprev, uprevprev, linsolve, + small_constant) + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Main step + σ = add_small_constant(σ, small_constant) + + u = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # statistics: 2 nsolve + + return u +end + +@muladd function perform_step_MPLM22!(u, P, P2, d, d2, dt, uprev, uprevprev, σ, linsolve, + small_constant) + # First σ approximation + @.. broadcast=false σ=uprev + small_constant + + # use lincomb! to handle cases in which d2 is nothing + P2 .= P + lincomb!(d2, 1, d) + + basic_patankar_step!(σ, uprev, P2, d2, σ, dt, linsolve) + + # Main step + @.. broadcast=false σ=σ + small_constant + + P2 .= P + lincomb!(d2, 1, d) + + basic_patankar_step!(u, uprevprev, P2, d2, σ, 2 * dt, linsolve) + + # statistics: 2 nsolve + + return nothing +end + +@muladd function perform_step!(integrator, cache::MPLM22oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, small_constant) = cache + + if integrator.u_modified + cache.step = 1 + end + + if cache.step <= 1 + + # increase step counter + cache.step += 1 + + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM22 + u, _, nf, ns = start_MPLM22_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + + integrator.stats.nf += nf + integrator.stats.nsolve += ns + else + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + u = perform_step_MPLM22_oop(P, d, dt, uprev, uprevprev, alg.linsolve, + small_constant) + integrator.stats.nsolve += 2 + end + + #TODO: Should be possible to use uprev2. But uprev2 is currently not updated. + cache.uprevprev = uprev + + integrator.u = u +end + +@muladd function perform_step!(integrator, cache::MPLM22Cache, repeat_step = false) + (; t, dt, u, uprev, uprev2, f, p) = integrator + (; uprevprev, small_constant, P, P2, D, D2, σ, linsolve) = cache + + if integrator.u_modified + cache.step = 1 + end + + if cache.step <= 1 + + # increase step counter + cache.step += 1 + + evaluate_pds!(P, D, f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM22 + nf, ns = start_MPLM22!(u, P, D, t, dt, uprev, σ, f, p, small_constant, linsolve) + + integrator.stats.nf += nf + integrator.stats.nsolve += ns + else + + # evaluate production matrix + evaluate_pds!(P, D, f, uprev, p, t) + integrator.stats.nf += 1 + + perform_step_MPLM22!(u, P, P2, D, D2, dt, uprev, uprevprev, σ, linsolve, + small_constant) + + integrator.stats.nsolve += 2 + end + + #TODO: Should be possible to use uprev2. But uprev2 is currently not updated. + uprevprev .= uprev +end + +#### MPLM33 ############################################################################ +@muladd function start_MPLM33_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + + # 1 macro step consists of 4 substeps + dts = dt / 4 + + ### first macro step ############################################################### + # substep 1 + u, t, nf, ns = start_MPLM22_oop(P, d, t, dts, uprev, f, p, small_constant, linsolve) + + # substeps 2 - 4 + for _ in 1:3 + uprevprev = uprev + uprev = u + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = perform_step_MPLM22_oop(P, d, dts, uprev, uprevprev, linsolve, small_constant) + t += dts + ns += 2 + end + + v = u + + ### second macro step ############################################################ + for _ in 1:4 + uprevprev = uprev + uprev = u + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = perform_step_MPLM22_oop(P, d, dts, uprev, uprevprev, linsolve, small_constant) + t += dts + ns += 2 + end + + return (v, u), t, nf, ns +end +@muladd function perform_step_MPLM33_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ, + small_constant) + P, P2, P3 = P_tup + d, d2, d3 = d_tup + uprev, uprevprev, uprev3 = u_tup + α1, α2, α3, β1, β2, β3 = αβ + + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Second σ approximation + σ = add_small_constant(σ, small_constant) + + σ = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # Main step + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3) + v = α1 * uprev + α2 * uprevprev + α3 * uprev3 + u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # statistics: 3 nsolve + + return u +end + +@muladd function perform_step!(integrator, cache::MPLM33oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, uprev3, P2, P3, d2, d3, αβ, small_constant) = cache + + #TODO: is this necessary? + if integrator.u_modified + cache.step = 1 + end + + if cache.step == 1 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM33 + v, t, nf, ns = start_MPLM33_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + integrator.stats.nf += nf + integrator.stats.nsolve += ns + + # u at time tspan[1] + dt + u = v[1] + + cache.uprevprev = uprev + + # we use uprev3 as temporary storage for the value of u needed in step 2. + cache.uprev3 = v[2] + + elseif cache.step == 2 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 2*dt (this was computed in step 1) + u = cache.uprev3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + else + + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + P_tup = (P, P2, P3) + d_tup = (d, d2, d3) + u_tup = (uprev, uprevprev, uprev3) + + u = perform_step_MPLM33_oop(P_tup, d_tup, dt, u_tup, alg.linsolve, αβ, + small_constant) + integrator.stats.nsolve += 3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + end + + integrator.u = u + + cache.P3 = cache.P2 + cache.P2 = P + cache.d3 = cache.d2 + cache.d2 = d +end + +#### MPLM43 ############################################################################ +@muladd function start_MPLM43_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + αβ33 = (0, 0, 1, 9 / 4, 0, 3 / 4) + + # 1 macro step consists of 4 substeps + dts = dt / 4 + + ### first macro step ############################################################### + # substeps 1 - 2 + v, t, nf, ns = start_MPLM33_oop(P, d, t, dts, uprev, f, p, small_constant, linsolve) + + uprevprev = uprev + P2 = P + d2 = d + + uprev = v[1] + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = v[2] + + # substeps 3-4 + for _ in 1:2 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P3 = P2 + P2 = P + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3) + d_tup = (d, d2, d3) + u_tup = (uprev, uprevprev, uprev3) + + u = perform_step_MPLM33_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ33, + small_constant) + t += dts + ns += 3 + end + + v1 = u + + ### second macro step ############################################################ + for _ in 1:4 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P3 = P2 + P2 = P + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3) + d_tup = (d, d2, d3) + u_tup = (uprev, uprevprev, uprev3) + + u = perform_step_MPLM33_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ33, + small_constant) + t += dts + ns += 3 + end + + v2 = u + + ### third macro step ############################################################ + for _ in 1:4 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P3 = P2 + P2 = P + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3) + d_tup = (d, d2, d3) + u_tup = (uprev, uprevprev, uprev3) + + u = perform_step_MPLM33_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ33, + small_constant) + t += dts + ns += 3 + end + + return (v1, v2, u), t, nf, ns +end +@muladd function perform_step_MPLM43_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ, + small_constant) + P, P2, P3, P4 = P_tup + d, d2, d3, d4 = d_tup + uprev, uprevprev, uprev3, uprev4 = u_tup + α1, α2, α3, α4, β1, β2, β3, β4 = αβ + + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Second σ approximation + σ = add_small_constant(σ, small_constant) + + σ = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # Main step + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3, β4, P4, d4) + v = α1 * uprev + α2 * uprevprev + α3 * uprev3 + α4 * uprev4 + u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # statistics: 3 nsolve + + return u +end + +@muladd function perform_step!(integrator, cache::MPLM43oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, uprev3, uprev4, P2, P3, P4, d2, d3, d4, αβ, small_constant) = cache + + #TODO: is this necessary? + if integrator.u_modified + cache.step = 1 + end + + if cache.step == 1 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM43 + v, _, nf, ns = start_MPLM43_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + integrator.stats.nf += nf + integrator.stats.nsolve += ns + + # u at time tspan[1] + dt + u = v[1] + + cache.uprevprev = uprev + + # we use uprev3 as temporary storage for the value of u needed in step 2. + cache.uprev3 = v[2] + # we use uprev4 as temporary storage for the value of u needed in step 3. + cache.uprev4 = v[3] + + elseif cache.step == 2 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 2*dt (this was computed in step 1) + u = cache.uprev3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 3 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 3*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 3*dt (this was computed in step 1) + u = cache.uprev4 + + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + else + + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + P_tup = (P, P2, P3, P4) + d_tup = (d, d2, d3, d4) + u_tup = (uprev, uprevprev, uprev3, uprev4) + + u = perform_step_MPLM43_oop(P_tup, d_tup, dt, u_tup, alg.linsolve, αβ, + small_constant) + integrator.stats.nsolve += 3 + + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + end + + integrator.u = u + + cache.P4 = P3 + cache.P3 = P2 + cache.P2 = P + cache.d4 = d3 + cache.d3 = d2 + cache.d2 = d +end + +#### MPLM54 ############################################################################ +#TODO Check nf and ns everywhere! +@muladd function start_MPLM54_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + αβ43 = (1 / 4, 0, 3 / 4, 0, 35 / 18, 1 / 3, 0, 2 / 9) + + # 1 macro step consists of 4 substeps + dts = dt / 4 + + ### first macro step ############################################################### + # substep 1 - 3 + v, t, nf, ns = start_MPLM43_oop(P, d, t, dts, uprev, f, p, small_constant, linsolve) + + uprev3 = uprev + P3 = P + d3 = d + + uprevprev = v[1] + P2, d2 = evaluate_pds(f, uprevprev, p, t) + nf += 1 + + uprev = v[2] + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = v[3] + + # substep 4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P4 = P3 + P3 = P2 + P2 = P + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4) + d_tup = (d, d2, d3, d4) + u_tup = (uprev, uprevprev, uprev3, uprev4) + + u = perform_step_MPLM43_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ43, + small_constant) + t += dts + ns += 3 + + v1 = u + + ### second macro step ############################################################ + for _ in 1:4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P4 = P3 + P3 = P2 + P2 = P + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4) + d_tup = (d, d2, d3, d4) + u_tup = (uprev, uprevprev, uprev3, uprev4) + + u = perform_step_MPLM43_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ43, + small_constant) + t += dts + ns += 3 + end + + v2 = u + + ### third macro step ############################################################ + for _ in 1:4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P4 = P3 + P3 = P2 + P2 = P + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4) + d_tup = (d, d2, d3, d4) + u_tup = (uprev, uprevprev, uprev3, uprev4) + + u = perform_step_MPLM43_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ43, + small_constant) + t += dts + ns += 3 + end + + v3 = u + + # fourth macro step + for _ in 1:4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P4 = P3 + P3 = P2 + P2 = P + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4) + d_tup = (d, d2, d3, d4) + u_tup = (uprev, uprevprev, uprev3, uprev4) + + u = perform_step_MPLM43_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ43, + small_constant) + t += dts + ns += 3 + end + + return (v1, v2, v3, u), t, nf, ns +end +@muladd function perform_step_MPLM54_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ, + small_constant) + P, P2, P3, P4, P5 = P_tup + d, d2, d3, d4, d5 = d_tup + uprev, uprevprev, uprev3, uprev4, uprev5 = u_tup + α1, α2, α3, α4, α5, β1, β2, β3, β4, β5 = αβ + + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Second σ approximation + σ = add_small_constant(σ, small_constant) + + σ = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # Third σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(35 / 18, P, d, 1 / 3, P2, d2, 2 / 9, P4, d4) + v = 1 / 4 * uprev + 3 / 4 * uprev3 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Main step + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3, β4, P4, d4, β5, P5, d5) + v = α1 * uprev + α2 * uprevprev + α3 * uprev3 + α4 * uprev4 + α5 * uprev5 + u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # statistics: 4 nsolve + + return u +end + +@muladd function perform_step!(integrator, cache::MPLM54oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, uprev3, uprev4, uprev5, P2, P3, P4, P5, d2, d3, d4, d5, αβ, small_constant) = cache + + #TODO: is this necessary? + if integrator.u_modified + cache.step = 1 + end + + if cache.step == 1 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM54 + v, _, nf, ns = start_MPLM54_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + integrator.stats.nf += nf + integrator.stats.nsolve += ns + + # u at time tspan[1] + dt + u = v[1] + + cache.uprevprev = uprev + + # we use uprev3 as temporary storage for the value of u needed in step 2. + cache.uprev3 = v[2] + # we use uprev4 as temporary storage for the value of u needed in step 3. + cache.uprev4 = v[3] + # we use uprev5 as temporary storage for the value of u needed in step 4. + cache.uprev5 = v[4] + + elseif cache.step == 2 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 2*dt (this was computed in step 1) + u = cache.uprev3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 3 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 3*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 3*dt (this was computed in step 1) + u = cache.uprev4 + + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 4 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 4*dt (this was computed in step 1) + u = cache.uprev5 + + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + else + # increase step count + cache.step += 1 + + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dt, u_tup, alg.linsolve, αβ, + small_constant) + integrator.stats.nsolve += 4 + + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + end + + integrator.u = u + + cache.P5 = P4 + cache.P4 = P3 + cache.P3 = P2 + cache.P2 = P + cache.d5 = d4 + cache.d4 = d3 + cache.d3 = d2 + cache.d2 = d +end +#### MPLM75 ############################################################################ +@muladd function start_MPLM75_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + αβ54 = (0, 0, 0, 0, 1, 225 / 96, 0, 50 / 96, 200 / 96, 5 / 96) + + # 1 macro steps consists of 4 substeps + dts = dt / 4 + + ### first macro step ############################################################### + # substeps 1 - 4 + v, t, nf, ns = start_MPLM54_oop(P, d, t, dts, uprev, f, p, small_constant, linsolve) + + uprev4 = uprev + P4 = P + d4 = d + + uprev3 = v[1] + P3, d3 = evaluate_pds(f, uprev3, p, t) + nf += 1 + + uprevprev = v[2] + P2, d2 = evaluate_pds(f, uprevprev, p, t) + nf += 1 + + uprev = v[3] + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = v[4] + + v1 = u + + ### second macro step ############################################################ + for _ in 1:4 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ54, + small_constant) + t += dts + ns += 4 + end + + v2 = u + + ### third macro step ############################################################ + for _ in 1:4 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ54, + small_constant) + t += dts + ns += 4 + end + + v3 = u + + ### fourth macro step ############################################################ + for _ in 1:4 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ54, + small_constant) + t += dts + ns += 4 + end + + v4 = u + + ### fifth macro step ############################################################ + for _ in 1:4 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ54, + small_constant) + t += dts + ns += 4 + end + + v5 = u + + ### sixth macro step ############################################################ + for _ in 1:4 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5) + d_tup = (d, d2, d3, d4, d5) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5) + + u = perform_step_MPLM54_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ54, + small_constant) + t += dts + ns += 4 + end + + return (v1, v2, v3, v4, v5, u), t, nf, ns +end + +@muladd function perform_step_MPLM75_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ, + small_constant) + P, P2, P3, P4, P5, P6, P7 = P_tup + d, d2, d3, d4, d5, d6, d7 = d_tup + uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7 = u_tup + α1, α2, α3, α4, α5, α6, α7, β1, β2, β3, β4, β5, β6, β7 = αβ + + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Second σ approximation + σ = add_small_constant(σ, small_constant) + + σ = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # Third σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(35 / 18, P, d, 1 / 3, P2, d2, 2 / 9, P4, d4) + v = 1 / 4 * uprev + 3 / 4 * uprev3 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Fourth σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(225 / 96, P, d, 50 / 96, P3, d3, 200 / 96, P4, d4, 5 / 96, P5, d5) + v = uprev5 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Main step + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3, β4, P4, d4, β5, P5, d5, β6, P6, + d6, β7, P7, d7) + v = α1 * uprev + α2 * uprevprev + α3 * uprev3 + α4 * uprev4 + α5 * uprev5 + + α6 * uprev6 + α7 * uprev7 + u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # statistics: 5 nsolve + + return u +end + +@muladd function perform_step!(integrator, cache::MPLM75oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7, P2, P3, P4, P5, P6, P7, d2, d3, d4, d5, d6, d7, αβ, small_constant) = cache + + #TODO: is this necessary? + if integrator.u_modified + cache.step = 1 + end + + if cache.step == 1 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM75 + v, _, nf, ns = start_MPLM75_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + integrator.stats.nf += nf + integrator.stats.nsolve += ns + + # u at time tspan[1] + dt + u = v[1] + + cache.uprevprev = uprev + + # we use uprev3 as temporary storage for the value of u needed in step 2. + cache.uprev3 = v[2] + # we use uprev4 as temporary storage for the value of u needed in step 3. + cache.uprev4 = v[3] + # we use uprev5 as temporary storage for the value of u needed in step 4. + cache.uprev5 = v[4] + # we use uprev6 as temporary storage for the value of u needed in step 5. + cache.uprev6 = v[5] + # we use uprev7 as temporary storage for the value of u needed in step 6. + cache.uprev7 = v[6] + + elseif cache.step == 2 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 2*dt (this was computed in step 1) + u = cache.uprev3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 3 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 3*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 3*dt (this was computed in step 1) + u = cache.uprev4 + + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 4 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 4*dt (this was computed in step 1) + u = cache.uprev5 + + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 5 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 5*dt (this was computed in step 1) + u = cache.uprev6 + + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 6 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 6*dt (this was computed in step 1) + u = cache.uprev7 + + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + else + # increase step count + cache.step += 1 + + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dt, u_tup, alg.linsolve, αβ, + small_constant) + integrator.stats.nsolve += 5 + + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + end + + integrator.u = u + + cache.P7 = P6 + cache.P6 = P5 + cache.P5 = P4 + cache.P4 = P3 + cache.P3 = P2 + cache.P2 = P + cache.d7 = d6 + cache.d6 = d5 + cache.d5 = d4 + cache.d4 = d3 + cache.d3 = d2 + cache.d2 = d +end + +#### MPLM106 ############################################################################ +@muladd function start_MPLM106_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve) + αβ75 = (0, 0, 0, 0, 0, 0, 1, 12 / 5, 0, 197 / 720, 701 / 360, 43 / 30, 107 / 360, + 467 / 720) + + # 1 macro step consists of 4 substeps + dts = dt / 4 + + ### 1.5 macro steps ############################################################### + v, t, nf, ns = start_MPLM75_oop(P, d, t, dts, uprev, f, p, small_constant, linsolve) + + uprev6 = uprev + P6 = P + d6 = d + + uprev5 = v[1] + P5, d5 = evaluate_pds(f, uprev5, p, t) + nf += 1 + + uprev4 = v[2] + P4, d4 = evaluate_pds(f, uprev4, p, t) + nf += 1 + + uprev3 = v[3] + P3, d3 = evaluate_pds(f, uprev3, p, t) + nf += 1 + + uprevprev = v[4] + P2, d2 = evaluate_pds(f, uprevprev, p, t) + nf += 1 + + uprev = v[5] + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + u = v[6] + + v1 = v[4] + + ### second macro step ############################################################ + # substeps 3 - 4 + for _ in 1:2 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v2 = u + + ### third macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v3 = u + + ### fourth macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v4 = u + + ### fifth macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v5 = u + + ### sixth macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v6 = u + + ### seventh macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v7 = u + + ### eighth macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + v8 = u + + ### ninth macro step ############################################################ + for _ in 1:4 + uprev7 = uprev6 + uprev6 = uprev5 + uprev5 = uprev4 + uprev4 = uprev3 + uprev3 = uprevprev + uprevprev = uprev + uprev = u + + P7 = P6 + P6 = P5 + P5 = P4 + P4 = P3 + P3 = P2 + P2 = P + d7 = d6 + d6 = d5 + d5 = d4 + d4 = d3 + d3 = d2 + d2 = d + + P, d = evaluate_pds(f, uprev, p, t) + nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7) + d_tup = (d, d2, d3, d4, d5, d6, d7) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7) + + u = perform_step_MPLM75_oop(P_tup, d_tup, dts, u_tup, linsolve, αβ75, + small_constant) + t += dts + ns += 5 + end + + return (v1, v2, v3, v4, v5, v6, v7, v8, u), t, nf, ns +end +@muladd function perform_step_MPLM106_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ, + small_constant) + P, P2, P3, P4, P5, P6, P7, P8, P9, P10 = P_tup + d, d2, d3, d4, d5, d6, d7, d8, d9, d10 = d_tup + uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7, uprev8, uprev9, uprev10 = u_tup + α1, α2, α3, α4, α5, α6, α7, α8, α9, α10, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10 = αβ + + # First σ approximation + σ = add_small_constant(uprev, small_constant) + + σ = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + # Second σ approximation + σ = add_small_constant(σ, small_constant) + + σ = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d) + + # Third σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(35 / 18, P, d, 1 / 3, P2, d2, 2 / 9, P4, d4) + v = 1 / 4 * uprev + 3 / 4 * uprev3 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Fourth σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(225 / 96, P, d, 50 / 96, P3, d3, 200 / 96, P4, d4, 5 / 96, P5, d5) + v = uprev5 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Fifth σ approximation + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(12 / 5, P, d, 197 / 720, P3, d3, 701 / 360, P4, d4, 43 / 30, P5, + d5, 107 / 360, P6, d6, 467 / 720, P7, d7) + v = uprev7 + σ = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # Main step + σ = add_small_constant(σ, small_constant) + + Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3, β4, P4, d4, β5, P5, d5, β6, P6, + d6, β7, P7, d7, β8, P8, d8, β9, P9, d9, β10, P10, d10) + v = α1 * uprev + α2 * uprevprev + α3 * uprev3 + α4 * uprev4 + α5 * uprev5 + + α6 * uprev6 + α7 * uprev7 + α8 * uprev8 + α9 * uprev9 + α10 * uprev10 + u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp) + + # statistics: 6 nsolve + + return u +end + +@muladd function perform_step!(integrator, cache::MPLM106oopCache, repeat_step = false) + (; alg, t, dt, uprev, uprev2, f, p) = integrator + (; uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7, uprev8, uprev9, uprev10, + P2, P3, P4, P5, P6, P7, P8, P9, P10, d2, d3, d4, d5, d6, d7, d8, d9, d10, + αβ, small_constant) = cache + + #TODO: is this necessary? + if integrator.u_modified + cache.step = 1 + end + + if cache.step == 1 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # compute initial values for MPLM106 + v, _, nf, ns = start_MPLM106_oop(P, d, t, dt, uprev, f, p, small_constant, + alg.linsolve) + integrator.stats.nf += nf + integrator.stats.nsolve += ns + + # u at time tspan[1] + dt + u = v[1] + + cache.uprevprev = uprev + + # we use uprev3 as temporary storage for the value of u needed in step 2. + cache.uprev3 = v[2] + # we use uprev4 as temporary storage for the value of u needed in step 3. + cache.uprev4 = v[3] + # we use uprev5 as temporary storage for the value of u needed in step 4. + cache.uprev5 = v[4] + # we use uprev6 as temporary storage for the value of u needed in step 5. + cache.uprev6 = v[5] + # we use uprev7 as temporary storage for the value of u needed in step 6. + cache.uprev7 = v[6] + # we use uprev8 as temporary storage for the value of u needed in step 7. + cache.uprev8 = v[7] + # we use uprev9 as temporary storage for the value of u needed in step 8. + cache.uprev9 = v[8] + # we use uprev10 as temporary storage for the value of u needed in step 9. + cache.uprev10 = v[9] + + elseif cache.step == 2 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 2*dt (this was computed in step 1) + u = cache.uprev3 + + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 3 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 3*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 3*dt (this was computed in step 1) + u = cache.uprev4 + + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 4 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 4*dt (this was computed in step 1) + u = cache.uprev5 + + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 5 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 5*dt (this was computed in step 1) + u = cache.uprev6 + + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 6 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 6*dt (this was computed in step 1) + u = cache.uprev7 + + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 7 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 7*dt (this was computed in step 1) + u = cache.uprev8 + + cache.uprev8 = uprev7 + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 8 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 8*dt (this was computed in step 1) + u = cache.uprev9 + + cache.uprev9 = uprev8 + cache.uprev8 = uprev7 + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + elseif cache.step == 9 + # increase step count + cache.step += 1 + + # evaluate production matrix at tspan[1] + 4*dt + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + # u at time tspan[1] + 9*dt (this was computed in step 1) + u = cache.uprev10 + + cache.uprev10 = uprev9 + cache.uprev9 = uprev8 + cache.uprev8 = uprev7 + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + else + + # evaluate production matrix + P, d = evaluate_pds(f, uprev, p, t) + integrator.stats.nf += 1 + + P_tup = (P, P2, P3, P4, P5, P6, P7, P8, P9, P10) + d_tup = (d, d2, d3, d4, d5, d6, d7, d8, d9, d10) + u_tup = (uprev, uprevprev, uprev3, uprev4, uprev5, uprev6, uprev7, uprev8, uprev9, + uprev10) + + u = perform_step_MPLM106_oop(P_tup, d_tup, dt, u_tup, alg.linsolve, αβ, + small_constant) + integrator.stats.nsolve += 6 + + cache.uprev10 = uprev9 + cache.uprev9 = uprev8 + cache.uprev8 = uprev7 + cache.uprev7 = uprev6 + cache.uprev6 = uprev5 + cache.uprev5 = uprev4 + cache.uprev4 = uprev3 + cache.uprev3 = uprevprev + cache.uprevprev = uprev + end + + integrator.u = u + + cache.P10 = P9 + cache.P9 = P8 + cache.P8 = P7 + cache.P7 = P6 + cache.P6 = P5 + cache.P5 = P4 + cache.P4 = P3 + cache.P3 = P2 + cache.P2 = P + cache.d10 = d9 + cache.d9 = d8 + cache.d8 = d7 + cache.d7 = d6 + cache.d6 = d5 + cache.d5 = d4 + cache.d4 = d3 + cache.d3 = d2 + cache.d2 = d +end diff --git a/src/mprk.jl b/src/mprk.jl index 61b73a068..9aaa285b0 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -483,18 +483,42 @@ end function initialize!(integrator, cache::MPEConstantCache) end +@muladd function perform_step_MPE_oop(P, d, dt, uprev, small_constant, linsolve) + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(uprev, small_constant) + + u = basic_patankar_step(uprev, P, σ, dt, linsolve, d) + + return u +end + +@muladd function perform_step_MPE!(u, P, d, dt, uprev, σ, small_constant, linsolve) + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + basic_patankar_step!(u, uprev, P, d, σ, dt, linsolve) + + return nothing +end + +@muladd function perform_step_MPE_conservative!(u, P, dt, uprev, σ, small_constant, + linsolve) + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=σ + small_constant + + basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve) + + return nothing +end + @muladd function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) (; alg, t, dt, uprev, f, p) = integrator (; small_constant) = cache - # evaluate production matrix and destruction vector P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 - # avoid division by zero due to zero Patankar weights - σ = add_small_constant(uprev, small_constant) - - u = basic_patankar_step(uprev, P, σ, dt, alg.linsolve, d) + u = perform_step_MPE_oop(P, d, dt, uprev, small_constant, alg.linsolve) integrator.stats.nsolve += 1 integrator.u = u diff --git a/test/runtests.jl b/test/runtests.jl index faf5862c8..8487543fe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2350,7 +2350,8 @@ end probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_bertolazzi, prob_pds_brusselator, prob_pds_npzd, - prob_pds_sir, prob_pds_stratreac) + prob_pds_sir, prob_pds_stratreac, + prob_pds_saceirqd, prob_pds_diffusion) @testset "$alg" for alg in algs @testset "$i" for (i, prob) in enumerate(probs) if prob == prob_pds_stratreac && alg == SSPMPRK22(0.5, 1.0) @@ -2378,7 +2379,8 @@ end #prob_pds_robertson not included probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_bertolazzi, prob_pds_brusselator, - prob_pds_npzd, prob_pds_sir, prob_pds_stratreac) + prob_pds_npzd, prob_pds_sir, prob_pds_stratreac, + prob_pds_saceirqd, prob_pds_diffusion) @testset "$alg" for alg in algs @testset "$prob" for prob in probs tspan = prob.tspan