Convection-diffusion equation

[Vítězslav Žabka]

Hi,
I'd like to solve a convection-diffusion equation simultaneously with
the flow equations. I suspect this isn't supported by Sailfish
currently. Are there any plans that such functionality will be added
soon? Or could you please give me some guidelines how to implement it
by myself?

Thanks,
Vita

[Michal Januszewski]

Hi Vita,

There are no immediate plans to add this functionality.

IIRC, to model the advection-diffusion equation one only needs a new
set of distributions which is relaxed using the standard equilibrium
distribution function and the velocity field from the fluid
simulation. If this is indeed what you want to do, it should be quite
simple and I'd be happy to provide some code pointers. If you want to
do something else, could you please provide a rough overview at the
level of the LB model so that I can understand better what needs to be
done? In any case, I'd be happy to work with you on getting this
implemented and subsequently added to our main repository.

Cheers,
Michal

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

[Vítězslav Žabka]

Thank you for your answer. You are correct. This is exactly what we want to do. We simulate 2D turbulent flow using the BGK model and LES-Smagorinsky subgrid model. I believe that the simplest way to model the advection-diffusion equation is to use four directions and the BGK approximation. We don't need any body force or boundary conditions other than full-way bounce-back.

It would be nice if the new code made it to the main repository.

Vita

[Michal Januszewski]

Hi Vita,

Just to clarify, the 2D turbulent flow would also be simulated using Sailfish?

If this is the case, the implementation could take the following form:

• in lb_single.py, you declare a new simulation class LBAdvectionDiffusion(LBFluidSim)
• in this class you need to override some methods to add a second lattice for the scalar field; this will follow the same pattern as lb_binary:LBBinaryFluidBase.  In particular, you will need to:
• override __init__ to add a second grid (you will need to add D2Q5 to sym.py if you want just 4 directions, or just use D2Q9, which would be simpler for a first version)
• override get_compute_kernels so that the new set of distributions is passed to the simulation; it might also make sense to use completely new kernels for this type of simulation 
• override initial_conditions to set the initial values of the distributions for the new lattice
• in templates/single_fluid.mako, you need to modify both kernels to optionally take the new sets of distributions as arguments, or define completely new kernels (you can start by copy-pasting the existing ones)
• in your new kernel / modified one, you will need to:
• declare a new set of distributions Dist d1 and call getDist to load it from global memory
• call get0thMoment to calculate the density for the scalar field
• modify bgk_args and bgk_args_decl to include this new density
• relaxation and propagation should be done automatically for you

Please let me know if you have any questions about this.

Thanks,

Michal

[Luca]

Hi Michal,

we have already been speaking about advection-diffusion equations.
Modifying the single_fluid.mako template as:


//------------------------------------------------------------------------

