alternative preconditioner

examples/stokes.cpp

@@ -0,0 +1,246 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make stokes
+//
+// Sample runs:
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void rhs_func(const Vector & x, Vector &f)
+{
+   f.SetSize(x.Size());
+   f = 0.0;
+}
+
+void u0_func(const Vector & x, Vector &u0)
+{
+   u0.SetSize(x.Size());
+   u0 = 0.0;
+   // u0[0] = -sin(M_PI * x[1]);
+   u0[0] = -x[1]*(1-x[1])*x[2]*(1-x[2]);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/inline-hex.mesh";
+   int sref = 0;
+   int pref = 0;
+   bool visualization = true;
+
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&sref, "-sr", "--refinements-serial",
+                  "Number of serial refinements");
+   args.AddOption(&pref, "-pr", "--refinements-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Mesh mesh(mesh_file, 1, 1);
+
+   Mesh mesh = Mesh::MakeCartesian3D(8,1,1,mfem::Element::HEXAHEDRON,8.0,1.0,1.0);
+   int dim = mesh.Dimension();
+
+   for (int l = 0; l < sref; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   for (int l = 0; l < pref; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   FiniteElementCollection *ufec = new H1_FECollection(2,dim);
+   FiniteElementCollection *pfec = new H1_FECollection(1,dim);
+
+   ParFiniteElementSpace ufes(&pmesh, ufec,dim,mfem::Ordering::byVDIM);
+   ParFiniteElementSpace pfes(&pmesh, pfec);
+   Array<ParFiniteElementSpace *> pfespaces;
+   pfespaces.Append(&ufes);
+   pfespaces.Append(&pfes);
+   HYPRE_BigInt usize = ufes.GlobalTrueVSize();
+   HYPRE_BigInt psize = pfes.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of velocity dofs: " << usize << endl;
+      cout << "Number of pressure dofs: " << psize << endl;
+   }
+
+   ParGridFunction u_gf(&ufes); u_gf = 0.0;
+   ParGridFunction p_gf(&pfes); p_gf = 0.0;
+   VectorFunctionCoefficient u0_cf(dim, u0_func);
+   Array<int> ess_bdr_in(pmesh.bdr_attributes.Max()); ess_bdr_in = 0;
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max()); ess_bdr = 1;
+   ess_bdr_in[2] = 1; // x = 1 plane
+   ess_bdr[2] = 0; // x = 1 plane
+   ess_bdr[4] = 0; // x = 0 plane
+   Vector vzero(dim); vzero = 0.0;
+   VectorConstantCoefficient zero_cf(vzero);
+   u_gf.ProjectBdrCoefficient(zero_cf,ess_bdr);
+   u_gf.ProjectBdrCoefficient(u0_cf,ess_bdr_in);
+
+   // essential dofs;
+   Array<int> ess_tdof_list;
+   Array<int> ess_tdof_listp;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> essbdr(pmesh.bdr_attributes.Max());
+      essbdr = 1; essbdr[4] = 0;
+      ufes.GetEssentialTrueDofs(essbdr, ess_tdof_list);
+      essbdr = 0; essbdr[4] = 1;
+      pfes.GetEssentialTrueDofs(essbdr, ess_tdof_listp);
+   }
+
+   ParBilinearForm prec_pq(&pfes);
+   ConstantCoefficient one(1.0);
+   prec_pq.AddDomainIntegrator(new MassIntegrator(one));
+   prec_pq.Assemble();
+   HypreParMatrix Mpq;
+   prec_pq.FormSystemMatrix(ess_tdof_listp,Mpq);
+
+   for (int i = 0; i < ess_tdof_listp.Size(); i++)
+   {
+      ess_tdof_listp[i] += ufes.GetTrueVSize();
+   }
+   ess_tdof_list.Append(ess_tdof_listp);
+
+   // linear and bilinear forms
+   ParLinearForm b_v(&ufes);
+   VectorFunctionCoefficient rhs_cf(dim, rhs_func);
+   b_v.AddDomainIntegrator(new VectorDomainLFIntegrator(rhs_cf));
+   b_v.Assemble();
+
+   ParLinearForm b_q(&pfes);
+   b_q.Assemble();
+
+
+   ConstantCoefficient negone(-1.0);
+   // (∇ u, ∇ v)
+   ParBilinearForm a_uv(&ufes);
+   a_uv.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
+
+   ParMixedBilinearForm a_pv(&pfes,&ufes);
+   a_pv.AddDomainIntegrator(new TransposeIntegrator(new VectorDivergenceIntegrator(
+                                                       negone)));
+
+   ParMixedBilinearForm a_uq(&ufes,&pfes);
+   a_uq.AddDomainIntegrator(new VectorDivergenceIntegrator(negone));
+
+   ParBilinearForm a_pq(&pfes);
+   ConstantCoefficient zero(0.0);
+   a_pq.AddDomainIntegrator(new MassIntegrator(zero));
+
+   ParBlockForm bform(pfespaces);
+   bform.SetBlock(&a_uv,0,0);
+   bform.SetBlock(&a_pv,0,1);
+   bform.SetBlock(&a_uq,1,0);
+   bform.SetBlock(&a_pq,1,1);
+
+   Array<int> tdof_offsets(3);
+   tdof_offsets[0] = 0;
+   tdof_offsets[1] = ufes.TrueVSize();
+   tdof_offsets[2] = pfes.TrueVSize();
+   tdof_offsets.PartialSum();
+   Array<int> dof_offsets(3);
+   dof_offsets[0] = 0;
+   dof_offsets[1] = ufes.GetVSize();
+   dof_offsets[2] = pfes.GetVSize();
+   dof_offsets.PartialSum();
+
+   Vector b(dof_offsets.Last());
+   b.SetVector(b_v,0);
+   b.SetVector(b_q,dof_offsets[1]);
+   Vector x(dof_offsets.Last());
+   x.SetVector(u_gf,0);
+   x.SetVector(p_gf,dof_offsets[1]);
+
+   OperatorPtr Ah;
+   Vector B,X;
+   bform.Assemble();
+   bform.FormLinearSystem(ess_tdof_list,x, b, Ah, X, B);
+
+   BlockOperator * A = Ah.As<BlockOperator>();
+
+   BlockDiagonalPreconditioner prec(tdof_offsets);
+   HypreBoomerAMG amg_v((HypreParMatrix&)A->GetBlock(0,0));
+   amg_v.SetSystemsOptions(dim);
+
+   HypreBoomerAMG amg_p(Mpq);
+
+   prec.SetDiagonalBlock(0,&amg_v);
+   prec.SetDiagonalBlock(1,&amg_p);
+
+
+   // ((HypreParMatrix*)&A->GetBlock(0,0))->Print("A00.mat");
+   // // ((HypreParMatrix*)&A->GetBlock(0,1))->Print("A01.mat");
+   // ((HypreParMatrix*)&A->GetBlock(1,0))->Print("A10.mat");
+   // // ((HypreParMatrix*)&A->GetBlock(1,1))->Print("A11.mat");
+   // return 0;
+
+
+   MINRESSolver solver(MPI_COMM_WORLD);
+   solver.SetRelTol(1e-12);
+   solver.SetMaxIter(60);
+   solver.SetPrintLevel(1);
+   solver.SetPreconditioner(prec);
+   solver.SetOperator(*A);
+   solver.Mult(B, X);
+
+
+   BlockVector Xb(X.GetData(),tdof_offsets);
+   u_gf.SetFromTrueDofs(Xb.GetBlock(0));
+   p_gf.SetFromTrueDofs(Xb.GetBlock(1));
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << u_gf << flush;
+
+      socketstream solp_sock(vishost, visport);
+      solp_sock << "parallel " << num_procs << " " << myid << "\n";
+      solp_sock.precision(8);
+      solp_sock << "solution\n" << pmesh << p_gf << flush;
+
+   }
+
+   delete ufec;
+   delete pfec;
+
+   return 0;
+}

