Merge pull request #219 from awslabs/sjg/estimator-libceed-dev

Use libCEED for element-wise error estimate integrals

palace/fem/integ/curlcurl.cpp

@@ -26,7 +26,7 @@ void CurlCurlIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                                   CeedElemRestriction geom_data_restr,
                                   CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/curlcurlmass.cpp

@@ -19,7 +19,7 @@ void CurlCurlMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr
                                       CeedElemRestriction geom_data_restr,
                                       CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.
@@ -59,8 +59,8 @@ void CurlCurlMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr
     info.trial_ops |= EvalMode::Weight;
   }
 
-  // Set up the coefficient and assemble.
-  auto ctx = PopulateCoefficientContext((dim < 3) ? 1 : dim, space_dim, Q, Q_mass);
+  // Set up the coefficient and assemble. Mass goes first.
+  auto ctx = PopulateCoefficientContext(space_dim, Q_mass, (dim < 3) ? 1 : dim, Q);
   AssembleCeedOperator(info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed,
                        trial_restr, test_restr, trial_basis, test_basis, geom_data,
                        geom_data_restr, op);

palace/fem/integ/diffusion.cpp

@@ -19,7 +19,7 @@ void DiffusionIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                                    CeedElemRestriction geom_data_restr,
                                    CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/diffusionmass.cpp

@@ -20,7 +20,7 @@ void DiffusionMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_rest
                                        CeedElemRestriction geom_data_restr,
                                        CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.
@@ -61,8 +61,8 @@ void DiffusionMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_rest
   info.trial_ops = EvalMode::Grad | EvalMode::Interp;
   info.test_ops = EvalMode::Grad | EvalMode::Interp;
 
-  // Set up the coefficient and assemble.
-  auto ctx = PopulateCoefficientContext(space_dim, 1, Q, Q_mass);
+  // Set up the coefficient and assemble. Mass goes first.
+  auto ctx = PopulateCoefficientContext(1, Q_mass, space_dim, Q);
   AssembleCeedOperator(info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed,
                        trial_restr, test_restr, trial_basis, test_basis, geom_data,
                        geom_data_restr, op);

palace/fem/integ/divdiv.cpp

@@ -18,7 +18,7 @@ void DivDivIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                                 CeedBasis test_basis, CeedVector geom_data,
                                 CeedElemRestriction geom_data_restr, CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/divdivmass.cpp

@@ -19,7 +19,7 @@ void DivDivMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                                     CeedElemRestriction geom_data_restr,
                                     CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.
@@ -60,8 +60,8 @@ void DivDivMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
   info.trial_ops = EvalMode::Div | EvalMode::Interp | EvalMode::Weight;
   info.test_ops = EvalMode::Div | EvalMode::Interp;
 
-  // Set up the coefficient and assemble.
-  auto ctx = PopulateCoefficientContext(1, space_dim, Q, Q_mass);
+  // Set up the coefficient and assemble. Mass goes first.
+  auto ctx = PopulateCoefficientContext(space_dim, Q_mass, 1, Q);
   AssembleCeedOperator(info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed,
                        trial_restr, test_restr, trial_basis, test_basis, geom_data,
                        geom_data_restr, op);

palace/fem/integ/grad.cpp

@@ -19,7 +19,7 @@ void GradientIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                                   CeedElemRestriction geom_data_restr,
                                   CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/mass.cpp

@@ -18,7 +18,7 @@ void MassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
                               CeedBasis test_basis, CeedVector geom_data,
                               CeedElemRestriction geom_data_restr, CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/mixedveccurl.cpp

@@ -27,7 +27,7 @@ void MixedVectorCurlIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_re
                                          CeedElemRestriction geom_data_restr,
                                          CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.
@@ -74,7 +74,7 @@ void MixedVectorWeakCurlIntegrator::Assemble(Ceed ceed, CeedElemRestriction tria
                                              CeedElemRestriction geom_data_restr,
                                              CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/mixedvecgrad.cpp

@@ -27,7 +27,7 @@ void MixedVectorGradientIntegrator::Assemble(Ceed ceed, CeedElemRestriction tria
                                              CeedElemRestriction geom_data_restr,
                                              CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.
@@ -132,7 +132,7 @@ void MixedVectorWeakDivergenceIntegrator::Assemble(
     CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
     CeedElemRestriction geom_data_restr, CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/integ/vecfemass.cpp

@@ -21,7 +21,7 @@ void VectorFEMassIntegrator::Assemble(Ceed ceed, CeedElemRestriction trial_restr
                                       CeedElemRestriction geom_data_restr,
                                       CeedOperator *op) const
 {
-  IntegratorInfo info;
+  CeedQFunctionInfo info;
   info.assemble_q_data = assemble_q_data;
 
   // Set up QFunctions.

palace/fem/libceed/basis.cpp

@@ -120,11 +120,10 @@ void InitMfemInterpolatorBasis(const mfem::FiniteElement &trial_fe,
   MFEM_VERIFY(trial_num_comp == test_num_comp && trial_num_comp == 1,
               "libCEED discrete linear operator requires same vdim = 1 for trial and test "
               "FE spaces!");
-  const int dim = trial_fe.GetDim();
   const int trial_P = trial_fe.GetDof();
   const int test_P = test_fe.GetDof();
-  mfem::DenseMatrix qX(dim, test_P), Gt(trial_P, test_P * dim), Bt;
-  mfem::Vector qW(test_P);
+  mfem::DenseMatrix Bt, Gt(trial_P, test_P);
+  mfem::Vector qX(test_P), qW(test_P);
   mfem::IsoparametricTransformation dummy;
   dummy.SetIdentityTransformation(trial_fe.GetGeomType());
   if (trial_fe.GetMapType() == test_fe.GetMapType())
@@ -159,9 +158,10 @@ void InitMfemInterpolatorBasis(const mfem::FiniteElement &trial_fe,
   qX = 0.0;
   qW = 0.0;
 
-  PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, GetCeedTopology(trial_fe.GetGeomType()),
-                                         trial_num_comp, trial_P, test_P, Bt.GetData(),
-                                         Gt.GetData(), qX.GetData(), qW.GetData(), basis));
+  // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
+  PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, CEED_TOPOLOGY_LINE, trial_num_comp, trial_P,
+                                         test_P, Bt.GetData(), Gt.GetData(), qX.GetData(),
+                                         qW.GetData(), basis));
 }
 
 }  // namespace

