CouetteSolver.h

// Copyright (C) 2015 Technische Universitaet Muenchen
// This file is part of the Mamico project. For conditions of distribution
// and use, please see the copyright notice in Mamico's main folder, or at
// www5.in.tum.de/mamico
#ifndef _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_
#define _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_
 
#include "tarch/la/Vector.h"
#include <cmath>
 
namespace coupling {
namespace solvers {
template <unsigned int dim> class AbstractCouetteSolver {
public:
  virtual ~AbstractCouetteSolver() {}
  virtual void advance(double dt) = 0;
  virtual tarch::la::Vector<dim, double> getVelocity(tarch::la::Vector<dim, double> pos) const = 0;
  virtual void setWallVelocity(const tarch::la::Vector<dim, double> wallVelocity) = 0;
};
 
template <unsigned int dim> class CouetteSolver;
} // namespace solvers
} // namespace coupling

template <unsigned int dim> class coupling::solvers::CouetteSolver : public coupling::solvers::AbstractCouetteSolver<dim> {
public:
  CouetteSolver(const double& channelheight, const double& wallVelocity, const double kinVisc)
      : AbstractCouetteSolver<dim>(), _channelheight(channelheight), _wallVelocity(wallVelocity), _kinVisc(kinVisc), _time(0.0) {}

  virtual ~CouetteSolver() {}

  virtual void advance(double dt) { _time += dt; }

  virtual tarch::la::Vector<dim, double> getVelocity(tarch::la::Vector<dim, double> pos) const {
    tarch::la::Vector<dim, double> v(0.0);
    v[0] = _wallVelocity * (1.0 - pos[dim - 1] / _channelheight);
    const double pi = 3.141592653589793238;
    double sum = 0.0;
    for (int k = 1; k < 30; k++) {
      sum += 1.0 / k * sin(k * pi * pos[dim - 1] / _channelheight) * exp(-k * k * pi * pi / (_channelheight * _channelheight) * _kinVisc * _time);
    }
    v[0] = v[0] - 2.0 * _wallVelocity / pi * sum;
    return v;
  }

  virtual void setWallVelocity(const tarch::la::Vector<dim, double> wallVelocity) { _wallVelocity = wallVelocity[0]; }
 
private:
  const double _channelheight;
  double _wallVelocity;
  const double _kinVisc;
  double _time;
};
 
#endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_