Fixed the interpolation precision bug and added more tests and work-p…

…recision benchmark generation codes

CMakeLists.txt

@@ -52,7 +52,15 @@ file(GLOB FLINT_SRC src/*.f90)
 add_library(${FLINTBIN} STATIC ${FLINT_SRC})
 target_include_directories(${FLINTBIN} INTERFACE ${FLINT_INSTALL_LIB_DIR})
 
-# The test executable 
+# The test executable for Work precision tests 
+set(TestFLINT_WP_SRC tests/test_WP.f90)
+set(TestFLINT_WPBIN TestFLINT_WP)
+add_executable(${TestFLINT_WPBIN} ${TestFLINT_WP_SRC})
+
+# link
+target_link_libraries(${TestFLINT_WPBIN} ${FLINTBIN})
+
+# The test executable for CR3BP test
 set(TestFLINT_SRC tests/test_CR3BP.f90)
 set(TestFLINTBIN TestFLINT_CR3BP)
 add_executable(${TestFLINTBIN} ${TestFLINT_SRC})
@@ -68,6 +76,13 @@ add_executable(${TestFLINT_LorenzBIN} ${TestFLINT_Lorenz_SRC})
 # link
 target_link_libraries(${TestFLINT_LorenzBIN} ${FLINTBIN})
 
+# The test executable for TBP test
+set(TestFLINT_SRC tests/test_TBP.f90)
+set(TestFLINTBIN TestFLINT_TBP)
+add_executable(${TestFLINTBIN} ${TestFLINT_SRC})
+
+# link
+target_link_libraries(${TestFLINTBIN} ${FLINTBIN})
 
 ###########################
 # Compile Options for FLINT

README.md

@@ -1,16 +1,17 @@
 # FLINT
 
+
 ## Fortran Library for numerical INTegration of differential equations
 
-[![Build Status](https://travis-ci.com/princemahajan/FLINT.svg?branch=master)](https://travis-ci.com/princemahajan/FLINT) 
 
 ### Author
 
 Bharat Mahajan
 
+
 ### Copyright
 
-Copyright 2021 Bharat Mahajan <br><br>
+Copyright 2024 Bharat Mahajan <br><br>
 The initial code was written by Bharat Mahajan at Odyssey Space Research LLC, Houston, TX as part of
 the work under contract no. 80JSC017D0001 with NASA-Johnson Space Center. 
 FLINT source code is licensed under the Apache License, Version 2.0 (the "License")
@@ -20,34 +21,47 @@ His original codes are available at http://www.unige.ch/~hairer/software.html.
 The coefficients for Verner65E and Verner98R methods were derived by Jim Verner, and 
 are available at http://people.math.sfu.ca/~jverner/.    
 
+
 ### Introduction
 
-FLINT is a modern object-oriented Fortran library that provides four adaptive step-size explicit Runge-Kutta (ERK) methods of order 5, 6, 8, and 9 along with dense-output and multiple event-detection support for each of the methods. The code is written such that any other ERK method can be implemented by including its coefficients with minimum changes required in the code. The DOP853 integrator is the default method chosen, and its implementation is hand-optimized specific to its coefficients. For other integrators, a generic routine for step-integration is implemented. This generic routine supports both FSAL and non-FSAL methods. Dense output is supported with delayed interpolation. When interpolation is enabled, FLINT computes the interpolation coefficients during the integration and stores them in internal memory. Thereafter, the interpolation method can be used any number of times to find the solution values at any user-specified grid within the initial and final conditions. Interpolation is much faster than integration, as the coefficients are all precomputed during the integration. Multiple event detection is supported for each integrator along with many features such as event root-finding, event step-size, event actions. In a nutshell, the features are:
+FLINT is a modern object-oriented Fortran library that provides four adaptive step-size explicit Runge-Kutta (ERK) methods of order 5, 6, 8, and 9 along with dense-output and multiple event-detection support for each of the methods. The code is written such that any other ERK method can be implemented by including its coefficients with minimum changes required in the code. The DOP54 and DOP853 integrators  step implementation is hand-optimized and for other integrators, a generic routine for step integration is implemented. This generic routine supports both FSAL and non-FSAL methods. Dense output is supported with delayed interpolation. When interpolation is enabled, FLINT computes the interpolation coefficients during the integration and stores them in internal memory. Thereafter, the interpolation method can be used any number of times to find the solution values at any user-specified grid within the initial and final conditions. Interpolation is much faster than integration, as the coefficients are all precomputed during the integration. Multiple event detection is supported for each integrator along with many features such as event root-finding, event step-size, event actions. In a nutshell, the features are:
                   
 + Modern object-oriented, thread-safe, and optimized Fortran code
-+ 4 Adaptive-step (fixed step-size also supported) ERK integrators: DOP54, DOP853, Verner98R, Verner65E
++ 4 Adaptive-step ERK integrators: DOP54, DOP853, Verner98R, Verner65E
++ Fixed step size for integration is supported as well as the ability to propagate for a single step and return control to the user 
 + Any other ERK method can be implemented by including their coefficients
 + Dense output with delayed interpolation (integrate once, interpolate as many times)
 + Multiple event-detection as well as finding location of events using root-finding  (Brent's algorithm) with static and dynamic event masking
 + Ability to set a maximum delay (referred to here as event step-size) after which events are guaranteed to be detected
 + Ability to restart the integration or change solution on the detection of events
 + Stiffness detection
 
+The following figure shows the multiple event detection capability of FLINT, in which the X-crossings in decreasing and Y-crossings in increasing direction of a three-body orbit are detected and reported to the user.
+
+
+![FINT's event-detection](media/FLINT_Events.png)
+
 
-### Performance benchmark against Julia's Differential Equation package
+### Performance Benchmarks
 
-The latest FLINT code (compiled using Intel Fortran compiler) is tested against Julia's DifferentialEquations package (https://diffeq.sciml.ai/stable/) and FLINT appears to be faster with and without event detection in the most cases as shown in the following screenshots. The Julia test code along with results are provided in the tests folder on FLINT's GitHub repository https://github.com/princemahajan/FLINT. 
+The latest FLINT code (compiled using Intel Fortran compiler ifort) is tested against Julia's (ver 1.10.0) DifferentialEquations package ver 7.12.0 (https://diffeq.sciml.ai/stable/) for different problems with and without event detection. The Julia test code along with results are provided in tests folder in FLINT's GitHub repository https://github.com/princemahajan/FLINT. More benchmarks are provided in media folder.
+
++ Work-Precision Data for Three-Body problem (No events)
+
+<img src="media/wp_sp_tr.jpg" alt="WP-data" width="800"/>
+
++ Three-Body problem propagation (w/ events)
 
-+ Three-Body problem propagation
 ![Julia Results](media/julia_screenshot_CR3BP.PNG)
 
-+ Lorenz equations integration
++ Lorenz equations integration (w/o events)
+
 ![Julia Results](media/julia_screenshot_Lorenz.PNG)
 
 
 ### Installation
 
-FLINT is tested with ifort (18.0.2) compiler from Intel Parallel Studio XE Composer for Windows 2016 integrated with Microsoft Visual Studio Community version 2017. Some testing is also done with MinGW-W64 gfortran (gcc 8.1.0). Doxyfile is provided for generating extensive API documentation using Doxygen. FLINT has no dependency on any other library. The CMakeLists file is provided to generate Visual Studio projects or makefiles on Windows and Linux using cmake. Additionally, it generates cmake config files to easily link FLINT using the find_package() command. The steps to link FLINT in cmake-based projects are:
+FLINT is tested with ifort (2021.11.1) and gfortran (11.4.0) compilers. Doxyfile is provided for generating extensive API documentation using Doxygen. FLINT has no dependency on any other library. The CMakeLists and a FPM (https://github.com/fortran-lang/fpm) toml files are provided to generate build files on Windows and Linux. Additionally, it generates cmake config files to easily link FLINT using the find_package() command when cmake is used. The steps to link FLINT in cmake-based projects are:
 
 + In cmake GUI or command-line, set FLINT_INSTALL_LIB_DIR to the desired directory, where the compiled library, modules, and cmake config files will be installed.
 + In cmake GUI or command-line, set FLINT_INSTALL_BIN_DIR to the desired directory, where the compiled test executables of FLINT will be installed.
@@ -61,7 +75,7 @@ FLINT is tested with ifort (18.0.2) compiler from Intel Parallel Studio XE Compo
 
 ### How to Use
 
-See the test program files, test.f90 and DiffEq.f90, in tests folder for a comparatively complex example problem that uses multiple events. Following are the steps in brief for a simpler problem:
+See the test codes in [tests](https://github.com/princemahajan/FLINT/tree/master/tests) folder. Following are the steps in brief for a simpler problem:
 
 1. Create a differential equation system class by providing differential equation function, events function (if any), and parameters (if any).
 
@@ -188,100 +202,8 @@ See the test program files, test.f90 and DiffEq.f90, in tests folder for a compa
     print *, Eventdata(8,:)     ! Event-ID number
 ```  
 
-For all the FLINT status codes and options supported by Init, Integrate, and Interpolate procedures along with the interfaces for user-supplied functions, see the FLINT_base module in FLINT_base.f90 file. 
-
-### Testing
-
-The following tests were performed using the very first verion of FLINT and they are not updated for the latest version. In the tests, an orbit is propagated for 4 orbital periods and the integration is repeated 5000 times. The tables give the total time for 5000 propagations and all other testing parameters are given for each integration. The relative tolerance is taken as 1e-11 and absolute tolerance as 1e-14 in all the cases except when explicitly mentioned otherwise. Interpolation is used to capture 1000 points per orbit period. Note that the delayed interpolation feature of FLINT is not used in these results.    
-
-The initial conditions in Cartesian coodinates used are as follows:
-- Two-Body circular Earth orbit  (Units: km, sec)
-    + GM: 398600.436233     
-    + Y0 = [6400.0,0.0,0.0, 0.0,5.58037857139,5.58037857139] 
-- Two-Body elliptic Earth orbit  (Units: km, sec)
-    + GM: 398600.436233     
-    + Y0 = [6400.0,0.0,0.0,0.0,7.69202528825512,7.69202528825512]
-- Arenstorf orbit in rotating frame  (Units: nondimensional)
-    + mass-ratio: 0.012277471
-    + time-period: 17.0652165601579625588917206249    
-    + Y0 = [0.994, 0.0, 0.0_wp, -2.00158510637908252240537862224]
-
-Note that JDOP853 and JDDEABM are modern fortran implementations of DOP853 and DDEABM by Jacob Williams, and are available at https://github.com/jacobwilliams. The plots are generated using https://github.com/jacobwilliams/pyplot-fortran.
-
-- Two-Body circular Earth orbit
-
-|Params\Integrator | DOP54   | DOP853 | Verner65E | Verner98R | JDOP853  | JDDEABM |
-|----------------- | ------- | ------ | --------- | --------- | -------  | --------------|
-|Closing Error     | 6.3e-6  | 4.3e-7 |   2.7e-6  |    1.1e-7 |  4.3e-7  |   3.1e-6         |  
-|IOM Error         | 3.1e-11 | 3.1e-11|   3.1e-11 |   3.1e-11 |  3.1e-11 |   3.9e-10        |       
-|Time (s)          |    2.42  | **0.42**   |   1.91    |    0.87   |  1.42    |   **0.62**  |         
-|Func Calls        |  7889   | 1886   |   6043    |    2218   |  1886    |   734           |
-|Accepted Steps    |   1302  | 157    |   749     |    131    |  157     |   365           |
-|Rejected Steps    |     15  | 0      |   7       |    8      |  0       |   NA           |
-
-- Two-Body circular Earth orbit with interpolation
-
-|Params\Integrator | DOP54   | DOP853 | Verner65E | Verner98R | JDOP853  | JDDEABM
-|----------------- | ------- | ------ | --------- | --------- | -------  | --------------
-|Closing Error     | 6.3e-6  | 4.3e-7 |   2.7e-6  |    1.1e-7 |  4.3e-7  |   3.1e-6           
-|IOM Error         | 3.1e-11 | 3.1e-11|   3.1e-11 |   3.1e-11 |  3.1e-11 |   3.9e-10               
-|Time (s)          |    6.3  | **4.3**    |   6.6     |    5.7    |  127.1   |   **5.1**           
-|Func Calls        |  7889   | 2357   |   6792    |    2873   |  2336    |   734           
-|Accepted Steps    |   1302  | 157    |   749     |    131    |  157     |   NA           
-|Rejected Steps    |     15  | 0      |   7       |    8      |  0       |   NA           
-
-- Two-Body elliptic Earth orbit  
-
-|Params\Integrator | DOP54   | DOP853  | Verner65E | Verner98R | JDOP853  | JDDEABM
-|----------------- | ------- | ------  | --------- | --------- | -------  | --------------
-|Closing Error     | 2.7     | 2.7     |   2.7     |    2.7    |  2.7     |   2.7           
-|IOM Error         | 6.7e-10 | 6.7e-10 |   6.7e-10 |   6.7e-10 |  6.7e-10 |   2.3e-9               
-|Time (s)          |    4.8    | **1.3**    |   3.6    |    2.43   |  3.9     |   **3.2**           
-|Func Calls        |  15409  | 6227    |   11520   |    6223   |  6227    |   3325           
-|Accepted Steps    |   2552  | 428     |   1431    |    296    |  428     |   1648           
-|Rejected Steps    |     19  | 99      |   10      |    99     |  99      |   NA         
-
-- Two-Body elliptic Earth orbit with interpolation
-
-|Params\Integrator | DOP54   | DOP853  | Verner65E | Verner98R | JDOP853  | JDDEABM
-|----------------- | ------- | ------  | --------- | --------- | -------  | --------------
-|Closing Error     | 2.7     | 2.7     |   2.7     |    2.7    |  2.7     |   2.7           
-|IOM Error         | 6.7e-10 | 6.7e-10 |   6.7e-10 |   6.7e-10 |  6.7e-10 |   2.3e-9               
-|Time (s)          |    10.7 | **6.0** |   10.5     |    8.6   |  130    |   **8.3**           
-|Func Calls        |  15409  | 7511    |   12951   |    7703   |  7397    |   3325           
-|Accepted Steps    |   2552  | 428     |   1431    |    296    |  428     |   4000           
-|Rejected Steps    |     19  | 99      |   10      |    99      |  99      |   NA         
-
-- Arenstorf orbit
-
-|Params\Integrator | DOP54   | DOP853  | Verner65E | Verner98R | JDOP853  | JDDEABM
-|----------------- | ------- | ------  | --------- | --------- | -------  | --------------
-|Closing Error     | 1.7     | 0.33    |   0.33    |    0.31   |  0.33    |   0.23           
-|IOM Error         | 7.2e-11 | 7.2e-11 |   7.2e-11 |   7.2e-11 |  7.2e-11 |   9.4e-11               
-|Time (s)          |    14.8 | **3.6**    |   10.7   |    6.3    |  12.9   |   **7.9**           
-|Func Calls        |  44847  | 15816   |   32337   |    15344  |  15816   |   8202           
-|Accepted Steps    |   7425  | 1074    |   4006    |    792    |  1074    |   4071           
-|Rejected Steps    |     59  | 266     |   41      |    178    |  266     |   NA           
-
-- Arenstorf orbit with interpolation
-
-|Params\Integrator | DOP54   | DOP853  | Verner65E | Verner98R | JDOP853  | JDDEABM
-|----------------- | ------- | ------  | --------- | --------- | -------  | --------------
-|Closing Error     | 1.7     | 0.33     |   0.33     |    0.31 |  0.33     |   0.23           
-|IOM Error         | 9.4e-11 | 7.1e-11 |   7.1e-11 |   7.1e-11 |  7.2e-11 |   9.4e-11               
-|Time (s)          |    29  | **10.8** |   26.4     |    16.4  |  136    |   **12.4**           
-|Func Calls        |  44847  | 19038   |   36343   |    19304  |  18315   |   9904           
-|Accepted Steps    |   7425  | 1074    |   4006    |    792    |  1074    |   4000           
-|Rejected Steps    |     59  | 266     |   41      |    178    |  266    |   NA           
-
-For relative tol=1e-9 and abs tol=1e-12, the Arenstorf orbit propagated for 4 periods by JDOP853, JDDEABM and FLINT's Verner98R are shown below. Verner98R diverges much slower than DOP853 and DDEABM for this orbit and tolerance values.
-
- ![JDOP853](media/JDOP853.png)
- ![JDDEABM](media/Jacob_DDEABM.png)
- ![FLINT's Verner98R](media/FLINT_Verner98R.png)
-
-The following figure shows the multiple event detection capability of FLINT, in which the X-crossings in decreasing and Y-crossings in increasing direction are detected and reported to the user.
-![FINT's event-detection](media/FLINT_Events.png)
+For all the FLINT status codes and options supported by Init, Integrate, and Interpolate procedures along with the interfaces for user-supplied functions, see the docs and/or [FLINT_base](https://github.com/princemahajan/FLINT/blob/master/src/FLINT_base.f90) module. 
+
 
 ### References
 

fpm.toml

@@ -1,10 +1,10 @@
 name = "FLINT"
 author = "Bharat Mahajan"
-copyright  = "Copyright 2021 Bharat Mahajan"
+copyright  = "Copyright 2024 Bharat Mahajan"
 license = "Apache-2.0"
 description = "Fortran Library for numerical INTegration of differential equations"
 homepage = "https://github.com/princemahajan/FLINT"
-keywords = ["ode", "runge-kutta"]
+keywords = ["ode", "runge-kutta", "ode-events"]
 
 [build]
 auto-executables = false
@@ -17,6 +17,11 @@ source-dir = "src"
 [install]
 library = true
 
+[[test]]
+name = "test_WP"
+source-dir = "tests"
+main = "test_WP.f90"
+
 [[test]]
 name = "test_CR3BP"
 source-dir = "tests"
@@ -26,3 +31,18 @@ main = "test_CR3BP.f90"
 name = "test_Lorenz"
 source-dir = "tests"
 main = "test_Lorenz.f90"
+
+[[test]]
+name = "test_TBP"
+source-dir = "tests"
+main = "test_TBP.f90"
+
+#[dev-dependencies]
+#rklib = { git="https://github.com/jacobwilliams/rklib.git" }
+#rklib = { path = "..\\rklib" }
+
+#[[test]]
+#name = "comp_TBP_rklib"
+#source-dir = "tests"
+#main = "comp_TBP_rklib.f90"
+

src/ButcherTableaus.f90

@@ -420,8 +420,12 @@ module ButcherTableaus
 
 !> Coefficients for 7th-order interpolation
 
+! TBD: this explicit definition is workaround as gfortran is 
+! giving error for including type-spec in array constructor
+integer :: i 
+
 !> Non-zero coefficients for non-zero stages
-integer, parameter :: DOP853_di_NZ(*) = [1,(i,integer:: i=6,16)]
+integer, parameter :: DOP853_di_NZ(*) = [ 1, (i, i=6,16)]
 integer, parameter :: DOP853_dj_NZ(*) = [4,5,6,7]
 
 !> dij, where i = 1 to s, j = 1 to pstar for a generic method
@@ -824,7 +828,7 @@ module ButcherTableaus
 !> Coefficients for 8th-order interpolation
 
 !> Non-zero coefficients for non-zero stages
-integer, parameter :: Verner98R_di_NZ(*) = [1,(i, integer:: i=8,15),(i, integer:: i=17,21)]
+integer, parameter :: Verner98R_di_NZ(*) = [1,(i, i=8,15),(i, i=17,21)]
 
 !> dij, where i = 1 to s, j = 1 to pstar for a generic method
 real(WP), parameter :: Verner98R_d(*,1:*) = reshape(&
@@ -1084,7 +1088,7 @@ module ButcherTableaus
 !> Coefficients for 5th-order interpolation
 
 !> Non-zero coefficients for non-zero stages
-integer, parameter :: Verner65E_di_NZ(*) = [1,(i, integer:: i=4,10)]
+integer, parameter :: Verner65E_di_NZ(*) = [1,(i, i=4,10)]
 
 !> dij, where i = 1 to s, j = 1 to pstar for a generic method
 real(WP), parameter :: Verner65E_d(*,1:*) = reshape(&

src/ERK.f90

@@ -222,10 +222,11 @@ end subroutine erk_destroy
 
 
 
-    module subroutine erk_stepint(me, X0, Y0, h, Y1, Yint12, FCalls, EstimateErr, Err, params)
+    module subroutine erk_stepint(me, X0, Y0, F0, h, Y1, Yint12, FCalls, EstimateErr, Err, params)
         class(ERK_class), intent(inout) :: me
         real(WP), intent(in) :: X0
         real(WP), dimension(me%pDiffEqSys%n), intent(in) :: Y0
+        real(WP), dimension(me%pDiffEqSys%n), intent(inout) :: F0
         real(WP), intent(in) :: h
         real(WP), dimension(size(Y0)), intent(out) :: Y1, Yint12
         integer, intent(out) :: FCalls