Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_fire.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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
9 // Public License as published by the Free Software Foundation; either
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 
22 #include <deal.II/base/logstream.h>
23 
24 #include <deal.II/lac/diagonal_matrix.h>
25 #include <deal.II/lac/solver.h>
26 
27 #include <functional>
28 
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
35 
90 template <typename VectorType = Vector<double>>
91 class SolverFIRE : public SolverBase<VectorType>
92 {
93 public:
98  {
104  explicit AdditionalData(const double initial_timestep = 0.1,
105  const double maximum_timestep = 1,
106  const double maximum_linfty_norm = 1);
107 
111  const double initial_timestep;
112 
116  const double maximum_timestep;
117 
121  const double maximum_linfty_norm;
122  };
123 
127  SolverFIRE(SolverControl & solver_control,
128  VectorMemory<VectorType> &vector_memory,
129  const AdditionalData & data = AdditionalData());
130 
135  SolverFIRE(SolverControl & solver_control,
136  const AdditionalData &data = AdditionalData());
137 
141  virtual ~SolverFIRE();
142 
152  template <typename PreconditionerType = DiagonalMatrix<VectorType>>
153  void
154  solve(const std::function<double(VectorType &, const VectorType &)> &compute,
155  VectorType & x,
156  const PreconditionerType &inverse_mass_matrix);
157 
163  template <typename MatrixType, typename PreconditionerType>
164  void
165  solve(const MatrixType & A,
166  VectorType & x,
167  const VectorType & b,
168  const PreconditionerType &preconditioner);
169 
170 protected:
177  virtual void
178  print_vectors(const unsigned int,
179  const VectorType &x,
180  const VectorType &v,
181  const VectorType &g) const;
182 
187 };
188 
191 /*------------------------- Implementation ----------------------------*/
192 
193 #ifndef DOXYGEN
194 
195 template <typename VectorType>
197  const double initial_timestep,
198  const double maximum_timestep,
199  const double maximum_linfty_norm)
200  : initial_timestep(initial_timestep)
201  , maximum_timestep(maximum_timestep)
202  , maximum_linfty_norm(maximum_linfty_norm)
203 {
204  AssertThrow(initial_timestep > 0. && maximum_timestep > 0. &&
205  maximum_linfty_norm > 0.,
206  ExcMessage("Expected positive values for initial_timestep, "
207  "maximum_timestep and maximum_linfty_norm but one "
208  "or more of the these values are not positive."));
209 }
210 
211 
212 
213 template <typename VectorType>
215  VectorMemory<VectorType> &vector_memory,
216  const AdditionalData & data)
217  : SolverBase<VectorType>(solver_control, vector_memory)
218  , additional_data(data)
219 {}
220 
221 
222 
223 template <typename VectorType>
225  const AdditionalData &data)
226  : SolverBase<VectorType>(solver_control)
227  , additional_data(data)
228 {}
229 
230 
231 
232 template <typename VectorType>
234 {}
235 
236 
237 
238 template <typename VectorType>
239 template <typename PreconditionerType>
240 void
242  const std::function<double(VectorType &, const VectorType &)> &compute,
243  VectorType & x,
244  const PreconditionerType &inverse_mass_matrix)
245 {
246  LogStream::Prefix prefix("FIRE");
247 
248  // FIRE algorithm constants
249  const double DELAYSTEP = 5;
250  const double TIMESTEP_GROW = 1.1;
251  const double TIMESTEP_SHRINK = 0.5;
252  const double ALPHA_0 = 0.1;
253  const double ALPHA_SHRINK = 0.99;
254 
255  using real_type = typename VectorType::real_type;
256 
257  typename VectorMemory<VectorType>::Pointer v(this->memory);
258  typename VectorMemory<VectorType>::Pointer g(this->memory);
259 
260  // Set velocities to zero but not gradients
261  // as we are going to compute them soon.
262  v->reinit(x, false);
263  g->reinit(x, true);
264 
265  // Refer to v and g with some readable names.
266  VectorType &velocities = *v;
267  VectorType &gradients = *g;
268 
269  // Update gradients for the new x.
270  compute(gradients, x);
271 
272  unsigned int iter = 0;
273 
275  conv = this->iteration_status(iter, gradients * gradients, x);
276  if (conv != SolverControl::iterate)
277  return;
278 
279  // Refer to additional data members with some readable names.
280  const auto &maximum_timestep = additional_data.maximum_timestep;
281  double timestep = additional_data.initial_timestep;
282 
283  // First scaling factor.
284  double alpha = ALPHA_0;
285 
286  unsigned int previous_iter_with_positive_v_dot_g = 0;
287 
288  while (conv == SolverControl::iterate)
289  {
290  ++iter;
291  // Euler integration step.
292  x.add(timestep, velocities); // x += dt * v
293  inverse_mass_matrix.vmult(gradients, gradients); // g = M^{-1} * g
294  velocities.add(-timestep, gradients); // v -= dt * h
295 
296  // Compute gradients for the new x.
297  compute(gradients, x);
298 
299  const real_type gradient_norm_squared = gradients * gradients;
300  conv = this->iteration_status(iter, gradient_norm_squared, x);
301  if (conv != SolverControl::iterate)
302  break;
303 
304  // v_dot_g = V * G
305  const real_type v_dot_g = velocities * gradients;
306 
307  if (v_dot_g < 0.)
308  {
309  const real_type velocities_norm_squared = velocities * velocities;
310 
311  // Check if we divide by zero in DEBUG mode.
312  Assert(gradient_norm_squared > 0., ExcInternalError());
313 
314  // beta = - alpha |V|/|G|
315  const real_type beta =
316  -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
317 
318  // V = (1-alpha) V + beta G.
319  velocities.sadd(1. - alpha, beta, gradients);
320 
321  if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
322  {
323  // Increase timestep and decrease alpha.
324  timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
325  alpha *= ALPHA_SHRINK;
326  }
327  }
328  else
329  {
330  // Decrease timestep, reset alpha and set V = 0.
331  previous_iter_with_positive_v_dot_g = iter;
332  timestep *= TIMESTEP_SHRINK;
333  alpha = ALPHA_0;
334  velocities = 0.;
335  }
336 
337  real_type vmax = velocities.linfty_norm();
338 
339  // Change timestep if any dof would move more than maximum_linfty_norm.
340  if (vmax > 0.)
341  {
342  const double minimal_timestep =
343  additional_data.maximum_linfty_norm / vmax;
344  if (minimal_timestep < timestep)
345  timestep = minimal_timestep;
346  }
347 
348  print_vectors(iter, x, velocities, gradients);
349 
350  } // While we need to iterate.
351 
352  // In the case of failure: throw exception.
353  if (conv != SolverControl::success)
354  AssertThrow(false,
355  SolverControl::NoConvergence(iter, gradients * gradients));
356 }
357 
358 
359 
360 template <typename VectorType>
361 template <typename MatrixType, typename PreconditionerType>
362 void
363 SolverFIRE<VectorType>::solve(const MatrixType & A,
364  VectorType & x,
365  const VectorType & b,
366  const PreconditionerType &preconditioner)
367 {
368  std::function<double(VectorType &, const VectorType &)> compute_func =
369  [&](VectorType &g, const VectorType &x) -> double {
370  // Residual of the quadratic form @f$ \frac{1}{2} xAx - xb @f$.
371  // G = b - Ax
372  A.residual(g, x, b);
373 
374  // Gradient G = Ax -b.
375  g *= -1.;
376 
377  // The quadratic form @f$\frac{1}{2} xAx - xb @f$.
378  return 0.5 * A.matrix_norm_square(x) - x * b;
379  };
380 
381  this->solve(compute_func, x, preconditioner);
382 }
383 
384 
385 
386 template <typename VectorType>
387 void
388 SolverFIRE<VectorType>::print_vectors(const unsigned int,
389  const VectorType &,
390  const VectorType &,
391  const VectorType &) const
392 {}
393 
394 
395 
396 #endif // DOXYGEN
397 
398 DEAL_II_NAMESPACE_CLOSE
399 
400 #endif
virtual void print_vectors(const unsigned int, const VectorType &x, const VectorType &v, const VectorType &g) const
Continue iteration.
AdditionalData(const double initial_timestep=0.1, const double maximum_timestep=1, const double maximum_linfty_norm=1)
const AdditionalData additional_data
Definition: solver_fire.h:186
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcMessage(std::string arg1)
Stop iteration, goal reached.
const double maximum_timestep
Definition: solver_fire.h:116
#define Assert(cond, exc)
Definition: exceptions.h:1407
void solve(const std::function< double(VectorType &, const VectorType &)> &compute, VectorType &x, const PreconditionerType &inverse_mass_matrix)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
const double initial_timestep
Definition: solver_fire.h:111
const double maximum_linfty_norm
Definition: solver_fire.h:121
SolverFIRE(SolverControl &solver_control, VectorMemory< VectorType > &vector_memory, const AdditionalData &data=AdditionalData())
virtual ~SolverFIRE()
static ::ExceptionBase & ExcInternalError()