LBMethod.org forums

Re: D3Q15 vs D3Q19 vs D3Q27

daralcan — Wed, 08 Sep 2010 18:04:37 +0200

Hi everyone,

Do you suggest any article about comparison of 3D lattice links?

Thank you

Daralcan

Re: Compiling Error on defineDynamics

hyoon — Tue, 07 Sep 2010 16:03:03 +0200

Hi,

Thanks for the examples.

In the Martys and Chen model, there are two potential terms for fluid-fluid and fluid-solid. Since we can assign a density value to the bounce back node, to have the Martys and Chen model, we just need to use two potential fns for fluid nodes without contact with bounceback node and fluid nodes with contact with bounceback nodes in lattice, respectively. I need to study the data process to understand the code structure better.

Hong

boussinesqThermal2d

giordy75 — Mon, 06 Sep 2010 18:09:24 +0200

Dear All.

Someone know how to change the boussinesqThermal2d solver to simulate the Natural convection in a square cavity (pag338, [rcpt.yousendit.com])?

Thanks in advance.
Stefano

Bed Elevation?

mathloverboy — Mon, 06 Sep 2010 08:55:02 +0200

I am trying to apply LBE to shallow water flow and I want to test it for a case where there is a small bump in a channel for the subcritical flow. Does anyone know how should I introduce the bed elevation in the model? I am using the sample FORTRAN code provided by Zhou at the end of his book.

Bounca back

r.imani — Sat, 04 Sep 2010 13:47:29 +0200

Hi dear all
I decided to simulate heat and fluid flow in a 2-d channel with a single square obstacle placed in it. at the first stage of my simulation I simulated the hydrodynamic using d2q9 with mid-plane bounce back but I have velocity on solid boundaries surely, this non-zero velocities can affect the heat transfer problem in the second stage of my simulation so i have 2 questions:
1- how much it will be accurate if we force the velocities on the solid boundaries to be zero right after simulation ?
2- are the half-way bounce back is significantly accurate than mid-plane bounce back?

thanks in advance

Re: About integrateProcessingFunctional

jlatt — Fri, 03 Sep 2010 17:20:45 +0200

Hi,

The effect of the Shan/Chen data processors (both the single- and the multi-component version) is to compute the density and the momentum for every component, and assign it to the lattice. In your example, the results of the single-component model will simply be overwritten by the results of the multi-component version, which is probably not what you want.

About integrateProcessingFunctional

Glary — Fri, 03 Sep 2010 17:11:43 +0200

Hi,

I would like to do some modify in the multicomponent2d.cpp code. I want to add the interparticle force for each component except the ShanChen Multicomponent interaction, which is already in the code. Thus I did it like below. Am I right? Can I use three interfrateProcessingFunctional in the same time? Thank you.

integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D (
                G11, new interparticlePotential::PsiShanChen93(rho0) ), lattice.getBoundingBox(), latticeOne, processorLevel );

integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D (
                G22, new interparticlePotential::PsiShanChen93(rho0) ),lattice.getBoundingBox(), latticeTwo, processorLevel );

integrateProcessingFunctional (new ShanChenMultiComponentProcessor2D(G12,constOmegaValues), Box2D(1,nx-2,0,ny-1),
            blockLattices, processorLevel );

Re: Few simple questions

malaspin — Fri, 03 Sep 2010 13:52:08 +0200

1) What do you mean by height of fluid. I guess that you implemented a shallow water equation or something similar. Otherwise you do not have any "height of watter " out of the box. If I am not mistaken in the shallow water equation the height of water is "equivalent" with the density in palabos. So initializing the density at the height you want is as simple as initializing the density.

2) You can use BoxProcessingFunctional2D instead. And then use integrateProcessingFunctional. But be aware that the DataProcessors integrate to the lattice with integrateProcessingFunctional are executed after the collision and streaming steps. This does not seem to be what you want since your DataProcessor is executed after collision.

