Merge pull request #27810 from oanaoana/ghosting

Ghosting boundaries

framework/doc/content/syntax/Mesh/splitting.md

@@ -51,9 +51,20 @@ It is noted that if there are mesh meta data generated by mesh generators, these
 will be written to a binary file under the generated directory that can be loaded when using
 split meshes.
 
+!alert note
+The mesh splitter commands do not work with a [distributed mesh](syntax/Mesh/index.md#replicated-and-distributed-mesh). You can only split starting from a replicated mesh.
+
 ## Using Split Meshes
 
-To use a mesh split configuration use the `--use-split` flag (which takes no arguments):
+There are two ways of referring to distributed meshes in MOOSE, as split or as distributed. A split mesh is a distributed mesh. A split mesh is read from a pre-split mesh file (from a [Checkpoint.md] or using the instructions above). A distributed mesh can be generated in parallel or loaded entirely replicated then distributed.
+
+### Loading a split mesh using [FileMeshGenerator.md]
+
+A split mesh, i.e. a `.cpr` file, may be loaded using a [FileMeshGenerator.md] if access to the mesh is needed before the full simulation case is set-up, notably for additional mesh generation with [mesh generators](syntax/Mesh/index.md#mesh-generators). In fact this approach is recommended; it leads to correct element ghosting whereas the approach we introduce in the following section may not.
+
+### Loading a split mesh using `--use-split`
+
+Alternatively, a previously split mesh can be loaded using the `--use-split` flag (which takes no arguments):
 
 ```
 $ mpiexec -n 42 moose-app-opt -i your_input.i --use-split
@@ -72,8 +83,5 @@ not use a file-based mesh, you will need to specify a split mesh file name using
 $ mpiexec -n 42 moose-app-opt -i your_input.i --use-split --split-file foo.cpr
 ```
 
-!alert note
-The mesh splitter commands do not work with DistributedMesh. You must only split with a ReplicatedMesh.
-
 !alert warning
-Using `--use-split` skips the `[Mesh]` block of the input file. Split mesh may be loaded using a [FileMeshGenerator.md] if further mesh generation is to be performed on the split mesh.
+The `--use-split` option skips entirely the `[Mesh]` block of the input file, and as such some data structure initializations required if further mesh operations are necessary at the simulation start-up stage. 

modules/heat_transfer/include/actions/ThermalContactAction.h

@@ -28,6 +28,7 @@ class ThermalContactAction : public Action
   virtual void addDiracKernels();
   virtual void addMaterials();
   virtual void addSecondaryFluxVector();
+  virtual void addRelationshipManagers(Moose::RelationshipManagerType input_rm) override;
 
   const bool _quadrature;
   const MooseEnum _order;
@@ -37,4 +38,5 @@ class ThermalContactAction : public Action
 
   /// Primary/Secondary boundary name pairs for thermal contact
   const std::vector<std::pair<BoundaryName, BoundaryName>> _boundary_pairs;
+  using Action::addRelationshipManagers;
 };

modules/heat_transfer/include/relationshipmanagers/GhostBoundary.h

@@ -0,0 +1,63 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+// Framework includes
+#include "AutomaticMortarGeneration.h"
+#include "RelationshipManager.h"
+
+// libMesh includes
+#include "libmesh/mesh_base.h"
+
+using libMesh::boundary_id_type;
+using libMesh::CouplingMatrix;
+using libMesh::Elem;
+using libMesh::GhostingFunctor;
+using libMesh::MeshBase;
+using libMesh::processor_id_type;
+
+/**
+ * GhostBoundary is used to ghost elements on a boundary. It is
+ * useful when non-local elements are required for geometric searches or when
+ * residual objects such as \p GapHeatTransfer require access into nonlocal entries
+ * in the solution vector corresponding to degrees of freedom from the boundary
+ * elements
+ */
+class GhostBoundary : public RelationshipManager
+{
+public:
+  GhostBoundary(const InputParameters &);
+
+  GhostBoundary(const GhostBoundary & other);
+
+  static InputParameters validParams();
+
+  virtual void operator()(const MeshBase::const_element_iterator & /*range_begin*/,
+                          const MeshBase::const_element_iterator & /*range_end*/,
+                          processor_id_type p,
+                          map_type & coupled_elements) override;
+
+  virtual std::unique_ptr<GhostingFunctor> clone() const override;
+
+  virtual void redistribute() override { this->mesh_reinit(); }
+
+  std::string getInfo() const override;
+
+  virtual bool operator>=(const RelationshipManager & other) const override;
+
+protected:
+  virtual void internalInitWithMesh(const MeshBase &) override;
+
+  /// The boundary for which we will ghost elements
+  const std::vector<BoundaryName> & _boundary_name;
+
+  /// null matrix for generating full variable coupling
+  const CouplingMatrix * const _null_mat = nullptr;
+};

modules/heat_transfer/src/actions/ThermalContactAction.C

@@ -11,6 +11,7 @@
 
 #include "AddVariableAction.h"
 #include "FEProblem.h"
+#include "GapHeatTransfer.h"
 #include "libmesh/string_to_enum.h"
 #include "GapConductance.h"
 #include "GapConductanceConstant.h"
@@ -393,3 +394,22 @@ ThermalContactAction::addSecondaryFluxVector()
   if (!_problem->isSNESMFReuseBaseSetbyUser())
     _problem->setSNESMFReuseBase(false, false);
 }
