Setup periodic boundary conditions in a scalable way

when mesh is generated in parallel

Refs idaholab#15501

framework/include/mesh/MooseMesh.h

@@ -630,6 +630,11 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf
    */
   void ghostGhostedBoundaries();
 
+  /**
+   *
+   */
+  void needGhostGhostedBoundaries(bool needghost) { _need_ghost_ghosted_boundaries = needghost; }
+
   /**
    * Getter for the patch_size parameter.
    */
@@ -1419,6 +1424,11 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf
 
   /// Set of elements ghosted by ghostGhostedBoundaries
   std::set<Elem *> _ghost_elems_from_ghost_boundaries;
+
+  /// Some mesh generator such as DistributedRectilinearMeshGenerator already
+  /// make everything ready. We do not need to gather all boundaries to
+  /// every single processor
+  bool _need_ghost_ghosted_boundaries;
 };
 
 /**

framework/src/mesh/MooseMesh.C

@@ -256,7 +256,8 @@ MooseMesh::MooseMesh(const InputParameters & parameters)
     _init_timer(registerTimedSection("init", 2)),
     _read_recovered_mesh_timer(registerTimedSection("readRecoveredMesh", 2)),
     _ghost_ghosted_boundaries_timer(registerTimedSection("GhostGhostedBoundaries", 3)),
-    _need_delete(false)
+    _need_delete(false),
+    _need_ghost_ghosted_boundaries(true)
 {
   if (isParamValid("ghosting_patch_size") && (_patch_update_strategy != Moose::Iteration))
     mooseError("Ghosting patch size parameter has to be set in the mesh block "
@@ -317,7 +318,8 @@ MooseMesh::MooseMesh(const MooseMesh & other_mesh)
     _init_timer(registerTimedSection("init", 2)),
     _read_recovered_mesh_timer(registerTimedSection("readRecoveredMesh", 2)),
     _ghost_ghosted_boundaries_timer(registerTimedSection("GhostGhostedBoundaries", 3)),
-    _need_delete(other_mesh._need_delete)
+    _need_delete(other_mesh._need_delete),
+    _need_ghost_ghosted_boundaries(other_mesh._need_ghost_ghosted_boundaries)
 {
   // Note: this calls BoundaryInfo::operator= without changing the
   // ownership semantics of either Mesh's BoundaryInfo object.
@@ -1569,6 +1571,17 @@ MooseMesh::detectPairedSidesets()
     }
   }
 
+  if (_use_distributed_mesh)
+    for (unsigned side_dim = 0; side_dim < dim; ++side_dim)
+    {
+      _communicator.set_union(minus_x_ids[side_dim]);
+      _communicator.set_union(plus_x_ids[side_dim]);
+      _communicator.set_union(minus_y_ids[side_dim]);
+      _communicator.set_union(plus_y_ids[side_dim]);
+      _communicator.set_union(minus_z_ids[side_dim]);
+      _communicator.set_union(plus_z_ids[side_dim]);
+    }
+
   for (unsigned side_dim = 0; side_dim < dim; ++side_dim)
   {
     // If unique pairings were found, fill up the _paired_boundary data
@@ -2559,7 +2572,7 @@ void
 MooseMesh::ghostGhostedBoundaries()
 {
   // No need to do this if using a serial mesh
-  if (!_use_distributed_mesh)
+  if (!_use_distributed_mesh || !_need_ghost_ghosted_boundaries)
     return;
 
   TIME_SECTION(_ghost_ghosted_boundaries_timer);

framework/src/meshgenerators/DistributedRectilinearMeshGenerator.C

@@ -134,9 +134,19 @@ DistributedRectilinearMeshGenerator::getNeighbors<Edge2>(const dof_id_type nx,
                                                          const dof_id_type /*j*/,
                                                          const dof_id_type /*k*/,
                                                          std::vector<dof_id_type> & neighbors,
