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};