Updates to Aitkens relaxation

Now completely handled from the sixDOFODE class

src/ODE/ODE/ODE.H

@@ -91,15 +91,6 @@ 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

@@ -54,9 +54,6 @@ void Foam::ODESolver::solve
     scalar& hEst
 )
 {
-    // Initialize ODE before solving
-    ode_.init();
-
     const label MAXSTP = 10000;
 
     scalar x = xStart;

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -76,10 +76,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
 {
     // Fix the total force in global coordinate system
     dimensionedVector fAbs =
-        // Force in global coordinate system
+        // External force
         force()
-        // Force in local coordinate system
-      + (R.T() & forceRelative())
         // Spring force in global coordinate system
       - (linSpringCoeffs() & xR)
         // Damping force in global coordinate system
@@ -100,10 +98,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
 {
     // External moment (torque) in local coordinate system
     dimensionedVector mRel =
-        // Moment in global coordinate system
+        // External moment
-        (dimensionedTensor("R", dimless, R) & moment())
+        (dimensionedTensor("R", dimless, R) & moment());
-        // Moment in local coordinate system
-      + momentRelative();
 
     // Note: constraints not implemented at the moment. They shall be
     // implemented in terms of Lagrange multipliers.

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -62,10 +62,8 @@ Foam::dimensionedVector Foam::quaternionSixDOF::A
 {
     // Fix the total force in global coordinate system
     dimensionedVector fAbs =
-        // Force in global coordinate system
+        // External force
         force()
-        // Force in local coordinate system
-      + (rotation.invR() & forceRelative())
         // Spring force in global coordinate system
       - (linSpringCoeffs() & xR)
         // Damping force in global coordinate system
@@ -96,7 +94,6 @@ Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot
             E(omega)
             // To relative
           + (rotation.R() & mAbs)
-          + momentRelative()
         );
 }
 

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -50,11 +50,7 @@ defineRunTimeSelectionTable(sixDOFODE, dictionary);
 
 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
 
-void Foam::sixDOFODE::aitkensRelaxation
+void Foam::sixDOFODE::updateRelaxFactors()
-(
-    const scalar min,
-    const scalar max
-)
 {
     // Calculate translational relax factor
     const scalar saveOldRelFacT = oldRelaxFactorT_;
@@ -71,17 +67,17 @@ void Foam::sixDOFODE::aitkensRelaxation
     }
     else
     {
-        relaxFactorT_ = min;
+        relaxFactorT_ = minRelaxFactor_;
     }
 
     // Bound the relaxation factor for stability
-    if (relaxFactorT_ > max)
+    if (relaxFactorT_ > maxRelaxFactor_)
     {
-        relaxFactorT_ = max;
+        relaxFactorT_ = maxRelaxFactor_;
     }
-    else if (relaxFactorT_ < min)
+    else if (relaxFactorT_ < minRelaxFactor_)
     {
-        relaxFactorT_ = min;
+        relaxFactorT_ = minRelaxFactor_;
     }
 
     // Calculate rotational relax factor
@@ -104,141 +100,35 @@ void Foam::sixDOFODE::aitkensRelaxation
     }
     else
     {
-        relaxFactorR_ = min;
+        relaxFactorR_ = minRelaxFactor_;
     }
 
     // Bound the relaxation factor for stability
-    if(relaxFactorR_ > max)
+    if (relaxFactorR_ > maxRelaxFactor_)
     {
-        relaxFactorR_ = max;
+        relaxFactorR_ = maxRelaxFactor_;
     }
-    else if(relaxFactorR_ < min)
+    else if (relaxFactorR_ < minRelaxFactor_)
     {
-        relaxFactorR_ = min;
+        relaxFactorR_ = minRelaxFactor_;
     }
 }
 
 
-// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+void Foam::sixDOFODE::relaxAcceleration()
-
-void Foam::sixDOFODE::setState(const sixDOFODE& sd)
 {
-    // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_.
+    if (mag(minRelaxFactor_ - maxRelaxFactor_) < SMALL)
-    // 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)
-:
-    ODE(),
-    dict_(io, *this),
-
-    mass_(dict_.lookup("mass")),
-    momentOfInertia_(dict_.lookup("momentOfInertia")),
-
-    Xequilibrium_(dict_.lookup("equilibriumPosition")),
-    linSpringCoeffs_(dict_.lookup("linearSpring")),
-    linDampingCoeffs_(dict_.lookup("linearDamping")),
-
-    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_(dict_.lookup("force")),
-    moment_(dict_.lookup("moment")),
-    forceRelative_(dict_.lookup("forceRelative")),
-    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_()
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sixDOFODE::~sixDOFODE()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::OutputControlDictionary<Foam::sixDOFODE>&
-Foam::sixDOFODE::dict() const
-{
-    return dict_;
-}
-
-
-void Foam::sixDOFODE::relaxAcceleration
-(
-    const scalar minRelFactor,
-    const scalar maxRelFactor
-)
-{
-    if (mag(minRelFactor - maxRelFactor) < SMALL)
     {
        // Fixed relaxation
-       relaxFactorT_ = minRelFactor;
+       relaxFactorT_ = minRelaxFactor_;
-       relaxFactorR_ = minRelFactor;
+       relaxFactorR_ = minRelaxFactor_;
     }
     else
     {
         // Use Aitkens relaxation
 
         // Update Aitkens relaxation factor
-        aitkensRelaxation(minRelFactor, maxRelFactor);
+        updateRelaxFactors();
 
         // Update non relaxed accelerations
         An_[1] = An_[2];
@@ -272,7 +162,7 @@ void Foam::sixDOFODE::relaxAcceleration
 }
 
 
-void Foam::sixDOFODE::init()
+void Foam::sixDOFODE::initState()
 {
     // Get time index
     const label timeIndex = dict().time().timeIndex();
@@ -299,6 +189,113 @@ void Foam::sixDOFODE::init()
 }
 
 
+// * * * * * * * * * * * 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_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::sixDOFODE(const IOobject& io)
+:
+    ODE(),
+    dict_(io, *this),
+
+    mass_(dict_.lookup("mass")),
+    momentOfInertia_(dict_.lookup("momentOfInertia")),
+
+    Xequilibrium_(dict_.lookup("equilibriumPosition")),
+    linSpringCoeffs_(dict_.lookup("linearSpring")),
+    linDampingCoeffs_(dict_.lookup("linearDamping")),
+
+    aitkensRelaxation_
+    (
+        dict_.lookupOrDefault<Switch>("useAitkensRelaxation", false)
+    ),
+    minRelaxFactor_(dict_.lookupOrDefault<scalar>("minRelaxFactor", 0.1)),
+    maxRelaxFactor_(dict_.lookupOrDefault<scalar>("maxRelaxFactor", 0.5)),
+
+    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_(dict_.lookup("force")),
+    moment_(dict_.lookup("moment")),
+
+    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_),
+
+    aitkensRelaxation_(sd.aitkensRelaxation_),
+    minRelaxFactor_(sd.minRelaxFactor_),
+    maxRelaxFactor_(sd.maxRelaxFactor_),
+
+    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_),
+
+    curTimeIndex_(-1),
+    oldStatePtr_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::~sixDOFODE()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::OutputControlDictionary<Foam::sixDOFODE>&
+Foam::sixDOFODE::dict() const
+{
+    return dict_;
+}
+
+
 bool Foam::sixDOFODE::writeData(Ostream& os) const
 {
     os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl;
@@ -312,14 +309,17 @@ bool Foam::sixDOFODE::writeData(Ostream& os) const
     os.writeKeyword("linearDamping") << tab << linDampingCoeffs_
         << token::END_STATEMENT << nl << nl;
 
+    os.writeKeyword("useAitkensRelaxation") << tab << aitkensRelaxation_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("minRelaxFactor") << tab << minRelaxFactor_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("maxRelaxFactor") << tab << maxRelaxFactor_
+        << token::END_STATEMENT << nl;
+
     os.writeKeyword("force") << tab << force_
         << token::END_STATEMENT << nl;
     os.writeKeyword("moment") << tab << moment_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("forceRelative") << tab << forceRelative_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("momentRelative") << tab << momentRelative_
-        << token::END_STATEMENT << endl;
 
     return os.good();
 }

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -48,6 +48,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "OutputControlDictionary.H"
 #include "runTimeSelectionTables.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +93,15 @@ class sixDOFODE
 
         // Aitkens relaxation data members
 
+            //- Switch to control whether to use Aitkens relaxation
+            Switch aitkensRelaxation_;
+
+            //- Minimum acceleration relaxation factor (default 0.1)
+            const scalar minRelaxFactor_;
+
+            //- Maximum acceleration relaxation factor (default 0.5)
+            const scalar maxRelaxFactor_;
+
             //- Translational relaxation factor
             scalar relaxFactorT_;
 
