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/twostep_time_solver.h"
20 :
21 : #include "libmesh/adjoint_refinement_estimator.h"
22 : #include "libmesh/diff_system.h"
23 : #include "libmesh/enum_norm_type.h"
24 : #include "libmesh/error_vector.h"
25 : #include "libmesh/euler_solver.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/numeric_vector.h"
28 : #include "libmesh/parameter_vector.h"
29 : #include "libmesh/sensitivity_data.h"
30 : #include "libmesh/solution_history.h"
31 :
32 : // C++ includes
33 : #include <memory>
34 :
35 : namespace libMesh
36 : {
37 :
38 :
39 :
40 420 : TwostepTimeSolver::TwostepTimeSolver (sys_type & s)
41 420 : : AdaptiveTimeSolver(s)
42 :
43 : {
44 : // We start with a reasonable time solver: implicit Euler
45 420 : core_time_solver = std::make_unique<EulerSolver>(s);
46 420 : }
47 :
48 :
49 :
50 792 : TwostepTimeSolver::~TwostepTimeSolver () = default;
51 :
52 :
53 :
54 3780 : void TwostepTimeSolver::solve()
55 : {
56 108 : libmesh_assert(core_time_solver.get());
57 :
58 : // All actual solution history operations are handled by the outer
59 : // solver, so the outer solver has to call advance_timestep and
60 : // call solution_history->store to store the initial conditions
61 3780 : if (first_solve)
62 : {
63 420 : advance_timestep();
64 420 : first_solve = false;
65 : }
66 :
67 : // We may have to repeat timesteps entirely if our error is bad
68 : // enough
69 108 : bool max_tolerance_met = false;
70 :
71 : // Calculating error values each time
72 3780 : Real single_norm(0.), double_norm(0.), error_norm(0.),
73 108 : relative_error(0.);
74 :
75 : // The loop below will change system time and deltat based on calculations.
76 : // We will need to save these for calculating the deltat for the next timestep
77 : // after the while loop has converged.
78 : Real old_time;
79 : Real old_deltat;
80 :
81 7560 : while (!max_tolerance_met)
82 : {
83 : // If we've been asked to reduce deltat if necessary, make sure
84 : // the core timesolver does so
85 3888 : core_time_solver->reduce_deltat_on_diffsolver_failure =
86 3780 : this->reduce_deltat_on_diffsolver_failure;
87 :
88 3780 : if (!quiet)
89 : {
90 108 : libMesh::out << "\n === Computing adaptive timestep === "
91 108 : << std::endl;
92 : }
93 :
94 : // Use the double-length timestep first (so the
95 : // old_nonlinear_solution won't have to change)
96 3780 : core_time_solver->solve();
97 :
98 : // Save a copy of the double-length nonlinear solution
99 : // and the old nonlinear solution
100 : std::unique_ptr<NumericVector<Number>> double_solution =
101 3888 : _system.solution->clone();
102 : std::unique_ptr<NumericVector<Number>> old_solution =
103 3780 : _system.get_vector("_old_nonlinear_solution").clone();
104 :
105 3888 : double_norm = calculate_norm(_system, *double_solution);
106 3780 : if (!quiet)
107 : {
108 108 : libMesh::out << "Double norm = " << double_norm << std::endl;
109 : }
110 :
111 : // Then reset the initial guess for our single-length calcs
112 3888 : *(_system.solution) = _system.get_vector("_old_nonlinear_solution");
113 :
114 : // Call two single-length timesteps
115 : // Be sure that the core_time_solver does not change the
116 : // timestep here. (This is unlikely because it just succeeded
117 : // with a timestep twice as large!)
118 : // FIXME: even if diffsolver failure is unlikely, we ought to
119 : // do *something* if it happens
120 3780 : core_time_solver->reduce_deltat_on_diffsolver_failure = 0;
121 :
122 3780 : old_time = _system.time;
123 3780 : old_deltat = _system.deltat;
124 3780 : _system.deltat *= 0.5;
125 :
126 : // Attempt the 'half timestep solve'
127 3780 : core_time_solver->solve();
128 :
129 : // If we successfully completed the solve, let the time solver know the deltat used
130 3780 : this->last_deltat = _system.deltat;
131 :
132 : // Increment system.time, and save the half solution to solution history
133 3780 : core_time_solver->advance_timestep();
134 :
135 3780 : core_time_solver->solve();
136 :
137 3888 : single_norm = calculate_norm(_system, *_system.solution);
138 3780 : if (!quiet)
139 : {
140 108 : libMesh::out << "Single norm = " << single_norm << std::endl;
141 : }
142 :
143 : // Reset the core_time_solver's reduce_deltat... value.
144 3888 : core_time_solver->reduce_deltat_on_diffsolver_failure =
145 3780 : this->reduce_deltat_on_diffsolver_failure;
146 :
147 : // Find the relative error
148 3888 : *double_solution -= *(_system.solution);
149 3888 : error_norm = calculate_norm(_system, *double_solution);
150 3888 : relative_error = error_norm / old_deltat /
151 3780 : std::max(double_norm, single_norm);
152 :
153 : // If the relative error makes no sense, we're done
154 3780 : if (!double_norm && !single_norm)
155 0 : return;
156 :
157 3780 : if (!quiet)
158 : {
159 108 : libMesh::out << "Error norm = " << error_norm << std::endl;
160 108 : libMesh::out << "Local relative error = "
161 3672 : << (error_norm /
162 216 : std::max(double_norm, single_norm))
163 108 : << std::endl;
164 108 : libMesh::out << "Global relative error = "
165 108 : << (error_norm / old_deltat /
166 216 : std::max(double_norm, single_norm))
167 108 : << std::endl;
168 108 : libMesh::out << "old delta t = " << old_deltat << std::endl;
169 : }
170 :
171 : // If our upper tolerance is negative, that means we want to set
172 : // it based on the first successful time step
173 3780 : if (this->upper_tolerance < 0)
174 0 : this->upper_tolerance = -this->upper_tolerance * relative_error;
175 :
176 : // If we haven't met our upper error tolerance, we'll have to
177 : // repeat this timestep entirely
178 3780 : if (this->upper_tolerance && relative_error > this->upper_tolerance)
179 : {
180 : // If we are saving solution histories, the core time solver
181 : // will save half solutions, and these solves can be attempted
182 : // repeatedly. If we failed to meet the tolerance, erase the
183 : // half solution from solution history.
184 0 : core_time_solver->get_solution_history().erase(_system.time);
185 :
186 : // We will be retrying this timestep with deltat/2, so restore
187 : // all the necessary state.
188 : // FIXME: this probably doesn't work with multistep methods
189 0 : _system.get_vector("_old_nonlinear_solution") = *old_solution;
190 0 : _system.time = old_time;
191 0 : _system.deltat = old_deltat;
192 :
193 : // Update to localize the old nonlinear solution
194 0 : core_time_solver->update();
195 :
196 : // Reset the initial guess for our next try
197 0 : *(_system.solution) =
198 0 : _system.get_vector("_old_nonlinear_solution");
199 :
200 : // Chop delta t in half
201 0 : _system.deltat /= 2.;
202 :
203 0 : if (!quiet)
204 : {
205 0 : libMesh::out << "Failed to meet upper error tolerance"
206 0 : << std::endl;
207 0 : libMesh::out << "Retrying with delta t = "
208 0 : << _system.deltat << std::endl;
209 : }
210 : }
211 : else
212 108 : max_tolerance_met = true;
213 :
214 3564 : }
215 :
216 : // We ended up taking two half steps of size system.deltat to
217 : // march our last time step.
218 3780 : this->last_deltat = _system.deltat;
219 3780 : this->completed_timestep_size = 2.0*_system.deltat;
220 :
221 : // TimeSolver::solve methods should leave system.time unchanged
222 3780 : _system.time = old_time;
223 :
224 : // Compare the relative error to the tolerance and adjust deltat
225 3780 : _system.deltat = old_deltat;
226 :
227 : // If our target tolerance is negative, that means we want to set
228 : // it based on the first successful time step
229 3780 : if (this->target_tolerance < 0)
230 0 : this->target_tolerance = -this->target_tolerance * relative_error;
231 :
232 : const Real global_shrink_or_growth_factor =
233 3780 : std::pow(this->target_tolerance / relative_error,
234 3780 : static_cast<Real>(1. / core_time_solver->error_order()));
235 :
236 : const Real local_shrink_or_growth_factor =
237 7560 : std::pow(this->target_tolerance /
238 3888 : (error_norm/std::max(double_norm, single_norm)),
239 3780 : static_cast<Real>(1. / (core_time_solver->error_order()+1.)));
240 :
241 3780 : if (!quiet)
242 : {
243 108 : libMesh::out << "The global growth/shrink factor is: "
244 108 : << global_shrink_or_growth_factor << std::endl;
245 108 : libMesh::out << "The local growth/shrink factor is: "
246 108 : << local_shrink_or_growth_factor << std::endl;
247 : }
248 :
249 : // The local s.o.g. factor is based on the expected **local**
250 : // truncation error for the timestepping method, the global
251 : // s.o.g. factor is based on the method's **global** truncation
252 : // error. You can shrink/grow the timestep to attempt to satisfy
253 : // either a global or local time-discretization error tolerance.
254 :
255 108 : Real shrink_or_growth_factor =
256 3780 : this->global_tolerance ? global_shrink_or_growth_factor :
257 : local_shrink_or_growth_factor;
258 :
259 3780 : if (this->max_growth && this->max_growth < shrink_or_growth_factor)
260 : {
261 3780 : if (!quiet && this->global_tolerance)
262 : {
263 108 : libMesh::out << "delta t is constrained by max_growth" << std::endl;
264 : }
265 3780 : shrink_or_growth_factor = this->max_growth;
266 : }
267 :
268 3780 : _system.deltat *= shrink_or_growth_factor;
269 :
270 : // Restrict deltat to max-allowable value if necessary
271 3780 : if ((this->max_deltat != 0.0) && (_system.deltat > this->max_deltat))
272 : {
273 0 : if (!quiet)
274 : {
275 0 : libMesh::out << "delta t is constrained by maximum-allowable delta t."
276 0 : << std::endl;
277 : }
278 0 : _system.deltat = this->max_deltat;
279 : }
280 :
281 : // Restrict deltat to min-allowable value if necessary
282 3780 : if ((this->min_deltat != 0.0) && (_system.deltat < this->min_deltat))
283 : {
284 0 : if (!quiet)
285 : {
286 0 : libMesh::out << "delta t is constrained by minimum-allowable delta t."
287 0 : << std::endl;
288 : }
289 0 : _system.deltat = this->min_deltat;
290 : }
291 :
292 3780 : if (!quiet)
293 : {
294 3780 : libMesh::out << "new delta t = " << _system.deltat << std::endl;
295 : }
296 : }
297 :
298 3780 : std::pair<unsigned int, Real> TwostepTimeSolver::adjoint_solve (const QoISet & qoi_indices)
299 : {
300 : // Take the first adjoint 'half timestep'
301 3780 : core_time_solver->adjoint_solve(qoi_indices);
302 :
303 : // We print the forward 'half solution' norms and we will do so for the adjoints if running in dbg.
304 : #ifdef DEBUG
305 296 : for(auto i : make_range(_system.n_qois()))
306 : {
307 376 : for(auto j : make_range(_system.n_vars()))
308 188 : libMesh::out<<"||Z_"<<"("<<_system.time<<";"<<i<<","<<j<<")||_H1: "<<_system.calculate_norm(_system.get_adjoint_solution(i), j,H1)<<std::endl;
309 : }
310 : #endif
311 :
312 : // Record the sub step deltat we used for the last adjoint solve.
313 3780 : last_deltat = _system.deltat;
314 :
315 : // Adjoint advance the timestep
316 3780 : core_time_solver->adjoint_advance_timestep();
317 :
318 : // We have to contend with the fact that the delta_t set by SolutionHistory will not be the
319 : // delta_t for the adjoint solve. At time t_i, the adjoint solve uses the same delta_t
320 : // as the primal solve, pulling the adjoint solution from t_i+1 to t_i.
321 : // FSH however sets delta_t to the value which takes us from t_i to t_i-1.
322 : // Therefore use the last_deltat for the solve and reset system delta_t after the solve.
323 3780 : Real temp_deltat = _system.deltat;
324 3780 : _system.deltat = last_deltat;
325 :
326 : // The second half timestep
327 3780 : std::pair<unsigned int, Real> full_adjoint_output = core_time_solver->adjoint_solve(qoi_indices);
328 :
329 : // Record the sub step deltat we used for the last adjoint solve and reset the system deltat to the
330 : // value set by SolutionHistory.
331 3780 : last_deltat = _system.deltat;
332 3780 : _system.deltat = temp_deltat;
333 :
334 : // Record the total size of the last timestep, for a 2StepTS, this is
335 : // simply twice the deltat for each sub(half) step.
336 3780 : this->completed_timestep_size = 2.0*_system.deltat;
337 :
338 3780 : return full_adjoint_output;
339 : }
340 :
341 2800 : void TwostepTimeSolver::integrate_qoi_timestep()
342 : {
343 : // Vectors to hold qoi contributions from the first and second half timesteps
344 2960 : std::vector<Number> qois_first_half(_system.n_qois(), 0.0);
345 2960 : std::vector<Number> qois_second_half(_system.n_qois(), 0.0);
346 :
347 : // First half contribution
348 2800 : core_time_solver->integrate_qoi_timestep();
349 :
350 8400 : for (auto j : make_range(_system.n_qois()))
351 : {
352 5600 : qois_first_half[j] = _system.get_qoi_value(j);
353 : }
354 :
355 : // Second half contribution
356 2800 : core_time_solver->integrate_qoi_timestep();
357 :
358 8400 : for (auto j : make_range(_system.n_qois()))
359 : {
360 5600 : qois_second_half[j] = _system.get_qoi_value(j);
361 : }
362 :
363 : // Zero out the system.qoi vector
364 8400 : for (auto j : make_range(_system.n_qois()))
365 : {
366 5600 : _system.set_qoi(j, 0.0);
367 : }
368 :
369 : // Add the contributions from the two halftimesteps to get the full QoI
370 : // contribution from this timestep
371 8400 : for (auto j : make_range(_system.n_qois()))
372 : {
373 5920 : _system.set_qoi(j, qois_first_half[j] + qois_second_half[j]);
374 : }
375 2800 : }
376 :
377 980 : void TwostepTimeSolver::integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities)
378 : {
379 : // We are using the trapezoidal rule to integrate each timestep, and then pooling the contributions here.
380 : // (f(t_j) + f(t_j+1/2))/2 (t_j+1/2 - t_j) + (f(t_j+1/2) + f(t_j+1))/2 (t_j+1 - t_j+1/2)
381 :
382 : // First half timestep
383 1008 : SensitivityData sensitivities_first_half(qois, _system, parameter_vector);
384 :
385 980 : core_time_solver->integrate_adjoint_sensitivity(qois, parameter_vector, sensitivities_first_half);
386 :
387 : // Second half timestep
388 1008 : SensitivityData sensitivities_second_half(qois, _system, parameter_vector);
389 :
390 980 : core_time_solver->integrate_adjoint_sensitivity(qois, parameter_vector, sensitivities_second_half);
391 :
392 : // Get the contributions for each sensitivity from this timestep
393 28 : const auto pv_size = parameter_vector.size();
394 2912 : for (auto i : make_range(qois.size(_system)))
395 1960 : for (auto j : make_range(pv_size))
396 1064 : sensitivities[i][j] = sensitivities_first_half[i][j] + sensitivities_second_half[i][j];
397 980 : }
398 :
399 : #ifdef LIBMESH_ENABLE_AMR
400 2800 : void TwostepTimeSolver::integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error)
401 : {
402 : // We use a numerical integration scheme consistent with the theta used for the timesolver.
403 :
404 : // Create first and second half error estimate vectors of the right size
405 2960 : std::vector<Number> qoi_error_estimates_first_half(_system.n_qois());
406 2960 : std::vector<Number> qoi_error_estimates_second_half(_system.n_qois());
407 :
408 : // First half timestep
409 160 : ErrorVector QoI_elementwise_error_first_half;
410 :
411 2800 : core_time_solver->integrate_adjoint_refinement_error_estimate(adjoint_refinement_error_estimator, QoI_elementwise_error_first_half);
412 :
413 : // Also get the first 'half step' spatially integrated errors for all the QoIs in the QoI set
414 8400 : for (auto j : make_range(_system.n_qois()))
415 : {
416 : // Skip this QoI if not in the QoI Set
417 5600 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
418 : {
419 5600 : qoi_error_estimates_first_half[j] = _system.get_qoi_error_estimate_value(j);
420 : }
421 : }
422 :
423 : // Second half timestep
424 160 : ErrorVector QoI_elementwise_error_second_half;
425 :
426 2800 : core_time_solver->integrate_adjoint_refinement_error_estimate(adjoint_refinement_error_estimator, QoI_elementwise_error_second_half);
427 :
428 : // Also get the second 'half step' spatially integrated errors for all the QoIs in the QoI set
429 8400 : for (auto j : make_range(_system.n_qois()))
430 : {
431 : // Skip this QoI if not in the QoI Set
432 5600 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
433 : {
434 5600 : qoi_error_estimates_second_half[j] = _system.get_qoi_error_estimate_value(j);
435 : }
436 : }
437 :
438 : // Error contribution from this timestep
439 2800 : for (auto i : index_range(QoI_elementwise_error))
440 0 : QoI_elementwise_error[i] = QoI_elementwise_error_first_half[i] + QoI_elementwise_error_second_half[i];
441 :
442 8400 : for (auto j : make_range(_system.n_qois()))
443 : {
444 : // Skip this QoI if not in the QoI Set
445 5600 : if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
446 : {
447 5920 : _system.set_qoi_error_estimate(j, qoi_error_estimates_first_half[j] + qoi_error_estimates_second_half[j]);
448 : }
449 : }
450 2800 : }
451 : #endif // LIBMESH_ENABLE_AMR
452 :
453 : } // namespace libMesh
|
