RZ example with Robin conditions

Temporarily this means disabling partial assembly for the RT mass term.

The PA functionality is implemented in PR #3539.

.gitignore

@@ -309,6 +309,9 @@ miniapps/solvers/sol.*
 
 miniapps/hdiv-linear-solver/darcy
 miniapps/hdiv-linear-solver/grad_div
+miniapps/hdiv-linear-solver/residual
+miniapps/hdiv-linear-solver/rz
+miniapps/hdiv-linear-solver/ParaView
 
 miniapps/parelag/MultilevelHcurlHdivSolver
 miniapps/parelag/*.mesh

fem/bilininteg_mass_pa.cpp

@@ -56,8 +56,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dim = mesh->Dimension();
    ne = fes.GetMesh()->GetNE();
    nq = ir->GetNPoints();
-   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
-                                    GeometricFactors::JACOBIANS, mt);
+   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS, mt);
    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
@@ -74,7 +73,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
       const bool const_c = coeff.Size() == 1;
       const bool by_val = map_type == FiniteElement::VALUE;
       const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D);
-      const auto J = Reshape(geom->J.Read(), Q1D,Q1D,2,2,NE);
+      const auto J = Reshape(geom->detJ.Read(), Q1D,Q1D,NE);
       const auto C = const_c ? Reshape(coeff.Read(), 1,1,1) :
                      Reshape(coeff.Read(), Q1D,Q1D,NE);
       auto v = Reshape(pa_data.Write(), Q1D,Q1D, NE);
@@ -84,11 +83,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
          {
             MFEM_FOREACH_THREAD(qy,y,Q1D)
             {
-               const double J11 = J(qx,qy,0,0,e);
-               const double J12 = J(qx,qy,1,0,e);
-               const double J21 = J(qx,qy,0,1,e);
-               const double J22 = J(qx,qy,1,1,e);
-               const double detJ = (J11*J22)-(J21*J12);
+               const double detJ = J(qx,qy,e);
                const double coeff = const_c ? C(0,0,0) : C(qx,qy,e);
                v(qx,qy,e) =  W(qx,qy) * coeff * (by_val ? detJ : 1.0/detJ);
             }
@@ -102,7 +97,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
       const bool const_c = coeff.Size() == 1;
       const bool by_val = map_type == FiniteElement::VALUE;
       const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D,Q1D);
-      const auto J = Reshape(geom->J.Read(), Q1D,Q1D,Q1D,3,3,NE);
+      const auto J = Reshape(geom->detJ.Read(), Q1D,Q1D,Q1D,NE);
       const auto C = const_c ? Reshape(coeff.Read(), 1,1,1,1) :
                      Reshape(coeff.Read(), Q1D,Q1D,Q1D,NE);
       auto v = Reshape(pa_data.Write(), Q1D,Q1D,Q1D,NE);
@@ -114,18 +109,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
             {
                MFEM_FOREACH_THREAD(qz,z,Q1D)
                {
-                  const double J11 = J(qx,qy,qz,0,0,e);
-                  const double J21 = J(qx,qy,qz,1,0,e);
-                  const double J31 = J(qx,qy,qz,2,0,e);
-                  const double J12 = J(qx,qy,qz,0,1,e);
-                  const double J22 = J(qx,qy,qz,1,1,e);
-                  const double J32 = J(qx,qy,qz,2,1,e);
-                  const double J13 = J(qx,qy,qz,0,2,e);
-                  const double J23 = J(qx,qy,qz,1,2,e);
-                  const double J33 = J(qx,qy,qz,2,2,e);
-                  const double detJ = J11 * (J22 * J33 - J32 * J23) -
-                  /* */               J21 * (J12 * J33 - J32 * J13) +
-                  /* */               J31 * (J12 * J23 - J22 * J13);
+                  const double detJ = J(qx,qy,qz,e);
                   const double coeff = const_c ? C(0,0,0,0) : C(qx,qy,qz,e);
                   v(qx,qy,qz,e) = W(qx,qy,qz) * coeff * (by_val ? detJ : 1.0/detJ);
                }

miniapps/hdiv-linear-solver/change_basis.cpp

@@ -63,7 +63,7 @@ void Transpose(const Array<double> &B, Array<double> &Bt)
 {
    const int n = sqrt(B.Size());
    Bt.SetSize(n*n);
-   for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) Bt[i+j*n] = B[j+i*n];
+   for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
 }
 
 ChangeOfBasis_L2::ChangeOfBasis_L2(FiniteElementSpace &fes)
