Has anybody experience with the topic? Is it generally possible to implement such a model with Palabos?

Thanks in advance for your help!]]>

Does anyone know if it is possible to implement a free surface field in the case of Octree multigrid? Have you got experience with the topic?

Best regards and thanks for the help!

Raffaele]]>

Language: C++T factor = param.rho * (param.dxFinest * param.dxFinest * param.dxFinest * param.dxFinest) / (param.dtFinest * param.dtFinest); Array<T,3> force = factor*boundaryCondition->getForceOnObject(); forces << (double) (iter * param.dtCoarsest) << " " << force[0] << " " << force[1] << " " << force[2] << std::endl;

My problem are now the values of the forces. They are quite big in comparison with the references (approximately 3x or 4x). I still have no free surface implemented in my code and i would therefore expect smaller values at least for the resistance. I also tried with different geometries, the problem remains the same.

Could the problem come from the conversion from lattice to physical units? Or, since i use a constant pressure as outflow condition, from the pressure reference at the outlet boundary? Would you suggest to use another method to compute the forces in my case?

Thanks in advance for your help!

Raffaele]]>

Thanks]]>

I am dipping my toes at multigrid modeling, and I am currently doing some exploration on the technique to understand it enough to model my problem. My goal is, at the end, to model a channel with a porous wall with a multigrid approach, where the duct is modeled with a coarser lattice and the porous wall with a grid with a much higher resolution.

This is a complicated problem, so I am first trying an easier problem to understand the best approach for this. Hence, I am developing a duct flow with an obstacle, where the grid is refined on the region around and past the obstacle, so the resolution is higher where the turbulent regime forms.

However, it seems that the defineDynamics function, which is used to parse bounce back dynamics to a region, cannot handle multi grid lattices. The following snippet is inside a loop which iterates at the iLevels of the multi grid structure:

Language: C++gridrefining.cpp:179:19: error: no matching function for call to ‘defineDynamics(plb::MultiGridLattice2D<double, plb::descriptors::D2Q9Descriptor>&, plb::Box2D, plb::Box2D&, plb::BounceBack<double, plb::descriptors::D2Q9Descriptor>*)’ defineDynamics(lattice, lattice.getBoundingBox(), ~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ obstacle, new plb::BounceBack<T,DESCRIPTOR>);

By inspecting the Palabos devel guide, I see that defineDynamics is defined only for objects of the type BlockLatticeXD or MultiBlockLatticeXD. Thus, I understand why the direct application of defineDynamics to a lattice of the MultiGridLattice2D type fails.

However, the Palabos users guide, at the Grid Refinement section, states that:

alabos simply generates for you three separate multi-blocks, each of which holds the data at a certain level of grid refinement. [...] Palabos then gathers the three multi-blocks into a single data structure (but, for convenience, you can still access them individually), called “multi-grid” structure

So it seems that the MultiGridLattice2D is a shortcut to parse the entire structure to domain functionals, instead of having to parse each underlying multiblock individually. I believe this is precisely the obstacle I reached.

With all this in mind, how to I correctly parse each sub-multiblock of a multigrid structure, as a MultiBlockLatticeXD object, to parse it to a domain functional that expects this kind of input? With this answered, how's the best way to, then, define an obstacle through the defineDynamics function?

Thanks in advance!]]>

It seems like I posted the same topic in a wrong place. Please pardon me leaving the same notes. The following is the post I made in LBM:algorithms section.

I don't know if the forum is no longer active or not, but I would like to ask a question if anyone can help me writing a data processor.

I would like to have a data processor that spits out a std::vector that can be accessed outside of the data processor. Here is the code I wrote.

Language: C++template<typename T, template<typename U> class Descriptor> class ComputeScatteredRho3D : public BoxProcessingFunctional3D_L<T,Descriptor> { public: ComputeScatteredRho3D(std::vector<std::vector<int> > &locationVector_, std::vector<T> &rhoVector_, plint len_) : locationVector(locationVector_), rhoVector(rhoVector_), len(len_) {} virtual void process(Box3D domain, BlockLattice3D<T,Descriptor> &lattice) { for (plint i=0; i<len; ++i) { plint iX = locationVector[i][0]; plint iY = locationVector[i][1]; plint iZ = locationVector[i][2]; rhoVector[i] = lattice.get(iX,iY,iZ).computeDensity(); } } void getRhoVector(std::vector<T> &vect) { for (plint i=0; i<len; ++i) { vect[i] = rhoVector[i]; } } virtual ComputeScatteredRho3D<T,Descriptor>* clone() const { return new ComputeScatteredRho3D<T,Descriptor>(*this); } void getTypeOfModification(std::vector<modif::ModifT>& modified) const { modified[0] = modif::nothing; } private: std::vector<std::vector<int> > locationVector; std::vector<T> rhoVector; plint len; }; void getScatteredDensity3D(MultiBlockLattice3D<T,AVDES>& lattice, vector<vector<int>>& locationVector, vector<T>& newVector) { plint len = locationVector.size(); ComputeScatteredRho3D<T,AVDES> *aaa = new ComputeScatteredRho3D<T,AVDES> (locationVector, newVector, len); applyProcessingFunctional(aaa, lattice.getBoundingBox(), lattice); aaa->getRhoVector(newVector); delete(aaa); }

