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           0 :   NieVelocityImposition(coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface, const unsigned int& outermostLayer,
      32             :                         const unsigned int& innermostLayer)
      33           0 :       : coupling::MomentumInsertion<LinkedCell, dim>(mdSolverInterface), _outermostLayer(outermostLayer), _innermostLayer(innermostLayer) {}
      34             : 
      35             :   /** @brief a simple destructor */
      36           0 :   virtual ~NieVelocityImposition() {}
      37             : 
      38             :   /** @brief momentum shall be inserted in every md time step, so this returns 1
      39             :    *  @returns the time step interval for momentum insertion, always 1 */
      40           0 :   unsigned int getTimeIntervalPerMomentumInsertion() const override { return 1; }
      41             : 
      42             :   /** @brief inserts momentum to a cell
      43             :    *  @param cell to the coupling cell will the momentum be inserted
      44             :    *  @param idx local linearised index for the
      45             :    * coupling cell */
      46           0 :   void insertMomentum(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 idx) const override {
      47             :     // nop if this is not an imposition cell
      48           0 :     if (!isInsideImpositionLayer(idx)) {
      49           0 :       return;
      50             :     }
      51             :     // set continuum velocity
      52           0 :     tarch::la::Vector<dim, double> continuumVelocity(cell.getMicroscopicMomentum());
      53             : 
      54           0 :     coupling::cellmappings::ComputeAvgForceAndVelocity<LinkedCell, dim> computeForceAndVelocity(
      55           0 :         coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
      56           0 :     cell.iterateConstCells(computeForceAndVelocity);
      57           0 :     const tarch::la::Vector<dim, double> avgVel(computeForceAndVelocity.getAvgVelocity());
      58           0 :     const tarch::la::Vector<dim, double> avgF(computeForceAndVelocity.getAvgForce());
      59             : 
      60           0 :     coupling::cellmappings::NieVelocityImpositionMapping<LinkedCell, dim> velocityImposition(continuumVelocity, avgVel, avgF,
      61           0 :                                                                                              coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
      62           0 :     cell.iterateCells(velocityImposition);
      63           0 :   }
      64             : 
      65             : private:
      66             :   /** returns true if the local cell at index currentLocalCouplingCellIndex is
      67             :    * inside the layer of imposition cells, given by outermostLayer and
      68             :    * innermostLayer. For, e.g., outermostLayer=2 and innermostLayer=3, the
      69             :    * layers for imposition are located in the 3rd and 4th strip of cells (we
      70             :    * start counting from cell layer=0 which corresponds to the outermost,
      71             :    * actually ghost-layer of cells which surrounds the MD domain).
      72             :    *  @brief based on the cell index, the function tells if the cell is inside
      73             :    * the imposition layer
      74             :    *  @param globalCellIndex global linearised index of a coupling
      75             :    * cell to check
      76             :    *  @returns a bool, that indicates if the given cell index is located in the
      77             :    * imposition layer (true) or not (false) */
      78           0 :   bool isInsideImpositionLayer(I01 globalCellIndex) const {
      79           0 :     bool inner = true;
      80           0 :     tarch::la::Vector<dim, unsigned int> globalIndexUnsigned{globalCellIndex.get()};
      81           0 :     for (unsigned int d = 0; d < dim; d++)
      82           0 :       inner = inner && (globalIndexUnsigned[d] > _innermostLayer && globalIndexUnsigned[d] < 1 + I09::numberCellsInDomain[d] - _innermostLayer);
      83             :     bool outer = false;
      84           0 :     for (unsigned int d = 0; d < dim; d++)
      85           0 :       outer = outer || (globalIndexUnsigned[d] < _outermostLayer || globalIndexUnsigned[d] > 1 + I09::numberCellsInDomain[d] - _outermostLayer);
      86           0 :     return !inner && !outer;
      87             :   }
      88             : 
      89             :   /** @brief the index of the outermost cell layer*/
      90             :   const unsigned int _outermostLayer;
      91             :   /** @brief the index of the innermost cell layer*/
      92             :   const unsigned int _innermostLayer;
      93             : };
      94             : 
      95             : #endif // _MOLECULARDYNAMICS_COUPLING_VELOCITYGRADIENTRELAXATION_H_