Re: access collision and relaxation time

[Michal Januszewski]

Hi Luca,

The code for the collision step is in
sailfish/templates/relaxation.mako:BGK_relaxate(), and the relaxation
time is calculated by LBMSim.get_tau() and passed to the templates as
the 'tau' symbol.

If you want to just play around with some changes, it might be easier
to first generate standard code for the geometry you want to run with
(--save_src=code.cu), then modify code.cu with your changes, and
finally run a simulation with the modified code (--use_src=code.cu).
This way you will not have to deal with the code of the multitude of
LB models that Sailfish supports and that is sometimes mixed into the
standard BGK path (you will see that in BGK_relaxate).

Hope this helps and please let us know if run into any problems,
Michal

On Wed, Nov 23, 2011 at 11:47, Luca <luca.be...@gmail.com> wrote:
> Hi all,
>
> I would like to access the collision step of BGK model as well as the
> relaxation time, but I am not locate them (say in a .mako template?).
> What I want to do is to modify relaxation time and the expression of
> the velocities for a 2D case (ux, uy). Consider that I am a beginner
> either in LBM and  Sailfish, therefore any help would be extremely
> appreciated...
>
> Thanks!
>

--
Michal Januszewski
http://people.gentoo.org/spock

[Luca]

Hi Michal,

thanks for the answer, I will try to go through it right now! While
looking for a way to handle this I also found that via
"sym.equilibrium_expr(eq, eq_vars)" it is possible to directly modify
the equilibrium. Do you think this can be also a solution? I just
would like to set a constant velocity profile for ux and uy..

Thanks a lot!
Luca

On 25 nov, 01:49, Michal Januszewski <mich...@gmail.com> wrote:
> Hi Luca,
>
> The code for the collision step is in
> sailfish/templates/relaxation.mako:BGK_relaxate(), and the relaxation
> time is calculated by LBMSim.get_tau() and passed to the templates as
> the 'tau' symbol.
>
> If you want to just play around with some changes, it might be easier
> to first generate standard code for the geometry you want to run with
> (--save_src=code.cu), then modify code.cu with your changes, and
> finally run a simulation with the modified code (--use_src=code.cu).
> This way you will not have to deal with the code of the multitude of
> LB models that Sailfish supports and that is sometimes mixed into the
> standard BGK path (you will see that in BGK_relaxate).
>
> Hope this helps and please let us know if run into any problems,
> Michal
>

[Michal Januszewski]

Hi Luca,

If you want to set a constant velocity profile, presumably at the
outlet or at the inlet, have you considered just using the velocity
boundary conditions?

Michal

--
Michal Januszewski
http://people.gentoo.org/spock

[Luca]

Hi Michal,

thanks, I would like to set a velocity profile within the domain, say
a gaussian for instance..

I initializate density by a distribution, and vx and vy with a known
profile. This velocity profile should keep constant during simulation
carrying density to steady state..

By the way, I am trying with your first suggestion editing the code.cu
in "perform the relaxation step in the BGK model given the density
rho, the velocity v and the distribution fi".. Do you think I am doing
correctly?

Thanks!
Luca

On 25 nov, 10:27, Michal Januszewski <mich...@gmail.com> wrote:
> Hi Luca,
>

[Michal Januszewski]

Hi Luca,

If you need this velocity profile to be specified inside the domain,
then the standard velocity boundary conditions might not work
correctly. You're probably right about just modifying the code
directly, but without seeing the changes you're making I can't tell
you if it's likely to work or not.

Michal

--
Michal Januszewski
http://people.gentoo.org/spock

[Luca]

Hi Michal,

thanks for your willingness in helping me. I have initialised density
by a distribution function in the domain. I have also initialised vx
and vy by a curve. Now, I would like to solve density field to steady
state while leaving constant in time the velocity profile I imposed as
initial condition. I am not using any boundary condition in the domain
(in this case might be fine).

I think the solution to my problem might be in this part of my
code.cu:

------------------------------------------------------------------------------------------