palace/fem/libceed/ceed.cpp

@@ -76,20 +76,22 @@ std::string Print()
   return std::string(ceed_resource);
 }
 
-void InitCeedVector(const mfem::Vector &v, Ceed ceed, CeedVector *cv)
+void InitCeedVector(const mfem::Vector &v, Ceed ceed, CeedVector *cv, bool init)
 {
   CeedMemType mem;
-  const CeedScalar *data;
-  PalaceCeedCall(ceed, CeedVectorCreate(ceed, v.Size(), cv));
   PalaceCeedCall(ceed, CeedGetPreferredMemType(ceed, &mem));
-  if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK) && mem == CEED_MEM_DEVICE)
+  if (!mfem::Device::Allows(mfem::Backend::DEVICE_MASK) && mem == CEED_MEM_DEVICE)
   {
-    data = v.Read();
+    mem = CEED_MEM_HOST;
+  }
+  const auto *data = v.Read(mem == CEED_MEM_DEVICE);
+  if (init)
+  {
+    PalaceCeedCall(ceed, CeedVectorCreate(ceed, v.Size(), cv));
   }
   else
   {
-    data = v.HostRead();
-    mem = CEED_MEM_HOST;
+    PalaceCeedCall(ceed, CeedVectorTakeArray(*cv, mem, nullptr));
   }
   PalaceCeedCall(
       ceed, CeedVectorSetArray(*cv, mem, CEED_USE_POINTER, const_cast<CeedScalar *>(data)));

palace/fem/libceed/ceed.hpp

@@ -54,8 +54,9 @@ void Finalize();
 // Get the configured libCEED backend.
 std::string Print();
 
-// Initialize a CeedVector from an mfem::Vector.
-void InitCeedVector(const mfem::Vector &v, Ceed ceed, CeedVector *cv);
+// Initialize a CeedVector from an mfem::Vector. When init is false, expects the CeedVector
+// has already been initialized and just sets the data pointer.
+void InitCeedVector(const mfem::Vector &v, Ceed ceed, CeedVector *cv, bool init = true);
 
 // Convert an MFEM geometry type to a libCEED one.
 CeedElemTopology GetCeedTopology(mfem::Geometry::Type geom);

palace/fem/libceed/coefficient.cpp

@@ -118,13 +118,12 @@ PopulateCoefficientContext(int dim, const MaterialPropertyCoefficient *Q, double
 }
 
 std::vector<CeedIntScalar>
-PopulateCoefficientContext(int dim, int dim_mass, const MaterialPropertyCoefficient *Q,
-                           const MaterialPropertyCoefficient *Q_mass, double a,
-                           double a_mass)
+PopulateCoefficientContext(int dim_mass, const MaterialPropertyCoefficient *Q_mass, int dim,
+                           const MaterialPropertyCoefficient *Q, double a_mass, double a)
 {
   // Mass coefficient comes first, then the other one for the QFunction.
-  auto ctx = PopulateCoefficientContext(dim, Q, a);
   auto ctx_mass = PopulateCoefficientContext(dim_mass, Q_mass, a_mass);
+  auto ctx = PopulateCoefficientContext(dim, Q, a);
   ctx_mass.insert(ctx_mass.end(), ctx.begin(), ctx.end());
   return ctx_mass;
 }

palace/fem/libceed/coefficient.hpp

@@ -20,9 +20,9 @@ std::vector<CeedIntScalar>
 PopulateCoefficientContext(int dim, const MaterialPropertyCoefficient *Q, double a = 1.0);
 
 std::vector<CeedIntScalar>
-PopulateCoefficientContext(int dim, int dim_mass, const MaterialPropertyCoefficient *Q,
-                           const MaterialPropertyCoefficient *Q_mass, double a = 1.0,
-                           double a_mass = 1.0);
+PopulateCoefficientContext(int dim_mass, const MaterialPropertyCoefficient *Q_mass, int dim,
+                           const MaterialPropertyCoefficient *Q, double a_mass = 1.0,
+                           double a = 1.0);
 
 }  // namespace ceed