ODESolver: added stepState sub-class to carry additional information needed during sub-cycling particularly for seulex

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -82,49 +82,72 @@ Foam::scalar Foam::ODESolver::normalizeError
 }
 
 
+void Foam::ODESolver::solve
+(
+    scalar& x,
+    scalarField& y,
+    stepState& step
+) const
+{
+    solve(x, y, step.dxTry);
+}
+
+
 void Foam::ODESolver::solve
 (
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
-    scalar& dxEst
+    scalar& dxTry
 ) const
 {
+    stepState step(dxTry);
     scalar x = xStart;
-    bool truncated = false;
 
     for (label nStep=0; nStep<maxSteps_; nStep++)
     {
-        // Store previous iteration dxEst
-        scalar dxEst0 = dxEst;
+        // Store previous iteration dxTry
+        scalar dxTry0 = step.dxTry;
 
-        // Check if this is a truncated step and set dxEst to integrate to xEnd
-        if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0)
+        step.reject = false;
+
+        // Check if this is a truncated step and set dxTry to integrate to xEnd
+        if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
         {
-            truncated = true;
-            dxEst = xEnd - x;
+            step.last = true;
+            step.dxTry = xEnd - x;
         }
 
-        // Integrate as far as possible up to dxEst
-        solve(x, y, dxEst);
+        // Integrate as far as possible up to step.dxTry
+        solve(x, y, step);
 
         // Check if reached xEnd
         if ((x - xEnd)*(xEnd - xStart) >= 0)
         {
-            if (nStep > 0 && truncated)
+            if (nStep > 0 && step.last)
             {
-                dxEst = dxEst0;
+                step.dxTry = dxTry0;
             }
 
+            dxTry = step.dxTry;
+
             return;
         }
+
+        step.first = false;
+
+        // If the step.dxTry was reject set step.prevReject
+        if (step.reject)
+        {
+            step.prevReject = true;
+        }
     }
 
     FatalErrorIn
     (
         "ODESolver::solve"
         "(const scalar xStart, const scalar xEnd,"
-        "scalarField& y, scalar& dxEst) const"
+        "scalarField& y, scalar& dxTry) const"
     )   << "Integration steps greater than maximum " << maxSteps_
         << exit(FatalError);
 }

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -93,6 +93,30 @@ public:
     //- Runtime type information
     TypeName("ODESolver");
 
+    class stepState
+    {
+        public:
+
+        const bool forward;
+        scalar dxTry;
+        scalar dxDid;
+        bool first;
+        bool last;
+        bool reject;
+        bool prevReject;
+
+        stepState(const scalar dx)
+        :
+            forward(dx > 0 ? true : false),
+            dxTry(dx),
+            dxDid(0),
+            first(true),
+            last(false),
+            reject(false),
+            prevReject(false)
+        {}
+    };
+
 
     // Declare run-time constructor selection table
 
@@ -156,7 +180,23 @@ public:
             scalar& x,
             scalarField& y,
             scalar& dxTry
-        ) const = 0;
+        ) const //= 0;
+        {
+            stepState step(dxTry);
+            solve(x, y, step);
+            dxTry = step.dxTry;
+        }
+
+        //- Solve the ODE system as far as possible upto dxTry
+        //  adjusting the step as necessary to provide a solution within
+        //  the specified tolerance.
+        //  Update the state and return an estimate for the next step in dxTry
+        virtual void solve
+        (
+            scalar& x,
+            scalarField& y,
+            stepState& step
+        ) const;
 
         //- Solve the ODE system from xStart to xEnd, update the state
         //  and return an estimate for the next step in dxTry