From 4b49b3188850fb33a90b6a7007f18e472d1cad10 Mon Sep 17 00:00:00 2001
From: Mikel <mikel.antonana@gmail.com>
Date: Wed, 7 Aug 2024 13:00:56 +0200
Subject: [PATCH] IRKGL16-SIMD

---
 test/runtests.jl | 77 ++++++++++++++++++++++++------------------------
 1 file changed, 39 insertions(+), 38 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 9112142..b9d2489 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,52 +1,52 @@
-using IRKGaussLegendre, ODEProblemLibrary, DiffEqDevTools, Test
+using .IRKGaussLegendre, ODEProblemLibrary, DiffEqDevTools, Test
 import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear
 
-function NbodyODE!(F, u, Gm, t)
+function NbodyODE!(F,u,Gm,t)
     N = length(Gm)
     for i in 1:N
-        for k in 1:3
-            F[k, i, 2] = 0
-        end
+       for k in 1:3
+           F[k, i, 2] = 0
+       end
     end
     for i in 1:N
-        xi = u[1, i, 1]
-        yi = u[2, i, 1]
-        zi = u[3, i, 1]
-        Gmi = Gm[i]
-        for j in (i + 1):N
-            xij = xi - u[1, j, 1]
-            yij = yi - u[2, j, 1]
-            zij = zi - u[3, j, 1]
-            Gmj = Gm[j]
-            dotij = (xij * xij + yij * yij + zij * zij)
-            auxij = 1 / (sqrt(dotij) * dotij)
-            Gmjauxij = Gmj * auxij
-            F[1, i, 2] -= Gmjauxij * xij
-            F[2, i, 2] -= Gmjauxij * yij
-            F[3, i, 2] -= Gmjauxij * zij
-            Gmiauxij = Gmi * auxij
-            F[1, j, 2] += Gmiauxij * xij
-            F[2, j, 2] += Gmiauxij * yij
-            F[3, j, 2] += Gmiauxij * zij
-        end
+       xi = u[1,i,1]
+       yi = u[2,i,1]
+       zi = u[3,i,1]
+       Gmi = Gm[i]
+       for j in i+1:N
+           xij = xi - u[1,j,1]
+           yij = yi - u[2,j,1]
+           zij = zi - u[3,j,1]
+           Gmj = Gm[j]
+           dotij = (xij*xij+yij*yij+zij*zij)
+           auxij = 1/(sqrt(dotij)*dotij)
+           Gmjauxij = Gmj*auxij
+           F[1,i,2] -= Gmjauxij*xij
+           F[2,i,2] -= Gmjauxij*yij
+           F[3,i,2] -= Gmjauxij*zij
+           Gmiauxij = Gmi*auxij
+           F[1,j,2] += Gmiauxij*xij
+           F[2,j,2] += Gmiauxij*yij
+           F[3,j,2] += Gmiauxij*zij
+       end
     end
     for i in 1:3, j in 1:N
-        F[i, j, 1] = u[i, j, 2]
+       F[i,j,1] = u[i,j,2]
     end
-    return nothing
+   return nothing
 end
 
 Gm = [5, 4, 3]
-N = length(Gm)
-q = [1, -1, 0, -2, -1, 0, 1, 3, 0]
-v = zeros(size(q))
-q0 = reshape(q, 3, :)
-v0 = reshape(v, 3, :)
-u0 = Array{Float64}(undef, 3, N, 2)
-u0[:, :, 1] = q0
-u0[:, :, 2] = v0
-tspan = (0.0, 63.0)
-prob = ODEProblem(NbodyODE!, u0, tspan, Gm);
+N=length(Gm)
+q=[1,-1,0,-2,-1,0,1,3,0]
+v=zeros(size(q))
+q0 = reshape(q,3,:)
+v0 = reshape(v,3,:)
+u0 = Array{Float64}(undef,3,N,2)
+u0[:,:,1] = q0
+u0[:,:,2] = v0
+tspan = (0.0,63.0)
+prob=ODEProblem(NbodyODE!,u0,tspan,Gm);
 sol1 = solve(prob, IRKGL16(), reltol = 1e-12, abstol = 1e-12)
 
 # Analytical tests
@@ -54,7 +54,8 @@ sol1 = solve(prob, IRKGL16(), reltol = 1e-12, abstol = 1e-12)
 sol = solve(prob_ode_2Dlinear, IRKGL16())
 @test sol.errors[:l2] < 1e-16
 
-dts = (1 // 2) .^ (4:-1:2)
+# dts = (1 // 2) .^ (4:-1:2)
+dts = BigFloat.( (1 // 2) .^ (4:-1:2))
 sim = test_convergence(dts, prob_ode_bigfloat2Dlinear, IRKGL16())
 @test abs(sim.𝒪est[:final] - 16) < 0.5