examples/stokes_ane.cpp

@@ -0,0 +1,260 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make stokes
+//
+// Sample runs:
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void rhs_func(const Vector & x, Vector &f)
+{
+   f.SetSize(x.Size());
+   f = 0.0;
+}
+
+void u0_func(const Vector & x, Vector &u0)
+{
+   u0.SetSize(x.Size());
+   u0 = 0.0;
+   double r0 = 2.393;
+   u0[0] = -(r0*r0 - x[1]*x[1] - x[2]*x[2]);
+   // u0[0] = +1;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "anu.msh";
+   int sref = 0;
+   int pref = 0;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&sref, "-sr", "--refinements-serial",
+                  "Number of serial refinements");
+   args.AddOption(&pref, "-pr", "--refinements-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   Mesh mesh(mesh_file, 1, 1);
+
+   int dim = mesh.Dimension();
+
+   for (int l = 0; l < sref; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   for (int l = 0; l < pref; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+
+   FiniteElementCollection *ufec = new H1_FECollection(2,dim);
+   FiniteElementCollection *pfec = new H1_FECollection(1,dim);
+
+   ParFiniteElementSpace ufes(&pmesh, ufec,dim,mfem::Ordering::byVDIM);
+   ParFiniteElementSpace pfes(&pmesh, pfec);
+   Array<ParFiniteElementSpace *> pfespaces;
+   pfespaces.Append(&ufes);
+   pfespaces.Append(&pfes);
+
+   HYPRE_BigInt usize = ufes.GlobalTrueVSize();
+   HYPRE_BigInt psize = pfes.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of velocity dofs: " << usize << endl;
+      cout << "Number of pressure dofs: " << psize << endl;
+   }
+
+   ParGridFunction u_gf(&ufes); u_gf = 0.0;
+   ParGridFunction p_gf(&pfes); p_gf = 0.0;
+   VectorFunctionCoefficient u0_cf(dim, u0_func);
+   Array<int> ess_bdr_in(pmesh.bdr_attributes.Max()); ess_bdr_in = 0;
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max()); ess_bdr = 1;
+   ess_bdr_in[1] = 1; // inflow
+   ess_bdr[0] = 0; // rest of bdr (except the outflow)
+   Vector vzero(dim); vzero = 0.0;
+   VectorConstantCoefficient zero_cf(vzero);
+   u_gf.ProjectBdrCoefficient(zero_cf,ess_bdr);
+   u_gf.ProjectBdrCoefficient(u0_cf,ess_bdr_in);
+
+   // essential dofs;
+   Array<int> ess_tdof_list;
+   Array<int> ess_tdof_listp;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> essbdr(pmesh.bdr_attributes.Max());
+      essbdr = 1; essbdr[2] = 0;
+      ufes.GetEssentialTrueDofs(essbdr, ess_tdof_list);
+      essbdr = 0; essbdr[2] = 1;
+      // pfes.GetEssentialTrueDofs(essbdr, ess_tdof_listp);
+   }
+   ConstantCoefficient one(1.0);
+
+   ParBilinearForm prec_pq(&pfes);
+   prec_pq.AddDomainIntegrator(new MassIntegrator(one));
+   prec_pq.Assemble();
+   HypreParMatrix Mpq;
+   prec_pq.FormSystemMatrix(ess_tdof_listp,Mpq);
+
+   for (int i = 0; i < ess_tdof_listp.Size(); i++)
+   {
+      ess_tdof_listp[i] += ufes.GetTrueVSize();
+   }
+   ess_tdof_list.Append(ess_tdof_listp);
+
+   // linear and bilinear forms
+   ParLinearForm b_v(&ufes);
+   VectorFunctionCoefficient rhs_cf(dim, rhs_func);
+   b_v.AddDomainIntegrator(new VectorDomainLFIntegrator(rhs_cf));
+   b_v.Assemble();
+
+   ParLinearForm b_q(&pfes);
+   b_q.Assemble();
+
+
+   ConstantCoefficient negone(-1.0);
+
+   // (∇ u, ∇ v)
+   ParBilinearForm a_uv(&ufes);
+   a_uv.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
+   ParMixedBilinearForm a_pv(&pfes,&ufes);
+   a_pv.AddDomainIntegrator(new TransposeIntegrator(new VectorDivergenceIntegrator(
+                                                       negone)));
+
+   ParMixedBilinearForm a_uq(&ufes,&pfes);
+   a_uq.AddDomainIntegrator(new VectorDivergenceIntegrator(negone));
+
+   ParBilinearForm a_pq(&pfes);
+
+   ParBlockForm bform(pfespaces);
+   bform.SetBlock(&a_uv,0,0);
+   bform.SetBlock(&a_pv,0,1);
+   bform.SetBlock(&a_uq,1,0);
+   bform.SetBlock(&a_pq,1,1);
+
+
+   Array<int> tdof_offsets(3);
+   tdof_offsets[0] = 0;
+   tdof_offsets[1] = ufes.TrueVSize();
+   tdof_offsets[2] = pfes.TrueVSize();
+   tdof_offsets.PartialSum();
+   Array<int> dof_offsets(3);
+   dof_offsets[0] = 0;
+   dof_offsets[1] = ufes.GetVSize();
+   dof_offsets[2] = pfes.GetVSize();
+   dof_offsets.PartialSum();
+
+   Vector b(dof_offsets.Last());
+   b.SetVector(b_v,0);
+   b.SetVector(b_q,dof_offsets[1]);
+   Vector x(dof_offsets.Last());
+   x.SetVector(u_gf,0);
+   x.SetVector(p_gf,dof_offsets[1]);
+
+   OperatorPtr Ah;
+   Vector B,X;
+   bform.Assemble();
+   bform.FormLinearSystem(ess_tdof_list,x, b, Ah, X, B);
+
+   BlockOperator * A = Ah.As<BlockOperator>();
+
+   BlockDiagonalPreconditioner prec(tdof_offsets);
+   HypreParMatrix A00 = (HypreParMatrix&)A->GetBlock(0,0);
+   HypreParMatrix A01 = (HypreParMatrix&)A->GetBlock(0,1);
+   HypreBoomerAMG amg_v(A00);
+   amg_v.SetSystemsOptions(dim);
+   prec.SetDiagonalBlock(0,&amg_v);
+
+   HypreParVector A00_diag(MPI_COMM_WORLD, A00.GetGlobalNumRows(),
+                                 A00.GetRowStarts());
+   A00.GetDiag(A00_diag);
+   HypreParMatrix S_tmp(A01);
+   S_tmp.InvScaleRows(A00_diag);
+   HypreParMatrix *S = ParMult(A01.Transpose(), &S_tmp, true);
+   HypreBoomerAMG amg_p(*S);
+   // HypreBoomerAMG amg_p(Mpq);
+   prec.SetDiagonalBlock(1,&amg_p);
+
+   MINRESSolver solver(MPI_COMM_WORLD);
+   solver.SetRelTol(1e-12);
+   solver.SetMaxIter(20000);
+   solver.SetPrintLevel(1);
+   solver.SetPreconditioner(prec);
+   solver.SetOperator(*A);
+   solver.Mult(B, X);
+
+   BlockVector Xb(X.GetData(),tdof_offsets);
+   u_gf.SetFromTrueDofs(Xb.GetBlock(0));
+   p_gf.SetFromTrueDofs(Xb.GetBlock(1));
+
+   ParaViewDataCollection * paraview_dc = nullptr;
+   bool paraview = true;
+   if (paraview)
+   {
+      paraview_dc = new ParaViewDataCollection(mesh_file, &pmesh);
+      paraview_dc->SetPrefixPath("ParaView");
+      paraview_dc->SetLevelsOfDetail(2);
+      paraview_dc->SetCycle(0);
+      paraview_dc->SetDataFormat(VTKFormat::BINARY);
+      paraview_dc->SetHighOrderOutput(true);
+      paraview_dc->SetTime(0.0); // set the time
+      paraview_dc->RegisterField("velocity",&u_gf);
+      paraview_dc->RegisterField("pressure",&p_gf);
+      paraview_dc->Save();
+   }
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << u_gf << flush;
+
+      socketstream solp_sock(vishost, visport);
+      solp_sock << "parallel " << num_procs << " " << myid << "\n";
+      solp_sock.precision(8);
+      solp_sock << "solution\n" << pmesh << p_gf << flush;
+
+   }
+
+   delete ufec;
+   delete pfec;
+
+   return 0;
+}

