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 "waveportoperator.hpp"
5 :
6 : #include <tuple>
7 : #include <fmt/ranges.h>
8 : #include "fem/bilinearform.hpp"
9 : #include "fem/coefficient.hpp"
10 : #include "fem/integrator.hpp"
11 : #include "linalg/arpack.hpp"
12 : #include "linalg/iterative.hpp"
13 : #include "linalg/mumps.hpp"
14 : #include "linalg/rap.hpp"
15 : #include "linalg/slepc.hpp"
16 : #include "linalg/solver.hpp"
17 : #include "linalg/strumpack.hpp"
18 : #include "linalg/superlu.hpp"
19 : #include "models/materialoperator.hpp"
20 : #include "utils/communication.hpp"
21 : #include "utils/geodata.hpp"
22 : #include "utils/iodata.hpp"
23 : #include "utils/timer.hpp"
24 :
25 : namespace palace
26 : {
27 :
28 : using namespace std::complex_literals;
29 :
30 : namespace
31 : {
32 :
33 0 : void GetEssentialTrueDofs(mfem::ParGridFunction &E0t, mfem::ParGridFunction &E0n,
34 : mfem::ParGridFunction &port_E0t, mfem::ParGridFunction &port_E0n,
35 : mfem::ParTransferMap &port_nd_transfer,
36 : mfem::ParTransferMap &port_h1_transfer,
37 : const mfem::Array<int> &dbc_attr,
38 : mfem::Array<int> &port_nd_dbc_tdof_list,
39 : mfem::Array<int> &port_h1_dbc_tdof_list)
40 : {
41 : auto &nd_fespace = *E0t.ParFESpace();
42 : auto &h1_fespace = *E0n.ParFESpace();
43 : auto &port_nd_fespace = *port_E0t.ParFESpace();
44 : auto &port_h1_fespace = *port_E0n.ParFESpace();
45 : const auto &mesh = *nd_fespace.GetParMesh();
46 :
47 : mfem::Array<int> dbc_marker, nd_dbc_tdof_list, h1_dbc_tdof_list;
48 0 : mesh::AttrToMarker(mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0, dbc_attr,
49 : dbc_marker);
50 0 : nd_fespace.GetEssentialTrueDofs(dbc_marker, nd_dbc_tdof_list);
51 0 : h1_fespace.GetEssentialTrueDofs(dbc_marker, h1_dbc_tdof_list);
52 :
53 0 : Vector tE0t(nd_fespace.GetTrueVSize()), tE0n(h1_fespace.GetTrueVSize());
54 : tE0t.UseDevice(true);
55 : tE0n.UseDevice(true);
56 0 : tE0t = 0.0;
57 0 : tE0n = 0.0;
58 0 : linalg::SetSubVector(tE0t, nd_dbc_tdof_list, 1.0);
59 0 : linalg::SetSubVector(tE0n, h1_dbc_tdof_list, 1.0);
60 0 : E0t.SetFromTrueDofs(tE0t);
61 0 : E0n.SetFromTrueDofs(tE0n);
62 0 : port_nd_transfer.Transfer(E0t, port_E0t);
63 0 : port_h1_transfer.Transfer(E0n, port_E0n);
64 :
65 0 : Vector port_tE0t(port_nd_fespace.GetTrueVSize()),
66 0 : port_tE0n(port_h1_fespace.GetTrueVSize());
67 : port_tE0t.UseDevice(true);
68 : port_tE0n.UseDevice(true);
69 0 : port_E0t.ParallelProject(port_tE0t);
70 0 : port_E0n.ParallelProject(port_tE0n);
71 : {
72 : const auto *h_port_tE0t = port_tE0t.HostRead();
73 : const auto *h_port_tE0n = port_tE0n.HostRead();
74 0 : for (int i = 0; i < port_tE0t.Size(); i++)
75 : {
76 0 : if (h_port_tE0t[i] != 0.0)
77 : {
78 0 : port_nd_dbc_tdof_list.Append(i);
79 : }
80 : }
81 0 : for (int i = 0; i < port_tE0n.Size(); i++)
82 : {
83 0 : if (h_port_tE0n[i] != 0.0)
84 : {
85 0 : port_h1_dbc_tdof_list.Append(i);
86 : }
87 : }
88 : }
89 0 : }
90 :
91 0 : void GetInitialSpace(const mfem::ParFiniteElementSpace &nd_fespace,
92 : const mfem::ParFiniteElementSpace &h1_fespace,
93 : const mfem::Array<int> &dbc_tdof_list, ComplexVector &v)
94 : {
95 : // Initial space which satisfies Dirichlet BCs.
96 0 : const int nd_size = nd_fespace.GetTrueVSize(), h1_size = h1_fespace.GetTrueVSize();
97 0 : v.SetSize(nd_size + h1_size);
98 0 : v.UseDevice(true);
99 0 : v = std::complex<double>(1.0, 0.0);
100 : // linalg::SetRandomReal(nd_fespace.GetComm(), v);
101 0 : linalg::SetSubVector(v, nd_size, nd_size + h1_size, 0.0);
102 0 : linalg::SetSubVector(v, dbc_tdof_list, 0.0);
103 0 : }
104 :
105 : using ComplexHypreParMatrix = std::tuple<std::unique_ptr<mfem::HypreParMatrix>,
106 : std::unique_ptr<mfem::HypreParMatrix>>;
107 : constexpr bool skip_zeros = false;
108 :
109 0 : ComplexHypreParMatrix GetAtt(const MaterialOperator &mat_op,
110 : const FiniteElementSpace &nd_fespace,
111 : const mfem::Vector &normal, double omega, double sigma)
112 : {
113 : // Stiffness matrix (shifted): Aₜₜ = (μ⁻¹ ∇ₜ x u, ∇ₜ x v) - ω² (ε u, v) - σ (μ⁻¹ u, v).
114 0 : MaterialPropertyCoefficient muinv_func(mat_op.GetBdrAttributeToMaterial(),
115 0 : mat_op.GetInvPermeability());
116 0 : muinv_func.NormalProjectedCoefficient(normal);
117 0 : MaterialPropertyCoefficient epsilon_func(mat_op.GetBdrAttributeToMaterial(),
118 0 : mat_op.GetPermittivityReal(), -omega * omega);
119 0 : epsilon_func.AddCoefficient(mat_op.GetBdrAttributeToMaterial(),
120 : mat_op.GetInvPermeability(), -sigma);
121 : BilinearForm attr(nd_fespace);
122 0 : attr.AddDomainIntegrator<CurlCurlMassIntegrator>(muinv_func, epsilon_func);
123 :
124 : // Contribution for loss tangent: ε -> ε * (1 - i tan(δ)).
125 0 : if (!mat_op.HasLossTangent())
126 : {
127 0 : return {ParOperator(attr.FullAssemble(skip_zeros), nd_fespace).StealParallelAssemble(),
128 : nullptr};
129 : }
130 : MaterialPropertyCoefficient negepstandelta_func(
131 0 : mat_op.GetBdrAttributeToMaterial(), mat_op.GetPermittivityImag(), -omega * omega);
132 : BilinearForm atti(nd_fespace);
133 0 : atti.AddDomainIntegrator<VectorFEMassIntegrator>(negepstandelta_func);
134 0 : return {ParOperator(attr.FullAssemble(skip_zeros), nd_fespace).StealParallelAssemble(),
135 0 : ParOperator(atti.FullAssemble(skip_zeros), nd_fespace).StealParallelAssemble()};
136 0 : }
137 :
138 0 : ComplexHypreParMatrix GetAtn(const MaterialOperator &mat_op,
139 : const FiniteElementSpace &nd_fespace,
140 : const FiniteElementSpace &h1_fespace)
141 : {
142 : // Coupling matrix: Aₜₙ = -(μ⁻¹ ∇ₜ u, v).
143 0 : MaterialPropertyCoefficient muinv_func(mat_op.GetBdrAttributeToMaterial(),
144 0 : mat_op.GetInvPermeability(), -1.0);
145 : BilinearForm atn(h1_fespace, nd_fespace);
146 0 : atn.AddDomainIntegrator<MixedVectorGradientIntegrator>(muinv_func);
147 0 : return {ParOperator(atn.FullAssemble(skip_zeros), h1_fespace, nd_fespace, false)
148 : .StealParallelAssemble(),
149 0 : nullptr};
150 0 : }
151 :
152 0 : ComplexHypreParMatrix GetAnt(const MaterialOperator &mat_op,
153 : const FiniteElementSpace &h1_fespace,
154 : const FiniteElementSpace &nd_fespace)
155 : {
156 : // Coupling matrix: Aₙₜ = -(ε u, ∇ₜ v).
157 0 : MaterialPropertyCoefficient epsilon_func(mat_op.GetBdrAttributeToMaterial(),
158 0 : mat_op.GetPermittivityReal(), 1.0);
159 :
160 : BilinearForm antr(nd_fespace, h1_fespace);
161 0 : antr.AddDomainIntegrator<MixedVectorWeakDivergenceIntegrator>(epsilon_func);
162 :
163 : // Contribution for loss tangent: ε -> ε * (1 - i tan(δ)).
164 0 : if (!mat_op.HasLossTangent())
165 : {
166 0 : return {ParOperator(antr.FullAssemble(skip_zeros), nd_fespace, h1_fespace, false)
167 : .StealParallelAssemble(),
168 : nullptr};
169 : }
170 0 : MaterialPropertyCoefficient negepstandelta_func(mat_op.GetBdrAttributeToMaterial(),
171 0 : mat_op.GetPermittivityImag(), 1.0);
172 : BilinearForm anti(nd_fespace, h1_fespace);
173 0 : anti.AddDomainIntegrator<MixedVectorWeakDivergenceIntegrator>(negepstandelta_func);
174 0 : return {ParOperator(antr.FullAssemble(skip_zeros), nd_fespace, h1_fespace, false)
175 : .StealParallelAssemble(),
176 0 : ParOperator(anti.FullAssemble(skip_zeros), nd_fespace, h1_fespace, false)
177 : .StealParallelAssemble()};
178 0 : }
179 :
180 0 : ComplexHypreParMatrix GetAnn(const MaterialOperator &mat_op,
181 : const FiniteElementSpace &h1_fespace,
182 : const mfem::Vector &normal)
183 : {
184 : // Mass matrix: Aₙₙ = -(ε u, v).
185 0 : MaterialPropertyCoefficient epsilon_func(mat_op.GetBdrAttributeToMaterial(),
186 0 : mat_op.GetPermittivityReal(), -1.0);
187 0 : epsilon_func.NormalProjectedCoefficient(normal);
188 : BilinearForm annr(h1_fespace);
189 0 : annr.AddDomainIntegrator<MassIntegrator>(epsilon_func);
190 :
191 : // Contribution for loss tangent: ε -> ε * (1 - i tan(δ)).
192 0 : if (!mat_op.HasLossTangent())
193 : {
194 0 : return {ParOperator(annr.FullAssemble(skip_zeros), h1_fespace).StealParallelAssemble(),
195 : nullptr};
196 : }
197 0 : MaterialPropertyCoefficient negepstandelta_func(mat_op.GetBdrAttributeToMaterial(),
198 0 : mat_op.GetPermittivityImag(), -1.0);
199 0 : negepstandelta_func.NormalProjectedCoefficient(normal);
200 : BilinearForm anni(h1_fespace);
201 0 : anni.AddDomainIntegrator<MassIntegrator>(negepstandelta_func);
202 0 : return {ParOperator(annr.FullAssemble(skip_zeros), h1_fespace).StealParallelAssemble(),
203 0 : ParOperator(anni.FullAssemble(skip_zeros), h1_fespace).StealParallelAssemble()};
204 0 : }
205 :
206 0 : ComplexHypreParMatrix GetBtt(const MaterialOperator &mat_op,
207 : const FiniteElementSpace &nd_fespace)
208 : {
209 : // Mass matrix: Bₜₜ = (μ⁻¹ u, v).
210 0 : MaterialPropertyCoefficient muinv_func(mat_op.GetBdrAttributeToMaterial(),
211 0 : mat_op.GetInvPermeability());
212 : BilinearForm btt(nd_fespace);
213 0 : btt.AddDomainIntegrator<VectorFEMassIntegrator>(muinv_func);
214 0 : return {ParOperator(btt.FullAssemble(skip_zeros), nd_fespace).StealParallelAssemble(),
215 0 : nullptr};
216 0 : }
217 :
218 : ComplexHypreParMatrix
219 0 : GetSystemMatrixA(const mfem::HypreParMatrix *Attr, const mfem::HypreParMatrix *Atti,
220 : const mfem::HypreParMatrix *Atnr, const mfem::HypreParMatrix *Atni,
221 : const mfem::HypreParMatrix *Antr, const mfem::HypreParMatrix *Anti,
222 : const mfem::HypreParMatrix *Annr, const mfem::HypreParMatrix *Anni,
223 : const mfem::Array<int> &dbc_tdof_list)
224 : {
225 : // Construct the 2x2 block matrices for the eigenvalue problem A e = λ B e.
226 : mfem::Array2D<const mfem::HypreParMatrix *> blocks(2, 2);
227 0 : blocks(0, 0) = Attr;
228 0 : blocks(0, 1) = Atnr;
229 0 : blocks(1, 0) = Antr;
230 0 : blocks(1, 1) = Annr;
231 0 : std::unique_ptr<mfem::HypreParMatrix> Ar(mfem::HypreParMatrixFromBlocks(blocks));
232 :
233 : std::unique_ptr<mfem::HypreParMatrix> Ai;
234 0 : if (Atti)
235 : {
236 0 : blocks(0, 0) = Atti;
237 0 : blocks(0, 1) = Atni;
238 0 : blocks(1, 0) = Anti;
239 0 : blocks(1, 1) = Anni;
240 0 : Ai.reset(mfem::HypreParMatrixFromBlocks(blocks));
241 : }
242 :
243 : // Eliminate boundary true dofs not associated with this wave port or constrained by
244 : // Dirichlet BCs.
245 0 : Ar->EliminateBC(dbc_tdof_list, Operator::DIAG_ONE);
246 0 : if (Ai)
247 : {
248 0 : Ai->EliminateBC(dbc_tdof_list, Operator::DIAG_ZERO);
249 : }
250 :
251 0 : return {std::move(Ar), std::move(Ai)};
252 : }
253 :
254 0 : ComplexHypreParMatrix GetSystemMatrixB(const mfem::HypreParMatrix *Bttr,
255 : const mfem::HypreParMatrix *Btti,
256 : const mfem::HypreParMatrix *Dnn,
257 : const mfem::Array<int> &dbc_tdof_list)
258 : {
259 : // Construct the 2x2 block matrices for the eigenvalue problem A e = λ B e.
260 : mfem::Array2D<const mfem::HypreParMatrix *> blocks(2, 2);
261 0 : blocks(0, 0) = Bttr;
262 0 : blocks(0, 1) = nullptr;
263 0 : blocks(1, 0) = nullptr;
264 0 : blocks(1, 1) = Dnn;
265 0 : std::unique_ptr<mfem::HypreParMatrix> Br(mfem::HypreParMatrixFromBlocks(blocks));
266 :
267 : std::unique_ptr<mfem::HypreParMatrix> Bi;
268 0 : if (Btti)
269 : {
270 0 : blocks(0, 0) = Btti;
271 0 : Bi.reset(mfem::HypreParMatrixFromBlocks(blocks));
272 : }
273 :
274 : // Eliminate boundary true dofs not associated with this wave port or constrained by
275 : // Dirichlet BCs.
276 0 : Br->EliminateBC(dbc_tdof_list, Operator::DIAG_ZERO);
277 0 : if (Bi)
278 : {
279 0 : Bi->EliminateBC(dbc_tdof_list, Operator::DIAG_ZERO);
280 : }
281 :
282 0 : return {std::move(Br), std::move(Bi)};
283 : }
284 :
285 0 : void Normalize(const GridFunction &S0t, GridFunction &E0t, GridFunction &E0n,
286 : mfem::LinearForm &sr, mfem::LinearForm &si)
287 : {
288 : // Normalize grid functions to a chosen polarization direction and unit power, |E x H⋆| ⋅
289 : // n, integrated over the port surface (+n is the direction of propagation). The n x H
290 : // coefficients are updated implicitly as the only store references to the Et, En grid
291 : // functions. We choose a (rather arbitrary) phase constraint to at least make results for
292 : // the same port consistent between frequencies/meshes.
293 :
294 : // |E x H⋆| ⋅ n = |E ⋅ (-n x H⋆)|. This also updates the n x H coefficients depending on
295 : // Et, En. Update linear forms for postprocessing too.
296 : std::complex<double> dot[2] = {
297 : {sr * S0t.Real(), si * S0t.Real()},
298 0 : {-(sr * E0t.Real()) - (si * E0t.Imag()), -(sr * E0t.Imag()) + (si * E0t.Real())}};
299 : Mpi::GlobalSum(2, dot, S0t.ParFESpace()->GetComm());
300 0 : auto scale = std::abs(dot[0]) / (dot[0] * std::sqrt(std::abs(dot[1])));
301 0 : ComplexVector::AXPBY(scale, E0t.Real(), E0t.Imag(), 0.0, E0t.Real(), E0t.Imag());
302 0 : ComplexVector::AXPBY(scale, E0n.Real(), E0n.Imag(), 0.0, E0n.Real(), E0n.Imag());
303 0 : ComplexVector::AXPBY(scale, sr, si, 0.0, sr, si);
304 :
305 : // This parallel communication is not required since wave port boundaries are true one-
306 : // sided boundaries.
307 : // E0t.Real().ExchangeFaceNbrData(); // Ready for parallel comm on shared faces for n x H
308 : // E0t.Imag().ExchangeFaceNbrData(); // coefficients evaluation
309 : // E0n.Real().ExchangeFaceNbrData();
310 : // E0n.Imag().ExchangeFaceNbrData();
311 0 : }
312 :
313 : // Helper for BdrSubmeshEVectorCoefficient and BdrSubmeshHVectorCoefficient.
314 : enum class ValueType
315 : {
316 : REAL,
317 : IMAG
318 : };
319 :
320 : // Return as a vector coefficient the boundary mode electric field.
321 : template <ValueType Type>
322 0 : class BdrSubmeshEVectorCoefficient : public mfem::VectorCoefficient
323 : {
324 : private:
325 : const GridFunction &Et, &En;
326 : const mfem::ParSubMesh &submesh;
327 : const std::unordered_map<int, int> &submesh_parent_elems;
328 : mfem::IsoparametricTransformation T_loc;
329 :
330 : public:
331 0 : BdrSubmeshEVectorCoefficient(const GridFunction &Et, const GridFunction &En,
332 : const mfem::ParSubMesh &submesh,
333 : const std::unordered_map<int, int> &submesh_parent_elems)
334 0 : : mfem::VectorCoefficient(Et.Real().VectorDim()), Et(Et), En(En), submesh(submesh),
335 0 : submesh_parent_elems(submesh_parent_elems)
336 : {
337 : }
338 :
339 0 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
340 : const mfem::IntegrationPoint &ip) override
341 : {
342 : // Always do the GridFunction evaluation in the submesh.
343 : mfem::ElementTransformation *T_submesh = nullptr;
344 0 : if (T.mesh == submesh.GetParent())
345 : {
346 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
347 : "BdrSubmeshEVectorCoefficient requires ElementType::BDR_ELEMENT when not "
348 : "used on a SubMesh!");
349 0 : auto it = submesh_parent_elems.find(T.ElementNo);
350 0 : if (it == submesh_parent_elems.end())
351 : {
352 : // Just return zero for a parent boundary element not in the submesh.
353 0 : V.SetSize(vdim);
354 0 : V = 0.0;
355 : return;
356 : }
357 : else
358 : {
359 0 : submesh.GetElementTransformation(it->second, &T_loc);
360 : T_loc.SetIntPoint(&ip);
361 0 : T_submesh = &T_loc;
362 : }
363 : }
364 0 : else if (T.mesh == &submesh)
365 : {
366 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::ELEMENT,
367 : "BdrSubmeshEVectorCoefficient requires ElementType::ELEMENT when used on "
368 : "a SubMesh!");
369 : T_submesh = &T;
370 : }
371 : else
372 : {
373 0 : MFEM_ABORT("Invalid mesh for BdrSubmeshEVectorCoefficient!");
374 : }
375 :
376 : // Compute Eₜ + n ⋅ Eₙ . The normal returned by GetNormal points out of the
377 : // computational domain, so we reverse it (direction of propagation is into the domain).
378 : double normal_data[3];
379 0 : mfem::Vector normal(normal_data, vdim);
380 0 : BdrGridFunctionCoefficient::GetNormal(*T_submesh, normal);
381 : if constexpr (Type == ValueType::REAL)
382 : {
383 0 : Et.Real().GetVectorValue(*T_submesh, ip, V);
384 0 : auto Vn = En.Real().GetValue(*T_submesh, ip);
385 0 : V.Add(-Vn, normal);
386 : }
387 : else
388 : {
389 0 : Et.Imag().GetVectorValue(*T_submesh, ip, V);
390 0 : auto Vn = En.Imag().GetValue(*T_submesh, ip);
391 0 : V.Add(-Vn, normal);
392 : }
393 : }
394 : };
395 :
396 : // Computes boundary mode n x H, where +n is the direction of wave propagation: n x H =
397 : // -1/(iωμ) (ikₙ Eₜ + ∇ₜ Eₙ), using the tangential and normal electric field component grid
398 : // functions evaluated on the (single-sided) boundary element.
399 : template <ValueType Type>
400 0 : class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
401 : {
402 : private:
403 : const GridFunction &Et, &En;
404 : const MaterialOperator &mat_op;
405 : const mfem::ParSubMesh &submesh;
406 : const std::unordered_map<int, int> &submesh_parent_elems;
407 : mfem::IsoparametricTransformation T_loc;
408 : std::complex<double> kn;
409 : double omega;
410 :
411 : public:
412 0 : BdrSubmeshHVectorCoefficient(const GridFunction &Et, const GridFunction &En,
413 : const MaterialOperator &mat_op,
414 : const mfem::ParSubMesh &submesh,
415 : const std::unordered_map<int, int> &submesh_parent_elems,
416 : std::complex<double> kn, double omega)
417 0 : : mfem::VectorCoefficient(Et.Real().VectorDim()), Et(Et), En(En), mat_op(mat_op),
418 0 : submesh(submesh), submesh_parent_elems(submesh_parent_elems), kn(kn), omega(omega)
419 : {
420 : }
421 :
422 0 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
423 : const mfem::IntegrationPoint &ip) override
424 : {
425 : // Always do the GridFunction evaluation in the submesh.
426 : mfem::ElementTransformation *T_submesh = nullptr;
427 0 : if (T.mesh == submesh.GetParent())
428 : {
429 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
430 : "BdrSubmeshHVectorCoefficient requires ElementType::BDR_ELEMENT when not "
431 : "used on a SubMesh!");
432 0 : auto it = submesh_parent_elems.find(T.ElementNo);
433 0 : if (it == submesh_parent_elems.end())
434 : {
435 : // Just return zero for a parent boundary element not in the submesh.
436 0 : V.SetSize(vdim);
437 0 : V = 0.0;
438 : return;
439 : }
440 : else
441 : {
442 0 : submesh.GetElementTransformation(it->second, &T_loc);
443 : T_loc.SetIntPoint(&ip);
444 0 : T_submesh = &T_loc;
445 : }
446 : }
447 0 : else if (T.mesh == &submesh)
448 : {
449 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::ELEMENT,
450 : "BdrSubmeshHVectorCoefficient requires ElementType::ELEMENT when used on "
451 : "a SubMesh!");
452 : T_submesh = &T;
453 : }
454 : else
455 : {
456 0 : MFEM_ABORT("Invalid mesh for BdrSubmeshHVectorCoefficient!");
457 : }
458 :
459 : // Get the attribute in the neighboring domain element of the parent mesh.
460 0 : int attr = [&T, this]()
461 : {
462 0 : int i = -1, o, iel1, iel2;
463 0 : if (T.mesh == submesh.GetParent())
464 : {
465 : MFEM_ASSERT(
466 : T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
467 : "BdrSubmeshHVectorCoefficient requires ElementType::BDR_ELEMENT when not "
468 : "used on a SubMesh!");
469 0 : T.mesh->GetBdrElementFace(T.ElementNo, &i, &o);
470 : }
471 0 : else if (T.mesh == &submesh)
472 : {
473 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::ELEMENT,
474 : "BdrSubmeshHVectorCoefficient requires ElementType::ELEMENT when used "
475 : "on a SubMesh!");
476 0 : submesh.GetParent()->GetBdrElementFace(submesh.GetParentElementIDMap()[T.ElementNo],
477 : &i, &o);
478 : }
479 : else
480 : {
481 0 : MFEM_ABORT("Invalid mesh for BdrSubmeshHVectorCoefficient!");
482 : }
483 0 : submesh.GetParent()->GetFaceElements(i, &iel1, &iel2);
484 0 : return submesh.GetParent()->GetAttribute(iel1);
485 0 : }();
486 :
487 : // Compute Re/Im{-1/i (ikₙ Eₜ + ∇ₜ Eₙ)} (t-gradient evaluated in boundary element).
488 : double U_data[3];
489 0 : mfem::Vector U(U_data, vdim);
490 : if constexpr (Type == ValueType::REAL)
491 : {
492 0 : Et.Real().GetVectorValue(*T_submesh, ip, U);
493 0 : U *= -kn.real();
494 :
495 : double dU_data[3];
496 0 : mfem::Vector dU(dU_data, vdim);
497 0 : En.Imag().GetGradient(*T_submesh, dU);
498 0 : U -= dU;
499 : }
500 : else
501 : {
502 0 : Et.Imag().GetVectorValue(*T_submesh, ip, U);
503 0 : U *= -kn.real();
504 :
505 : double dU_data[3];
506 0 : mfem::Vector dU(dU_data, vdim);
507 0 : En.Real().GetGradient(*T_submesh, dU);
508 0 : U += dU;
509 : }
510 :
511 : // Scale by 1/(ωμ) with μ evaluated in the neighboring element.
512 0 : V.SetSize(U.Size());
513 0 : mat_op.GetInvPermeability(attr).Mult(U, V);
514 0 : V *= (1.0 / omega);
515 : }
516 : };
517 :
518 : } // namespace
519 :
520 0 : WavePortData::WavePortData(const config::WavePortData &data,
521 : const config::SolverData &solver, const MaterialOperator &mat_op,
522 : mfem::ParFiniteElementSpace &nd_fespace,
523 : mfem::ParFiniteElementSpace &h1_fespace,
524 0 : const mfem::Array<int> &dbc_attr)
525 0 : : mat_op(mat_op), excitation(data.excitation), active(data.active)
526 : {
527 0 : mode_idx = data.mode_idx;
528 0 : d_offset = data.d_offset;
529 : kn0 = 0.0;
530 0 : omega0 = 0.0;
531 :
532 : // Construct the SubMesh.
533 0 : MFEM_VERIFY(!data.attributes.empty(), "Wave port boundary found with no attributes!");
534 : const auto &mesh = *nd_fespace.GetParMesh();
535 0 : attr_list.Append(data.attributes.data(), data.attributes.size());
536 0 : port_mesh = std::make_unique<Mesh>(std::make_unique<mfem::ParSubMesh>(
537 0 : mfem::ParSubMesh::CreateFromBoundary(mesh, attr_list)));
538 0 : port_normal = mesh::GetSurfaceNormal(*port_mesh);
539 :
540 0 : port_nd_fec = std::make_unique<mfem::ND_FECollection>(nd_fespace.GetMaxElementOrder(),
541 0 : port_mesh->Dimension());
542 0 : port_h1_fec = std::make_unique<mfem::H1_FECollection>(h1_fespace.GetMaxElementOrder(),
543 0 : port_mesh->Dimension());
544 0 : port_nd_fespace = std::make_unique<FiniteElementSpace>(*port_mesh, port_nd_fec.get());
545 0 : port_h1_fespace = std::make_unique<FiniteElementSpace>(*port_mesh, port_h1_fec.get());
546 :
547 0 : GridFunction E0t(nd_fespace), E0n(h1_fespace);
548 0 : port_E0t = std::make_unique<GridFunction>(*port_nd_fespace, true);
549 0 : port_E0n = std::make_unique<GridFunction>(*port_h1_fespace, true);
550 0 : port_E = std::make_unique<GridFunction>(*port_nd_fespace, true);
551 :
552 0 : port_nd_transfer = std::make_unique<mfem::ParTransferMap>(
553 0 : mfem::ParSubMesh::CreateTransferMap(E0t.Real(), port_E0t->Real()));
554 0 : port_h1_transfer = std::make_unique<mfem::ParTransferMap>(
555 0 : mfem::ParSubMesh::CreateTransferMap(E0n.Real(), port_E0n->Real()));
556 :
557 : // Construct mapping from parent (boundary) element indices to submesh (domain)
558 : // elements.
559 : {
560 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
561 : const mfem::Array<int> &parent_elems = port_submesh.GetParentElementIDMap();
562 0 : for (int i = 0; i < parent_elems.Size(); i++)
563 : {
564 0 : submesh_parent_elems[parent_elems[i]] = i;
565 : }
566 : }
567 :
568 : // Extract Dirichlet BC true dofs for the port FE spaces.
569 : {
570 : mfem::Array<int> port_nd_dbc_tdof_list, port_h1_dbc_tdof_list;
571 0 : GetEssentialTrueDofs(E0t.Real(), E0n.Real(), port_E0t->Real(), port_E0n->Real(),
572 : *port_nd_transfer, *port_h1_transfer, dbc_attr,
573 : port_nd_dbc_tdof_list, port_h1_dbc_tdof_list);
574 : int nd_tdof_offset = port_nd_fespace->GetTrueVSize();
575 0 : port_dbc_tdof_list.Reserve(port_nd_dbc_tdof_list.Size() + port_h1_dbc_tdof_list.Size());
576 0 : for (auto tdof : port_nd_dbc_tdof_list)
577 : {
578 0 : port_dbc_tdof_list.Append(tdof);
579 : }
580 0 : for (auto tdof : port_h1_dbc_tdof_list)
581 : {
582 0 : port_dbc_tdof_list.Append(tdof + nd_tdof_offset);
583 : }
584 : }
585 :
586 : // Create vector for initial space for eigenvalue solves and eigenmode solution.
587 0 : GetInitialSpace(*port_nd_fespace, *port_h1_fespace, port_dbc_tdof_list, v0);
588 0 : e0.SetSize(port_nd_fespace->GetTrueVSize() + port_h1_fespace->GetTrueVSize());
589 0 : e0.UseDevice(true);
590 :
591 : // The operators for the generalized eigenvalue problem are:
592 : // [Aₜₜ Aₜₙ] [eₜ] = -kₙ² [Bₜₜ 0ₜₙ] [eₜ]
593 : // [Aₙₜ Aₙₙ] [eₙ] [0ₙₜ 0ₙₙ] [eₙ]
594 : // for the wave port of the given index. The transformed variables are related to the true
595 : // field by Eₜ = eₜ and Eₙ = eₙ / ikₙ. We will actually solve the shift-and-inverse
596 : // problem (A - σ B)⁻¹ B e = λ e, with λ = 1 / (-kₙ² - σ).
597 : // Reference: Vardapetyan and Demkowicz, Full-wave analysis of dielectric waveguides at a
598 : // given frequency, Math. Comput. (2003).
599 : // See also: Halla and Monk, On the analysis of waveguide modes in an electromagnetic
600 : // transmission line, arXiv:2302.11994 (2023).
601 0 : const double c_min = mat_op.GetLightSpeedMax().Min();
602 0 : MFEM_VERIFY(c_min > 0.0 && c_min < mfem::infinity(),
603 : "Invalid material speed of light detected in WavePortOperator!");
604 0 : mu_eps_max = 1.0 / (c_min * c_min) * 1.1; // Add a safety factor for maximum
605 : // propagation constant possible
606 0 : std::tie(Atnr, Atni) = GetAtn(mat_op, *port_nd_fespace, *port_h1_fespace);
607 0 : std::tie(Antr, Anti) = GetAnt(mat_op, *port_h1_fespace, *port_nd_fespace);
608 0 : std::tie(Annr, Anni) = GetAnn(mat_op, *port_h1_fespace, port_normal);
609 : {
610 : // The HypreParMatrix constructor from a SparseMatrix on each process does not copy
611 : // the SparseMatrix data, but that's OK since this Dnn is copied in the block system
612 : // matrix construction.
613 : Vector d(port_h1_fespace->GetTrueVSize());
614 : d.UseDevice(false); // SparseMatrix constructor uses Vector on host
615 0 : d = 0.0;
616 0 : mfem::SparseMatrix diag(d);
617 : auto Dnn = std::make_unique<mfem::HypreParMatrix>(
618 0 : port_h1_fespace->GetComm(), port_h1_fespace->Get().GlobalTrueVSize(),
619 0 : port_h1_fespace->Get().GetTrueDofOffsets(), &diag);
620 0 : auto [Bttr, Btti] = GetBtt(mat_op, *port_nd_fespace);
621 0 : auto [Br, Bi] = GetSystemMatrixB(Bttr.get(), Btti.get(), Dnn.get(), port_dbc_tdof_list);
622 0 : opB = std::make_unique<ComplexWrapperOperator>(std::move(Br), std::move(Bi));
623 0 : }
624 :
625 : // Configure a communicator for the processes which have elements for this port.
626 : MPI_Comm comm = nd_fespace.GetComm();
627 0 : int color = (port_nd_fespace->GetVSize() > 0 || port_h1_fespace->GetVSize() > 0)
628 0 : ? 0
629 : : MPI_UNDEFINED;
630 0 : MPI_Comm_split(comm, color, Mpi::Rank(comm), &port_comm);
631 0 : MFEM_VERIFY((color == 0 && port_comm != MPI_COMM_NULL) ||
632 : (color == MPI_UNDEFINED && port_comm == MPI_COMM_NULL),
633 : "Unexpected error splitting communicator for wave port boundaries!");
634 0 : port_root = (color == MPI_UNDEFINED) ? Mpi::Size(comm) : Mpi::Rank(comm);
635 0 : Mpi::GlobalMin(1, &port_root, comm);
636 0 : MFEM_VERIFY(port_root < Mpi::Size(comm), "No root process found for port!");
637 :
638 : // Configure the eigenvalue problem solver. As for the full 3D case, the system matrices
639 : // are in general complex and symmetric. We supply the operators to the solver in
640 : // shift-inverted form and handle the back-transformation externally.
641 0 : if (port_comm != MPI_COMM_NULL)
642 : {
643 : // Define the linear solver to be used for solving systems associated with the
644 : // generalized eigenvalue problem.
645 0 : auto gmres = std::make_unique<GmresSolver<ComplexOperator>>(port_comm, data.verbose);
646 0 : gmres->SetInitialGuess(false);
647 0 : gmres->SetRelTol(data.ksp_tol);
648 0 : gmres->SetMaxIter(data.ksp_max_its);
649 : gmres->SetRestartDim(data.ksp_max_its);
650 : // gmres->SetPrecSide(PreconditionerSide::RIGHT);
651 :
652 0 : LinearSolver pc_type = solver.linear.type;
653 0 : if (pc_type == LinearSolver::SUPERLU)
654 : {
655 : #if !defined(MFEM_USE_SUPERLU)
656 0 : MFEM_ABORT("Solver was not built with SuperLU_DIST support, please choose a "
657 : "different solver!");
658 : #endif
659 : }
660 0 : else if (pc_type == LinearSolver::STRUMPACK || pc_type == LinearSolver::STRUMPACK_MP)
661 : {
662 : #if !defined(MFEM_USE_STRUMPACK)
663 : MFEM_ABORT("Solver was not built with STRUMPACK support, please choose a "
664 : "different solver!");
665 : #endif
666 : }
667 0 : else if (pc_type == LinearSolver::MUMPS)
668 : {
669 : #if !defined(MFEM_USE_MUMPS)
670 0 : MFEM_ABORT("Solver was not built with MUMPS support, please choose a "
671 : "different solver!");
672 : #endif
673 : }
674 : else // Default choice
675 : {
676 : #if defined(MFEM_USE_SUPERLU)
677 : pc_type = LinearSolver::SUPERLU;
678 : #elif defined(MFEM_USE_STRUMPACK)
679 0 : pc_type = LinearSolver::STRUMPACK;
680 : #elif defined(MFEM_USE_MUMPS)
681 : pc_type = LinearSolver::MUMPS;
682 : #else
683 : #error "Wave port solver requires building with SuperLU_DIST, STRUMPACK, or MUMPS!"
684 : #endif
685 : }
686 : auto pc = std::make_unique<MfemWrapperSolver<ComplexOperator>>(
687 0 : [&]() -> std::unique_ptr<mfem::Solver>
688 : {
689 0 : if (pc_type == LinearSolver::SUPERLU)
690 : {
691 : #if defined(MFEM_USE_SUPERLU)
692 : auto slu = std::make_unique<SuperLUSolver>(
693 : port_comm, SymbolicFactorization::DEFAULT, false, data.verbose - 1);
694 : // slu->GetSolver().SetColumnPermutation(mfem::superlu::MMD_AT_PLUS_A);
695 : return slu;
696 : #endif
697 : }
698 0 : else if (pc_type == LinearSolver::STRUMPACK)
699 : {
700 : #if defined(MFEM_USE_STRUMPACK)
701 : auto strumpack = std::make_unique<StrumpackSolver>(
702 0 : port_comm, SymbolicFactorization::DEFAULT, SparseCompression::NONE, 0.0, 0,
703 0 : 0, data.verbose - 1);
704 : // strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::AMD);
705 : return strumpack;
706 : #endif
707 : }
708 : else if (pc_type == LinearSolver::MUMPS)
709 : {
710 : #if defined(MFEM_USE_MUMPS)
711 : auto mumps = std::make_unique<MumpsSolver>(
712 : port_comm, mfem::MUMPSSolver::UNSYMMETRIC, SymbolicFactorization::DEFAULT,
713 : 0.0, data.verbose - 1);
714 : // mumps->SetReorderingStrategy(mfem::MUMPSSolver::AMD);
715 : return mumps;
716 : #endif
717 : }
718 : return {};
719 0 : }());
720 : pc->SetSaveAssembled(false);
721 0 : ksp = std::make_unique<ComplexKspSolver>(std::move(gmres), std::move(pc));
722 :
723 : // Define the eigenvalue solver.
724 0 : constexpr int print = 0;
725 0 : EigenSolverBackend type = data.eigen_solver;
726 0 : if (type == EigenSolverBackend::SLEPC)
727 : {
728 : #if !defined(PALACE_WITH_SLEPC)
729 : MFEM_ABORT("Solver was not built with SLEPc support, please choose a "
730 : "different solver!");
731 : #endif
732 : }
733 0 : else if (type == EigenSolverBackend::ARPACK)
734 : {
735 : #if !defined(PALACE_WITH_ARPACK)
736 0 : MFEM_ABORT("Solver was not built with ARPACK support, please choose a "
737 : "different solver!");
738 : #endif
739 : }
740 : else // Default choice
741 : {
742 : #if defined(PALACE_WITH_SLEPC)
743 : type = EigenSolverBackend::SLEPC;
744 : #elif defined(PALACE_WITH_ARPACK)
745 : type = EigenSolverBackend::ARPACK;
746 : #else
747 : #error "Wave port solver requires building with ARPACK or SLEPc!"
748 : #endif
749 : }
750 : if (type == EigenSolverBackend::ARPACK)
751 : {
752 : #if defined(PALACE_WITH_ARPACK)
753 : eigen = std::make_unique<arpack::ArpackEPSSolver>(port_comm, print);
754 : #endif
755 : }
756 : else // EigenSolverBackend::SLEPC
757 : {
758 : #if defined(PALACE_WITH_SLEPC)
759 0 : auto slepc = std::make_unique<slepc::SlepcEPSSolver>(port_comm, print);
760 0 : slepc->SetType(slepc::SlepcEigenvalueSolver::Type::KRYLOVSCHUR);
761 0 : slepc->SetProblemType(slepc::SlepcEigenvalueSolver::ProblemType::GEN_NON_HERMITIAN);
762 : eigen = std::move(slepc);
763 : #endif
764 : }
765 0 : eigen->SetNumModes(mode_idx, std::max(2 * mode_idx + 1, 5));
766 0 : eigen->SetTol(data.eig_tol);
767 0 : eigen->SetLinearSolver(*ksp);
768 :
769 : // We want to ignore evanescent modes (kₙ with large imaginary component). The
770 : // eigenvalue 1 / (-kₙ² - σ) of the shifted problem will be a large-magnitude positive
771 : // real number for an eigenvalue kₙ² with real part close to but not above the cutoff σ.
772 0 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::LARGEST_REAL);
773 : }
774 :
775 : // Configure port mode sign convention: 1ᵀ Re{-n x H} >= 0 on the "upper-right quadrant"
776 : // of the wave port boundary, in order to deal with symmetry effectively.
777 : {
778 : Vector bbmin, bbmax;
779 0 : mesh::GetAxisAlignedBoundingBox(*port_mesh, bbmin, bbmax);
780 : const int dim = port_mesh->SpaceDimension();
781 :
782 : double la = 0.0, lb = 0.0;
783 : int da = -1, db = -1;
784 0 : for (int d = 0; d < dim; d++)
785 : {
786 0 : double diff = bbmax(d) - bbmin(d);
787 0 : if (diff > la)
788 : {
789 : lb = la;
790 : la = diff;
791 : db = da;
792 : da = d;
793 : }
794 0 : else if (diff > lb)
795 : {
796 : lb = diff;
797 : db = d;
798 : }
799 : }
800 0 : MFEM_VERIFY(da >= 0 && db >= 0 && da != db,
801 : "Unexpected wave port geometry for normalization!");
802 0 : double ca = 0.5 * (bbmax[da] + bbmin[da]), cb = 0.5 * (bbmax[db] + bbmin[db]);
803 :
804 0 : auto TDirection = [da, db, ca, cb, dim](const Vector &x, Vector &f)
805 : {
806 : MFEM_ASSERT(x.Size() == dim,
807 : "Invalid dimension mismatch for wave port mode normalization!");
808 0 : f.SetSize(dim);
809 0 : if (x[da] >= ca && x[db] >= cb)
810 : {
811 0 : f = 1.0;
812 : }
813 : else
814 : {
815 0 : f = 0.0;
816 : }
817 0 : };
818 0 : mfem::VectorFunctionCoefficient tfunc(dim, TDirection);
819 0 : port_S0t = std::make_unique<GridFunction>(*port_nd_fespace);
820 0 : port_S0t->Real().ProjectCoefficient(tfunc);
821 0 : }
822 0 : }
823 :
824 0 : WavePortData::~WavePortData()
825 : {
826 : // Free the solvers before the communicator on which they are based.
827 : ksp.reset();
828 : eigen.reset();
829 0 : if (port_comm != MPI_COMM_NULL)
830 : {
831 0 : MPI_Comm_free(&port_comm);
832 : }
833 0 : }
834 :
835 0 : void WavePortData::Initialize(double omega)
836 : {
837 0 : if (omega == omega0)
838 : {
839 0 : return;
840 : }
841 :
842 : // Construct matrices and solve the generalized eigenvalue problem for the desired wave
843 : // port mode. The B matrix is operating frequency-independent and has already been
844 : // constructed.
845 : std::unique_ptr<ComplexOperator> opA;
846 0 : const double sigma = -omega * omega * mu_eps_max;
847 : {
848 0 : auto [Attr, Atti] = GetAtt(mat_op, *port_nd_fespace, port_normal, omega, sigma);
849 : auto [Ar, Ai] =
850 : GetSystemMatrixA(Attr.get(), Atti.get(), Atnr.get(), Atni.get(), Antr.get(),
851 0 : Anti.get(), Annr.get(), Anni.get(), port_dbc_tdof_list);
852 0 : opA = std::make_unique<ComplexWrapperOperator>(std::move(Ar), std::move(Ai));
853 : }
854 :
855 : // Configure and solve the (inverse) eigenvalue problem for the desired boundary mode.
856 : // Linear solves are preconditioned with the real part of the system matrix (ignore loss
857 : // tangent).
858 0 : std::complex<double> lambda;
859 0 : if (port_comm != MPI_COMM_NULL)
860 : {
861 0 : ComplexWrapperOperator opP(opA->Real(), nullptr); // Non-owning constructor
862 0 : ksp->SetOperators(*opA, opP);
863 0 : eigen->SetOperators(*opB, *opA, EigenvalueSolver::ScaleType::NONE);
864 0 : eigen->SetInitialSpace(v0);
865 0 : int num_conv = eigen->Solve();
866 0 : MFEM_VERIFY(num_conv >= mode_idx, "Wave port eigensolver did not converge!");
867 0 : lambda = eigen->GetEigenvalue(mode_idx - 1);
868 : // Mpi::Print(port_comm, " ... Wave port eigensolver error = {} (bkwd), {} (abs)\n",
869 : // eigen->GetError(mode_idx - 1, EigenvalueSolver::ErrorType::BACKWARD),
870 : // eigen->GetError(mode_idx - 1, EigenvalueSolver::ErrorType::ABSOLUTE));
871 0 : }
872 0 : Mpi::Broadcast(1, &lambda, port_root, port_mesh->GetComm());
873 :
874 : // Extract the eigenmode solution and postprocess. The extracted eigenvalue is λ =
875 : // 1 / (-kₙ² - σ).
876 0 : kn0 = std::sqrt(-sigma - 1.0 / lambda);
877 0 : omega0 = omega;
878 :
879 : // Separate the computed field out into eₜ and eₙ and and transform back to true
880 : // electric field variables: Eₜ = eₜ and Eₙ = eₙ / ikₙ.
881 : {
882 0 : if (port_comm != MPI_COMM_NULL)
883 : {
884 0 : eigen->GetEigenvector(mode_idx - 1, e0);
885 0 : linalg::NormalizePhase(port_comm, e0);
886 : }
887 : else
888 : {
889 : MFEM_ASSERT(e0.Size() == 0,
890 : "Unexpected non-empty port FE space in wave port boundary mode solve!");
891 : }
892 0 : e0.Real().Read(); // Ensure memory is allocated on device before aliasing
893 0 : e0.Imag().Read();
894 : Vector e0tr(e0.Real(), 0, port_nd_fespace->GetTrueVSize());
895 : Vector e0nr(e0.Real(), port_nd_fespace->GetTrueVSize(),
896 : port_h1_fespace->GetTrueVSize());
897 : Vector e0ti(e0.Imag(), 0, port_nd_fespace->GetTrueVSize());
898 : Vector e0ni(e0.Imag(), port_nd_fespace->GetTrueVSize(),
899 : port_h1_fespace->GetTrueVSize());
900 : e0tr.UseDevice(true);
901 : e0nr.UseDevice(true);
902 : e0ti.UseDevice(true);
903 : e0ni.UseDevice(true);
904 0 : ComplexVector::AXPBY(1.0 / (1i * kn0), e0nr, e0ni, 0.0, e0nr, e0ni);
905 0 : port_E0t->Real().SetFromTrueDofs(e0tr); // Parallel distribute
906 0 : port_E0t->Imag().SetFromTrueDofs(e0ti);
907 0 : port_E0n->Real().SetFromTrueDofs(e0nr);
908 0 : port_E0n->Imag().SetFromTrueDofs(e0ni);
909 : }
910 :
911 : // Configure the linear forms for computing S-parameters (projection of the field onto the
912 : // port mode). Normalize the mode for a chosen polarization direction and unit power,
913 : // |E x H⋆| ⋅ n, integrated over the port surface (+n is the direction of propagation).
914 : {
915 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
916 : BdrSubmeshHVectorCoefficient<ValueType::REAL> port_nxH0r_func(
917 0 : *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
918 : BdrSubmeshHVectorCoefficient<ValueType::IMAG> port_nxH0i_func(
919 0 : *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
920 : {
921 0 : port_sr = std::make_unique<mfem::LinearForm>(&port_nd_fespace->Get());
922 0 : port_sr->AddDomainIntegrator(new VectorFEDomainLFIntegrator(port_nxH0r_func));
923 0 : port_sr->UseFastAssembly(false);
924 0 : port_sr->UseDevice(false);
925 0 : port_sr->Assemble();
926 0 : port_sr->UseDevice(true);
927 : }
928 : {
929 0 : port_si = std::make_unique<mfem::LinearForm>(&port_nd_fespace->Get());
930 0 : port_si->AddDomainIntegrator(new VectorFEDomainLFIntegrator(port_nxH0i_func));
931 0 : port_si->UseFastAssembly(false);
932 0 : port_si->UseDevice(false);
933 0 : port_si->Assemble();
934 0 : port_si->UseDevice(true);
935 : }
936 0 : Normalize(*port_S0t, *port_E0t, *port_E0n, *port_sr, *port_si);
937 : }
938 : }
939 :
940 : std::unique_ptr<mfem::VectorCoefficient>
941 0 : WavePortData::GetModeExcitationCoefficientReal() const
942 : {
943 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
944 : return std::make_unique<
945 0 : RestrictedVectorCoefficient<BdrSubmeshHVectorCoefficient<ValueType::REAL>>>(
946 0 : attr_list, *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0,
947 0 : omega0);
948 : }
949 :
950 : std::unique_ptr<mfem::VectorCoefficient>
951 0 : WavePortData::GetModeExcitationCoefficientImag() const
952 : {
953 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
954 : return std::make_unique<
955 0 : RestrictedVectorCoefficient<BdrSubmeshHVectorCoefficient<ValueType::IMAG>>>(
956 0 : attr_list, *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0,
957 0 : omega0);
958 : }
959 :
960 0 : std::unique_ptr<mfem::VectorCoefficient> WavePortData::GetModeFieldCoefficientReal() const
961 : {
962 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
963 : return std::make_unique<
964 0 : RestrictedVectorCoefficient<BdrSubmeshEVectorCoefficient<ValueType::REAL>>>(
965 0 : attr_list, *port_E0t, *port_E0n, port_submesh, submesh_parent_elems);
966 : }
967 :
968 0 : std::unique_ptr<mfem::VectorCoefficient> WavePortData::GetModeFieldCoefficientImag() const
969 : {
970 : const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
971 : return std::make_unique<
972 0 : RestrictedVectorCoefficient<BdrSubmeshEVectorCoefficient<ValueType::IMAG>>>(
973 0 : attr_list, *port_E0t, *port_E0n, port_submesh, submesh_parent_elems);
974 : }
975 :
976 0 : double WavePortData::GetExcitationPower() const
977 : {
978 : // The computed port modes are normalized such that the power integrated over the port is
979 : // 1: ∫ (E_inc x H_inc⋆) ⋅ n dS = 1.
980 0 : return HasExcitation() ? 1.0 : 0.0;
981 : }
982 :
983 0 : std::complex<double> WavePortData::GetPower(GridFunction &E, GridFunction &B) const
984 : {
985 : // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface using
986 : // the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation (into the
987 : // domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an outward normal,
988 : // so we multiply by -1. The linear form is reconstructed from scratch each time due to
989 : // changing H.
990 0 : MFEM_VERIFY(E.HasImag() && B.HasImag(),
991 : "Wave ports expect complex-valued E and B fields in port power "
992 : "calculation!");
993 : auto &nd_fespace = *E.ParFESpace();
994 : const auto &mesh = *nd_fespace.GetParMesh();
995 0 : BdrSurfaceCurrentVectorCoefficient nxHr_func(B.Real(), mat_op);
996 0 : BdrSurfaceCurrentVectorCoefficient nxHi_func(B.Imag(), mat_op);
997 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
998 0 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
999 0 : std::complex<double> dot;
1000 : {
1001 0 : mfem::LinearForm pr(&nd_fespace);
1002 0 : pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(nxHr_func), attr_marker);
1003 0 : pr.UseFastAssembly(false);
1004 : pr.UseDevice(false);
1005 0 : pr.Assemble();
1006 : pr.UseDevice(true);
1007 0 : dot = -(pr * E.Real()) - 1i * (pr * E.Imag());
1008 0 : }
1009 : {
1010 0 : mfem::LinearForm pi(&nd_fespace);
1011 0 : pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(nxHi_func), attr_marker);
1012 0 : pi.UseFastAssembly(false);
1013 : pi.UseDevice(false);
1014 0 : pi.Assemble();
1015 : pi.UseDevice(true);
1016 0 : dot += -(pi * E.Imag()) + 1i * (pi * E.Real());
1017 0 : }
1018 : Mpi::GlobalSum(1, &dot, nd_fespace.GetComm());
1019 0 : return dot;
1020 : }
1021 :
1022 0 : std::complex<double> WavePortData::GetSParameter(GridFunction &E) const
1023 : {
1024 : // Compute port S-parameter, or the projection of the field onto the port mode:
1025 : // (E x H_inc⋆) ⋅ n = E ⋅ (-n x H_inc⋆), integrated over the port surface.
1026 0 : MFEM_VERIFY(E.HasImag(),
1027 : "Wave ports expect complex-valued E and B fields in port S-parameter "
1028 : "calculation!");
1029 0 : port_nd_transfer->Transfer(E.Real(), port_E->Real());
1030 0 : port_nd_transfer->Transfer(E.Imag(), port_E->Imag());
1031 0 : std::complex<double> dot(-((*port_sr) * port_E->Real()) - ((*port_si) * port_E->Imag()),
1032 0 : -((*port_sr) * port_E->Imag()) + ((*port_si) * port_E->Real()));
1033 : Mpi::GlobalSum(1, &dot, port_nd_fespace->GetComm());
1034 0 : return dot;
1035 : }
1036 :
1037 17 : WavePortOperator::WavePortOperator(const IoData &iodata, const MaterialOperator &mat_op,
1038 : mfem::ParFiniteElementSpace &nd_fespace,
1039 17 : mfem::ParFiniteElementSpace &h1_fespace)
1040 17 : : suppress_output(false),
1041 17 : fc(iodata.units.Dimensionalize<Units::ValueType::FREQUENCY>(1.0)),
1042 17 : kc(1.0 / iodata.units.Dimensionalize<Units::ValueType::LENGTH>(1.0))
1043 : {
1044 : // Set up wave port boundary conditions.
1045 17 : MFEM_VERIFY(nd_fespace.GetParMesh() == h1_fespace.GetParMesh(),
1046 : "Mesh mismatch in WavePortOperator FE spaces!");
1047 17 : SetUpBoundaryProperties(iodata, mat_op, nd_fespace, h1_fespace);
1048 17 : PrintBoundaryInfo(iodata, *nd_fespace.GetParMesh());
1049 17 : }
1050 :
1051 17 : void WavePortOperator::SetUpBoundaryProperties(const IoData &iodata,
1052 : const MaterialOperator &mat_op,
1053 : mfem::ParFiniteElementSpace &nd_fespace,
1054 : mfem::ParFiniteElementSpace &h1_fespace)
1055 : {
1056 : // Check that wave port boundary attributes have been specified correctly.
1057 : const auto &mesh = *nd_fespace.GetParMesh();
1058 17 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
1059 17 : if (!iodata.boundaries.waveport.empty())
1060 : {
1061 : mfem::Array<int> bdr_attr_marker(bdr_attr_max), port_marker(bdr_attr_max);
1062 : bdr_attr_marker = 0;
1063 : port_marker = 0;
1064 0 : for (auto attr : mesh.bdr_attributes)
1065 : {
1066 0 : bdr_attr_marker[attr - 1] = 1;
1067 : }
1068 0 : for (const auto &[idx, data] : iodata.boundaries.waveport)
1069 : {
1070 0 : for (auto attr : data.attributes)
1071 : {
1072 0 : MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
1073 : "Port boundary attribute tags must be non-negative and correspond to "
1074 : "boundaries in the mesh!");
1075 0 : MFEM_VERIFY(bdr_attr_marker[attr - 1],
1076 : "Unknown port boundary attribute " << attr << "!");
1077 0 : MFEM_VERIFY(!data.active || !port_marker[attr - 1],
1078 : "Boundary attribute is assigned to more than one wave port!");
1079 0 : port_marker[attr - 1] = 1;
1080 : }
1081 : }
1082 : }
1083 :
1084 : // List of all boundaries which will be marked as essential for the purposes of computing
1085 : // wave port modes. This includes all PEC surfaces, but may also include others like when
1086 : // a kinetic inductance or other BC is applied for the 3D simulation but should be
1087 : // considered as PEC for the 2D problem. In addition, we mark as Dirichlet boundaries all
1088 : // wave ports other than the wave port being currently considered, in case two wave ports
1089 : // are touching and share one or more edges.
1090 : mfem::Array<int> dbc_bcs;
1091 17 : dbc_bcs.Reserve(static_cast<int>(iodata.boundaries.pec.attributes.size() +
1092 17 : iodata.boundaries.auxpec.attributes.size() +
1093 : iodata.boundaries.conductivity.size()));
1094 26 : for (auto attr : iodata.boundaries.pec.attributes)
1095 : {
1096 9 : if (attr <= 0 || attr > bdr_attr_max)
1097 : {
1098 0 : continue; // Can just ignore if wrong
1099 : }
1100 9 : dbc_bcs.Append(attr);
1101 : }
1102 17 : for (auto attr : iodata.boundaries.auxpec.attributes)
1103 : {
1104 0 : if (attr <= 0 || attr > bdr_attr_max)
1105 : {
1106 0 : continue; // Can just ignore if wrong
1107 : }
1108 0 : dbc_bcs.Append(attr);
1109 : }
1110 17 : for (const auto &data : iodata.boundaries.conductivity)
1111 : {
1112 0 : for (auto attr : data.attributes)
1113 : {
1114 0 : if (attr <= 0 || attr > bdr_attr_max)
1115 : {
1116 0 : continue; // Can just ignore if wrong
1117 : }
1118 0 : dbc_bcs.Append(attr);
1119 : }
1120 : }
1121 : // If user accidentally specifies a surface as both "PEC" and "WavePortPEC", this is fine
1122 : // so allow for duplicates in the attribute list.
1123 : dbc_bcs.Sort();
1124 17 : dbc_bcs.Unique();
1125 :
1126 : // Set up wave port data structures.
1127 17 : for (const auto &[idx, data] : iodata.boundaries.waveport)
1128 : {
1129 0 : mfem::Array<int> port_dbc_bcs(dbc_bcs);
1130 0 : for (const auto &[other_idx, other_data] : iodata.boundaries.waveport)
1131 : {
1132 0 : if (other_idx == idx || !other_data.active)
1133 : {
1134 0 : continue;
1135 : }
1136 0 : for (auto attr : other_data.attributes)
1137 : {
1138 0 : if (std::binary_search(data.attributes.begin(), data.attributes.end(), attr))
1139 : {
1140 0 : continue;
1141 : }
1142 0 : port_dbc_bcs.Append(attr);
1143 : }
1144 : }
1145 : port_dbc_bcs.Sort();
1146 0 : port_dbc_bcs.Unique();
1147 0 : ports.try_emplace(idx, data, iodata.solver, mat_op, nd_fespace, h1_fespace,
1148 : port_dbc_bcs);
1149 : }
1150 17 : MFEM_VERIFY(
1151 : ports.empty() || iodata.problem.type == ProblemType::DRIVEN ||
1152 : iodata.problem.type == ProblemType::EIGENMODE,
1153 : "Wave port boundaries are only available for frequency domain driven simulations!");
1154 17 : }
1155 :
1156 17 : void WavePortOperator::PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh)
1157 : {
1158 17 : if (ports.empty())
1159 : {
1160 17 : return;
1161 : }
1162 : fmt::memory_buffer buf{}; // Output buffer & buffer append lambda for cleaner code
1163 0 : auto to = [&buf](auto fmt, auto &&...args)
1164 0 : { fmt::format_to(std::back_inserter(buf), fmt, std::forward<decltype(args)>(args)...); };
1165 :
1166 : // Print out BC info for all active port attributes.
1167 0 : for (const auto &[idx, data] : ports)
1168 : {
1169 0 : if (!data.active)
1170 : {
1171 0 : continue;
1172 : }
1173 0 : for (auto attr : data.GetAttrList())
1174 : {
1175 0 : to(" {:d}: Index = {:d}, mode = {:d}, d = {:.3e} m, n = ({:+.1f})\n", attr, idx,
1176 0 : data.mode_idx,
1177 0 : iodata.units.Dimensionalize<Units::ValueType::LENGTH>(data.d_offset),
1178 0 : fmt::join(data.port_normal, ","));
1179 : }
1180 : }
1181 0 : if (buf.size() > 0)
1182 : {
1183 0 : Mpi::Print("\nConfiguring Robin impedance BC for wave ports at attributes:\n");
1184 0 : Mpi::Print("{}", fmt::to_string(buf));
1185 : buf.clear();
1186 : }
1187 :
1188 : // Print some information for excited wave ports.
1189 0 : for (const auto &[idx, data] : ports)
1190 : {
1191 0 : if (!data.HasExcitation())
1192 : {
1193 0 : continue;
1194 : }
1195 0 : for (auto attr : data.GetAttrList())
1196 : {
1197 0 : to(" {:d}: Index = {:d}\n", attr, idx);
1198 : }
1199 : }
1200 0 : if (buf.size() > 0)
1201 : {
1202 0 : Mpi::Print("\nConfiguring wave port excitation source term at attributes:\n");
1203 0 : Mpi::Print("{}", fmt::to_string(buf));
1204 : }
1205 : }
1206 :
1207 0 : const WavePortData &WavePortOperator::GetPort(int idx) const
1208 : {
1209 : auto it = ports.find(idx);
1210 0 : MFEM_VERIFY(it != ports.end(), "Unknown wave port index requested!");
1211 0 : return it->second;
1212 : }
1213 :
1214 17 : mfem::Array<int> WavePortOperator::GetAttrList() const
1215 : {
1216 : mfem::Array<int> attr_list;
1217 17 : for (const auto &[idx, data] : ports)
1218 : {
1219 0 : if (!data.active)
1220 : {
1221 0 : continue;
1222 : }
1223 : attr_list.Append(data.GetAttrList());
1224 : }
1225 17 : return attr_list;
1226 : }
1227 :
1228 0 : void WavePortOperator::Initialize(double omega)
1229 : {
1230 : bool init = false, first = true;
1231 0 : for (const auto &[idx, data] : ports)
1232 : {
1233 0 : init = init || (data.omega0 != omega);
1234 0 : first = first && (data.omega0 == 0.0);
1235 : }
1236 0 : if (!init)
1237 : {
1238 0 : return;
1239 : }
1240 0 : BlockTimer bt(Timer::WAVE_PORT);
1241 0 : if (!suppress_output)
1242 : {
1243 0 : Mpi::Print(
1244 : "\nCalculating boundary modes at wave ports for ω/2π = {:.3e} GHz ({:.3e})\n",
1245 0 : omega * fc, omega);
1246 : }
1247 0 : for (auto &[idx, data] : ports)
1248 : {
1249 0 : data.Initialize(omega);
1250 0 : if (!suppress_output)
1251 : {
1252 0 : if (first)
1253 : {
1254 0 : Mpi::Print(" Number of global unknowns for port {:d}:\n"
1255 : " H1: {:d}, ND: {:d}\n",
1256 0 : idx, data.GlobalTrueH1Size(), data.GlobalTrueNDSize());
1257 : }
1258 0 : Mpi::Print(" Port {:d}, mode {:d}: kₙ = {:.3e}{:+.3e}i m⁻¹\n", idx, data.mode_idx,
1259 0 : data.kn0.real() * kc, data.kn0.imag() * kc);
1260 : }
1261 : }
1262 0 : }
1263 :
1264 0 : void WavePortOperator::AddExtraSystemBdrCoefficients(double omega,
1265 : MaterialPropertyCoefficient &fbr,
1266 : MaterialPropertyCoefficient &fbi)
1267 : {
1268 : // Add wave port boundaries to the bilinear form. This looks a lot like the lumped port
1269 : // boundary, except the iω / Z_s coefficient goes to ikₙ / μ where kₙ is specific to the
1270 : // port mode at the given operating frequency (note only the real part of the propagation
1271 : // constant contributes).
1272 0 : Initialize(omega);
1273 0 : for (const auto &[idx, data] : ports)
1274 : {
1275 0 : if (!data.active)
1276 : {
1277 0 : continue;
1278 : }
1279 0 : const MaterialOperator &mat_op = data.mat_op;
1280 0 : MaterialPropertyCoefficient muinv_func(mat_op.GetBdrAttributeToMaterial(),
1281 0 : mat_op.GetInvPermeability());
1282 0 : muinv_func.RestrictCoefficient(mat_op.GetCeedBdrAttributes(data.GetAttrList()));
1283 : // fbr.AddCoefficient(muinv_func.GetAttributeToMaterial(),
1284 : // muinv_func.GetMaterialProperties(),
1285 : // -data.kn0.imag());
1286 0 : fbi.AddCoefficient(muinv_func.GetAttributeToMaterial(),
1287 : muinv_func.GetMaterialProperties(), data.kn0.real());
1288 0 : }
1289 0 : }
1290 :
1291 0 : void WavePortOperator::AddExcitationBdrCoefficients(int excitation_idx, double omega,
1292 : SumVectorCoefficient &fbr,
1293 : SumVectorCoefficient &fbi)
1294 : {
1295 : // Re/Im{-U_inc} = Re/Im{+2 (-iω) n x H_inc}, which is a function of E_inc as computed by
1296 : // the modal solution (stored as a grid function and coefficient during initialization).
1297 0 : Initialize(omega);
1298 0 : for (const auto &[idx, data] : ports)
1299 : {
1300 0 : if (data.excitation != excitation_idx)
1301 : {
1302 0 : continue;
1303 : }
1304 0 : fbr.AddCoefficient(data.GetModeExcitationCoefficientImag(), 2.0 * omega);
1305 0 : fbi.AddCoefficient(data.GetModeExcitationCoefficientReal(), -2.0 * omega);
1306 : }
1307 0 : }
1308 :
1309 : } // namespace palace
|
