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_