1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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/diff_solver.h"
22 #include "libmesh/diff_system.h"
23 #include "libmesh/dof_map.h"
24 #include "libmesh/numeric_vector.h"
25 #include "libmesh/parameter_vector.h"
26 #include "libmesh/sensitivity_data.h"
27 #include "libmesh/solution_history.h"
28
29 namespace libMesh
30 {
31
32
33
UnsteadySolver(sys_type & s)34 UnsteadySolver::UnsteadySolver (sys_type & s)
35 : TimeSolver(s),
36 old_local_nonlinear_solution (NumericVector<Number>::build(s.comm())),
37 first_solve (true),
38 first_adjoint_step (true)
39 {
40 }
41
42
43
~UnsteadySolver()44 UnsteadySolver::~UnsteadySolver ()
45 {
46 }
47
48
49
init()50 void UnsteadySolver::init ()
51 {
52 TimeSolver::init();
53
54 _system.add_vector("_old_nonlinear_solution");
55 }
56
57
58
init_data()59 void UnsteadySolver::init_data()
60 {
61 TimeSolver::init_data();
62
63 #ifdef LIBMESH_ENABLE_GHOSTED
64 old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
65 _system.get_dof_map().get_send_list(), false,
66 GHOSTED);
67 #else
68 old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
69 #endif
70 }
71
72
73
reinit()74 void UnsteadySolver::reinit ()
75 {
76 TimeSolver::reinit();
77
78 #ifdef LIBMESH_ENABLE_GHOSTED
79 old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
80 _system.get_dof_map().get_send_list(), false,
81 GHOSTED);
82 #else
83 old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL);
84 #endif
85
86 // localize the old solution
87 NumericVector<Number> & old_nonlinear_soln =
88 _system.get_vector("_old_nonlinear_solution");
89
90 old_nonlinear_soln.localize
91 (*old_local_nonlinear_solution,
92 _system.get_dof_map().get_send_list());
93 }
94
95
96
solve()97 void UnsteadySolver::solve ()
98 {
99 if (first_solve)
100 {
101 advance_timestep();
102 first_solve = false;
103 }
104
105 unsigned int solve_result = _diff_solver->solve();
106
107 // If we requested the UnsteadySolver to attempt reducing dt after a
108 // failed DiffSolver solve, check the results of the solve now.
109 if (reduce_deltat_on_diffsolver_failure)
110 {
111 bool backtracking_failed =
112 solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
113
114 bool max_iterations =
115 solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
116
117 if (backtracking_failed || max_iterations)
118 {
119 // Cut timestep in half
120 for (unsigned int nr=0; nr<reduce_deltat_on_diffsolver_failure; ++nr)
121 {
122 _system.deltat *= 0.5;
123 libMesh::out << "Newton backtracking failed. Trying with smaller timestep, dt="
124 << _system.deltat << std::endl;
125
126 solve_result = _diff_solver->solve();
127
128 // Check solve results with reduced timestep
129 bool backtracking_still_failed =
130 solve_result & DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
131
132 bool backtracking_max_iterations =
133 solve_result & DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
134
135 if (!backtracking_still_failed && !backtracking_max_iterations)
136 {
137 if (!quiet)
138 libMesh::out << "Reduced dt solve succeeded." << std::endl;
139 return;
140 }
141 }
142
143 // If we made it here, we still couldn't converge the solve after
144 // reducing deltat
145 libMesh::out << "DiffSolver::solve() did not succeed after "
146 << reduce_deltat_on_diffsolver_failure
147 << " attempts." << std::endl;
148 libmesh_convergence_failure();
149
150 } // end if (backtracking_failed || max_iterations)
151 } // end if (reduce_deltat_on_diffsolver_failure)
152 }
153
154
155
advance_timestep()156 void UnsteadySolver::advance_timestep ()
157 {
158 // The first access of advance_timestep happens via solve, not user code
159 // It is used here to store any initial conditions data
160 if (!first_solve)
161 {
162 // We call advance_timestep in user code after solve, so any solutions
163 // we will be storing will be for the next time instance
164 _system.time += _system.deltat;
165 }
166 else
167 {
168 // We are here because of a call to advance_timestep that happens
169 // via solve, the very first solve. All we are doing here is storing
170 // the initial condition. The actual solution computed via this solve
171 // will be stored when we call advance_timestep in the user's timestep loop
172 first_solve = false;
173 }
174
175 // If the user has attached a memory or file solution history object
176 // to the solver, this will store the current solution indexed with
177 // the current time
178 solution_history->store(false, _system.time);
179
180 NumericVector<Number> & old_nonlinear_soln =
181 _system.get_vector("_old_nonlinear_solution");
182 NumericVector<Number> & nonlinear_solution =
183 *(_system.solution);
184
185 old_nonlinear_soln = nonlinear_solution;
186
187 old_nonlinear_soln.localize
188 (*old_local_nonlinear_solution,
189 _system.get_dof_map().get_send_list());
190 }
191
adjoint_solve(const QoISet & qoi_indices)192 std::pair<unsigned int, Real> UnsteadySolver::adjoint_solve(const QoISet & qoi_indices)
193 {
194 return _system.ImplicitSystem::adjoint_solve(qoi_indices);
195 }
196
adjoint_advance_timestep()197 void UnsteadySolver::adjoint_advance_timestep ()
198 {
199 // All calls to adjoint_advance_timestep are made in the user's
200 // code. This first call is made immediately after the adjoint initial conditions
201 // are set. This is in the user code outside the adjoint time loop.
202 if (!first_adjoint_step)
203 {
204 // The adjoint system has been solved. We need to store the adjoint solution and
205 // load the primal solutions for the next time instance (t - delta_ti).
206 _system.time -= _system.deltat;
207 }
208 else
209 {
210 // The first adjoint step simply saves the given adjoint initial condition
211 // So there is a store, but no solve, no actual timestep, so no need to change system time
212 first_adjoint_step = false;
213 }
214
215 // Retrieve the primal solution vectors at this new (or for
216 // first_adjoint_step, initial) time instance. These provide the
217 // data to solve the adjoint problem for the next time instance.
218 solution_history->retrieve(true, _system.time);
219
220 // Call the store function to store the adjoint we have computed (or
221 // for first_adjoint_step, the adjoint initial condition) in this
222 // time step for the time instance.
223 solution_history->store(true, _system.time);
224
225 // Dont forget to localize the old_nonlinear_solution !
226 _system.get_vector("_old_nonlinear_solution").localize
227 (*old_local_nonlinear_solution,
228 _system.get_dof_map().get_send_list());
229 }
230
retrieve_timestep()231 void UnsteadySolver::retrieve_timestep()
232 {
233 // Retrieve all the stored vectors at the current time
234 solution_history->retrieve(false, _system.time);
235
236 // Dont forget to localize the old_nonlinear_solution !
237 _system.get_vector("_old_nonlinear_solution").localize
238 (*old_local_nonlinear_solution,
239 _system.get_dof_map().get_send_list());
240 }
241
integrate_adjoint_sensitivity(const QoISet & qois,const ParameterVector & parameter_vector,SensitivityData & sensitivities)242 void UnsteadySolver::integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities)
243 {
244 // We are using the midpoint rule to integrate each timestep
245 // (f(t_j) + f(t_j+1))/2 (t_j+1 - t_j)
246
247 // Get t_j
248 Real time_left = _system.time;
249
250 // Left side sensitivities to hold f(t_j)
251 SensitivityData sensitivities_left(qois, _system, parameter_vector);
252
253 // Get f(t_j)
254 _system.adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivities_left);
255
256 // Advance to t_j+1
257 _system.time = _system.time + _system.deltat;
258
259 // Get t_j+1
260 Real time_right = _system.time;
261
262 // Right side sensitivities f(t_j+1)
263 SensitivityData sensitivities_right(qois, _system, parameter_vector);
264
265 // Remove the sensitivity rhs vector from system since we did not write it to file and it cannot be retrieved
266 _system.remove_vector("sensitivity_rhs0");
267
268 // Retrieve the primal and adjoint solutions at the current timestep
269 retrieve_timestep();
270
271 // Get f(t_j+1)
272 _system.adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivities_right);
273
274 // Remove the sensitivity rhs vector from system since we did not write it to file and it cannot be retrieved
275 _system.remove_vector("sensitivity_rhs0");
276
277 // Get the contributions for each sensitivity from this timestep
278 for(unsigned int i = 0; i != qois.size(_system); i++)
279 for(unsigned int j = 0; j != parameter_vector.size(); j++)
280 sensitivities[i][j] = ( (sensitivities_left[i][j] + sensitivities_right[i][j])/2. )*(time_right - time_left);
281
282 }
283
284
old_nonlinear_solution(const dof_id_type global_dof_number)285 Number UnsteadySolver::old_nonlinear_solution(const dof_id_type global_dof_number)
286 const
287 {
288 libmesh_assert_less (global_dof_number, _system.get_dof_map().n_dofs());
289 libmesh_assert_less (global_dof_number, old_local_nonlinear_solution->size());
290
291 return (*old_local_nonlinear_solution)(global_dof_number);
292 }
293
294
295
du(const SystemNorm & norm)296 Real UnsteadySolver::du(const SystemNorm & norm) const
297 {
298
299 std::unique_ptr<NumericVector<Number>> solution_copy =
300 _system.solution->clone();
301
302 solution_copy->add(-1., _system.get_vector("_old_nonlinear_solution"));
303
304 solution_copy->close();
305
306 return _system.calculate_norm(*solution_copy, norm);
307 }
308
309 } // namespace libMesh
310