formatted code

SKopecz · Jan 29, 2024 · 8ef244f · 8ef244f
1 parent 6930aff
commit 8ef244f
Show file tree

Hide file tree

Showing 5 changed files with 45 additions and 30 deletions.
diff --git a/examples/04_example_problemlibrary.jl b/examples/04_example_problemlibrary.jl
@@ -18,9 +18,9 @@ include("utilities.jl")
 f2(t, x, y) = (t, x + y)
 f3(t, x, y, z) = (t, x + y + z)
 f4(t, u1, u2, u3, u4) = (t, u1 + u2 + u3 + u4)
-f_npzd(t, u1, u2, u3, u4) = (t, 0.66*(u1 + u2 + u3 + u4))
+f_npzd(t, u1, u2, u3, u4) = (t, 0.66 * (u1 + u2 + u3 + u4))
 f6(t, u1, u2, u3, u4, u5, u6) = (t, u1 + u2 + u3 + u4 + u5 + u6)
-f_brusselator(t, u1, u2, u3, u4, u5, u6) = (t, 0.55*(u1 + u2 + u3 + u4 + u5 + u6))
+f_brusselator(t, u1, u2, u3, u4, u5, u6) = (t, 0.55 * (u1 + u2 + u3 + u4 + u5 + u6))
 
 ## linear model ##########################################################
 sol_linmod_MPE = solve(prob_pds_linmod, MPE(), dt = 0.2);
@@ -31,11 +31,13 @@ plot!(sol_linmod_MPE, idxs = (f2, 0, 1, 2))
 
 # convergence order
 # error based on analytic solution
-convergence_tab_plot(prob_pds_linmod, [MPE(), Euler()]; dts = 0.5 .^ (3:18), analytic = true, order_plot = true)
+convergence_tab_plot(prob_pds_linmod, [MPE(), Euler()]; dts = 0.5 .^ (3:18),
+                     analytic = true, order_plot = true)
 savefig("figs/error_linmod_analytic.svg")
 # error based on reference solution
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14);
-convergence_tab_plot(prob_pds_linmod, [MPE(), Euler()], test_setup; dts = 0.5 .^ (3:18), order_plot = true)
+convergence_tab_plot(prob_pds_linmod, [MPE(), Euler()], test_setup; dts = 0.5 .^ (3:18),
+                     order_plot = true)
 savefig("figs/error_linmod_reference.svg")
 
 ## nonlinear model ########################################################
@@ -49,7 +51,8 @@ plot!(sol_nonlinmod_MPE, idxs = (f3, 0, 1, 2, 3))
 
 # convergence order 
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
-convergence_tab_plot(prob_pds_nonlinmod, [MPE(), Euler()], test_setup; dts = 0.5 .^ (3:17), order_plot = true)
+convergence_tab_plot(prob_pds_nonlinmod, [MPE(), Euler()], test_setup; dts = 0.5 .^ (3:17),
+                     order_plot = true)
 
 ## robertson problem ######################################################
 sol_robertson = solve(prob_pds_robertson, Rosenbrock23());
@@ -68,11 +71,13 @@ sol_brusselator_MPE = solve(prob_pds_brusselator, MPE(), dt = 0.25);
 # plot
 plot(sol_brusselator, legend = :outerright)
 myplot!(sol_brusselator_MPE, "MPE")
-plot!(sol_brusselator_MPE, idxs = (f_brusselator, 0, 1, 2, 3, 4, 5, 6), label="f_brusselator")
+plot!(sol_brusselator_MPE, idxs = (f_brusselator, 0, 1, 2, 3, 4, 5, 6),
+      label = "f_brusselator")
 
 # convergence order 
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
-convergence_tab_plot(prob_pds_brusselator, [MPE()], test_setup; dts = 0.5 .^ (3:17), order_plot = true)
+convergence_tab_plot(prob_pds_brusselator, [MPE()], test_setup; dts = 0.5 .^ (3:17),
+                     order_plot = true)
 
 ## SIR model ##############################################################
 sol_sir = solve(prob_pds_sir, Tsit5());
@@ -86,41 +91,43 @@ plot!(sol_sir_MPE, idxs = (f3, 0, 1, 2, 3), label = "f3")
 p2 = plot(sol_sir)
 myplot!(sol_sir_Euler, "Euler")
 plot!(sol_sir_Euler, idxs = (f3, 0, 1, 2, 3), label = "f3")
