Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/euler_solver.h"
20 :
21 : #include "libmesh/adjoint_refinement_estimator.h"
22 : #include "libmesh/diff_system.h"
23 : #include "libmesh/error_vector.h"
24 : #include "libmesh/int_range.h"
25 : #include "libmesh/numeric_vector.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 :
32 163 : EulerSolver::EulerSolver (sys_type & s)
33 163 : : FirstOrderUnsteadySolver(s), theta(1.)
34 : {
35 163 : }
36 :
37 :
38 :
39 118 : EulerSolver::~EulerSolver () = default;
40 :
41 :
42 :
43 648 : Real EulerSolver::error_order() const
44 : {
45 648 : if (theta == 0.5)
46 0 : return 2.;
47 216 : return 1.;
48 : }
49 :
50 :
51 :
52 :
53 5288520 : bool EulerSolver::element_residual (bool request_jacobian,
54 : DiffContext & context)
55 : {
56 5288520 : bool compute_second_order_eqns = this->_system.get_physics()->have_second_order_vars();
57 :
58 6992120 : return this->_general_residual(request_jacobian,
59 : context,
60 : &DifferentiablePhysics::mass_residual,
61 : &DifferentiablePhysics::damping_residual,
62 : &DifferentiablePhysics::_eulerian_time_deriv,
63 : &DifferentiablePhysics::element_constraint,
64 : &DiffContext::elem_reinit,
65 5288520 : compute_second_order_eqns);
66 : }
67 :
68 :
69 :
70 512304 : bool EulerSolver::side_residual (bool request_jacobian,
71 : DiffContext & context)
72 : {
73 669312 : return this->_general_residual(request_jacobian,
74 : context,
75 : &DifferentiablePhysics::side_mass_residual,
76 : &DifferentiablePhysics::side_damping_residual,
77 : &DifferentiablePhysics::side_time_derivative,
78 : &DifferentiablePhysics::side_constraint,
79 : &DiffContext::elem_side_reinit,
80 512304 : false);
81 : }
82 :
83 :
84 :
85 0 : bool EulerSolver::nonlocal_residual (bool request_jacobian,
86 : DiffContext & context)
87 : {
88 0 : bool compute_second_order_eqns = this->_system.have_second_order_scalar_vars();
89 :
90 0 : return this->_general_residual(request_jacobian,
91 : context,
92 : &DifferentiablePhysics::nonlocal_mass_residual,
93 : &DifferentiablePhysics::nonlocal_damping_residual,
94 : &DifferentiablePhysics::nonlocal_time_derivative,
95 : &DifferentiablePhysics::nonlocal_constraint,
96 : &DiffContext::nonlocal_reinit,
97 0 : compute_second_order_eqns);
98 : }
99 :
100 :
101 :
102 5800824 : bool EulerSolver::_general_residual (bool request_jacobian,
103 : DiffContext & context,
104 : ResFuncType mass,
105 : ResFuncType damping,
106 : ResFuncType time_deriv,
107 : ResFuncType constraint,
108 : ReinitFuncType reinit_func,
109 : bool compute_second_order_eqns)
110 : {
111 1860608 : unsigned int n_dofs = context.get_elem_solution().size();
112 :
113 : // We might need to save the old jacobian in case one of our physics
114 : // terms later is unable to update it analytically.
115 9522040 : DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
116 5800824 : if (request_jacobian)
117 1101980 : old_elem_jacobian.swap(context.get_elem_jacobian());
118 :
119 : // Local nonlinear solution at old timestep
120 5800824 : DenseVector<Number> old_elem_solution(n_dofs);
121 50101096 : for (unsigned int i=0; i != n_dofs; ++i)
122 44300272 : old_elem_solution(i) =
123 57305236 : old_nonlinear_solution(context.get_dof_indices()[i]);
124 :
125 : // Local time derivative of solution
126 1860608 : context.get_elem_solution_rate() = context.get_elem_solution();
127 3940216 : context.get_elem_solution_rate() -= old_elem_solution;
128 5800824 : context.elem_solution_rate_derivative = 1 / _system.deltat;
129 1860608 : context.get_elem_solution_rate() *=
130 2079608 : context.elem_solution_rate_derivative;
131 :
132 : // If we are asked to compute residuals for second order variables,
133 : // we also populate the acceleration part so the user can use that.
134 5800824 : if (compute_second_order_eqns)
135 106216 : this->prepare_accel(context);
136 :
137 : // Local nonlinear solution at time t_theta
138 1860608 : DenseVector<Number> theta_solution(context.get_elem_solution());
139 5800824 : theta_solution *= theta;
140 5800824 : theta_solution.add(1. - theta, old_elem_solution);
141 :
142 5800824 : context.elem_solution_derivative = theta;
143 5800824 : context.fixed_solution_derivative = theta;
144 :
145 : // If a fixed solution is requested, we'll use theta_solution
146 5800824 : if (_system.use_fixed_solution)
147 0 : context.get_elem_fixed_solution() = theta_solution;
148 :
149 : // Move theta_->elem_, elem_->theta_
150 1860608 : context.get_elem_solution().swap(theta_solution);
151 :
152 : // Move the mesh into place first if necessary, set t = t_{\theta}
153 5800824 : (context.*reinit_func)(theta);
154 :
155 : // Get the time derivative at t_theta
156 : bool jacobian_computed =
157 9741040 : (_system.get_physics()->*time_deriv)(request_jacobian, context);
158 :
159 9741040 : jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
160 : jacobian_computed;
161 :
162 : // If we have second-order variables, we need to get damping terms
163 : // and the velocity equations
164 5800824 : if (compute_second_order_eqns)
165 : {
166 182380 : jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
167 : jacobian_computed;
168 :
169 106216 : jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
170 : jacobian_computed;
171 : }
172 :
173 : // Restore the elem position if necessary, set t = t_{n+1}
174 5800824 : (context.*reinit_func)(1);
175 :
176 : // Move elem_->elem_, theta_->theta_
177 1860608 : context.get_elem_solution().swap(theta_solution);
178 5800824 : context.elem_solution_derivative = 1;
179 :
180 : // Add the constraint term
181 9741040 : jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
182 : jacobian_computed;
183 :
184 : // Add back (or restore) the old jacobian
185 5800824 : if (request_jacobian)
186 : {
187 1636620 : if (jacobian_computed)
188 1101980 : context.get_elem_jacobian() += old_elem_jacobian;
189 : else
190 0 : context.get_elem_jacobian().swap(old_elem_jacobian);
191 : }
192 :
193 7661432 : return jacobian_computed;
194 2079608 : }
195 :
196 720 : void EulerSolver::integrate_qoi_timestep()
197 : {
198 : // We use a numerical integration scheme consistent with the theta used for the timesolver.
199 :
200 : // Zero out the system.qoi vector
201 2160 : for (auto j : make_range(_system.n_qois()))
202 : {
203 1440 : _system.set_qoi(j, 0.0);
204 : }
205 :
206 : // Left and right side contributions
207 1200 : std::vector<Number> left_contribution(_system.n_qois(), 0.0);
208 240 : Number time_left = 0.0;
209 1200 : std::vector<Number> right_contribution(_system.n_qois(), 0.0);
210 240 : Number time_right = 0.0;
211 :
212 720 : time_left = _system.time;
213 :
214 : // Base class assumes a direct steady evaluation
215 960 : this->_system.assemble_qoi();
216 :
217 : // Also get the spatially integrated errors for all the QoIs in the QoI set
218 2160 : for (auto j : make_range(_system.n_qois()))
219 : {
220 1440 : left_contribution[j] = _system.get_qoi_value(j);
221 : }
222 :
223 : // Advance to t_j+1
224 720 : _system.time = _system.time + _system.deltat;
225 :
226 240 : time_right = _system.time;
227 :
228 : // Load the solution at the next timestep
229 720 : retrieve_timestep();
230 :
231 : // Zero out the system.qoi vector
232 2160 : for (auto j : make_range(_system.n_qois()))
233 : {
234 1440 : _system.set_qoi(j, 0.0);
235 : }
236 :
237 : // Base class assumes a direct steady evaluation
238 960 : this->_system.assemble_qoi();
239 :
240 2160 : for(auto j : make_range(_system.n_qois()))
241 : {
242 1440 : right_contribution[j] = _system.get_qoi_value(j);
243 : }
244 :
245 : // Combine the left and right side contributions as per the specified theta
246 : // theta = 0.5 (Crank-Nicholson) gives the trapezoidal rule.
247 2160 : for (auto j : make_range(_system.n_qois()))
248 : {
249 2400 : _system.set_qoi(j, ( ((1.0 - theta)*left_contribution[j]) + (theta*right_contribution[j]) )*(time_right - time_left));
250 : }
251 720 : }
252 :
253 : #ifdef LIBMESH_ENABLE_AMR
254 720 : void EulerSolver::integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error)
255 : {
256 : // Currently, we only support this functionality when Backward-Euler time integration is used.
257 720 : if (theta != 1.0)
258 0 : libmesh_not_implemented();
259 :
260 : // There are two possibilities regarding the integration rule we need to use for time integration.
261 : // If we have a instantaneous QoI, then we need to use a left sided Riemann sum, otherwise the trapezoidal rule for temporally smooth QoIs.
262 :
263 : // Create left and right error estimate vectors of the right size
264 1200 : std::vector<Number> qoi_error_estimates_left(_system.n_qois());
265 1200 : std::vector<Number> qoi_error_estimates_right(_system.n_qois());
266 :
267 : // Get t_j
268 720 : Real time_left = _system.time;
269 :
270 : // Get f(t_j)
271 480 : ErrorVector QoI_elementwise_error_left;
272 :
273 : // If we are at the very initial step, the error contribution is zero,
274 : // otherwise the old ajoint vector has been filled and we are the left end
275 : // of a subsequent timestep or sub-timestep
276 960 : if(old_adjoints[0] != nullptr)
277 : {
278 : // For evaluating the residual, we need to use the deltat that was used
279 : // to get us to this solution, so we save the current deltat as next_step_deltat
280 : // and set _system.deltat to the last completed deltat.
281 672 : next_step_deltat = _system.deltat;
282 672 : _system.deltat = last_step_deltat;
283 :
284 : // The adjoint error estimate expression for a backwards facing step
285 : // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
286 2016 : for (auto j : make_range(_system.n_qois()))
287 : {
288 : // Swap for residual weighting
289 1344 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
290 : }
291 :
292 672 : _system.update();
293 :
294 : // The residual has to be evaluated at the last time
295 672 : _system.time = _system.time - last_step_deltat;
296 :
297 672 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_left);
298 :
299 : // Shift the time back
300 672 : _system.time = _system.time + last_step_deltat;
301 :
302 : // Swap back the current and old adjoints
303 2016 : for (auto j : make_range(_system.n_qois()))
304 : {
305 1344 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
306 : }
307 :
308 : // Set the system deltat back to what it should be to march to the next time
309 672 : _system.deltat = next_step_deltat;
310 :
311 : }
312 : else
313 : {
314 48 : for (auto i : index_range(QoI_elementwise_error))
315 0 : QoI_elementwise_error_left[i] = 0.0;
316 : }
317 :
318 : // Also get the left side contributions for the spatially integrated errors for all the QoIs in the QoI set
319 2160 : for (auto j : make_range(_system.n_qois()))
320 : {
321 : // Skip this QoI if not in the QoI Set
322 1440 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
323 : {
324 : // If we are at the initial time, the error contribution is zero
325 1920 : if(std::abs(_system.time) > TOLERANCE*sqrt(TOLERANCE))
326 : {
327 1792 : qoi_error_estimates_left[j] = adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j);
328 : }
329 : else
330 : {
331 128 : qoi_error_estimates_left[j] = 0.0;
332 : }
333 : }
334 : }
335 :
336 : // Advance to t_j+1
337 720 : _system.time = _system.time + _system.deltat;
338 :
339 : // Get t_j+1
340 240 : Real time_right = _system.time;
341 :
342 : // We will need to use the last step deltat for the weighted residual evaluation
343 720 : last_step_deltat = _system.deltat;
344 :
345 : // The adjoint error estimate expression for a backwards facing step
346 : // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
347 2160 : for (auto j : make_range(_system.n_qois()))
348 : {
349 1920 : old_adjoints[j] = _system.get_adjoint_solution(j).clone();
350 : }
351 :
352 : // Retrieve the state and adjoint vectors for the next time instant
353 720 : retrieve_timestep();
354 :
355 : // Swap for residual weighting
356 2160 : for (auto j : make_range(_system.n_qois()))
357 : {
358 1440 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
359 : }
360 :
361 : // Swap out the deltats as we did for the left side
362 720 : next_step_deltat = _system.deltat;
363 720 : _system.deltat = last_step_deltat;
364 :
365 : // Get f(t_j+1)
366 480 : ErrorVector QoI_elementwise_error_right;
367 :
368 720 : _system.update();
369 :
370 : // The residual has to be evaluated at the last time
371 720 : _system.time = _system.time - last_step_deltat;
372 :
373 720 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_right);
374 :
375 : // Shift the time back
376 720 : _system.time = _system.time + last_step_deltat;
377 :
378 : // Set the system deltat back to what it needs to be able to march to the next time
379 720 : _system.deltat = next_step_deltat;
380 :
381 : // Swap back now that the residual weighting is done
382 2160 : for (auto j : make_range(_system.n_qois()))
383 : {
384 1440 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
385 : }
386 :
387 : // Also get the right side contributions for the spatially integrated errors for all the QoIs in the QoI set
388 2160 : for (auto j : make_range(_system.n_qois()))
389 : {
390 : // Skip this QoI if not in the QoI Set
391 1440 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
392 : {
393 1920 : qoi_error_estimates_right[j] = adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j);
394 : }
395 : }
396 :
397 : // Error contribution from this timestep
398 720 : for (auto i : index_range(QoI_elementwise_error))
399 0 : QoI_elementwise_error[i] = float(((QoI_elementwise_error_right[i] + QoI_elementwise_error_left[i])/2)
400 0 : * (time_right - time_left));
401 :
402 : // QoI set spatially integrated errors contribution from this timestep
403 2160 : for (auto j : make_range(_system.n_qois()))
404 : {
405 : // Skip this QoI if not in the QoI Set
406 1440 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
407 : {
408 2400 : _system.set_qoi_error_estimate(j, ( (1.0 - theta)*qoi_error_estimates_left[j] + theta*qoi_error_estimates_right[j] )*last_step_deltat);
409 : }
410 : }
411 :
412 720 : }
413 : #endif // LIBMESH_ENABLE_AMR
414 :
415 : } // namespace libMesh
|
