From 25472673d968f55d431992a4a015b507208238ff Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 13 Jun 2018 17:43:36 +0100
Subject: [PATCH] ENH: prescribedRotation: fix testcase. Fixes #872.

---
 .../prescribedRotation/prescribedRotation.C   | 16 +++++--
 .../rigidBodyMeshMotion/rigidBodyMeshMotion.C | 48 +++++++++++--------
 .../rigidBodyMeshMotionSolver.C               |  3 +-
 .../overInterDyMFoam/boatAndPropeller/Allrun  |  2 +-
 .../boatAndPropeller/constant/dynamicMeshDict |  4 ++
 5 files changed, 46 insertions(+), 27 deletions(-)

diff --git a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
index 483e61c29c4..a2f0f19cfa7 100644
--- a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
+++ b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
@@ -122,19 +122,27 @@ void Foam::RBD::restraints::prescribedRotation::restrain
     vector omega = model_.v(model_.master(bodyID_)).w();
     scalar Inertia = mag(model_.I(model_.master(bodyID_)).Ic());
 
-
     // from the definition of the angular momentum:
     //  moment = Inertia*ddt(omega)
-    // with time step 1:
-    //  moment = Inertia*(omega_wanted - omega)
+    const scalar relax = 0.5;
+
     vector moment
     (
-        (Inertia * (omegaSet_.value(model_.time().value()) - omega) & a) * a
+        (
+            relax
+          * Inertia
+          * (omegaSet_.value(model_.time().value()) - omega)
+          / model_.time().deltaTValue()
+          & a
+        )
+      * a
     );
 
     if (model_.debug)
     {
         Info<< " angle  " << theta*sign(a & axis_) << endl
+            << " omega  " << omega << endl
+            << " wanted " << omegaSet_.value(model_.time().value()) << endl
             << " moment " << moment << endl
             << " oldDir " << oldDir << endl
             << " newDir " << newDir << endl
diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
index 93f2fbffded..d69537ff5e7 100644
--- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
+++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
@@ -263,32 +263,37 @@ void Foam::rigidBodyMeshMotion::solve()
     }
     else
     {
-        Field<spatialVector> fx(model_.nBodies(), Zero);
+        label nIter(coeffDict().lookupOrDefault("nIter", 1));
 
-        forAll(bodyMeshes_, bi)
+        for (label i=0; i<nIter; i++)
         {
-            const label bodyID = bodyMeshes_[bi].bodyID_;
+            Field<spatialVector> fx(model_.nBodies(), Zero);
 
-            dictionary forcesDict;
-            forcesDict.add("type", functionObjects::forces::typeName);
-            forcesDict.add("patches", bodyMeshes_[bi].patches_);
-            forcesDict.add("rhoInf", rhoInf_);
-            forcesDict.add("rho", rhoName_);
-            forcesDict.add("CofR", vector::zero);
+            forAll(bodyMeshes_, bi)
+            {
+                const label bodyID = bodyMeshes_[bi].bodyID_;
 
-            functionObjects::forces f("forces", db(), forcesDict);
-            f.calcForcesMoment();
+                dictionary forcesDict;
+                forcesDict.add("type", functionObjects::forces::typeName);
+                forcesDict.add("patches", bodyMeshes_[bi].patches_);
+                forcesDict.add("rhoInf", rhoInf_);
+                forcesDict.add("rho", rhoName_);
+                forcesDict.add("CofR", vector::zero);
 
-            fx[bodyID] = ramp*spatialVector(f.momentEff(), f.forceEff());
-        }
+                functionObjects::forces f("forces", db(), forcesDict);
+                f.calcForcesMoment();
 
-        model_.solve
-        (
-            t.value(),
-            t.deltaTValue(),
-            scalarField(model_.nDoF(), Zero),
-            fx
-        );
+                fx[bodyID] = ramp*spatialVector(f.momentEff(), f.forceEff());
+            }
+
+            model_.solve
+            (
+                t.value(),
+                t.deltaTValue(),
+                scalarField(model_.nDoF(), Zero),
+                fx
+            );
+        }
     }
 
     if (Pstream::master() && model_.report())
@@ -354,7 +359,8 @@ bool Foam::rigidBodyMeshMotion::writeObject
     );
 
     model_.state().write(dict);
-    return dict.regIOobject::write();
+    // Force ascii writing
+    return dict.regIOobject::writeObject(IOstream::ASCII, ver, cmp, valid);
 }
 
 
diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
index 56ac81c80ec..abf03128124 100644
--- a/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
+++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
@@ -296,7 +296,8 @@ bool Foam::rigidBodyMeshMotionSolver::writeObject
     );
 
     model_.state().write(dict);
-    return dict.regIOobject::write();
+    // Force ascii writing
+    return dict.regIOobject::writeObject(IOstream::ASCII, ver, cmp, valid);
 }
 
 
diff --git a/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/Allrun b/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/Allrun
index d5979700f56..787d72a9dce 100755
--- a/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/Allrun
+++ b/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/Allrun
@@ -5,6 +5,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 ./Allrun.pre
 
 runApplication decomposePar
-runParallel overInterDyMFoam
+runParallel $(getApplication)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/constant/dynamicMeshDict b/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/constant/dynamicMeshDict
index de53a7d4241..4142814ded2 100644
--- a/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/constant/dynamicMeshDict
+++ b/tutorials/multiphase/overInterDyMFoam/boatAndPropeller/constant/dynamicMeshDict
@@ -34,6 +34,10 @@ motionSolver        rigidBodyMotion;
 
     accelerationRelaxation 1.0;
 
+    //- prescribedRotation requires some sub-iterations to converge faster
+    //  to desired value.
+    nIter   3;
+
     bodies
     {
         hull
-- 
GitLab