You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: src/mpdec.jl
+18-14Lines changed: 18 additions & 14 deletions
Original file line number
Diff line number
Diff line change
@@ -4,15 +4,15 @@
4
4
A family of arbitrary order modified Patankar-Runge-Kutta algorithms for
5
5
production-destruction systems. Each member of this family is an adaptive, one-step method which is
6
6
Kth order accurate, unconditionally positivity-preserving, and linearly
7
-
implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10.
7
+
implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10.
8
8
Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default.
9
9
These methods support adaptive time stepping, using the numerical solution obtained with one correction step less as a lower-order approximation to estimate the error.
10
10
The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems and
11
11
further investigated in Torlo, Öffner and Ranocha (2022).
12
12
13
13
For nonconservative production–destruction systems we use a straight forward extension
14
14
analogous to [`MPE`](@ref).
15
-
A general discussion of DeC schemes applied to non-autonomous differential equations
15
+
A general discussion of DeC schemes applied to non-autonomous differential equations
16
16
and using general integration nodes is given by Ong and Spiteri (2020).
17
17
18
18
The MPDeC methods require the special structure of a
@@ -51,7 +51,7 @@ end
51
51
52
52
functionsmall_constant_function_MPDeC(type)
53
53
if type == Float64
54
-
# small_constant is chosen such that
54
+
# small_constant is chosen such that
55
55
# the testset "Zero initial values" passes.
56
56
small_constant =1e-300
57
57
else
@@ -159,7 +159,7 @@ function get_constant_parameters(alg::MPDeC, type)
159
159
else
160
160
error("MPDeC requires 2 ≤ K ≤ 10.")
161
161
end
162
-
else# alg.nodes === :gausslobatto
162
+
else# alg.nodes === :gausslobatto
163
163
if alg.M ==1
164
164
nodes = [0, oneType]
165
165
theta = [0 oneType/2; 0 oneType/2]
@@ -247,12 +247,12 @@ function build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, theta,
247
247
small_constant)
248
248
N =length(uprev)
249
249
if f isa PDSFunction
250
-
# Additional destruction terms
250
+
# Additional destruction terms
251
251
Mmat, rhs =_build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
252
252
theta,
253
253
small_constant, f.d)
254
254
else
255
-
# No additional destruction terms
255
+
# No additional destruction terms
256
256
Mmat, rhs =_build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
257
257
theta,
258
258
small_constant)
@@ -479,14 +479,14 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS
479
479
fill!(tmp, zero(eltype(tmp)))
480
480
481
481
if dt_th ≥0
482
-
for j in1:n # run through columns of P
482
+
for j in1:n # run through columns of P
483
483
for idx_P innzrange(P, j) # run through rows of P
484
484
i = P_rows[idx_P]
485
485
dt_th_P = dt_th * P_vals[idx_P]
486
486
if i != j
487
487
for idx_M innzrange(M, j)
488
488
if M_rows[idx_M] == i
489
-
M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij
489
+
M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij
490
490
break
491
491
end
492
492
end
@@ -513,9 +513,9 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS
513
513
end
514
514
end
515
515
else# dt ≤ 0
516
-
for j in1:n # j is column index
516
+
for j in1:n # j is column index
517
517
for idx_P innzrange(P, j)
518
-
i = P_rows[idx_P] # i is row index
518
+
i = P_rows[idx_P] # i is row index
519
519
dt_th_P = dt_th * P_vals[idx_P]
520
520
if i != j
521
521
for idx_M innzrange(M, i)
@@ -635,13 +635,13 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
635
635
P2 =p_prototype(u, f) # stores the linear system matrix
636
636
ifissparse(P2)
637
637
# We need to ensure that evaluating the production function does
638
-
# not alter the sparsity pattern given by the production matrix prototype
638
+
# not alter the sparsity pattern given by the production matrix prototype
639
639
f.p(P2, uprev, p, t)
640
640
@assert P.rowval == P2.rowval&&P.colptr == P2.colptr "Evaluation of the production terms must not alter the sparsity pattern given by the prototype."
641
641
642
642
if alg.K >2
643
643
# Negative weights of MPDeC(K) , K > 2 require
644
-
# a symmetric sparsity pattern of the linear system matrix P2
644
+
# a symmetric sparsity pattern of the linear system matrix P2
645
645
P2 = P + P'
646
646
end
647
647
end
@@ -657,7 +657,9 @@ function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
657
657
# not be altered, since it is needed to compute the adaptive time step
0 commit comments