Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_MODELS_POST_OPERATOR_HPP
5 : #define PALACE_MODELS_POST_OPERATOR_HPP
6 :
7 : #include <complex>
8 : #include <map>
9 : #include <memory>
10 : #include <optional>
11 : #include <type_traits>
12 : #include <vector>
13 : #include <mfem.hpp>
14 : #include "fem/gridfunction.hpp"
15 : #include "fem/interpolator.hpp"
16 : #include "linalg/operator.hpp"
17 : #include "linalg/vector.hpp"
18 : #include "models/domainpostoperator.hpp"
19 : #include "models/lumpedportoperator.hpp"
20 : #include "models/postoperatorcsv.hpp"
21 : #include "models/surfacepostoperator.hpp"
22 : #include "utils/configfile.hpp"
23 : #include "utils/filesystem.hpp"
24 : #include "utils/units.hpp"
25 :
26 : namespace palace
27 : {
28 :
29 : class CurlCurlOperator;
30 : class ErrorIndicator;
31 : class IoData;
32 : class LaplaceOperator;
33 : class MaterialOperator;
34 : class SpaceOperator;
35 : class SurfaceCurrentOperator;
36 : class WavePortOperator;
37 :
38 : // Statically specify if solver uses real or complex fields.
39 :
40 : template <ProblemType solver_t>
41 : constexpr bool HasComplexGridFunction()
42 : {
43 : return solver_t == ProblemType::DRIVEN || solver_t == ProblemType::EIGENMODE;
44 : }
45 :
46 : // Statically specify what fields a solver uses
47 : // TODO(C++20): Change these to inline consteval and use with requires.
48 :
49 : template <ProblemType solver_t>
50 : constexpr bool HasVGridFunction()
51 : {
52 : return solver_t == ProblemType::ELECTROSTATIC;
53 : }
54 :
55 : template <ProblemType solver_t>
56 : constexpr bool HasAGridFunction()
57 : {
58 : return solver_t == ProblemType::MAGNETOSTATIC;
59 : }
60 :
61 : template <ProblemType solver_t>
62 : constexpr bool HasEGridFunction()
63 : {
64 : return solver_t != ProblemType::MAGNETOSTATIC;
65 : }
66 :
67 : template <ProblemType solver_t>
68 : constexpr bool HasBGridFunction()
69 : {
70 : return solver_t != ProblemType::ELECTROSTATIC;
71 : }
72 :
73 : //
74 : // A class to handle solution postprocessing for all solvers.
75 : //
76 : template <ProblemType solver_t>
77 : class PostOperator
78 : {
79 : protected:
80 : // Pointer to operator handling discretization and FEM space appropriate to solver. It
81 : // also contains the reference to all domains, boundary conditions, etc. needed for
82 : // measurement and printing.
83 : // TODO(C++20): Use std::reference_wrapper with incomplete types.
84 : fem_op_t<solver_t> *fem_op;
85 :
86 : // Unit converter from IOData to scale mesh and measurements. Lightweight class so it is
87 : // cheap to copy, rather than keep another reference to IOData.
88 : Units units;
89 :
90 : // Base post-op output directory.
91 : fs::path post_dir;
92 :
93 : // Fields: Electric, Magnetic, Scalar Potential, Vector Potential.
94 : std::unique_ptr<GridFunction> E, B, V, A;
95 :
96 : // Field output format control flags.
97 : bool enable_paraview_output = false;
98 : bool enable_gridfunction_output = false;
99 :
100 : // How many / which fields to output.
101 : int output_delta_post = 0; // printing rate (TRANSIENT)
102 : int output_n_post = 0; // max printing (OTHER SOLVERS)
103 : std::vector<std::size_t> output_save_indices = {}; // explicit saves
104 :
105 : // Whether any output formats were specified.
106 0 : bool AnyOutputFormats() const
107 : {
108 15 : return enable_paraview_output || enable_gridfunction_output;
109 : }
110 0 : bool AnythingToSave() const
111 : {
112 30 : return (output_delta_post > 0) || (output_n_post > 0) || !output_save_indices.empty();
113 : }
114 :
115 : // Whether any fields should be written at all.
116 0 : bool ShouldWriteFields() const { return AnyOutputFormats() && AnythingToSave(); }
117 :
118 : // Whether any fields should be written for this step.
119 30 : bool ShouldWriteFields(std::size_t step) const
120 : {
121 0 : return AnyOutputFormats() &&
122 30 : ((output_delta_post > 0 && step % output_delta_post == 0) ||
123 30 : (output_n_post > 0 && step < output_n_post) ||
124 6 : std::binary_search(output_save_indices.cbegin(), output_save_indices.cend(),
125 30 : step));
126 : }
127 :
128 : // Whether fields should be written for a particular output format (at a given step).
129 0 : bool ShouldWriteParaviewFields() const
130 : {
131 15 : return enable_paraview_output && AnythingToSave();
132 : }
133 0 : bool ShouldWriteParaviewFields(std::size_t step) const
134 : {
135 15 : return enable_paraview_output && ShouldWriteFields(step);
136 : }
137 0 : bool ShouldWriteGridFunctionFields() const
138 : {
139 0 : return enable_gridfunction_output && AnythingToSave();
140 : }
141 0 : bool ShouldWriteGridFunctionFields(std::size_t step) const
142 : {
143 15 : return enable_gridfunction_output && ShouldWriteFields(step);
144 : }
145 :
146 : // ParaView data collection: writing fields to disk for visualization.
147 : // This is an optional, since ParaViewDataCollection has no default (empty) ctor,
148 : // and we only want initialize it if ShouldWriteParaviewFields() returns true.
149 : std::optional<mfem::ParaViewDataCollection> paraview, paraview_bdr;
150 :
151 : // MFEM grid function output details.
152 : std::string gridfunction_output_dir;
153 : const std::size_t pad_digits_default = 6;
154 :
155 : // Measurements of field solution for ParaView files (full domain or surfaces).
156 :
157 : // Poynting Coefficient, Electric Boundary Field (re+im), Magnetic Boundary Field (re+im),
158 : // Vector Potential Boundary Field, Surface Current (re+im).
159 : std::unique_ptr<mfem::VectorCoefficient> S, E_sr, E_si, B_sr, B_si, A_s, J_sr, J_si;
160 :
161 : // Electric Energy Density, Magnetic Energy Density, Scalar Potential Boundary Field,
162 : // Surface Charge (re+im).
163 : std::unique_ptr<mfem::Coefficient> U_e, U_m, V_s, Q_sr, Q_si;
164 :
165 : // Wave port boundary mode field postprocessing.
166 0 : struct WavePortFieldData
167 : {
168 : std::unique_ptr<mfem::VectorCoefficient> E0r, E0i;
169 : };
170 : std::map<int, WavePortFieldData> port_E0;
171 :
172 : // Setup coefficients for field postprocessing.
173 : void SetupFieldCoefficients();
174 :
175 : // Initialize Paraview, register all fields to write.
176 : void InitializeParaviewDataCollection(const fs::path &sub_folder_name = "");
177 :
178 : public:
179 : // Public overload for the driven solver only, that takes in an excitation index and
180 : // sets the correct sub_folder_name path for the primary function above.
181 : template <ProblemType U = solver_t>
182 : auto InitializeParaviewDataCollection(int ex_idx)
183 : -> std::enable_if_t<U == ProblemType::DRIVEN, void>;
184 :
185 : protected:
186 : // Write to disk the E- and B-fields extracted from the solution vectors. Note that
187 : // fields are not redimensionalized, to do so one needs to compute: B <= B * (μ₀ H₀), E
188 : // <= E * (Z₀ H₀), V <= V * (Z₀ H₀ L₀), etc.
189 : void WriteParaviewFields(double time, int step);
190 : void WriteParaviewFieldsFinal(const ErrorIndicator *indicator = nullptr);
191 : void WriteMFEMGridFunctions(double time, int step);
192 : void WriteMFEMGridFunctionsFinal(const ErrorIndicator *indicator = nullptr);
193 :
194 : // CSV Measure & Print.
195 :
196 : // PostOperatorCSV<solver_t> is a class that contains csv tables and printers of
197 : // measurements. Conceptually, its members could be a part of this class, like the
198 : // ParaView fields and functions above. It has been separated out for code readability. To
199 : // achieve this, it is has a pointer back to its "parent" PostOperator class and is a
200 : // friend class so it can access the private measurement_cache and references of the
201 : // system from fem_op.
202 : friend PostOperatorCSV<solver_t>;
203 :
204 : PostOperatorCSV<solver_t> post_op_csv;
205 :
206 : // Helper classes that actually do some measurements that will be saved to csv files.
207 :
208 : DomainPostOperator dom_post_op; // Energy in bulk
209 : SurfacePostOperator surf_post_op; // Dielectric Interface Energy, Flux, and FarField
210 : mutable InterpolationOperator interp_op; // E & B fields: mutates during measure
211 :
212 : mutable Measurement measurement_cache;
213 :
214 : // Individual measurements to fill the cache/workspace. Measurements functions are not
215 : // constrained by solver type in the signature since they are private member functions.
216 : // They dispatch on solver type within the function itself using `if constexpr`, and do
217 : // nothing if the measurement is not solver appropriate.
218 : void MeasureDomainFieldEnergy() const;
219 : void MeasureLumpedPorts() const;
220 : void MeasureWavePorts() const;
221 : void MeasureLumpedPortsEig() const; // Depends: DomainFieldEnergy, LumpedPorts
222 : void MeasureSParameter() const; // Depends: LumpedPorts, WavePorts
223 : void MeasureSurfaceFlux() const;
224 : void MeasureFarField() const;
225 : void MeasureInterfaceEFieldEnergy() const; // Depends: LumpedPorts
226 : void MeasureProbes() const;
227 :
228 : // Helper function called by all solvers. Has to ensure correct call order to deal with
229 : // dependant measurements.
230 15 : void MeasureAllImpl() const
231 : {
232 15 : MeasureDomainFieldEnergy();
233 9 : MeasureLumpedPorts();
234 3 : MeasureWavePorts();
235 3 : MeasureLumpedPortsEig();
236 3 : MeasureSParameter();
237 15 : MeasureSurfaceFlux();
238 12 : MeasureInterfaceEFieldEnergy();
239 15 : MeasureProbes();
240 6 : MeasureFarField();
241 15 : }
242 :
243 : // Setting grid functions.
244 : //
245 : // Populate the grid function solutions for the E- and B-field using the solution vectors
246 : // on the true dofs. For the real-valued overload, the electric scalar potential can be
247 : // specified too for electrostatic simulations. The output mesh and fields are
248 : // non-dimensionalized consistently (B ~ E (L₀ ω₀ E₀⁻¹)).
249 : //
250 : // These functions are private helper functions. We want to enforce that a caller passes
251 : // the appropriate ones as part of the MeasureAndPrintAll interface, rather than do a
252 : // runtime check to see that they have been set.
253 : //
254 : // TODO(C++20): Switch SFINAE to requires.
255 :
256 : template <ProblemType U = solver_t>
257 6 : auto SetEGridFunction(const ComplexVector &e, bool exchange_face_nbr_data = true)
258 : -> std::enable_if_t<HasEGridFunction<U>() && HasComplexGridFunction<U>(), void>
259 : {
260 6 : E->Real().SetFromTrueDofs(e.Real()); // Parallel distribute
261 6 : E->Imag().SetFromTrueDofs(e.Imag());
262 6 : if (exchange_face_nbr_data)
263 : {
264 6 : E->Real().ExchangeFaceNbrData(); // Ready for parallel comm on shared faces
265 6 : E->Imag().ExchangeFaceNbrData();
266 : }
267 6 : }
268 :
269 : template <ProblemType U = solver_t>
270 : auto SetEGridFunction(const Vector &e, bool exchange_face_nbr_data = true)
271 : -> std::enable_if_t<HasEGridFunction<U>() && !HasComplexGridFunction<U>(), void>
272 : {
273 6 : E->Real().SetFromTrueDofs(e);
274 : if (exchange_face_nbr_data)
275 : {
276 6 : E->Real().ExchangeFaceNbrData();
277 : }
278 : }
279 :
280 : template <ProblemType U = solver_t>
281 6 : auto SetBGridFunction(const ComplexVector &b, bool exchange_face_nbr_data = true)
282 : -> std::enable_if_t<HasBGridFunction<U>() && HasComplexGridFunction<U>(), void>
283 : {
284 6 : B->Real().SetFromTrueDofs(b.Real()); // Parallel distribute
285 6 : B->Imag().SetFromTrueDofs(b.Imag());
286 6 : if (exchange_face_nbr_data)
287 : {
288 6 : B->Real().ExchangeFaceNbrData(); // Ready for parallel comm on shared faces
289 6 : B->Imag().ExchangeFaceNbrData();
290 : }
291 6 : }
292 :
293 : template <ProblemType U = solver_t>
294 : auto SetBGridFunction(const Vector &b, bool exchange_face_nbr_data = true)
295 : -> std::enable_if_t<HasBGridFunction<U>() && !HasComplexGridFunction<U>(), void>
296 : {
297 6 : B->Real().SetFromTrueDofs(b);
298 : if (exchange_face_nbr_data)
299 : {
300 6 : B->Real().ExchangeFaceNbrData();
301 : }
302 : }
303 :
304 : template <ProblemType U = solver_t>
305 : auto SetVGridFunction(const Vector &v, bool exchange_face_nbr_data = true)
306 : -> std::enable_if_t<HasVGridFunction<U>() && !HasComplexGridFunction<U>(), void>
307 : {
308 3 : V->Real().SetFromTrueDofs(v);
309 : if (exchange_face_nbr_data)
310 : {
311 3 : V->Real().ExchangeFaceNbrData();
312 : }
313 : }
314 :
315 : template <ProblemType U = solver_t>
316 : auto SetAGridFunction(const Vector &a, bool exchange_face_nbr_data = true)
317 : -> std::enable_if_t<HasAGridFunction<U>() && !HasComplexGridFunction<U>(), void>
318 : {
319 3 : A->Real().SetFromTrueDofs(a);
320 : if (exchange_face_nbr_data)
321 : {
322 3 : A->Real().ExchangeFaceNbrData();
323 : }
324 : }
325 :
326 : public:
327 : explicit PostOperator(const IoData &iodata, fem_op_t<solver_t> &fem_op);
328 :
329 : // MeasureAndPrintAll is the primary public interface of this class. It is specialized by
330 : // solver type, since each solver has different fields and extra data required. These
331 : // functions all:
332 : // 1) Set the GridFunctions which have to be passed as part of the call.
333 : // 2) Perform all measurements and populate measurement_cache with temporary results. This
334 : // cache structure exists since measurements have dependencies; we may use some
335 : // measurement results in later measurements.
336 : // 3) Pass the measurement cache to the csv printer which will add the appropriate
337 : // rows/cols to the csv tables and print to file.
338 : // 4) Trigger ParaView field computation and save.
339 : //
340 : // The functions return the total domain energy which is the only thing needed in the
341 : // solver to normalize the error indicator. If more measurements were needed by the solver
342 : // loop, we could imagine passing a small struct (like Measurement above or some sub-set
343 : // therefore).
344 : //
345 : // The measure functions will also do logging of (some) measurements to stdout.
346 : //
347 : // TODO(C++20): Upgrade SFINAE to C++20 concepts to simplify static selection since we can
348 : // just write `MeasureAndPrintAll(...) requires (solver_t == Type::A)` without extra
349 : // template.
350 :
351 : template <ProblemType U = solver_t>
352 : auto MeasureAndPrintAll(int ex_idx, int step, const ComplexVector &e,
353 : const ComplexVector &b, std::complex<double> omega)
354 : -> std::enable_if_t<U == ProblemType::DRIVEN, double>;
355 :
356 : template <ProblemType U = solver_t>
357 : auto MeasureAndPrintAll(int step, const ComplexVector &e, const ComplexVector &b,
358 : std::complex<double> omega, double error_abs, double error_bkwd,
359 : int num_conv)
360 : -> std::enable_if_t<U == ProblemType::EIGENMODE, double>;
361 :
362 : template <ProblemType U = solver_t>
363 : auto MeasureAndPrintAll(int step, const Vector &v, const Vector &e, int idx)
364 : -> std::enable_if_t<U == ProblemType::ELECTROSTATIC, double>;
365 :
366 : template <ProblemType U = solver_t>
367 : auto MeasureAndPrintAll(int step, const Vector &a, const Vector &b, int idx)
368 : -> std::enable_if_t<U == ProblemType::MAGNETOSTATIC, double>;
369 :
370 : template <ProblemType U = solver_t>
371 : auto MeasureAndPrintAll(int step, const Vector &e, const Vector &b, double t,
372 : double J_coef)
373 : -> std::enable_if_t<U == ProblemType::TRANSIENT, double>;
374 :
375 : // Write error indicator into ParaView file and print summary statistics to csv. Should be
376 : // called once at the end of the solver loop.
377 : void MeasureFinalize(const ErrorIndicator &indicator);
378 :
379 : // Measurement of the domain energy without printing. This is needed during the driven
380 : // simulation with PROM. There samples are taken and we need the total domain energy for
381 : // the error indicator, but no other measurement / printing should be done.
382 : //
383 : // TODO(C++20): SFINAE to requires.
384 : template <ProblemType U = solver_t>
385 : auto MeasureDomainFieldEnergyOnly(const ComplexVector &e, const ComplexVector &b)
386 : -> std::enable_if_t<U == ProblemType::DRIVEN, double>;
387 :
388 : // Access grid functions for field solutions. Note that these are NOT const functions. The
389 : // electrostatics / magnetostatics solver do measurements of the capacitance/ inductance
390 : // matrix globally at the end of all solves. This is done in the solver class, but uses
391 : // the GridFunctions in this (PostOp) class as already allocated scratch workspace.
392 : //
393 : // Future: Consider moving those cap/ind measurements into this class and MeasureFinalize?
394 : // Would need to store vector of V,A.
395 : //
396 : // TODO(C++20): Switch SFINAE to requires.
397 : template <ProblemType U = solver_t>
398 : auto GetEGridFunction() -> std::enable_if_t<HasEGridFunction<U>(), decltype(*E) &>
399 : {
400 : return *E;
401 : }
402 :
403 : template <ProblemType U = solver_t>
404 : auto GetBGridFunction() -> std::enable_if_t<HasBGridFunction<U>(), decltype(*B) &>
405 : {
406 : return *B;
407 : }
408 :
409 : template <ProblemType U = solver_t>
410 : auto GetVGridFunction() -> std::enable_if_t<HasVGridFunction<U>(), decltype(*V) &>
411 : {
412 : return *V;
413 : }
414 :
415 : template <ProblemType U = solver_t>
416 : auto GetAGridFunction() -> std::enable_if_t<HasAGridFunction<U>(), decltype(*A) &>
417 : {
418 : return *A;
419 : }
420 :
421 : // Access to number of padding digits.
422 15 : constexpr auto GetPadDigitsDefault() const { return pad_digits_default; }
423 :
424 : // Access to domain postprocessing objects. Use in electrostatic & magnetostatic matrix
425 : // measurement (see above).
426 0 : const auto &GetDomainPostOp() const { return dom_post_op; }
427 :
428 : // Expose MPI communicator from fem_op for electrostatic & magnetostatic matrix processing
429 : // (see above).
430 0 : auto GetComm() const { return fem_op->GetComm(); }
431 : };
432 :
433 : } // namespace palace
434 :
435 : #endif // PALACE_MODELS_POST_OPERATOR_HPP
|
