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 "magnetostaticsolver.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "fem/errorindicator.hpp"
8 : #include "fem/mesh.hpp"
9 : #include "linalg/errorestimator.hpp"
10 : #include "linalg/ksp.hpp"
11 : #include "linalg/operator.hpp"
12 : #include "models/curlcurloperator.hpp"
13 : #include "models/postoperator.hpp"
14 : #include "models/surfacecurrentoperator.hpp"
15 : #include "utils/communication.hpp"
16 : #include "utils/iodata.hpp"
17 : #include "utils/timer.hpp"
18 :
19 : namespace palace
20 : {
21 :
22 : std::pair<ErrorIndicator, long long int>
23 0 : MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
24 : {
25 : // Construct the system matrix defining the linear operator. Dirichlet boundaries are
26 : // handled eliminating the rows and columns of the system matrix for the corresponding
27 : // dofs.
28 0 : BlockTimer bt0(Timer::CONSTRUCT);
29 0 : CurlCurlOperator curlcurl_op(iodata, mesh);
30 0 : auto K = curlcurl_op.GetStiffnessMatrix();
31 : const auto &Curl = curlcurl_op.GetCurlMatrix();
32 0 : SaveMetadata(curlcurl_op.GetNDSpaces());
33 :
34 : // Set up the linear solver.
35 0 : KspSolver ksp(iodata, curlcurl_op.GetNDSpaces(), &curlcurl_op.GetH1Spaces());
36 0 : ksp.SetOperators(*K, *K);
37 :
38 : // Terminal indices are the set of boundaries over which to compute the inductance matrix.
39 0 : PostOperator<ProblemType::MAGNETOSTATIC> post_op(iodata, curlcurl_op);
40 0 : int n_step = static_cast<int>(curlcurl_op.GetSurfaceCurrentOp().Size());
41 0 : MFEM_VERIFY(n_step > 0,
42 : "No surface current boundaries specified for magnetostatic simulation!");
43 :
44 : // Source term and solution vector storage.
45 : Vector RHS(Curl.Width()), B(Curl.Height());
46 0 : std::vector<Vector> A(n_step);
47 0 : std::vector<double> I_inc(n_step);
48 :
49 : // Initialize structures for storing and reducing the results of error estimation.
50 : CurlFluxErrorEstimator estimator(
51 : curlcurl_op.GetMaterialOp(), curlcurl_op.GetRTSpace(), curlcurl_op.GetNDSpaces(),
52 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
53 0 : iodata.solver.linear.estimator_mg);
54 : ErrorIndicator indicator;
55 :
56 : // Main loop over current source boundaries.
57 0 : Mpi::Print("\nComputing magnetostatic fields for {:d} source {}\n", n_step,
58 0 : (n_step > 1) ? "boundaries" : "boundary");
59 : int step = 0;
60 : auto t0 = Timer::Now();
61 0 : for (const auto &[idx, data] : curlcurl_op.GetSurfaceCurrentOp())
62 : {
63 0 : Mpi::Print("\nIt {:d}/{:d}: Index = {:d} (elapsed time = {:.2e} s)\n", step + 1, n_step,
64 0 : idx, Timer::Duration(Timer::Now() - t0).count());
65 :
66 : // Form and solve the linear system for a prescribed current on the specified source.
67 0 : Mpi::Print("\n");
68 0 : A[step].SetSize(RHS.Size());
69 0 : A[step].UseDevice(true);
70 0 : A[step] = 0.0;
71 0 : curlcurl_op.GetExcitationVector(idx, RHS);
72 0 : ksp.Mult(RHS, A[step]);
73 :
74 : // Start Post-processing.
75 0 : BlockTimer bt2(Timer::POSTPRO);
76 0 : Mpi::Print(" Sol. ||A|| = {:.6e} (||RHS|| = {:.6e})\n",
77 0 : linalg::Norml2(curlcurl_op.GetComm(), A[step]),
78 0 : linalg::Norml2(curlcurl_op.GetComm(), RHS));
79 :
80 : // Compute B = ∇ x A on the true dofs.
81 0 : Curl.Mult(A[step], B);
82 :
83 : // Save excitation current for inductance matrix calculation.
84 0 : I_inc[step] = data.GetExcitationCurrent();
85 :
86 : // Measurement and printing.
87 0 : auto total_domain_energy = post_op.MeasureAndPrintAll(step, A[step], B, idx);
88 :
89 : // Calculate and record the error indicators.
90 0 : Mpi::Print(" Updating solution error estimates\n");
91 0 : estimator.AddErrorIndicator(B, total_domain_energy, indicator);
92 :
93 : // Next source.
94 : step++;
95 0 : }
96 :
97 : // Postprocess the inductance matrix from the computed field solutions.
98 0 : BlockTimer bt1(Timer::POSTPRO);
99 0 : SaveMetadata(ksp);
100 0 : PostprocessTerminals(post_op, curlcurl_op.GetSurfaceCurrentOp(), A, I_inc);
101 0 : post_op.MeasureFinalize(indicator);
102 0 : return {indicator, curlcurl_op.GlobalTrueVSize()};
103 0 : }
104 :
105 0 : void MagnetostaticSolver::PostprocessTerminals(
106 : PostOperator<ProblemType::MAGNETOSTATIC> &post_op,
107 : const SurfaceCurrentOperator &surf_j_op, const std::vector<Vector> &A,
108 : const std::vector<double> &I_inc) const
109 : {
110 : // Postprocess the Maxwell inductance matrix. See p. 97 of the COMSOL AC/DC Module manual
111 : // for the associated formulas based on the magnetic field energy based on a current
112 : // excitation for each port. Alternatively, we could compute the resulting loop fluxes to
113 : // get M directly as:
114 : // Φ_i = ∫ B ⋅ n_j dS
115 : // and M_ij = Φ_i/I_j. The energy formulation avoids having to locally integrate B =
116 : // ∇ x A.
117 0 : mfem::DenseMatrix M(A.size()), Mm(A.size());
118 0 : for (int i = 0; i < M.Height(); i++)
119 : {
120 : // Diagonal: Mᵢᵢ = 2 Uₘ(Aᵢ) / Iᵢ² = (Aᵢᵀ K Aᵢ) / Iᵢ²
121 : auto &A_gf = post_op.GetAGridFunction().Real();
122 0 : auto &H_gf = post_op.GetDomainPostOp().H;
123 0 : A_gf.SetFromTrueDofs(A[i]);
124 0 : post_op.GetDomainPostOp().M_mag->Mult(A_gf, H_gf);
125 0 : M(i, i) = Mm(i, i) =
126 0 : linalg::Dot<Vector>(post_op.GetComm(), A_gf, H_gf) / (I_inc[i] * I_inc[i]);
127 :
128 : // Off-diagonals: Mᵢⱼ = Uₘ(Aᵢ + Aⱼ) / (Iᵢ Iⱼ) - 1/2 (Iᵢ/Iⱼ Mᵢᵢ + Iⱼ/Iᵢ Mⱼⱼ)
129 : // = (Aⱼᵀ K Aᵢ) / (Iᵢ Iⱼ)
130 0 : for (int j = i + 1; j < M.Width(); j++)
131 : {
132 0 : A_gf.SetFromTrueDofs(A[j]);
133 0 : M(i, j) = linalg::Dot<Vector>(post_op.GetComm(), A_gf, H_gf) / (I_inc[i] * I_inc[j]);
134 0 : Mm(i, j) = -M(i, j);
135 0 : Mm(i, i) -= Mm(i, j);
136 : }
137 :
138 : // Copy lower triangle from already computed upper triangle.
139 0 : for (int j = 0; j < i; j++)
140 : {
141 0 : M(i, j) = M(j, i);
142 0 : Mm(i, j) = Mm(j, i);
143 0 : Mm(i, i) -= Mm(i, j);
144 : }
145 : }
146 0 : mfem::DenseMatrix Minv(M);
147 0 : Minv.Invert(); // In-place, uses LAPACK (when available) and should be cheap
148 :
149 : // Only root writes to disk (every process has full matrices).
150 0 : if (!root)
151 : {
152 : return;
153 : }
154 : using fmt::format;
155 :
156 : // Write inductance matrix data.
157 0 : auto PrintMatrix = [&surf_j_op, this](const std::string &file, const std::string &name,
158 : const std::string &unit,
159 : const mfem::DenseMatrix &mat, double scale)
160 : {
161 0 : TableWithCSVFile output(post_dir / file);
162 0 : output.table.insert(Column("i", "i", 0, 0, 2, ""));
163 : int j = 0;
164 0 : for (const auto &[idx2, data2] : surf_j_op)
165 : {
166 0 : output.table.insert(format("i2{}", idx2), format("{}[i][{}] {}", name, idx2, unit));
167 : // Use the fact that iterator over i and j is the same span.
168 0 : output.table["i"] << idx2;
169 :
170 0 : auto &col = output.table[format("i2{}", idx2)];
171 0 : for (std::size_t i = 0; i < surf_j_op.Size(); i++)
172 : {
173 0 : col << mat(i, j) * scale;
174 : }
175 0 : j++;
176 : }
177 0 : output.WriteFullTableTrunc();
178 0 : };
179 0 : const double H = iodata.units.GetScaleFactor<Units::ValueType::INDUCTANCE>();
180 0 : PrintMatrix("terminal-M.csv", "M", "(H)", M, H);
181 0 : PrintMatrix("terminal-Minv.csv", "M⁻¹", "(1/H)", Minv, 1.0 / H);
182 0 : PrintMatrix("terminal-Mm.csv", "M_m", "(H)", Mm, H);
183 :
184 : // Also write out a file with source current excitations.
185 : {
186 0 : TableWithCSVFile terminal_I(post_dir / "terminal-I.csv");
187 0 : terminal_I.table.insert(Column("i", "i", 0, 0, 2, ""));
188 0 : terminal_I.table.insert("Iinc", "I_inc[i] (A)");
189 : int i = 0;
190 0 : for (const auto &[idx, data] : surf_j_op)
191 : {
192 0 : terminal_I.table["i"] << double(idx);
193 0 : terminal_I.table["Iinc"] << iodata.units.Dimensionalize<Units::ValueType::CURRENT>(
194 0 : I_inc[i]);
195 0 : i++;
196 : }
197 0 : terminal_I.WriteFullTableTrunc();
198 : }
199 0 : }
200 :
201 : } // namespace palace
|
