DESCRIPTOR D2Q9Descriptor

[chris]

Hello,

I would like to know whether the example “bstep2d.cpp” uses the functions in lbHelpers.h or lbHelpersD2Q9.h. If the functions in lbHelpers.h are used, then how can I use the functions in lbHelpersD2Q9?

And “bstep2d.cpp” uses “zouHe boundary conditions”? If not, how do I know the information of the boundary condition that is used for the example?

For using a new Descriptor, which files should I prepare? Up to now, I have found that latticeDescriptors and lbHelpers(lbHelpersD2Q9) are related. Is there any other files to prepare?

Best
Chris

[Markus Mohrhard]

Hello Chris,

bstep2d uses a BGKDynamics with a D2Q9Descriptor which in the end calls into lbHelpersD2Q9.h.

prepareLattice in bstep2d.cpp tells you which boundaries are used in the example. The variable bulkDynamics corresponds to the BGK fluid dynamics, bounce back is obvious and only the velocity and the pressure boundary require some more work.

By default the example uses local boundaries (273 createLocalBoundaryCondition2D<T,DESCRIPTOR>( sBoundaryCondition ); ) which means that we get a RegularizedBoundaryManager2D and as a Result CombinedRLBdynamics (with MixinDynamics as BGKdynamics). The only difference between the pressure and the velocity boundary in that case is the used boundary momenta.

Which descriptor are you currently missing. Depending on the type of descriptor you might just need to add a new descriptor. There are only a few cases where putting work into the lbHelpers and the corresponding template specializations is ever required.

Regards,
Markus

[chris]

Hello, Markus,

Thank you so much for your explanation.

I can see that bstep2d.cpp uses the class BGKdynamics and an instance of the class is bulkDynamics. And the class BGKdynamics is declared in dynamics.h. But I do not see that bulkDynamics uses functions in lbHelpersD2Q9.h but not in lbHelpers.h.

I would like to use D2Q9 lattice with a different f_i^eq. Then, it seems lbHelpersD2Q9.h and latticeDescriptors.hh should be modified. Do I need to modify any other files, too?

Best
Chris

[Markus Mohrhard]

Hey Chris,

lbHelpers.h and lbHelpersD2Q9.h belong together. lbHelpersD2Q9.h contains the template specialization for the lbDynamicsHelpers class for the D2Q9DescriptorBase (and therefore for the D2Q9Descriptor). The call to lbHelpers<D2Q9Descriptor>::some_method() is forwarded to lbDynamicsHelper<D2Q9Descriptor::BaseDescriptor>::some_method. The compiler now has the generic version for that method in lbHelpers.h but an explicit template specialization in lbHelpersD2Q9.h that fits as well. As the C++ standard requires the compiler to select the explicit template specialization the code will use the lbHelpersD2Q9.h version in the end.

My suggestion for using a different equilibrium implementation would be to define a new descriptor (keep in mind that you also need to change the base descriptor) and then implement an explicit template specialization for this new descriptor (new descriptor in latticeDescriptors.h[h] and new template specialization in a new lbHelpersXY.h file (keep in mind that the new file needs to be included in lbHelpers.h)). That way you can keep your changes separate and can easily switch back to the classical equilibrium calculation.

Regards,
Markus

[chris]

Hi, Markus,

Thank you again for your helpful explanation.

I have tried to modify D2Q9 and obtained a result which is not perfect.

The followings are what I modified:

Adding codes for D2Q9M in latticeDescriptors.h

Code:

/// D2Q9M lattice
template <typename T> struct D2Q9MDescriptorBase {
typedef D2Q9MDescriptorBase<T> BaseDescriptor;
enum { d = 2, q = 9 }; ///< number of dimensions/distr. functions
static const int vicinity; ///< size of neighborhood
static const int c[q][d]; ///< lattice directions
static const int opposite[q]; ///< opposite entry
static const T t[q]; ///< lattice weights
static const T invCs2; ///< inverse square of speed of sound
};

template <typename T> struct D2Q9MDescriptor
: public D2Q9MDescriptorBase<T>, public NoExternalFieldBase {
};

