Merge pull request #1 from DominicDirkx/AbamIntegrator

Abam/Bs integrator

Tudat/Astrodynamics/BasicAstrodynamics/convertMeanToEccentricAnomalies.h

@@ -271,7 +271,6 @@ ScalarType convertMeanAnomalyToEccentricAnomaly(
             eccentricAnomaly = bisectionRootfinder->execute( rootFunction, initialGuess );
         }
     }
-
     //  Eccentricity is invalid: eccentricity < 0.0 or eccentricity >= 1.0.
     else
     {

Tudat/Astrodynamics/Propagators/UnitTests/unitTestForwardsBackwardsIntegration.cpp

@@ -69,7 +69,7 @@ boost::shared_ptr< IntegratorSettings< TimeType > > getIntegrationSettings(
     }
     else if( integratorCase < 5 )
     {
-        RungeKuttaCoefficients::CoefficientSets coefficientSet;
+        RungeKuttaCoefficients::CoefficientSets coefficientSet = RungeKuttaCoefficients::undefinedCoefficientSet;
         if( integratorCase == 1 )
         {
             coefficientSet = RungeKuttaCoefficients::rungeKuttaFehlberg45;

Tudat/Astrodynamics/Propagators/dynamicsStateDerivativeModel.h

@@ -68,7 +68,8 @@ class DynamicsStateDerivativeModel
                 const std::vector< IntegratedStateType > ) > environmentUpdateFunction,
             const boost::shared_ptr< VariationalEquations > variationalEquations =
             boost::shared_ptr< VariationalEquations >( ) ):
-        environmentUpdateFunction_( environmentUpdateFunction ), variationalEquations_( variationalEquations )
+        environmentUpdateFunction_( environmentUpdateFunction ), variationalEquations_( variationalEquations ),
+        functionEvaluationCounter_( 0 )
     {
         std::vector< IntegratedStateType > stateTypeList;
         totalStateSize_ = 0;
@@ -219,6 +220,8 @@ class DynamicsStateDerivativeModel
                         stateDerivative_.block( 0, 0, totalStateSize_, variationalEquations_->getNumberOfParameterValues( ) )  );
         }
 
+        functionEvaluationCounter_++;
+
         return stateDerivative_;
     }
 
@@ -446,6 +449,28 @@ class DynamicsStateDerivativeModel
         return stateTypeStartIndex_;
     }
 
+    //! Function to retrieve number of calls to the computeStateDerivative function
+    /*!
+     * Function to retrieve number of calls to the computeStateDerivative function since object creation/last call to
+     * resetFunctionEvaluationCounter function
+     * \return Number of calls to the computeStateDerivative function since object creation/last call to
+     * resetFunctionEvaluationCounter function
+     */
+    int getNumberOfFunctionEvaluations( )
+    {
+        return functionEvaluationCounter_;
+    }
+
+    //! Function to resetr the number of calls to the computeStateDerivative function to zero.
+    /*!
+     * Function to resetr the number of calls to the computeStateDerivative function to zero.  Typically called before any
+     * start of numerical integration of dynamics (automatically by DynamicsSimulator)
+     */
+    void resetFunctionEvaluationCounter( )
+    {
+        functionEvaluationCounter_ = 0;
+    }
+
 private:
 
     //! Function to convert the to the conventional form in the global frame per dynamics type.
@@ -490,19 +515,13 @@ class DynamicsStateDerivativeModel
                 // Get state block indices of current state derivative model
                 currentIndices = stateIndices_.at( stateDerivativeModelsIterator_->first ).at( i );
 
-//                std::cout << "Pre-converted state: " << state.block( currentIndices.first, startColumn, currentIndices.second, 1 ).transpose( ) << std::endl;
-
                 // Set current block in split state (in global form)
                 stateDerivativeModelsIterator_->second.at( i )->convertCurrentStateToGlobalRepresentation(
                             state.block( currentIndices.first, startColumn, currentIndices.second, 1 ), time,
                             currentStatesPerTypeInConventionalRepresentation_.at(
                                 stateDerivativeModelsIterator_->first ).block(
                                 currentStateTypeSize, 0, currentIndices.second, 1 ) );
 
