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 "drivensolver.hpp"
5 :
6 : #include <complex>
7 : #include <mfem.hpp>
8 : #include "fem/errorindicator.hpp"
9 : #include "fem/mesh.hpp"
10 : #include "linalg/errorestimator.hpp"
11 : #include "linalg/floquetcorrection.hpp"
12 : #include "linalg/ksp.hpp"
13 : #include "linalg/operator.hpp"
14 : #include "linalg/vector.hpp"
15 : #include "models/lumpedportoperator.hpp"
16 : #include "models/portexcitations.hpp"
17 : #include "models/postoperator.hpp"
18 : #include "models/romoperator.hpp"
19 : #include "models/spaceoperator.hpp"
20 : #include "models/surfacecurrentoperator.hpp"
21 : #include "models/waveportoperator.hpp"
22 : #include "utils/communication.hpp"
23 : #include "utils/iodata.hpp"
24 : #include "utils/prettyprint.hpp"
25 : #include "utils/timer.hpp"
26 :
27 : namespace palace
28 : {
29 :
30 : using namespace std::complex_literals;
31 :
32 : std::pair<ErrorIndicator, long long int>
33 0 : DrivenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
34 : {
35 : // Set up the spatial discretization and frequency sweep.
36 0 : BlockTimer bt0(Timer::CONSTRUCT);
37 0 : SpaceOperator space_op(iodata, mesh);
38 : const auto &port_excitations = space_op.GetPortExcitations();
39 0 : SaveMetadata(port_excitations);
40 :
41 0 : const auto &omega_sample = iodata.solver.driven.sample_f;
42 :
43 0 : bool adaptive = (iodata.solver.driven.adaptive_tol > 0.0);
44 0 : if (adaptive && omega_sample.size() <= iodata.solver.driven.prom_indices.size())
45 : {
46 0 : Mpi::Warning("Adaptive frequency sweep requires > {} total frequency samples!\n"
47 : "Reverting to uniform sweep!\n",
48 0 : iodata.solver.driven.prom_indices.size());
49 : adaptive = false;
50 : }
51 0 : SaveMetadata(space_op.GetNDSpaces());
52 0 : Mpi::Print("\nComputing {}frequency response for:\n{}", adaptive ? "adaptive fast " : "",
53 0 : port_excitations.FmtLog());
54 :
55 0 : auto restart = iodata.solver.driven.restart;
56 0 : if (restart != 1)
57 : {
58 0 : int max_iter = omega_sample.size() * space_op.GetPortExcitations().Size();
59 0 : MFEM_VERIFY(
60 : restart - 1 < max_iter,
61 : fmt::format("\"Restart\" ({}) is greater than the number of total samples ({})!",
62 : restart, max_iter));
63 :
64 0 : Mpi::Print("\nRestarting from solve {}", iodata.solver.driven.restart);
65 : }
66 :
67 : // Main frequency sweep loop.
68 0 : return {adaptive ? SweepAdaptive(space_op) : SweepUniform(space_op),
69 0 : space_op.GlobalTrueVSize()};
70 0 : }
71 :
72 0 : ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op) const
73 : {
74 : const auto &port_excitations = space_op.GetPortExcitations();
75 0 : const auto &omega_sample = iodata.solver.driven.sample_f;
76 :
77 : // Initialize postprocessing for measurement and printers.
78 : // Initialize write directory with default path; will be changed for multi-excitations.
79 0 : PostOperator<ProblemType::DRIVEN> post_op(iodata, space_op);
80 :
81 : // Construct the system matrices defining the linear operator. PEC boundaries are handled
82 : // simply by setting diagonal entries of the system matrix for the corresponding dofs.
83 : // Because the Dirichlet BC is always homogeneous, no special elimination is required on
84 : // the RHS. Assemble the linear system for the initial frequency (so we can call
85 : // KspSolver::SetOperators). Compute everything at the first frequency step.
86 0 : auto K = space_op.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
87 0 : auto C = space_op.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
88 0 : auto M = space_op.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
89 : const auto &Curl = space_op.GetCurlMatrix();
90 :
91 : // Set up the linear solver.
92 : // The operators are constructed for each frequency step and used to initialize the ksp.
93 0 : ComplexKspSolver ksp(iodata, space_op.GetNDSpaces(), &space_op.GetH1Spaces());
94 :
95 : // Set up RHS vector for the incident field at port boundaries, and the vector for the
96 : // first frequency step.
97 0 : ComplexVector RHS(Curl.Width()), E(Curl.Width()), B(Curl.Height());
98 0 : RHS.UseDevice(true);
99 0 : E.UseDevice(true);
100 0 : B.UseDevice(true);
101 : E = 0.0;
102 : B = 0.0;
103 :
104 : // Initialize structures for storing and reducing the results of error estimation.
105 : TimeDependentFluxErrorEstimator<ComplexVector> estimator(
106 : space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
107 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
108 0 : iodata.solver.linear.estimator_mg);
109 : ErrorIndicator indicator;
110 :
111 : // If using Floquet BCs, a correction term (kp x E) needs to be added to the B field.
112 0 : std::unique_ptr<FloquetCorrSolver<ComplexVector>> floquet_corr;
113 0 : if (space_op.GetMaterialOp().HasWaveVector())
114 : {
115 0 : floquet_corr = std::make_unique<FloquetCorrSolver<ComplexVector>>(
116 : space_op.GetMaterialOp(), space_op.GetNDSpace(), space_op.GetRTSpace(),
117 0 : iodata.solver.linear.tol, iodata.solver.linear.max_it, 0);
118 : }
119 :
120 : // Main excitation and frequency loop.
121 : auto t0 = Timer::Now();
122 0 : int excitation_counter = 0;
123 : const int excitation_restart_counter =
124 0 : ((iodata.solver.driven.restart - 1) / omega_sample.size()) + 1;
125 0 : const int freq_restart_idx = (iodata.solver.driven.restart - 1) % omega_sample.size();
126 0 : for (const auto &[excitation_idx, excitation_spec] : port_excitations)
127 : {
128 0 : if (++excitation_counter < excitation_restart_counter)
129 : {
130 0 : continue;
131 : }
132 0 : if (port_excitations.Size() > 1)
133 : {
134 0 : Mpi::Print("\nSweeping excitation index {:d} ({:d}/{:d}):\n", excitation_idx,
135 0 : excitation_counter, port_excitations.Size());
136 : }
137 : // Switch paraview subfolders: one for each excitation, if nr_excitations > 1.
138 0 : post_op.InitializeParaviewDataCollection(excitation_idx);
139 :
140 : // Frequency loop.
141 : for (std::size_t omega_i =
142 0 : ((excitation_counter == excitation_restart_counter) ? freq_restart_idx : 0);
143 0 : omega_i < omega_sample.size(); omega_i++)
144 : {
145 0 : auto omega = omega_sample[omega_i];
146 : // Assemble frequency dependent matrices and initialize operators in linear
147 : // solver.
148 0 : auto A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
149 0 : auto A = space_op.GetSystemMatrix(1.0 + 0.0i, 1i * omega, -omega * omega + 0.0i,
150 0 : K.get(), C.get(), M.get(), A2.get());
151 : auto P = space_op.GetPreconditionerMatrix<ComplexOperator>(
152 0 : 1.0 + 0.0i, 1i * omega, -omega * omega + 0.0i, omega);
153 0 : ksp.SetOperators(*A, *P);
154 :
155 0 : Mpi::Print(
156 : "\nIt {:d}/{:d}: ω/2π = {:.3e} GHz (total elapsed time = {:.2e} s{})\n",
157 0 : omega_i + 1, omega_sample.size(),
158 0 : iodata.units.Dimensionalize<Units::ValueType::FREQUENCY>(omega),
159 0 : Timer::Duration(Timer::Now() - t0).count(),
160 : (port_excitations.Size() > 1)
161 0 : ? fmt::format(", solve {:d}/{:d}",
162 0 : 1 + omega_i + (excitation_counter - 1) * omega_sample.size(),
163 0 : omega_sample.size() * port_excitations.Size())
164 0 : : "");
165 :
166 : // Solve linear system.
167 0 : space_op.GetExcitationVector(excitation_idx, omega, RHS);
168 0 : Mpi::Print("\n");
169 0 : ksp.Mult(RHS, E);
170 :
171 : // Start Post-processing.
172 0 : BlockTimer bt0(Timer::POSTPRO);
173 0 : Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n",
174 0 : linalg::Norml2(space_op.GetComm(), E),
175 0 : linalg::Norml2(space_op.GetComm(), RHS));
176 :
177 : // Compute B = -1/(iω) ∇ x E on the true dofs.
178 0 : Curl.Mult(E.Real(), B.Real());
179 0 : Curl.Mult(E.Imag(), B.Imag());
180 0 : B *= -1.0 / (1i * omega);
181 0 : if (space_op.GetMaterialOp().HasWaveVector())
182 : {
183 : // Calculate B field correction for Floquet BCs.
184 : // B = -1/(iω) ∇ x E + 1/ω kp x E
185 0 : floquet_corr->AddMult(E, B, 1.0 / omega);
186 : }
187 :
188 : auto total_domain_energy =
189 0 : post_op.MeasureAndPrintAll(excitation_idx, int(omega_i), E, B, omega);
190 :
191 : // Calculate and record the error indicators.
192 0 : Mpi::Print(" Updating solution error estimates\n");
193 0 : estimator.AddErrorIndicator(E, B, total_domain_energy, indicator);
194 0 : }
195 :
196 : // Final postprocessing & printing.
197 0 : BlockTimer bt0(Timer::POSTPRO);
198 0 : SaveMetadata(ksp);
199 0 : }
200 0 : post_op.MeasureFinalize(indicator);
201 0 : return indicator;
202 0 : }
203 :
204 0 : ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
205 : {
206 : const auto &port_excitations = space_op.GetPortExcitations();
207 0 : const auto &omega_sample = iodata.solver.driven.sample_f;
208 :
209 : // Initialize postprocessing for measurement and printers.
210 : // Initialize write directory with default path; will be changed for multi-excitations.
211 0 : PostOperator<ProblemType::DRIVEN> post_op(iodata, space_op);
212 :
213 : // Configure PROM parameters if not specified.
214 0 : double offline_tol = iodata.solver.driven.adaptive_tol;
215 0 : int convergence_memory = iodata.solver.driven.adaptive_memory;
216 0 : int max_size_per_excitation = iodata.solver.driven.adaptive_max_size;
217 0 : int nprom_indices = static_cast<int>(iodata.solver.driven.prom_indices.size());
218 0 : MFEM_VERIFY(max_size_per_excitation <= 0 || max_size_per_excitation >= nprom_indices,
219 : "Adaptive frequency sweep must sample at least " << nprom_indices
220 : << " frequency points!");
221 : // Maximum size — no more than nr steps needed.
222 0 : max_size_per_excitation =
223 0 : std::min(max_size_per_excitation, static_cast<int>(omega_sample.size()));
224 :
225 : // Allocate negative curl matrix for postprocessing the B-field and vectors for the
226 : // high-dimensional field solution.
227 : const auto &Curl = space_op.GetCurlMatrix();
228 0 : ComplexVector E(Curl.Width()), Eh(Curl.Width()), B(Curl.Height());
229 0 : E.UseDevice(true);
230 0 : Eh.UseDevice(true);
231 0 : B.UseDevice(true);
232 : E = 0.0;
233 : Eh = 0.0;
234 : B = 0.0;
235 :
236 : // Initialize structures for storing and reducing the results of error estimation.
237 : TimeDependentFluxErrorEstimator<ComplexVector> estimator(
238 : space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
239 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
240 0 : iodata.solver.linear.estimator_mg);
241 : ErrorIndicator indicator;
242 :
243 : // If using Floquet BCs, a correction term (kp x E) needs to be added to the B field.
244 0 : std::unique_ptr<FloquetCorrSolver<ComplexVector>> floquet_corr;
245 0 : if (space_op.GetMaterialOp().HasWaveVector())
246 : {
247 0 : floquet_corr = std::make_unique<FloquetCorrSolver<ComplexVector>>(
248 : space_op.GetMaterialOp(), space_op.GetNDSpace(), space_op.GetRTSpace(),
249 0 : iodata.solver.linear.tol, iodata.solver.linear.max_it, 0);
250 : }
251 :
252 : // Configure the PROM operator which performs the parameter space sampling and basis
253 : // construction during the offline phase as well as the PROM solution during the online
254 : // phase.
255 : auto t0 = Timer::Now();
256 0 : const double unit_GHz = iodata.units.Dimensionalize<Units::ValueType::FREQUENCY>(1.0);
257 0 : Mpi::Print("\nBeginning PROM construction offline phase:\n"
258 : " {:d} points for frequency sweep over [{:.3e}, {:.3e}] GHz\n",
259 0 : omega_sample.size(), omega_sample.front() * unit_GHz,
260 0 : omega_sample.back() * unit_GHz);
261 0 : RomOperator prom_op(iodata, space_op, max_size_per_excitation);
262 : space_op.GetWavePortOp().SetSuppressOutput(true);
263 :
264 : // Initialize the basis with samples from the top and bottom of the frequency
265 : // range of interest. Each call for an HDM solution adds the frequency sample to P_S and
266 : // removes it from P \ P_S. Timing for the HDM construction and solve is handled inside
267 : // of the RomOperator.
268 0 : auto UpdatePROM = [&](int excitation_idx, double omega)
269 : {
270 : // Add the HDM solution to the PROM reduced basis.
271 0 : prom_op.UpdatePROM(E);
272 0 : prom_op.UpdateMRI(excitation_idx, omega, E);
273 :
274 : // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
275 : // PostOperator for energy postprocessing and error estimation.
276 0 : BlockTimer bt0(Timer::POSTPRO);
277 0 : Curl.Mult(E.Real(), B.Real());
278 0 : Curl.Mult(E.Imag(), B.Imag());
279 0 : B *= -1.0 / (1i * omega);
280 0 : if (space_op.GetMaterialOp().HasWaveVector())
281 : {
282 : // Calculate B field correction for Floquet BCs.
283 : // B = -1/(iω) ∇ x E + 1/ω kp x E
284 0 : floquet_corr->AddMult(E, B, 1.0 / omega);
285 : }
286 :
287 : // Measure domain energies for the error indicator only. Don't exchange face_nbr_data,
288 : // unless printing paraview fields.
289 0 : auto total_domain_energy = post_op.MeasureDomainFieldEnergyOnly(E, B);
290 0 : estimator.AddErrorIndicator(E, B, total_domain_energy, indicator);
291 0 : };
292 :
293 : // Loop excitations to add to PROM.
294 : //
295 : // Restart should not really be used for adaptive sweeps, but must work. Construct PROM in
296 : // the same way same regardless of restart for consistency. Don't shift excitation start.
297 0 : int excitation_counter = 0;
298 0 : for (const auto &[excitation_idx, excitation_spec] : port_excitations)
299 : {
300 0 : if (port_excitations.Size() > 1)
301 : {
302 0 : Mpi::Print("\nAdding excitation index {:d} ({:d}/{:d}):\n", excitation_idx,
303 0 : ++excitation_counter, port_excitations.Size());
304 : }
305 0 : prom_op.SetExcitationIndex(excitation_idx); // Pre-compute RHS1
306 :
307 : // Initialize PROM with explicit HDM samples, record the estimate but do not act on it.
308 : std::vector<double> max_errors;
309 0 : for (auto i : iodata.solver.driven.prom_indices)
310 : {
311 0 : auto omega = omega_sample[i];
312 0 : prom_op.SolveHDM(excitation_idx, omega, E);
313 0 : prom_op.SolvePROM(excitation_idx, omega, Eh);
314 0 : linalg::AXPY(-1.0, E, Eh);
315 0 : max_errors.push_back(linalg::Norml2(space_op.GetComm(), Eh) /
316 0 : linalg::Norml2(space_op.GetComm(), E));
317 0 : UpdatePROM(excitation_idx, omega);
318 : }
319 : // The estimates associated to the end points are assumed inaccurate.
320 0 : max_errors[0] = std::numeric_limits<double>::infinity();
321 0 : max_errors[1] = std::numeric_limits<double>::infinity();
322 : int memory = std::distance(max_errors.rbegin(),
323 0 : std::find_if(max_errors.rbegin(), max_errors.rend(),
324 0 : [=](auto x) { return x > offline_tol; }));
325 :
326 : // Greedy procedure for basis construction (offline phase). Basis is initialized with
327 : // solutions at frequency sweep endpoints and explicit sample frequencies.
328 0 : int it = static_cast<int>(max_errors.size());
329 0 : for (int it0 = it; it < max_size_per_excitation && memory < convergence_memory; it++)
330 : {
331 : // Compute the location of the maximum error in parameter domain (bounded by the
332 : // previous samples).
333 0 : double omega_star = prom_op.FindMaxError(excitation_idx)[0];
334 :
335 : // Sample HDM and add solution to basis.
336 0 : prom_op.SolveHDM(excitation_idx, omega_star, E);
337 0 : prom_op.SolvePROM(excitation_idx, omega_star, Eh);
338 0 : linalg::AXPY(-1.0, E, Eh);
339 0 : max_errors.push_back(linalg::Norml2(space_op.GetComm(), Eh) /
340 0 : linalg::Norml2(space_op.GetComm(), E));
341 0 : memory = max_errors.back() < offline_tol ? memory + 1 : 0;
342 :
343 0 : Mpi::Print("\nGreedy iteration {:d} (n = {:d}): ω* = {:.3e} GHz ({:.3e}), error = "
344 : "{:.3e}{}\n",
345 0 : it - it0 + 1, prom_op.GetReducedDimension(), omega_star * unit_GHz,
346 : omega_star, max_errors.back(),
347 : (memory == 0)
348 0 : ? ""
349 0 : : fmt::format(", memory = {:d}/{:d}", memory, convergence_memory));
350 0 : UpdatePROM(excitation_idx, omega_star);
351 : }
352 0 : Mpi::Print("\nAdaptive sampling{} {:d} frequency samples:\n"
353 : " n = {:d}, error = {:.3e}, tol = {:.3e}, memory = {:d}/{:d}\n",
354 0 : (it == max_size_per_excitation) ? " reached maximum" : " converged with", it,
355 0 : prom_op.GetReducedDimension(), max_errors.back(), offline_tol, memory,
356 : convergence_memory);
357 0 : utils::PrettyPrint(prom_op.GetSamplePoints(excitation_idx), unit_GHz,
358 : " Sampled frequencies (GHz):");
359 0 : utils::PrettyPrint(max_errors, 1.0, " Sample errors:");
360 : }
361 :
362 0 : Mpi::Print(" Total offline phase elapsed time: {:.2e} s\n",
363 0 : Timer::Duration(Timer::Now() - t0).count()); // Timing on root
364 :
365 : // XX TODO: Add output of eigenvalue estimates from the PROM system (and nonlinear EVP
366 : // in the general case with wave ports, etc.?)
367 :
368 : // Main fast frequency sweep loop (online phase).
369 0 : Mpi::Print("\nBeginning fast frequency sweep online phase\n");
370 : space_op.GetWavePortOp().SetSuppressOutput(false); // Disable output suppression
371 0 : for (const auto &[excitation_idx, excitation_spec] : port_excitations)
372 : {
373 0 : if (port_excitations.Size() > 1)
374 : {
375 0 : Mpi::Print("\nSweeping excitation index {:d} ({:d}/{:d}):\n", excitation_idx,
376 0 : excitation_counter, port_excitations.Size());
377 : }
378 : // Switch paraview subfolders: one for each excitation, if nr_excitations > 1.
379 0 : post_op.InitializeParaviewDataCollection(excitation_idx);
380 :
381 : // Frequency loop.
382 0 : for (std::size_t omega_i = 0; omega_i < omega_sample.size(); omega_i++)
383 : {
384 0 : auto omega = omega_sample[omega_i];
385 0 : Mpi::Print("\nIt {:d}/{:d}: ω/2π = {:.3e} GHz (total elapsed time = {:.2e} s)\n",
386 0 : omega_i + 1, omega_sample.size(),
387 0 : iodata.units.Dimensionalize<Units::ValueType::FREQUENCY>(omega),
388 0 : Timer::Duration(Timer::Now() - t0).count());
389 :
390 : // Assemble and solve the PROM linear system.
391 0 : prom_op.SolvePROM(excitation_idx, omega, E);
392 0 : Mpi::Print("\n");
393 :
394 : // Start Post-processing.
395 0 : BlockTimer bt0(Timer::POSTPRO);
396 0 : Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(space_op.GetComm(), E));
397 :
398 : // Compute B = -1/(iω) ∇ x E on the true dofs.
399 0 : Curl.Mult(E.Real(), B.Real());
400 0 : Curl.Mult(E.Imag(), B.Imag());
401 0 : B *= -1.0 / (1i * omega);
402 0 : if (space_op.GetMaterialOp().HasWaveVector())
403 : {
404 : // Calculate B field correction for Floquet BCs.
405 : // B = -1/(iω) ∇ x E + 1/ω kp x E
406 0 : floquet_corr->AddMult(E, B, 1.0 / omega);
407 : }
408 0 : post_op.MeasureAndPrintAll(excitation_idx, int(omega_i), E, B, omega);
409 0 : }
410 :
411 : // Final postprocessing & printing: no change to indicator since these are in PROM.
412 0 : BlockTimer bt0(Timer::POSTPRO);
413 0 : SaveMetadata(prom_op.GetLinearSolver());
414 0 : }
415 0 : post_op.MeasureFinalize(indicator);
416 0 : return indicator;
417 0 : }
418 :
419 : } // namespace palace
|
