15#ifndef dealii_solver_fire_h
16#define dealii_solver_fire_h
89template <
typename VectorType = Vector<
double>>
146 template <
typename PreconditionerType = DiagonalMatrix<VectorType>>
148 solve(
const std::function<
double(VectorType &,
const VectorType &)> &compute,
150 const PreconditionerType &inverse_mass_matrix);
157 template <
typename MatrixType,
typename PreconditionerType>
162 const PreconditionerType &preconditioner);
175 const VectorType &g)
const;
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)
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."));
207template <
typename VectorType>
210 const AdditionalData &data)
211 :
SolverBase<VectorType>(solver_control, vector_memory)
212 , additional_data(data)
217template <
typename VectorType>
219 const AdditionalData &data)
221 , additional_data(data)
226template <
typename VectorType>
227template <
typename PreconditionerType>
230 const std::function<
double(VectorType &,
const VectorType &)> &compute,
232 const PreconditionerType &inverse_mass_matrix)
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;
243 using real_type =
typename VectorType::real_type;
254 VectorType &velocities = *v;
258 compute(gradients, x);
260 unsigned int iter = 0;
263 conv = this->iteration_status(iter, gradients * gradients, x);
268 const auto &maximum_timestep = additional_data.maximum_timestep;
269 double timestep = additional_data.initial_timestep;
272 double alpha = ALPHA_0;
274 unsigned int previous_iter_with_positive_v_dot_g = 0;
280 x.add(timestep, velocities);
281 inverse_mass_matrix.vmult(gradients, gradients);
282 velocities.add(-timestep, gradients);
285 compute(gradients, x);
288 conv = this->iteration_status(iter, gradient_norm_squared, x);
293 const real_type v_dot_g = velocities *
gradients;
297 const real_type velocities_norm_squared = velocities * velocities;
303 const real_type beta =
304 -alpha *
std::sqrt(velocities_norm_squared / gradient_norm_squared);
307 velocities.sadd(1. - alpha, beta, gradients);
309 if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
312 timestep =
std::min(timestep * TIMESTEP_GROW, maximum_timestep);
313 alpha *= ALPHA_SHRINK;
319 previous_iter_with_positive_v_dot_g = iter;
320 timestep *= TIMESTEP_SHRINK;
325 real_type vmax = velocities.linfty_norm();
330 const double minimal_timestep =
331 additional_data.maximum_linfty_norm / vmax;
332 if (minimal_timestep < timestep)
333 timestep = minimal_timestep;
336 print_vectors(iter, x, velocities, gradients);
348template <
typename VectorType>
349template <
typename MatrixType,
typename PreconditionerType>
354 const PreconditionerType &preconditioner)
356 std::function<double(VectorType &,
const VectorType &)> compute_func =
357 [&](VectorType &g,
const VectorType &x) ->
double {
366 return 0.5 * A.matrix_norm_square(x) - x *
b;
369 this->solve(compute_func, x, preconditioner);
374template <
typename VectorType>
379 const VectorType &)
const
@ 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
#define DEAL_II_NAMESPACE_CLOSE
#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_timestep
const double initial_timestep
const double maximum_linfty_norm
AdditionalData(const double initial_timestep=0.1, const double maximum_timestep=1, const double maximum_linfty_norm=1)