Merge pull request #272 from awslabs/sjg/eigen-error-estimator-memory

Release memory associated with libCEED operator assembly

palace/drivers/eigensolver.cpp

@@ -250,6 +250,13 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   ksp->SetOperators(*A, *P);
   eigen->SetLinearSolver(*ksp);
 
+  // Initialize structures for storing and reducing the results of error estimation.
+  TimeDependentFluxErrorEstimator<ComplexVector> estimator(
+      space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
+  ErrorIndicator indicator;
+
   // Eigenvalue problem solve.
   BlockTimer bt1(Timer::EPS);
   Mpi::Print("\n");
@@ -267,11 +274,6 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 
   // Calculate and record the error indicators, and postprocess the results.
   Mpi::Print("\nComputing solution error estimates and performing postprocessing\n");
-  TimeDependentFluxErrorEstimator<ComplexVector> estimator(
-      space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
-      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
-      iodata.solver.linear.estimator_mg);
-  ErrorIndicator indicator;
   if (!KM)
   {
     // Normalize the finalized eigenvectors with respect to mass matrix (unit electric field

palace/fem/bilinearform.cpp

@@ -75,7 +75,7 @@ BilinearForm::PartialAssemble(const FiniteElementSpace &trial_fespace,
           integ->SetMapTypes(trial_map_type, test_map_type);
           integ->Assemble(ceed, trial_restr, test_restr, trial_basis, test_basis,
                           data.geom_data, data.geom_data_restr, &sub_op);
-          op->AddOper(sub_op);  // Sub-operator owned by ceed::Operator
+          op->AddSubOperator(sub_op);  // Sub-operator owned by ceed::Operator
         }
       }
       else if (mfem::Geometry::Dimension[geom] == mesh.Dimension() - 1 &&
@@ -95,7 +95,7 @@ BilinearForm::PartialAssemble(const FiniteElementSpace &trial_fespace,
           integ->SetMapTypes(trial_map_type, test_map_type);
           integ->Assemble(ceed, trial_restr, test_restr, trial_basis, test_basis,
                           data.geom_data, data.geom_data_restr, &sub_op);
-          op->AddOper(sub_op);  // Sub-operator owned by ceed::Operator
+          op->AddSubOperator(sub_op);  // Sub-operator owned by ceed::Operator
         }
       }
     }
@@ -181,13 +181,14 @@ BilinearForm::Assemble(const FiniteElementSpaceHierarchy &fespaces, bool skip_ze
     }
   }
 
-  // Construct the final operators using full or partial assemble as needed. Force the
-  // coarse-level operator to be fully assembled always.
+  // Construct the final operators using full or partial assemble as needed. We do not
+  // force the coarse-level operator to be fully assembled always, it will be only assembled
+  // as needed for parallel assembly.
   std::vector<std::unique_ptr<Operator>> ops;
   ops.reserve(fespaces.GetNumLevels() - l0);
   for (std::size_t l = l0; l < fespaces.GetNumLevels(); l++)
   {
-    if (l == 0 || UseFullAssembly(fespaces.GetFESpaceAtLevel(l), pa_order_threshold))
+    if (UseFullAssembly(fespaces.GetFESpaceAtLevel(l), pa_order_threshold))
     {
       ops.push_back(FullAssemble(*pa_ops[l - l0], skip_zeros));
     }
@@ -242,7 +243,7 @@ std::unique_ptr<ceed::Operator> DiscreteLinearOperator::PartialAssemble() const
         {
           CeedOperator sub_op, sub_op_t;
           interp->Assemble(ceed, trial_restr, test_restr, interp_basis, &sub_op, &sub_op_t);
-          op->AddOper(sub_op, sub_op_t);  // Sub-operator owned by ceed::Operator
+          op->AddSubOperator(sub_op, sub_op_t);  // Sub-operator owned by ceed::Operator
         }
 
         // Basis is owned by the operator.

palace/fem/libceed/operator.cpp

@@ -57,7 +57,7 @@ Operator::~Operator()
   }
 }
 
-void Operator::AddOper(CeedOperator sub_op, CeedOperator sub_op_t)
+void Operator::AddSubOperator(CeedOperator sub_op, CeedOperator sub_op_t)
 {
   // This should be called from within a OpenMP parallel region.
   const int id = utils::GetThreadNum();
@@ -100,6 +100,19 @@ void Operator::Finalize()
   }
 }
 