fem/CMakeLists.txt

@@ -268,7 +268,8 @@ if (MFEM_USE_MPI)
     pgridfunc.cpp
     plinearform.cpp
     pnonlinearform.cpp
-    prestriction.cpp)
+    prestriction.cpp
+    pblockform.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the
   # headers added all the time.
   list(APPEND HDRS
@@ -277,7 +278,8 @@ if (MFEM_USE_MPI)
     pgridfunc.hpp
     plinearform.hpp
     pnonlinearform.hpp
-    prestriction.hpp)
+    prestriction.hpp
+    pblockform.hpp)
 endif()
 
 convert_filenames_to_full_paths(SRCS)

fem/fem.hpp

@@ -55,6 +55,7 @@
 #include "plinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pnonlinearform.hpp"
+#include "pblockform.hpp"
 #endif
 
 #ifdef MFEM_USE_SIDRE

fem/pblockform.cpp

@@ -0,0 +1,241 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// Implementation of class ParBlockForm
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+void ParBlockForm::FillEssTdofLists(const Array<int> & ess_tdof_list)
+{
+   int j;
+   for (int i = 0; i<ess_tdof_list.Size(); i++)
+   {
+      int tdof = ess_tdof_list[i];
+      for (j = 0; j < nblocks; j++)
+      {
+         if (tdof_offsets[j+1] > tdof) { break; }
+      }
+      ess_tdofs[j]->Append(tdof-tdof_offsets[j]);
+   }
+}
+
+void ParBlockForm::BuildProlongation()
+{
+   P = new BlockOperator(dof_offsets, tdof_offsets);
+   R = new BlockMatrix(tdof_offsets, dof_offsets);
+   P->owns_blocks = 0;
+   R->owns_blocks = 0;
+
+   for (int i = 0; i<nblocks; i++)
+   {
+      HypreParMatrix * P_ = pfes[i]->Dof_TrueDof_Matrix();
+      P->SetBlock(i,i,P_);
+      const SparseMatrix * R_ = pfes[i]->GetRestrictionMatrix();
+      R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
+   }
+}
+
+void ParBlockForm::ParallelAssemble(BlockMatrix *m)
+{
+   if (!P) { BuildProlongation(); }
+
+   p_mat = new BlockOperator(tdof_offsets);
+   p_mat_e = new BlockOperator(tdof_offsets);
+   p_mat->owns_blocks = 1;
+   p_mat_e->owns_blocks = 1;
+   HypreParMatrix * A = nullptr;
+   HypreParMatrix * PtAP = nullptr;
+   for (int i = 0; i<nblocks; i++)
+   {
+      HypreParMatrix * Pi = (HypreParMatrix*)(&P->GetBlock(i,i));
+      HypreParMatrix * Pit = Pi->Transpose();
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (m->IsZeroBlock(i,j)) { continue; }
+         if (i == j)
+         {
+            // Make block diagonal square hypre matrix
+            A = new HypreParMatrix(pfes[i]->GetComm(), pfes[i]->GlobalVSize(),
+                                   pfes[i]->GetDofOffsets(),&m->GetBlock(i,i));
+            PtAP = RAP(A,Pi);
+            delete A;
+            p_mat_e->SetBlock(i,i,PtAP->EliminateRowsCols(*ess_tdofs[i]));
+         }
+         else
+         {
+            HypreParMatrix * Pj = (HypreParMatrix*)(&P->GetBlock(j,j));
+            A = new HypreParMatrix(pfes[i]->GetComm(), pfes[i]->GlobalVSize(),
+                                   pfes[j]->GlobalVSize(), pfes[i]->GetDofOffsets(),
+                                   pfes[j]->GetDofOffsets(), &m->GetBlock(i,j));
+            HypreParMatrix * APj = ParMult(A, Pj,true);
+            delete A;
+            PtAP = ParMult(Pit,APj,true);
+            delete APj;
+            p_mat_e->SetBlock(i,j,PtAP->EliminateCols(*ess_tdofs[j]));
+            PtAP->EliminateRows(*ess_tdofs[i]);
+         }
+         p_mat->SetBlock(i,j,PtAP);
+      }
+      delete Pit;
+   }
+}
+
+
+ParBlockForm::ParBlockForm(const Array<ParFiniteElementSpace*> pfes_ ): pfes(
+      pfes_)
+{
+   nblocks = pfes.Size();
+   bforms.SetSize(nblocks,nblocks);
+   mforms.SetSize(nblocks,nblocks);
+   ess_tdofs.SetSize(nblocks);
+   dof_offsets.Append(0);
+   tdof_offsets.Append(0);
+   for (int i = 0; i<nblocks; i++)
+   {
+      dof_offsets.Append(pfes[i]->GetVSize());
+      tdof_offsets.Append(pfes[i]->TrueVSize());
+      ess_tdofs[i] = new Array<int>();
+      for (int j = 0; j<nblocks; j++)
+      {
+         bforms(i,j) = nullptr;
+         mforms(i,j) = nullptr;
+      }
+   }
+   dof_offsets.PartialSum();
+   tdof_offsets.PartialSum();
+}
+
+void ParBlockForm::SetBlock(ParBilinearForm * bform, int row_idx, int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   bforms(row_idx,col_idx) = bform;
+}
+void ParBlockForm::SetBlock(ParMixedBilinearForm * mform, int row_idx,
+                            int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   mforms(row_idx,col_idx) = mform;
+}
+
+/// Assemble the local matrix
+void ParBlockForm::Assemble(int skip_zeros)
+{
+   sp_mat = new BlockMatrix(dof_offsets);
+   for (int i = 0; i<nblocks; i++)
+   {
+      int h = dof_offsets[i+1]-dof_offsets[i];
+      for (int j = 0; j<nblocks; j++)
+      {
+         int w = dof_offsets[j+1]-dof_offsets[j];
+         if (bforms(i,j))
+         {
+            bforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = bforms(i,j)->Height(), "inconsistent height of bilinear form");
+            MFEM_VERIFY(w = bforms(i,j)->Width(), "inconsistent width of bilinear form");
+            sp_mat->SetBlock(i,j,&bforms(i,j)->SpMat());
+         }
+         else if (mforms(i,j))
+         {
+            mforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = mforms(i,j)->Height(),
+                        "inconsistent height of MixedBilinear form");
+            MFEM_VERIFY(w = mforms(i,j)->Width(),
+                        "inconsistent width of Mixedbilinear form");
+            sp_mat->SetBlock(i,j,&mforms(i,j)->SpMat());
+         }
+         else
+         {
+            sp_mat->SetBlock(i,j,nullptr);
+         }
+      }
+   }
+}
+
+void ParBlockForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
+                                    Vector & b,
+                                    OperatorHandle &A, Vector &X,
+                                    Vector &B, int copy_interior)
+{
+   FormSystemMatrix(ess_tdof_list, A);
+   B.SetSize(P->Width());
+   P->MultTranspose(b,B);
+   X.SetSize(R->Height());
+   R->Mult(x,X);
+
+   // eliminate tdof in RHS
+   // B -= Ae*X
+   Vector tmp(B.Size());
+   p_mat_e->Mult(X,tmp);
+   B-=tmp;
+
+   for (int j = 0; j<nblocks; j++)
+   {
+      if (!ess_tdofs[j]->Size()) { continue; }
+      for (int i = 0; i < ess_tdofs[j]->Size(); i++)
+      {
+         int tdof = (*ess_tdofs[j])[i];
+         int gdof = tdof + tdof_offsets[j];
+         B(gdof) = X(gdof); // diagonal policy in always one in parallel
+      }
+   }
+   if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+}
+
+void ParBlockForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
+                                    OperatorHandle &A)
+{
+   FillEssTdofLists(ess_tdof_list);
+   if (sp_mat)
+   {
+      sp_mat->Finalize();
+      ParallelAssemble(sp_mat);
+      delete sp_mat;
+      sp_mat = nullptr;
+   }
+   A.Reset(p_mat,false);
+}
+
+void ParBlockForm::RecoverFEMSolution(const Vector &X, Vector &x)
+{
+   x.SetSize(P->Height());
+   P->Mult(X, x);
+}
+
+ParBlockForm::~ParBlockForm()
+{
+   delete p_mat_e;
+   p_mat_e = nullptr;
+   delete p_mat;
+   p_mat = nullptr;
+   for (int i = 0; i<nblocks; i++)
+   {
+      delete ess_tdofs[i];
+   }
+   delete P;
+   delete R;
+}
+
+};
+
+#endif

