diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..6fddca0 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +version: 2 +updates: + # Maintain dependencies for GitHub Actions + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..0e2b08f --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,64 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - "1" + - "1.9" + # - "nightly" + os: + - ubuntu-latest + # - macos-latest + # - windows-latest + arch: + - x86 + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v2 + with: + files: lcov.info + docs: + name: Documentation + runs-on: ubuntu-latest + permissions: + contents: write + statuses: write + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: '1' + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + - run: | + julia --project=docs -e ' + using Documenter: DocMeta, doctest + using DormandPrince + DocMeta.setdocmeta!(DormandPrince, :DocTestSetup, :(using DormandPrince); recursive=true) + doctest(DormandPrince)' diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..cba9134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..464b9ea --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,44 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == '{{{TRIGGER}}}' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: {{{TOKEN}}} + {{#SSH}} + ssh: {{{SSH}}} + {{/SSH}} + {{#SSH_PASSWORD}} + ssh_password: {{{SSH_PASSWORD}}} + {{/SSH_PASSWORD}} + {{#CHANGELOG}} + changelog: {{{CHANGELOG}}} + {{/CHANGELOG}} + {{#CHANGELOG_IGNORE}} + changelog_ignore: {{{CHANGELOG_IGNORE}}} + {{/CHANGELOG_IGNORE}} + {{#GPG}} + gpg: {{{GPG}}} + {{/GPG}} + {{#GPG_PASSWORD}} + gpg_password: {{{GPG_PASSWORD}}} + {{/GPG_PASSWORD}} + {{#REGISTRY}} + registry: {{{REGISTRY}}} + {{/REGISTRY}} + {{#BRANCHES}} + branches: {{{BRANCHES}}} + {{/BRANCHES}} + {{#DISPATCH}} + dispatch: {{{DISPATCH}}} + {{/DISPATCH}} + {{#DISPATCH_DELAY}} + dispatch_delay: {{{DISPATCH_DELAY}}} + {{/DISPATCH_DELAY}} diff --git a/.gitignore b/.gitignore index 2e528c7..eff5f6d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,11 @@ -.vscode/ -.CondaPkg/ +/.vscode +**/**/**.cov + +/Manifest.toml +/docs/build/ +/docs/Manifest.toml +/docs/src/assets/main.css +/docs/src/assets/indigo.css + .DS_Store -Manifest.toml \ No newline at end of file +.CondaPkg diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..0735f0b --- /dev/null +++ b/LICENSE @@ -0,0 +1,13 @@ +Copyright 2023 John Long + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/Project.toml b/Project.toml index 2b33c23..c273fd0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,9 @@ +name = "DormandPrince" +uuid = "5e45e72d-22b8-4dd0-9c8b-f96714864bcd" +authors = ["John Long",] +version = "0.1.0" + [deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1" -DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" -JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" + +[compat] +julia = "1.9" diff --git a/README.md b/README.md new file mode 100644 index 0000000..7a97971 --- /dev/null +++ b/README.md @@ -0,0 +1,25 @@ +# DormandPrince + +Julia-native implementation of the Dormand-Prince 5th order solver + +## Installation + +

+DormandPrince is a   + + + Julia Language + +   package. To install DormandPrince, + please open + Julia's interactive session (known as REPL) and press ] + key in the REPL to use the package mode, then type the following command +