3) You cannot extract dx and dt directly from the lattice since it contains no information about the "physical" parameters of your flow. But if you want to use these parameters in a DataProcessor you can just pass them as members of the DataProcessor and initialize them with the constructor.

4) There is no way to disable stream . The only way to do it is as you mention to enable it only on certain parts of the domain.

5) On the boundaries of the domain, if the periodicity is "off" then by default the distribution functions that are incoming into the domain are doing a "half-way" bounceback. If periodicity is "on" then the f_i are juste streamed periodically.

I hope it helps.

Re: Compiling Error on defineDynamics

jlatt — Fri, 03 Sep 2010 11:25:17 +0200

Hi,

Are you referring to the single-component multi-phase Shan/Chen model in which the interaction potential is psi(rho)=rho? In this case it is already implemented in Palabos. Have a look at the example program codesByTopic/shanChenMultiPhase/segregation2D.cpp (the example is 2D, but the same works in 3D). This example uses the ShanChenSingleComponentProcessor, which takes as one of the arguments a function for the interaction potential. In the example, the chosen potential is interparticlePotential::PsiShanChen94(psi0,rho0). You can replace this by interparticlePotential::PsiIsRho() to have the Martys and Chen model.

Have a look at the file multiPhysics/interparticlePotential.h to get a list of implemented potentials, and to see how to implement new potentials.

Yes, please do share your thoughts and results with us, we are looking forward to gathering more knowledge about these models.

Re: about relaxation time and kinematic viscosity

Holzdolph — Thu, 02 Sep 2010 16:24:48 +0200

So far I have only used the simple BGK-dynamics model with only one component, so it may be not exactly the same in your case.
But the principle idea should be the same, so If you take a look at tau:

tau = nu_lb / c + 1/2

with nu_lb = dt/dx^2 * nu_d and nu_d = 1/Re

tau = 3 * dt/dx^2 * 1/Re + 1/2

Now if you fix the values for tau (=1), and Re (or nu_d which comes from nu_ph), choose a value for dx (or the number of cells N, with dx=1/N), you are left with picking the right timestep dt, and nu_ph has still the value you wanted it to have.
Of course this usually means longer calculation time, but with the values you need.

nu_lb ... viscosity in LB-units
nu_d.... dimensionless viscosity
nu_ph... physical visc.
Re... Reynoldsnumber
dt.... timestep LB
dx... grid spacing LB

hope that helps a bit^^

Re: Simulating a canal with sliding gate

makhan — Thu, 02 Sep 2010 13:12:39 +0200

Ah sorry, I didn't update the information. I don't get the segmentation fault anymore. I had it when I tried to use one gate may times, or tried to delete it. I think palabos is taking care of all the pointers that I pass as arguments, so it's ok.
I already switched from dynamics to data processors though.

I'm trying to write a data processor that'd connect two separate lattices, but with no success. Don't know how to turn off boundary condition completely, it always gets bounced back.
I decided this topic is way too messy and created another - [www.lbmethod.org]

Re: Simulating a canal with sliding gate

malaspin — Thu, 02 Sep 2010 12:14:41 +0200

Hello,

You never delete your gleft and gright pointers which may be a problem. A you sure it is not related with your segmentation fault?

Re: Compiling Error on defineDynamics

hyoon — Thu, 02 Sep 2010 05:51:46 +0200

Hi Jonas,

Thanks for correcting the problem. Type works and you're right -- no reason for double type of geometry.

I have one more question related to my current work.
In coming weeks, I plan to update Shan-Chen model by adding fluid-surface forces developed by Martys and Chen (1996).
This method has been already used for many cases and demonstrated well in Sukop and Thorne's textbook.

I plan to dig out how the current Shan-Chen model is coded and follow the current data processing.
I wonder if you have a good example or some tips for this.

Hopefully I can share my work here later.

Hong

Re: stuck with own data processor

Holzdolph — Wed, 01 Sep 2010 12:14:31 +0200

