Rebalance

.gitignore

@@ -0,0 +1,3 @@
+/user.mk
+*.o
+/laghos

laghos.cpp

@@ -125,6 +125,9 @@ int main(int argc, char *argv[])
    double blast_energy = 0.25;
    double blast_position[] = {0.0, 0.0, 0.0};
 
+   bool enable_nc = true;
+   bool enable_rebalance = true;
+
    OptionsParser args(argc, argv);
    args.AddOption(&dim, "-dim", "--dimension", "Dimension of the problem.");
    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
@@ -192,6 +195,13 @@ int main(int argc, char *argv[])
                   "Enable figure of merit output.");
    args.AddOption(&gpu_aware_mpi, "-gam", "--gpu-aware-mpi", "-no-gam",
                   "--no-gpu-aware-mpi", "Enable GPU aware MPI communications.");
+   args.AddOption(&enable_nc, "-nc", "--nonconforming", "-no-nc",
+                  "--conforming",
+                  "Use non-conforming meshes. Requires a 2D or 3D mesh.");
+   args.AddOption(&enable_rebalance, "-b", "--balance", "-no-b",
+                  "--no-rebalance",
+                  "Perform a rebalance after parallel refinement. Only enabled \n\t"
+                  "for non-conforming meshes with Metis partitioning.");
    args.AddOption(&dev, "-dev", "--dev", "GPU device to use.");
    args.Parse();
    if (!args.Good())
@@ -259,6 +269,15 @@ int main(int argc, char *argv[])
       }
    }
 
+   if (enable_nc && dim > 1)
+   {
+      if (Mpi::Root())
+      {
+         cout << "Using non-conforming mesh." << endl;
+      }
+      mesh->EnsureNCMesh();
+   }
+
    // Refine the mesh in serial to increase the resolution.
    for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
    const int mesh_NE = mesh->GetNE();
@@ -400,6 +419,12 @@ int main(int argc, char *argv[])
    // Refine the mesh further in parallel to increase the resolution.
    for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }
 
+   if (!cartesian_partitioning && enable_nc && dim > 1)
+   {
+      if (myid == 0) { cout << "Rebalancing mesh" << endl; }
+      pmesh->Rebalance();
+   }
+
    int NE = pmesh->GetNE(), ne_min, ne_max;
    MPI_Reduce(&NE, &ne_min, 1, MPI_INT, MPI_MIN, 0, pmesh->GetComm());
    MPI_Reduce(&NE, &ne_max, 1, MPI_INT, MPI_MAX, 0, pmesh->GetComm());
@@ -452,8 +477,8 @@ int main(int argc, char *argv[])
          return 3;
    }
 