@@ -135,7 +135,8 @@ ChangeOfBasis_RT::ChangeOfBasis_RT(FiniteElementSpace &fes)
    const int cb_type = rt_fec->GetClosedBasisType();
    const int ob_type = rt_fec->GetOpenBasisType();
 
-   no_op = (cb_type == BasisType::GaussLobatto && ob_type == BasisType::IntegratedGLL);
+   no_op = (cb_type == BasisType::GaussLobatto &&
+            ob_type == BasisType::IntegratedGLL);
    if (no_op) { return; }
 
    const int pp1 = p + 1;
@@ -313,17 +314,17 @@ void ChangeOfBasis_RT::Mult(const Vector &x, Vector &y, Mode mode) const
 
    const Operator *P = fes.GetProlongationMatrix();
 
-   if (P)
+   if (IsIdentityProlongation(P))
+   {
+      x_l.MakeRef(const_cast<Vector&>(x), 0, fes.GetVSize());
+      y_l.MakeRef(y, 0, fes.GetVSize());
+   }
+   else
    {
       x_l.SetSize(fes.GetVSize());
       y_l.SetSize(fes.GetVSize());
       P->Mult(x, x_l);
    }
-   else
-   {
-      x_l.MakeRef(const_cast<Vector&>(x), 0);
-      y_l.MakeRef(y, 0);
-   }
 
    x_e.SetSize(elem_restr->Height());
    y_e.SetSize(elem_restr->Height());

miniapps/hdiv-linear-solver/darcy.cpp

@@ -53,6 +53,36 @@ using namespace mfem;
 
 ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
 
+class RobinCoefficient : public Coefficient
+{
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      double xdata[3];
+      Vector xvec(xdata, 3);
+      T.Transform(ip, xvec);
+      const int dim = xvec.Size();
+
+      Vector n(dim);
+      CalcOrtho(T.Jacobian(), n);
+      n /= n.Norml2();
+
+      const double p_val = u(xvec);
+      const double x = pi*xvec[0];
+      const double y = pi*xvec[1];
+
+      if (dim == 2)
+      {
+         const double u_val = -pi*(n[0]*cos(x)*sin(y) + n[1]*sin(x)*cos(y));
+         return p_val - u_val;
+      }
+      else
+      {
+         MFEM_ABORT("Not implemented");
+      }
+      return 0.0;
+   }
+};
+
 int main(int argc, char *argv[])
 {
    Mpi::Init(argc, argv);
@@ -85,7 +115,7 @@ int main(int argc, char *argv[])
    MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
 
    const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
-   const int mt = FiniteElement::INTEGRAL;
+   const int mt = FiniteElement::VALUE;
    RT_FECollection fec_rt(order-1, dim, b1, b2);
    L2_FECollection fec_l2(order-1, dim, b2, mt);
    ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
@@ -103,6 +133,8 @@ int main(int argc, char *argv[])
 
    // f is the RHS, u is the exact solution
    FunctionCoefficient f_coeff(f(alpha)), u_coeff(u);
+   // Coefficient to enforce Robin boundary condition
+   RobinCoefficient bc_coeff;
 
    // Assemble the right-hand side for the scalar (L2) unknown.
    ParLinearForm b_l2(&fes_l2);
@@ -113,7 +145,7 @@ int main(int argc, char *argv[])
    // Enforce Dirichlet boundary conditions on the scalar unknown by adding
    // the boundary term to the flux equation.
    ParLinearForm b_rt(&fes_rt);
-   b_rt.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(u_coeff));
+   b_rt.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(bc_coeff));
    b_rt.UseFastAssembly(true);
    b_rt.Assemble();
 
@@ -141,7 +173,7 @@ int main(int argc, char *argv[])
    ConstantCoefficient alpha_coeff(alpha);
    const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
    HdivSaddlePointSolver saddle_point_solver(
-      mesh, fes_rt, fes_l2, alpha_coeff, one, ess_rt_dofs, solver_mode);
+      mesh, fes_rt, fes_l2, alpha_coeff, one, one, ess_rt_dofs, solver_mode);
 
    const Array<int> &offsets = saddle_point_solver.GetOffsets();
    BlockVector X_block(offsets), B_block(offsets);
