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