-plot(p1,p2)
+plot(p1, p2)
 
 # convergence order 
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
-convergence_tab_plot(prob_pds_sir, [MPE(), Euler()], test_setup; dts = 0.5 .^ (1:15), order_plot = true)
+convergence_tab_plot(prob_pds_sir, [MPE(), Euler()], test_setup; dts = 0.5 .^ (1:15),
+                     order_plot = true)
 
 ## bertolazzi problem #####################################################
 sol_bertolazzi = solve(prob_pds_bertolazzi, TRBDF2());
 sol_bertolazzi_MPE = solve(prob_pds_bertolazzi, MPE(), dt = 0.01);
 
 # plot 
 plot(sol_bertolazzi, legend = :right)
-myplot!(sol_bertolazzi_MPE,"MPE")
+myplot!(sol_bertolazzi_MPE, "MPE")
 ylims!((-0.5, 3.5))
 plot!(sol_bertolazzi_MPE, idxs = (f3, 0, 1, 2, 3))
 
 # convergence order 
 test_setup = Dict(:alg => Rosenbrock23(), :reltol => 1e-8, :abstol => 1e-8)
-convergence_tab_plot(prob_pds_bertolazzi, [MPE(), ImplicitEuler()], test_setup; dts = 0.5 .^ (10:15), order_plot = true)
+convergence_tab_plot(prob_pds_bertolazzi, [MPE(), ImplicitEuler()], test_setup;
+                     dts = 0.5 .^ (10:15), order_plot = true)
 
 ### npzd problem ##########################################################
 sol_npzd = solve(prob_pds_npzd, Rosenbrock23());
 sol_npzd_MPE = solve(prob_pds_npzd, MPE(), dt = 0.1);
 
 # plot
 plot(sol_npzd)
-myplot!(sol_npzd_MPE,"MPE")
+myplot!(sol_npzd_MPE, "MPE")
 plot!(sol_npzd_MPE, idxs = (f_npzd, 0, 1, 2, 3, 4), label = "f_npzd")
-plot!(legend=:bottomright)
+plot!(legend = :bottomright)
 
 # convergence order 
 # error should take all time steps into account, not only the final time! 
 test_setup = Dict(:alg => Rosenbrock23(), :reltol => 1e-14, :abstol => 1e-14)
-convergence_tab_plot(prob_pds_npzd, [MPE(), ImplicitEuler()], test_setup; dts = 0.5 .^ (5:17), order_plot = true)
-
+convergence_tab_plot(prob_pds_npzd, [MPE(), ImplicitEuler()], test_setup;
+                     dts = 0.5 .^ (5:17), order_plot = true)
 
 ### strat reac problem ####################################################
 sol_stratreac = solve(prob_pds_stratreac, TRBDF2(autodiff = false));
@@ -162,12 +169,13 @@ p6 = plot(sol_stratreac, idxs = (0, 6), xticks = [tspan[1], tspan[2]], legend =
 #plot!(sol_stratreac_MPE.t, tmp[6, :])
 ylims!((1.08e9, 1.1e9))
 
-p7 = plot(sol_stratreac, idxs = (g1, 0, 1, 2, 3, 4, 5, 6), xticks = [tspan[1], tspan[2]], legend = :outertop)
+p7 = plot(sol_stratreac, idxs = (g1, 0, 1, 2, 3, 4, 5, 6), xticks = [tspan[1], tspan[2]],
+          legend = :outertop)
 #plot!(sol_stratreac_MPE.t,
 #      vec(mapslices((x) -> abs(x[1] + x[2] + 3 * x[3] + 2 * x[4] + x[5] + 2 * x[6] - linear_invariant_1) /
 #                           linear_invariant_1, tmp, dims = 1)))
 p8 = plot(sol_stratreac, idxs = (g2, 0, 1, 2, 3, 4, 5, 6),
-                 xticks = [tspan[1], tspan[2]], legend = :outertop)
+          xticks = [tspan[1], tspan[2]], legend = :outertop)
 #plot!(sol_stratreac_MPE.t,
 #      vec(mapslices((x) -> abs(x[5] + x[6] - linear_invariant_2) / linear_invariant_2, tmp,
 #                    dims = 1)))

