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_