time_dependent.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2018 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;

  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);


  TimeStepBase(const TimeStepBase &) = delete;

  TimeStepBase &
  operator=(const TimeStepBase &) = delete;

  // 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 = -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 ---------------------------*/