99#include < mfem.hpp>
1010#include " linalg/operator.hpp"
1111#include " linalg/rap.hpp"
12+ #include " utils/iodata.hpp"
1213
1314namespace palace ::utils
1415{
@@ -17,61 +18,118 @@ namespace palace::utils
1718// Methods for constructing hierarchies of finite element spaces for geometric multigrid.
1819//
1920
21+ // Helper function for getting the order of the finite element space underlying a bilinear
22+ // form.
23+ inline auto GetMaxElementOrder (mfem::BilinearForm &a)
24+ {
25+ return a.FESpace ()->GetMaxElementOrder ();
26+ }
27+
28+ // Helper function for getting the order of the finite element space underlying a mixed
29+ // bilinear form.
30+ inline auto GetMaxElementOrder (mfem::MixedBilinearForm &a)
31+ {
32+ return std::max (a.TestFESpace ()->GetMaxElementOrder (),
33+ a.TrialFESpace ()->GetMaxElementOrder ());
34+ }
35+
36+ // Assemble a bilinear or mixed bilinear form. If the order is lower than the specified
37+ // threshold, the operator is assembled as a sparse matrix.
38+ template <typename BilinearForm>
39+ inline std::unique_ptr<Operator>
40+ AssembleOperator (std::unique_ptr<BilinearForm> &&a, bool mfem_pa_support,
41+ int pa_order_threshold, int skip_zeros = 1 )
42+ {
43+ mfem::AssemblyLevel assembly_level =
44+ (mfem::DeviceCanUseCeed () ||
45+ (mfem_pa_support && GetMaxElementOrder (*a) >= pa_order_threshold))
46+ ? mfem::AssemblyLevel::PARTIAL
47+ : mfem::AssemblyLevel::LEGACY;
48+ a->SetAssemblyLevel (assembly_level);
49+ a->Assemble (skip_zeros);
50+ a->Finalize (skip_zeros);
51+ if (assembly_level == mfem::AssemblyLevel::LEGACY ||
52+ (assembly_level == mfem::AssemblyLevel::PARTIAL &&
53+ GetMaxElementOrder (*a) < pa_order_threshold &&
54+ std::is_base_of<mfem::BilinearForm, BilinearForm>::value))
55+ {
56+ // libCEED full assembly does not support mixed forms.
57+ #ifdef MFEM_USE_CEED
58+ mfem::SparseMatrix *spm =
59+ a->HasExt () ? mfem::ceed::CeedOperatorFullAssemble (*a) : a->LoseMat ();
60+ #else
61+ mfem::SparseMatrix *spm = a->LoseMat ();
62+ #endif
63+ MFEM_VERIFY (spm, " Missing assembled sparse matrix!"
64+ return std::unique_ptr<Operator>(spm);
65+ }
66+ else
67+ {
68+ return std::move (a);
69+ }
70+ }
71+
2072// Construct sequence of FECollection objects.
2173template <typename FECollection>
22- std::vector<std::unique_ptr<FECollection>> ConstructFECollections (bool pc_pmg, bool pc_lor,
23- int p, int dim)
74+ std::vector<std::unique_ptr<FECollection>> inline ConstructFECollections (
75+ int p, int dim, int mg_max_levels,
76+ config::LinearSolverData::MultigridCoarsenType mg_coarsen_type, bool mat_lor)
2477{
2578 // If the solver will use a LOR preconditioner, we need to construct with a specific basis
2679 // type.
27- constexpr int pmin = (std::is_same<FECollection, mfem::H1_FECollection>::value ||
28- std::is_same<FECollection, mfem::ND_FECollection>::value)
80+ constexpr int pmin = (std::is_base_of< mfem::H1_FECollection, FECollection >::value ||
81+ std::is_base_of< mfem::ND_FECollection, FECollection >::value)
2982 ? 1
3083 : 0 ;
3184 MFEM_VERIFY (p >= pmin, " FE space order must not be less than " " !"
3285 int b1 = mfem::BasisType::GaussLobatto, b2 = mfem::BasisType::GaussLegendre;
33- if (pc_lor )
86+ if (mat_lor )
3487 {
3588 b2 = mfem::BasisType::IntegratedGLL;
3689 }
90+
91+ // Construct the p-multigrid hierarchy, first finest to coarsest and then reverse the
92+ // order.
3793 std::vector<std::unique_ptr<FECollection>> fecs;
38- if (pc_pmg)
39- {
40- fecs.reserve (p);
41- for (int o = pmin; o <= p; o++)
42- {
43- if constexpr (std::is_same<FECollection, mfem::ND_FECollection>::value ||
44- std::is_same<FECollection, mfem::RT_FECollection>::value)
45- {
46- fecs.push_back (std::make_unique<FECollection>(o, dim, b1, b2));
47- }
48- else
49- {
50- fecs.push_back (std::make_unique<FECollection>(o, dim, b1));
51- }
52- }
53- }
54- else
94+ for (int l = 0 ; l < std::max (1 , mg_max_levels); l++)
5595 {
56- fecs.reserve (1 );
57- if constexpr (std::is_same<FECollection, mfem::ND_FECollection>::value ||
58- std::is_same<FECollection, mfem::RT_FECollection>::value)
96+ if constexpr (std::is_base_of<mfem::ND_FECollection, FECollection>::value ||
97+ std::is_base_of<mfem::RT_FECollection, FECollection>::value)
5998 {
6099 fecs.push_back (std::make_unique<FECollection>(p, dim, b1, b2));
61100 }
62101 else
63102 {
64103 fecs.push_back (std::make_unique<FECollection>(p, dim, b1));
104+ MFEM_CONTRACT_VAR (b2);
105+ }
106+ if (p == pmin)
107+ {
108+ break ;
109+ }
110+ switch (mg_coarsen_type)
111+ {
112+ case config::LinearSolverData::MultigridCoarsenType::LINEAR:
113+ p--;
114+ break ;
115+ case config::LinearSolverData::MultigridCoarsenType::LOGARITHMIC:
116+ p = (p + pmin) / 2 ;
117+ break ;
118+ case config::LinearSolverData::MultigridCoarsenType::INVALID:
119+ MFEM_ABORT (" Invalid coarsening type for p-multigrid levels!"
120+ break ;
65121 }
66122 }
123+ std::reverse (fecs.begin (), fecs.end ());
67124 return fecs;
68125}
69126
70127// Construct a hierarchy of finite element spaces given a sequence of meshes and
71128// finite element collections. Dirichlet boundary conditions are additionally
72129// marked.
73130template <typename FECollection>
74- mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy (
131+ inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy (
132+ int mg_max_levels, bool mg_legacy_transfer, int pa_order_threshold,
75133 const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
76134 const std::vector<std::unique_ptr<FECollection>> &fecs,
77135 const mfem::Array<int > *dbc_marker = nullptr ,
@@ -80,17 +138,18 @@ mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
80138 MFEM_VERIFY (!mesh.empty () && !fecs.empty () &&
81139 (!dbc_tdof_lists || dbc_tdof_lists->empty ()),
82140 " Empty mesh or FE collection for FE space construction!"
83- auto *fespace = new mfem::ParFiniteElementSpace (mesh[0 ].get (), fecs[0 ].get ());
141+ int coarse_mesh_l =
142+ std::max (0 , static_cast <int >(mesh.size () + fecs.size ()) - 1 - mg_max_levels);
143+ auto *fespace = new mfem::ParFiniteElementSpace (mesh[coarse_mesh_l].get (), fecs[0 ].get ());
84144 if (dbc_marker && dbc_tdof_lists)
85145 {
86146 fespace->GetEssentialTrueDofs (*dbc_marker, dbc_tdof_lists->emplace_back ());
87147 }
88- mfem::ParFiniteElementSpaceHierarchy fespaces (mesh[0 ].get (), fespace, false , true );
89-
90- // XX TODO: LibCEED transfer operators!
148+ mfem::ParFiniteElementSpaceHierarchy fespaces (mesh[coarse_mesh_l].get (), fespace, false ,
149+ true );
91150
92151 // h-refinement
93- for (std::size_t l = 1 ; l < mesh.size (); l++)
152+ for (std::size_t l = coarse_mesh_l + 1 ; l < mesh.size (); l++)
94153 {
95154 fespace = new mfem::ParFiniteElementSpace (mesh[l].get (), fecs[0 ].get ());
96155 if (dbc_marker && dbc_tdof_lists)
@@ -111,31 +170,27 @@ mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
111170 {
112171 fespace->GetEssentialTrueDofs (*dbc_marker, dbc_tdof_lists->emplace_back ());
113172 }
114- auto *P = new ParOperator (
115- std::make_unique<mfem::TransferOperator>(fespaces.GetFinestFESpace (), *fespace),
116- fespaces.GetFinestFESpace (), *fespace, true );
173+ ParOperator *P;
174+ if (!mg_legacy_transfer && mfem::DeviceCanUseCeed ())
175+ {
176+ // Partial and full assembly for this operator is only available with libCEED backend.
177+ auto p = std::make_unique<mfem::DiscreteLinearOperator>(&fespaces.GetFinestFESpace (),
178+ fespace);
179+ p->AddDomainInterpolator (new mfem::IdentityInterpolator);
180+ P = new ParOperator (AssembleOperator (std::move (p), false , pa_order_threshold),
181+ fespaces.GetFinestFESpace (), *fespace, true );
182+ }
183+ else
184+ {
185+ P = new ParOperator (
186+ std::make_unique<mfem::TransferOperator>(fespaces.GetFinestFESpace (), *fespace),
187+ fespaces.GetFinestFESpace (), *fespace, true );
188+ }
117189 fespaces.AddLevel (mesh.back ().get (), fespace, P, false , true , true );
118190 }
119191 return fespaces;
120192}
121193
122- // Construct a single-level finite element space hierarchy from a single mesh and
123- // finite element collection. Unnecessary to pass the dirichlet boundary
124- // conditions as they need not be incorporated in any inter-space projectors.
125- template <typename FECollection>
126- mfem::ParFiniteElementSpaceHierarchy
127- ConstructFiniteElementSpaceHierarchy (mfem::ParMesh &mesh, const FECollection &fec,
128- const mfem::Array<int > *dbc_marker = nullptr ,
129- mfem::Array<int > *dbc_tdof_list = nullptr )
130- {
131- auto *fespace = new mfem::ParFiniteElementSpace (&mesh, &fec);
132- if (dbc_marker && dbc_tdof_list)
133- {
134- fespace->GetEssentialTrueDofs (*dbc_marker, *dbc_tdof_list);
135- }
136- return mfem::ParFiniteElementSpaceHierarchy (&mesh, fespace, false , true );
137- }
138-
139194} // namespace palace::utils
140195
141196#endif // PALACE_FEM_MULTIGRID_HPP
0 commit comments