-//                std::cout << "Converted state: " << currentStatesPerTypeInConventionalRepresentation_.at(
-//                               stateDerivativeModelsIterator_->first ).block(
-//                               currentStateTypeSize, 0, currentIndices.second, 1 ).transpose( ) << std::endl;
-
                 currentStateTypeSize += currentIndices.second;
             }
         }
@@ -558,6 +577,9 @@ class DynamicsStateDerivativeModel
     //! convertCurrentStateToGlobalRepresentationPerType
     std::unordered_map< IntegratedStateType, Eigen::Matrix< StateScalarType, Eigen::Dynamic, 1 > >
     currentStatesPerTypeInConventionalRepresentation_;
+
+    //! Variable to keep track of the number of calls to the computeStateDerivative function
+    int functionEvaluationCounter_ = 0;
 };
 
 //! Function to retrieve a single given acceleration model from a list of models

Tudat/JsonInterface/Mathematics/integrator.h

@@ -27,7 +27,8 @@ static std::map< AvailableIntegrators, std::string > integratorTypes =
     { rungeKutta4, "rungeKutta4" },
     { euler, "euler" },
     { rungeKuttaVariableStepSize, "rungeKuttaVariableStepSize" },    
-    { adamsBashforthMoulton, "adamsBashforthMoulton" }
+    { adamsBashforthMoulton, "adamsBashforthMoulton" },
+    { bulirschStoer, "bulirschStoer" },
 };
 
 //! `AvailableIntegrators` not supported by `json_interface`.
@@ -127,9 +128,32 @@ void to_json( nlohmann::json& jsonObject, const boost::shared_ptr< IntegratorSet
         jsonObject[ K::maximumStepSize ] = adamsBashforthMoultonSettings->maximumStepSize_;
         jsonObject[ K::relativeErrorTolerance ] = adamsBashforthMoultonSettings->relativeErrorTolerance_;
         jsonObject[ K::absoluteErrorTolerance ] = adamsBashforthMoultonSettings->absoluteErrorTolerance_;
+        jsonObject[ K::minimumOrder ] = adamsBashforthMoultonSettings->maximumOrder_;
+        jsonObject[ K::maximumOrder ] = adamsBashforthMoultonSettings->maximumOrder_;
         jsonObject[ K::bandwidth ] = adamsBashforthMoultonSettings->bandwidth_;
         return;
     }
+    case bulirschStoer:
+    {
+        boost::shared_ptr< BulirschStoerIntegratorSettings< TimeType > > bulirschStoerSettings =
+                boost::dynamic_pointer_cast< BulirschStoerIntegratorSettings< TimeType > >( integratorSettings );
+        assertNonNullPointer( bulirschStoerSettings );
+        jsonObject[ K::initialStepSize ] = bulirschStoerSettings->initialTimeStep_;
+        jsonObject[ K::minimumStepSize ] = bulirschStoerSettings->minimumStepSize_;
+        jsonObject[ K::maximumStepSize ] = bulirschStoerSettings->maximumStepSize_;
+        jsonObject[ K::relativeErrorTolerance ] = bulirschStoerSettings->relativeErrorTolerance_;
+        jsonObject[ K::absoluteErrorTolerance ] = bulirschStoerSettings->absoluteErrorTolerance_;
+        jsonObject[ K::extrapolationSequence ] =  bulirschStoerSettings->extrapolationSequence_;
+        jsonObject[ K::maximumNumberOfSteps ] =  bulirschStoerSettings->maximumNumberOfSteps_;
+        jsonObject[ K::safetyFactorForNextStepSize ] =
+                bulirschStoerSettings->safetyFactorForNextStepSize_;
+        jsonObject[ K::maximumFactorIncreaseForNextStepSize ] =
+                bulirschStoerSettings->maximumFactorIncreaseForNextStepSize_;
+        jsonObject[ K::minimumFactorDecreaseForNextStepSize ] =
+                bulirschStoerSettings->minimumFactorDecreaseForNextStepSize_;
+
+        return;
+    }
     default:
         handleUnimplementedEnumValue( integratorType, integratorTypes, unsupportedIntegratorTypes );
     }
