Skip to content

Commit

Permalink
Examples for PDS and MPRK updated (need to learn how to use VSCode ...)
Browse files Browse the repository at this point in the history
  • Loading branch information
@SKopecz
SKopecz committed Oct 31, 2023
1 parent b68c91a commit 65ad63e
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 44 deletions.
16 changes: 9 additions & 7 deletions examples/01_example_proddest.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,12 +89,12 @@ linmod_PDS_ip = ProdDestODEProblem(linmodP!, linmodD!, u0, tspan)
linmod_ConsPDS_ip = ConsProdDestODEProblem(linmodP!, u0, tspan)

# solutions
sol_linmod_f_op = solve(linmod_f_op, Tsit5())
sol_linmod_f_ip = solve(linmod_f_ip,Tsit5())
sol_linmod_PDS_op = solve(linmod_PDS_op, Tsit5())
sol_linmod_PDS_ip = solve(linmod_PDS_ip,Tsit5())
sol_linmod_ConsPDS_op = solve(linmod_ConsPDS_op, Tsit5())
sol_linmod_ConsPDS_ip = solve(linmod_ConsPDS_ip,Tsit5())
sol_linmod_f_op = solve(linmod_f_op, Tsit5());
sol_linmod_f_ip = solve(linmod_f_ip,Tsit5());
sol_linmod_PDS_op = solve(linmod_PDS_op, Tsit5());
sol_linmod_PDS_ip = solve(linmod_PDS_ip,Tsit5());
sol_linmod_ConsPDS_op = solve(linmod_ConsPDS_op, Tsit5());
sol_linmod_ConsPDS_ip = solve(linmod_ConsPDS_ip,Tsit5());

# check equality of solutions
@assert sol_linmod_f_op.t ≈ sol_linmod_f_ip.t ≈ sol_linmod_PDS_op.t ≈ sol_linmod_PDS_ip.t ≈ sol_linmod_ConsPDS_op.t ≈ sol_linmod_ConsPDS_ip.t
Expand All @@ -103,7 +103,9 @@ sol_linmod_ConsPDS_ip = solve(linmod_ConsPDS_ip,Tsit5())
# check that we really do not use too many additional allocations for in-place implementations
alloc1 = @allocated(solve(linmod_f_ip, Tsit5()))
alloc2 = @allocated(solve(linmod_PDS_ip, Tsit5()))
alloc3 = @allocated(solve(linmod_ConsPDS_ip, Tsit5()))
@assert 0.95 < alloc1/alloc2 < 1.05
@assert 0.95 < alloc1/alloc3 < 1.05

##########################################################################################################################
### Example 2: Lotka-Volterra ############################################################################################
Expand Down Expand Up @@ -222,7 +224,7 @@ function fdupwindD!(D,u,p,t)
end
# problem with dense matrices
fdupwind_PDS_dense = ProdDestODEProblem(fdupwindP!, fdupwindD!, u0, tspan);
# proboem with sparse matrices
# problem with sparse matrices
p_prototype = spdiagm(-1 => ones(eltype(u0),N-1), N-1 => ones(eltype(u0),1))
d_prototype = zero(u0)
PD_sparse = ProdDestFunction(fdupwindP!,fdupwindD!;p_prototype = p_prototype, d_prototype = d_prototype);
Expand Down
56 changes: 19 additions & 37 deletions examples/03_example_mprk22.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,9 +40,7 @@ end

# out-of-place implementation
linmodP(u,p,t) = [0.0 p[2]u[2]; p[1]*u[1] 0.0]
linmodD(u,p,t) = [0.0; 0.0]
linmodPD = ProdDestFunction(linmodP,linmodD; analytic=f_analytic)
linmod_op = ProdDestODEProblem(linmodPD, u0, tspan, p)
linmod_op = ConsProdDestODEProblem(linmodP, u0, tspan, p; analytic=f_analytic)

# solutions with constant equispaced time steps
dt0 = 0.25
Expand All @@ -69,19 +67,13 @@ function linmodP!(P,u,p,t)
P[2, 1] = 5.0*u[1]
return nothing
end
function linmodD!(D,u,p,t)
D .= 0
return nothing
end
PD_ip = ProdDestFunction(linmodP!,linmodD!,p_prototype=zeros(2,2), d_prototype=zeros(2,1),analytic=f_analytic)
#BUG: linmod_ip cannot be solved if prototypes are not given.
linmod_ip = ProdDestODEProblem(PD_ip, u0, tspan, p)
linmod_ip = ConsProdDestODEProblem(linmodP!, u0, tspan, p; analytic=f_analytic)

# solutions with constant equispaced time steps
dt0 = 0.25
sol_linmod_ip_MPE = solve(linmod_ip, MPE(), dt=dt0)
sol_linmod_ip_MPRK22_1 = solve(linmod_ip, MPRK22(1.0), dt=dt0, adaptive=false)
sol_linmod_ip_MPRK22_½ = solve(linmod_ip, MPRK22(0.5), dt=dt0, adaptive=false)
sol_linmod_ip_MPE = solve(linmod_ip, MPE(), dt=dt0);
sol_linmod_ip_MPRK22_1 = solve(linmod_ip, MPRK22(1.0), dt=dt0, adaptive=false);
sol_linmod_ip_MPRK22_½ = solve(linmod_ip, MPRK22(0.5), dt=dt0, adaptive=false);

# plots
using Plots
Expand All @@ -104,14 +96,13 @@ tspan = (0.0, 30.0)

