Fix ptscotch test

.gitignore

@@ -110,6 +110,36 @@ build
 framework/contrib/libtorch
 test.cache
 vgcore*
+test/moose_test.yaml
+modules/chemical_reactions/chemical_reactions.yaml
+modules/combined/combined.yaml
+modules/contact/contact.yaml
+modules/electromagnetics/electromagnetics.yaml
+modules/external_petsc_solver/external_petsc_solver.yaml
+modules/fluid_properties/fluid_properties.yaml
+modules/fsi/fsi.yaml
+modules/functional_expansion_tools/functional_expansion_tools.yaml
+modules/geochemistry/geochemistry.yaml
+modules/heat_transfer/heat_transfer.yaml
+modules/level_set/level_set.yaml
+modules/misc/misc.yaml
+modules/module_loader/module_loader.yaml
+modules/navier_stokes/navier_stokes.yaml
+modules/optimization/optimization.yaml
+modules/peridynamics/peridynamics.yaml
+modules/phase_field/phase_field.yaml
+modules/porous_flow/porous_flow.yaml
+modules/ray_tracing/ray_tracing.yaml
+modules/rdg/rdg.yaml
+modules/reactor/reactor.yaml
+modules/richards/richards.yaml
+modules/scalar_transport/scalar_transport.yaml
+modules/solid_mechanics/solid_mechanics.yaml
+modules/solid_properties/solid_properties.yaml
+modules/stochastic_tools/stochastic_tools.yaml
+modules/subchannel/subchannel.yaml
+modules/thermal_hydraulics/thermal_hydraulics.yaml
+modules/xfem/xfem.yaml
 
 # Allow certain files in gold directories
 !**/gold/**/*.e

framework/include/partitioner/BlockWeightedPartitioner.h

@@ -26,7 +26,7 @@ class BlockWeightedPartitioner : public PetscExternalPartitioner
 
   virtual std::unique_ptr<Partitioner> clone() const override;
 
-  virtual dof_id_type computeElementWeight(Elem & elm) override;
+  virtual dof_id_type computeElementWeight(Elem & elm) const override;
 
   /**
    * Fills _blocks_to_weights before performing the partition

framework/include/partitioner/PetscExternalPartitioner.h

@@ -28,17 +28,17 @@ class PetscExternalPartitioner : public MoosePartitioner
 
   virtual std::unique_ptr<Partitioner> clone() const override;
 
-  virtual dof_id_type computeElementWeight(Elem & elm);
+  virtual dof_id_type computeElementWeight(Elem & elm) const;
 
-  virtual dof_id_type computeSideWeight(Elem & elem, unsigned int side);
+  virtual dof_id_type computeSideWeight(Elem & elem, unsigned int side) const;
 
   using Partitioner::partition;
 
   virtual void partition(MeshBase & mesh, const unsigned int n) override;
 
-  bool applySideWeight() { return _apply_side_weight; }
+  bool applySideWeight() const { return _apply_side_weight; }
 
-  bool applyElementEeight() { return _apply_element_weight; }
+  bool applyElementWeight() const { return _apply_element_weight; }
 
   static void partitionGraph(const Parallel::Communicator & comm,
                              const std::vector<std::vector<dof_id_type>> & graph,
@@ -52,7 +52,7 @@ class PetscExternalPartitioner : public MoosePartitioner
   /**
    * Called immediately before partitioning
    */
-  virtual void initialize(MeshBase & /* mesh */){};
+  virtual void initialize(MeshBase & /* mesh */) {}
 
 protected:
   virtual void _do_partition(MeshBase & mesh, const unsigned int n) override;

framework/src/loops/FlagElementsThread.C

@@ -85,11 +85,11 @@ FlagElementsThread::onElement(const Elem * elem)
       dof_value = current_local_solution(dof_number);
     }
 
-    // round() is a C99 function, it is not located in the std:: namespace.
-    marker_value = static_cast<Marker::MarkerValue>(round(dof_value));
+    const auto int_value = std::lround(dof_value);
+    marker_value = static_cast<Marker::MarkerValue>(int_value);
 
     // Make sure we aren't masking an issue in the Marker system by rounding its values.
-    if (std::abs(marker_value - dof_value) > TOLERANCE * TOLERANCE)
+    if (std::abs(int_value - dof_value) > TOLERANCE * TOLERANCE)
       mooseError("Invalid Marker value detected: ", dof_value);
   }
 

framework/src/meshgenerators/DistributedRectilinearMeshGenerator.C

