Updates to ODESolver and sixDOFODE classes

Updates enable automatic handling of multiple calls to ODESolve::solve within a
single time step

src/ODE/ODE/ODE.H

@@ -91,6 +91,15 @@ public:
 
         //- Update ODE after the solution, advancing by delta
         virtual void update(const scalar delta) = 0;
+
+        //- Member function for initialising the ODE state. Required for some
+        //  ODE's (specifically sixDOFODE's) which might be called multiple
+        //  times per single time-step. Called at the beginning of
+        //  ODESolver::solve and does nothing by default
+        virtual void init()
+        {
+            // Does nothing by default
+        }
 };
 
 

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -52,8 +52,11 @@ void Foam::ODESolver::solve
     const scalar xEnd,
     const scalar eps,
     scalar& hEst
-) const
+)
 {
+    // Initialize ODE before solving
+    ode_.init();
+
     const label MAXSTP = 10000;
 
     scalar x = xStart;

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -53,7 +53,7 @@ class ODESolver
 
 protected:
 
-    // Private data
+    // Protected data
 
         //- Reference to ODE
         ODE& ode_;
@@ -62,6 +62,8 @@ protected:
         mutable scalarField dydx_;
 
 
+private:
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -126,13 +128,13 @@ public:
         ) const = 0;
 
 
-        virtual void solve
+        void solve
         (
             const scalar xStart,
             const scalar xEnd,
             const scalar eps,
             scalar& hEst
-        ) const;
+        );
 };
 
 

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -233,6 +233,39 @@ Foam::vector Foam::geometricSixDOF::dexpMap
 }
 
 
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to geometricSixDOF&
+    const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = gsd.Xrel_;
+    U_ = gsd.U_;
+    Uaverage_ = gsd.Uaverage_;
+    rotation_ = gsd.rotation_;
+    omega_ = gsd.omega_;
+    omegaAverage_ = gsd.omegaAverage_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = gsd.coeffs_;
+
+    fixedSurge_ = gsd.fixedSurge_;
+    fixedSway_ = gsd.fixedSway_;
+    fixedHeave_ = gsd.fixedHeave_;
+    fixedRoll_ = gsd.fixedRoll_;
+    fixedPitch_ = gsd.fixedPitch_;
+    fixedYaw_ = gsd.fixedYaw_;
+    referentMotionConstraints_ = gsd.referentMotionConstraints_;
+    referentRotation_ = gsd.referentRotation_;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
@@ -310,18 +343,7 @@ Foam::geometricSixDOF::geometricSixDOF
     const geometricSixDOF& gsd
 )
 :
