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             : 
       6             : template <class LinkedCell, unsigned int dim>
       7           0 : coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::TransferStrategy4NieCoupling(
       8             :     coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface, unsigned int numberMDSteps, double shiftTimestep,
       9             :     tarch::la::Vector<2 * dim, bool> massFluxBoundary)
      10           0 :     : coupling::transferstrategies::TransferStrategy<LinkedCell, dim>(mdSolverInterface), _massMapping(mdSolverInterface), _momentumMapping(mdSolverInterface),
      11           0 :       _oldSolution(NULL), _newSolution(NULL), _numberMDSteps(numberMDSteps), _shiftTimestep(shiftTimestep), _numberLocalCells(I02::linearNumberCellsInDomain),
      12           0 :       _excessMass(NULL), _massFluxBoundary(massFluxBoundary) {
      13             :   // allocate arrays for old and new solution (required for linear time
      14             :   // interpolation) and init. them with zero vectors
      15           0 :   _oldSolution = new tarch::la::Vector<dim, double>[_numberLocalCells];
      16           0 :   _newSolution = new tarch::la::Vector<dim, double>[_numberLocalCells];
      17           0 :   _excessMass = new double[_numberLocalCells];
      18           0 :   if ((_oldSolution == NULL) || (_newSolution == NULL) || (_excessMass == NULL)) {
      19           0 :     std::cout << "ERROR "
      20             :                  "coupling::transferstrategies::TransferStrategy4NieCoupling::"
      21             :                  "TransferStrategy4NieCoupling(...): _ptr==NULL!"
      22           0 :               << std::endl;
      23           0 :     exit(EXIT_FAILURE);
      24             :   }
      25           0 :   for (unsigned int i = 0; i < _numberLocalCells; i++) {
      26           0 :     _oldSolution[i] = tarch::la::Vector<dim, double>(0.0);
      27           0 :     _newSolution[i] = tarch::la::Vector<dim, double>(0.0);
      28           0 :     _excessMass[i] = 0.0;
      29             :   }
      30           0 : }
      31             : 
      32           0 : template <class LinkedCell, unsigned int dim> coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::~TransferStrategy4NieCoupling() {
      33             :   // delete arrays
      34           0 :   if (_oldSolution != NULL) {
      35           0 :     delete[] _oldSolution;
      36           0 :     _oldSolution = NULL;
      37             :   }
      38           0 :   if (_newSolution != NULL) {
      39           0 :     delete[] _newSolution;
      40           0 :     _newSolution = NULL;
      41             :   }
      42           0 :   if (_excessMass != NULL) {
      43           0 :     delete[] _excessMass;
      44           0 :     _excessMass = NULL;
      45             :   }
      46           0 : }
      47             : 
      48             : template <class LinkedCell, unsigned int dim>
      49           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::beginProcessInnerCouplingCellsBeforeReceivingMacroscopicSolverData() {
      50             :   // reset time counter (always counts from 0 to _numberMDSteps)
      51           0 :   _timestepCounter = 0;
      52             :   //_totalMass=0;
      53             :   //_cellCount=0;
      54             :   // copy new solution to old solution
      55           0 :   for (unsigned int i = 0; i < _numberLocalCells; i++) {
      56           0 :     _oldSolution[i] = _newSolution[i];
      57             :   }
      58           0 : }
      59             : 
      60             : template <class LinkedCell, unsigned int dim>
      61           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellBeforeReceivingMacroscopicSolverData(
      62             :     coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
      63             :   // backup old mass that could not be inserted into cells
      64           0 :   _excessMass[index.get()] = cell.getMicroscopicMass();
      65             : 
      66             :   // TODO: remove temporary debug output
      67             :   // if(_excessMass[index.get()] > 0) std::cout << "TransferStrategy4NieCoupling:
      68             :   // _excessMass[" << index << "]: " << _excessMass[index.get()] << std::endl;
      69             :   // TODO: make threshold customizable
      70             :   // const unsigned int threshold = 3;
      71             :   /*if(_excessMass[index.get()] > threshold) {
      72             :         std::cout << "TransferStrategy4NieCoupling: ERROR: Excess mass above
      73             :   threshold (= " << threshold << "): _excessMass[" << index << "] = " <<
      74             :   _excessMass[index.get()] << std::endl; exit(EXIT_FAILURE);
      75             :   }*/
      76             : 
      77             :   // reset cell-mass (would be overwritten anyway; but just to be on safe side)
      78           0 :   cell.setMicroscopicMass(0.0);
      79           0 : }
      80             : 
      81             : template <class LinkedCell, unsigned int dim>
      82           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellAfterReceivingMacroscopicSolverData(
      83             :     coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
      84             : 
      85             :   // convert momentum to velocity values for cont->MD transfer
      86           0 :   if (cell.getMicroscopicMass() == 0.0) {
      87           0 :     cell.setMicroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
      88             :   } else {
      89           0 :     cell.setMicroscopicMomentum((1.0 / cell.getMicroscopicMass()) * cell.getMicroscopicMomentum());
      90             :   }
      91             : 
      92             :   // write velocity value to new solution array and interpolate first solution
      93             :   // (= value of old solution)
      94           0 :   _newSolution[index.get()] = cell.getMicroscopicMomentum();
      95           0 :   cell.setMicroscopicMomentum((1.0 - _shiftTimestep) * _oldSolution[index.get()] + _shiftTimestep * _newSolution[index.get()]);
      96             : 
      97           0 :   const tarch::la::Vector<dim, double> velocity = (0.5 - _shiftTimestep) * _oldSolution[index.get()] + (_shiftTimestep + 0.5) * _newSolution[index.get()];
      98             : 
      99             :   // for mass transfer: compute mass flux if required and add excess mass
     100           0 :   const double massFlux = computeMassFlux(cell.getMicroscopicMass(), velocity, index) + _excessMass[index.get()];
     101           0 :   cell.setMicroscopicMass(massFlux);
     102             : 
     103             :   // reset macroscopic solver values (from averaging)
     104           0 :   cell.setMacroscopicMass(0.0);
     105           0 :   cell.setMacroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
     106           0 : }
     107             : 
     108             : // template<class LinkedCell,unsigned int dim>
     109             : // void
     110             : // coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell,dim>::
     111             : // beginProcessInnerCouplingCellsBeforeSendingMDSolverData(){
     112             : // }
     113             : 
     114             : template <class LinkedCell, unsigned int dim>
     115           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellBeforeSendingMDSolverData(
     116             :     coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
     117             : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
     118             :   std::cout << "processInnerCouplingCellBeforeSendingMDSolverData(): Data "
     119             :                "before averaging from cell "
     120             :             << index << "Mass: " << cell.getMacroscopicMass() << " , momentum: " << cell.getMacroscopicMomentum() << std::endl;
     121             : #endif
     122             : 
     123             :   // average sampled data from MD
     124           0 :   if (_timestepCounter != 0) {
     125           0 :     const double mass = cell.getMacroscopicMass() / ((double)_timestepCounter);
     126           0 :     const tarch::la::Vector<dim, double> momentum = (1.0 / ((double)_timestepCounter)) * cell.getMacroscopicMomentum();
     127           0 :     cell.setMacroscopicMass(mass);
     128           0 :     cell.setMacroscopicMomentum(momentum);
     129             :   } else {
     130           0 :     cell.setMacroscopicMass(0.0);
     131           0 :     cell.setMacroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
     132             :   }
     133             : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
     134             :   std::cout << "Time step counter=" << _timestepCounter << std::endl;
     135             :   std::cout << "processInnerCouplingCellBeforeSendingMDSolverData(): Data "
     136             :                "from cell "
     137             :             << index << "Mass: " << cell.getMacroscopicMass() << " , momentum: " << cell.getMacroscopicMomentum() << std::endl;
     138             : #endif
     139             : 
     140             :   // cell.iterateConstCells(_massMapping);
     141             :   // _totalMass += _massMapping.getMass();
     142             :   // _totalMass +=  cell.getMicroscopicMass();
     143           0 : }
     144             : 
     145             : // template<class LinkedCell,unsigned int dim>
     146             : // void
     147             : // coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell,dim>::
     148             : // endProcessInnerCouplingCellsBeforeSendingMDSolverData(){
     149             : // std::cout << "Avg massFluxBoundary MD density: " <<
     150             : // _totalMass/(_cellCount*15.625) << std::endl;
     151             : //}
     152             : 
     153             : template <class LinkedCell, unsigned int dim>
     154           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::beginProcessInnerCouplingCellsAfterMDTimestep() {
     155             :   // increment time step counter
     156           0 :   _timestepCounter++;
     157           0 : }
     158             : 
     159             : template <class LinkedCell, unsigned int dim>
     160           0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellAfterMDTimestep(
     161             :     coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
     162             :   // sampling: add up mass and momentum in macroscopic buffers
     163           0 :   cell.iterateConstCells(_massMapping);
     164           0 :   const double mass = _massMapping.getMass();
     165             :   //_totalMass += mass;
     166             :   //_cellCount++;
     167           0 :   cell.iterateConstCells(_momentumMapping);
     168           0 :   const tarch::la::Vector<dim, double> momentum(_momentumMapping.getMomentum());
     169           0 :   cell.addMacroscopicMass(mass);
     170           0 :   cell.addMacroscopicMomentum(momentum);
     171             : 
     172             :   // interpolation of velocity value from old and new solution
     173             :   // we shift the factor by _shiftTimestep time steps and inter- or extrapolate
     174             :   // velocities. A CFD solver will deliver values at t and t+dt and we want to
     175             :   // advance the simulation potentially differently
     176           0 :   const double factor = _shiftTimestep + ((double)_timestepCounter) / ((double)_numberMDSteps);
     177             : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
     178             :   if (factor < _shiftTimestep || factor > 1.0 + _shiftTimestep) {
     179             :     std::cout << "ERROR "
     180             :                  "coupling::transferstrategies::TransferStrategy4NieCoupling<"
     181             :                  "LinkedCell,dim>::processInnerCouplingCellAfterMDTimestep(."
     182             :                  ".): factor out of range!"
     183             :               << std::endl;
     184             :     exit(EXIT_FAILURE);
     185             :   }
     186             : #endif
     187           0 :   const tarch::la::Vector<dim, double> velocity = (1.0 - factor) * _oldSolution[index.get()] + factor * _newSolution[index.get()];
     188           0 :   cell.setMicroscopicMomentum(velocity);
     189           0 : }
     190             : 
     191             : template <class LinkedCell, unsigned int dim>
     192           0 : double coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::computeMassFlux(const double& mass,
     193             :                                                                                                     const tarch::la::Vector<dim, double>& velocity,
     194             :                                                                                                     const I01 index) {
     195             :   // meshsize and volume of cell
     196           0 :   const tarch::la::Vector<dim, double> meshsize(IDXS.getCouplingCellSize());
     197             :   double volume = 1.0;
     198           0 :   for (unsigned int d = 0; d < dim; d++) {
     199           0 :     volume = volume * meshsize[d];
     200             :   }
     201             : 
     202             :   // mass density of cell
     203           0 :   const double density = mass / volume;
     204             :   // coupling time interval
     205           0 :   const double couplingTimeInterval = coupling::transferstrategies::TransferStrategy<LinkedCell, dim>::_mdSolverInterface->getDt() * _numberMDSteps;
     206             : 
     207           0 :   double massFlux = 0.0;
     208             : 
     209             :   // compute normal vector for this cell, checking whether there is mass flux or
     210             :   // not
     211           0 :   tarch::la::Vector<dim, int> n(0);
     212           0 :   for (unsigned int d = 0; d < dim; d++) {
     213             :     // compute normal vector (pointing inwards MD domain)
     214           0 :     if ((index.get()[d] == 1) && _massFluxBoundary[2 * d]) {
     215           0 :       n[d] = 1; /*_totalMass+=mass;_cellCount++;*/
     216           0 :     } else if ((index.get()[d] == (int)I09::numberCellsInDomain[d]) && _massFluxBoundary[2 * d + 1]) {
     217           0 :       n[d] = -1; /*_totalMass+=mass;_cellCount++;*/
     218             :     }
     219             :   }
     220             : 
     221             :   // compute mass flux for each boundary
     222           0 :   for (unsigned int d = 0; d < dim; d++) {
     223             :     double surface = 1.0;
     224           0 :     for (unsigned int e = 0; e < d; e++) {
     225           0 :       surface = surface * meshsize[e];
     226             :     }
     227           0 :     for (unsigned int e = d + 1; e < dim; e++) {
     228           0 :       surface = surface * meshsize[e];
     229             :     }
     230             : 
     231           0 :     massFlux += n[d] * surface * density * velocity[d] * couplingTimeInterval;
     232             :   }
     233           0 :   return massFlux;
     234             : }