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_KINETICENERGYCONTROLLER_H_
       6             : #define _MOLECULARDYNAMICS_COUPLING_KINETICENERGYCONTROLLER_H_
       7             : 
       8             : #include "coupling/cell-mappings/ComputeKineticEnergyMapping.h"
       9             : #include "coupling/cell-mappings/ComputeMassMapping.h"
      10             : #include "coupling/cell-mappings/ComputeMomentumMapping.h"
      11             : #include "coupling/cell-mappings/ComputeTemperatureMapping.h"
      12             : #include "coupling/cell-mappings/SetKineticEnergyMapping.h"
      13             : #include "coupling/cell-mappings/SetTemperatureMapping.h"
      14             : #include "coupling/datastructures/CouplingCell.h"
      15             : #include "coupling/datastructures/CouplingCellWithLinkedCells.h"
      16             : 
      17             : namespace coupling {
      18             : template <class LinkedCell, unsigned int dim> class KineticEnergyController;
      19             : }
      20             : 
      21             : /** @brief controles and regulates the kinetic energy of the MD system.
      22             :  *  @author Philipp Neumann
      23             :  *  @tparam LinkedCell the LinkedCell class is given by the implementation of
      24             :  * linked cells in the molecular dynamics simulation
      25             :  *  @tparam dim  refers to the spacial dimension of the simulation, can be 1, 2,
      26             :  * or 3*/
      27             : template <class LinkedCell, unsigned int dim> class coupling::KineticEnergyController {
      28             : public:
      29             :   /** @brief a simple constructor
      30             :    *  @param mdSolverInterface interface to the md solver */
      31           4 :   KineticEnergyController(coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface) : _mdSolverInterface(mdSolverInterface) {}
      32             :   /** @brief a simple destructor*/
      33           0 :   ~KineticEnergyController() {}
      34             : 
      35             :   /** @brief computes and returns the kinetic energy within a coupling cell
      36             :    *  @param cell coupling cell to compute the kinetic energy for
      37             :    *  @returns the kinetic energy in the cell */
      38             :   double computeKineticEnergy(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell) const {
      39             :     coupling::cellmappings::ComputeKineticEnergyMapping<LinkedCell, dim> computeKineticEnergyMapping(_mdSolverInterface);
      40             :     cell.iterateConstCells(computeKineticEnergyMapping);
      41             :     return computeKineticEnergyMapping.getKineticEnergy();
      42             :   }
      43             : 
      44             :   /** @brief computes and returns the temperature within a coupling cell
      45             :    *  @param cell coupling cell to compute the temperature in
      46             :    *  @returns the temperature in the cell */
      47           0 :   double computeTemperature(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell) const {
      48           0 :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      49           0 :     cell.iterateConstCells(computeMomentumMapping);
      50           0 :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      51             : 
      52           0 :     coupling::cellmappings::ComputeTemperatureMapping<LinkedCell, dim> computeTemperatureMapping(meanVelocity, _mdSolverInterface);
      53           0 :     cell.iterateConstCells(computeTemperatureMapping);
      54           0 :     return computeTemperatureMapping.getTemperature();
      55           0 :   }
      56             : 
      57             :   /** Therefore, the mean velocity is computed first. Afterwards the deviation
      58             :    * from the mean velocity are rescaled such that momentum is conserved.
      59             :    *  @brief sets the kinetic energy within a coupling cell to the input
      60             :    * value.
      61             :    *  @param cell the coupling cell to set the kinetic energy in
      62             :    *  @param kineticEnergy the value the kinetic energy shall be set to*/
      63             :   void setKineticEnergy(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const double& kineticEnergy) const {
      64             :     // determine mass, momentum and old kinetic energy
      65             :     coupling::cellmappings::ComputeMassMapping<LinkedCell, dim> computeMassMapping(_mdSolverInterface);
      66             :     cell.iterateConstCells(computeMassMapping);
      67             :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      68             :     cell.iterateConstCells(computeMomentumMapping);
      69             :     coupling::cellmappings::ComputeKineticEnergyMapping<LinkedCell, dim> computeKineticEnergyMapping(_mdSolverInterface);
      70             :     cell.iterateConstCells(computeKineticEnergyMapping);
      71             :     // set new kinetic energy
      72             :     unsigned int numberParticles = computeMassMapping.getNumberOfParticles();
      73             :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      74             :     double oldKineticEnergy = computeKineticEnergyMapping.getKineticEnergy();
      75             :     coupling::cellmappings::SetKineticEnergyMapping<LinkedCell, dim> setKineticEnergyMapping(oldKineticEnergy, kineticEnergy, numberParticles, meanVelocity);
      76             :     cell.iterateCells(setKineticEnergyMapping);
      77             :   }
      78             : 
      79             :   /** Here we just scale the deviation from the mean velocity accordingly, i.e.
      80             :    * we set: v_molecule = v_mean +
      81             :    * sqrt(temperature/current_temperature)*(v_molecule-v_mean)
      82             :    *  @brief sets the temperature within the cell to 'temperature'.
      83             :    *  @param cell the coupling cell the temperature shall be applied in
      84             :    *  @param temperature the value the temperature shall be set at*/
      85           0 :   void setTemperature(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const double& temperature) const {
      86           0 :     coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> computeMomentumMapping(_mdSolverInterface);
      87           0 :     cell.iterateConstCells(computeMomentumMapping);
      88           0 :     tarch::la::Vector<dim, double> meanVelocity = computeMomentumMapping.getMeanVelocity();
      89             : 
      90           0 :     coupling::cellmappings::ComputeTemperatureMapping<LinkedCell, dim> computeTemperatureMapping(meanVelocity, _mdSolverInterface);
      91           0 :     cell.iterateConstCells(computeTemperatureMapping);
      92           0 :     double currentTemperature = computeTemperatureMapping.getTemperature();
      93             : 
      94           0 :     coupling::cellmappings::SetTemperatureMapping<LinkedCell, dim> setTemperatureMapping(currentTemperature, temperature, meanVelocity, _mdSolverInterface);
      95           0 :     cell.iterateCells(setTemperatureMapping);
      96           0 :   }
      97             : 
      98             : private:
      99             :   /** interface of the md solver*/
     100             :   coupling::interface::MDSolverInterface<LinkedCell, dim>* const _mdSolverInterface;
     101             : };
     102             : #endif // _MOLECULARDYNAMICS_COUPLING_KINETICENERGYCONTROLLER_H_