-    sixDOFODE
+    sixDOFODE(name, gsd),
-    (
-        IOobject
-        (
-            name,
-            gsd.dict().instance(),
-            gsd.dict().local(),
-            gsd.dict().db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
 
     Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
     U_(gsd.U_.name(), gsd.U_),
@@ -399,37 +421,6 @@ const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
 }
 
 
-void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
-{
-    // First set the state in base class subobject
-    sixDOFODE::setState(sd);
-
-    // Cast sixDOFODE& to geometricSixDOF&
-    const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
-
-    // Set state variables for this class
-    Xrel_ = gsd.Xrel_;
-    U_ = gsd.U_;
-    Uaverage_ = gsd.Uaverage_;
-    rotation_ = gsd.rotation_;
-    omega_ = gsd.omega_;
-    omegaAverage_ = gsd.omegaAverage_;
-
-    // Copy ODE coefficients: this carries actual ODE state
-    // HJ, 23/Mar/2015
-    coeffs_ = gsd.coeffs_;
-
-    fixedSurge_ = gsd.fixedSurge_;
-    fixedSway_ = gsd.fixedSway_;
-    fixedHeave_ = gsd.fixedHeave_;
-    fixedRoll_ = gsd.fixedRoll_;
-    fixedPitch_ = gsd.fixedPitch_;
-    fixedYaw_ = gsd.fixedYaw_;
-    referentMotionConstraints_ = gsd.referentMotionConstraints_;
-    referentRotation_ = gsd.referentRotation_;
-}
-
-
 const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
 {
     return omegaAverage_;

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H

@@ -182,6 +182,16 @@ class geometricSixDOF
             vector dexpMap(const vector& rotInc, const vector& omega) const;
 
 
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
 public:
 
     // Run-time type information
@@ -246,12 +256,6 @@ public:
                 virtual const dimensionedVector& Uaverage() const;
 
 
-            // Non-access control for setting state variables
-
-                //- Set ODE parameters from another ODE
-                virtual void setState(const sixDOFODE&);
-
-
             // Average motion per time-step
 
                 //- Return average rotational velocity in relative coordinate

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -190,6 +190,40 @@ void Foam::quaternionSixDOF::constrainTranslation(vector& vec) const
 }
 
 
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to quaternionSixDOF&
+    const quaternionSixDOF& qsd = refCast<const quaternionSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = qsd.Xrel_;
+    U_ = qsd.U_;
+    Uaverage_ = qsd.Uaverage_;
+    rotation_ = qsd.rotation_;
+    omega_ = qsd.omega_;
+    omegaAverage_ = qsd.omegaAverage_;
+    omegaAverageAbsolute_ = qsd.omegaAverageAbsolute_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = qsd.coeffs_;
+
+    fixedSurge_ = qsd.fixedSurge_;
+    fixedSway_ = qsd.fixedSway_;
+    fixedHeave_ = qsd.fixedHeave_;
+    fixedRoll_ = qsd.fixedRoll_;
+    fixedPitch_ = qsd.fixedPitch_;
+    fixedYaw_ = qsd.fixedYaw_;
+    referentMotionConstraints_ = qsd.referentMotionConstraints_;
+    referentRotation_ = qsd.referentRotation_;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
@@ -260,18 +294,7 @@ Foam::quaternionSixDOF::quaternionSixDOF
     const quaternionSixDOF& qsd
 )
 :