template <typename T> struct ForcedD2Q9MDescriptor
: public D2Q9MDescriptorBase<T>, public Force2dDescriptorBase {
};

Adding codes for D2Q9M in latticeDescriptors.hh

Code:

// D2Q9M /////////////////////////////////////////////////////////

template<typename T>
const int D2Q9MDescriptorBase<T>::vicinity = 1;

template<typename T>
const int D2Q9MDescriptorBase<T>::c
[D2Q9MDescriptorBase<T>::q][D2Q9MDescriptorBase<T>::d] = {
{ 0, 0},
{-1, 1}, {-1, 0}, {-1,-1}, { 0,-1},
{ 1,-1}, { 1, 0}, { 1, 1}, { 0, 1}
};

template<typename T>
const int D2Q9MDescriptorBase<T>::opposite[D2Q9MDescriptorBase<T>::q] = {
0, 5, 6, 7, 8, 1, 2, 3, 4
};

template<typename T>
const T D2Q9MDescriptorBase<T>::t[D2Q9MDescriptorBase<T>::q] = {
(T)4/(T)9, (T)1/(T)36, (T)1/(T)9, (T)1/(T)36, (T)1/(T)9,
(T)1/(T)36, (T)1/(T)9, (T)1/(T)36, (T)1/(T)9
};

template<typename T>
const T D2Q9MDescriptorBase<T>::invCs2 = (T)3;

Adding codes for D2Q9M in latticeDescriptors.cpp

Code:

template struct D2Q9MDescriptorBase<double>;
template struct D2Q9MDescriptor<double>;

Adding a file of lbHelpersD2Q9M.h
From lbHelpersD2Q9.h, we have obtained lbHelpersD2Q9M.h by changing
descriptors::D2Q9DescriptorsBase —> descriptors::D2Q9MDescriptorsBase
descriptors::ForcedD2Q9Descriptor —> descriptors::ForcedD2Q9MDescriptor
and modifying the followings:

Code:

static T equilibrium( int iPop, T rho, const T u[2], T uSqr )
{
typedef descriptors::D2Q9MDescriptorBase<T> L;
T c_u = L::c[iPop][0]*u[0] + L::c[iPop][1]*u[1];
T ccuu = L::c[iPop][0]*u[0]*L::c[iPop][0]*u[0]
+ L::c[iPop][1]*u[1]*L::c[iPop][1]*u[1];
return rho * L::t[iPop] * (
1. + 3.*c_u + 4.5*c_u*c_u – 1.5*uSqr
+ 81/4.*ccuu )
– L::t[iPop];
}

static T incEquilibrium( int iPop, const T j[2], const T jSqr, const T pressure )
{
typedef descriptors::D2Q9MDescriptorBase<T> L;
T c_j = L::c[iPop][0]*j[0] + L::c[iPop][1]*j[1];
T ccjj = L::c[iPop][0]*j[0]*L::c[iPop][0]*j[0]
+ L::c[iPop][1]*j[1]*L::c[iPop][1]*j[1];
return L::t[iPop] * (
3.*pressure + 3.*c_j + 4.5*c_j*c_j – 1.5*jSqr
+ 81/4.*ccjj )
– L::t[iPop];
}

Code:

static T bgkCollision ( CellBase<T,descriptors::D2Q9MDescriptorBase<T> >& cell, T const& rho, const T u[2], T const& omega)
{

… … … … …
cell[7] *= (T)1-omega;
cell[7] += omega*(cf_ + rho_*(ux_ + uy_ + uxPySqr_ – uSqr_ + (T)81/(T)4*uxSqr*uySqr));
cell[1] *= (T)1-omega;
cell[1] += omega*(cf_ + rho_*(-ux_ + uy_ + uxMySqr_ – uSqr_+ (T)81/(T)4*uxSqr*uySqr));
cell[3] *= (T)1-omega;
cell[3] += omega*(cf_ + rho_*(-ux_ – uy_ + uxPySqr_ – uSqr_+ (T)81/(T)4*uxSqr*uySqr));
cell[5] *= (T)1-omega;
cell[5] += omega*(cf_ + rho_*(ux_ – uy_ + uxMySqr_ – uSqr_+ (T)81/(T)4*uxSqr*uySqr))
return uxSqr + uySqr; }

