Fix a case when mfront would raise an exception and leave PetSc vectors in an unwanted state

MathLib/LinAlg/LinAlg.h

@@ -17,6 +17,8 @@
 
 #ifdef USE_PETSC
 #include <petscsystypes.h>
+
+#include "PETSc/PETScVector.h"
 #endif
 
 namespace MathLib
@@ -204,6 +206,14 @@ void linearSysNormalize(PETScMatrix const& A, PETScMatrix& new_A,
 void finalizeAssembly(PETScMatrix& A);
 void finalizeAssembly(PETScVector& x);
 
+// Reduce operations for interprocess communications while using Petsc
+static inline int reduceMin(int val)
+{
+    int result;
+    MPI_Allreduce(&val, &result, 1, MPI_INTEGER, MPI_MIN, PETSC_COMM_WORLD);
+    return result;
+}
+
 }  // namespace LinAlg
 }  // namespace MathLib
 
@@ -276,6 +286,12 @@ void linearSysNormalize(EigenMatrix const& A, EigenMatrix& new_A,
 void finalizeAssembly(EigenMatrix& x);
 void finalizeAssembly(EigenVector& A);
 
+// Reduce operations for interprocess communications while using Petsc
+static inline int reduceMin(int val)
+{
+    return val;
+}
+
 }  // namespace LinAlg
 
 }  // namespace MathLib

NumLib/ODESolver/NonlinearSolver.cpp

@@ -356,6 +356,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
         BaseLib::RunTime time_assembly;
         time_assembly.start();
+        int ok = 1, g_ok = 0;
         try
         {
             sys.assemble(x, x_prev, process_id);
@@ -366,8 +367,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
                 e.what());
             error_norms_met = false;
             iteration = _maxiter;
-            break;
+            ok = 0;
         }
+        g_ok = MathLib::LinAlg::reduceMin(ok);
+        if (g_ok == 0)
+            break;
         sys.getResidual(*x[process_id], *x_prev[process_id], res);
         sys.getJacobian(J);
         INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());

NumLib/ODESolver/TimeDiscretizedODESystem.cpp

@@ -79,11 +79,18 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
 
     _b->setZero();
     _Jac->setZero();
-
-    _ode.preAssemble(t, dt, x_curr);
-    _ode.assembleWithJacobian(t, dt, x_new_timestep, x_prev, process_id, *_b,
-                              *_Jac);
-
+    try
+    {
+        _ode.preAssemble(t, dt, x_curr);
+        _ode.assembleWithJacobian(t, dt, x_new_timestep, x_prev, process_id,
+                                  *_b, *_Jac);
+    }
+    catch (AssemblyException const&)
+    {
+        MathLib::LinAlg::finalizeAssembly(*_b);
+        MathLib::LinAlg::finalizeAssembly(*_Jac);
+        throw;
+    }
     MathLib::LinAlg::finalizeAssembly(*_b);
     MathLib::LinAlg::finalizeAssembly(*_Jac);
 }

ProcessLib/SmallDeformation/SmallDeformationProcess.cpp

@@ -14,6 +14,7 @@
 
 #include "MeshLib/Utils/IntegrationPointWriter.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
+#include "NumLib/Exceptions.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
 #include "ProcessLib/Output/CellAverageAlgorithm.h"
 #include "ProcessLib/Process.h"
@@ -162,11 +163,20 @@ void SmallDeformationProcess<DisplacementDim>::
     std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
         _local_to_global_index_map.get()};
     // Call global assembler for each local assembly item.
-    GlobalExecutor::executeSelectedMemberDereferenced(
-        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x, x_prev,
-        process_id, &b, &Jac);
-
+    try
+    {
+        GlobalExecutor::executeSelectedMemberDereferenced(
+            _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+            _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x,
+            x_prev, process_id, &b, &Jac);
+    }
+    catch (NumLib::AssemblyException const&)
+    {
+        transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
+                                          *_nodal_forces,
+                                          std::negate<double>());
+        throw;
+    }
     transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
                                       *_nodal_forces, std::negate<double>());