Merge branch 'main' into ralberd/material_qoi

examples/arch_bc/ArchBc.py

@@ -156,10 +156,10 @@ def run(self):
 
             disp -= maxDisp/steps
             p = Objective.param_index_update(p, 0, disp)
-            Uu = EqSolver.nonlinear_equation_solve(objective,
-                                                   Uu,
-                                                   p,
-                                                   self.trSettings)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective,
+                                                                  Uu,
+                                                                  p,
+                                                                  self.trSettings)
 
             self.write_output(Uu, p, i)
 

examples/arch_bc/ArchTraction.py

@@ -172,7 +172,7 @@ def run(self):
 
             force += maxForce/steps
             p = Objective.param_index_update(p, 0, force)
-            Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
 
             self.write_output(Uu, p, i)
 
@@ -184,7 +184,7 @@ def run(self):
 
                 force -= maxForce/steps
                 p = Objective.param_index_update(p, 0, force)
-                Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
+                Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
 
                 self.write_output(Uu, p, i)
 

examples/configurational_forces/CrackConfigurationalForce.py

@@ -132,8 +132,7 @@ def run():
         appliedK += 1.0
         p = Objective.param_index_update(p, 0, appliedK)
 
-        Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)
-
+        Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)
 
         Wu = np.zeros_like(Uu)
         _, configForcesOnUnknowns = sensit(Wu, Uu, p, mesh)

examples/eigenvalue_sensitivity/nonlinear_elasticity_eigenvalue_sensitivity.ipynb

@@ -414,7 +414,7 @@
     "\n",
     "force = 2.0\n",
     "p = Objective.param_index_update(p, 0, force) # put current total force in parameter set\n",
-    "Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings, useWarmStart=False)\n"
+    "Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings, useWarmStart=False)\n"
    ]
   },
   {
@@ -718,7 +718,7 @@
     "df = 1e-4\n",
     "forcePlus = force + df\n",
     "p = Objective.param_index_update(p, 0, forcePlus)\n",
-    "UuPlus = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings, useWarmStart=False)"
+    "UuPlus, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings, useWarmStart=False)"
    ]
   },
   {

examples/euler_buckle/EulerBuckle.py

@@ -154,7 +154,7 @@ def run(self):
             print('LOAD STEP ', i)
             force += maxForce/steps
             p = Objective.param_index_update(p, 0, force)
-            Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
 
             self.write_output(Uu, p, i)
 
@@ -165,7 +165,7 @@ def run(self):
                 print('LOAD STEP ', i)
                 force -= maxForce/steps
                 p = Objective.param_index_update(p, 0, force)
-                Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
+                Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings, solver_algorithm=solver)
 
                 self.write_output(Uu, p, i)
 

examples/hemisphere_cap/hemisphere_plastic_disp_control.py

@@ -156,8 +156,7 @@ def run(self):
 
             disp += maxDisp/steps
             p = Objective.param_index_update(p, 0, disp)
-            Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
-
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
 
             U = self.dofManager.create_field(Uu, self.get_ubcs(p))
             ivs = self.bvpFuncs.\
@@ -175,7 +174,7 @@ def run(self):
 
                 disp -= maxDisp/steps
                 p = Objective.param_index_update(p, 0, disp)
-                Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
+                Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
 
                 U = self.dofManager.create_field(Uu, self.get_ubcs(p))
                 ivs = self.bvpFuncs.\

examples/hemisphere_cap/hemisphere_plastic_load_control.py

@@ -197,7 +197,7 @@ def run(self):
 
             force += maxForce/steps
             p = Objective.param_index_update(p, 0, force)
-            Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
 
             U = self.dofManager.create_field(Uu, self.get_ubcs(p))
             ivs = self.bvpFuncs.compute_updated_internal_variables(U, p[1])
@@ -214,7 +214,7 @@ def run(self):
 
                 force -= maxForce/steps
                 p = Objective.param_index_update(p, 0, force)
-                Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
+                Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, self.trSettings)
 
                 U = self.dofManager.create_field(Uu, self.get_ubcs(p))
                 ivs = self.bvpFuncs.compute_updated_internal_variables(U, p[1])

examples/inverse/BuckleInverse.py

@@ -221,7 +221,7 @@ def run(self):
             endDispOrTraction = maxDispOrTraction*(i+1)/N
 
             p = Objective.param_index_update(p, 0, np.array([endDispOrTraction, 0.0]))
-            self.Uu = EquationSolver.nonlinear_equation_solve(self.objective, self.Uu, p, settings)
+            self.Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(self.objective, self.Uu, p, settings)
 
             # post process
             print('step ', i)

examples/multi_material_tutorial/optimism_tutorial.ipynb

@@ -801,7 +801,7 @@
     "p = Objective.param_index_update(p, 0, appliedDisp)\n",
     "\n",
     "# Find unknown displacements by minimizing the objective\n",
-    "Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)\n",
+    "Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)\n",
     "\n"
    ]
   },

examples/multi_material_tutorial/uniaxial_two_material.py

