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 "timeoperator.hpp"
5 :
6 : #include <limits>
7 : #include <vector>
8 : #include "linalg/iterative.hpp"
9 : #include "linalg/jacobi.hpp"
10 : #include "linalg/solver.hpp"
11 : #include "models/portexcitations.hpp"
12 : #include "models/spaceoperator.hpp"
13 : #include "utils/communication.hpp"
14 : #include "utils/iodata.hpp"
15 :
16 : namespace palace
17 : {
18 :
19 : namespace
20 : {
21 :
22 0 : class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
23 : {
24 : public:
25 : // MPI communicator.
26 : MPI_Comm comm;
27 :
28 : // System matrices and excitation RHS.
29 : std::unique_ptr<Operator> K, M, C;
30 : Vector NegJ;
31 :
32 : // Time dependence of current pulse for excitation: -J'(t) = -g'(t) J. This function
33 : // returns g'(t).
34 : std::function<double(double)> dJ_coef;
35 :
36 : // Internal objects for solution of linear systems during time stepping.
37 : double dt_, saved_gamma;
38 : std::unique_ptr<KspSolver> kspM, kspA;
39 : std::unique_ptr<Operator> A, B;
40 : mutable Vector RHS;
41 : int size_E, size_B;
42 :
43 : const Operator &Curl;
44 :
45 : // Bindings to SpaceOperator functions to get the system matrix and preconditioner, and
46 : // construct the linear solver.
47 : std::function<void(double dt)> ConfigureLinearSolver;
48 :
49 : public:
50 0 : TimeDependentFirstOrderOperator(const IoData &iodata, SpaceOperator &space_op,
51 : std::function<double(double)> dJ_coef, double t0,
52 : mfem::TimeDependentOperator::Type type)
53 0 : : mfem::TimeDependentOperator(2 * space_op.GetNDSpace().GetTrueVSize() +
54 : space_op.GetRTSpace().GetTrueVSize(),
55 : t0, type),
56 0 : comm(space_op.GetComm()), dJ_coef(dJ_coef),
57 0 : size_E(space_op.GetNDSpace().GetTrueVSize()),
58 0 : size_B(space_op.GetRTSpace().GetTrueVSize()), Curl(space_op.GetCurlMatrix())
59 : {
60 : // Construct the system matrices defining the linear operator. PEC boundaries are
61 : // handled simply by setting diagonal entries of the mass matrix for the corresponding
62 : // dofs. Because the Dirichlet BC is always homogeneous, no special elimination is
63 : // required on the RHS. Diagonal entries are set in M (so M is non-singular).
64 0 : K = space_op.GetStiffnessMatrix<Operator>(Operator::DIAG_ZERO);
65 0 : C = space_op.GetDampingMatrix<Operator>(Operator::DIAG_ZERO);
66 0 : M = space_op.GetMassMatrix<Operator>(Operator::DIAG_ONE);
67 :
68 : // Already asserted that only that time dependant solver only has a single excitation.
69 : auto excitation_helper = space_op.GetPortExcitations();
70 0 : auto excitation_idx = excitation_helper.excitations.begin()->first;
71 : // Set up RHS vector for the current source term: -g'(t) J, where g(t) handles the time
72 : // dependence.
73 0 : space_op.GetExcitationVector(excitation_idx, NegJ);
74 0 : RHS.SetSize(2 * size_E + size_B);
75 : RHS.UseDevice(true);
76 :
77 : // Set up linear solvers.
78 : {
79 0 : auto pcg = std::make_unique<CgSolver<Operator>>(comm, 0);
80 0 : pcg->SetInitialGuess(0);
81 0 : pcg->SetRelTol(iodata.solver.linear.tol);
82 : pcg->SetAbsTol(std::numeric_limits<double>::epsilon());
83 0 : pcg->SetMaxIter(iodata.solver.linear.max_it);
84 0 : auto jac = std::make_unique<JacobiSmoother<Operator>>(comm);
85 0 : kspM = std::make_unique<KspSolver>(std::move(pcg), std::move(jac));
86 0 : kspM->SetOperators(*M, *M);
87 : }
88 : {
89 : // For explicit schemes, recommended to just use cheaper preconditioners. Otherwise,
90 : // use AMS or a direct solver. The system matrix is formed as a sequence of matrix
91 : // vector products, and is only assembled for preconditioning.
92 0 : ConfigureLinearSolver = [this, &iodata, &space_op](double dt)
93 : {
94 : // Configure the system matrix and also the matrix (matrices) from which the
95 : // preconditioner will be constructed.
96 0 : A = space_op.GetSystemMatrix(dt * dt, dt, 1.0, K.get(), C.get(), M.get());
97 0 : B = space_op.GetPreconditionerMatrix<Operator>(dt * dt, dt, 1.0, 0.0);
98 :
99 : // Configure the solver.
100 0 : if (!kspA)
101 : {
102 0 : kspA = std::make_unique<KspSolver>(iodata, space_op.GetNDSpaces(),
103 0 : &space_op.GetH1Spaces());
104 : }
105 0 : kspA->SetOperators(*A, *B);
106 0 : };
107 : }
108 0 : }
109 :
110 : // Form the RHS for the first-order ODE system.
111 0 : void FormRHS(const Vector &u, Vector &rhs) const
112 : {
113 : Vector u1, u2, u3, rhs1, rhs2, rhs3;
114 : u1.UseDevice(true);
115 : u2.UseDevice(true);
116 : u3.UseDevice(true);
117 : rhs1.UseDevice(true);
118 : rhs2.UseDevice(true);
119 : rhs3.UseDevice(true);
120 0 : u.Read();
121 0 : u1.MakeRef(const_cast<Vector &>(u), 0, size_E);
122 0 : u2.MakeRef(const_cast<Vector &>(u), size_E, size_E);
123 0 : u3.MakeRef(const_cast<Vector &>(u), 2 * size_E, size_B);
124 0 : rhs.ReadWrite();
125 0 : rhs1.MakeRef(rhs, 0, size_E);
126 0 : rhs2.MakeRef(rhs, size_E, size_E);
127 0 : rhs3.MakeRef(rhs, 2 * size_E, size_B);
128 :
129 : // u1 = Edot, u2 = E, u3 = B
130 : // rhs1 = -(K * u2 + C * u1) - J(t)
131 : // rhs2 = u1
132 : // rhs3 = -curl u2
133 0 : K->Mult(u2, rhs1);
134 0 : if (C)
135 : {
136 0 : C->AddMult(u1, rhs1, 1.0);
137 : }
138 0 : linalg::AXPBYPCZ(-1.0, rhs1, dJ_coef(t), NegJ, 0.0, rhs1);
139 :
140 0 : rhs2 = u1;
141 :
142 0 : Curl.Mult(u2, rhs3);
143 0 : rhs3 *= -1;
144 0 : }
145 :
146 : // Solve M du = rhs
147 : // |M 0 0| |du1| = |-(K * u2 + C * u1) - J(t) |
148 : // |0 I 0| |du2| | u1 |
149 : // |0 0 I| |du3| = |-curl u2 |
150 0 : void Mult(const Vector &u, Vector &du) const override
151 : {
152 0 : if (kspM->NumTotalMult() == 0)
153 : {
154 : // Operators have already been set in constructor.
155 0 : du = 0.0;
156 : }
157 0 : FormRHS(u, RHS);
158 :
159 : Vector du1, du2, du3, RHS1, RHS2, RHS3;
160 : du1.UseDevice(true);
161 : du2.UseDevice(true);
162 : du3.UseDevice(true);
163 : RHS1.UseDevice(true);
164 : RHS2.UseDevice(true);
165 : RHS3.UseDevice(true);
166 0 : du.ReadWrite();
167 0 : du1.MakeRef(du, 0, size_E);
168 0 : du2.MakeRef(du, size_E, size_E);
169 0 : du3.MakeRef(du, 2 * size_E, size_B);
170 : RHS.ReadWrite();
171 0 : RHS1.MakeRef(RHS, 0, size_E);
172 0 : RHS2.MakeRef(RHS, size_E, size_E);
173 0 : RHS3.MakeRef(RHS, 2 * size_E, size_B);
174 :
175 0 : kspM->Mult(RHS1, du1);
176 0 : du2 = RHS2;
177 0 : du3 = RHS3;
178 0 : }
179 :
180 0 : void ImplicitSolve(double dt, const Vector &u, Vector &k) override
181 : {
182 : // Solve: M k = f(u + dt k, t)
183 : // Use block elimination to avoid solving a 3n x 3n linear system.
184 0 : if (!kspA || dt != dt_)
185 : {
186 : // Configure the linear solver, including the system matrix and also the matrix
187 : // (matrices) from which the preconditioner will be constructed.
188 0 : ConfigureLinearSolver(dt);
189 0 : dt_ = dt;
190 0 : k = 0.0;
191 : }
192 0 : Mpi::Print("\n");
193 0 : FormRHS(u, RHS);
194 :
195 : Vector k1, k2, k3, RHS1, RHS2, RHS3;
196 : k1.UseDevice(true);
197 : k2.UseDevice(true);
198 : k3.UseDevice(true);
199 : RHS1.UseDevice(true);
200 : RHS2.UseDevice(true);
201 : RHS3.UseDevice(true);
202 0 : k.ReadWrite();
203 0 : k1.MakeRef(k, 0, size_E);
204 0 : k2.MakeRef(k, size_E, size_E);
205 0 : k3.MakeRef(k, 2 * size_E, size_B);
206 : RHS.ReadWrite();
207 0 : RHS1.MakeRef(RHS, 0, size_E);
208 0 : RHS2.MakeRef(RHS, size_E, size_E);
209 0 : RHS3.MakeRef(RHS, 2 * size_E, size_B);
210 :
211 : // A k1 = RHS1 - dt K RHS2
212 0 : K->AddMult(RHS2, RHS1, -dt);
213 0 : kspA->Mult(RHS1, k1);
214 :
215 : // k2 = rhs2 + dt k1
216 0 : linalg::AXPBYPCZ(1.0, RHS2, dt, k1, 0.0, k2);
217 :
218 : // k3 = rhs3 - dt curl k2
219 0 : k3 = RHS3;
220 0 : Curl.AddMult(k2, RHS3, -dt);
221 0 : }
222 :
223 0 : void ExplicitMult(const Vector &u, Vector &v) const override { Mult(u, v); }
224 :
225 : // Setup A = M - gamma J = M + gamma C + gamma^2 K
226 0 : int SUNImplicitSetup(const Vector &y, const Vector &fy, int jok, int *jcur,
227 : double gamma) override
228 : {
229 : // Update Jacobian matrix.
230 0 : if (!kspA || gamma != saved_gamma)
231 : {
232 0 : ConfigureLinearSolver(gamma);
233 : }
234 :
235 : // Indicate Jacobian was updated.
236 0 : *jcur = 1;
237 :
238 : // Save gamma for use in solve.
239 0 : saved_gamma = gamma;
240 :
241 0 : return 0;
242 : }
243 :
244 : // Solve (Mass - dt Jacobian) x = Mass b
245 0 : int SUNImplicitSolve(const Vector &b, Vector &x, double tol) override
246 : {
247 : Vector b1, b2, b3, x1, x2, x3, RHS1;
248 : b1.UseDevice(true);
249 : b2.UseDevice(true);
250 : b3.UseDevice(true);
251 : x1.UseDevice(true);
252 : x2.UseDevice(true);
253 : x3.UseDevice(true);
254 : RHS1.UseDevice(true);
255 0 : b.Read();
256 0 : b1.MakeRef(const_cast<Vector &>(b), 0, size_E);
257 0 : b2.MakeRef(const_cast<Vector &>(b), size_E, size_E);
258 0 : b3.MakeRef(const_cast<Vector &>(b), 2 * size_E, size_B);
259 0 : x.ReadWrite();
260 0 : x1.MakeRef(x, 0, size_E);
261 0 : x2.MakeRef(x, size_E, size_E);
262 0 : x3.MakeRef(x, 2 * size_E, size_B);
263 : RHS.ReadWrite();
264 0 : RHS1.MakeRef(RHS, 0, size_E);
265 :
266 : // A x1 = M b1 - dt K b2
267 0 : M->Mult(b1, RHS1);
268 0 : K->AddMult(b2, RHS1, -saved_gamma);
269 0 : kspA->Mult(RHS1, x1);
270 :
271 : // x2 = b2 + dt x1
272 0 : linalg::AXPBYPCZ(1.0, b2, saved_gamma, x1, 0.0, x2);
273 :
274 : // x3 = b3 - dt curl x2
275 0 : x3 = b3;
276 0 : Curl.AddMult(x2, x3, -saved_gamma);
277 :
278 0 : return 0;
279 : }
280 : };
281 :
282 : } // namespace
283 :
284 0 : TimeOperator::TimeOperator(const IoData &iodata, SpaceOperator &space_op,
285 0 : std::function<double(double)> dJ_coef)
286 0 : : rel_tol(iodata.solver.transient.rel_tol), abs_tol(iodata.solver.transient.abs_tol),
287 0 : order(iodata.solver.transient.order)
288 : {
289 : auto excitation_helper = space_op.GetPortExcitations();
290 : // Should have already asserted that time dependant solver only has a single excitation.
291 0 : MFEM_VERIFY(excitation_helper.Size() == 1,
292 : fmt::format("Transient evolution currently only allows for a single "
293 : "excitation, received {}",
294 : excitation_helper.Size()));
295 :
296 : // Get sizes.
297 : int size_E = space_op.GetNDSpace().GetTrueVSize();
298 : int size_B = space_op.GetRTSpace().GetTrueVSize();
299 :
300 : // Allocate space for solution vectors.
301 0 : sol.SetSize(2 * size_E + size_B);
302 : sol.UseDevice(true);
303 : E.UseDevice(true);
304 : B.UseDevice(true);
305 : sol.ReadWrite();
306 : E.MakeRef(sol, size_E, size_E);
307 : B.MakeRef(sol, 2 * size_E, size_B);
308 :
309 : // Create ODE solver for 1st-order IVP.
310 0 : mfem::TimeDependentOperator::Type type = mfem::TimeDependentOperator::IMPLICIT;
311 0 : op = std::make_unique<TimeDependentFirstOrderOperator>(iodata, space_op, dJ_coef, 0.0,
312 : type);
313 0 : switch (iodata.solver.transient.type)
314 : {
315 0 : case TimeSteppingScheme::GEN_ALPHA:
316 : {
317 0 : constexpr double rho_inf = 1.0;
318 0 : use_mfem_integrator = true;
319 0 : ode = std::make_unique<mfem::GeneralizedAlphaSolver>(rho_inf);
320 : }
321 0 : break;
322 0 : case TimeSteppingScheme::RUNGE_KUTTA:
323 : {
324 0 : constexpr int gamma_opt = 2;
325 0 : use_mfem_integrator = true;
326 0 : ode = std::make_unique<mfem::SDIRK23Solver>(gamma_opt);
327 : }
328 0 : break;
329 0 : case TimeSteppingScheme::ARKODE:
330 : {
331 : #if defined(MFEM_USE_SUNDIALS)
332 : // SUNDIALS ARKODE solver.
333 : std::unique_ptr<mfem::ARKStepSolver> arkode;
334 0 : arkode = std::make_unique<mfem::ARKStepSolver>(space_op.GetComm(),
335 0 : mfem::ARKStepSolver::IMPLICIT);
336 : // Initialize ARKODE.
337 0 : arkode->Init(*op);
338 : // Use implicit setup/solve defined in SUNImplicit*.
339 0 : arkode->UseMFEMLinearSolver();
340 : // Implicit solve is linear and J is not time-dependent.
341 0 : ARKodeSetLinear(arkode->GetMem(), 0);
342 : // Relative and absolute tolerances.
343 0 : arkode->SetSStolerances(rel_tol, abs_tol);
344 : // Set the order of the RK scheme.
345 0 : ARKodeSetOrder(arkode->GetMem(), order);
346 : // Set the ODE solver to ARKODE.
347 : ode = std::move(arkode);
348 : #else
349 : MFEM_ABORT("Solver was not built with SUNDIALS support, please choose a "
350 : "different transient solver type!");
351 : #endif
352 : }
353 : break;
354 0 : case TimeSteppingScheme::CVODE:
355 : {
356 : #if defined(MFEM_USE_SUNDIALS)
357 : // SUNDIALS CVODE solver.
358 : std::unique_ptr<mfem::CVODESolver> cvode;
359 0 : cvode = std::make_unique<mfem::CVODESolver>(space_op.GetComm(), CV_BDF);
360 : // Initialize CVODE.
361 0 : cvode->Init(*op);
362 : // Relative and absolute tolerances for time step control.
363 0 : cvode->SetSStolerances(rel_tol, abs_tol);
364 : // Use implicit setup/solve defined in SUNImplicit*.
365 0 : cvode->UseMFEMLinearSolver();
366 : // Set the max order of the multistep scheme.
367 : // CV_BDF can go up to 5, but >= 3 is not unconditionally stable.
368 0 : cvode->SetMaxOrder(order);
369 : // Set the max number of steps allowed in one CVODE step() call.
370 0 : cvode->SetMaxNSteps(10000);
371 : // Set the ODE solver to CVODE.
372 : ode = std::move(cvode);
373 : #else
374 : MFEM_ABORT("Solver was not built with SUNDIALS support, please choose a "
375 : "different transient solver type!");
376 : #endif
377 : }
378 : break;
379 : }
380 0 : }
381 :
382 0 : const KspSolver &TimeOperator::GetLinearSolver() const
383 : {
384 0 : const auto &first_order = dynamic_cast<const TimeDependentFirstOrderOperator &>(*op);
385 0 : MFEM_VERIFY(first_order.kspA,
386 : "No linear solver for time-dependent operator has been constructed!\n");
387 0 : return *first_order.kspA;
388 : }
389 :
390 0 : void TimeOperator::Init()
391 : {
392 : // Always use zero initial conditions.
393 0 : sol = 0.0;
394 0 : if (use_mfem_integrator)
395 : {
396 0 : ode->Init(*op);
397 : }
398 0 : }
399 :
400 0 : void TimeOperator::Step(double &t, double &dt)
401 : {
402 0 : double dt_input = dt;
403 0 : ode->Step(sol, t, dt);
404 : // Ensure user-specified dt does not change.
405 0 : dt = dt_input;
406 0 : }
407 :
408 0 : void TimeOperator::PrintStats()
409 : {
410 : #if defined(MFEM_USE_SUNDIALS)
411 0 : if (mfem::ARKStepSolver *arkode = dynamic_cast<mfem::ARKStepSolver *>(ode.get()))
412 : {
413 : long int expsteps, accsteps, step_attempts, nfe_evals, nfi_evals, nlinsetups, netfails;
414 0 : ARKStepGetTimestepperStats(arkode->GetMem(), &expsteps, &accsteps, &step_attempts,
415 : &nfe_evals, &nfi_evals, &nlinsetups, &netfails);
416 :
417 : long int nniters;
418 0 : ARKodeGetNumNonlinSolvIters(arkode->GetMem(), &nniters);
419 :
420 0 : Mpi::Print("\nARKODE time-stepper statistics\n");
421 0 : Mpi::Print(" Stability-limited steps: {:d}\n", expsteps);
422 0 : Mpi::Print(" Accuracy-limited steps: {:d}\n", accsteps);
423 0 : Mpi::Print(" Calls to explicit RHS function: {:d}\n", nfe_evals);
424 0 : Mpi::Print(" Calls to implicit RHS function: {:d}\n", nfi_evals);
425 0 : Mpi::Print(" Calls to linear solver setup function: {:d}\n", nlinsetups);
426 0 : Mpi::Print(" Calls to linear solver solve function: {:d}\n", nniters);
427 0 : Mpi::Print(" Number of error test failures: {:d}\n", netfails);
428 : }
429 0 : else if (mfem::CVODESolver *cvode = dynamic_cast<mfem::CVODESolver *>(ode.get()))
430 : {
431 : long int nsteps, nfevals, nlinsetups, netfails;
432 : int qlast, qcur;
433 : double hinused, hlast, hcur, tcur;
434 :
435 : // Get integrator stats.
436 0 : CVodeGetIntegratorStats(cvode->GetMem(), &nsteps, &nfevals, &nlinsetups, &netfails,
437 : &qlast, &qcur, &hinused, &hlast, &hcur, &tcur);
438 0 : Mpi::Print("\n CVODE time-stepper statistics\n");
439 0 : Mpi::Print(" Number of steps: {:d}\n", nsteps);
440 0 : Mpi::Print(" Calls to RHS function: {:d}\n", nfevals);
441 0 : Mpi::Print(" Calls to linear solver setup function: {:d}\n", nlinsetups);
442 0 : Mpi::Print(" Number of error test failures: {:d}\n", netfails);
443 0 : Mpi::Print("\n");
444 : }
445 : #endif
446 0 : }
447 :
448 : } // namespace palace
|