@@ -1086,6 +1086,10 @@ DistributedRectilinearMeshGenerator::buildCube(UnstructuredMesh & mesh,
     for (auto neighbor : neighbors)
       if (neighbor != Elem::invalid_id)
         row.push_back(neighbor);
+    // PETSc documentation requires that rows be sorted. We've seen this matter for
+    // ptscotch, e.g. it crashes when the rows are not sorted even though its developers say in
+    // theory this shouldn't matter but in practice they admit it might.
+    std::sort(row.begin(), row.end());
   }
 
   // Partition the distributed graph

framework/src/partitioner/BlockWeightedPartitioner.C

@@ -67,13 +67,14 @@ BlockWeightedPartitioner::initialize(MeshBase & mesh)
 
   // Setup the block -> weight map for use in computeElementWeight
   _blocks_to_weights.reserve(_weights.size());
-  for (MooseIndex(block_ids.size()) i = 0; i < block_ids.size(); i++)
+  for (const auto i : index_range(block_ids))
     _blocks_to_weights[block_ids[i]] = _weights[i];
 }
 
 dof_id_type
-BlockWeightedPartitioner::computeElementWeight(Elem & elem)
+BlockWeightedPartitioner::computeElementWeight(Elem & elem) const
 {
-  mooseAssert(_blocks_to_weights.count(elem.subdomain_id()), "Missing weight for block");
-  return _blocks_to_weights[elem.subdomain_id()];
+  auto it = _blocks_to_weights.find(elem.subdomain_id());
+  mooseAssert(it != _blocks_to_weights.end(), "Missing weight for block");
+  return it->second;
 }

framework/src/partitioner/PetscExternalPartitioner.C

