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_NIEVELOCITYIMPOSITION_H_
       6             : #define _MOLECULARDYNAMICS_COUPLING_NIEVELOCITYIMPOSITION_H_
       7             : 
       8             : #include "coupling/MomentumInsertion.h"
       9             : #include "coupling/cell-mappings/ComputeAvgForceAndVelocity.h"
      10             : #include "coupling/cell-mappings/NieVelocityImpositionMapping.h"
      11             : #include "coupling/interface/MDSolverInterface.h"
      12             : 
      13             : namespace coupling {
      14             : template <class LinkedCell, unsigned int dim> class NieVelocityImposition;
      15             : }
      16             : 
      17             : /** @brief Velocity imposition scheme following the respective paper by Nie et
      18             :  * al., J. Fluid. Mech. 500, 2004
      19             :  *  @author Philipp Neumann
      20             :  *  @tparam LinkedCell the LinkedCell class is given by the implementation of
      21             :  * linked cells in the molecular dynamics simulation
      22             :  *  @tparam dim  refers to the spacial dimension of the simulation, can be 1, 2,
      23             :  * or 3
      24             :  */
      25             : template <class LinkedCell, unsigned int dim> class coupling::NieVelocityImposition : public coupling::MomentumInsertion<LinkedCell, dim> {
      26             : public:
      27             :   /** @brief a simple constructor
      28             :    *  @param mdSolverInterface interface to the md solver
      29             :    *  @param outermostLayer the index of the outermost cell layer
      30             :    *  @param innermostLayer the index of the innermost cell layer
      31             :    *  @param impositionEnabled whether imposition is allowed at the specified domain face */
      32           0 :   NieVelocityImposition(coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface, const unsigned int& outermostLayer,
      33             :                         const unsigned int& innermostLayer, const tarch::la::Vector<2 * dim, bool> impositionEnabled)
      34           0 :       : coupling::MomentumInsertion<LinkedCell, dim>(mdSolverInterface), _outermostLayer(outermostLayer), _innermostLayer(innermostLayer),
      35           0 :         _impositionEnabled(impositionEnabled), _enableInnerImposition(false) {}
      36             : 
      37             :   /** @brief a simple destructor */
      38           0 :   virtual ~NieVelocityImposition() {}
      39             : 
      40             :   /** @brief momentum shall be inserted in every md time step, so this returns 1
      41             :    *  @returns the time step interval for momentum insertion, always 1 */
      42           0 :   unsigned int getTimeIntervalPerMomentumInsertion() const override { return 1; }
      43             : 
      44             :   /** @brief inserts momentum to a cell
      45             :    *  @param cell to the coupling cell will the momentum be inserted
      46             :    *  @param idx local linearised index for the
      47             :    * coupling cell */
      48           0 :   void insertMomentum(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 idx) const override {
      49             :     // nop if this is not an imposition cell
      50           0 :     if (!isInsideImpositionLayer(idx)) {
      51           0 :       return;
      52             :     }
      53             :     // set continuum velocity
      54           0 :     tarch::la::Vector<dim, double> continuumVelocity(cell.getMicroscopicMomentum());
      55             : 
      56           0 :     coupling::cellmappings::ComputeAvgForceAndVelocity<LinkedCell, dim> computeForceAndVelocity(
      57           0 :         coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
      58           0 :     cell.iterateConstCells(computeForceAndVelocity);
      59           0 :     const tarch::la::Vector<dim, double> avgVel(computeForceAndVelocity.getAvgVelocity());
      60           0 :     const tarch::la::Vector<dim, double> avgF(computeForceAndVelocity.getAvgForce());
      61             : 
      62           0 :     coupling::cellmappings::NieVelocityImpositionMapping<LinkedCell, dim> velocityImposition(continuumVelocity, avgVel, avgF,
      63           0 :                                                                                              coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
      64           0 :     cell.iterateCells(velocityImposition);
      65           0 :   }
      66             : 
      67           0 :   void setInnerImposition(bool enable) override { _enableInnerImposition = enable; }
      68             : 
      69             : private:
      70             :   /** returns true if the local cell at index currentLocalCouplingCellIndex is
      71             :    * inside the layer of imposition cells, given by outermostLayer and
      72             :    * innermostLayer. For, e.g., outermostLayer=2 and innermostLayer=3, the
      73             :    * layers for imposition are located in the 3rd and 4th strip of cells (we
      74             :    * start counting from cell layer=0 which corresponds to the outermost,
      75             :    * actually ghost-layer of cells which surrounds the MD domain).
      76             :    * Also checks if imposition is enabled in that direction.
      77             :    *  @brief based on the cell index, the function tells if the cell is inside
      78             :    * the imposition layer
      79             :    *  @param globalCellIndex global linearised index of a coupling
      80             :    * cell to check
      81             :    *  @returns a bool, that indicates if the given cell index is located in the
      82             :    * imposition layer (true) or not (false) */
      83           0 :   bool isInsideImpositionLayer(I01 globalCellIndex) const {
      84           0 :     bool inner = true;
      85           0 :     bool outer = false;
      86           0 :     tarch::la::Vector<dim, unsigned int> globalIndexUnsigned{globalCellIndex.get()};
      87           0 :     for (unsigned int d = 0; d < dim; d++) {
      88           0 :       if (_impositionEnabled[d * 2]) { // negative direction
      89           0 :         inner = inner && (globalIndexUnsigned[d] > _innermostLayer);
      90           0 :         outer = outer || (globalIndexUnsigned[d] < _outermostLayer);
      91             :       }
      92           0 :       if (_impositionEnabled[d * 2 + 1]) { // positive direction
      93           0 :         inner = inner && (globalIndexUnsigned[d] < 1 + I09::numberCellsInDomain[d] - _innermostLayer);
      94           0 :         outer = outer || (globalIndexUnsigned[d] > 1 + I09::numberCellsInDomain[d] - _outermostLayer);
      95             :       }
      96             :     }
      97           0 :     return (_enableInnerImposition || !inner) && !outer;
      98             :   }
      99             : 
     100             :   /** @brief the index of the outermost cell layer*/
     101             :   const unsigned int _outermostLayer;
     102             :   /** @brief the index of the innermost cell layer*/
     103             :   const unsigned int _innermostLayer;
     104             :   const tarch::la::Vector<2 * dim, bool> _impositionEnabled;
     105             :   bool _enableInnerImposition;
     106             : };
     107             : 
     108             : #endif // _MOLECULARDYNAMICS_COUPLING_VELOCITYGRADIENTRELAXATION_H_