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/unsteady_solver.h"
20 :
21 : #include "libmesh/adjoint_refinement_estimator.h"
22 : #include "libmesh/diff_solver.h"
23 : #include "libmesh/diff_system.h"
24 : #include "libmesh/dof_map.h"
25 : #include "libmesh/error_vector.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 : namespace libMesh
33 : {
34 :
35 :
36 :
37 3099 : UnsteadySolver::UnsteadySolver (sys_type & s)
38 : : TimeSolver(s),
39 6110 : old_local_nonlinear_solution (NumericVector<Number>::build(s.comm()).release()),
40 2923 : first_solve (true),
41 6286 : first_adjoint_step (true)
42 : {
43 3099 : old_adjoints.resize(s.n_qois());
44 :
45 : // Set the old adjoint pointers to nullptrs
46 : // We will use this nullness to skip the initial time instant,
47 : // when there is no older adjoint.
48 5899 : for(auto j : make_range(s.n_qois()))
49 : {
50 2800 : old_adjoints[j] = nullptr;
51 : }
52 3099 : }
53 :
54 :
55 :
56 5846 : UnsteadySolver::~UnsteadySolver () = default;
57 :
58 :
59 :
60 2259 : void UnsteadySolver::init ()
61 : {
62 2259 : TimeSolver::init();
63 :
64 2259 : _system.add_vector("_old_nonlinear_solution");
65 2259 : }
66 :
67 :
68 840 : void UnsteadySolver::init_adjoints ()
69 : {
70 840 : TimeSolver::init_adjoints();
71 :
72 : // Add old adjoint solutions
73 : // To keep the number of vectors consistent between the primal and adjoint
74 : // time loops, we will also add the adjoint rhs vector during initialization
75 2240 : for(auto i : make_range(_system.n_qois()))
76 : {
77 1440 : std::string old_adjoint_solution_name = "_old_adjoint_solution";
78 1400 : old_adjoint_solution_name+= std::to_string(i);
79 1440 : _system.add_vector(old_adjoint_solution_name, false, GHOSTED);
80 :
81 1440 : std::string adjoint_rhs_name = "adjoint_rhs";
82 1360 : adjoint_rhs_name+= std::to_string(i);
83 1440 : _system.add_vector(adjoint_rhs_name, false, GHOSTED);
84 : }
85 :
86 840 : }
87 :
88 :
89 2259 : void UnsteadySolver::init_data()
90 : {
91 2259 : TimeSolver::init_data();
92 :
93 : #ifdef LIBMESH_ENABLE_GHOSTED
94 2323 : old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
95 2259 : _system.get_dof_map().get_send_list(), false,
96 128 : GHOSTED);
97 : #else
98 : old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
99 : #endif
100 2259 : }
101 :
102 :
103 :
104 1120 : void UnsteadySolver::reinit ()
105 : {
106 1120 : TimeSolver::reinit();
107 :
108 : #ifdef LIBMESH_ENABLE_GHOSTED
109 1152 : old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
110 1120 : _system.get_dof_map().get_send_list(), false,
111 64 : GHOSTED);
112 : #else
113 : old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
114 : #endif
115 :
116 : // localize the old solution
117 : NumericVector<Number> & old_nonlinear_soln =
118 1120 : _system.get_vector("_old_nonlinear_solution");
119 :
120 : old_nonlinear_soln.localize
121 1120 : (*old_local_nonlinear_solution,
122 1120 : _system.get_dof_map().get_send_list());
123 1120 : }
124 :
125 :
126 :
127 27465 : void UnsteadySolver::solve ()
128 : {
129 27465 : if (first_solve)
130 : {
131 1839 : advance_timestep();
132 1839 : first_solve = false;
133 : }
134 :
135 27465 : unsigned int solve_result = _diff_solver->solve();
136 :
137 : // If we requested the UnsteadySolver to attempt reducing dt after a
138 : // failed DiffSolver solve, check the results of the solve now.
139 27465 : if (reduce_deltat_on_diffsolver_failure)
140 : {
141 0 : bool backtracking_failed =
142 0 : solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
143 :
144 0 : bool max_iterations =
145 0 : solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
146 :
147 0 : if (backtracking_failed || max_iterations)
148 : {
149 : // Cut timestep in half
150 0 : for (unsigned int nr=0; nr<reduce_deltat_on_diffsolver_failure; ++nr)
151 : {
152 0 : _system.deltat *= 0.5;
153 0 : libMesh::out << "Newton backtracking failed. Trying with smaller timestep, dt="
154 0 : << _system.deltat << std::endl;
155 :
156 0 : solve_result = _diff_solver->solve();
157 :
158 : // Check solve results with reduced timestep
159 0 : bool backtracking_still_failed =
160 0 : solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
161 :
162 0 : bool backtracking_max_iterations =
163 0 : solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
164 :
165 0 : if (!backtracking_still_failed && !backtracking_max_iterations)
166 : {
167 : // Set the successful deltat as the last deltat
168 0 : last_deltat = _system.deltat;
169 :
170 0 : if (!quiet)
171 0 : libMesh::out << "Reduced dt solve succeeded." << std::endl;
172 0 : return;
173 : }
174 : }
175 :
176 : // If we made it here, we still couldn't converge the solve after
177 : // reducing deltat
178 0 : libMesh::out << "DiffSolver::solve() did not succeed after "
179 0 : << reduce_deltat_on_diffsolver_failure
180 0 : << " attempts." << std::endl;
181 0 : libmesh_convergence_failure();
182 :
183 : } // end if (backtracking_failed || max_iterations)
184 : } // end if (reduce_deltat_on_diffsolver_failure)
185 :
186 : // Set the successful deltat as the last deltat
187 27465 : last_deltat = _system.deltat;
188 : }
189 :
190 :
191 :
192 21744 : void UnsteadySolver::advance_timestep ()
193 : {
194 : // The first access of advance_timestep happens via solve, not user code
195 : // It is used here to store any initial conditions data
196 21744 : if (!first_solve)
197 : {
198 : // We call advance_timestep in user code after solve, so any solutions
199 : // we will be storing will be for the next time instance
200 19905 : _system.time += _system.deltat;
201 : }
202 : else
203 : {
204 : // We are here because of a call to advance_timestep that happens
205 : // via solve, the very first solve. All we are doing here is storing
206 : // the initial condition. The actual solution computed via this solve
207 : // will be stored when we call advance_timestep in the user's timestep loop
208 1839 : first_solve = false;
209 : }
210 :
211 : // If the user has attached a memory or file solution history object
212 : // to the solver, this will store the current solution indexed with
213 : // the current time
214 21744 : solution_history->store(false, _system.time);
215 :
216 : NumericVector<Number> & old_nonlinear_soln =
217 21744 : _system.get_vector("_old_nonlinear_solution");
218 : NumericVector<Number> & nonlinear_solution =
219 21744 : *(_system.solution);
220 :
221 21744 : old_nonlinear_soln = nonlinear_solution;
222 :
223 : old_nonlinear_soln.localize
224 21744 : (*old_local_nonlinear_solution,
225 21744 : _system.get_dof_map().get_send_list());
226 21744 : }
227 :
228 11760 : std::pair<unsigned int, Real> UnsteadySolver::adjoint_solve(const QoISet & qoi_indices)
229 : {
230 11760 : std::pair<unsigned int, Real> adjoint_output = _system.ImplicitSystem::adjoint_solve(qoi_indices);
231 :
232 : // Record the deltat we used for this adjoint timestep. This was determined completely
233 : // by SolutionHistory::retrieve methods. The adjoint_solve methods should never change deltat.
234 11760 : last_deltat = _system.deltat;
235 :
236 11760 : return adjoint_output;
237 : }
238 :
239 8400 : void UnsteadySolver::adjoint_advance_timestep ()
240 : {
241 : // Call the store function to store the adjoint we have computed (or
242 : // for first_adjoint_step, the adjoint initial condition) in this
243 : // time step for the time instance.
244 8400 : solution_history->store(true, _system.time);
245 :
246 : // Before moving to the next time instant, copy over the current adjoint solutions into _old_adjoint_solutions
247 22680 : for(auto i : make_range(_system.n_qois()))
248 : {
249 14688 : std::string old_adjoint_solution_name = "_old_adjoint_solution";
250 13872 : old_adjoint_solution_name+= std::to_string(i);
251 14688 : NumericVector<Number> & old_adjoint_solution_i = _system.get_vector(old_adjoint_solution_name);
252 14280 : NumericVector<Number> & adjoint_solution_i = _system.get_adjoint_solution(i);
253 14280 : old_adjoint_solution_i = adjoint_solution_i;
254 : }
255 :
256 8400 : _system.time -= _system.deltat;
257 :
258 : // Retrieve the primal solution vectors at this new (or for
259 : // first_adjoint_step, initial) time instance. These provide the
260 : // data to solve the adjoint problem for the next time instance.
261 8400 : solution_history->retrieve(true, _system.time);
262 :
263 : // Dont forget to localize the old_nonlinear_solution !
264 8400 : _system.get_vector("_old_nonlinear_solution").localize
265 8400 : (*old_local_nonlinear_solution,
266 8400 : _system.get_dof_map().get_send_list());
267 8400 : }
268 :
269 21560 : void UnsteadySolver::retrieve_timestep()
270 : {
271 : // Retrieve all the stored vectors at the current time
272 21560 : solution_history->retrieve(false, _system.time);
273 :
274 : // Dont forget to localize the old_nonlinear_solution !
275 21560 : _system.get_vector("_old_nonlinear_solution").localize
276 21560 : (*old_local_nonlinear_solution,
277 21560 : _system.get_dof_map().get_send_list());
278 21560 : }
279 :
280 3360 : void UnsteadySolver::integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities)
281 : {
282 : // CURRENTLY using the trapezoidal rule to integrate each timestep
283 : // (f(t_j) + f(t_j+1))/2 (t_j+1 - t_j)
284 : // Fix me: This function needs to be moved to the EulerSolver classes like the
285 : // other integrate_timestep functions, and use an integration rule consistent with
286 : // the theta method used for the time integration.
287 :
288 : // Get t_j
289 3360 : Real time_left = _system.time;
290 :
291 : // Left side sensitivities to hold f(t_j)
292 3456 : SensitivityData sensitivities_left(qois, _system, parameter_vector);
293 :
294 : // Get f(t_j)
295 3360 : _system.adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivities_left);
296 :
297 : // Advance to t_j+1
298 3360 : _system.time = _system.time + _system.deltat;
299 :
300 : // Get t_j+1
301 96 : Real time_right = _system.time;
302 :
303 : // Right side sensitivities f(t_j+1)
304 3456 : SensitivityData sensitivities_right(qois, _system, parameter_vector);
305 :
306 : // Remove the sensitivity rhs vector from system since we did not write it to file and it cannot be retrieved
307 3360 : _system.remove_vector("sensitivity_rhs0");
308 :
309 : // Retrieve the primal and adjoint solutions at the current timestep
310 3360 : retrieve_timestep();
311 :
312 : // Get f(t_j+1)
313 3360 : _system.adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivities_right);
314 :
315 : // Remove the sensitivity rhs vector from system since we did not write it to file and it cannot be retrieved
316 3360 : _system.remove_vector("sensitivity_rhs0");
317 :
318 : // Get the contributions for each sensitivity from this timestep
319 96 : const auto pv_size = parameter_vector.size();
320 9888 : for (auto i : make_range(qois.size(_system)))
321 6720 : for (auto j : make_range(pv_size))
322 3648 : sensitivities[i][j] = ( (sensitivities_left[i][j] + sensitivities_right[i][j])/2. )*(time_right - time_left);
323 3360 : }
324 :
325 0 : void UnsteadySolver::integrate_qoi_timestep()
326 : {
327 0 : libmesh_not_implemented();
328 : }
329 :
330 : #ifdef LIBMESH_ENABLE_AMR
331 0 : void UnsteadySolver::integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & /*adjoint_refinement_error_estimator*/, ErrorVector & /*QoI_elementwise_error*/)
332 : {
333 0 : libmesh_not_implemented();
334 : }
335 : #endif // LIBMESH_ENABLE_AMR
336 :
337 176326912 : Number UnsteadySolver::old_nonlinear_solution(const dof_id_type global_dof_number)
338 : const
339 : {
340 15400200 : libmesh_assert_less (global_dof_number, _system.get_dof_map().n_dofs());
341 15400200 : libmesh_assert_less (global_dof_number, old_local_nonlinear_solution->size());
342 :
343 176326912 : return (*old_local_nonlinear_solution)(global_dof_number);
344 : }
345 :
346 :
347 :
348 0 : Real UnsteadySolver::du(const SystemNorm & norm) const
349 : {
350 :
351 : std::unique_ptr<NumericVector<Number>> solution_copy =
352 0 : _system.solution->clone();
353 :
354 0 : solution_copy->add(-1., _system.get_vector("_old_nonlinear_solution"));
355 :
356 0 : solution_copy->close();
357 :
358 0 : return _system.calculate_norm(*solution_copy, norm);
359 0 : }
360 :
361 0 : void UnsteadySolver::update()
362 : {
363 : // Dont forget to localize the old_nonlinear_solution !
364 0 : _system.get_vector("_old_nonlinear_solution").localize
365 0 : (*old_local_nonlinear_solution,
366 0 : _system.get_dof_map().get_send_list());
367 0 : }
368 :
369 : } // namespace libMesh
|
