Reference documentation for deal.II version 9.4.0
solver_fire.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2020 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_solver_fire_h
17#define dealii_solver_fire_h
18
19
20#include <deal.II/base/config.h>
21
23
25#include <deal.II/lac/solver.h>
26
27#include <functional>
28
29
31
32
35
88template <typename VectorType = Vector<double>>
89class SolverFIRE : public SolverBase<VectorType>
90{
91public:
96 {
102 explicit AdditionalData(const double initial_timestep = 0.1,
103 const double maximum_timestep = 1,
104 const double maximum_linfty_norm = 1);
105
109 const double initial_timestep;
110
114 const double maximum_timestep;
115
120 };
121
125 SolverFIRE(SolverControl & solver_control,
126 VectorMemory<VectorType> &vector_memory,
128
133 SolverFIRE(SolverControl & solver_control,
135
145 template <typename PreconditionerType = DiagonalMatrix<VectorType>>
146 void
147 solve(const std::function<double(VectorType &, const VectorType &)> &compute,
148 VectorType & x,
149 const PreconditionerType &inverse_mass_matrix);
150
156 template <typename MatrixType, typename PreconditionerType>
157 void
158 solve(const MatrixType & A,
159 VectorType & x,
160 const VectorType & b,
161 const PreconditionerType &preconditioner);
162
163protected:
170 virtual void
171 print_vectors(const unsigned int,
172 const VectorType &x,
173 const VectorType &v,
174 const VectorType &g) const;
175
180};
181
184/*------------------------- Implementation ----------------------------*/
185
186#ifndef DOXYGEN
187
188template <typename VectorType>
190 const double initial_timestep,
191 const double maximum_timestep,
192 const double maximum_linfty_norm)
193 : initial_timestep(initial_timestep)
194 , maximum_timestep(maximum_timestep)
195 , maximum_linfty_norm(maximum_linfty_norm)
196{
197 AssertThrow(initial_timestep > 0. && maximum_timestep > 0. &&
198 maximum_linfty_norm > 0.,
199 ExcMessage("Expected positive values for initial_timestep, "
200 "maximum_timestep and maximum_linfty_norm but one "
201 "or more of the these values are not positive."));
202}
203
204
205
206template <typename VectorType>
208 VectorMemory<VectorType> &vector_memory,
210 : SolverBase<VectorType>(solver_control, vector_memory)
212{}
213
214
215
216template <typename VectorType>
219 : SolverBase<VectorType>(solver_control)
221{}
222
223
224
225template <typename VectorType>
226template <typename PreconditionerType>
227void
229 const std::function<double(VectorType &, const VectorType &)> &compute,
230 VectorType & x,
231 const PreconditionerType &inverse_mass_matrix)
232{
233 LogStream::Prefix prefix("FIRE");
234
235 // FIRE algorithm constants
236 const double DELAYSTEP = 5;
237 const double TIMESTEP_GROW = 1.1;
238 const double TIMESTEP_SHRINK = 0.5;
239 const double ALPHA_0 = 0.1;
240 const double ALPHA_SHRINK = 0.99;
241
242 using real_type = typename VectorType::real_type;
243
244 typename VectorMemory<VectorType>::Pointer v(this->memory);
245 typename VectorMemory<VectorType>::Pointer g(this->memory);
246
247 // Set velocities to zero but not gradients
248 // as we are going to compute them soon.
249 v->reinit(x, false);
250 g->reinit(x, true);
251
252 // Refer to v and g with some readable names.
253 VectorType &velocities = *v;
255
256 // Update gradients for the new x.
258
259 unsigned int iter = 0;
260
263 if (conv != SolverControl::iterate)
264 return;
265
267 const auto &maximum_timestep = additional_data.maximum_timestep;
269
270 // First scaling factor.
271 double alpha = ALPHA_0;
272
273 unsigned int previous_iter_with_positive_v_dot_g = 0;
274
275 while (conv == SolverControl::iterate)
276 {
277 ++iter;
278 // Euler integration step.
279 x.add(timestep, velocities); // x += dt * v
282
283 // Compute gradients for the new x.
285
287 conv = this->iteration_status(iter, gradient_norm_squared, x);
288 if (conv != SolverControl::iterate)
289 break;
290
291 // v_dot_g = V * G
292 const real_type v_dot_g = velocities * gradients;
293
294 if (v_dot_g < 0.)
295 {
296 const real_type velocities_norm_squared = velocities * velocities;
297
298 // Check if we divide by zero in DEBUG mode.
300
301 // beta = - alpha |V|/|G|
302 const real_type beta =
303 -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
304
305 // V = (1-alpha) V + beta G.
307
308 if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
309 {
310 // Increase timestep and decrease alpha.
311 timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
312 alpha *= ALPHA_SHRINK;
313 }
314 }
315 else
316 {
317 // Decrease timestep, reset alpha and set V = 0.
318 previous_iter_with_positive_v_dot_g = iter;
319 timestep *= TIMESTEP_SHRINK;
320 alpha = ALPHA_0;
321 velocities = 0.;
322 }
323
324 real_type vmax = velocities.linfty_norm();
325
326 // Change timestep if any dof would move more than maximum_linfty_norm.
327 if (vmax > 0.)
328 {
329 const double minimal_timestep =
331 if (minimal_timestep < timestep)
332 timestep = minimal_timestep;
333 }
334
336
337 } // While we need to iterate.
338
339 // In the case of failure: throw exception.
340 if (conv != SolverControl::success)
341 AssertThrow(false,
343}
344
345
346
347template <typename VectorType>
348template <typename MatrixType, typename PreconditionerType>
349void
350SolverFIRE<VectorType>::solve(const MatrixType & A,
351 VectorType & x,
352 const VectorType & b,
353 const PreconditionerType &preconditioner)
354{
355 std::function<double(VectorType &, const VectorType &)> compute_func =
356 [&](VectorType &g, const VectorType &x) -> double {
357 // Residual of the quadratic form @f$\frac{1}{2} xAx - xb @f$.
358 // G = b - Ax
359 A.residual(g, x, b);
360
361 // Gradient G = Ax -b.
362 g *= -1.;
363
364 // The quadratic form @f$\frac{1}{2} xAx - xb @f$.
365 return 0.5 * A.matrix_norm_square(x) - x * b;
366 };
367
368 this->solve(compute_func, x, preconditioner);
369}
370
371
372
373template <typename VectorType>
374void
376 const VectorType &,
377 const VectorType &,
378 const VectorType &) const
379{}
380
381
382
383#endif // DOXYGEN
384
386
387#endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
Definition: solver_fire.h:179
virtual void print_vectors(const unsigned int, const VectorType &x, const VectorType &v, const VectorType &g) const
void solve(const std::function< double(VectorType &, const VectorType &)> &compute, VectorType &x, const PreconditionerType &inverse_mass_matrix)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
static const char A
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const double maximum_timestep
Definition: solver_fire.h:114
const double initial_timestep
Definition: solver_fire.h:109
const double maximum_linfty_norm
Definition: solver_fire.h:119
AdditionalData(const double initial_timestep=0.1, const double maximum_timestep=1, const double maximum_linfty_norm=1)