locationVector is a vector storing location information of a domain.

I want to have an access to rhoVector outside of the data processor using a function getRhoVector. However, it seems like the rhoVector is destroyed after the function applyProcessingFunctional.

Can anyone give me advice on solving this problem? Any help would be greatly appreciated.

Best,

hjung]]>

I am running the simulation to compute the permeability of a 3D Porous Medium. My sample dimensions are 1232x1197x1197. It seems that it is using a considerable amount of memory as the job submitted to the cluster is rejected. I downsized the sample to 50% and now it is running. I would like to know how to calculate the amount of memory the simulation uses knowing the size of the sample.

Thanks in advance

jleon]]>

Thanks a lot！]]>

In these problems, compressibility effects are present even at low speeds. In such a way that the density is not constant, therefore the pressure gradient can not be represented as a force body.

On the other hand, the relaxation time is not related to the viscosity (by means of Reynolds' number), but to the Knudsen's number, additionally there is a slip condition of the fluid with respect to the wall.

Therefore, I modified a file in Palabos (unit.h), so that the relaxation time is obtained by means of Knudsen's number. Already verify this change, with benchmark problems. Observing how to work with the concept of pressure gradient in the example permeability.cpp, I implemented it to work with the Poiseuille problem and a pressure gradient, but it does not work.

Thanks to everyone for reading, I appreciate any contribution. It can be very short, that is, "read this document" or "looking at the source code that I append".

Thank you very much,

Mauricio

]]>Language: C++

