Line data    Source code

       1             : // Copyright (C) 2015 Technische Universitaet Muenchen
       2             : // This file is part of the Mamico project. For conditions of distribution
       3             : // and use, please see the copyright notice in Mamico's main folder, or at
       4             : // www5.in.tum.de/mamico
       5             : #ifndef _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_
       6             : #define _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_
       7             : 
       8             : #include "tarch/la/Vector.h"
       9             : #include <cmath>
      10             : 
      11             : namespace coupling {
      12             : /** @brief all numerical solvers are defined in the namespace, and their
      13             :  * interfaces  */
      14             : namespace solvers {
      15             : /** @brief interface for continuum/macro fluid solvers for the Couette scenario
      16             :  *  @author Philipp Neumann
      17             :  *  @tparam dim  refers to the spacial dimension of the simulation, can be 1, 2,
      18             :  * or 3*/
      19         116 : template <unsigned int dim> class AbstractCouetteSolver {
      20             : public:
      21             :   /** @brief a dummy destructor
      22             :    *  @todo Why a dummy destructor buut not a constructor */
      23         116 :   virtual ~AbstractCouetteSolver() {}
      24             :   /** @brief advances the solver in time
      25             :    *  @param dt the solver will be advanced from the current time t to t+dt */
      26             :   virtual void advance(double dt) = 0;
      27             :   /** @brief returns the current velocity at the given position
      28             :    *  @param pos position to return the velocity for
      29             :    *  @returns a velocity vector */
      30             :   virtual tarch::la::Vector<dim, double> getVelocity(tarch::la::Vector<dim, double> pos) const = 0;
      31             :   /** @brief changes the velocity at the moving for, refers to Couette scenario
      32             :    *  @param wallVelocity value to which the veloctiy will be set */
      33             :   virtual void setWallVelocity(const tarch::la::Vector<dim, double> wallVelocity) = 0;
      34             : };
      35             : 
      36             : template <unsigned int dim> class CouetteSolver;
      37             : } // namespace solvers
      38             : } // namespace coupling
      39             : 
      40             : /** In our scenario, the lower wall is accelerated and the upper wall stands
      41             :  * still. The lower wall is located at zero height.
      42             :  *  @brief implements an analytic Couette flow solver.
      43             :  *  @author Philipp Neumann
      44             :  *  @tparam dim  refers to the spacial dimension of the simulation, can be 1, 2,
      45             :  * or 3 */
      46             : template <unsigned int dim> class coupling::solvers::CouetteSolver : public coupling::solvers::AbstractCouetteSolver<dim> {
      47             : public:
      48             :   /** @brief a simple constructor
      49             :    *  @param channelheight the height and width of the channel in y and z
      50             :    * direction
      51             :    *  @param wallVelocity the velocity at the moving wall
      52             :    *  @param kinVisc the kinematic viscosity of the fluid  */
      53             :   CouetteSolver(const double& channelheight, const double& wallVelocity, const double kinVisc)
      54             :       : AbstractCouetteSolver<dim>(), _channelheight(channelheight), _wallVelocity(wallVelocity), _kinVisc(kinVisc), _time(0.0) {}
      55             : 
      56             :   /** @brief a dummy destructor */
      57             :   virtual ~CouetteSolver() {}
      58             : 
      59             :   /** @brief advances one time step dt in time
      60             :    *  @param dt size of the time step to advance */
      61             :   virtual void advance(double dt) { _time += dt; }
      62             : 
      63             :   /** for the first entry of the vector is the analytic solution of the Couette
      64             :    * scenario returned, given by: u(z,t)= V(1-z/H) - 2V/pi*sum_k=1^infty
      65             :    * 1/k*sin(k*pi*z/H)*exp(-k^2 pi^2/H^2 * nu*t) The other two entries are 0
      66             :    *  @brief returns the velocity vector at a certain channel position
      67             :    *  @param pos the position to return the velocity for
      68             :    *  @returns a velocity vector */
      69             :   virtual tarch::la::Vector<dim, double> getVelocity(tarch::la::Vector<dim, double> pos) const {
      70             :     tarch::la::Vector<dim, double> v(0.0);
      71             :     v[0] = _wallVelocity * (1.0 - pos[dim - 1] / _channelheight);
      72             :     const double pi = 3.141592653589793238;
      73             :     double sum = 0.0;
      74             :     for (int k = 1; k < 30; k++) {
      75             :       sum += 1.0 / k * sin(k * pi * pos[dim - 1] / _channelheight) * exp(-k * k * pi * pi / (_channelheight * _channelheight) * _kinVisc * _time);
      76             :     }
      77             :     v[0] = v[0] - 2.0 * _wallVelocity / pi * sum;
      78             :     return v;
      79             :   }
      80             : 
      81             :   /** @brief changes the velocity at the moving for, refers to Couette scenario
      82             :    *  @param wallVelocity value to which the veloctiy will be set */
      83             :   virtual void setWallVelocity(const tarch::la::Vector<dim, double> wallVelocity) { _wallVelocity = wallVelocity[0]; }
      84             : 
      85             : private:
      86             :   /** @brief height of couette channel */
      87             :   const double _channelheight;
      88             :   /** @brief velocity of moving wall */
      89             :   double _wallVelocity;
      90             :   /** @brief kinematic viscosity */
      91             :   const double _kinVisc;
      92             :   /** @brief current time */
      93             :   double _time;
      94             : };
      95             : 
      96             : #endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_COUETTESOLVER_H_