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 "lumpedportoperator.hpp"
5 :
6 : #include <fmt/ranges.h>
7 : #include "fem/coefficient.hpp"
8 : #include "fem/gridfunction.hpp"
9 : #include "fem/integrator.hpp"
10 : #include "models/materialoperator.hpp"
11 : #include "utils/communication.hpp"
12 : #include "utils/geodata.hpp"
13 : #include "utils/iodata.hpp"
14 :
15 : namespace palace
16 : {
17 :
18 : using namespace std::complex_literals;
19 :
20 30 : LumpedPortData::LumpedPortData(const config::LumpedPortData &data,
21 30 : const MaterialOperator &mat_op, const mfem::ParMesh &mesh)
22 30 : : mat_op(mat_op), excitation(data.excitation), active(data.active)
23 : {
24 : // Check inputs. Only one of the circuit or per square properties should be specified
25 : // for the port boundary.
26 30 : bool has_circ = (std::abs(data.R) + std::abs(data.L) + std::abs(data.C) > 0.0);
27 30 : bool has_surf = (std::abs(data.Rs) + std::abs(data.Ls) + std::abs(data.Cs) > 0.0);
28 30 : MFEM_VERIFY(has_circ || has_surf,
29 : "Lumped port boundary has no R/L/C or Rs/Ls/Cs defined, needs "
30 : "at least one!");
31 30 : MFEM_VERIFY(!(has_circ && has_surf),
32 : "Lumped port boundary has both R/L/C and Rs/Ls/Cs defined, "
33 : "should only use one!");
34 :
35 30 : if (HasExcitation())
36 : {
37 16 : if (has_circ)
38 : {
39 16 : MFEM_VERIFY(data.R > 0.0, "Excited lumped port must have nonzero resistance!");
40 16 : MFEM_VERIFY(data.C == 0.0 && data.L == 0.0,
41 : "Lumped port excitations do not support nonzero reactance!");
42 : }
43 : else
44 : {
45 0 : MFEM_VERIFY(data.Rs > 0.0, "Excited lumped port must have nonzero resistance!");
46 0 : MFEM_VERIFY(data.Cs == 0.0 && data.Ls == 0.0,
47 : "Lumped port excitations do not support nonzero reactance!");
48 : }
49 : }
50 :
51 : // Construct the port elements allowing for a possible multielement lumped port.
52 60 : for (const auto &elem : data.elements)
53 : {
54 : mfem::Array<int> attr_list;
55 30 : attr_list.Append(elem.attributes.data(), elem.attributes.size());
56 30 : switch (elem.coordinate_system)
57 : {
58 0 : case CoordinateSystem::CYLINDRICAL:
59 0 : elems.push_back(
60 0 : std::make_unique<CoaxialElementData>(elem.direction, attr_list, mesh));
61 0 : break;
62 30 : case CoordinateSystem::CARTESIAN:
63 60 : elems.push_back(
64 30 : std::make_unique<UniformElementData>(elem.direction, attr_list, mesh));
65 30 : break;
66 : }
67 : }
68 :
69 : // Populate the property data for the lumped port.
70 30 : if (std::abs(data.Rs) + std::abs(data.Ls) + std::abs(data.Cs) == 0.0)
71 : {
72 30 : R = data.R;
73 30 : L = data.L;
74 30 : C = data.C;
75 : }
76 : else
77 : {
78 : // If defined by surface properties, need to compute circuit properties for the
79 : // multielement port.
80 : double ooR = 0.0, ooL = 0.0;
81 0 : R = L = C = 0.0;
82 0 : for (const auto &elem : elems)
83 : {
84 0 : const double sq = elem->GetGeometryWidth() / elem->GetGeometryLength();
85 0 : if (std::abs(data.Rs) > 0.0)
86 : {
87 0 : ooR += sq / data.Rs;
88 : }
89 0 : if (std::abs(data.Ls) > 0.0)
90 : {
91 0 : ooL += sq / data.Ls;
92 : }
93 0 : if (std::abs(data.Cs) > 0.0)
94 : {
95 0 : C += sq * data.Cs;
96 : }
97 : }
98 0 : if (std::abs(ooR) > 0.0)
99 : {
100 0 : R = 1.0 / ooR;
101 : }
102 0 : if (std::abs(ooL) > 0.0)
103 : {
104 0 : L = 1.0 / ooL;
105 : }
106 : }
107 30 : }
108 :
109 : std::complex<double>
110 3 : LumpedPortData::GetCharacteristicImpedance(double omega,
111 : LumpedPortData::Branch branch) const
112 : {
113 3 : MFEM_VERIFY((L == 0.0 && C == 0.0) || branch == Branch::R || omega > 0.0,
114 : "Lumped port with nonzero reactance requires frequency in order to define "
115 : "characteristic impedance!");
116 3 : std::complex<double> Y = 0.0;
117 3 : if (std::abs(R) > 0.0 && (branch == Branch::TOTAL || branch == Branch::R))
118 : {
119 3 : Y += 1.0 / R;
120 : }
121 3 : if (std::abs(L) > 0.0 && (branch == Branch::TOTAL || branch == Branch::L))
122 : {
123 : Y += 1.0 / (1i * omega * L);
124 : }
125 3 : if (std::abs(C) > 0.0 && (branch == Branch::TOTAL || branch == Branch::C))
126 : {
127 : Y += 1i * omega * C;
128 : }
129 3 : MFEM_VERIFY(std::abs(Y) > 0.0,
130 : "Characteristic impedance requested for lumped port with zero admittance!")
131 3 : return 1.0 / Y;
132 : }
133 :
134 2 : double LumpedPortData::GetExcitationPower() const
135 : {
136 : // The lumped port excitation is normalized such that the power integrated over the port
137 : // is 1: ∫ (E_inc x H_inc) ⋅ n dS = 1.
138 2 : return HasExcitation() ? 1.0 : 0.0;
139 : }
140 :
141 4 : double LumpedPortData::GetExcitationVoltage() const
142 : {
143 : // Incident voltage should be the same across all elements of an excited lumped port.
144 4 : if (HasExcitation())
145 : {
146 : double V_inc = 0.0;
147 8 : for (const auto &elem : elems)
148 : {
149 4 : const double Rs = R * GetToSquare(*elem);
150 4 : const double E_inc = std::sqrt(
151 4 : Rs / (elem->GetGeometryWidth() * elem->GetGeometryLength() * elems.size()));
152 4 : V_inc += E_inc * elem->GetGeometryLength() / elems.size();
153 : }
154 : return V_inc;
155 : }
156 : else
157 : {
158 : return 0.0;
159 : }
160 : }
161 :
162 9 : void LumpedPortData::InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespace) const
163 : {
164 : const auto &mesh = *nd_fespace.GetParMesh();
165 : mfem::Array<int> attr_marker;
166 9 : if (!s || !v)
167 : {
168 : mfem::Array<int> attr_list;
169 12 : for (const auto &elem : elems)
170 : {
171 : attr_list.Append(elem->GetAttrList());
172 : }
173 6 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
174 : mesh::AttrToMarker(bdr_attr_max, attr_list, attr_marker);
175 : }
176 :
177 : // The port S-parameter, or the projection of the field onto the port mode, is computed
178 : // as: (E x H_inc) ⋅ n = E ⋅ (E_inc / Z_s), integrated over the port surface.
179 9 : if (!s)
180 : {
181 : SumVectorCoefficient fb(mesh.SpaceDimension());
182 12 : for (const auto &elem : elems)
183 : {
184 6 : const double Rs = R * GetToSquare(*elem);
185 : const double Hinc = (std::abs(Rs) > 0.0)
186 6 : ? 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
187 6 : elem->GetGeometryLength() * elems.size())
188 : : 0.0;
189 12 : fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
190 : }
191 12 : s = std::make_unique<mfem::LinearForm>(&nd_fespace);
192 6 : s->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
193 : s->UseFastAssembly(false);
194 6 : s->UseDevice(false);
195 6 : s->Assemble();
196 6 : s->UseDevice(true);
197 : }
198 :
199 : // The voltage across a port is computed using the electric field solution.
200 : // We have:
201 : // V = ∫ E ⋅ l̂ dl = 1/w ∫ E ⋅ l̂ dS (for rectangular ports)
202 : // or,
203 : // V = 1/(2π) ∫ E ⋅ r̂ / r dS (for coaxial ports).
204 : // We compute the surface integral via an inner product between the linear form with the
205 : // averaging function as a vector coefficient and the solution expansion coefficients.
206 9 : if (!v)
207 : {
208 : SumVectorCoefficient fb(mesh.SpaceDimension());
209 12 : for (const auto &elem : elems)
210 : {
211 6 : fb.AddCoefficient(
212 12 : elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * elems.size())));
213 : }
214 12 : v = std::make_unique<mfem::LinearForm>(&nd_fespace);
215 6 : v->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
216 6 : v->UseFastAssembly(false);
217 6 : v->UseDevice(false);
218 6 : v->Assemble();
219 6 : v->UseDevice(true);
220 : }
221 9 : }
222 :
223 6 : std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B) const
224 : {
225 : // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface using
226 : // the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation (into the
227 : // domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an outward normal,
228 : // so we multiply by -1. The linear form is reconstructed from scratch each time due to
229 : // changing H.
230 6 : MFEM_VERIFY((E.HasImag() && B.HasImag()) || (!E.HasImag() && !B.HasImag()),
231 : "Mismatch between real- and complex-valued E and B fields in port power "
232 : "calculation!");
233 : const bool has_imag = E.HasImag();
234 : auto &nd_fespace = *E.ParFESpace();
235 : const auto &mesh = *nd_fespace.GetParMesh();
236 : SumVectorCoefficient fbr(mesh.SpaceDimension()), fbi(mesh.SpaceDimension());
237 : mfem::Array<int> attr_list;
238 12 : for (const auto &elem : elems)
239 : {
240 6 : fbr.AddCoefficient(
241 6 : std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
242 : elem->GetAttrList(), B.Real(), mat_op));
243 6 : if (has_imag)
244 : {
245 3 : fbi.AddCoefficient(
246 6 : std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
247 : elem->GetAttrList(), B.Imag(), mat_op));
248 : }
249 : attr_list.Append(elem->GetAttrList());
250 : }
251 6 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
252 6 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
253 6 : std::complex<double> dot;
254 : {
255 6 : mfem::LinearForm pr(&nd_fespace);
256 6 : pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbr), attr_marker);
257 6 : pr.UseFastAssembly(false);
258 : pr.UseDevice(false);
259 6 : pr.Assemble();
260 : pr.UseDevice(true);
261 9 : dot = -(pr * E.Real()) + (has_imag ? -1i * (pr * E.Imag()) : 0.0);
262 6 : }
263 6 : if (has_imag)
264 : {
265 3 : mfem::LinearForm pi(&nd_fespace);
266 3 : pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi), attr_marker);
267 3 : pi.UseFastAssembly(false);
268 : pi.UseDevice(false);
269 3 : pi.Assemble();
270 : pi.UseDevice(true);
271 3 : dot += -(pi * E.Imag()) + 1i * (pi * E.Real());
272 : Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
273 3 : return dot;
274 3 : }
275 : else
276 : {
277 3 : double rdot = dot.real();
278 : Mpi::GlobalSum(1, &rdot, E.ParFESpace()->GetComm());
279 3 : return rdot;
280 : }
281 : }
282 :
283 3 : std::complex<double> LumpedPortData::GetSParameter(GridFunction &E) const
284 : {
285 : // Compute port S-parameter, or the projection of the field onto the port mode.
286 3 : InitializeLinearForms(*E.ParFESpace());
287 3 : std::complex<double> dot((*s) * E.Real(), 0.0);
288 3 : if (E.HasImag())
289 : {
290 3 : dot.imag((*s) * E.Imag());
291 : }
292 : Mpi::GlobalSum(1, &dot, E.GetComm());
293 3 : return dot;
294 : }
295 :
296 6 : std::complex<double> LumpedPortData::GetVoltage(GridFunction &E) const
297 : {
298 : // Compute the average voltage across the port.
299 6 : InitializeLinearForms(*E.ParFESpace());
300 6 : std::complex<double> dot((*v) * E.Real(), 0.0);
301 6 : if (E.HasImag())
302 : {
303 3 : dot.imag((*v) * E.Imag());
304 : }
305 : Mpi::GlobalSum(1, &dot, E.GetComm());
306 6 : return dot;
307 : }
308 :
309 17 : LumpedPortOperator::LumpedPortOperator(const IoData &iodata, const MaterialOperator &mat_op,
310 : const mfem::ParMesh &mesh)
311 : {
312 : // Set up lumped port boundary conditions.
313 17 : SetUpBoundaryProperties(iodata, mat_op, mesh);
314 17 : PrintBoundaryInfo(iodata, mesh);
315 17 : }
316 :
317 17 : void LumpedPortOperator::SetUpBoundaryProperties(const IoData &iodata,
318 : const MaterialOperator &mat_op,
319 : const mfem::ParMesh &mesh)
320 : {
321 : // Check that lumped port boundary attributes have been specified correctly.
322 17 : if (!iodata.boundaries.lumpedport.empty())
323 : {
324 14 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
325 : mfem::Array<int> bdr_attr_marker(bdr_attr_max), port_marker(bdr_attr_max);
326 : bdr_attr_marker = 0;
327 : port_marker = 0;
328 242 : for (auto attr : mesh.bdr_attributes)
329 : {
330 228 : bdr_attr_marker[attr - 1] = 1;
331 : }
332 44 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
333 : {
334 60 : for (const auto &elem : data.elements)
335 : {
336 108 : for (auto attr : elem.attributes)
337 : {
338 78 : MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
339 : "Port boundary attribute tags must be non-negative and correspond to "
340 : "boundaries in the mesh!");
341 78 : MFEM_VERIFY(bdr_attr_marker[attr - 1],
342 : "Unknown port boundary attribute " << attr << "!");
343 78 : MFEM_VERIFY(!data.active || !port_marker[attr - 1],
344 : "Boundary attribute is assigned to more than one lumped port!");
345 78 : port_marker[attr - 1] = 1;
346 : }
347 : }
348 : }
349 : }
350 :
351 : // Set up lumped port data structures.
352 47 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
353 : {
354 30 : ports.try_emplace(idx, data, mat_op, mesh);
355 : }
356 17 : }
357 :
358 17 : void LumpedPortOperator::PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh)
359 : {
360 17 : if (ports.empty())
361 : {
362 3 : return;
363 : }
364 :
365 : fmt::memory_buffer buf{}; // Output buffer & buffer append lambda for cleaner code
366 434 : auto to = [](auto &buf, auto fmt, auto &&...args)
367 434 : { fmt::format_to(std::back_inserter(buf), fmt, std::forward<decltype(args)>(args)...); };
368 : using VT = Units::ValueType;
369 :
370 : // Print out BC info for all port attributes, for both active and inactive ports.
371 14 : to(buf, "\nConfiguring Robin impedance BC for lumped ports at attributes:\n");
372 44 : for (const auto &[idx, data] : ports)
373 : {
374 60 : for (const auto &elem : data.elems)
375 : {
376 108 : for (auto attr : elem->GetAttrList())
377 : {
378 78 : to(buf, " {:d}:", attr);
379 78 : if (std::abs(data.R) > 0.0)
380 : {
381 54 : double Rs = data.R * data.GetToSquare(*elem);
382 54 : to(buf, " Rs = {:.3e} Ω/sq,", iodata.units.Dimensionalize<VT::IMPEDANCE>(Rs));
383 : }
384 78 : if (std::abs(data.L) > 0.0)
385 : {
386 24 : double Ls = data.L * data.GetToSquare(*elem);
387 24 : to(buf, " Ls = {:.3e} H/sq,", iodata.units.Dimensionalize<VT::INDUCTANCE>(Ls));
388 : }
389 78 : if (std::abs(data.C) > 0.0)
390 : {
391 24 : double Cs = data.C / data.GetToSquare(*elem);
392 24 : to(buf, " Cs = {:.3e} F/sq,", iodata.units.Dimensionalize<VT::CAPACITANCE>(Cs));
393 : }
394 156 : to(buf, " n = ({:+.1f})\n", fmt::join(mesh::GetSurfaceNormal(mesh, attr), ","));
395 : }
396 : }
397 : }
398 :
399 : // Print out port info for all active ports.
400 : fmt::memory_buffer buf_a{};
401 44 : for (const auto &[idx, data] : ports)
402 : {
403 30 : if (!data.active)
404 : {
405 0 : continue;
406 : }
407 30 : to(buf_a, " Index = {:d}: ", idx);
408 30 : if (std::abs(data.R) > 0.0)
409 : {
410 22 : to(buf_a, "R = {:.3e} Ω,", iodata.units.Dimensionalize<VT::IMPEDANCE>(data.R));
411 : }
412 30 : if (std::abs(data.L) > 0.0)
413 : {
414 8 : to(buf_a, "L = {:.3e} H,", iodata.units.Dimensionalize<VT::INDUCTANCE>(data.L));
415 : }
416 30 : if (std::abs(data.C) > 0.0)
417 : {
418 8 : to(buf_a, "C = {:.3e} F,", iodata.units.Dimensionalize<VT::CAPACITANCE>(data.C));
419 : }
420 30 : buf_a.resize(buf_a.size() - 1); // Remove last ","
421 30 : to(buf_a, "\n");
422 : }
423 14 : if (buf_a.size() > 0)
424 : {
425 14 : to(buf, "\nConfiguring lumped port circuit properties:\n");
426 : buf.append(buf_a);
427 : buf_a.clear();
428 : }
429 :
430 : // Print some information for excited lumped ports.
431 44 : for (const auto &[idx, data] : ports)
432 : {
433 30 : if (!data.HasExcitation())
434 : {
435 14 : continue;
436 : }
437 :
438 32 : for (const auto &elem : data.elems)
439 : {
440 52 : for (auto attr : elem->GetAttrList())
441 : {
442 36 : to(buf_a, " {:d}: Index = {:d}\n", attr, idx);
443 : }
444 : }
445 : }
446 14 : if (buf_a.size() > 0)
447 : {
448 14 : to(buf, "\nConfiguring lumped port excitation source term at attributes:\n");
449 : buf.append(buf_a);
450 : }
451 :
452 28 : Mpi::Print("{}", fmt::to_string(buf));
453 : }
454 :
455 3 : const LumpedPortData &LumpedPortOperator::GetPort(int idx) const
456 : {
457 : auto it = ports.find(idx);
458 3 : MFEM_VERIFY(it != ports.end(), "Unknown lumped port index requested!");
459 3 : return it->second;
460 : }
461 :
462 17 : mfem::Array<int> LumpedPortOperator::GetAttrList() const
463 : {
464 : mfem::Array<int> attr_list;
465 47 : for (const auto &[idx, data] : ports)
466 : {
467 30 : if (!data.active)
468 : {
469 0 : continue;
470 : }
471 60 : for (const auto &elem : data.elems)
472 : {
473 : attr_list.Append(elem->GetAttrList());
474 : }
475 : }
476 17 : return attr_list;
477 : }
478 :
479 17 : mfem::Array<int> LumpedPortOperator::GetRsAttrList() const
480 : {
481 : mfem::Array<int> attr_list;
482 47 : for (const auto &[idx, data] : ports)
483 : {
484 30 : if (!data.active)
485 : {
486 0 : continue;
487 : }
488 30 : if (std::abs(data.R) > 0.0)
489 : {
490 44 : for (const auto &elem : data.elems)
491 : {
492 : attr_list.Append(elem->GetAttrList());
493 : }
494 : }
495 : }
496 17 : return attr_list;
497 : }
498 :
499 17 : mfem::Array<int> LumpedPortOperator::GetLsAttrList() const
500 : {
501 : mfem::Array<int> attr_list;
502 47 : for (const auto &[idx, data] : ports)
503 : {
504 30 : if (!data.active)
505 : {
506 0 : continue;
507 : }
508 30 : if (std::abs(data.L) > 0.0)
509 : {
510 16 : for (const auto &elem : data.elems)
511 : {
512 : attr_list.Append(elem->GetAttrList());
513 : }
514 : }
515 : }
516 17 : return attr_list;
517 : }
518 :
519 0 : mfem::Array<int> LumpedPortOperator::GetCsAttrList() const
520 : {
521 : mfem::Array<int> attr_list;
522 0 : for (const auto &[idx, data] : ports)
523 : {
524 0 : if (!data.active)
525 : {
526 0 : continue;
527 : }
528 0 : if (std::abs(data.C) > 0.0)
529 : {
530 0 : for (const auto &elem : data.elems)
531 : {
532 : attr_list.Append(elem->GetAttrList());
533 : }
534 : }
535 : }
536 0 : return attr_list;
537 : }
538 :
539 0 : void LumpedPortOperator::AddStiffnessBdrCoefficients(double coeff,
540 : MaterialPropertyCoefficient &fb)
541 : {
542 : // Add lumped inductor boundaries to the bilinear form.
543 0 : for (const auto &[idx, data] : ports)
544 : {
545 0 : if (!data.active)
546 : {
547 0 : continue;
548 : }
549 0 : if (std::abs(data.L) > 0.0)
550 : {
551 0 : for (const auto &elem : data.elems)
552 : {
553 0 : const double Ls = data.L * data.GetToSquare(*elem);
554 0 : fb.AddMaterialProperty(data.mat_op.GetCeedBdrAttributes(elem->GetAttrList()),
555 0 : coeff / Ls);
556 : }
557 : }
558 : }
559 0 : }
560 :
561 0 : void LumpedPortOperator::AddDampingBdrCoefficients(double coeff,
562 : MaterialPropertyCoefficient &fb)
563 : {
564 : // Add lumped resistor boundaries to the bilinear form.
565 0 : for (const auto &[idx, data] : ports)
566 : {
567 0 : if (!data.active)
568 : {
569 0 : continue;
570 : }
571 0 : if (std::abs(data.R) > 0.0)
572 : {
573 0 : for (const auto &elem : data.elems)
574 : {
575 0 : const double Rs = data.R * data.GetToSquare(*elem);
576 0 : fb.AddMaterialProperty(data.mat_op.GetCeedBdrAttributes(elem->GetAttrList()),
577 0 : coeff / Rs);
578 : }
579 : }
580 : }
581 0 : }
582 :
583 0 : void LumpedPortOperator::AddMassBdrCoefficients(double coeff,
584 : MaterialPropertyCoefficient &fb)
585 : {
586 : // Add lumped capacitance boundaries to the bilinear form.
587 0 : for (const auto &[idx, data] : ports)
588 : {
589 0 : if (!data.active)
590 : {
591 0 : continue;
592 : }
593 0 : if (std::abs(data.C) > 0.0)
594 : {
595 0 : for (const auto &elem : data.elems)
596 : {
597 0 : const double Cs = data.C / data.GetToSquare(*elem);
598 0 : fb.AddMaterialProperty(data.mat_op.GetCeedBdrAttributes(elem->GetAttrList()),
599 0 : coeff * Cs);
600 : }
601 : }
602 : }
603 0 : }
604 :
605 0 : void LumpedPortOperator::AddExcitationBdrCoefficients(int excitation_idx,
606 : SumVectorCoefficient &fb)
607 : {
608 : // Construct the RHS source term for lumped port boundaries, which looks like -U_inc =
609 : // +2 iω/Z_s E_inc for a port boundary with an incident field E_inc. The chosen incident
610 : // field magnitude corresponds to a unit incident power over the full port boundary. See
611 : // p. 49 and p. 82 of the COMSOL RF Module manual for more detail.
612 : // Note: The real RHS returned here does not yet have the factor of (iω) included, so
613 : // works for time domain simulations requiring RHS -U_inc(t).
614 0 : for (const auto &[idx, data] : ports)
615 : {
616 0 : if (data.excitation != excitation_idx)
617 : {
618 0 : continue;
619 : }
620 0 : MFEM_VERIFY(std::abs(data.R) > 0.0,
621 : "Unexpected zero resistance in excited lumped port!");
622 0 : for (const auto &elem : data.elems)
623 : {
624 0 : const double Rs = data.R * data.GetToSquare(*elem);
625 0 : const double Hinc = 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
626 0 : elem->GetGeometryLength() * data.elems.size());
627 0 : fb.AddCoefficient(elem->GetModeCoefficient(2.0 * Hinc));
628 : }
629 : }
630 0 : }
631 :
632 : } // namespace palace
|