#include"palabos3D.h" #include"palabos3D.hh" #include <cstdlib> #include <vector> #include <cmath> #include <iostream> #include <sstream> #include <fstream> #include <iomanip> using namespace plb; using namespace std; using namespace plb::descriptors; typedef double T; #define DESCRIPTOR D3Q27Descriptor #define DYNAMICS BGKdynamics<T, DESCRIPTOR>(parameters.getOmega()) // // Aca definimos el gradiente de presion class PressureGradient { public: PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_) { } void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const { velocity.resetToZero(); density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX; } private: T deltaP; plint nx; }; void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, IncomprFlowParam<T> const& parameters, OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition, T deltaP) { const plint nx = parameters.getNx(); const plint ny = parameters.getNy(); const plint nz = parameters.getNz(); // // Microchannel';s Geometry //================================================================== Box3D bottom(0, nx-1, 0, 0, 0, nz-1); Box3D top (0, nx-1, ny-1, ny-1, 0, nz-1); Box3D right (0, nx-1, 0, ny-1, nz-1, nz-1); Box3D left (0, nx-1, 0, ny-1, 0, 0); Box3D outlet(nx-1, nx-1, 1, ny-2, 1, nz-2); Box3D inlet (0, 0, 1, ny-2, 1, nz-2); // //boudary==================================================== / boundaryCondition->addPressureBoundary0N(inlet, lattice); setBoundaryDensity(lattice, inlet, (T) 1.); // initializeAtEquilibrium(lattice, inlet, PressureGradient(deltaP, nx) // ); boundaryCondition->addPressureBoundary0P(outlet, lattice); setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2); initializeAtEquilibrium(lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) ); lattice.initialize(); delete boundaryCondition; } template<class BlockLatticeT> void writeVTK(BlockLatticeT& lattice, IncomprFlowParam<T> const& parameters, plint iter){ T dx = parameters.getDeltaX(); T dt = parameters.getDeltaT(); VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx); vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt); vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.); } T computeVmedia ( MultiBlockLattice3D<T,DESCRIPTOR>& lattice, IncomprFlowParam<T> const& parameters , Box3D domain ) { T meanU = computeAverage ( *computeVelocityNorm (lattice, domain) ); T meanUphys = (parameters.getDeltaX() / parameters.getDeltaT()) * meanU; return meanUphys; } int main(int argc, char **argv){ plbInit(&argc,&argv); global::directories().setOutputDir("./tmp/"); IncomprFlowParam<T> parameters( (T) 0.63746, //Velocity at physical units (T) 0.03, //Velocity LB (T) 0.025, //Knudsen';s number (T) 1e-6, // caractheristc Lenght (T) 16., //Resolution (T) 18, //times lenght x (T) 1, //times lenght y (T) 3 //times lenght z ); T gradP = (T)1; const T maxT = (T)39e-6; const T imSave = (T)maxT/(T)10; T deltaP = (T)gradP * parameters.getDeltaX(); writeLogFile(parameters, "Parameters_Microchannel"); plint nx = parameters.getNx(); plint ny = parameters.getNy(); plint nz = parameters.getNz(); MultiBlockLattice3D<T, DESCRIPTOR> lattice( nx,ny,nz, new DYNAMICS); util::ValueTracer<T> convergen(parameters.getPhysicalU(), parameters.getPhysicalLength(), 1e-5); lattice.periodicity().toggleAll(false); OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>(); channelSetup(lattice, parameters, *boundaryCondition, deltaP); plb_ofstream UPvsT("./tmp/VvsT_320.dat"); Box3D profileSectionFinal(nx/2, nx/2, ny/2, ny/2,0, nz - 1);//A lo largo del eje Z Box3D profileSectionUP1(nx/2, nx/2, ny/2, ny/2, nz/2, nz/2); plint counter = 0; T unit_velo = parameters.getDeltaX() / parameters.getDeltaT() / parameters.getLatticeU(); for(plint iT=0; iT*parameters.getDeltaT()<1.01*maxT; ++iT){ counter += 1; double TimE = iT*parameters.getDeltaT(); lattice.collideAndStream(); convergencia.takeValue(getStoredAverageEnergy(lattice),true); if(convergen.hasConverged()){ pcout << "Simulation respect to kinetic energy mean" << endl; pcout << "Simulation';s Time =" << TimE << endl; break; } UPvsT << setprecision(16) << TimE << " " << *multiply(*computeVelocityNorm (lattice, profileSectionUP1),unit_velo) << endl; } T Vmedia = computeVmedia(lattice, parameters, profileSectionFinal); pcout <<"mean velocity ="<<" "<< Vmedia << endl; return 0; }

When I set const pressure value to one side of a box,and other side were set to freeslip .there will be velocity in the box.But the four edges of inlet( the pressureboundary ) have velocity=0,which will make the model error after several steps.

If you meet the samilar problem,we can have a communication.

dongxueyang

email: dxy159753@sina.com]]>

I am modeling a flow with const velocity.In a box, there is one inlet, one outlet and four sidewall. About the boundarycondition, I set like this:

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

Box3D inlet = Box3D(0 ,nx-1,0 , 0 , 0 , nz-1);

Box3D outlet = Box3D(0 ,nx-1,ny-1, ny-1, 0 , nz-1);

Box3D front = Box3D(nx-1,nx-1,1 , ny-2, 0 , nz-1);

Box3D back = Box3D(0 ,0 ,1 , ny-2, 0 , nz-1);

Box3D top = Box3D(1 ,nx-2,1 , ny-2, nz-1, nz-1);

Box3D bottom = Box3D(1 ,nx-2,1 , ny-2, 0, 0 );

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, inlet, boundary::dirichlet );

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, outlet, boundary::outflow );

T vy = getInletVelocity(param.inletVelocityLB, 0, param.startIter);

setBoundaryVelocity(*lattice, inlet, Array<T,3>( (T) 0,vy, (T) 0));

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, bottom, boundary::freeslip );

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, front, boundary::freeslip );

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, back, boundary::freeslip );

boundaryCondition->setVelocityConditionOnBlockBoundaries ( *lattice, top, boundary::freeslip );

initializeAtEquilibrium( *lattice, lattice->getBoundingBox(), param.rhoLoLB, Array<T,3>( (T) 0,vx, (T) 0));

lattice->initialize();

lattice->toggleInternalStatistics(false);

The speed at the boundary advances from the entrance to the exit. After reaching the exit, the speed increases and in turn advances toward the entrance. After several cycles, the overall speed becomes larger until the calculation is wrong.

is there anyone can give me any advice?]]>