+ +```julia +pkg> add DormandPrince +``` + +## License + +Apache License 2.0 diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml new file mode 100644 index 0000000..541dee7 --- /dev/null +++ b/benchmarks/Project.toml @@ -0,0 +1,3 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" diff --git a/benchmarks/dp5_profile_calls.jl b/benchmarks/dp5_profile_calls.jl deleted file mode 100644 index 7a10aee..0000000 --- a/benchmarks/dp5_profile_calls.jl +++ /dev/null @@ -1,31 +0,0 @@ -include("../solver.jl") -include("../types.jl") - -# fcn(x,y,f) -## f is the end result of the derivative itself -## x is the dependent variable -## y is y(x) -function fcn(x,y,f) - f[1] = -15.0*y[1] -end - -@kwdef mutable struct CounterFunc{F} - fun::F - counter::Int64 = 0 -end - -function (f::CounterFunc)(x...) - f.counter += 1 - return f.fun(x...) -end - -f = CounterFunc(fun = fcn) - -solver = DP5Solver( - f, - 0.0, # initial x - [1.0] # initial y at x -) - -dopri5(solver, 0.5) - diff --git a/benchmarks/ordinarydiffeq_profile_calls.jl b/benchmarks/ordinarydiffeq_profile_calls.jl deleted file mode 100644 index c06ccde..0000000 --- a/benchmarks/ordinarydiffeq_profile_calls.jl +++ /dev/null @@ -1,26 +0,0 @@ -using DifferentialEquations -call_count = 0 - -function f(u, p, t) - global call_count += 1 - return -15.0 * u -end - -#=- -@kwdef mutable struct CounterFunc{F} - fun::F - counter::Int64 = 0 -end - -function (f::CounterFunc)(x...) - f.counter += 1 - return f.fun(x...) -end - -counted_f = CounterFunc(fun = f) -=# - -u0 = 1 -tspan = (0.0, 0.5) -prob = ODEProblem(f, u0, tspan) -sol = solve(prob, DP5(), reltol = 1e-10, abstol = 1e-10) \ No newline at end of file diff --git a/benchmarks/rabi-diffeq.jl b/benchmarks/rabi-diffeq.jl new file mode 100644 index 0000000..1c2de1a --- /dev/null +++ b/benchmarks/rabi-diffeq.jl @@ -0,0 +1,19 @@ +using DifferentialEquations +using BenchmarkTools + + +function f(du, u, p, t) + g(t) = 2.2*2π*sin(2π*t) + + du[1] = -1im * g(t)/2 * u[2] + du[2] = -1im * g(t)/2 * u[1] +end + +u0 = ComplexF64[1.0, 0.0] +tspan = (0.0, 2π) +prob = ODEProblem(f, u0, tspan) +# get precompilation out of the way +sol = solve(prob, DP5(), reltol=1e-10, abstol=1e-10) + +# terminate benchmark after maximum of 5 minutes +@benchmark solve(prob, DP5(), reltol=1e-10, abstol=1e-10) samples=10000 evals=5 seconds=300 \ No newline at end of file diff --git a/benchmarks/rabi-dormand-prince.jl b/benchmarks/rabi-dormand-prince.jl new file mode 100644 index 0000000..097d72c --- /dev/null +++ b/benchmarks/rabi-dormand-prince.jl @@ -0,0 +1,19 @@ +using BenchmarkTools +using DormandPrince:DP5Solver, dopri5 + +function fcn(x, y, f) + g(x) = 2.2*2π*sin(2π*x) + + f[1] = -1im * g(x)/2 * y[2] + f[2] = -1im * g(x)/2 * y[1] +end + +solver = DP5Solver( + fcn, + 0.0, + ComplexF64[1.0, 0.0] +) + +dopri5(solver, 2π) + +@benchmark dopri5(clean_solver, 2π) setup=(clean_solver = DP5Solver(fcn, 0.0, ComplexF64[1.0, 0.0])) samples=10000 evals=5 seconds=500 \ No newline at end of file diff --git a/benchmarks/single_atom_bloqade.jl b/benchmarks/single_atom_bloqade.jl deleted file mode 100644 index 078ce85..0000000 --- a/benchmarks/single_atom_bloqade.jl +++ /dev/null @@ -1,19 +0,0 @@ -using Bloqade -using BenchmarkTools - -function generate_clean_problem() - atoms = generate_sites(ChainLattice(), 1, scale = 5.74) - h = rydberg_h(atoms; Ω = 20 * 2π, Δ = 0) - reg = zero_state(1) - - prob = SchrodingerProblem(reg, 1.6, h; algo=DP8()) - - return prob -end - -emulate!(generate_clean_problem()) - -# b = @benchmark emulate!(clean_problem) samples=100 evals=1 seconds=172800 setup=(clean_problem=SchrodingerProblem(deepcopy(reg), 1.6, h)) -b = @benchmark emulate!(clean_problem) samples=10000 evals=5 seconds=172800 setup=(clean_problem=generate_clean_problem()) - -# statevec(prob.reg) \ No newline at end of file diff --git a/benchmarks/single_atom_dp5.jl b/benchmarks/single_atom_dp5.jl deleted file mode 100644 index 36c1641..0000000 --- a/benchmarks/single_atom_dp5.jl +++ /dev/null @@ -1,43 +0,0 @@ -include("../solver.jl") -include("../types.jl") -using Bloqade -using BloqadeExpr: Hamiltonian -using BenchmarkTools - -function generate_clean_solver() - nsites = 1; - atoms = generate_sites(ChainLattice(), nsites, scale = 5.74) - h = rydberg_h(atoms; Ω = 20 * 2π, Δ = 0) - reg = zero_state(1) - # stop short of creating the SchrodingerProblem, we want the SchrodignerEquation - # which we can than trivially wrap with the fcn(n,x,y,f) that DP5 accepts - - # Generate SchrodingerEquation - state = statevec(reg) - space = YaoSubspaceArrayReg.space(reg) - T = real(eltype(state)) - T = isreal(h) ? T : Complex{T} - eq = SchrodingerEquation(h, Hamiltonian(T, h, space)) - - function fcn(x, y, f) - eq(f, y, nothing, x) - end - - y = statevec(reg) - - solver = DP5Solver( - fcn, - 0.0, - y - ) - - return solver -end - -# invoke eq via eq(dstate, state, p, t::Number) - -# trigger precompilation -dopri5(generate_clean_solver(), 1.6) - -# b = @benchmark dopri5(clean_solver, 1.6) samples=100 evals=1 seconds=172800 setup=(clean_solver = DP5Solver(fcn, 0.0, deepcopy(y))) -@benchmark dopri5(clean_solver, 1.6) samples=10000 evals=5 seconds=172800 setup=(clean_solver = generate_clean_solver()) \ No newline at end of file diff --git a/benchmarks/struct_profiler.jl b/benchmarks/struct_profiler.jl deleted file mode 100644 index 8c46cc9..0000000 --- a/benchmarks/struct_profiler.jl +++ /dev/null @@ -1,12 +0,0 @@ -# recommended by Phillip Weinberg -# for profiling the number of calls to the ODE function - -@kwdef mutable struct CounterFunc{F} - fun::F - counter::Int64 = 0 -end - -function (f::CounterFunc)(x...) - f.counter += 1 - return f.fun(x...) -end \ No newline at end of file diff --git a/classic_dp5/dp5.jl b/classic_dp5/dp5.jl deleted file mode 100644 index dd50eb8..0000000 --- a/classic_dp5/dp5.jl +++ /dev/null @@ -1,514 +0,0 @@ -# include("dp5_types.jl") -using Base.Iterators:repeated - -function dopri5( - n, - fcn, - x, # mutate - y, # mutate - xend, - rtol, - atol, - options, - work, -) - - arret = false - - ###### nmax - maximal number of steps - if options.maximum_allowed_steps < 0 - if options.print_error_messages - println("Maximum Allowed steps cannot be negative") - end - arret = true - end - - nmax = options.maximum_allowed_steps - - ###### nstiff - parameters for stiffness detection - nstiff = options.stiffness_test_activation_step - - if nstiff < 0 - nstiff = nmax + 10 - end - - ###### uround - smallest number satisfying 1.0 + uround > 1.0 - - uround = options.uround - if (uround <= 1e-35) || (uround >= 1.0) - if options.print_error_messages - println("WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:", work[1]) - end - arret = true - end - - ####### safety factor - safe = options.safety_factor - - if (safe >= 1.0) || (safe <= 1e-4) - if options.print_error_messages - println("CURIOUS INPUT FOR SAFETY FACTOR WORK[2]=", work[2]) - end - arret = true - end - - ###### fac1, fac2 - parameters for step size selection - fac1 = options.step_size_selection_one - fac2 = options.step_size_selection_two - - ###### beta - for step control stabilization - - beta = options.beta - if beta > 0.2 - if options.print_error_messages - println("CURIOUS INPUT FOR BETA: ", beta) - end - arret = true - end - - ####### maximal step size - - if work[6] == 0.0 - hmax = xend-x - else - hmax = options.maximal_step_size - end - - ####### initial step size - h = options.initial_step_size - - ####### prepare entry-points for arrays in work - iey1 = 1 - iek1 = iey1 + n - iek2 = iek1 + n - iek3 = iek2 + n - iek4 = iek3 + n - iek5 = iek4 + n - iek6 = iek5 + n - ieys = iek6 + n - - ####### total storage requirement - istore = ieys - if istore > length(work) - if options.print_error_messages - println("INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=", istore) - end - arret = true - end - - ###### When a fail has occurred, return with idid = -1 - if arret - return DP5Report(x, -1, 0, 0, 0, 0) - end - - # indices to work start at the starting locations but for a view we need - dp5_report = dopcor( - n, fcn, x, y, xend, hmax, h, rtol, atol, nmax, uround, nstiff, - safe, beta, fac1, fac2, - #view(work, iey1:iey1+n-1), - view(work, iek1:iek1+n-1), - view(work, iek2:iek2+n-1), - view(work, iek3:iek3+n-1), - view(work, iek4:iek4+n-1), - view(work, iek5:iek5+n-1), - view(work, iek6:iek6+n-1), - #view(work, ieys:ieys+n-1) - ) - - return dp5_report - -end - -# return idid, x, StepsEvalReport -function dopcor( - n, - fcn, - x, - y, - xend, - hmax, - h, - rtol, - atol, - nmax, - uround, - nstiff, - safe, - beta, - fac1, - fac2, - #y1, - k1, - k2, - k3, - k4, - k5, - k6, - #ysti -) - ##### Initializations - - coeffs = cdopri() - c2 = coeffs[1] - c3 = coeffs[2] - c4 = coeffs[3] - c5 = coeffs[4] - a21 = coeffs[5] - a31 = coeffs[6] - a32 = coeffs[7] - a41 = coeffs[8] - a42 = coeffs[9] - a43 = coeffs[10] - a51 = coeffs[11] - a52 = coeffs[12] - a53 = coeffs[13] - a54 = coeffs[14] - a61 = coeffs[15] - a62 = coeffs[16] - a63 = coeffs[17] - a64 = coeffs[18] - a65 = coeffs[19] - a71 = coeffs[20] - a73 = coeffs[21] - a74 = coeffs[22] - a75 = coeffs[23] - a76 = coeffs[24] - e1 = coeffs[25] - e3 = coeffs[26] - e4 = coeffs[27] - e5 = coeffs[28] - e6 = coeffs[29] - e7 = coeffs[30] - - facold = 1e-4 - expo1 = 0.20-beta*0.75 - facc1 = 1.0/fac1 - facc2 = 1.0/fac2 - posneg = sign(1.0, xend-x) - - ###### Initial Preparations - nfcn = 0 - nstep = 0 - naccpt = 0 - nrejct = 0 - - #atoli = atol[1] # works with integers or arrays - #rtoli = rtol[1] # works with integers or arrays - atol_iter = atol isa Number ? repeated(atol) : atol - rtol_iter = rtol isa Number ? repeated(rtol) : rtol - - last = false - hlamb = 0.0 - iasti = 0 - nonsti = 0 - fcn(n, x, y, k1) - hmax = abs(hmax) - iord = 5 - - if h == 0.0 - h = hinit(n, fcn, x, y, xend, posneg, k1, k2, k3, iord, hmax, atol, rtol) - end - - nfcn += 2 - reject = false - - - ###### Basic Integration Step - xph = 0.0 # need to put this here because FORTRAN seems to be looser with scoping - while true - if nstep > nmax - # GOTO 78 - println(" MORE THAN NMAX = ", nmax, " STEPS ARE NEEDED") - return DP5Report(x, -2, 0, 0, 0, 0) - end - - if (0.10 * abs(h)) <= abs(x)*uround - # GOTO 77 - println("STEP SIZE TOO SMALL, H = ", h) - return DP5Report(x, -3, 0, 0, 0, 0) - end - - if ((x+1.01*h-xend)*posneg) > 0.0 - h = xend-x - last = true - end - - nstep += 1 - - ####### First 6 stages, just set to equality and should work bc everything is vector (no need for loops) - # 22 - y1 = y + h * a21 * k1 - fcn(n, x + c2 * h, y1, k2) - # 23 - y1 = y + h * ( a31 * k1 + a32 * k2) - fcn(n, x + c3 * h, y1, k3) - # 24 - y1 = y + h * (a41 * k1 + a42 * k2 + a43 * k3) - fcn(n, x + c4 * h, y1, k4) - # 25 - y1 = y + h * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4) - fcn(n, x+c5*h, y1, k5) - # 26 - ysti = y + h * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5) - xph = x + h - fcn(n, xph, ysti, k6) - # 27 - y1 = y + h * (a71 * k1 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6) - fcn(n, xph, y1, k2) - # 28 - k4 = h * (e1 * k1 + e3 * k3 + e4 * k4 + e5 * k5 + e6 * k6 + e7 * k2) - - nfcn += 6 - - ###### Error Estimation - err = 0.0 - - #= - if itol == 0 # do some kind of reduction here - for i in range(1, n) - sk = atoli + rtoli*max(abs(y[i]), abs(y1[i])) - err += abs(k4[i]/sk)^2 - end - else - for i in range(1, n) - sk = atoli[i] + rtoli[i]*max(abs(y[i]), abs(y1[i])) - err += abs(k4[i]/sk)^2 - end - end - =# - - err = mapreduce(+, atol_iter, rtol_iter, k4, y, ysti) do atoli, rtoli, k4i, yi, ystii - sk = atoli + rtoli*max(abs(yi), abs(ystii)) - abs(k4i/sk)^2 - end - - err = sqrt(err/n) - - ###### Computation of hnew - fac11 = err^expo1 - ###### Lund-Stabilization - fac = fac11/(facold^beta) - ###### we require fac1 <= hnew/h <= fac2 - fac = max(facc2, min(facc1, fac/safe)) # facc1, facc2, fac must be Float64 - hnew = h/fac - if err <= 1.0 - ###### Step is accepted - facold = max(err, 1e-4) - naccpt += 1 - ###### Stiffness Detection - if (mod(naccpt, nstiff) == 0) || (iasti > 0) - stnum = 0.0 - stden = 0.0 - - #= - for i in range(1, n) - stnum += abs(k2[i]-k6[i])^2 # added "abs" per Phillip's advice - stden += abs(y1[i]-ysti[i])^2 - end - =# - - stnum, stden = mapreduce(.+, k2, k6, y1, ysti) do k2i, k6i, y1i, ystii - stnum = abs(k2i-k6i)^2 - stden = abs(y1i-ystii)^2 - stnum, stden - end - - if stden > 0.0 - hlamb = h*sqrt(stnum/stden) - end - - - if hlamb > 3.25 - iasti += 1 - if iasti == 15 - println("THE PROBLEM SEEMS TO BECOME STIFF AT X = ", x) - end - else - nonsti += 1 - if nonsti == 6 - iasti = 0 - end - end - end - - # 44 - #= - for i in range(1, n) - k1[i] = k2[i] - y[i] = y1[i] - end - =# - copyto!(k1, k2) - copyto!(y, y1) - - x = xph - - ###### Normal Exit - if last - h = hnew - return DP5Report( - x, - 1, - nfcn, - nstep, - naccpt, - nrejct - ) - end - - if(abs(hnew) > hmax) - hnew = posneg*hmax - end - - if reject - hnew = posneg*min(abs(hnew), abs(h)) - end - reject = false - else - ###### Step is rejected - hnew = h/min(facc1, fac11/safe) - reject = true - if naccpt > 1 - nrejct += 1 - end - last = false - - end - h = hnew - end - -end - -function hinit( - n, - fcn, - x, - y, - xend, - posneg, - f0, - f1, - y1, - iord, - hmax, - atol, - rtol, -) - #= - Compute a first guess for explicit euler as - h = 0.01 * norm (y0) / norm (f0) - the increment for explicit euler is small - compared to the solution - =# - dnf = 0.0 - dny = 0.0 - #atoli = atol[1] - #rtoli = rtol[1] - - # [dnf, dny] = mapreduce(+, atol, rtol, f0, y; init=[0.0, 0.0]) do (atoli, rtoli, f0i, yi) - # sk = atoli + rtoli*abs(yi) - # dnf = abs(f0i/sk)^2 - # dny = abs(yi/sk)^2 - # [dnf, dny] - # end - - - atol_iter = atol isa Number ? repeated(atol) : atol - rtol_iter = rtol isa Number ? repeated(rtol) : rtol - - dnf, dny = mapreduce(.+, atol_iter, rtol_iter, f0, y) do atoli, rtoli, f0i, yi - sk = atoli + rtoli*abs(yi) - dnf = abs(f0i/sk)^2 - dny = abs(yi/sk)^2 - dnf, dny - end - - - - # problem with comparing ComplexF64 with Float64 - # take abs of dnf - if (dnf <= 1.0e-10) || (dny <= 1.0e-10) - h = 1.0e-6 - else - h = 0.01*sqrt(dny/dnf) - end - h = min(h, hmax) - h = sign(h, posneg) - - ###### Perform an explicit step - #= - for i in range(1, n) - y1[i] = y[i] + h*f0[i] - end - =# - y1 = y + h*f0 - fcn(n, x+h, y1, f1) - ###### Estimate the second derivative of the solution - der2 = 0.0 - - der2 = mapreduce(+, atol_iter, rtol_iter, f1, f0, y) do atoli, rtoli, f1i, f0i, yi - sk = atoli + rtoli*abs(yi) - ((f1i-f0i)/sk)^2 - end - - der2 = sqrt(der2)/h - - ##### Step size is computed such that - ##### H**IORD * MAX ( NORM(F0), NORM(F1), DER2 ) = 0.01 - der12 = max(abs(der2), sqrt(dnf)) - if der12 <= 1e-15 - h1 = max(1.0e-6, abs(h)*1.0e-3) - else - h1 = (0.01/der12)^(1.0/iord) - end - h = min(100*abs(h), h1, hmax) - return sign(h, posneg) - -end - -function cdopri() - - return [ - 0.2, #C2 - 0.3, #C3 - 0.8, #C4 - 8.0/9.0, #C5 - 0.2, #A21 - 3.0/40.0, #A31 - 9.0/40.0, #A32 - 44.0/45.0, #A4 - -56.0/15.0, #A42 - 32.0/9.0, #A43 - 19372.0/6561.0, #A52 - -25360.0/2187.0, # - 64448.0/6561.0, - -212.0/729.0, - 9017.0/3168.0, - -355.0/33.0, - 46732.0/5247.0, - 49.0/176.0, - -5103.0/18656.0, - 35.0/384.0, - 500.0/1113.0, - 125.0/192.0, - -2187.0/6784.0, - 11.0/84.0, - 71.0/57600.0, - -71.0/16695.0, - 71.0/1920.0, - -17253.0/339200.0, - 22.0/525.0, - -1.0/40.0 - ] -end - -# FORTRAN sign - returns the value of A with the sign of B -function sign(a,b) - if b >= 0 - sign = 1.0 - else - sign = -1.0 - end - - return a*sign -end \ No newline at end of file diff --git a/classic_dp5/test_dp5.jl b/classic_dp5/test_dp5.jl deleted file mode 100644 index 657fedd..0000000 --- a/classic_dp5/test_dp5.jl +++ /dev/null @@ -1,24 +0,0 @@ -include("dp5.jl") -include("dp5_types.jl") - -# FCN must follow format FCN(N, X, Y, F, RPAR, IPAR) - -function fcn(n, x, y, f) - f[1] = ℯ^x - f[2] = ℯ^x -end - -y = [1.0, 1.0] -work = zeros(16) - -dopri5( - 2, - fcn, - 0.0, - y, - 1.0, - 1e-10, - 1e-10, - DP5Options(), - work,# work, should be 8*N -) \ No newline at end of file diff --git a/classic_dp5/test_schrodinger_eq.jl b/classic_dp5/test_schrodinger_eq.jl deleted file mode 100644 index 8781438..0000000 --- a/classic_dp5/test_schrodinger_eq.jl +++ /dev/null @@ -1,41 +0,0 @@ -using Bloqade -using BloqadeExpr: Hamiltonian -include("dp5.jl") -include("dp5_types.jl") - -nsites = 1; -atoms = generate_sites(ChainLattice(), nsites, scale = 5.74) -h = rydberg_h(atoms; Ω = 11 * 2π, Δ = 0) -reg = zero_state(1) -# stop short of creating the SchrodingerProblem, we want the SchrodignerEquation -# which we can than trivially wrap with the fcn(n,x,y,f) that DP5 accepts - -# Generate SchrodingerEquation -state = statevec(reg) -space = YaoSubspaceArrayReg.space(reg) -T = real(eltype(state)) -T = isreal(h) ? T : Complex{T} -eq = SchrodingerEquation(h, Hamiltonian(T, h, space)) -# invoke eq via eq(dstate, state, p, t::Number) - -function fcn(n, x, y, f) - eq(f, y, nothing, x) -end - -# y = [zeros(ComplexF64, 2)] # initial y_0 at x = 0`` -y = statevec(reg) -# work = [zeros(ComplexF64, 2) for i in range(1,8)] -work = zeros(ComplexF64, 16) - -dopri5( - 2, - fcn, - 0.0, - y, - 1.6, - 1e-10, - 1e-10, - DP5Options(), - work,# work, should be 8*N -) - diff --git a/classic_dp5/test_stateful_dp5.jl b/classic_dp5/test_stateful_dp5.jl deleted file mode 100644 index 821bc91..0000000 --- a/classic_dp5/test_stateful_dp5.jl +++ /dev/null @@ -1,38 +0,0 @@ -include("dp5.jl") - - -# differential equation to solve for -function fcn(n, x, y, f) - f[1] = y[1] - x^2 + 1 -end - -# initial y -y = [0.5] - -# intermediate times to save -times = [1.0, 2.0, 3.0, 4.0] -# place to save values -values = [] -# initial x -x = 0.0 - -for t in times - xend = t - - dopri5( - 1, - fcn, - x, - y, - xend, - 1e-10, - 1e-10, - 0, - DP5Options(), - zeros(8),# work, should be 8*N + 5*NRDENS+21 = 8*2 + 5*0 + 21 - ) - - push!(values, y[1]) - global x = xend # set next initial x to last x - -end \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..9245e0d --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,7 @@ +[deps] +DocThemeIndigo = "8bac0ac5-51bf-41f9-885e-2bf1ac2bec5f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DormandPrince = "5e45e72d-22b8-4dd0-9c8b-f96714864bcd" + +[compat] +DocThemeIndigo = "0.1.1" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..023ce5f --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,21 @@ +using Documenter +using DormandPrince +using DocThemeIndigo + +indigo = DocThemeIndigo.install(Configurations) + +makedocs(; + modules = [DormandPrince], + format = Documenter.HTML( + prettyurls = !("local" in ARGS), + canonical="https://John Long.github.io/DormandPrince.jl", + assets=String[indigo], + ), + pages = [ + "Home" => "index.md", + ], + repo = "https://github.com/John Long/DormandPrince.jl", + sitename = "DormandPrince.jl", +) + +deploydocs(; repo = "https://github.com/John Long/DormandPrince.jl") diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..eaac378 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,14 @@ +```@meta +CurrentModule = DormandPrince +``` + +# DormandPrince + +Documentation for [DormandPrince](https://github.com/John Long/DormandPrince.jl). + +```@index +``` + +```@autodocs +Modules = [DormandPrince] +``` diff --git a/integrate_test.jl b/integrate_test.jl deleted file mode 100644 index a82e7c6..0000000 --- a/integrate_test.jl +++ /dev/null @@ -1,45 +0,0 @@ -include("integrate.jl") - -using Bloqade -using BloqadeExpr: Hamiltonian - -nsites = 1; -atoms = generate_sites(ChainLattice(), nsites, scale = 5.74) -h = rydberg_h(atoms; Ω = 11 * 2π, Δ = 0) -reg = zero_state(1) - -state = statevec(reg) -space = YaoSubspaceArrayReg.space(reg) -T = real(eltype(state)) -T = isreal(h) ? T : Complex{T} -eq = SchrodingerEquation(h, Hamiltonian(T, h, space)) - -function fcn(x, y, f) - eq(f, y, nothing, x) -end - -y = statevec(reg) - -solver = DP5Solver( - fcn, - 0.0, - y -) - -iterator = integrate(solver, [1.8, 1.9, 2.0]) - - -for (t, y) in iterator - println(t) - println(y) -end - - -#= -function simple_callback(time, state) - println("State ", state, " at time ", time) - return copy(state) -end - -integrate(simple_callback, solver, [1.8, 1.9, 2.0]) -=# \ No newline at end of file diff --git a/schrodinger_eq_reference.jl b/schrodinger_eq_reference.jl deleted file mode 100644 index e7de9d2..0000000 --- a/schrodinger_eq_reference.jl +++ /dev/null @@ -1,10 +0,0 @@ -using Bloqade - -atoms = generate_sites(ChainLattice(), 1, scale = 5.74) -h = rydberg_h(atoms; Ω = 20 * 2π, Δ = 0) -reg = zero_state(1) - -prob = SchrodingerProblem(reg, 1.6, h; algo=DP5()) -emulate!(prob) - -statevec(prob.reg) \ No newline at end of file diff --git a/solver_test.jl b/solver_test.jl deleted file mode 100644 index d18562c..0000000 --- a/solver_test.jl +++ /dev/null @@ -1,44 +0,0 @@ -include("solver.jl") -include("types.jl") -using Bloqade -using BloqadeExpr: Hamiltonian -using JET - -nsites = 1; -atoms = generate_sites(ChainLattice(), nsites, scale = 5.74) -h = rydberg_h(atoms; Ω = 20 * 2π, Δ = 0) -reg = zero_state(1) -# stop short of creating the SchrodingerProblem, we want the SchrodignerEquation -# which we can than trivially wrap with the fcn(n,x,y,f) that DP5 accepts - -# Generate SchrodingerEquation -state = statevec(reg) -space = YaoSubspaceArrayReg.space(reg) -T = real(eltype(state)) -T = isreal(h) ? T : Complex{T} -eq = SchrodingerEquation(h, Hamiltonian(T, h, space)) -# invoke eq via eq(dstate, state, p, t::Number) - -function fcn(x, y, f) - #println(y) - #println(x) - #println(f) - eq(f, y, nothing, x) -end - -y = statevec(reg) - -solver = DP5Solver( - fcn, - 0.0, - y -) - -#@report_opt dopri5(solver, 1.6) -# @code_warntype dopri5(solver, 1.6) -dopri5(solver, 1.6) - -# narrow down to Dopcor -## dopcor(solver, xend, hmax, h) -## @report_opt dopcor(solver, 1.6, 1.6, 0.0) - diff --git a/src/DormandPrince.jl b/src/DormandPrince.jl new file mode 100644 index 0000000..832347f --- /dev/null +++ b/src/DormandPrince.jl @@ -0,0 +1,10 @@ +module DormandPrince +using Base.Iterators:repeated, Repeated + +include("types.jl") +include("solver.jl") +include("integrate.jl") +include("checks.jl") +include("helpers.jl") + +end # DormandPrince diff --git a/checks.jl b/src/checks.jl similarity index 97% rename from checks.jl rename to src/checks.jl index ef87401..95465de 100644 --- a/checks.jl +++ b/src/checks.jl @@ -1,4 +1,4 @@ -include("types.jl") +# include("types.jl") # make enums for every error type and return that instead of printing errors function check_max_allowed_steps(options::DP5Options{T}) where T diff --git a/helpers.jl b/src/helpers.jl similarity index 88% rename from helpers.jl rename to src/helpers.jl index 47937cf..73a76cd 100644 --- a/helpers.jl +++ b/src/helpers.jl @@ -1,14 +1,3 @@ -# FORTRAN sign - returns the value of A with the sign of B -function sign(a,b) - if b >= 0 - sign = 1.0 - else - sign = -1.0 - end - - return a*sign -end - function do_step!(solver, h) # define constants @@ -43,27 +32,27 @@ function do_step!(solver, h) e6=22.0/525.0 e7=-1.0/40.0 - ####### First 6 stages, just set to equality and should work bc everything is vector (no need for loops) - # 22 + ####### First 6 stages + solver.y1 .= solver.y .+ h .* a21 .* solver.k1 solver.f(solver.vars.x + c2 * h, solver.y1, solver.k2) - # 23 + solver.y1 .= solver.y .+ h .* (a31 .* solver.k1 .+ a32 .* solver.k2) solver.f(solver.vars.x + c3 * h, solver.y1, solver.k3) - # 24 + solver.y1 .= solver.y .+ h .* (a41 .* solver.k1 .+ a42 .* solver.k2 .+ a43 .* solver.k3) solver.f(solver.vars.x + c4 * h, solver.y1, solver.k4) - # 25 + solver.y1 .= solver.y .+ h .* (a51 .* solver.k1 .+ a52 .* solver.k2 .+ a53 .* solver.k3 .+ a54 .* solver.k4) solver.f(solver.vars.x + c5*h, solver.y1, solver.k5) - # 26 + solver.ysti .= solver.y .+ h .* (a61 .* solver.k1 .+ a62 .* solver.k2 .+ a63 .* solver.k3 .+ a64 .* solver.k4 .+ a65 .* solver.k5) xph = solver.vars.x + h solver.f(xph, solver.ysti, solver.k6) - # 27 + solver.y1 .= solver.y .+ h .* (a71 .* solver.k1 .+ a73 .* solver.k3 .+ a74 .* solver.k4 .+ a75 .* solver.k5 .+ a76 .* solver.k6) solver.f(xph, solver.y1, solver.k2) - # 28 + solver.k4 .= h .* (e1 .* solver.k1 .+ e3 .* solver.k3 .+ e4 .* solver.k4 .+ e5 .* solver.k5 .+ e6 .* solver.k6 .+ e7 .* solver.k2) end @@ -115,7 +104,6 @@ function stiffness_detection!(solver, naccpt, h) if solver.vars.hlamb > 3.25 solver.vars.iasti += 1 if solver.vars.iasti == 15 - # turn this into a debug statement @debug "The problem seems to become stiff at $x" end else @@ -134,9 +122,7 @@ function euler_first_guess(solver, hmax, posneg) abs(f0i/sk)^2, abs(yi/sk)^2 # dnf, dny end - - # problem with comparing ComplexF64 with Float64 - # take abs of dnf + if (dnf <= 1.0e-10) || (dny <= 1.0e-10) h = 1.0e-6 else diff --git a/integrate.jl b/src/integrate.jl similarity index 66% rename from integrate.jl rename to src/integrate.jl index dd6e468..3066b8c 100644 --- a/integrate.jl +++ b/src/integrate.jl @@ -1,47 +1,3 @@ -#= -# low level interface: - -solver = DP5Solver(fcn,0,y0) -# low level interface -integrate(solver, 0.1) -# do stuff here -integrate(solver, 0.2) -# do stuff here -integrate(solver, 0.3) - -# iterator interface -# TODO: Check how to dispatch to Iterator of a specific type objects -iterator = integrate(solver, times) - -for t,y in iterator: - #do stuff -end - -# callback for each time in times -integrate(callback, solver, times) -> Vector - -results = integrate(solver, times) do t, state - # do stuff with t and state -end - - - -# 99% of use cases: -function integrate(callback, solver::DP5Solver{StateVec, T}, times::AbstractVector{T}; sort_times::Bool = true) where {StateVec, T} - times = sort_times ? sorted(collect(times)) : times - - result = [] - for time in times - integrate(solver, time) - push!(result, callback(time, solver.y)) - end -end -=# - -# integrate(solver, times) should return an iterator - -include("solver.jl") - mutable struct DP5Iterator{T <: Real} solver::DP5Solver times::AbstractVector{T} @@ -81,9 +37,10 @@ function Base.iterate(dp5_iterator::DP5Iterator, state) end # 3 modes of operation for integrate -# 1. integrate(solver, time) -> state +# 1. integrate(solver, time) -> state (modify solver object in place) # 2. integrate(solver, times) -> iterator # 3. integrate(callback, solver, times) -> vector of states with callback applied + integrate(solver::DP5Solver, time::Real) = dopri5(solver, time) integrate(solver::DP5Solver, times::AbstractVector{T}) where {T <: Real} = DP5Iterator(solver, times) diff --git a/solver.jl b/src/solver.jl similarity index 95% rename from solver.jl rename to src/solver.jl index a3b8bd0..5540843 100644 --- a/solver.jl +++ b/src/solver.jl @@ -1,7 +1,7 @@ -using Base.Iterators:repeated -include("types.jl") -include("checks.jl") -include("helpers.jl") +# using Base.Iterators:repeated +#include("types.jl") +#include("checks.jl") +#include("helpers.jl") function dopri5( solver, @@ -57,6 +57,9 @@ function dopri5( # update with final h solver.vars.h = h + # reset the necessary vars + solver.vars.last = false + return dp5_report end @@ -195,7 +198,8 @@ function hinit( ###### Perform an explicit step #y1 = y + h*f0 #fcn(n, x+h, y1, f1) - copyto!(solver.y1, solver.y + h*solver.k1) + # copyto!(solver.y1, solver.y + h*solver.k1) + solver.y1 .= solver.y .+ h .*solver.k1 solver.f(solver.vars.x + h, solver.k3, solver.k2) ###### Estimate the second derivative of the solution diff --git a/types.jl b/src/types.jl similarity index 91% rename from types.jl rename to src/types.jl index 7222b6c..1b25f39 100644 --- a/types.jl +++ b/src/types.jl @@ -1,4 +1,4 @@ -using Base.Iterators:repeated, Repeated +# using Base.Iterators:repeated, Repeated @enum Idid begin COMPUTATION_SUCCESSFUL = 1 @@ -20,10 +20,10 @@ struct DP5Report{T <: Real} checks::Checks idid::Idid - num_func_evals::Int64 - num_computed_steps::Int64 - num_accepted_steps::Int64 - num_rejected_steps::Int64 + num_func_evals::Int + num_computed_steps::Int + num_accepted_steps::Int + num_rejected_steps::Int end @kwdef struct DP5Options{T <: Real} @@ -37,9 +37,9 @@ end initial_step_size::T = 0.0 # originally in iwork[1] - iwork[4] - maximum_allowed_steps::Int64 = 100000 + maximum_allowed_steps::Int = 100000 print_error_messages::Bool = true - stiffness_test_activation_step::Int64 = 1000 + stiffness_test_activation_step::Int = 1000 # should be either vector or repeated for type atol::Union{T, Vector{T}} = 1e-10 @@ -69,8 +69,8 @@ end x::T = zero(T) h::T = zero(T) facold::T = 1e-4 - iasti::Int64 = 0 - nonsti::Int64 = 0 + iasti::Int = 0 + nonsti::Int = 0 hlamb::T = zero(T) last::Bool = false end diff --git a/test/Manifest.toml b/test/Manifest.toml new file mode 100644 index 0000000..99e7320 --- /dev/null +++ b/test/Manifest.toml @@ -0,0 +1,59 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.9.3" +manifest_format = "2.0" +project_hash = "70c6548fc0267b7c924ca6e56c4af9fd2abca604" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.5+0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.21+4" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.8.0+0" diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..a2cd2f8 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,3 @@ +[deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/exact_evol_helpers.jl b/test/exact_evol_helpers.jl new file mode 100644 index 0000000..c92f5a5 --- /dev/null +++ b/test/exact_evol_helpers.jl @@ -0,0 +1,16 @@ +function evolution_operator(t::Float64) + ϕ = 2.2 * sin(π * t)^2 + U = zeros(ComplexF64, 2,2) + U[1,1] = 1 / sqrt(2) + U[2,1] = 1 / sqrt(2) + U[2,2] = 1 / sqrt(2) + U[1,2] = -1 / sqrt(2) + + U * diagm(exp.([-im*ϕ, im*ϕ])) * U' +end + +function solution(t) + U = evolution_operator(t) + return U * [1.0, 0.0] + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..0c6dda0 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,85 @@ +using Test +using LinearAlgebra +using DormandPrince: DP5Solver, dopri5, integrate + +include("exact_evol_helpers.jl") + + +function fcn(x, y, f) + g(x) = 2.2*2π*sin(2π*x) + + f[1] = -1im * g(x)/2 * y[2] + f[2] = -1im * g(x)/2 * y[1] +end + +# standalone solver test +@testset "Integration Test" begin + solver = DP5Solver( + fcn, + 0.0, + ComplexF64[1.0, 0.0] + ) + + dopri5(solver, 2π) + + @test solver.y ≈ solution(2π) +end + +# Test integrate() +@testset "Integration Interface Test" begin + + # test iterator generation + @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)) + 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 + + @test dp5_values ≈ exact_values + end +end + +