Skip to content

Commit d2577da

Browse files
ChrisRackauckas-ClaudeChrisRackauckas
ChrisRackauckas-Claude
and
ChrisRackauckas
authored
ExponentialRK: remove per-step allocations from EPIRK/Exp4/EXPRB53s3 steppers (#3928)
The in-place perform_step! implementations of Exp4, EPIRK4s3A/B, EPIRK5s3, EXPRB53s3, and EPIRK5P1/P2 allocated every step: - literal snapshot-time arrays ([dt/3, 2dt/3, dt], [g11, g21, g31], ...) constructed per step for phiv_timestep!; now a ts vector preallocated in the cache and filled in place - coefficient vectors for 3-column linear combinations built per step (mul!(rtmp, K, [-7/300, 97/150, -37/300])); now fused broadcasts over the columns of K - unfused broadcast right-hand sides that materialized N-sized temporaries (B[:, 4] .= (32 / dt^2) * rtmp and friends, and two broadcast-plus-broadcast sums); now fully fused with dots The N-sized temporaries in particular scale with the PDE size. Together with ExponentialUtilities#254/#255 this takes per-step allocations on the Kuramoto-Sivashinsky pseudospectral benchmark (N=128, non-allocating nonlinearity) from 0.6-1.2 MB/step to 11-23 KB/step for the EPIRK family; solutions are identical. Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
1 parent 4b3cafd commit d2577da

3 files changed

Lines changed: 70 additions & 42 deletions

File tree

‎lib/OrdinaryDiffEqExponentialRK/Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "OrdinaryDiffEqExponentialRK"
22
uuid = "e0540318-69ee-4070-8777-9e2de6de23de"
33
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4-
version = "2.1.0"
4+
version = "2.1.1"
55

66
[deps]
77
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"

‎lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl‎

Lines changed: 24 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -443,7 +443,7 @@ for (Alg, Cache) in [
443443
end
444444
end
445445

446-
@cache struct Exp4Cache{uType, rateType, JCType, FType, matType, JType, KsType} <:
446+
@cache struct Exp4Cache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
447447
ExpRKCache
448448
u::uType
449449
uprev::uType
@@ -458,6 +458,7 @@ end
458458
J::JType
459459
B::matType
460460
KsCache::KsType
461+
ts::tsType
461462
end
462463
function alg_cache(
463464
alg::Exp4, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -489,10 +490,11 @@ function alg_cache(
489490
# Allocate caches for phiv_timestep
490491
maxiter = min(alg.m, n)
491492
KsCache = _phiv_timestep_caches(u, maxiter, 1)
492-
return Exp4Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache)
493+
ts = Vector{typeof(dt)}(undef, 3)
494+
return Exp4Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache, ts)
493495
end
494496

495-
@cache struct EPIRK4s3ACache{uType, rateType, JCType, FType, matType, JType, KsType} <:
497+
@cache struct EPIRK4s3ACache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
496498
ExpRKCache
497499
u::uType
498500
uprev::uType
@@ -507,6 +509,7 @@ end
507509
J::JType
508510
B::matType
509511
KsCache::KsType
512+
ts::tsType
510513
end
511514
function alg_cache(
512515
alg::EPIRK4s3A, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -537,10 +540,11 @@ function alg_cache(
537540
# Allocate caches for phiv_timestep
538541
maxiter = min(alg.m, n)
539542
KsCache = _phiv_timestep_caches(u, maxiter, 4)
540-
return EPIRK4s3ACache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache)
543+
ts = Vector{typeof(dt)}(undef, 2)
544+
return EPIRK4s3ACache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache, ts)
541545
end
542546

543-
@cache struct EPIRK4s3BCache{uType, rateType, JCType, FType, matType, JType, KsType} <:
547+
@cache struct EPIRK4s3BCache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
544548
ExpRKCache
545549
u::uType
546550
uprev::uType
@@ -555,6 +559,7 @@ end
555559
J::JType
556560
B::matType
557561
KsCache::KsType
562+
ts::tsType
558563
end
559564
function alg_cache(
560565
alg::EPIRK4s3B, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -585,7 +590,8 @@ function alg_cache(
585590
# Allocate caches for phiv_timestep
586591
maxiter = min(alg.m, n)
587592
KsCache = _phiv_timestep_caches(u, maxiter, 4)
588-
return EPIRK4s3BCache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache)
593+
ts = Vector{typeof(dt)}(undef, 2)
594+
return EPIRK4s3BCache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache, ts)
589595
end
590596

591597
@cache struct EPIRK5s3Cache{uType, rateType, JCType, FType, matType, JType, KsType} <:
@@ -635,7 +641,7 @@ function alg_cache(
635641
return EPIRK5s3Cache(u, uprev, tmp, dz, k, rtmp, rtmp2, du1, jac_config, uf, J, B, KsCache)
636642
end
637643

638-
@cache struct EXPRB53s3Cache{uType, rateType, JCType, FType, matType, JType, KsType} <:
644+
@cache struct EXPRB53s3Cache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
639645
ExpRKCache
640646
u::uType
641647
uprev::uType
@@ -650,6 +656,7 @@ end
650656
J::JType
651657
B::matType
652658
KsCache::KsType
659+
ts::tsType
653660
end
654661
function alg_cache(
655662
alg::EXPRB53s3, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -680,10 +687,11 @@ function alg_cache(
680687
# Allocate caches for phiv_timestep
681688
maxiter = min(alg.m, n)
682689
KsCache = _phiv_timestep_caches(u, maxiter, 4)
683-
return EXPRB53s3Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache)
690+
ts = Vector{typeof(dt)}(undef, 2)
691+
return EXPRB53s3Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache, ts)
684692
end
685693

686-
@cache struct EPIRK5P1Cache{uType, rateType, JCType, FType, matType, JType, KsType} <:
694+
@cache struct EPIRK5P1Cache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
687695
ExpRKCache
688696
u::uType
689697
uprev::uType
@@ -698,6 +706,7 @@ end
698706
J::JType
699707
B::matType
700708
KsCache::KsType
709+
ts::tsType
701710
end
702711
function alg_cache(
703712
alg::EPIRK5P1, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -728,10 +737,11 @@ function alg_cache(
728737
# Allocate caches for phiv_timestep
729738
maxiter = min(alg.m, n)
730739
KsCache = _phiv_timestep_caches(u, maxiter, 3)
731-
return EPIRK5P1Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache)
740+
ts = Vector{typeof(dt)}(undef, 3)
741+
return EPIRK5P1Cache(u, uprev, tmp, dz, rtmp, rtmp2, du1, jac_config, uf, K, J, B, KsCache, ts)
732742
end
733743

734-
@cache struct EPIRK5P2Cache{uType, rateType, JCType, FType, matType, JType, KsType} <:
744+
@cache struct EPIRK5P2Cache{uType, rateType, JCType, FType, matType, JType, KsType, tsType} <:
735745
ExpRKCache
736746
u::uType
737747
uprev::uType
@@ -747,6 +757,7 @@ end
747757
J::JType
748758
B::matType
749759
KsCache::KsType
760+
ts::tsType
750761
end
751762
function alg_cache(
752763
alg::EPIRK5P2, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -777,7 +788,8 @@ function alg_cache(
777788
# Allocate caches for phiv_timestep
778789
maxiter = min(alg.m, n)
779790
KsCache = _phiv_timestep_caches(u, maxiter, 3)
780-
return EPIRK5P2Cache(u, uprev, tmp, dz, rtmp, rtmp2, dR, du1, jac_config, uf, K, J, B, KsCache)
791+
ts = Vector{typeof(dt)}(undef, 3)
792+
return EPIRK5P2Cache(u, uprev, tmp, dz, rtmp, rtmp2, dR, du1, jac_config, uf, K, J, B, KsCache, ts)
781793
end
782794

783795
####################################

‎lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_perform_step.jl‎

Lines changed: 45 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -801,11 +801,13 @@ end
801801

802802
function perform_step!(integrator, cache::Exp4Cache, repeat_step = false)
803803
(; t, dt, uprev, u, f, p) = integrator
804-
(; tmp, rtmp, rtmp2, K, J, B, KsCache) = cache
804+
(; tmp, rtmp, rtmp2, K, J, B, KsCache, ts) = cache
805805
calc_J!(J, integrator, cache)
806806
alg = unwrap_alg(integrator, true)
807807
f0 = integrator.fsalfirst # f(u0) is fsaled
808-
ts = [dt / 3, 2dt / 3, dt]
808+
ts[1] = dt / 3
809+
ts[2] = 2dt / 3
810+
ts[3] = dt
809811
kwargs = (
810812
tol = integrator.opts.reltol, iop = alg.iop,
811813
opnorm = integrator.opts.internalopnorm,
@@ -819,29 +821,31 @@ function perform_step!(integrator, cache::Exp4Cache, repeat_step = false)
819821
@inbounds for i in 1:3
820822
K[:, i] ./= ts[i]
821823
end
822-
mul!(rtmp, K, [-7 / 300, 97 / 150, -37 / 300]) # rtmp is now w4
824+
@views @.. broadcast = false rtmp = (-7 / 300) * K[:, 1] + (97 / 150) * K[:, 2] +
825+
(-37 / 300) * K[:, 3] # rtmp is now w4
823826
@muladd @.. broadcast = false tmp = uprev + dt * rtmp # tmp is now u4
824827
mul!(rtmp2, J, rtmp)
825828
f(rtmp, tmp, p, t + dt) # TODO: what should be the time?
826829
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
827830
@muladd @.. broadcast = false @view(B[:, 2]) = rtmp - f0 - dt * rtmp2 # B[:,2] is now d4
828831
# Partially update entities that use k1, k2, k3
829-
mul!(rtmp, K, [59 / 300, -7 / 75, 269 / 300]) # rtmp is now w7
832+
@views @.. broadcast = false rtmp = (59 / 300) * K[:, 1] + (-7 / 75) * K[:, 2] +
833+
(269 / 300) * K[:, 3] # rtmp is now w7
830834
@muladd @.. broadcast = false u = uprev + dt * @view(K[:, 3])
831835
# Krylov for the first remainder d4
832836
phiv_timestep!(K, ts, J, B; kwargs...)
833837
@inbounds for i in 1:3
834838
K[:, i] ./= ts[i]
835839
end
836-
mul!(rtmp2, K, [2 / 3, 2 / 3, 2 / 3])
840+
@views @.. broadcast = false rtmp2 = (2 / 3) * (K[:, 1] + K[:, 2] + K[:, 3])
837841
rtmp .+= rtmp2 # w7 fully updated
838842
@muladd @.. broadcast = false tmp = uprev + dt * rtmp # tmp is now u7
839843
mul!(rtmp2, J, rtmp)
840844
f(rtmp, tmp, p, t + dt) # TODO: what should be the time?
841845
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
842846
@muladd @.. broadcast = false @view(B[:, 2]) = rtmp - f0 - dt * rtmp2 # B[:,2] is now d7
843847
# Partially update entities that use k4, k5, k6
844-
mul!(rtmp, K, [1.0, -4 / 3, 1.0])
848+
@views @.. broadcast = false rtmp = K[:, 1] + (-4 / 3) * K[:, 2] + K[:, 3]
845849
axpy!(dt, rtmp, u)
846850
# Krylov for the second remainder d7
847851
k7 = @view(K[:, 1])
@@ -892,7 +896,7 @@ end
892896

893897
function perform_step!(integrator, cache::EPIRK4s3ACache, repeat_step = false)
894898
(; t, dt, uprev, u, f, p) = integrator
895-
(; tmp, rtmp, rtmp2, K, J, B, KsCache) = cache
899+
(; tmp, rtmp, rtmp2, K, J, B, KsCache, ts) = cache
896900
calc_J!(J, integrator, cache)
897901
alg = unwrap_alg(integrator, true)
898902
f0 = integrator.fsalfirst # f(u0) is fsaled
@@ -905,21 +909,23 @@ function perform_step!(integrator, cache::EPIRK4s3ACache, repeat_step = false)
905909

906910
# Compute U2 and U3 vertically
907911
B[:, 2] .= f0
908-
phiv_timestep!(K, [dt / 2, 2dt / 3], J, @view(B[:, 1:2]); kwargs...)
912+
ts[1] = dt / 2
913+
ts[2] = 2dt / 3
914+
phiv_timestep!(K, ts, J, @view(B[:, 1:2]); kwargs...)
909915
## U2 and R2
910916
@.. broadcast = false tmp = uprev + @view(K[:, 1]) # tmp is now U2
911917
f(rtmp, tmp, p, t + dt / 2)
912918
mul!(rtmp2, J, @view(K[:, 1]))
913919
@.. broadcast = false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
914-
B[:, 4] .= (32 / dt^2) * rtmp
915-
B[:, 5] .= (-144 / dt^3) * rtmp
920+
B[:, 4] .= (32 / dt^2) .* rtmp
921+
B[:, 5] .= (-144 / dt^3) .* rtmp
916922
## U3 and R3
917923
@.. broadcast = false tmp = uprev + @view(K[:, 2]) # tmp is now U3
918924
f(rtmp, tmp, p, t + 2dt / 3)
919925
mul!(rtmp2, J, @view(K[:, 2]))
920926
@.. broadcast = false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
921-
B[:, 4] .-= (13.5 / dt^2) * rtmp
922-
B[:, 5] .+= (81 / dt^3) * rtmp
927+
B[:, 4] .-= (13.5 / dt^2) .* rtmp
928+
B[:, 5] .+= (81 / dt^3) .* rtmp
923929
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2)
924930

925931
# Update u
@@ -972,7 +978,7 @@ end
972978

973979
function perform_step!(integrator, cache::EPIRK4s3BCache, repeat_step = false)
974980
(; t, dt, uprev, u, f, p) = integrator
975-
(; tmp, rtmp, rtmp2, K, J, B, KsCache) = cache
981+
(; tmp, rtmp, rtmp2, K, J, B, KsCache, ts) = cache
976982
calc_J!(J, integrator, cache)
977983
alg = unwrap_alg(integrator, true)
978984
f0 = integrator.fsalfirst # f(u0) is fsaled
@@ -986,23 +992,25 @@ function perform_step!(integrator, cache::EPIRK4s3BCache, repeat_step = false)
986992
# Compute U2 and U3 vertically
987993
fill!(@view(B[:, 2]), zero(eltype(B)))
988994
B[:, 3] .= f0
989-
phiv_timestep!(K, [dt / 2, 3dt / 4], J, @view(B[:, 1:3]); kwargs...)
995+
ts[1] = dt / 2
996+
ts[2] = 3dt / 4
997+
phiv_timestep!(K, ts, J, @view(B[:, 1:3]); kwargs...)
990998
K[:, 1] .*= 8 / (3 * dt)
991999
K[:, 2] .*= 16 / (9 * dt)
9921000
## U2 and R2
9931001
@.. broadcast = false tmp = uprev + @view(K[:, 1]) # tmp is now U2
9941002
f(rtmp, tmp, p, t + dt / 2)
9951003
mul!(rtmp2, J, @view(K[:, 1]))
9961004
@.. broadcast = false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
997-
B[:, 4] .= (54 / dt^2) * rtmp
998-
B[:, 5] .= (-324 / dt^3) * rtmp
1005+
B[:, 4] .= (54 / dt^2) .* rtmp
1006+
B[:, 5] .= (-324 / dt^3) .* rtmp
9991007
## U3 and R3
10001008
@.. broadcast = false tmp = uprev + @view(K[:, 2]) # tmp is now U3
10011009
f(rtmp, tmp, p, t + 3dt / 4)
10021010
mul!(rtmp2, J, @view(K[:, 2]))
10031011
@.. broadcast = false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
1004-
B[:, 4] .-= (16 / dt^2) * rtmp
1005-
B[:, 5] .+= (144 / dt^3) * rtmp
1012+
B[:, 4] .-= (16 / dt^2) .* rtmp
1013+
B[:, 5] .+= (144 / dt^3) .* rtmp
10061014
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2)
10071015

10081016
# Update u
@@ -1091,7 +1099,7 @@ function perform_step!(integrator, cache::EPIRK5s3Cache, repeat_step = false)
10911099
# Compute U3 horizontally
10921100
B[:, 2] .= (53 / 5) .* f0
10931101
B[:, 3] .= (-648 / (5 * dt)) .* f0
1094-
B[:, 4] .= (2916 / (5 * dt^2)) .* f0 + (32065 / (1152 * dt^2)) .* rtmp
1102+
B[:, 4] .= (2916 / (5 * dt^2)) .* f0 .+ (32065 / (1152 * dt^2)) .* rtmp
10951103
phiv_timestep!(k, 4dt / 9, J, @view(B[:, 1:4]); kwargs...)
10961104
## Update B matrix using R2
10971105
B[:, 2] .= f0
@@ -1161,7 +1169,7 @@ end
11611169

11621170
function perform_step!(integrator, cache::EXPRB53s3Cache, repeat_step = false)
11631171
(; t, dt, uprev, u, f, p) = integrator
1164-
(; tmp, rtmp, rtmp2, K, J, B, KsCache) = cache
1172+
(; tmp, rtmp, rtmp2, K, J, B, KsCache, ts) = cache
11651173
calc_J!(J, integrator, cache)
11661174
alg = unwrap_alg(integrator, true)
11671175
f0 = integrator.fsalfirst # f(u0) is fsaled
@@ -1174,7 +1182,9 @@ function perform_step!(integrator, cache::EXPRB53s3Cache, repeat_step = false)
11741182

11751183
# Compute the first group for U2 and U3
11761184
B[:, 2] .= f0
1177-
phiv_timestep!(K, [dt / 2, 9dt / 10], J, @view(B[:, 1:2]); kwargs...)
1185+
ts[1] = dt / 2
1186+
ts[2] = 9dt / 10
1187+
phiv_timestep!(K, ts, J, @view(B[:, 1:2]); kwargs...)
11781188
## U2 and R2
11791189
@.. broadcast = false tmp = uprev + @view(K[:, 1]) # tmp is now U2
11801190
f(rtmp, tmp, p, t + dt / 2)
@@ -1186,21 +1196,21 @@ function perform_step!(integrator, cache::EXPRB53s3Cache, repeat_step = false)
11861196
# Compute the second group for U3
11871197
fill!(@view(B[:, 2]), zero(eltype(B)))
11881198
B[:, 4] .= rtmp
1189-
phiv_timestep!(K, [dt / 2, 9dt / 10], J, @view(B[:, 1:4]); kwargs...)
1199+
phiv_timestep!(K, ts, J, @view(B[:, 1:4]); kwargs...)
11901200
## Update B using R2
11911201
B[:, 2] .= f0
11921202
B[:, 4] .= (18 / dt^2) .* rtmp
11931203
B[:, 5] .= (-60 / dt^3) .* rtmp
11941204
## U3 and R3
1195-
@views tmp .+= 216 / (25 * dt^2) .* K[:, 1] + 8 / dt^2 .* K[:, 2] # tmp is now U3
1205+
@views tmp .+= 216 / (25 * dt^2) .* K[:, 1] .+ 8 / dt^2 .* K[:, 2] # tmp is now U3
11961206
f(rtmp, tmp, p, t + 9dt / 10)
11971207
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
11981208
tmp .-= uprev
11991209
mul!(rtmp2, J, tmp)
12001210
@.. broadcast = false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
12011211
## Update B using R3
1202-
B[:, 4] .-= (250 / (81 * dt^2)) * rtmp
1203-
B[:, 5] .+= (500 / (27 * dt^3)) * rtmp
1212+
B[:, 4] .-= (250 / (81 * dt^2)) .* rtmp
1213+
B[:, 5] .+= (500 / (27 * dt^3)) .* rtmp
12041214

12051215
# Update u
12061216
du = @view(K[:, 1])
@@ -1268,7 +1278,7 @@ end
12681278

12691279
function perform_step!(integrator, cache::EPIRK5P1Cache, repeat_step = false)
12701280
(; t, dt, uprev, u, f, p) = integrator
1271-
(; tmp, rtmp, rtmp2, K, J, B, KsCache) = cache
1281+
(; tmp, rtmp, rtmp2, K, J, B, KsCache, ts) = cache
12721282
calc_J!(J, integrator, cache)
12731283
alg = unwrap_alg(integrator, true)
12741284
f0 = integrator.fsalfirst # f(u0) is fsaled
@@ -1291,7 +1301,10 @@ function perform_step!(integrator, cache::EPIRK5P1Cache, repeat_step = false)
12911301

12921302
# Compute the first column (f0)
12931303
B[:, 2] .= f0
1294-
phiv_timestep!(K, [g11, g21, g31], J, @view(B[:, 1:2]); kwargs...)
1304+
ts[1] = g11
1305+
ts[2] = g21
1306+
ts[3] = g31
1307+
phiv_timestep!(K, ts, J, @view(B[:, 1:2]); kwargs...)
12951308
## U1 and R1
12961309
@.. broadcast = false tmp = uprev + @view(K[:, 1]) # tmp is now U1
12971310
f(rtmp, tmp, p, t + g11)
@@ -1384,7 +1397,7 @@ end
13841397

13851398
function perform_step!(integrator, cache::EPIRK5P2Cache, repeat_step = false)
13861399
(; t, dt, uprev, u, f, p) = integrator
1387-
(; tmp, rtmp, rtmp2, dR, K, J, B, KsCache) = cache
1400+
(; tmp, rtmp, rtmp2, dR, K, J, B, KsCache, ts) = cache
13881401
calc_J!(J, integrator, cache)
13891402
alg = unwrap_alg(integrator, true)
13901403
f0 = integrator.fsalfirst # f(u0) is fsaled
@@ -1409,7 +1422,10 @@ function perform_step!(integrator, cache::EPIRK5P2Cache, repeat_step = false)
14091422

14101423
# Compute the first column (f0)
14111424
B[:, 2] .= f0
1412-
phiv_timestep!(K, [g11, g21, g31], J, @view(B[:, 1:2]); kwargs...)
1425+
ts[1] = g11
1426+
ts[2] = g21
1427+
ts[3] = g31
1428+
phiv_timestep!(K, ts, J, @view(B[:, 1:2]); kwargs...)
14131429
## U1 and R1
14141430
@.. broadcast = false tmp = uprev + @view(K[:, 1]) # tmp is now U1
14151431
f(rtmp, tmp, p, t + g11)

0 commit comments

Comments
 (0)