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 "transientsolver.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "fem/errorindicator.hpp"
8 : #include "fem/mesh.hpp"
9 : #include "linalg/errorestimator.hpp"
10 : #include "linalg/vector.hpp"
11 : #include "models/lumpedportoperator.hpp"
12 : #include "models/portexcitations.hpp"
13 : #include "models/postoperator.hpp"
14 : #include "models/spaceoperator.hpp"
15 : #include "models/surfacecurrentoperator.hpp"
16 : #include "models/timeoperator.hpp"
17 : #include "utils/communication.hpp"
18 : #include "utils/excitations.hpp"
19 : #include "utils/iodata.hpp"
20 : #include "utils/timer.hpp"
21 :
22 : namespace palace
23 : {
24 :
25 : std::pair<ErrorIndicator, long long int>
26 0 : TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
27 : {
28 : // Set up the spatial discretization and time integrators for the E and B fields.
29 0 : BlockTimer bt0(Timer::CONSTRUCT);
30 0 : std::function<double(double)> J_coef = GetTimeExcitation(false);
31 0 : std::function<double(double)> dJdt_coef = GetTimeExcitation(true);
32 0 : SpaceOperator space_op(iodata, mesh);
33 0 : TimeOperator time_op(iodata, space_op, dJdt_coef);
34 :
35 0 : double delta_t = iodata.solver.transient.delta_t;
36 0 : int n_step = config::GetNumSteps(0.0, iodata.solver.transient.max_t, delta_t);
37 0 : SaveMetadata(space_op.GetNDSpaces());
38 :
39 : // Time stepping is uniform in the time domain. Index sets are for computing things like
40 : // port voltages and currents in postprocessing.
41 0 : PostOperator<ProblemType::TRANSIENT> post_op(iodata, space_op);
42 :
43 : // Transient solver only supports a single excitation, this is check in SpaceOperator.
44 0 : Mpi::Print("\nComputing transient response for:\n{}",
45 0 : space_op.GetPortExcitations().FmtLog());
46 :
47 : // Initialize structures for storing and reducing the results of error estimation.
48 : TimeDependentFluxErrorEstimator<Vector> estimator(
49 : space_op.GetMaterialOp(), space_op.GetNDSpaces(), space_op.GetRTSpaces(),
50 0 : iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
51 0 : iodata.solver.linear.estimator_mg);
52 : ErrorIndicator indicator;
53 :
54 : // Main time integration loop.
55 0 : double t = -delta_t;
56 : auto t0 = Timer::Now();
57 0 : for (int step = 0; step < n_step; step++)
58 : {
59 0 : const double ts = iodata.units.Dimensionalize<Units::ValueType::TIME>(t + delta_t);
60 0 : Mpi::Print("\nIt {:d}/{:d}: t = {:e} ns (elapsed time = {:.2e} s)\n", step, n_step - 1,
61 0 : ts, Timer::Duration(Timer::Now() - t0).count());
62 :
63 : // Single time step t -> t + dt.
64 0 : BlockTimer bt1(Timer::TS);
65 0 : if (step == 0)
66 : {
67 0 : Mpi::Print("\n");
68 0 : t += delta_t;
69 0 : time_op.Init(); // Initial conditions
70 : }
71 : else
72 : {
73 0 : time_op.Step(t, delta_t); // Advances t internally
74 : }
75 :
76 : // Postprocess for the time step.
77 0 : BlockTimer bt2(Timer::POSTPRO);
78 : const Vector &E = time_op.GetE();
79 : const Vector &B = time_op.GetB();
80 0 : Mpi::Print(" Sol. ||E|| = {:.6e}, ||B|| = {:.6e}\n",
81 0 : linalg::Norml2(space_op.GetComm(), E),
82 0 : linalg::Norml2(space_op.GetComm(), B));
83 :
84 0 : auto total_domain_energy = post_op.MeasureAndPrintAll(step, E, B, t, J_coef(t));
85 :
86 : // Calculate and record the error indicators.
87 0 : Mpi::Print(" Updating solution error estimates\n");
88 0 : estimator.AddErrorIndicator(E, B, total_domain_energy, indicator);
89 0 : }
90 : // Final postprocessing & printing.
91 0 : BlockTimer bt1(Timer::POSTPRO);
92 0 : time_op.PrintStats();
93 0 : SaveMetadata(time_op.GetLinearSolver());
94 0 : post_op.MeasureFinalize(indicator);
95 0 : return {indicator, space_op.GlobalTrueVSize()};
96 0 : }
97 :
98 0 : std::function<double(double)> TransientSolver::GetTimeExcitation(bool dot) const
99 : {
100 : using namespace excitations;
101 : using F = std::function<double(double)>;
102 0 : const config::TransientSolverData &data = iodata.solver.transient;
103 : const Excitation &type = data.excitation;
104 0 : if (type == Excitation::SINUSOIDAL || type == Excitation::MOD_GAUSSIAN)
105 : {
106 0 : MFEM_VERIFY(data.pulse_f > 0.0,
107 : "Excitation frequency is missing for transient simulation!");
108 : }
109 0 : if (type == Excitation::GAUSSIAN || type == Excitation::DIFF_GAUSSIAN ||
110 0 : type == Excitation::MOD_GAUSSIAN || type == Excitation::SMOOTH_STEP)
111 : {
112 0 : MFEM_VERIFY(data.pulse_tau > 0.0,
113 : "Excitation width is missing for transient simulation!");
114 : }
115 : const double delay = (type == Excitation::GAUSSIAN || type == Excitation::DIFF_GAUSSIAN ||
116 : type == Excitation::MOD_GAUSSIAN)
117 0 : ? 4.5 * data.pulse_tau
118 : : 0.0;
119 0 : switch (type)
120 : {
121 0 : case Excitation::SINUSOIDAL:
122 0 : if (dot)
123 : {
124 0 : return F{[=](double t) { return dpulse_sinusoidal(t, data.pulse_f, delay); }};
125 : }
126 : else
127 : {
128 0 : return F{[=](double t) { return pulse_sinusoidal(t, data.pulse_f, delay); }};
129 : }
130 : break;
131 0 : case Excitation::GAUSSIAN:
132 0 : if (dot)
133 : {
134 0 : return F{[=](double t) { return dpulse_gaussian(t, data.pulse_tau, delay); }};
135 : }
136 : else
137 : {
138 0 : return F{[=](double t) { return pulse_gaussian(t, data.pulse_tau, delay); }};
139 : }
140 : break;
141 0 : case Excitation::DIFF_GAUSSIAN:
142 0 : if (dot)
143 : {
144 0 : return F{[=](double t) { return dpulse_gaussian_diff(t, data.pulse_tau, delay); }};
145 : }
146 : else
147 : {
148 0 : return F{[=](double t) { return pulse_gaussian_diff(t, data.pulse_tau, delay); }};
149 : }
150 : break;
151 0 : case Excitation::MOD_GAUSSIAN:
152 0 : if (dot)
153 : {
154 0 : return F{[=](double t)
155 0 : { return dpulse_gaussian_mod(t, data.pulse_f, data.pulse_tau, delay); }};
156 : }
157 : else
158 : {
159 0 : return F{[=](double t)
160 0 : { return pulse_gaussian_mod(t, data.pulse_f, data.pulse_tau, delay); }};
161 : }
162 : break;
163 0 : case Excitation::RAMP_STEP:
164 0 : if (dot)
165 : {
166 0 : return F{[=](double t) { return dpulse_ramp(t, data.pulse_tau, delay); }};
167 : }
168 : else
169 : {
170 0 : return F{[=](double t) { return pulse_ramp(t, data.pulse_tau, delay); }};
171 : }
172 : break;
173 0 : case Excitation::SMOOTH_STEP:
174 0 : if (dot)
175 : {
176 0 : return F{[=](double t) { return dpulse_smootherstep(t, data.pulse_tau, delay); }};
177 : }
178 : else
179 : {
180 0 : return F{[=](double t) { return pulse_smootherstep(t, data.pulse_tau, delay); }};
181 : }
182 : break;
183 : }
184 : return F{};
185 : }
186 :
187 : } // namespace palace
|