@@ -144,7 +144,7 @@ def assemble_sparse_preconditioner_matrix(Uu, p):
 p = Objective.param_index_update(p, 0, appliedDisp)
 
 # Find unknown displacements by minimizing the objective
-Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)
+Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)
 
 # update the state variables in the new equilibrium configuration
 U = create_field(Uu, p)

examples/surfing_fracture/SurfingFracture.py

@@ -462,11 +462,11 @@ def run_alternating_min(self):
                 print("------------------------------")
 
                 print("Minimizing disp: objective = ", dispObjective.get_value(UuDisp))
-                UuDisp = EqSolver.nonlinear_equation_solve(dispObjective,
-                                                           UuDisp,
-                                                           p,
-                                                           dispSettings,
-                                                           useWarmStart=False)
+                UuDisp, solverSuccess = EqSolver.nonlinear_equation_solve(dispObjective,
+                                                                          UuDisp,
+                                                                          p,
+                                                                          dispSettings,
+                                                                          useWarmStart=False)
 
                 Uu = Uu.at[self.dispIds].set(UuDisp)
                 print("Minimized disp: objective = ", dispObjective.get_value(UuDisp))

examples/uniaxial/Uniaxial.py

@@ -167,7 +167,7 @@ def run(self):
             xDisp = i/N*maxDisp
             p = Objective.param_index_update(p, 0, xDisp)
 
-            Uu = EqSolver.nonlinear_equation_solve(objective, Uu, p, settings)
+            Uu, solverSuccess = EqSolver.nonlinear_equation_solve(objective, Uu, p, settings)
 
             state = self.mechanicsFunctions.\
                 compute_updated_internal_variables(self.create_field(Uu, p), p[1])

examples/uniaxial_dynamic/UniaxialDynamic.py

@@ -208,11 +208,11 @@ def run(self):
             #print('Predicted Vx = \n', self.create_field(Vu,p)[:,0])
             p = Objective.param_index_update(p, 5, UuPredicted)
 
-            Uu = EquationSolver.nonlinear_equation_solve(objective,
-                                                         Uu,
-                                                         p,
-                                                         solverSettings,
-                                                         useWarmStart=False)
+            Uu, solverSuccess = EquationSolver.nonlinear_equation_solve(objective,
+                                                                        Uu,
+                                                                        p,
+                                                                        solverSettings,
+                                                                        useWarmStart=False)
 
             UuCorrection = Uu - UuPredicted
             #print('Corrected Ux = \n', self.create_field(UuCorrection,p)[:,0])

optimism/AlSolver.py