+void Operator::DestroyAssemblyData() const
+{
+  PalacePragmaOmp(parallel if (op.size() > 1))
+  {
+    const int id = utils::GetThreadNum();
+    MFEM_ASSERT(static_cast<std::size_t>(id) < op.size(),
+                "Out of bounds access for thread number " << id << "!");
+    Ceed ceed;
+    PalaceCeedCallBackend(CeedOperatorGetCeed(op[id], &ceed));
+    PalaceCeedCall(ceed, CeedOperatorAssemblyDataStrip(op[id]));
+  }
+}
+
 void Operator::AssembleDiagonal(Vector &diag) const
 {
   Ceed ceed;
@@ -125,6 +138,7 @@ void Operator::AssembleDiagonal(Vector &diag) const
     PalaceCeedCall(
         ceed, CeedOperatorLinearAssembleAddDiagonal(op[id], v[id], CEED_REQUEST_IMMEDIATE));
     PalaceCeedCall(ceed, CeedVectorTakeArray(v[id], mem, nullptr));
+    PalaceCeedCall(ceed, CeedOperatorAssemblyDataStrip(op[id]));
   }
 }
 
@@ -467,6 +481,7 @@ std::unique_ptr<hypre::HypreCSRMatrix> CeedOperatorFullAssemble(const Operator &
       CeedVector vals;
       CeedMemType mem;
       CeedOperatorAssembleCOO(ceed, op[id], skip_zeros, &nnz, &rows, &cols, &vals, &mem);
+      PalaceCeedCall(ceed, CeedOperatorAssemblyDataStrip(op[id]));
 
       // Convert COO to CSR (on each thread). The COO memory is free'd internally.
       loc_mat[id] =
@@ -527,6 +542,7 @@ std::unique_ptr<Operator> CeedOperatorCoarsen(const Operator &op_fine,
     PalaceCeedCall(ceed, CeedOperatorMultigridLevelCreate(op_fine, nullptr, restr_coarse,
                                                           basis_coarse, op_coarse, nullptr,
                                                           nullptr));
+    PalaceCeedCall(ceed, CeedOperatorAssemblyDataStrip(*op_coarse));
   };
 
   // Initialize the coarse operator.
@@ -558,7 +574,7 @@ std::unique_ptr<Operator> CeedOperatorCoarsen(const Operator &op_fine,
     {
       CeedOperator sub_op_coarse;
       SingleOperatorCoarsen(ceed, sub_ops_fine[k], &sub_op_coarse);
-      op_coarse->AddOper(sub_op_coarse);  // Sub-operator owned by ceed::Operator
+      op_coarse->AddSubOperator(sub_op_coarse);  // Sub-operator owned by ceed::Operator
     }
   }
 

palace/fem/libceed/operator.hpp

@@ -44,10 +44,12 @@ class Operator : public palace::Operator
   CeedOperator operator[](std::size_t i) const { return op[i]; }
   auto Size() const { return op.size(); }
 
-  void AddOper(CeedOperator sub_op, CeedOperator sub_op_t = nullptr);
+  void AddSubOperator(CeedOperator sub_op, CeedOperator sub_op_t = nullptr);
 
   void Finalize();
 
+  void DestroyAssemblyData() const;
+
   void SetDofMultiplicity(Vector &&mult) { dof_multiplicity = std::move(mult); }
 
   void AssembleDiagonal(Vector &diag) const override;

palace/linalg/errorestimator.cpp

@@ -345,7 +345,7 @@ GradFluxErrorEstimator<VecType>::GradFluxErrorEstimator(
           info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed, E_gf_vec,
           D_gf_vec, nd_restr, rt_restr, nd_basis, rt_basis, mesh_elem_restr, data.geom_data,
           data.geom_data_restr, &sub_op);
-      integ_op.AddOper(sub_op);  // Sub-operator owned by ceed::Operator
+      integ_op.AddSubOperator(sub_op);  // Sub-operator owned by ceed::Operator
 
       // Element restriction and passive input vectors are owned by the operator.
       PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&mesh_elem_restr));
@@ -458,7 +458,7 @@ CurlFluxErrorEstimator<VecType>::CurlFluxErrorEstimator(
           info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed, B_gf_vec,
           H_gf_vec, rt_restr, nd_restr, rt_basis, nd_basis, mesh_elem_restr, data.geom_data,
           data.geom_data_restr, &sub_op);
-      integ_op.AddOper(sub_op);  // Sub-operator owned by ceed::Operator
+      integ_op.AddSubOperator(sub_op);  // Sub-operator owned by ceed::Operator
 
       // Element restriction and passive input vectors are owned by the operator.
       PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&mesh_elem_restr));

palace/linalg/operator.hpp

@@ -14,12 +14,12 @@
 namespace palace
 {
 
-using Operator = mfem::Operator;
-
 //
 // Functionality extending mfem::Operator from MFEM.
 //
 
+using Operator = mfem::Operator;
+
 // Abstract base class for complex-valued operators.
 class ComplexOperator
 {

palace/linalg/vector.hpp

@@ -12,13 +12,13 @@
 namespace palace
 {
 
-using Vector = mfem::Vector;
-
 //
 // Functionality extending mfem::Vector from MFEM, including basic functions for parallel
 // vectors distributed across MPI processes.
 //
 
+using Vector = mfem::Vector;
+
 // A complex-valued vector represented as two real vectors, one for each component.
 class ComplexVector
 {