mrtDynamics.h

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 MRT_DYNAMICS_H
#define MRT_DYNAMICS_H

#include "core/globalDefs.h"
#include "basicDynamics/isoThermalDynamics.h"

namespace plb {

template<typename T, template<typename U> class Descriptor>
class MRTparam {
public:
    typedef T MatT[Descriptor<T>::q][Descriptor<T>::q];
public:
    MRTparam();
    virtual ~MRTparam() { }
    MRTparam(T omega_);
    MRTparam(T omega_, T lambda_);
    MRTparam(std::vector<T> s_);
    MRTparam(HierarchicUnserializer& unserializer);
    void serialize(HierarchicSerializer& serializer) const;
    void unserialize(HierarchicUnserializer& unserializer);
    MRTparam<T,Descriptor>* clone() const;
    virtual void setParameter(plint whichParameter, T value);
    virtual T getParameter(plint whichParameter) const;
    virtual T getOmega() const;
    virtual void setOmega(T omega_);
    T getLambda() const;
    void setLambda(T lambda_);
    T getOmegaQ() const;
    void setOmegaQ(T q_);
    T getOmegaEpsilon() const;
    void setOmegaEpsilon(T epsilon_);
    std::vector<T> const& getS() const;
    MatT& getInvM();
private:
    void iniRelaxationVector();
    void precompute_invM_S();
private:
    T invM_S[Descriptor<T>::q][Descriptor<T>::q]; // relaxation time matrix.
    std::vector<T> s;
};

template<typename T, template<typename U> class Descriptor>
class MRTparamList {
public:
    void set(plint id, MRTparam<T,Descriptor> const& param);
    MRTparam<T,Descriptor>& get(plint id);
    MRTparam<T,Descriptor> const& get(plint id) const;
private:
    std::vector<MRTparam<T,Descriptor> > parameters;
};

template<typename T, template<typename U> class Descriptor>
MRTparamList<T,Descriptor>& mrtParam();

template<typename T, template<typename U> class Descriptor>
class MRTdynamics : public IsoThermalBulkDynamics<T,Descriptor> {
public:
/* *************** Construction / Destruction ************************ */
    MRTdynamics(plint externalParam_);
    MRTdynamics(MRTparam<T,Descriptor>* param_);
    MRTdynamics(HierarchicUnserializer& unserializer);
    MRTdynamics(MRTdynamics<T,Descriptor> const& rhs);
    MRTdynamics<T,Descriptor>& operator=(MRTdynamics<T,Descriptor> const& rhs);
    ~MRTdynamics<T,Descriptor>();
    void swap(MRTdynamics<T,Descriptor>& rhs);
    
    virtual MRTdynamics<T,Descriptor>* clone() const;

    virtual int getId() const;

    virtual void serialize(HierarchicSerializer& serializer) const;
    virtual void unserialize(HierarchicUnserializer& unserializer);

/* *************** Collision and Equilibrium ************************* */

    virtual void collide(Cell<T,Descriptor>& cell,
                         BlockStatistics& statistics_);

    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const;
                                 
    MRTparam<T,Descriptor> const& getMrtParameter() const;
private:
    MRTparam<T,Descriptor>* param;
    plint externalParam;
    static int id;
};

template<typename T, template<typename U> class Descriptor>
class ExternalVelocityMRTdynamics : public IsoThermalBulkDynamics<T,Descriptor> {
public:
/* *************** Construction / Destruction ************************ */
    ExternalVelocityMRTdynamics(plint externalParam_);
    ExternalVelocityMRTdynamics(MRTparam<T,Descriptor>* param_);
    ExternalVelocityMRTdynamics(HierarchicUnserializer& unserializer);
    ExternalVelocityMRTdynamics(ExternalVelocityMRTdynamics<T,Descriptor> const& rhs);
    ExternalVelocityMRTdynamics<T,Descriptor>& operator=(ExternalVelocityMRTdynamics<T,Descriptor> const& rhs);
    ~ExternalVelocityMRTdynamics<T,Descriptor>();
    void swap(ExternalVelocityMRTdynamics<T,Descriptor>& rhs);
    
    virtual ExternalVelocityMRTdynamics<T,Descriptor>* clone() const;

    virtual int getId() const;

    virtual void serialize(HierarchicSerializer& serializer) const;
    virtual void unserialize(HierarchicUnserializer& unserializer);

/* *************** Collision and Equilibrium ************************* */

    virtual void collide(Cell<T,Descriptor>& cell,
                         BlockStatistics& statistics_);

    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const;
                                 
/* *************** Moments ******************************************* */

    // The function computeDensity() is not overridden, and the default
    //   implementation is kept, for two reasons. First, it is equivalent to access
    //   rho from the external scalar or to recompute it from the bulk (this
    //   is not the case for the velocity, in a Shan/Chen multicomponent model).
    //   Second, the Shan/Chen data-processor needs computeDensity() to be
    //   default implemented, because it uses this function to treat walls with
    //   a virtual-density mechanism.

