calculation of momentum flux PI.

Palabos General
#1

Is the momentun flux PI, which is second-order moments of the non-equilibrium distribution function:

PI_ij = Sigma_a e_ai e_aj *[f_a - f_a^{eq}].

From momentTemplates2D.h, there are a few functions to compute the PiNeq


static void compute_PiNeq(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,2> const& j, Array<T,3>& PiNeq, T invRho)
{
    typedef SymmetricTensorImpl<T,2> S;

    T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
    partial_rho(f, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);

    PiNeq[S::xx] = lineX_P1+lineX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
    PiNeq[S::yy] = lineY_P1+lineY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
    PiNeq[S::xy] = -f[1] + f[3] - f[5] + f[7] - invRho*j[0]*j[1];
}



static void compute_rhoBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,2>& j, Array<T,3>& PiNeq)
{
    typedef SymmetricTensorImpl<T,2> S;

    T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
    partial_rho(f, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);

    rhoBar = lineX_P1 + lineX_0 + lineX_M1;
    j[0]   = lineX_P1 - lineX_M1;
    j[1]   = lineY_P1 - lineY_M1;
    T invRho = Descriptor::invRho(rhoBar);
    PiNeq[S::xx] = lineX_P1+lineX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
    PiNeq[S::yy] = lineY_P1+lineY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
    PiNeq[S::xy] = -f[1] + f[3] - f[5] + f[7] - invRho*j[0]*j[1];
}.

It seems that the PiNeq does not match the formula of PI_ij = Sigma_a e_ai e_aj *[f_a - f_a^{eq}]. Could you please provide a reference for the calculation of PiNeq in compute_rhoBar_j_PiNeq.

Plus, in function compute_rhoBar_j_PiNeq, the j seems not to represent the /
/ Momentum is equal to density times velocity (component-wise)


j[iD] = rho * u[iD];

However, in file smagorinskyDynamics.hh,


template<typename T, template<typename U> class Descriptor>
void SmagorinskyBGKdynamics<T,Descriptor>::collide (
        Cell<T,Descriptor>& cell,
        BlockStatistics& statistics )
{
    T rhoBar;
    Array<T,Descriptor<T>::d> j;
    Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
    momentTemplates<T,Descriptor>::compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq);
    T omega = SmagoOperations<T,Descriptor>::computeOmega (
            omega0, preFactor, rhoBar, PiNeq );
    T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, omega);
    if (cell.takesStatistics()) {
        gatherStatistics(statistics, rhoBar, uSqr);
    }
}

the j in bgk_ma2_collision(cell, rhoBar, j, omega) should be the momentum
// Momentum is equal to density times velocity (component-wise)


j[iD] = rho * u[iD];

but its value is based on the j cacluated by compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq), which is not the momentum.

Hopefully, Dr. Latt or Dr. Malaspinas could provide me some help on this or reference.

Cheers.

#2

Dear Palabos Developers,

It seems that my class template inherited from IsoThermalBulkDynamics does not work properly as I checked the three files:

[ul]
[1] momentTemplates.h
[2] momentTemplates2D.h
[3] momentTemplates3D.h
[/ul]

And found out that all the helper functions have been overloaded. The implemetions in momentTemplates.h seems to be straightforward, however, the implementations in momentTemplates2D.h and momentTemplates3D.h appears to be a little bit of confusing. Could you please provide a reference for that?

Now my questions is that why the functions in my class template calling these helper functions do not work properly and they seem to be calling the implemetations 2D or 3D specialization. I am a little bit rusty on C++.

Plus, according the limited documetion on Extending without modifying, I mimic the implemention of BGKdynamics and create two new files myNewModel.h and myNewModel.hh that were put in the same directory as the project code, containing a class MyNewDynamics which inherits from IsoThermalBulkDynamics and defines the same three methods as BGKdynamics, adapted to my case.

#3

Hi,

The mechanism involved in momentTemplates.h and momentTemplatesXD.h is not function overloading, but template specialization.