@@ -157,8 +189,8 @@ int main(int argc, char *argv[])
    if (Mpi::Root())
    {
       cout << "Done.\nIterations: "
-            << saddle_point_solver.GetNumIterations()
-            << "\nElapsed: " << tic_toc.RealTime() << endl;
+           << saddle_point_solver.GetNumIterations()
+           << "\nElapsed: " << tic_toc.RealTime() << endl;
    }
 
    ParGridFunction x(&fes_l2);

miniapps/hdiv-linear-solver/hdiv_linear_solver.cpp

@@ -25,7 +25,8 @@ void Reciprocal(Vector &x)
 }
 
 /// Return a new HypreParMatrix with given diagonal entries
-HypreParMatrix *MakeDiagonalMatrix(Vector &diag, const ParFiniteElementSpace &fes)
+HypreParMatrix *MakeDiagonalMatrix(Vector &diag,
+                                   const ParFiniteElementSpace &fes)
 {
    const int n = diag.Size();
 
@@ -63,8 +64,8 @@ const IntegrationRule &GetMassIntRule(FiniteElementSpace &fes_l2)
 
 HdivSaddlePointSolver::HdivSaddlePointSolver(
    ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
-   Coefficient &L_coeff_, Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_,
-   Mode mode_)
+   Coefficient &L_coeff_, Coefficient &R_coeff_, Coefficient &B_coeff_,
+   const Array<int> &ess_rt_dofs_, Mode mode_)
    : minres(mesh.GetComm()),
      order(fes_rt_.GetMaxElementOrder()),
      fec_l2(order - 1, mesh.Dimension(), b2, mt),
@@ -74,14 +75,18 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
      ess_rt_dofs(ess_rt_dofs_),
      basis_l2(fes_l2_),
      basis_rt(fes_rt_),
+     convert_map_type(fes_l2_.GetFE(0)->GetMapType() == FiniteElement::VALUE),
      mass_l2(&fes_l2),
      mass_rt(&fes_rt),
      L_coeff(L_coeff_),
      R_coeff(R_coeff_),
+     B_coeff(B_coeff_),
      mode(mode_),
      qs(mesh, GetMassIntRule(fes_l2)),
-     qf(qs),
-     l2_qf_coeff(qf)
+     W_coeff_qf(qs),
+     W_mix_coeff_qf(qs),
+     W_coeff(W_coeff_qf),
+     W_mix_coeff(W_mix_coeff_qf)
 {
    // If the user gives zero L coefficient, switch mode to DARCY_ZERO
    auto *L_const_coeff = dynamic_cast<ConstantCoefficient*>(&L_coeff);
@@ -93,11 +98,12 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
                   "Mode::GRAD_DIV incompatible with zero coefficient.");
    }
 
-   mass_l2.AddDomainIntegrator(new MassIntegrator(l2_qf_coeff));
+   mass_l2.AddDomainIntegrator(new MassIntegrator(W_coeff));
    mass_l2.SetAssemblyLevel(AssemblyLevel::PARTIAL);
 
    mass_rt.AddDomainIntegrator(new VectorFEMassIntegrator(&R_coeff));
-   mass_rt.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   mass_rt.AddBoundaryIntegrator(new MassIntegrator(B_coeff));
+   // mass_rt.SetAssemblyLevel(AssemblyLevel::PARTIAL);
 
    D.reset(FormDiscreteDivergenceMatrix(fes_rt, fes_l2, ess_rt_dofs));
    Dt.reset(D->Transpose());
@@ -125,7 +131,19 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
    if (mode == Mode::DARCY && !zero_l2_block)
    {
       ParBilinearForm mass_l2_unweighted(&fes_l2);
-      mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator);
+      QuadratureFunction det_J_qf(qs);
+      QuadratureFunctionCoefficient det_J_coeff(det_J_qf);
+      if (convert_map_type)
+      {
+         const auto flags = GeometricFactors::DETERMINANTS;
+         auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
+         det_J_qf = geom->detJ;
+         mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator(det_J_coeff));
+      }
+      else
+      {
+         mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator);
+      }
       mass_l2_unweighted.SetAssemblyLevel(AssemblyLevel::PARTIAL);
       mass_l2_unweighted.Assemble();
       const int n_l2 = fes_l2.GetTrueVSize();
@@ -136,45 +154,83 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
    Setup();
 }
 