# out-of-place implementation
nonlinmodP(u,p,t) = [0.0 0.0 0.0; u[2]*u[1]/(u[1]+1.0) 0.0 0.0; 0.0 0.3*u[2] 0.0]
nonlinmodD(u,p,t) = [0.0; 0.0; 0.0]
nonlinmod_op = ProdDestODEProblem(nonlinmodP, nonlinmodD, u0, tspan, p)
nonlinmod_op = ConsProdDestODEProblem(nonlinmodP, u0, tspan, p)

# solutions with constant equispaced time steps
dt0 = 1.0
sol_nonlinmod_op_Tsit5 = solve(nonlinmod_op, Tsit5(), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_op_MPRK22_1 = solve(nonlinmod_op, MPRK22(1.0), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_op_MPRK22_½ = solve(nonlinmod_op, MPRK22(0.5), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_op_Tsit5 = solve(nonlinmod_op, Tsit5(), dt=dt0, abstol=1e-2, reltol=1e-3);
sol_nonlinmod_op_MPRK22_1 = solve(nonlinmod_op, MPRK22(1.0), dt=dt0, abstol=1e-2, reltol=1e-3);
sol_nonlinmod_op_MPRK22_½ = solve(nonlinmod_op, MPRK22(0.5), dt=dt0, abstol=1e-2, reltol=1e-3);

# plots
using Plots
Expand All @@ -127,17 +118,13 @@ function nonlinmodP!(P,u,p,t)
P[3, 2] = 0.3*u[2]
return nothing
end
function nonlinmodD!(D,u,p,t)
D .= 0
return nothing
end
nonlinmod_ip = ProdDestODEProblem(nonlinmodP!, nonlinmodD!, u0, tspan, p)
nonlinmod_ip = ConsProdDestODEProblem(nonlinmodP!, u0, tspan, p)

# solutions with constant equispaced time steps
dt0 = 1.0
sol_nonlinmod_ip_Tsit5 = solve(nonlinmod_ip, Tsit5(), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_ip_MPRK22_1 = solve(nonlinmod_ip, MPRK22(1.0), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_ip_MPRK22_½ = solve(nonlinmod_ip, MPRK22(0.5), dt=dt0, abstol=1e-2, reltol=1e-3)
sol_nonlinmod_ip_Tsit5 = solve(nonlinmod_ip, Tsit5(), dt=dt0, abstol=1e-2, reltol=1e-3);
sol_nonlinmod_ip_MPRK22_1 = solve(nonlinmod_ip, MPRK22(1.0), dt=dt0, abstol=1e-2, reltol=1e-3);
sol_nonlinmod_ip_MPRK22_½ = solve(nonlinmod_ip, MPRK22(0.5), dt=dt0, abstol=1e-2, reltol=1e-3);

# plots
using Plots
Expand All @@ -157,12 +144,11 @@ p = [1e4,4e-2, 3e7]

# out-of-place
robertsonP(u,p,t) = [0.0 1e4*u[2]*u[3] 0.0; 4e-2*u[1] 0.0 0.0; 0.0 3e7*u[2]^2 0.0]
robertsonD(u,p,t) = [0.0; 0.0; 0.0]
robertson_op = ProdDestODEProblem(robertsonP, robertsonD, u0, tspan, p)
robertson_op = ConsProdDestODEProblem(robertsonP, u0, tspan, p)

# Test constant time step size
sol_robertson_op_RB23 = solve(robertson_op, Rosenbrock23(), abstol=1e-3, reltol=1e-2)
sol_robertson_op_MPRK22_1 = solve(robertson_op, MPRK22(1.1), abstol=1e-3, reltol=1e-2)
sol_robertson_op_RB23 = solve(robertson_op, Rosenbrock23(), abstol=1e-3, reltol=1e-2);
sol_robertson_op_MPRK22_1 = solve(robertson_op, MPRK22(1.1), abstol=1e-3, reltol=1e-2);

# plot solution
p1 = plot(sol_robertson_op_RB23[2:end],idxs = [(0, 1), ((x,y)-> (x,1e4.*y) , 0, 2), (0, 3)], color=palette(:default)[1:3]',legend=:right, xaxis=:log)
Expand All @@ -181,15 +167,11 @@ function robertsonP!(P,u,p,t)
P[3, 2] = 3e7*u[2]^2
return nothing
end
function robertsonD!(D,u,p,t)
D .= 0
return nothing
end
robertson_ip = ProdDestODEProblem(robertsonP!, robertsonD!, u0, tspan, p)
robertson_ip = ConsProdDestODEProblem(robertsonP!, u0, tspan, p)

# Test constant time step size
sol_robertson_ip_RB23 = solve(robertson_ip, Rosenbrock23(autodiff=false), abstol=1e-3, reltol=1e-2)
sol_robertson_ip_MPRK22_1 = solve(robertson_ip, MPRK22(1.1), abstol=1e-3, reltol=1e-2)
sol_robertson_ip_RB23 = solve(robertson_ip, Rosenbrock23(autodiff=false), abstol=1e-3, reltol=1e-2);
sol_robertson_ip_MPRK22_1 = solve(robertson_ip, MPRK22(1.1), abstol=1e-3, reltol=1e-2);

# plot solution
p1 = plot(sol_robertson_op_RB23[2:end],idxs = [(0, 1), ((x,y)-> (x,1e4.*y) , 0, 2), (0, 3)], color=palette(:default)[1:3]',legend=:right, xaxis=:log)
Expand Down

0 comments on commit 65ad63e

Please sign in to comment.