-   const HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
-   const HYPRE_Int glob_size_h1 = H1FESpace.GlobalTrueVSize();
+   const HYPRE_BigInt glob_size_l2 = L2FESpace.GlobalTrueVSize();
+   const HYPRE_BigInt glob_size_h1 = H1FESpace.GlobalTrueVSize();
    if (Mpi::Root())
    {
       cout << "Number of kinematic (position, velocity) dofs: "

laghos_solver.cpp

@@ -21,6 +21,9 @@
 
 #ifdef MFEM_USE_MPI
 
+// for benchmark timing purposes; for a regular run this can be a no-op
+#define LAGHOS_DEVICE_SYNC MFEM_DEVICE_SYNC
+
 namespace mfem
 {
 
@@ -335,8 +338,10 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S,
 
    if (p_assembly)
    {
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Start();
       ForcePA->Mult(one, rhs);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Stop();
       rhs.Neg();
 
@@ -365,8 +370,10 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S,
          H1c.GetRestrictionMatrix()->Mult(dvc_gf, X);
          VMassPA->SetEssentialTrueDofs(c_tdofs[c]);
          VMassPA->EliminateRHS(B);
+         LAGHOS_DEVICE_SYNC;
          timer.sw_cgH1.Start();
          CG_VMass.Mult(B, X);
+         LAGHOS_DEVICE_SYNC;
          timer.sw_cgH1.Stop();
          timer.H1iter += CG_VMass.GetNumIterations();
          if (Pconf) { Pconf->Mult(X, dvc_gf); }
@@ -378,8 +385,10 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S,
    }
    else
    {
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Start();
       Force.Mult(one, rhs);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Stop();
       rhs.Neg();
 
@@ -402,8 +411,10 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S,
       cg.SetAbsTol(0.0);
       cg.SetMaxIter(cg_max_iter);
       cg.SetPrintLevel(-1);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_cgH1.Start();
       cg.Mult(B, X);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_cgH1.Stop();
       timer.H1iter += cg.GetNumIterations();
       Mv.RecoverFEMSolution(X, rhs, dv);
@@ -438,12 +449,16 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v,
    Array<int> l2dofs;
    if (p_assembly)
    {
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Start();
       ForcePA->MultTranspose(v, e_rhs);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Stop();
       if (e_source) { e_rhs += *e_source; }
+      LAGHOS_DEVICE_SYNC;
       timer.sw_cgL2.Start();
       CG_EMass.Mult(e_rhs, de);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_cgL2.Stop();
       const HYPRE_Int cg_num_iter = CG_EMass.GetNumIterations();
       timer.L2iter += (cg_num_iter==0) ? 1 : cg_num_iter;
@@ -453,17 +468,21 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v,
    }
    else // not p_assembly
    {
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Start();
       Force.MultTranspose(v, e_rhs);
+      LAGHOS_DEVICE_SYNC;
       timer.sw_force.Stop();
       if (e_source) { e_rhs += *e_source; }
       Vector loc_rhs(l2dofs_cnt), loc_de(l2dofs_cnt);
       for (int e = 0; e < NE; e++)
       {
          L2.GetElementDofs(e, l2dofs);
          e_rhs.GetSubVector(l2dofs, loc_rhs);
+         LAGHOS_DEVICE_SYNC;
          timer.sw_cgL2.Start();
          Me_inv(e).Mult(loc_rhs, loc_de);
+         LAGHOS_DEVICE_SYNC;
          timer.sw_cgL2.Stop();
          timer.L2iter += 1;
          de.SetSubVector(l2dofs, loc_de);
@@ -663,11 +682,11 @@ void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps,
    my_rt[4] = my_rt[0] + my_rt[2] + my_rt[3];
    MPI_Reduce(my_rt, T, 5, MPI_DOUBLE, MPI_MAX, 0, com);
 
-   HYPRE_Int mydata[3], alldata[3];
-   mydata[0] = timer.L2dof * timer.L2iter;
+   HYPRE_BigInt mydata[3], alldata[3];
+   mydata[0] = static_cast<HYPRE_BigInt>(timer.L2dof) * static_cast<HYPRE_BigInt>(timer.L2iter);
    mydata[1] = timer.quad_tstep;
    mydata[2] = NE;
-   MPI_Reduce(mydata, alldata, 3, HYPRE_MPI_INT, MPI_SUM, 0, com);
+   MPI_Reduce(mydata, alldata, 3, HYPRE_MPI_BIG_INT, MPI_SUM, 0, com);
 
    if (IamRoot)
    {
@@ -752,6 +771,7 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
    if (dim > 1 && p_assembly) { return qupdate->UpdateQuadratureData(S, qdata); }
 
    // This code is only for the 1D/FA mode
+   LAGHOS_DEVICE_SYNC;
    timer.sw_qdata.Start();
    const int nqp = ir.GetNPoints();
    ParGridFunction x, v, e;
@@ -914,6 +934,7 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
    delete [] p_b;
    delete [] cs_b;
    delete [] Jpr_b;
+   LAGHOS_DEVICE_SYNC;
    timer.sw_qdata.Stop();
    timer.quad_tstep += NE;
 }
@@ -1287,6 +1308,7 @@ void QKernel(const int NE, const int NQ,
 
 void QUpdate::UpdateQuadratureData(const Vector &S, QuadratureData &qdata)
 {
+   LAGHOS_DEVICE_SYNC;
    timer->sw_qdata.Start();
    Vector* S_p = const_cast<Vector*>(&S);
    const int H1_size = H1.GetVSize();
@@ -1335,6 +1357,7 @@ void QUpdate::UpdateQuadratureData(const Vector &S, QuadratureData &qdata)
                qdata.rho0DetJ0w, q_e, q_dv,
                qdata.Jac0inv, q_dt_est, qdata.stressJinvT);
    qdata.dt_est = q_dt_est.Min();
+   LAGHOS_DEVICE_SYNC;
    timer->sw_qdata.Stop();
    timer->quad_tstep += NE;
 }
@@ -1343,8 +1366,10 @@ void LagrangianHydroOperator::AssembleForceMatrix() const
 {
    if (forcemat_is_assembled || p_assembly) { return; }
    Force = 0.0;
+   LAGHOS_DEVICE_SYNC;
    timer.sw_force.Start();
    Force.Assemble();
+   LAGHOS_DEVICE_SYNC;
    timer.sw_force.Stop();
    forcemat_is_assembled = true;
 }

laghos_solver.hpp

@@ -103,10 +103,10 @@ class LagrangianHydroOperator : public TimeDependentOperator
    // FE spaces local and global sizes
    const int H1Vsize;
    const int H1TVSize;
-   const HYPRE_Int H1GTVSize;
+   const HYPRE_BigInt H1GTVSize;
    const int L2Vsize;
    const int L2TVSize;
-   const HYPRE_Int L2GTVSize;
+   const HYPRE_BigInt L2GTVSize;
    Array<int> block_offsets;
    // Reference to the current mesh configuration.
    mutable ParGridFunction x_gf;

makefile

@@ -55,6 +55,8 @@ make style
 
 endef
 
+-include user.mk
+
 NPROC = $(shell getconf _NPROCESSORS_ONLN)
 GOALS = help clean distclean style setup mfem metis hypre
 
@@ -76,8 +78,8 @@ TEST_MK = $(MFEM_TEST_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 ifeq (,$(filter $(GOALS),$(MAKECMDGOALS)))
    -include $(CONFIG_MK)
-   ifneq ($(realpath $(MFEM_DIR)),$(MFEM_SOURCE_DIR))
-      ifneq ($(realpath $(MFEM_DIR)),$(MFEM_INSTALL_DIR))
+   ifneq ($(realpath $(MFEM_DIR)),$(realpath $(MFEM_SOURCE_DIR)))
+      ifneq ($(realpath $(MFEM_DIR)),$(realpath $(MFEM_INSTALL_DIR)))
          MFEM_BUILD_DIR := $(MFEM_DIR)
          override MFEM_DIR := $(MFEM_SOURCE_DIR)
       endif