Skip to content

Commit

Permalink
Construct coarse level multigrid operators without duplicating quadra…
Browse files Browse the repository at this point in the history 
…ture data
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Jan 9, 2024
1 parent 2cfeaf9 commit eebe04d
Show file tree
Hide file tree
Showing 6 changed files with 177 additions and 126 deletions.
14 changes: 8 additions & 6 deletions palace/linalg/divfree.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,27 +22,29 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, FiniteElementSpace
const std::vector<mfem::Array<int>> &h1_bdr_tdof_lists,
double tol, int max_it, int print)
{
constexpr bool skip_zeros = false;
MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
mat_op.GetPermittivityReal());
{
constexpr bool skip_zeros = false;
BilinearForm m(h1_fespaces.GetFinestFESpace());
m.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
// m.AssembleQuadratureData();
auto m_vec = m.Assemble(h1_fespaces, skip_zeros);
auto M_mg = std::make_unique<MultigridOperator>(h1_fespaces.GetNumLevels());
for (std::size_t l = 0; l < h1_fespaces.GetNumLevels(); l++)
{
// Force coarse level operator to be fully assembled always.
const auto &h1_fespace_l = h1_fespaces.GetFESpaceAtLevel(l);
BilinearForm m(h1_fespace_l);
m.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
auto M_l = std::make_unique<ParOperator>(m.Assemble(skip_zeros), h1_fespace_l);
auto M_l = std::make_unique<ParOperator>(std::move(m_vec[l]), h1_fespace_l);
M_l->SetEssentialTrueDofs(h1_bdr_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE);
M_mg->AddOperator(std::move(M_l));
}
M = std::move(M_mg);
}
{
// Weak divergence operator is always partially assembled.
BilinearForm weakdiv(nd_fespace, h1_fespaces.GetFinestFESpace());
weakdiv.AddDomainIntegrator<MixedVectorWeakDivergenceIntegrator>(epsilon_func);
WeakDiv = std::make_unique<ParOperator>(weakdiv.Assemble(skip_zeros), nd_fespace,
WeakDiv = std::make_unique<ParOperator>(weakdiv.PartialAssemble(), nd_fespace,
h1_fespaces.GetFinestFESpace(), false);
}
Grad = &h1_fespaces.GetFinestFESpace().GetDiscreteInterpolator();
Expand Down
20 changes: 9 additions & 11 deletions palace/linalg/hcurl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,25 +32,23 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver(
mat_op.GetInvPermeability());
MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
mat_op.GetPermittivityReal());
BilinearForm a(nd_fespaces.GetFinestFESpace()), a_aux(h1_fespaces.GetFinestFESpace());
a.AddDomainIntegrator<CurlCurlMassIntegrator>(muinv_func, epsilon_func);
a_aux.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
// a.AssembleQuadratureData();
// a_aux.AssembleQuadratureData();
auto a_vec = a.Assemble(nd_fespaces, skip_zeros);
auto a_aux_vec = a_aux.Assemble(h1_fespaces, skip_zeros);
auto A_mg = std::make_unique<MultigridOperator>(n_levels);
for (bool aux : {false, true})
{
for (std::size_t l = 0; l < n_levels; l++)
{
// Force coarse level operator to be fully assembled always.
const auto &fespace_l =
aux ? h1_fespaces.GetFESpaceAtLevel(l) : nd_fespaces.GetFESpaceAtLevel(l);
const auto &dbc_tdof_lists_l = aux ? h1_dbc_tdof_lists[l] : nd_dbc_tdof_lists[l];
BilinearForm a(fespace_l);
if (aux)
{
a.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
}
else
{
a.AddDomainIntegrator<CurlCurlMassIntegrator>(muinv_func, epsilon_func);
}
auto A_l = std::make_unique<ParOperator>(a.Assemble(skip_zeros), fespace_l);
auto A_l = std::make_unique<ParOperator>(std::move(aux ? a_aux_vec[l] : a_vec[l]),
fespace_l);
A_l->SetEssentialTrueDofs(dbc_tdof_lists_l, Operator::DiagonalPolicy::DIAG_ONE);
if (aux)
{
Expand Down
26 changes: 13 additions & 13 deletions palace/models/curlcurloperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,28 +154,26 @@ void PrintHeader(const mfem::ParFiniteElementSpace &h1_fespace,

std::unique_ptr<Operator> CurlCurlOperator::GetStiffnessMatrix()
{
// When partially assembled, the coarse operators can reuse the fine operator quadrature
// data if the spaces correspond to the same mesh.
PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);

constexpr bool skip_zeros = false;
MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
mat_op.GetInvPermeability());
BilinearForm k(GetNDSpace());
k.AddDomainIntegrator<CurlCurlIntegrator>(muinv_func);
k.AssembleQuadratureData();
auto k_vec = k.Assemble(GetNDSpaces(), skip_zeros);
auto K = std::make_unique<MultigridOperator>(GetNDSpaces().GetNumLevels());
for (std::size_t l = 0; l < GetNDSpaces().GetNumLevels(); l++)
{
// Force coarse level operator to be fully assembled always.
const auto &nd_fespace_l = GetNDSpaces().GetFESpaceAtLevel(l);
if (print_hdr)
{
Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l,
nd_fespace_l.GetMaxElementOrder(), nd_fespace_l.GlobalTrueVSize());
}
constexpr bool skip_zeros = false;
MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
mat_op.GetInvPermeability());
BilinearForm k(nd_fespace_l);
k.AddDomainIntegrator<CurlCurlIntegrator>(muinv_func);
auto K_l = std::make_unique<ParOperator>(
(l > 0) ? k.Assemble(skip_zeros) : k.FullAssemble(skip_zeros), nd_fespace_l);
if (print_hdr)
{
if (const auto *k_spm =
dynamic_cast<const mfem::SparseMatrix *>(&K_l->LocalOperator()))
if (const auto *k_spm = dynamic_cast<const mfem::SparseMatrix *>(k_vec[l].get()))
{
HYPRE_BigInt nnz = k_spm->NumNonZeroElems();
Mpi::GlobalSum(1, &nnz, nd_fespace_l.GetComm());
Expand All @@ -186,9 +184,11 @@ std::unique_ptr<Operator> CurlCurlOperator::GetStiffnessMatrix()
Mpi::Print("\n");
}
}
auto K_l = std::make_unique<ParOperator>(std::move(k_vec[l]), nd_fespace_l);
K_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE);
K->AddOperator(std::move(K_l));
}