    virtual void computeVelocity( Cell<T,Descriptor> const& cell,
                                  Array<T,Descriptor<T>::d>& u ) const;

    virtual T computeRhoBar(Cell<T,Descriptor> const& cell) const;


    virtual void computeRhoBarJ(Cell<T,Descriptor> const& cell,
                                T& rhoBar, Array<T,Descriptor<T>::d>& j) const;
                                 
public:
    MRTparam<T,Descriptor> const &getMrtParameter() const {
        return *param;
    }
private:
    MRTparam<T,Descriptor>* param;
    plint externalParam;
    static int id;
};


template<typename T, template<typename U> class Descriptor>
class IncMRTdynamics : public IsoThermalBulkDynamics<T,Descriptor> {
public:
/* *************** Construction / Destruction ************************ */
    IncMRTdynamics(plint externalParam_);
    IncMRTdynamics(MRTparam<T,Descriptor>* param_);
    IncMRTdynamics(HierarchicUnserializer& unserializer);
    IncMRTdynamics(IncMRTdynamics<T,Descriptor> const& rhs);
    IncMRTdynamics<T,Descriptor>& operator=(IncMRTdynamics<T,Descriptor> const& rhs);
    ~IncMRTdynamics<T,Descriptor>();
    void swap(IncMRTdynamics<T,Descriptor>& rhs);
    
    virtual IncMRTdynamics<T,Descriptor>* clone() const;

    virtual int getId() const;

    virtual void serialize(HierarchicSerializer& serializer) const;
    virtual void unserialize(HierarchicUnserializer& unserializer);

/* *************** Collision and Equilibrium ************************* */

    virtual void collide(Cell<T,Descriptor>& cell,
                         BlockStatistics& statistics_);

    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const;
                                 
    virtual bool velIsJ() const;

/* *************** Macroscopic variables ***************************** */
    
    virtual void computeVelocity( Cell<T,Descriptor> const& cell,
                                  Array<T,Descriptor<T>::d>& u ) const;
public:
    MRTparam<T,Descriptor> const &getMrtParameter() const {
        return *param;
    }
private:
    MRTparam<T,Descriptor>* param;
    plint externalParam;
    static int id;
};



template<typename T, template<typename U> class Descriptor>
class ExternalVelocityIncMRTdynamics : public IsoThermalBulkDynamics<T,Descriptor> {
public:
/* *************** Construction / Destruction ************************ */
    ExternalVelocityIncMRTdynamics(plint externalParam_);
    ExternalVelocityIncMRTdynamics(MRTparam<T,Descriptor>* param_);
    ExternalVelocityIncMRTdynamics(HierarchicUnserializer& unserializer);
    ExternalVelocityIncMRTdynamics(ExternalVelocityIncMRTdynamics<T,Descriptor> const& rhs);
    ExternalVelocityIncMRTdynamics<T,Descriptor>& operator=(ExternalVelocityIncMRTdynamics<T,Descriptor> const& rhs);
    ~ExternalVelocityIncMRTdynamics<T,Descriptor>();
    void swap(ExternalVelocityIncMRTdynamics<T,Descriptor>& rhs);
    
    virtual ExternalVelocityIncMRTdynamics<T,Descriptor>* clone() const;

    virtual int getId() const;

    virtual void serialize(HierarchicSerializer& serializer) const;
    virtual void unserialize(HierarchicUnserializer& unserializer);

/* *************** Collision and Equilibrium ************************* */

    virtual void collide(Cell<T,Descriptor>& cell,
                         BlockStatistics& statistics_);

    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const;
                                 
    virtual bool velIsJ() const;
    
/* *************** Moments ******************************************* */

    // The function computeDensity() is not overridden, and the default
    //   implementation is kept, for two reasons. First, it is equivalent to access
    //   rho from the external scalar or to recompute it from the bulk (this
    //   is not the case for the velocity, in a Shan/Chen multicomponent model).
    //   Second, the Shan/Chen data-processor needs computeDensity() to be
    //   default implemented, because it uses this function to treat walls with
    //   a virtual-density mechanism.

    virtual void computeVelocity( Cell<T,Descriptor> const& cell,
                                  Array<T,Descriptor<T>::d>& u ) const;

    virtual T computeRhoBar(Cell<T,Descriptor> const& cell) const;


    virtual void computeRhoBarJ(Cell<T,Descriptor> const& cell,
                                T& rhoBar, Array<T,Descriptor<T>::d>& j) const;
public:
    MRTparam<T,Descriptor> const &getMrtParameter() const {
        return *param;
    }
private:
    MRTparam<T,Descriptor>* param;
    plint externalParam;
    static int id;
};

}  // namespace plb

#endif  // MRT_DYNAMICS_H