//
// Performs the relaxation step in the BGK model given the density
rho,
// the velocity v and the distribution fi.
//
__device__ inline void BGK_relaxate(float rho, float *iv0, Dist *d0,
int node_type, int ncode)
{
Dist feq0;
float v0[2];

// Body force acceleration.

v0[0] = iv0[0] + 0.00000000000000000000e+00f;
v0[1] = iv0[1] + 0.00000000000000000000e+00f;
feq0.fC = 4*rho/9 + 4*rho*(-3*v0[0]*v0[0]/2 - 3*v0[1]*v0[1]/2)/9;
feq0.fE = rho/9 + rho*(3*v0[0]*(1 + v0[0]) - 3*v0[1]*v0[1]/2)/9;
feq0.fN = rho/9 + rho*(3*v0[1]*(1 + v0[1]) - 3*v0[0]*v0[0]/2)/9;
feq0.fW = rho/9 + rho*(-3*v0[0]*(1 - v0[0]) - 3*v0[1]*v0[1]/2)/9;
feq0.fS = rho/9 + rho*(-3*v0[1]*(1 - v0[1]) - 3*v0[0]*v0[0]/2)/9;
feq0.fNE = rho/36 + rho*(3*v0[0]*(1 + v0[0]) + 3*v0[1]*(1 + v0[1] +
3*v0[0]))/36;
feq0.fNW = rho/36 + rho*(-3*v0[0]*(1 - v0[0]) + 3*v0[1]*(1 + v0[1] -
3*v0[0]))/36;
feq0.fSW = rho/36 + rho*(-3*v0[0]*(1 - v0[0]) - 3*v0[1]*(1 - v0[1] -
3*v0[0]))/36;
feq0.fSE = rho/36 + rho*(-3*v0[1]*(1 - v0[1] + 3*v0[0]) + 3*v0[0]*(1 +
v0[0]))/36;
d0->fC += (feq0.fC - d0->fC) / tau0;
d0->fE += (feq0.fE - d0->fE) / tau0;
d0->fN += (feq0.fN - d0->fN) / tau0;
d0->fW += (feq0.fW - d0->fW) / tau0;
d0->fS += (feq0.fS - d0->fS) / tau0;
d0->fNE += (feq0.fNE - d0->fNE) / tau0;
d0->fNW += (feq0.fNW - d0->fNW) / tau0;
d0->fSW += (feq0.fSW - d0->fSW) / tau0;
d0->fSE += (feq0.fSE - d0->fSE) / tau0;

// Body force acceleration.

// FIXME: This should be moved to postcollision boundary conditions.
iv0[0] += 0.00000000000000000000e+00f; iv0[1] +=
0.00000000000000000000e+00f;}

---------------------------------------------------------------------------------------------------

My idea should be to not to update velocities but keep them constant,
reading always the initial ones I set. Sorry for the roughness of my
knowledge, maybe the answer is straightforward but I am having some
problems to figure it out as I am just a beginner in the field...

Thanks!
Luca

On 27 Nov, 21:09, Michal Januszewski <mich...@gmail.com> wrote:
> Hi Luca,
>

[Luca]

... or maybe when computing moments:

---------------------------------------------------------------------------------------

// Compute the 0th moment of the distributions, i.e. density.
__device__ inline void compute_0th_moment(Dist *fi, float *out)
{
*out = fi->fC + fi->fE + fi->fN + fi->fNE + fi->fNW + fi->fS + fi->fSE
+ fi->fSW + fi->fW;
}

// Compute the 1st moments of the distributions, i.e. momentum.
__device__ inline void compute_1st_moment(Dist *fi, float *out, int
add, float factor)
{
if (add)
{
out[0] += factor * (fi->fE + fi->fNE + fi->fSE - fi->fNW - fi->fSW -
fi->fW); //orig +=
out[1] += factor * (fi->fN + fi->fNE + fi->fNW - fi->fS - fi->fSE - fi-
>fSW); //orig +=
}
else
{
out[0] = factor * (fi->fE + fi->fNE + fi->fSE - fi->fNW - fi->fSW - fi-
>fW);
out[1] = factor * (fi->fN + fi->fNE + fi->fNW - fi->fS - fi->fSE - fi-
>fSW);
}
}

// Compute the 2nd moments of the distributions. Order of components
is:
// 2D: xx, xy, yy
// 3D: xx, xy, xz, yy, yz, zz
__device__ inline void compute_2nd_moment(Dist *fi, float *out)
{
out[0] = fi->fE + fi->fNE + fi->fNW + fi->fSE + fi->fSW + fi->fW;
out[1] = fi->fNE + fi->fSW - fi->fNW - fi->fSE;
out[2] = fi->fN + fi->fNE + fi->fNW + fi->fS + fi->fSE + fi->fSW;
}