In momentTemplates.h you find generic code, that works for any LB lattice (D2Q9, D3Q19, D3Q27, etc.). In momentTemplates3D.h for example, the same code is written again for a specific 3D lattice, like D3Q19, and is optimized for this case. This is why it looks more confusing to you, because the various optimizations (loop unrolling, aggregation of common arithmetic terms, etc.) leave the code in a more complex shape.

The effect of the code is the same in momentTemplates.h and momentTemplates3D.h, though. That’s why I emphasize that this has nothing to do with function overloading. It’s just an efficient specialization of the code for specific lattices. If you want to read the code as a matter of documentation, it is always better to read the generic one, in momentTemplates.h.

To write your own model, you did the right thing, inheriting from IsoThermalBulkDynamics. As you implement the collision step, you can take the same path as the Palabos developers, and work with templates that are specialized (and optimized) in specific cases. To start with, I suggest however that you skip this technically advanced step, and insert the generic version of your collision model right there, without referring to external templates.

Cheers,
Jonas

#4

Dear Dr. Latt,

Thank you very much for your reply. I will refresh my knowledge on template specilization.

Could you please give a little bit more explanation about the difference in Palabos?

Now that you metioned that “The effect of the code is the same in momentTemplates.h and momentTemplates3D.h, though.”, but in my class template when I call


get_rhoBar_j();

the j seems not to be the supposed momentum. I therefore reckon the called get_rhoBar_j() does not come form momentumTemplates.h.

Plus, what is the generic version of collision? In my class template, I mimic the GuoExternalForceBGKdynamics to add another force with another my own ExternalForceDynamics which inherits from IsoThermalBulkDynamics.

I noted that the collision method in GuoExternalForceBGKdynamics


    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 uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

, insteads of


    Array<T,Descriptor<T>::d> j;
    momentTemplates<T,Descriptor>::get_rhoBar_j(cell, rhoBar, j);
    T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

used in BGKdynamics. I thought I could use


    momentTemplates<T,Descriptor>::get_rhoBar_j(cell, rhoBar, j);

to get the momentun j, but it seems not be to correct.

Dear Dr. Latt, if you could explain this to me, properly I will understand a little bit more about the art of Palabos.

Thank you very much again.

jonas@flowkit Wrote:

As you

implement the collision step, you can take the
same path as the Palabos developers, and work with
templates that are specialized (and optimized) in
specific cases. To start with, I suggest however
that you skip this technically advanced step, and
insert the generic version of your collision model
right there, without referring to external
templates.

Cheers,
Jonas

#5

Hi,

You are confusing the words “momentum” (which is velocity times density) and “moment” (the macroscopic variables are moments of the particle distributions).

Cheers,
Jonas

#6

Hi,

From GuoExternalBGKdynamics, I thought j is momentum as appeared in the code

 
  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 uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

, however, in BGKdynamics, you mean the j is moments rather than momentum?
I thought


    T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

applies to collision that appears in different dynamics with j.

jonas@flowkit Wrote:

Hi,

You are confusing the words “momentum” (which is
velocity times density) and “moment” (the
macroscopic variables are moments of the particle
distributions).

Cheers,
Jonas

#7

Dear Dr. Latt,

in bgk_ma2_collision,


static T bgk_ma2_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T omega) {
    T invRho = Descriptor::invRho(rhoBar);
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        f[iPop] *= (T)1-omega;
        f[iPop] += omega * dynamicsTemplatesImpl<T,Descriptor>::bgk_ma2_equilibrium (
                                iPop, rhoBar, invRho, j, jSqr );
    }
    return jSqr*invRho*invRho;
}

,
my understanding is that j is momentum rather than moments of particle distribution function. I hope made a sense about the j.

Cheers.

#8

Dear Dr. Latt,

Yes, you are right.

from


static void get_j(Array<T,Descriptor::q> const& f, Array<T,Descriptor::d>& j ) {
    for (int iD=0; iD < Descriptor::d; ++iD) {
        j[iD] = f[0]*Descriptor::c[0][iD];
    }
    for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
        for (int iD=0; iD < Descriptor::d; ++iD) {
            j[iD] += f[iPop]*Descriptor::c[iPop][iD];
        }
    }
}

