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 "iodata.hpp"
5 :
6 : #include <charconv>
7 : #include <fstream>
8 : #include <functional>
9 : #include <iostream>
10 : #include <regex>
11 : #include <sstream>
12 : #include <stack>
13 : #include <string>
14 : #include <string_view>
15 : #include <mfem.hpp>
16 : #include <nlohmann/json.hpp>
17 : #include "fem/bilinearform.hpp"
18 : #include "fem/integrator.hpp"
19 : #include "utils/communication.hpp"
20 : #include "utils/geodata.hpp"
21 :
22 : namespace palace
23 : {
24 :
25 26 : std::stringstream PreprocessFile(const char *filename)
26 : {
27 : // Read configuration file into memory and return as a stringstream.
28 : std::string file;
29 : {
30 26 : std::ifstream fi(filename);
31 26 : std::stringstream buf;
32 26 : if (!fi.is_open())
33 : {
34 0 : MFEM_ABORT("Unable to open configuration file \"" << filename << "\"!");
35 : }
36 26 : buf << fi.rdbuf();
37 26 : fi.close();
38 26 : file = buf.str();
39 26 : }
40 :
41 : // Strip C and C++ style comments (//, /* */) using regex. Correctly handles comments
42 : // within strings and escaped comment markers (see tinyurl.com/2s3n8dkr). An alternative
43 : // for the middle line is: R"(|(\/\*([^*]|[\r\n]|(\*+([^*\/]|[\r\n])))*\*+\/))", but this
44 : // seems to sometimes lead to issues with std::regex_replace for long files.
45 : {
46 : std::regex rgx(R"((([\"'])(?:(?=(\\?))\3.)*?\2))"
47 : R"(|(\/\*(.|[\r\n])*?\*\/))"
48 26 : R"(|(\/\/.*))");
49 26 : file = std::regex_replace(file, rgx, "$1");
50 26 : }
51 :
52 : // Also strip whitespace.
53 : {
54 : std::regex rgx(R"((([\"'])(?:(?=(\\?))\3.)*?\2))"
55 26 : R"(|(\s+))");
56 26 : file = std::regex_replace(file, rgx, "$1");
57 26 : }
58 :
59 : // Also strip erroneous trailing commas.
60 : {
61 : std::regex rgx(R"((([\"'])(?:(?=(\\?))\3.)*?\2))"
62 26 : R"(|,+(?=\s*?[\}\]]))");
63 26 : file = std::regex_replace(file, rgx, "$1");
64 26 : }
65 :
66 : // Perform integer range expansion for arrays ([a - b, c] = [a-b,c] =
67 : // [a,a+1,...,b-1,b,c]). The whole file is now one line and arrays have no spaces after
68 : // whitespace stripping.
69 26 : std::stringstream output;
70 1114 : auto RangeExpand = [](std::string_view str) -> std::string
71 : {
72 : // Handle the given string which is only numeric with possible hyphens.
73 1114 : if (str.empty())
74 : {
75 0 : return "";
76 : }
77 : int num;
78 1114 : auto [ptr, ec] = std::from_chars(str.data(), str.data() + str.length(), num);
79 1114 : MFEM_VERIFY(
80 : ec == std::errc(),
81 : "Invalid integer conversion in range expansion"
82 : << (ec == std::errc::result_out_of_range ? " (integer out of range)!" : "!"));
83 1114 : if (ptr == str.data() + str.length())
84 : {
85 1114 : return std::string(str);
86 : }
87 : // Range specified, expand the bounds.
88 : int num2;
89 0 : auto [ptr2, ec2] = std::from_chars(ptr + 1, str.data() + str.length(), num2);
90 0 : MFEM_VERIFY(
91 : ec2 == std::errc(),
92 : "Invalid integer conversion in range expansion"
93 : << (ec2 == std::errc::result_out_of_range ? " (integer out of range)!" : "!"));
94 : std::string rng;
95 0 : while (num < num2)
96 : {
97 0 : rng += std::to_string(num++) + ",";
98 : }
99 0 : rng += std::to_string(num);
100 0 : return rng;
101 : };
102 : {
103 26 : const std::string range_vals = "-0123456789,";
104 : auto start = file.begin();
105 : bool inside = false;
106 84362 : for (auto it = start; it != file.end(); ++it)
107 : {
108 84336 : if (inside)
109 : {
110 3304 : if (*it == ']')
111 : {
112 : // Apply integer range expansion (as needed) to the array, which includes only
113 : // digits, commas, and '-'. Exclude the outer square brackets.
114 1066 : std::string_view str(file.data() + (start - file.cbegin() + 1), it - start - 1);
115 : std::size_t s = 0, pos;
116 1066 : output << '[';
117 1114 : while ((pos = str.find(',', s)) != std::string::npos)
118 : {
119 96 : output << RangeExpand(str.substr(s, pos - s)) << ',';
120 48 : s = pos + 1;
121 : }
122 3198 : output << RangeExpand(str.substr(s)) << ']';
123 : start = it + 1;
124 : inside = false;
125 : }
126 2238 : else if (*it == '[')
127 : {
128 0 : output << std::string(start, it);
129 : start = it;
130 : }
131 2238 : else if (range_vals.find(*it) == std::string::npos)
132 : {
133 1998 : output << std::string(start, it);
134 : start = it;
135 : inside = false;
136 : }
137 : }
138 81032 : else if (*it == '[')
139 : {
140 5196 : output << std::string(start, it);
141 : start = it;
142 : inside = true;
143 : }
144 : }
145 78 : output << std::string(start, file.end());
146 : }
147 26 : return output;
148 0 : }
149 :
150 : using json = nlohmann::json;
151 :
152 8 : IoData::IoData(const char *filename, bool print) : units(1.0, 1.0), init(false)
153 : {
154 : // Open configuration file and preprocess: strip whitespace, comments, and expand integer
155 : // ranges.
156 8 : std::stringstream buffer = PreprocessFile(filename);
157 :
158 : // Parse the configuration file. Use a callback function to detect and throw errors for
159 : // duplicate keys.
160 : json config;
161 : std::stack<std::set<json>> parse_stack;
162 : json::parser_callback_t check_duplicate_keys =
163 1932 : [&](int, json::parse_event_t event, json &parsed)
164 : {
165 1932 : switch (event)
166 : {
167 192 : case json::parse_event_t::object_start:
168 192 : parse_stack.push(std::set<json>());
169 192 : break;
170 192 : case json::parse_event_t::object_end:
171 192 : parse_stack.pop();
172 : break;
173 650 : case json::parse_event_t::key:
174 : {
175 650 : const auto result = parse_stack.top().insert(parsed);
176 650 : if (!result.second)
177 : {
178 0 : MFEM_ABORT("Error parsing configuration file!\nDuplicate key "
179 : << parsed << " was already seen in this object!");
180 : return false;
181 : }
182 : }
183 650 : break;
184 : default:
185 : break;
186 : }
187 : return true;
188 : };
189 : try
190 : {
191 8 : config = json::parse(buffer, check_duplicate_keys);
192 : }
193 0 : catch (json::parse_error &e)
194 : {
195 0 : MFEM_ABORT("Error parsing configuration file!\n " << e.what());
196 0 : }
197 8 : if (print)
198 : {
199 0 : Mpi::Print("\n{}\n", config.dump(2));
200 : }
201 :
202 : // Set up configuration option data structures.
203 8 : problem.SetUp(config);
204 8 : model.SetUp(config);
205 8 : domains.SetUp(config);
206 8 : boundaries.SetUp(config);
207 8 : solver.SetUp(config);
208 :
209 : // Cleanup and error checking.
210 : config.erase("Problem");
211 : config.erase("Model");
212 : config.erase("Domains");
213 : config.erase("Boundaries");
214 : config.erase("Solver");
215 8 : MFEM_VERIFY(config.empty(), "Found an unsupported configuration file section!\n"
216 : << config.dump(2));
217 :
218 : // Check compatibility of configuration file and problem type.
219 8 : CheckConfiguration();
220 8 : }
221 :
222 8 : void IoData::CheckConfiguration()
223 : {
224 : // Check that the provided domain and boundary objects are all supported by the requested
225 : // problem type.
226 8 : if (problem.type == ProblemType::DRIVEN)
227 : {
228 : // No unsupported domain or boundary objects for frequency domain driven simulations.
229 : }
230 : else if (problem.type == ProblemType::EIGENMODE)
231 : {
232 : // No unsupported domain or boundary objects for frequency domain driven simulations.
233 : }
234 : else if (problem.type == ProblemType::ELECTROSTATIC)
235 : {
236 0 : if (!boundaries.farfield.empty())
237 : {
238 0 : Mpi::Warning(
239 : "Electrostatic problem type does not support absorbing boundary conditions!\n");
240 : }
241 0 : if (!boundaries.conductivity.empty())
242 : {
243 0 : Mpi::Warning("Electrostatic problem type does not support surface conductivity "
244 : "boundary conditions!\n");
245 : }
246 0 : if (!boundaries.impedance.empty())
247 : {
248 0 : Mpi::Warning("Electrostatic problem type does not support surface impedance boundary "
249 : "conditions!\n");
250 : }
251 0 : if (!boundaries.auxpec.empty() || !boundaries.waveport.empty())
252 : {
253 0 : Mpi::Warning(
254 : "Electrostatic problem type does not support wave port boundary conditions!\n");
255 : }
256 0 : if (!boundaries.current.empty())
257 : {
258 0 : Mpi::Warning(
259 : "Electrostatic problem type does not support surface current excitation!\n");
260 : }
261 : }
262 : else if (problem.type == ProblemType::MAGNETOSTATIC)
263 : {
264 0 : if (!boundaries.farfield.empty())
265 : {
266 0 : Mpi::Warning(
267 : "Magnetostatic problem type does not support absorbing boundary conditions!\n");
268 : }
269 0 : if (!boundaries.conductivity.empty())
270 : {
271 0 : Mpi::Warning("Magnetostatic problem type does not support surface conductivity "
272 : "boundary conditions!\n");
273 : }
274 0 : if (!boundaries.impedance.empty())
275 : {
276 0 : Mpi::Warning("Magnetostatic problem type does not support surface impedance boundary "
277 : "conditions!\n");
278 : }
279 0 : if (!boundaries.lumpedport.empty())
280 : {
281 0 : Mpi::Warning(
282 : "Magnetostatic problem type does not support lumped port boundary conditions!\n");
283 : }
284 0 : if (!boundaries.auxpec.empty() || !boundaries.waveport.empty())
285 : {
286 0 : Mpi::Warning(
287 : "Magnetostatic problem type does not support wave port boundary conditions!\n");
288 : }
289 0 : if (!boundaries.postpro.dielectric.empty())
290 : {
291 0 : Mpi::Warning("Magnetostatic problem type does not support boundary interface "
292 : "dielectric loss postprocessing!\n");
293 : }
294 : }
295 : else if (problem.type == ProblemType::TRANSIENT)
296 : {
297 0 : if (!boundaries.conductivity.empty())
298 : {
299 0 : Mpi::Warning("Transient problem type does not support surface conductivity boundary "
300 : "conditions!\n");
301 : }
302 0 : if (!boundaries.auxpec.empty() || !boundaries.waveport.empty())
303 : {
304 0 : Mpi::Warning(
305 : "Transient problem type does not support wave port boundary conditions!\n");
306 : }
307 0 : if (!boundaries.farfield.empty() && boundaries.farfield.order > 1)
308 : {
309 0 : Mpi::Warning("Transient problem type does not support absorbing boundary conditions "
310 : "with order > 1!\n");
311 : }
312 : }
313 :
314 : // Resolve default values in configuration file.
315 8 : if (solver.linear.type == LinearSolver::DEFAULT)
316 : {
317 8 : if (problem.type == ProblemType::ELECTROSTATIC)
318 : {
319 0 : solver.linear.type = LinearSolver::BOOMER_AMG;
320 : }
321 8 : else if (problem.type == ProblemType::MAGNETOSTATIC ||
322 : problem.type == ProblemType::TRANSIENT)
323 : {
324 0 : solver.linear.type = LinearSolver::AMS;
325 : }
326 : else
327 : {
328 : // Prefer sparse direct solver for frequency domain problems if available.
329 : #if defined(MFEM_USE_SUPERLU)
330 : solver.linear.type = LinearSolver::SUPERLU;
331 : #elif defined(MFEM_USE_STRUMPACK)
332 8 : solver.linear.type = LinearSolver::STRUMPACK;
333 : #elif defined(MFEM_USE_MUMPS)
334 : solver.linear.type = LinearSolver::MUMPS;
335 : #else
336 : solver.linear.type = LinearSolver::AMS;
337 : #endif
338 : }
339 : }
340 8 : if (solver.linear.krylov_solver == KrylovSolver::DEFAULT)
341 : {
342 : // Problems with SPD operators use CG by default, else GMRES.
343 0 : if (problem.type == ProblemType::ELECTROSTATIC ||
344 0 : problem.type == ProblemType::MAGNETOSTATIC ||
345 : problem.type == ProblemType::TRANSIENT)
346 : {
347 0 : solver.linear.krylov_solver = KrylovSolver::CG;
348 : }
349 : else
350 : {
351 0 : solver.linear.krylov_solver = KrylovSolver::GMRES;
352 : }
353 : }
354 8 : if (solver.linear.max_size < 0)
355 : {
356 8 : solver.linear.max_size = solver.linear.max_it;
357 : }
358 8 : if (solver.linear.initial_guess < 0)
359 : {
360 8 : if ((problem.type == ProblemType::DRIVEN && solver.driven.adaptive_tol <= 0.0) ||
361 0 : problem.type == ProblemType::TRANSIENT ||
362 0 : problem.type == ProblemType::ELECTROSTATIC ||
363 : problem.type == ProblemType::MAGNETOSTATIC)
364 : {
365 : // Default true only driven simulations without adaptive frequency sweep, transient
366 : // simulations, electrostatics, or magnetostatics.
367 8 : solver.linear.initial_guess = 1;
368 : }
369 : else
370 : {
371 0 : solver.linear.initial_guess = 0;
372 : }
373 : }
374 8 : if (solver.linear.pc_mat_shifted < 0)
375 : {
376 8 : if (problem.type == ProblemType::DRIVEN && solver.linear.type == LinearSolver::AMS)
377 : {
378 : // Default true only driven simulations using AMS (false for most cases).
379 0 : solver.linear.pc_mat_shifted = 1;
380 : }
381 : else
382 : {
383 8 : solver.linear.pc_mat_shifted = 0;
384 : }
385 : }
386 8 : if (solver.linear.mg_smooth_aux < 0)
387 : {
388 8 : if (problem.type == ProblemType::ELECTROSTATIC ||
389 : problem.type == ProblemType::MAGNETOSTATIC)
390 : {
391 : // Disable auxiliary space smoothing using distributive relaxation by default for
392 : // problems which don't need it.
393 0 : solver.linear.mg_smooth_aux = 0;
394 : }
395 : else
396 : {
397 8 : solver.linear.mg_smooth_aux = 1;
398 : }
399 : }
400 8 : if (solver.linear.mg_smooth_order < 0)
401 : {
402 8 : solver.linear.mg_smooth_order = std::max(2 * solver.order, 4);
403 : }
404 8 : if (solver.linear.ams_singular_op < 0)
405 : {
406 8 : solver.linear.ams_singular_op = (problem.type == ProblemType::MAGNETOSTATIC);
407 : }
408 8 : if (solver.linear.amg_agg_coarsen < 0)
409 : {
410 8 : solver.linear.amg_agg_coarsen = (problem.type == ProblemType::ELECTROSTATIC ||
411 8 : problem.type == ProblemType::MAGNETOSTATIC ||
412 : problem.type == ProblemType::TRANSIENT);
413 : }
414 :
415 : // Configure settings for quadrature rules and partial assembly.
416 8 : BilinearForm::pa_order_threshold = solver.pa_order_threshold;
417 8 : fem::DefaultIntegrationOrder::p_trial = solver.order;
418 8 : fem::DefaultIntegrationOrder::q_order_jac = solver.q_order_jac;
419 8 : fem::DefaultIntegrationOrder::q_order_extra_pk = solver.q_order_extra;
420 8 : fem::DefaultIntegrationOrder::q_order_extra_qk = solver.q_order_extra;
421 8 : }
422 :
423 : namespace
424 : {
425 :
426 : template <std::size_t N>
427 : constexpr config::SymmetricMatrixData<N> &operator/=(config::SymmetricMatrixData<N> &data,
428 : double s)
429 : {
430 104 : for (auto &x : data.s)
431 : {
432 78 : x /= s;
433 : }
434 : return data;
435 : }
436 :
437 : } // namespace
438 :
439 26 : void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh)
440 : {
441 : // Nondimensionalization of the equations is based on a given length Lc in [m], typically
442 : // the largest domain dimension. Configuration file lengths and the mesh coordinates are
443 : // provided with units of model.L0 x [m].
444 26 : MFEM_VERIFY(!init, "NondimensionalizeInputs should only be called once!");
445 26 : init = true;
446 :
447 : // Calculate the reference length and time. A user specified model.Lc is in mesh length
448 : // units.
449 26 : if (model.Lc <= 0.0)
450 : {
451 : mfem::Vector bbmin, bbmax;
452 26 : mesh::GetAxisAlignedBoundingBox(mesh, bbmin, bbmax);
453 26 : bbmax -= bbmin;
454 26 : model.Lc = *std::max_element(bbmax.begin(), bbmax.end());
455 : }
456 : // Define units now mesh length set. Note: In model field Lc is measured in units of L0.
457 26 : units = Units(model.L0, model.Lc * model.L0);
458 :
459 : // Mesh refinement parameters.
460 : auto DivideLengthScale = [Lc0 = units.GetMeshLengthRelativeScale()](double val)
461 96 : { return val / Lc0; };
462 26 : for (auto &box : model.refinement.GetBoxes())
463 : {
464 : std::transform(box.bbmin.begin(), box.bbmin.end(), box.bbmin.begin(),
465 : DivideLengthScale);
466 : std::transform(box.bbmax.begin(), box.bbmax.end(), box.bbmax.begin(),
467 : DivideLengthScale);
468 : }
469 26 : for (auto &sphere : model.refinement.GetSpheres())
470 : {
471 0 : sphere.r /= units.GetMeshLengthRelativeScale();
472 : std::transform(sphere.center.begin(), sphere.center.end(), sphere.center.begin(),
473 : DivideLengthScale);
474 : }
475 :
476 : // Materials: conductivity and London penetration depth.
477 52 : for (auto &data : domains.materials)
478 : {
479 : data.sigma /= units.GetScaleFactor<Units::ValueType::CONDUCTIVITY>();
480 26 : data.lambda_L /= units.GetMeshLengthRelativeScale();
481 : }
482 :
483 : // Probe location coordinates.
484 42 : for (auto &[idx, data] : domains.postpro.probe)
485 : {
486 : std::transform(data.center.begin(), data.center.end(), data.center.begin(),
487 : DivideLengthScale);
488 : }
489 :
490 : // Finite conductivity boundaries.
491 26 : for (auto &data : boundaries.conductivity)
492 : {
493 0 : data.sigma /= units.GetScaleFactor<Units::ValueType::CONDUCTIVITY>();
494 0 : data.h /= units.GetMeshLengthRelativeScale();
495 : }
496 :
497 : // Impedance boundaries and lumped ports.
498 26 : for (auto &data : boundaries.impedance)
499 : {
500 0 : data.Rs /= units.GetScaleFactor<Units::ValueType::IMPEDANCE>();
501 0 : data.Ls /= units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
502 0 : data.Cs /= units.GetScaleFactor<Units::ValueType::CAPACITANCE>();
503 : }
504 50 : for (auto &[idx, data] : boundaries.lumpedport)
505 : {
506 24 : data.R /= units.GetScaleFactor<Units::ValueType::IMPEDANCE>();
507 24 : data.L /= units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
508 24 : data.C /= units.GetScaleFactor<Units::ValueType::CAPACITANCE>();
509 24 : data.Rs /= units.GetScaleFactor<Units::ValueType::IMPEDANCE>();
510 24 : data.Ls /= units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
511 24 : data.Cs /= units.GetScaleFactor<Units::ValueType::CAPACITANCE>();
512 : }
513 :
514 : // Floquet periodic boundaries.
515 104 : for (auto &k : boundaries.periodic.wave_vector)
516 : {
517 78 : k *= units.GetMeshLengthRelativeScale();
518 : }
519 :
520 : // Wave port offset distance.
521 26 : for (auto &[idx, data] : boundaries.waveport)
522 : {
523 0 : data.d_offset /= units.GetMeshLengthRelativeScale();
524 : }
525 :
526 : // Center coordinates for surface flux.
527 42 : for (auto &[idx, data] : boundaries.postpro.flux)
528 : {
529 : std::transform(data.center.begin(), data.center.end(), data.center.begin(),
530 : DivideLengthScale);
531 : }
532 :
533 : // Dielectric interface thickness.
534 50 : for (auto &[idx, data] : boundaries.postpro.dielectric)
535 : {
536 24 : data.t /= units.GetMeshLengthRelativeScale();
537 : }
538 :
539 : // For eigenmode simulations:
540 26 : solver.eigenmode.target /= units.GetScaleFactor<Units::ValueType::FREQUENCY>();
541 26 : solver.eigenmode.target_upper /= units.GetScaleFactor<Units::ValueType::FREQUENCY>();
542 :
543 : // For driven simulations:
544 74 : for (auto &f : solver.driven.sample_f)
545 48 : f /= units.GetScaleFactor<Units::ValueType::FREQUENCY>();
546 :
547 : // For transient simulations:
548 26 : solver.transient.pulse_f /= units.GetScaleFactor<Units::ValueType::FREQUENCY>();
549 26 : solver.transient.pulse_tau /= units.GetScaleFactor<Units::ValueType::TIME>();
550 26 : solver.transient.max_t /= units.GetScaleFactor<Units::ValueType::TIME>();
551 26 : solver.transient.delta_t /= units.GetScaleFactor<Units::ValueType::TIME>();
552 :
553 : // Scale mesh vertices for correct nondimensionalization.
554 26 : mesh::NondimensionalizeMesh(mesh, units.GetMeshLengthRelativeScale());
555 :
556 : // Print some information.
557 26 : Mpi::Print(mesh.GetComm(),
558 : "\nCharacteristic length and time scales:\n L₀ = {:.3e} m, t₀ = {:.3e} ns\n",
559 26 : units.GetScaleFactor<Units::ValueType::LENGTH>(),
560 26 : units.GetScaleFactor<Units::ValueType::TIME>());
561 26 : }
562 :
563 : } // namespace palace
|