@@ -116,18 +126,12 @@ class sixDOFODE
 
         // External forces and moments
 
-            //- Force driving the motion in inertial coord. sys.
+            //- Force driving the motion in global (inertial) coord. sys.
             dimensionedVector force_;
 
-            //- Moment driving the motion in inertial coord. sys.
+            //- Moment driving the motion in global (inertial) coord. sys.
             dimensionedVector moment_;
 
-            //- Force driving the motion in relative coord. sys.
-            dimensionedVector forceRelative_;
-
-            //- Moment driving the motion in relative coord. sys.
-            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
@@ -157,7 +161,18 @@ class sixDOFODE
         // Aitkens relaxation helper function
 
             //- Update Aitkens relaxation parameters
-            void aitkensRelaxation(const scalar min, const scalar max);
+            void updateRelaxFactors();
+
+            //- Relax the force (acceleration) using Aitkens or fixed relaxation
+            void relaxAcceleration();
+
+
+        // Helper function for controlling multiple ODE solver calls per
+        // time-step
+
+            //- Initialise ODE before setting external forces/moments and
+            //  solving
+            virtual void initState();
 
 
 protected:
@@ -245,41 +260,17 @@ public:
 
         // Access to forces and moments
 
-            //- Return force in inertial coord. sys.
+            //- Return force in global (inertial) coord. sys.
             inline const dimensionedVector& force() const;
 
-            //- Return access to force in inertial coord. sys.
+            //- Return moment in global (inertial) coord. sys.
-            inline dimensionedVector& force();
-
-            //- Return moment in inertial coord. sys.
             inline const dimensionedVector& moment() const;
 
-            //- Return access to moment in inertial coord. sys.
+            //- Set external force and moment
-            inline dimensionedVector& moment();
+            inline void setExternalForceAndMoment
-
-            //- Return force in relative coord. sys.
-            inline const dimensionedVector& forceRelative() const;
-
-            //- Return access to force in relative coord. sys.
-            inline dimensionedVector& forceRelative();
-
-            //- Return moment in relative coord. sys.
-            inline const dimensionedVector& momentRelative() const;
-
-            //- Return access to moment in relative coord. sys.
-            inline dimensionedVector& momentRelative();
-
-            //- Return total force in inertial coord. sys.
-            inline dimensionedVector forceTotal() const;
-
-            //- Return total moment in inertial coord. sys.
-            inline dimensionedVector momentTotal() const;
-
-            //- Relax the force (acceleration) using Aitkens or fixed relaxation
-            void relaxAcceleration
             (
-                const scalar minRelFactor,
+                const vector& externalForce,
-                const scalar maxRelFactor
+                const vector& externalMoment
             );
 
 
@@ -358,10 +349,6 @@ 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

@@ -91,57 +91,30 @@ const Foam::dimensionedVector& Foam::sixDOFODE::force() const
 }
 
 
-Foam::dimensionedVector& Foam::sixDOFODE::force()
-{
-    return force_;
-}
-
-
 const Foam::dimensionedVector& Foam::sixDOFODE::moment() const
 {
     return moment_;
 }
 
 
-Foam::dimensionedVector& Foam::sixDOFODE::moment()
+void Foam::sixDOFODE::setExternalForceAndMoment
+(
+    const vector& externalForce,
+    const vector& externalMoment
+)
 {
-    return moment_;
+    // Initialise the state before setting external forces and moments
-}
+    initState();
 
+    // Set forces and moments
+    force_ = externalForce;
+    moment_ = externalMoment;
 
-const Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() const
+    // Relax acceleration if the Aitkens relaxation is turned on
-{
+    if (aitkensRelaxation_)
-    return forceRelative_;
+    {
-}
+        relaxAcceleration();
-
+    }
-
-Foam::dimensionedVector& Foam::sixDOFODE::forceRelative()
-{
-    return forceRelative_;
-}
-
-
-const Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() const
-{
-    return momentRelative_;
-}
-
-
-Foam::dimensionedVector& Foam::sixDOFODE::momentRelative()
-{
-    return momentRelative_;
-}
-
-
-Foam::dimensionedVector Foam::sixDOFODE::forceTotal() const
-{
-    return force_ + (this->toAbsolute() & forceRelative_);
-}
-
-
-Foam::dimensionedVector Foam::sixDOFODE::momentTotal() const
-{
-    return moment_ + (this->toAbsolute() & momentRelative_);
 }