Yes I´ve started with that tutorial to write my functional. And I think the functional i wrote is OK, the only problem that remains is to get my data :)

That parallel programming is still something I´ve to get used to:
In the example with the average energy, you use some BlockStatistics class to calculate the energy (for each process) and then somehow sum it up.....
and since i only use one simple variable inside my class to store the value i want, of course it only works for a single process, since it gets split up in parallel mode....

edit:
So, I cant say that I completely understand how that BlockStatistics class is working, but it´s doing what I want right now.
Here´s the final code:

class MyFunctional :
    public ReductiveBoxProcessingFunctional3D_S
{
public:
    MyFunctional();
    virtual void process(Box3D domain, ScalarField3D& matrix);
    virtual MyFunctional* clone() const;
    plint getArea() const;
private:
    plint sumEnergyId;
    plint area;
};

MyFunctional::MyFunctional() : area(this->getStatistics().subscribeSum()) { 

}

void MyFunctional::process (Box3D domain, ScalarField3D& matrix )

{
  plint subarea=0;
  BlockStatistics& statistics = this->getStatistics();
    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 (matrix.get(iX, iY, iZ)) subarea++;
            }
	}
    }
statistics.gatherSum(area, subarea);
}

MyFunctional* MyFunctional::clone() const
{
    return new MyFunctional(*this);
}

plint MyFunctional::getArea() const {
  return this->getStatistics().getSum(area);
  
}


plint computeArea(MultiScalarField3D& matrix, Box3D& domain)  {
  MyFunctional*  functional;
  functional = new MyFunctional();
  applyProcessingFunctional(*functional, domain, matrix);
  return functional->getArea();
}

Re: Compiling Error on defineDynamics

jlatt — Wed, 01 Sep 2010 11:28:02 +0200

Hello,

Yes, we made a little change it the user's interface, sorry for that (the change made sense, though, as you will see in the following). As of version 0.7, the geometry argument should no longer be of type MultiScalarField2D, but of type MultiScalarField2D.

Of course, "int" makes more sense than "double," because "geometry" is a flag matrix. In old versions of Palabos, it was unfortunately not possible to mix data types, which means, it was not possible to bring a MultiBlockLattice2D in contact with a MultiScalarField2D if T1 is different from T2. The workaround was to have the flag matrix to be of type "double". This was a bit awkward and has been fixed: from now on flag matrices are simply of type int.

Re: stuck with own data processor

jlatt — Wed, 01 Sep 2010 11:16:51 +0200

Hi,

The philosophy of Palabos is that, while the end-user works with MultiBlocks, the data-processors are applied individually to each AtomicBlock inside the MultiBlock. That's how parallelism works in Palabos, and it is explained in the user's guide. There is even a tutorial that investigates this point in detail.

Re: stuck with own data processor

Holzdolph — Wed, 01 Sep 2010 10:33:20 +0200

Hi,

you were right, taking the existing code of a functional and changing it is the best way, I finally got it to work :)

Although, two problems still remain (were one of them I´m sure I´ll solve on my own), about the other one I´m not so sure:

The thing is, the ReductiveBoxProcessingFunctional3D_S´s virtual process function requires a ScalarField, and not a MultiScalarField. And I didnt find any other ProcessingFunctional which has a MultiScalarField in its process function.

So if I use a ScalarField, the code compiles and works - but not in parallel mode. Which is the point of the DataProcessor in the first place, to work in mpi-mode :(

Thats my Functional:

class MyFunctional :
    public ReductiveBoxProcessingFunctional3D_S
{
public:
    MyFunctional();
    virtual void process(Box3D domain, ScalarField3D& matrix);
    virtual MyFunctional* clone() const;
    plint getArea() const;
private:
    plint sumEnergyId;
    plint area;
};

MyFunctional::MyFunctional() { 
area = 0;
}

void MyFunctional::process (Box3D domain, ScalarField3D& matrix )

{
    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 (matrix.get(iX, iY, iZ)) area++;
                pcout << matrix.get(iX, iY, iZ);
            }
	}
    pcout << endl;
    }
