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 "postoperatorcsv.hpp"
5 :
6 : #include <mfem.hpp>
7 :
8 : #include "models/curlcurloperator.hpp"
9 : #include "models/laplaceoperator.hpp"
10 : #include "models/materialoperator.hpp"
11 : #include "models/postoperator.hpp"
12 : #include "models/spaceoperator.hpp"
13 : #include "utils/iodata.hpp"
14 : #include "utils/timer.hpp"
15 :
16 : namespace palace
17 : {
18 :
19 : // static
20 11 : Measurement Measurement::Dimensionalize(const Units &units,
21 : const Measurement &nondim_measurement_cache)
22 : {
23 11 : Measurement measurement_cache;
24 : measurement_cache.freq =
25 11 : units.Dimensionalize<Units::ValueType::FREQUENCY>(nondim_measurement_cache.freq);
26 11 : measurement_cache.ex_idx = nondim_measurement_cache.ex_idx; // NONE
27 11 : measurement_cache.Jcoeff_excitation = nondim_measurement_cache.Jcoeff_excitation; // NONE
28 11 : measurement_cache.eigenmode_Q = nondim_measurement_cache.eigenmode_Q; // NONE
29 11 : measurement_cache.error_abs = nondim_measurement_cache.error_abs; // NONE
30 11 : measurement_cache.error_bkwd = nondim_measurement_cache.error_bkwd; // NONE
31 :
32 11 : measurement_cache.domain_E_field_energy_all =
33 : units.Dimensionalize<Units::ValueType::ENERGY>(
34 11 : nondim_measurement_cache.domain_E_field_energy_all);
35 11 : measurement_cache.domain_H_field_energy_all =
36 : units.Dimensionalize<Units::ValueType::ENERGY>(
37 11 : nondim_measurement_cache.domain_H_field_energy_all);
38 16 : for (const auto &e : nondim_measurement_cache.domain_E_field_energy_i)
39 : {
40 5 : measurement_cache.domain_E_field_energy_i.emplace_back(Measurement::DomainData{
41 5 : e.idx, units.Dimensionalize<Units::ValueType::ENERGY>(e.energy),
42 5 : e.participation_ratio});
43 : }
44 16 : for (const auto &e : nondim_measurement_cache.domain_H_field_energy_i)
45 : {
46 5 : measurement_cache.domain_H_field_energy_i.emplace_back(Measurement::DomainData{
47 5 : e.idx, units.Dimensionalize<Units::ValueType::ENERGY>(e.energy),
48 5 : e.participation_ratio});
49 : }
50 11 : measurement_cache.lumped_port_capacitor_energy =
51 : units.Dimensionalize<Units::ValueType::ENERGY>(
52 11 : nondim_measurement_cache.lumped_port_capacitor_energy);
53 11 : measurement_cache.lumped_port_inductor_energy =
54 : units.Dimensionalize<Units::ValueType::ENERGY>(
55 11 : nondim_measurement_cache.lumped_port_inductor_energy);
56 :
57 : auto dimensionalize_port_post_data =
58 22 : [&units](const std::map<int, Measurement::PortPostData> &nondim)
59 : {
60 : std::map<int, Measurement::PortPostData> dim;
61 34 : for (const auto &[k, data] : nondim)
62 : {
63 12 : dim[k] = Measurement::PortPostData();
64 12 : dim[k].P = units.Dimensionalize<Units::ValueType::POWER>(data.P);
65 12 : dim[k].V = units.Dimensionalize<Units::ValueType::VOLTAGE>(data.V),
66 12 : dim[k].I = units.Dimensionalize<Units::ValueType::CURRENT>(data.I),
67 12 : dim[k].I_RLC = {units.Dimensionalize<Units::ValueType::CURRENT>(data.I_RLC[0]),
68 : units.Dimensionalize<Units::ValueType::CURRENT>(data.I_RLC[1]),
69 : units.Dimensionalize<Units::ValueType::CURRENT>(data.I_RLC[2])};
70 12 : dim[k].S = data.S; // NONE
71 :
72 12 : dim[k].inductor_energy =
73 12 : units.Dimensionalize<Units::ValueType::ENERGY>(data.inductor_energy);
74 12 : dim[k].capacitor_energy =
75 12 : units.Dimensionalize<Units::ValueType::ENERGY>(data.capacitor_energy);
76 :
77 12 : dim[k].mode_port_kappa =
78 12 : units.Dimensionalize<Units::ValueType::FREQUENCY>(data.mode_port_kappa);
79 12 : dim[k].quality_factor = data.quality_factor; // NONE
80 12 : dim[k].inductive_energy_participation = data.inductive_energy_participation; // NONE
81 : }
82 22 : return dim;
83 11 : };
84 : measurement_cache.lumped_port_vi =
85 11 : dimensionalize_port_post_data(nondim_measurement_cache.lumped_port_vi);
86 : measurement_cache.wave_port_vi =
87 11 : dimensionalize_port_post_data(nondim_measurement_cache.wave_port_vi);
88 :
89 22 : measurement_cache.probe_E_field = units.Dimensionalize<Units::ValueType::FIELD_E>(
90 11 : nondim_measurement_cache.probe_E_field);
91 11 : measurement_cache.probe_B_field = units.Dimensionalize<Units::ValueType::FIELD_B>(
92 11 : nondim_measurement_cache.probe_B_field);
93 :
94 16 : for (const auto &data : nondim_measurement_cache.surface_flux_i)
95 : {
96 5 : auto &flux = measurement_cache.surface_flux_i.emplace_back(data);
97 5 : if (data.type == SurfaceFlux::ELECTRIC)
98 : {
99 : flux.Phi *= units.GetScaleFactor<Units::ValueType::CAPACITANCE>();
100 : flux.Phi *= units.GetScaleFactor<Units::ValueType::VOLTAGE>();
101 : }
102 3 : else if (data.type == SurfaceFlux::MAGNETIC)
103 : {
104 : flux.Phi *= units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
105 : flux.Phi *= units.GetScaleFactor<Units::ValueType::CURRENT>();
106 : }
107 3 : else if (data.type == SurfaceFlux::POWER)
108 : {
109 : flux.Phi *= units.GetScaleFactor<Units::ValueType::POWER>();
110 : }
111 : }
112 :
113 16 : for (const auto &data : nondim_measurement_cache.interface_eps_i)
114 : {
115 5 : auto &eps = measurement_cache.interface_eps_i.emplace_back(data);
116 5 : eps.energy = units.Dimensionalize<Units::ValueType::ENERGY>(data.energy);
117 : }
118 :
119 : measurement_cache.farfield.thetaphis =
120 11 : nondim_measurement_cache.farfield.thetaphis; // NONE
121 22 : measurement_cache.farfield.E_field = units.Nondimensionalize<Units::ValueType::FIELD_E>(
122 11 : nondim_measurement_cache.farfield.E_field);
123 :
124 11 : return measurement_cache;
125 0 : }
126 :
127 : // static.
128 1 : Measurement Measurement::Nondimensionalize(const Units &units,
129 : const Measurement &dim_measurement_cache)
130 : {
131 1 : Measurement measurement_cache;
132 : measurement_cache.freq =
133 1 : units.Nondimensionalize<Units::ValueType::FREQUENCY>(dim_measurement_cache.freq);
134 1 : measurement_cache.ex_idx = dim_measurement_cache.ex_idx; // NONE
135 1 : measurement_cache.Jcoeff_excitation = dim_measurement_cache.Jcoeff_excitation; // NONE
136 1 : measurement_cache.eigenmode_Q = dim_measurement_cache.eigenmode_Q; // NONE
137 1 : measurement_cache.error_abs = dim_measurement_cache.error_abs; // NONE
138 1 : measurement_cache.error_bkwd = dim_measurement_cache.error_bkwd; // NONE
139 :
140 1 : measurement_cache.domain_E_field_energy_all =
141 : units.Nondimensionalize<Units::ValueType::ENERGY>(
142 1 : dim_measurement_cache.domain_E_field_energy_all);
143 1 : measurement_cache.domain_H_field_energy_all =
144 : units.Nondimensionalize<Units::ValueType::ENERGY>(
145 1 : dim_measurement_cache.domain_H_field_energy_all);
146 6 : for (const auto &e : dim_measurement_cache.domain_E_field_energy_i)
147 : {
148 5 : measurement_cache.domain_E_field_energy_i.emplace_back(Measurement::DomainData{
149 5 : e.idx, units.Nondimensionalize<Units::ValueType::ENERGY>(e.energy),
150 5 : e.participation_ratio});
151 : }
152 6 : for (const auto &e : dim_measurement_cache.domain_H_field_energy_i)
153 : {
154 5 : measurement_cache.domain_H_field_energy_i.emplace_back(Measurement::DomainData{
155 5 : e.idx, units.Nondimensionalize<Units::ValueType::ENERGY>(e.energy),
156 5 : e.participation_ratio});
157 : }
158 1 : measurement_cache.lumped_port_capacitor_energy =
159 : units.Nondimensionalize<Units::ValueType::ENERGY>(
160 1 : dim_measurement_cache.lumped_port_capacitor_energy);
161 1 : measurement_cache.lumped_port_inductor_energy =
162 : units.Nondimensionalize<Units::ValueType::ENERGY>(
163 1 : dim_measurement_cache.lumped_port_inductor_energy);
164 :
165 : auto dimensionalize_port_post_data =
166 2 : [&units](const std::map<int, Measurement::PortPostData> &nondim)
167 : {
168 : std::map<int, Measurement::PortPostData> dim;
169 10 : for (const auto &[k, data] : nondim)
170 : {
171 8 : dim[k] = Measurement::PortPostData();
172 8 : dim[k].P = units.Nondimensionalize<Units::ValueType::POWER>(data.P);
173 8 : dim[k].V = units.Nondimensionalize<Units::ValueType::VOLTAGE>(data.V),
174 8 : dim[k].I = units.Nondimensionalize<Units::ValueType::CURRENT>(data.I),
175 8 : dim[k].I_RLC = {units.Nondimensionalize<Units::ValueType::CURRENT>(data.I_RLC[0]),
176 : units.Nondimensionalize<Units::ValueType::CURRENT>(data.I_RLC[1]),
177 : units.Nondimensionalize<Units::ValueType::CURRENT>(data.I_RLC[2])};
178 8 : dim[k].S = data.S; // NONE
179 :
180 8 : dim[k].inductor_energy =
181 8 : units.Nondimensionalize<Units::ValueType::ENERGY>(data.inductor_energy);
182 8 : dim[k].capacitor_energy =
183 8 : units.Nondimensionalize<Units::ValueType::ENERGY>(data.capacitor_energy);
184 :
185 8 : dim[k].mode_port_kappa =
186 8 : units.Nondimensionalize<Units::ValueType::FREQUENCY>(data.mode_port_kappa);
187 8 : dim[k].quality_factor = data.quality_factor; // NONE
188 8 : dim[k].inductive_energy_participation = data.inductive_energy_participation; // NONE
189 : }
190 2 : return dim;
191 1 : };
192 : measurement_cache.lumped_port_vi =
193 1 : dimensionalize_port_post_data(dim_measurement_cache.lumped_port_vi);
194 : measurement_cache.wave_port_vi =
195 1 : dimensionalize_port_post_data(dim_measurement_cache.wave_port_vi);
196 :
197 2 : measurement_cache.probe_E_field = units.Nondimensionalize<Units::ValueType::FIELD_E>(
198 1 : dim_measurement_cache.probe_E_field);
199 1 : measurement_cache.probe_B_field = units.Nondimensionalize<Units::ValueType::FIELD_B>(
200 1 : dim_measurement_cache.probe_B_field);
201 :
202 6 : for (const auto &data : dim_measurement_cache.surface_flux_i)
203 : {
204 5 : auto &flux = measurement_cache.surface_flux_i.emplace_back(data);
205 5 : if (data.type == SurfaceFlux::ELECTRIC)
206 : {
207 : flux.Phi /= units.GetScaleFactor<Units::ValueType::CAPACITANCE>();
208 : flux.Phi /= units.GetScaleFactor<Units::ValueType::VOLTAGE>();
209 : }
210 3 : else if (data.type == SurfaceFlux::MAGNETIC)
211 : {
212 : flux.Phi /= units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
213 : flux.Phi /= units.GetScaleFactor<Units::ValueType::CURRENT>();
214 : }
215 3 : else if (data.type == SurfaceFlux::POWER)
216 : {
217 : flux.Phi /= units.GetScaleFactor<Units::ValueType::POWER>();
218 : }
219 : }
220 :
221 6 : for (const auto &data : dim_measurement_cache.interface_eps_i)
222 : {
223 5 : auto &eps = measurement_cache.interface_eps_i.emplace_back(data);
224 5 : eps.energy = units.Nondimensionalize<Units::ValueType::ENERGY>(data.energy);
225 : }
226 :
227 1 : measurement_cache.farfield.thetaphis = dim_measurement_cache.farfield.thetaphis; // NONE
228 2 : measurement_cache.farfield.E_field = units.Nondimensionalize<Units::ValueType::FIELD_E>(
229 1 : dim_measurement_cache.farfield.E_field);
230 :
231 1 : return measurement_cache;
232 0 : }
233 :
234 : namespace
235 : {
236 :
237 : // TODO(C++20): Do constexpr with string.
238 0 : std::string DimLabel(int i)
239 : {
240 0 : switch (i)
241 : {
242 : // Note: Zero-based indexing here.
243 : case 0:
244 0 : return "x";
245 : case 1:
246 0 : return "y";
247 : case 2:
248 0 : return "z";
249 : default:
250 : return fmt::format("d{}", i);
251 : }
252 : }
253 :
254 : // TODO(C++20): Do constexpr with string.
255 : std::string LabelIndexCol(const ProblemType solver_t)
256 : {
257 : switch (solver_t)
258 : {
259 : case ProblemType::DRIVEN:
260 76 : return "f (GHz)";
261 : case ProblemType::EIGENMODE:
262 4 : return "m";
263 : case ProblemType::ELECTROSTATIC:
264 : case ProblemType::MAGNETOSTATIC:
265 8 : return "i";
266 : case ProblemType::TRANSIENT:
267 12 : return "t (ns)";
268 : default:
269 : return "unknown";
270 : }
271 : }
272 : int PrecIndexCol(const ProblemType solver_t)
273 : {
274 : switch (solver_t)
275 : {
276 : case ProblemType::DRIVEN:
277 : case ProblemType::TRANSIENT:
278 : return 8;
279 : case ProblemType::EIGENMODE:
280 : case ProblemType::ELECTROSTATIC:
281 : case ProblemType::MAGNETOSTATIC:
282 : return 2;
283 : default:
284 : return 8;
285 : }
286 : }
287 :
288 : // Index checking when adding to a new excitation block: When adding data to data_col with
289 : // index idx, checks that idx matches what is already written in the corresponding row of
290 : // idx_col. Adds a new idx row to idx_col if needed.
291 20 : void CheckAppendIndex(Column &idx_col, double idx_value, size_t m_idx_row)
292 : {
293 20 : if (m_idx_row == idx_col.n_rows())
294 : {
295 20 : idx_col << idx_value;
296 : }
297 : else
298 : {
299 0 : auto current_idx = idx_col.data.at(m_idx_row);
300 0 : MFEM_VERIFY(idx_value == current_idx,
301 : fmt::format("Writing data table at incorrect index. Data has index {} "
302 : "while table is at {}",
303 : idx_value, current_idx));
304 : }
305 20 : }
306 :
307 : } // namespace
308 :
309 16 : std::vector<std::size_t> _impl::table_expected_filling(std::size_t m_idx_row,
310 : std::size_t ex_idx_i,
311 : std::size_t nr_rows,
312 : std::size_t nr_col_blocks)
313 : {
314 : // Expected column group filling pattern. Include leading index (freq, ...)
315 16 : std::vector<std::size_t> filling_pattern(nr_col_blocks + 1, 0);
316 16 : filling_pattern.at(0) = (ex_idx_i == 0) ? m_idx_row : nr_rows; // index column
317 23 : for (std::size_t i = 1; i < ex_idx_i + 1; i++)
318 : {
319 7 : filling_pattern.at(i) = nr_rows;
320 : }
321 16 : filling_pattern.at(ex_idx_i + 1) = m_idx_row;
322 16 : return filling_pattern;
323 : }
324 :
325 : template <ProblemType solver_t>
326 44 : void PostOperatorCSV<solver_t>::MoveTableValidateReload(TableWithCSVFile &t_csv_base,
327 : Table &&t_ref)
328 : {
329 : // For non-driven solvers or driven with default restart, no table was loaded.
330 44 : if (!reload_table)
331 : {
332 30 : t_csv_base.table = std::move(t_ref);
333 30 : return;
334 : }
335 :
336 : // At this point we have a non-default restart. We need to verify that (a) the structure
337 : // of the table is valid, (b) the cursor location matches the expected restart location.
338 :
339 : auto file = t_csv_base.get_csv_filepath();
340 : Table &t_base = t_csv_base.table;
341 17 : MFEM_VERIFY(!t_base.empty(),
342 : fmt::format("The data table loaded from path {} was empty, but the "
343 : "simulation expected a restart with existing data!",
344 : file));
345 :
346 : auto err_msg = fmt::format("The results table loaded from path {} contains pre-existing "
347 : "data, but it doest not match the "
348 : "expected table structure.",
349 : file);
350 13 : if (t_base.n_cols() != t_ref.n_cols())
351 : {
352 4 : MFEM_ABORT(fmt::format("{} [Mismatched number of columns: expected {}, got {}.]",
353 : err_msg, t_base.n_cols(), t_ref.n_cols()))
354 : }
355 : std::vector<std::size_t> base_ex_idx_nrows;
356 : long current_ex_idx_v = std::numeric_limits<long>::max();
357 141 : for (std::size_t i = 0; i < t_base.n_cols(); i++)
358 : {
359 : auto &t_base_i = t_base[i];
360 : auto &t_ref_i = t_ref[i];
361 :
362 131 : if (t_base_i.header_text != t_ref_i.header_text)
363 : {
364 4 : MFEM_ABORT(fmt::format("{} [Mismatched column header: expected {}, got {}.]", err_msg,
365 : t_base_i.header_text, t_ref_i.header_text))
366 : }
367 : // Since we cannot parse the column name (internal label) from the printed csv file or
368 : // other options from csv file, we will over-write them from t_ref.
369 130 : t_base_i.name = t_ref_i.name;
370 130 : t_base_i.column_group_idx = t_ref_i.column_group_idx;
371 130 : t_base_i.min_left_padding = t_ref_i.min_left_padding;
372 130 : t_base_i.float_precision = t_ref_i.float_precision;
373 : t_base_i.fmt_sign = t_ref_i.fmt_sign;
374 130 : t_base_i.print_as_int = t_ref_i.print_as_int;
375 :
376 : // Check that columns in same group have the same row number. Assumes that column groups
377 : // are contiguous. If no error, save row number to compare to expected pattern.
378 130 : if (t_base_i.column_group_idx != current_ex_idx_v)
379 : {
380 : current_ex_idx_v = t_base_i.column_group_idx;
381 36 : base_ex_idx_nrows.push_back(t_base_i.n_rows());
382 : }
383 : else
384 : {
385 100 : if (t_base_i.n_rows() != base_ex_idx_nrows.back())
386 : {
387 4 : MFEM_ABORT(fmt::format("{} [Mismatched rows in excitation {}.]", err_msg,
388 : current_ex_idx_v))
389 : }
390 : }
391 : }
392 : // Match expected column group pattern.
393 10 : auto expected_ex_idx_nrows = _impl::table_expected_filling(
394 : row_i, ex_idx_i, nr_expected_measurement_rows, ex_idx_v_all.size());
395 :
396 : // Copy over other options from reference table, since we don't recover them on load.
397 10 : t_csv_base.table.col_options = t_ref.col_options;
398 :
399 22 : MFEM_VERIFY(base_ex_idx_nrows == expected_ex_idx_nrows,
400 : fmt::format("{} [Specified restart position is incompatible with reloaded "
401 : "file. Row filling by excitation expected {}, got {}]",
402 : err_msg, expected_ex_idx_nrows, base_ex_idx_nrows))
403 :
404 : // Don't check index column (frequency) values or size. Size should match with sizing from
405 : // cursor from printer below. Values will be checked as new frequencies are written.
406 : }
407 :
408 : template <ProblemType solver_t>
409 10 : void PostOperatorCSV<solver_t>::InitializeDomainE(const DomainPostOperator &dom_post_op)
410 : {
411 : using fmt::format;
412 30 : domain_E = TableWithCSVFile(post_dir / "domain-E.csv", reload_table);
413 :
414 10 : Table t; // Define table locally first due to potential reload.
415 10 : auto nr_expected_measurement_cols =
416 10 : 1 + ex_idx_v_all.size() * 4 * (1 + dom_post_op.M_i.size());
417 10 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
418 20 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
419 20 : for (const auto ex_idx : ex_idx_v_all)
420 : {
421 10 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
422 :
423 20 : t.insert(format("Ee_{}", ex_idx), format("E_elec{} (J)", ex_label), ex_idx);
424 20 : t.insert(format("Em_{}", ex_idx), format("E_mag{} (J)", ex_label), ex_idx);
425 20 : t.insert(format("Ec_{}", ex_idx), format("E_cap{} (J)", ex_label), ex_idx);
426 20 : t.insert(format("Ei_{}", ex_idx), format("E_ind{} (J)", ex_label), ex_idx);
427 :
428 10 : for (const auto &[idx, data] : dom_post_op.M_i)
429 : {
430 0 : t.insert(format("Ee_{}_{}", idx, ex_idx), format("E_elec[{}]{} (J)", idx, ex_label),
431 : ex_idx);
432 0 : t.insert(format("pe_{}_{}", idx, ex_idx), format("p_elec[{}]{}", idx, ex_label),
433 : ex_idx);
434 0 : t.insert(format("Em_{}_{}", idx, ex_idx), format("E_mag[{}]{} (J)", idx, ex_label),
435 : ex_idx);
436 0 : t.insert(format("pm_{}_{}", idx, ex_idx), format("p_mag[{}]{}", idx, ex_label),
437 : ex_idx);
438 : }
439 : }
440 10 : MoveTableValidateReload(*domain_E, std::move(t));
441 : // No longer WriteFullTableTrunc here. We want to potentially reload and check all
442 : // measurement tables, which may have existing values, before overwriting anything. Just
443 : // write on first measurement.
444 10 : }
445 :
446 : template <ProblemType solver_t>
447 10 : void PostOperatorCSV<solver_t>::PrintDomainE()
448 : {
449 10 : if (!domain_E) // trivial check: always written and we are always on root
450 : {
451 : return;
452 : }
453 : using fmt::format;
454 10 : CheckAppendIndex(domain_E->table["idx"], row_idx_v, row_i);
455 20 : domain_E->table[format("Ee_{}", m_ex_idx)] << measurement_cache.domain_E_field_energy_all;
456 20 : domain_E->table[format("Em_{}", m_ex_idx)] << measurement_cache.domain_H_field_energy_all;
457 10 : domain_E->table[format("Ec_{}", m_ex_idx)]
458 20 : << measurement_cache.lumped_port_capacitor_energy;
459 10 : domain_E->table[format("Ei_{}", m_ex_idx)]
460 20 : << measurement_cache.lumped_port_inductor_energy;
461 10 : for (const auto &data : measurement_cache.domain_E_field_energy_i)
462 : {
463 0 : domain_E->table[format("Ee_{}_{}", data.idx, m_ex_idx)] << data.energy;
464 0 : domain_E->table[format("pe_{}_{}", data.idx, m_ex_idx)] << data.participation_ratio;
465 : }
466 10 : for (const auto &data : measurement_cache.domain_H_field_energy_i)
467 : {
468 0 : domain_E->table[format("Em_{}_{}", data.idx, m_ex_idx)] << data.energy;
469 0 : domain_E->table[format("pm_{}_{}", data.idx, m_ex_idx)] << data.participation_ratio;
470 : }
471 10 : domain_E->WriteFullTableTrunc();
472 : }
473 :
474 : template <ProblemType solver_t>
475 10 : void PostOperatorCSV<solver_t>::InitializeSurfaceF(const SurfacePostOperator &surf_post_op)
476 : {
477 10 : if (!(surf_post_op.flux_surfs.size() > 0))
478 : {
479 10 : return;
480 : }
481 : using fmt::format;
482 0 : surface_F = TableWithCSVFile(post_dir / "surface-F.csv", reload_table);
483 :
484 0 : Table t; // Define table locally first due to potential reload.
485 0 : auto nr_expected_measurement_cols = 1 + ex_idx_v_all.size() *
486 0 : (HasComplexGridFunction<solver_t>() ? 2 : 1) *
487 : surf_post_op.flux_surfs.size();
488 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
489 0 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
490 0 : for (const auto ex_idx : ex_idx_v_all)
491 : {
492 0 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
493 0 : for (const auto &[idx, data] : surf_post_op.flux_surfs)
494 : {
495 0 : switch (data.type)
496 : {
497 : case SurfaceFlux::ELECTRIC:
498 : if (HasComplexGridFunction<solver_t>())
499 : {
500 0 : t.insert(format("F_{}_{}_re", idx, ex_idx),
501 0 : format("Re{{Φ_elec[{}]{}}} (C)", idx, ex_label), ex_idx);
502 0 : t.insert(format("F_{}_{}_im", idx, ex_idx),
503 0 : format("Im{{Φ_elec[{}]{}}} (C)", idx, ex_label), ex_idx);
504 : }
505 : else
506 : {
507 0 : t.insert(format("F_{}_{}_re", idx, ex_idx),
508 0 : format("Φ_elec[{}]{} (C)", idx, ex_label), ex_idx);
509 : }
510 0 : break;
511 : case SurfaceFlux::MAGNETIC:
512 : if (HasComplexGridFunction<solver_t>())
513 : {
514 0 : t.insert(format("F_{}_{}_re", idx, ex_idx),
515 0 : format("Re{{Φ_mag[{}]{}}} (Wb)", idx, ex_label), ex_idx);
516 0 : t.insert(format("F_{}_{}_im", idx, ex_idx),
517 0 : format("Im{{Φ_mag[{}]{}}} (Wb)", idx, ex_label), ex_idx);
518 : }
519 : else
520 : {
521 0 : t.insert(format("F_{}_{}_re", idx, ex_idx),
522 0 : format("Φ_mag[{}]{} (Wb)", idx, ex_label), ex_idx);
523 : }
524 0 : break;
525 : case SurfaceFlux::POWER:
526 0 : t.insert(format("F_{}_{}_re", idx, ex_idx),
527 0 : format("Φ_pow[{}]{} (W)", idx, ex_label), ex_idx);
528 0 : break;
529 : }
530 : }
531 : }
532 0 : MoveTableValidateReload(*surface_F, std::move(t));
533 0 : }
534 :
535 : template <ProblemType solver_t>
536 10 : void PostOperatorCSV<solver_t>::PrintSurfaceF()
537 : {
538 10 : if (!surface_F)
539 : {
540 : return;
541 : }
542 : using fmt::format;
543 0 : CheckAppendIndex(surface_F->table["idx"], row_idx_v, row_i);
544 0 : for (const auto &data : measurement_cache.surface_flux_i)
545 : {
546 0 : surface_F->table[format("F_{}_{}_re", data.idx, m_ex_idx)] << data.Phi.real();
547 0 : if (HasComplexGridFunction<solver_t>() &&
548 0 : (data.type == SurfaceFlux::ELECTRIC || data.type == SurfaceFlux::MAGNETIC))
549 : {
550 0 : surface_F->table[format("F_{}_{}_im", data.idx, m_ex_idx)] << data.Phi.imag();
551 : }
552 : }
553 0 : surface_F->WriteFullTableTrunc();
554 : }
555 :
556 : template <ProblemType solver_t>
557 10 : void PostOperatorCSV<solver_t>::InitializeSurfaceQ(const SurfacePostOperator &surf_post_op)
558 : {
559 10 : if (!(surf_post_op.eps_surfs.size() > 0))
560 : {
561 10 : return;
562 : }
563 : using fmt::format;
564 0 : surface_Q = TableWithCSVFile(post_dir / "surface-Q.csv", reload_table);
565 :
566 0 : Table t; // Define table locally first due to potential reload.
567 0 : auto nr_expected_measurement_cols =
568 0 : 1 + ex_idx_v_all.size() * (2 * surf_post_op.eps_surfs.size());
569 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
570 0 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
571 0 : for (const auto ex_idx : ex_idx_v_all)
572 : {
573 0 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
574 0 : for (const auto &[idx, data] : surf_post_op.eps_surfs)
575 : {
576 0 : t.insert(format("p_{}_{}", idx, ex_idx), format("p_surf[{}]{}", idx, ex_label),
577 : ex_idx);
578 0 : t.insert(format("Q_{}_{}", idx, ex_idx), format("Q_surf[{}]{}", idx, ex_label),
579 : ex_idx);
580 : }
581 : }
582 0 : MoveTableValidateReload(*surface_Q, std::move(t));
583 0 : }
584 :
585 : template <ProblemType solver_t>
586 10 : void PostOperatorCSV<solver_t>::PrintSurfaceQ()
587 : {
588 10 : if (!surface_Q)
589 : {
590 : return;
591 : }
592 : using fmt::format;
593 0 : CheckAppendIndex(surface_Q->table["idx"], row_idx_v, row_i);
594 :
595 0 : for (const auto &data : measurement_cache.interface_eps_i)
596 : {
597 0 : surface_Q->table[format("p_{}_{}", data.idx, m_ex_idx)] << data.energy_participation;
598 0 : surface_Q->table[format("Q_{}_{}", data.idx, m_ex_idx)] << data.quality_factor;
599 : }
600 0 : surface_Q->WriteFullTableTrunc();
601 : }
602 :
603 : template <ProblemType solver_t>
604 : template <ProblemType U>
605 4 : auto PostOperatorCSV<solver_t>::InitializeFarFieldE(const SurfacePostOperator &surf_post_op)
606 : -> std::enable_if_t<U == ProblemType::DRIVEN || U == ProblemType::EIGENMODE, void>
607 : {
608 4 : if (!(surf_post_op.farfield.size() > 0))
609 : {
610 4 : return;
611 : }
612 : using fmt::format;
613 :
614 0 : farfield_E = TableWithCSVFile(post_dir / "farfield-rE.csv");
615 :
616 0 : Table t; // Define table locally first due to potential reload.
617 :
618 : int v_dim = surf_post_op.GetVDim();
619 0 : int scale_col = 2 * v_dim; // Real + Imag components
620 0 : int nr_expected_measurement_cols = 3 + scale_col; // freq, theta, phi
621 0 : int nr_expected_measurement_rows = surf_post_op.farfield.size();
622 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
623 : if constexpr (U == ProblemType::EIGENMODE)
624 : {
625 0 : t.insert("idx", "m", -1, 0, PrecIndexCol(solver_t), "");
626 0 : t.insert("f_re", "f_re (GHz)");
627 0 : t.insert("f_im", "f_im (GHz)");
628 : }
629 : else
630 : {
631 0 : t.insert("idx", "f (GHz)", -1, 0, PrecIndexCol(solver_t), "");
632 : }
633 0 : t.insert(Column("theta", "theta (deg.)", 0, PrecIndexCol(solver_t), {}, ""));
634 0 : t.insert(Column("phi", "phi (deg.)", 0, PrecIndexCol(solver_t), {}, ""));
635 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
636 : {
637 0 : t.insert(format("rE{}_re", i_dim), format("r*Re{{E_{}}} (V)", DimLabel(i_dim)));
638 0 : t.insert(format("rE{}_im", i_dim), format("r*Im{{E_{}}} (V)", DimLabel(i_dim)));
639 : }
640 :
641 0 : MoveTableValidateReload(*farfield_E, std::move(t));
642 0 : }
643 :
644 : template <ProblemType solver_t>
645 : template <ProblemType U>
646 4 : auto PostOperatorCSV<solver_t>::PrintFarFieldE(const SurfacePostOperator &surf_post_op)
647 : -> std::enable_if_t<U == ProblemType::DRIVEN || U == ProblemType::EIGENMODE, void>
648 : {
649 4 : if (!farfield_E)
650 : {
651 : return;
652 : }
653 : using fmt::format;
654 : int v_dim = surf_post_op.GetVDim();
655 0 : for (size_t i = 0; i < measurement_cache.farfield.thetaphis.size(); i++)
656 : {
657 0 : farfield_E->table["idx"] << row_idx_v;
658 : if constexpr (U == ProblemType::EIGENMODE)
659 : {
660 0 : farfield_E->table["f_re"] << measurement_cache.freq.real();
661 0 : farfield_E->table["f_im"] << measurement_cache.freq.imag();
662 : }
663 : const auto &[theta, phi] = measurement_cache.farfield.thetaphis[i];
664 : const auto &E_field = measurement_cache.farfield.E_field[i];
665 :
666 : // Print as degrees instead of radians.
667 0 : farfield_E->table["theta"] << 180 / M_PI * theta;
668 0 : farfield_E->table["phi"] << 180 / M_PI * phi;
669 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
670 : {
671 0 : farfield_E->table[format("rE{}_re", i_dim)] << E_field[i_dim].real();
672 0 : farfield_E->table[format("rE{}_im", i_dim)] << E_field[i_dim].imag();
673 : }
674 : }
675 0 : farfield_E->WriteFullTableTrunc();
676 : }
677 :
678 : template <ProblemType solver_t>
679 8 : void PostOperatorCSV<solver_t>::InitializeProbeE(const InterpolationOperator &interp_op)
680 : {
681 8 : if (!(interp_op.GetProbes().size() > 0) || !HasEGridFunction<solver_t>())
682 : {
683 8 : return;
684 : }
685 : using fmt::format;
686 0 : probe_E = TableWithCSVFile(post_dir / "probe-E.csv", reload_table);
687 :
688 0 : Table t; // Define table locally first due to potential reload.
689 : auto v_dim = interp_op.GetVDim();
690 0 : int scale_col = (HasComplexGridFunction<solver_t>() ? 2 : 1) * v_dim;
691 0 : auto nr_expected_measurement_cols =
692 0 : 1 + ex_idx_v_all.size() * scale_col * interp_op.GetProbes().size();
693 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
694 0 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
695 0 : for (const auto ex_idx : ex_idx_v_all)
696 : {
697 0 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
698 0 : for (const auto &idx : interp_op.GetProbes())
699 : {
700 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
701 : {
702 : if constexpr (HasComplexGridFunction<solver_t>())
703 : {
704 0 : t.insert(format("E{}_{}_{}_re", i_dim, idx, ex_idx),
705 0 : format("Re{{E_{}[{}]{}}} (V/m)", DimLabel(i_dim), idx, ex_label),
706 : ex_idx);
707 0 : t.insert(format("E{}_{}_{}_im", i_dim, idx, ex_idx),
708 0 : format("Im{{E_{}[{}]{}}} (V/m)", DimLabel(i_dim), idx, ex_label),
709 : ex_idx);
710 : }
711 : else
712 : {
713 0 : t.insert(format("E{}_{}_{}_re", i_dim, idx, ex_idx),
714 0 : format("E_{}[{}]{} (V/m)", DimLabel(i_dim), idx, ex_label), ex_idx);
715 : }
716 : }
717 : }
718 : }
719 0 : MoveTableValidateReload(*probe_E, std::move(t));
720 0 : }
721 :
722 : template <ProblemType solver_t>
723 10 : void PostOperatorCSV<solver_t>::PrintProbeE(const InterpolationOperator &interp_op)
724 : {
725 10 : if (!probe_E)
726 : {
727 10 : return;
728 : }
729 : using fmt::format;
730 : auto v_dim = interp_op.GetVDim();
731 0 : auto probe_field = measurement_cache.probe_E_field;
732 0 : MFEM_VERIFY(probe_field.size() == v_dim * interp_op.GetProbes().size(),
733 : format("Size mismatch: expect vector field to have size {} * {} = {}; got {}",
734 : v_dim, interp_op.GetProbes().size(),
735 : v_dim * interp_op.GetProbes().size(), probe_field.size()))
736 :
737 0 : CheckAppendIndex(probe_E->table["idx"], row_idx_v, row_i);
738 :
739 : size_t i = 0;
740 0 : for (const auto &idx : interp_op.GetProbes())
741 : {
742 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
743 : {
744 0 : auto val = probe_field[i * v_dim + i_dim];
745 0 : probe_E->table[format("E{}_{}_{}_re", i_dim, idx, m_ex_idx)] << val.real();
746 : if (HasComplexGridFunction<solver_t>())
747 : {
748 0 : probe_E->table[format("E{}_{}_{}_im", i_dim, idx, m_ex_idx)] << val.imag();
749 : }
750 : }
751 0 : i++;
752 : }
753 0 : probe_E->WriteFullTableTrunc();
754 : }
755 :
756 : template <ProblemType solver_t>
757 8 : void PostOperatorCSV<solver_t>::InitializeProbeB(const InterpolationOperator &interp_op)
758 : {
759 8 : if (!(interp_op.GetProbes().size() > 0) || !HasBGridFunction<solver_t>())
760 : {
761 8 : return;
762 : }
763 : using fmt::format;
764 0 : probe_B = TableWithCSVFile(post_dir / "probe-B.csv", reload_table);
765 0 : Table t; // Define table locally first due to potential reload.
766 : auto v_dim = interp_op.GetVDim();
767 0 : int scale_col = (HasComplexGridFunction<solver_t>() ? 2 : 1) * v_dim;
768 0 : auto nr_expected_measurement_cols =
769 0 : 1 + ex_idx_v_all.size() * scale_col * interp_op.GetProbes().size();
770 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
771 0 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
772 0 : for (const auto ex_idx : ex_idx_v_all)
773 : {
774 0 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
775 0 : for (const auto &idx : interp_op.GetProbes())
776 : {
777 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
778 : {
779 : if (HasComplexGridFunction<solver_t>())
780 : {
781 0 : t.insert(format("B{}_{}_{}_re", i_dim, idx, ex_idx),
782 0 : format("Re{{B_{}[{}]{}}} (Wb/m²)", DimLabel(i_dim), idx, ex_label),
783 : ex_idx);
784 0 : t.insert(format("B{}_{}_{}_im", i_dim, idx, ex_idx),
785 0 : format("Im{{B_{}[{}]{}}} (Wb/m²)", DimLabel(i_dim), idx, ex_label),
786 : ex_idx);
787 : }
788 : else
789 : {
790 0 : t.insert(format("B{}_{}_{}_re", i_dim, idx, ex_idx),
791 0 : format("B_{}[{}]{} (Wb/m²)", DimLabel(i_dim), idx, ex_label), ex_idx);
792 : }
793 : }
794 : }
795 : }
796 0 : MoveTableValidateReload(*probe_B, std::move(t));
797 0 : }
798 :
799 : template <ProblemType solver_t>
800 10 : void PostOperatorCSV<solver_t>::PrintProbeB(const InterpolationOperator &interp_op)
801 : {
802 10 : if (!probe_B)
803 : {
804 10 : return;
805 : }
806 : using fmt::format;
807 :
808 : auto v_dim = interp_op.GetVDim();
809 0 : auto probe_field = measurement_cache.probe_B_field;
810 0 : MFEM_VERIFY(probe_field.size() == v_dim * interp_op.GetProbes().size(),
811 : format("Size mismatch: expect vector field to have size {} * {} = {}; got {}",
812 : v_dim, interp_op.GetProbes().size(),
813 : v_dim * interp_op.GetProbes().size(), probe_field.size()))
814 :
815 0 : CheckAppendIndex(probe_B->table["idx"], row_idx_v, row_i);
816 :
817 : size_t i = 0;
818 0 : for (const auto &idx : interp_op.GetProbes())
819 : {
820 0 : for (int i_dim = 0; i_dim < v_dim; i_dim++)
821 : {
822 0 : auto val = probe_field[i * v_dim + i_dim];
823 0 : probe_B->table[format("B{}_{}_{}_re", i_dim, idx, m_ex_idx)] << val.real();
824 : if (HasComplexGridFunction<solver_t>())
825 : {
826 0 : probe_B->table[format("B{}_{}_{}_im", i_dim, idx, m_ex_idx)] << val.imag();
827 : }
828 : }
829 0 : i++;
830 : }
831 0 : probe_B->WriteFullTableTrunc();
832 : }
833 :
834 : template <ProblemType solver_t>
835 : template <ProblemType U>
836 4 : auto PostOperatorCSV<solver_t>::InitializeSurfaceI(const SurfaceCurrentOperator &surf_j_op)
837 : -> std::enable_if_t<U == ProblemType::DRIVEN || U == ProblemType::TRANSIENT, void>
838 : {
839 4 : if (!(surf_j_op.Size() > 0))
840 : {
841 4 : return;
842 : }
843 : using fmt::format;
844 0 : surface_I = TableWithCSVFile(post_dir / "surface-I.csv", reload_table);
845 :
846 0 : Table t; // Define table locally first due to potential reload.
847 0 : auto nr_expected_measurement_cols = 1 + ex_idx_v_all.size() * surf_j_op.Size();
848 0 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
849 0 : t.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
850 0 : for (const auto ex_idx : ex_idx_v_all)
851 : {
852 0 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
853 0 : for (const auto &[idx, data] : surf_j_op)
854 : {
855 0 : t.insert(format("I_{}_{}", idx, ex_idx), format("I_inc[{}]{} (A)", idx, ex_label),
856 : ex_idx);
857 : }
858 : }
859 0 : MoveTableValidateReload(*surface_I, std::move(t));
860 0 : }
861 :
862 : template <ProblemType solver_t>
863 : template <ProblemType U>
864 4 : auto PostOperatorCSV<solver_t>::PrintSurfaceI(const SurfaceCurrentOperator &surf_j_op,
865 : const Units &units)
866 : -> std::enable_if_t<U == ProblemType::DRIVEN || U == ProblemType::TRANSIENT, void>
867 : {
868 4 : if (!surface_I)
869 : {
870 : return;
871 : }
872 : using fmt::format;
873 0 : CheckAppendIndex(surface_I->table["idx"], row_idx_v, row_i);
874 0 : for (const auto &[idx, data] : surf_j_op)
875 : {
876 0 : auto I_inc_raw = data.GetExcitationCurrent() * measurement_cache.Jcoeff_excitation;
877 : auto I_inc = units.Dimensionalize<Units::ValueType::CURRENT>(I_inc_raw);
878 0 : surface_I->table[format("I_{}_{}", idx, m_ex_idx)] << I_inc;
879 : }
880 0 : surface_I->WriteFullTableTrunc();
881 : }
882 :
883 : template <ProblemType solver_t>
884 : template <ProblemType U>
885 22 : auto PostOperatorCSV<solver_t>::InitializePortVI(const SpaceOperator &fem_op)
886 : -> std::enable_if_t<U == ProblemType::EIGENMODE || U == ProblemType::DRIVEN ||
887 : U == ProblemType::TRANSIENT,
888 : void>
889 : {
890 22 : if (!(fem_op.GetLumpedPortOp().Size() > 0))
891 : {
892 2 : return;
893 : }
894 : using fmt::format;
895 : // Currently only works for lumped ports.
896 : const auto &lumped_port_op = fem_op.GetLumpedPortOp();
897 60 : port_V = TableWithCSVFile(post_dir / "port-V.csv", reload_table);
898 60 : port_I = TableWithCSVFile(post_dir / "port-I.csv", reload_table);
899 :
900 20 : Table tV; // Define table locally first due to potential reload.
901 20 : Table tI;
902 :
903 20 : auto nr_expected_measurement_cols = 1 + ex_idx_v_all.size() * lumped_port_op.Size();
904 20 : tV.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
905 20 : tI.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
906 :
907 20 : tV.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
908 40 : tI.insert("idx", LabelIndexCol(solver_t), -1, 0, PrecIndexCol(solver_t), "");
909 46 : for (const auto ex_idx : ex_idx_v_all)
910 : {
911 26 : std::string ex_label = HasSingleExIdx() ? "" : format("[{}]", ex_idx);
912 :
913 : // Print incident signal, if solver supports excitation on ports.
914 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::TRANSIENT)
915 : {
916 26 : auto ex_spec = fem_op.GetPortExcitations().excitations.at(ex_idx);
917 52 : for (const auto &idx : ex_spec.lumped_port)
918 : {
919 52 : tV.insert(format("inc{}_{}", idx, ex_idx), format("V_inc[{}]{} (V)", idx, ex_label),
920 : ex_idx);
921 52 : tI.insert(format("inc{}_{}", idx, ex_idx), format("I_inc[{}]{} (A)", idx, ex_label),
922 : ex_idx);
923 : }
924 26 : }
925 96 : for (const auto &[idx, data] : lumped_port_op)
926 : {
927 : if constexpr (HasComplexGridFunction<solver_t>())
928 : {
929 136 : tV.insert(format("re{}_{}", idx, ex_idx),
930 68 : format("Re{{V[{}]{}}} (V)", idx, ex_label), ex_idx);
931 136 : tV.insert(format("im{}_{}", idx, ex_idx),
932 68 : format("Im{{V[{}]{}}} (V)", idx, ex_label), ex_idx);
933 136 : tI.insert(format("re{}_{}", idx, ex_idx),
934 68 : format("Re{{I[{}]{}}} (A)", idx, ex_label), ex_idx);
935 136 : tI.insert(format("im{}_{}", idx, ex_idx),
936 68 : format("Im{{I[{}]{}}} (A)", idx, ex_label), ex_idx);
937 : }
938 : else
939 : {
940 4 : tV.insert(format("re{}_{}", idx, ex_idx), format("V[{}]{} (V)", idx, ex_label),
941 : ex_idx);
942 4 : tI.insert(format("re{}_{}", idx, ex_idx), format("I[{}]{} (A)", idx, ex_label),
943 : ex_idx);
944 : }
945 : }
946 : }
947 20 : MoveTableValidateReload(*port_V, std::move(tV));
948 12 : MoveTableValidateReload(*port_I, std::move(tI));
949 28 : }
950 :
951 : template <ProblemType solver_t>
952 : template <ProblemType U>
953 6 : auto PostOperatorCSV<solver_t>::PrintPortVI(const LumpedPortOperator &lumped_port_op,
954 : const Units &units)
955 : -> std::enable_if_t<U == ProblemType::EIGENMODE || U == ProblemType::DRIVEN ||
956 : U == ProblemType::TRANSIENT,
957 : void>
958 : {
959 6 : if (!port_V) // no need to recheck port_I
960 : {
961 : return;
962 : }
963 : using fmt::format;
964 : // Currently only works for lumped ports.
965 : // Postprocess the frequency domain lumped port voltages and currents (complex magnitude
966 : // = sqrt(2) * RMS).
967 :
968 4 : CheckAppendIndex(port_V->table["idx"], row_idx_v, row_i);
969 4 : CheckAppendIndex(port_I->table["idx"], row_idx_v, row_i);
970 :
971 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::TRANSIENT)
972 : {
973 8 : for (const auto &[idx, data] : lumped_port_op)
974 : {
975 4 : if (data.excitation == m_ex_idx)
976 : {
977 4 : auto Jcoeff = measurement_cache.Jcoeff_excitation;
978 4 : double V_inc = data.GetExcitationVoltage() * Jcoeff;
979 : double I_inc = (std::abs(V_inc) > 0.0)
980 4 : ? data.GetExcitationPower() * Jcoeff * Jcoeff / V_inc
981 : : 0.0;
982 :
983 4 : port_V->table[format("inc{}_{}", idx, m_ex_idx)]
984 8 : << units.Dimensionalize<Units::ValueType::VOLTAGE>(V_inc);
985 4 : port_I->table[format("inc{}_{}", idx, m_ex_idx)]
986 8 : << units.Dimensionalize<Units::ValueType::CURRENT>(I_inc);
987 : }
988 : }
989 : }
990 :
991 8 : for (const auto &[idx, data] : measurement_cache.lumped_port_vi)
992 : {
993 8 : port_V->table[fmt::format("re{}_{}", idx, m_ex_idx)] << data.V.real();
994 8 : port_I->table[fmt::format("re{}_{}", idx, m_ex_idx)] << data.I.real();
995 :
996 : if constexpr (HasComplexGridFunction<solver_t>())
997 : {
998 4 : port_V->table[fmt::format("im{}_{}", idx, m_ex_idx)] << data.V.imag();
999 4 : port_I->table[fmt::format("im{}_{}", idx, m_ex_idx)] << data.I.imag();
1000 : }
1001 : }
1002 4 : port_V->WriteFullTableTrunc();
1003 4 : port_I->WriteFullTableTrunc();
1004 : }
1005 :
1006 : template <ProblemType solver_t>
1007 : template <ProblemType U>
1008 2 : auto PostOperatorCSV<solver_t>::InitializePortS(const SpaceOperator &fem_op)
1009 : -> std::enable_if_t<U == ProblemType::DRIVEN, void>
1010 : {
1011 2 : if (!fem_op.GetPortExcitations().IsMultipleSimple() ||
1012 2 : !((fem_op.GetLumpedPortOp().Size() > 0) xor (fem_op.GetWavePortOp().Size() > 0)))
1013 : {
1014 0 : return;
1015 : }
1016 : using fmt::format;
1017 6 : port_S = TableWithCSVFile(post_dir / "port-S.csv", reload_table);
1018 :
1019 2 : Table t; // Define table locally first due to potential reload.
1020 :
1021 2 : auto nr_ports = fem_op.GetLumpedPortOp().Size() + fem_op.GetWavePortOp().Size();
1022 :
1023 2 : auto nr_expected_measurement_cols = 1 + ex_idx_v_all.size() * nr_ports;
1024 2 : t.reserve(nr_expected_measurement_rows, nr_expected_measurement_cols);
1025 2 : t.insert("idx", "f (GHz)", -1, 0, PrecIndexCol(solver_t), "");
1026 :
1027 4 : for (const auto ex_idx : ex_idx_v_all)
1028 : {
1029 : // TODO(C++20): Combine identical loops with ranges + projection.
1030 4 : for (const auto &[o_idx, data] : fem_op.GetLumpedPortOp())
1031 : {
1032 4 : t.insert(format("abs_{}_{}", o_idx, ex_idx),
1033 2 : format("|S[{}][{}]| (dB)", o_idx, ex_idx), ex_idx);
1034 4 : t.insert(format("arg_{}_{}", o_idx, ex_idx),
1035 2 : format("arg(S[{}][{}]) (deg.)", o_idx, ex_idx), ex_idx);
1036 : }
1037 2 : for (const auto &[o_idx, data] : fem_op.GetWavePortOp())
1038 : {
1039 0 : t.insert(format("abs_{}_{}", o_idx, ex_idx),
1040 0 : format("|S[{}][{}]| (dB)", o_idx, ex_idx), ex_idx);
1041 0 : t.insert(format("arg_{}_{}", o_idx, ex_idx),
1042 0 : format("arg(S[{}][{}]) (deg.)", o_idx, ex_idx), ex_idx);
1043 : }
1044 : }
1045 2 : MoveTableValidateReload(*port_S, std::move(t));
1046 2 : }
1047 :
1048 : template <ProblemType solver_t>
1049 : template <ProblemType U>
1050 2 : auto PostOperatorCSV<solver_t>::PrintPortS()
1051 : -> std::enable_if_t<U == ProblemType::DRIVEN, void>
1052 : {
1053 2 : if (!port_S)
1054 : {
1055 : return;
1056 : }
1057 : using fmt::format;
1058 2 : CheckAppendIndex(port_S->table["idx"], row_idx_v, row_i);
1059 4 : for (const auto &[idx, data] : measurement_cache.lumped_port_vi)
1060 : {
1061 4 : port_S->table[format("abs_{}_{}", idx, m_ex_idx)] << Measurement::Magnitude(data.S);
1062 4 : port_S->table[format("arg_{}_{}", idx, m_ex_idx)] << Measurement::Phase(data.S);
1063 : }
1064 2 : for (const auto &[idx, data] : measurement_cache.wave_port_vi)
1065 : {
1066 0 : port_S->table[format("abs_{}_{}", idx, m_ex_idx)] << Measurement::Magnitude(data.S);
1067 0 : port_S->table[format("arg_{}_{}", idx, m_ex_idx)] << Measurement::Phase(data.S);
1068 : }
1069 2 : port_S->WriteFullTableTrunc();
1070 : }
1071 :
1072 : template <ProblemType solver_t>
1073 : template <ProblemType U>
1074 2 : auto PostOperatorCSV<solver_t>::InitializeEig()
1075 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1076 : {
1077 : using fmt::format;
1078 6 : eig = TableWithCSVFile(post_dir / "eig.csv");
1079 2 : eig->table.reserve(nr_expected_measurement_rows, 6);
1080 2 : eig->table.insert("idx", "m", -1, 0, PrecIndexCol(solver_t), "");
1081 2 : eig->table.insert("f_re", "Re{f} (GHz)");
1082 2 : eig->table.insert("f_im", "Im{f} (GHz)");
1083 2 : eig->table.insert("q", "Q");
1084 2 : eig->table.insert("err_back", "Error (Bkwd.)");
1085 2 : eig->table.insert("err_abs", "Error (Abs.)");
1086 2 : eig->WriteFullTableTrunc();
1087 2 : }
1088 :
1089 : template <ProblemType solver_t>
1090 : template <ProblemType U>
1091 2 : auto PostOperatorCSV<solver_t>::PrintEig()
1092 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1093 : {
1094 2 : if (!eig) // trivial check
1095 : {
1096 : return;
1097 : }
1098 2 : eig->table["idx"] << row_idx_v;
1099 4 : eig->table["f_re"] << measurement_cache.freq.real();
1100 4 : eig->table["f_im"] << measurement_cache.freq.imag();
1101 2 : eig->table["q"] << measurement_cache.eigenmode_Q;
1102 2 : eig->table["err_back"] << measurement_cache.error_bkwd;
1103 2 : eig->table["err_abs"] << measurement_cache.error_abs;
1104 2 : eig->WriteFullTableTrunc();
1105 : }
1106 :
1107 : template <ProblemType solver_t>
1108 : template <ProblemType U>
1109 2 : auto PostOperatorCSV<solver_t>::InitializeEigPortEPR(
1110 : const LumpedPortOperator &lumped_port_op)
1111 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1112 : {
1113 : // TODO(C++20): Make this a filtered iterator in LumpedPortOp.
1114 2 : for (const auto &[idx, data] : lumped_port_op)
1115 : {
1116 0 : if (std::abs(data.L) > 0.0)
1117 : {
1118 0 : ports_with_L.push_back(idx);
1119 : }
1120 : }
1121 2 : if (ports_with_L.empty())
1122 : {
1123 : return;
1124 : }
1125 : using fmt::format;
1126 0 : port_EPR = TableWithCSVFile(post_dir / "port-EPR.csv");
1127 0 : port_EPR->table.reserve(nr_expected_measurement_rows, 1 + ports_with_L.size());
1128 0 : port_EPR->table.insert("idx", "m", -1, 0, PrecIndexCol(solver_t), "");
1129 0 : for (const auto idx : ports_with_L)
1130 : {
1131 0 : port_EPR->table.insert(format("p_{}", idx), format("p[{}]", idx));
1132 : }
1133 0 : port_EPR->WriteFullTableTrunc();
1134 : }
1135 :
1136 : template <ProblemType solver_t>
1137 : template <ProblemType U>
1138 2 : auto PostOperatorCSV<solver_t>::PrintEigPortEPR()
1139 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1140 : {
1141 2 : if (!port_EPR)
1142 : {
1143 : return;
1144 : }
1145 : using fmt::format;
1146 0 : port_EPR->table["idx"] << row_idx_v;
1147 0 : for (const auto idx : ports_with_L)
1148 : {
1149 0 : auto vi = measurement_cache.lumped_port_vi.at(idx);
1150 0 : port_EPR->table[format("p_{}", idx)] << vi.inductive_energy_participation;
1151 : }
1152 0 : port_EPR->WriteFullTableTrunc();
1153 : }
1154 :
1155 : template <ProblemType solver_t>
1156 : template <ProblemType U>
1157 2 : auto PostOperatorCSV<solver_t>::InitializeEigPortQ(const LumpedPortOperator &lumped_port_op)
1158 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1159 : {
1160 : // TODO(C++20): Make this a filtered iterator in LumpedPortOp.
1161 2 : for (const auto &[idx, data] : lumped_port_op)
1162 : {
1163 0 : if (std::abs(data.R) > 0.0)
1164 : {
1165 0 : ports_with_R.push_back(idx);
1166 : }
1167 : }
1168 2 : if (ports_with_R.empty())
1169 : {
1170 : return;
1171 : }
1172 : using fmt::format;
1173 0 : port_Q = TableWithCSVFile(post_dir / "port-Q.csv");
1174 0 : port_Q->table.reserve(nr_expected_measurement_rows, 1 + ports_with_R.size());
1175 0 : port_Q->table.insert("idx", "m", -1, 0, PrecIndexCol(solver_t), "");
1176 0 : for (const auto idx : ports_with_R)
1177 : {
1178 0 : port_Q->table.insert(format("Ql_{}", idx), format("Q_ext[{}]", idx));
1179 0 : port_Q->table.insert(format("Kl_{}", idx), format("κ_ext[{}] (GHz)", idx));
1180 : }
1181 0 : port_Q->WriteFullTableTrunc();
1182 : }
1183 :
1184 : template <ProblemType solver_t>
1185 : template <ProblemType U>
1186 2 : auto PostOperatorCSV<solver_t>::PrintEigPortQ()
1187 : -> std::enable_if_t<U == ProblemType::EIGENMODE, void>
1188 : {
1189 2 : if (!port_Q)
1190 : {
1191 : return;
1192 : }
1193 : using fmt::format;
1194 0 : port_Q->table["idx"] << row_idx_v;
1195 0 : for (const auto idx : ports_with_R)
1196 : {
1197 0 : auto vi = measurement_cache.lumped_port_vi.at(idx);
1198 0 : port_Q->table[format("Ql_{}", idx)] << vi.quality_factor;
1199 0 : port_Q->table[format("Kl_{}", idx)] << vi.mode_port_kappa;
1200 : }
1201 0 : port_Q->WriteFullTableTrunc();
1202 : }
1203 :
1204 : template <ProblemType solver_t>
1205 0 : void PostOperatorCSV<solver_t>::PrintErrorIndicator(
1206 : bool is_root, const ErrorIndicator::SummaryStatistics &indicator_stats)
1207 : {
1208 0 : if (!is_root)
1209 : {
1210 0 : return;
1211 : }
1212 :
1213 0 : TableWithCSVFile error_indicator(post_dir / "error-indicators.csv");
1214 0 : error_indicator.table.reserve(1, 4);
1215 :
1216 0 : error_indicator.table.insert(Column("norm", "Norm") << indicator_stats.norm);
1217 0 : error_indicator.table.insert(Column("min", "Minimum") << indicator_stats.min);
1218 0 : error_indicator.table.insert(Column("max", "Maximum") << indicator_stats.max);
1219 0 : error_indicator.table.insert(Column("mean", "Mean") << indicator_stats.mean);
1220 :
1221 0 : error_indicator.WriteFullTableTrunc();
1222 : }
1223 :
1224 : template <ProblemType solver_t>
1225 15 : void PostOperatorCSV<solver_t>::InitializeCSVDataCollection(
1226 : const PostOperator<solver_t> &post_op)
1227 : {
1228 15 : if (!Mpi::Root(post_op.fem_op->GetComm()))
1229 : {
1230 : return;
1231 : }
1232 10 : InitializeDomainE(post_op.dom_post_op);
1233 10 : InitializeSurfaceF(post_op.surf_post_op);
1234 10 : InitializeSurfaceQ(post_op.surf_post_op);
1235 :
1236 : #if defined(MFEM_USE_GSLIB)
1237 8 : InitializeProbeE(post_op.interp_op);
1238 8 : InitializeProbeB(post_op.interp_op);
1239 : #endif
1240 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::TRANSIENT)
1241 : {
1242 4 : InitializeSurfaceI(post_op.fem_op->GetSurfaceCurrentOp());
1243 : }
1244 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::EIGENMODE ||
1245 : solver_t == ProblemType::TRANSIENT)
1246 : {
1247 6 : InitializePortVI(*post_op.fem_op);
1248 : }
1249 : if constexpr (solver_t == ProblemType::DRIVEN)
1250 : {
1251 2 : InitializePortS(*post_op.fem_op);
1252 : }
1253 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::EIGENMODE)
1254 : {
1255 4 : InitializeFarFieldE(post_op.surf_post_op);
1256 : }
1257 : if constexpr (solver_t == ProblemType::EIGENMODE)
1258 : {
1259 2 : InitializeEig();
1260 2 : InitializeEigPortEPR(post_op.fem_op->GetLumpedPortOp());
1261 2 : InitializeEigPortQ(post_op.fem_op->GetLumpedPortOp());
1262 : }
1263 : }
1264 :
1265 : template <ProblemType solver_t>
1266 15 : void PostOperatorCSV<solver_t>::PrintAllCSVData(
1267 : const PostOperator<solver_t> &post_op, const Measurement &non_dim_measurement_cache,
1268 : double idx_value_dimensionful, int step)
1269 : {
1270 15 : if (!Mpi::Root(post_op.fem_op->GetComm()))
1271 : {
1272 : return;
1273 : }
1274 10 : row_idx_v = idx_value_dimensionful;
1275 10 : row_i = step;
1276 :
1277 : // PostOperator acts on a nondimensional measurement cache, we write a dimensional cache.
1278 10 : measurement_cache = Measurement::Dimensionalize(post_op.units, non_dim_measurement_cache);
1279 :
1280 10 : PrintDomainE();
1281 10 : PrintSurfaceF();
1282 10 : PrintSurfaceQ();
1283 :
1284 : #if defined(MFEM_USE_GSLIB)
1285 10 : PrintProbeE(post_op.interp_op);
1286 10 : PrintProbeB(post_op.interp_op);
1287 : #endif
1288 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::TRANSIENT)
1289 : {
1290 4 : PrintSurfaceI(post_op.fem_op->GetSurfaceCurrentOp(), post_op.units);
1291 : }
1292 :
1293 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::EIGENMODE ||
1294 : solver_t == ProblemType::TRANSIENT)
1295 : {
1296 6 : PrintPortVI(post_op.fem_op->GetLumpedPortOp(), post_op.units);
1297 : }
1298 : if constexpr (solver_t == ProblemType::DRIVEN)
1299 : {
1300 2 : PrintPortS();
1301 : }
1302 : if constexpr (solver_t == ProblemType::EIGENMODE || solver_t == ProblemType::DRIVEN)
1303 : {
1304 4 : PrintFarFieldE(post_op.surf_post_op);
1305 : }
1306 : if constexpr (solver_t == ProblemType::EIGENMODE)
1307 : {
1308 2 : PrintEig();
1309 2 : PrintEigPortEPR();
1310 2 : PrintEigPortQ();
1311 : }
1312 : }
1313 :
1314 : template <ProblemType solver_t>
1315 31 : PostOperatorCSV<solver_t>::PostOperatorCSV(const IoData &iodata,
1316 31 : const fem_op_t<solver_t> &fem_op)
1317 : {
1318 31 : if (!Mpi::Root(fem_op.GetComm()))
1319 : {
1320 : return;
1321 : }
1322 :
1323 26 : post_dir = iodata.problem.output;
1324 :
1325 : // Initialize multi-excitation column group index. Only driven or transient support
1326 : // excitations; for other solvers this is default to a single idx=0.
1327 : if constexpr (solver_t == ProblemType::DRIVEN || solver_t == ProblemType::TRANSIENT)
1328 : {
1329 : auto excitation_helper = fem_op.GetPortExcitations();
1330 : ex_idx_v_all.clear();
1331 20 : ex_idx_v_all.reserve(excitation_helper.Size());
1332 20 : std::transform(excitation_helper.begin(), excitation_helper.end(),
1333 : std::back_inserter(ex_idx_v_all),
1334 26 : [](const auto &pair) { return pair.first; });
1335 : // Default to the first excitation.
1336 20 : ex_idx_i = 0;
1337 20 : m_ex_idx = ex_idx_v_all.front();
1338 : }
1339 :
1340 : // Driven solver: can have non-trivial restart.
1341 : if constexpr (solver_t == ProblemType::DRIVEN)
1342 : {
1343 18 : nr_expected_measurement_rows = iodata.solver.driven.sample_f.size();
1344 18 : reload_table = (iodata.solver.driven.restart != 1);
1345 :
1346 18 : row_i = std::size_t(iodata.solver.driven.restart - 1) % nr_expected_measurement_rows;
1347 18 : ex_idx_i = std::size_t(iodata.solver.driven.restart - 1) / nr_expected_measurement_rows;
1348 18 : m_ex_idx = ex_idx_v_all.at(ex_idx_i);
1349 : }
1350 :
1351 : // Non-driven solver: get nr_expected_measurement_rows to reserve table space.
1352 : if (solver_t == ProblemType::EIGENMODE)
1353 : {
1354 2 : nr_expected_measurement_rows = iodata.solver.eigenmode.n;
1355 : }
1356 : else if (solver_t == ProblemType::ELECTROSTATIC)
1357 : {
1358 2 : nr_expected_measurement_rows = iodata.solver.electrostatic.n_post;
1359 : }
1360 : else if (solver_t == ProblemType::MAGNETOSTATIC)
1361 : {
1362 2 : nr_expected_measurement_rows = iodata.solver.magnetostatic.n_post;
1363 : }
1364 : else if (solver_t == ProblemType::TRANSIENT)
1365 : {
1366 : // Estimate number for fixed (linear) stepping.
1367 2 : nr_expected_measurement_rows =
1368 2 : std::size_t(iodata.solver.transient.max_t / iodata.solver.transient.delta_t) + 1;
1369 : }
1370 0 : }
1371 :
1372 : // Explicit template instantiation.
1373 : template class PostOperatorCSV<ProblemType::DRIVEN>;
1374 : template class PostOperatorCSV<ProblemType::EIGENMODE>;
1375 : template class PostOperatorCSV<ProblemType::ELECTROSTATIC>;
1376 : template class PostOperatorCSV<ProblemType::MAGNETOSTATIC>;
1377 : template class PostOperatorCSV<ProblemType::TRANSIENT>;
1378 :
1379 : // Function explict needed testing since everywhere it's through PostOperator.
1380 : // TODO(C++20): with requires, we won't need a second template.
1381 :
1382 : template auto PostOperatorCSV<ProblemType::DRIVEN>::InitializePortVI<ProblemType::DRIVEN>(
1383 : const SpaceOperator &fem_op) -> void;
1384 :
1385 : } // namespace palace
|