+HdivSaddlePointSolver::HdivSaddlePointSolver(
+   ParMesh &mesh_, ParFiniteElementSpace &fes_rt_,
+   ParFiniteElementSpace &fes_l2_, Coefficient &L_coeff_, Coefficient &R_coeff_,
+   const Array<int> &ess_rt_dofs_, Mode mode_)
+   : HdivSaddlePointSolver(mesh_, fes_rt_, fes_l2_, L_coeff_, R_coeff_, zero,
+                           ess_rt_dofs_, mode_)
+{ }
+
 HdivSaddlePointSolver::HdivSaddlePointSolver(
    ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
    Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_)
-   : HdivSaddlePointSolver(mesh, fes_rt_, fes_l2_, zero, R_coeff_, ess_rt_dofs_, Mode::DARCY)
+   : HdivSaddlePointSolver(mesh, fes_rt_, fes_l2_, zero, R_coeff_, zero,
+                           ess_rt_dofs_, Mode::DARCY)
 { }
 
 void HdivSaddlePointSolver::Setup()
 {
-   if (!zero_l2_block) { L_coeff.Project(qf); }
+   const auto flags = GeometricFactors::DETERMINANTS;
+   auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
 
-   // Reassemble the L2 mass diagonal with the new coefficient
-   mass_l2.Assemble();
-   mass_l2.AssembleDiagonal(L_diag);
-   mass_l2.FormSystemMatrix(empty, L);
+   if (!zero_l2_block) { L_coeff.Project(W_coeff_qf); }
+   // In "grad-div mode", the transformation matrix is scaled by the coefficient
+   // of the mass and divergence matrices.
+   // In "Darcy mode", the transformation matrix is unweighted.
+   if (mode == Mode::GRAD_DIV) { W_mix_coeff_qf = W_coeff_qf; }
+   else { W_mix_coeff_qf = 1.0; }
 
-   if (mode == Mode::GRAD_DIV)
+   // The transformation matrix has to be "mixed" value and integral map type,
+   // which means that the coefficient has to be scaled like the Jacobian
+   // determinant.
+   if (convert_map_type)
    {
-      L_inv.reset(new DGMassInverse(fes_l2, l2_qf_coeff));
-      A_11 = L_inv;
-      Reciprocal(L_diag);
+      const int n = W_mix_coeff_qf.Size();
+      const double *d_detJ = geom->detJ.Read();
+      double *d_w_mix = W_mix_coeff_qf.ReadWrite();
+      double *d_w = W_coeff_qf.ReadWrite();
+      const bool zero_l2 = zero_l2_block;
+      MFEM_FORALL(i, n,
+      {
+         const double detJ = d_detJ[i];
+         if (!zero_l2) { d_w[i] *= detJ*detJ; }
+         d_w_mix[i] *= detJ;
+      });
+   }
+
+   L_inv.reset(new DGMassInverse(fes_l2, W_mix_coeff));
+
+   if (zero_l2_block)
+   {
+      A_11.reset();
    }
    else
    {
-      L_inv.reset(new DGMassInverse(fes_l2));
-      if (zero_l2_block)
+      mass_l2.Assemble();
+      mass_l2.AssembleDiagonal(L_diag);
+      mass_l2.FormSystemMatrix(empty, L);
+
+      A_11.reset(new RAPOperator(*L_inv, *L, *L_inv));
+
+      if (mode == GRAD_DIV)
       {
-         A_11.reset();
+         L_diag_unweighted.SetSize(L_diag.Size());
+
+         BilinearForm mass_l2_mix(&fes_l2);
+         mass_l2_mix.AddDomainIntegrator(new MassIntegrator(W_mix_coeff));
+         mass_l2_mix.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+         mass_l2_mix.Assemble();
+         mass_l2_mix.AssembleDiagonal(L_diag_unweighted);
       }
-      else
+
+      const double *d_L_diag_unweighted = L_diag_unweighted.Read();
+      double *d_L_diag = L_diag.ReadWrite();
+      MFEM_FORALL(i, L_diag.Size(),
       {
-         A_11.reset(new RAPOperator(*L_inv, *L, *L_inv));
-         const double *d_L_diag_unweighted = L_diag_unweighted.Read();
-         double *d_L_diag = L_diag.ReadWrite();
-         MFEM_FORALL(i, L_diag.Size(),
-         {
-            const double d = d_L_diag_unweighted[i];
-            d_L_diag[i] /= d*d;
-         });
-      }
+         const double d = d_L_diag_unweighted[i];
+         d_L_diag[i] /= d*d;
+      });
    }
 
    // Reassmble the RT mass operator with the new coefficient
