time_dependent.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_time_dependent_h
#  define dealii_time_dependent_h
 
 
/*----------------------------   time-dependent.h ---------------------------*/
 
 
#  include <deal.II/base/config.h>
 
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/smartpointer.h>
#  include <deal.II/base/subscriptor.h>
 
#  include <utility>
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
class TimeStepBase;
template <typename number>
class Vector;
template <int dim, int spacedim>
class Triangulation;
 
class TimeDependent
{
public:
  struct TimeSteppingData
  {
    TimeSteppingData(const unsigned int look_ahead,
                     const unsigned int look_back);
 
    const unsigned int look_ahead;
 
    const unsigned int look_back;
  };
 
  enum Direction
  {
    forward,
    backward
  };
 
  TimeDependent(const TimeSteppingData &data_primal,
                const TimeSteppingData &data_dual,
                const TimeSteppingData &data_postprocess);
 
 
  virtual ~TimeDependent();
 
  void
  insert_timestep(const TimeStepBase *position, TimeStepBase *new_timestep);
 
  void
  add_timestep(TimeStepBase *new_timestep);
 
  void
  delete_timestep(const unsigned int position);
 
  void
  solve_primal_problem();
 
  void
  solve_dual_problem();
 
  void
  postprocess();
 
  template <typename InitFunctionObject, typename LoopFunctionObject>
  void
  do_loop(InitFunctionObject      init_function,
          LoopFunctionObject      loop_function,
          const TimeSteppingData &timestepping_data,
          const Direction         direction);
 
 
  virtual void
  start_sweep(const unsigned int sweep_no);
 
  virtual void
  end_sweep();
 
  std::size_t
  memory_consumption() const;
 
  DeclExceptionMsg(ExcInvalidPosition,
                   "You cannot insert a time step at the specified position.");
 
protected:
  std::vector<SmartPointer<TimeStepBase, TimeDependent>> timesteps;
 
  unsigned int sweep_no;
 
  const TimeSteppingData timestepping_data_primal;
 
  const TimeSteppingData timestepping_data_dual;
 
  const TimeSteppingData timestepping_data_postprocess;
 
private:
  void
  end_sweep(const unsigned int begin_timestep, const unsigned int end_timestep);
};
 
 
 
class TimeStepBase : public Subscriptor
{
public:
  enum SolutionState
  {
    primal_problem = 0x0,
    dual_problem = 0x1,
    postprocess = 0x2
  };
 
  TimeStepBase(const double time);
 
  virtual ~TimeStepBase() override = default;
 
  TimeStepBase(const TimeStepBase &) = delete;
 
  TimeStepBase &
  operator=(const TimeStepBase &) = delete;
 
  virtual void
  wake_up(const unsigned int);
 
  virtual void
  sleep(const unsigned int);
 
  virtual void
  start_sweep();
 
  virtual void
  end_sweep();
 
  virtual void
  init_for_primal_problem();
 
  virtual void
  init_for_dual_problem();
 
  virtual void
  init_for_postprocessing();
 
  virtual void
  solve_primal_problem() = 0;
 
  virtual void
  solve_dual_problem();
 
  virtual void
  postprocess_timestep();
 
  double
  get_time() const;
 
  unsigned int
  get_timestep_no() const;
 
  double
  get_backward_timestep() const;
 
  double
  get_forward_timestep() const;
 
  virtual std::size_t
  memory_consumption() const;
 
protected:
  const TimeStepBase *previous_timestep;
 
  const TimeStepBase *next_timestep;
 
  unsigned int sweep_no;
 
  unsigned int timestep_no;
 
  const double time;
 
  unsigned int next_action;
 
private:
  void
  set_previous_timestep(const TimeStepBase *previous);
 
  void
  set_next_timestep(const TimeStepBase *next);
 
  void
  set_timestep_no(const unsigned int step_no);
 
  void
  set_sweep_no(const unsigned int sweep_no);
 
  // make the manager object a friend
  friend class TimeDependent;
};
 
 
 
namespace TimeStepBase_Tria_Flags
{
  template <int dim>
  struct Flags
  {
    Flags();
 
