16 #ifndef dealii_solver_fire_h 17 #define dealii_solver_fire_h 20 #include <deal.II/base/config.h> 22 #include <deal.II/base/logstream.h> 24 #include <deal.II/lac/diagonal_matrix.h> 25 #include <deal.II/lac/solver.h> 30 DEAL_II_NAMESPACE_OPEN
90 template <
typename VectorType = Vector<
double>>
151 template <
typename PreconditionerType = DiagonalMatrix<VectorType>>
153 solve(
const std::function<
double(VectorType &,
const VectorType &)> &compute,
155 const PreconditionerType &inverse_mass_matrix);
162 template <
typename MatrixType,
typename PreconditionerType>
164 solve(
const MatrixType & A,
166 const VectorType & b,
167 const PreconditionerType &preconditioner);
180 const VectorType &g)
const;
194 template <
typename VectorType>
199 : initial_timestep(initial_timestep)
200 , maximum_timestep(maximum_timestep)
201 , maximum_linfty_norm(maximum_linfty_norm)
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."));
212 template <
typename VectorType>
222 template <
typename VectorType>
231 template <
typename VectorType>
237 template <
typename VectorType>
238 template <
typename PreconditionerType>
241 const std::function<
double(VectorType &,
const VectorType &)> &compute,
243 const PreconditionerType &inverse_mass_matrix)
248 const double DELAYSTEP = 5;
249 const double TIMESTEP_GROW = 1.1;
250 const double TIMESTEP_SHRINK = 0.5;
251 const double ALPHA_0 = 0.1;
252 const double ALPHA_SHRINK = 0.99;
254 using real_type =
typename VectorType::real_type;
265 VectorType &velocities = *v;
266 VectorType &gradients = *g;
269 compute(gradients, x);
271 unsigned int iter = 0;
283 double alpha = ALPHA_0;
285 unsigned int previous_iter_with_positive_v_dot_g = 0;
291 x.add(timestep, velocities);
292 inverse_mass_matrix.vmult(gradients, gradients);
293 velocities.add(-timestep, gradients);
296 compute(gradients, x);
298 const real_type gradient_norm_squared = gradients * gradients;
304 const real_type v_dot_g = velocities * gradients;
308 const real_type velocities_norm_squared = velocities * velocities;
314 const real_type beta =
315 -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
318 velocities.sadd(1. - alpha, beta, gradients);
320 if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
323 timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
324 alpha *= ALPHA_SHRINK;
330 previous_iter_with_positive_v_dot_g = iter;
331 timestep *= TIMESTEP_SHRINK;
336 real_type vmax = velocities.linfty_norm();
341 const double minimal_timestep =
343 if (minimal_timestep < timestep)
344 timestep = minimal_timestep;
359 template <
typename VectorType>
360 template <
typename MatrixType,
typename PreconditionerType>
364 const VectorType & b,
365 const PreconditionerType &preconditioner)
367 std::function<double(VectorType &, const VectorType &)> compute_func =
368 [&](VectorType &g,
const VectorType &x) ->
double {
377 return 0.5 * A.matrix_norm_square(x) - x * b;
380 this->
solve(compute_func, x, preconditioner);
385 template <
typename VectorType>
390 const VectorType &)
const 397 DEAL_II_NAMESPACE_CLOSE
AdditionalData(const double initial_timestep=0.1, const double maximum_timestep=1, const double maximum_linfty_norm=1)
const AdditionalData additional_data
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcMessage(std::string arg1)
Stop iteration, goal reached.
const double maximum_timestep
#define Assert(cond, exc)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType ¤t_iterate), StateCombiner > iteration_status
void solve(const std::function< double(VectorType &, const VectorType &)> &compute, VectorType &x, const PreconditionerType &inverse_mass_matrix)
virtual void print_vectors(const unsigned int, const VectorType &x, const VectorType &v, const VectorType &g) const
VectorMemory< VectorType > & memory
const double initial_timestep
const double maximum_linfty_norm
SolverFIRE(SolverControl &solver_control, VectorMemory< VectorType > &vector_memory, const AdditionalData &data=AdditionalData())
static::ExceptionBase & ExcInternalError()