diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index 8f12cba7e..bc5445288 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -250,6 +250,13 @@ EigenSolver::Solve(const std::vector> &mesh) const ksp->SetOperators(*A, *P); eigen->SetLinearSolver(*ksp); + // Initialize structures for storing and reducing the results of error estimation. + TimeDependentFluxErrorEstimator 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> &mesh) const // Calculate and record the error indicators, and postprocess the results. Mpi::Print("\nComputing solution error estimates and performing postprocessing\n"); - TimeDependentFluxErrorEstimator 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 diff --git a/palace/fem/bilinearform.cpp b/palace/fem/bilinearform.cpp index 7fca63cf3..e589e1013 100644 --- a/palace/fem/bilinearform.cpp +++ b/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> 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 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. diff --git a/palace/fem/libceed/operator.cpp b/palace/fem/libceed/operator.cpp index 8fdc162bd..2fba1f7eb 100644 --- a/palace/fem/libceed/operator.cpp +++ b/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(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 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 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 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 } } diff --git a/palace/fem/libceed/operator.hpp b/palace/fem/libceed/operator.hpp index 72848fcea..b3d40d1f7 100644 --- a/palace/fem/libceed/operator.hpp +++ b/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; diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index bbf99fdc4..479818155 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -345,7 +345,7 @@ GradFluxErrorEstimator::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::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)); diff --git a/palace/linalg/operator.hpp b/palace/linalg/operator.hpp index f0512a769..f32cacda2 100644 --- a/palace/linalg/operator.hpp +++ b/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 { diff --git a/palace/linalg/vector.hpp b/palace/linalg/vector.hpp index 5e9248797..0163b48a9 100644 --- a/palace/linalg/vector.hpp +++ b/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 {