[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

[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

[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

[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

[chris]

Thank you, Mathias!

I have found the functions “incBgkCollision” and “constRhoBgkCollision”. I suppose that the former is related to the incEquilibrium [Zou et al., 1995]. Is this correct?

Then, what is the related work to the latter?

Best,
Chris

[chris]

Your equation f_OpenLB = f – f^eq means Cell[iPop] = f_iPop – t_iPop? or it is something else?

In openLB, I have found two expressions to calculate pressure.
1) p = rho / invCs2 which can be reduced to p = rho * Cs2
2) p = (rho-1) * Cs2.
Do you mean 1) is the usual pressure and 2) is the special pressure only for the codes of openLB?

Best
Chris

[chris]

Thank you again for your help!

In the line 353 of dynamics.hh, the pressure is calculated as p = rho / invCs2, which is calculated as I expect.
However, in the line 450 of util.h, the lattice pressure is calculated by p = (rho-1) * Cs2.
The difference between rho and rho-1 comes from the difference between pressure and lattice pressure?
Then, is there any reason that the lattice pressure is defined by using the factor of rho-1 instead of rho?

Best,
Chris

[Author]

Posts