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 "postoperator.hpp"
5 :
6 : #include <algorithm>
7 : #include <string>
8 : #include "fem/coefficient.hpp"
9 : #include "fem/errorindicator.hpp"
10 : #include "models/curlcurloperator.hpp"
11 : #include "models/laplaceoperator.hpp"
12 : #include "models/materialoperator.hpp"
13 : #include "models/spaceoperator.hpp"
14 : #include "models/surfacecurrentoperator.hpp"
15 : #include "models/waveportoperator.hpp"
16 : #include "utils/communication.hpp"
17 : #include "utils/geodata.hpp"
18 : #include "utils/iodata.hpp"
19 : #include "utils/timer.hpp"
20 :
21 : namespace palace
22 : {
23 :
24 : namespace
25 : {
26 :
27 : std::string OutputFolderName(const ProblemType solver_t)
28 : {
29 : switch (solver_t)
30 : {
31 : case ProblemType::DRIVEN:
32 9 : return "driven";
33 : case ProblemType::EIGENMODE:
34 6 : return "eigenmode";
35 : case ProblemType::ELECTROSTATIC:
36 6 : return "electrostatic";
37 : case ProblemType::MAGNETOSTATIC:
38 6 : return "magnetostatic";
39 : case ProblemType::TRANSIENT:
40 6 : return "transient";
41 : default:
42 : return "unknown";
43 : }
44 : }
45 :
46 : } // namespace
47 :
48 : template <ProblemType solver_t>
49 15 : PostOperator<solver_t>::PostOperator(const IoData &iodata, fem_op_t<solver_t> &fem_op_)
50 15 : : fem_op(&fem_op_), units(iodata.units), post_dir(iodata.problem.output),
51 15 : post_op_csv(iodata, fem_op_),
52 : // dom_post_op does not have a default ctor so specialize via immediate lambda.
53 15 : dom_post_op(std::move(
54 : [&iodata, &fem_op_]()
55 : {
56 : if constexpr (solver_t == ProblemType::ELECTROSTATIC)
57 : {
58 : return DomainPostOperator(iodata, fem_op_.GetMaterialOp(),
59 3 : fem_op_.GetH1Space());
60 : }
61 : else if constexpr (solver_t == ProblemType::MAGNETOSTATIC)
62 : {
63 : return DomainPostOperator(iodata, fem_op_.GetMaterialOp(),
64 3 : fem_op_.GetNDSpace());
65 : }
66 : else
67 : {
68 : return DomainPostOperator(iodata, fem_op_.GetMaterialOp(), fem_op_.GetNDSpace(),
69 9 : fem_op_.GetRTSpace());
70 : }
71 15 : }())),
72 15 : surf_post_op(iodata, fem_op->GetMaterialOp(), fem_op->GetH1Space(),
73 15 : fem_op->GetNDSpace()),
74 30 : interp_op(iodata, fem_op->GetNDSpace())
75 : {
76 : // Define primary grid-functions.
77 : if constexpr (HasVGridFunction<solver_t>())
78 : {
79 3 : V = std::make_unique<GridFunction>(fem_op->GetH1Space());
80 : }
81 : if constexpr (HasAGridFunction<solver_t>())
82 : {
83 3 : A = std::make_unique<GridFunction>(fem_op->GetNDSpace());
84 : }
85 : if constexpr (HasEGridFunction<solver_t>())
86 : {
87 24 : E = std::make_unique<GridFunction>(fem_op->GetNDSpace(),
88 12 : HasComplexGridFunction<solver_t>());
89 : }
90 : if constexpr (HasBGridFunction<solver_t>())
91 : {
92 24 : B = std::make_unique<GridFunction>(fem_op->GetRTSpace(),
93 12 : HasComplexGridFunction<solver_t>());
94 : }
95 :
96 : // Add wave port boundary mode postprocessing, if available.
97 : if constexpr (std::is_same_v<fem_op_t<solver_t>, SpaceOperator>)
98 : {
99 9 : for (const auto &[idx, data] : fem_op->GetWavePortOp())
100 : {
101 0 : auto ret = port_E0.emplace(idx, WavePortFieldData());
102 0 : ret.first->second.E0r = data.GetModeFieldCoefficientReal();
103 0 : ret.first->second.E0i = data.GetModeFieldCoefficientImag();
104 : }
105 : }
106 :
107 : // Prepare for saving fields
108 15 : enable_paraview_output = iodata.problem.output_formats.paraview;
109 15 : enable_gridfunction_output = iodata.problem.output_formats.gridfunction;
110 : if (solver_t == ProblemType::DRIVEN)
111 : {
112 3 : output_save_indices = iodata.solver.driven.save_indices;
113 : }
114 : else if (solver_t == ProblemType::EIGENMODE)
115 : {
116 3 : output_n_post = iodata.solver.eigenmode.n_post;
117 : }
118 : else if (solver_t == ProblemType::ELECTROSTATIC)
119 : {
120 3 : output_n_post = iodata.solver.electrostatic.n_post;
121 : }
122 : else if (solver_t == ProblemType::MAGNETOSTATIC)
123 : {
124 3 : output_n_post = iodata.solver.magnetostatic.n_post;
125 : }
126 : else if (solver_t == ProblemType::TRANSIENT)
127 : {
128 3 : output_delta_post = iodata.solver.transient.delta_post;
129 : }
130 :
131 : gridfunction_output_dir =
132 45 : (post_dir / "gridfunction" / OutputFolderName(solver_t)).string();
133 :
134 15 : SetupFieldCoefficients();
135 15 : InitializeParaviewDataCollection();
136 :
137 : // Initialize CSV files for measurements.
138 15 : post_op_csv.InitializeCSVDataCollection(*this);
139 15 : }
140 :
141 : template <ProblemType solver_t>
142 : template <ProblemType U>
143 0 : auto PostOperator<solver_t>::InitializeParaviewDataCollection(int ex_idx)
144 : -> std::enable_if_t<U == ProblemType::DRIVEN, void>
145 : {
146 0 : fs::path sub_folder_name = "";
147 0 : auto nr_excitations = fem_op->GetPortExcitations().Size();
148 0 : if ((nr_excitations > 1) && (ex_idx > 0))
149 : {
150 0 : int spacing = 1 + int(std::log10(nr_excitations));
151 0 : sub_folder_name = fmt::format(FMT_STRING("excitation_{:0>{}}"), ex_idx, spacing);
152 : }
153 0 : InitializeParaviewDataCollection(sub_folder_name);
154 0 : }
155 :
156 : template <ProblemType solver_t>
157 15 : void PostOperator<solver_t>::SetupFieldCoefficients()
158 : {
159 : // We currently don't use the dependent grid functions apart from saving fields, so only
160 : // initialize if needed.
161 : if (!ShouldWriteFields())
162 : {
163 : return;
164 : }
165 :
166 : // Set-up grid-functions for the paraview output / measurement.
167 : if constexpr (HasVGridFunction<solver_t>())
168 : {
169 3 : V_s = std::make_unique<BdrFieldCoefficient>(V->Real());
170 : }
171 :
172 : if constexpr (HasAGridFunction<solver_t>())
173 : {
174 3 : A_s = std::make_unique<BdrFieldVectorCoefficient>(A->Real());
175 : }
176 :
177 : if constexpr (HasEGridFunction<solver_t>())
178 : {
179 : // Electric Energy Density.
180 24 : U_e = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>>(
181 12 : *E, fem_op->GetMaterialOp());
182 :
183 : // Electric Boundary Field & Surface Charge.
184 24 : E_sr = std::make_unique<BdrFieldVectorCoefficient>(E->Real());
185 12 : Q_sr = std::make_unique<BdrSurfaceFluxCoefficient<SurfaceFlux::ELECTRIC>>(
186 24 : &E->Real(), nullptr, fem_op->GetMaterialOp(), true, mfem::Vector());
187 :
188 : if constexpr (HasComplexGridFunction<solver_t>())
189 : {
190 12 : E_si = std::make_unique<BdrFieldVectorCoefficient>(E->Imag());
191 6 : Q_si = std::make_unique<BdrSurfaceFluxCoefficient<SurfaceFlux::ELECTRIC>>(
192 12 : &E->Imag(), nullptr, fem_op->GetMaterialOp(), true, mfem::Vector());
193 : }
194 : }
195 :
196 : if constexpr (HasBGridFunction<solver_t>())
197 : {
198 : // Magnetic Energy Density.
199 24 : U_m = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>>(
200 12 : *B, fem_op->GetMaterialOp());
201 :
202 : // Magnetic Boundary Field & Surface Current.
203 12 : B_sr = std::make_unique<BdrFieldVectorCoefficient>(B->Real());
204 21 : J_sr = std::make_unique<BdrSurfaceCurrentVectorCoefficient>(B->Real(),
205 12 : fem_op->GetMaterialOp());
206 :
207 : if constexpr (HasComplexGridFunction<solver_t>())
208 : {
209 6 : B_si = std::make_unique<BdrFieldVectorCoefficient>(B->Imag());
210 6 : J_si = std::make_unique<BdrSurfaceCurrentVectorCoefficient>(B->Imag(),
211 6 : fem_op->GetMaterialOp());
212 : }
213 : }
214 :
215 : if constexpr (HasEGridFunction<solver_t>() && HasBGridFunction<solver_t>())
216 : {
217 : // Poynting Vector.
218 18 : S = std::make_unique<PoyntingVectorCoefficient>(*E, *B, fem_op->GetMaterialOp());
219 : }
220 : }
221 :
222 : template <ProblemType solver_t>
223 15 : void PostOperator<solver_t>::InitializeParaviewDataCollection(
224 : const fs::path &sub_folder_name)
225 : {
226 : if (!ShouldWriteParaviewFields())
227 : {
228 0 : return;
229 : }
230 30 : fs::path paraview_dir_v = post_dir / "paraview" / OutputFolderName(solver_t);
231 15 : fs::path paraview_dir_b =
232 45 : post_dir / "paraview" / fmt::format("{}_boundary", OutputFolderName(solver_t));
233 15 : if (!sub_folder_name.empty())
234 : {
235 0 : paraview_dir_v /= sub_folder_name;
236 0 : paraview_dir_b /= sub_folder_name;
237 : }
238 : // Set up postprocessing for output to disk.
239 30 : paraview = {paraview_dir_v, &fem_op->GetNDSpace().GetParMesh()};
240 30 : paraview_bdr = {paraview_dir_b, &fem_op->GetNDSpace().GetParMesh()};
241 :
242 : const mfem::VTKFormat format = mfem::VTKFormat::BINARY32;
243 : #if defined(MFEM_USE_ZLIB)
244 : const int compress = -1; // Default compression level
245 : #else
246 : const int compress = 0;
247 : #endif
248 : const bool use_ho = true;
249 : const int refine_ho = HasEGridFunction<solver_t>()
250 15 : ? E->ParFESpace()->GetMaxElementOrder()
251 3 : : B->ParFESpace()->GetMaxElementOrder();
252 :
253 : // Output mesh coordinate units same as input.
254 : paraview->SetCycle(-1);
255 15 : paraview->SetDataFormat(format);
256 15 : paraview->SetCompressionLevel(compress);
257 15 : paraview->SetHighOrderOutput(use_ho);
258 15 : paraview->SetLevelsOfDetail(refine_ho);
259 :
260 15 : paraview_bdr->SetBoundaryOutput(true);
261 : paraview_bdr->SetCycle(-1);
262 15 : paraview_bdr->SetDataFormat(format);
263 15 : paraview_bdr->SetCompressionLevel(compress);
264 15 : paraview_bdr->SetHighOrderOutput(use_ho);
265 15 : paraview_bdr->SetLevelsOfDetail(refine_ho);
266 :
267 : // Output fields @ phase = 0 and π/2 for frequency domain (rather than, for example,
268 : // peak phasors or magnitude = sqrt(2) * RMS). Also output fields evaluated on mesh
269 : // boundaries. For internal boundary surfaces, this takes the field evaluated in the
270 : // neighboring element with the larger dielectric permittivity or magnetic
271 : // permeability.
272 15 : if (E)
273 : {
274 : if (HasComplexGridFunction<solver_t>())
275 : {
276 6 : paraview->RegisterField("E_real", &E->Real());
277 12 : paraview->RegisterField("E_imag", &E->Imag());
278 12 : paraview_bdr->RegisterVCoeffField("E_real", E_sr.get());
279 12 : paraview_bdr->RegisterVCoeffField("E_imag", E_si.get());
280 : }
281 : else
282 : {
283 12 : paraview->RegisterField("E", &E->Real());
284 12 : paraview_bdr->RegisterVCoeffField("E", E_sr.get());
285 : }
286 : }
287 15 : if (B)
288 : {
289 : if (HasComplexGridFunction<solver_t>())
290 : {
291 6 : paraview->RegisterField("B_real", &B->Real());
292 12 : paraview->RegisterField("B_imag", &B->Imag());
293 12 : paraview_bdr->RegisterVCoeffField("B_real", B_sr.get());
294 12 : paraview_bdr->RegisterVCoeffField("B_imag", B_si.get());
295 : }
296 : else
297 : {
298 12 : paraview->RegisterField("B", &B->Real());
299 12 : paraview_bdr->RegisterVCoeffField("B", B_sr.get());
300 : }
301 : }
302 15 : if (V)
303 : {
304 6 : paraview->RegisterField("V", &V->Real());
305 6 : paraview_bdr->RegisterCoeffField("V", V_s.get());
306 : }
307 15 : if (A)
308 : {
309 6 : paraview->RegisterField("A", &A->Real());
310 6 : paraview_bdr->RegisterVCoeffField("A", A_s.get());
311 : }
312 :
313 : // Extract energy density field for electric field energy 1/2 Dᴴ E or magnetic field
314 : // energy 1/2 Hᴴ B. Also Poynting vector S = E x H⋆.
315 15 : if (U_e)
316 : {
317 24 : paraview->RegisterCoeffField("U_e", U_e.get());
318 24 : paraview_bdr->RegisterCoeffField("U_e", U_e.get());
319 : }
320 15 : if (U_m)
321 : {
322 24 : paraview->RegisterCoeffField("U_m", U_m.get());
323 24 : paraview_bdr->RegisterCoeffField("U_m", U_m.get());
324 : }
325 15 : if (S)
326 : {
327 18 : paraview->RegisterVCoeffField("S", S.get());
328 18 : paraview_bdr->RegisterVCoeffField("S", S.get());
329 : }
330 :
331 : // Extract surface charge from normally discontinuous ND E-field. Also extract surface
332 : // currents from tangentially discontinuous RT B-field The surface charge and surface
333 : // currents are single-valued at internal boundaries.
334 15 : if (Q_sr)
335 : {
336 : if (HasComplexGridFunction<solver_t>())
337 : {
338 12 : paraview_bdr->RegisterCoeffField("Q_s_real", Q_sr.get());
339 12 : paraview_bdr->RegisterCoeffField("Q_s_imag", Q_si.get());
340 : }
341 : else
342 : {
343 12 : paraview_bdr->RegisterCoeffField("Q_s", Q_sr.get());
344 : }
345 : }
346 15 : if (J_sr)
347 : {
348 : if (HasComplexGridFunction<solver_t>())
349 : {
350 12 : paraview_bdr->RegisterVCoeffField("J_s_real", J_sr.get());
351 12 : paraview_bdr->RegisterVCoeffField("J_s_imag", J_si.get());
352 : }
353 : else
354 : {
355 12 : paraview_bdr->RegisterVCoeffField("J_s", J_sr.get());
356 : }
357 : }
358 :
359 : // Add wave port boundary mode postprocessing when available.
360 15 : for (const auto &[idx, data] : port_E0)
361 : {
362 0 : paraview_bdr->RegisterVCoeffField(fmt::format("E0_{}_real", idx), data.E0r.get());
363 0 : paraview_bdr->RegisterVCoeffField(fmt::format("E0_{}_imag", idx), data.E0i.get());
364 : }
365 15 : }
366 :
367 : namespace
368 : {
369 :
370 : template <typename T>
371 60 : void ScaleGridFunctions(double L, int dim, bool imag, T &E, T &B, T &V, T &A)
372 : {
373 : // For fields on H(curl) and H(div) spaces, we "undo" the effect of redimensionalizing
374 : // the mesh which would carry into the fields during the mapping from reference to
375 : // physical space through the element Jacobians. No transformation for V is needed (H1
376 : // interpolation). Because the coefficients are always evaluating E, B in neighboring
377 : // elements, the Jacobian scaling is the same for the domain and boundary data
378 : // collections (instead of being different for B due to the dim - 1 evaluation). Wave
379 : // port fields also do not require rescaling since their submesh object where they are
380 : // evaluated remains nondimensionalized.
381 60 : if (E)
382 : {
383 : // Piola transform: J^-T
384 48 : E->Real() *= L;
385 48 : E->Real().FaceNbrData() *= L;
386 48 : if (imag)
387 : {
388 24 : E->Imag() *= L;
389 24 : E->Imag().FaceNbrData() *= L;
390 : }
391 : }
392 60 : if (B)
393 : {
394 : // Piola transform: J / |J|
395 48 : const auto Ld = std::pow(L, dim - 1);
396 48 : B->Real() *= Ld;
397 48 : B->Real().FaceNbrData() *= Ld;
398 48 : if (imag)
399 : {
400 24 : B->Imag() *= Ld;
401 24 : B->Imag().FaceNbrData() *= Ld;
402 : }
403 : }
404 60 : if (A)
405 : {
406 : // Piola transform: J^-T
407 12 : A->Real() *= L;
408 12 : A->Real().FaceNbrData() *= L;
409 : }
410 60 : }
411 :
412 : } // namespace
413 :
414 : template <ProblemType solver_t>
415 15 : void PostOperator<solver_t>::WriteParaviewFields(double time, int step)
416 : {
417 15 : BlockTimer bt(Timer::POSTPRO_PARAVIEW);
418 :
419 : auto mesh_Lc0 = units.GetMeshLengthRelativeScale();
420 :
421 : // Given the electric field and magnetic flux density, write the fields to disk for
422 : // visualization. Write the mesh coordinates in the same units as originally input.
423 15 : mfem::ParMesh &mesh = E ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
424 15 : mesh::DimensionalizeMesh(mesh, mesh_Lc0);
425 15 : ScaleGridFunctions(mesh_Lc0, mesh.Dimension(), HasComplexGridFunction<solver_t>(), E, B,
426 15 : V, A);
427 : paraview->SetCycle(step);
428 : paraview->SetTime(time);
429 : paraview_bdr->SetCycle(step);
430 : paraview_bdr->SetTime(time);
431 15 : paraview->Save();
432 15 : paraview_bdr->Save();
433 15 : mesh::NondimensionalizeMesh(mesh, mesh_Lc0);
434 15 : ScaleGridFunctions(1.0 / mesh_Lc0, mesh.Dimension(), HasComplexGridFunction<solver_t>(),
435 : E, B, V, A);
436 15 : Mpi::Barrier(fem_op->GetComm());
437 15 : }
438 :
439 : template <ProblemType solver_t>
440 0 : void PostOperator<solver_t>::WriteParaviewFieldsFinal(const ErrorIndicator *indicator)
441 : {
442 0 : BlockTimer bt(Timer::POSTPRO_PARAVIEW);
443 :
444 : auto mesh_Lc0 = units.GetMeshLengthRelativeScale();
445 :
446 : // Write the mesh partitioning and (optionally) error indicators at the final step. No
447 : // need for these to be parallel objects, since the data is local to each process and
448 : // there isn't a need to ever access the element neighbors. We set the time to some
449 : // non-used value to make the step identifiable within the data collection.
450 0 : mfem::ParMesh &mesh = E ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
451 0 : mesh::DimensionalizeMesh(mesh, mesh_Lc0);
452 0 : paraview->SetCycle(paraview->GetCycle() + 1);
453 0 : if (paraview->GetTime() < 1.0)
454 : {
455 : paraview->SetTime(99.0);
456 : }
457 : else
458 : {
459 : // 1 -> 99, 10 -> 999, etc.
460 0 : paraview->SetTime(
461 0 : std::pow(10.0, 2.0 + static_cast<int>(std::log10(paraview->GetTime()))) - 1.0);
462 : }
463 : auto field_map = paraview->GetFieldMap(); // Copy, so can reregister later
464 0 : for (const auto &[name, gf] : field_map)
465 : {
466 0 : paraview->DeregisterField(name);
467 : }
468 : auto coeff_field_map = paraview->GetCoeffFieldMap();
469 0 : for (const auto &[name, gf] : coeff_field_map)
470 : {
471 : paraview->DeregisterCoeffField(name);
472 : }
473 : auto vcoeff_field_map = paraview->GetVCoeffFieldMap();
474 0 : for (const auto &[name, gf] : vcoeff_field_map)
475 : {
476 : paraview->DeregisterVCoeffField(name);
477 : }
478 0 : mfem::L2_FECollection pwconst_fec(0, mesh.Dimension());
479 0 : mfem::FiniteElementSpace pwconst_fespace(&mesh, &pwconst_fec);
480 : std::unique_ptr<mfem::GridFunction> rank, eta;
481 : {
482 0 : rank = std::make_unique<mfem::GridFunction>(&pwconst_fespace);
483 0 : *rank = mesh.GetMyRank() + 1;
484 0 : paraview->RegisterField("Rank", rank.get());
485 : }
486 0 : if (indicator)
487 : {
488 0 : eta = std::make_unique<mfem::GridFunction>(&pwconst_fespace);
489 0 : MFEM_VERIFY(eta->Size() == indicator->Local().Size(),
490 : "Size mismatch for provided ErrorIndicator for postprocessing!");
491 0 : *eta = indicator->Local();
492 0 : paraview->RegisterField("Indicator", eta.get());
493 : }
494 0 : paraview->Save();
495 : if (rank)
496 : {
497 0 : paraview->DeregisterField("Rank");
498 : }
499 0 : if (eta)
500 : {
501 0 : paraview->DeregisterField("Indicator");
502 : }
503 0 : for (const auto &[name, gf] : field_map)
504 : {
505 0 : paraview->RegisterField(name, gf);
506 : }
507 0 : for (const auto &[name, gf] : coeff_field_map)
508 : {
509 0 : paraview->RegisterCoeffField(name, gf);
510 : }
511 0 : for (const auto &[name, gf] : vcoeff_field_map)
512 : {
513 0 : paraview->RegisterVCoeffField(name, gf);
514 : }
515 0 : mesh::NondimensionalizeMesh(mesh, mesh_Lc0);
516 0 : Mpi::Barrier(fem_op->GetComm());
517 0 : }
518 :
519 : template <ProblemType solver_t>
520 15 : void PostOperator<solver_t>::WriteMFEMGridFunctions(double time, int step)
521 : {
522 15 : BlockTimer bt(Timer::POSTPRO_GRIDFUNCTION);
523 :
524 : // Create output directory if it doesn't exist.
525 15 : if (Mpi::Root(fem_op->GetComm()))
526 : {
527 10 : fs::create_directories(gridfunction_output_dir);
528 : }
529 :
530 : auto mesh_Lc0 = units.GetMeshLengthRelativeScale();
531 :
532 : // Given the electric field and magnetic flux density, write the fields to disk for
533 : // visualization. Write the mesh coordinates in the same units as originally input.
534 15 : mfem::ParMesh &mesh = E ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
535 15 : mesh::DimensionalizeMesh(mesh, mesh_Lc0);
536 15 : ScaleGridFunctions(mesh_Lc0, mesh.Dimension(), HasComplexGridFunction<solver_t>(), E, B,
537 15 : V, A);
538 :
539 : // Create grid function for vector coefficients.
540 15 : mfem::ParFiniteElementSpace &fespace = E ? *E->ParFESpace() : *B->ParFESpace();
541 15 : mfem::ParGridFunction gridfunc_vector(&fespace);
542 :
543 : // Create grid function for scalar coefficients.
544 15 : mfem::H1_FECollection h1_fec(fespace.GetMaxElementOrder(), mesh.Dimension());
545 15 : mfem::ParFiniteElementSpace h1_fespace(&mesh, &h1_fec);
546 15 : mfem::ParGridFunction gridfunc_scalar(&h1_fespace);
547 :
548 15 : const int local_rank = mesh.GetMyRank();
549 :
550 90 : auto write_grid_function = [&](const auto &gridfunc, const std::string &name)
551 : {
552 75 : auto path = fs::path(gridfunction_output_dir) /
553 150 : fmt::format("{}_{:0{}d}.gf.{:0{}d}", name, step, pad_digits_default,
554 : local_rank, pad_digits_default);
555 : std::ofstream file(path);
556 75 : gridfunc.Save(file);
557 75 : };
558 :
559 : // Write grid functions using MFEM's built-in Save method.
560 : // Use 6-digit padding to match MFEM's pad_digits_default.
561 : if constexpr (HasEGridFunction<solver_t>())
562 : {
563 12 : if (E)
564 : {
565 : if constexpr (HasComplexGridFunction<solver_t>())
566 : {
567 : // Write real and imaginary parts separately.
568 6 : write_grid_function(E->Real(), "E_real");
569 12 : write_grid_function(E->Imag(), "E_imag");
570 : }
571 : else
572 : {
573 : // Write real part only.
574 12 : write_grid_function(E->Real(), "E");
575 : }
576 : }
577 : }
578 :
579 : if constexpr (HasBGridFunction<solver_t>())
580 : {
581 12 : if (B)
582 : {
583 : if constexpr (HasComplexGridFunction<solver_t>())
584 : {
585 : // Write real and imaginary parts separately.
586 6 : write_grid_function(B->Real(), "B_real");
587 12 : write_grid_function(B->Imag(), "B_imag");
588 : }
589 : else
590 : {
591 : // Write real part only.
592 12 : write_grid_function(B->Real(), "B");
593 : }
594 : }
595 : }
596 :
597 : if constexpr (HasVGridFunction<solver_t>())
598 : {
599 3 : if (V)
600 : {
601 6 : write_grid_function(V->Real(), "V");
602 : }
603 : }
604 :
605 : if constexpr (HasAGridFunction<solver_t>())
606 : {
607 3 : if (A)
608 : {
609 6 : write_grid_function(A->Real(), "A");
610 : }
611 : }
612 :
613 15 : if (U_e)
614 : {
615 : gridfunc_scalar = 0.0;
616 12 : gridfunc_scalar.ProjectCoefficient(*U_e.get());
617 24 : write_grid_function(gridfunc_scalar, "U_e");
618 : }
619 :
620 15 : if (U_m)
621 : {
622 : gridfunc_scalar = 0.0;
623 12 : gridfunc_scalar.ProjectCoefficient(*U_m.get());
624 24 : write_grid_function(gridfunc_scalar, "U_m");
625 : }
626 :
627 15 : if (S)
628 : {
629 : gridfunc_vector = 0.0;
630 9 : gridfunc_vector.ProjectCoefficient(*S.get());
631 18 : write_grid_function(gridfunc_vector, "S");
632 : }
633 :
634 15 : mesh::NondimensionalizeMesh(mesh, mesh_Lc0);
635 15 : ScaleGridFunctions(1.0 / mesh_Lc0, mesh.Dimension(), HasComplexGridFunction<solver_t>(),
636 : E, B, V, A);
637 15 : Mpi::Barrier(fem_op->GetComm());
638 15 : }
639 :
640 : template <ProblemType solver_t>
641 0 : void PostOperator<solver_t>::WriteMFEMGridFunctionsFinal(const ErrorIndicator *indicator)
642 : {
643 0 : BlockTimer bt(Timer::POSTPRO_GRIDFUNCTION);
644 :
645 : auto mesh_Lc0 = units.GetMeshLengthRelativeScale();
646 :
647 : // Write the mesh partitioning and (optionally) error indicators at the final step.
648 : // Write the mesh coordinates in the same units as originally input.
649 0 : mfem::ParMesh &mesh = E ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
650 0 : mesh::DimensionalizeMesh(mesh, mesh_Lc0);
651 :
652 : // Create output directory if it doesn't exist.
653 0 : if (Mpi::Root(fem_op->GetComm()))
654 : {
655 0 : fs::create_directories(gridfunction_output_dir);
656 : }
657 :
658 : // Create piecewise constant finite element space for rank and error indicator.
659 0 : mfem::L2_FECollection pwconst_fec(0, mesh.Dimension());
660 0 : mfem::FiniteElementSpace pwconst_fespace(&mesh, &pwconst_fec);
661 :
662 0 : const int local_rank = mesh.GetMyRank();
663 :
664 0 : auto write_grid_function = [&](const auto &gridfunc, const std::string &name)
665 : {
666 0 : auto path = fs::path(gridfunction_output_dir) /
667 0 : fmt::format("{}.gf.{:0{}d}", name, local_rank, pad_digits_default);
668 : std::ofstream file(path);
669 0 : gridfunc.Save(file);
670 0 : };
671 :
672 : // Write mesh partitioning (rank information).
673 : {
674 0 : mfem::GridFunction rank(&pwconst_fespace);
675 0 : rank = local_rank + 1;
676 0 : write_grid_function(rank, "rank");
677 0 : }
678 :
679 : // Write error indicator if provided.
680 0 : if (indicator)
681 : {
682 0 : mfem::GridFunction eta(&pwconst_fespace);
683 0 : MFEM_VERIFY(eta.Size() == indicator->Local().Size(),
684 : "Size mismatch for provided ErrorIndicator for postprocessing!");
685 0 : eta = indicator->Local();
686 0 : write_grid_function(eta, "indicator");
687 0 : }
688 :
689 : // Save ParMesh files; necessary to visualize grid functions.
690 0 : fs::path mesh_filename = fs::path(gridfunction_output_dir) / "mesh";
691 0 : mesh.Save(mesh_filename);
692 :
693 0 : mesh::NondimensionalizeMesh(mesh, mesh_Lc0);
694 0 : Mpi::Barrier(fem_op->GetComm());
695 0 : }
696 :
697 : // Measurements.
698 :
699 : template <ProblemType solver_t>
700 15 : void PostOperator<solver_t>::MeasureDomainFieldEnergy() const
701 : {
702 : measurement_cache.domain_E_field_energy_i.clear();
703 : measurement_cache.domain_H_field_energy_i.clear();
704 :
705 15 : measurement_cache.domain_E_field_energy_i.reserve(dom_post_op.M_i.size());
706 15 : measurement_cache.domain_H_field_energy_i.reserve(dom_post_op.M_i.size());
707 :
708 : if constexpr (HasEGridFunction<solver_t>())
709 : {
710 : // Use V if it has it rather than E.
711 12 : auto &field = V ? *V : *E;
712 12 : auto energy = dom_post_op.GetElectricFieldEnergy(field);
713 12 : measurement_cache.domain_E_field_energy_all = energy;
714 :
715 12 : for (const auto &[idx, data] : dom_post_op.M_i)
716 : {
717 0 : auto energy_i = dom_post_op.GetDomainElectricFieldEnergy(idx, field);
718 0 : auto participation_ratio = std::abs(energy_i) > 0.0 ? energy_i / energy : 0.0;
719 0 : measurement_cache.domain_E_field_energy_i.emplace_back(
720 : Measurement::DomainData{idx, energy_i, participation_ratio});
721 : }
722 : }
723 : else
724 : {
725 : // Magnetic field only.
726 3 : measurement_cache.domain_E_field_energy_all = 0.0;
727 3 : for (const auto &[idx, data] : dom_post_op.M_i)
728 : {
729 0 : measurement_cache.domain_E_field_energy_i.emplace_back(
730 : Measurement::DomainData{idx, 0.0, 0.0});
731 : }
732 : }
733 :
734 : if (HasBGridFunction<solver_t>())
735 : {
736 12 : auto &field = A ? *A : *B;
737 12 : auto energy = dom_post_op.GetMagneticFieldEnergy(field);
738 12 : measurement_cache.domain_H_field_energy_all = energy;
739 :
740 12 : for (const auto &[idx, data] : dom_post_op.M_i)
741 : {
742 0 : auto energy_i = dom_post_op.GetDomainMagneticFieldEnergy(idx, field);
743 0 : auto participation_ratio = std::abs(energy) > 0.0 ? energy_i / energy : 0.0;
744 0 : measurement_cache.domain_H_field_energy_i.emplace_back(
745 : Measurement::DomainData{idx, energy_i, participation_ratio});
746 : }
747 : }
748 : else
749 : {
750 : // Electric field only.
751 3 : measurement_cache.domain_H_field_energy_all = 0.0;
752 3 : for (const auto &[idx, data] : dom_post_op.M_i)
753 : {
754 0 : measurement_cache.domain_H_field_energy_i.emplace_back(
755 : Measurement::DomainData{idx, 0.0, 0.0});
756 : }
757 : }
758 :
759 : // Log Domain Energy.
760 9 : const auto domain_E = units.Dimensionalize<Units::ValueType::ENERGY>(
761 : measurement_cache.domain_E_field_energy_all);
762 9 : const auto domain_H = units.Dimensionalize<Units::ValueType::ENERGY>(
763 : measurement_cache.domain_H_field_energy_all);
764 : if constexpr (HasEGridFunction<solver_t>() && !HasBGridFunction<solver_t>())
765 : {
766 3 : Mpi::Print(" Field energy E = {:.3e} J\n", domain_E);
767 : }
768 : else if constexpr (!HasEGridFunction<solver_t>() && HasBGridFunction<solver_t>())
769 : {
770 3 : Mpi::Print(" Field energy H = {:.3e} J\n", domain_H);
771 : }
772 : else if constexpr (solver_t != ProblemType::EIGENMODE)
773 : {
774 6 : Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", domain_E, domain_H,
775 6 : domain_E + domain_H);
776 : }
777 15 : }
778 :
779 : template <ProblemType solver_t>
780 9 : void PostOperator<solver_t>::MeasureLumpedPorts() const
781 : {
782 : measurement_cache.lumped_port_vi.clear();
783 15 : measurement_cache.lumped_port_inductor_energy = 0.0;
784 15 : measurement_cache.lumped_port_capacitor_energy = 0.0;
785 :
786 : if constexpr (solver_t == ProblemType::EIGENMODE || solver_t == ProblemType::DRIVEN ||
787 : solver_t == ProblemType::TRANSIENT)
788 : {
789 15 : for (const auto &[idx, data] : fem_op->GetLumpedPortOp())
790 : {
791 6 : auto &vi = measurement_cache.lumped_port_vi[idx];
792 6 : vi.P = data.GetPower(*E, *B);
793 6 : vi.V = data.GetVoltage(*E);
794 : if constexpr (solver_t == ProblemType::EIGENMODE || solver_t == ProblemType::DRIVEN)
795 : {
796 : // Compute current from the port impedance, separate contributions for R, L, C
797 : // branches.
798 : // Get value and make real: Matches current behaviour (even for eigensolver!).
799 3 : MFEM_VERIFY(
800 : measurement_cache.freq.real() > 0.0,
801 : "Frequency domain lumped port postprocessing requires nonzero frequency!");
802 3 : vi.I_RLC[0] =
803 3 : (std::abs(data.R) > 0.0)
804 3 : ? vi.V / data.GetCharacteristicImpedance(measurement_cache.freq.real(),
805 : LumpedPortData::Branch::R)
806 : : 0.0;
807 3 : vi.I_RLC[1] =
808 3 : (std::abs(data.L) > 0.0)
809 3 : ? vi.V / data.GetCharacteristicImpedance(measurement_cache.freq.real(),
810 : LumpedPortData::Branch::L)
811 : : 0.0;
812 3 : vi.I_RLC[2] =
813 3 : (std::abs(data.C) > 0.0)
814 6 : ? vi.V / data.GetCharacteristicImpedance(measurement_cache.freq.real(),
815 : LumpedPortData::Branch::C)
816 : : 0.0;
817 3 : vi.I = std::accumulate(vi.I_RLC.begin(), vi.I_RLC.end(),
818 : std::complex<double>{0.0, 0.0});
819 3 : vi.S = data.GetSParameter(*E);
820 :
821 : // Add contribution due to all inductive lumped boundaries in the model:
822 : // E_ind = ∑_j 1/2 L_j I_mj².
823 3 : if (std::abs(data.L) > 0.0)
824 : {
825 0 : std::complex<double> I_mj = vi.I_RLC[1];
826 0 : vi.inductor_energy = 0.5 * std::abs(data.L) * std::real(I_mj * std::conj(I_mj));
827 0 : measurement_cache.lumped_port_inductor_energy += vi.inductor_energy;
828 : }
829 :
830 : // Add contribution due to all capacitive lumped boundaries in the model:
831 : // E_cap = ∑_j 1/2 C_j V_mj².
832 3 : if (std::abs(data.C) > 0.0)
833 : {
834 0 : std::complex<double> V_mj = vi.V;
835 0 : vi.capacitor_energy = 0.5 * std::abs(data.C) * std::real(V_mj * std::conj(V_mj));
836 0 : measurement_cache.lumped_port_capacitor_energy += vi.capacitor_energy;
837 : }
838 : }
839 : else
840 : {
841 : // Compute current from P = V I^* since there is no frequency & characteristic
842 : // impedance of the lumped element.
843 6 : vi.I = (std::abs(vi.V) > 0.0) ? std::conj(vi.P / vi.V) : 0.0;
844 : }
845 : }
846 : }
847 9 : }
848 :
849 : template <ProblemType solver_t>
850 3 : void PostOperator<solver_t>::MeasureLumpedPortsEig() const
851 : {
852 : // Depends on MeasureLumpedPorts.
853 : if constexpr (solver_t == ProblemType::EIGENMODE)
854 : {
855 : auto freq_re = measurement_cache.freq.real();
856 3 : auto energy_electric_all = measurement_cache.domain_E_field_energy_all +
857 3 : measurement_cache.lumped_port_capacitor_energy;
858 3 : for (const auto &[idx, data] : fem_op->GetLumpedPortOp())
859 : {
860 : // Get previously computed data: should never fail as defined by MeasureLumpedPorts.
861 0 : auto &vi = measurement_cache.lumped_port_vi.at(idx);
862 :
863 : // Resistive Lumped Ports:
864 : // Compute participation ratio of external ports (given as any port boundary with
865 : // nonzero resistance). Currently no reactance of the ports is supported. The κ of
866 : // the port follows from:
867 : // κ_mj = 1/2 R_j I_mj² / E_m
868 : // from which the mode coupling quality factor is computed as:
869 : // Q_mj = ω_m / κ_mj.
870 0 : if (std::abs(data.R) > 0.0)
871 : {
872 0 : std::complex<double> I_mj = vi.I_RLC[0];
873 : // Power = 1/2 R_j I_mj².
874 : // Note conventions: mean(I²) = (I_r² + I_i²) / 2;
875 0 : auto resistor_power = 0.5 * std::abs(data.R) * std::real(I_mj * std::conj(I_mj));
876 0 : vi.mode_port_kappa =
877 0 : std::copysign(resistor_power / energy_electric_all, I_mj.real());
878 0 : vi.quality_factor = (vi.mode_port_kappa == 0.0)
879 0 : ? mfem::infinity()
880 : : freq_re / std::abs(vi.mode_port_kappa);
881 : }
882 :
883 : // Inductive Lumped Ports:
884 : // Compute energy-participation ratio of junction given by index idx for the field
885 : // mode. We first get the port line voltage, and use lumped port circuit impedance to
886 : // get peak current through the inductor: I_mj = V_mj / Z_mj, Z_mj = i ω_m L_j. E_m
887 : // is the total energy in mode m: E_m = E_elec + E_cap = E_mag + E_ind. The signed EPR
888 : // for a lumped inductive element is computed as:
889 : // p_mj = 1/2 L_j I_mj² / E_m.
890 : // An element with no assigned inductance will be treated as having zero admittance
891 : // and thus zero current.
892 0 : if (std::abs(data.L) > 0.0)
893 : {
894 0 : std::complex<double> I_mj = vi.I_RLC[1];
895 0 : vi.inductive_energy_participation =
896 0 : std::copysign(vi.inductor_energy / energy_electric_all, I_mj.real());
897 : }
898 : }
899 : }
900 3 : }
901 :
902 : template <ProblemType solver_t>
903 3 : void PostOperator<solver_t>::MeasureWavePorts() const
904 : {
905 : measurement_cache.wave_port_vi.clear();
906 :
907 : if constexpr (solver_t == ProblemType::DRIVEN)
908 : {
909 3 : for (const auto &[idx, data] : fem_op->GetWavePortOp())
910 : {
911 : // Get value and make real: Matches current behaviour.
912 : auto freq_re = measurement_cache.freq.real(); // TODO: Fix
913 0 : MFEM_VERIFY(freq_re > 0.0,
914 : "Frequency domain wave port postprocessing requires nonzero frequency!");
915 0 : auto &vi = measurement_cache.wave_port_vi[idx];
916 0 : vi.P = data.GetPower(*E, *B);
917 0 : vi.S = data.GetSParameter(*E);
918 : // vi.V = vi.I[0] = vi.I[1] = vi.I[2] = 0.0; // Not yet implemented
919 : // // (Z = V² / P, I = V / Z)
920 : }
921 : }
922 3 : }
923 :
924 : template <ProblemType solver_t>
925 3 : void PostOperator<solver_t>::MeasureSParameter() const
926 : {
927 : // Depends on LumpedPorts, WavePorts.
928 : if constexpr (solver_t == ProblemType::DRIVEN)
929 : {
930 : using fmt::format;
931 : using std::complex_literals::operator""i;
932 :
933 : // Don't measure S-Matrix unless there is only one excitation per port. Also, we current
934 : // don't support mixing wave and lumped ports, because we need to fix consistent
935 : // conventions / de-embedding.
936 3 : if (!fem_op->GetPortExcitations().IsMultipleSimple() ||
937 3 : !((fem_op->GetLumpedPortOp().Size() > 0) xor (fem_op->GetWavePortOp().Size() > 0)))
938 : {
939 : return;
940 : }
941 :
942 : // Assumes that for single driving port the excitation index is equal to the port index.
943 3 : auto drive_port_idx = measurement_cache.ex_idx;
944 :
945 : // Currently S-Parameters are not calculated for mixed lumped & wave ports, so don't
946 : // combine output iterators.
947 6 : for (const auto &[idx, data] : fem_op->GetLumpedPortOp())
948 : {
949 : // Get previously computed data: should never fail as defined by MeasureLumpedPorts.
950 3 : auto &vi = measurement_cache.lumped_port_vi.at(idx);
951 :
952 3 : const LumpedPortData &src_data = fem_op->GetLumpedPortOp().GetPort(drive_port_idx);
953 3 : if (idx == drive_port_idx)
954 : {
955 3 : vi.S.real(vi.S.real() - 1.0);
956 : }
957 : // Generalized S-parameters if the ports are resistive (avoids divide-by-zero).
958 3 : if (std::abs(data.R) > 0.0)
959 : {
960 3 : vi.S *= std::sqrt(src_data.R / data.R);
961 : }
962 :
963 3 : Mpi::Print(" {0} = {1:+.3e}{2:+.3e}i, |{0}| = {3:+.3e}, arg({0}) = {4:+.3e}\n",
964 3 : format("S[{}][{}]", idx, drive_port_idx), vi.S.real(), vi.S.imag(),
965 6 : Measurement::Magnitude(vi.S), Measurement::Phase(vi.S));
966 : }
967 3 : for (const auto &[idx, data] : fem_op->GetWavePortOp())
968 : {
969 : // Get previously computed data: should never fail as defined by MeasureWavePorts.
970 0 : auto &vi = measurement_cache.wave_port_vi.at(idx);
971 :
972 : // Wave port modes are not normalized to a characteristic impedance so no generalized
973 : // S-parameters are available.
974 0 : const WavePortData &src_data = fem_op->GetWavePortOp().GetPort(drive_port_idx);
975 0 : if (idx == drive_port_idx)
976 : {
977 0 : vi.S.real(vi.S.real() - 1.0);
978 : }
979 : // Port de-embedding: S_demb = S exp(ikₙᵢ dᵢ) exp(ikₙⱼ dⱼ) (distance offset is default
980 : // 0 unless specified).
981 : vi.S *= std::exp(1i * src_data.kn0 * src_data.d_offset);
982 : vi.S *= std::exp(1i * data.kn0 * data.d_offset);
983 :
984 0 : Mpi::Print(" {0} = {1:+.3e}{2:+.3e}i, |{0}| = {3:+.3e}, arg({0}) = {4:+.3e}\n",
985 0 : format("S[{}][{}]", idx, drive_port_idx), vi.S.real(), vi.S.imag(),
986 0 : Measurement::Magnitude(vi.S), Measurement::Phase(vi.S));
987 : }
988 : }
989 0 : }
990 :
991 : template <ProblemType solver_t>
992 15 : void PostOperator<solver_t>::MeasureSurfaceFlux() const
993 : {
994 : // Compute the flux through a surface as Φ_j = ∫ F ⋅ n_j dS, with F = B, F = ε D, or F =
995 : // E x H. The special coefficient is used to avoid issues evaluating MFEM GridFunctions
996 : // which are discontinuous at interior boundary elements.
997 : measurement_cache.surface_flux_i.clear();
998 15 : measurement_cache.surface_flux_i.reserve(surf_post_op.flux_surfs.size());
999 15 : for (const auto &[idx, data] : surf_post_op.flux_surfs)
1000 : {
1001 0 : measurement_cache.surface_flux_i.emplace_back(Measurement::FluxData{
1002 0 : idx, surf_post_op.GetSurfaceFlux(idx, E.get(), B.get()), data.type});
1003 : }
1004 15 : }
1005 :
1006 : template <ProblemType solver_t>
1007 6 : void PostOperator<solver_t>::MeasureFarField() const
1008 : {
1009 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::EIGENMODE)
1010 : {
1011 6 : measurement_cache.farfield.thetaphis = surf_post_op.farfield.thetaphis;
1012 :
1013 : // NOTE: measurement_cache.freq is omega (it has a factor of 2pi).
1014 6 : measurement_cache.farfield.E_field = surf_post_op.GetFarFieldrE(
1015 : measurement_cache.farfield.thetaphis, *E, *B, measurement_cache.freq.real(),
1016 6 : measurement_cache.freq.imag());
1017 : }
1018 6 : }
1019 :
1020 : template <ProblemType solver_t>
1021 12 : void PostOperator<solver_t>::MeasureInterfaceEFieldEnergy() const
1022 : {
1023 : // Depends on Lumped Port Energy since this is used in normalization of participation
1024 : // ratio.
1025 :
1026 : // Compute the surface dielectric participation ratio and associated quality factor for
1027 : // the material interface given by index idx. We have:
1028 : // 1/Q_mj = p_mj tan(δ)_j
1029 : // with:
1030 : // p_mj = 1/2 t_j Re{∫_{Γ_j} (ε_j E_m)ᴴ E_m dS} / (E_elec + E_cap).
1031 : measurement_cache.interface_eps_i.clear();
1032 : if constexpr (HasEGridFunction<solver_t>())
1033 : {
1034 : // Domain and port energies must have been measured first. E_cap returns zero if the
1035 : // solver does not support lumped ports.
1036 : //
1037 : // TODO: Should this not include other types of energy too (surface impedance case)?
1038 12 : auto energy_electric_all = measurement_cache.domain_E_field_energy_all +
1039 12 : measurement_cache.lumped_port_capacitor_energy;
1040 :
1041 12 : measurement_cache.interface_eps_i.reserve(surf_post_op.eps_surfs.size());
1042 12 : for (const auto &[idx, data] : surf_post_op.eps_surfs)
1043 : {
1044 0 : auto energy = surf_post_op.GetInterfaceElectricFieldEnergy(idx, *E);
1045 :
1046 0 : auto energy_participation_p = energy / energy_electric_all;
1047 0 : auto loss_tangent_delta = surf_post_op.GetInterfaceLossTangent(idx);
1048 0 : auto quality_factor_Q = (energy_participation_p == 0.0 || loss_tangent_delta == 0.0)
1049 0 : ? mfem::infinity()
1050 0 : : 1.0 / (loss_tangent_delta * energy_participation_p);
1051 :
1052 0 : measurement_cache.interface_eps_i.emplace_back(Measurement::InterfaceData{
1053 : idx, energy, loss_tangent_delta, energy_participation_p, quality_factor_Q});
1054 : }
1055 : }
1056 12 : }
1057 :
1058 : template <ProblemType solver_t>
1059 15 : void PostOperator<solver_t>::MeasureProbes() const
1060 : {
1061 : measurement_cache.probe_E_field.clear();
1062 : measurement_cache.probe_B_field.clear();
1063 :
1064 : #if defined(MFEM_USE_GSLIB)
1065 : if constexpr (HasEGridFunction<solver_t>())
1066 : {
1067 12 : if (interp_op.GetProbes().size() > 0)
1068 : {
1069 0 : measurement_cache.probe_E_field = interp_op.ProbeField(*E);
1070 : }
1071 : }
1072 : if constexpr (HasBGridFunction<solver_t>())
1073 : {
1074 12 : if (interp_op.GetProbes().size() > 0)
1075 : {
1076 0 : measurement_cache.probe_B_field = interp_op.ProbeField(*B);
1077 : }
1078 : }
1079 : #endif
1080 15 : }
1081 :
1082 : using fmt::format;
1083 :
1084 : template <ProblemType solver_t>
1085 : template <ProblemType U>
1086 3 : auto PostOperator<solver_t>::MeasureAndPrintAll(int ex_idx, int step,
1087 : const ComplexVector &e,
1088 : const ComplexVector &b,
1089 : std::complex<double> omega)
1090 : -> std::enable_if_t<U == ProblemType::DRIVEN, double>
1091 : {
1092 3 : BlockTimer bt0(Timer::POSTPRO);
1093 3 : SetEGridFunction(e);
1094 3 : SetBGridFunction(b);
1095 :
1096 3 : measurement_cache = {};
1097 3 : measurement_cache.freq = omega;
1098 3 : measurement_cache.ex_idx = ex_idx;
1099 3 : MeasureAllImpl();
1100 :
1101 : omega = units.Dimensionalize<Units::ValueType::FREQUENCY>(omega);
1102 3 : post_op_csv.PrintAllCSVData(*this, measurement_cache, omega.real(), step, ex_idx);
1103 6 : if (ShouldWriteParaviewFields(step))
1104 : {
1105 3 : Mpi::Print("\n");
1106 3 : auto ind = 1 + std::distance(output_save_indices.begin(),
1107 : std::lower_bound(output_save_indices.begin(),
1108 : output_save_indices.end(), step));
1109 3 : WriteParaviewFields(omega.real(), ind);
1110 6 : Mpi::Print(" Wrote fields to disk (Paraview) at step {:d}\n", step + 1);
1111 : }
1112 3 : if (ShouldWriteGridFunctionFields(step))
1113 : {
1114 3 : Mpi::Print("\n");
1115 3 : auto ind = 1 + std::distance(output_save_indices.begin(),
1116 : std::lower_bound(output_save_indices.begin(),
1117 : output_save_indices.end(), step));
1118 3 : WriteMFEMGridFunctions(omega.real(), ind);
1119 6 : Mpi::Print(" Wrote fields to disk (grid function) at step {:d}\n", step + 1);
1120 : }
1121 3 : return measurement_cache.domain_E_field_energy_all +
1122 3 : measurement_cache.domain_H_field_energy_all;
1123 6 : }
1124 :
1125 : template <ProblemType solver_t>
1126 : template <ProblemType U>
1127 3 : auto PostOperator<solver_t>::MeasureAndPrintAll(int step, const ComplexVector &e,
1128 : const ComplexVector &b,
1129 : std::complex<double> omega,
1130 : double error_abs, double error_bkwd,
1131 : int num_conv)
1132 : -> std::enable_if_t<U == ProblemType::EIGENMODE, double>
1133 : {
1134 3 : BlockTimer bt0(Timer::POSTPRO);
1135 3 : SetEGridFunction(e);
1136 3 : SetBGridFunction(b);
1137 :
1138 3 : measurement_cache = {};
1139 3 : measurement_cache.freq = omega;
1140 3 : measurement_cache.eigenmode_Q =
1141 3 : (omega == 0.0) ? mfem::infinity() : 0.5 * std::abs(omega) / std::abs(omega.imag());
1142 3 : measurement_cache.error_abs = error_abs;
1143 3 : measurement_cache.error_bkwd = error_bkwd;
1144 :
1145 : // Mini pretty-print table of eig summaries: always print with header since other
1146 : // measurements may log their results.
1147 3 : if (Mpi::Root(fem_op->GetComm()))
1148 : {
1149 2 : Table table;
1150 2 : int idx_pad = 1 + static_cast<int>(std::log10(num_conv));
1151 4 : table.col_options = {6, 6};
1152 6 : table.insert(Column("idx", "m", idx_pad, {}, {}, "") << step + 1);
1153 6 : table.insert(Column("f_re", "Re{f} (GHz)")
1154 : << units.Dimensionalize<Units::ValueType::FREQUENCY>(omega.real()));
1155 6 : table.insert(Column("f_im", "Im{f} (GHz)")
1156 : << units.Dimensionalize<Units::ValueType::FREQUENCY>(omega.imag()));
1157 4 : table.insert(Column("q", "Q") << measurement_cache.eigenmode_Q);
1158 6 : table.insert(Column("err_back", "Error (Bkwd.)") << error_bkwd);
1159 6 : table.insert(Column("err_abs", "Error (Abs.)") << error_abs);
1160 2 : table[0].print_as_int = true;
1161 4 : Mpi::Print("{}", (step == 0) ? table.format_table() : table.format_row(0));
1162 2 : }
1163 3 : MeasureAllImpl();
1164 :
1165 3 : int print_idx = step + 1;
1166 3 : post_op_csv.PrintAllCSVData(*this, measurement_cache, print_idx, step);
1167 6 : if (ShouldWriteParaviewFields(step))
1168 : {
1169 3 : WriteParaviewFields(step, print_idx);
1170 3 : Mpi::Print(" Wrote mode {:d} to disk (Paraview)\n", print_idx);
1171 : }
1172 3 : if (ShouldWriteGridFunctionFields(step))
1173 : {
1174 3 : WriteMFEMGridFunctions(step, print_idx);
1175 3 : Mpi::Print(" Wrote mode {:d} to disk (grid function)\n", print_idx);
1176 : }
1177 3 : return measurement_cache.domain_E_field_energy_all +
1178 3 : measurement_cache.domain_H_field_energy_all;
1179 8 : }
1180 :
1181 : template <ProblemType solver_t>
1182 : template <ProblemType U>
1183 3 : auto PostOperator<solver_t>::MeasureAndPrintAll(int step, const Vector &v, const Vector &e,
1184 : int idx)
1185 : -> std::enable_if_t<U == ProblemType::ELECTROSTATIC, double>
1186 : {
1187 3 : BlockTimer bt0(Timer::POSTPRO);
1188 : SetVGridFunction(v);
1189 : SetEGridFunction(e);
1190 :
1191 3 : measurement_cache = {};
1192 3 : MeasureAllImpl();
1193 :
1194 3 : int print_idx = step + 1;
1195 3 : post_op_csv.PrintAllCSVData(*this, measurement_cache, print_idx, step);
1196 6 : if (ShouldWriteParaviewFields(step))
1197 : {
1198 3 : Mpi::Print("\n");
1199 3 : WriteParaviewFields(step, idx);
1200 3 : Mpi::Print(" Wrote fields to disk (Paraview) for source {:d}\n", idx);
1201 : }
1202 3 : if (ShouldWriteGridFunctionFields(step))
1203 : {
1204 3 : Mpi::Print("\n");
1205 3 : WriteMFEMGridFunctions(step, idx);
1206 3 : Mpi::Print(" Wrote fields to disk (grid function) for source {:d}\n", idx);
1207 : }
1208 3 : return measurement_cache.domain_E_field_energy_all +
1209 3 : measurement_cache.domain_H_field_energy_all;
1210 6 : }
1211 : template <ProblemType solver_t>
1212 : template <ProblemType U>
1213 3 : auto PostOperator<solver_t>::MeasureAndPrintAll(int step, const Vector &a, const Vector &b,
1214 : int idx)
1215 : -> std::enable_if_t<U == ProblemType::MAGNETOSTATIC, double>
1216 : {
1217 3 : BlockTimer bt0(Timer::POSTPRO);
1218 : SetAGridFunction(a);
1219 : SetBGridFunction(b);
1220 :
1221 3 : measurement_cache = {};
1222 3 : MeasureAllImpl();
1223 :
1224 3 : int print_idx = step + 1;
1225 3 : post_op_csv.PrintAllCSVData(*this, measurement_cache, print_idx, step);
1226 6 : if (ShouldWriteParaviewFields(step))
1227 : {
1228 3 : Mpi::Print("\n");
1229 3 : WriteParaviewFields(step, idx);
1230 3 : Mpi::Print(" Wrote fields to disk (Paraview) for source {:d}\n", idx);
1231 : }
1232 3 : if (ShouldWriteGridFunctionFields(step))
1233 : {
1234 3 : Mpi::Print("\n");
1235 3 : WriteMFEMGridFunctions(step, idx);
1236 3 : Mpi::Print(" Wrote fields to disk (grid function) for source {:d}\n", idx);
1237 : }
1238 3 : return measurement_cache.domain_E_field_energy_all +
1239 3 : measurement_cache.domain_H_field_energy_all;
1240 6 : }
1241 :
1242 : template <ProblemType solver_t>
1243 : template <ProblemType U>
1244 3 : auto PostOperator<solver_t>::MeasureAndPrintAll(int step, const Vector &e, const Vector &b,
1245 : double time, double J_coef)
1246 : -> std::enable_if_t<U == ProblemType::TRANSIENT, double>
1247 : {
1248 3 : BlockTimer bt0(Timer::POSTPRO);
1249 : SetEGridFunction(e);
1250 : SetBGridFunction(b);
1251 :
1252 3 : measurement_cache = {};
1253 3 : measurement_cache.Jcoeff_excitation = J_coef;
1254 3 : MeasureAllImpl();
1255 :
1256 : // Time must be converted before passing into csv due to the shared PrintAllCSVData
1257 : // method.
1258 : time = units.Dimensionalize<Units::ValueType::TIME>(time);
1259 3 : post_op_csv.PrintAllCSVData(*this, measurement_cache, time, step);
1260 6 : if (ShouldWriteParaviewFields(step))
1261 : {
1262 3 : Mpi::Print("\n");
1263 3 : WriteParaviewFields(double(step) / output_delta_post, time);
1264 6 : Mpi::Print(" Wrote fields to disk (Paraview) at step {:d}\n", step + 1);
1265 : }
1266 3 : if (ShouldWriteGridFunctionFields(step))
1267 : {
1268 3 : Mpi::Print("\n");
1269 3 : WriteMFEMGridFunctions(double(step) / output_delta_post, time);
1270 6 : Mpi::Print(" Wrote fields to disk (grid function) at step {:d}\n", step + 1);
1271 : }
1272 3 : return measurement_cache.domain_E_field_energy_all +
1273 3 : measurement_cache.domain_H_field_energy_all;
1274 6 : }
1275 :
1276 : template <ProblemType solver_t>
1277 0 : void PostOperator<solver_t>::MeasureFinalize(const ErrorIndicator &indicator)
1278 : {
1279 0 : BlockTimer bt0(Timer::POSTPRO);
1280 0 : auto indicator_stats = indicator.GetSummaryStatistics(fem_op->GetComm());
1281 0 : post_op_csv.PrintErrorIndicator(Mpi::Root(fem_op->GetComm()), indicator_stats);
1282 : if (ShouldWriteParaviewFields())
1283 : {
1284 0 : WriteParaviewFieldsFinal(&indicator);
1285 : }
1286 : if (ShouldWriteGridFunctionFields())
1287 : {
1288 0 : WriteMFEMGridFunctionsFinal(&indicator);
1289 : }
1290 0 : }
1291 :
1292 : template <ProblemType solver_t>
1293 : template <ProblemType U>
1294 0 : auto PostOperator<solver_t>::MeasureDomainFieldEnergyOnly(const ComplexVector &e,
1295 : const ComplexVector &b)
1296 : -> std::enable_if_t<U == ProblemType::DRIVEN, double>
1297 : {
1298 0 : SetEGridFunction(e);
1299 0 : SetBGridFunction(b);
1300 0 : MeasureDomainFieldEnergy();
1301 0 : Mpi::Barrier(fem_op->GetComm());
1302 :
1303 : // Return total domain energy for normalizing error indicator.
1304 0 : return measurement_cache.domain_E_field_energy_all +
1305 0 : measurement_cache.domain_H_field_energy_all;
1306 : }
1307 :
1308 : // Explict template instantiation.
1309 :
1310 : template class PostOperator<ProblemType::DRIVEN>;
1311 : template class PostOperator<ProblemType::EIGENMODE>;
1312 : template class PostOperator<ProblemType::ELECTROSTATIC>;
1313 : template class PostOperator<ProblemType::MAGNETOSTATIC>;
1314 : template class PostOperator<ProblemType::TRANSIENT>;
1315 :
1316 : // Function explict instantiation.
1317 : // TODO(C++20): with requires, we won't need a second template.
1318 :
1319 : template auto PostOperator<ProblemType::DRIVEN>::MeasureAndPrintAll<ProblemType::DRIVEN>(
1320 : int ex_idx, int step, const ComplexVector &e, const ComplexVector &b,
1321 : std::complex<double> omega) -> double;
1322 :
1323 : template auto
1324 : PostOperator<ProblemType::EIGENMODE>::MeasureAndPrintAll<ProblemType::EIGENMODE>(
1325 : int step, const ComplexVector &e, const ComplexVector &b, std::complex<double> omega,
1326 : double error_abs, double error_bkwd, int num_conv) -> double;
1327 :
1328 : template auto
1329 : PostOperator<ProblemType::ELECTROSTATIC>::MeasureAndPrintAll<ProblemType::ELECTROSTATIC>(
1330 : int step, const Vector &v, const Vector &e, int idx) -> double;
1331 :
1332 : template auto
1333 : PostOperator<ProblemType::MAGNETOSTATIC>::MeasureAndPrintAll<ProblemType::MAGNETOSTATIC>(
1334 : int step, const Vector &a, const Vector &b, int idx) -> double;
1335 :
1336 : template auto
1337 : PostOperator<ProblemType::TRANSIENT>::MeasureAndPrintAll<ProblemType::TRANSIENT>(
1338 : int step, const Vector &e, const Vector &b, double t, double J_coef) -> double;
1339 :
1340 : template auto
1341 : PostOperator<ProblemType::DRIVEN>::MeasureDomainFieldEnergyOnly<ProblemType::DRIVEN>(
1342 : const ComplexVector &e, const ComplexVector &b) -> double;
1343 :
1344 : template auto
1345 : PostOperator<ProblemType::DRIVEN>::InitializeParaviewDataCollection(int ex_idx) -> void;
1346 :
1347 : } // namespace palace
|