+
+void
+ThermalContactAction::addRelationshipManagers(Moose::RelationshipManagerType input_rm_type)
+{
+  if (!_quadrature)
+    return;
+
+  for (const auto & contact_pair : _boundary_pairs)
+  {
+    const auto & object_name = getParam<std::string>("type");
+    auto params = _factory.getValidParams(object_name);
+    params.applyParameters(parameters());
+
+    params.set<BoundaryName>("paired_boundary") = contact_pair.first;
+    params.set<bool>("use_displaced_mesh") = true;
+    params.set<std::vector<BoundaryName>>("boundary") = {contact_pair.second};
+    addRelationshipManagers(input_rm_type, params);
+  }
+}

modules/heat_transfer/src/bcs/GapHeatTransfer.C

@@ -16,6 +16,7 @@
 #include "MooseVariable.h"
 #include "PenetrationLocator.h"
 #include "SystemBase.h"
+#include "GhostBoundary.h"
 
 #include "libmesh/string_to_enum.h"
 
@@ -78,6 +79,16 @@ GapHeatTransfer::validParams()
   params.addCoupledVar("gap_distance", "Distance across the gap");
   params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
 
+  params.addRelationshipManager(
+      "GhostBoundary",
+      Moose::RelationshipManagerType::GEOMETRIC,
+      [](const InputParameters & obj_params, InputParameters & rm_params)
+      {
+        auto & boundary = rm_params.set<std::vector<BoundaryName>>("boundary");
+        boundary = obj_params.get<std::vector<BoundaryName>>("boundary");
+        boundary.push_back(obj_params.get<BoundaryName>("paired_boundary"));
+      });
+
   return params;
 }
 

modules/heat_transfer/src/relationshipmanagers/GhostBoundary.C

