Merge branch 'bugfix-petsc-slope-case' into 'master'

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

See merge request ogs/ogs!5096

BaseLib/MPI.h

@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace BaseLib::MPI
+{
+
+#ifdef USE_PETSC
+// 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;
+}
+#else
+// Reduce operations for interprocess communications without using Petsc
+static inline int reduceMin(int val)
+{
+    return val;
+}
+#endif
+
+}  // namespace BaseLib::MPI

MathLib/LinAlg/LinAlg.h

@@ -17,6 +17,8 @@
 
 #ifdef USE_PETSC
 #include <petscsystypes.h>
+
+#include "PETSc/PETScVector.h"
 #endif
 
 namespace MathLib

NumLib/ODESolver/NonlinearSolver.cpp

@@ -15,6 +15,7 @@
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "ConvergenceCriterion.h"
 #include "MathLib/LinAlg/LinAlg.h"
@@ -356,6 +357,8 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
         BaseLib::RunTime time_assembly;
         time_assembly.start();
+        int mpi_rank_assembly_ok = 1;
+        int mpi_all_assembly_ok = 0;
         try
         {
             sys.assemble(x, x_prev, process_id);
@@ -366,6 +369,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
                 e.what());
             error_norms_met = false;
             iteration = _maxiter;
+            mpi_rank_assembly_ok = 0;
+        }
+        mpi_all_assembly_ok = BaseLib::MPI::reduceMin(mpi_rank_assembly_ok);
+        if (mpi_all_assembly_ok == 0)
+        {
             break;
         }
         sys.getResidual(*x[process_id], *x_prev[process_id], res);

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>());
 

ProcessLib/SmallDeformation/Tests.cmake

@@ -112,6 +112,9 @@ if (OGS_USE_MFRONT)
         PROPERTIES PASS_REGULAR_EXPRESSION "The nonlinear solver failed in time step #.* at t = 5.37.* s for process #0"
     )
     OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloan/load_test_mc.prj)
+    if (OGS_USE_MPI)
+        OgsTest(WRAPPER mpirun -np 4 PROJECTFILE Mechanics/MohrCoulombAbboSloan/PetscMpi/slope_hexa.prj)
+    endif()
     OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_1e0_47.prj)
     OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_1e0_47.prj)
     OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_1e0_47.prj)

Tests/Data/Mechanics/MohrCoulombAbboSloan/PetscMpi/README

@@ -0,0 +1,8 @@
+Instruction to generate this test case:
+
+gmsh -3 -format msh2 hexa_slope.geo
+msh2vtu -r hexa_slope.msh
+partmesh -i hexa_slope_domain.vtu --ogs2metis
+partmesh -n 4 -m -i hexa_slope_domain.vtu -- hexa_slope_physical_group*.vtu
+mpirun -n 4 ./soft/build/release-petsc/bin/ogs slope_hexa.prj
+

Tests/Data/Mechanics/MohrCoulombAbboSloan/PetscMpi/hexa_slope.geo

@@ -0,0 +1,72 @@
+lc = 3.0;
+nh = 8;
+//+
+Point(1) = {0, 0, 0, lc};
+//+
+Point(2) = {18, 0, 0, lc};
+//+
+Point(3) = {18, 2, 0, lc};
+//+
+Point(4) = {15, 2, 0, lc};
+//+
+Point(5) = {6, 6.5, 0, lc};
+//+
+Point(6) = {0, 6.5, 0, lc};
+//+
+Point(7) = {15, 0, 0, lc};
+//+
+Point(8) = {6, 0, 0, lc};
+
+//+
+Line(1) = {1, 8};
+//+
+Line(2) = {8, 7};
+//+
+Line(3) = {7, 2};
+//+
+Line(4) = {2, 3};
+//+
+Line(5) = {3, 4};
+//+
+Line(6) = {4, 5};
+//+
+Line(7) = {5, 6};
+//+
+Line(8) = {6, 1};
+//+
+Line(9) = {5, 8};
+//+
+Line(10) = {4, 7};
+//+
+Curve Loop(1) = {8, 1, -9, 7};
+//+
+Plane Surface(1) = {1};
+//+
+Curve Loop(2) = {9, 2, -10, 6};
+//+
+Plane Surface(2) = {2};
+//+
+Curve Loop(3) = {10, 3, 4, 5};
+//+
+Plane Surface(3) = {3};
+
+Extrude {0, 0, 20} {
+	Surface{1,2,3};
+}
+
+Transfinite Line "*" = nh;
+Transfinite Line {17,18,22,26,44,48,66,70} = 3*nh;
+Transfinite Surface "*";
+Recombine Surface "*";
+Transfinite Volume "*";
+//+
+//+
+Physical Surface("cote", 77) = {32, 54, 76, 1, 2, 3};
+//+
+Physical Surface("bas", 78) = {23, 45, 67};
+//+
+Physical Surface("derriere", 79) = {19};
+//+
+Physical Surface("devant", 80) = {71};
+//+
+Physical Volume("bloc", 81) = {2, 1, 3};