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 "surfacepostoperator.hpp"
5 :
6 : #include <complex>
7 : #include <set>
8 : #include "fem/gridfunction.hpp"
9 : #include "fem/integrator.hpp"
10 : #include "linalg/vector.hpp"
11 : #include "models/materialoperator.hpp"
12 : #include "models/strattonchu.hpp"
13 : #include "utils/communication.hpp"
14 : #include "utils/geodata.hpp"
15 : #include "utils/iodata.hpp"
16 : #include "utils/prettyprint.hpp"
17 : #include "utils/timer.hpp"
18 :
19 : namespace palace
20 : {
21 :
22 : namespace
23 : {
24 :
25 : template <typename T>
26 17 : mfem::Array<int> SetUpBoundaryProperties(const T &data,
27 : const mfem::Array<int> &bdr_attr_marker)
28 : {
29 : mfem::Array<int> attr_list;
30 17 : attr_list.Reserve(static_cast<int>(data.attributes.size()));
31 : std::set<int> bdr_warn_list;
32 29 : for (auto attr : data.attributes)
33 : {
34 : // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
35 : // "Boundary postprocessing attribute tags must be non-negative and "
36 : // "correspond to attributes in the mesh!");
37 : // MFEM_VERIFY(bdr_attr_marker[attr - 1],
38 : // "Unknown boundary postprocessing attribute " << attr << "!");
39 12 : if (attr <= 0 || attr > bdr_attr_marker.Size() || !bdr_attr_marker[attr - 1])
40 : {
41 : bdr_warn_list.insert(attr);
42 : }
43 : else
44 : {
45 12 : attr_list.Append(attr);
46 : }
47 : }
48 17 : if (!bdr_warn_list.empty())
49 : {
50 0 : Mpi::Print("\n");
51 0 : Mpi::Warning(
52 : "Unknown boundary postprocessing attributes!\nSolver will just ignore them!");
53 0 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
54 0 : Mpi::Print("\n");
55 : }
56 17 : return attr_list;
57 : }
58 :
59 : } // namespace
60 :
61 0 : SurfacePostOperator::SurfaceFluxData::SurfaceFluxData(
62 : const config::SurfaceFluxData &data, const mfem::ParMesh &mesh,
63 0 : const mfem::Array<int> &bdr_attr_marker)
64 : {
65 : // Store boundary attributes for this postprocessing boundary.
66 0 : attr_list = SetUpBoundaryProperties(data, bdr_attr_marker);
67 :
68 : // Store the type of flux.
69 0 : switch (data.type)
70 : {
71 0 : case SurfaceFlux::ELECTRIC:
72 0 : type = SurfaceFlux::ELECTRIC;
73 0 : break;
74 0 : case SurfaceFlux::MAGNETIC:
75 0 : type = SurfaceFlux::MAGNETIC;
76 0 : break;
77 0 : case SurfaceFlux::POWER:
78 0 : type = SurfaceFlux::POWER;
79 0 : break;
80 : }
81 :
82 : // Store information about the global direction for orientation. Note the true boundary
83 : // normal is used in calculating the flux, this is just used to determine the sign.
84 0 : two_sided = data.two_sided;
85 0 : if (!two_sided)
86 : {
87 0 : center.SetSize(mesh.SpaceDimension());
88 0 : if (data.no_center)
89 : {
90 : // Compute the center as the bounding box centroid for all boundary elements making up
91 : // this postprocessing boundary.
92 : mfem::Vector bbmin, bbmax;
93 0 : mesh::GetAxisAlignedBoundingBox(
94 0 : mesh, mesh::AttrToMarker(bdr_attr_marker.Size(), attr_list), true, bbmin, bbmax);
95 0 : for (int d = 0; d < mesh.SpaceDimension(); d++)
96 : {
97 0 : center(d) = 0.5 * (bbmin(d) + bbmax(d));
98 : }
99 : }
100 : else
101 : {
102 : std::copy(data.center.begin(), data.center.end(), center.begin());
103 : }
104 : }
105 0 : }
106 :
107 : std::unique_ptr<mfem::Coefficient>
108 0 : SurfacePostOperator::SurfaceFluxData::GetCoefficient(const mfem::ParGridFunction *E,
109 : const mfem::ParGridFunction *B,
110 : const MaterialOperator &mat_op) const
111 : {
112 0 : switch (type)
113 : {
114 0 : case SurfaceFlux::ELECTRIC:
115 : return std::make_unique<
116 0 : RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFlux::ELECTRIC>>>(
117 0 : attr_list, E, nullptr, mat_op, two_sided, center);
118 0 : case SurfaceFlux::MAGNETIC:
119 : return std::make_unique<
120 0 : RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFlux::MAGNETIC>>>(
121 0 : attr_list, nullptr, B, mat_op, two_sided, center);
122 0 : case SurfaceFlux::POWER:
123 : return std::make_unique<
124 0 : RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFlux::POWER>>>(
125 0 : attr_list, E, B, mat_op, two_sided, center);
126 : }
127 : return {};
128 : }
129 :
130 0 : SurfacePostOperator::InterfaceDielectricData::InterfaceDielectricData(
131 : const config::InterfaceDielectricData &data, const mfem::ParMesh &mesh,
132 0 : const mfem::Array<int> &bdr_attr_marker)
133 : {
134 : // Store boundary attributes for this postprocessing boundary.
135 0 : attr_list = SetUpBoundaryProperties(data, bdr_attr_marker);
136 :
137 : // Calculate surface dielectric loss according to the formulas from J. Wenner et al.,
138 : // Surface loss simulations of superconducting coplanar waveguide resonators, Appl. Phys.
139 : // Lett. (2011). If only a general layer permittivity is specified and not any special
140 : // metal-air (MA), metal-substrate (MS), or substrate-air (SA) permittivity, compute the
141 : // numerator of the participation ratio according to the regular formula
142 : // p * E_elec = 1/2 t Re{∫ (ε E)ᴴ E_m dS} .
143 0 : switch (data.type)
144 : {
145 0 : case InterfaceDielectric::DEFAULT:
146 0 : type = InterfaceDielectric::DEFAULT;
147 0 : break;
148 0 : case InterfaceDielectric::MA:
149 0 : type = InterfaceDielectric::MA;
150 0 : break;
151 0 : case InterfaceDielectric::MS:
152 0 : type = InterfaceDielectric::MS;
153 0 : break;
154 0 : case InterfaceDielectric::SA:
155 0 : type = InterfaceDielectric::SA;
156 0 : break;
157 : }
158 0 : t = data.t;
159 0 : epsilon = data.epsilon_r;
160 0 : tandelta = data.tandelta;
161 0 : }
162 :
163 : std::unique_ptr<mfem::Coefficient>
164 0 : SurfacePostOperator::InterfaceDielectricData::GetCoefficient(
165 : const GridFunction &E, const MaterialOperator &mat_op) const
166 : {
167 0 : switch (type)
168 : {
169 0 : case InterfaceDielectric::DEFAULT:
170 : return std::make_unique<RestrictedCoefficient<
171 0 : InterfaceDielectricCoefficient<InterfaceDielectric::DEFAULT>>>(
172 0 : attr_list, E, mat_op, t, epsilon);
173 0 : case InterfaceDielectric::MA:
174 : return std::make_unique<
175 0 : RestrictedCoefficient<InterfaceDielectricCoefficient<InterfaceDielectric::MA>>>(
176 0 : attr_list, E, mat_op, t, epsilon);
177 0 : case InterfaceDielectric::MS:
178 : return std::make_unique<
179 0 : RestrictedCoefficient<InterfaceDielectricCoefficient<InterfaceDielectric::MS>>>(
180 0 : attr_list, E, mat_op, t, epsilon);
181 0 : case InterfaceDielectric::SA:
182 : return std::make_unique<
183 0 : RestrictedCoefficient<InterfaceDielectricCoefficient<InterfaceDielectric::SA>>>(
184 0 : attr_list, E, mat_op, t, epsilon);
185 : }
186 : return {}; // For compiler warning
187 : }
188 :
189 17 : SurfacePostOperator::FarFieldData::FarFieldData(const config::FarFieldPostData &data,
190 : const mfem::ParMesh &mesh,
191 17 : const mfem::Array<int> &bdr_attr_marker)
192 17 : : thetaphis(data.thetaphis)
193 : {
194 : // Store boundary attributes for this postprocessing boundary.
195 17 : attr_list = SetUpBoundaryProperties(data, bdr_attr_marker);
196 17 : }
197 :
198 18 : SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
199 : const MaterialOperator &mat_op,
200 : mfem::ParFiniteElementSpace &h1_fespace,
201 18 : mfem::ParFiniteElementSpace &nd_fespace)
202 18 : : mat_op(mat_op), h1_fespace(h1_fespace), nd_fespace(nd_fespace)
203 : {
204 : // Check that boundary attributes have been specified correctly.
205 : const auto &mesh = *h1_fespace.GetParMesh();
206 18 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
207 : mfem::Array<int> bdr_attr_marker;
208 18 : if (!iodata.boundaries.postpro.flux.empty() ||
209 36 : !iodata.boundaries.postpro.dielectric.empty() ||
210 : !iodata.boundaries.postpro.farfield.empty())
211 : {
212 : bdr_attr_marker.SetSize(bdr_attr_max);
213 : bdr_attr_marker = 0;
214 21 : for (auto attr : mesh.bdr_attributes)
215 : {
216 18 : bdr_attr_marker[attr - 1] = 1;
217 : }
218 : }
219 :
220 : // Surface flux postprocessing.
221 18 : for (const auto &[idx, data] : iodata.boundaries.postpro.flux)
222 : {
223 0 : MFEM_VERIFY(iodata.problem.type != ProblemType::ELECTROSTATIC ||
224 : data.type == SurfaceFlux::ELECTRIC,
225 : "Magnetic field or power surface flux postprocessing are not available "
226 : "for electrostatic problems!");
227 0 : MFEM_VERIFY(iodata.problem.type != ProblemType::MAGNETOSTATIC ||
228 : data.type == SurfaceFlux::MAGNETIC,
229 : "Electric field or power surface flux postprocessing are not available "
230 : "for magnetostatic problems!");
231 0 : flux_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
232 : }
233 :
234 : // Interface dielectric postprocessing.
235 18 : MFEM_VERIFY(iodata.boundaries.postpro.dielectric.empty() ||
236 : iodata.problem.type != ProblemType::MAGNETOSTATIC,
237 : "Interface dielectric loss postprocessing is not available for "
238 : "magnetostatic problems!");
239 18 : for (const auto &[idx, data] : iodata.boundaries.postpro.dielectric)
240 : {
241 0 : eps_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
242 : }
243 :
244 : // FarField postprocessing.
245 18 : MFEM_VERIFY(iodata.boundaries.postpro.farfield.empty() ||
246 : iodata.problem.type == ProblemType::DRIVEN ||
247 : iodata.problem.type == ProblemType::EIGENMODE,
248 : "Far-field extraction is only available for driven and eigenmode problems!");
249 :
250 : // Check that we don't have anisotropic materials.
251 18 : if (!iodata.boundaries.postpro.farfield.empty())
252 : {
253 : const auto &mesh = *nd_fespace.GetParMesh();
254 3 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
255 : mfem::Array<int> bdr_attr_marker =
256 3 : mesh::AttrToMarker(bdr_attr_max, iodata.boundaries.postpro.farfield.attributes);
257 :
258 : std::set<int> domain_attrs;
259 :
260 9651 : for (int i = 0; i < mesh.GetNBE(); i++)
261 : {
262 9648 : if (bdr_attr_marker[mesh.GetBdrAttribute(i) - 1])
263 : {
264 : int elem_id, _face_id;
265 9608 : mesh.GetBdrElementAdjacentElement(i, elem_id, _face_id);
266 9608 : if (elem_id >= 0)
267 : {
268 9609 : domain_attrs.insert(mesh.GetAttribute(elem_id));
269 : }
270 : }
271 : }
272 :
273 5 : for (int attr : domain_attrs)
274 : {
275 7 : MFEM_VERIFY(mat_op.IsIsotropic(attr),
276 : "FarField requires isotropic materials, but attribute " +
277 : std::to_string(attr) + " is not.");
278 : }
279 : }
280 :
281 18 : farfield = FarFieldData(iodata.boundaries.postpro.farfield, *nd_fespace.GetParMesh(),
282 17 : bdr_attr_marker);
283 18 : }
284 :
285 0 : std::complex<double> SurfacePostOperator::GetSurfaceFlux(int idx, const GridFunction *E,
286 : const GridFunction *B) const
287 : {
288 : // For complex-valued fields, output the separate real and imaginary parts for the time-
289 : // harmonic quantity. For power flux (Poynting vector), output only the stationary real
290 : // part and not the part which has double the frequency.
291 : auto it = flux_surfs.find(idx);
292 0 : MFEM_VERIFY(it != flux_surfs.end(),
293 : "Unknown surface flux postprocessing index requested!");
294 0 : const bool has_imag = (E) ? E->HasImag() : B->HasImag();
295 0 : const auto &mesh = *h1_fespace.GetParMesh();
296 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
297 0 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
298 : auto f =
299 0 : it->second.GetCoefficient(E ? &E->Real() : nullptr, B ? &B->Real() : nullptr, mat_op);
300 0 : std::complex<double> dot(GetLocalSurfaceIntegral(*f, attr_marker), 0.0);
301 0 : if (has_imag)
302 : {
303 0 : f = it->second.GetCoefficient(E ? &E->Imag() : nullptr, B ? &B->Imag() : nullptr,
304 : mat_op);
305 0 : double doti = GetLocalSurfaceIntegral(*f, attr_marker);
306 0 : if (it->second.type == SurfaceFlux::POWER)
307 : {
308 : dot += doti;
309 : }
310 : else
311 : {
312 : dot.imag(doti);
313 : }
314 : }
315 0 : Mpi::GlobalSum(1, &dot, (E) ? E->GetComm() : B->GetComm());
316 0 : return dot;
317 : }
318 :
319 0 : double SurfacePostOperator::GetInterfaceLossTangent(int idx) const
320 : {
321 : auto it = eps_surfs.find(idx);
322 0 : MFEM_VERIFY(it != eps_surfs.end(),
323 : "Unknown interface dielectric postprocessing index requested!");
324 0 : return it->second.tandelta;
325 : }
326 :
327 0 : double SurfacePostOperator::GetInterfaceElectricFieldEnergy(int idx,
328 : const GridFunction &E) const
329 : {
330 : auto it = eps_surfs.find(idx);
331 0 : MFEM_VERIFY(it != eps_surfs.end(),
332 : "Unknown interface dielectric postprocessing index requested!");
333 0 : const auto &mesh = *h1_fespace.GetParMesh();
334 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
335 0 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
336 0 : auto f = it->second.GetCoefficient(E, mat_op);
337 0 : double dot = GetLocalSurfaceIntegral(*f, attr_marker);
338 : Mpi::GlobalSum(1, &dot, E.GetComm());
339 0 : return dot;
340 : }
341 :
342 : double
343 0 : SurfacePostOperator::GetLocalSurfaceIntegral(mfem::Coefficient &f,
344 : const mfem::Array<int> &attr_marker) const
345 : {
346 : // Integrate the coefficient over the boundary attributes making up this surface index.
347 0 : mfem::LinearForm s(&h1_fespace);
348 0 : s.AddBoundaryIntegrator(new BoundaryLFIntegrator(f),
349 : const_cast<mfem::Array<int> &>(attr_marker));
350 0 : s.UseFastAssembly(false);
351 : s.UseDevice(false);
352 0 : s.Assemble();
353 : s.UseDevice(true);
354 0 : return linalg::LocalSum(s);
355 0 : }
356 :
357 8 : std::vector<std::array<std::complex<double>, 3>> SurfacePostOperator::GetFarFieldrE(
358 : const std::vector<std::pair<double, double>> &theta_phi_pairs, const GridFunction &E,
359 : const GridFunction &B, double omega_re, double omega_im) const
360 : {
361 8 : if (theta_phi_pairs.empty())
362 6 : return {};
363 2 : BlockTimer bt0(Timer::POSTPRO_FARFIELD);
364 :
365 : // Compute target unit vectors from the given theta and phis.
366 : std::vector<std::array<double, 3>> r_naughts;
367 2 : r_naughts.reserve(theta_phi_pairs.size());
368 :
369 2 : r_naughts.reserve(theta_phi_pairs.size());
370 2050 : for (const auto &[theta, phi] : theta_phi_pairs)
371 : {
372 2048 : r_naughts.emplace_back(std::array<double, 3>{
373 2048 : std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta)});
374 : }
375 :
376 2 : const auto &mesh = *nd_fespace.GetParMesh();
377 2 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
378 2 : mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, farfield.attr_list);
379 :
380 : // Integrate. Each MPI process computes its contribution and we will reduce
381 : // everything at the end. We make them std::vector<std::array<double, 3>>
382 : // because we want a very simple memory layout so that we can reduce
383 : // everything with two MPI calls.
384 2 : std::vector<std::array<double, 3>> integrals_r(theta_phi_pairs.size());
385 2 : std::vector<std::array<double, 3>> integrals_i(theta_phi_pairs.size());
386 :
387 9602 : for (int i = 0; i < mesh.GetNBE(); i++)
388 : {
389 9600 : if (!attr_marker[mesh.GetBdrAttribute(i) - 1])
390 0 : continue;
391 :
392 9600 : auto *T = const_cast<mfem::ParMesh &>(mesh).GetBdrElementTransformation(i);
393 9600 : const auto *fe = nd_fespace.GetBE(i);
394 : const auto *ir =
395 9600 : &mfem::IntRules.Get(fe->GetGeomType(), fem::DefaultIntegrationOrder::Get(*T));
396 :
397 9600 : AddStrattonChuIntegrandAtElement(E, B, mat_op, omega_re, omega_im, r_naughts, *T, *ir,
398 : integrals_r, integrals_i);
399 : }
400 :
401 : double *data_r_ptr = integrals_r.data()->data();
402 : double *data_i_ptr = integrals_i.data()->data();
403 2 : size_t total_elements = integrals_r.size() * 3;
404 2 : Mpi::GlobalSum(total_elements, data_i_ptr, E.GetComm());
405 : Mpi::GlobalSum(total_elements, data_r_ptr, E.GetComm());
406 :
407 : // Finally, we apply cross product to reduced integrals and package the result
408 : // in a neatly accessible vector of arrays of complex numbers.
409 2 : std::vector<std::array<std::complex<double>, 3>> result(theta_phi_pairs.size());
410 2 : StaticVector<3> tmp_r, tmp_i;
411 2050 : for (size_t k = 0; k < theta_phi_pairs.size(); k++)
412 : {
413 2048 : linalg::Cross3(r_naughts[k], integrals_r[k], tmp_r);
414 2048 : linalg::Cross3(r_naughts[k], integrals_i[k], tmp_i);
415 8192 : for (size_t d = 0; d < 3; d++)
416 : {
417 6144 : result[k][d] = std::complex<double>{tmp_r[d], tmp_i[d]};
418 : }
419 : }
420 : return result;
421 2 : }
422 :
423 : } // namespace palace
|