pcout << area << endl;
}

MyFunctional* MyFunctional::clone() const
{
    return new MyFunctional(*this);
}

plint MyFunctional::getArea() const {
  return area;
}

plint computeArea(MultiScalarField3D& matrix, Box3D& domain)  {
  MyFunctional  functional;
  applyProcessingFunctional(functional, domain, matrix);
  return functional.getArea();
}

(note: those pcouts are just for testing)

the ReductiveBoxProcessingFunctional3D_S in reductiveDataProcessingFunctional3D.h looks like:

/

// Easy instantiation of boxed data processor for a single scalar field
template
struct ReductiveBoxProcessingFunctional3D_S : public PlainReductiveBoxProcessingFunctional3D {
    virtual void process(Box3D domain, ScalarField3D& field) =0;
    /// Invoke parent-method "processGenericBlocks" through a type-cast
    virtual void processGenericBlocks(Box3D domain, std::vector atomicBlocks);
};

Can I just make another one, with a MultiScalarField (and does it work with an AtomicBlock or shouldn´t it be a MultiBlock). Well thats the question, and I´ll try it, but I think it´s not going to be that simple ?

Re: Few simple questions

makhan — Wed, 01 Sep 2010 10:22:09 +0200

5) As I define boundary condition using data processor, what's the proper way to disable any bounce back that occurs? Right now, I'm running the processor after collision phase and then apply streaming, but judging by the results, some default boundary condition is applied to these cells anyway (I defined BounceBack dynamics at the top and bottom wall, but I can see the water being bounced back also at inlet and outlet).

Compiling Error on defineDynamics

hyoon — Tue, 31 Aug 2010 20:03:35 +0200

Hello,

I tried to use ShanChen model to simulate single component liquid-vapor flow.
I modifed the sample code - segregation2D and tutorial code - permeability.
I defined lattice and geometry as in permeability.
I compiled the code using version 6 and the latest version.
With version 6, it was ok, but with the latest version (I complied on Cygwin),
I've got the error message as below:

error: no matching function for call to 'defineDynamics(plb::MultiBlockLattice2D&,
plb::MultiScalarField2D&, plb::BounceBack*, int)

Here is my code:
defineDynamics(lattice, geometry, new BounceBack(65.), 1);

the same error on defineDynamics(lattice, geometry, new NoDynamics(), 2);

It was complied correctly with permeability.cpp where similar call was used.

Any help will be appreciated.

Thanks,

Hong

Few simple questions

makhan — Tue, 31 Aug 2010 11:46:38 +0200

My other topic grew too large already so I'll start fresh here. Especially since I switched from dynamics to data processors.
Will number my questions to make it easier to answer.

1) What's the proper way to initialize the height of fluid? InitializeAtEquilibrium with rho parameter seems like initializing the density and I want the density throughout the whole lattice to be constant. Want to initialize and monitor the water level.

2) I'm using OneCellIndexedFunctional2D. What's the proper way to run it in each iteration between collision and streaming phase? Now I'm using lattice.collide(), applyIndexed(lattice, box, new Functional()), lattice.stream(). Is there a way to reach the same result with single function call? Avoid creating new instance in each loop?

3) I can get dx and dt using parameters.getDeltaX() and parameters.getDeltaT(), where parameters are of type IncomprFlowParam. Can I (how) extract parameters from lattice (Multi)BlockLattice2D?

4) Is there a way to permanently disable streaming between two columns? Or do I have to use stream with Box2D parameter defining the domain?

Re: Non newtonian flow

Timm — Tue, 31 Aug 2010 10:47:09 +0200

Hello,

there are a few references. I think that some of them are also discussed here. Otherwise you should search for 'lattice boltzmann non-newtonian' in Google Scholar. You will find plenty of literature.

Timm

Non newtonian flow

piyush — Tue, 31 Aug 2010 09:52:45 +0200