Adding codes in lbHelpers.h
#include “lbHelpersD2Q9M.h”

But the result is not perfect. It seems there is a problem on the right end: the abnormal values propagate to the left direction. (On the left end, it seems there is no problem: a flow is generated.)

Below is the console output and we can see that around t=4, the values become nan.

Code:

[LatticeStatistics] step=9600; t=3.2; uMax=0.0243807; avEnergy=7.62342e-06; avRho=1.00168
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.186842; meanVelocity[m/s]=0.127392
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=2.08237; meanPressure[Pa]=1.4198
[Timer] step=9900; percent=8.25; passedTime=15.886; remTime=176.672; MLUPs=55.277
[LatticeStatistics] step=9900; t=3.3; uMax=0.0132551; avEnergy=8.30207e-06; avRho=1.00198
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.130843; meanVelocity[m/s]=0.089211
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=4.22094; meanPressure[Pa]=2.87791
[Timer] step=10200; percent=8.5; passedTime=16.311; remTime=175.583; MLUPs=67.2148
[LatticeStatistics] step=10200; t=3.4; uMax=0.0331086; avEnergy=8.79137e-06; avRho=1.00255
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.115719; meanVelocity[m/s]=0.0788991
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=5.19487; meanPressure[Pa]=3.54196
[Timer] step=10500; percent=8.75; passedTime=16.833; remTime=175.544; MLUPs=54.9587
[LatticeStatistics] step=10500; t=3.5; uMax=0.0332306; avEnergy=9.43505e-06; avRho=1.00304
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.186313; meanVelocity[m/s]=0.127031
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=3.7029; meanPressure[Pa]=2.52471
[Timer] step=10800; percent=9; passedTime=17.261; remTime=174.528; MLUPs=66.9007
[LatticeStatistics] step=10800; t=3.6; uMax=0.0252073; avEnergy=9.85379e-06; avRho=1.00357
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.163435; meanVelocity[m/s]=0.111433
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=4.9358; meanPressure[Pa]=3.36532
[Timer] step=11100; percent=9.25; passedTime=17.781; remTime=174.446; MLUPs=55.1705
[LatticeStatistics] step=11100; t=3.7; uMax=0.0247429; avEnergy=1.02304e-05; avRho=1.00417
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.15393; meanVelocity[m/s]=0.104952
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=5.7863; meanPressure[Pa]=3.9452
[Timer] step=11400; percent=9.5; passedTime=18.19; remTime=173.284; MLUPs=69.8378
[LatticeStatistics] step=11400; t=3.8; uMax=0.0414893; avEnergy=1.03922e-05; avRho=1.00489
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.155016; meanVelocity[m/s]=0.105693
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=6.37754; meanPressure[Pa]=4.34832
[Timer] step=11700; percent=9.75; passedTime=18.703; remTime=173.123; MLUPs=55.9248
[LatticeStatistics] step=11700; t=3.9; uMax=0.0510242; avEnergy=1.11023e-05; avRho=1.0057
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.139137; meanVelocity[m/s]=0.0948658
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=7.46674; meanPressure[Pa]=5.09096
[Timer] step=12000; percent=10; passedTime=19.125; remTime=172.125; MLUPs=67.8519
[LatticeStatistics] step=12000; t=4; uMax=0.0502752; avEnergy=1.26411e-05; avRho=1.00665
[getResults] Checkpointing the system at t=12000
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=0.0973408; meanVelocity[m/s]=0.0663687
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=9.33927; meanPressure[Pa]=6.36769
[Timer] step=12300; percent=10.25; passedTime=19.694; remTime=172.443; MLUPs=50.4111
[LatticeStatistics] step=12300; t=4.1; uMax=0.055468; avEnergy=1.5976e-05; avRho=1.00776
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=-73.3333; meanVelocity[m/s]=-50
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=1.3899e+25; meanPressure[Pa]=9.47662e+24
[Timer] step=12600; percent=10.5; passedTime=20.111; remTime=171.422; MLUPs=68.5012
[LatticeStatistics] step=12600; t=4.2; uMax=1.41421; avEnergy=-nan; avRho=-nan
“./tmp/imageData/data/heatMapl2physVelocity.p”, line 20: warning: matrix contains missing or undefined values
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=-nan; meanVelocity[m/s]=-nan
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=-nan; meanPressure[Pa]=nan
[Timer] step=12900; percent=10.75; passedTime=20.641; remTime=171.368; MLUPs=54.1276
[LatticeStatistics] step=12900; t=4.3; uMax=1.05969; avEnergy=-nan; avRho=-nan
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=-nan; meanVelocity[m/s]=-nan
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=-nan; meanPressure[Pa]=nan
[Timer] step=13200; percent=11; passedTime=21.07; remTime=170.475; MLUPs=66.7448
[LatticeStatistics] step=13200; t=4.4; uMax=1.05974; avEnergy=-nan; avRho=-nan
“./tmp/imageData/data/heatMapl2physVelocity.p”, line 20: warning: matrix contains missing or undefined values
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=-nan; meanVelocity[m/s]=-nan
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=-nan; meanPressure[Pa]=nan
[Timer] step=13500; percent=11.25; passedTime=21.594; remTime=170.353; MLUPs=54.6441
[LatticeStatistics] step=13500; t=4.5; uMax=1.05976; avEnergy=-nan; avRho=-nan
[SuperPlaneIntegralFluxVelocity2D] regionSize[m]=1.46667; flowRate[m^2/s]=-nan; meanVelocity[m/s]=-nan
[SuperPlaneIntegralFluxPressure2D] regionSize[m]=1.46667; force[N]=-nan; meanPressure[Pa]=nan
[Timer] step=13800; percent=11.5; passedTime=22.019; remTime=169.451; MLUPs=67.2148

