externalForceMrtDynamics.hh

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* Orestis Malaspinas contributed this code. */

#ifndef EXTERNAL_FORCE_MRT_DYNAMICS_HH
#define EXTERNAL_FORCE_MRT_DYNAMICS_HH

#include "externalForceMrtDynamics.h"
#include "latticeBoltzmann/mrtTemplates.h"
#include "latticeBoltzmann/dynamicsTemplates.h"
#include "latticeBoltzmann/momentTemplates.h"
#include "core/latticeStatistics.h"
#include <algorithm>
#include <limits>

namespace plb {
    
/* *************** Class ShanChenExternalForceMRTdynamics *********************************************** */

template<typename T, template<typename U> class Descriptor>
int GuoExternalForceMRTdynamics<T,Descriptor>::id =
    meta::registerGeneralDynamics<T,Descriptor,GuoExternalForceMRTdynamics<T,Descriptor> >("MRT_ExternalForce_Guo");

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>::GuoExternalForceMRTdynamics(plint externalParam_)
    : ExternalForceDynamics<T,Descriptor>(mrtParam<T,Descriptor>().get(externalParam_).getOmega()),
      param(0),
      externalParam(externalParam_)
{ } 

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>::GuoExternalForceMRTdynamics(MRTparam<T,Descriptor>* param_)
    : ExternalForceDynamics<T,Descriptor>(param_->getOmega()),
      param(param_),
      externalParam(-1)
{ }

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>::GuoExternalForceMRTdynamics(HierarchicUnserializer& unserializer) 
    : ExternalForceDynamics<T,Descriptor>(T())
{
    unserialize(unserializer);
}

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>::GuoExternalForceMRTdynamics(GuoExternalForceMRTdynamics<T,Descriptor> const& rhs)
    : ExternalForceDynamics<T,Descriptor>(rhs),
      param(rhs.param ? rhs.param->clone() : 0),
      externalParam(rhs.externalParam)
{ }

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>&
    GuoExternalForceMRTdynamics<T,Descriptor>::operator=(GuoExternalForceMRTdynamics<T,Descriptor> const& rhs)
{
    GuoExternalForceMRTdynamics<T,Descriptor>(rhs).swap(*this);
    return *this;
}

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>::~GuoExternalForceMRTdynamics()
{
    delete param;
}

template<typename T, template<typename U> class Descriptor>
void GuoExternalForceMRTdynamics<T,Descriptor>::swap(GuoExternalForceMRTdynamics<T,Descriptor>& rhs)
{
    std::swap(param, rhs.param);
    std::swap(externalParam, rhs.externalParam);
    T tmpOmega = this->getOmega();
    this->setOmega(rhs.getOmega());
    rhs.setOmega(tmpOmega);
}

template<typename T, template<typename U> class Descriptor>
GuoExternalForceMRTdynamics<T,Descriptor>* GuoExternalForceMRTdynamics<T,Descriptor>::clone() const {
    return new GuoExternalForceMRTdynamics<T,Descriptor>(*this);
}
 
template<typename T, template<typename U> class Descriptor>
int GuoExternalForceMRTdynamics<T,Descriptor>::getId() const {
    return id;
}

template<typename T, template<typename U> class Descriptor>
void GuoExternalForceMRTdynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
{
    int useExternalParamFlag = param ? 0:1;
    serializer.addValue(useExternalParamFlag);
    if (param) {
        param->serialize(serializer);
    }
    else {
        serializer.addValue(externalParam);
    }
    ExternalForceDynamics<T,Descriptor>::serialize(serializer);
}

template<typename T, template<typename U> class Descriptor>
void GuoExternalForceMRTdynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
{
    int useExternalParamFlag = unserializer.readValue<int>();
    if (useExternalParamFlag) {
        externalParam = unserializer.readValue<plint>();
        param = 0;
    }
    else {
        delete param;
        param = new MRTparam<T,Descriptor>(unserializer);
        externalParam = -1;
    }
    ExternalForceDynamics<T,Descriptor>::unserialize(unserializer);
}

template<typename T, template<typename U> class Descriptor>
void GuoExternalForceMRTdynamics<T,Descriptor>::collide (
        Cell<T,Descriptor>& cell, BlockStatistics& statistics )
{
    typedef mrtTemplates<T,Descriptor> mrtTemp;
    
    T rhoBar = momentTemplates<T,Descriptor>::get_rhoBar(cell);
    MRTparam<T,Descriptor>* parameter = param ? param : &(mrtParam<T,Descriptor>().get(externalParam));
    
    Array<T,Descriptor<T>::d> u, j;
    this->computeVelocity(cell, u);
    T rho = Descriptor<T>::fullRho(rhoBar);
    for (plint iD = 0; iD < Descriptor<T>::d; ++iD) {
        j[iD] = rho * u[iD];
    }
    
    T jSqr = mrtTemp::mrtCollisionWithForce(cell, rhoBar, j, parameter->getInvM(), (T)1);
    
    if (cell.takesStatistics()) {
        T invRho = Descriptor<T>::invRho(rhoBar);
        gatherStatistics(statistics, rhoBar, jSqr*invRho*invRho);
//         gatherStatistics(statistics, rhoBar, uSqr);
    }
}

template<typename T, template<typename U> class Descriptor>
T GuoExternalForceMRTdynamics<T,Descriptor>::computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                                T jSqr, T thetaBar) const
{
    T invRho = Descriptor<T>::invRho(rhoBar);
    return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
}

template<typename T, template<typename U> class Descriptor>
MRTparam<T,Descriptor> const& GuoExternalForceMRTdynamics<T,Descriptor>::getMrtParameter() const {
    MRTparam<T,Descriptor>* parameter = param ? param : &(mrtParam<T,Descriptor>().get(externalParam));
    return *parameter;
}

}
#endif