timestepping.cc

template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme {
public:
  TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
                     VectorType const &_u_old,
                     Dune::BitSetVector<dim> const &_dirichletNodes,
                     FunctionType const &_dirichletFunction, double _time,
                     double _tau)
      : ell(_ell),
        A(_A),
        u_old(_u_old),
        dirichletNodes(_dirichletNodes),
        dirichletFunction(_dirichletFunction),
        time(_time),
        tau(_tau) {}

  virtual ~TimeSteppingScheme() {}

  void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                     MatrixType &problem_A) const = 0;

  void virtual extractSolution(VectorType const &problem_iterate,
                               VectorType &solution) const = 0;

  void virtual extractVelocity(VectorType const &problem_iterate,
                               VectorType &velocity) const = 0;

protected:
  VectorType const &ell;
  MatrixType const &A;
  VectorType const &u_old;
  Dune::BitSetVector<dim> const &dirichletNodes;
  FunctionType const &dirichletFunction;
  double const time;
  double const tau;
};

template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitEuler
    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
  ImplicitEuler(
      VectorType const &_ell, MatrixType const &_A, VectorType const &_u_old,
      VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
      Dune::BitSetVector<dim> const &_dirichletNodes,
      FunctionType const &_dirichletFunction, double _time, double _tau)
      : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
            _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
        u_old_old_ptr(_u_old_old_ptr) {}

  void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                     MatrixType &problem_A) const {
    problem_rhs = this->ell;
    this->A.mmv(this->u_old, problem_rhs);

    problem_A = this->A;
    problem_A *= this->tau;

    // Use the old velocity as an initial iterate if possible
    if (u_old_old_ptr) {
      problem_iterate = this->u_old;
      problem_iterate -= *u_old_old_ptr;
      problem_iterate /= this->tau;
    } else
      problem_iterate = 0.0;

    for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
      if (this->dirichletNodes[i].count() == dim) {
        problem_iterate[i] = 0;
        this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
      } else if (this->dirichletNodes[i][1])
        problem_iterate[i][1] = 0; // Y direction prescribed
  }

  void virtual extractSolution(VectorType const &problem_iterate,
                               VectorType &solution) const {
    solution = this->u_old;
    solution.axpy(this->tau, problem_iterate);
  }

  void virtual extractVelocity(VectorType const &problem_iterate,
                               VectorType &velocity) const {
    velocity = problem_iterate;
  }

private:
  VectorType const *u_old_old_ptr;
};

template <class VectorType, class MatrixType, class FunctionType, int dim>
class ImplicitTwoStep
    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
  ImplicitTwoStep(VectorType const &_ell, MatrixType const &_A,
                  VectorType const &_u_old, VectorType const *_u_old_old_ptr,
                  Dune::BitSetVector<dim> const &_dirichletNodes,
                  FunctionType const &_dirichletFunction, double _time,
                  double _tau)
      : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
            _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
        u_old_old_ptr(_u_old_old_ptr) {}

  void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
                     MatrixType &problem_A) const {
    problem_rhs = this->ell;
    this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
    this->A.usmv(+1.0 / 3.0, *u_old_old_ptr, problem_rhs);

    problem_A = this->A;
    problem_A *= 2.0 / 3.0 * this->tau;

    // Use the old velocity as an initial iterate
    problem_iterate = this->u_old;
    problem_iterate -= *u_old_old_ptr;
    problem_iterate /= this->tau;

    for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
      if (this->dirichletNodes[i].count() == dim) {
        problem_iterate[i] = 0;
        this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
      } else if (this->dirichletNodes[i][1])
        problem_iterate[i][1] = 0; // Y direction prescribed
  }

  void virtual extractSolution(VectorType const &problem_iterate,
                               VectorType &solution) const {
    solution = problem_iterate;
    solution *= this->tau;
    solution.axpy(2, this->u_old);
    solution.axpy(-.5, *u_old_old_ptr);
    solution *= 2.0 / 3.0;

    // Check if we split correctly
    {
      VectorType test = problem_iterate;
      test *= this->tau;
      test.axpy(-1.5, solution);
      test.axpy(+2, this->u_old);
      test.axpy(-.5, *u_old_old_ptr);
      assert(test.two_norm() < 1e-10);
    }
  }

  void virtual extractVelocity(VectorType const &problem_iterate,
                               VectorType &velocity) const {
    velocity = problem_iterate;
  }

private:
  VectorType const *u_old_old_ptr;
};