-    sixDOFODE
+    sixDOFODE(name, qsd),
-    (
-        IOobject
-        (
-            name,
-            qsd.dict().instance(),
-            qsd.dict().local(),
-            qsd.dict().db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
 
     Xrel_(qsd.Xrel_.name(), qsd.Xrel_),
     U_(qsd.U_.name(), qsd.U_),
@@ -354,38 +377,6 @@ const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const
 }
 
 
-void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
-{
-    // First set the state in base class subobject
-    sixDOFODE::setState(sd);
-
-    // Cast sixDOFODE& to quaternionSixDOF&
-    const quaternionSixDOF& qsd = refCast<const quaternionSixDOF>(sd);
-
-    // Set state variables for this class
-    Xrel_ = qsd.Xrel_;
-    U_ = qsd.U_;
-    Uaverage_ = qsd.Uaverage_;
-    rotation_ = qsd.rotation_;
-    omega_ = qsd.omega_;
-    omegaAverage_ = qsd.omegaAverage_;
-    omegaAverageAbsolute_ = qsd.omegaAverageAbsolute_;
-
-    // Copy ODE coefficients: this carries actual ODE state
-    // HJ, 23/Mar/2015
-    coeffs_ = qsd.coeffs_;
-
-    fixedSurge_ = qsd.fixedSurge_;
-    fixedSway_ = qsd.fixedSway_;
-    fixedHeave_ = qsd.fixedHeave_;
-    fixedRoll_ = qsd.fixedRoll_;
-    fixedPitch_ = qsd.fixedPitch_;
-    fixedYaw_ = qsd.fixedYaw_;
-    referentMotionConstraints_ = qsd.referentMotionConstraints_;
-    referentRotation_ = qsd.referentRotation_;
-}
-
-
 const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
 {
     return omegaAverage_;

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H

@@ -158,6 +158,16 @@ class quaternionSixDOF
             void constrainTranslation(vector& vec) const;
 
 
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
 public:
 
     // Run-time type information
@@ -211,12 +221,6 @@ public:
                 virtual const dimensionedVector& Uaverage() const;
 
 
-            // Non-access control for setting state variables
-
-                //- Set ODE parameters from another ODE
-                virtual void setState(const sixDOFODE&);
-
-
             // Average motion per time-step
 
                 //- Return average rotational velocity in relative coordinate

src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C

@@ -27,7 +27,7 @@ Class
 
 Description
     Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
-    equations
+    equations with necessary interface for two-way coupling with CFD solvers.
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -26,7 +26,7 @@ Class
 
 Description
     Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
-    equations
+    equations with necessary interface for two-way coupling with CFD solvers.
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
@@ -119,6 +119,27 @@ void Foam::sixDOFODE::aitkensRelaxation
 }
 
 
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::sixDOFODE::setState(const sixDOFODE& sd)
+{
+    // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_.
+    // In case of multiple updates, overwriting Aitkens relaxation parameters
+    // would invalidate the underrelaxation.  IG, 5/May/2016
+    mass_ = sd.mass_;
+    momentOfInertia_ = sd.momentOfInertia_;
+
+    Xequilibrium_ = sd.Xequilibrium_;
+    linSpringCoeffs_ = sd.linSpringCoeffs_;
+    linDampingCoeffs_ = sd.linDampingCoeffs_;
+
+    force_ = sd.force_;
+    moment_ = sd.moment_;
+    forceRelative_ = sd.forceRelative_;
+    momentRelative_ = sd.momentRelative_;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sixDOFODE::sixDOFODE(const IOobject& io)
@@ -146,7 +167,42 @@ Foam::sixDOFODE::sixDOFODE(const IOobject& io)
     force_(dict_.lookup("force")),
     moment_(dict_.lookup("moment")),
     forceRelative_(dict_.lookup("forceRelative")),
-    momentRelative_(dict_.lookup("momentRelative"))
+    momentRelative_(dict_.lookup("momentRelative")),
+
+    curTimeIndex_(-1),
+    oldStatePtr_()
+{}
+
+
+Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd)
+:
+    ODE(),
+    dict_(sd.dict_),
+
+    mass_(sd.mass_),
+    momentOfInertia_(sd.momentOfInertia_),
+
+    Xequilibrium_(sd.Xequilibrium_),
+    linSpringCoeffs_(sd.linSpringCoeffs_),
+    linDampingCoeffs_(sd.linDampingCoeffs_),
+
+    relaxFactorT_(1.0),
+    relaxFactorR_(1.0),
+    oldRelaxFactorT_(1.0),
+    oldRelaxFactorR_(1.0),
+
+    A_(3, vector::zero),
+    OmegaDot_(3, vector::zero),
+    An_(3, vector::zero),
+    OmegaDotn_(3, vector::zero),
+
+    force_(sd.force_),
+    moment_(sd.moment_),
+    forceRelative_(sd.forceRelative_),
+    momentRelative_(sd.momentRelative_),
+
+    curTimeIndex_(-1),
+    oldStatePtr_()
 {}
 
 
@@ -216,23 +272,30 @@ void Foam::sixDOFODE::relaxAcceleration
 }
 
 