    Flags(const bool         delete_and_rebuild_tria,
          const unsigned int wakeup_level_to_build_grid,
          const unsigned int sleep_level_to_delete_grid);
 
    const bool delete_and_rebuild_tria;
 
    const unsigned int wakeup_level_to_build_grid;
 
    const unsigned int sleep_level_to_delete_grid;
  };
 
 
 
  template <int dim>
  struct RefinementFlags
  {
    using CorrectionRelaxations =
      std::vector<std::vector<std::pair<unsigned int, double>>>;
 
    static CorrectionRelaxations default_correction_relaxations;
 
    RefinementFlags(const unsigned int max_refinement_level        = 0,
                    const unsigned int first_sweep_with_correction = 0,
                    const unsigned int min_cells_for_correction    = 0,
                    const double       cell_number_corridor_top    = (1 << dim),
                    const double       cell_number_corridor_bottom = 1,
                    const CorrectionRelaxations &correction_relaxations =
                      CorrectionRelaxations(),
                    const unsigned int cell_number_correction_steps  = 0,
                    const bool         mirror_flags_to_previous_grid = false,
                    const bool         adapt_grids                   = false);
 
    const unsigned int max_refinement_level;
 
    const unsigned int first_sweep_with_correction;
 
 
    const unsigned int min_cells_for_correction;
 
    const double cell_number_corridor_top;
 
    const double cell_number_corridor_bottom;
 
    const std::vector<std::vector<std::pair<unsigned int, double>>>
      correction_relaxations;
 
    const unsigned int cell_number_correction_steps;
 
    const bool mirror_flags_to_previous_grid;
 
    const bool adapt_grids;
 
    DeclException1(ExcInvalidValue,
                   double,
                   << "The value " << arg1
                   << " for the cell number corridor does not fulfill "
                      "its natural requirements.");
  };
 
 
 
  template <int dim>
  struct RefinementData
  {
    RefinementData(const double refinement_threshold,
                   const double coarsening_threshold = 0);
 
    const double refinement_threshold;
 
    const double coarsening_threshold;
 
    DeclException1(ExcInvalidValue,
                   double,
                   << "The value " << arg1
                   << " for the cell refinement thresholds does not fulfill "
                      "its natural requirements.");
  };
} // namespace TimeStepBase_Tria_Flags
 
 
 
template <int dim>
class TimeStepBase_Tria : public TimeStepBase
{
public:
  using Flags = typename TimeStepBase_Tria_Flags::Flags<dim>;
  using RefinementFlags =
    typename TimeStepBase_Tria_Flags::RefinementFlags<dim>;
  using RefinementData = typename TimeStepBase_Tria_Flags::RefinementData<dim>;
 
 
  enum SolutionState
  {
    grid_refinement = 0x1000
  };
 
 
  TimeStepBase_Tria();
 
  TimeStepBase_Tria(
    const double                   time,
    const Triangulation<dim, dim> &coarse_grid,
    const Flags &                  flags,
    const RefinementFlags &        refinement_flags = RefinementFlags());
 
  virtual ~TimeStepBase_Tria() override;
 
  virtual void
  wake_up(const unsigned int wakeup_level) override;
 
  virtual void
  sleep(const unsigned int) override;
 
  void
  refine_grid(const RefinementData data);
 
  virtual void
  init_for_refinement();
 
  virtual void
  get_tria_refinement_criteria(Vector<float> &criteria) const = 0;
 
  void
  save_refine_flags();
 
  virtual std::size_t
  memory_consumption() const override;
 
  DeclExceptionMsg(ExcGridNotDeleted,
                   "When calling restore_grid(), you must have previously "
                   "deleted the triangulation.");
 
protected:
  SmartPointer<Triangulation<dim, dim>, TimeStepBase_Tria<dim>> tria;
 
  SmartPointer<const Triangulation<dim, dim>, TimeStepBase_Tria<dim>>
    coarse_grid;
 
  const Flags flags;
 