print_hdr = false;
return K;
}
Expand Down
26 changes: 13 additions & 13 deletions palace/models/laplaceoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,28 +175,26 @@ void PrintHeader(const mfem::ParFiniteElementSpace &h1_fespace,

std::unique_ptr<Operator> LaplaceOperator::GetStiffnessMatrix()
{
// When partially assembled, the coarse operators can reuse the fine operator quadrature
// data if the spaces correspond to the same mesh.
PrintHeader(GetH1Space(), GetNDSpace(), print_hdr);

constexpr bool skip_zeros = false;
MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
mat_op.GetPermittivityReal());
BilinearForm k(GetH1Space());
k.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
k.AssembleQuadratureData();
auto k_vec = k.Assemble(GetH1Spaces(), skip_zeros);
auto K = std::make_unique<MultigridOperator>(GetH1Spaces().GetNumLevels());
for (std::size_t l = 0; l < GetH1Spaces().GetNumLevels(); l++)
{
// Force coarse level operator to be fully assembled always.
const auto &h1_fespace_l = GetH1Spaces().GetFESpaceAtLevel(l);
if (print_hdr)
{
Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l,
h1_fespace_l.GetMaxElementOrder(), h1_fespace_l.GlobalTrueVSize());
}
constexpr bool skip_zeros = false;
MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
mat_op.GetPermittivityReal());
BilinearForm k(h1_fespace_l);
k.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
auto K_l = std::make_unique<ParOperator>(
(l > 0) ? k.Assemble(skip_zeros) : k.FullAssemble(skip_zeros), h1_fespace_l);
if (print_hdr)
{
if (const auto *k_spm =
dynamic_cast<const mfem::SparseMatrix *>(&K_l->LocalOperator()))
if (const auto *k_spm = dynamic_cast<const mfem::SparseMatrix *>(k_vec[l].get()))
{
HYPRE_BigInt nnz = k_spm->NumNonZeroElems();
Mpi::GlobalSum(1, &nnz, h1_fespace_l.GetComm());
Expand All @@ -207,9 +205,11 @@ std::unique_ptr<Operator> LaplaceOperator::GetStiffnessMatrix()
Mpi::Print("\n");
}
}
auto K_l = std::make_unique<ParOperator>(std::move(k_vec[l]), h1_fespace_l);
K_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE);
K->AddOperator(std::move(K_l));
}

print_hdr = false;
return K;
}
Expand Down
Loading

0 comments on commit eebe04d

Please sign in to comment.