Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_fire_h
16#define dealii_solver_fire_h
17
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/lac/solver.h>
26
27#include <functional>
28
29
31
32
90template <typename VectorType = Vector<double>>
92class SolverFIRE : public SolverBase<VectorType>
93{
94public:
99 {
105 explicit AdditionalData(const double initial_timestep = 0.1,
106 const double maximum_timestep = 1,
107 const double maximum_linfty_norm = 1);
108
112 const double initial_timestep;
113
117 const double maximum_timestep;
118
123 };
124
128 SolverFIRE(SolverControl &solver_control,
129 VectorMemory<VectorType> &vector_memory,
130 const AdditionalData &data = AdditionalData());
131
136 SolverFIRE(SolverControl &solver_control,
137 const AdditionalData &data = AdditionalData());
138
148 template <typename PreconditionerType = DiagonalMatrix<VectorType>>
149 void
150 solve(const std::function<double(VectorType &, const VectorType &)> &compute,
151 VectorType &x,
152 const PreconditionerType &inverse_mass_matrix);
153
159 template <typename MatrixType, typename PreconditionerType>
163 void solve(const MatrixType &A,
164 VectorType &x,
165 const VectorType &b,
166 const PreconditionerType &preconditioner);
167
168protected:
175 virtual void
176 print_vectors(const unsigned int,
177 const VectorType &x,
178 const VectorType &v,
179 const VectorType &g) const;
180
184 const AdditionalData additional_data;
185};
186
189/*------------------------- Implementation ----------------------------*/
190
191#ifndef DOXYGEN
192
193template <typename VectorType>
196 const double initial_timestep,
197 const double maximum_timestep,
198 const double maximum_linfty_norm)
199 : initial_timestep(initial_timestep)
200 , maximum_timestep(maximum_timestep)
201 , maximum_linfty_norm(maximum_linfty_norm)
202{
203 AssertThrow(initial_timestep > 0. && maximum_timestep > 0. &&
204 maximum_linfty_norm > 0.,
205 ExcMessage("Expected positive values for initial_timestep, "
206 "maximum_timestep and maximum_linfty_norm but one "
207 "or more of the these values are not positive."));
208}
209
210
211
212template <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
223template <typename VectorType>
226 const AdditionalData &data)
227 : SolverBase<VectorType>(solver_control)
228 , additional_data(data)
229{}
230
231
232
233template <typename VectorType>
235template <typename PreconditionerType>
237 const std::function<double(VectorType &, const VectorType &)> &compute,
238 VectorType &x,
239 const PreconditionerType &inverse_mass_matrix)
240{
241 LogStream::Prefix prefix("FIRE");
242
243 // FIRE algorithm constants
244 const double DELAYSTEP = 5;
245 const double TIMESTEP_GROW = 1.1;
246 const double TIMESTEP_SHRINK = 0.5;
247 const double ALPHA_0 = 0.1;
248 const double ALPHA_SHRINK = 0.99;
249
250 using real_type = typename VectorType::real_type;
251
252 typename VectorMemory<VectorType>::Pointer v(this->memory);
253 typename VectorMemory<VectorType>::Pointer g(this->memory);
254
255 // Set velocities to zero but not gradients
256 // as we are going to compute them soon.
257 v->reinit(x, false);
258 g->reinit(x, true);
259
260 // Refer to v and g with some readable names.
261 VectorType &velocities = *v;
262 VectorType &gradients = *g;
263
264 // Update gradients for the new x.
265 compute(gradients, x);
266
267 unsigned int iter = 0;
268
270 conv = this->iteration_status(iter, gradients * gradients, x);
271 if (conv != SolverControl::iterate)
272 return;
273
274 // Refer to additional data members with some readable names.
275 const auto &maximum_timestep = additional_data.maximum_timestep;
276 double timestep = additional_data.initial_timestep;
277
278 // First scaling factor.
279 double alpha = ALPHA_0;
280
281 unsigned int previous_iter_with_positive_v_dot_g = 0;
282
283 while (conv == SolverControl::iterate)
284 {
285 ++iter;
286 // Euler integration step.
287 x.add(timestep, velocities); // x += dt * v
288 inverse_mass_matrix.vmult(gradients, gradients); // g = M^{-1} * g
289 velocities.add(-timestep, gradients); // v -= dt * h
290
291 // Compute gradients for the new x.
292 compute(gradients, x);
293
294 const real_type gradient_norm_squared = gradients * gradients;
295 conv = this->iteration_status(iter, gradient_norm_squared, x);
296 if (conv != SolverControl::iterate)
297 break;
298
299 // v_dot_g = V * G
300 const real_type v_dot_g = velocities * gradients;
301
302 if (v_dot_g < 0.)
303 {
304 const real_type velocities_norm_squared = velocities * velocities;
305
306 // Check if we divide by zero in DEBUG mode.
307 Assert(gradient_norm_squared > 0., ExcInternalError());
308
309 // beta = - alpha |V|/|G|
310 const real_type beta =
311 -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
312
313 // V = (1-alpha) V + beta G.
314 velocities.sadd(1. - alpha, beta, gradients);
315
316 if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
317 {
318 // Increase timestep and decrease alpha.
319 timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
320 alpha *= ALPHA_SHRINK;
321 }
322 }
323 else
324 {
325 // Decrease timestep, reset alpha and set V = 0.
326 previous_iter_with_positive_v_dot_g = iter;
327 timestep *= TIMESTEP_SHRINK;
328 alpha = ALPHA_0;
329 velocities = 0.;
330 }
331
332 real_type vmax = velocities.linfty_norm();
333
334 // Change timestep if any dof would move more than maximum_linfty_norm.
335 if (vmax > 0.)
336 {
337 const double minimal_timestep =
338 additional_data.maximum_linfty_norm / vmax;
339 if (minimal_timestep < timestep)
340 timestep = minimal_timestep;
341 }
342
343 print_vectors(iter, x, velocities, gradients);
344
345 } // While we need to iterate.
346
347 // In the case of failure: throw exception.
348 if (conv != SolverControl::success)
349 AssertThrow(false,
350 SolverControl::NoConvergence(iter, gradients * gradients));
351}
352
353
354
355template <typename VectorType>
357template <typename MatrixType, typename PreconditionerType>
361void SolverFIRE<VectorType>::solve(const MatrixType &A,
362 VectorType &x,
363 const VectorType &b,
364 const PreconditionerType &preconditioner)
365{
366 std::function<double(VectorType &, const VectorType &)> compute_func =
367 [&](VectorType &g, const VectorType &x) -> double {
368 // Residual of the quadratic form @f$ \frac{1}{2} xAx - xb @f$.
369 // G = b - Ax
370 A.residual(g, x, b);
371
372 // Gradient G = Ax -b.
373 g *= -1.;
374
375 // The quadratic form @f$\frac{1}{2} xAx - xb @f$.
376 return 0.5 * A.matrix_norm_square(x) - x * b;
377 };
378
379 this->solve(compute_func, x, preconditioner);
380}
381
382
383
384template <typename VectorType>
386void SolverFIRE<VectorType>::print_vectors(const unsigned int,
387 const VectorType &,
388 const VectorType &,
389 const VectorType &) const
390{}
391
392
393
394#endif // DOXYGEN
395
397
398#endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
SolverFIRE(SolverControl &solver_control, VectorMemory< VectorType > &vector_memory, const AdditionalData &data=AdditionalData())
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
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)