@@ -121,6 +121,9 @@ PetscExternalPartitioner::_do_partition(MeshBase & mesh, const unsigned int n_pa
   // Call libmesh to build the dual graph of mesh
   build_graph(mesh);
   num_local_elems = _dual_graph.size();
+  // PETSc requires that each row be sorted
+  for (auto & row : _dual_graph)
+    std::sort(row.begin(), row.end());
 
   elem_weights.clear();
   if (_apply_element_weight)
@@ -149,11 +152,13 @@ PetscExternalPartitioner::_do_partition(MeshBase & mesh, const unsigned int n_pa
 
   local_elem_id = 0;
   nj = 0;
+  std::map<dof_id_type, dof_id_type> global_index_to_weight;
   for (auto & row : _dual_graph)
   {
     mooseAssert(local_elem_id < static_cast<dof_id_type>(_local_id_to_elem.size()),
                 "Local element id " << local_elem_id << " is not smaller than "
                                     << _local_id_to_elem.size());
+    global_index_to_weight.clear();
     auto elem = _local_id_to_elem[local_elem_id];
     unsigned int n_neighbors = 0;
 
@@ -165,9 +170,8 @@ PetscExternalPartitioner::_do_partition(MeshBase & mesh, const unsigned int n_pa
       if (neighbor != nullptr && neighbor->active())
       {
         if (_apply_side_weight)
-          side_weights[nj] = computeSideWeight(*elem, side);
-
-        nj++;
+          global_index_to_weight.emplace(libmesh_map_find(_global_index_by_pid_map, neighbor->id()),
+                                         computeSideWeight(*elem, side));
         n_neighbors++;
       }
 
@@ -176,6 +180,12 @@ PetscExternalPartitioner::_do_partition(MeshBase & mesh, const unsigned int n_pa
     if (n_neighbors != row.size())
       mooseError(
           "Cannot construct dual graph correctly since the number of neighbors is inconsistent");
+    if (_apply_side_weight)
+    {
+      mooseAssert(global_index_to_weight.size() == row.size(), "These must match");
+      for (const auto [_, weight] : global_index_to_weight)
+        side_weights[nj++] = weight;
+    }
 
     local_elem_id++;
   }
@@ -231,7 +241,7 @@ PetscExternalPartitioner::partitionGraph(const Parallel::Communicator & comm,
   // Fill up adjacency
   i = 0;
   for (auto & row : graph)
-    for (auto elem : row)
+    for (const auto elem : row)
       adjncy[i++] = elem;
 
   // If there are no neighbors at all, no side weights should be proivded
@@ -258,16 +268,15 @@ PetscExternalPartitioner::partitionGraph(const Parallel::Communicator & comm,
       MatCreateMPIAdj(comm.get(), num_local_elems, num_elems, xadj, adjncy, values, &dual));
 
   LibmeshPetscCallA(comm.get(), MatPartitioningCreate(comm.get(), &part));
-#if !PETSC_VERSION_LESS_THAN(3, 12, 3)
-  LibmeshPetscCallA(comm.get(), MatPartitioningSetUseEdgeWeights(part, PETSC_TRUE));
-#endif
+  if (values)
+    // This should only be set to true if the adjacency matrix has valid edge weights (which
+    // correspond to \p values)
+    LibmeshPetscCallA(comm.get(), MatPartitioningSetUseEdgeWeights(part, PETSC_TRUE));
   LibmeshPetscCallA(comm.get(), MatPartitioningSetAdjacency(part, dual));
 
   if (!num_local_elems)
-  {
     mooseAssert(!elem_weights.size(),
                 "No element weights should be provided since there are no elements at all");
-  }
 
   // Handle element weights
   if (elem_weights.size())
@@ -286,12 +295,8 @@ PetscExternalPartitioner::partitionGraph(const Parallel::Communicator & comm,
   }
 
   LibmeshPetscCallA(comm.get(), MatPartitioningSetNParts(part, num_parts));
-#if PETSC_VERSION_LESS_THAN(3, 9, 2)
   mooseAssert(part_package != "party", "PETSc-3.9.3 or higher is required for using party");
-#endif
-#if PETSC_VERSION_LESS_THAN(3, 9, 0)
   mooseAssert(part_package != "chaco", "PETSc-3.9.0 or higher is required for using chaco");
-#endif
   LibmeshPetscCallA(comm.get(), MatPartitioningSetType(part, part_package.c_str()));
   if (part_package == "hierarch")
     LibmeshPetscCallA(comm.get(),
@@ -313,13 +318,13 @@ PetscExternalPartitioner::partitionGraph(const Parallel::Communicator & comm,
 }
 
 dof_id_type
-PetscExternalPartitioner::computeElementWeight(Elem & /*elem*/)
+PetscExternalPartitioner::computeElementWeight(Elem & /*elem*/) const
 {
   return 1;
 }
 
 dof_id_type
-PetscExternalPartitioner::computeSideWeight(Elem & /*elem*/, unsigned int /*side*/)
+PetscExternalPartitioner::computeSideWeight(Elem & /*elem*/, unsigned int /*side*/) const
 {
   return 1;
 }

framework/src/preconditioners/StaticCondensationFieldSplitPreconditioner.C

@@ -67,7 +67,7 @@ StaticCondensationFieldSplitPreconditioner::setupDM()
   // Create and set up the DM that will consume the split options and deal with block matrices.
   auto & petsc_solver = cast_ref<PetscLinearSolver<Number> &>(scSysMat().reduced_system_solver());
   auto ksp = petsc_solver.ksp();
-  // if there exists a DMMoose object, not to recreate a new one
+  // if there exists a DMMoose object, then we do not need to recreate one
   LibmeshPetscCall(KSPGetDM(ksp, &dm));
   if (dm)
   {
@@ -77,8 +77,12 @@ StaticCondensationFieldSplitPreconditioner::setupDM()
   }
   createMooseDM(&dm);
   LibmeshPetscCall(KSPSetDM(ksp, dm));
-  // We set the operators ourselves. We do not want the DM to generate the operators
+// We set the operators ourselves. We do not want the DM to generate the operators
+#if PETSC_VERSION_LESS_THAN(3, 25, 0)
   LibmeshPetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
+#else
+  LibmeshPetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_OPERATOR, PETSC_FALSE));
+#endif
   LibmeshPetscCall(DMDestroy(&dm));
 }
 

framework/src/vectorpostprocessors/WorkBalance.C

@@ -106,9 +106,8 @@ public:
       _local_partition_hardware_id_surface_area(0),
       _this_pid(_mesh.processor_id()) // Get this once because it is expensive
   {
-    // This is required because dynamic_pointer_cast() requires an l-value
-    auto partitioner = mesh.getMesh().partitioner()->clone();
-    _petsc_partitioner = dynamic_pointer_cast<PetscExternalPartitioner>(partitioner);
+    _petsc_partitioner =
+        dynamic_cast<PetscExternalPartitioner *>(mesh.getMesh().partitioner().get());
   }
 
   WBElementLoop(WBElementLoop & x, Threads::split split)
@@ -122,14 +121,9 @@ public:
       _local_partition_surface_area(0),
       _local_num_partition_hardware_id_sides(0),
       _local_partition_hardware_id_surface_area(0),
-      _this_pid(x._this_pid)
+      _this_pid(x._this_pid),
+      _petsc_partitioner(x._petsc_partitioner)
   {
-    if (x._petsc_partitioner)
-    {
-      // This is required because dynamic_pointer_cast() requires an l-value
-      auto partitioner = x._petsc_partitioner->clone();
-      _petsc_partitioner = dynamic_pointer_cast<PetscExternalPartitioner>(partitioner);
-    }
   }
 
   virtual ~WBElementLoop() {}
@@ -146,7 +140,7 @@ public:
 
   virtual void onElement(const Elem * elem) override
   {
-    if (_petsc_partitioner && _petsc_partitioner->applyElementEeight())
+    if (_petsc_partitioner && _petsc_partitioner->applyElementWeight())
     {
       // We should change partitioner interface to take const
       // But at this point let us keep API intact
@@ -230,7 +224,7 @@ public:
 
   libMesh::ElemSideBuilder _elem_side_builder;
 
-  std::unique_ptr<PetscExternalPartitioner> _petsc_partitioner;
+  const PetscExternalPartitioner * _petsc_partitioner;
 
 private:
   bool shouldComputeInternalSide(const Elem & /*elem*/, const Elem & /*neighbor*/) const override

modules/contact/src/linesearches/PetscContactLineSearch.C

@@ -39,7 +39,6 @@ PetscContactLineSearch::lineSearch()
   Vec X, F, Y, W, G, W1;
   SNESLineSearch line_search;
   PetscReal fnorm, xnorm, ynorm, gnorm;
-  PetscBool domainerror;
   PetscReal ksp_rtol, ksp_abstol, ksp_dtol;
   PetscInt ksp_maxits;
   KSP ksp;
@@ -73,9 +72,17 @@ PetscContactLineSearch::lineSearch()
    * residual evaluation */
   _current_contact_state.clear();
   LibmeshPetscCall(SNESComputeFunction(snes, W, F));
+#if PETSC_VERSION_LESS_THAN(3, 25, 0)
+  PetscBool domainerror;
   LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
   if (domainerror)
     LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
+#else
+  SNESConvergedReason reason;
+  LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
+  if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
+    LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
+#endif
 
   LibmeshPetscCall(VecNorm(F, NORM_2, &fnorm));
   std::set<dof_id_type> contact_state_stored = _current_contact_state;
@@ -101,9 +108,16 @@ PetscContactLineSearch::lineSearch()
     LibmeshPetscCall(VecWAXPY(W1, -_contact_lambda, Y, X));
 
     LibmeshPetscCall(SNESComputeFunction(snes, W1, G));
+#if PETSC_VERSION_LESS_THAN(3, 25, 0)
     LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
     if (domainerror)
       LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
+#else
+    LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
+    if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
+      LibmeshPetscCall(
+          SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
+#endif
 
     LibmeshPetscCall(VecNorm(G, NORM_2, &gnorm));
     if (gnorm < fnorm)
@@ -130,10 +144,16 @@ PetscContactLineSearch::lineSearch()
   if (changed_w || changed_y)
   {
     LibmeshPetscCall(SNESComputeFunction(snes, W, F));
+#if PETSC_VERSION_LESS_THAN(3, 25, 0)
     LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
     if (domainerror)
       LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
-
+#else
+    LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
+    if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
+      LibmeshPetscCall(
+          SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
+#endif
     contact_state_stored.swap(_current_contact_state);
     _current_contact_state.clear();
     printContactInfo(contact_state_stored);

modules/contact/src/linesearches/PetscProjectSolutionOntoBounds.C

@@ -52,7 +52,6 @@ PetscProjectSolutionOntoBounds::lineSearch()
   Vec X, F, Y, W, G;
   SNESLineSearch line_search;
   PetscReal fnorm, xnorm, ynorm;
-  PetscBool domainerror;
   SNES snes = _solver->snes();
 
   LibmeshPetscCall(SNESGetLineSearch(snes, &line_search));
@@ -134,9 +133,17 @@ PetscProjectSolutionOntoBounds::lineSearch()
   LibmeshPetscCall(VecAssemblyEnd(W));
 
   LibmeshPetscCall(SNESComputeFunction(snes, W, F));
+#if PETSC_VERSION_LESS_THAN(3, 25, 0)
+  PetscBool domainerror;
   LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
   if (domainerror)
     LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
+#else
+  SNESConvergedReason reason;
+  LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
+  if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
+    LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
+#endif
 
   LibmeshPetscCall(VecNorm(F, NORM_2, &fnorm));