solver is fixed.

QuEraComputing · Dec 22, 2023 · 1a1af14 · 1a1af14
1 parent 13db0bc
commit 1a1af14
Show file tree

Hide file tree

Showing 3 changed files with 41 additions and 37 deletions.
diff --git a/src/dp8/helpers.jl b/src/dp8/helpers.jl
@@ -123,8 +123,6 @@ function error_estimation(solver::DP8Solver{T}, h::T) where T
         sk = atoli + rtoli*max(abs(yi), abs(k5i))
         erri2 = k4i - bhh1*k1i - bhh2*k9i - bhh3*k3i
         erri1 = er1*k1i + er6*k6i + er7*k7i + er8*k8i + er9*k9i + er10*k10i + er11*k2i + er12*k3i
-        println("sk = $sk, erri1 = $erri1, erri2 = $erri2")
-
         ((abs(erri1)/sk)^2, (abs(erri2)/sk)^2)
     end
 
@@ -134,8 +132,6 @@ function error_estimation(solver::DP8Solver{T}, h::T) where T
     den0 = err1 + 0.01*err2
     den0 = den0 <= 0.0 ? 1.0 : den0
     err = abs(h)*err1*sqrt(1.0/(length(solver.y) * den0))
-    println("err = $err, err1 = $err1, err2 = $err2")
-    println()
     return err 
 end
 

diff --git a/src/dp8/solver.jl b/src/dp8/solver.jl
@@ -64,10 +64,10 @@ function  DormandPrince.integrate(
 end
 
 function dop853(
-    solver, # contains f, x, y, k1, k2, k3, k4, k5, k6, y1, ysti, options
-    xend, 
-    hmax,
-    h, 
+    solver::DP8Solver{T}, # contains f, x, y, k1, k2, k3, k4, k5, k6, y1, ysti, options
+    xend::T, 
+    hmax::T,
+    h::T, 
 ) where T
     ##### Initializations
     # replace sign with Julia-native Base.sign
@@ -97,7 +97,6 @@ function dop853(
 
     ###### Basic Integration Step
     for _ in 1:solver.options.maximum_allowed_steps
-        println("x = $(solver.vars.x), h = $h, y = $(solver.y)")
         # if nstep > solver.options.maximum_allowed_steps
         #     # GOTO 78
         #     # println(" MORE THAN NMAX = ", solver.options.maximum_allowed_steps, " STEPS ARE NEEDED")
@@ -134,7 +133,6 @@ function dop853(
         ###### we require fac1 <= hnew/h <= fac2
         fac = max(solver.consts.facc2, min(solver.consts.facc1, fac/solver.options.safety_factor)) # facc1, facc2, fac must be Float64 
         hnew = h/fac
-        println("hnew = $hnew, h = $h, fac = $fac, fac11 = $fac11, err = $err")
         if err <= 1.0 
             ###### Step is accepted
             solver.vars.facold = max(err, 1e-4)

diff --git a/test/interface.jl b/test/interface.jl
@@ -29,7 +29,7 @@ end
 # standalone solver test
 @testset "Integration Test" begin
 
-    for SolverType in [DP5Solver]
+    for SolverType in [DP5Solver, DP8Solver]
         solver = SolverType(
             fcn,
             0.0,
@@ -50,52 +50,62 @@ end
     @testset "Iterator Interface" begin
         times = [0.1, 0.5, 1.1]
         exact_values = []
-        dp5_values = []
 
         # get exact values
         for t in times
             push!(exact_values, solution(t))
         end
 
         # use iterator to get exact values
-        solver = DP5Solver(
-            fcn,
-            0.0,
-            ComplexF64[1.0, 0.0]
-        )
-
-        iter = integrate(solver, times)
-
-        for (t,y) in iter
-            push!(dp5_values, copy(y))
+
+        for SolverType in [DP5Solver, DP8Solver]
+            values = []
+            solver = SolverType(
+                fcn,
+                0.0,
+                ComplexF64[1.0, 0.0]
+            )
+
+            iter = integrate(solver, times)
+
+            for (t,y) in iter
+                push!(values, copy(y))
+            end
+
+            @test values ≈ exact_values
         end
 
-        @test dp5_values ≈ exact_values
+
     end
 
     @testset "Callback Interface" begin
         times = [0.1, 0.5, 1.1]
-        callback_times = []
         exact_values = []
-        dp5_values = []
-
+
         # get exact values
         for t in times
             push!(exact_values, solution(t))
         end
 
-        # use iterator to get exact values
-        solver = DP5Solver(
-            fcn,
-            0.0,
-            ComplexF64[1.0, 0.0]
-        )
 
-        integrate(solver, times) do t, y
-            push!(callback_times, t)
-            push!(dp5_values, copy(y))
-        end
+        for SolverType in [DP5Solver, DP8Solver]
+            values = []
+            callback_times = []
+
+            # use iterator to get exact values
+            solver = SolverType(
+                fcn,
+                0.0,
+                ComplexF64[1.0, 0.0]
+            )
 
-        @test dp5_values ≈ exact_values
+            integrate(solver, times) do t, y
+                push!(callback_times, t)
+                push!(values, copy(y))
+            end
+
+            @test values ≈ exact_values
+            @test callback_times == times
+        end
     end
 end