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 "eigensolver.hpp"
5 :
6 : #include <complex>
7 : #include <mfem.hpp>
8 : #include "fem/errorindicator.hpp"
9 : #include "fem/mesh.hpp"
10 : #include "linalg/arpack.hpp"
11 : #include "linalg/divfree.hpp"
12 : #include "linalg/errorestimator.hpp"
13 : #include "linalg/floquetcorrection.hpp"
14 : #include "linalg/ksp.hpp"
15 : #include "linalg/nleps.hpp"
16 : #include "linalg/operator.hpp"
17 : #include "linalg/rap.hpp"
18 : #include "linalg/slepc.hpp"
19 : #include "linalg/vector.hpp"
20 : #include "models/lumpedportoperator.hpp"
21 : #include "models/postoperator.hpp"
22 : #include "models/spaceoperator.hpp"
23 : #include "utils/communication.hpp"
24 : #include "utils/iodata.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 : EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
34 : {
35 : // Construct and extract the system matrices defining the eigenvalue problem. The diagonal
36 : // values for the mass matrix PEC dof shift the Dirichlet eigenvalues out of the
37 : // computational range. The damping matrix may be nullptr.
38 0 : BlockTimer bt0(Timer::CONSTRUCT);
39 0 : SpaceOperator space_op(iodata, mesh);
40 0 : auto K = space_op.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
41 0 : auto C = space_op.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
42 0 : auto M = space_op.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
43 :
44 : // Check if there are nonlinear terms and, if so, setup interpolation operator.
45 : auto funcA2 = [&space_op](double omega) -> std::unique_ptr<ComplexOperator>
46 0 : { return space_op.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO); };
47 : auto funcP = [&space_op](std::complex<double> a0, std::complex<double> a1,
48 : std::complex<double> a2,
49 : double omega) -> std::unique_ptr<ComplexOperator>
50 0 : { return space_op.GetPreconditionerMatrix<ComplexOperator>(a0, a1, a2, omega); };
51 0 : const double target = iodata.solver.eigenmode.target;
52 : auto A2 = funcA2(target);
53 : bool has_A2 = (A2 != nullptr);
54 :
55 : // Extend K, C, M operators with interpolated A2 operator.
56 : // K' = K + A2_0, C' = C + A2_1, M' = M + A2_2
57 : std::unique_ptr<ComplexOperator> Kp, Cp, Mp;
58 : std::unique_ptr<Interpolation> interp_op;
59 : std::unique_ptr<ComplexOperator> A2_0, A2_1, A2_2;
60 0 : NonlinearEigenSolver nonlinear_type = iodata.solver.eigenmode.nonlinear_type;
61 0 : if (has_A2 && nonlinear_type == NonlinearEigenSolver::HYBRID)
62 : {
63 0 : const double target_max = iodata.solver.eigenmode.target_upper;
64 0 : interp_op = std::make_unique<NewtonInterpolationOperator>(funcA2, A2->Width());
65 0 : interp_op->Interpolate(1i * target, 1i * target_max);
66 0 : A2_0 = interp_op->GetInterpolationOperator(0);
67 0 : A2_1 = interp_op->GetInterpolationOperator(1);
68 0 : A2_2 = interp_op->GetInterpolationOperator(2);
69 0 : Kp = BuildParSumOperator({1.0 + 0i, 1.0 + 0i}, {K.get(), A2_0.get()});
70 0 : Cp = BuildParSumOperator({1.0 + 0i, 1.0 + 0i}, {C.get(), A2_1.get()});
71 0 : Mp = BuildParSumOperator({1.0 + 0i, 1.0 + 0i}, {M.get(), A2_2.get()});
72 : }
73 :
74 : const auto &Curl = space_op.GetCurlMatrix();
75 0 : SaveMetadata(space_op.GetNDSpaces());
76 :
77 : // Configure objects for postprocessing.
78 0 : PostOperator<ProblemType::EIGENMODE> post_op(iodata, space_op);
79 0 : ComplexVector E(Curl.Width()), B(Curl.Height());
80 0 : E.UseDevice(true);
81 0 : B.UseDevice(true);
82 :
83 : // Define and configure the eigensolver to solve the eigenvalue problem:
84 : // (K + λ C + λ² M) u = 0 or K u = -λ² M u
85 : // with λ = iω. In general, the system matrices are complex and symmetric.
86 0 : std::unique_ptr<EigenvalueSolver> eigen;
87 0 : EigenSolverBackend type = iodata.solver.eigenmode.type;
88 :
89 : #if defined(PALACE_WITH_ARPACK) && defined(PALACE_WITH_SLEPC)
90 : if (type == EigenSolverBackend::DEFAULT)
91 : {
92 : type = EigenSolverBackend::SLEPC;
93 : }
94 : #elif defined(PALACE_WITH_ARPACK)
95 : if (type == EigenSolverBackend::SLEPC)
96 : {
97 : Mpi::Warning("SLEPc eigensolver not available, using ARPACK!\n");
98 : }
99 : type = EigenSolverBackend::ARPACK;
100 : if (nonlinear_type == NonlinearEigenSolver::SLP)
101 : {
102 : Mpi::Warning("SLP nonlinear eigensolver not available, using Hybrid!\n");
103 : }
104 : nonlinear_type = NonlinearEigenSolver::HYBRID;
105 : #elif defined(PALACE_WITH_SLEPC)
106 0 : if (type == EigenSolverBackend::ARPACK)
107 : {
108 0 : Mpi::Warning("ARPACK eigensolver not available, using SLEPc!\n");
109 : }
110 : type = EigenSolverBackend::SLEPC;
111 : #else
112 : #error "Eigenmode solver requires building with ARPACK or SLEPc!"
113 : #endif
114 : if (type == EigenSolverBackend::ARPACK)
115 : {
116 : #if defined(PALACE_WITH_ARPACK)
117 : Mpi::Print("\nConfiguring ARPACK eigenvalue solver:\n");
118 : if (C || has_A2)
119 : {
120 : eigen = std::make_unique<arpack::ArpackPEPSolver>(space_op.GetComm(),
121 : iodata.problem.verbose);
122 : }
123 : else
124 : {
125 : eigen = std::make_unique<arpack::ArpackEPSSolver>(space_op.GetComm(),
126 : iodata.problem.verbose);
127 : }
128 : #endif
129 : }
130 : else // EigenSolverBackend::SLEPC
131 : {
132 : #if defined(PALACE_WITH_SLEPC)
133 0 : Mpi::Print("\nConfiguring SLEPc eigenvalue solver:\n");
134 : std::unique_ptr<slepc::SlepcEigenvalueSolver> slepc;
135 0 : if (nonlinear_type == NonlinearEigenSolver::SLP)
136 : {
137 0 : slepc = std::make_unique<slepc::SlepcNEPSolver>(space_op.GetComm(),
138 0 : iodata.problem.verbose);
139 0 : slepc->SetType(slepc::SlepcEigenvalueSolver::Type::SLP);
140 0 : slepc->SetProblemType(slepc::SlepcEigenvalueSolver::ProblemType::GENERAL);
141 : }
142 : else
143 : {
144 0 : if (C || has_A2)
145 : {
146 0 : if (!iodata.solver.eigenmode.pep_linear)
147 : {
148 0 : slepc = std::make_unique<slepc::SlepcPEPSolver>(space_op.GetComm(),
149 0 : iodata.problem.verbose);
150 0 : slepc->SetType(slepc::SlepcEigenvalueSolver::Type::TOAR);
151 : }
152 : else
153 : {
154 0 : slepc = std::make_unique<slepc::SlepcPEPLinearSolver>(space_op.GetComm(),
155 0 : iodata.problem.verbose);
156 0 : slepc->SetType(slepc::SlepcEigenvalueSolver::Type::KRYLOVSCHUR);
157 : }
158 : }
159 : else
160 : {
161 0 : slepc = std::make_unique<slepc::SlepcEPSSolver>(space_op.GetComm(),
162 0 : iodata.problem.verbose);
163 0 : slepc->SetType(slepc::SlepcEigenvalueSolver::Type::KRYLOVSCHUR);
164 : }
165 0 : slepc->SetProblemType(slepc::SlepcEigenvalueSolver::ProblemType::GEN_NON_HERMITIAN);
166 : }
167 0 : slepc->SetOrthogonalization(iodata.solver.linear.gs_orthog == Orthogonalization::MGS,
168 0 : iodata.solver.linear.gs_orthog == Orthogonalization::CGS2);
169 : eigen = std::move(slepc);
170 : #endif
171 : }
172 0 : EigenvalueSolver::ScaleType scale = iodata.solver.eigenmode.scale
173 0 : ? EigenvalueSolver::ScaleType::NORM_2
174 : : EigenvalueSolver::ScaleType::NONE;
175 0 : if (nonlinear_type == NonlinearEigenSolver::SLP)
176 : {
177 0 : eigen->SetOperators(*K, *C, *M, EigenvalueSolver::ScaleType::NONE);
178 0 : eigen->SetExtraSystemMatrix(funcA2);
179 0 : eigen->SetPreconditionerUpdate(funcP);
180 : }
181 : else
182 : {
183 0 : if (has_A2)
184 : {
185 0 : eigen->SetOperators(*Kp, *Cp, *Mp, scale);
186 : }
187 0 : else if (C)
188 : {
189 0 : eigen->SetOperators(*K, *C, *M, scale);
190 : }
191 : else
192 : {
193 0 : eigen->SetOperators(*K, *M, scale);
194 : }
195 : }
196 0 : eigen->SetNumModes(iodata.solver.eigenmode.n, iodata.solver.eigenmode.max_size);
197 : const double tol = (has_A2 && nonlinear_type == NonlinearEigenSolver::HYBRID)
198 0 : ? iodata.solver.eigenmode.linear_tol
199 0 : : iodata.solver.eigenmode.tol;
200 0 : eigen->SetTol(tol);
201 0 : eigen->SetMaxIter(iodata.solver.eigenmode.max_it);
202 0 : Mpi::Print(" Scaling γ = {:.3e}, δ = {:.3e}\n", eigen->GetScalingGamma(),
203 0 : eigen->GetScalingDelta());
204 :
205 : // If desired, use an M-inner product for orthogonalizing the eigenvalue subspace. The
206 : // constructed matrix just references the real SPD part of the mass matrix (no copy is
207 : // performed). Boundary conditions don't need to be eliminated here.
208 : std::unique_ptr<Operator> KM;
209 0 : if (iodata.solver.eigenmode.mass_orthog)
210 : {
211 0 : Mpi::Print(" Basis uses M-inner product\n");
212 0 : KM = space_op.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get());
213 0 : eigen->SetBMat(*KM);
214 :
215 : // Mpi::Print(" Basis uses (K + M)-inner product\n");
216 : // KM = space_op.GetInnerProductMatrix(1.0, 1.0, K.get(), M.get());
217 : // eigen->SetBMat(*KM);
218 : }
219 :
220 : // Construct a divergence-free projector so the eigenvalue solve is performed in the space
221 : // orthogonal to the zero eigenvalues of the stiffness matrix.
222 0 : std::unique_ptr<DivFreeSolver<ComplexVector>> divfree;
223 0 : if (iodata.solver.linear.divfree_max_it > 0 &&
224 0 : !space_op.GetMaterialOp().HasWaveVector() &&
225 : !space_op.GetMaterialOp().HasLondonDepth())
226 : {
227 0 : Mpi::Print(" Configuring divergence-free projection\n");
228 0 : constexpr int divfree_verbose = 0;
229 0 : divfree = std::make_unique<DivFreeSolver<ComplexVector>>(
230 : space_op.GetMaterialOp(), space_op.GetNDSpace(), space_op.GetH1Spaces(),
231 0 : space_op.GetAuxBdrTDofLists(), iodata.solver.linear.divfree_tol,
232 0 : iodata.solver.linear.divfree_max_it, divfree_verbose);
233 0 : eigen->SetDivFreeProjector(*divfree);
234 : }
235 :
236 : // If using Floquet BCs, a correction term (kp x E) needs to be added to the B field.
237 0 : std::unique_ptr<FloquetCorrSolver<ComplexVector>> floquet_corr;
238 0 : if (space_op.GetMaterialOp().HasWaveVector())
239 : {
240 0 : floquet_corr = std::make_unique<FloquetCorrSolver<ComplexVector>>(
241 : space_op.GetMaterialOp(), space_op.GetNDSpace(), space_op.GetRTSpace(),
242 0 : iodata.solver.linear.tol, iodata.solver.linear.max_it, 0);
243 : }
244 :
245 : // Set up the initial space for the eigenvalue solve. Satisfies boundary conditions and is
246 : // projected appropriately.
247 0 : if (iodata.solver.eigenmode.init_v0)
248 : {
249 0 : ComplexVector v0;
250 0 : if (iodata.solver.eigenmode.init_v0_const)
251 : {
252 0 : Mpi::Print(" Using constant starting vector\n");
253 0 : space_op.GetConstantInitialVector(v0);
254 : }
255 : else
256 : {
257 0 : Mpi::Print(" Using random starting vector\n");
258 0 : space_op.GetRandomInitialVector(v0);
259 : }
260 0 : if (divfree)
261 : {
262 0 : divfree->Mult(v0);
263 : }
264 0 : eigen->SetInitialSpace(v0); // Copies the vector
265 :
266 : // Debug
267 : // const auto &Grad = space_op.GetGradMatrix();
268 : // ComplexVector r0(Grad->Width());
269 : // r0.UseDevice(true);
270 : // Grad.MultTranspose(v0.Real(), r0.Real());
271 : // Grad.MultTranspose(v0.Imag(), r0.Imag());
272 : // r0.Print();
273 : }
274 :
275 : // Configure the shift-and-invert strategy is employed to solve for the eigenvalues
276 : // closest to the specified target, σ.
277 : {
278 : const double f_target =
279 0 : iodata.units.Dimensionalize<Units::ValueType::FREQUENCY>(target);
280 0 : Mpi::Print(" Shift-and-invert σ = {:.3e} GHz ({:.3e})\n", f_target, target);
281 : }
282 0 : if (C || has_A2 || nonlinear_type == NonlinearEigenSolver::SLP)
283 : {
284 : // Search for eigenvalues closest to λ = iσ.
285 0 : eigen->SetShiftInvert(1i * target);
286 : if (type == EigenSolverBackend::ARPACK)
287 : {
288 : // ARPACK searches based on eigenvalues of the transformed problem. The eigenvalue
289 : // 1 / (λ - σ) will be a large-magnitude negative imaginary number for an eigenvalue
290 : // λ with frequency close to but not below the target σ.
291 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::SMALLEST_IMAGINARY);
292 : }
293 0 : else if (nonlinear_type == NonlinearEigenSolver::SLP)
294 : {
295 0 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::TARGET_MAGNITUDE);
296 : }
297 : else
298 : {
299 0 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::TARGET_IMAGINARY);
300 : }
301 : }
302 : else
303 : {
304 : // Linear EVP has eigenvalues μ = -λ² = ω². Search for eigenvalues closest to μ = σ².
305 0 : eigen->SetShiftInvert(target * target);
306 : if (type == EigenSolverBackend::ARPACK)
307 : {
308 : // ARPACK searches based on eigenvalues of the transformed problem. 1 / (μ - σ²)
309 : // will be a large-magnitude positive real number for an eigenvalue μ with frequency
310 : // close to but below the target σ².
311 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::LARGEST_REAL);
312 : }
313 : else
314 : {
315 0 : eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::TARGET_REAL);
316 : }
317 : }
318 :
319 : // Set up the linear solver required for solving systems involving the shifted operator
320 : // (K - σ² M) or P(iσ) = (K + iσ C - σ² M) during the eigenvalue solve. The
321 : // preconditioner for complex linear systems is constructed from a real approximation
322 : // to the complex system matrix.
323 0 : auto A = space_op.GetSystemMatrix(1.0 + 0.0i, 1i * target, -target * target + 0.0i,
324 0 : K.get(), C.get(), M.get(), A2.get());
325 : auto P = space_op.GetPreconditionerMatrix<ComplexOperator>(
326 0 : 1.0 + 0.0i, 1i * target, -target * target + 0.0i, target);
327 : auto ksp = std::make_unique<ComplexKspSolver>(iodata, space_op.GetNDSpaces(),
328 0 : &space_op.GetH1Spaces());
329 0 : ksp->SetOperators(*A, *P);
330 0 : eigen->SetLinearSolver(*ksp);
331 :
332 : // Initialize structures for storing and reducing the results of error estimation.
333 : TimeDependentFluxErrorEstimator<ComplexVector> estimator(
334 : space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
335 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
336 0 : iodata.solver.linear.estimator_mg);
337 : ErrorIndicator indicator;
338 :
339 : // Eigenvalue problem solve.
340 0 : BlockTimer bt1(Timer::EPS);
341 0 : Mpi::Print("\n");
342 0 : int num_conv = eigen->Solve();
343 : {
344 0 : std::complex<double> lambda = (num_conv > 0) ? eigen->GetEigenvalue(0) : 0.0;
345 0 : Mpi::Print(" Found {:d} converged eigenvalue{}{}\n", num_conv,
346 0 : (num_conv > 1) ? "s" : "",
347 : (num_conv > 0)
348 0 : ? fmt::format(" (first = {:.3e}{:+.3e}i)", lambda.real(), lambda.imag())
349 0 : : "");
350 : }
351 :
352 0 : if (has_A2 && nonlinear_type == NonlinearEigenSolver::HYBRID)
353 : {
354 0 : Mpi::Print("\n Refining eigenvalues with Quasi-Newton solver\n");
355 0 : auto qn = std::make_unique<QuasiNewtonSolver>(space_op.GetComm(), std::move(eigen),
356 0 : num_conv, iodata.problem.verbose,
357 0 : iodata.solver.eigenmode.refine_nonlinear);
358 0 : qn->SetTol(iodata.solver.eigenmode.tol);
359 0 : qn->SetMaxIter(iodata.solver.eigenmode.max_it);
360 0 : if (C)
361 : {
362 0 : qn->SetOperators(*K, *C, *M, EigenvalueSolver::ScaleType::NONE);
363 : }
364 : else
365 : {
366 0 : qn->SetOperators(*K, *M, EigenvalueSolver::ScaleType::NONE);
367 : }
368 0 : qn->SetExtraSystemMatrix(funcA2);
369 0 : qn->SetPreconditionerUpdate(funcP);
370 0 : qn->SetNumModes(iodata.solver.eigenmode.n, iodata.solver.eigenmode.max_size);
371 0 : qn->SetPreconditionerLag(iodata.solver.eigenmode.preconditioner_lag,
372 0 : iodata.solver.eigenmode.preconditioner_lag_tol);
373 0 : qn->SetMaxRestart(iodata.solver.eigenmode.max_restart);
374 0 : qn->SetLinearSolver(*ksp);
375 0 : qn->SetShiftInvert(1i * target);
376 : eigen = std::move(qn);
377 :
378 : // Suppress wave port output during nonlinear eigensolver iterations.
379 : space_op.GetWavePortOp().SetSuppressOutput(true);
380 0 : num_conv = eigen->Solve();
381 : space_op.GetWavePortOp().SetSuppressOutput(false);
382 : }
383 :
384 0 : BlockTimer bt2(Timer::POSTPRO);
385 0 : SaveMetadata(*ksp);
386 :
387 : // Calculate and record the error indicators, and postprocess the results.
388 0 : Mpi::Print("\nComputing solution error estimates and performing postprocessing\n");
389 0 : if (!KM)
390 : {
391 : // Normalize the finalized eigenvectors with respect to mass matrix (unit electric field
392 : // energy) even if they are not computed to be orthogonal with respect to it.
393 0 : KM = space_op.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get());
394 0 : eigen->SetBMat(*KM);
395 0 : eigen->RescaleEigenvectors(num_conv);
396 : }
397 0 : Mpi::Print("\n");
398 :
399 0 : for (int i = 0; i < num_conv; i++)
400 : {
401 : // Get the eigenvalue and relative error.
402 0 : std::complex<double> omega = eigen->GetEigenvalue(i);
403 0 : double error_bkwd = eigen->GetError(i, EigenvalueSolver::ErrorType::BACKWARD);
404 0 : double error_abs = eigen->GetError(i, EigenvalueSolver::ErrorType::ABSOLUTE);
405 0 : if (!C && !has_A2)
406 : {
407 : // Linear EVP has eigenvalue μ = -λ² = ω².
408 : omega = std::sqrt(omega);
409 : }
410 : else
411 : {
412 : // Quadratic EVP solves for eigenvalue λ = iω.
413 : omega /= 1i;
414 : }
415 :
416 : // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
417 : // PostOperator for all postprocessing operations.
418 0 : eigen->GetEigenvector(i, E);
419 :
420 0 : linalg::NormalizePhase(space_op.GetComm(), E);
421 :
422 0 : Curl.Mult(E.Real(), B.Real());
423 0 : Curl.Mult(E.Imag(), B.Imag());
424 0 : B *= -1.0 / (1i * omega);
425 0 : if (space_op.GetMaterialOp().HasWaveVector())
426 : {
427 : // Calculate B field correction for Floquet BCs.
428 : // B = -1/(iω) ∇ x E + 1/ω kp x E.
429 0 : floquet_corr->AddMult(E, B, 1.0 / omega);
430 : }
431 :
432 : auto total_domain_energy =
433 0 : post_op.MeasureAndPrintAll(i, E, B, omega, error_abs, error_bkwd, num_conv);
434 :
435 : // Calculate and record the error indicators.
436 0 : if (i < iodata.solver.eigenmode.n)
437 : {
438 0 : estimator.AddErrorIndicator(E, B, total_domain_energy, indicator);
439 : }
440 :
441 : // Final write: Different condition than end of loop (i = num_conv - 1).
442 0 : if (i == iodata.solver.eigenmode.n - 1)
443 : {
444 0 : post_op.MeasureFinalize(indicator);
445 : }
446 : }
447 0 : MFEM_VERIFY(num_conv >= iodata.solver.eigenmode.n, "Eigenmode solve only found "
448 : << num_conv << " modes when "
449 : << iodata.solver.eigenmode.n
450 : << " were requested!");
451 0 : return {indicator, space_op.GlobalTrueVSize()};
452 0 : }
453 :
454 : } // namespace palace
|
