completed restructuring and code refactor

src/DormandPrince.jl

@@ -1,10 +1,12 @@
 module DormandPrince
-using Base.Iterators:repeated, Repeated
 
-include("types.jl")
-include("solver.jl")
-include("integrate.jl")
-include("checks.jl")
-include("helpers.jl")
+
+# internal imports
+include("dp5_core/mod.jl")
+include("interface.jl")
+
+# export Interface
+export DP5Solver, integrate
+
 
 end # DormandPrince

src/dp5_core/mod.jl

@@ -0,0 +1,14 @@
+module DP5Core
+
+# external imports
+using Base.Iterators:repeated, Repeated
+
+include("types.jl")
+include("helpers.jl")
+include("checks.jl")
+include("solver.jl")
+
+export DP5Solver, dopri5
+
+
+end

src/solver.jl → src/dp5_core/solver.jl

@@ -94,18 +94,23 @@ function dopcor(
     nfcn += 2
     reject = false
 
+    idid = LARGER_NMAX_NEEDED
+
     ###### Basic Integration Step
-    while true
-        if nstep > solver.options.maximum_allowed_steps
-            # GOTO 78
-            # println(" MORE THAN NMAX = ", solver.options.maximum_allowed_steps, " STEPS ARE NEEDED")
-            return h, DP5Report(solver.vars.x, CHECKS_SUCCESSFUL, LARGER_NMAX_NEEDED , 0, 0, 0, 0)
-        end
+    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, DP5Report(solver.vars.x, INPUT_CHECKS_SUCCESSFUL, LARGER_NMAX_NEEDED , 0, 0, 0, 0)
+        # end
 
         if (0.10 * abs(h)) <= abs(solver.vars.x)*solver.options.uround 
             # GOTO 77
             # println("STEP SIZE TOO SMALL, H = ", h)
-            return h, DP5Report(solver.vars.x, CHECKS_SUCCESSFUL, STEP_SIZE_BECOMES_TOO_SMALL, 0, 0, 0, 0)
+            # return h, DP5Report(solver.vars.x, INPUT_CHECKS_SUCCESSFUL, STEP_SIZE_BECOMES_TOO_SMALL, 0, 0, 0, 0)
+
+            idid = STEP_SIZE_BECOMES_TOO_SMALL
+            break
         end
 
         if ((solver.vars.x+1.01*h-xend)*posneg) > 0.0
@@ -144,15 +149,17 @@ function dopcor(
             ###### Normal Exit
             if solver.vars.last 
                 h = hnew
-                return h, DP5Report(
-                    solver.vars.x, 
-                    CHECKS_SUCCESSFUL,
-                    COMPUTATION_SUCCESSFUL, 
-                    nfcn, 
-                    nstep, 
-                    naccpt, 
-                    nrejct
-                )
+                idid = COMPUTATION_SUCCESSFUL
+                break
+                # return h, DP5Report(
+                #     solver.vars.x, 
+                #     INPUT_CHECKS_SUCCESSFUL,
+                #     COMPUTATION_SUCCESSFUL, 
+                #     nfcn, 
+                #     nstep, 
+                #     naccpt, 
+                #     nrejct
+                # )
             end
 
             if(abs(hnew) > hmax)
@@ -176,6 +183,16 @@ function dopcor(
         h = hnew
     end
 
+    return h, DP5Report(
+        solver.vars.x, 
+        INPUT_CHECKS_SUCCESSFUL,
+        idid, 
+        nfcn, 
+        nstep, 
+        naccpt, 
+        nrejct
+    )
+
 end
 
 function hinit(

src/types.jl → src/dp5_core/types.jl

@@ -8,7 +8,7 @@
 end
 
 @enum Checks begin
-    CHECKS_SUCCESSFUL
+    INPUT_CHECKS_SUCCESSFUL
     MAX_ALLOWED_STEPS_NEGATIVE
     UNSUPPORTED_UROUND
     CURIOUS_BETA 
@@ -92,15 +92,21 @@ struct DP5Solver{StateType <: AbstractVector, T <: Real, F}
     consts::DP5Consts{T}
     vars::DP5Vars{T}
 
-    function DP5Solver(f::F, x::T, y::StateType; kw...) where {StateType <: AbstractVector, T<:Real, F}
-        k1 = copy(y)
-        k2 = copy(y)
-        k3 = copy(y)
-        k4 = copy(y)
-        k5 = copy(y)
-        k6 = copy(y)
-        y1 = copy(y)
-        ysti = copy(y)
+    function DP5Solver(
+        f::F, 
+        x::T, 
+        y::StateType,
+        k1::StateType,
+        k2::StateType,
+        k3::StateType,
+        k4::StateType,
+        k5::StateType,
+        k6::StateType,
+        y1::StateType,
+        ysti::StateType; kw...) where {StateType <: AbstractVector, T<:Real, F}
+
+        #TODO: check if y, k1, k2, k3, k4, k5, k6, y1, ysti have the same length
+
         options = DP5Options{T}(;kw...)
         consts = DP5Consts(options)
         vars = DP5Vars{T}(;x=x, h=options.initial_step_size)
@@ -109,16 +115,17 @@ struct DP5Solver{StateType <: AbstractVector, T <: Real, F}
     end
 end
 
-
-# DP5Options() = DP5Options(
-#     2.3e-16, #uround
-#     0.9,     # safety_factor
-#     0.2,     # step_size_selection_one
-#     10.0,    # step_size_selection_two
-#     0.04,    # beta
-#     0.0,     # maximal step size - default to 0.0, later set to xend - x
-#     0.0,     # initial step size - default to 0.0, trigger hinit later
-#     100000,  # maximum number of allowed steps 
-#     true,    # whether or not error messages should be printed
-#     1000 # stiffness test activated after step J * this number
-# )
+function DP5Solver(
+    f::F, 
+    x::T, 
+    y::StateType; kw...) where {StateType <: AbstractVector, T<:Real, F}
+    k1 = copy(y)
+    k2 = copy(y)
+    k3 = copy(y)
+    k4 = copy(y)
+    k5 = copy(y)
+    k6 = copy(y)
+    y1 = copy(y)
+    ysti = copy(y)
+    DP5Solver(f, x, y, k1, k2, k3, k4, k5, k6, y1, ysti;kw...)
+end

src/integrate.jl → src/interface.jl

@@ -1,3 +1,5 @@
+using DormandPrince.DP5Core: DP5Solver, dopri5
+
 mutable struct DP5Iterator{T <: Real}
     solver::DP5Solver
     times::AbstractVector{T}

test/runtests.jl

@@ -1,6 +1,7 @@
 using Test
 using LinearAlgebra
 using DormandPrince: DP5Solver, dopri5, integrate
+# test error paths of dopcor, bypass checks if necessary 
 
 include("exact_evol_helpers.jl")
 
@@ -20,7 +21,7 @@ end
         ComplexF64[1.0, 0.0]
     )
 
-    dopri5(solver, 2π)
+    integrate(solver, 2π)
 
     @test solver.y ≈ solution(2π)
 end