@@ -0,0 +1,111 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+// App includes
+#include "GhostBoundary.h"
+#include "Executioner.h"
+#include "FEProblemBase.h"
+#include "MooseApp.h"
+
+// libMesh includes
+#include "libmesh/elem.h"
+#include "libmesh/mesh_base.h"
+#include "libmesh/boundary_info.h"
+
+registerMooseObject("HeatTransferApp", GhostBoundary);
+
+using namespace libMesh;
+
+InputParameters
+GhostBoundary::validParams()
+{
+  InputParameters params = RelationshipManager::validParams();
+  params.addRequiredParam<std::vector<BoundaryName>>("boundary",
+                                                     "The name of the primary boundary sideset.");
+  params.addClassDescription("This class constructs a relationship manager system' "
+                             "to communicate ghost elements on a boundary.");
+  return params;
+}
+
+GhostBoundary::GhostBoundary(const InputParameters & params)
+  : RelationshipManager(params), _boundary_name(getParam<std::vector<BoundaryName>>("boundary"))
+{
+}
+
+GhostBoundary::GhostBoundary(const GhostBoundary & other)
+  : RelationshipManager(other), _boundary_name(other._boundary_name)
+{
+}
+
+void
+GhostBoundary::internalInitWithMesh(const MeshBase &)
+{
+}
+
+std::string
+GhostBoundary::getInfo() const
+{
+  std::ostringstream oss;
+  oss << "GhostBoundary";
+  return oss.str();
+}
+
+void
+GhostBoundary::operator()(const MeshBase::const_element_iterator & /*range_begin*/,
+                          const MeshBase::const_element_iterator & /*range_end*/,
+                          const processor_id_type p,
+                          map_type & coupled_elements)
+{
+  // We ask the user to pass boundary names instead of ids to our constraint object.  However, We
+  // are unable to get the boundary ids from boundary names until we've attached the MeshBase object
+  // to the MooseMesh
+  const bool generating_mesh = !_moose_mesh->getMeshPtr();
+  const auto boundary_ids = generating_mesh ? std::vector<BoundaryID>{Moose::INVALID_BOUNDARY_ID}
+                                            : _moose_mesh->getBoundaryIDs(_boundary_name);
+
+  for (const Elem * const elem : _mesh->active_element_ptr_range())
+  {
+    if (generating_mesh)
+    { // We are still generating the mesh, so it's possible we don't even have the right boundary
+      // ids created yet! So we actually ghost all boundary elements and all lower dimensional
+      // elements who have parents on a boundary
+      if (elem->on_boundary())
+        coupled_elements.insert(std::make_pair(elem, _null_mat));
+    }
+    else
+    {
+      // We've finished generating our mesh so we can be selective and only ghost elements lying on
+      // our boundary
+      const BoundaryInfo & binfo = _mesh->get_boundary_info();
+      for (auto side : elem->side_index_range())
+        for (auto boundary_id : boundary_ids)
+          if ((elem->processor_id() != p) && (binfo.has_boundary_id(elem, side, boundary_id)))
+          {
+            coupled_elements.insert(std::make_pair(elem, _null_mat));
+            goto countBreak;
+          }
+    countBreak:;
+    }
+  }
+}
+
+bool
+GhostBoundary::operator>=(const RelationshipManager & other) const
+{
+  if (auto asoi = dynamic_cast<const GhostBoundary *>(&other))
+    if (_boundary_name == asoi->_boundary_name && baseGreaterEqual(*asoi))
+      return true;
+  return false;
+}
+
+std::unique_ptr<GhostingFunctor>
+GhostBoundary::clone() const
+{
+  return _app.getFactory().copyConstruct(*this);
+}

modules/heat_transfer/test/tests/heat_conduction/2d_quadrature_gap_heat_transfer/perfect_split.i

@@ -0,0 +1,64 @@
+[Mesh]
+  [fmg]
+    type = FileMeshGenerator
+    file = 'perfect.cpa.gz'
+  []
+  parallel_type = distributed
+[]
+
+[Variables]
+  [temp]
+  []
+[]
+
+[Kernels]
+  [hc]
+    type = HeatConduction
+    variable = temp
+  []
+[]
+
+[BCs]
+  [left]
+    type = DirichletBC
+    variable = temp
+    boundary = leftleft
+    value = 300
+  []
+  [right]
+    type = DirichletBC
+    variable = temp
+    boundary = rightright
+    value = 400
+  []
+[]
+
+[ThermalContact]
+  [left_to_right]
+    secondary = leftright
+    quadrature = true
+    primary = rightleft
+    emissivity_primary = 0
+    emissivity_secondary = 0
+    variable = temp
+    type = GapHeatTransfer
+  []
+[]
+
+[Materials]
+  [hcm]
+    type = HeatConductionMaterial
+    block = 'left right'
+    specific_heat = 1
+    thermal_conductivity = 1
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = 'PJFNK'
+[]
+
+[Outputs]
+  exodus = true
+[]

modules/heat_transfer/test/tests/heat_conduction/2d_quadrature_gap_heat_transfer/tests

@@ -76,4 +76,30 @@
     design = "/actions/ThermalContactAction.md"
     issues = "#13043"
   [../]
+
+  [perfect_prereq]
+    type = 'CheckFiles'
+    prereq = perfect
+    input = 'perfect.i'
+    cli_args = '--split-mesh=2'
+    mesh_mode = 'REPLICATED'
+    check_files = 'perfect.cpa.gz/2/header.gz perfect.cpa.gz/2/split-2-0.gz perfect.cpa.gz/2/split-2-1.gz'
+    recover = false
+    requirement = "The system shall generate a parallel mesh split across 2 processes."
+    design = "Mesh/splitting.md"
+    issues = "#27203"
+  []
+
+  [perfect_split_mesh]
+    type = 'Exodiff'
+    prereq = perfect_prereq
+    input = 'perfect_split.i'
+    exodiff = 'perfect_out.e'
+    cli_args = 'Outputs/file_base=perfect_out'
+    requirement = "The system shall read in parallel a mesh split across 2 processes."
+    design = "Mesh/splitting.md"
+    min_parallel = 2
+    max_parallel = 2
+    issues = "#27203"
+  []
 []