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 "electrostaticsolver.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/laplaceoperator.hpp"
13 : #include "models/postoperator.hpp"
14 : #include "utils/communication.hpp"
15 : #include "utils/iodata.hpp"
16 : #include "utils/timer.hpp"
17 :
18 : namespace palace
19 : {
20 :
21 : std::pair<ErrorIndicator, long long int>
22 0 : ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
23 : {
24 : // Construct the system matrix defining the linear operator. Dirichlet boundaries are
25 : // handled eliminating the rows and columns of the system matrix for the corresponding
26 : // dofs. The eliminated matrix is stored in order to construct the RHS vector for nonzero
27 : // prescribed BC values.
28 0 : BlockTimer bt0(Timer::CONSTRUCT);
29 0 : LaplaceOperator laplace_op(iodata, mesh);
30 0 : auto K = laplace_op.GetStiffnessMatrix();
31 : const auto &Grad = laplace_op.GetGradMatrix();
32 0 : SaveMetadata(laplace_op.GetH1Spaces());
33 :
34 : // Set up the linear solver.
35 0 : KspSolver ksp(iodata, laplace_op.GetH1Spaces());
36 0 : ksp.SetOperators(*K, *K);
37 :
38 : // Terminal indices are the set of boundaries over which to compute the capacitance
39 : // matrix. Terminal boundaries are aliases for ports.
40 0 : PostOperator<ProblemType::ELECTROSTATIC> post_op(iodata, laplace_op);
41 0 : int n_step = static_cast<int>(laplace_op.GetSources().size());
42 0 : MFEM_VERIFY(n_step > 0, "No terminal boundaries specified for electrostatic simulation!");
43 :
44 : // Right-hand side term and solution vector storage.
45 : Vector RHS(Grad.Width()), E(Grad.Height());
46 0 : std::vector<Vector> V(n_step);
47 :
48 : // Initialize structures for storing and reducing the results of error estimation.
49 : GradFluxErrorEstimator estimator(
50 : laplace_op.GetMaterialOp(), laplace_op.GetNDSpace(), laplace_op.GetRTSpaces(),
51 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
52 0 : iodata.solver.linear.estimator_mg);
53 : ErrorIndicator indicator;
54 :
55 : // Main loop over terminal boundaries.
56 0 : Mpi::Print("\nComputing electrostatic fields for {:d} terminal {}\n", n_step,
57 0 : (n_step > 1) ? "boundaries" : "boundary");
58 : int step = 0;
59 : auto t0 = Timer::Now();
60 0 : for (const auto &[idx, data] : laplace_op.GetSources())
61 : {
62 0 : Mpi::Print("\nIt {:d}/{:d}: Index = {:d} (elapsed time = {:.2e} s)\n", step + 1, n_step,
63 0 : idx, Timer::Duration(Timer::Now() - t0).count());
64 :
65 : // Form and solve the linear system for a prescribed nonzero voltage on the specified
66 : // terminal.
67 0 : Mpi::Print("\n");
68 0 : laplace_op.GetExcitationVector(idx, *K, V[step], RHS);
69 0 : ksp.Mult(RHS, V[step]);
70 :
71 : // Start Post-processing.
72 0 : BlockTimer bt2(Timer::POSTPRO);
73 0 : Mpi::Print(" Sol. ||V|| = {:.6e} (||RHS|| = {:.6e})\n",
74 0 : linalg::Norml2(laplace_op.GetComm(), V[step]),
75 0 : linalg::Norml2(laplace_op.GetComm(), RHS));
76 :
77 : // Compute E = -∇V on the true dofs.
78 0 : E = 0.0;
79 0 : Grad.AddMult(V[step], E, -1.0);
80 :
81 : // Measurement and printing.
82 0 : auto total_domain_energy = post_op.MeasureAndPrintAll(step, V[step], E, idx);
83 :
84 : // Calculate and record the error indicators.
85 0 : Mpi::Print(" Updating solution error estimates\n");
86 0 : estimator.AddErrorIndicator(E, total_domain_energy, indicator);
87 :
88 : // Next terminal.
89 : step++;
90 0 : }
91 :
92 : // Postprocess the capacitance matrix from the computed field solutions.
93 0 : BlockTimer bt1(Timer::POSTPRO);
94 0 : SaveMetadata(ksp);
95 0 : PostprocessTerminals(post_op, laplace_op.GetSources(), V);
96 0 : post_op.MeasureFinalize(indicator);
97 0 : return {indicator, laplace_op.GlobalTrueVSize()};
98 0 : }
99 :
100 0 : void ElectrostaticSolver::PostprocessTerminals(
101 : PostOperator<ProblemType::ELECTROSTATIC> &post_op,
102 : const std::map<int, mfem::Array<int>> &terminal_sources,
103 : const std::vector<Vector> &V) const
104 : {
105 : // Postprocess the Maxwell capacitance matrix. See p. 97 of the COMSOL AC/DC Module manual
106 : // for the associated formulas based on the electric field energy based on a unit voltage
107 : // excitation for each terminal. Alternatively, we could compute the resulting terminal
108 : // charges from the prescribed voltage to get C directly as:
109 : // Q_i = ∫ ρ dV = ∫ ∇ ⋅ (ε E) dV = ∫ (ε E) ⋅ n dS
110 : // and C_ij = Q_i/V_j. The energy formulation avoids having to locally integrate E = -∇V.
111 0 : mfem::DenseMatrix C(V.size()), Cm(V.size());
112 0 : for (int i = 0; i < C.Height(); i++)
113 : {
114 : // Diagonal: Cᵢᵢ = 2 Uₑ(Vᵢ) / Vᵢ² = (Vᵢᵀ K Vᵢ) / Vᵢ² (with ∀i, Vᵢ = 1)
115 : auto &V_gf = post_op.GetVGridFunction().Real();
116 0 : auto &D_gf = post_op.GetDomainPostOp().D;
117 0 : V_gf.SetFromTrueDofs(V[i]);
118 0 : post_op.GetDomainPostOp().M_elec->Mult(V_gf, D_gf);
119 0 : C(i, i) = Cm(i, i) = linalg::Dot<Vector>(post_op.GetComm(), V_gf, D_gf);
120 :
121 : // Off-diagonals: Cᵢⱼ = Uₑ(Vᵢ + Vⱼ) / (Vᵢ Vⱼ) - 1/2 (Vᵢ/Vⱼ Cᵢᵢ + Vⱼ/Vᵢ Cⱼⱼ)
122 : // = (Vⱼᵀ K Vᵢ) / (Vᵢ Vⱼ)
123 0 : for (int j = i + 1; j < C.Width(); j++)
124 : {
125 0 : V_gf.SetFromTrueDofs(V[j]);
126 0 : C(i, j) = linalg::Dot<Vector>(post_op.GetComm(), V_gf, D_gf);
127 0 : Cm(i, j) = -C(i, j);
128 0 : Cm(i, i) -= Cm(i, j);
129 : }
130 :
131 : // Copy lower triangle from already computed upper triangle.
132 0 : for (int j = 0; j < i; j++)
133 : {
134 0 : C(i, j) = C(j, i);
135 0 : Cm(i, j) = Cm(j, i);
136 0 : Cm(i, i) -= Cm(i, j);
137 : }
138 : }
139 0 : mfem::DenseMatrix Cinv(C);
140 0 : Cinv.Invert(); // In-place, uses LAPACK (when available) and should be cheap
141 :
142 : // Only root writes to disk (every process has full matrices).
143 0 : if (!root)
144 : {
145 : return;
146 : }
147 : using VT = Units::ValueType;
148 : using fmt::format;
149 :
150 : // Write capacitance matrix data.
151 0 : auto PrintMatrix = [&terminal_sources, this](const std::string &file,
152 : const std::string &name,
153 : const std::string &unit,
154 : const mfem::DenseMatrix &mat, double scale)
155 : {
156 0 : TableWithCSVFile output(post_dir / file);
157 0 : output.table.insert(Column("i", "i", 0, 0, 2, ""));
158 : int j = 0;
159 0 : for (const auto &[idx2, data2] : terminal_sources)
160 : {
161 0 : output.table.insert(format("i2{}", idx2), format("{}[i][{}] {}", name, idx2, unit));
162 : // Use the fact that iterator over i and j is the same span.
163 0 : output.table["i"] << idx2;
164 :
165 0 : auto &col = output.table[format("i2{}", idx2)];
166 0 : for (std::size_t i = 0; i < terminal_sources.size(); i++)
167 : {
168 0 : col << mat(i, j) * scale;
169 : }
170 0 : j++;
171 : }
172 0 : output.WriteFullTableTrunc();
173 0 : };
174 0 : const double F = iodata.units.Dimensionalize<VT::CAPACITANCE>(1.0);
175 0 : PrintMatrix("terminal-C.csv", "C", "(F)", C, F);
176 0 : PrintMatrix("terminal-Cinv.csv", "C⁻¹", "(1/F)", Cinv, 1.0 / F);
177 0 : PrintMatrix("terminal-Cm.csv", "C_m", "(F)", Cm, F);
178 :
179 : // Also write out a file with terminal voltage excitations.
180 : {
181 0 : TableWithCSVFile terminal_V(post_dir / "terminal-V.csv");
182 0 : terminal_V.table.insert(Column("i", "i", 0, 0, 2, ""));
183 0 : terminal_V.table.insert("Vinc", "V_inc[i] (V)");
184 0 : for (const auto &[idx, data] : terminal_sources)
185 : {
186 0 : terminal_V.table["i"] << double(idx);
187 0 : terminal_V.table["Vinc"] << iodata.units.Dimensionalize<VT::VOLTAGE>(1.0);
188 : }
189 0 : terminal_V.WriteFullTableTrunc();
190 : }
191 0 : }
192 :
193 : } // namespace palace
|
