fix issues with aliasing, finish up tests and export algorithm

src/OrdinaryDiffEq.jl

@@ -403,7 +403,7 @@ export SplitEuler
 
 export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent,
        Nystrom5VelocityIndependent,
-       IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7
+       IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, RKN4
 
 export ROCK2, ROCK4, RKC, IRKC, ESERK4, ESERK5, SERK2
 

src/perform_step/rkn_perform_step.jl

@@ -1838,7 +1838,7 @@ end
 @muladd function perform_step!(integrator, cache::RKN4ConstantCache, repeat_step = false)
     @unpack t, dt, f, p = integrator
     duprev, uprev = integrator.uprev.x
-
+    u, du = integrator.u.x
     #define dt values
     halfdt = dt/2
     dtsq = dt^2
@@ -1851,6 +1851,7 @@ end
     #perform operations to find k values
     k₁ = integrator.fsalfirst.x[1]
     ku = uprev + halfdt * duprev + eightdtsq * k₁
+    print(ku)
     kdu = duprev + halfdt * k₁
 
     k₂ = f.f1(kdu, ku, p, ttmp)
@@ -1876,7 +1877,7 @@ end
 @muladd function perform_step!(integrator, cache::RKN4Cache, repeat_step = false)
     @unpack t, dt, f, p = integrator
     duprev, uprev = integrator.uprev.x
-
+    du, u = integrator.u.x
     @unpack tmp, fsalfirst, k₂, k₃, k = cache
     kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
 
@@ -1891,20 +1892,20 @@ end
 
     #perform operations to find k values
     k₁ = integrator.fsalfirst.x[1]
-    ku = uprev + halfdt * duprev + eightdtsq * k₁
-    kdu = duprev + halfdt * k₁
+    @.. broadcast=false ku = uprev + halfdt * duprev + eightdtsq * k₁
+    @.. broadcast=false kdu = duprev + halfdt * k₁
 
     f.f1(k₂, kdu, ku, p, ttmp)
-    ku = uprev + dt * duprev + halfdtsq * k₂
-    kdu = duprev + dt * k₂
+    @.. broadcast=false ku = uprev + dt * duprev + halfdtsq * k₂
+    @.. broadcast=false kdu = duprev + dt * k₂
 
     f.f1(k₃, kdu, ku, p, t + dt)
-    ku = uprev + dt * duprev + eightdtsq * k₃
-    kdu = duprev + dt * k₃
+    @.. broadcast=false ku = uprev + dt * duprev + eightdtsq * k₃
+    @.. broadcast=false kdu = duprev + dt * k₃
 
     #perform final calculations to determine new y and y'.
-    u = uprev + sixthdtsq* (1*k₁ + 2*k₂ + 0*k₃) + dt * duprev
-    du = duprev + sixthdt * (1*k₁ + 4*k₂ + 1*k₃)
+    @.. broadcast=false u = uprev + sixthdtsq* (1*k₁ + 2*k₂ + 0*k₃) + dt * duprev
+    @.. broadcast=false du = duprev + sixthdt * (1*k₁ + 4*k₂ + 1*k₃)
 
     f.f1(k.x[1], du, u, p, t + dt)
     f.f2(k.x[2], du, u, p, t + dt)

test/algconvergence/partitioned_methods_tests.jl

@@ -111,9 +111,11 @@ sim = test_convergence(dts, prob, KahanLi8(), dense_errors = true)
 
 sol = solve(prob, Nystrom4(), dt = 1 / 1000)
 
+sol = solve(prob, RKN4(), dt = 1/1000)
+
 # Nyström method
 dts = 1 .// 2 .^ (9:-1:6)
-sim = test_convergence(dts, prob, OrdinaryDiffEq.RKN4(), dense_errors = true)
+sim = test_convergence(dts, prob, RKN4(), dense_errors = true)
 @test sim.𝒪est[:l2]≈4 rtol=1e-1
 @test sim.𝒪est[:L2]≈4 rtol=1e-1
 sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true)