@@ -185,6 +241,7 @@ void HdivSaddlePointSolver::Setup()
    // Form the updated approximate Schur complement
    mass_rt.AssembleDiagonal(R_diag);
 
+   // Update the mass RT diagonal for essential DOFs
    {
       const int *d_I = ess_rt_dofs.Read();
       double *d_R_diag = R_diag.ReadWrite();
@@ -230,7 +287,7 @@ void HdivSaddlePointSolver::Setup()
 void HdivSaddlePointSolver::EliminateBC(Vector &b) const
 {
    const int n_ess_dofs = ess_rt_dofs.Size();
-   if (n_ess_dofs == 0) { return; }
+   if (fes_l2.GetParMesh()->ReduceInt(n_ess_dofs) == 0) { return; }
 
    const int n_l2 = offsets[1];
    const int n_rt = offsets[2]-offsets[1];
@@ -238,7 +295,7 @@ void HdivSaddlePointSolver::EliminateBC(Vector &b) const
    Vector bF(b, n_l2, n_rt);
 
    // SetBC must be called first
-   MFEM_VERIFY(x_bc.Size() == n_rt, "BCs not set");
+   MFEM_VERIFY(x_bc.Size() == n_rt || n_ess_dofs == 0, "BCs not set");
 
    // Create a vector z that has the BC values at essential DOFs, zero elsewhere
    z.SetSize(n_rt);
@@ -282,6 +339,7 @@ void HdivSaddlePointSolver::EliminateBC(Vector &b) const
 
 void HdivSaddlePointSolver::Mult(const Vector &b, Vector &x) const
 {
+   w.SetSize(fes_l2.GetTrueVSize());
    b_prime.SetSize(b.Size());
    x_prime.SetSize(x.Size());
 
@@ -315,6 +373,7 @@ void HdivSaddlePointSolver::Mult(const Vector &b, Vector &x) const
    Vector xE(x, offsets[0], offsets[1]-offsets[0]);
    Vector xF(x, offsets[1], offsets[2]-offsets[1]);
 
+   z.SetSize(bE.Size()); // Size of z may have changed in EliminateBC
    L_inv->Mult(xE_prime, z);
 
    basis_l2.Mult(z, xE);

miniapps/hdiv-linear-solver/hdiv_linear_solver.hpp

@@ -55,6 +55,9 @@ private:
    ChangeOfBasis_L2 basis_l2;
    ChangeOfBasis_RT basis_rt;
 
+   /// Whether conversion from map type VALUE to INTEGRAL is required.
+   const bool convert_map_type;
+
    ParBilinearForm mass_l2, mass_rt;
 
    // Components needed for the block operator
@@ -79,13 +82,13 @@ private:
    /// The block-diagonal preconditioner.
    std::unique_ptr<BlockDiagonalPreconditioner> D_prec;
 
-   Coefficient &L_coeff, &R_coeff;
+   Coefficient &L_coeff, &R_coeff, &B_coeff;
 
    const Mode mode;
    bool zero_l2_block = false;
    QuadratureSpace qs;
-   QuadratureFunction qf;
-   QuadratureFunctionCoefficient l2_qf_coeff;
+   QuadratureFunction W_coeff_qf, W_mix_coeff_qf;
+   QuadratureFunctionCoefficient W_coeff, W_mix_coeff;
 
    ConstantCoefficient zero = ConstantCoefficient(0.0);
 
@@ -130,12 +133,22 @@ public:
    /// ParFiniteElementSpace&, Coefficient&, const Array<int>&) the zero-block
    /// HdivSaddlePointSolver constructor@endlink).
    HdivSaddlePointSolver(ParMesh &mesh_,
-                               ParFiniteElementSpace &fes_rt_,
-                               ParFiniteElementSpace &fes_l2_,
-                               Coefficient &L_coeff_,
-                               Coefficient &R_coeff_,
-                               const Array<int> &ess_rt_dofs_,
-                               Mode mode_);
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &L_coeff_,
+                         Coefficient &R_coeff_,
+                         Coefficient &B_coeff_,
+                         const Array<int> &ess_rt_dofs_,
+                         Mode mode_);
+
+   /// Same as the main constructor, but with B_coeff set to zero.
+   HdivSaddlePointSolver(ParMesh &mesh_,
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &L_coeff_,
+                         Coefficient &R_coeff_,
+                         const Array<int> &ess_rt_dofs_,
+                         Mode mode_);
 
    /// @brief Creates a linear solver for the case when the L2 diagonal block is
    /// zero (for Darcy problems).
@@ -145,10 +158,10 @@ public:
    /// ParFiniteElementSpace&, ParFiniteElementSpace&, Coefficient &,
    /// Coefficient&, const Array<int>&, Mode) the primary constructor@endlink.
    HdivSaddlePointSolver(ParMesh &mesh_,
-                               ParFiniteElementSpace &fes_rt_,
-                               ParFiniteElementSpace &fes_l2_,
-                               Coefficient &R_coeff_,
-                               const Array<int> &ess_rt_dofs_);
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &R_coeff_,
+                         const Array<int> &ess_rt_dofs_);
 
    /// @brief Build the linear operator and solver. Must be called when the
    /// coefficients change.

miniapps/hdiv-linear-solver/makefile

@@ -26,7 +26,7 @@ HDIV_HEADERS = hdiv_linear_solver.hpp discrete_divergence.hpp change_basis.hpp \
 HDIV_SRC = hdiv_linear_solver.cpp change_basis.cpp discrete_divergence.cpp
 HDIV_OBJ = $(HDIV_SRC:.cpp=.o)
 
-MINIAPPS = grad_div darcy
+MINIAPPS = grad_div darcy residual rz
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/hdiv-linear-solver/residual.cpp

@@ -0,0 +1,205 @@
+#include "mfem.hpp"
+#include <iostream>
+
+#include "hdiv_linear_solver.hpp"
+#include "discrete_divergence.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+
+double f(const Vector &xvec);
+double g(const Vector &xvec);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "../../data/star.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+   bool mt_value = true;
+   bool darcy = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&mt_value, "-val", "--value", "-int", "--integral",
+                  "Map type integral or value.");
+   args.AddOption(&darcy, "-da", "--darcy", "-g", "--grad-div",
+                  "Grad-div or Darcy problem");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   const int mt = mt_value ? FiniteElement::VALUE : FiniteElement::INTEGRAL;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   L2_FECollection fec_l2(order-1, dim, b2, mt);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+   ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+   HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
+   HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
+
+   if (Mpi::Root())
+   {
+      cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
+   }
+
+   Array<int> ess_rt_dofs;
+
+   FunctionCoefficient a_coeff(f);
+   FunctionCoefficient b_coeff(g);
+   ConstantCoefficient one(1.0);
+
+   Coefficient &div_coeff = darcy ? (Coefficient&)one : (Coefficient&)a_coeff;
+
+   // Solve the system with the saddle-point solver
+   const auto solver_mode = darcy ? HdivSaddlePointSolver::Mode::DARCY
+                            : HdivSaddlePointSolver::Mode::GRAD_DIV;
+   HdivSaddlePointSolver saddle_point_solver(
+      mesh, fes_rt, fes_l2, a_coeff, b_coeff, ess_rt_dofs, solver_mode);
+
+   const Array<int> &offsets = saddle_point_solver.GetOffsets();
+   BlockVector X_block(offsets), B_block(offsets);
+
+   saddle_point_solver.GetMINRES().SetAbsTol(1e-18);
+   saddle_point_solver.GetMINRES().SetRelTol(1e-20);
+   saddle_point_solver.GetMINRES().SetPrintLevel(
+      IterativeSolver::PrintLevel().FirstAndLast());
+   X_block = 0.0;
+   B_block.Randomize(1);
+   B_block.GetBlock(0) = 0.0;
+   if (Mpi::Root()) { std::cout << "Saddle point solver... " << std::endl; }
+   saddle_point_solver.Mult(B_block, X_block);
+
+   // Form the matrix-based system
+   ParBilinearForm w(&fes_l2);
+   w.AddDomainIntegrator(new MassIntegrator(a_coeff));
+   w.Assemble();
+   w.Finalize();
+   std::unique_ptr<HypreParMatrix> W(w.ParallelAssemble());
+
+   ParMixedBilinearForm b(&fes_rt, &fes_l2);
+   b.AddDomainIntegrator(new VectorFEDivergenceIntegrator(div_coeff));
+   b.Assemble();
+   b.Finalize();
+   std::unique_ptr<HypreParMatrix> B(b.ParallelAssemble());
+   std::unique_ptr<HypreParMatrix> Bt(B->Transpose());
+
+   ParBilinearForm m(&fes_rt);
+   m.AddDomainIntegrator(new VectorFEMassIntegrator(b_coeff));
+   m.Assemble();
+   m.Finalize();
+   std::unique_ptr<HypreParMatrix> M(m.ParallelAssemble());
+
+   BlockOperator A(offsets);
+   A.SetBlock(0, 0, W.get());
+   A.SetBlock(0, 1, B.get());
+   A.SetBlock(1, 0, Bt.get());
+   A.SetBlock(1, 1, M.get(), -1.0);
+
+   // Compute the residual
+   BlockVector Y_block(offsets);
+   A.Mult(X_block, Y_block);
+   Y_block -= B_block;
+
+   auto nrm2 = [](const Vector &x)
+   {
+      return sqrt(InnerProduct(MPI_COMM_WORLD, x, x));
+   };
+
+   const double resnorm1 = nrm2(Y_block)/nrm2(B_block);
+   if (Mpi::Root()) { std::cout << "Linear residual norm: " << resnorm1 << "\n\n"; }
+
+   // Solve the system with a matrix-based solver (see ex5p)
+   HypreParVector Md(MPI_COMM_WORLD, M->GetGlobalNumRows(),
+                     M->GetRowStarts());
+   M->GetDiag(Md);
+   std::unique_ptr<HypreParMatrix> MinvBt(B->Transpose());
+   MinvBt->InvScaleRows(Md);
+   std::unique_ptr<HypreParMatrix> S(ParMult(B.get(), MinvBt.get()));
+
+   HypreDiagScale M_inv(*M);
+   HypreBoomerAMG S_inv(*S);
+   S_inv.SetPrintLevel(0);
+
+   BlockDiagonalPreconditioner D(offsets);
+   D.SetDiagonalBlock(0, &S_inv);
+   D.SetDiagonalBlock(1, &M_inv);
+
+   X_block = 0.0;
+   MINRESSolver minres(MPI_COMM_WORLD);
+   minres.SetAbsTol(1e-18);
+   minres.SetRelTol(1e-20);
+   minres.SetMaxIter(500);
+   minres.SetOperator(A);
+   minres.SetPreconditioner(D);
+   minres.SetPrintLevel(IterativeSolver::PrintLevel().FirstAndLast());
+   if (Mpi::Root()) { std::cout << "Matrix-based solver... " << std::endl; }
+   minres.Mult(B_block, X_block);
+
+   A.Mult(X_block, Y_block);
+   Y_block -= B_block;
+
+   const double resnorm2 = nrm2(Y_block)/nrm2(B_block);
+   if (Mpi::Root()) { std::cout << "Linear residual norm: " << resnorm2 << "\n\n"; }
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}
+
+double f(const Vector &xvec)
+{
+   const int dim = xvec.Size();
+   const double x = xvec[0], y = xvec[1];
+   if (dim == 2)
+   {
+      return 2*(2.0 + sin(x)*sin(y));
+   }
+   else // dim == 3
+   {
+      const double z = xvec[2];
+      return 3*(2.0 + sin(x)*sin(y)*sin(z));
+   }
+}
+
+double g(const Vector &xvec)
+{
+   const int dim = xvec.Size();
+   const double x = xvec[0], y = xvec[1];
+   if (dim == 2)
+   {
+      return 2*(2.0 + cos(x)*cos(y));
+   }
+   else // dim == 3
+   {
+      const double z = xvec[2];
+      return 3*(2.0 + cos(x)*cos(y)*cos(z));
+   }
+}