@@ -42,20 +42,21 @@ def get_settings(penalty_scaling=4.,
 
 
 def solve_sub_step(alObjective, x, ncpErrorOld, alSettings, subSettings, sub_problem_solver, sub_problem_callback=None):
-    x = sub_problem_solver(alObjective, x, subSettings, sub_problem_callback)
+    x, solverSuccess = sub_problem_solver(alObjective, x, subSettings, sub_problem_callback)
 
     c = alObjective.constraint(x)
     kappa = alObjective.kappa
     alObjective.lam = np.maximum(alObjective.lam-kappa*c, 0.0)
 
     ncpError = np.abs( alObjective.ncp(x) )
-    poorProgress = ncpError > alSettings.target_constraint_decrease_factor * ncpErrorOld
+    # check if each constraint is making progress, or if they are already small relative to the specificed constraint tolerances
+    poorProgress = ncpError > np.maximum(alSettings.target_constraint_decrease_factor * ncpErrorOld, 10 * alSettings.tol / np.sqrt(len(ncpError)))
 
-    if (np.any(poorProgress)):
+    if (np.any(poorProgress) and solverSuccess):
         print('Poor progress on ncp detected, increasing some penalty parameters')
         alObjective.kappa = kappa.at[poorProgress].set(alSettings.penalty_scaling*kappa[poorProgress])
 
-    return x, ncpError
+    return x, ncpError, solverSuccess
 
 
 def linear_update(alObjective, x, rhs_func, alSettings):
@@ -168,13 +169,13 @@ def augmented_lagrange_solve(alObjective, x, p, alSettings, subSettings,
             alObjective.update_precond(x)
 
         if not alSettings.use_newton_only:
-            x, ncpError = solve_sub_step(alObjective,
-                                         x,
-                                         ncpError,
-                                         alSettings,
-                                         settings,
-                                         sub_problem_solver,
-                                         sub_problem_callback)
+            x, ncpError, solverSuccess = solve_sub_step(alObjective,
+                                                        x,
+                                                        ncpError,
+                                                        alSettings,
+                                                        settings,
+                                                        sub_problem_solver,
+                                                        sub_problem_callback)
 
             forceErrorNorm = norm(alObjective.gradient(x))
             print('force error = ', forceErrorNorm)
@@ -187,6 +188,5 @@ def augmented_lagrange_solve(alObjective, x, p, alSettings, subSettings,
                 return x
 
     raise NameError('Loadstep failed to converge in', maxAlIters, 'iterations.')
-
     return x
 

optimism/EquationSolver.py

@@ -278,7 +278,7 @@ def trust_region_least_squares_solve(objective, x, settings):
             print_banner(oy, modelObjective, cgIters, trSize, stepType, True)
             if settings.check_stability:
                 objective.check_stability(y)
-            return y
+            return y, True
 
         modelImprove = o - modelObjective
         realImprove = o - oy
@@ -300,11 +300,10 @@ def trust_region_least_squares_solve(objective, x, settings):
             o = oy
             J.update( objective.hessian(x) )
 
-
     print("Reached the maximum number of trust region iterations or trust region is too small.")
     if settings.check_stability:
         objective.check_stability(x)
-    return x
+    return x, False
 
 
 def is_converged(objective, x, realO, modelO, realRes, modelRes, cgIters, trSize, settings):
@@ -369,7 +368,7 @@ def trust_region_minimize(objective, x, settings, callback=None):
     # this could potentially return an unstable solution
     if is_converged(objective, x, 0.0, 0.0, g, g, 0, trSize, settings):
         if callback: callback(x, objective)
-        return x
+        return x, True
 
     cumulativeCgIters=0
 
@@ -427,7 +426,7 @@ def trust_region_minimize(objective, x, settings, callback=None):
             if is_converged(objective, y, realObjective, modelObjective,
                             gy, g + Jd, cgIters, trSizeUsed, settings):
                 if callback: callback(y, objective)
-                return y
+                return y, True
 
             modelImprove = -modelObjective
             realImprove = -realObjective
@@ -490,14 +489,14 @@ def trust_region_minimize(objective, x, settings, callback=None):
                 else:
                     print("The trust region is still too small.  Accepting, but be careful.")
                     if callback: callback(x, objective)
-                    return x
+                    return x, False
 
     print("Reached the maximum number of trust region iterations.")
     if settings.check_stability:
         objective.check_stability(x)
 
         if callback: callback(x, objective)
-    return x
+    return x, False
 
 
 def newton(objective, x, settings, callback=None):
@@ -546,12 +545,12 @@ def apply_a(v):
         o = objective.value(x)
         if is_converged(objective, x, o, o,
                         g, g, 1, np.inf, settings):
-            return x
+            return x, True
 
     print("Reached the maximum number of newton iterations.")
     if settings.check_stability:
         objective.check_stability(x)
-    return x
+    return x, False
 
 
 # This solve uses a different interface.  There is no globalization yet
@@ -579,7 +578,7 @@ def newton_solve(objective_func, solution, maxSteps=1):
             with Timer(name="LinearSolve"):
                 solution = solution - np.reshape(np.linalg.solve(hess,gradients), gradShape)
 
-    return solution
+    return solution, True
 
 
 def nonlinear_equation_solve(objective, x0, p, settings,
@@ -603,7 +602,7 @@ def nonlinear_equation_solve(objective, x0, p, settings,
     if updatePrecond:
         objective.update_precond(xBar0)
 
-    xBar = solver_algorithm(objective, xBar0, settings, callback=callback)
+    xBar, solverSuccess = solver_algorithm(objective, xBar0, settings, callback=callback)
 
-    return objective.invScaling * xBar
+    return objective.invScaling * xBar, solverSuccess
 

optimism/TrustRegionSPG.py

@@ -342,7 +342,7 @@ def bound_constrained_trust_region_minimize(objective, x, bounds, settings, call
     if is_converged(objective, x, 0.0, 0.0, prevOptimality, prevOptimality, 0,
                     trSize, settings):
         if callback: callback(x, objective)
-        return x
+        return x, True
 
     # Set guess for first cauchy point step length
     gHg = g@objective.hessian_vec(x, g)
@@ -391,7 +391,7 @@ def bound_constrained_trust_region_minimize(objective, x, bounds, settings, call
                         realOptimality, modelOptimality, spgIters, trSizeUsed,
                         settings):
             if callback: callback(y, objective)
-            return y
+            return y, True
 
         modelImprove = -modelObjective
         realImprove = -realObjective
@@ -445,14 +445,14 @@ def bound_constrained_trust_region_minimize(objective, x, bounds, settings, call
             else:
                 print("The trust region is still too small.  Accepting, but be careful.")
                 if callback: callback(x, objective)
-                return x
+                return x, False
 
     print("Reached the maximum number of trust region iterations.")
     if settings.check_stability:
         objective.check_stability(x)
 
         if callback: callback(x, objective)
-    return x
+    return x, False
 
 
 def solve(objective, x0, p, lowerBounds, upperBounds, settings, callback=None,
@@ -477,8 +477,8 @@ def solve(objective, x0, p, lowerBounds, upperBounds, settings, callback=None,
 
     bounds = np.column_stack((lBar, uBar))
 
-    xBar = bound_constrained_trust_region_minimize(objective, xBar0, bounds, settings,
+    xBar, solverSuccess = bound_constrained_trust_region_minimize(objective, xBar0, bounds, settings,
                                                    callback=callback)
 
-    return objective.invScaling * xBar
+    return objective.invScaling * xBar, solverSuccess