@@ -192,23 +216,50 @@ void from_json( const nlohmann::json& jsonObject, boost::shared_ptr< IntegratorS
     case adamsBashforthMoulton:
     {
         AdamsBashforthMoultonSettings< TimeType > defaults(
-                    integratorType, 0.0, 0.0, 0.0, 0.0 );
+                    0.0, 0.0, 0.0, 0.0 );
 
         integratorSettings = boost::make_shared< AdamsBashforthMoultonSettings< TimeType > >(
-                    integratorType,
                     initialTime,
                     getValue< TimeType >( jsonObject, K::initialStepSize ),
                     getValue< TimeType >( jsonObject, K::minimumStepSize ),
                     getValue< TimeType >( jsonObject, K::maximumStepSize ),
                     getValue( jsonObject, K::relativeErrorTolerance, defaults.relativeErrorTolerance_ ),
                     getValue( jsonObject, K::absoluteErrorTolerance, defaults.absoluteErrorTolerance_ ),
+                    getValue( jsonObject, K::minimumOrder, defaults.minimumOrder_ ),
+                    getValue( jsonObject, K::maximumOrder, defaults.maximumOrder_ ),
                     getValue( jsonObject, K::saveFrequency, defaults.saveFrequency_ ),
                     getValue( jsonObject, K::assessPropagationTerminationConditionDuringIntegrationSubsteps,
                               defaults.assessPropagationTerminationConditionDuringIntegrationSubsteps_ ),
                     getValue( jsonObject, K::bandwidth,
                               defaults.bandwidth_ ) );
         return;
     }
+    case bulirschStoer:
+    {
+        BulirschStoerIntegratorSettings< TimeType > defaults(
+                    0.0, 0.0, bulirsch_stoer_sequence, 6, std::numeric_limits< double >::epsilon( ),
+                    std::numeric_limits< double >::infinity( ) );
+
+        integratorSettings = boost::make_shared< BulirschStoerIntegratorSettings< TimeType > >(
+                    initialTime,
+                    getValue< TimeType >( jsonObject, K::initialStepSize ),
+                    getValue( jsonObject, K::extrapolationSequence, defaults.extrapolationSequence_ ),
+                    getValue( jsonObject, K::maximumNumberOfSteps, defaults.maximumNumberOfSteps_ ),
+                    getValue< TimeType >( jsonObject, K::minimumStepSize ),
+                    getValue< TimeType >( jsonObject, K::maximumStepSize ),
+                    getValue( jsonObject, K::relativeErrorTolerance, defaults.relativeErrorTolerance_ ),
+                    getValue( jsonObject, K::absoluteErrorTolerance, defaults.absoluteErrorTolerance_ ),
+                    getValue( jsonObject, K::saveFrequency, defaults.saveFrequency_ ),
+                    getValue( jsonObject, K::assessPropagationTerminationConditionDuringIntegrationSubsteps,
+                              defaults.assessPropagationTerminationConditionDuringIntegrationSubsteps_ ),
+                    getValue( jsonObject, K::safetyFactorForNextStepSize,
+                              defaults.safetyFactorForNextStepSize_ ),
+                    getValue( jsonObject, K::maximumFactorIncreaseForNextStepSize,
+                              defaults.maximumFactorIncreaseForNextStepSize_ ),
+                    getValue( jsonObject, K::minimumFactorDecreaseForNextStepSize,
+                              defaults.minimumFactorDecreaseForNextStepSize_ ) );
+        return;
+    }
     default:
         handleUnimplementedEnumValue( integratorType, integratorTypes, unsupportedIntegratorTypes );
     }

Tudat/JsonInterface/Support/keys.cpp

@@ -330,7 +330,10 @@ const std::string Keys::Integrator::safetyFactorForNextStepSize = "safetyFactorF
 const std::string Keys::Integrator::maximumFactorIncreaseForNextStepSize = "maximumFactorIncreaseForNextStepSize";
 const std::string Keys::Integrator::minimumFactorDecreaseForNextStepSize = "minimumFactorDecreaseForNextStepSize";
 const std::string Keys::Integrator::bandwidth = "bandwidth";
