Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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
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 
25 #include <deal.II/lac/solver.h>
26 
27 #include <functional>
28 
29 
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 
147  template <typename PreconditionerType = DiagonalMatrix<VectorType>>
148  void
149  solve(const std::function<double(VectorType &, const VectorType &)> &compute,
150  VectorType & x,
151  const PreconditionerType &inverse_mass_matrix);
152 
158  template <typename MatrixType, typename PreconditionerType>
159  void
160  solve(const MatrixType & A,
161  VectorType & x,
162  const VectorType & b,
163  const PreconditionerType &preconditioner);
164 
165 protected:
172  virtual void
173  print_vectors(const unsigned int,
174  const VectorType &x,
175  const VectorType &v,
176  const VectorType &g) const;
177 
182 };
183 
186 /*------------------------- Implementation ----------------------------*/
187 
188 #ifndef DOXYGEN
189 
190 template <typename VectorType>
192  const double initial_timestep,
193  const double maximum_timestep,
194  const double maximum_linfty_norm)
195  : initial_timestep(initial_timestep)
196  , maximum_timestep(maximum_timestep)
197  , maximum_linfty_norm(maximum_linfty_norm)
198 {
199  AssertThrow(initial_timestep > 0. && maximum_timestep > 0. &&
200  maximum_linfty_norm > 0.,
201  ExcMessage("Expected positive values for initial_timestep, "
202  "maximum_timestep and maximum_linfty_norm but one "
203  "or more of the these values are not positive."));
204 }
205 
206 
207 
208 template <typename VectorType>
210  VectorMemory<VectorType> &vector_memory,
211  const AdditionalData & data)
212  : SolverBase<VectorType>(solver_control, vector_memory)
213  , additional_data(data)
214 {}
215 
216 
217 
218 template <typename VectorType>
220  const AdditionalData &data)
221  : SolverBase<VectorType>(solver_control)
222  , additional_data(data)
223 {}
224 
225 
226 
227 template <typename VectorType>
228 template <typename PreconditionerType>
229 void
231  const std::function<double(VectorType &, const VectorType &)> &compute,
232  VectorType & x,
233  const PreconditionerType &inverse_mass_matrix)
234 {
235  LogStream::Prefix prefix("FIRE");
236 
237  // FIRE algorithm constants
238  const double DELAYSTEP = 5;
239  const double TIMESTEP_GROW = 1.1;
240  const double TIMESTEP_SHRINK = 0.5;
241  const double ALPHA_0 = 0.1;
242  const double ALPHA_SHRINK = 0.99;
243 
244  using real_type = typename VectorType::real_type;
245 
246  typename VectorMemory<VectorType>::Pointer v(this->memory);
247  typename VectorMemory<VectorType>::Pointer g(this->memory);
248 
249  // Set velocities to zero but not gradients
250  // as we are going to compute them soon.
251  v->reinit(x, false);
252  g->reinit(x, true);
253 
254  // Refer to v and g with some readable names.
255  VectorType &velocities = *v;
256  VectorType &gradients = *g;
257 
258  // Update gradients for the new x.
259  compute(gradients, x);
260 
261  unsigned int iter = 0;
262 
264  conv = this->iteration_status(iter, gradients * gradients, x);
265  if (conv != SolverControl::iterate)
266  return;
267 
268  // Refer to additional data members with some readable names.
269  const auto &maximum_timestep = additional_data.maximum_timestep;
270  double timestep = additional_data.initial_timestep;
271 
272  // First scaling factor.
273  double alpha = ALPHA_0;
274 
275  unsigned int previous_iter_with_positive_v_dot_g = 0;
276 
277  while (conv == SolverControl::iterate)
278  {
279  ++iter;
280  // Euler integration step.
281  x.add(timestep, velocities); // x += dt * v
282  inverse_mass_matrix.vmult(gradients, gradients); // g = M^{-1} * g
283  velocities.add(-timestep, gradients); // v -= dt * h
284 
285  // Compute gradients for the new x.
286  compute(gradients, x);
287 
288  const real_type gradient_norm_squared = gradients * gradients;
289  conv = this->iteration_status(iter, gradient_norm_squared, x);
290  if (conv != SolverControl::iterate)
291  break;
292 
293  // v_dot_g = V * G
294  const real_type v_dot_g = velocities * gradients;
295 
296  if (v_dot_g < 0.)
297  {
298  const real_type velocities_norm_squared = velocities * velocities;
299 
300  // Check if we divide by zero in DEBUG mode.
301  Assert(gradient_norm_squared > 0., ExcInternalError());
302 
303  // beta = - alpha |V|/|G|
304  const real_type beta =
305  -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
306 
307  // V = (1-alpha) V + beta G.
308  velocities.sadd(1. - alpha, beta, gradients);
309 
310  if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
311  {
312  // Increase timestep and decrease alpha.
313  timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
314  alpha *= ALPHA_SHRINK;
315  }
316  }
317  else
318  {
319  // Decrease timestep, reset alpha and set V = 0.
320  previous_iter_with_positive_v_dot_g = iter;
321  timestep *= TIMESTEP_SHRINK;
322  alpha = ALPHA_0;
323  velocities = 0.;
324  }
325 
326  real_type vmax = velocities.linfty_norm();
327 
328  // Change timestep if any dof would move more than maximum_linfty_norm.
329  if (vmax > 0.)
330  {
331  const double minimal_timestep =
333  if (minimal_timestep < timestep)
334  timestep = minimal_timestep;
335  }
336 
337  print_vectors(iter, x, velocities, gradients);
338 
339  } // While we need to iterate.
340 
341  // In the case of failure: throw exception.
342  if (conv != SolverControl::success)
343  AssertThrow(false,
344  SolverControl::NoConvergence(iter, gradients * gradients));
345 }
346 
347 
348 
349 template <typename VectorType>
350 template <typename MatrixType, typename PreconditionerType>
351 void
352 SolverFIRE<VectorType>::solve(const MatrixType & A,
353  VectorType & x,
354  const VectorType & b,
355  const PreconditionerType &preconditioner)
356 {
357  std::function<double(VectorType &, const VectorType &)> compute_func =
358  [&](VectorType &g, const VectorType &x) -> double {
359  // Residual of the quadratic form @f$ \frac{1}{2} xAx - xb @f$.
360  // G = b - Ax
361  A.residual(g, x, b);
362 
363  // Gradient G = Ax -b.
364  g *= -1.;
365 
366  // The quadratic form @f$\frac{1}{2} xAx - xb @f$.
367  return 0.5 * A.matrix_norm_square(x) - x * b;
368  };
369 
370  this->solve(compute_func, x, preconditioner);
371 }
372 
373 
374 
375 template <typename VectorType>
376 void
377 SolverFIRE<VectorType>::print_vectors(const unsigned int,
378  const VectorType &,
379  const VectorType &,
380  const VectorType &) const
381 {}
382 
383 
384 
385 #endif // DOXYGEN
386 
388 
389 #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)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:463
const AdditionalData additional_data
Definition: solver_fire.h:181
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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)
static const char A
const double initial_timestep
Definition: solver_fire.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
T min(const T &t, const MPI_Comm &mpi_communicator)
const double maximum_linfty_norm
Definition: solver_fire.h:121
SolverFIRE(SolverControl &solver_control, VectorMemory< VectorType > &vector_memory, const AdditionalData &data=AdditionalData())
VectorMemory< VectorType > & memory
Definition: solver.h:412
static ::ExceptionBase & ExcInternalError()