adding error code for when step size is NaN

QuEraComputing · Jan 2, 2024 · b2d36ca · b2d36ca
1 parent ca081ff
commit b2d36ca
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 17 deletions.
diff --git a/src/dp5/solver.jl b/src/dp5/solver.jl
@@ -98,11 +98,10 @@ function dopri5(
 
     ###### Basic Integration Step
     for _ in 1:solver.options.maximum_allowed_steps
-        # if nstep > solver.options.maximum_allowed_steps
-        #     # GOTO 78
-        #     # println(" MORE THAN NMAX = ", solver.options.maximum_allowed_steps, " STEPS ARE NEEDED")
-        #     return h, Report(solver.vars.x, DormandPrince.INPUT_CHECKS_SUCCESSFUL, DormandPrince.LARGER_NMAX_NEEDED , 0, 0, 0, 0)
-        # end
+        if isnan(h)
+            idid = DormandPrince.STEP_SIZE_BECOMES_NAN
+            break
+        end
 
         if (0.10 * abs(h)) <= abs(solver.vars.x)*solver.options.uround 
             # GOTO 77

diff --git a/src/dp8/solver.jl b/src/dp8/solver.jl
@@ -97,12 +97,11 @@ function dop853(
 
     ###### Basic Integration Step
     for _ in 1:solver.options.maximum_allowed_steps
-        # if nstep > solver.options.maximum_allowed_steps
-        #     # GOTO 78
-        #     # println(" MORE THAN NMAX = ", solver.options.maximum_allowed_steps, " STEPS ARE NEEDED")
-        #     return h, Report(solver.vars.x, DormandPrince.INPUT_CHECKS_SUCCESSFUL, DormandPrince.LARGER_NMAX_NEEDED , 0, 0, 0, 0)
-        # end
-
+        if isnan(h)
+            idid = DormandPrince.STEP_SIZE_BECOMES_NAN
+            break
+        end
+
         if (0.10 * abs(h)) <= abs(solver.vars.x)*solver.options.uround 
             # GOTO 77
             # println("STEP SIZE TOO SMALL, H = ", h)

diff --git a/src/types.jl b/src/types.jl
@@ -6,6 +6,7 @@ abstract type AbstractDPSolver{T <: Real, StateType <: AbstractVector, F} end
     INPUT_NOT_CONSISTENT = -1 # use for check failures in the beginning of  core_integrator call
     LARGER_NMAX_NEEDED = -2
     STEP_SIZE_BECOMES_TOO_SMALL = -3
+    STEP_SIZE_BECOMES_NAN = -4
 end
 
 @enum Checks begin

diff --git a/test/debug.jl b/test/debug.jl
@@ -1,5 +1,5 @@
 using LinearAlgebra
-using DormandPrince: DP5Solver, DP8Solver, integrate
+using DormandPrince: DP8Solver, DP5Solver, integrate
 
 function evolution_operator(t::Float64)
     ϕ = 2.2 * sin(π * t)^2
@@ -33,13 +33,20 @@ function run()
         ComplexF64[1.0, 0.0]
     )
 
-    report = integrate(solver, 1.0)
+    report = integrate(solver, 2π)
 
-    println(report.num_rejected_steps)
-    println(norm(solver.y))
+    if !(solver.y ≈ solution(2π))
+        println("failed, $(solver.y) != $(solution(2π)), $(report.idid)")
+    end
+    @assert solver.y ≈ solution(2π)
+    # println(report.num_rejected_steps)
+    # println(norm(solver.y))
     solver.vars.h
-    @assert solver.y ≈ solution(1.0)
+    solver.y
 
 end
 
-run()
+for i in 1:1000
+    run()
+    println("-------------------------------------------------------------------------------")
+end