How to set time-dependent omega

Palabos General
#1

In the expamle code rayleighTaylor2D.cpp, the omega is a constant value, which is only component-denpendent. I want to set the omega as the time-denpendent values. I find the constructor as below,.


/// With this constructor, space- and time-dependent values of the
///   relaxation parameters omega are accounted for.
    ShanChenMultiComponentProcessor2D(T G_);

Also I find the command


virtual void setOmega(T omega_)=0

but I am not sure how to use this. Is there anyone know how to do that? Thank you. Another question is that I want to output the distribution fuction on each time step. wheter am I right?



pcout << "density fluid one = "
                  << *computeDensity(heavyFluid);
#2

Hello,

The simplest way of varying the relaxation parameter omega is to tell the Shan/Chen data processor that omega is not constant; the data processor will then read the value of omega from the dynamics object, on every cell, when it is needed. To do this, you need (as you pointed out) to construct the Shan/Chen data processor through the standard constructor:


new ShanChenMultiComponentProcessor2D(G);

After this, you can change omega as often as needed by accessing the cells of the lattice directly. A data processor which does this for you will be provided in the next Palabos release, but is not implemented in the current release. In the mean time, you can copy-paste the following code into your program:


namespace plb {

template<typename T, template<typename U> class Descriptor>
class AssignOmegaFunctional2D : public BoxProcessingFunctional2D_L<T,Descriptor>
{
public:
    AssignOmegaFunctional2D(T omega);
    virtual void process(Box2D domain, BlockLattice2D<T,Descriptor>& lattice);
    virtual AssignOmegaFunctional2D<T,Descriptor>* clone() const;
    virtual BlockDomain::DomainT appliesTo() const;
    virtual void getModificationPattern(std::vector<bool>& isWritten) const;
private:
    T omega;
};

/* ************* Class AssignOmegaFunctional2D ******************* */

template<typename T, template<typename U> class Descriptor>
AssignOmegaFunctional2D<T,Descriptor>::AssignOmegaFunctional2D(T omega_)
    : omega(omega_)
{ }

template<typename T, template<typename U> class Descriptor>
void AssignOmegaFunctional2D<T,Descriptor>::process (
        Box2D domain, BlockLattice2D<T,Descriptor>& lattice )
{
    for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
        for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
            lattice.get(iX,iY).getDynamics().setOmega(omega);
        }
    }
}

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

template<typename T, template<typename U> class Descriptor>
BlockDomain::DomainT AssignOmegaFunctional2D<T,Descriptor>::appliesTo() const
{
    // Omega needs to be set on envelope nodes as well, because the dynamics object
    //   is being modified.
    return BlockDomain::bulkAndEnvelope;
}

template<typename T, template<typename U> class Descriptor>
void AssignOmegaFunctional2D<T,Descriptor>::getModificationPattern (
        std::vector<bool>& isWritten) const
{
    isWritten[0] = true;
}

template<typename T, template<class U> class Descriptor>
void setOmega(MultiBlockLattice2D<T,Descriptor>& lattice, Box2D domain, T omega) {
    applyProcessingFunctional(new AssignOmegaFunctional2D<T,Descriptor>(omega), domain, lattice);
}

}

Then, you can modify the relaxation parameter at any time by invoking something like


setOmega(lattice, lattice.getBoundingBox(), omega);

Hope this works. If not, just let us know.

#3

Hi jlatt,

Thank you for your reply. I have tried it as you said, but it does not work. perhaps I make a mistake. Could you please help me to check?
Firstly, I copy-paste the code you listed in last , and then to construct the Shan/Chen data processor through the standard constructor, in the main I add some sentense as below


int main (int argc, char *argv[ ]) {
         ...
         ...
 integrateProcessingFunctional (
            new ShanChenMultiComponentProcessor2D<T,DESCRIPTOR>(G),
            Box2D(0,nx-1,0,ny-1),
            blockLattices,
            processorLevel );

   rayleighTaylorSetup(FluidOne, FluidTwo, force);  
   
  FluidOne.collideAndStream();
  FluidTwo.collideAndStream();

        if (iT%1==0) {
            for (int iY=0; iY<ny; ++iY) {
                for (int iX=0; iX<nx; ++iX) {
            if (FluidOne.get(iX,iY).computeDensity()>20 && outsideFluid.get(iX,iY).computeDensity()>20){
                         T omega1 = 0.5;
                         T omega2 = 1.2;
             }

             insideFluid.get(iX,iY).getDynamics().setOmega(omega1);
             outsideFluid.get(iX,iY).getDynamics().setOmega(omega2);
		
            }
                }
        }

Whether am I right? I appreciate your help. Thank you.

#4

Hi @jlatt,

Is there a way of setting two blobs of fluid1 in a container full of fluid2, where each of these blobs of fluid1 has a different viscosity (independent of their spatial location) ?

Thanks a lot!

#5

Do you have any hints @orestis ?

#6

Hi jlatt,

I am implementing a similar version to this problem. Though I am not implementing a multi-phase flow. I would just like to change the local omega of a standard implementation of a D3Q19 LBM fluid. I have written a data processor following your example provided above. I am making sure that omega is being set on the envelop along with the bulk i.e. BlockDomain::bulkAndEnvelope, and making sure that the getTypeofModification … { modified[0]= dynamicVariables} along with the getModificationPattern {isWritten[0] = True} …

The behavior I am observing is that omega is staying constant at the equilibrium value (say 0.54) only on the envelopes while everywhere else in the bulk it is getting set to the value I want (say 0.6). Is there something that I am missing when writing my data processor? I have attached my code below.

template <typename T, template <typename U> class Descriptor>
class MacroGetViscSetOmega : public plb::BoxProcessingFunctional3D_L<T, Descriptor>
{
public:
    plint level;
    MacroGetViscSetOmega(plint level) {
        this->level = level;
    }
    virtual void process(plb::Box3D domain, plb::BlockLattice3D<T, Descriptor>& lattice) {
    for (plint iX=domain.x0; iX<=domain.x1; ++iX){
        for (plint iY=domain.y0; iY<=domain.y1; ++iY){
            for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
            if(lattice.get(iX,iY,iZ).getDynamics().isBoundary())
              {continue;}
            else{
            lattice.get(iX,iY,iZ).getDynamics().setOmega(0.6);
            }
          }
          }
        }
    }
    virtual void getTypeOfModification(std::vector<plb::modif::ModifT> &modified) const
    {
        modified[0] = modif::dynamicVariables;
    }
    virtual plb::BlockDomain::DomainT appliesTo() const
    {
        return BlockDomain::bulkAndEnvelope;
    }
    virtual MacroGetViscSetOmega<T, Descriptor> *clone() const
    {
        return new MacroGetViscSetOmega<T, Descriptor>(*this);
    }
    virtual void getModificationPattern(std::vector<bool>& isWritten) const
    {
    isWritten[0] = true;
    }
};