Hi
I'm a ChemE undergrad, final year. As a project I'm modifying a D3Q19 LBM implementation to include the modellng of non newtonian fluids, power law to start with.
As I understand it, what I'd essentially need to change is the collision function, which currently implements the standard BGK model.
Is there a well known formulation of such a thing? Some references would be great..

Thanks
Piyush

Re: time dependent external force

Glary — Mon, 30 Aug 2010 13:49:30 +0200

Hi jlatt,

Thank you for your quick reply and advise. I will upgrade to the new version. If I find further questions in the new version, I will come back. Thank you very much for your help.

Re: time dependent external force

jlatt — Mon, 30 Aug 2010 13:43:52 +0200

Hi,

In the version 0.6 you are using, there exists unfortunately no convenient way of changing the external force in a given domain. The only solution in this version is to write your own data processor. This is not difficult, as you can take example on other data processors in the multi-component code, but it is definitely not convenient.

I much more recommend that you upgrade to version 0.7. While this upgrade may require a few minor adaptations in your end-user application, this should be quick, especially if you get inspiration from the example programs.

In version 0.7, you can change the exernal force at any time by calling something like

Array externalForce(0.,0.001);
setExternalVector(lattice,lattice.getBoundingBox(),DESCRIPTOR::ExternalField::forceBeginsAt,force);

Re: Simulating a canal with sliding gate

makhan — Mon, 30 Aug 2010 12:37:40 +0200

Ok, I learned that I can't access coordinates from boundary condition dynamics, so I switched to data processors.

I wrote a simple OneCellIndexed functional that does nothing and am trying to run it with each iteration, but I keep getting segmentation fault. Is there any method that allows me to define an indexed functional that'd be run with each iteration? (something like integrateProcessingFunctional?) Or should I create a new functional object for each iteration?

Another thing I'd like to do is to disable streaming between 2 columns (in this case nx/2 and nx/2+1). Can it be done? Or should I define 2 separate lattices and apply functional to the end of one and beginning of another?

#include "palabos2D.h"
#include "palabos2D.hh"
#include 
#include 

using namespace plb;
using namespace std;

typedef double T;

#define DESCRIPTOR descriptors::D2Q9Descriptor

template class Descriptor>
class Gate : public OneCellIndexedFunctional2D {
public:
    Gate(plint ny_, bool left_)
        : ny(ny_),
          left(left_)
    { }
    Gate* clone() const {
        return new Gate(*this);
    }
    virtual void execute(plint iX, plint iY, Cell& cell) const {
//	T h = 0.;
//	for (plint i=0; i<9; i++)
//		h += cell[i];
//	cout << "Height in (" << iX << "," << iY << ") = " << h << endl;
    }
private:
    plint ny;
    bool left;
};

Gate *gleft, *gright;


void canalSetup( MultiBlockLattice2D& water,
                          T rho0,
                          T force )
{
    int nx = water.getNx();
    int ny = water.getNy();
    
    // Bounce-back on bottom wall
    defineDynamics(water, Box2D(0,nx-1, 0,0), new BounceBack(rho0) );
    // Bounce-back on top wall
    defineDynamics(water, Box2D(0,nx-1, ny-1,ny-1), new BounceBack(rho0) );



	
    // External force
//    setExternalVector(water, water.getBoundingBox(),DESCRIPTOR::ExternalField::forceBeginsAt, Array(0.,-force));

    water.initialize();
}

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");
    
    const T omega = 1.0;
    const int nx   = 1200;
    const int ny   = 400;
    T force        = 0.15/(T)ny;
    T rho0 = 0.;
    const int maxIter  = 16000;
    const int saveIter = 100;

    MultiBlockLattice2D water (nx,ny, new BGKdynamics(omega) );

    water.periodicity().toggleAll(false);

    canalSetup(water, rho0, force);
	
    pcout << "Starting simulation" << endl;
    for (int iT=0; iT imageWriter("leeloo");
			imageWriter.writeScaledGif(createFileName("u", iT, 6), *computeVelocityNorm(water));
        }

        water.collide();
        gleft = new Gate(ny,true);
        gright = new Gate(ny,false);
	applyIndexed(water, Box2D(nx/2, nx/2, 0, ny-1), gleft);
	applyIndexed(water, Box2D(nx/2+1, nx/2+1, 0, ny-1), gright);
        water.stream();
    }
}

