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 "spaceoperator.hpp"
5 :
6 : #include <set>
7 : #include <type_traits>
8 : #include "fem/bilinearform.hpp"
9 : #include "fem/coefficient.hpp"
10 : #include "fem/integrator.hpp"
11 : #include "fem/mesh.hpp"
12 : #include "fem/multigrid.hpp"
13 : #include "linalg/hypre.hpp"
14 : #include "linalg/rap.hpp"
15 : #include "utils/communication.hpp"
16 : #include "utils/geodata.hpp"
17 : #include "utils/iodata.hpp"
18 : #include "utils/prettyprint.hpp"
19 :
20 : namespace palace
21 : {
22 :
23 : using namespace std::complex_literals;
24 :
25 17 : SpaceOperator::SpaceOperator(const IoData &iodata,
26 17 : const std::vector<std::unique_ptr<Mesh>> &mesh)
27 17 : : pc_mat_real(iodata.solver.linear.pc_mat_real),
28 17 : pc_mat_shifted(iodata.solver.linear.pc_mat_shifted), print_hdr(true),
29 17 : print_prec_hdr(true), dbc_attr(SetUpBoundaryProperties(iodata, *mesh.back())),
30 17 : nd_fecs(fem::ConstructFECollections<mfem::ND_FECollection>(
31 17 : iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels,
32 17 : iodata.solver.linear.mg_coarsening, false)),
33 17 : h1_fecs(fem::ConstructFECollections<mfem::H1_FECollection>(
34 17 : iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels,
35 17 : iodata.solver.linear.mg_coarsening, false)),
36 0 : rt_fecs(fem::ConstructFECollections<mfem::RT_FECollection>(
37 17 : iodata.solver.order - 1, mesh.back()->Dimension(),
38 17 : iodata.solver.linear.estimator_mg ? iodata.solver.linear.mg_max_levels : 1,
39 17 : iodata.solver.linear.mg_coarsening, false)),
40 17 : nd_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
41 17 : iodata.solver.linear.mg_max_levels, mesh, nd_fecs, &dbc_attr, &nd_dbc_tdof_lists)),
42 17 : h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
43 17 : iodata.solver.linear.mg_max_levels, mesh, h1_fecs, &dbc_attr, &h1_dbc_tdof_lists)),
44 17 : rt_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::RT_FECollection>(
45 17 : iodata.solver.linear.estimator_mg ? iodata.solver.linear.mg_max_levels : 1, mesh,
46 17 : rt_fecs)),
47 17 : mat_op(iodata, *mesh.back()), farfield_op(iodata, mat_op, *mesh.back()),
48 17 : surf_sigma_op(iodata, mat_op, *mesh.back()), surf_z_op(iodata, mat_op, *mesh.back()),
49 17 : lumped_port_op(iodata, mat_op, *mesh.back()),
50 17 : wave_port_op(iodata, mat_op, GetNDSpace(), GetH1Space()),
51 17 : surf_j_op(iodata, *mesh.back()),
52 17 : port_excitation_helper(lumped_port_op, wave_port_op, surf_j_op)
53 : {
54 : // Check Excitations.
55 17 : if (iodata.problem.type == ProblemType::DRIVEN)
56 : {
57 11 : MFEM_VERIFY(!port_excitation_helper.Empty(),
58 : "Driven problems must specify at least one excitation!");
59 : }
60 6 : else if (iodata.problem.type == ProblemType::EIGENMODE)
61 : {
62 3 : MFEM_VERIFY(port_excitation_helper.Empty(),
63 : "Eigenmode problems must not specify any excitation!");
64 : }
65 3 : else if (iodata.problem.type == ProblemType::TRANSIENT)
66 : {
67 3 : MFEM_VERIFY(
68 : port_excitation_helper.Size() == 1,
69 : "Transient problems currently only support a single excitation per simulation!");
70 : }
71 : else
72 : {
73 0 : MFEM_ABORT("Internal Error: Solver type incompatible with SpaceOperator.");
74 : }
75 :
76 : // Finalize setup.
77 17 : CheckBoundaryProperties();
78 :
79 : // Print essential BC information.
80 17 : if (dbc_attr.Size())
81 : {
82 9 : Mpi::Print("\nConfiguring Dirichlet PEC BC at attributes:\n");
83 18 : utils::PrettyPrint(dbc_attr);
84 : }
85 17 : }
86 :
87 17 : mfem::Array<int> SpaceOperator::SetUpBoundaryProperties(const IoData &iodata,
88 : const mfem::ParMesh &mesh)
89 : {
90 : // Check that boundary attributes have been specified correctly.
91 17 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
92 : mfem::Array<int> bdr_attr_marker;
93 17 : if (!iodata.boundaries.pec.empty())
94 : {
95 : bdr_attr_marker.SetSize(bdr_attr_max);
96 : bdr_attr_marker = 0;
97 63 : for (auto attr : mesh.bdr_attributes)
98 : {
99 54 : bdr_attr_marker[attr - 1] = 1;
100 : }
101 : std::set<int> bdr_warn_list;
102 18 : for (auto attr : iodata.boundaries.pec.attributes)
103 : {
104 : // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
105 : // "PEC boundary attribute tags must be non-negative and correspond to "
106 : // "attributes in the mesh!");
107 : // MFEM_VERIFY(bdr_attr_marker[attr - 1],
108 : // "Unknown PEC boundary attribute " << attr << "!");
109 9 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
110 : {
111 : bdr_warn_list.insert(attr);
112 : }
113 : }
114 9 : if (!bdr_warn_list.empty())
115 : {
116 0 : Mpi::Print("\n");
117 0 : Mpi::Warning("Unknown PEC boundary attributes!\nSolver will just ignore them!");
118 0 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
119 0 : Mpi::Print("\n");
120 : }
121 : }
122 :
123 : // Mark selected boundary attributes from the mesh as essential (Dirichlet).
124 : mfem::Array<int> dbc_bcs;
125 17 : dbc_bcs.Reserve(static_cast<int>(iodata.boundaries.pec.attributes.size()));
126 26 : for (auto attr : iodata.boundaries.pec.attributes)
127 : {
128 9 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
129 : {
130 0 : continue; // Can just ignore if wrong
131 : }
132 9 : dbc_bcs.Append(attr);
133 : }
134 17 : return dbc_bcs;
135 : }
136 :
137 17 : void SpaceOperator::CheckBoundaryProperties()
138 : {
139 : // Mark selected boundary attributes from the mesh as having some Dirichlet, Neumann, or
140 : // mixed BC applied.
141 : const mfem::ParMesh &mesh = GetMesh();
142 17 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
143 17 : const auto dbc_marker = mesh::AttrToMarker(bdr_attr_max, dbc_attr);
144 17 : const auto farfield_marker = mesh::AttrToMarker(bdr_attr_max, farfield_op.GetAttrList());
145 : const auto surf_sigma_marker =
146 17 : mesh::AttrToMarker(bdr_attr_max, surf_sigma_op.GetAttrList());
147 17 : const auto surf_z_Rs_marker = mesh::AttrToMarker(bdr_attr_max, surf_z_op.GetRsAttrList());
148 17 : const auto surf_z_Ls_marker = mesh::AttrToMarker(bdr_attr_max, surf_z_op.GetLsAttrList());
149 : const auto lumped_port_Rs_marker =
150 17 : mesh::AttrToMarker(bdr_attr_max, lumped_port_op.GetRsAttrList());
151 : const auto lumped_port_Ls_marker =
152 17 : mesh::AttrToMarker(bdr_attr_max, lumped_port_op.GetLsAttrList());
153 : const auto wave_port_marker =
154 34 : mesh::AttrToMarker(bdr_attr_max, wave_port_op.GetAttrList());
155 : mfem::Array<int> aux_bdr_marker(dbc_marker.Size());
156 263 : for (int i = 0; i < dbc_marker.Size(); i++)
157 : {
158 246 : aux_bdr_marker[i] =
159 237 : (dbc_marker[i] || farfield_marker[i] || surf_sigma_marker[i] ||
160 237 : surf_z_Rs_marker[i] || surf_z_Ls_marker[i] || lumped_port_Rs_marker[i] ||
161 516 : lumped_port_Ls_marker[i] || wave_port_marker[i]);
162 246 : if (aux_bdr_marker[i])
163 : {
164 87 : aux_bdr_attr.Append(i + 1);
165 : }
166 : }
167 : // aux_bdr_marker = 1; // Mark all boundaries (including material interfaces
168 : // // added during mesh preprocessing)
169 : // // As tested, this does not eliminate all DC modes!
170 34 : for (std::size_t l = 0; l < GetH1Spaces().GetNumLevels(); l++)
171 : {
172 17 : GetH1Spaces().GetFESpaceAtLevel(l).Get().GetEssentialTrueDofs(
173 17 : aux_bdr_marker, aux_bdr_tdof_lists.emplace_back());
174 : }
175 :
176 : // A final check that no boundary attribute is assigned multiple boundary conditions.
177 17 : const auto surf_z_marker = mesh::AttrToMarker(bdr_attr_max, surf_z_op.GetAttrList());
178 : const auto lumped_port_marker =
179 17 : mesh::AttrToMarker(bdr_attr_max, lumped_port_op.GetAttrList());
180 17 : const auto surf_j_marker = mesh::AttrToMarker(bdr_attr_max, surf_j_op.GetAttrList());
181 263 : for (int i = 0; i < dbc_marker.Size(); i++)
182 : {
183 246 : MFEM_VERIFY(dbc_marker[i] + farfield_marker[i] + surf_sigma_marker[i] +
184 : surf_z_marker[i] + lumped_port_marker[i] + wave_port_marker[i] +
185 : surf_j_marker[i] <=
186 : 1,
187 : "Boundary attributes should not be specified with multiple BC!");
188 : }
189 17 : }
190 :
191 : namespace
192 : {
193 :
194 0 : void PrintHeader(const mfem::ParFiniteElementSpace &h1_fespace,
195 : const mfem::ParFiniteElementSpace &nd_fespace,
196 : const mfem::ParFiniteElementSpace &rt_fespace, bool &print_hdr)
197 : {
198 0 : if (print_hdr)
199 : {
200 0 : Mpi::Print("\nAssembling system matrices, number of global unknowns:\n"
201 : " H1 (p = {:d}): {:d}, ND (p = {:d}): {:d}, RT (p = {:d}): {:d}\n Operator "
202 : "assembly level: {}\n",
203 0 : h1_fespace.GetMaxElementOrder(), h1_fespace.GlobalTrueVSize(),
204 0 : nd_fespace.GetMaxElementOrder(), nd_fespace.GlobalTrueVSize(),
205 0 : rt_fespace.GetMaxElementOrder(), rt_fespace.GlobalTrueVSize(),
206 0 : (nd_fespace.GetMaxElementOrder() >= BilinearForm::pa_order_threshold)
207 0 : ? "Partial"
208 : : "Full");
209 :
210 : const auto &mesh = *nd_fespace.GetParMesh();
211 0 : const auto geom_types = mesh::CheckElements(mesh).GetGeomTypes();
212 0 : Mpi::Print(" Mesh geometries:\n");
213 0 : for (auto geom : geom_types)
214 : {
215 0 : const auto *fe = nd_fespace.FEColl()->FiniteElementForGeometry(geom);
216 0 : MFEM_VERIFY(fe, "MFEM does not support ND spaces on geometry = "
217 : << mfem::Geometry::Name[geom] << "!");
218 0 : const int q_order = fem::DefaultIntegrationOrder::Get(mesh, geom);
219 0 : Mpi::Print(" {}: P = {:d}, Q = {:d} (quadrature order = {:d}){}\n",
220 0 : mfem::Geometry::Name[geom], fe->GetDof(),
221 0 : mfem::IntRules.Get(geom, q_order).GetNPoints(), q_order,
222 0 : (geom == geom_types.back()) ? "" : ",");
223 : }
224 : }
225 0 : print_hdr = false;
226 0 : }
227 :
228 0 : void AddIntegrators(BilinearForm &a, const MaterialPropertyCoefficient *df,
229 : const MaterialPropertyCoefficient *f,
230 : const MaterialPropertyCoefficient *dfb,
231 : const MaterialPropertyCoefficient *fb,
232 : const MaterialPropertyCoefficient *fp, bool assemble_q_data = false)
233 : {
234 0 : if (df && !df->empty() && f && !f->empty())
235 : {
236 0 : a.AddDomainIntegrator<CurlCurlMassIntegrator>(*df, *f);
237 : }
238 : else
239 : {
240 0 : if (df && !df->empty())
241 : {
242 0 : a.AddDomainIntegrator<CurlCurlIntegrator>(*df);
243 : }
244 0 : if (f && !f->empty())
245 : {
246 0 : a.AddDomainIntegrator<VectorFEMassIntegrator>(*f);
247 : }
248 : }
249 0 : if (dfb && !dfb->empty() && fb && !fb->empty())
250 : {
251 0 : a.AddBoundaryIntegrator<CurlCurlMassIntegrator>(*dfb, *fb);
252 : }
253 : else
254 : {
255 0 : if (dfb && !dfb->empty())
256 : {
257 0 : a.AddBoundaryIntegrator<CurlCurlIntegrator>(*dfb);
258 : }
259 0 : if (fb && !fb->empty())
260 : {
261 0 : a.AddBoundaryIntegrator<VectorFEMassIntegrator>(*fb);
262 : }
263 : }
264 0 : if (fp && !fp->empty())
265 : {
266 0 : a.AddDomainIntegrator<MixedVectorWeakCurlIntegrator>(*fp);
267 0 : a.AddDomainIntegrator<MixedVectorCurlIntegrator>(*fp, true);
268 : }
269 0 : if (assemble_q_data)
270 : {
271 0 : a.AssembleQuadratureData();
272 : }
273 0 : }
274 :
275 0 : void AddAuxIntegrators(BilinearForm &a, const MaterialPropertyCoefficient *f,
276 : const MaterialPropertyCoefficient *fb, bool assemble_q_data = false)
277 : {
278 0 : if (f && !f->empty())
279 : {
280 0 : a.AddDomainIntegrator<DiffusionIntegrator>(*f);
281 : }
282 0 : if (fb && !fb->empty())
283 : {
284 0 : a.AddBoundaryIntegrator<DiffusionIntegrator>(*fb);
285 : }
286 0 : if (assemble_q_data)
287 : {
288 0 : a.AssembleQuadratureData();
289 : }
290 0 : }
291 :
292 0 : auto AssembleOperator(const FiniteElementSpace &fespace,
293 : const MaterialPropertyCoefficient *df,
294 : const MaterialPropertyCoefficient *f,
295 : const MaterialPropertyCoefficient *dfb,
296 : const MaterialPropertyCoefficient *fb,
297 : const MaterialPropertyCoefficient *fp, bool skip_zeros = false,
298 : bool assemble_q_data = false)
299 : {
300 : BilinearForm a(fespace);
301 0 : AddIntegrators(a, df, f, dfb, fb, fp, assemble_q_data);
302 0 : return a.Assemble(skip_zeros);
303 : }
304 :
305 0 : auto AssembleOperators(const FiniteElementSpaceHierarchy &fespaces,
306 : const MaterialPropertyCoefficient *df,
307 : const MaterialPropertyCoefficient *f,
308 : const MaterialPropertyCoefficient *dfb,
309 : const MaterialPropertyCoefficient *fb,
310 : const MaterialPropertyCoefficient *fp, bool skip_zeros = false,
311 : bool assemble_q_data = false, std::size_t l0 = 0)
312 : {
313 : BilinearForm a(fespaces.GetFinestFESpace());
314 0 : AddIntegrators(a, df, f, dfb, fb, fp, assemble_q_data);
315 0 : return a.Assemble(fespaces, skip_zeros, l0);
316 : }
317 :
318 0 : auto AssembleAuxOperators(const FiniteElementSpaceHierarchy &fespaces,
319 : const MaterialPropertyCoefficient *f,
320 : const MaterialPropertyCoefficient *fb, bool skip_zeros = false,
321 : bool assemble_q_data = false, std::size_t l0 = 0)
322 : {
323 : BilinearForm a(fespaces.GetFinestFESpace());
324 0 : AddAuxIntegrators(a, f, fb, assemble_q_data);
325 0 : return a.Assemble(fespaces, skip_zeros, l0);
326 : }
327 :
328 : } // namespace
329 :
330 : template <typename OperType>
331 : std::unique_ptr<OperType>
332 0 : SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy)
333 : {
334 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
335 0 : MaterialPropertyCoefficient df(mat_op.MaxCeedAttribute()), f(mat_op.MaxCeedAttribute()),
336 0 : fb(mat_op.MaxCeedBdrAttribute()), fc(mat_op.MaxCeedAttribute());
337 0 : AddStiffnessCoefficients(1.0, df, f);
338 0 : AddStiffnessBdrCoefficients(1.0, fb);
339 0 : AddRealPeriodicCoefficients(1.0, f);
340 0 : AddImagPeriodicCoefficients(1.0, fc);
341 0 : int empty[2] = {(df.empty() && f.empty() && fb.empty()), (fc.empty())};
342 : Mpi::GlobalMin(2, empty, GetComm());
343 0 : if (empty[0] && empty[1])
344 : {
345 : return {};
346 : }
347 : constexpr bool skip_zeros = false;
348 0 : std::unique_ptr<Operator> kr, ki;
349 0 : if (!empty[0])
350 : {
351 0 : kr = AssembleOperator(GetNDSpace(), &df, &f, nullptr, &fb, nullptr, skip_zeros);
352 : }
353 0 : if (!empty[1])
354 : {
355 : ki =
356 0 : AssembleOperator(GetNDSpace(), nullptr, nullptr, nullptr, nullptr, &fc, skip_zeros);
357 : }
358 : if constexpr (std::is_same<OperType, ComplexOperator>::value)
359 : {
360 0 : auto K =
361 : std::make_unique<ComplexParOperator>(std::move(kr), std::move(ki), GetNDSpace());
362 0 : K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
363 : return K;
364 : }
365 : else
366 : {
367 0 : MFEM_VERIFY(!ki, "Unexpected imaginary part in GetStiffnessMatrix<Operator>!");
368 0 : auto K = std::make_unique<ParOperator>(std::move(kr), GetNDSpace());
369 0 : K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
370 : return K;
371 : }
372 0 : }
373 :
374 : template <typename OperType>
375 : std::unique_ptr<OperType>
376 0 : SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy diag_policy)
377 : {
378 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
379 0 : MaterialPropertyCoefficient f(mat_op.MaxCeedAttribute()),
380 0 : fb(mat_op.MaxCeedBdrAttribute());
381 0 : AddDampingCoefficients(1.0, f);
382 0 : AddDampingBdrCoefficients(1.0, fb);
383 0 : int empty = (f.empty() && fb.empty());
384 : Mpi::GlobalMin(1, &empty, GetComm());
385 0 : if (empty)
386 : {
387 : return {};
388 : }
389 : constexpr bool skip_zeros = false;
390 0 : auto c = AssembleOperator(GetNDSpace(), nullptr, &f, nullptr, &fb, nullptr, skip_zeros);
391 : if constexpr (std::is_same<OperType, ComplexOperator>::value)
392 : {
393 0 : auto C = std::make_unique<ComplexParOperator>(std::move(c), nullptr, GetNDSpace());
394 0 : C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
395 : return C;
396 : }
397 : else
398 : {
399 0 : auto C = std::make_unique<ParOperator>(std::move(c), GetNDSpace());
400 0 : C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
401 : return C;
402 : }
403 0 : }
404 :
405 : template <typename OperType>
406 0 : std::unique_ptr<OperType> SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy diag_policy)
407 : {
408 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
409 0 : MaterialPropertyCoefficient fr(mat_op.MaxCeedAttribute()), fi(mat_op.MaxCeedAttribute()),
410 0 : fbr(mat_op.MaxCeedBdrAttribute()), fbi(mat_op.MaxCeedBdrAttribute());
411 0 : AddRealMassCoefficients(1.0, fr);
412 0 : AddRealMassBdrCoefficients(1.0, fbr);
413 : if constexpr (std::is_same<OperType, ComplexOperator>::value)
414 : {
415 0 : AddImagMassCoefficients(1.0, fi);
416 : }
417 0 : int empty[2] = {(fr.empty() && fbr.empty()), (fi.empty() && fbi.empty())};
418 : Mpi::GlobalMin(2, empty, GetComm());
419 0 : if (empty[0] && empty[1])
420 : {
421 : return {};
422 : }
423 : constexpr bool skip_zeros = false;
424 0 : std::unique_ptr<Operator> mr, mi;
425 0 : if (!empty[0])
426 : {
427 0 : mr = AssembleOperator(GetNDSpace(), nullptr, &fr, nullptr, &fbr, nullptr, skip_zeros);
428 : }
429 0 : if (!empty[1])
430 : {
431 0 : mi = AssembleOperator(GetNDSpace(), nullptr, &fi, nullptr, &fbi, nullptr, skip_zeros);
432 : }
433 : if constexpr (std::is_same<OperType, ComplexOperator>::value)
434 : {
435 0 : auto M =
436 : std::make_unique<ComplexParOperator>(std::move(mr), std::move(mi), GetNDSpace());
437 0 : M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
438 : return M;
439 : }
440 : else
441 : {
442 0 : auto M = std::make_unique<ParOperator>(std::move(mr), GetNDSpace());
443 0 : M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
444 : return M;
445 : }
446 0 : }
447 :
448 : template <typename OperType>
449 : std::unique_ptr<OperType>
450 0 : SpaceOperator::GetExtraSystemMatrix(double omega, Operator::DiagonalPolicy diag_policy)
451 : {
452 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
453 0 : MaterialPropertyCoefficient dfbr(mat_op.MaxCeedBdrAttribute()),
454 0 : dfbi(mat_op.MaxCeedBdrAttribute()), fbr(mat_op.MaxCeedBdrAttribute()),
455 0 : fbi(mat_op.MaxCeedBdrAttribute());
456 0 : AddExtraSystemBdrCoefficients(omega, dfbr, dfbi, fbr, fbi);
457 0 : int empty[2] = {(dfbr.empty() && fbr.empty()), (dfbi.empty() && fbi.empty())};
458 : Mpi::GlobalMin(2, empty, GetComm());
459 0 : if (empty[0] && empty[1])
460 : {
461 : return {};
462 : }
463 : constexpr bool skip_zeros = false;
464 0 : std::unique_ptr<Operator> ar, ai;
465 0 : if (!empty[0])
466 : {
467 0 : ar = AssembleOperator(GetNDSpace(), nullptr, nullptr, &dfbr, &fbr, nullptr, skip_zeros);
468 : }
469 0 : if (!empty[1])
470 : {
471 0 : ai = AssembleOperator(GetNDSpace(), nullptr, nullptr, &dfbi, &fbi, nullptr, skip_zeros);
472 : }
473 : if constexpr (std::is_same<OperType, ComplexOperator>::value)
474 : {
475 0 : auto A =
476 : std::make_unique<ComplexParOperator>(std::move(ar), std::move(ai), GetNDSpace());
477 0 : A->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
478 : return A;
479 : }
480 : else
481 : {
482 0 : MFEM_VERIFY(!ai, "Unexpected imaginary part in GetExtraSystemMatrix<Operator>!");
483 0 : auto A = std::make_unique<ParOperator>(std::move(ar), GetNDSpace());
484 0 : A->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy);
485 : return A;
486 : }
487 0 : }
488 :
489 : template <typename OperType, typename ScalarType>
490 : std::unique_ptr<OperType>
491 0 : SpaceOperator::GetSystemMatrix(ScalarType a0, ScalarType a1, ScalarType a2,
492 : const OperType *K, const OperType *C, const OperType *M,
493 : const OperType *A2)
494 : {
495 0 : return BuildParSumOperator({a0, a1, a2, ScalarType{1}}, {K, C, M, A2});
496 : }
497 :
498 0 : std::unique_ptr<Operator> SpaceOperator::GetInnerProductMatrix(double a0, double a2,
499 : const ComplexOperator *K,
500 : const ComplexOperator *M)
501 : {
502 0 : const auto *PtAP_K = (K) ? dynamic_cast<const ComplexParOperator *>(K) : nullptr;
503 0 : const auto *PtAP_M = (M) ? dynamic_cast<const ComplexParOperator *>(M) : nullptr;
504 0 : return BuildParSumOperator(
505 0 : {a0, a2}, {PtAP_K ? PtAP_K->Real() : nullptr, PtAP_M ? PtAP_M->Real() : nullptr});
506 : }
507 :
508 : namespace
509 : {
510 :
511 : template <typename OperType>
512 : auto BuildLevelParOperator(std::unique_ptr<Operator> &&br, std::unique_ptr<Operator> &&bi,
513 : const FiniteElementSpace &fespace);
514 :
515 : template <>
516 0 : auto BuildLevelParOperator<Operator>(std::unique_ptr<Operator> &&br,
517 : std::unique_ptr<Operator> &&bi,
518 : const FiniteElementSpace &fespace)
519 : {
520 0 : MFEM_VERIFY(
521 : !bi,
522 : "Should not be constructing a real-valued ParOperator with non-zero imaginary part!");
523 0 : return std::make_unique<ParOperator>(std::move(br), fespace);
524 : }
525 :
526 : template <>
527 : auto BuildLevelParOperator<ComplexOperator>(std::unique_ptr<Operator> &&br,
528 : std::unique_ptr<Operator> &&bi,
529 : const FiniteElementSpace &fespace)
530 : {
531 0 : return std::make_unique<ComplexParOperator>(std::move(br), std::move(bi), fespace);
532 : }
533 :
534 : } // namespace
535 :
536 0 : void SpaceOperator::AssemblePreconditioner(
537 : std::complex<double> a0, std::complex<double> a1, std::complex<double> a2, double a3,
538 : std::vector<std::unique_ptr<Operator>> &br_vec,
539 : std::vector<std::unique_ptr<Operator>> &br_aux_vec,
540 : std::vector<std::unique_ptr<Operator>> &bi_vec,
541 : std::vector<std::unique_ptr<Operator>> &bi_aux_vec)
542 : {
543 : constexpr bool skip_zeros = false, assemble_q_data = false;
544 0 : MaterialPropertyCoefficient dfr(mat_op.MaxCeedAttribute()),
545 0 : dfi(mat_op.MaxCeedAttribute()), fr(mat_op.MaxCeedAttribute()),
546 0 : fi(mat_op.MaxCeedAttribute()), dfbr(mat_op.MaxCeedBdrAttribute()),
547 0 : dfbi(mat_op.MaxCeedBdrAttribute()), fbr(mat_op.MaxCeedBdrAttribute()),
548 0 : fbi(mat_op.MaxCeedBdrAttribute()), fpi(mat_op.MaxCeedAttribute()),
549 0 : fpr(mat_op.MaxCeedAttribute());
550 0 : AddStiffnessCoefficients(a0.real(), dfr, fr);
551 0 : AddStiffnessCoefficients(a0.imag(), dfi, fi);
552 0 : AddStiffnessBdrCoefficients(a0.real(), fbr);
553 0 : AddStiffnessBdrCoefficients(a0.imag(), fbi);
554 0 : AddDampingCoefficients(a1.real(), fr);
555 0 : AddDampingCoefficients(a1.imag(), fi);
556 0 : AddDampingBdrCoefficients(a1.real(), fbr);
557 0 : AddDampingBdrCoefficients(a1.imag(), fbi);
558 0 : AddRealMassCoefficients(pc_mat_shifted ? std::abs(a2.real()) : a2.real(), fr);
559 0 : AddRealMassCoefficients(a2.imag(), fi);
560 0 : AddRealMassBdrCoefficients(pc_mat_shifted ? std::abs(a2.real()) : a2.real(), fbr);
561 0 : AddRealMassBdrCoefficients(a2.imag(), fbi);
562 0 : AddImagMassCoefficients(a2.real(), fi);
563 0 : AddImagMassCoefficients(-a2.imag(), fr);
564 0 : AddExtraSystemBdrCoefficients(a3, dfbr, dfbi, fbr, fbi);
565 0 : AddRealPeriodicCoefficients(a0.real(), fr);
566 0 : AddRealPeriodicCoefficients(a0.imag(), fi);
567 0 : AddImagPeriodicCoefficients(a0.real(), fpi);
568 0 : AddImagPeriodicCoefficients(-a0.imag(), fpr);
569 : int empty[2] = {
570 0 : (dfr.empty() && fr.empty() && dfbr.empty() && fbr.empty() && fpr.empty()),
571 0 : (dfi.empty() && fi.empty() && dfbi.empty() && fbi.empty() && fpi.empty())};
572 : Mpi::GlobalMin(2, empty, GetComm());
573 0 : if (!empty[0])
574 : {
575 0 : br_vec = AssembleOperators(GetNDSpaces(), &dfr, &fr, &dfbr, &fbr, &fpr, skip_zeros,
576 0 : assemble_q_data);
577 : br_aux_vec =
578 0 : AssembleAuxOperators(GetH1Spaces(), &fr, &fbr, skip_zeros, assemble_q_data);
579 : }
580 0 : if (!empty[1])
581 : {
582 0 : bi_vec = AssembleOperators(GetNDSpaces(), &dfi, &fi, &dfbi, &fbi, &fpi, skip_zeros,
583 0 : assemble_q_data);
584 : bi_aux_vec =
585 0 : AssembleAuxOperators(GetH1Spaces(), &fi, &fbi, skip_zeros, assemble_q_data);
586 : }
587 0 : }
588 :
589 0 : void SpaceOperator::AssemblePreconditioner(
590 : std::complex<double> a0, std::complex<double> a1, std::complex<double> a2, double a3,
591 : std::vector<std::unique_ptr<Operator>> &br_vec,
592 : std::vector<std::unique_ptr<Operator>> &br_aux_vec)
593 : {
594 : constexpr bool skip_zeros = false, assemble_q_data = false;
595 0 : MaterialPropertyCoefficient dfr(mat_op.MaxCeedAttribute()), fr(mat_op.MaxCeedAttribute()),
596 0 : dfbr(mat_op.MaxCeedBdrAttribute()), fbr(mat_op.MaxCeedBdrAttribute());
597 0 : AddStiffnessCoefficients(a0.real(), dfr, fr);
598 0 : AddStiffnessBdrCoefficients(a0.real(), fbr);
599 0 : AddDampingCoefficients(a1.imag(), fr);
600 0 : AddDampingBdrCoefficients(a1.imag(), fbr);
601 0 : AddAbsMassCoefficients(pc_mat_shifted ? std::abs(a2.real()) : a2.real(), fr);
602 0 : AddRealMassBdrCoefficients(pc_mat_shifted ? std::abs(a2.real()) : a2.real(), fbr);
603 0 : AddExtraSystemBdrCoefficients(a3, dfbr, dfbr, fbr, fbr);
604 0 : AddRealPeriodicCoefficients(a0.real(), fr);
605 0 : int empty = (dfr.empty() && fr.empty() && dfbr.empty() && fbr.empty());
606 : Mpi::GlobalMin(1, &empty, GetComm());
607 0 : if (!empty)
608 : {
609 0 : br_vec = AssembleOperators(GetNDSpaces(), &dfr, &fr, &dfbr, &fbr, nullptr, skip_zeros,
610 0 : assemble_q_data);
611 : br_aux_vec =
612 0 : AssembleAuxOperators(GetH1Spaces(), &fr, &fbr, skip_zeros, assemble_q_data);
613 : }
614 0 : }
615 :
616 0 : void SpaceOperator::AssemblePreconditioner(
617 : double a0, double a1, double a2, double a3,
618 : std::vector<std::unique_ptr<Operator>> &br_vec,
619 : std::vector<std::unique_ptr<Operator>> &br_aux_vec)
620 : {
621 : constexpr bool skip_zeros = false, assemble_q_data = false;
622 0 : MaterialPropertyCoefficient dfr(mat_op.MaxCeedAttribute()), fr(mat_op.MaxCeedAttribute()),
623 0 : dfbr(mat_op.MaxCeedBdrAttribute()), fbr(mat_op.MaxCeedBdrAttribute());
624 0 : AddStiffnessCoefficients(a0, dfr, fr);
625 0 : AddStiffnessBdrCoefficients(a0, fbr);
626 0 : AddDampingCoefficients(a1, fr);
627 0 : AddDampingBdrCoefficients(a1, fbr);
628 0 : AddAbsMassCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fr);
629 0 : AddRealMassBdrCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fbr);
630 0 : AddExtraSystemBdrCoefficients(a3, dfbr, dfbr, fbr, fbr);
631 0 : AddRealPeriodicCoefficients(a0, fr);
632 0 : int empty = (dfr.empty() && fr.empty() && dfbr.empty() && fbr.empty());
633 : Mpi::GlobalMin(1, &empty, GetComm());
634 0 : if (!empty)
635 : {
636 0 : br_vec = AssembleOperators(GetNDSpaces(), &dfr, &fr, &dfbr, &fbr, nullptr, skip_zeros,
637 0 : assemble_q_data);
638 : br_aux_vec =
639 0 : AssembleAuxOperators(GetH1Spaces(), &fr, &fbr, skip_zeros, assemble_q_data);
640 : }
641 0 : }
642 :
643 : template <typename OperType, typename ScalarType>
644 0 : std::unique_ptr<OperType> SpaceOperator::GetPreconditionerMatrix(ScalarType a0,
645 : ScalarType a1,
646 : ScalarType a2, double a3)
647 : {
648 : // When partially assembled, the coarse operators can reuse the fine operator quadrature
649 : // data if the spaces correspond to the same mesh. When appropriate, we build the
650 : // preconditioner on all levels based on the actual complex-valued system matrix. The
651 : // coarse operator is always fully assembled.
652 0 : if (print_prec_hdr)
653 : {
654 0 : Mpi::Print("\nAssembling multigrid hierarchy:\n");
655 : }
656 0 : MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(),
657 : "Multigrid hierarchy mismatch for auxiliary space preconditioning!");
658 :
659 0 : const auto n_levels = GetNDSpaces().GetNumLevels();
660 0 : std::vector<std::unique_ptr<Operator>> br_vec(n_levels), bi_vec(n_levels),
661 0 : br_aux_vec(n_levels), bi_aux_vec(n_levels);
662 0 : if (std::is_same<OperType, ComplexOperator>::value && !pc_mat_real)
663 : {
664 0 : AssemblePreconditioner(a0, a1, a2, a3, br_vec, br_aux_vec, bi_vec, bi_aux_vec);
665 : }
666 : else
667 : {
668 0 : AssemblePreconditioner(a0, a1, a2, a3, br_vec, br_aux_vec);
669 : }
670 :
671 0 : auto B = std::make_unique<BaseMultigridOperator<OperType>>(n_levels);
672 0 : for (bool aux : {false, true})
673 : {
674 0 : for (std::size_t l = 0; l < n_levels; l++)
675 : {
676 0 : const auto &fespace_l =
677 : aux ? GetH1Spaces().GetFESpaceAtLevel(l) : GetNDSpaces().GetFESpaceAtLevel(l);
678 0 : const auto &dbc_tdof_lists_l = aux ? h1_dbc_tdof_lists[l] : nd_dbc_tdof_lists[l];
679 0 : auto &br_l = aux ? br_aux_vec[l] : br_vec[l];
680 0 : auto &bi_l = aux ? bi_aux_vec[l] : bi_vec[l];
681 0 : if (print_prec_hdr)
682 : {
683 0 : Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, aux ? " (auxiliary)" : "",
684 0 : fespace_l.GetMaxElementOrder(), fespace_l.GlobalTrueVSize());
685 0 : const auto *b_spm = dynamic_cast<const hypre::HypreCSRMatrix *>(br_l.get());
686 0 : if (!b_spm)
687 : {
688 0 : b_spm = dynamic_cast<const hypre::HypreCSRMatrix *>(bi_l.get());
689 : }
690 0 : if (b_spm)
691 : {
692 0 : HYPRE_BigInt nnz = b_spm->NNZ();
693 : Mpi::GlobalSum(1, &nnz, fespace_l.GetComm());
694 0 : Mpi::Print(", {:d} NNZ\n", nnz);
695 : }
696 : else
697 : {
698 0 : Mpi::Print("\n");
699 : }
700 : }
701 0 : auto B_l =
702 : BuildLevelParOperator<OperType>(std::move(br_l), std::move(bi_l), fespace_l);
703 0 : B_l->SetEssentialTrueDofs(dbc_tdof_lists_l, Operator::DiagonalPolicy::DIAG_ONE);
704 0 : if (aux)
705 : {
706 0 : B->AddAuxiliaryOperator(std::move(B_l));
707 : }
708 : else
709 : {
710 0 : B->AddOperator(std::move(B_l));
711 : }
712 : }
713 : }
714 :
715 0 : print_prec_hdr = false;
716 0 : return B;
717 0 : }
718 :
719 0 : void SpaceOperator::AddStiffnessCoefficients(double coeff, MaterialPropertyCoefficient &df,
720 : MaterialPropertyCoefficient &f)
721 : {
722 : // Contribution from material permeability.
723 0 : df.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetInvPermeability(), coeff);
724 :
725 : // Contribution for London superconductors.
726 0 : if (mat_op.HasLondonDepth())
727 : {
728 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetInvLondonDepth(), coeff);
729 : }
730 0 : }
731 :
732 0 : void SpaceOperator::AddStiffnessBdrCoefficients(double coeff,
733 : MaterialPropertyCoefficient &fb)
734 : {
735 : // Robin BC contributions due to surface impedance and lumped ports (inductance).
736 0 : surf_z_op.AddStiffnessBdrCoefficients(coeff, fb);
737 0 : lumped_port_op.AddStiffnessBdrCoefficients(coeff, fb);
738 0 : }
739 :
740 0 : void SpaceOperator::AddDampingCoefficients(double coeff, MaterialPropertyCoefficient &f)
741 : {
742 : // Contribution for domain conductivity.
743 0 : if (mat_op.HasConductivity())
744 : {
745 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetConductivity(), coeff);
746 : }
747 0 : }
748 :
749 0 : void SpaceOperator::AddDampingBdrCoefficients(double coeff, MaterialPropertyCoefficient &fb)
750 : {
751 : // Robin BC contributions due to surface impedance, lumped ports, and absorbing
752 : // boundaries (resistance).
753 0 : farfield_op.AddDampingBdrCoefficients(coeff, fb);
754 0 : surf_z_op.AddDampingBdrCoefficients(coeff, fb);
755 0 : lumped_port_op.AddDampingBdrCoefficients(coeff, fb);
756 0 : }
757 :
758 0 : void SpaceOperator::AddRealMassCoefficients(double coeff, MaterialPropertyCoefficient &f)
759 : {
760 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetPermittivityReal(), coeff);
761 0 : }
762 :
763 0 : void SpaceOperator::AddRealMassBdrCoefficients(double coeff,
764 : MaterialPropertyCoefficient &fb)
765 : {
766 : // Robin BC contributions due to surface impedance and lumped ports (capacitance).
767 0 : surf_z_op.AddMassBdrCoefficients(coeff, fb);
768 0 : lumped_port_op.AddMassBdrCoefficients(coeff, fb);
769 0 : }
770 :
771 0 : void SpaceOperator::AddImagMassCoefficients(double coeff, MaterialPropertyCoefficient &f)
772 : {
773 : // Contribution for loss tangent: ε -> ε * (1 - i tan(δ)).
774 0 : if (mat_op.HasLossTangent())
775 : {
776 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetPermittivityImag(), coeff);
777 : }
778 0 : }
779 :
780 0 : void SpaceOperator::AddAbsMassCoefficients(double coeff, MaterialPropertyCoefficient &f)
781 : {
782 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetPermittivityAbs(), coeff);
783 0 : }
784 :
785 0 : void SpaceOperator::AddExtraSystemBdrCoefficients(double omega,
786 : MaterialPropertyCoefficient &dfbr,
787 : MaterialPropertyCoefficient &dfbi,
788 : MaterialPropertyCoefficient &fbr,
789 : MaterialPropertyCoefficient &fbi)
790 : {
791 : // Contribution for second-order farfield boundaries and finite conductivity boundaries.
792 0 : farfield_op.AddExtraSystemBdrCoefficients(omega, dfbr, dfbi);
793 0 : surf_sigma_op.AddExtraSystemBdrCoefficients(omega, fbr, fbi);
794 :
795 : // Contribution for numeric wave ports.
796 0 : wave_port_op.AddExtraSystemBdrCoefficients(omega, fbr, fbi);
797 0 : }
798 :
799 0 : void SpaceOperator::AddRealPeriodicCoefficients(double coeff,
800 : MaterialPropertyCoefficient &f)
801 : {
802 : // Floquet periodicity contributions.
803 0 : if (mat_op.HasWaveVector())
804 : {
805 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetFloquetMass(), coeff);
806 : }
807 0 : }
808 :
809 0 : void SpaceOperator::AddImagPeriodicCoefficients(double coeff,
810 : MaterialPropertyCoefficient &f)
811 : {
812 : // Floquet periodicity contributions.
813 0 : if (mat_op.HasWaveVector())
814 : {
815 0 : f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetFloquetCurl(), coeff);
816 : }
817 0 : }
818 :
819 0 : bool SpaceOperator::GetExcitationVector(int excitation_idx, Vector &RHS)
820 : {
821 : // Time domain excitation vector.
822 0 : RHS.SetSize(GetNDSpace().GetTrueVSize());
823 0 : RHS.UseDevice(true);
824 0 : RHS = 0.0;
825 0 : bool nnz = AddExcitationVector1Internal(excitation_idx, RHS);
826 0 : linalg::SetSubVector(RHS, nd_dbc_tdof_lists.back(), 0.0);
827 0 : return nnz;
828 : }
829 :
830 0 : bool SpaceOperator::GetExcitationVector(int excitation_idx, double omega,
831 : ComplexVector &RHS)
832 : {
833 : // Frequency domain excitation vector: RHS = iω RHS1 + RHS2(ω).
834 0 : RHS.SetSize(GetNDSpace().GetTrueVSize());
835 0 : RHS.UseDevice(true);
836 : RHS = 0.0;
837 0 : bool nnz1 = AddExcitationVector1Internal(excitation_idx, RHS.Real());
838 0 : RHS *= 1i * omega;
839 0 : bool nnz2 = AddExcitationVector2Internal(excitation_idx, omega, RHS);
840 0 : linalg::SetSubVector(RHS, nd_dbc_tdof_lists.back(), 0.0);
841 0 : return nnz1 || nnz2;
842 : }
843 :
844 0 : bool SpaceOperator::GetExcitationVector1(int excitation_idx, ComplexVector &RHS1)
845 : {
846 : // Assemble the frequency domain excitation term with linear frequency dependence
847 : // (coefficient iω, see GetExcitationVector above, is accounted for later).
848 0 : RHS1.SetSize(GetNDSpace().GetTrueVSize());
849 0 : RHS1.UseDevice(true);
850 : RHS1 = 0.0;
851 0 : bool nnz1 = AddExcitationVector1Internal(excitation_idx, RHS1.Real());
852 0 : linalg::SetSubVector(RHS1.Real(), nd_dbc_tdof_lists.back(), 0.0);
853 0 : return nnz1;
854 : }
855 :
856 0 : bool SpaceOperator::GetExcitationVector2(int excitation_idx, double omega,
857 : ComplexVector &RHS2)
858 : {
859 0 : RHS2.SetSize(GetNDSpace().GetTrueVSize());
860 0 : RHS2.UseDevice(true);
861 : RHS2 = 0.0;
862 0 : bool nnz2 = AddExcitationVector2Internal(excitation_idx, omega, RHS2);
863 0 : linalg::SetSubVector(RHS2, nd_dbc_tdof_lists.back(), 0.0);
864 0 : return nnz2;
865 : }
866 :
867 0 : bool SpaceOperator::AddExcitationVector1Internal(int excitation_idx, Vector &RHS1)
868 : {
869 : // Assemble the time domain excitation -g'(t) J or frequency domain excitation -iω J.
870 : // The g'(t) or iω factors are not accounted for here, they is accounted for in the time
871 : // integration or frequency sweep later.
872 0 : MFEM_VERIFY(RHS1.Size() == GetNDSpace().GetTrueVSize(),
873 : "Invalid T-vector size for AddExcitationVector1Internal!");
874 : SumVectorCoefficient fb(GetMesh().SpaceDimension());
875 0 : lumped_port_op.AddExcitationBdrCoefficients(excitation_idx, fb);
876 0 : surf_j_op.AddExcitationBdrCoefficients(fb); // No excitation_idx — currently in all
877 0 : int empty = (fb.empty());
878 : Mpi::GlobalMin(1, &empty, GetComm());
879 0 : if (empty)
880 : {
881 : return false;
882 : }
883 0 : mfem::LinearForm rhs1(&GetNDSpace().Get());
884 0 : rhs1.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb));
885 0 : rhs1.UseFastAssembly(false);
886 : rhs1.UseDevice(false);
887 0 : rhs1.Assemble();
888 : rhs1.UseDevice(true);
889 0 : GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs1, RHS1);
890 : return true;
891 0 : }
892 :
893 0 : bool SpaceOperator::AddExcitationVector2Internal(int excitation_idx, double omega,
894 : ComplexVector &RHS2)
895 : {
896 : // Assemble the contribution of wave ports to the frequency domain excitation term at the
897 : // specified frequency.
898 0 : MFEM_VERIFY(RHS2.Size() == GetNDSpace().GetTrueVSize(),
899 : "Invalid T-vector size for AddExcitationVector2Internal!");
900 : SumVectorCoefficient fbr(GetMesh().SpaceDimension()), fbi(GetMesh().SpaceDimension());
901 0 : wave_port_op.AddExcitationBdrCoefficients(excitation_idx, omega, fbr, fbi);
902 0 : int empty = (fbr.empty() && fbi.empty());
903 : Mpi::GlobalMin(1, &empty, GetComm());
904 0 : if (empty)
905 : {
906 : return false;
907 : }
908 : {
909 0 : mfem::LinearForm rhs2(&GetNDSpace().Get());
910 0 : rhs2.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbr));
911 0 : rhs2.UseFastAssembly(false);
912 : rhs2.UseDevice(false);
913 0 : rhs2.Assemble();
914 : rhs2.UseDevice(true);
915 0 : GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs2, RHS2.Real());
916 0 : }
917 : {
918 0 : mfem::LinearForm rhs2(&GetNDSpace().Get());
919 0 : rhs2.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi));
920 0 : rhs2.UseFastAssembly(false);
921 : rhs2.UseDevice(false);
922 0 : rhs2.Assemble();
923 : rhs2.UseDevice(true);
924 0 : GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs2, RHS2.Imag());
925 0 : }
926 0 : return true;
927 : }
928 :
929 0 : void SpaceOperator::GetConstantInitialVector(ComplexVector &v)
930 : {
931 0 : v.SetSize(GetNDSpace().GetTrueVSize());
932 0 : v.UseDevice(true);
933 : v = 1.0;
934 0 : linalg::SetSubVector(v.Real(), nd_dbc_tdof_lists.back(), 0.0);
935 0 : }
936 :
937 0 : void SpaceOperator::GetRandomInitialVector(ComplexVector &v)
938 : {
939 0 : v.SetSize(GetNDSpace().GetTrueVSize());
940 0 : v.UseDevice(true);
941 0 : linalg::SetRandom(GetNDSpace().GetComm(), v);
942 0 : linalg::SetSubVector(v, nd_dbc_tdof_lists.back(), 0.0);
943 0 : }
944 :
945 : template std::unique_ptr<Operator>
946 : SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy);
947 : template std::unique_ptr<ComplexOperator>
948 : SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy);
949 :
950 : template std::unique_ptr<Operator>
951 : SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy);
952 : template std::unique_ptr<ComplexOperator>
953 : SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy);
954 :
955 : template std::unique_ptr<Operator> SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy);
956 : template std::unique_ptr<ComplexOperator>
957 : SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy);
958 :
959 : template std::unique_ptr<Operator>
960 : SpaceOperator::GetExtraSystemMatrix(double, Operator::DiagonalPolicy);
961 : template std::unique_ptr<ComplexOperator>
962 : SpaceOperator::GetExtraSystemMatrix(double, Operator::DiagonalPolicy);
963 :
964 : template std::unique_ptr<Operator>
965 : SpaceOperator::GetSystemMatrix<Operator, double>(double, double, double, const Operator *,
966 : const Operator *, const Operator *,
967 : const Operator *);
968 : template std::unique_ptr<ComplexOperator>
969 : SpaceOperator::GetSystemMatrix<ComplexOperator, std::complex<double>>(
970 : std::complex<double>, std::complex<double>, std::complex<double>,
971 : const ComplexOperator *, const ComplexOperator *, const ComplexOperator *,
972 : const ComplexOperator *);
973 :
974 : template std::unique_ptr<Operator>
975 : SpaceOperator::GetPreconditionerMatrix<Operator, double>(double, double, double, double);
976 : template std::unique_ptr<ComplexOperator>
977 : SpaceOperator::GetPreconditionerMatrix<ComplexOperator, std::complex<double>>(
978 : std::complex<double>, std::complex<double>, std::complex<double>, double);
979 :
980 : } // namespace palace
|
