Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local includes
20 : #include "libmesh/euler2_solver.h"
21 :
22 : #include "libmesh/adjoint_refinement_estimator.h"
23 : #include "libmesh/diff_system.h"
24 : #include "libmesh/error_vector.h"
25 : #include "libmesh/int_range.h"
26 : #include "libmesh/numeric_vector.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 426 : Euler2Solver::Euler2Solver (sys_type & s)
34 426 : : FirstOrderUnsteadySolver(s), theta(1.)
35 : {
36 426 : }
37 :
38 :
39 :
40 804 : Euler2Solver::~Euler2Solver () = default;
41 :
42 :
43 :
44 0 : Real Euler2Solver::error_order() const
45 : {
46 0 : if (theta == 0.5)
47 0 : return 2.;
48 0 : return 1.;
49 : }
50 :
51 :
52 :
53 :
54 327632 : bool Euler2Solver::element_residual (bool request_jacobian,
55 : DiffContext & context)
56 : {
57 327632 : bool compute_second_order_eqns = !this->_system.get_physics()->get_second_order_vars().empty();
58 :
59 357744 : return this->_general_residual(request_jacobian,
60 : context,
61 : &DifferentiablePhysics::mass_residual,
62 : &DifferentiablePhysics::damping_residual,
63 : &DifferentiablePhysics::_eulerian_time_deriv,
64 : &DifferentiablePhysics::element_constraint,
65 : &DiffContext::elem_reinit,
66 327632 : compute_second_order_eqns);
67 : }
68 :
69 :
70 :
71 185600 : bool Euler2Solver::side_residual (bool request_jacobian,
72 : DiffContext & context)
73 : {
74 201856 : return this->_general_residual(request_jacobian,
75 : context,
76 : &DifferentiablePhysics::side_mass_residual,
77 : &DifferentiablePhysics::side_damping_residual,
78 : &DifferentiablePhysics::side_time_derivative,
79 : &DifferentiablePhysics::side_constraint,
80 : &DiffContext::elem_side_reinit,
81 185600 : false);
82 : }
83 :
84 :
85 :
86 0 : bool Euler2Solver::nonlocal_residual (bool request_jacobian,
87 : DiffContext & context)
88 : {
89 0 : bool compute_second_order_eqns = this->_system.have_second_order_scalar_vars();
90 :
91 0 : return this->_general_residual(request_jacobian,
92 : context,
93 : &DifferentiablePhysics::nonlocal_mass_residual,
94 : &DifferentiablePhysics::nonlocal_damping_residual,
95 : &DifferentiablePhysics::nonlocal_time_derivative,
96 : &DifferentiablePhysics::nonlocal_constraint,
97 : &DiffContext::nonlocal_reinit,
98 0 : compute_second_order_eqns);
99 : }
100 :
101 :
102 :
103 513232 : bool Euler2Solver::_general_residual (bool request_jacobian,
104 : DiffContext & context,
105 : ResFuncType mass,
106 : ResFuncType damping,
107 : ResFuncType time_deriv,
108 : ResFuncType constraint,
109 : ReinitFuncType reinit_func,
110 : bool compute_second_order_eqns)
111 : {
112 46368 : unsigned int n_dofs = context.get_elem_solution().size();
113 :
114 : // Local nonlinear solution at old timestep
115 513232 : DenseVector<Number> old_elem_solution(n_dofs);
116 21556672 : for (unsigned int i=0; i != n_dofs; ++i)
117 21043440 : old_elem_solution(i) =
118 22903444 : old_nonlinear_solution(context.get_dof_indices()[i]);
119 :
120 : // Local time derivative of solution
121 46368 : context.get_elem_solution_rate() = context.get_elem_solution();
122 466864 : context.get_elem_solution_rate() -= old_elem_solution;
123 513232 : context.elem_solution_rate_derivative = 1 / _system.deltat;
124 46368 : context.get_elem_solution_rate() *=
125 71448 : context.elem_solution_rate_derivative;
126 :
127 : // Our first evaluations are at the final elem_solution
128 513232 : context.elem_solution_derivative = 1.0;
129 :
130 : // If a fixed solution is requested, we'll use the elem_solution
131 : // at the new timestep
132 : // FIXME - should this be the theta solution instead?
133 513232 : if (_system.use_fixed_solution)
134 0 : context.get_elem_fixed_solution() = context.get_elem_solution();
135 :
136 513232 : context.fixed_solution_derivative = 1.0;
137 :
138 : // We need to save the old jacobian and old residual since we'll be
139 : // multiplying some of the new contributions by theta or 1-theta
140 605968 : DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
141 513232 : DenseVector<Number> old_elem_residual(n_dofs);
142 46368 : old_elem_residual.swap(context.get_elem_residual());
143 513232 : if (request_jacobian)
144 233432 : old_elem_jacobian.swap(context.get_elem_jacobian());
145 :
146 : // Local time derivative of solution
147 46368 : context.get_elem_solution_rate() = context.get_elem_solution();
148 466864 : context.get_elem_solution_rate() -= old_elem_solution;
149 513232 : context.elem_solution_rate_derivative = 1 / _system.deltat;
150 46368 : context.get_elem_solution_rate() *=
151 71448 : context.elem_solution_rate_derivative;
152 :
153 : // If we are asked to compute residuals for second order variables,
154 : // we also populate the acceleration part so the user can use that.
155 513232 : if (compute_second_order_eqns)
156 326912 : this->prepare_accel(context);
157 :
158 : // Move the mesh into place first if necessary, set t = t_{n+1}
159 513232 : (context.*reinit_func)(1.);
160 :
161 : // First, evaluate time derivative at the new timestep.
162 : // The element should already be in the proper place
163 : // even for a moving mesh problem.
164 : bool jacobian_computed =
165 980096 : (_system.get_physics()->*time_deriv)(request_jacobian, context);
166 :
167 : // Next, evaluate the mass residual at the new timestep
168 :
169 980096 : jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
170 : jacobian_computed;
171 :
172 : // If we have second-order variables, we need to get damping terms
173 : // and the velocity equations
174 513232 : if (compute_second_order_eqns)
175 : {
176 623772 : jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
177 : jacobian_computed;
178 :
179 326912 : jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
180 : jacobian_computed;
181 : }
182 :
183 : // Add the constraint term
184 980096 : jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
185 : jacobian_computed;
186 :
187 : // The new solution's contribution is scaled by theta
188 513232 : context.get_elem_residual() *= theta;
189 513232 : context.get_elem_jacobian() *= theta;
190 :
191 : // Save the new solution's term
192 605968 : DenseMatrix<Number> elem_jacobian_newterm(n_dofs, n_dofs);
193 466864 : DenseVector<Number> elem_residual_newterm(n_dofs);
194 46368 : elem_residual_newterm.swap(context.get_elem_residual());
195 466864 : if (request_jacobian)
196 233432 : elem_jacobian_newterm.swap(context.get_elem_jacobian());
197 :
198 : // Add the time-dependent term for the old solution
199 :
200 : // Make sure elem_solution is set up for elem_reinit to use
201 : // Move elem_->old_, old_->elem_
202 46368 : context.get_elem_solution().swap(old_elem_solution);
203 513232 : context.elem_solution_derivative = 0.0;
204 :
205 : // Move the mesh into place if necessary, set t = t_{n}
206 513232 : (context.*reinit_func)(0.);
207 :
208 46368 : jacobian_computed =
209 980096 : (_system.get_physics()->*time_deriv)(jacobian_computed, context) &&
210 : jacobian_computed;
211 :
212 : // Add the mass residual term for the old solution
213 :
214 : // Evaluating the mass residual at both old and new timesteps will be
215 : // redundant in most problems but may be necessary for time accuracy
216 : // or stability in moving mesh problems or problems with user-overridden
217 : // mass_residual functions
218 :
219 46368 : jacobian_computed =
220 980096 : (_system.get_physics()->*mass)(jacobian_computed, context) &&
221 : jacobian_computed;
222 :
223 : // If we have second-order variables, we need to get damping terms
224 : // and the velocity equations
225 513232 : if (compute_second_order_eqns)
226 : {
227 623772 : jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
228 : jacobian_computed;
229 :
230 326912 : jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
231 : jacobian_computed;
232 : }
233 :
234 : // The old solution's contribution is scaled by (1-theta)
235 513232 : context.get_elem_residual() *= (1-theta);
236 513232 : context.get_elem_jacobian() *= (1-theta);
237 :
238 : // Restore the elem_solution
239 : // Move elem_->elem_, old_->old_
240 46368 : context.get_elem_solution().swap(old_elem_solution);
241 513232 : context.elem_solution_derivative = 1;
242 :
243 : // Restore the elem position if necessary, set t = t_{n+1}
244 513232 : (context.*reinit_func)(1.);
245 :
246 : // Add back (or restore) the old residual/jacobian
247 466864 : context.get_elem_residual() += old_elem_residual;
248 513232 : if (request_jacobian)
249 : {
250 256616 : if (jacobian_computed)
251 233432 : context.get_elem_jacobian() += old_elem_jacobian;
252 : else
253 0 : context.get_elem_jacobian().swap(old_elem_jacobian);
254 : }
255 :
256 : // Add the saved new-solution terms
257 466864 : context.get_elem_residual() += elem_residual_newterm;
258 513232 : if (jacobian_computed)
259 233432 : context.get_elem_jacobian() += elem_jacobian_newterm;
260 :
261 559600 : return jacobian_computed;
262 840992 : }
263 :
264 0 : void Euler2Solver::integrate_qoi_timestep()
265 : {
266 : // We use a numerical integration scheme consistent with the theta used for the timesolver.
267 :
268 : // Zero out the system.qoi vector
269 0 : for (auto j : make_range(_system.n_qois()))
270 : {
271 0 : _system.set_qoi(j, 0.0);
272 : }
273 :
274 : // Left and right side contributions
275 0 : std::vector<Number> left_contribution(_system.n_qois(), 0.0);
276 0 : Number time_left = 0.0;
277 0 : std::vector<Number> right_contribution(_system.n_qois(), 0.0);
278 0 : Number time_right = 0.0;
279 :
280 0 : time_left = _system.time;
281 :
282 : // Base class assumes a direct steady evaluation
283 0 : this->_system.assemble_qoi();
284 :
285 : // Also get the spatially integrated errors for all the QoIs in the QoI set
286 0 : for (auto j : make_range(_system.n_qois()))
287 : {
288 0 : left_contribution[j] = _system.get_qoi_value(j);
289 : }
290 :
291 : // Advance to t_j+1
292 0 : _system.time = _system.time + _system.deltat;
293 :
294 0 : time_right = _system.time;
295 :
296 : // Load the solution at the next timestep
297 0 : retrieve_timestep();
298 :
299 : // Zero out the system.qoi vector
300 0 : for (auto j : make_range(_system.n_qois()))
301 : {
302 0 : _system.set_qoi(j, 0.0);
303 : }
304 :
305 : // Base class assumes a direct steady evaluation
306 0 : this->_system.assemble_qoi();
307 :
308 0 : for(auto j : make_range(_system.n_qois()))
309 : {
310 0 : right_contribution[j] = _system.get_qoi_value(j);
311 : }
312 :
313 : // Combine the left and right side contributions as per the specified theta
314 : // theta = 0.5 (Crank-Nicholson) gives the trapezoidal rule.
315 0 : for (auto j : make_range(_system.n_qois()))
316 : {
317 0 : _system.set_qoi(j, ( ( ((1.0 - theta)*left_contribution[j]) + (theta*right_contribution[j]) )/2.0 )*(time_right - time_left));
318 : }
319 0 : }
320 :
321 : #ifdef LIBMESH_ENABLE_AMR
322 0 : void Euler2Solver::integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error)
323 : {
324 : // Currently, we only support this functionality when Backward-Euler time integration is used.
325 0 : if (theta != 1.0)
326 0 : libmesh_not_implemented();
327 :
328 : // There are two possibilities regarding the integration rule we need to use for time integration.
329 : // If we have a instantaneous QoI, then we need to use a left sided Riemann sum, otherwise the trapezoidal rule for temporally smooth QoIs.
330 :
331 : // Create left and right error estimate vectors of the right size
332 0 : std::vector<Number> qoi_error_estimates_left(_system.n_qois());
333 0 : std::vector<Number> qoi_error_estimates_right(_system.n_qois());
334 :
335 : // Get t_j
336 0 : Real time_left = _system.time;
337 :
338 : // Get f(t_j)
339 0 : ErrorVector QoI_elementwise_error_left;
340 :
341 : // If we are at the very initial step, the error contribution is zero,
342 : // otherwise the old ajoint vector has been filled and we are the left end
343 : // of a subsequent timestep or sub-timestep
344 0 : if(old_adjoints[0] != nullptr)
345 : {
346 : // For evaluating the residual, we need to use the deltat that was used
347 : // to get us to this solution, so we save the current deltat as next_step_deltat
348 : // and set _system.deltat to the last completed deltat.
349 0 : next_step_deltat = _system.deltat;
350 0 : _system.deltat = last_step_deltat;
351 :
352 : // The adjoint error estimate expression for a backwards facing step
353 : // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
354 0 : for (auto j : make_range(_system.n_qois()))
355 : {
356 : // Swap for residual weighting
357 0 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
358 : }
359 :
360 0 : _system.update();
361 :
362 : // The residual has to be evaluated at the last time
363 0 : _system.time = _system.time - last_step_deltat;
364 :
365 0 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_left);
366 :
367 : // Shift the time back
368 0 : _system.time = _system.time + last_step_deltat;
369 :
370 : // Swap back the current and old adjoints
371 0 : for (auto j : make_range(_system.n_qois()))
372 : {
373 0 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
374 : }
375 :
376 : // Set the system deltat back to what it should be to march to the next time
377 0 : _system.deltat = next_step_deltat;
378 :
379 : }
380 : else
381 : {
382 0 : for(auto i : index_range(QoI_elementwise_error))
383 0 : QoI_elementwise_error_left[i] = 0.0;
384 : }
385 :
386 : // Also get the left side contributions for the spatially integrated errors for all the QoIs in the QoI set
387 0 : for (auto j : make_range(_system.n_qois()))
388 : {
389 : // Skip this QoI if not in the QoI Set
390 0 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
391 : {
392 : // If we are at the initial time, the error contribution is zero
393 0 : if(std::abs(_system.time) > TOLERANCE*sqrt(TOLERANCE))
394 : {
395 0 : qoi_error_estimates_left[j] = adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j);
396 : }
397 : else
398 : {
399 0 : qoi_error_estimates_left[j] = 0.0;
400 : }
401 : }
402 : }
403 :
404 : // Advance to t_j+1
405 0 : _system.time = _system.time + _system.deltat;
406 :
407 : // Get t_j+1
408 0 : Real time_right = _system.time;
409 :
410 : // We will need to use the last step deltat for the weighted residual evaluation
411 0 : last_step_deltat = _system.deltat;
412 :
413 : // The adjoint error estimate expression for a backwards facing step
414 : // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
415 0 : for (auto j : make_range(_system.n_qois()))
416 : {
417 0 : old_adjoints[j] = _system.get_adjoint_solution(j).clone();
418 : }
419 :
420 : // Retrieve the state and adjoint vectors for the next time instant
421 0 : retrieve_timestep();
422 :
423 : // Swap for residual weighting
424 0 : for (auto j : make_range(_system.n_qois()))
425 : {
426 0 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
427 : }
428 :
429 : // Swap out the deltats as we did for the left side
430 0 : next_step_deltat = _system.deltat;
431 0 : _system.deltat = last_step_deltat;
432 :
433 : // Get f(t_j+1)
434 0 : ErrorVector QoI_elementwise_error_right;
435 :
436 0 : _system.update();
437 :
438 : // The residual has to be evaluated at the last time
439 0 : _system.time = _system.time - last_step_deltat;
440 :
441 0 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_right);
442 :
443 : // Shift the time back
444 0 : _system.time = _system.time + last_step_deltat;
445 :
446 : // Set the system deltat back to what it needs to be able to march to the next time
447 0 : _system.deltat = next_step_deltat;
448 :
449 : // Swap back now that the residual weighting is done
450 0 : for (auto j : make_range(_system.n_qois()))
451 : {
452 0 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
453 : }
454 :
455 : // Also get the right side contributions for the spatially integrated errors for all the QoIs in the QoI set
456 0 : for (auto j : make_range(_system.n_qois()))
457 : {
458 : // Skip this QoI if not in the QoI Set
459 0 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
460 : {
461 0 : qoi_error_estimates_right[j] = adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j);
462 : }
463 : }
464 :
465 : // Error contribution from this timestep
466 0 : for (auto i : index_range(QoI_elementwise_error))
467 0 : QoI_elementwise_error[i] = float(((QoI_elementwise_error_right[i] + QoI_elementwise_error_left[i])/2)
468 0 : * (time_right - time_left));
469 :
470 : // QoI set spatially integrated errors contribution from this timestep
471 0 : for (auto j : make_range(_system.n_qois()))
472 : {
473 : // Skip this QoI if not in the QoI Set
474 0 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
475 : {
476 0 : _system.set_qoi_error_estimate(j, ( (1.0 - theta)*qoi_error_estimates_left[j] + theta*qoi_error_estimates_right[j] )*last_step_deltat);
477 : }
478 : }
479 :
480 0 : }
481 : #endif // LIBMESH_ENABLE_AMR
482 :
483 : } // namespace libMesh
|