${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')}

//------------------------------------------------------------------------

actually works. Nevertheless accessing BGK equilibrium can be useful,
say in my case for example to get rid of second order terms for this
kind of problems. I have read a post where you explain how to deal
with introducing a D2Q5 lattice and a new simulation class for adv-
diff in lb_single.py. By now I am just interested in commenting out
2nd order, can you please point me right in the templates?

Thanks a lot!
Luca

[Michal Januszewski]

Hi Luca,

I think your use case was a little different, as you had a constant
(?) velocity field, while Vita wants to have the fluid simulation
running alongside the advection-diffusion solver. I'm glad the simple
approach worked for you though.

For BGK equilibrium, if you want to change the formula, you can just modify:

https://github.com/sailfish-team/sailfish/blob/master/sailfish/sym.py#L491

If you want to do custom calculations on the device (i.e. modify the
relaxation step), you can put your code somewhere around:

https://github.com/sailfish-team/sailfish/blob/master/sailfish/templates/relaxation.mako#L316

(the BGK equilibrium is available in feq0, which is local Dist struct
on the device).

Hope that helps,
Michal

[Luca]

Hi Michal,

thanks, I have found the equilibrium. Despite I more or less got the
solution with v0.2, I am now stuck with v0.3 that I have on a
workstation. I wanted just to ask you a couple of things and if you
think what I have done until now could be correct. I rapidly sum up
what I want to do/have done:

1. Fix relaxation time: https://github.com/sailfish-team/sailfish/blob/master/sailfish/sym.py#L1022

def relaxation_time(viscosity):
return my_tau #(6.0 * viscosity + 1.0) / 2.0

2. Impose my own dx and dt: in v0.2 I did it respectively in geo.py
(line 118) and lbm.py (line 121), namely:

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

@property
def dx(self):
return my_length/min(self.shape) #1.0/min(self.shape)

-----

@property
def dt(self):
return my_expression*self.geo.dx**2 #self.geo.dx**2

(obviously in my_expression goes tau -> self.get_tau() )

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

3. Modify equilibrium as you suggested above, to correct speed of
sound from 1/sqrt(3) to my_c/sqrt(3) (where my_c is my_dx/my_dt),
then:

https://github.com/sailfish-team/sailfish/blob/master/sailfish/sym.py#L506


for i, ei in enumerate(grid.basis):
t = (grid.weights[i] * (
(rho + rho0 * poly_factorize(
3*ei.dot(grid.v) +
Rational(9, 2) * (ei.dot(grid.v))**2 -
Rational(3, 2) * grid.v.dot(grid.v)))))

becomes (cutting 2-nd order):

my_c = my_dx/my_dt
my_Cs = my_c/(3**0.5)
c_over_Cs2 = my_c/(my_Cs**2)


for i, ei in enumerate(grid.basis):
t = (grid.weights[i] * ((rho + rho0 * c_over_Cs2 *
ei.dot(grid.v))))

4. Modify single_fluid.mako as said above to never update velocity
field and keep that given in sim.vx[:] and sim.vy[:]

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

I think this steps may be sufficient to resolve an advection-diffusion
equation on a double periodic 2D domain (applying no BCs). In
conclusion my problem is that I cannot find definition of dx and dt in
v0.3. Could you please point me there?

Do you think what I have done could be correct?

Thanks a lot!
Luca

[Michal Januszewski]

Hi Luca,

On Wed, Jun 13, 2012 at 4:08 PM, Luca <luca.be...@gmail.com> wrote:
> Hi Michal,
>
> thanks, I have found the equilibrium. Despite I more or less got the
> solution with v0.2, I am now stuck with v0.3 that I have on a
> workstation. I wanted just to ask you a couple of things and if you
> think what I have done until now could be correct. I rapidly sum up
> what I want to do/have done:
>
> 1. Fix relaxation time: https://github.com/sailfish-team/sailfish/blob/master/sailfish/sym.py#L1022
>
> def relaxation_time(viscosity):
>    return my_tau #(6.0 * viscosity + 1.0) / 2.0

Wouldn't it be easier to just provide a corresponding viscosity as a
default option in your simulation? (see e.g.
https://github.com/sailfish-team/sailfish/blob/master/examples/poiseuille.py#L100)

> 2. Impose my own dx and dt: in v0.2 I did it respectively in geo.py
> (line 118) and lbm.py (line 121), namely:
>
> ------------------------
>
> @property
> def dx(self):
>    return my_length/min(self.shape) #1.0/min(self.shape)
>
> -----
>
> @property
> def dt(self):
>    return my_expression*self.geo.dx**2 #self.geo.dx**2
>
> (obviously in my_expression goes tau -> self.get_tau() )

This is fine, but note that these two properties were there for
informative purposes only -- redefining them will not change anything
about the dynamics of the simulation.

You can't find them because they are not there :) The simplest thing
would be to just change equilibrium function directly. You can define
dx and dt locally in bgk_equilibrium if this makes the code cleaner or
is helpful in some other way.

> Do you think what I have done could be correct?

Yeah, I don't see anything wrong with it off the top of my head.

Cheers,

[Michal Januszewski]

Hi Vita,

Replies inline.

On Wed, Aug 8, 2012 at 11:33 PM, Vítězslav Žabka <zab...@gmail.com> wrote:

Hi Michal,

I have been trying to implement the solution of an advection-diffusion equation according to your guidelines, but I'm stuck now. Could you please help me? The flow is also simulated using Sailfish. I'm using the previous version of Sailfish (v0.2 or so).

Please be advised that the code you write for this might require non-trivial porting if you later decide to switch to the current version.

 

What I did is:

• I created a new simulation class LBAdvectionDiffusion derived from lb_single.FluidLBMSim in the same manner as BinaryFluidBase is created. I replaced 'phi' with 'c' (concentration is the unknown of the AD equation) and added ctx['tau_c'] to represent tau in the A-D equation.

Please make sure to also modify the relaxation code in the templates to use this new symbol (tau_c). 

• I decided to use the same grid as for the flow (D2Q9), so I didn't care of adding any new grid at all. Do I need to add another D2Q9 grid, or one grid is enough? And will it work if the LES subgrid is used for the flow and no subgrid model for AD?
  

You don't need to build a new grid class for this, but you should have two sets of distributions in your simulation, one to represent the fluid density/velocity and another one for the concentration, much like the binary fluid free energy model does. Using LES on the first lattice might require some minor hacking, but shouldn't be a problem.

• I added a new distribution dist2 and overrode get_compute_kernels (as it is done in BinaryFluidBase).

Cool, that's exactly what I meant above. 

• I created ad.mako file based on single_fluid.mako and changed the kernels, so that they accept also the new dist2 distributions (similarly to binary_fluid.mako). I'd like to set the diffusivity in AD somehow. Is it sufficient to use tau computed from viscosity for the flow and tau_c (depending on the diffusivity) for AD?

It should be sufficient.  If the diffusivity is constant, you could just calculate tau_c on the host (it looks like you are already doing that).  If you need access to the raw diffusivity value in your kernel, you can just export it in the context dictionary like tau_c.

I can't override initial_conditions (in v0.2 there is no such method). I set the initial value of 'c' in overriden LBMGeo2D.init_fields. Is it what you meant?

Yeah, init_fields is the legacy branch equivalent of initial_conditions.

 

Hope that helps,

Michal

I'm sorry for all the questions. I am new to LBM (I use FEM instead).
Thank you,
Vita

[Vítězslav Žabka]

Hi Everyone,

I managed to make the advection-diffusion simulation work. I'm posting the files. They are based on the legacy branch of Sailfish.

advection_diffusion.py -> sailfish/
advection_diffusion.mako -> sailfish/templates/
lbm_ldc_ad.py -> examples/

Please feel free to contact me if you have any idea to improve the simulation.

Thank you, Michal, for your help.

Regards,
Vita

[Marcin Kostur]

Dear Vita,

It is great.
We will try to take your code and see if it could be integrated with v0.3.
Can you provide any validation example for AD problem? With analytical
solutions.

the best

Marcin

--
Pełnomocnik Rektora
ds. zastosowania nowoczesnych komputerowych
metod kształcenia dla regionalnych
kadr innowacyjnej gospodarki
tel. +48 32 359 1360
http://zft.us.edu.pl/kostur
http://icse.us.edu.pl
http://twing.us.edu.pl
---------------------------

[Alexandr Kuzmin]

Hi guys,

I can send you c++ validation codes for advection-diffusion if you are interested.

Cheers,
Alex

[Vítězslav Žabka]

Dear Marcin,

Unfortunately, I have no example with analytical solutions.

Regards,
Vita