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 "configfile.hpp"
5 :
6 : #include <algorithm>
7 : #include <iterator>
8 : #include <string_view>
9 : #include <fmt/format.h>
10 : #include <fmt/ranges.h>
11 : #include <mfem.hpp>
12 : #include <nlohmann/json.hpp>
13 :
14 : // This is similar to NLOHMANN_JSON_SERIALIZE_ENUM, but results in an error if an enum
15 : // value corresponding to the string cannot be found. Also adds an overload for stream
16 : // printing enum values.
17 : #define PALACE_JSON_SERIALIZE_ENUM(ENUM_TYPE, ...) \
18 : template <typename BasicJsonType> \
19 : inline void to_json(BasicJsonType &j, const ENUM_TYPE &e) \
20 : { \
21 : static_assert(std::is_enum<ENUM_TYPE>::value, #ENUM_TYPE " must be an enum!"); \
22 : static const std::pair<ENUM_TYPE, BasicJsonType> m[] = __VA_ARGS__; \
23 : auto it = std::find_if(std::begin(m), std::end(m), \
24 : [e](const std::pair<ENUM_TYPE, BasicJsonType> &ej_pair) \
25 : { return ej_pair.first == e; }); \
26 : MFEM_VERIFY(it != std::end(m), \
27 : "Invalid value for " << #ENUM_TYPE " given when parsing to JSON!"); \
28 : j = it->second; \
29 : } \
30 : template <typename BasicJsonType> \
31 : inline void from_json(const BasicJsonType &j, ENUM_TYPE &e) \
32 : { \
33 : static_assert(std::is_enum<ENUM_TYPE>::value, #ENUM_TYPE " must be an enum!"); \
34 : static const std::pair<ENUM_TYPE, BasicJsonType> m[] = __VA_ARGS__; \
35 : auto it = std::find_if(std::begin(m), std::end(m), \
36 : [j](const std::pair<ENUM_TYPE, BasicJsonType> &ej_pair) \
37 : { return ej_pair.second == j; }); \
38 : MFEM_VERIFY(it != std::end(m), \
39 : "Invalid value (" << j << ") for " \
40 : << #ENUM_TYPE \
41 : " given in the configuration file when parsing from JSON!"); \
42 : e = it->first; \
43 : } \
44 : std::ostream &operator<<(std::ostream &os, const ENUM_TYPE &e) \
45 : { \
46 : static const std::pair<ENUM_TYPE, const char *> m[] = __VA_ARGS__; \
47 : os << std::find_if(std::begin(m), std::end(m), \
48 : [e](const std::pair<ENUM_TYPE, const char *> &ej_pair) \
49 : { return ej_pair.first == e; }) \
50 : ->second; \
51 : return os; \
52 : }
53 :
54 : using json = nlohmann::json;
55 : namespace palace
56 : {
57 : // Helpers for converting enums specified in labels.hpp. Must be done in palace scope rather
58 : // than palace::config scope to ensure argument-dependent-lookup succeeds in json.
59 :
60 : // Helper for converting string keys to enum for CoordinateSystem.
61 0 : PALACE_JSON_SERIALIZE_ENUM(CoordinateSystem,
62 : {{CoordinateSystem::CARTESIAN, "Cartesian"},
63 : {CoordinateSystem::CYLINDRICAL, "Cylindrical"}})
64 :
65 : // Helper for converting string keys to enum for ProblemType.
66 66 : PALACE_JSON_SERIALIZE_ENUM(ProblemType, {{ProblemType::DRIVEN, "Driven"},
67 : {ProblemType::EIGENMODE, "Eigenmode"},
68 : {ProblemType::ELECTROSTATIC, "Electrostatic"},
69 : {ProblemType::MAGNETOSTATIC, "Magnetostatic"},
70 : {ProblemType::TRANSIENT, "Transient"}})
71 :
72 : // Helper for converting string keys to enum for EigenSolverBackend.
73 0 : PALACE_JSON_SERIALIZE_ENUM(EigenSolverBackend, {{EigenSolverBackend::DEFAULT, "Default"},
74 : {EigenSolverBackend::SLEPC, "SLEPc"},
75 : {EigenSolverBackend::ARPACK, "ARPACK"}})
76 :
77 : // Helper for converting string keys to enum for EigenSolverBackend.
78 0 : PALACE_JSON_SERIALIZE_ENUM(NonlinearEigenSolver, {{NonlinearEigenSolver::HYBRID, "Hybrid"},
79 : {NonlinearEigenSolver::SLP, "SLP"}})
80 :
81 : // Helper for converting string keys to enum for SurfaceFlux.
82 146 : PALACE_JSON_SERIALIZE_ENUM(SurfaceFlux, {{SurfaceFlux::ELECTRIC, "Electric"},
83 : {SurfaceFlux::MAGNETIC, "Magnetic"},
84 : {SurfaceFlux::POWER, "Power"}})
85 :
86 : // Helper for converting string keys to enum for InterfaceDielectric.
87 242 : PALACE_JSON_SERIALIZE_ENUM(InterfaceDielectric, {{InterfaceDielectric::DEFAULT, "Default"},
88 : {InterfaceDielectric::MA, "MA"},
89 : {InterfaceDielectric::MS, "MS"},
90 : {InterfaceDielectric::SA, "SA"}})
91 :
92 : // Helper for converting string keys to enum for FrequencySampling.
93 125 : PALACE_JSON_SERIALIZE_ENUM(FrequencySampling, {{FrequencySampling::DEFAULT, "Default"},
94 : {FrequencySampling::LINEAR, "Linear"},
95 : {FrequencySampling::LOG, "Log"},
96 : {FrequencySampling::POINT, "Point"}})
97 :
98 : // Helper for converting string keys to enum for TimeSteppingScheme and Excitation.
99 0 : PALACE_JSON_SERIALIZE_ENUM(TimeSteppingScheme,
100 : {{TimeSteppingScheme::DEFAULT, "Default"},
101 : {TimeSteppingScheme::GEN_ALPHA, "GeneralizedAlpha"},
102 : {TimeSteppingScheme::RUNGE_KUTTA, "RungeKutta"},
103 : {TimeSteppingScheme::CVODE, "CVODE"},
104 : {TimeSteppingScheme::ARKODE, "ARKODE"}})
105 0 : PALACE_JSON_SERIALIZE_ENUM(Excitation,
106 : {{Excitation::SINUSOIDAL, "Sinusoidal"},
107 : {Excitation::GAUSSIAN, "Gaussian"},
108 : {Excitation::DIFF_GAUSSIAN, "DifferentiatedGaussian"},
109 : {Excitation::MOD_GAUSSIAN, "ModulatedGaussian"},
110 : {Excitation::RAMP_STEP, "Ramp"},
111 : {Excitation::SMOOTH_STEP, "SmoothStep"}})
112 :
113 : // Helper for converting string keys to enum for LinearSolver, KrylovSolver, and
114 : // MultigridCoarsening
115 66 : PALACE_JSON_SERIALIZE_ENUM(LinearSolver, {{LinearSolver::DEFAULT, "Default"},
116 : {LinearSolver::AMS, "AMS"},
117 : {LinearSolver::BOOMER_AMG, "BoomerAMG"},
118 : {LinearSolver::MUMPS, "MUMPS"},
119 : {LinearSolver::SUPERLU, "SuperLU"},
120 : {LinearSolver::STRUMPACK, "STRUMPACK"},
121 : {LinearSolver::STRUMPACK_MP, "STRUMPACK-MP"},
122 : {LinearSolver::JACOBI, "Jacobi"}})
123 90 : PALACE_JSON_SERIALIZE_ENUM(KrylovSolver, {{KrylovSolver::DEFAULT, "Default"},
124 : {KrylovSolver::CG, "CG"},
125 : {KrylovSolver::MINRES, "MINRES"},
126 : {KrylovSolver::GMRES, "GMRES"},
127 : {KrylovSolver::FGMRES, "FGMRES"},
128 : {KrylovSolver::BICGSTAB, "BiCGSTAB"}})
129 0 : PALACE_JSON_SERIALIZE_ENUM(MultigridCoarsening,
130 : {{MultigridCoarsening::LINEAR, "Linear"},
131 : {MultigridCoarsening::LOGARITHMIC, "Logarithmic"}})
132 :
133 : // Helpers for converting string keys to enum for PreconditionerSide, SymbolicFactorization,
134 : // SparseCompression, and Orthogonalization.
135 0 : PALACE_JSON_SERIALIZE_ENUM(PreconditionerSide, {{PreconditionerSide::DEFAULT, "Default"},
136 : {PreconditionerSide::RIGHT, "Right"},
137 : {PreconditionerSide::LEFT, "Left"}})
138 0 : PALACE_JSON_SERIALIZE_ENUM(SymbolicFactorization,
139 : {{SymbolicFactorization::DEFAULT, "Default"},
140 : {SymbolicFactorization::METIS, "METIS"},
141 : {SymbolicFactorization::PARMETIS, "ParMETIS"},
142 : {SymbolicFactorization::SCOTCH, "Scotch"},
143 : {SymbolicFactorization::PTSCOTCH, "PTScotch"},
144 : {SymbolicFactorization::PORD, "PORD"},
145 : {SymbolicFactorization::AMD, "AMD"},
146 : {SymbolicFactorization::RCM, "RCM"}})
147 0 : PALACE_JSON_SERIALIZE_ENUM(SparseCompression,
148 : {{SparseCompression::NONE, "None"},
149 : {SparseCompression::BLR, "BLR"},
150 : {SparseCompression::HSS, "HSS"},
151 : {SparseCompression::HODLR, "HODLR"},
152 : {SparseCompression::ZFP, "ZFP"},
153 : {SparseCompression::BLR_HODLR, "BLR-HODLR"},
154 : {SparseCompression::ZFP_BLR_HODLR, "ZFP-BLR-HODLR"}})
155 0 : PALACE_JSON_SERIALIZE_ENUM(Orthogonalization, {{Orthogonalization::MGS, "MGS"},
156 : {Orthogonalization::CGS, "CGS"},
157 : {Orthogonalization::CGS2, "CGS2"}})
158 :
159 : // Helpers for converting string keys to enum for Device.
160 66 : PALACE_JSON_SERIALIZE_ENUM(Device, {{Device::CPU, "CPU"},
161 : {Device::GPU, "GPU"},
162 : {Device::DEBUG, "Debug"}})
163 : } // namespace palace
164 :
165 : namespace palace::config
166 : {
167 :
168 : namespace
169 : {
170 :
171 149 : int AtIndex(json::iterator &port_it, std::string_view errmsg_parent)
172 : {
173 149 : MFEM_VERIFY(
174 : port_it->find("Index") != port_it->end(),
175 : fmt::format("Missing {} \"Index\" in the configuration file!", errmsg_parent));
176 149 : auto index = port_it->at("Index").get<int>();
177 155 : MFEM_VERIFY(index > 0, fmt::format("The {} \"Index\" should be an integer > 0; got {}",
178 : errmsg_parent, index));
179 147 : return index;
180 : }
181 :
182 : template <std::size_t N>
183 32 : void ParseSymmetricMatrixData(json &mat, const std::string &name,
184 : SymmetricMatrixData<N> &data)
185 : {
186 32 : auto it = mat.find(name);
187 32 : if (it != mat.end() && it->is_array())
188 : {
189 : // Attempt to parse as an array.
190 0 : data.s = it->get<std::array<double, N>>();
191 : }
192 : else
193 : {
194 : // Fall back to scalar parsing with default.
195 32 : double s = mat.value(name, data.s[0]);
196 : data.s.fill(s);
197 : }
198 32 : data.v = mat.value("MaterialAxes", data.v);
199 32 : }
200 :
201 : // Helper function for extracting element data from the configuration file, either from a
202 : // provided keyword argument of from a specified vector. In extracting the direction various
203 : // checks are performed for validity of the input combinations.
204 63 : void ParseElementData(json &elem, const std::string &name, bool required,
205 : internal::ElementData &data)
206 : {
207 189 : data.attributes = elem.at("Attributes").get<std::vector<int>>(); // Required
208 63 : std::sort(data.attributes.begin(), data.attributes.end());
209 63 : auto it = elem.find(name);
210 63 : if (it != elem.end() && it->is_array())
211 : {
212 : // Attempt to parse as an array.
213 0 : data.direction = it->get<std::array<double, 3>>();
214 0 : data.coordinate_system = elem.value("CoordinateSystem", data.coordinate_system);
215 : }
216 : else
217 : {
218 : // Fall back to parsing as a string (value is optional).
219 63 : MFEM_VERIFY(elem.find("CoordinateSystem") == elem.end(),
220 : "Cannot specify \"CoordinateSystem\" when specifying a direction or side "
221 : "using a string in the configuration file!");
222 : std::string direction;
223 126 : direction = elem.value(name, direction);
224 189 : for (auto &c : direction)
225 : {
226 126 : c = std::tolower(c);
227 : }
228 : const auto xpos = direction.find("x");
229 : const auto ypos = direction.find("y");
230 : const auto zpos = direction.find("z");
231 : const auto rpos = direction.find("r");
232 : const bool xfound = xpos != std::string::npos;
233 : const bool yfound = ypos != std::string::npos;
234 : const bool zfound = zpos != std::string::npos;
235 : const bool rfound = rpos != std::string::npos;
236 63 : if (xfound)
237 : {
238 19 : MFEM_VERIFY(direction.length() == 1 || direction[xpos - 1] == '-' ||
239 : direction[xpos - 1] == '+',
240 : "Missing required sign specification on \"X\" for \""
241 : << name << "\" in the configuration file!");
242 19 : MFEM_VERIFY(!yfound && !zfound && !rfound,
243 : "\"X\" cannot be combined with \"Y\", \"Z\", or \"R\" for \""
244 : << name << "\" in the configuration file!");
245 19 : data.direction[0] =
246 19 : (direction.length() == 1 || direction[xpos - 1] == '+') ? 1.0 : -1.0;
247 19 : data.coordinate_system = CoordinateSystem::CARTESIAN;
248 : }
249 63 : if (yfound)
250 : {
251 44 : MFEM_VERIFY(direction.length() == 1 || direction[ypos - 1] == '-' ||
252 : direction[ypos - 1] == '+',
253 : "Missing required sign specification on \"Y\" for \""
254 : << name << "\" in the configuration file!");
255 44 : MFEM_VERIFY(!xfound && !zfound && !rfound,
256 : "\"Y\" cannot be combined with \"X\", \"Z\", or \"R\" for \""
257 : << name << "\" in the configuration file!");
258 44 : data.direction[1] =
259 44 : direction.length() == 1 || direction[ypos - 1] == '+' ? 1.0 : -1.0;
260 44 : data.coordinate_system = CoordinateSystem::CARTESIAN;
261 : }
262 63 : if (zfound)
263 : {
264 0 : MFEM_VERIFY(direction.length() == 1 || direction[zpos - 1] == '-' ||
265 : direction[zpos - 1] == '+',
266 : "Missing required sign specification on \"Z\" for \""
267 : << name << "\" in the configuration file!");
268 0 : MFEM_VERIFY(!xfound && !yfound && !rfound,
269 : "\"Z\" cannot be combined with \"X\", \"Y\", or \"R\" for \""
270 : << name << "\" in the configuration file!");
271 0 : data.direction[2] =
272 0 : direction.length() == 1 || direction[zpos - 1] == '+' ? 1.0 : -1.0;
273 0 : data.coordinate_system = CoordinateSystem::CARTESIAN;
274 : }
275 63 : if (rfound)
276 : {
277 0 : MFEM_VERIFY(direction.length() == 1 || direction[rpos - 1] == '-' ||
278 : direction[rpos - 1] == '+',
279 : "Missing required sign specification on \"R\" for \""
280 : << name << "\" in the configuration file!");
281 0 : MFEM_VERIFY(!xfound && !yfound && !zfound,
282 : "\"R\" cannot be combined with \"X\", \"Y\", or \"Z\" for \""
283 : << name << "\" in the configuration file!");
284 0 : data.direction[0] =
285 0 : direction.length() == 1 || direction[rpos - 1] == '+' ? 1.0 : -1.0;
286 0 : data.direction[1] = 0.0;
287 0 : data.direction[2] = 0.0;
288 0 : data.coordinate_system = CoordinateSystem::CYLINDRICAL;
289 : }
290 : }
291 63 : MFEM_VERIFY(data.coordinate_system != CoordinateSystem::CYLINDRICAL ||
292 : (data.direction[1] == 0.0 && data.direction[2] == 0.0),
293 : "Parsing azimuthal and longitudinal directions for cylindrical coordinate "
294 : "system directions from the configuration file is not currently supported!");
295 63 : MFEM_VERIFY(
296 : !required || data.direction[0] != 0.0 || data.direction[1] != 0.0 ||
297 : data.direction[2] != 0.0,
298 : "Missing \"" << name
299 : << "\" for an object which requires it in the configuration file!");
300 63 : }
301 :
302 : template <typename T>
303 : std::ostream &operator<<(std::ostream &os, const std::vector<T> &data)
304 : {
305 : os << fmt::format("{}", fmt::join(data, " "));
306 : return os;
307 : }
308 :
309 : template <typename T, std::size_t N>
310 : std::ostream &operator<<(std::ostream &os, const std::array<T, N> &data)
311 : {
312 : os << fmt::format("{}", fmt::join(data, " "));
313 : return os;
314 : }
315 :
316 : template <std::size_t N>
317 : std::ostream &operator<<(std::ostream &os, const SymmetricMatrixData<N> &data)
318 : {
319 : os << "s: " << data.s;
320 : int j = 0;
321 : for (const auto &x : data.v)
322 : {
323 : os << ", v" << j++ << ": " << x;
324 : }
325 : return os;
326 : }
327 :
328 : constexpr bool JSON_DEBUG = false;
329 :
330 : } // namespace
331 :
332 8 : void ProblemData::SetUp(json &config)
333 : {
334 8 : auto problem = config.find("Problem");
335 8 : MFEM_VERIFY(problem != config.end(),
336 : "\"Problem\" must be specified in the configuration file!");
337 8 : MFEM_VERIFY(problem->find("Type") != problem->end(),
338 : "Missing config[\"Problem\"][\"Type\"] in the configuration file!");
339 8 : type = problem->at("Type"); // Required
340 8 : verbose = problem->value("Verbose", verbose);
341 8 : output = problem->value("Output", output);
342 :
343 : // Parse output formats.
344 8 : auto output_formats_it = problem->find("OutputFormats");
345 8 : if (output_formats_it != problem->end())
346 : {
347 0 : output_formats.paraview = output_formats_it->value("Paraview", output_formats.paraview);
348 0 : output_formats.gridfunction =
349 0 : output_formats_it->value("GridFunction", output_formats.gridfunction);
350 : }
351 :
352 : // Check for provided solver configuration data (not required for electrostatics or
353 : // magnetostatics since defaults can be used for every option).
354 8 : auto solver = config.find("Solver");
355 8 : if (type == ProblemType::DRIVEN)
356 : {
357 8 : MFEM_VERIFY(solver->find("Driven") != solver->end(),
358 : "config[\"Problem\"][\"Type\"] == \"Driven\" should be accompanied by a "
359 : "config[\"Solver\"][\"Driven\"] configuration!");
360 : }
361 : else if (type == ProblemType::EIGENMODE)
362 : {
363 0 : MFEM_VERIFY(solver->find("Eigenmode") != solver->end(),
364 : "config[\"Problem\"][\"Type\"] == \"Eigenmode\" should be accompanied by a "
365 : "config[\"Solver\"][\"Eigenmode\"] configuration!");
366 : }
367 : else if (type == ProblemType::ELECTROSTATIC)
368 : {
369 : // MFEM_VERIFY(
370 : // solver->find("Electrostatic") != solver->end(),
371 : // "config[\"Problem\"][\"Type\"] == \"Electrostatic\" should be accompanied by a "
372 : // "config[\"Solver\"][\"Electrostatic\"] configuration!");
373 : }
374 : else if (type == ProblemType::MAGNETOSTATIC)
375 : {
376 : // MFEM_VERIFY(
377 : // solver->find("Magnetostatic") != solver->end(),
378 : // "config[\"Problem\"][\"Type\"] == \"Magnetostatic\" should be accompanied by a "
379 : // "config[\"Solver\"][\"Magnetostatic\"] configuration!");
380 : }
381 : else if (type == ProblemType::TRANSIENT)
382 : {
383 0 : MFEM_VERIFY(solver->find("Transient") != solver->end(),
384 : "config[\"Problem\"][\"Type\"] == \"Transient\" should be accompanied by a "
385 : "config[\"Solver\"][\"Transient\"] configuration!");
386 : }
387 :
388 : // Cleanup
389 8 : problem->erase("Type");
390 8 : problem->erase("Verbose");
391 8 : problem->erase("Output");
392 8 : problem->erase("OutputFormats");
393 16 : MFEM_VERIFY(problem->empty(),
394 : "Found an unsupported configuration file keyword under \"Problem\"!\n"
395 : << problem->dump(2));
396 :
397 : // Debug
398 : if constexpr (JSON_DEBUG)
399 : {
400 : std::cout << "Type: " << type << '\n';
401 : std::cout << "Verbose: " << verbose << '\n';
402 : std::cout << "Output: " << output << '\n';
403 : std::cout << "OutputFormats.Paraview: " << output_formats.paraview << '\n';
404 : std::cout << "OutputFormats.GridFunction: " << output_formats.gridfunction << '\n';
405 : }
406 8 : }
407 :
408 8 : void RefinementData::SetUp(json &model)
409 : {
410 8 : auto refinement = model.find("Refinement");
411 8 : if (refinement == model.end())
412 : {
413 0 : return;
414 : }
415 :
416 : // Options for AMR.
417 8 : tol = refinement->value("Tol", tol);
418 8 : max_it = refinement->value("MaxIts", max_it);
419 8 : max_size = refinement->value("MaxSize", max_size);
420 8 : nonconformal = refinement->value("Nonconformal", nonconformal);
421 8 : max_nc_levels = refinement->value("MaxNCLevels", max_nc_levels);
422 8 : update_fraction = refinement->value("UpdateFraction", update_fraction);
423 8 : maximum_imbalance = refinement->value("MaximumImbalance", maximum_imbalance);
424 8 : save_adapt_iterations = refinement->value("SaveAdaptIterations", save_adapt_iterations);
425 8 : save_adapt_mesh = refinement->value("SaveAdaptMesh", save_adapt_mesh);
426 8 : MFEM_VERIFY(tol > 0.0, "config[\"Refinement\"][\"Tol\"] must be strictly positive!");
427 8 : MFEM_VERIFY(max_it >= 0, "config[\"Refinement\"][\"MaxIts\"] must be non-negative!");
428 8 : MFEM_VERIFY(max_size >= 0, "config[\"Refinement\"][\"MaxSize\"] must be non-negative!");
429 8 : MFEM_VERIFY(max_nc_levels >= 0,
430 : "config[\"Refinement\"][\"MaxNCLevels\"] must be non-negative!");
431 8 : MFEM_VERIFY(update_fraction > 0 && update_fraction < 1,
432 : "config[\"Refinement\"][\"UpdateFraction\" must be in (0,1)!");
433 8 : MFEM_VERIFY(
434 : maximum_imbalance >= 1,
435 : "config[\"Refinement\"][\"MaximumImbalance\"] must be greater than or equal to 1!");
436 :
437 : // Options for a priori refinement.
438 8 : uniform_ref_levels = refinement->value("UniformLevels", uniform_ref_levels);
439 8 : ser_uniform_ref_levels = refinement->value("SerialUniformLevels", ser_uniform_ref_levels);
440 8 : MFEM_VERIFY(uniform_ref_levels >= 0 && ser_uniform_ref_levels >= 0,
441 : "Number of uniform mesh refinement levels must be non-negative!");
442 8 : auto boxes = refinement->find("Boxes");
443 8 : if (boxes != refinement->end())
444 : {
445 0 : MFEM_VERIFY(boxes->is_array(), "config[\"Refinement\"][\"Boxes\"] should specify an "
446 : "array in the configuration file!");
447 0 : for (auto it = boxes->begin(); it != boxes->end(); ++it)
448 : {
449 0 : MFEM_VERIFY(
450 : it->find("Levels") != it->end(),
451 : "Missing \"Boxes\" refinement region \"Levels\" in the configuration file!");
452 0 : auto bbmin = it->find("BoundingBoxMin");
453 0 : auto bbmax = it->find("BoundingBoxMax");
454 0 : MFEM_VERIFY(bbmin != it->end() && bbmax != it->end(),
455 : "Missing \"Boxes\" refinement region \"BoundingBoxMin/Max\" in the "
456 : "configuration file!");
457 0 : MFEM_VERIFY(bbmin->is_array() && bbmin->is_array(),
458 : "config[\"Refinement\"][\"Boxes\"][\"BoundingBoxMin/Max\"] should "
459 : "specify an array in the configuration file!");
460 0 : BoxRefinementData &data = box_list.emplace_back();
461 0 : data.ref_levels = it->at("Levels"); // Required
462 0 : data.bbmin = bbmin->get<std::array<double, 3>>(); // Required
463 0 : data.bbmax = bbmax->get<std::array<double, 3>>(); // Required
464 :
465 : // Cleanup
466 0 : it->erase("Levels");
467 0 : it->erase("BoundingBoxMin");
468 0 : it->erase("BoundingBoxMax");
469 0 : MFEM_VERIFY(it->empty(), "Found an unsupported configuration file keyword under "
470 : "config[\"Refinement\"][\"Boxes\"]!\n"
471 : << it->dump(2));
472 :
473 : // Debug
474 : if constexpr (JSON_DEBUG)
475 : {
476 : std::cout << "Levels: " << data.ref_levels << '\n';
477 : std::cout << "BoundingBoxMin: " << data.bbmin << '\n';
478 : std::cout << "BoundingBoxMax: " << data.bbmax << '\n';
479 : }
480 : }
481 : }
482 8 : auto spheres = refinement->find("Spheres");
483 8 : if (spheres != refinement->end())
484 : {
485 0 : MFEM_VERIFY(spheres->is_array(), "config[\"Refinement\"][\"Spheres\"] should specify "
486 : "an array in the configuration file!");
487 0 : for (auto it = spheres->begin(); it != spheres->end(); ++it)
488 : {
489 0 : MFEM_VERIFY(
490 : it->find("Levels") != it->end(),
491 : "Missing \"Spheres\" refinement region \"Levels\" in the configuration file!");
492 0 : auto ctr = it->find("Center");
493 0 : MFEM_VERIFY(ctr != it->end() && it->find("Radius") != it->end(),
494 : "Missing \"Spheres\" refinement region \"Center\" or \"Radius\" in "
495 : "configuration file!");
496 0 : MFEM_VERIFY(ctr->is_array(),
497 : "config[\"Refinement\"][\"Spheres\"][\"Center\"] should specify "
498 : "an array in the configuration file!");
499 0 : SphereRefinementData &data = sphere_list.emplace_back();
500 0 : data.ref_levels = it->at("Levels"); // Required
501 0 : data.r = it->at("Radius"); // Required
502 0 : data.center = ctr->get<std::array<double, 3>>(); // Required
503 :
504 : // Cleanup
505 0 : it->erase("Levels");
506 0 : it->erase("Radius");
507 0 : it->erase("Center");
508 0 : MFEM_VERIFY(it->empty(), "Found an unsupported configuration file keyword under "
509 : "config[\"Refinement\"][\"Spheres\"]!\n"
510 : << it->dump(2));
511 :
512 : // Debug
513 : if constexpr (JSON_DEBUG)
514 : {
515 : std::cout << "Levels: " << data.ref_levels << '\n';
516 : std::cout << "Radius: " << data.r << '\n';
517 : std::cout << "Center: " << data.center << '\n';
518 : }
519 : }
520 : }
521 :
522 : // Cleanup
523 8 : refinement->erase("Tol");
524 8 : refinement->erase("MaxIts");
525 8 : refinement->erase("MaxSize");
526 8 : refinement->erase("Nonconformal");
527 8 : refinement->erase("MaxNCLevels");
528 8 : refinement->erase("UpdateFraction");
529 8 : refinement->erase("MaximumImbalance");
530 8 : refinement->erase("SaveAdaptIterations");
531 8 : refinement->erase("SaveAdaptMesh");
532 8 : refinement->erase("UniformLevels");
533 8 : refinement->erase("SerialUniformLevels");
534 8 : refinement->erase("Boxes");
535 8 : refinement->erase("Spheres");
536 16 : MFEM_VERIFY(refinement->empty(),
537 : "Found an unsupported configuration file keyword under \"Refinement\"!\n"
538 : << refinement->dump(2));
539 :
540 : // Debug
541 : if constexpr (JSON_DEBUG)
542 : {
543 : std::cout << "Tol: " << tol << '\n';
544 : std::cout << "MaxIts: " << max_it << '\n';
545 : std::cout << "MaxSize: " << max_size << '\n';
546 : std::cout << "Nonconformal: " << nonconformal << '\n';
547 : std::cout << "MaxNCLevels: " << max_nc_levels << '\n';
548 : std::cout << "UpdateFraction: " << update_fraction << '\n';
549 : std::cout << "MaximumImbalance: " << maximum_imbalance << '\n';
550 : std::cout << "SaveAdaptIterations: " << save_adapt_iterations << '\n';
551 : std::cout << "SaveAdaptMesh: " << save_adapt_mesh << '\n';
552 : std::cout << "UniformLevels: " << uniform_ref_levels << '\n';
553 : std::cout << "SerialUniformLevels: " << ser_uniform_ref_levels << '\n';
554 : }
555 : }
556 :
557 8 : void ModelData::SetUp(json &config)
558 : {
559 8 : auto model = config.find("Model");
560 8 : MFEM_VERIFY(model != config.end(),
561 : "\"Model\" must be specified in the configuration file!");
562 8 : MFEM_VERIFY(model->find("Mesh") != model->end(),
563 : "Missing config[\"Model\"][\"Mesh\"] file in the configuration file!");
564 16 : mesh = model->at("Mesh"); // Required
565 8 : L0 = model->value("L0", L0);
566 8 : Lc = model->value("Lc", Lc);
567 8 : remove_curvature = model->value("RemoveCurvature", remove_curvature);
568 8 : make_simplex = model->value("MakeSimplex", make_simplex);
569 8 : make_hex = model->value("MakeHexahedral", make_hex);
570 8 : reorder_elements = model->value("ReorderElements", reorder_elements);
571 8 : clean_unused_elements = model->value("CleanUnusedElements", clean_unused_elements);
572 8 : crack_bdr_elements = model->value("CrackInternalBoundaryElements", crack_bdr_elements);
573 8 : refine_crack_elements = model->value("RefineCrackElements", refine_crack_elements);
574 8 : crack_displ_factor = model->value("CrackDisplacementFactor", crack_displ_factor);
575 8 : add_bdr_elements = model->value("AddInterfaceBoundaryElements", add_bdr_elements);
576 8 : export_prerefined_mesh = model->value("ExportPrerefinedMesh", export_prerefined_mesh);
577 8 : reorient_tet_mesh = model->value("ReorientTetMesh", reorient_tet_mesh);
578 8 : partitioning = model->value("Partitioning", partitioning);
579 8 : refinement.SetUp(*model);
580 :
581 : // Cleanup
582 8 : model->erase("Mesh");
583 8 : model->erase("L0");
584 8 : model->erase("Lc");
585 8 : model->erase("RemoveCurvature");
586 8 : model->erase("MakeSimplex");
587 8 : model->erase("MakeHexahedral");
588 8 : model->erase("ReorderElements");
589 8 : model->erase("CleanUnusedElements");
590 8 : model->erase("CrackInternalBoundaryElements");
591 8 : model->erase("RefineCrackElements");
592 8 : model->erase("CrackDisplacementFactor");
593 8 : model->erase("AddInterfaceBoundaryElements");
594 8 : model->erase("ExportPrerefinedMesh");
595 8 : model->erase("ReorientTetMesh");
596 8 : model->erase("Partitioning");
597 8 : model->erase("Refinement");
598 16 : MFEM_VERIFY(model->empty(),
599 : "Found an unsupported configuration file keyword under \"Model\"!\n"
600 : << model->dump(2));
601 :
602 : // Debug
603 : if constexpr (JSON_DEBUG)
604 : {
605 : std::cout << "Mesh: " << mesh << '\n';
606 : std::cout << "L0: " << L0 << '\n';
607 : std::cout << "Lc: " << Lc << '\n';
608 : std::cout << "RemoveCurvature: " << remove_curvature << '\n';
609 : std::cout << "MakeSimplex: " << make_simplex << '\n';
610 : std::cout << "MakeHexahedral: " << make_hex << '\n';
611 : std::cout << "ReorderElements: " << reorder_elements << '\n';
612 : std::cout << "CleanUnusedElements: " << clean_unused_elements << '\n';
613 : std::cout << "CrackInternalBoundaryElements: " << crack_bdr_elements << '\n';
614 : std::cout << "RefineCrackElements: " << refine_crack_elements << '\n';
615 : std::cout << "CrackDisplacementFactor: " << crack_displ_factor << '\n';
616 : std::cout << "AddInterfaceBoundaryElements: " << add_bdr_elements << '\n';
617 : std::cout << "ExportPrerefinedMesh: " << export_prerefined_mesh << '\n';
618 : std::cout << "ReorientTetMesh: " << reorient_tet_mesh << '\n';
619 : std::cout << "Partitioning: " << partitioning << '\n';
620 : }
621 8 : }
622 :
623 8 : void DomainMaterialData::SetUp(json &domains)
624 : {
625 8 : auto materials = domains.find("Materials");
626 8 : MFEM_VERIFY(materials != domains.end() && materials->is_array(),
627 : "\"Materials\" must be specified as an array in the configuration file!");
628 16 : for (auto it = materials->begin(); it != materials->end(); ++it)
629 : {
630 8 : MFEM_VERIFY(
631 : it->find("Attributes") != it->end(),
632 : "Missing \"Attributes\" list for \"Materials\" domain in the configuration file!");
633 : MaterialData &data = emplace_back();
634 24 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
635 8 : std::sort(data.attributes.begin(), data.attributes.end());
636 8 : ParseSymmetricMatrixData(*it, "Permeability", data.mu_r);
637 8 : ParseSymmetricMatrixData(*it, "Permittivity", data.epsilon_r);
638 8 : ParseSymmetricMatrixData(*it, "LossTan", data.tandelta);
639 8 : ParseSymmetricMatrixData(*it, "Conductivity", data.sigma);
640 8 : data.lambda_L = it->value("LondonDepth", data.lambda_L);
641 :
642 : // Cleanup
643 8 : it->erase("Attributes");
644 8 : it->erase("Permeability");
645 8 : it->erase("Permittivity");
646 8 : it->erase("LossTan");
647 8 : it->erase("Conductivity");
648 8 : it->erase("MaterialAxes");
649 8 : it->erase("LondonDepth");
650 16 : MFEM_VERIFY(it->empty(),
651 : "Found an unsupported configuration file keyword under \"Materials\"!\n"
652 : << it->dump(2));
653 :
654 : // Debug
655 : if constexpr (JSON_DEBUG)
656 : {
657 : std::cout << "Attributes: " << data.attributes << '\n';
658 : std::cout << "Permeability: " << data.mu_r << '\n';
659 : std::cout << "Permittivity: " << data.epsilon_r << '\n';
660 : std::cout << "LossTan: " << data.tandelta << '\n';
661 : std::cout << "Conductivity: " << data.sigma << '\n';
662 : std::cout << "LondonDepth: " << data.lambda_L << '\n';
663 : }
664 : }
665 8 : }
666 :
667 8 : void DomainEnergyPostData::SetUp(json &postpro)
668 : {
669 8 : auto energy = postpro.find("Energy");
670 8 : if (energy == postpro.end())
671 : {
672 0 : return;
673 : }
674 8 : MFEM_VERIFY(energy->is_array(),
675 : "\"Energy\" should specify an array in the configuration file!");
676 16 : for (auto it = energy->begin(); it != energy->end(); ++it)
677 : {
678 8 : auto index = AtIndex(it, "\"Energy\" domain");
679 8 : MFEM_VERIFY(
680 : it->find("Attributes") != it->end(),
681 : "Missing \"Attributes\" list for \"Energy\" domain in the configuration file!");
682 8 : auto ret = mapdata.insert(std::make_pair(index, DomainEnergyData()));
683 8 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"Energy\" domains "
684 : "in the configuration file!");
685 : auto &data = ret.first->second;
686 24 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
687 8 : std::sort(data.attributes.begin(), data.attributes.end());
688 :
689 : // Cleanup
690 8 : it->erase("Index");
691 8 : it->erase("Attributes");
692 16 : MFEM_VERIFY(it->empty(),
693 : "Found an unsupported configuration file keyword under \"Energy\"!\n"
694 : << it->dump(2));
695 :
696 : // Debug
697 : if constexpr (JSON_DEBUG)
698 : {
699 : std::cout << "Index: " << ret.first->first << '\n';
700 : std::cout << "Attributes: " << data.attributes << '\n';
701 : }
702 : }
703 : }
704 :
705 8 : void ProbePostData::SetUp(json &postpro)
706 : {
707 8 : auto probe = postpro.find("Probe");
708 8 : if (probe == postpro.end())
709 : {
710 0 : return;
711 : }
712 8 : MFEM_VERIFY(probe->is_array(),
713 : "\"Probe\" should specify an array in the configuration file!");
714 24 : for (auto it = probe->begin(); it != probe->end(); ++it)
715 : {
716 16 : auto index = AtIndex(it, "\"Probe\" point");
717 16 : auto ctr = it->find("Center");
718 16 : MFEM_VERIFY(ctr != it->end() && ctr->is_array(),
719 : "Missing \"Probe\" point \"Center\" or \"Center\" should specify an array "
720 : "in the configuration file!");
721 48 : auto ret = mapdata.insert(std::make_pair(index, ProbeData()));
722 16 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"Probe\" points in "
723 : "the configuration file!");
724 : auto &data = ret.first->second;
725 16 : data.center = ctr->get<std::array<double, 3>>(); // Required
726 :
727 : // Cleanup
728 16 : it->erase("Index");
729 16 : it->erase("Center");
730 32 : MFEM_VERIFY(it->empty(),
731 : "Found an unsupported configuration file keyword under \"Probe\"!\n"
732 : << it->dump(2));
733 :
734 : // Debug
735 : if constexpr (JSON_DEBUG)
736 : {
737 : std::cout << "Index: " << ret.first->first << '\n';
738 : std::cout << "Center: " << data.center << '\n';
739 : }
740 : }
741 : }
742 :
743 8 : void DomainPostData::SetUp(json &domains)
744 : {
745 8 : auto postpro = domains.find("Postprocessing");
746 8 : if (postpro == domains.end())
747 : {
748 0 : return;
749 : }
750 8 : energy.SetUp(*postpro);
751 8 : probe.SetUp(*postpro);
752 :
753 : // Store all unique postprocessing domain attributes.
754 16 : for (const auto &[idx, data] : energy)
755 : {
756 8 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
757 : }
758 8 : std::sort(attributes.begin(), attributes.end());
759 8 : attributes.erase(std::unique(attributes.begin(), attributes.end()), attributes.end());
760 : attributes.shrink_to_fit();
761 :
762 : // Cleanup
763 8 : postpro->erase("Energy");
764 8 : postpro->erase("Probe");
765 16 : MFEM_VERIFY(postpro->empty(),
766 : "Found an unsupported configuration file keyword under \"Postprocessing\"!\n"
767 : << postpro->dump(2));
768 : }
769 :
770 8 : void DomainData::SetUp(json &config)
771 : {
772 8 : auto domains = config.find("Domains");
773 8 : MFEM_VERIFY(domains != config.end(),
774 : "\"Domains\" must be specified in the configuration file!");
775 8 : materials.SetUp(*domains);
776 8 : postpro.SetUp(*domains);
777 :
778 : // Store all unique domain attributes.
779 16 : for (const auto &data : materials)
780 : {
781 8 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
782 : }
783 8 : std::sort(attributes.begin(), attributes.end());
784 8 : attributes.erase(std::unique(attributes.begin(), attributes.end()), attributes.end());
785 : attributes.shrink_to_fit();
786 16 : for (const auto &attr : postpro.attributes)
787 : {
788 8 : MFEM_VERIFY(std::lower_bound(attributes.begin(), attributes.end(), attr) !=
789 : attributes.end(),
790 : "Domain postprocessing can only be enabled on domains which have a "
791 : "corresponding \"Materials\" entry!");
792 : }
793 :
794 : // Cleanup
795 8 : domains->erase("Materials");
796 8 : domains->erase("Postprocessing");
797 16 : MFEM_VERIFY(domains->empty(),
798 : "Found an unsupported configuration file keyword under \"Domains\"!\n"
799 : << domains->dump(2));
800 8 : }
801 :
802 39 : void PecBoundaryData::SetUp(json &boundaries)
803 : {
804 39 : auto pec = boundaries.find("PEC");
805 39 : auto ground = boundaries.find("Ground");
806 39 : if (pec == boundaries.end() && ground == boundaries.end())
807 : {
808 37 : return;
809 : }
810 2 : if (pec == boundaries.end())
811 : {
812 : pec = ground;
813 : }
814 2 : else if (ground == boundaries.end()) // Do nothing
815 : {
816 : }
817 : else
818 : {
819 0 : MFEM_ABORT(
820 : "Configuration file should not specify both \"PEC\" and \"Ground\" boundaries!");
821 : }
822 2 : MFEM_VERIFY(
823 : pec->find("Attributes") != pec->end(),
824 : "Missing \"Attributes\" list for \"PEC\" boundary in the configuration file!");
825 6 : attributes = pec->at("Attributes").get<std::vector<int>>(); // Required
826 2 : std::sort(attributes.begin(), attributes.end());
827 :
828 : // Cleanup
829 2 : pec->erase("Attributes");
830 4 : MFEM_VERIFY(pec->empty(),
831 : "Found an unsupported configuration file keyword under \"PEC\"!\n"
832 : << pec->dump(2));
833 :
834 : // Debug
835 : if constexpr (JSON_DEBUG)
836 : {
837 : std::cout << "PEC:" << attributes << '\n';
838 : }
839 : }
840 :
841 39 : void PmcBoundaryData::SetUp(json &boundaries)
842 : {
843 39 : auto pmc = boundaries.find("PMC");
844 39 : auto zeroq = boundaries.find("ZeroCharge");
845 39 : if (pmc == boundaries.end() && zeroq == boundaries.end())
846 : {
847 39 : return;
848 : }
849 0 : if (pmc == boundaries.end())
850 : {
851 : pmc = zeroq;
852 : }
853 0 : else if (zeroq == boundaries.end()) // Do nothing
854 : {
855 : }
856 : else
857 : {
858 0 : MFEM_ABORT("Configuration file should not specify both \"PMC\" and \"ZeroCharge\" "
859 : "boundaries!");
860 : }
861 0 : MFEM_VERIFY(
862 : pmc->find("Attributes") != pmc->end(),
863 : "Missing \"Attributes\" list for \"PMC\" boundary in the configuration file!");
864 0 : attributes = pmc->at("Attributes").get<std::vector<int>>(); // Required
865 0 : std::sort(attributes.begin(), attributes.end());
866 :
867 : // Cleanup
868 0 : pmc->erase("Attributes");
869 0 : MFEM_VERIFY(pmc->empty(),
870 : "Found an unsupported configuration file keyword under \"PMC\"!\n"
871 : << pmc->dump(2));
872 :
873 : // Debug
874 : if constexpr (JSON_DEBUG)
875 : {
876 : std::cout << "PMC:" << attributes << '\n';
877 : }
878 : }
879 :
880 39 : void WavePortPecBoundaryData::SetUp(json &boundaries)
881 : {
882 39 : auto pec = boundaries.find("WavePortPEC");
883 39 : if (pec == boundaries.end())
884 : {
885 39 : return;
886 : }
887 0 : MFEM_VERIFY(pec->find("Attributes") != pec->end(),
888 : "Missing \"Attributes\" list for \"WavePortPEC\" boundary in the "
889 : "configuration file!");
890 0 : attributes = pec->at("Attributes").get<std::vector<int>>(); // Required
891 0 : std::sort(attributes.begin(), attributes.end());
892 :
893 : // Cleanup
894 0 : pec->erase("Attributes");
895 0 : MFEM_VERIFY(pec->empty(),
896 : "Found an unsupported configuration file keyword under \"WavePortPEC\"!\n"
897 : << pec->dump(2));
898 :
899 : // Debug
900 : if constexpr (JSON_DEBUG)
901 : {
902 : std::cout << "WavePortPEC:" << attributes << '\n';
903 : }
904 : }
905 :
906 39 : void FarfieldBoundaryData::SetUp(json &boundaries)
907 : {
908 39 : auto absorbing = boundaries.find("Absorbing");
909 39 : if (absorbing == boundaries.end())
910 : {
911 37 : return;
912 : }
913 2 : MFEM_VERIFY(
914 : absorbing->find("Attributes") != absorbing->end(),
915 : "Missing \"Attributes\" list for \"Absorbing\" boundary in the configuration file!");
916 6 : attributes = absorbing->at("Attributes").get<std::vector<int>>(); // Required
917 2 : std::sort(attributes.begin(), attributes.end());
918 2 : order = absorbing->value("Order", order);
919 2 : MFEM_VERIFY(order == 1 || order == 2,
920 : "Supported values for config[\"Absorbing\"][\"Order\"] are 1 or 2!");
921 :
922 : // Cleanup
923 2 : absorbing->erase("Attributes");
924 2 : absorbing->erase("Order");
925 4 : MFEM_VERIFY(absorbing->empty(),
926 : "Found an unsupported configuration file keyword under \"Absorbing\"!\n"
927 : << absorbing->dump(2));
928 :
929 : // Debug
930 : if constexpr (JSON_DEBUG)
931 : {
932 : std::cout << "Absorbing:" << attributes << '\n';
933 : std::cout << "Order: " << order << '\n';
934 : }
935 : }
936 :
937 39 : void ConductivityBoundaryData::SetUp(json &boundaries)
938 : {
939 39 : auto conductivity = boundaries.find("Conductivity");
940 39 : if (conductivity == boundaries.end())
941 : {
942 39 : return;
943 : }
944 0 : MFEM_VERIFY(conductivity->is_array(),
945 : "\"Conductivity\" should specify an array in the configuration file!");
946 0 : for (auto it = conductivity->begin(); it != conductivity->end(); ++it)
947 : {
948 0 : MFEM_VERIFY(it->find("Attributes") != it->end(),
949 : "Missing \"Attributes\" list for \"Conductivity\" boundary in the "
950 : "configuration file!");
951 0 : MFEM_VERIFY(
952 : it->find("Conductivity") != it->end(),
953 : "Missing \"Conductivity\" boundary \"Conductivity\" in the configuration file!");
954 : ConductivityData &data = emplace_back();
955 0 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
956 0 : std::sort(data.attributes.begin(), data.attributes.end());
957 0 : data.sigma = it->at("Conductivity"); // Required
958 0 : data.mu_r = it->value("Permeability", data.mu_r);
959 0 : data.h = it->value("Thickness", data.h);
960 0 : data.external = it->value("External", data.external);
961 :
962 : // Cleanup
963 0 : it->erase("Attributes");
964 0 : it->erase("Conductivity");
965 0 : it->erase("Permeability");
966 0 : it->erase("Thickness");
967 0 : it->erase("External");
968 0 : MFEM_VERIFY(it->empty(),
969 : "Found an unsupported configuration file keyword under \"Conductivity\"!\n"
970 : << it->dump(2));
971 :
972 : // Debug
973 : if constexpr (JSON_DEBUG)
974 : {
975 : std::cout << "Attributes: " << data.attributes << '\n';
976 : std::cout << "Conductivity: " << data.sigma << '\n';
977 : std::cout << "Permeability: " << data.mu_r << '\n';
978 : std::cout << "Thickness: " << data.h << '\n';
979 : std::cout << "External: " << data.external << '\n';
980 : }
981 : }
982 : }
983 :
984 39 : void ImpedanceBoundaryData::SetUp(json &boundaries)
985 : {
986 39 : auto impedance = boundaries.find("Impedance");
987 39 : if (impedance == boundaries.end())
988 : {
989 39 : return;
990 : }
991 0 : MFEM_VERIFY(impedance->is_array(),
992 : "\"Impedance\" should specify an array in the configuration file!");
993 0 : for (auto it = impedance->begin(); it != impedance->end(); ++it)
994 : {
995 0 : MFEM_VERIFY(it->find("Attributes") != it->end(),
996 : "Missing \"Attributes\" list for \"Impedance\" boundary in the "
997 : "configuration file!");
998 : ImpedanceData &data = emplace_back();
999 0 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
1000 0 : std::sort(data.attributes.begin(), data.attributes.end());
1001 0 : data.Rs = it->value("Rs", data.Rs);
1002 0 : data.Ls = it->value("Ls", data.Ls);
1003 0 : data.Cs = it->value("Cs", data.Cs);
1004 :
1005 : // Cleanup
1006 0 : it->erase("Attributes");
1007 0 : it->erase("Rs");
1008 0 : it->erase("Ls");
1009 0 : it->erase("Cs");
1010 0 : MFEM_VERIFY(it->empty(),
1011 : "Found an unsupported configuration file keyword under \"Impedance\"!\n"
1012 : << it->dump(2));
1013 :
1014 : // Debug
1015 : if constexpr (JSON_DEBUG)
1016 : {
1017 : std::cout << "Attributes: " << data.attributes << '\n';
1018 : std::cout << "Rs: " << data.Rs << '\n';
1019 : std::cout << "Ls: " << data.Ls << '\n';
1020 : std::cout << "Cs: " << data.Cs << '\n';
1021 : }
1022 : }
1023 : }
1024 :
1025 81 : int ParsePortExcitation(json::iterator &port_it, int default_excitation)
1026 : {
1027 81 : auto it_excitation = port_it->find("Excitation");
1028 81 : if (it_excitation == port_it->end())
1029 : {
1030 : // Keep default; don't set input flag.
1031 : return default_excitation;
1032 : }
1033 51 : else if (it_excitation->is_boolean())
1034 : {
1035 58 : return int(it_excitation->get<bool>()); // 0 false; 1 true
1036 : }
1037 22 : else if (it_excitation->is_number_unsigned())
1038 : {
1039 40 : return it_excitation->get<int>();
1040 : }
1041 : else
1042 : {
1043 10 : MFEM_ABORT(fmt::format("\"Excitation\" on port index {:d} could not be parsed "
1044 : "as a bool or unsigned (non-negative) integer; got {}",
1045 : int(port_it->at("Index")), it_excitation->dump(2)));
1046 : }
1047 : }
1048 :
1049 39 : void LumpedPortBoundaryData::SetUp(json &boundaries)
1050 : {
1051 39 : auto port = boundaries.find("LumpedPort");
1052 39 : auto terminal = boundaries.find("Terminal");
1053 39 : if (port == boundaries.end() && terminal == boundaries.end())
1054 : {
1055 1 : return;
1056 : }
1057 :
1058 38 : if (port == boundaries.end())
1059 : {
1060 : port = terminal;
1061 : }
1062 38 : else if (terminal == boundaries.end()) // Do nothing
1063 : {
1064 : }
1065 : else
1066 : {
1067 0 : MFEM_ABORT("Configuration file should not specify both \"LumpedPort\" and \"Terminal\" "
1068 : "boundaries!");
1069 : }
1070 :
1071 76 : std::string label = (terminal != boundaries.end()) ? "\"Terminal\"" : "\"LumpedPort\"";
1072 38 : MFEM_VERIFY(port->is_array(),
1073 : label << " should specify an array in the configuration file!");
1074 99 : for (auto it = port->begin(); it != port->end(); ++it)
1075 : {
1076 64 : auto index = AtIndex(it, label);
1077 63 : auto ret = mapdata.insert(std::make_pair(index, LumpedPortData()));
1078 66 : MFEM_VERIFY(ret.second, fmt::format("Repeated \"Index\" found when processing {} "
1079 : "boundaries in the configuration file!",
1080 : label));
1081 : auto &data = ret.first->second;
1082 62 : data.R = it->value("R", data.R);
1083 62 : data.L = it->value("L", data.L);
1084 62 : data.C = it->value("C", data.C);
1085 62 : data.Rs = it->value("Rs", data.Rs);
1086 62 : data.Ls = it->value("Ls", data.Ls);
1087 62 : data.Cs = it->value("Cs", data.Cs);
1088 :
1089 62 : data.excitation = ParsePortExcitation(it, data.excitation);
1090 61 : data.active = it->value("Active", data.active);
1091 61 : if (it->find("Attributes") != it->end())
1092 : {
1093 59 : MFEM_VERIFY(
1094 : it->find("Elements") == it->end(),
1095 : fmt::format(
1096 : "Cannot specify both top-level \"Attributes\" list and \"Elements\" for "
1097 : "{} in the configuration file!",
1098 : label));
1099 59 : auto &elem = data.elements.emplace_back();
1100 118 : ParseElementData(*it, "Direction", terminal == boundaries.end(), elem);
1101 : }
1102 : else
1103 : {
1104 2 : auto elements = it->find("Elements");
1105 2 : MFEM_VERIFY(elements != it->end(),
1106 : fmt::format("Missing top-level \"Attributes\" list or \"Elements\" for "
1107 : "{} in the configuration file!",
1108 : label));
1109 6 : for (auto elem_it = elements->begin(); elem_it != elements->end(); ++elem_it)
1110 : {
1111 4 : MFEM_VERIFY(elem_it->find("Attributes") != elem_it->end(),
1112 : fmt::format("Missing \"Attributes\" list for {} element in "
1113 : "the configuration file!",
1114 : label));
1115 4 : auto &elem = data.elements.emplace_back();
1116 7 : ParseElementData(*elem_it, "Direction", terminal == boundaries.end(), elem);
1117 :
1118 : // Cleanup
1119 4 : elem_it->erase("Attributes");
1120 4 : elem_it->erase("Direction");
1121 4 : elem_it->erase("CoordinateSystem");
1122 8 : MFEM_VERIFY(elem_it->empty(), fmt::format("Found an unsupported configuration file "
1123 : "keyword under {} element!\n{}",
1124 : label, elem_it->dump(2)));
1125 : }
1126 : }
1127 :
1128 : // Cleanup
1129 61 : it->erase("Index");
1130 61 : it->erase("R");
1131 61 : it->erase("L");
1132 61 : it->erase("C");
1133 61 : it->erase("Rs");
1134 61 : it->erase("Ls");
1135 61 : it->erase("Cs");
1136 61 : it->erase("Excitation");
1137 61 : it->erase("Active");
1138 61 : it->erase("Attributes");
1139 61 : it->erase("Direction");
1140 61 : it->erase("CoordinateSystem");
1141 61 : it->erase("Elements");
1142 122 : MFEM_VERIFY(it->empty(),
1143 : fmt::format("Found an unsupported configuration file keyword under {}!\n{}",
1144 : label, it->dump(2)));
1145 :
1146 : // Debug
1147 : if constexpr (JSON_DEBUG)
1148 : {
1149 : std::cout << "Index: " << ret.first->first << '\n';
1150 : std::cout << "R: " << data.R << '\n';
1151 : std::cout << "L: " << data.L << '\n';
1152 : std::cout << "C: " << data.C << '\n';
1153 : std::cout << "Rs: " << data.Rs << '\n';
1154 : std::cout << "Ls: " << data.Ls << '\n';
1155 : std::cout << "Cs: " << data.Cs << '\n';
1156 : std::cout << "Excitation: " << data.excitation << '\n';
1157 : std::cout << "Active: " << data.active << '\n';
1158 : for (const auto &elem : data.elements)
1159 : {
1160 : std::cout << "Attributes: " << elem.attributes << '\n';
1161 : std::cout << "Direction: " << elem.direction << '\n';
1162 : std::cout << "CoordinateSystem: " << elem.coordinate_system << '\n';
1163 : }
1164 : }
1165 : }
1166 : }
1167 :
1168 36 : void PeriodicBoundaryData::SetUp(json &boundaries)
1169 : {
1170 36 : auto periodic = boundaries.find("Periodic");
1171 36 : if (periodic == boundaries.end())
1172 : {
1173 35 : return;
1174 : }
1175 1 : auto floquet = periodic->find("FloquetWaveVector");
1176 1 : if (floquet != periodic->end())
1177 : {
1178 0 : MFEM_VERIFY(floquet->is_array(),
1179 : "\"FloquetWaveVector\" should specify an array in the configuration file!");
1180 0 : wave_vector = floquet->get<std::array<double, 3>>();
1181 : }
1182 :
1183 : // Debug
1184 : if constexpr (JSON_DEBUG)
1185 : {
1186 : std::cout << "FloquetWaveVector: " << wave_vector << '\n';
1187 : }
1188 :
1189 1 : auto pairs = periodic->find("BoundaryPairs");
1190 1 : MFEM_VERIFY(pairs->is_array(),
1191 : "\"BoundaryPairs\" should specify an array in the configuration file!");
1192 2 : for (auto it = pairs->begin(); it != pairs->end(); ++it)
1193 : {
1194 1 : MFEM_VERIFY(it->find("DonorAttributes") != it->end(),
1195 : "Missing \"DonorAttributes\" list for \"Periodic\" boundary in the "
1196 : "configuration file!");
1197 1 : MFEM_VERIFY(it->find("ReceiverAttributes") != it->end(),
1198 : "Missing \"ReceiverAttributes\" list for \"Periodic\" boundary in the "
1199 : "configuration file!");
1200 :
1201 1 : PeriodicData &data = boundary_pairs.emplace_back();
1202 2 : data.donor_attributes = it->at("DonorAttributes").get<std::vector<int>>(); // Required
1203 : data.receiver_attributes =
1204 2 : it->at("ReceiverAttributes").get<std::vector<int>>(); // Required
1205 1 : auto translation = it->find("Translation");
1206 1 : if (translation != it->end())
1207 : {
1208 0 : MFEM_VERIFY(translation->is_array(),
1209 : "\"Translation\" should specify an array in the configuration file!");
1210 0 : std::array<double, 3> translation_array = translation->get<std::array<double, 3>>();
1211 0 : for (int i = 0; i < 3; i++)
1212 : {
1213 0 : data.affine_transform[i * 4 + i] = 1.0;
1214 0 : data.affine_transform[i * 4 + 3] = translation_array[i];
1215 : }
1216 0 : data.affine_transform[3 * 4 + 3] = 1.0;
1217 : }
1218 1 : auto transformation = it->find("AffineTransformation");
1219 1 : if (transformation != it->end())
1220 : {
1221 0 : MFEM_VERIFY(
1222 : transformation->is_array(),
1223 : "\"AffineTransformation\" should specify an array in the configuration file!");
1224 0 : data.affine_transform = transformation->get<std::array<double, 16>>();
1225 : }
1226 :
1227 : // Cleanup
1228 1 : it->erase("DonorAttributes");
1229 1 : it->erase("ReceiverAttributes");
1230 1 : it->erase("Translation");
1231 1 : it->erase("AffineTransformation");
1232 2 : MFEM_VERIFY(it->empty(),
1233 : "Found an unsupported configuration file keyword under \"Periodic\"!\n"
1234 : << it->dump(2));
1235 :
1236 : // Debug
1237 : if constexpr (JSON_DEBUG)
1238 : {
1239 : std::cout << "DonorAttributes: " << data.donor_attributes << '\n';
1240 : std::cout << "ReceiverAttributes: " << data.receiver_attributes << '\n';
1241 : std::cout << "AffineTransformation: " << data.affine_transform << '\n';
1242 : }
1243 : }
1244 : }
1245 :
1246 36 : void WavePortBoundaryData::SetUp(json &boundaries)
1247 : {
1248 36 : auto port = boundaries.find("WavePort");
1249 36 : if (port == boundaries.end())
1250 : {
1251 24 : return;
1252 : }
1253 12 : MFEM_VERIFY(port->is_array(),
1254 : "\"WavePort\" should specify an array in the configuration file!");
1255 :
1256 30 : for (auto it = port->begin(); it != port->end(); ++it)
1257 : {
1258 21 : MFEM_VERIFY(
1259 : it->find("Attributes") != it->end(),
1260 : "Missing \"Attributes\" list for \"WavePort\" boundary in the configuration file!");
1261 21 : auto index = AtIndex(it, "\"WavePort\" boundary");
1262 20 : auto ret = mapdata.insert(std::make_pair(index, WavePortData()));
1263 23 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"WavePort\" "
1264 : "boundaries in the configuration file!");
1265 : auto &data = ret.first->second;
1266 57 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
1267 19 : std::sort(data.attributes.begin(), data.attributes.end());
1268 19 : data.mode_idx = it->value("Mode", data.mode_idx);
1269 19 : MFEM_VERIFY(data.mode_idx > 0,
1270 : "\"WavePort\" boundary \"Mode\" must be positive (1-based)!");
1271 19 : data.d_offset = it->value("Offset", data.d_offset);
1272 19 : data.eigen_solver = it->value("SolverType", data.eigen_solver);
1273 :
1274 19 : data.excitation = ParsePortExcitation(it, data.excitation);
1275 18 : data.active = it->value("Active", data.active);
1276 18 : data.ksp_max_its = it->value("MaxIts", data.ksp_max_its);
1277 18 : data.ksp_tol = it->value("KSPTol", data.ksp_tol);
1278 18 : data.eig_tol = it->value("EigenTol", data.eig_tol);
1279 18 : data.verbose = it->value("Verbose", data.verbose);
1280 :
1281 : // Cleanup
1282 18 : it->erase("Index");
1283 18 : it->erase("Attributes");
1284 18 : it->erase("Mode");
1285 18 : it->erase("Offset");
1286 18 : it->erase("SolverType");
1287 18 : it->erase("Excitation");
1288 18 : it->erase("Active");
1289 18 : it->erase("MaxIts");
1290 18 : it->erase("KSPTol");
1291 18 : it->erase("EigenTol");
1292 18 : it->erase("Verbose");
1293 36 : MFEM_VERIFY(it->empty(),
1294 : "Found an unsupported configuration file keyword under \"WavePort\"!\n"
1295 : << it->dump(2));
1296 :
1297 : // Debug
1298 : if constexpr (JSON_DEBUG)
1299 : {
1300 : std::cout << "Index: " << ret.first->first << '\n';
1301 : std::cout << "Attributes: " << data.attributes << '\n';
1302 : std::cout << "Mode: " << data.mode_idx << '\n';
1303 : std::cout << "Offset: " << data.d_offset << '\n';
1304 : std::cout << "SolverType: " << data.eigen_solver << '\n';
1305 : std::cout << "Excitation: " << data.excitation << '\n';
1306 : std::cout << "Active: " << data.active << '\n';
1307 : std::cout << "MaxIts: " << data.ksp_max_its << '\n';
1308 : std::cout << "KSPTol: " << data.ksp_tol << '\n';
1309 : std::cout << "EigenTol: " << data.eig_tol << '\n';
1310 : std::cout << "Verbose: " << data.verbose << '\n';
1311 : }
1312 : }
1313 : }
1314 :
1315 33 : void SurfaceCurrentBoundaryData::SetUp(json &boundaries)
1316 : {
1317 33 : auto source = boundaries.find("SurfaceCurrent");
1318 33 : if (source == boundaries.end())
1319 : {
1320 33 : return;
1321 : }
1322 0 : MFEM_VERIFY(source->is_array(),
1323 : "\"SurfaceCurrent\" should specify an array in the configuration file!");
1324 0 : for (auto it = source->begin(); it != source->end(); ++it)
1325 : {
1326 0 : auto index = AtIndex(it, "\"SurfaceCurrent\" source");
1327 0 : auto ret = mapdata.insert(std::make_pair(index, SurfaceCurrentData()));
1328 0 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"SurfaceCurrent\" "
1329 : "boundaries in the configuration file!");
1330 : auto &data = ret.first->second;
1331 0 : if (it->find("Attributes") != it->end())
1332 : {
1333 0 : MFEM_VERIFY(it->find("Elements") == it->end(),
1334 : "Cannot specify both top-level \"Attributes\" list and \"Elements\" for "
1335 : "\"SurfaceCurrent\" boundary in the configuration file!");
1336 0 : auto &elem = data.elements.emplace_back();
1337 0 : ParseElementData(*it, "Direction", true, elem);
1338 : }
1339 : else
1340 : {
1341 0 : auto elements = it->find("Elements");
1342 0 : MFEM_VERIFY(
1343 : elements != it->end(),
1344 : "Missing top-level \"Attributes\" list or \"Elements\" for \"SurfaceCurrent\" "
1345 : "boundary in the configuration file!");
1346 0 : for (auto elem_it = elements->begin(); elem_it != elements->end(); ++elem_it)
1347 : {
1348 0 : MFEM_VERIFY(
1349 : elem_it->find("Attributes") != elem_it->end(),
1350 : "Missing \"Attributes\" list for \"SurfaceCurrent\" boundary element in "
1351 : "configuration file!");
1352 0 : auto &elem = data.elements.emplace_back();
1353 0 : ParseElementData(*elem_it, "Direction", true, elem);
1354 :
1355 : // Cleanup
1356 0 : elem_it->erase("Attributes");
1357 0 : elem_it->erase("Direction");
1358 0 : elem_it->erase("CoordinateSystem");
1359 0 : MFEM_VERIFY(elem_it->empty(), "Found an unsupported configuration file keyword "
1360 : "under \"SurfaceCurrent\" boundary element!\n"
1361 : << elem_it->dump(2));
1362 : }
1363 : }
1364 :
1365 : // Cleanup
1366 0 : it->erase("Index");
1367 0 : it->erase("Attributes");
1368 0 : it->erase("Direction");
1369 0 : it->erase("Elements");
1370 0 : it->erase("CoordinateSystem");
1371 0 : MFEM_VERIFY(
1372 : it->empty(),
1373 : "Found an unsupported configuration file keyword under \"SurfaceCurrent\"!\n"
1374 : << it->dump(2));
1375 :
1376 : // Debug
1377 : if constexpr (JSON_DEBUG)
1378 : {
1379 : std::cout << "Index: " << ret.first->first << '\n';
1380 : for (const auto &elem : data.elements)
1381 : {
1382 : std::cout << "Attributes: " << elem.attributes << '\n';
1383 : std::cout << "Direction: " << elem.direction << '\n';
1384 : std::cout << "CoordinateSystem: " << elem.coordinate_system << '\n';
1385 : }
1386 : }
1387 : }
1388 : }
1389 :
1390 8 : void SurfaceFluxPostData::SetUp(json &postpro)
1391 : {
1392 8 : auto flux = postpro.find("SurfaceFlux");
1393 8 : if (flux == postpro.end())
1394 : {
1395 0 : return;
1396 : }
1397 8 : MFEM_VERIFY(flux->is_array(),
1398 : "\"SurfaceFlux\" should specify an array in the configuration file!");
1399 24 : for (auto it = flux->begin(); it != flux->end(); ++it)
1400 : {
1401 16 : auto index = AtIndex(it, "\"SurfaceFlux\" boundary");
1402 16 : MFEM_VERIFY(it->find("Attributes") != it->end() && it->find("Type") != it->end(),
1403 : "Missing \"Attributes\" list or \"Type\" for \"SurfaceFlux\" boundary "
1404 : "in the configuration file!");
1405 16 : auto ret = mapdata.insert(std::make_pair(index, SurfaceFluxData()));
1406 16 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"SurfaceFlux\" "
1407 : "boundaries in the configuration file!");
1408 : auto &data = ret.first->second;
1409 48 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
1410 16 : std::sort(data.attributes.begin(), data.attributes.end());
1411 16 : data.type = it->at("Type"); // Required
1412 16 : data.two_sided = it->value("TwoSided", data.two_sided);
1413 16 : auto ctr = it->find("Center");
1414 16 : if (ctr != it->end())
1415 : {
1416 0 : MFEM_VERIFY(ctr->is_array(),
1417 : "\"Center\" should specify an array in the configuration file!");
1418 0 : data.center = ctr->get<std::array<double, 3>>();
1419 0 : data.no_center = false;
1420 : }
1421 :
1422 : // Cleanup
1423 16 : it->erase("Index");
1424 16 : it->erase("Attributes");
1425 16 : it->erase("Type");
1426 16 : it->erase("TwoSided");
1427 16 : it->erase("Center");
1428 32 : MFEM_VERIFY(it->empty(),
1429 : "Found an unsupported configuration file keyword under \"SurfaceFlux\"!\n"
1430 : << it->dump(2));
1431 :
1432 : // Debug
1433 : if constexpr (JSON_DEBUG)
1434 : {
1435 : std::cout << "Index: " << ret.first->first << '\n';
1436 : std::cout << "Attributes: " << data.attributes << '\n';
1437 : std::cout << "Type: " << data.type << '\n';
1438 : std::cout << "TwoSided: " << data.two_sided << '\n';
1439 : std::cout << "Center: " << data.center << '\n';
1440 : }
1441 : }
1442 : }
1443 :
1444 8 : void InterfaceDielectricPostData::SetUp(json &postpro)
1445 : {
1446 8 : auto dielectric = postpro.find("Dielectric");
1447 8 : if (dielectric == postpro.end())
1448 : {
1449 0 : return;
1450 : }
1451 8 : MFEM_VERIFY(dielectric->is_array(),
1452 : "\"Dielectric\" should specify an array in the configuration file!");
1453 32 : for (auto it = dielectric->begin(); it != dielectric->end(); ++it)
1454 : {
1455 24 : auto index = AtIndex(it, "\"Dielectric\" boundary");
1456 24 : MFEM_VERIFY(it->find("Attributes") != it->end() && it->find("Thickness") != it->end() &&
1457 : it->find("Permittivity") != it->end(),
1458 : "Missing \"Dielectric\" boundary \"Attributes\" list, \"Thickness\", or "
1459 : "\"Permittivity\" in the configuration file!");
1460 24 : auto ret = mapdata.insert(std::make_pair(index, InterfaceDielectricData()));
1461 24 : MFEM_VERIFY(ret.second, "Repeated \"Index\" found when processing \"Dielectric\" "
1462 : "boundaries in the configuration file!");
1463 : auto &data = ret.first->second;
1464 72 : data.attributes = it->at("Attributes").get<std::vector<int>>(); // Required
1465 24 : std::sort(data.attributes.begin(), data.attributes.end());
1466 24 : data.type = it->value("Type", data.type);
1467 24 : data.t = it->at("Thickness"); // Required
1468 24 : data.epsilon_r = it->at("Permittivity"); // Required
1469 24 : data.tandelta = it->value("LossTan", data.tandelta);
1470 :
1471 : // Cleanup
1472 24 : it->erase("Index");
1473 24 : it->erase("Attributes");
1474 24 : it->erase("Type");
1475 24 : it->erase("Thickness");
1476 24 : it->erase("Permittivity");
1477 24 : it->erase("LossTan");
1478 48 : MFEM_VERIFY(it->empty(),
1479 : "Found an unsupported configuration file keyword under \"Dielectric\"!\n"
1480 : << it->dump(2));
1481 :
1482 : // Debug
1483 : if constexpr (JSON_DEBUG)
1484 : {
1485 : std::cout << "Index: " << ret.first->first << '\n';
1486 : std::cout << "Attributes: " << data.attributes << '\n';
1487 : std::cout << "Type: " << data.type << '\n';
1488 : std::cout << "Thickness: " << data.t << '\n';
1489 : std::cout << "Permittivity: " << data.epsilon_r << '\n';
1490 : std::cout << "LossTan: " << data.tandelta << '\n';
1491 : }
1492 : }
1493 : }
1494 25 : void FarFieldPostData::SetUp(json &postpro)
1495 : {
1496 25 : auto farfield = postpro.find("FarField");
1497 25 : if (farfield == postpro.end())
1498 : {
1499 9 : return;
1500 : }
1501 :
1502 16 : MFEM_VERIFY(farfield->find("Attributes") != farfield->end(),
1503 : "Missing \"Attributes\" list for \"FarField\" postprocessing in the "
1504 : "configuration file!");
1505 :
1506 48 : attributes = farfield->at("Attributes").get<std::vector<int>>(); // Required
1507 16 : std::sort(attributes.begin(), attributes.end());
1508 :
1509 : // Generate NSample points with the following properties:
1510 : // - If NSample >= 2, the generated points are precisely NSample, otherwise NSample = 2.
1511 : // - The poles, the equator, and the XZ plane are always included.
1512 : // - The points are almost uniformly on a sphere, with a small bias due to satisfying the
1513 : // previous condition.
1514 : // - The points are on rings of constant theta.
1515 :
1516 16 : auto nsample_json = farfield->find("NSample");
1517 : int nsample = 0;
1518 16 : if (nsample_json != farfield->end())
1519 : {
1520 10 : nsample = nsample_json->get<int>();
1521 10 : if (nsample > 0)
1522 : {
1523 : // Always include poles.
1524 9 : thetaphis.emplace_back(0.0, 0.0); // North pole.
1525 9 : thetaphis.emplace_back(M_PI, 0.0); // South pole.
1526 :
1527 9 : if (nsample > 2)
1528 : {
1529 7 : int remaining = nsample - 2;
1530 :
1531 : // Distribute all remaining points across rings with number weighted by the
1532 : // local circumference.
1533 7 : int n_theta = std::max(1, static_cast<int>(std::sqrt(remaining)));
1534 7 : n_theta = std::min(n_theta, remaining); // Can't have more rings than points.
1535 :
1536 7 : std::vector<int> points_per_level(n_theta);
1537 7 : std::vector<double> sin_theta_values(n_theta);
1538 : double total_sin_theta = 0.0;
1539 :
1540 : // Calculate sin(theta) for each ring and total (sin(theta) is proportional to the
1541 : // circumference).
1542 278 : for (int i = 0; i < n_theta; ++i)
1543 : {
1544 271 : double theta = std::acos(1.0 - 2.0 * (i + 1) / (n_theta + 1.0));
1545 271 : sin_theta_values[i] = std::sin(theta);
1546 271 : total_sin_theta += sin_theta_values[i];
1547 : }
1548 :
1549 : // Distribute points proportional to sin(theta).
1550 : int assigned_points = 0;
1551 271 : for (int i = 0; i < n_theta - 1; ++i)
1552 : {
1553 264 : points_per_level[i] =
1554 264 : static_cast<int>(remaining * sin_theta_values[i] / total_sin_theta + 0.5);
1555 264 : assigned_points += points_per_level[i];
1556 : }
1557 : // Assign remaining points to last ring to ensure exact total.
1558 7 : points_per_level[n_theta - 1] = remaining - assigned_points;
1559 :
1560 278 : for (int i = 1; i <= n_theta; ++i)
1561 : {
1562 : // Ensure equator and XZ plane inclusion.
1563 271 : bool is_equator = (i == (n_theta + 1) / 2);
1564 271 : double theta = is_equator ? M_PI / 2 : std::acos(1.0 - 2.0 * i / (n_theta + 1.0));
1565 271 : int points_in_level = points_per_level[i - 1];
1566 :
1567 65143 : for (int j = 0; j < points_in_level; ++j)
1568 : {
1569 64872 : double phi = 2.0 * M_PI * j / points_in_level;
1570 :
1571 : // Force XZ plane points (phi = 0 or π).
1572 64872 : if (j == 0)
1573 : {
1574 271 : phi = 0.0;
1575 : }
1576 64601 : else if (j == points_in_level / 2)
1577 : {
1578 271 : phi = M_PI;
1579 : }
1580 :
1581 64872 : thetaphis.emplace_back(theta, phi);
1582 : }
1583 : }
1584 : }
1585 :
1586 : if (nsample > 2)
1587 : MFEM_ASSERT(thetaphis.size() == nsample,
1588 : "Sampled number of points is not NSample!");
1589 : }
1590 : }
1591 :
1592 16 : auto thetaphis_json = farfield->find("ThetaPhis");
1593 16 : if (thetaphis_json != farfield->end())
1594 : {
1595 6 : MFEM_VERIFY(thetaphis_json->is_array(),
1596 : "\"ThetaPhis\" should specify an array in the configuration file!");
1597 :
1598 : // JSON does not support the notion of pair, so we read the theta and phis as vectors
1599 : // of vectors, and then cast them to vectors of pairs.
1600 : //
1601 : // Convert to radians in the process.
1602 6 : auto vec_of_vec = thetaphis_json->get<std::vector<std::vector<double>>>();
1603 20 : for (const auto &vec : vec_of_vec)
1604 : {
1605 14 : thetaphis.emplace_back(vec[0] * M_PI / 180, vec[1] * M_PI / 180);
1606 : }
1607 6 : }
1608 :
1609 : // Remove duplicate entries with numerical tolerance.
1610 : constexpr double tol = 1e-6;
1611 16 : std::sort(thetaphis.begin(), thetaphis.end());
1612 16 : auto it = std::unique(thetaphis.begin(), thetaphis.end(),
1613 64890 : [tol](const auto &a, const auto &b)
1614 : {
1615 : // At poles (theta ≈ 0 or π), phi is irrelevant.
1616 64890 : if ((std::abs(a.first) < tol || std::abs(a.first - M_PI) < tol) &&
1617 14 : (std::abs(b.first) < tol || std::abs(b.first - M_PI) < tol))
1618 : {
1619 5 : return std::abs(a.first - b.first) < tol;
1620 : }
1621 :
1622 : // Check direct match.
1623 64885 : if (std::abs(a.first - b.first) < tol)
1624 : {
1625 64603 : double phi_diff = std::abs(a.second - b.second);
1626 64603 : return phi_diff < tol || std::abs(phi_diff - 2.0 * M_PI) < tol;
1627 : }
1628 :
1629 : // Check theta periodicity: (θ, φ) ≡ (π-θ, φ+π).
1630 282 : if (std::abs(a.first - (M_PI - b.first)) < tol)
1631 : {
1632 1 : double phi_diff = std::abs(a.second - (b.second + M_PI));
1633 1 : if (phi_diff > M_PI)
1634 1 : phi_diff = 2.0 * M_PI - phi_diff;
1635 1 : return phi_diff < tol;
1636 : }
1637 :
1638 : return false;
1639 : });
1640 16 : thetaphis.erase(it, thetaphis.end());
1641 :
1642 16 : if (thetaphis.empty())
1643 : {
1644 8 : MFEM_WARNING("No target points specified under farfield \"FarField\"!\n");
1645 : }
1646 :
1647 : // Cleanup
1648 16 : farfield->erase("Attributes");
1649 16 : farfield->erase("NSample");
1650 16 : farfield->erase("ThetaPhis");
1651 32 : MFEM_VERIFY(farfield->empty(),
1652 : "Found an unsupported configuration file keyword under \"FarField\"!\n"
1653 : << farfield->dump(2));
1654 :
1655 : // Debug
1656 : if constexpr (JSON_DEBUG)
1657 : {
1658 : std::cout << "Attributes: " << attributes << '\n';
1659 : std::cout << "NSample: " << nsample << '\n';
1660 : std::cout << "ThetaPhis: " << thetaphis << '\n';
1661 : }
1662 : }
1663 33 : void BoundaryPostData::SetUp(json &boundaries)
1664 : {
1665 33 : auto postpro = boundaries.find("Postprocessing");
1666 33 : if (postpro == boundaries.end())
1667 : {
1668 25 : return;
1669 : }
1670 8 : flux.SetUp(*postpro);
1671 8 : dielectric.SetUp(*postpro);
1672 8 : farfield.SetUp(*postpro);
1673 :
1674 : // Store all unique postprocessing boundary attributes.
1675 24 : for (const auto &[idx, data] : flux)
1676 : {
1677 16 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
1678 : }
1679 32 : for (const auto &[idx, data] : dielectric)
1680 : {
1681 24 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
1682 : }
1683 :
1684 8 : attributes.insert(attributes.end(), farfield.attributes.begin(),
1685 : farfield.attributes.end());
1686 :
1687 8 : std::sort(attributes.begin(), attributes.end());
1688 8 : attributes.erase(std::unique(attributes.begin(), attributes.end()), attributes.end());
1689 : attributes.shrink_to_fit();
1690 :
1691 : // Cleanup
1692 8 : postpro->erase("SurfaceFlux");
1693 8 : postpro->erase("Dielectric");
1694 8 : postpro->erase("FarField");
1695 16 : MFEM_VERIFY(postpro->empty(),
1696 : "Found an unsupported configuration file keyword under \"Postprocessing\"!\n"
1697 : << postpro->dump(2));
1698 : }
1699 :
1700 39 : void BoundaryData::SetUp(json &config)
1701 : {
1702 39 : auto boundaries = config.find("Boundaries");
1703 39 : MFEM_VERIFY(boundaries != config.end(),
1704 : "\"Boundaries\" must be specified in the configuration file!");
1705 39 : pec.SetUp(*boundaries);
1706 39 : pmc.SetUp(*boundaries);
1707 39 : auxpec.SetUp(*boundaries);
1708 39 : farfield.SetUp(*boundaries);
1709 39 : conductivity.SetUp(*boundaries);
1710 39 : impedance.SetUp(*boundaries);
1711 39 : lumpedport.SetUp(*boundaries);
1712 36 : periodic.SetUp(*boundaries);
1713 36 : waveport.SetUp(*boundaries);
1714 33 : current.SetUp(*boundaries);
1715 33 : postpro.SetUp(*boundaries);
1716 :
1717 : // Ensure unique indexing of lumpedport, waveport, current.
1718 : {
1719 : std::map<int, std::string> index_map;
1720 : std::map<int, std::vector<int>> excitation_map;
1721 36 : const std::string lumpedport_str = "\"LumpedPort\"";
1722 36 : const std::string waveport_str = "WavePort";
1723 36 : const std::string current_str = "SurfaceCurrent";
1724 :
1725 90 : for (const auto &data : lumpedport)
1726 : {
1727 60 : auto result = index_map.insert({data.first, lumpedport_str});
1728 57 : MFEM_VERIFY(result.second, "Duplicate \"Index\": " << data.first << " in "
1729 : << index_map[data.first] << "!");
1730 57 : excitation_map[data.second.excitation].emplace_back(data.first);
1731 : }
1732 49 : for (const auto &data : waveport)
1733 : {
1734 20 : auto result = index_map.insert({data.first, waveport_str});
1735 22 : MFEM_VERIFY(result.second, "Duplicate \"Index\": " << data.first << " in "
1736 : << index_map[data.first] << "!");
1737 16 : excitation_map[data.second.excitation].emplace_back(data.first);
1738 : }
1739 32 : for (const auto &data : current)
1740 : {
1741 3 : auto result = index_map.insert({data.first, current_str});
1742 0 : MFEM_VERIFY(result.second, "Duplicate \"Index\": " << data.first << " in "
1743 : << index_map[data.first] << "!");
1744 : }
1745 : // Typical usecase: If each excitation is simple, S-parameters will be calculated.
1746 : // If there were multiple excitations specified, check their indices match the
1747 : // port indices. If there was only one, assign it.
1748 32 : excitation_map.erase(0); // zeroth index is unexcited.
1749 : bool calc_s_params = std::all_of(excitation_map.begin(), excitation_map.end(),
1750 : [](const auto &x) { return x.second.size() == 1; });
1751 32 : if (calc_s_params && !excitation_map.empty())
1752 : {
1753 : // If there's one excitation, needs to be 1 (set with bool) or the port index.
1754 : const auto &ext1 = *excitation_map.begin();
1755 42 : MFEM_VERIFY(
1756 : (excitation_map.size() == 1 &&
1757 : (ext1.first == 1 || ext1.second[0] == ext1.first)) ||
1758 : std::all_of(excitation_map.begin(), excitation_map.end(),
1759 : [](const auto &x) { return x.first == x.second[0]; }),
1760 : "\"Excitation\" must match \"Index\" for single ports to avoid ambiguity!");
1761 :
1762 71 : for (auto &[port_idx, lp] : lumpedport)
1763 : {
1764 45 : if (lp.excitation == 1)
1765 : {
1766 25 : lp.excitation = port_idx;
1767 : }
1768 : }
1769 32 : for (auto &[port_idx, wp] : waveport)
1770 : {
1771 6 : if (wp.excitation == 1)
1772 : {
1773 1 : wp.excitation = port_idx;
1774 : }
1775 : }
1776 : }
1777 : }
1778 :
1779 : // Store all unique boundary attributes.
1780 30 : attributes.insert(attributes.end(), pec.attributes.begin(), pec.attributes.end());
1781 30 : attributes.insert(attributes.end(), pmc.attributes.begin(), pmc.attributes.end());
1782 30 : attributes.insert(attributes.end(), auxpec.attributes.begin(), auxpec.attributes.end());
1783 30 : attributes.insert(attributes.end(), farfield.attributes.begin(),
1784 : farfield.attributes.end());
1785 30 : for (const auto &data : conductivity)
1786 : {
1787 0 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
1788 : }
1789 30 : for (const auto &data : impedance)
1790 : {
1791 0 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
1792 : }
1793 83 : for (const auto &[idx, data] : lumpedport)
1794 : {
1795 108 : for (const auto &elem : data.elements)
1796 : {
1797 55 : attributes.insert(attributes.end(), elem.attributes.begin(), elem.attributes.end());
1798 : }
1799 : }
1800 44 : for (const auto &[idx, data] : waveport)
1801 : {
1802 14 : attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
1803 : }
1804 30 : for (const auto &[idx, data] : current)
1805 : {
1806 0 : for (const auto &elem : data.elements)
1807 : {
1808 0 : attributes.insert(attributes.end(), elem.attributes.begin(), elem.attributes.end());
1809 : }
1810 : }
1811 30 : std::sort(attributes.begin(), attributes.end());
1812 30 : attributes.erase(std::unique(attributes.begin(), attributes.end()), attributes.end());
1813 : attributes.shrink_to_fit();
1814 :
1815 : // Cleanup
1816 30 : boundaries->erase("PEC");
1817 30 : boundaries->erase("PMC");
1818 30 : boundaries->erase("WavePortPEC");
1819 30 : boundaries->erase("Absorbing");
1820 30 : boundaries->erase("Conductivity");
1821 30 : boundaries->erase("Impedance");
1822 30 : boundaries->erase("LumpedPort");
1823 30 : boundaries->erase("Periodic");
1824 30 : boundaries->erase("FloquetWaveVector");
1825 30 : boundaries->erase("WavePort");
1826 30 : boundaries->erase("SurfaceCurrent");
1827 30 : boundaries->erase("Ground");
1828 30 : boundaries->erase("ZeroCharge");
1829 30 : boundaries->erase("Terminal");
1830 30 : boundaries->erase("Postprocessing");
1831 60 : MFEM_VERIFY(boundaries->empty(),
1832 : "Found an unsupported configuration file keyword under \"Boundaries\"!\n"
1833 : << boundaries->dump(2));
1834 30 : }
1835 :
1836 11 : std::vector<double> ConstructLinearRange(double start, double end, double delta)
1837 : {
1838 11 : auto n_step = GetNumSteps(start, end, delta);
1839 11 : std::vector<double> f(n_step);
1840 : std::iota(f.begin(), f.end(), 0);
1841 81 : std::for_each(f.begin(), f.end(), [=](double &x) { x = start + x * delta; });
1842 11 : return f;
1843 : }
1844 5 : std::vector<double> ConstructLinearRange(double start, double end, int n_sample)
1845 : {
1846 5 : std::vector<double> f(n_sample);
1847 44 : for (int i = 0; i < n_sample; i++)
1848 : {
1849 39 : f[i] = start + (double(i) / (n_sample - 1)) * (end - start);
1850 : }
1851 5 : return f;
1852 : }
1853 2 : std::vector<double> ConstructLogRange(double start, double end, int n_sample)
1854 : {
1855 5 : MFEM_VERIFY(start > 0 && end > 0,
1856 : "\"Type\": \"Log\" only valid for non-zero start and end!");
1857 1 : std::vector<double> f(n_sample);
1858 1 : double log_start = std::log10(start);
1859 1 : double log_end = std::log10(end);
1860 6 : for (int i = 0; i < n_sample; i++)
1861 : {
1862 5 : double log_val = log_start + (double(i) / (n_sample - 1)) * (log_end - log_start);
1863 5 : f[i] = std::pow(10.0, log_val);
1864 : }
1865 1 : return f;
1866 : }
1867 :
1868 : // Helper to find entry closest to x in vec, up to tol. If no match, returns end.
1869 4 : auto FindNearestValue(const std::vector<double> &vec, double x, double tol)
1870 : {
1871 : // Find the first element not less than x.
1872 : auto it = std::lower_bound(vec.begin(), vec.end(), x);
1873 : // Check if we found an exact match or a close enough value.
1874 4 : if (it != vec.end() && std::abs(*it - x) <= tol)
1875 : {
1876 3 : return it;
1877 : }
1878 : // If we're not at the beginning, check the previous element too.
1879 1 : if (it != vec.begin())
1880 : {
1881 : auto prev = std::prev(it);
1882 1 : if (std::abs(*prev - x) <= tol)
1883 : {
1884 0 : return prev;
1885 : }
1886 : }
1887 : // No value within tol found
1888 : return vec.end();
1889 : }
1890 :
1891 20 : void DrivenSolverData::SetUp(json &solver)
1892 : {
1893 20 : auto driven = solver.find("Driven");
1894 20 : if (driven == solver.end())
1895 : {
1896 0 : return;
1897 : }
1898 :
1899 20 : restart = driven->value("Restart", restart);
1900 20 : adaptive_tol = driven->value("AdaptiveTol", adaptive_tol);
1901 20 : adaptive_max_size = driven->value("AdaptiveMaxSamples", adaptive_max_size);
1902 20 : adaptive_memory = driven->value("AdaptiveConvergenceMemory", adaptive_memory);
1903 :
1904 20 : MFEM_VERIFY(!(restart != 1 && adaptive_tol > 0.0),
1905 : "\"Restart\" is incompatible with adaptive frequency sweep!");
1906 :
1907 : std::vector<double> save_f, prom_f; // samples to be saved to paraview and added to prom
1908 : // Backwards compatible top level interface.
1909 22 : if (driven->find("MinFreq") != driven->end() &&
1910 22 : driven->find("MaxFreq") != driven->end() && driven->find("FreqStep") != driven->end())
1911 : {
1912 2 : double min_f = driven->at("MinFreq"); // Required
1913 2 : double max_f = driven->at("MaxFreq"); // Required
1914 2 : double delta_f = driven->at("FreqStep"); // Required
1915 2 : sample_f = ConstructLinearRange(min_f, max_f, delta_f);
1916 2 : if (int save_step = driven->value("SaveStep", 0); save_step > 0)
1917 : {
1918 14 : for (std::size_t n = 0; n < sample_f.size(); n += save_step)
1919 : {
1920 12 : save_f.emplace_back(sample_f[n]);
1921 : }
1922 : }
1923 : }
1924 20 : if (auto freq_samples = driven->find("Samples"); freq_samples != driven->end())
1925 : {
1926 72 : for (auto &r : *freq_samples)
1927 : {
1928 40 : auto type = r.value("Type", r.find("Freq") != r.end() ? FrequencySampling::POINT
1929 22 : : FrequencySampling::DEFAULT);
1930 68 : auto f = [&]()
1931 : {
1932 22 : switch (type)
1933 : {
1934 15 : case FrequencySampling::LINEAR:
1935 : {
1936 15 : auto min_f = r.at("MinFreq");
1937 14 : auto max_f = r.at("MaxFreq");
1938 14 : auto delta_f = r.value("FreqStep", 0.0);
1939 14 : auto n_sample = r.value("NSample", 0);
1940 14 : MFEM_VERIFY(delta_f > 0 ^ n_sample > 0,
1941 : "Only one of \"FreqStep\" or \"NSample\" can be specified for "
1942 : "\"Type\": \"Linear\"!");
1943 14 : if (delta_f > 0)
1944 : {
1945 9 : return ConstructLinearRange(min_f, max_f, delta_f);
1946 : }
1947 5 : if (n_sample > 0)
1948 : {
1949 5 : return ConstructLinearRange(min_f, max_f, n_sample);
1950 : }
1951 : }
1952 : case FrequencySampling::LOG:
1953 : {
1954 3 : auto min_f = r.at("MinFreq");
1955 3 : auto max_f = r.at("MaxFreq");
1956 3 : auto n_sample = r.at("NSample");
1957 2 : return ConstructLogRange(min_f, max_f, n_sample);
1958 : }
1959 4 : case FrequencySampling::POINT:
1960 4 : return r.at("Freq").get<std::vector<double>>();
1961 : }
1962 0 : return std::vector<double>{};
1963 22 : }();
1964 18 : sample_f.insert(sample_f.end(), f.begin(), f.end());
1965 :
1966 18 : if (auto save_step = r.value("SaveStep", 0); save_step > 0)
1967 : {
1968 52 : for (std::size_t n = 0; n < f.size(); n += save_step)
1969 : {
1970 42 : save_f.emplace_back(f[n]);
1971 : }
1972 : }
1973 18 : if (auto prom_sample = r.value("AddToPROM", false); prom_sample)
1974 : {
1975 3 : if (adaptive_tol == 0)
1976 : {
1977 0 : MFEM_WARNING("Ignoring \"AddToPROM\" for non-adaptive simulation!");
1978 : }
1979 3 : prom_f.insert(prom_f.end(), f.begin(), f.end());
1980 : }
1981 :
1982 : // Debug
1983 : if constexpr (JSON_DEBUG)
1984 : {
1985 : std::cout << "Type: " << type << '\n';
1986 : std::cout << "Freq: " << f << '\n';
1987 : std::cout << "FreqStep: " << r.value("FreqStep", 0.0) << '\n';
1988 : std::cout << "NSample: " << r.value("NSample", 0.0) << '\n';
1989 : std::cout << "SaveStep: " << r.value("SaveStep", 0) << '\n';
1990 : std::cout << "AddToPROM: " << r.value("AddToPROM", false) << '\n';
1991 : }
1992 :
1993 : // Cleanup
1994 : r.erase("Type");
1995 : r.erase("MinFreq");
1996 : r.erase("MaxFreq");
1997 : r.erase("FreqStep");
1998 : r.erase("NSample");
1999 : r.erase("Freq");
2000 : r.erase("SaveStep");
2001 : r.erase("AddToPROM");
2002 18 : MFEM_VERIFY(r.empty(),
2003 : "Found an unsupported configuration file keyword in \"Samples\"!\n"
2004 : << r.dump(2));
2005 : }
2006 : }
2007 :
2008 : // Deduplicate all samples, and find indices of save and prom samples.
2009 : constexpr double delta_eps = 1.0e-9; // Precision in frequency comparisons (Hz)
2010 165 : auto equal_f = [=](auto x, auto y) { return std::abs(x - y) < delta_eps; };
2011 46 : auto deduplicate = [&equal_f](auto &f)
2012 : {
2013 46 : std::sort(f.begin(), f.end());
2014 46 : f.erase(std::unique(f.begin(), f.end(), equal_f), f.end());
2015 62 : };
2016 :
2017 : // Enforce explicit saves exactly match the sample frequencies.
2018 16 : deduplicate(sample_f);
2019 38 : auto explicit_save_f = driven->value("Save", std::vector<double>());
2020 19 : for (auto &f : explicit_save_f)
2021 : {
2022 4 : auto it = FindNearestValue(sample_f, f, delta_eps);
2023 8 : MFEM_VERIFY(it != sample_f.end(),
2024 : "Entry " << f << " in \"Save\" must be an explicitly sampled frequency!");
2025 3 : f = *it;
2026 : }
2027 15 : save_f.insert(save_f.end(), explicit_save_f.begin(), explicit_save_f.end());
2028 15 : deduplicate(save_f);
2029 15 : deduplicate(prom_f);
2030 :
2031 : // Given the matched ordering, and values are assigned by copying, can do a
2032 : // paired-iterator scan.
2033 53 : for (auto it_sample = sample_f.begin(), it_save = save_f.begin(); it_save != save_f.end();
2034 : ++it_save)
2035 : {
2036 92 : while (*it_sample != *it_save) // safe because save samples is a subset of samples
2037 : {
2038 : ++it_sample;
2039 54 : MFEM_VERIFY(it_sample != sample_f.end(),
2040 : "Save frequency " << *it_save << " not found in sample frequencies!");
2041 : }
2042 38 : save_indices.emplace_back(std::distance(sample_f.begin(), it_sample));
2043 : }
2044 : // PROM sampling always begins with the minimum and maximum frequencies. Exclude them from
2045 : // extra samples. Can use equality comparison given no floating point operations have been
2046 : // done.
2047 30 : prom_indices = {0, sample_f.size() - 1};
2048 15 : if (prom_f.size() > 0 && prom_f.back() == sample_f.back())
2049 : {
2050 : prom_f.pop_back();
2051 : }
2052 15 : if (prom_f.size() > 0 && prom_f.front() == sample_f.front())
2053 : {
2054 : prom_f.erase(prom_f.begin(), std::next(prom_f.begin()));
2055 : }
2056 21 : for (auto it_sample = sample_f.begin(), it_prom = prom_f.begin(); it_prom != prom_f.end();
2057 : ++it_prom)
2058 : {
2059 18 : while (*it_sample != *it_prom) // safe because prom samples is a subset of samples
2060 : {
2061 : ++it_sample;
2062 12 : MFEM_VERIFY(it_sample != sample_f.end(), "PROM sample frequency "
2063 : << *it_prom
2064 : << " not found in sample frequencies!");
2065 : }
2066 8 : prom_indices.emplace_back(std::distance(sample_f.begin(), it_sample));
2067 : }
2068 :
2069 18 : MFEM_VERIFY(!sample_f.empty(), "No sample frequency samples specified in \"Driven\"!");
2070 :
2071 : // Debug
2072 : if constexpr (JSON_DEBUG)
2073 : {
2074 : std::cout << "MinFreq: " << driven->value("MinFreq", 0.0) << '\n';
2075 : std::cout << "MaxFreq: " << driven->value("MaxFreq", 0.0) << '\n';
2076 : std::cout << "FreqStep: " << driven->value("FreqStep", 0) << '\n';
2077 : std::cout << "Samples: " << sample_f << '\n';
2078 : std::cout << "SaveSamples: " << save_f << '\n';
2079 : std::cout << "PROMSamples: " << prom_f << '\n';
2080 : std::cout << "SaveIndices: " << save_indices << '\n';
2081 : std::cout << "PromIndices: " << prom_indices << '\n';
2082 : std::cout << "Restart: " << restart << '\n';
2083 : std::cout << "AdaptiveTol: " << adaptive_tol << '\n';
2084 : std::cout << "AdaptiveMaxSamples: " << adaptive_max_size << '\n';
2085 : std::cout << "AdaptiveConvergenceMemory: " << adaptive_memory << '\n';
2086 : }
2087 :
2088 : // Cleanup
2089 14 : driven->erase("MinFreq");
2090 14 : driven->erase("MaxFreq");
2091 14 : driven->erase("FreqStep");
2092 14 : driven->erase("Samples");
2093 14 : driven->erase("Save");
2094 14 : driven->erase("SaveStep");
2095 14 : driven->erase("Restart");
2096 14 : driven->erase("AdaptiveTol");
2097 14 : driven->erase("AdaptiveMaxSamples");
2098 14 : driven->erase("AdaptiveConvergenceMemory");
2099 28 : MFEM_VERIFY(driven->empty(),
2100 : "Found an unsupported configuration file keyword under \"Driven\"!\n"
2101 : << driven->dump(2));
2102 : }
2103 :
2104 8 : void EigenSolverData::SetUp(json &solver)
2105 : {
2106 8 : auto eigenmode = solver.find("Eigenmode");
2107 8 : if (eigenmode == solver.end())
2108 : {
2109 8 : return;
2110 : }
2111 0 : MFEM_VERIFY(eigenmode->find("Target") != eigenmode->end() ||
2112 : solver.find("Driven") != solver.end(),
2113 : "Missing \"Eigenmode\" solver \"Target\" in the configuration file!");
2114 0 : target = eigenmode->value("Target", target); // Required (only for eigenmode simulations)
2115 0 : tol = eigenmode->value("Tol", tol);
2116 0 : max_it = eigenmode->value("MaxIts", max_it);
2117 0 : max_size = eigenmode->value("MaxSize", max_size);
2118 0 : n = eigenmode->value("N", n);
2119 0 : n_post = eigenmode->value("Save", n_post);
2120 0 : type = eigenmode->value("Type", type);
2121 0 : pep_linear = eigenmode->value("PEPLinear", pep_linear);
2122 0 : scale = eigenmode->value("Scaling", scale);
2123 0 : init_v0 = eigenmode->value("StartVector", init_v0);
2124 0 : init_v0_const = eigenmode->value("StartVectorConstant", init_v0_const);
2125 0 : mass_orthog = eigenmode->value("MassOrthogonal", mass_orthog);
2126 0 : nonlinear_type = eigenmode->value("NonlinearType", nonlinear_type);
2127 0 : refine_nonlinear = eigenmode->value("RefineNonlinear", refine_nonlinear);
2128 0 : linear_tol = eigenmode->value("LinearTol", linear_tol);
2129 0 : target_upper = eigenmode->value("TargetUpper", target_upper);
2130 0 : preconditioner_lag = eigenmode->value("PreconditionerLag", preconditioner_lag);
2131 0 : preconditioner_lag_tol = eigenmode->value("PreconditionerLagTol", preconditioner_lag_tol);
2132 0 : max_restart = eigenmode->value("MaxRestart", max_restart);
2133 :
2134 0 : target_upper = (target_upper < 0) ? 3 * target : target_upper; // default = 3 * target
2135 0 : MFEM_VERIFY(target > 0.0, "config[\"Eigenmode\"][\"Target\"] must be strictly positive!");
2136 0 : MFEM_VERIFY(target_upper > target, "config[\"Eigenmode\"][\"TargetUpper\"] must be "
2137 : "greater than config[\"Eigenmode\"][\"Target\"]!");
2138 0 : MFEM_VERIFY(preconditioner_lag >= 0,
2139 : "config[\"Eigenmode\"][\"PreconditionerLag\"] must be non-negative!");
2140 0 : MFEM_VERIFY(preconditioner_lag_tol >= 0,
2141 : "config[\"Eigenmode\"][\"PreconditionerLagTol\"] must be non-negative!");
2142 0 : MFEM_VERIFY(max_restart >= 0,
2143 : "config[\"Eigenmode\"][\"MaxRestart\"] must be non-negative!");
2144 :
2145 0 : MFEM_VERIFY(n > 0, "\"N\" must be greater than 0!");
2146 :
2147 : // Cleanup
2148 0 : eigenmode->erase("Target");
2149 0 : eigenmode->erase("Tol");
2150 0 : eigenmode->erase("MaxIts");
2151 0 : eigenmode->erase("MaxSize");
2152 0 : eigenmode->erase("N");
2153 0 : eigenmode->erase("Save");
2154 0 : eigenmode->erase("Type");
2155 0 : eigenmode->erase("PEPLinear");
2156 0 : eigenmode->erase("Scaling");
2157 0 : eigenmode->erase("StartVector");
2158 0 : eigenmode->erase("StartVectorConstant");
2159 0 : eigenmode->erase("MassOrthogonal");
2160 0 : eigenmode->erase("NonlinearType");
2161 0 : eigenmode->erase("RefineNonlinear");
2162 0 : eigenmode->erase("LinearTol");
2163 0 : eigenmode->erase("TargetUpper");
2164 0 : eigenmode->erase("PreconditionerLag");
2165 0 : eigenmode->erase("PreconditionerLagTol");
2166 0 : eigenmode->erase("MaxRestart");
2167 0 : MFEM_VERIFY(eigenmode->empty(),
2168 : "Found an unsupported configuration file keyword under \"Eigenmode\"!\n"
2169 : << eigenmode->dump(2));
2170 :
2171 : // Debug
2172 : if constexpr (JSON_DEBUG)
2173 : {
2174 : std::cout << "Target: " << target << '\n';
2175 : std::cout << "Tol: " << tol << '\n';
2176 : std::cout << "MaxIts: " << max_it << '\n';
2177 : std::cout << "MaxSize: " << max_size << '\n';
2178 : std::cout << "N: " << n << '\n';
2179 : std::cout << "Save: " << n_post << '\n';
2180 : std::cout << "Type: " << type << '\n';
2181 : std::cout << "PEPLinear: " << pep_linear << '\n';
2182 : std::cout << "Scaling: " << scale << '\n';
2183 : std::cout << "StartVector: " << init_v0 << '\n';
2184 : std::cout << "StartVectorConstant: " << init_v0_const << '\n';
2185 : std::cout << "MassOrthogonal: " << mass_orthog << '\n';
2186 : std::cout << "NonlinearType: " << nonlinear_type << '\n';
2187 : std::cout << "RefineNonlinear: " << refine_nonlinear << '\n';
2188 : std::cout << "LinearTol: " << linear_tol << '\n';
2189 : std::cout << "TargetUpper: " << target_upper << '\n';
2190 : std::cout << "PreconditionerLag: " << preconditioner_lag << '\n';
2191 : std::cout << "PreconditionerLagTol: " << preconditioner_lag_tol << '\n';
2192 : std::cout << "MaxRestart: " << max_restart << '\n';
2193 : }
2194 : }
2195 :
2196 8 : void ElectrostaticSolverData::SetUp(json &solver)
2197 : {
2198 8 : auto electrostatic = solver.find("Electrostatic");
2199 8 : if (electrostatic == solver.end())
2200 : {
2201 8 : return;
2202 : }
2203 0 : n_post = electrostatic->value("Save", n_post);
2204 :
2205 : // Cleanup
2206 0 : electrostatic->erase("Save");
2207 0 : MFEM_VERIFY(electrostatic->empty(),
2208 : "Found an unsupported configuration file keyword under \"Electrostatic\"!\n"
2209 : << electrostatic->dump(2));
2210 :
2211 : // Debug
2212 : if constexpr (JSON_DEBUG)
2213 : {
2214 : std::cout << "Save: " << n_post << '\n';
2215 : }
2216 : }
2217 :
2218 8 : void MagnetostaticSolverData::SetUp(json &solver)
2219 : {
2220 8 : auto magnetostatic = solver.find("Magnetostatic");
2221 8 : if (magnetostatic == solver.end())
2222 : {
2223 8 : return;
2224 : }
2225 0 : n_post = magnetostatic->value("Save", n_post);
2226 :
2227 : // Cleanup
2228 0 : magnetostatic->erase("Save");
2229 0 : MFEM_VERIFY(magnetostatic->empty(),
2230 : "Found an unsupported configuration file keyword under \"Magnetostatic\"!\n"
2231 : << magnetostatic->dump(2));
2232 :
2233 : // Debug
2234 : if constexpr (JSON_DEBUG)
2235 : {
2236 : std::cout << "Save: " << n_post << '\n';
2237 : }
2238 : }
2239 :
2240 8 : void TransientSolverData::SetUp(json &solver)
2241 : {
2242 8 : auto transient = solver.find("Transient");
2243 8 : if (transient == solver.end())
2244 : {
2245 8 : return;
2246 : }
2247 0 : MFEM_VERIFY(
2248 : transient->find("Excitation") != transient->end(),
2249 : "Missing \"Transient\" solver \"Excitation\" type in the configuration file!");
2250 0 : MFEM_VERIFY(transient->find("MaxTime") != transient->end() &&
2251 : transient->find("TimeStep") != transient->end(),
2252 : "Missing \"Transient\" solver \"MaxTime\" or \"TimeStep\" in the "
2253 : "configuration file!");
2254 0 : type = transient->value("Type", type);
2255 0 : excitation = transient->at("Excitation"); // Required
2256 0 : pulse_f = transient->value("ExcitationFreq", pulse_f);
2257 0 : pulse_tau = transient->value("ExcitationWidth", pulse_tau);
2258 0 : max_t = transient->at("MaxTime"); // Required
2259 0 : delta_t = transient->at("TimeStep"); // Required
2260 0 : delta_post = transient->value("SaveStep", delta_post);
2261 0 : order = transient->value("Order", order);
2262 0 : rel_tol = transient->value("RelTol", rel_tol);
2263 0 : abs_tol = transient->value("AbsTol", abs_tol);
2264 0 : MFEM_VERIFY(delta_t > 0, "\"TimeStep\" must be greater than 0.0!");
2265 :
2266 0 : if (type == TimeSteppingScheme::GEN_ALPHA || type == TimeSteppingScheme::RUNGE_KUTTA)
2267 : {
2268 0 : if (transient->contains("Order"))
2269 : {
2270 0 : MFEM_WARNING("GeneralizedAlpha and RungeKutta transient solvers do not use "
2271 : "config[\"Transient\"][\"Order\"]!");
2272 : }
2273 0 : if (transient->contains("RelTol") || transient->contains("AbsTol"))
2274 : {
2275 0 : MFEM_WARNING(
2276 : "GeneralizedAlpha and RungeKutta transient solvers do not use\n"
2277 : "config[\"Transient\"][\"RelTol\"] and config[\"Transient\"][\"AbsTol\"]!");
2278 : }
2279 : }
2280 : else
2281 : {
2282 0 : MFEM_VERIFY(rel_tol > 0,
2283 : "config[\"Transient\"][\"RelTol\"] must be strictly positive!");
2284 0 : MFEM_VERIFY(abs_tol > 0,
2285 : "config[\"Transient\"][\"AbsTol\"] must be strictly positive!");
2286 0 : MFEM_VERIFY(order >= 2 && order <= 5,
2287 : "config[\"Transient\"][\"Order\"] must be between 2 and 5!");
2288 : }
2289 :
2290 : // Cleanup
2291 0 : transient->erase("Type");
2292 0 : transient->erase("Excitation");
2293 0 : transient->erase("ExcitationFreq");
2294 0 : transient->erase("ExcitationWidth");
2295 0 : transient->erase("MaxTime");
2296 0 : transient->erase("TimeStep");
2297 0 : transient->erase("SaveStep");
2298 0 : transient->erase("Order");
2299 0 : transient->erase("RelTol");
2300 0 : transient->erase("AbsTol");
2301 0 : MFEM_VERIFY(transient->empty(),
2302 : "Found an unsupported configuration file keyword under \"Transient\"!\n"
2303 : << transient->dump(2));
2304 :
2305 : // Debug
2306 : if constexpr (JSON_DEBUG)
2307 : {
2308 : std::cout << "Type: " << type << '\n';
2309 : std::cout << "Excitation: " << excitation << '\n';
2310 : std::cout << "ExcitationFreq: " << pulse_f << '\n';
2311 : std::cout << "ExcitationWidth: " << pulse_tau << '\n';
2312 : std::cout << "MaxTime: " << max_t << '\n';
2313 : std::cout << "TimeStep: " << delta_t << '\n';
2314 : std::cout << "SaveStep: " << delta_post << '\n';
2315 : std::cout << "Order: " << order << '\n';
2316 : std::cout << "RelTol: " << rel_tol << '\n';
2317 : std::cout << "AbsTol: " << abs_tol << '\n';
2318 : }
2319 : }
2320 :
2321 8 : void LinearSolverData::SetUp(json &solver)
2322 : {
2323 8 : auto linear = solver.find("Linear");
2324 8 : if (linear == solver.end())
2325 : {
2326 0 : return;
2327 : }
2328 8 : type = linear->value("Type", type);
2329 8 : krylov_solver = linear->value("KSPType", krylov_solver);
2330 8 : tol = linear->value("Tol", tol);
2331 8 : max_it = linear->value("MaxIts", max_it);
2332 8 : max_size = linear->value("MaxSize", max_size);
2333 8 : initial_guess = linear->value("InitialGuess", initial_guess);
2334 :
2335 : // Options related to multigrid.
2336 8 : mg_max_levels = linear->value("MGMaxLevels", mg_max_levels);
2337 8 : mg_coarsening = linear->value("MGCoarsenType", mg_coarsening);
2338 8 : mg_use_mesh = linear->value("MGUseMesh", mg_use_mesh);
2339 8 : mg_cycle_it = linear->value("MGCycleIts", mg_cycle_it);
2340 8 : mg_smooth_aux = linear->value("MGAuxiliarySmoother", mg_smooth_aux);
2341 8 : mg_smooth_it = linear->value("MGSmoothIts", mg_smooth_it);
2342 8 : mg_smooth_order = linear->value("MGSmoothOrder", mg_smooth_order);
2343 8 : mg_smooth_sf_max = linear->value("MGSmoothEigScaleMax", mg_smooth_sf_max);
2344 8 : mg_smooth_sf_min = linear->value("MGSmoothEigScaleMin", mg_smooth_sf_min);
2345 8 : mg_smooth_cheby_4th = linear->value("MGSmoothChebyshev4th", mg_smooth_cheby_4th);
2346 :
2347 : // Preconditioner-specific options.
2348 8 : pc_mat_real = linear->value("PCMatReal", pc_mat_real);
2349 8 : pc_mat_shifted = linear->value("PCMatShifted", pc_mat_shifted);
2350 8 : complex_coarse_solve = linear->value("ComplexCoarseSolve", complex_coarse_solve);
2351 8 : pc_side = linear->value("PCSide", pc_side);
2352 8 : sym_factorization = linear->value("ColumnOrdering", sym_factorization);
2353 8 : strumpack_compression_type =
2354 8 : linear->value("STRUMPACKCompressionType", strumpack_compression_type);
2355 8 : strumpack_lr_tol = linear->value("STRUMPACKCompressionTol", strumpack_lr_tol);
2356 8 : strumpack_lossy_precision =
2357 8 : linear->value("STRUMPACKLossyPrecision", strumpack_lossy_precision);
2358 8 : strumpack_butterfly_l = linear->value("STRUMPACKButterflyLevels", strumpack_butterfly_l);
2359 8 : superlu_3d = linear->value("SuperLU3DCommunicator", superlu_3d);
2360 8 : ams_vector_interp = linear->value("AMSVectorInterpolation", ams_vector_interp);
2361 8 : ams_singular_op = linear->value("AMSSingularOperator", ams_singular_op);
2362 8 : amg_agg_coarsen = linear->value("AMGAggressiveCoarsening", amg_agg_coarsen);
2363 :
2364 : // Other linear solver options.
2365 8 : divfree_tol = linear->value("DivFreeTol", divfree_tol);
2366 8 : divfree_max_it = linear->value("DivFreeMaxIts", divfree_max_it);
2367 8 : estimator_tol = linear->value("EstimatorTol", estimator_tol);
2368 8 : estimator_max_it = linear->value("EstimatorMaxIts", estimator_max_it);
2369 8 : estimator_mg = linear->value("EstimatorMG", estimator_mg);
2370 8 : gs_orthog = linear->value("GSOrthogonalization", gs_orthog);
2371 :
2372 : // Cleanup
2373 8 : linear->erase("Type");
2374 8 : linear->erase("KSPType");
2375 8 : linear->erase("Tol");
2376 8 : linear->erase("MaxIts");
2377 8 : linear->erase("MaxSize");
2378 8 : linear->erase("InitialGuess");
2379 :
2380 8 : linear->erase("MGMaxLevels");
2381 8 : linear->erase("MGCoarsenType");
2382 8 : linear->erase("MGUseMesh");
2383 8 : linear->erase("MGCycleIts");
2384 8 : linear->erase("MGAuxiliarySmoother");
2385 8 : linear->erase("MGSmoothIts");
2386 8 : linear->erase("MGSmoothOrder");
2387 8 : linear->erase("MGSmoothEigScaleMax");
2388 8 : linear->erase("MGSmoothEigScaleMin");
2389 8 : linear->erase("MGSmoothChebyshev4th");
2390 :
2391 8 : linear->erase("PCMatReal");
2392 8 : linear->erase("PCMatShifted");
2393 8 : linear->erase("ComplexCoarseSolve");
2394 8 : linear->erase("PCSide");
2395 8 : linear->erase("ColumnOrdering");
2396 8 : linear->erase("STRUMPACKCompressionType");
2397 8 : linear->erase("STRUMPACKCompressionTol");
2398 8 : linear->erase("STRUMPACKLossyPrecision");
2399 8 : linear->erase("STRUMPACKButterflyLevels");
2400 8 : linear->erase("SuperLU3DCommunicator");
2401 8 : linear->erase("AMSVectorInterpolation");
2402 8 : linear->erase("AMSSingularOperator");
2403 8 : linear->erase("AMGAggressiveCoarsening");
2404 :
2405 8 : linear->erase("DivFreeTol");
2406 8 : linear->erase("DivFreeMaxIts");
2407 8 : linear->erase("EstimatorTol");
2408 8 : linear->erase("EstimatorMaxIts");
2409 8 : linear->erase("EstimatorMG");
2410 8 : linear->erase("GSOrthogonalization");
2411 16 : MFEM_VERIFY(linear->empty(),
2412 : "Found an unsupported configuration file keyword under \"Linear\"!\n"
2413 : << linear->dump(2));
2414 :
2415 : // Debug
2416 : if constexpr (JSON_DEBUG)
2417 : {
2418 : std::cout << "Type: " << type << '\n';
2419 : std::cout << "KSPType: " << krylov_solver << '\n';
2420 : std::cout << "Tol: " << tol << '\n';
2421 : std::cout << "MaxIts: " << max_it << '\n';
2422 : std::cout << "MaxSize: " << max_size << '\n';
2423 : std::cout << "InitialGuess: " << initial_guess << '\n';
2424 :
2425 : std::cout << "MGMaxLevels: " << mg_max_levels << '\n';
2426 : std::cout << "MGCoarsenType: " << mg_coarsening << '\n';
2427 : std::cout << "MGUseMesh: " << mg_use_mesh << '\n';
2428 : std::cout << "MGCycleIts: " << mg_cycle_it << '\n';
2429 : std::cout << "MGAuxiliarySmoother: " << mg_smooth_aux << '\n';
2430 : std::cout << "MGSmoothIts: " << mg_smooth_it << '\n';
2431 : std::cout << "MGSmoothOrder: " << mg_smooth_order << '\n';
2432 : std::cout << "MGSmoothEigScaleMax: " << mg_smooth_sf_max << '\n';
2433 : std::cout << "MGSmoothEigScaleMin: " << mg_smooth_sf_min << '\n';
2434 : std::cout << "MGSmoothChebyshev4th: " << mg_smooth_cheby_4th << '\n';
2435 :
2436 : std::cout << "PCMatReal: " << pc_mat_real << '\n';
2437 : std::cout << "PCMatShifted: " << pc_mat_shifted << '\n';
2438 : std::cout << "ComplexCoarseSolve: " << complex_coarse_solve << '\n';
2439 : std::cout << "PCSide: " << pc_side << '\n';
2440 : std::cout << "ColumnOrdering: " << sym_factorization << '\n';
2441 : std::cout << "STRUMPACKCompressionType: " << strumpack_compression_type << '\n';
2442 : std::cout << "STRUMPACKCompressionTol: " << strumpack_lr_tol << '\n';
2443 : std::cout << "STRUMPACKLossyPrecision: " << strumpack_lossy_precision << '\n';
2444 : std::cout << "STRUMPACKButterflyLevels: " << strumpack_butterfly_l << '\n';
2445 : std::cout << "SuperLU3DCommunicator: " << superlu_3d << '\n';
2446 : std::cout << "AMSVectorInterpolation: " << ams_vector_interp << '\n';
2447 : std::cout << "AMSSingularOperator: " << ams_singular_op << '\n';
2448 : std::cout << "AMGAggressiveCoarsening: " << amg_agg_coarsen << '\n';
2449 :
2450 : std::cout << "DivFreeTol: " << divfree_tol << '\n';
2451 : std::cout << "DivFreeMaxIts: " << divfree_max_it << '\n';
2452 : std::cout << "EstimatorTol: " << estimator_tol << '\n';
2453 : std::cout << "EstimatorMaxIts: " << estimator_max_it << '\n';
2454 : std::cout << "EstimatorMG: " << estimator_mg << '\n';
2455 : std::cout << "GSOrthogonalization: " << gs_orthog << '\n';
2456 : }
2457 : }
2458 :
2459 8 : void SolverData::SetUp(json &config)
2460 : {
2461 8 : auto solver = config.find("Solver");
2462 8 : if (solver == config.end())
2463 : {
2464 0 : return;
2465 : }
2466 8 : order = solver->value("Order", order);
2467 8 : pa_order_threshold = solver->value("PartialAssemblyOrder", pa_order_threshold);
2468 8 : q_order_jac = solver->value("QuadratureOrderJacobian", q_order_jac);
2469 8 : q_order_extra = solver->value("QuadratureOrderExtra", q_order_extra);
2470 8 : device = solver->value("Device", device);
2471 8 : ceed_backend = solver->value("Backend", ceed_backend);
2472 :
2473 8 : driven.SetUp(*solver);
2474 8 : eigenmode.SetUp(*solver);
2475 8 : electrostatic.SetUp(*solver);
2476 8 : magnetostatic.SetUp(*solver);
2477 8 : transient.SetUp(*solver);
2478 8 : linear.SetUp(*solver);
2479 :
2480 : // Cleanup
2481 8 : solver->erase("Order");
2482 8 : solver->erase("PartialAssemblyOrder");
2483 8 : solver->erase("QuadratureOrderJacobian");
2484 8 : solver->erase("QuadratureOrderExtra");
2485 8 : solver->erase("Device");
2486 8 : solver->erase("Backend");
2487 :
2488 8 : solver->erase("Driven");
2489 8 : solver->erase("Eigenmode");
2490 8 : solver->erase("Electrostatic");
2491 8 : solver->erase("Magnetostatic");
2492 8 : solver->erase("Transient");
2493 8 : solver->erase("Linear");
2494 16 : MFEM_VERIFY(solver->empty(),
2495 : "Found an unsupported configuration file keyword under \"Solver\"!\n"
2496 : << solver->dump(2));
2497 :
2498 : // Debug
2499 : if constexpr (JSON_DEBUG)
2500 : {
2501 : std::cout << "Order: " << order << '\n';
2502 : std::cout << "PartialAssemblyOrder: " << pa_order_threshold << '\n';
2503 : std::cout << "QuadratureOrderJacobian: " << q_order_jac << '\n';
2504 : std::cout << "QuadratureOrderExtra: " << q_order_extra << '\n';
2505 : std::cout << "Device: " << device << '\n';
2506 : std::cout << "Backend: " << ceed_backend << '\n';
2507 : }
2508 : }
2509 :
2510 11 : int GetNumSteps(double start, double end, double delta)
2511 : {
2512 11 : if (end < start)
2513 : {
2514 : return 1;
2515 : }
2516 : constexpr double delta_eps = 1.0e-9; // 9 digits of precision comparing endpoint
2517 11 : double dn = std::abs(end - start) / std::abs(delta);
2518 11 : int n_step = 1 + static_cast<int>(dn);
2519 11 : double dfinal = start + n_step * delta;
2520 11 : return n_step + ((delta < 0.0 && dfinal - end > -delta_eps * end) ||
2521 11 : (delta > 0.0 && dfinal - end < delta_eps * end));
2522 : }
2523 :
2524 : } // namespace palace::config
|
