Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
solver_fire.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2022 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
23
25#include <deal.II/lac/solver.h>
26
27#include <functional>
28
29
31
32
90template <typename VectorType = Vector<double>>
91class SolverFIRE : public SolverBase<VectorType>
92{
93public:
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
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
165protected:
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
190template <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
208template <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
218template <typename VectorType>
220 const AdditionalData &data)
221 : SolverBase<VectorType>(solver_control)
222 , additional_data(data)
223{}
224
225
226
227template <typename VectorType>
228template <typename PreconditionerType>
229void
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 =
332 additional_data.maximum_linfty_norm / vmax;
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
349template <typename VectorType>
350template <typename MatrixType, typename PreconditionerType>
351void
352SolverFIRE<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
375template <typename VectorType>
376void
378 const VectorType &,
379 const VectorType &,
380 const VectorType &) const
381{}
382
383
384
385#endif // DOXYGEN
386
388
389#endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
SolverFIRE(SolverControl &solver_control, VectorMemory< VectorType > &vector_memory, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
virtual void print_vectors(const unsigned int, const VectorType &x, const VectorType &v, const VectorType &g) const
SolverFIRE(SolverControl &solver_control, const AdditionalData &data=AdditionalData())
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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_linfty_norm
AdditionalData(const double initial_timestep=0.1, const double maximum_timestep=1, const double maximum_linfty_norm=1)