// Compute the 1st moments of the distributions and divide it by the 0-
th moment
// i.e. compute velocity.
__device__ inline void compute_1st_div_0th(Dist *fi, float *out, float
zero)
{
out[0] = (fi->fE + fi->fNE + fi->fSE - fi->fNW - fi->fSW - fi->fW)/
zero;
out[1] = (fi->fN + fi->fNE + fi->fNW - fi->fS - fi->fSE - fi->fSW)/
zero;
}

__device__ inline void compute_macro_quant(Dist *fi, float *rho, float
*v)
{
compute_0th_moment(fi, rho);
compute_1st_div_0th(fi, v, *rho);
}

__device__ inline void get0thMoment(Dist *fi, int node_type, int
orientation, float *out)
{
compute_0th_moment(fi, out);
}

//
// Get macroscopic density rho and velocity v given a distribution fi,
and
// the node class node_type.
//
__device__ inline void getMacro(Dist *fi, int ncode, int node_type,
int orientation, float *rho, float *v0)
{
if (isFluidOrWallNode(node_type) || isSlipNode(node_type) ||
orientation == 0)
{
compute_macro_quant(fi, rho, v0);
if (isWallNode(node_type))
{
}
}
else if (isVelocityNode(node_type))
{
compute_macro_quant(fi, rho, v0);
}
else if (isPressureNode(node_type))
{
compute_macro_quant(fi, rho, v0);
}
}

Luca

[Michal Januszewski]

Hi Luca,

If you want to keep the velocity profile constant in the whole domain,
then are you sure that your system is modeled by the Navier-Stokes
equation? It seems to me that perhaps the advection-diffusion
equation might be a better match. If so, you should look up a
consistent LB model for this equation. It would then have to be added
to Sailfish by introducing a new equilibrium function. I'd be happy
to help you with that.

If this is not what you want, then overriding the contents of v0 in
getMacro should provide a constant velocity field, towards which the
model will then be relaxed by the standard relaxation procedure which
uses v0 to compute the equilibrium distribution.

Cheers,
Michal

--
Michal Januszewski
http://people.gentoo.org/spock

[Alexandr Kuzmin]

Dear all,

If you fix velocity you indeed model the advection-diffusion equation and you don't need to change equilibrium  - it will work like that.

Cheers,

Alex

[Michal Januszewski]

Hi Alex,

Thanks for clarifying that! In this case, the getMacro approach
mentioned earlier should work. An alternative could be to apply the
following patch:

diff --git a/sailfish/templates/single_fluid.mako
b/sailfish/templates/single_fluid.mako
index 1bd1fa4..430d839 100644
--- a/sailfish/templates/single_fluid.mako
+++ b/sailfish/templates/single_fluid.mako
@@ -168,17 +168,16 @@ ${kernel} void CollideAndPropagate(
%endif

precollisionBoundaryConditions(&d0, ncode, type, orientation, &g0m0, v);
+
+ v[0] = ovx[0];
+ v[1] = ovy[1];
+
${relaxate(bgk_args)}
postcollisionBoundaryConditions(&d0, ncode, type, orientation,
&g0m0, v, gi, dist_out);

// only save the macroscopic quantities if requested to do so
if (save_macro == 1) {
orho[gi] = g0m0;
- ovx[gi] = v[0];
- ovy[gi] = v[1];
- %if dim == 3:
- ovz[gi] = v[2];
- %endif
}

${propagate('dist_out', 'd0')}

This will:cause:
- the velocity field to never be updated,
- the initial velocity field to always be used for relaxation.

Hope this helps,
Michal

[Luca]

Hi all,

Michal, it is indeed an advection-diffusion equation, but as Alex
pointed out, there is no need to change equilibrium, I got it simply
fixing velocity (similarly to your suggestion) and I can confirm that
it works.

Thank you all, cheers!
Luca

> > On Sun, Dec 4, 2011 at 5:53 PM, Michal Januszewski <mich...@gmail.com>
> > wrote:
>
> >> Hi Luca,
>

> >> If you want to keep the velocity profile constant in the whole domain,
> >> then are you sure that your system is modeled by the Navier-Stokes
> >> equation?  It seems to me that perhaps the advection-diffusion
> >> equation might be a better match.  If so, you should look up a
> >> consistent LB model for this equation.  It would then have to be added
> >> to Sailfish by introducing a new equilibrium function.  I'd be happy
> >> to help you with that.
>
> >> If this is not what you want, then overriding the contents of v0 in
> >> getMacro should provide a constant velocity field, towards which the
> >> model will then be relaxed by the standard relaxation procedure which
> >> uses v0 to compute the equilibrium distribution.
>
> >> Cheers,
> >> Michal
>