From 8d9172f9f164e1d0eed68ca04db2bf21c623f120 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 06:43:23 +0530 Subject: [PATCH 01/28] Added OrdinaryDiffEqSymplecticRK --- lib/OrdinaryDiffEqSymplecticRK/Project.toml | 4 + .../src/OrdinaryDiffEqSymplecticRK.jl | 16 ++ .../src/alg_utils.jl | 17 ++ .../src/algorithms.jl | 224 ++++++++++++++++++ .../src}/symplectic_caches.jl | 2 +- .../src}/symplectic_perform_step.jl | 2 +- .../src}/symplectic_tableaus.jl | 2 +- .../test/runtests.jl | 3 + .../test}/symplectic_tests.jl | 2 +- src/OrdinaryDiffEq.jl | 13 +- src/alg_utils.jl | 18 -- src/algorithms.jl | 212 ----------------- test/runtests.jl | 5 +- 13 files changed, 278 insertions(+), 242 deletions(-) create mode 100644 lib/OrdinaryDiffEqSymplecticRK/Project.toml create mode 100644 lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl create mode 100644 lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl create mode 100644 lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl rename {src/caches => lib/OrdinaryDiffEqSymplecticRK/src}/symplectic_caches.jl (99%) rename {src/perform_step => lib/OrdinaryDiffEqSymplecticRK/src}/symplectic_perform_step.jl (99%) rename {src/tableaus => lib/OrdinaryDiffEqSymplecticRK/src}/symplectic_tableaus.jl (99%) create mode 100644 lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl rename {test/algconvergence => lib/OrdinaryDiffEqSymplecticRK/test}/symplectic_tests.jl (99%) diff --git a/lib/OrdinaryDiffEqSymplecticRK/Project.toml b/lib/OrdinaryDiffEqSymplecticRK/Project.toml new file mode 100644 index 0000000000..3b304389da --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/Project.toml @@ -0,0 +1,4 @@ +name = "OrdinaryDiffEqSymplecticRK" +uuid = "fa646aed-7ef9-47eb-84c4-9443fc8cbfa8" +authors = ["ParamThakkar123 "] +version = "0.1.0" diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl new file mode 100644 index 0000000000..4bc650a671 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -0,0 +1,16 @@ +module OrdinaryDiffEqSymplecticRK + +using OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache + + +include("algorithms.jl") +include("alg_utils.jl") +include("symplectic_caches.jl") +include("symplectic_perform_step.jl") +include("symplectic_tableaus.jl") + +export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, + McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, + CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 + +end diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl b/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl new file mode 100644 index 0000000000..ee6e0bee37 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl @@ -0,0 +1,17 @@ +alg_order(alg::SymplecticEuler) = 1 +alg_order(alg::VelocityVerlet) = 2 +alg_order(alg::VerletLeapfrog) = 2 +alg_order(alg::PseudoVerletLeapfrog) = 2 +alg_order(alg::McAte2) = 2 +alg_order(alg::Ruth3) = 3 +alg_order(alg::McAte3) = 3 +alg_order(alg::McAte4) = 4 +alg_order(alg::CandyRoz4) = 4 +alg_order(alg::CalvoSanz4) = 4 +alg_order(alg::McAte42) = 4 +alg_order(alg::McAte5) = 5 +alg_order(alg::Yoshida6) = 6 +alg_order(alg::KahanLi6) = 6 +alg_order(alg::McAte8) = 8 +alg_order(alg::KahanLi8) = 8 +alg_order(alg::SofSpa10) = 10 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl b/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl new file mode 100644 index 0000000000..22f0f9258b --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl @@ -0,0 +1,224 @@ +struct SymplecticEuler <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct VelocityVerlet <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct VerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct PseudoVerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte2 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{ruth1983canonical, +title={A canonical integration technique}, +author={Ruth, Ronald D}, +journal={IEEE Trans. Nucl. Sci.}, +volume={30}, +number={CERN-LEP-TH-83-14}, +pages={2669--2671}, +year={1983} +} +""" + +struct Ruth3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{candy1991symplectic, +title={A symplectic integration algorithm for separable Hamiltonian functions}, +author={Candy, J and Rozmus, W}, +journal={Journal of Computational Physics}, +volume={92}, +number={1}, +pages={230--256}, +year={1991}, +publisher={Elsevier} +} +""" + +struct CandyRoz4 <: OrdinaryDiffEqPartitionedAlgorithm end + +struct McAte4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{sanz1993symplectic, +title={Symplectic numerical methods for Hamiltonian problems}, +author={Sanz-Serna, Jes{\'u}s Maria and Calvo, Mari-Paz}, +journal={International Journal of Modern Physics C}, +volume={4}, +number={02}, +pages={385--392}, +year={1993}, +publisher={World Scientific} +} +""" + +struct CalvoSanz4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte42 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte5 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{yoshida1990construction, +title={Construction of higher order symplectic integrators}, +author={Yoshida, Haruo}, +journal={Physics letters A}, +volume={150}, +number={5-7}, +pages={262--268}, +year={1990}, +publisher={Elsevier} +} +""" + +struct Yoshida6 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{kahan1997composition, +title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, +author={Kahan, William and Li, Ren-Cang}, +journal={Mathematics of computation}, +volume={66}, +number={219}, +pages={1089--1099}, +year={1997} +} +""" + +struct KahanLi6 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1995numerical, +title={On the numerical integration of ordinary differential equations by symmetric composition methods}, +author={McLachlan, Robert I}, +journal={SIAM Journal on Scientific Computing}, +volume={16}, +number={1}, +pages={151--168}, +year={1995}, +publisher={SIAM} +} +""" + +struct McAte8 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{kahan1997composition, +title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, +author={Kahan, William and Li, Ren-Cang}, +journal={Mathematics of computation}, +volume={66}, +number={219}, +pages={1089--1099}, +year={1997} +} +""" + +struct KahanLi8 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{sofroniou2005derivation, +title={Derivation of symmetric composition constants for symmetric integrators}, +author={Sofroniou, Mark and Spaletta, Giulia}, +journal={Optimization Methods and Software}, +volume={20}, +number={4-5}, +pages={597--613}, +year={2005}, +publisher={Taylor \\& Francis} +} +""" +struct SofSpa10 <: OrdinaryDiffEqPartitionedAlgorithm end \ No newline at end of file diff --git a/src/caches/symplectic_caches.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl similarity index 99% rename from src/caches/symplectic_caches.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl index 71c55340dc..483fa6996b 100644 --- a/src/caches/symplectic_caches.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl @@ -416,4 +416,4 @@ function alg_cache(alg::SofSpa10, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} SofSpa10ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) -end +end \ No newline at end of file diff --git a/src/perform_step/symplectic_perform_step.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl similarity index 99% rename from src/perform_step/symplectic_perform_step.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl index 413f1882a2..f309e50098 100644 --- a/src/perform_step/symplectic_perform_step.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl @@ -2022,4 +2022,4 @@ end integrator.stats.nf += 36 integrator.stats.nf2 += 36 store_symp_state!(integrator, cache, kdu, ku) -end +end \ No newline at end of file diff --git a/src/tableaus/symplectic_tableaus.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl similarity index 99% rename from src/tableaus/symplectic_tableaus.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl index 98676979c3..c7bee07b0c 100644 --- a/src/tableaus/symplectic_tableaus.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl @@ -804,4 +804,4 @@ function SofSpa10ConstantCache(T::Type, T2::Type) b19, b20, b21, b22, b23, b24, b25, b26, b27, b28, b29, b30, b31, b32, b33, b34, b35, b36) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl b/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl new file mode 100644 index 0000000000..231d0ee343 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl @@ -0,0 +1,3 @@ +using SafeTestsets + +@time @safetestset "Synplectic Tests" include("symplectic_tests.jl") \ No newline at end of file diff --git a/test/algconvergence/symplectic_tests.jl b/lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl similarity index 99% rename from test/algconvergence/symplectic_tests.jl rename to lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl index 489f86a9c5..d47e6312d9 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl @@ -99,4 +99,4 @@ ref = solve( else @test maximum(ref[4, :] - sol[4, :]) < 3e-3 end -end +end \ No newline at end of file diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index a688dacd6a..86b4fafc04 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -154,7 +154,6 @@ include("caches/firk_caches.jl") include("caches/kencarp_kvaerno_caches.jl") include("caches/linear_caches.jl") include("caches/linear_nonlinear_caches.jl") -include("caches/symplectic_caches.jl") include("caches/rosenbrock_caches.jl") include("caches/rkn_caches.jl") include("caches/adams_bashforth_moulton_caches.jl") @@ -167,7 +166,6 @@ include("caches/qprk_caches.jl") include("tableaus/low_order_rk_tableaus.jl") include("tableaus/high_order_rk_tableaus.jl") -include("tableaus/symplectic_tableaus.jl") include("tableaus/verner_tableaus.jl") include("tableaus/rosenbrock_tableaus.jl") include("tableaus/sdirk_tableaus.jl") @@ -184,7 +182,6 @@ include("initialize_dae.jl") include("wrappers.jl") include("perform_step/fixed_timestep_perform_step.jl") -include("perform_step/symplectic_perform_step.jl") include("perform_step/rkn_perform_step.jl") include("perform_step/split_perform_step.jl") include("perform_step/linear_perform_step.jl") @@ -266,6 +263,12 @@ include("../lib/OrdinaryDiffEqFeagin/src/OrdinaryDiffEqFeagin.jl") using ..OrdinaryDiffEqFeagin export Feagin10, Feagin12, Feagin14 +include("../lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl") +using ..OrdinaryDiffEqSymplecticRK +export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, + McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, + CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 + import PrecompileTools PrecompileTools.@compile_workload begin @@ -426,10 +429,6 @@ export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3B, EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43 -export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, - McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, - CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 - export SHLDDRK52, SHLDDRK_2N export SplitEuler diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4467a10eb1..878b6d9f32 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -421,24 +421,6 @@ alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 alg_order(alg::FRK65) = 6 -alg_order(alg::SymplecticEuler) = 1 -alg_order(alg::VelocityVerlet) = 2 -alg_order(alg::VerletLeapfrog) = 2 -alg_order(alg::PseudoVerletLeapfrog) = 2 -alg_order(alg::McAte2) = 2 -alg_order(alg::Ruth3) = 3 -alg_order(alg::McAte3) = 3 -alg_order(alg::McAte4) = 4 -alg_order(alg::CandyRoz4) = 4 -alg_order(alg::CalvoSanz4) = 4 -alg_order(alg::McAte42) = 4 -alg_order(alg::McAte5) = 5 -alg_order(alg::Yoshida6) = 6 -alg_order(alg::KahanLi6) = 6 -alg_order(alg::McAte8) = 8 -alg_order(alg::KahanLi8) = 8 -alg_order(alg::SofSpa10) = 10 - alg_order(alg::IRKN3) = 3 alg_order(alg::Nystrom4) = 4 alg_order(alg::FineRKN4) = 4 diff --git a/src/algorithms.jl b/src/algorithms.jl index 973ebb43d4..4152a7057b 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -307,218 +307,6 @@ end ################################################################################ -# Symplectic methods - -struct SymplecticEuler <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct VelocityVerlet <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct VerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct PseudoVerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte2 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{ruth1983canonical, -title={A canonical integration technique}, -author={Ruth, Ronald D}, -journal={IEEE Trans. Nucl. Sci.}, -volume={30}, -number={CERN-LEP-TH-83-14}, -pages={2669--2671}, -year={1983} -} -""" -struct Ruth3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{candy1991symplectic, -title={A symplectic integration algorithm for separable Hamiltonian functions}, -author={Candy, J and Rozmus, W}, -journal={Journal of Computational Physics}, -volume={92}, -number={1}, -pages={230--256}, -year={1991}, -publisher={Elsevier} -} -""" -struct CandyRoz4 <: OrdinaryDiffEqPartitionedAlgorithm end -struct McAte4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{sanz1993symplectic, -title={Symplectic numerical methods for Hamiltonian problems}, -author={Sanz-Serna, Jes{\'u}s Maria and Calvo, Mari-Paz}, -journal={International Journal of Modern Physics C}, -volume={4}, -number={02}, -pages={385--392}, -year={1993}, -publisher={World Scientific} -} -""" -struct CalvoSanz4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte42 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte5 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{yoshida1990construction, -title={Construction of higher order symplectic integrators}, -author={Yoshida, Haruo}, -journal={Physics letters A}, -volume={150}, -number={5-7}, -pages={262--268}, -year={1990}, -publisher={Elsevier} -} -""" -struct Yoshida6 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{kahan1997composition, -title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, -author={Kahan, William and Li, Ren-Cang}, -journal={Mathematics of computation}, -volume={66}, -number={219}, -pages={1089--1099}, -year={1997} -} -""" -struct KahanLi6 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1995numerical, -title={On the numerical integration of ordinary differential equations by symmetric composition methods}, -author={McLachlan, Robert I}, -journal={SIAM Journal on Scientific Computing}, -volume={16}, -number={1}, -pages={151--168}, -year={1995}, -publisher={SIAM} -} -""" -struct McAte8 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{kahan1997composition, -title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, -author={Kahan, William and Li, Ren-Cang}, -journal={Mathematics of computation}, -volume={66}, -number={219}, -pages={1089--1099}, -year={1997} -} -""" -struct KahanLi8 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{sofroniou2005derivation, -title={Derivation of symmetric composition constants for symmetric integrators}, -author={Sofroniou, Mark and Spaletta, Giulia}, -journal={Optimization Methods and Software}, -volume={20}, -number={4-5}, -pages={597--613}, -year={2005}, -publisher={Taylor \\& Francis} -} -""" -struct SofSpa10 <: OrdinaryDiffEqPartitionedAlgorithm end - # Nyström methods """ diff --git a/test/runtests.jl b/test/runtests.jl index 14aa5a0703..54b6bddff6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -199,10 +199,13 @@ end @time @safetestset "FIRK Tests" include("algconvergence/ode_firk_tests.jl") @time @safetestset "Linear-Nonlinear Methods Tests" include("algconvergence/linear_nonlinear_convergence_tests.jl") @time @safetestset "Linear-Nonlinear Krylov Methods Tests" include("algconvergence/linear_nonlinear_krylov_tests.jl") - @time @safetestset "Symplectic Tests" include("algconvergence/symplectic_tests.jl") @time @safetestset "Quadruple precision Runge-Kutta Tests" include("algconvergence/ode_quadruple_precision_tests.jl") end + if !is_APPVEYOR && GROUP == "Symplectic" + @time @safetestset "Symplectic Tests" include("../lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl") + end + if !is_APPVEYOR && GROUP == "Extrapolation" @time @safetestset "Extrapolation Tests" include("../lib/OrdinaryDiffEqExtrapolation/test/runtests.jl") end From 6422d1b0c50a7c8ab7e19b0cfa12d941a1d45c5c Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 06:47:30 +0530 Subject: [PATCH 02/28] Cache --- .../src/OrdinaryDiffEqSymplecticRK.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index 4bc650a671..50c526934b 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -1,6 +1,7 @@ module OrdinaryDiffEqSymplecticRK -using OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache +import OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + @cache include("algorithms.jl") From 229e29cd8a4da052562922769fb81e0f217254cb Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 06:52:27 +0530 Subject: [PATCH 03/28] muladd --- .../src/OrdinaryDiffEqSymplecticRK.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index 50c526934b..c97d964922 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -1,7 +1,8 @@ module OrdinaryDiffEqSymplecticRK import OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @cache + @cache +import MuladdMacro: @muladd include("algorithms.jl") From f7ff2dabd900269891d351c2dc3d5f8f4eaa5e38 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 06:55:04 +0530 Subject: [PATCH 04/28] Unpack --- .../src/OrdinaryDiffEqSymplecticRK.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index c97d964922..cae0ae65a5 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -1,7 +1,7 @@ module OrdinaryDiffEqSymplecticRK import OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @cache + @cache, @unpack import MuladdMacro: @muladd From 10534f88001909e75c07608d8d6775038678a0c4 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 06:58:57 +0530 Subject: [PATCH 05/28] Changes --- lib/OrdinaryDiffEqSymplecticRK/Project.toml | 17 ++++++++++++++++- .../src/OrdinaryDiffEqSymplecticRK.jl | 14 +++++++++++--- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/Project.toml b/lib/OrdinaryDiffEqSymplecticRK/Project.toml index 3b304389da..d25f5cfd73 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/Project.toml +++ b/lib/OrdinaryDiffEqSymplecticRK/Project.toml @@ -1,4 +1,19 @@ name = "OrdinaryDiffEqSymplecticRK" uuid = "fa646aed-7ef9-47eb-84c4-9443fc8cbfa8" authors = ["ParamThakkar123 "] -version = "0.1.0" +version = "1.0.0" + +[deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +julia = "1.10" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index cae0ae65a5..ee8d74bfa5 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -1,8 +1,16 @@ module OrdinaryDiffEqSymplecticRK -import OrdinaryDiffEq: OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @cache, @unpack -import MuladdMacro: @muladd +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, + explicit_rk_docstring, trivial_limiter!, + _ode_interpolant!, _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools include("algorithms.jl") From 0149a3a244b00e563f99c7d83f1ca6f43153550e Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 07:01:44 +0530 Subject: [PATCH 06/28] OrdinaryDiffEqPartitionedAlgorithm --- .../src/OrdinaryDiffEqSymplecticRK.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index ee8d74bfa5..473144bfc4 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -5,7 +5,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, calculate_residuals, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqPartitionedAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, du_alias_or_new, explicit_rk_docstring, trivial_limiter!, From 9b849e5a82b4ea1bd22cbd5a7ae95769a7708a7f Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Mon, 1 Jul 2024 07:06:17 +0530 Subject: [PATCH 07/28] Symplectic2ConstantCache --- .../src/OrdinaryDiffEqSymplecticRK.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl index 473144bfc4..7fc1144e47 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -16,8 +16,8 @@ using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools include("algorithms.jl") include("alg_utils.jl") include("symplectic_caches.jl") -include("symplectic_perform_step.jl") include("symplectic_tableaus.jl") +include("symplectic_perform_step.jl") export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, From a53466e71fc201eab1385a0dbcfdc56465ad1d6a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 1 Jul 2024 09:36:39 -0400 Subject: [PATCH 08/28] Update CI.yml --- .github/workflows/CI.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index bc7b5a436d..97989c7464 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,6 +22,7 @@ jobs: - AlgConvergence_II - AlgConvergence_III - Downstream + - Symplectic - Extrapolation - StabilizedRK - StabilizedIRK From 217e5da75c12a6e8ce11eb97d8e2cc2d3e0c7d74 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 11:25:52 +0530 Subject: [PATCH 09/28] RKN Methods --- lib/OrdinaryDiffEqRKN/Project.toml | 23 ++ .../src/OrdinaryDiffEqRKN.jl | 26 ++ lib/OrdinaryDiffEqRKN/src/alg_utils.jl | 20 + lib/OrdinaryDiffEqRKN/src/algorithms.jl | 341 +++++++++++++++++ .../OrdinaryDiffEqRKN/src}/rkn_caches.jl | 2 +- .../src}/rkn_perform_step.jl | 2 +- .../OrdinaryDiffEqRKN/src}/rkn_tableaus.jl | 2 +- lib/OrdinaryDiffEqRKN/test/runtests.jl | 0 src/OrdinaryDiffEq.jl | 15 +- src/alg_utils.jl | 20 - src/algorithms.jl | 346 ------------------ test/runtests.jl | 6 + 12 files changed, 426 insertions(+), 377 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/Project.toml create mode 100644 lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl create mode 100644 lib/OrdinaryDiffEqRKN/src/alg_utils.jl create mode 100644 lib/OrdinaryDiffEqRKN/src/algorithms.jl rename {src/caches => lib/OrdinaryDiffEqRKN/src}/rkn_caches.jl (99%) rename {src/perform_step => lib/OrdinaryDiffEqRKN/src}/rkn_perform_step.jl (99%) rename {src/tableaus => lib/OrdinaryDiffEqRKN/src}/rkn_tableaus.jl (99%) create mode 100644 lib/OrdinaryDiffEqRKN/test/runtests.jl diff --git a/lib/OrdinaryDiffEqRKN/Project.toml b/lib/OrdinaryDiffEqRKN/Project.toml new file mode 100644 index 0000000000..05d567ea4f --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/Project.toml @@ -0,0 +1,23 @@ +name = "OrdinaryDiffEqRKN" +uuid = "af6ede74-add8-4cfd-b1df-9a4dbb109d7a" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[compat] +julia = "1.10" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl new file mode 100644 index 0000000000..bf08dc8b9d --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -0,0 +1,26 @@ +module OrdinaryDiffEqRKN + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, + trivial_limiter!, _ode_interpolant!, _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools + +include("algorithms.jl") +include("alg_utils.jl") +include("rkn_tableaus.jl") +include("rkn_caches.jl") +include("rkn_perform_step.jl") + +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + +end diff --git a/lib/OrdinaryDiffEqRKN/src/alg_utils.jl b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl new file mode 100644 index 0000000000..8f22b97e9c --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl @@ -0,0 +1,20 @@ +alg_extrapolates(alg::IRKN3) = true +alg_extrapolates(alg::IRKN4) = true + +alg_order(alg::IRKN3) = 3 +alg_order(alg::Nystrom4) = 4 +alg_order(alg::FineRKN4) = 4 +alg_order(alg::FineRKN5) = 5 +alg_order(alg::Nystrom4VelocityIndependent) = 4 +alg_order(alg::IRKN4) = 4 +alg_order(alg::Nystrom5VelocityIndependent) = 5 +alg_order(alg::DPRKN4) = 4 +alg_order(alg::DPRKN5) = 5 +alg_order(alg::DPRKN6) = 6 +alg_order(alg::DPRKN6FM) = 6 +alg_order(alg::DPRKN8) = 8 +alg_order(alg::DPRKN12) = 12 +alg_order(alg::ERKN4) = 4 +alg_order(alg::ERKN5) = 5 +alg_order(alg::ERKN7) = 7 +alg_order(alg::RKN4) = 4 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/algorithms.jl b/lib/OrdinaryDiffEqRKN/src/algorithms.jl new file mode 100644 index 0000000000..f1cf14fd93 --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/algorithms.jl @@ -0,0 +1,341 @@ +""" + IRKN3 + +Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not depend on the first derivative. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom4 + +A 4th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. + +In case the ODE Problem is not dependent on the first derivative consider using +[`Nystrom4VelocityIndependent`](@ref) to increase performance. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + FineRKN4() + +A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +In particular, this method allows the acceleration equation to depend on the velocity. + +## References + +``` +@article{fine1987low, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +``` +""" +struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + FineRKN5() + +A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +In particular, this method allows the acceleration equation to depend on the velocity. + +## References + +``` +@article{fine1987low, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +``` +""" +struct FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + Nystrom4VelocityIdependent + +A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). + +More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. + +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + IRKN4 + +Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not be dependent on the first derivative. + +Recommended for smooth problems with expensive functions to evaluate. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom5VelocityIndependent + +A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + DPRKN4 + +4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN5 + +5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. + +## References + +@article{Bettis1973ARN, +title={A Runge-Kutta Nystrom algorithm}, +author={Dale G. Bettis}, +journal={Celestial mechanics}, +year={1973}, +volume={8}, +pages={229-233}, +publisher={Springer} +} +""" +struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6 + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. + +## References + +@article{dormand1987runge, +title={Runge-kutta-nystrom triples}, +author={Dormand, JR and Prince, PJ}, +journal={Computers \\& Mathematics with Applications}, +volume={13}, +number={12}, +pages={937--949}, +year={1987}, +publisher={Elsevier} +} +""" +struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6FM + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain +when only the main solution points are considered. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN8 + +8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with +[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN12 + +12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Most efficient when high accuracy is needed. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN4 + +Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba2017embedded, +title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, +author={Demba, MA and Senu, N and Ismail, F}, +journal={Applied Mathematical Sciences}, +volume={11}, +number={17}, +pages={819--838}, +year={2017} +} +""" +struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN5 + +Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba20165, +title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, +author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, +journal={Mathematical and Computational Applications}, +volume={21}, +number={4}, +pages={46}, +year={2016}, +publisher={Multidisciplinary Digital Publishing Institute} +} +""" +struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN7 + +Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic Problems. + +## References + +@article{SimosOnHO, +title={On high order Runge-Kutta-Nystr{\"o}m pairs}, +author={Theodore E. Simos and Ch. Tsitouras}, +journal={J. Comput. Appl. Math.}, +volume={400}, +pages={113753} +} +""" +struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" +3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. + +Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. + +!!! warn + + This method is only fourth order for these systems, the method is second order otherwise! + +## References + +@article{MONTIJANO2024115533, +title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, +author = {J.I. Montijano and L. Rández and M. Calvo}, +journal = {Journal of Computational and Applied Mathematics}, +volume = {438}, +pages = {115533}, +year = {2024}, +} +""" +struct RKN4 <: OrdinaryDiffEqAlgorithm end \ No newline at end of file diff --git a/src/caches/rkn_caches.jl b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl similarity index 99% rename from src/caches/rkn_caches.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_caches.jl index 02014deef1..9119c135f3 100644 --- a/src/caches/rkn_caches.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl @@ -712,4 +712,4 @@ function alg_cache(alg::RKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} RKN4ConstantCache() -end +end \ No newline at end of file diff --git a/src/perform_step/rkn_perform_step.jl b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl similarity index 99% rename from src/perform_step/rkn_perform_step.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl index 70077f9268..1bacc80b93 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl @@ -1907,4 +1907,4 @@ end integrator.stats.nf += 2 integrator.stats.nf2 += 1 -end +end \ No newline at end of file diff --git a/src/tableaus/rkn_tableaus.jl b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl similarity index 99% rename from src/tableaus/rkn_tableaus.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl index c9b50e05ff..63dc672144 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl @@ -2469,4 +2469,4 @@ function DPRKN12ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloa convert(T, 0.015662325288859434), convert(T, -0.029155189630077964), convert(T, 0.025)) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/test/runtests.jl b/lib/OrdinaryDiffEqRKN/test/runtests.jl new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 86b4fafc04..e1a0319b57 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -155,7 +155,6 @@ include("caches/kencarp_kvaerno_caches.jl") include("caches/linear_caches.jl") include("caches/linear_nonlinear_caches.jl") include("caches/rosenbrock_caches.jl") -include("caches/rkn_caches.jl") include("caches/adams_bashforth_moulton_caches.jl") include("caches/nordsieck_caches.jl") include("caches/bdf_caches.jl") @@ -170,7 +169,6 @@ include("tableaus/verner_tableaus.jl") include("tableaus/rosenbrock_tableaus.jl") include("tableaus/sdirk_tableaus.jl") include("tableaus/firk_tableaus.jl") -include("tableaus/rkn_tableaus.jl") include("tableaus/qprk_tableaus.jl") include("integrators/type.jl") @@ -182,7 +180,6 @@ include("initialize_dae.jl") include("wrappers.jl") include("perform_step/fixed_timestep_perform_step.jl") -include("perform_step/rkn_perform_step.jl") include("perform_step/split_perform_step.jl") include("perform_step/linear_perform_step.jl") include("perform_step/exponential_rk_perform_step.jl") @@ -269,6 +266,13 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 +include("../lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl") +using ..OrdinaryDiffEqRKN +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + import PrecompileTools PrecompileTools.@compile_workload begin @@ -433,11 +437,6 @@ export SHLDDRK52, SHLDDRK_2N export SplitEuler -export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, - Nystrom5VelocityIndependent, - IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, - RKN4 - export AB3, AB4, AB5, ABM32, ABM43, ABM54 export VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 878b6d9f32..ff78c80801 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -355,8 +355,6 @@ alg_extrapolates(alg::DImplicitEuler) = true alg_extrapolates(alg::DABDF2) = true alg_extrapolates(alg::Trapezoid) = true alg_extrapolates(alg::SDIRK22) = true -alg_extrapolates(alg::IRKN4) = true -alg_extrapolates(alg::IRKN3) = true alg_extrapolates(alg::ABDF2) = true alg_extrapolates(alg::SBDF) = true alg_extrapolates(alg::MEBDF2) = true @@ -421,24 +419,6 @@ alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 alg_order(alg::FRK65) = 6 -alg_order(alg::IRKN3) = 3 -alg_order(alg::Nystrom4) = 4 -alg_order(alg::FineRKN4) = 4 -alg_order(alg::FineRKN5) = 5 -alg_order(alg::Nystrom4VelocityIndependent) = 4 -alg_order(alg::IRKN4) = 4 -alg_order(alg::Nystrom5VelocityIndependent) = 5 -alg_order(alg::DPRKN4) = 4 -alg_order(alg::DPRKN5) = 5 -alg_order(alg::DPRKN6) = 6 -alg_order(alg::DPRKN6FM) = 6 -alg_order(alg::DPRKN8) = 8 -alg_order(alg::DPRKN12) = 12 -alg_order(alg::ERKN4) = 4 -alg_order(alg::ERKN5) = 5 -alg_order(alg::ERKN7) = 7 -alg_order(alg::RKN4) = 4 - alg_order(alg::Midpoint) = 2 alg_order(alg::RK4) = 4 diff --git a/src/algorithms.jl b/src/algorithms.jl index 4152a7057b..807d425419 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -307,352 +307,6 @@ end ################################################################################ -# Nyström methods - -""" - IRKN3 - -Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not depend on the first derivative. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom4 - -A 4th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. - -In case the ODE Problem is not dependent on the first derivative consider using -[`Nystrom4VelocityIndependent`](@ref) to increase performance. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - FineRKN4() - -A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. -In particular, this method allows the acceleration equation to depend on the velocity. - -## References - -``` -@article{fine1987low, - title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, - author={Fine, Jerry Michael}, - journal={Computing}, - volume={38}, - number={4}, - pages={281--297}, - year={1987}, - publisher={Springer} -} -``` -""" -struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - FineRKN5() - -A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. -In particular, this method allows the acceleration equation to depend on the velocity. - -## References - -``` -@article{fine1987low, - title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, - author={Fine, Jerry Michael}, - journal={Computing}, - volume={38}, - number={4}, - pages={281--297}, - year={1987}, - publisher={Springer} -} -``` -""" -struct FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - Nystrom4VelocityIdependent - -A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). - -More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. - -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - IRKN4 - -Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not be dependent on the first derivative. - -Recommended for smooth problems with expensive functions to evaluate. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom5VelocityIndependent - -A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - DPRKN4 - -4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN5 - -5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. - -## References - -@article{Bettis1973ARN, -title={A Runge-Kutta Nystrom algorithm}, -author={Dale G. Bettis}, -journal={Celestial mechanics}, -year={1973}, -volume={8}, -pages={229-233}, -publisher={Springer} -} -""" -struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6 - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. - -## References - -@article{dormand1987runge, -title={Runge-kutta-nystrom triples}, -author={Dormand, JR and Prince, PJ}, -journal={Computers \\& Mathematics with Applications}, -volume={13}, -number={12}, -pages={937--949}, -year={1987}, -publisher={Elsevier} -} -""" -struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6FM - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain -when only the main solution points are considered. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN8 - -8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with -[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN12 - -12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Most efficient when high accuracy is needed. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN4 - -Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba2017embedded, -title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, -author={Demba, MA and Senu, N and Ismail, F}, -journal={Applied Mathematical Sciences}, -volume={11}, -number={17}, -pages={819--838}, -year={2017} -} -""" -struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN5 - -Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba20165, -title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, -author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, -journal={Mathematical and Computational Applications}, -volume={21}, -number={4}, -pages={46}, -year={2016}, -publisher={Multidisciplinary Digital Publishing Institute} -} -""" -struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN7 - -Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic Problems. - -## References - -@article{SimosOnHO, -title={On high order Runge-Kutta-Nystr{\"o}m pairs}, -author={Theodore E. Simos and Ch. Tsitouras}, -journal={J. Comput. Appl. Math.}, -volume={400}, -pages={113753} -} -""" -struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" -3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. - -Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. - -!!! warn - - This method is only fourth order for these systems, the method is second order otherwise! - -## References - -@article{MONTIJANO2024115533, -title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, -author = {J.I. Montijano and L. Rández and M. Calvo}, -journal = {Journal of Computational and Applied Mathematics}, -volume = {438}, -pages = {115533}, -year = {2024}, -} -""" -struct RKN4 <: OrdinaryDiffEqAlgorithm end - -################################################################################ - # Adams Bashforth and Adams moulton methods """ diff --git a/test/runtests.jl b/test/runtests.jl index 54b6bddff6..b5a0fd333e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,12 @@ function activate_feagin() Pkg.instantiate() end +function activate_symplectic_rk() + Pkg.activate("../lib/OrdinaryDiffEqSymplecticRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin From 0a213b283f0705072e11e11748120cb83506bd5b Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:13:36 +0530 Subject: [PATCH 10/28] RKN fix --- src/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 807d425419..7b7ae3ce6e 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -2289,7 +2289,7 @@ function PDIRK44(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{tru end ### Algorithm Groups -const MultistepAlgorithms = Union{IRKN3, IRKN4, +const MultistepAlgorithms = Union{ ABDF2, AB3, AB4, AB5, ABM32, ABM43, ABM54} From ddc95b07040dd92e0557b45fcd558a7f0e9856e0 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:30:32 +0530 Subject: [PATCH 11/28] interpolants.jl --- .../src/OrdinaryDiffEqRKN.jl | 1 + lib/OrdinaryDiffEqRKN/src/interpolants.jl | 110 +++++++++++++++++ src/dense/interpolants.jl | 111 ------------------ 3 files changed, 111 insertions(+), 111 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/src/interpolants.jl diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index bf08dc8b9d..6bb8989a25 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -16,6 +16,7 @@ include("algorithms.jl") include("alg_utils.jl") include("rkn_tableaus.jl") include("rkn_caches.jl") +include("interpolants.jl") include("rkn_perform_step.jl") export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl new file mode 100644 index 0000000000..61aea3a42d --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -0,0 +1,110 @@ +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), + uprev + + dt * Θ * + (duprev + + dt * Θ * (b1Θ * k1 + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +if idxs <= halfsize + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) +else + idxs = idxs - halfsize + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) +end +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +@inbounds @.. broadcast=false out.x[2]=uprev + + dt * Θ * + (duprev + + dt * Θ * + (b1Θ * k1 + + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) +@inbounds @.. broadcast=false out.x[1]=duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) +#for i in eachindex(out.x[1]) +# out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + +# b3Θ*k3[i] + +# b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) +# out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + +# bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) +#end +out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +isfirsthalf = idxs .<= halfsize +secondhalf = idxs .> halfsize +firstidxs = idxs[isfirsthalf] +secondidxs_shifted = idxs[secondhalf] +secondidxs = secondidxs_shifted .- halfsize + +@views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + + dt * Θ * + (duprev[secondidxs] + + dt * Θ * + (b1Θ * k1[secondidxs] + + b3Θ * k3[secondidxs] + + b4Θ * k4[secondidxs] + + b5Θ * k5[secondidxs] + + b6Θ * k6[secondidxs])) +@views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + + dt * Θ * + (bp1Θ * k1[firstidxs] + + bp3Θ * k3[firstidxs] + + bp4Θ * k4[firstidxs] + + bp5Θ * k5[firstidxs] + + bp6Θ * k6[firstidxs]) +out +end \ No newline at end of file diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 7356b8865d..6c20b7401e 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2756,114 +2756,3 @@ end duprev, uprev = y₀.x dtsq = dt^2 end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), - uprev + - dt * Θ * - (duprev + - dt * Θ * (b1Θ * k1 + b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - if idxs <= halfsize - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) - else - idxs = idxs - halfsize - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) - end -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - @inbounds @.. broadcast=false out.x[2]=uprev + - dt * Θ * - (duprev + - dt * Θ * - (b1Θ * k1 + - b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) - @inbounds @.. broadcast=false out.x[1]=duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) - #for i in eachindex(out.x[1]) - # out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + - # b3Θ*k3[i] + - # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) - # out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + - # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) - #end - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - isfirsthalf = idxs .<= halfsize - secondhalf = idxs .> halfsize - firstidxs = idxs[isfirsthalf] - secondidxs_shifted = idxs[secondhalf] - secondidxs = secondidxs_shifted .- halfsize - - @views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + - dt * Θ * - (duprev[secondidxs] + - dt * Θ * - (b1Θ * k1[secondidxs] + - b3Θ * k3[secondidxs] + - b4Θ * k4[secondidxs] + - b5Θ * k5[secondidxs] + - b6Θ * k6[secondidxs])) - @views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + - dt * Θ * - (bp1Θ * k1[firstidxs] + - bp3Θ * k3[firstidxs] + - bp4Θ * k4[firstidxs] + - bp5Θ * k5[firstidxs] + - bp6Θ * k6[firstidxs]) - out -end From 1ebfc0edd273bdda462adac50768ee059288d79c Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:37:28 +0530 Subject: [PATCH 12/28] interp_func.jl --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 1 + lib/OrdinaryDiffEqRKN/src/interp_func.jl | 7 +++++++ src/interp_func.jl | 7 ------- 3 files changed, 8 insertions(+), 7 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/src/interp_func.jl diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 6bb8989a25..66d37b62dd 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -16,6 +16,7 @@ include("algorithms.jl") include("alg_utils.jl") include("rkn_tableaus.jl") include("rkn_caches.jl") +include("interp_func.jl") include("interpolants.jl") include("rkn_perform_step.jl") diff --git a/lib/OrdinaryDiffEqRKN/src/interp_func.jl b/lib/OrdinaryDiffEqRKN/src/interp_func.jl new file mode 100644 index 0000000000..f5b2c03a2b --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interp_func.jl @@ -0,0 +1,7 @@ +function DiffEqBase.interp_summary(::Type{cacheType}, + dense::Bool) where { + cacheType <: + Union{DPRKN6ConstantCache, + DPRKN6Cache}} +dense ? "specialized 6th order interpolation" : "1st order linear" +end \ No newline at end of file diff --git a/src/interp_func.jl b/src/interp_func.jl index 14c95b7aba..4265219b7d 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -148,13 +148,6 @@ function DiffEqBase.interp_summary(::Type{cacheType}, caches = fieldtype(cacheType, :caches) join([DiffEqBase.interp_summary(ct, dense) for ct in fieldtypes(caches)], ", ") end -function DiffEqBase.interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{DPRKN6ConstantCache, - DPRKN6Cache}} - dense ? "specialized 6th order interpolation" : "1st order linear" -end function (interp::InterpolationData)(tvals, idxs, deriv, p, continuity::Symbol = :left) ode_interpolation(tvals, interp, idxs, deriv, p, continuity) From 0e13020a2da0d81807e71c97d8808a85bae65abb Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:43:28 +0530 Subject: [PATCH 13/28] OrdinaryDiffEqPartitionedAlgorithm --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 66d37b62dd..588e86157c 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -6,6 +6,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, du_alias_or_new, @@ -25,4 +26,4 @@ export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, RKN4 -end +end \ No newline at end of file From 3323dfbddeb7f64d751d1b30dbf36c36a68b7546 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:48:31 +0530 Subject: [PATCH 14/28] Changes --- lib/OrdinaryDiffEqRKN/src/interpolants.jl | 23 +++++++++++++++++++++ src/dense/interpolants.jl | 25 +---------------------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl index 61aea3a42d..bad93428b4 100644 --- a/lib/OrdinaryDiffEqRKN/src/interpolants.jl +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -1,3 +1,26 @@ +@def dprkn6pre0 begin + @dprkn6unpack + b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) + b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) + b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) + b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) + b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) + + bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) + bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) + bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) + bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) + bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) + + kk1, kk2, kk3 = k + k1, k2 = kk1.x + k3, k4 = kk2.x + k5, k6 = kk3.x + + duprev, uprev = y₀.x + dtsq = dt^2 +end + @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 6c20b7401e..275b70deaf 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2732,27 +2732,4 @@ end else @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache end -end - -@def dprkn6pre0 begin - @dprkn6unpack - b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) - b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) - b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) - b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) - b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) - - bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) - bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) - bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) - bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) - bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) - - kk1, kk2, kk3 = k - k1, k2 = kk1.x - k3, k4 = kk2.x - k5, k6 = kk3.x - - duprev, uprev = y₀.x - dtsq = dt^2 -end +end \ No newline at end of file From f91395dcc71b324ec743f5e5a5fa56b93600f997 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:51:34 +0530 Subject: [PATCH 15/28] @def --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 588e86157c..ca97625d07 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -12,6 +12,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, constvalue, _unwrap_val, du_alias_or_new, trivial_limiter!, _ode_interpolant!, _ode_addsteps! using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools +using DiffEqBase: @def include("algorithms.jl") include("alg_utils.jl") From 5984917ed050e6103b7e041565e3a8e8c46239eb Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:56:13 +0530 Subject: [PATCH 16/28] dprkn6unpack --- lib/OrdinaryDiffEqRKN/src/interpolants.jl | 8 ++++++++ src/dense/interpolants.jl | 9 +-------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl index bad93428b4..4e919f9f47 100644 --- a/lib/OrdinaryDiffEqRKN/src/interpolants.jl +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -1,3 +1,11 @@ +@def dprkn6unpack begin + if cache isa OrdinaryDiffEqMutableCache + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab + else + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache + end +end + @def dprkn6pre0 begin @dprkn6unpack b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 275b70deaf..ba99d9f4aa 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2725,11 +2725,4 @@ end end """ -""" -@def dprkn6unpack begin - if cache isa OrdinaryDiffEqMutableCache - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab - else - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache - end -end \ No newline at end of file +""" \ No newline at end of file From 8b1a6e9d4136cd270800d54cc0630fb80e155cb5 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 13:03:49 +0530 Subject: [PATCH 17/28] Tight Loop Macros --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index ca97625d07..b8823c9cd2 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -12,7 +12,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, constvalue, _unwrap_val, du_alias_or_new, trivial_limiter!, _ode_interpolant!, _ode_addsteps! using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools -using DiffEqBase: @def +using DiffEqBase: @def, @tight_loop_macros include("algorithms.jl") include("alg_utils.jl") From 97456ecd05ae9a28efbd7130e6380733484c8e0c Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 11:25:52 +0530 Subject: [PATCH 18/28] RKN Methods --- lib/OrdinaryDiffEqRKN/Project.toml | 23 ++ .../src/OrdinaryDiffEqRKN.jl | 26 ++ lib/OrdinaryDiffEqRKN/src/alg_utils.jl | 20 + lib/OrdinaryDiffEqRKN/src/algorithms.jl | 341 +++++++++++++++++ .../OrdinaryDiffEqRKN/src}/rkn_caches.jl | 2 +- .../src}/rkn_perform_step.jl | 2 +- .../OrdinaryDiffEqRKN/src}/rkn_tableaus.jl | 2 +- lib/OrdinaryDiffEqRKN/test/runtests.jl | 0 src/OrdinaryDiffEq.jl | 15 +- src/alg_utils.jl | 20 - src/algorithms.jl | 346 ------------------ test/runtests.jl | 6 + 12 files changed, 426 insertions(+), 377 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/Project.toml create mode 100644 lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl create mode 100644 lib/OrdinaryDiffEqRKN/src/alg_utils.jl create mode 100644 lib/OrdinaryDiffEqRKN/src/algorithms.jl rename {src/caches => lib/OrdinaryDiffEqRKN/src}/rkn_caches.jl (99%) rename {src/perform_step => lib/OrdinaryDiffEqRKN/src}/rkn_perform_step.jl (99%) rename {src/tableaus => lib/OrdinaryDiffEqRKN/src}/rkn_tableaus.jl (99%) create mode 100644 lib/OrdinaryDiffEqRKN/test/runtests.jl diff --git a/lib/OrdinaryDiffEqRKN/Project.toml b/lib/OrdinaryDiffEqRKN/Project.toml new file mode 100644 index 0000000000..05d567ea4f --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/Project.toml @@ -0,0 +1,23 @@ +name = "OrdinaryDiffEqRKN" +uuid = "af6ede74-add8-4cfd-b1df-9a4dbb109d7a" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[compat] +julia = "1.10" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl new file mode 100644 index 0000000000..bf08dc8b9d --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -0,0 +1,26 @@ +module OrdinaryDiffEqRKN + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, + trivial_limiter!, _ode_interpolant!, _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools + +include("algorithms.jl") +include("alg_utils.jl") +include("rkn_tableaus.jl") +include("rkn_caches.jl") +include("rkn_perform_step.jl") + +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + +end diff --git a/lib/OrdinaryDiffEqRKN/src/alg_utils.jl b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl new file mode 100644 index 0000000000..8f22b97e9c --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl @@ -0,0 +1,20 @@ +alg_extrapolates(alg::IRKN3) = true +alg_extrapolates(alg::IRKN4) = true + +alg_order(alg::IRKN3) = 3 +alg_order(alg::Nystrom4) = 4 +alg_order(alg::FineRKN4) = 4 +alg_order(alg::FineRKN5) = 5 +alg_order(alg::Nystrom4VelocityIndependent) = 4 +alg_order(alg::IRKN4) = 4 +alg_order(alg::Nystrom5VelocityIndependent) = 5 +alg_order(alg::DPRKN4) = 4 +alg_order(alg::DPRKN5) = 5 +alg_order(alg::DPRKN6) = 6 +alg_order(alg::DPRKN6FM) = 6 +alg_order(alg::DPRKN8) = 8 +alg_order(alg::DPRKN12) = 12 +alg_order(alg::ERKN4) = 4 +alg_order(alg::ERKN5) = 5 +alg_order(alg::ERKN7) = 7 +alg_order(alg::RKN4) = 4 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/algorithms.jl b/lib/OrdinaryDiffEqRKN/src/algorithms.jl new file mode 100644 index 0000000000..f1cf14fd93 --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/algorithms.jl @@ -0,0 +1,341 @@ +""" + IRKN3 + +Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not depend on the first derivative. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom4 + +A 4th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. + +In case the ODE Problem is not dependent on the first derivative consider using +[`Nystrom4VelocityIndependent`](@ref) to increase performance. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + FineRKN4() + +A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +In particular, this method allows the acceleration equation to depend on the velocity. + +## References + +``` +@article{fine1987low, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +``` +""" +struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + FineRKN5() + +A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +In particular, this method allows the acceleration equation to depend on the velocity. + +## References + +``` +@article{fine1987low, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +``` +""" +struct FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + Nystrom4VelocityIdependent + +A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). + +More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. + +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + IRKN4 + +Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not be dependent on the first derivative. + +Recommended for smooth problems with expensive functions to evaluate. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom5VelocityIndependent + +A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + DPRKN4 + +4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN5 + +5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. + +## References + +@article{Bettis1973ARN, +title={A Runge-Kutta Nystrom algorithm}, +author={Dale G. Bettis}, +journal={Celestial mechanics}, +year={1973}, +volume={8}, +pages={229-233}, +publisher={Springer} +} +""" +struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6 + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. + +## References + +@article{dormand1987runge, +title={Runge-kutta-nystrom triples}, +author={Dormand, JR and Prince, PJ}, +journal={Computers \\& Mathematics with Applications}, +volume={13}, +number={12}, +pages={937--949}, +year={1987}, +publisher={Elsevier} +} +""" +struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6FM + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain +when only the main solution points are considered. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN8 + +8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with +[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN12 + +12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Most efficient when high accuracy is needed. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN4 + +Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba2017embedded, +title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, +author={Demba, MA and Senu, N and Ismail, F}, +journal={Applied Mathematical Sciences}, +volume={11}, +number={17}, +pages={819--838}, +year={2017} +} +""" +struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN5 + +Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba20165, +title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, +author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, +journal={Mathematical and Computational Applications}, +volume={21}, +number={4}, +pages={46}, +year={2016}, +publisher={Multidisciplinary Digital Publishing Institute} +} +""" +struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN7 + +Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic Problems. + +## References + +@article{SimosOnHO, +title={On high order Runge-Kutta-Nystr{\"o}m pairs}, +author={Theodore E. Simos and Ch. Tsitouras}, +journal={J. Comput. Appl. Math.}, +volume={400}, +pages={113753} +} +""" +struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" +3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. + +Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. + +!!! warn + + This method is only fourth order for these systems, the method is second order otherwise! + +## References + +@article{MONTIJANO2024115533, +title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, +author = {J.I. Montijano and L. Rández and M. Calvo}, +journal = {Journal of Computational and Applied Mathematics}, +volume = {438}, +pages = {115533}, +year = {2024}, +} +""" +struct RKN4 <: OrdinaryDiffEqAlgorithm end \ No newline at end of file diff --git a/src/caches/rkn_caches.jl b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl similarity index 99% rename from src/caches/rkn_caches.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_caches.jl index 02014deef1..9119c135f3 100644 --- a/src/caches/rkn_caches.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl @@ -712,4 +712,4 @@ function alg_cache(alg::RKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} RKN4ConstantCache() -end +end \ No newline at end of file diff --git a/src/perform_step/rkn_perform_step.jl b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl similarity index 99% rename from src/perform_step/rkn_perform_step.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl index 70077f9268..1bacc80b93 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl @@ -1907,4 +1907,4 @@ end integrator.stats.nf += 2 integrator.stats.nf2 += 1 -end +end \ No newline at end of file diff --git a/src/tableaus/rkn_tableaus.jl b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl similarity index 99% rename from src/tableaus/rkn_tableaus.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl index c9b50e05ff..63dc672144 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl @@ -2469,4 +2469,4 @@ function DPRKN12ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloa convert(T, 0.015662325288859434), convert(T, -0.029155189630077964), convert(T, 0.025)) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/test/runtests.jl b/lib/OrdinaryDiffEqRKN/test/runtests.jl new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 86b4fafc04..e1a0319b57 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -155,7 +155,6 @@ include("caches/kencarp_kvaerno_caches.jl") include("caches/linear_caches.jl") include("caches/linear_nonlinear_caches.jl") include("caches/rosenbrock_caches.jl") -include("caches/rkn_caches.jl") include("caches/adams_bashforth_moulton_caches.jl") include("caches/nordsieck_caches.jl") include("caches/bdf_caches.jl") @@ -170,7 +169,6 @@ include("tableaus/verner_tableaus.jl") include("tableaus/rosenbrock_tableaus.jl") include("tableaus/sdirk_tableaus.jl") include("tableaus/firk_tableaus.jl") -include("tableaus/rkn_tableaus.jl") include("tableaus/qprk_tableaus.jl") include("integrators/type.jl") @@ -182,7 +180,6 @@ include("initialize_dae.jl") include("wrappers.jl") include("perform_step/fixed_timestep_perform_step.jl") -include("perform_step/rkn_perform_step.jl") include("perform_step/split_perform_step.jl") include("perform_step/linear_perform_step.jl") include("perform_step/exponential_rk_perform_step.jl") @@ -269,6 +266,13 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 +include("../lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl") +using ..OrdinaryDiffEqRKN +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + import PrecompileTools PrecompileTools.@compile_workload begin @@ -433,11 +437,6 @@ export SHLDDRK52, SHLDDRK_2N export SplitEuler -export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, - Nystrom5VelocityIndependent, - IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, - RKN4 - export AB3, AB4, AB5, ABM32, ABM43, ABM54 export VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 878b6d9f32..ff78c80801 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -355,8 +355,6 @@ alg_extrapolates(alg::DImplicitEuler) = true alg_extrapolates(alg::DABDF2) = true alg_extrapolates(alg::Trapezoid) = true alg_extrapolates(alg::SDIRK22) = true -alg_extrapolates(alg::IRKN4) = true -alg_extrapolates(alg::IRKN3) = true alg_extrapolates(alg::ABDF2) = true alg_extrapolates(alg::SBDF) = true alg_extrapolates(alg::MEBDF2) = true @@ -421,24 +419,6 @@ alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 alg_order(alg::FRK65) = 6 -alg_order(alg::IRKN3) = 3 -alg_order(alg::Nystrom4) = 4 -alg_order(alg::FineRKN4) = 4 -alg_order(alg::FineRKN5) = 5 -alg_order(alg::Nystrom4VelocityIndependent) = 4 -alg_order(alg::IRKN4) = 4 -alg_order(alg::Nystrom5VelocityIndependent) = 5 -alg_order(alg::DPRKN4) = 4 -alg_order(alg::DPRKN5) = 5 -alg_order(alg::DPRKN6) = 6 -alg_order(alg::DPRKN6FM) = 6 -alg_order(alg::DPRKN8) = 8 -alg_order(alg::DPRKN12) = 12 -alg_order(alg::ERKN4) = 4 -alg_order(alg::ERKN5) = 5 -alg_order(alg::ERKN7) = 7 -alg_order(alg::RKN4) = 4 - alg_order(alg::Midpoint) = 2 alg_order(alg::RK4) = 4 diff --git a/src/algorithms.jl b/src/algorithms.jl index 4152a7057b..807d425419 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -307,352 +307,6 @@ end ################################################################################ -# Nyström methods - -""" - IRKN3 - -Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not depend on the first derivative. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom4 - -A 4th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. - -In case the ODE Problem is not dependent on the first derivative consider using -[`Nystrom4VelocityIndependent`](@ref) to increase performance. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - FineRKN4() - -A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. -In particular, this method allows the acceleration equation to depend on the velocity. - -## References - -``` -@article{fine1987low, - title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, - author={Fine, Jerry Michael}, - journal={Computing}, - volume={38}, - number={4}, - pages={281--297}, - year={1987}, - publisher={Springer} -} -``` -""" -struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - FineRKN5() - -A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. -In particular, this method allows the acceleration equation to depend on the velocity. - -## References - -``` -@article{fine1987low, - title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, - author={Fine, Jerry Michael}, - journal={Computing}, - volume={38}, - number={4}, - pages={281--297}, - year={1987}, - publisher={Springer} -} -``` -""" -struct FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - Nystrom4VelocityIdependent - -A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). - -More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. - -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - IRKN4 - -Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not be dependent on the first derivative. - -Recommended for smooth problems with expensive functions to evaluate. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom5VelocityIndependent - -A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - DPRKN4 - -4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN5 - -5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. - -## References - -@article{Bettis1973ARN, -title={A Runge-Kutta Nystrom algorithm}, -author={Dale G. Bettis}, -journal={Celestial mechanics}, -year={1973}, -volume={8}, -pages={229-233}, -publisher={Springer} -} -""" -struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6 - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. - -## References - -@article{dormand1987runge, -title={Runge-kutta-nystrom triples}, -author={Dormand, JR and Prince, PJ}, -journal={Computers \\& Mathematics with Applications}, -volume={13}, -number={12}, -pages={937--949}, -year={1987}, -publisher={Elsevier} -} -""" -struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6FM - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain -when only the main solution points are considered. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN8 - -8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with -[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN12 - -12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Most efficient when high accuracy is needed. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN4 - -Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba2017embedded, -title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, -author={Demba, MA and Senu, N and Ismail, F}, -journal={Applied Mathematical Sciences}, -volume={11}, -number={17}, -pages={819--838}, -year={2017} -} -""" -struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN5 - -Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba20165, -title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, -author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, -journal={Mathematical and Computational Applications}, -volume={21}, -number={4}, -pages={46}, -year={2016}, -publisher={Multidisciplinary Digital Publishing Institute} -} -""" -struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN7 - -Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic Problems. - -## References - -@article{SimosOnHO, -title={On high order Runge-Kutta-Nystr{\"o}m pairs}, -author={Theodore E. Simos and Ch. Tsitouras}, -journal={J. Comput. Appl. Math.}, -volume={400}, -pages={113753} -} -""" -struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" -3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. - -Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. - -!!! warn - - This method is only fourth order for these systems, the method is second order otherwise! - -## References - -@article{MONTIJANO2024115533, -title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, -author = {J.I. Montijano and L. Rández and M. Calvo}, -journal = {Journal of Computational and Applied Mathematics}, -volume = {438}, -pages = {115533}, -year = {2024}, -} -""" -struct RKN4 <: OrdinaryDiffEqAlgorithm end - -################################################################################ - # Adams Bashforth and Adams moulton methods """ diff --git a/test/runtests.jl b/test/runtests.jl index 54b6bddff6..b5a0fd333e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,12 @@ function activate_feagin() Pkg.instantiate() end +function activate_symplectic_rk() + Pkg.activate("../lib/OrdinaryDiffEqSymplecticRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin From fde1a1ded83ee33ce4ba7fa24072e881cb94558e Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:13:36 +0530 Subject: [PATCH 19/28] RKN fix --- src/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 807d425419..7b7ae3ce6e 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -2289,7 +2289,7 @@ function PDIRK44(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{tru end ### Algorithm Groups -const MultistepAlgorithms = Union{IRKN3, IRKN4, +const MultistepAlgorithms = Union{ ABDF2, AB3, AB4, AB5, ABM32, ABM43, ABM54} From 773407ea09b7cab4c750342f40cf0be602017081 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:30:32 +0530 Subject: [PATCH 20/28] interpolants.jl --- .../src/OrdinaryDiffEqRKN.jl | 1 + lib/OrdinaryDiffEqRKN/src/interpolants.jl | 110 +++++++++++++++++ src/dense/interpolants.jl | 111 ------------------ 3 files changed, 111 insertions(+), 111 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/src/interpolants.jl diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index bf08dc8b9d..6bb8989a25 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -16,6 +16,7 @@ include("algorithms.jl") include("alg_utils.jl") include("rkn_tableaus.jl") include("rkn_caches.jl") +include("interpolants.jl") include("rkn_perform_step.jl") export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl new file mode 100644 index 0000000000..61aea3a42d --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -0,0 +1,110 @@ +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), + uprev + + dt * Θ * + (duprev + + dt * Θ * (b1Θ * k1 + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +if idxs <= halfsize + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) +else + idxs = idxs - halfsize + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) +end +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +@inbounds @.. broadcast=false out.x[2]=uprev + + dt * Θ * + (duprev + + dt * Θ * + (b1Θ * k1 + + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) +@inbounds @.. broadcast=false out.x[1]=duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) +#for i in eachindex(out.x[1]) +# out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + +# b3Θ*k3[i] + +# b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) +# out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + +# bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) +#end +out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +isfirsthalf = idxs .<= halfsize +secondhalf = idxs .> halfsize +firstidxs = idxs[isfirsthalf] +secondidxs_shifted = idxs[secondhalf] +secondidxs = secondidxs_shifted .- halfsize + +@views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + + dt * Θ * + (duprev[secondidxs] + + dt * Θ * + (b1Θ * k1[secondidxs] + + b3Θ * k3[secondidxs] + + b4Θ * k4[secondidxs] + + b5Θ * k5[secondidxs] + + b6Θ * k6[secondidxs])) +@views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + + dt * Θ * + (bp1Θ * k1[firstidxs] + + bp3Θ * k3[firstidxs] + + bp4Θ * k4[firstidxs] + + bp5Θ * k5[firstidxs] + + bp6Θ * k6[firstidxs]) +out +end \ No newline at end of file diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 7356b8865d..6c20b7401e 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2756,114 +2756,3 @@ end duprev, uprev = y₀.x dtsq = dt^2 end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), - uprev + - dt * Θ * - (duprev + - dt * Θ * (b1Θ * k1 + b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - if idxs <= halfsize - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) - else - idxs = idxs - halfsize - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) - end -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - @inbounds @.. broadcast=false out.x[2]=uprev + - dt * Θ * - (duprev + - dt * Θ * - (b1Θ * k1 + - b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) - @inbounds @.. broadcast=false out.x[1]=duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) - #for i in eachindex(out.x[1]) - # out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + - # b3Θ*k3[i] + - # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) - # out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + - # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) - #end - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - isfirsthalf = idxs .<= halfsize - secondhalf = idxs .> halfsize - firstidxs = idxs[isfirsthalf] - secondidxs_shifted = idxs[secondhalf] - secondidxs = secondidxs_shifted .- halfsize - - @views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + - dt * Θ * - (duprev[secondidxs] + - dt * Θ * - (b1Θ * k1[secondidxs] + - b3Θ * k3[secondidxs] + - b4Θ * k4[secondidxs] + - b5Θ * k5[secondidxs] + - b6Θ * k6[secondidxs])) - @views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + - dt * Θ * - (bp1Θ * k1[firstidxs] + - bp3Θ * k3[firstidxs] + - bp4Θ * k4[firstidxs] + - bp5Θ * k5[firstidxs] + - bp6Θ * k6[firstidxs]) - out -end From 5feeb2d6cbf60114428c184ae18ccfe0b78a9479 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:37:28 +0530 Subject: [PATCH 21/28] interp_func.jl --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 1 + lib/OrdinaryDiffEqRKN/src/interp_func.jl | 7 +++++++ src/interp_func.jl | 7 ------- 3 files changed, 8 insertions(+), 7 deletions(-) create mode 100644 lib/OrdinaryDiffEqRKN/src/interp_func.jl diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 6bb8989a25..66d37b62dd 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -16,6 +16,7 @@ include("algorithms.jl") include("alg_utils.jl") include("rkn_tableaus.jl") include("rkn_caches.jl") +include("interp_func.jl") include("interpolants.jl") include("rkn_perform_step.jl") diff --git a/lib/OrdinaryDiffEqRKN/src/interp_func.jl b/lib/OrdinaryDiffEqRKN/src/interp_func.jl new file mode 100644 index 0000000000..f5b2c03a2b --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interp_func.jl @@ -0,0 +1,7 @@ +function DiffEqBase.interp_summary(::Type{cacheType}, + dense::Bool) where { + cacheType <: + Union{DPRKN6ConstantCache, + DPRKN6Cache}} +dense ? "specialized 6th order interpolation" : "1st order linear" +end \ No newline at end of file diff --git a/src/interp_func.jl b/src/interp_func.jl index 14c95b7aba..4265219b7d 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -148,13 +148,6 @@ function DiffEqBase.interp_summary(::Type{cacheType}, caches = fieldtype(cacheType, :caches) join([DiffEqBase.interp_summary(ct, dense) for ct in fieldtypes(caches)], ", ") end -function DiffEqBase.interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{DPRKN6ConstantCache, - DPRKN6Cache}} - dense ? "specialized 6th order interpolation" : "1st order linear" -end function (interp::InterpolationData)(tvals, idxs, deriv, p, continuity::Symbol = :left) ode_interpolation(tvals, interp, idxs, deriv, p, continuity) From ba6ce4c5e3d6eb9df94b610d99d8dbaa6599aa71 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:43:28 +0530 Subject: [PATCH 22/28] OrdinaryDiffEqPartitionedAlgorithm --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 66d37b62dd..588e86157c 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -6,6 +6,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, du_alias_or_new, @@ -25,4 +26,4 @@ export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, RKN4 -end +end \ No newline at end of file From ffeb107157942048b9b28fc0ba50d26e16ed3134 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:48:31 +0530 Subject: [PATCH 23/28] Changes --- lib/OrdinaryDiffEqRKN/src/interpolants.jl | 23 +++++++++++++++++++++ src/dense/interpolants.jl | 25 +---------------------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl index 61aea3a42d..bad93428b4 100644 --- a/lib/OrdinaryDiffEqRKN/src/interpolants.jl +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -1,3 +1,26 @@ +@def dprkn6pre0 begin + @dprkn6unpack + b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) + b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) + b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) + b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) + b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) + + bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) + bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) + bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) + bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) + bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) + + kk1, kk2, kk3 = k + k1, k2 = kk1.x + k3, k4 = kk2.x + k5, k6 = kk3.x + + duprev, uprev = y₀.x + dtsq = dt^2 +end + @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 6c20b7401e..275b70deaf 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2732,27 +2732,4 @@ end else @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache end -end - -@def dprkn6pre0 begin - @dprkn6unpack - b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) - b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) - b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) - b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) - b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) - - bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) - bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) - bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) - bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) - bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) - - kk1, kk2, kk3 = k - k1, k2 = kk1.x - k3, k4 = kk2.x - k5, k6 = kk3.x - - duprev, uprev = y₀.x - dtsq = dt^2 -end +end \ No newline at end of file From f4ff7875fbd79350f73d72c1817acdc28f8b30f5 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:51:34 +0530 Subject: [PATCH 24/28] @def --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index 588e86157c..ca97625d07 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -12,6 +12,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, constvalue, _unwrap_val, du_alias_or_new, trivial_limiter!, _ode_interpolant!, _ode_addsteps! using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools +using DiffEqBase: @def include("algorithms.jl") include("alg_utils.jl") From a797132e9409802c5f12005b6ed09d16e1c8e671 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 12:56:13 +0530 Subject: [PATCH 25/28] dprkn6unpack --- lib/OrdinaryDiffEqRKN/src/interpolants.jl | 8 ++++++++ src/dense/interpolants.jl | 9 +-------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl index bad93428b4..4e919f9f47 100644 --- a/lib/OrdinaryDiffEqRKN/src/interpolants.jl +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -1,3 +1,11 @@ +@def dprkn6unpack begin + if cache isa OrdinaryDiffEqMutableCache + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab + else + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache + end +end + @def dprkn6pre0 begin @dprkn6unpack b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 275b70deaf..ba99d9f4aa 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2725,11 +2725,4 @@ end end """ -""" -@def dprkn6unpack begin - if cache isa OrdinaryDiffEqMutableCache - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab - else - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache - end -end \ No newline at end of file +""" \ No newline at end of file From 2d0d491d8bafbafb27b4006a6670d54d0082836f Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 2 Jul 2024 13:03:49 +0530 Subject: [PATCH 26/28] Tight Loop Macros --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index ca97625d07..b8823c9cd2 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -12,7 +12,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, constvalue, _unwrap_val, du_alias_or_new, trivial_limiter!, _ode_interpolant!, _ode_addsteps! using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools -using DiffEqBase: @def +using DiffEqBase: @def, @tight_loop_macros include("algorithms.jl") include("alg_utils.jl") From cc3b6c583174bd867b8d5d73c24fd21df4fbefd5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 6 Jul 2024 02:19:39 -0400 Subject: [PATCH 27/28] Update lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index b8823c9cd2..d1b0e197d3 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -2,7 +2,7 @@ module OrdinaryDiffEqRKN import OrdinaryDiffEq: alg_order, calculate_residuals!, initialize!, perform_step!, @unpack, unwrap_alg, - calculate_residuals, + calculate_residuals, alg_extrapolates, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptivePartitionedAlgorithm, From 96ffdfb760cb712818420e4289cc990ad45bd026 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 6 Jul 2024 14:38:32 -0400 Subject: [PATCH 28/28] Update lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl --- lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl index d1b0e197d3..7ddfbc2c8e 100644 --- a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -9,7 +9,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, OrdinaryDiffEqPartitionedAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, - constvalue, _unwrap_val, du_alias_or_new, + constvalue, _unwrap_val, du_alias_or_new, _ode_interpolant, trivial_limiter!, _ode_interpolant!, _ode_addsteps! using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools using DiffEqBase: @def, @tight_loop_macros