diff --git a/examples/utilities.jl b/examples/utilities.jl
@@ -2,7 +2,8 @@ using Plots
 using DiffEqDevTools: test_convergence, analyticless_test_convergence
 using PrettyTables: pretty_table
 
-function convergence_tab_plot(prob, algs, test_setup = nothing; dts = 0.5 .^ (1:10), order_plot = false, analytic = false)
+function convergence_tab_plot(prob, algs, test_setup = nothing; dts = 0.5 .^ (1:10),
+                              order_plot = false, analytic = false)
     sims = Array{Any}(undef, length(algs))
     for i in eachindex(algs)
         #convergence order
@@ -50,7 +51,7 @@ function _myplot(plotf, sol, name = "", analytic = false)
     N = length(sol.u[1])
     if analytic == true
         plotf(sol, color = palette(:default)[1:(2 * N)]', legend = :right,
-             plot_analytic = true)
+              plot_analytic = true)
     else
         plotf(sol, color = palette(:default)[1:N]', legend = :right, plot_analytic = false)
     end

diff --git a/src/PDSProblemLibrary.jl b/src/PDSProblemLibrary.jl
@@ -9,7 +9,8 @@ function f_linmod_analytic(u0, p, t)
     return ((u₁⁰ + u₂⁰) * [b; a] + exp(-c * t) * (a * u₁⁰ - b * u₂⁰) * [1; -1]) / c
 end
 u0_linmod = @SVector [0.9, 0.1]
-prob_pds_linmod = ConservativePDSProblem(P_linmod, u0_linmod, (0.0, 2.0), analytic = f_linmod_analytic)
+prob_pds_linmod = ConservativePDSProblem(P_linmod, u0_linmod, (0.0, 2.0),
+                                         analytic = f_linmod_analytic)
 
 # nonlinear model problem
 function P_nonlinmod(u, p, t)
@@ -120,7 +121,7 @@ function P_stratreac(u, p, t)
 end
 function d_stratreac(u, p, t)
     O1D, O, O3, O2, NO, NO2 = u
-    
+
     k2 = 8.018e-17
     k11 = 1.0e-8
 

diff --git a/src/mprk.jl b/src/mprk.jl
@@ -1,6 +1,7 @@
-add_small_constant(v, small_constant) = v .+ small_constant    
-add_small_constant(v::SVector{N,T}, small_constant::T) where {N,T} = v + SVector{N,T}(ntuple(i -> small_constant, N))
-
+add_small_constant(v, small_constant) = v .+ small_constant
+function add_small_constant(v::SVector{N, T}, small_constant::T) where {N, T}
+    v + SVector{N, T}(ntuple(i -> small_constant, N))
+end
 
 function build_mprk_matrix(P, sigma, dt)
     #=
@@ -148,7 +149,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
              linsolve_tmp, linsolve, weight)
 end
 
-struct MPEConstantCache{T} <: OrdinaryDiffEqConstantCache 
+struct MPEConstantCache{T} <: OrdinaryDiffEqConstantCache
     small_constant::T
 end
 
@@ -191,7 +192,7 @@ function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false)
     sol = solve(linprob, alg.linsolve, Pl = alg.precs[1], Pr = alg.precs[2],
                 alias_A = false, alias_b = false,
                 assumptions = LinearSolve.OperatorAssumptions(true))
-    u = sol.u    
+    u = sol.u
 
     k = f(u, p, t + dt) # For the interpolation, needs k at the updated point
     integrator.stats.nf += 1

diff --git a/src/proddest.jl b/src/proddest.jl
@@ -60,8 +60,10 @@ function Base.getproperty(obj::PDSFunction, sym::Symbol)
         return nothing
     elseif sym === :sparsity
         return nothing
-    elseif sym === :sys 
-        return SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing, nothing, nothing)        
+    elseif sym === :sys
+        return SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing,
+                                                                                nothing,
+                                                                                nothing)
     else # fallback to getfield
         return getfield(obj, sym)
     end
@@ -210,8 +212,10 @@ function Base.getproperty(obj::ConservativePDSFunction, sym::Symbol)
         return nothing
     elseif sym === :sparsity
         return nothing
-    elseif sym === :sys 
-        return SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing, nothing, nothing)
+    elseif sym === :sys
+        return SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing,
+                                                                                nothing,
+                                                                                nothing)
     else # fallback to getfield
         return getfield(obj, sym)
     end