,
j is the moments from particle distribution fucntion f[iPop] for the macrovariable velocity.

But rho = 1. So j equals to the momentum of density rho times velocity.

I am still confusing about the difference within the collision,


  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 uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

from GuoExternalForceBGKdynamics and


    momentTemplates<T,Descriptor>::get_rhoBar_j(cell, rhoBar, j);
    T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());

.

What’s the difference betweent the two j?

Thank you very much for your time.

jonas@flowkit Wrote:

Hi,

You are confusing the words “momentum” (which is
velocity times density) and “moment” (the
macroscopic variables are moments of the particle
distributions).

Cheers,
Jonas

#9

Dear Dr. Latt,

Properly you feel a little bit disappointed about my knowledge on LBM.

From the limited documentation of Palabos, I found out the following part:

Velocity and Density
Velocity and mass density are the two most recurrent macroscopic variables in lattice Boltzmann. It is therefore
surprising that variables with names like rho or u rarely occur inside the Palabos code (nothing prevents you from
perusing them in end-user code, though). Instead, internal Palabos structures most often work with the variables
rhoBar and j. The variable j is nothing else than the first-order velocity moment which, for most models, is
identical with the fluid momentum: j = rho * u.

So, I mistook the momentum and moments and regarded them as two equal quantities. When I ignore the step (I thought rho = 1. )




    for (plint iD = 0; iD < Descriptor<T>::d; ++iD)
    {
        j[iD] = rho * u[iD];
    }

my results seems to be very weird.

However, when I add the above step, everything seems to work. What caused this problem (round-off error)?

Could I get back the original question? How about the momentum flux?

Thanks.

jonas@flowkit Wrote:

Hi,

You are confusing the words “momentum” (which is
velocity times density) and “moment” (the
macroscopic variables are moments of the particle
distributions).

Cheers,
Jonas

#10

Dear Dr. Latt,

from Palabos online documentation about velocity and density,

Velocity and Density

Velocity and mass density are the two most recurrent macroscopic variables in lattice Boltzmann. It is therefore surprising that variables with names like rho or u rarely occur inside the Palabos code (nothing prevents you from perusing them in end-user code, though). Instead, internal Palabos structures most often work with the variables rhoBar and j. The variable j is nothing else than the first-order velocity moment which, for most models, is identical with the fluid momentum: j = rho * u. The definition of the variable rhoBar on the other hand can be customized. By default, it is defined as rhoBar=rho-1 and is used to improve the numerical accuracy of the method. Indeed, as the density is often close to 1, you gain significant digits in the representation of floating point variables by representing a quantity fluctuating around 0 rather than 1.

Without getting into details, here are a few rules for switching between the rho / u and the rhoBar / j representation:


// Use custom definitions in lattice descriptor to compute rho from rhoBar

rho = Descriptor<T>::fullRho(rhoBar);

// Use custom definitions in lattice descriptor to compute rhoBar from rho
rhoBar = Descriptor<T>::rhoBar(rho);

// Momentum is equal to density times velocity (component-wise)
j[iD] = rho * u[iD];

// Velocity is equal to inverse-density times momentum (component-wise)
u[iD] = 1./rho * j[iD];

// Inverse-density can be computed right away from rhoBar. Depending on
//   the custom definitions in the descriptor, this is substantially more
//   efficient than computing 1./rho, because the Taylor expansion of the
//   inverse-function is truncated at an appropriate position.
u[iD] = Descriptor<T>::invRho(rhoBar) * j[iD];

,

It was mentioned that the variable j is nothing else than the first-order velocity moment which, for most models, is identical with the fluid momentum: j = rho * u…

So I was confused about j.

Could you go a little more detail about “Without getting into details, here are a few rules for switching between the rho / u and the rhoBar / j representation:”?

Cheers.

palabosAlpha.

jonas@flowkit Wrote:

Hi,

You are confusing the words “momentum” (which is
velocity times density) and “moment” (the
macroscopic variables are moments of the particle
distributions).

Cheers,
Jonas