From 9c2d37399b4800f95399dff0e93bae69adf7dcb6 Mon Sep 17 00:00:00 2001 From: John Long Date: Fri, 15 Dec 2023 11:08:28 -0500 Subject: [PATCH] completed restructuring and code refactor --- src/DormandPrince.jl | 14 ++++---- src/{ => dp5_core}/checks.jl | 0 src/{ => dp5_core}/helpers.jl | 0 src/dp5_core/mod.jl | 14 ++++++++ src/{ => dp5_core}/solver.jl | 49 ++++++++++++++++++--------- src/{ => dp5_core}/types.jl | 53 +++++++++++++++++------------- src/{integrate.jl => interface.jl} | 2 ++ test/runtests.jl | 3 +- 8 files changed, 89 insertions(+), 46 deletions(-) rename src/{ => dp5_core}/checks.jl (100%) rename src/{ => dp5_core}/helpers.jl (100%) create mode 100644 src/dp5_core/mod.jl rename src/{ => dp5_core}/solver.jl (83%) rename src/{ => dp5_core}/types.jl (76%) rename src/{integrate.jl => interface.jl} (97%) diff --git a/src/DormandPrince.jl b/src/DormandPrince.jl index 832347f..e404696 100644 --- a/src/DormandPrince.jl +++ b/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 diff --git a/src/checks.jl b/src/dp5_core/checks.jl similarity index 100% rename from src/checks.jl rename to src/dp5_core/checks.jl diff --git a/src/helpers.jl b/src/dp5_core/helpers.jl similarity index 100% rename from src/helpers.jl rename to src/dp5_core/helpers.jl diff --git a/src/dp5_core/mod.jl b/src/dp5_core/mod.jl new file mode 100644 index 0000000..8d67b76 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/solver.jl b/src/dp5_core/solver.jl similarity index 83% rename from src/solver.jl rename to src/dp5_core/solver.jl index 5540843..9b293ca 100644 --- a/src/solver.jl +++ b/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( diff --git a/src/types.jl b/src/dp5_core/types.jl similarity index 76% rename from src/types.jl rename to src/dp5_core/types.jl index 1b25f39..f6f3742 100644 --- a/src/types.jl +++ b/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 -# ) \ No newline at end of file +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 diff --git a/src/integrate.jl b/src/interface.jl similarity index 97% rename from src/integrate.jl rename to src/interface.jl index 3066b8c..e318b29 100644 --- a/src/integrate.jl +++ b/src/interface.jl @@ -1,3 +1,5 @@ +using DormandPrince.DP5Core: DP5Solver, dopri5 + mutable struct DP5Iterator{T <: Real} solver::DP5Solver times::AbstractVector{T} diff --git a/test/runtests.jl b/test/runtests.jl index 0c6dda0..bc91a95 100644 --- a/test/runtests.jl +++ b/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