From 7e0eee0f941080150329b4767247ca019b2a4d52 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 1 Mar 2017 15:38:38 +0100 Subject: [PATCH] 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 | 9 ++ src/ODE/ODESolvers/ODESolver/ODESolver.C | 5 +- src/ODE/ODESolvers/ODESolver/ODESolver.H | 8 +- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 77 +++++++-------- .../sixDOF/geometricSixDOF/geometricSixDOF.H | 16 ++-- .../quaternionSixDOF/quaternionSixDOF.C | 79 +++++++-------- .../quaternionSixDOF/quaternionSixDOF.H | 16 ++-- src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C | 2 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 95 +++++++++++++++---- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 45 ++++++--- src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H | 2 +- 11 files changed, 221 insertions(+), 133 deletions(-) diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H index 6b86f9107..1b9aba269 100644 --- a/src/ODE/ODE/ODE.H +++ b/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 + } }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 39218d7ff..ccfc18bda 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/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; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index aed2b49dc..ab11a06fa 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/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; + ); }; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index f70ad4908..6cdd54b23 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/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(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 - ( - IOobject - ( - name, - gsd.dict().instance(), - gsd.dict().local(), - gsd.dict().db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ) - ), + sixDOFODE(name, gsd), 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(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_; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H index 64a26e9e5..ee5e660e4 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H +++ b/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 diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 45143c562..3e7ac9d57 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/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(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 - ( - IOobject - ( - name, - qsd.dict().instance(), - qsd.dict().local(), - qsd.dict().db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ) - ), + sixDOFODE(name, qsd), 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(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_; diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H index 9eed2941b..59a72899b 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H +++ b/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 diff --git a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C index 30e0e8eca..3fe6d9c7c 100644 --- a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C +++ b/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. diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index e0c70b045..485108c7e 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/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 - // relaxFactorOld_. In case of multiple updates, overwriting Aitkens - // relaxation parameters would invalidate the underrelaxation. - // IG, 5/May/2016 - mass_ = sd.mass_; - momentOfInertia_ = sd.momentOfInertia_; + // Get time index + const label timeIndex = dict().time().timeIndex(); - Xequilibrium_ = sd.Xequilibrium_; - linSpringCoeffs_ = sd.linSpringCoeffs_; - linDampingCoeffs_ = sd.linDampingCoeffs_; + if (curTimeIndex_ < timeIndex) + { + // First call in this time index, store data + oldStatePtr_ = this->clone(dict().name() + "_0"); - force_ = sd.force_; - moment_ = sd.moment_; - forceRelative_ = sd.forceRelative_; - momentRelative_ = sd.momentRelative_; + Info<< "First 6DOF solution within a time step, storing old data..." + << endl; + } + 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; } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index 70ab14f38..d88adc759 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/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 - rotations in quaternions or rotation tensors. + equations with necessary interface for two-way coupling with CFD solvers. 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 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 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 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 diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index b45deb964..01b9dfad5 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H +++ b/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.