-void Foam::sixDOFODE::setState(const sixDOFODE& sd)
+void Foam::sixDOFODE::init()
 {
-    // Set state does not copy AList_, AOld_, relaxFactor_ and
+    // Get time index
-    // relaxFactorOld_. In case of multiple updates, overwriting Aitkens
+    const label timeIndex = dict().time().timeIndex();
-    // relaxation parameters would invalidate the underrelaxation.
-    // IG, 5/May/2016
-    mass_ = sd.mass_;
-    momentOfInertia_ = sd.momentOfInertia_;
 
-    Xequilibrium_ = sd.Xequilibrium_;
+    if (curTimeIndex_ < timeIndex)
-    linSpringCoeffs_ = sd.linSpringCoeffs_;
+    {
-    linDampingCoeffs_ = sd.linDampingCoeffs_;
+        // First call in this time index, store data
+        oldStatePtr_ = this->clone(dict().name() + "_0");
 
-    force_ = sd.force_;
+        Info<< "First 6DOF solution within a time step, storing old data..."
-    moment_ = sd.moment_;
+            << endl;
-    forceRelative_ = sd.forceRelative_;
+    }
-    momentRelative_ = sd.momentRelative_;
+    else
+    {
+        // Multiple calls in this time index, reset this data
+        this->setState(oldStatePtr_());
+
+        Info<< "Repeated 6DOF solution within a time step, restoring data..."
+            << endl;
+    }
+
+    // Update local time index
+    curTimeIndex_ = timeIndex;
 }
 
 

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -26,8 +26,7 @@ Class
 
 Description
     Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
-    equations with a fat interface in order to provide representation of
+    equations with necessary interface for two-way coupling with CFD solvers.
-    rotations in quaternions or rotation tensors.
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
@@ -55,9 +54,6 @@ SourceFiles
 namespace Foam
 {
 
-class finiteRotation;
-class HamiltonRodriguezRot;
-
 /*---------------------------------------------------------------------------*\
                          Class sixDOFODE Declaration
 \*---------------------------------------------------------------------------*/
@@ -118,7 +114,7 @@ class sixDOFODE
             List<vector> OmegaDotn_;
 
 
-        // External forces
+        // External forces and moments
 
             //- Force driving the motion in inertial coord. sys.
             dimensionedVector force_;
@@ -133,6 +129,20 @@ class sixDOFODE
             dimensionedVector momentRelative_;
 
 
+        // Private data used to control multiple ODE updates per time step
+        // Note: Before solving the ODE from the top level, we will store the
+        // previous state if this is the first update in a given time step. The
+        // state will be kept as the pointer to sixDOFODE, which is not optimal
+        // because we are keeping track of some unnecessary data. We could use
+        // more elegant and efficient code design. VV, 1/Mar/2017.
+
+            //- Local time index
+            label curTimeIndex_;
+
+            //- Pointer to the sixDOFODE object carrying old state
+            autoPtr<sixDOFODE> oldStatePtr_;
+
+
     // Private Member Functions
 
         // Copy control
@@ -150,6 +160,16 @@ class sixDOFODE
             void aitkensRelaxation(const scalar min, const scalar max);
 
 
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
 public:
 
     //- Run-time type information
@@ -173,6 +193,9 @@ public:
         //- Construct from dictionary
         sixDOFODE(const IOobject& io);
 
+        //- Copy construct given new name
+        sixDOFODE(const word& name, const sixDOFODE& sd);
+
         //- Return a clone, changing the name
         virtual autoPtr<sixDOFODE> clone(const word& name) const = 0;
 
@@ -284,12 +307,6 @@ public:
                 virtual const dimensionedVector& Uaverage() const = 0;
 
 
-            // Non-access control for setting state variables
-
-                //- Set ODE parameters from another ODE
-                virtual void setState(const sixDOFODE&);
-
-
             // Average motion per time-step
 
                 //- Return average rotational velocity in relative coordinate
@@ -341,6 +358,10 @@ public:
                 //- Update ODE after the solution, advancing by delta
                 virtual void update(const scalar delta) = 0;
 
+                //- Initialise ODE before solving, handling multiple calls per
+                //  single time step
+                virtual void init();
+
 
         // Write control
 

src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H

@@ -26,7 +26,7 @@ Class
 
 Description
     Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
-    equations
+    equations with necessary interface for two-way coupling with CFD solvers.
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.