-
+const std::string Keys::Integrator::extrapolationSequence = "extrapolationSequence";
+const std::string Keys::Integrator::maximumNumberOfSteps = "maximumNumberOfSteps";
+const std::string Keys::Integrator::maximumOrder = "maximumOrder";
+const std::string Keys::Integrator::minimumOrder = "minimumOrder";
 
 //  Interpolation
 

Tudat/JsonInterface/Support/keys.h

@@ -351,6 +351,12 @@ struct Keys
         static const std::string maximumFactorIncreaseForNextStepSize;
         static const std::string minimumFactorDecreaseForNextStepSize;
         static const std::string bandwidth;
+        static const std::string extrapolationSequence;
+        static const std::string maximumNumberOfSteps;
+        static const std::string minimumOrder;
+        static const std::string maximumOrder;
+
+
     };
 
     struct Interpolation

Tudat/JsonInterface/UnitTests/inputs/unitTestIntegrator/bulirschStoer.json

@@ -0,0 +1,10 @@
+{
+  "initialStepSize": 1.4,
+  "minimumStepSize": 0.4,
+  "maximumStepSize": 2.4,
+  "maximumNumberOfSteps": 8,
+  "relativeErrorTolerance": 0.0001,
+  "absoluteErrorTolerance": 0.01,
+  "type": "bulirschStoer",
+  "initialTime": -0.3
+}

Tudat/JsonInterface/UnitTests/inputs/unitTestIntegrator/types.json

@@ -2,5 +2,6 @@
   "euler",
   "rungeKutta4",
   "rungeKuttaVariableStepSize",
-  "adamsBashforthMoulton"
+  "adamsBashforthMoulton",
+  "bulirschStoer"
 ]

Tudat/JsonInterface/UnitTests/unitTestIntegrator.cpp

@@ -152,17 +152,48 @@ BOOST_AUTO_TEST_CASE( test_json_integrator_adamsBashforthMoulton )
     const double relativeErrorTolerance = 1.0E-4;
     const double absoluteErrorTolerance = 1.0E-2;
     const double bandwidth = 200;
+    const int minimumOrder = 6;
+    const int maximumOrder = 11;
     const boost::shared_ptr< IntegratorSettings< double > > manualSettings =
-            boost::make_shared< AdamsBashforthMoultonSettings< double > >( integratorType,
-                                                                                initialTime,
-                                                                                initialStepSize,
-                                                                                minimumStepSize,
-                                                                                maximumStepSize,
-                                                                                relativeErrorTolerance,
-                                                                                absoluteErrorTolerance,
-                                                                                1,
-                                                                                false,
-                                                                                bandwidth );
+            boost::make_shared< AdamsBashforthMoultonSettings< double > >(  initialTime,
+                                                                initialStepSize,
+                                                                minimumStepSize,
+                                                                maximumStepSize,
+                                                                relativeErrorTolerance,
+                                                                absoluteErrorTolerance,
+                                                                minimumOrder,
+                                                                maximumOrder,
+                                                                1,
+                                                                false,
+                                                                bandwidth );
+
+    // Compare
+    BOOST_CHECK_EQUAL_JSON( fromFileSettings, manualSettings );
+}
+
+// Test 6: bulirschStoer
+BOOST_AUTO_TEST_CASE( test_json_integrator_bulirschStoer )
+{
+    using namespace tudat::numerical_integrators;
+    using namespace tudat::json_interface;
+
+    // Create IntegratorSettings from JSON file
+    const boost::shared_ptr< IntegratorSettings< double > > fromFileSettings =
+            parseJSONFile< boost::shared_ptr< IntegratorSettings< double > > >( INPUT( "bulirschStoer" ) );
+
+    // Create IntegratorSettings manually
+    const double initialTime = -0.3;
+    const double initialStepSize = 1.4;
+    const double minimumStepSize = 0.4;
+    const double maximumStepSize = 2.4;
+    const double relativeErrorTolerance = 1.0E-4;
+    const double absoluteErrorTolerance = 1.0E-2;
+    const int maximumNumberOfSteps = 8;
+
+    const boost::shared_ptr< IntegratorSettings< double > > manualSettings =
+            boost::make_shared< BulirschStoerIntegratorSettings< double > >(
+                initialTime, initialStepSize, bulirsch_stoer_sequence, maximumNumberOfSteps,
+                minimumStepSize, maximumStepSize, relativeErrorTolerance, absoluteErrorTolerance );
 
     // Compare
     BOOST_CHECK_EQUAL_JSON( fromFileSettings, manualSettings );

Tudat/JsonInterface/UnitTests/unitTestSupport.h

@@ -82,8 +82,6 @@ void checkConsistentEnum( const std::string& filename,
         }
     }
 
