Fix issues with nonlinear problems when using MPhys wrapper

tacs/TACS.pyx

@@ -1266,6 +1266,13 @@ cdef class KSM:
             self.ptr.decref()
         return
 
+    def setOperators(self, Mat mat, Pc pc):
+        """
+        Set the matrix and preconditioner operators
+        """
+        self.ptr.setOperators(mat.ptr, pc.ptr)
+        return
+
     def solve(self, Vec b, Vec x, int zero_guess=1):
         """
         Try to solve the linear system using GMRES.

tacs/mphys/solver.py

@@ -145,10 +145,12 @@ def solve_nonlinear(self, inputs, outputs):
             self.sp.writeSolution(baseName=f"{self.sp.name}-failed")
             # TODO: In future we could add something here to distinguish between fatal failures and those that could be recovered from
             self.sp.zeroVariables()
-            raise om.AnalysisError("TACS solver did not converge")
 
         self.sp.getVariables(states=outputs[self.states_name])
 
+        if not hasConverged:
+            raise om.AnalysisError("TACS solver did not converge")
+
     def solve_linear(self, d_outputs, d_residuals, mode):
         if mode == "fwd":
             if self.check_partials:

tacs/problems/static.py

@@ -31,6 +31,11 @@ class StaticProblem(TACSProblem):
     defaultOptions = {
         "outputDir": [str, "./", "Output directory for F5 file writer."],
         # Solution Options
+        "resetBeforeSolve": [
+            bool,
+            False,
+            "Reset the states before every solve, this can be useful for avoiding issues when testing derivatives against finite-difference/complex-step",
+        ],
         "linearSolver": [
             str,
             "GMRES",
@@ -354,6 +359,16 @@ def _createVariables(self):
                 freq=opt("monitorFrequency"),
             )
 
+        # Pass new matrix and preconditioner to nonlinear solver linear solvers
+        if self.nonlinearSolver is not None:
+            self.nonlinearSolver.linearSolver.setOperators(self.K, self.PC)
+            try:
+                self.nonlinearSolver.innerSolver.linearSolver.setOperators(
+                    self.K, self.PC
+                )
+            except AttributeError:
+                pass
+
         # Linear solver factor flag
         self._jacobianUpdateRequired = True
         self._preconditionerUpdateRequired = True
@@ -400,6 +415,7 @@ def setOption(self, name, value):
                 "printtiming",
                 "numbersolutions",
                 "outputdir",
+                "printLevel",
             ]:
                 createVariables = False
 
@@ -841,6 +857,8 @@ def solve(self, Fext=None):
         self.callCounter += 1
 
         # Set problem vars to assembler
+        if self.getOption("resetBeforeSolve"):
+            self.zeroVariables()
         self._updateAssemblerVars()
 
         # Check if we need to initialize

tacs/solvers/base.py

@@ -105,7 +105,8 @@ def hasConverged(self) -> bool:
     def fatalFailure(self) -> bool:
         """Whether the solver has failed, set as a property rather than an attribute so that it is read-only
 
-        Note that a fatalFailure is not the same as not converging, this flag is meant to reflect that there has been a fatal failure in the solver which requires a full reset
+        Note that a fatalFailure is not the same as not converging, this flag is meant to reflect that there has been a
+        fatal failure in the solver which requires a full reset
         """
         return self._fatalFailure
 

tacs/solvers/continuation.py

@@ -279,10 +279,7 @@ def solve(
         USE_PREDICTOR = self.getOption("UsePredictor")
         # PREDICTOR_USE_DERIVATIVE = self.getOption("predictorUseDerivative")
 
-        self.initializeSolve()
-
-        if u0 is not None:
-            self.stateVec.copyValues(u0)
+        self.initializeSolve(u0=u0)
 
         # Compute the internal and external forcing vectors at the current point
         self.computeForceVectors()
@@ -421,7 +418,8 @@ def continuationcallBack(solver, u, res, monitorVars):
                     loadStepDirection = 1
                 # If the inner solve failed then we'll reduce the step size and try again, unless we've hit the increment or step size limits
                 elif not isLastIncrement and stepSize > MIN_STEP:
-                    self.setStateFunc(self.incStartState)
+                    self.stateVec.copyValues(self.incStartState)
+                    self.setStateFunc(self.stateVec)
                     currentLambda -= stepSize * loadStepDirection
                     self.setLambdaFunc(currentLambda)
                     stepSize *= STEP_RETRACT_FACTOR

tests/constitutive_tests/test_composite_shell_constitutive.py

@@ -92,7 +92,7 @@ def setUp(self):
         self.ply_thicknesses = np.array([ply_thickness] * nplies, dtype=self.dtype)
         self.ply_angles = np.array([0.0, -45.0, 90.0], dtype=self.dtype) * DEG2RAD
         # Distance between mid-plane and reference plane (nodes) of shell
-        self.tOffset = 0.5 # This places the nodes at the top ply
+        self.tOffset = 0.5  # This places the nodes at the top ply
 
         # Seed random number generator in tacs for consistent test results
         elements.SeedRandomGenerator(0)
@@ -102,7 +102,10 @@ def test_constitutive_density(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveDensity(
                     con,
@@ -122,7 +125,10 @@ def test_constitutive_specific_heat(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveSpecificHeat(
                     con,
@@ -142,7 +148,10 @@ def test_constitutive_heat_flux(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveHeatFlux(
                     con,
@@ -162,7 +171,10 @@ def test_constitutive_stress(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveStress(
                     con,
@@ -182,7 +194,10 @@ def test_constitutive_thermal_strain(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveThermalStrain(
                     con,
@@ -202,7 +217,10 @@ def test_constitutive_failure(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveFailure(
                     con,
@@ -221,7 +239,10 @@ def test_constitutive_failure_strain_sens(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
+                    layup,
+                    self.ply_thicknesses,
+                    self.ply_angles,
+                    tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveFailureStrainSens(
                     con,

tests/constitutive_tests/test_iso_shell_constitutive.py

@@ -46,7 +46,9 @@ def setUp(self):
         )
 
         # Create stiffness (need class)
-        self.con = constitutive.IsoShellConstitutive(self.props, t=1.0, tNum=0, tOffset=-0.5)
+        self.con = constitutive.IsoShellConstitutive(
+            self.props, t=1.0, tNum=0, tOffset=-0.5
+        )
 
         # Seed random number generator in tacs for consistent test results
         elements.SeedRandomGenerator(0)