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 : #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 1827 : EulerSolver::EulerSolver (sys_type & s)
33 1827 : : FirstOrderUnsteadySolver(s), theta(1.)
34 : {
35 1827 : }
36 :
37 :
38 :
39 3446 : EulerSolver::~EulerSolver () = default;
40 :
41 :
42 :
43 7560 : Real EulerSolver::error_order() const
44 : {
45 7560 : if (theta == 0.5)
46 0 : return 2.;
47 216 : return 1.;
48 : }
49 :
50 :
51 :
52 :
53 18934400 : bool EulerSolver::element_residual (bool request_jacobian,
54 : DiffContext & context)
55 : {
56 18934400 : bool compute_second_order_eqns = this->_system.get_physics()->have_second_order_vars();
57 :
58 20638000 : 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 18934400 : compute_second_order_eqns);
66 : }
67 :
68 :
69 :
70 1774032 : bool EulerSolver::side_residual (bool request_jacobian,
71 : DiffContext & context)
72 : {
73 1931040 : 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 1774032 : 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 20708432 : 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 24429648 : DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
116 20708432 : if (request_jacobian)
117 5373912 : old_elem_jacobian.swap(context.get_elem_jacobian());
118 :
119 : // Local nonlinear solution at old timestep
120 20708432 : DenseVector<Number> old_elem_solution(n_dofs);
121 169594752 : for (unsigned int i=0; i != n_dofs; ++i)
122 148886320 : old_elem_solution(i) =
123 161891284 : 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 18847824 : context.get_elem_solution_rate() -= old_elem_solution;
128 20708432 : 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 20708432 : if (compute_second_order_eqns)
135 326912 : this->prepare_accel(context);
136 :
137 : // Local nonlinear solution at time t_theta
138 1860608 : DenseVector<Number> theta_solution(context.get_elem_solution());
139 20708432 : theta_solution *= theta;
140 20708432 : theta_solution.add(1. - theta, old_elem_solution);
141 :
142 20708432 : context.elem_solution_derivative = theta;
143 20708432 : context.fixed_solution_derivative = theta;
144 :
145 : // If a fixed solution is requested, we'll use theta_solution
146 20708432 : 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 20708432 : (context.*reinit_func)(theta);
154 :
155 : // Get the time derivative at t_theta
156 : bool jacobian_computed =
157 39556256 : (_system.get_physics()->*time_deriv)(request_jacobian, context);
158 :
159 39556256 : 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 20708432 : if (compute_second_order_eqns)
165 : {
166 623772 : jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
167 : jacobian_computed;
168 :
169 326912 : 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 20708432 : (context.*reinit_func)(1);
175 :
176 : // Move elem_->elem_, theta_->theta_
177 1860608 : context.get_elem_solution().swap(theta_solution);
178 20708432 : context.elem_solution_derivative = 1;
179 :
180 : // Add the constraint term
181 39556256 : jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
182 : jacobian_computed;
183 :
184 : // Add back (or restore) the old jacobian
185 20708432 : if (request_jacobian)
186 : {
187 5908552 : if (jacobian_computed)
188 5373912 : context.get_elem_jacobian() += old_elem_jacobian;
189 : else
190 0 : context.get_elem_jacobian().swap(old_elem_jacobian);
191 : }
192 :
193 22569040 : return jacobian_computed;
194 16987216 : }
195 :
196 8400 : 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 25200 : for (auto j : make_range(_system.n_qois()))
202 : {
203 16800 : _system.set_qoi(j, 0.0);
204 : }
205 :
206 : // Left and right side contributions
207 8880 : std::vector<Number> left_contribution(_system.n_qois(), 0.0);
208 240 : Number time_left = 0.0;
209 8880 : std::vector<Number> right_contribution(_system.n_qois(), 0.0);
210 240 : Number time_right = 0.0;
211 :
212 8400 : time_left = _system.time;
213 :
214 : // Base class assumes a direct steady evaluation
215 8640 : this->_system.assemble_qoi();
216 :
217 : // Also get the spatially integrated errors for all the QoIs in the QoI set
218 25200 : for (auto j : make_range(_system.n_qois()))
219 : {
220 16800 : left_contribution[j] = _system.get_qoi_value(j);
221 : }
222 :
223 : // Advance to t_j+1
224 8400 : _system.time = _system.time + _system.deltat;
225 :
226 240 : time_right = _system.time;
227 :
228 : // Load the solution at the next timestep
229 8400 : retrieve_timestep();
230 :
231 : // Zero out the system.qoi vector
232 25200 : for (auto j : make_range(_system.n_qois()))
233 : {
234 16800 : _system.set_qoi(j, 0.0);
235 : }
236 :
237 : // Base class assumes a direct steady evaluation
238 8640 : this->_system.assemble_qoi();
239 :
240 25200 : for(auto j : make_range(_system.n_qois()))
241 : {
242 16800 : 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 25200 : for (auto j : make_range(_system.n_qois()))
248 : {
249 17760 : _system.set_qoi(j, ( ((1.0 - theta)*left_contribution[j]) + (theta*right_contribution[j]) )*(time_right - time_left));
250 : }
251 8400 : }
252 :
253 : #ifdef LIBMESH_ENABLE_AMR
254 8400 : 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 8400 : 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 8880 : std::vector<Number> qoi_error_estimates_left(_system.n_qois());
265 8880 : std::vector<Number> qoi_error_estimates_right(_system.n_qois());
266 :
267 : // Get t_j
268 8400 : 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 8640 : 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 7840 : next_step_deltat = _system.deltat;
282 7840 : _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 23520 : for (auto j : make_range(_system.n_qois()))
287 : {
288 : // Swap for residual weighting
289 15680 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
290 : }
291 :
292 7840 : _system.update();
293 :
294 : // The residual has to be evaluated at the last time
295 7840 : _system.time = _system.time - last_step_deltat;
296 :
297 7840 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_left);
298 :
299 : // Shift the time back
300 7840 : _system.time = _system.time + last_step_deltat;
301 :
302 : // Swap back the current and old adjoints
303 23520 : for (auto j : make_range(_system.n_qois()))
304 : {
305 15680 : _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 7840 : _system.deltat = next_step_deltat;
310 :
311 : }
312 : else
313 : {
314 560 : 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 25200 : for (auto j : make_range(_system.n_qois()))
320 : {
321 : // Skip this QoI if not in the QoI Set
322 16800 : 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 17280 : if(std::abs(_system.time) > TOLERANCE*sqrt(TOLERANCE))
326 : {
327 16128 : qoi_error_estimates_left[j] = adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j);
328 : }
329 : else
330 : {
331 1152 : qoi_error_estimates_left[j] = 0.0;
332 : }
333 : }
334 : }
335 :
336 : // Advance to t_j+1
337 8400 : _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 8400 : 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 25200 : for (auto j : make_range(_system.n_qois()))
348 : {
349 17280 : old_adjoints[j] = _system.get_adjoint_solution(j).clone();
350 : }
351 :
352 : // Retrieve the state and adjoint vectors for the next time instant
353 8400 : retrieve_timestep();
354 :
355 : // Swap for residual weighting
356 25200 : for (auto j : make_range(_system.n_qois()))
357 : {
358 16800 : _system.get_adjoint_solution(j).swap(*old_adjoints[j]);
359 : }
360 :
361 : // Swap out the deltats as we did for the left side
362 8400 : next_step_deltat = _system.deltat;
363 8400 : _system.deltat = last_step_deltat;
364 :
365 : // Get f(t_j+1)
366 480 : ErrorVector QoI_elementwise_error_right;
367 :
368 8400 : _system.update();
369 :
370 : // The residual has to be evaluated at the last time
371 8400 : _system.time = _system.time - last_step_deltat;
372 :
373 8400 : adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error_right);
374 :
375 : // Shift the time back
376 8400 : _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 8400 : _system.deltat = next_step_deltat;
380 :
381 : // Swap back now that the residual weighting is done
382 25200 : for (auto j : make_range(_system.n_qois()))
383 : {
384 16800 : _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 25200 : for (auto j : make_range(_system.n_qois()))
389 : {
390 : // Skip this QoI if not in the QoI Set
391 16800 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
392 : {
393 17280 : 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 8400 : 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 25200 : for (auto j : make_range(_system.n_qois()))
404 : {
405 : // Skip this QoI if not in the QoI Set
406 16800 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
407 : {
408 17760 : _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 8400 : }
413 : #endif // LIBMESH_ENABLE_AMR
414 :
415 : } // namespace libMesh
|
