Implement matrix assembly and factorization persistent through entire simulation

src/chrono/timestepper/ChTimestepperHHT.cpp

@@ -34,7 +34,8 @@ ChTimestepperHHT::ChTimestepperHHT(ChIntegrableIIorder* intgr)
       h_min(1e-10),
       h(1e6),
       num_successful_steps(0),
-      modified_Newton(true) {
+      modified_Newton(EVERY_STEP),
+      call_setup(true) {
     SetAlpha(-0.2);  // default: some dissipation
 }
 
@@ -50,6 +51,7 @@ void ChTimestepperHHT::SetAlpha(double val) {
 
 // Performs a step of HHT (generalized alpha) implicit for II order systems
 void ChTimestepperHHT::Advance(const double dt) {
+
     // Downcast
     ChIntegrableIIorder* integrable2 = (ChIntegrableIIorder*)this->integrable;
 
@@ -97,14 +99,9 @@ void ChTimestepperHHT::Advance(const double dt) {
         h = std::min(h, dt);
     }
 
-    // Monitor flags controlling whther or not the Newton matrix must be updated.
-    // If using modified Newton, a matrix update occurs:
-    //   - at the beginning of a step
-    //   - on a stepsize decrease
-    //   - if the Newton iteration does not converge with an out-of-date matrix
-    // Otherwise, the matrix is updated at each iteration.
+    // TODO: Flag overwritten after every step inside Increment(). initializing it here does not seem useful
     matrix_is_current = false;
-    call_setup = true;
+    bool converged = false; // Taken outside the while loop so each sub-step knows if the previous step converged. TODO: remove this comment after review
 
     // Loop until reaching final time
     while (true) {
@@ -113,11 +110,24 @@ void ChTimestepperHHT::Advance(const double dt) {
         // Newton for state at T+h
         Da_nrm_hist.fill(0.0);
         Dl_nrm_hist.fill(0.0);
-        bool converged = false;
         unsigned int it;
 
         for (it = 0; it < maxiters; it++) {
-            if (verbose && modified_Newton && call_setup)
+            // Monitor flags controlling whether or not the Newton matrix must be updated.
+            // If using ModifiedNewton::UNMODIFIED, a matrix update occurs:
+            //   - at every iteration
+            // If using ModifiedNewton::EVERY_STEP, a matrix update occurs:
+            //   - at the beginning of a step
+            //   - on a stepsize decrease
+            //   - if the Newton iteration does not converge with an out-of-date matrix
+            // If using ModifiedNewton::CONSTANT, a matrix update occurs:
+            //   - only at the beginning of the very first step
+            // If using ModifiedNewton::AUTOMATIC, a matrix update occurs:
+            //   - when appropriate. TODO: implement
+            // Otherwise, the matrix is updated at each iteration.
+            call_setup = SetupRequiredForIteration(it, converged);
+
+            if (verbose && (modified_Newton != ModifiedNewton::UNMODIFIED) && call_setup)
                 std::cout << " HHT call Setup." << std::endl;
 
             // Solve linear system and increment state
@@ -130,9 +140,6 @@ void ChTimestepperHHT::Advance(const double dt) {
                 numsetups++;
             }
 
-            // If using modified Newton, do not call Setup again
-            call_setup = !modified_Newton;
-
             // Check convergence
             converged = CheckConvergence(it);
             if (converged)
@@ -178,7 +185,7 @@ void ChTimestepperHHT::Advance(const double dt) {
                     std::cout << " HHT re-attempt step with updated matrix." << std::endl;
                 }
 
-                call_setup = true;
+                call_setup = true; // TODO: update this if this code block ever get uncommented
             */
 
         } else if (!step_control) {
@@ -217,9 +224,6 @@ void ChTimestepperHHT::Advance(const double dt) {
                     std::cerr << " HHT at minimum stepsize. Exiting..." << std::endl;
                 throw std::runtime_error("HHT: Reached minimum allowable step size.");
             }
-
-            // force a matrix re-evaluation (due to change in stepsize)
-            call_setup = true;
         }
 
         if (T >= tfinal) {
@@ -356,6 +360,30 @@ void ChTimestepperHHT::CalcErrorWeights(const ChVectorDynamic<>& x, double rtol,
     ewt = (rtol * x.cwiseAbs() + atol).cwiseInverse();
 }
 
+bool ChTimestepperHHT::SetupRequiredForIteration(int iteration, bool previous_substep_converged) {
+    bool setup;
+    switch (modified_Newton) {
+        case UNMODIFIED:
+            setup = true;
+            break;
+        case EVERY_STEP:
+            // Force a matrix re-evaluation (due to change in stepsize) if previous sub-step did not converge
+            // Do not setup again at start of sub-step if previous sub-step converged
+            setup = (iteration == 0 && !previous_substep_converged);
+            break;
+        case CONSTANT:
+            // Only call setup on very first iteration of very first run
+            // call_setup(true) in constructor, or reset to true if modified Newton changed to another style
+            setup = (call_setup && iteration == 0);
+            break;
+        case AUTOMATIC:
+            // TODO: implement. function potentially needs more arguments to determine automatic update schedule
+            setup = true; // TEMP
+            break;
+    }
+    return setup;
+}
+
 void ChTimestepperHHT::ArchiveOut(ChArchiveOut& archive) {
     // version number
     archive.VersionWrite<ChTimestepperHHT>();

src/chrono/timestepper/ChTimestepperHHT.h

@@ -70,12 +70,26 @@ class ChApi ChTimestepperHHT : public ChTimestepperIIorder, public ChImplicitIte
     /// Default: 0.5.
     void SetStepDecreaseFactor(double factor) { step_decrease_factor = factor; }
 
+    /// Available styles of modified Newton methods
+    enum ModifiedNewton : int {
+        UNMODIFIED, // Classical Newton. Jacobian updated at every iteration
+        EVERY_STEP, // Jacobian updated at every step
+        CONSTANT,   // Jacobian never updated
+        AUTOMATIC   // Jaobian updated when appropriate. TODO: implement
+    };
+
     /// Enable/disable modified Newton.
     /// If enabled, the Newton matrix is evaluated, assembled, and factorized only once
     /// per step or if the Newton iteration does not converge with an out-of-date matrix.
     /// If disabled, the Newton matrix is evaluated at every iteration of the nonlinear solver.
     /// Default: true.
-    void SetModifiedNewton(bool enable) { modified_Newton = enable; }
+    void SetModifiedNewton(ModifiedNewton style) { modified_Newton = style; }
+    void SetModifiedNewton(bool enable) { // Overload for backward compatibility with boolean implementation
+        if (enable)
+            modified_Newton = ModifiedNewton::EVERY_STEP;
+        else
+            modified_Newton = ModifiedNewton::UNMODIFIED;
+    }
 
     /// Perform an integration timestep, by advancing the state by the specified time step.
     virtual void Advance(const double dt) override;
@@ -96,6 +110,7 @@ class ChApi ChTimestepperHHT : public ChTimestepperIIorder, public ChImplicitIte
     void Increment(ChIntegrableIIorder* integrable2);
     bool CheckConvergence(int it);
     void CalcErrorWeights(const ChVectorDynamic<>& x, double rtol, double atol, ChVectorDynamic<>& ewt);
+    bool SetupRequiredForIteration(int iteration, bool previous_substep_converged);
 
   private:
     double alpha;  ///< HHT method parameter:  -1/3 <= alpha <= 0
@@ -110,7 +125,7 @@ class ChApi ChTimestepperHHT : public ChTimestepperIIorder, public ChImplicitIte
     ChVectorDynamic<> Lnew;  ///< current estimate of Lagrange multipliers
     ChVectorDynamic<> R;     ///< residual of nonlinear system (dynamics portion)
     ChVectorDynamic<> Rold;  ///< residual terms depending on previous state
-    ChVectorDynamic<> Qc;    ///< residual of nonlinear system (constranints portion)
+    ChVectorDynamic<> Qc;    ///< residual of nonlinear system (constraints portion)
 
     std::array<double, 3> Da_nrm_hist;  ///< last 3 update norms
     std::array<double, 3> Dl_nrm_hist;  ///< last 3 update norms
@@ -125,7 +140,7 @@ class ChApi ChTimestepperHHT : public ChTimestepperIIorder, public ChImplicitIte
     double h;                           ///< internal stepsize
     unsigned int num_successful_steps;  ///< number of successful steps
 
-    bool modified_Newton;    ///< use modified Newton?
+    ModifiedNewton modified_Newton;    ///< style of modified Newton?
     bool matrix_is_current;  ///< is the Newton matrix up-to-date?
     bool call_setup;         ///< should the solver's Setup function be called?