Poiseuille2D: Interpolated and local values(velocity and pressure)

[achodankar]

Hello Developers,
I unable to track how the values for the velocity and pressure are set on the boundaries and the whole domain. I recognize the code lines for it.

// Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( superGeometry, 1, rho, u );
sLattice.iniEquilibrium( superGeometry, 1, rho, u );
sLattice.defineRhoU( superGeometry, 2, rho, u );
sLattice.iniEquilibrium( superGeometry, 2, rho, u );
sLattice.defineRhoU( superGeometry, 3, rho, u );
sLattice.iniEquilibrium( superGeometry, 3, rho, u );
sLattice.defineRhoU( superGeometry, 4, rho, u );
sLattice.iniEquilibrium( superGeometry, 4, rho, u );

For interpolated(I can track back to the setInterpolatedPressureBoundary2D.h & setInterpolatedPressureBoundary2D.hh files for pressure, and similarly, setInterpolatedVelocityBoundary2D.h & setInterpolatedVelocityBoundary2D.hh files for velocity)

setInterpolatedVelocityBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 4);

For local(I can track back to the setlocalPressureBoundary2D.h & setlocalPressureBoundary2D.hh files for pressure, and similarly, setlocalVelocityBoundary2D.h & setlocalVelocityBoundary2D.hh files for velocity)
setLocalVelocityBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 3);
setLocalPressureBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 4);

For pressure, I do see this code snippet, but I don’t understand how it is tied back to above code to set the pressure values

From util.h

483 template <typename T, typename DESCRIPTOR>
484 T densityFromPressure( T latticePressure )
485 {
486 // rho = p / c_s^2 + 1
487 return latticePressure * descriptors::invCs2<T,DESCRIPTOR>() + 1.0;
488 }

I do come across code to initialise the velocity, momenta, and other parameters in momentonboundaries.hh file, but no mention of pressure values. Is stress function used to set pressure values? Here, I fail to understand how it gets tied back to the interpolated and local pressure and velocity values file.

I would really appreciate your feedback on this matter. A code mapping for it would be of great help to understand the code structure. The map next to the code is not detailed to understand it.

Thank you.

Yours sincerely,

Abhijeet C.

[mathias]

You may navigate within the code easily using doxygen. In LBM, you set the pressure by setting rho — the transformation formular is the equation of state.

[achodankar]

Hello Mathias,
I see in the util.h file mentioning the equation of state, which can be employed to set the pressure values. However, in the code, I don’t know understand how the pressure values are set using the equation of set. There are functions to initialise the density, velocity, and momenta, but nothing for pressure. I have been using Doxygen for the same. The two searches closest to pressure are: PressureBM and pressureFromDensity. I would be grateful if you can direct me how the pressureFromDensity function is employed while setting the local and interpolated pressure values. In the pressureBM file, they compute density, velocity, momenta, and stress, but no mention of pressure. I have seen the same pattern while browsing through other code files. I would appreciate any help on this matter.

Thank you.

Yours sincerely,

Abhijeet C.

[mathias]

You can use the unitConverter to get rho for pressure and initialize this rho. Best Mathias

[achodankar]

Hello Mathias,
I did follow the procedure as you suggested. The density value is already set in the unit converter.

SuperLatticeDensity2D<T,DESCRIPTOR> densityF(sLattice);
AnalyticalFfromSuperF2D<T> intpolatedensity_lat(densityF, true );

intpolatedensity_lat( density_lat,point );

//compute lattice pressure
pressure_lat[0]=(density_lat[0])/descriptors::invCs2<T,DESCRIPTOR>();

//Compute physical pressure
pressure_phys[0] = (pressure_lat[0] * converter.getConversionFactorPressure())+converter.getCharPhysPressure();

If I set characteristic pressure as 101325 Pa, the lattice pressure is 0.3333 and the physical pressure is 101349 Pa.
However, when I extract the pressure values from the usual method, the value of physical pressure is 101325 Pa, which is not correct. I tried the code to extract lattice pressure from one of the previous posts in the forum but was unable to extract. What would the syntax be? The commented lines compile without any errors. I am unable to extract those values. Also, why do we subtract 1. Based on equation of state: P = rho * cs^2 (in the blocklatticepressurephys code file also they have subtracted one).

Here is the code snippet:

SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure_phys(sLattice, converter);
AnalyticalFfromSuperF2D<T> intpolatepressure_phys( pressure_phys, true );
intpolatepressure_phys( numericalpressure_phys,point );

//Store the lattice Pressure
// SuperLatticeDensity2D<T,DESCRIPTOR> densityF(sLattice);
// auto latticePressureF = (functor_dsl::lift(densityF) – 1.0) / descriptors::invCs2<T,DESCRIPTOR>();

I would be grateful for your suggestions on this issue

Thank you.

Yours sincerely,

Abhijeet C.

[achodankar]

Hello Mathias,
Once that 1 is removed from the formula inputted in the code file, the pressure values do match properly. I would be grateful if you can point if I am doing anything wrong. Thank you once again for your responses.

Thank you.

Yours sincerely,

Abhijeet C.

[jjessberger]

Hello Abhijeet,

I think you are doing everything right. The constant offset of 1.0 arises because the populations are stored with such an offset (for the reason of numerical accuracy). On the other hand, you may use functors like SuperLatticePhysPressure2D which apply this offset automatically.

Yours,
Julius

[achodankar]

Hello Julius,
Actually, I removed the offset one from the SuperLatticePhysPressure2D code file. That’s how my results are matching now. Before that, the results were differing and incorrect values. Will removing the offset make a difference ? I don’t understand how it is related to numerical accuracy. Thank you for your time once again.

Thank you.

Yours sincerely,

Abhijeet C.

[jjessberger]

Hello Abhijeet,

The density offset is introduced in order to increase machine precision.

If you need to remove the offset in the pressure computation, it has probably been introduced in a wrong way. How did you define your characteristic pressure?

Yours,
Julius J.

[achodankar]

Hello Julius,
I defined the characteristic pressure as 101325 Pa using the unit converter.

Here is the code snippet:

UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
int {N}, // resolution: number of voxels per charPhysL
(T) lat_relax_time,//0.8 // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) char_phys_length, // charPhysLength: reference length of simulation geometry
(T) char_phys_velocity, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) phys_kinematic_viscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) phys_density, // physDensity: physical density in __kg / m^3__
(T) 101325. //charphysPressure
);
(latticePhysPressure2D.hh) Here is the code snippet where I removed the offset for pressure:

template <typename T, typename DESCRIPTOR>
bool BlockLatticePhysPressure2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
{
// lattice pressure = c_s^2 ( rho -1 )
T latticePressure = ( this->_blockLattice.get( input[0], input[1] ).computeRho()) / descriptors::invCs2<T,DESCRIPTOR>();
output[0] = this->_converter.getPhysPressure(latticePressure);

return true;
}

Thank you.

Yours sincerely,

Abhijeet C.

[jjessberger]

Hello Abhijeet,

Unfortunately, I cannot reproduce this problem: I entered the charPhysPressure as you did, did not modify any offset and got the correct result (tested on olb release 1.4, poiseuille2d, both forced and nonForced mode).

Yours,
Julius

[achodankar]

Hello Julius,
Thank you for your prompt response again. The source code has the offset. The formula for lattice pressure in the code is given as: p= cs2 * (rho-1). You had mentioned earlier the offset is for the reason of numerical accuracy. But, I am unable to comprehend how the pressure values will match.

template <typename T, typename DESCRIPTOR>
bool BlockLatticePhysPressure2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
{
// lattice pressure = c_s^2 ( rho -1 )
T latticePressure = ( this->_blockLattice.get( input[0], input[1] ).computeRho()-1) / descriptors::invCs2<T,DESCRIPTOR>();
output[0] = this->_converter.getPhysPressure(latticePressure);

return true;
}

There is no offset in the density in the source code:

template <typename T, typename DESCRIPTOR>
bool BlockLatticeDensity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
{
output[0] = this->_blockLattice.get( input[0], input[1] ).computeRho();
return true;
}

Thank you once again for your patience and help. I really appreciate it.

Thank you.

Yours sincerely,

Abhijeet C.

[jjessberger]

I’m sorry I think I mixed things up. The offset that we are talking about is not because of the storage of the lattice populations (this happens deeper in the code). Instead, it simply seems to be since one wants to correspond density=1 to pressure=0 as a default. Consequently, one has to modify the offset when a different char. phys. pressure is defined. The UnitConverter does this.

[achodankar]

Hello Julius,
I understand it well now.

Thank you.

Yours sincerely,

Abhijeet C.

[Author]

Posts