-    std::cout << "JSON file: " << filename << std::endl;
-
     // Check that values and supportedValues are equivalent
     const std::vector< Enum > values = parseJSONFile< std::vector< Enum > >( filename );
     BOOST_CHECK( containsAllOf( values, supportedValues ) && containsAllOf( supportedValues, values ) );

Tudat/Mathematics/NumericalIntegrators/CMakeLists.txt

@@ -11,13 +11,15 @@
 # Add source files.
 set(NUMERICALINTEGRATORS_SOURCES
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/rungeKuttaCoefficients.cpp"
+  "${SRCROOT}${NUMERICALINTEGRATORSDIR}/bulirschStoerVariableStepsizeIntegrator.cpp"
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/UnitTests/numericalIntegratorTests.cpp"
 )
 
 # Add header files.
 set(NUMERICALINTEGRATORS_HEADERS 
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/adamsBashforthMoultonIntegrator.h"
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/createNumericalIntegrator.h"
+  "${SRCROOT}${NUMERICALINTEGRATORSDIR}/bulirschStoerVariableStepsizeIntegrator.h"
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/euler.h"
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/numericalIntegrator.h"
   "${SRCROOT}${NUMERICALINTEGRATORSDIR}/reinitializableNumericalIntegrator.h"
@@ -74,3 +76,8 @@ target_link_libraries(test_RungeKuttaFehlberg78Integrator tudat_numerical_integr
 add_executable(test_RungeKutta87DormandPrinceIntegrator "${SRCROOT}${NUMERICALINTEGRATORSDIR}/UnitTests/unitTestRungeKutta87DormandPrinceIntegrator.cpp")
 setup_custom_test_program(test_RungeKutta87DormandPrinceIntegrator "${SRCROOT}${NUMERICALINTEGRATORSDIR}")
 target_link_libraries(test_RungeKutta87DormandPrinceIntegrator tudat_numerical_integrators tudat_input_output ${Boost_LIBRARIES})
+
+add_executable(test_BulirschStoerVariableStepSizeIntegrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators/UnitTests/unitTestBulirschStoerVariableStepSizeIntegrator.cpp")
+setup_custom_test_program(test_BulirschStoerVariableStepSizeIntegrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators")
+target_link_libraries(test_BulirschStoerVariableStepSizeIntegrator tudat_numerical_integrators tudat_input_output ${Boost_LIBRARIES})
+

Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h

@@ -37,6 +37,15 @@ static inline Eigen::VectorXd computeZeroStateDerivative( const double time,
     return Eigen::VectorXd::Zero( state.rows( ) );
 }
 
+static inline Eigen::VectorXd computeConstantStateDerivative( const double time,
+                                                              const Eigen::VectorXd& state )
+{
+    Eigen::VectorXd stateDerivative = Eigen::VectorXd::Zero( state.rows( ) );
+    stateDerivative( 0 ) = 1.0;
+    //std::cout<<"State der: "<<time<<" "<<stateDerivative.transpose( )<<" "<<std::endl<<state.transpose( )<<std::endl;
+    return stateDerivative;
+}
+
 //! Compute analytical state for zero-state derivative.
 /*!
  * Computes analytical state for zero-state derivative function. This function always returns the