I trying to compile the tutorial_2 on a macOS, but I have a issue. My mac has phyton3.x and the Palabos uses previous phyton versions to compile. The error is

"scons: *** SCons version 2.1.0 does not run under Python version 3.6.0. Python 3.0 and later are not yet supported".

Anyone knows how to resolve this issue?.

I did not find the right solution from the internet.

References:-

[www.palabos.org]

B2B Online Video

Thanks]]>

Thank you again for your answers. I will try some further tests.

With best,

Song]]>

1) I do not think so, as the name suggests it is used for domains that are bounded. If you need something like that, I would suggest to use a "normal" BoxProcessingFunctionalXXX and then check the nodetype inside the dataprocessor (via getDynamics() ) and apply different routines depending on the type of node. However, this does not seem to be the optimal solution.

2) Yes, you can of course still use the multiblock data structure.

Regards,

Knut]]>

I really appreciate your answers which are helpful. I am going to have a try.

Here I have another two questions:

- [*] Does the boundedBoxProcessingFunctional work for the periodic boundary condition? I tried to write a functional inherits from boundedBoxProcessingFunctional. However, when I toggle on the periodicity of the boundary, there seem to be no effects. Does this mean that I have to write another functional inheriting from BoxProcessingFunctional with the bulk version process alone?

[*] Thank you for telling me the parameter "vicinity". I was wondering if my non-local operation is global, would it still be possible to use the multiblock data structure?

With best,

Song]]>

To the first point: There is way to create a data processors that handles boundaries separately. You have to use a boundedBoxProcessingFunctional for that. Therein, depending on the dimension of the problem, you have specify separate methods for every corner/edge and the bulk of the domain.

The best way to find out how this works is to look at the computeGradient() Function, which is located in

src/dataprocessors/dataAnalysisWrapper2D

with the corresponding Functional in

src/dataprocessors/dataAnalysisFunctional2D

To the second point: In the lattice Descriptors (D2Q9, D3Q19), there is a parameter called 'vicinity' which specifies 'how local' you have to be (for D2Q9 and D3Q19 it is set to 1). If you want to access two cells further away, you have change it. I have had a similar problem so you might also check out

[www.palabos.org]

for that.

Hope this helps.

Regards,

Knut]]>

I am interested in implementing some simple spatial filters and trying to write a functional for it. I am now facing some problems and was wondering if anyone has any idea about it:

- As descript in the documentation:

QuoteRule 0 of data processors: A data processor must always be written in such a way that executing the data processor on a given domain has the same effect as splitting the domain into two sub-domains, and then executing the data processor consecutively on each of these sub-domains.

Does it mean, there is no way to treat the boundary separately?

- Also, as written in the documentation:

QuoteOn nearest-neighbor lattices (D2Q9, D3Q19, etc.), you can be non-local by one cell but no more

If so, what can I do if I want to access two or even more cells further neighbors? Does it mean this is no way to write a functional for such operation?

I will very appreciate your answers,

With best,

Song]]>

I am trying to modify the aneurysm show case to run a sloshing simulation in a tank (I have the STL file). I started off by removing the inlet outlet parts of the code, since the tank won't have any inlet and outlet in my case. I was also able to successfully load the STL file. I already have a function which calculates the force at a time instant and have changed the BGKdynamics to GuoExternalForceSmagorinskyBGKdynamics. I am stuck in how to initialize the fluid inside the tank. At start time, there should be 20% water in the tank. Now I could do this like it has been done in damBreak case but in this case I also have baffles in the tank and owing to the large number of baffles, it's not possible for me to manually write a condition for X,Y and Z coordinates to initialize fluid. Is it possible to have a general piece of code which can initialize the fluid inside such a domain without manually specifying the locations of the obstacles (as has been done in the damBreak case)?

Thank you in advance

Vishwesh]]>

This system of equations contain two nonlinear equation with two main unknowns. In the one of these two equatios, we have a third order differential term that is first order temporal and second order spatial.

Now I do not know that how can I expand the LB equation to cover this term.

I will be so grateful if you help me.

Regards]]>

I'm having trouble imposing boundary conditions on a voxelized domain similar to the aneurysm example. In that example a velocity inlet boundary condition is used, while in the other two outlets, a constant pressure Neumann boundary condition is imposed.

I am trying to make something similar, but with a constant pressure Neumann boundary condition on the inlet and a time dependant Neumann boundary condition on the outlet, so that the flow oscillates.