-                                                         const bool /*corner*/)
+                                                         const bool corner)
 
 {
+  if (corner)
+  {
+    // The elements on the other side are useful for
+    // periodic boundary conditions
+    neighbors[0] = (i - 1 + nx) % nx;
+    neighbors[1] = (i + 1 + nx) % nx;
+
+    return;
+  }
+
   neighbors[0] = i - 1;
   neighbors[1] = i + 1;
 
@@ -169,26 +179,17 @@ DistributedRectilinearMeshGenerator::getGhostNeighbors<Edge2>(const dof_id_type
                                                               const MeshBase & mesh,
                                                               std::set<dof_id_type> & ghost_elems)
 {
-  auto & boundary_info = mesh.get_boundary_info();
-
   std::vector<dof_id_type> neighbors(2);
 
   for (auto elem_ptr : mesh.element_ptr_range())
   {
-    for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
-    {
-      // No current neighbor
-      if (!elem_ptr->neighbor_ptr(s))
-      {
-        // Not on a boundary
-        if (!boundary_info.n_boundary_ids(elem_ptr, s))
-        {
-          getNeighbors<Edge2>(nx, 0, 0, elem_ptr->id(), 0, 0, neighbors, false);
+    auto elem_id = elem_ptr->id();
 
-          ghost_elems.insert(neighbors[s]);
-        }
-      }
-    }
+    getNeighbors<Edge2>(nx, 0, 0, elem_id, 0, 0, neighbors, true);
+
+    for (auto neighbor : neighbors)
+      if (neighbor != Elem::invalid_id && !mesh.query_elem_ptr(neighbor))
+        ghost_elems.insert(neighbor);
   }
 }
 
@@ -337,11 +338,13 @@ DistributedRectilinearMeshGenerator::getNeighbors<Quad4>(const dof_id_type nx,
   if (corner)
   {
     // libMesh dof_id_type looks like unsigned int
+    // We add one layer of point neighbors
+    // If the elems are on on boundary, the elems on the other side
+    // are included in case there are some periodic boundary conditions in simulation
     unsigned int nnb = 0;
     for (unsigned int ii = 0; ii <= 2; ii++)
       for (unsigned int jj = 0; jj <= 2; jj++)
-        if ((i + ii >= 1) && (i + ii <= nx) && (j + jj >= 1) && (j + jj <= ny))
-          neighbors[nnb++] = elemId<Quad4>(nx, 0, i + ii - 1, j + jj - 1, 0);
+        neighbors[nnb++] = elemId<Quad4>(nx, 0, (i + ii - 1 + nx) % nx, (j + jj - 1 + ny) % ny, 0);
 
     return;
   }
@@ -600,13 +603,14 @@ DistributedRectilinearMeshGenerator::getNeighbors<Hex8>(const dof_id_type nx,
 
   if (corner)
   {
+    // We collect one layer of point neighbors
+    // When the elements are on a boundary, we add the elements on the other side
     unsigned int nnb = 0;
     for (unsigned int ii = 0; ii <= 2; ii++)
       for (unsigned int jj = 0; jj <= 2; jj++)
         for (unsigned int kk = 0; kk <= 2; kk++)
-          if ((i + ii >= 1) && (i + ii <= nx) && (j + jj >= 1) && (j + jj <= ny) && (k + kk >= 1) &&
-              (k + kk <= nz))
-            neighbors[nnb++] = elemId<Hex8>(nx, ny, i + ii - 1, j + jj - 1, k + kk - 1);
+          neighbors[nnb++] = elemId<Hex8>(
+              nx, ny, (i + ii - 1 + nx) % nx, (j + jj - 1 + ny) % ny, (k + kk - 1 + nz) % nz);
 
     return;
   }
@@ -1015,6 +1019,11 @@ DistributedRectilinearMeshGenerator::buildCube(UnstructuredMesh & mesh,
   bool old_allow_renumbering = mesh.allow_renumbering();
   mesh.allow_renumbering(false);
   mesh.allow_find_neighbors(false);
+  // No, I do not want to remove anything
+  // Periodic boundary conditions will be added late
+  // We should stop libMesh from deleting the valuable boundary elements
+  // paired with the local elements
+  mesh.allow_remote_element_removal(false);
   mesh.prepare_for_use();
 
   // But we'll want to at least find neighbors after any future mesh changes!
@@ -1030,6 +1039,11 @@ DistributedRectilinearMeshGenerator::generate()
 {
   // DistributedRectilinearMeshGenerator always generates a distributed mesh
   _mesh->setParallelType(MooseMesh::ParallelType::DISTRIBUTED);
+  // We will setup boundaries accordingly
+  // We do not need to call ghostGhostedBoundaries in which allgather_packed_range
+  // is unscalable.
+  // ghostGhostedBoundaries will gather all boundaries to every single processor
+  _mesh->needGhostGhostedBoundaries(false);
   auto mesh = _mesh->buildMeshBaseObject(MooseMesh::ParallelType::DISTRIBUTED);
 
   MooseEnum elem_type_enum = getParam<MooseEnum>("elem_type");

test/tests/bcs/dmg_periodic/dmg_periodic_bc.i

@@ -0,0 +1,85 @@
+[Mesh]
+  [dmg]
+    type = DistributedRectilinearMeshGenerator
+    dim = 2
+    nx = 40
+    ny = 40
+    nz = 0
+    xmax = 40
+    ymax = 40
+    zmax = 0
+  []
+  // We need this for recovering
+  parallel_type = distributed
+[]
+
+[Variables]
+  [./u]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[AuxVariables]
+  [./periodic_dist]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+  [./pid]
+    order = CONSTANT
+    family = monomial
+  []
+[]
+
+[AuxKernels]
+  [./pidaux]
+    type = ProcessorIDAux
+    variable = pid
+  [../]
+[]
+
+[Kernels]
+  [./diff]
+    type = Diffusion
+    variable = u
+  [../]
+
+  [./forcing]
+    type = GaussContForcing
+    variable = u
+  [../]
+
+  [./dot]
+    type = TimeDerivative
+    variable = u
+  [../]
+[]
+
+[AuxKernels]
+  [./periodic_dist]
+    type = PeriodicDistanceAux
+    variable = periodic_dist
+    point = '4 6 0'
+  [../]
+[]
+
+[BCs]
+  [./Periodic]
+    [./all]
+      variable = u
+      auto_direction = 'x y'
+    [../]
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  dt = 1
+  num_steps = 20
+  solve_type = NEWTON
+[]
+
+[Outputs]
+  execute_on = 'timestep_end'
+  exodus = true
+[]

test/tests/bcs/dmg_periodic/tests

@@ -0,0 +1,38 @@
+[Tests]
+  issues = '#15501'
+  design = 'AddPeriodicBCAction.md  DistributedRectilinearMeshGenerator.md'
+
+  [./2d]
+    type = 'Exodiff'
+    input = 'dmg_periodic_bc.i'
+    exodiff = 'dmg_periodic_bc_out.e'
+    min_parallel = 2
+    cli_args = 'Outputs/hide="pid"'
+    abs_zero = 1e-6
+    rel_err = 1e-5
+    requirement = "MOOSE shall support periodic boundary conditions in a parallel generated mesh."
+  [../]
+
+  [./3d]
+    type = 'Exodiff'
+    input = 'dmg_periodic_bc.i'
+    exodiff = 'dmg_periodic_bc_out_3d.e'
+    min_parallel = 2
+    cli_args = 'Outputs/hide="pid" Outputs/file_base=dmg_periodic_bc_out_3d Mesh/dmg/dim=3 Mesh/dmg/nx=10 Mesh/dmg/ny=10 Mesh/dmg/nz=10 Mesh/dmg/zmax=40 BCs/Periodic/all/auto_direction="x y z" Kernels/forcing/x_center=6'
+    abs_zero = 1e-6
+    rel_err = 1e-5
+    requirement = "MOOSE shall support periodic boundary conditions in a parallel generated 3D mesh."
+  [../]
+
+  [./1d]
+    type = 'Exodiff'
+    input = 'dmg_periodic_bc.i'
+    exodiff = 'dmg_periodic_bc_out_1d.e'
+    min_parallel = 2
+    cli_args = 'Outputs/hide="pid" Outputs/file_base=dmg_periodic_bc_out_1d Mesh/dmg/dim=1 Mesh/dmg/nx=100 Mesh/dmg/ny=0 Mesh/dmg/ymax=0 BCs/Periodic/all/auto_direction="x" AuxKernels/periodic_dist/point="1 0 0" Kernels/forcing/y_center=0  Kernels/forcing/x_center=0 '
+    abs_zero = 1e-6
+    rel_err = 1e-5
+    requirement = "MOOSE shall support periodic boundary conditions in a parallel generated 1D mesh."
+  [../]
+
+[]