  const RefinementFlags refinement_flags;
 
private:
  std::vector<std::vector<bool>> refine_flags;
 
  std::vector<std::vector<bool>> coarsen_flags;
 
  void
  restore_grid();
};
 
 
 
/*--------------------------- template functions ----------------------------*/
 
template <typename InitFunctionObject, typename LoopFunctionObject>
void
TimeDependent::do_loop(InitFunctionObject      init_function,
                       LoopFunctionObject      loop_function,
                       const TimeSteppingData &timestepping_data,
                       const Direction         direction)
{
  // the following functions looks quite
  // disrupted due to the recurring switches
  // for forward and backward running loops.
  //
  // I chose to switch at every place where
  // it is needed, since it is so easy
  // to overlook something when you change
  // some code at one place when it needs
  // to be changed at a second place, here
  // for the other direction, also.
 
  const unsigned int n_timesteps = timesteps.size();
 
  // initialize the time steps for
  // a round of this loop
  for (unsigned int step = 0; step < n_timesteps; ++step)
    switch (direction)
      {
        case forward:
          init_function((&*timesteps[step]));
          break;
        case backward:
          init_function((&*timesteps[n_timesteps - step - 1]));
          break;
      };
 
 
  // wake up the first few time levels
  for (int step = -static_cast<int>(timestepping_data.look_ahead); step < 0;
       ++step)
    for (int look_ahead = 0;
         look_ahead <= static_cast<int>(timestepping_data.look_ahead);
         ++look_ahead)
      switch (direction)
        {
          case forward:
            if (step + look_ahead >= 0)
              timesteps[step + look_ahead]->wake_up(look_ahead);
            break;
          case backward:
            if (n_timesteps - (step + look_ahead) < n_timesteps)
              timesteps[n_timesteps - (step + look_ahead)]->wake_up(look_ahead);
            break;
        };
 
 
  for (unsigned int step = 0; step < n_timesteps; ++step)
    {
      // first thing: wake up the
      // timesteps ahead as necessary
      for (unsigned int look_ahead = 0;
           look_ahead <= timestepping_data.look_ahead;
           ++look_ahead)
        switch (direction)
          {
            case forward:
              if (step + look_ahead < n_timesteps)
                timesteps[step + look_ahead]->wake_up(look_ahead);
              break;
            case backward:
              if (n_timesteps > (step + look_ahead))
                timesteps[n_timesteps - (step + look_ahead) - 1]->wake_up(
                  look_ahead);
              break;
          };
 
 
      // actually do the work
      switch (direction)
        {
          case forward:
            loop_function((&*timesteps[step]));
            break;
          case backward:
            loop_function((&*timesteps[n_timesteps - step - 1]));
            break;
        };
 
      // let the timesteps behind sleep
      for (unsigned int look_back = 0; look_back <= timestepping_data.look_back;
           ++look_back)
        switch (direction)
          {
            case forward:
              if (step >= look_back)
                timesteps[step - look_back]->sleep(look_back);
              break;
            case backward:
              if (n_timesteps - (step - look_back) <= n_timesteps)
                timesteps[n_timesteps - (step - look_back) - 1]->sleep(
                  look_back);
              break;
          }
    }
 
  // make the last few timesteps sleep
  for (int step = n_timesteps;
       step < static_cast<int>(n_timesteps + timestepping_data.look_back);
       ++step)
    for (int look_back = 0;
         look_back <= static_cast<int>(timestepping_data.look_back);
         ++look_back)
      switch (direction)
        {
          case forward:
            if ((step - look_back >= 0) &&
                (step - look_back < static_cast<int>(n_timesteps)))
              timesteps[step - look_back]->sleep(look_back);
            break;
          case backward:
            if ((step - look_back >= 0) &&
                (step - look_back < static_cast<int>(n_timesteps)))
              timesteps[n_timesteps - (step - look_back) - 1]->sleep(look_back);
            break;
        };
}
 
DEAL_II_NAMESPACE_CLOSE
 
/*----------------------------   time-dependent.h ---------------------------*/
#endif
/*----------------------------   time-dependent.h ---------------------------*/