From 523b16b418e4f8fc8cc53c9720e95640e27025ec Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 16 Apr 2024 18:09:34 -0400
Subject: [PATCH] fix issues with aliasing, finish up tests and export
 algorithm

---
 src/OrdinaryDiffEq.jl                         |  2 +-
 src/perform_step/rkn_perform_step.jl          | 21 ++++++++++---------
 .../partitioned_methods_tests.jl              |  4 +++-
 3 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 96c9ce2499..5ca5ccf1e0 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/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
 
diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl
index b112677f26..64fd0677de 100644
--- a/src/perform_step/rkn_perform_step.jl
+++ b/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)
diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl
index e7eb7a663e..c75d4f7e86 100644
--- a/test/algconvergence/partitioned_methods_tests.jl
+++ b/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)