time dependent external force

Glary — Sat, 28 Aug 2010 16:30:04 +0200

Hi,

I would like to add a time dependent ExternalForce in the multicomponent code in palabos-0.6r2. (I know there is a new version, but I am sorry what I used is the old version.) I found an opption about set External Force in the core/dynamic.hh file. But I didnot used it successfully. Please give me a help. The command I used to set the time-dependent external force is as below

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");
.........

     // Main loop over time iterations.
    for (int iT=0; iT 20){
       Array externalforce (0.,0.001);
       setExternalField (insideFluid, Box2D(0, nx-1, 0, ny-1), 1, &externalforce[0]);
     }

the error message is
rayleighTaylor2D.cpp:322: error: `setExternalField' undeclared (first use this function)

How to deal with this problem? I really appreciate for your help! By the way what is the meaning of size in

void Dynamics::setExternalField (
        Cell& cell, plint pos, plint size, const T* ext)

Thank you in advance.

Re: constant density at inlet

hyoon — Fri, 27 Aug 2010 19:39:00 +0200

Hi Cotrf,

I wonder if you figured out what problem was. I will try to work on wetting porous media this weekend.

Please let me know if you have any progress. I will update mine after the weekend if I will get any meaningful resutls.

Thanks,

Hong

Re: scaling in CUDA !!

n0thing — Thu, 26 Aug 2010 14:36:34 +0200

Hi Alex,

I am trying to implement an efficient LBM simulation in OpenCL for ATI GPUs. Can you please send the code to me too?
My mail ID - priyadarshi@darshan3d.com

Thanks
Priyadarshi

Re: Simulating a canal with sliding gate

makhan — Thu, 26 Aug 2010 14:13:13 +0200

I copied ZouHe* files, changed the names, classes and added necessary includes, but to avoid confusion, lets just assume I'm editing ZouHe classes.

How can I pass an argument to the ZouHeVelocityDynamics constructor? In the program, I create OnLatticeBoundaryCondition2D object using createZouHeBoundaryCondition2D function.
createZouHeBoundaryCondition2D is defined in zouHeBoundary2D.hh file as "return new BoundaryConditionInstantiator2D"
that's where I get lost, because BoundaryConditionInstantiator2D class is defined in boundaryInstantiator2D.h and it doesn't contain any constructor.
And the only predecessor is OnLatticeBoundaryCondition2D defined in file boundaryCondition2D.h which doesn't have any constructor either.

Obviously I missed something, cause the program works and constructor body is evaluated. Could you help me find out what?

I'll look into the global definitions for now, but I'll need to pass some arguments to the constructor anyway.

---edit---
I just encountered another problem. To calculate new f values, I need to know the height difference on both sides of the gate. The easiest way would be to store all these values in the VelocityDynamics class, but completePopulations function is constant, so the access is read-only.
Even if I used the global variable, I would need to know the exact coordinates of the cell. Can I access these informations?

Because both sides of the gate are on different lattices, the best way would be to store pointer to the other side of the gate in each of them, get information on the actual cell position and access appropriate cell in the other lattice to read the distribution values. Is it possible?

Or is it not possible and I should use data processor instead?

Re: how to modify multiComponent3d

motti — Thu, 26 Aug 2010 11:56:59 +0200

Hi all,

I found the Units.h in the directory CORE.
and I changed the LatticeNu with changing LatticeU,Re,resolution.
and Run the code. then the sphpere was going to get smaller and Vanishing at last.
it happens both low latticeNu and High latticeNu.so,
I'd like to control the behavior.
and one more things G Variable is what?
give me any idea.

cheers