Also, I'm having trouble on the geometry importing step. I can't seem to find a way to make the domain fit to the size of my drawing. Also like in the aneurysm case, where the domain envelops the 3D mesh almost perfectly. When I thy to follow a similar procedure, my domain doesn't match the size of my drawing.

Thank you all.]]>

I am also interested in the implementation of dynamic LES. The dynamicSmagorinsky dynamics or functionals in the "complexDynamics" folder seems to be a framework only instead of an implementation.

I am trying to follow the paper by Premnath et al. (2009) “Dynamic subgridscale modeling of turbulent flows using lattice-Boltzmann method”, would anyone provide any advice on the implementation in Palabos? I will very appreciate it.

With best,

Song]]>

I've tried to use the following line of code but it didn't compile.

lattice.get(someX,someY,someZ).computeDensity(someVariable)

What I'm trying to do with the above code is accessing the cell and sum up the distribution functions of the cell and spit out the density value to "someVariable".

Could anyone tell me what is wrong with this approach?

In addition, I found out that I can get the density of cells over the entire domain by simply doing this

*computeDensity(lattice)

which can be useful later for me with one condition.

I need to store the computed values in a vector to be used for a third-party library but I have no idea how to convert MultiScalarFieldXD to a vector.

Any help would be greatly appreciated

Thanks in advance,

hjung]]>

How do you find out the G and force value? Is there any formula or relationship? Also, how do we find delta t in the current code?

I am trying to write a piece of code which can apply a periodic external force on the liquids. Now, I have written the function to do it but the value is dependent on time but in the code given in the showCases (MultiComponent2d), there is no mention of delta t. I can understand it is because of some lattice unit conversion but can someone explain it to me?

Vishwesh]]>

I am trying to validate the multicomponent2D- rayleigh taylor code given in palabos show cases using this paper. There are several queries that I have regarding this:

1. The paper mentions:

Quote

The initial perturbation of the interface is a sinusoidal function, with amplitude A = 0.1W and wavelength λ = 2π/W

When I refer to the current code, the interface is defined using the following piece of code:

Language: C++// Initialize top layer. applyIndexed(heavyFluid, Box2D(0, nx-1, 0, ny-1), new TwoLayerInitializer<T,DESCRIPTOR>(ny, true) ); // Initialize bottom layer. applyIndexed(lightFluid, Box2D(0, nx-1, 0, ny-1), new TwoLayerInitializer<T,DESCRIPTOR>(ny, false) );

Considering that I do know the equation for the interface which will be a function of x, how can I initialize the interface according to this?

2. Also, the paper does not mention about perturbation due to density variation so I changed the densityFluctuation to 0. Will it cause any problems?

3. The computational domain in the paper has been described to be W = 100 and height = 4W = 400. Also, the grid has also been defined. But, in the current example, the domain size is defined by nx and ny and changing them is the only way to change the domain as well as the mesh. Any other way to resolve it?

4. I referred to a paper on Shan-Chen model which described G to be a quantity which defines the interaction between the fluids. Also, according to the paper, G value cannot be calculated using a formula. If this is the case then how did we arrive to G=0.15

5. Why is gravity defined to be 0.45/ny? Has this been done to make the quantity dimensionless? If that's the case, how did we arrive to 0.45?

I would be really grateful for any help I can get regarding this matter.

Thanks in advance

Vishwesh]]>

According with the xml file you provided, I think that you are based on the aneurysm example. Are you using a provided example (like aneurysm.cpp) and just changing the stl file, or are you using your own palabos implementation?

In my case, the problem was that I wasn't setting the inlets and outlets properly. I figured it out by plotting the profiles at these boundaries in the first iteration. I solved it by following the aneurysm example. Let me know if that's your case so I can try to help you.

Regards,

Emilio]]>

That's really helpful! Thank you very much!

Cheers,

Eric]]>

[github.com]]]>

I was trying to understand how external geometry can be used as an input on Palabos and I was using the externalFlowAroundObstacle as a starting point for this. I wanted to simulate a circular duct, similar to a poiseuille flow example but with an imported geometry. So drew the duct on a CAD software and gave it a length slightly lesser than the simulation domain, if I try to have it with the same size or bigger than the domain, I get the following error messages:

Errors

------

Caught exception in process 0. Palabos generic exception:

List of error messages (not chronological):

1) Error treating the geometry in the Guo off-lattice model.

2) Guo off-lattice model could not find an intersection with a triangle.

But once I had the duct half outside the domain by accident and the simulation ran fine, but I can't understand how.

I wanted to do this in order to have the boundary condition completely inside the duct.

Best regards,

Fernando]]>