fem/pblockform.hpp

@@ -0,0 +1,82 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PBLOCKFORM
+#define MFEM_PBLOCKFORM
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include <mpi.h>
+#include "pfespace.hpp"
+#include "pgridfunc.hpp"
+#include "pbilinearform.hpp"
+
+namespace mfem
+{
+// square block forms
+class ParBlockForm
+{
+
+private:
+   // BilinearForms
+   Array2D<ParBilinearForm * > bforms;
+   Array2D<ParMixedBilinearForm * > mforms;
+   Array<ParFiniteElementSpace *> pfes;
+   Array<int> dof_offsets;
+   Array<int> tdof_offsets;
+   int nblocks;
+
+   // ess_tdof list for each space
+   Array<Array<int> *> ess_tdofs;
+
+   // Block operator of HypreParMatrix
+   BlockOperator * P = nullptr; // Block Prolongation
+   BlockMatrix * R = nullptr; // Block Restriction
+
+   BlockMatrix * sp_mat = nullptr;
+   // Block operator of HypreParMatrix
+   BlockOperator * p_mat = nullptr;
+   BlockOperator * p_mat_e = nullptr;
+
+   void FillEssTdofLists(const Array<int> & ess_tdof_list);
+   void BuildProlongation();
+   void ParallelAssemble(BlockMatrix *m);
+
+public:
+
+   ParBlockForm(const Array<ParFiniteElementSpace*> pfes_ );
+
+   void SetBlock(ParBilinearForm * bform, int row_idx, int col_idx);
+   void SetBlock(ParMixedBilinearForm * mform, int row_idx, int col_idx);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector & b,
+                         OperatorHandle &A, Vector &X,
+                         Vector &B, int copy_interior = 0);
+
+   void FormSystemMatrix(const Array<int> &ess_tdof_list,
+                         OperatorHandle &A);
+
+   void RecoverFEMSolution(const Vector &X, Vector &x);
+
+   /// Destroys bilinear form.
+   ~ParBlockForm();
+};
+
+}
+
+#endif // MFEM_USE_MPI
+
+#endif