Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "bilinearform.hpp"
5 :
6 : #include "fem/fespace.hpp"
7 : #include "fem/libceed/basis.hpp"
8 : #include "fem/libceed/ceed.hpp"
9 : #include "fem/mesh.hpp"
10 : #include "utils/omp.hpp"
11 :
12 : namespace palace
13 : {
14 :
15 0 : void BilinearForm::AssembleQuadratureData()
16 : {
17 0 : for (auto &integ : domain_integs)
18 : {
19 : integ->AssembleQuadratureData();
20 : }
21 0 : for (auto &integ : boundary_integs)
22 : {
23 : integ->AssembleQuadratureData();
24 : }
25 0 : }
26 :
27 : std::unique_ptr<ceed::Operator>
28 10788 : BilinearForm::PartialAssemble(const FiniteElementSpace &trial_fespace,
29 : const FiniteElementSpace &test_fespace) const
30 : {
31 10788 : MFEM_VERIFY(&trial_fespace.GetMesh() == &test_fespace.GetMesh(),
32 : "Trial and test finite element spaces must correspond to the same mesh!");
33 : const auto &mesh = trial_fespace.GetMesh();
34 :
35 : // Initialize the operator.
36 10788 : std::unique_ptr<ceed::Operator> op;
37 10788 : if (&trial_fespace == &test_fespace)
38 : {
39 11856 : op = std::make_unique<ceed::SymmetricOperator>(test_fespace.GetVSize(),
40 11856 : trial_fespace.GetVSize());
41 : }
42 : else
43 : {
44 : op =
45 9720 : std::make_unique<ceed::Operator>(test_fespace.GetVSize(), trial_fespace.GetVSize());
46 : }
47 :
48 : // Assemble the libCEED operator in parallel, each thread builds a composite operator.
49 : // This should work fine if some threads create an empty operator (no elements or boundary
50 : // elements).
51 10788 : PalacePragmaOmp(parallel if (ceed::internal::NumCeeds() > 1))
52 : {
53 : Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
54 : for (const auto &[geom, data] : mesh.GetCeedGeomFactorData(ceed))
55 : {
56 : const auto trial_map_type =
57 : trial_fespace.GetFEColl().GetMapType(mfem::Geometry::Dimension[geom]);
58 : const auto test_map_type =
59 : test_fespace.GetFEColl().GetMapType(mfem::Geometry::Dimension[geom]);
60 :
61 : if (mfem::Geometry::Dimension[geom] == mesh.Dimension() && !domain_integs.empty())
62 : {
63 : // Assemble domain integrators on this element geometry type.
64 : CeedElemRestriction trial_restr =
65 : trial_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
66 : CeedElemRestriction test_restr =
67 : test_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
68 : CeedBasis trial_basis = trial_fespace.GetCeedBasis(ceed, geom);
69 : CeedBasis test_basis = test_fespace.GetCeedBasis(ceed, geom);
70 :
71 : for (const auto &integ : domain_integs)
72 : {
73 : CeedOperator sub_op;
74 : integ->SetMapTypes(trial_map_type, test_map_type);
75 : integ->Assemble(ceed, trial_restr, test_restr, trial_basis, test_basis,
76 : data.geom_data, data.geom_data_restr, &sub_op);
77 : op->AddSubOperator(sub_op); // Sub-operator owned by ceed::Operator
78 : }
79 : }
80 : else if (mfem::Geometry::Dimension[geom] == mesh.Dimension() - 1 &&
81 : !boundary_integs.empty())
82 : {
83 : // Assemble boundary integrators on this element geometry type.
84 : CeedElemRestriction trial_restr =
85 : trial_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
86 : CeedElemRestriction test_restr =
87 : test_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
88 : CeedBasis trial_basis = trial_fespace.GetCeedBasis(ceed, geom);
89 : CeedBasis test_basis = test_fespace.GetCeedBasis(ceed, geom);
90 :
91 : for (const auto &integ : boundary_integs)
92 : {
93 : CeedOperator sub_op;
94 : integ->SetMapTypes(trial_map_type, test_map_type);
95 : integ->Assemble(ceed, trial_restr, test_restr, trial_basis, test_basis,
96 : data.geom_data, data.geom_data_restr, &sub_op);
97 : op->AddSubOperator(sub_op); // Sub-operator owned by ceed::Operator
98 : }
99 : }
100 : }
101 : }
102 :
103 : // Finalize the operator (call CeedOperatorCheckReady).
104 10788 : op->Finalize();
105 :
106 10788 : return op;
107 : }
108 :
109 10950 : std::unique_ptr<hypre::HypreCSRMatrix> BilinearForm::FullAssemble(const ceed::Operator &op,
110 : bool skip_zeros, bool set)
111 : {
112 10950 : return ceed::CeedOperatorFullAssemble(op, skip_zeros, set);
113 : }
114 :
115 : namespace
116 : {
117 :
118 18 : bool UseFullAssembly(const FiniteElementSpace &trial_fespace,
119 : const FiniteElementSpace &test_fespace, int pa_order_threshold)
120 : {
121 : // Returns order such that the miniumum for all element types is 1. MFEM's
122 : // RT_FECollection actually already returns order + 1 for GetOrder() for historical
123 : // reasons.
124 : const auto &trial_fec = trial_fespace.GetFEColl();
125 : const auto &test_fec = test_fespace.GetFEColl();
126 : int max_order = std::max(
127 18 : dynamic_cast<const mfem::L2_FECollection *>(&trial_fec) ? trial_fec.GetOrder() + 1
128 : : trial_fec.GetOrder(),
129 18 : dynamic_cast<const mfem::L2_FECollection *>(&test_fec) ? test_fec.GetOrder() + 1
130 : : test_fec.GetOrder());
131 18 : return (max_order < pa_order_threshold);
132 : }
133 :
134 : bool UseFullAssembly(const FiniteElementSpace &fespace, int pa_order_threshold)
135 : {
136 0 : return UseFullAssembly(fespace, fespace, pa_order_threshold);
137 : }
138 :
139 : } // namespace
140 :
141 18 : std::unique_ptr<Operator> BilinearForm::Assemble(bool skip_zeros) const
142 : {
143 18 : if (UseFullAssembly(trial_fespace, test_fespace, pa_order_threshold))
144 : {
145 0 : return FullAssemble(skip_zeros);
146 : }
147 : else
148 : {
149 18 : return PartialAssemble();
150 : }
151 : }
152 :
153 : std::vector<std::unique_ptr<Operator>>
154 0 : BilinearForm::Assemble(const FiniteElementSpaceHierarchy &fespaces, bool skip_zeros,
155 : std::size_t l0) const
156 : {
157 : // Only available for square operators (same test and trial spaces).
158 0 : MFEM_VERIFY(&trial_fespace == &test_fespace &&
159 : &fespaces.GetFinestFESpace() == &trial_fespace,
160 : "Assembly on a FiniteElementSpaceHierarchy should have the same BilinearForm "
161 : "spaces and fine space of the hierarchy!");
162 :
163 : // First partially assemble all of the operators.
164 0 : MFEM_VERIFY(l0 < fespaces.GetNumLevels(),
165 : "No levels available for operator coarsening (l0 = " << l0 << ")!");
166 : std::vector<std::unique_ptr<ceed::Operator>> pa_ops;
167 0 : pa_ops.reserve(fespaces.GetNumLevels() - l0);
168 0 : for (std::size_t l = l0; l < fespaces.GetNumLevels(); l++)
169 : {
170 0 : if (l > l0 && &fespaces.GetFESpaceAtLevel(l).GetMesh() ==
171 0 : &fespaces.GetFESpaceAtLevel(l - 1).GetMesh())
172 : {
173 : pa_ops.push_back(
174 0 : ceed::CeedOperatorCoarsen(*pa_ops.back(), fespaces.GetFESpaceAtLevel(l)));
175 : }
176 : else
177 : {
178 : pa_ops.push_back(
179 0 : PartialAssemble(fespaces.GetFESpaceAtLevel(l), fespaces.GetFESpaceAtLevel(l)));
180 : }
181 : }
182 :
183 : // Construct the final operators using full or partial assemble as needed. We do not
184 : // force the coarse-level operator to be fully assembled always, it will be only assembled
185 : // as needed for parallel assembly.
186 : std::vector<std::unique_ptr<Operator>> ops;
187 0 : ops.reserve(fespaces.GetNumLevels() - l0);
188 0 : for (std::size_t l = l0; l < fespaces.GetNumLevels(); l++)
189 : {
190 0 : if (UseFullAssembly(fespaces.GetFESpaceAtLevel(l), pa_order_threshold))
191 : {
192 0 : ops.push_back(FullAssemble(*pa_ops[l - l0], skip_zeros));
193 : }
194 : else
195 : {
196 0 : ops.push_back(std::move(pa_ops[l - l0]));
197 : }
198 : }
199 :
200 0 : return ops;
201 0 : }
202 :
203 558 : std::unique_ptr<ceed::Operator> DiscreteLinearOperator::PartialAssemble() const
204 : {
205 558 : MFEM_VERIFY(&trial_fespace.GetMesh() == &test_fespace.GetMesh(),
206 : "Trial and test finite element spaces must correspond to the same mesh!");
207 : const auto &mesh = trial_fespace.GetMesh();
208 :
209 : // Initialize the operator.
210 : auto op =
211 558 : std::make_unique<ceed::Operator>(test_fespace.GetVSize(), trial_fespace.GetVSize());
212 :
213 : // Assemble the libCEED operator in parallel, each thread builds a composite operator.
214 : // This should work fine if some threads create an empty operator (no elements or boundary
215 : // elements).
216 558 : PalacePragmaOmp(parallel if (ceed::internal::NumCeeds() > 1))
217 : {
218 : Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
219 : for (const auto &[geom, data] : mesh.GetCeedGeomFactorData(ceed))
220 : {
221 : if (mfem::Geometry::Dimension[geom] == mesh.Dimension() && !domain_interps.empty())
222 : {
223 : // Assemble domain interpolators on this element geometry type.
224 : CeedElemRestriction trial_restr =
225 : trial_fespace.GetInterpCeedElemRestriction(ceed, geom, data.indices);
226 : CeedElemRestriction test_restr =
227 : test_fespace.GetInterpRangeCeedElemRestriction(ceed, geom, data.indices);
228 :
229 : // Construct the interpolator basis.
230 : CeedBasis interp_basis;
231 : const mfem::FiniteElement &trial_fe =
232 : *trial_fespace.GetFEColl().FiniteElementForGeometry(geom);
233 : const mfem::FiniteElement &test_fe =
234 : *test_fespace.GetFEColl().FiniteElementForGeometry(geom);
235 : const int trial_vdim = trial_fespace.GetVDim();
236 : const int test_vdim = test_fespace.GetVDim();
237 : ceed::InitInterpolatorBasis(trial_fe, test_fe, trial_vdim, test_vdim, ceed,
238 : &interp_basis);
239 :
240 : for (const auto &interp : domain_interps)
241 : {
242 : CeedOperator sub_op, sub_op_t;
243 : interp->Assemble(ceed, trial_restr, test_restr, interp_basis, &sub_op, &sub_op_t);
244 : op->AddSubOperator(sub_op, sub_op_t); // Sub-operator owned by ceed::Operator
245 : }
246 :
247 : // Basis is owned by the operator.
248 : PalaceCeedCall(ceed, CeedBasisDestroy(&interp_basis));
249 : }
250 : }
251 : }
252 :
253 : // Finalize the operator (call CeedOperatorCheckReady).
254 558 : op->Finalize();
255 :
256 : // Construct dof multiplicity vector for scaling to account for dofs shared between
257 : // elements (on host, then copy to device).
258 558 : Vector test_multiplicity(test_fespace.GetVSize());
259 558 : test_multiplicity = 0.0;
260 : auto *h_mult = test_multiplicity.HostReadWrite();
261 558 : PalacePragmaOmp(parallel)
262 : {
263 : mfem::Array<int> dofs;
264 : mfem::DofTransformation dof_trans;
265 : PalacePragmaOmp(for schedule(static))
266 : for (int i = 0; i < test_fespace.GetMesh().GetNE(); i++)
267 : {
268 : test_fespace.Get().GetElementVDofs(i, dofs, dof_trans);
269 : for (int j = 0; j < dofs.Size(); j++)
270 : {
271 : const int k = dofs[j];
272 : PalacePragmaOmp(atomic update)
273 : h_mult[(k >= 0) ? k : -1 - k] += 1.0;
274 : }
275 : }
276 : }
277 : test_multiplicity.UseDevice(true);
278 558 : test_multiplicity.Reciprocal();
279 : op->SetDofMultiplicity(std::move(test_multiplicity));
280 :
281 558 : return op;
282 : }
283 :
284 : } // namespace palace
|