Which parts should be again modified?

Best
Chris

[Markus Mohrhard]

Hey Chris,

I’m sorry but I can honestly not help you as I have no idea what you are trying to implement. I assume that you have a paper that you are trying to reproduce but you need to carefully check the results of the paper.

The lattice descriptor part looks fine.

While implementing you need to be a bit careful as our f_i and therefore the f_i^eq is shifted by t_i. This is usually the only part that one needs to pay special attention to while implementing but our rho and u calculation should already deal correctly with the shift (see e.g. adding 1 back to rho).

Kind regards,
Markus

[chris]

Hi, Markus,

Thank you again.

I am just modifying the f_i^eq of the standard D2Q9 to f_i^eq + 81/4 * u_x ^2 * u_y ^2 * e_x ^2 * e_y ^2. The extra adding term is one of the higher-order term in u.

I wonder is there anything to modify in lbHelpers or in boundary conditions or elsewhere.

Best regards,
Chris

[Markus Mohrhard]

Hey Chris,

normally nothing should be required. I’d start by writing a small test for one single grid point and make sure that the results are correct. I’m sorry but as I have never looked into higher order solutions I can not provide you with more guidance.

Regards,
Markus

[chris]

Hey Markus and lbTeam,

Thank you a lot.

bstep2d.cpp calls UnitConverter<T, DESCRPITOR>, however, unitConverter.cpp does contain the only D2Q9Descriptor for 2D as the following:

Code:

template class UnitConverter<double,descriptors::D2Q9Descriptor>;
template class UnitConverter<double,descriptors::D3Q19Descriptor>;
template class UnitConverter<double,descriptors::AdvectionDiffusionD3Q7Descriptor>;

template UnitConverter<double,descriptors::D2Q9Descriptor>* createUnitConverter(const olb::XMLreader&);

I would like to know what happens when we use

Code:

#define DESCRIPTOR D2Q9_M_Descriptor

in bstep2d.cpp.

Also in firstOrderLbHelpers2d.h, there exists only “descriptors::D2Q9Descriptor”.

Best,
Chris

[Author]

Posts