miniapps/hdiv-linear-solver/rz.cpp

@@ -0,0 +1,207 @@
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+
+#include "discrete_divergence.hpp"
+#include "hdiv_linear_solver.hpp"
+
+#include "../solvers/lor_mms.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+
+double one_over_r(const Vector &xvec)
+{
+   // xvec = [z, r]
+   const double r = xvec[1];
+   return r == 0.0 ? 0.0 : 1.0/r;
+}
+
+double f_rz(const Vector &xvec)
+{
+   const double z = xvec[0];
+   const double r = xvec[1];
+
+   // alpha is the coefficient in the equation -Delta(u) + alpha*u = f
+   const double alpha = 1.0;
+   const double f = -cos(z)*(4*sin(r) + 5*r*cos(r) - (2 + alpha)*r*r*sin(r));
+
+   // scale integral by r because of coordinate transformation
+   return r*f;
+}
+
+double u_rz(const Vector &xvec)
+{
+   const double z = xvec[0];
+   const double r = xvec[1];
+
+   return r*r*sin(r)*cos(z);
+}
+
+class RobinCoefficient : public Coefficient
+{
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      double xdata[3];
+      Vector xvec(xdata, 3);
+      T.Transform(ip, xvec);
+      const int dim = xvec.Size();
+
+      Vector n(dim);
+      CalcOrtho(T.Jacobian(), n);
+      n /= n.Norml2();
+
+      const double p_val = u_rz(xvec);
+      const double z = xvec[0];
+      const double r = xvec[1];
+
+      if (dim == 2)
+      {
+         const double dpdz = -r*r*sin(r)*sin(z);
+         const double dpdr = r*cos(z)*(r*cos(r) + 2*sin(r));
+         const double u_val = n[0]*dpdz + n[1]*dpdr;
+         return p_val + u_val;
+      }
+      else
+      {
+         MFEM_ABORT("Not implemented");
+      }
+      return 0.0;
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "rz.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   const int mt = FiniteElement::VALUE;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   L2_FECollection fec_l2(order-1, dim, b2, mt);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+   ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+   HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
+   HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
+
+   if (Mpi::Root())
+   {
+      cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
+   }
+
+   Array<int> ess_rt_dofs; // empty
+
+   // f is the RHS, u is the exact solution
+   FunctionCoefficient f_coeff(f_rz), u_coeff(u_rz);
+
+   // Assemble the right-hand side for the scalar (L2) unknown.
+   ParLinearForm b_l2(&fes_l2);
+   // f_coeff has to include the r scaling for the coordinate transformation
+   b_l2.AddDomainIntegrator(new DomainLFIntegrator(f_coeff));
+   b_l2.UseFastAssembly(true);
+   b_l2.Assemble();
+
+   // Coefficient to enforce Robin boundary condition
+   RobinCoefficient bc_coeff;
+
+   // Enforce Robin boundary conditions by adding the boundary term to the flux
+   // equation.
+   ParLinearForm b_rt(&fes_rt);
+   b_rt.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(bc_coeff));
+   b_rt.UseFastAssembly(true);
+   b_rt.Assemble();
+
+   if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
+   tic_toc.Clear(); tic_toc.Start();
+
+   // Have to scale the RT mass matrix by (1/r)
+   FunctionCoefficient r_recip_coeff(one_over_r);
+   // Have to scale the L2 mass matrix by r
+   FunctionCoefficient r_coeff([](const Vector &xvec) { return xvec[1]; });
+
+   const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
+   HdivSaddlePointSolver saddle_point_solver(
+      mesh, fes_rt, fes_l2, r_coeff, r_recip_coeff, r_recip_coeff, ess_rt_dofs, solver_mode);
+
+   const Array<int> &offsets = saddle_point_solver.GetOffsets();
+   BlockVector X_block(offsets), B_block(offsets);
+
+   b_l2.ParallelAssemble(B_block.GetBlock(0));
+   b_rt.ParallelAssemble(B_block.GetBlock(1));
+   B_block.SyncFromBlocks();
+
+   X_block = 0.0;
+   saddle_point_solver.Mult(B_block, X_block);
+   X_block.SyncToBlocks();
+
+   if (Mpi::Root())
+   {
+      cout << "Done.\nIterations: "
+           << saddle_point_solver.GetNumIterations()
+           << "\nElapsed: " << tic_toc.RealTime() << endl;
+   }
+
+   ParGridFunction x(&fes_l2);
+   x.SetFromTrueDofs(X_block.GetBlock(0));
+   ParGridFunction flux(&fes_rt);
+   flux.SetFromTrueDofs(X_block.GetBlock(1));
+
+   const double error = x.ComputeL2Error(u_coeff);
+   if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+
+   ParGridFunction u_ex(&fes_l2), er(&fes_l2);
+   u_ex.ProjectCoefficient(u_coeff);
+   er = x;
+   er -= u_ex;
+
+   ParaViewDataCollection pv("RZ", &mesh);
+   pv.SetPrefixPath("ParaView");
+   pv.SetHighOrderOutput(true);
+   pv.SetLevelsOfDetail(order + 1);
+   pv.RegisterField("u", &x);
+   pv.RegisterField("flux", &flux);
+   pv.RegisterField("exact", &u_ex);
+   pv.RegisterField("error", &er);
+   pv.SetCycle(0);
+   pv.SetTime(0);
+   pv.Save();
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}

miniapps/hdiv-linear-solver/rz.mesh

@@ -0,0 +1,35 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+0 0
+12 0
+12 7.5
+0 7.5