How to create large parallel Mesh

[ilsunshine]

Enviroment: mfem version:4.6; metis version: 5.1.0; hypre version: 2.22.0
Problem: I'm stuck on how to construct a large-scale mesh within a fixed memory limit (without considering the time it takes to build the mesh).

Problem details：When building parallel algorithms using MFEM, I used ParMesh to construct parallel meshes. Among the seven constructors provided in the ParMesh documentation,
one is a default constructor, four require an already constructed ParMesh as input, one requires reading an already created parallel mesh file, and only one function can be used without relying on an existing parallel mesh:
mfem::ParMesh::ParMesh ( MPI_Comm comm,Mesh & mesh,int * partitioning_ = NULL,int part_method = 1). However, this approach requires creating a mesh on each process, which can lead to memory explosion when the mesh size and the number of processes increase. I couldn't find any constructors in the mesh library that return a mesh object (e.g., Mesh::MakeCartesian2D). Therefore, I'm stuck on how to construct a large-scale mesh given a fixed amount of memory (without considering the time it takes to build the mesh).

The help I need：Tell me if there are other ways to relax the constraints for building a larger-scale ParMesh (with a rectangular or cuboid region). Alternatively, can the following code be optimized in terms of memory usage?

Below is the minimal reproducible code for my problem. As the number of processes increases, memory consumption increases rapidly. At 24 processes, it exceeds 192GB of memory, causing errors in my operating environment.

#include <mfem.hpp>
#include <mpi.h>
int main(int argc, char *argv[]) {
MPI_Init(&argc, &argv); 

int Dim = 2;  // Or 3 for 3D case
int nx = 4000, ny = 4000;  // Large values triggering the error
double sx = 1.0, sy = 1.0;
int nz = 4000;  // For 3D
double sz = 1.0;
int rank;

MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
mfem::Mesh *mesh=nullptr;
int NE;
if (Dim == 2) {
mesh = new mfem::Mesh(
    mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, sx, sy));
} else if (Dim == 3) {
    mesh = new mfem::Mesh(
        mfem::Mesh::MakeCartesian3D(nx, ny, nz, mfem::Element::TETRAHEDRON, sx, sy, sz, false));
}

NE=mesh->GetNE();
int *par_partition = nullptr;
if (rank == 0) {
    par_partition = mesh->GeneratePartitioning(mfem::Mpi::WorldSize());
}
else
{
    par_partition = new int[NE];
}
MPI_Bcast(par_partition, NE, MPI_INT, 0, MPI_COMM_WORLD);
ParMesh *parmesh = new ParMesh(comm, *mesh, par_partition);
delete mesh;
delete parmesh;
delete[] par_partition;

MPI_Finalize(); 
return 0;
}

[ilsunshine]

This record documents a feasible method obtained from issue #5185. We can use mfem::MeshPartitioner and mfem::MeshPart which requires MFEM version 4.7 or later to extract a portion of the mesh and save that. Below is the example code provided by #5185.

#include <mfem.hpp>

#include <fstream>

int main(int argc, char *argv[])
{
   int Dim = 2;              // Or 3 for 3D case
   int nx = 4000, ny = 4000;
   double sx = 1.0, sy = 1.0;
   int nz = 4000; // For 3D
   double sz = 1.0;
   int nparts = 1; // must match how many ranks you intend to read in from

   mfem::OptionsParser args(argc, argv);
   args.AddOption(&nparts, "-n", "--nparts", "Number of ranks to partition to.");
   args.ParseCheck();

   
   mfem::Mesh *mesh;
   if (Dim == 2)
   {
      mesh = new mfem::Mesh(mfem::Mesh::MakeCartesian2D(
         nx, ny, mfem::Element::TRIANGLE, true, sx, sy));
   }
   else if (Dim == 3)
   {
      mesh = new mfem::Mesh(mfem::Mesh::MakeCartesian3D(
         nx, ny, nz, mfem::Element::TETRAHEDRON, sx, sy, sz, false));
   }

   int *par_partition = mesh->GeneratePartitioning(nparts);

   mfem::MeshPartitioner partitioner(*mesh, nparts, par_partition);
   mfem::MeshPart mesh_part;
   for (int i = 0; i < nparts; i++)
   {
      partitioner.ExtractPart(i, mesh_part);
      std::ofstream omesh(mfem::MakeParFilename("my-mesh.", i));
      mesh_part.Print(omesh);
   }

   // Cleanup
   delete mesh;
   delete[] par_partition;
   return 0;
}