Microflow 3D  v1.0
Collision.cpp
Go to the documentation of this file.
1 // ==============================================================================================
2 // Microflow 3D, http://www.microflow.pwr.edu.pl/
3 // Created by Roman Szafran on 29.05.19.
4 // Modified by Roman Szafran, last on 07.09.19.
5 // Copyright (c) 2019 Wroclaw University of Science and Technology.
6 // Distributed under the Apache License, Version 2.0. You may obtain a copy of the License at
7 // http://www.apache.org/licenses/LICENSE-2.0 or see accompanying file license.txt.
8 // Redistributions of source code must retain the above copyright and license notice.
9 // ==============================================================================================
10 
11 #include "Collision.h"
12 
13 #define f(x) (pNode->FQ19[x])
14 #define U (Vector4.x)
15 #define V (Vector4.y)
16 #define W (Vector4.z)
17 #define RHO (Vector4.rho)
18 #define DirectionVectorComponent_Cx(x) (m_LatticeParameters_Ptr->DirectionVectorComponent_Cx[x])
19 #define DirectionVectorComponent_Cy(x) (m_LatticeParameters_Ptr->DirectionVectorComponent_Cy[x])
20 #define DirectionVectorComponent_Cz(x) (m_LatticeParameters_Ptr->DirectionVectorComponent_Cz[x])
21 #define DirectionWeight_W(x) (m_LatticeParameters_Ptr->DirectionWeight_W[x])
22 
23 void MF::Solver_CPU::Collision::SetCollisionPointerToFunc(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
24  for (auto & Thread_Ptr : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
25  switch (Thread_Ptr->m_NodeType) {
26  case 0: {break;} // for solid node
27  case 61: { // for BB2 61 node
28  Thread_Ptr->m_pColideFunc = MF::Solver_CPU::Collision::Coll_BounceBackForNode61;
29  Thread_Ptr.get()->m_DoCollision = false;
30  break;
31  }
32 
33  default:
34  case 1: { //All other nodes except BB2 61
35  if (m_CaseParameters_Ptr->CollisionModel_KT == MF::Solver_CPU::ModelTypeClass::MRT) { //MRT collision model
36  if (m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {
37  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_MRT_Incompr; //Incompressible MRT collision model
38  Thread_Ptr.get()->m_DoCollision = true;
39  }
40  else {//QuasiCompressible MRT
41  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_MRT_Qcompr;
42  Thread_Ptr.get()->m_DoCollision = true;
43  }
44  }
45  else if (m_CaseParameters_Ptr->CollisionModel_KT == MF::Solver_CPU::ModelTypeClass::MRT2) { //MRT collision model
46  if (m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {
47  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_MRT2_Incompr; //Incompressible MRT collision model
48  Thread_Ptr.get()->m_DoCollision = true;
49  }
50  else {//QuasiCompressible MRT
51  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_MRT2_Qcompr;
52  Thread_Ptr.get()->m_DoCollision = true;
53  }
54  }
55  else if (m_CaseParameters_Ptr->CollisionModel_KT == MF::Solver_CPU::ModelTypeClass::FBGK) //FBGK
56  {
57  if (m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {
58  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_FBGK_Incompr; //Incompressible FBGK
59  Thread_Ptr.get()->m_DoCollision = true;
60  }
61  else {//QuasiCompressible FBGK
62  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_FBGK_Qcompr;
63  Thread_Ptr.get()->m_DoCollision = true;
64  }
65  }
66  else //BGK
67  {
68  if (m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {//Incompressible BGK
69  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_BGK_Incompr;
70  Thread_Ptr.get()->m_DoCollision = true;
71  }
72  else {//QuasiCompressible BGK
73  Thread_Ptr.get()->m_pColideFunc = MF::Solver_CPU::Collision::Coll_BGK_Qcompr;
74  Thread_Ptr.get()->m_DoCollision = true;
75  }
76  }
77  break;
78  }
79  }
80  }
81 }
82 
83 //-----------------------------------------------------------------------------------------------------------------------------
84 // BGK ------------------------------------------------------------------------------------------------------------------------
85 //-----------------------------------------------------------------------------------------------------------------------------
86 
87 inline double MF::Solver_CPU::Collision::Feq_Qcompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
88  const double csq = m_LatticeParameters_Ptr->LatticeConstant_CSQ;
89  const double cx=DirectionVectorComponent_Cx(k);
90  const double cy=DirectionVectorComponent_Cy(k);
91  const double cz=DirectionVectorComponent_Cz(k);
92  const double wg=DirectionWeight_W(k);
93  const double cu=cx*U+cy*V+cz*W; // ck * u
94  const double U2=U*U+V*V+W*W; //u*u
95  return (wg*RHO*(1.0 + 1.0 / csq * cu + 1.0 / (2.0 * csq * csq) * cu * cu - 1.0 / (2.0 * csq) * U2));
96 }
97 
98 inline double MF::Solver_CPU::Collision::Feq_Incompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
99  const double csq = m_LatticeParameters_Ptr->LatticeConstant_CSQ;
100  const double cx=DirectionVectorComponent_Cx(k);
101  const double cy=DirectionVectorComponent_Cy(k);
102  const double cz=DirectionVectorComponent_Cz(k);
103  const double wg=DirectionWeight_W(k);
104  const double cu=cx*U+cy*V+cz*W; // ck * u
105  const double U2=U*U+V*V+W*W; //u*u
106  return (wg*(RHO + 1.0 / csq * cu + 1.0 / (2.0 * csq * csq) * cu * cu - 1.0 / (2.0 * csq) * U2));
107 }
108 
109 void MF::Solver_CPU::Collision::Coll_BGK_Qcompr(MF::Database::Node * pNode) {
110  double FEQ;
111  const double TAU = m_CaseParameters_Ptr->Tau;
112  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
113  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
114  for(unsigned char k=0; k < MFQ19; k++) {
115  FEQ = Feq_Qcompr(Vector4, k);
116  f(k) = f(k)-(f(k) - FEQ) / TAU;
117  }
118  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
119  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
120  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
121  }
122 }
123 
124 void MF::Solver_CPU::Collision::Coll_BGK_Incompr(MF::Database::Node * pNode) {
125  double FEQ;
126  const double TAU = m_CaseParameters_Ptr->Tau;
127  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
128  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
129  for(unsigned char k=0; k < MFQ19; k++) {
130  FEQ = Feq_Incompr(Vector4, k);
131  f(k) = f(k)- (f(k) - FEQ) / TAU;
132  }
133  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
134  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
135  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
136  }
137 }
138 
139 //-----------------------------------------------------------------------------------------------------------------------------
140 // BGK with external force F --------------------------------------------------------------------------------------------------
141 //-----------------------------------------------------------------------------------------------------------------------------
142 
143 inline double MF::Solver_CPU::Collision::Feq_Qcompr_F(const MF::Database::Vec4<double> &Vector4, const unsigned char &k, const MF::Database::Vec3<double> & XYZForce_LB) {
144  const double its_csq = m_LatticeParameters_Ptr->LatticeConstant_CSQ;
145  const double TAU = m_CaseParameters_Ptr->Tau;
146  const double ax = XYZForce_LB.x / RHO;
147  const double ay = XYZForce_LB.y / RHO;
148  const double az = XYZForce_LB.z / RHO;
149 
150  const double cx=DirectionVectorComponent_Cx(k);
151  const double cy=DirectionVectorComponent_Cy(k);
152  const double cz=DirectionVectorComponent_Cz(k);
153  const double wg=DirectionWeight_W(k);
154  const double cu=cx*U+cy*V+cz*W; // ck * u
155  const double U2=U*U+V*V+W*W; //u*u
156  const double cvx = cx-U;
157  const double cvy = cy-V;
158  const double cvz = cz-W;
159  const double cva = cvx*ax+cvy*ay+cvz*az;
160  const double part = 1.0/(2.0*its_csq)*cva;
161  const double feq = (wg*RHO*(1.0+1.0/its_csq*cu+1.0/(2.0*its_csq*its_csq)*cu*cu-1.0/(2.0*its_csq)*U2)) * (1.0-part*(1.0-2.0*TAU));
162  return (feq);
163 }
164 
165 inline double MF::Solver_CPU::Collision::Feq_Incompr_F(const MF::Database::Vec4<double> &Vector4, const unsigned char &k, const MF::Database::Vec3<double> &XYZForce_LB) {
166  const double its_csq = m_LatticeParameters_Ptr->LatticeConstant_CSQ;
167  const double TAU = m_CaseParameters_Ptr->Tau;
168  const double ax = XYZForce_LB.x;
169  const double ay = XYZForce_LB.y;
170  const double az = XYZForce_LB.z;
171 
172  const double cx=DirectionVectorComponent_Cx(k);
173  const double cy=DirectionVectorComponent_Cy(k);
174  const double cz=DirectionVectorComponent_Cz(k);
175  const double wg=DirectionWeight_W(k);
176  const double cu=cx*U+cy*V+cz*W; // ck * u
177  const double U2=U*U+V*V+W*W; //u*u
178  const double cvx = cx-U;
179  const double cvy = cy-V;
180  const double cvz = cz-W;
181  const double cva = cvx*ax+cvy*ay+cvz*az;
182  const double part = 1.0/(2.0*its_csq)*cva;
183  const double feq = (wg*(RHO+1.0/its_csq*cu+1.0/(2.0*its_csq*its_csq)*cu*cu-1.0/(2.0*its_csq)*U2)) * (1.0-part*(1.0-2.0*TAU));
184  return (feq);
185 }
186 
187 void MF::Solver_CPU::Collision::Coll_FBGK_Qcompr(MF::Database::Node * pNode) {
188  double FEQ;
189  const double TAU = m_CaseParameters_Ptr->Tau;
190  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
191  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
192  for(unsigned char k=0; k < MFQ19; k++) {
193  FEQ=Feq_Qcompr_F(Vector4, k, Thread_Ptr->m_XYZForce_LB);
194  f(k)= f(k)-(f(k)-FEQ)/TAU;
195  }
196  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
197  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
198  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
199  }
200 }
201 
202 void MF::Solver_CPU::Collision::Coll_FBGK_Incompr(MF::Database::Node * pNode) {
203  double FEQ;
204  const double TAU = m_CaseParameters_Ptr->Tau;
205  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
206  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
207  for(unsigned char k=0; k < MFQ19; k++) {
208  FEQ=Feq_Incompr_F(Vector4, k, Thread_Ptr->m_XYZForce_LB);
209  f(k)= f(k)- (f(k) - FEQ)/TAU;
210  }
211  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
212  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
213  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
214  }
215 }
216 
217 //-----------------------------------------------------------------------------------------------------------------------------
218 // MRT ------------------------------------------------------------------------------------------------------------------------
219 //-----------------------------------------------------------------------------------------------------------------------------
220 
221 void MF::Solver_CPU::Collision::Coll_MRT_Qcompr(MF::Database::Node * pNode) {
222  double MEQ;
223  double A[MFQ19];
224  const double TAU = m_CaseParameters_Ptr->Tau;
225  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
226  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
227 
228 //Frequencies of relaxation for the MRT model, according to D.D'Humieres et al 2002
229  const double s[MFQ19] = {
230  0.0, //rho (conserved)
231  1.19, //results from the formula for the bulk viscosity = 1.19
232  1.4,
233  0.0, // rho*ux (conserved)
234  1.2,
235  0.0, // rho*uy (conserved)
236  1.2, // =s[4]
237  0.0, // rho*uz (conserved)
238  1.2, // =s[4]
239  1.0/TAU, //results from the formula for kinematic viscosity
240  1.4,
241  1.0/TAU, //the same as s[9]
242  1.4,
243  1.0/TAU, //the same as s[9]
244  1.0/TAU, //the same as s[13]
245  1.0/TAU, //the same as s[13]
246  1.2, //1.98;
247  1.2, //1.98;
248  1.2 //1.98;
249  };
250 
251 //Matrix of coefficients for converting velocities to the momentum space and vice versa
252  const double M[MFQ19][MFQ19] = {
253  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
254  {-30, -11, -11, -11, -11, -11, -11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8},
255  {12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
256  {0, 1, 0, -1, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
257  {0, -4, 0, 4, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
258  {0, 0, 1, 0, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
259  {0, 0, -4, 0, 4, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
260  {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
261  {0, 0, 0, 0, 0, -4, 4, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
262  {0, 2, -1, 2, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
263  {0, -4, 2, -4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
264  {0, 0, 1, 0, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
265  {0, 0, -2, 0, -2, 2, 2, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
266  {0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0},
267  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1},
268  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0},
269  {0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, -1, 1, 1, 0, 0, 0, 0},
270  {0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 1, 1, -1, -1},
271  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, 1, 1, -1}
272  };
273 
274  const double d[MFQ19]={19.0, 2394.0, 252.0, 10.0, 40.0, 10.0, 40.0, 10.0, 40.0, 36.0, 72.0, 12.0, 24.0, 4.0, 4.0, 4.0, 8.0, 8.0, 8.0}; //D=M*MT
275 
276 //Transformation from the velocity to momentum space
277  for(unsigned char k=0; k < MFQ19; k++) {
278  A[k]=0.0;
279  for(unsigned char n=0; n < MFQ19; n++)
280  A[k]+=f(n)*M[k][n]; //M
281  }
282 //Relaxation in the momentum space
283  for(unsigned char k=0; k < MFQ19; k++) {
284  if(s[k]!=0.0) {
285  MEQ=Meq_Qcompr(Vector4, k);
286  A[k]=A[k]-s[k]*(A[k]-MEQ);
287  }
288  A[k]/=d[k];
289  }
290 //Transformation of results back to the velocity space
291  for(unsigned char k=0; k < MFQ19; k++) {
292  f(k)=0.0;
293  for(unsigned char n=0; n < MFQ19; n++)
294  f(k)+=A[n]*M[n][k]; //MT (M transposed)
295  }
296  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
297  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
298  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
299  }
300 }
301 
302 void MF::Solver_CPU::Collision::Coll_MRT_Incompr(MF::Database::Node * pNode) {
303  double MEQ;
304  double A[MFQ19];
305  const double TAU = m_CaseParameters_Ptr->Tau;
306  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
307  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
308 
309 //Frequencies of relaxation for the MRT model, according to D.D'Humieres et al 2002
310  const double s[MFQ19] = {
311  0.0, //rho (conserved)
312  1.19, //results from the formula for the bulk viscosity = 1.19
313  1.4,
314  0.0, // rho*ux (conserved)
315  1.2,
316  0.0, // rho*uy (conserved)
317  1.2, // =s[4]
318  0.0, // rho*uz (conserved)
319  1.2, // =s[4]
320  1.0/TAU, //results from the formula for kinematic viscosity
321  1.4,
322  1.0/TAU, //=s[9]
323  1.4,
324  1.0/TAU, //the same as s[9]
325  1.0/TAU, //the same as s[13]
326  1.0/TAU, //the same as s[13]
327  1.2, //1.98;
328  1.2, //1.98;
329  1.2 //1.98;;
330  };
331 
332 //Matrix of coefficients for converting velocities to the momentum space and vice versa
333  const double M[MFQ19][MFQ19] = {
334  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
335  {-30, -11, -11, -11, -11, -11, -11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8},
336  {12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
337  {0, 1, 0, -1, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
338  {0, -4, 0, 4, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
339  {0, 0, 1, 0, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
340  {0, 0, -4, 0, 4, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
341  {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
342  {0, 0, 0, 0, 0, -4, 4, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
343  {0, 2, -1, 2, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
344  {0, -4, 2, -4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
345  {0, 0, 1, 0, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
346  {0, 0, -2, 0, -2, 2, 2, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
347  {0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0},
348  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1},
349  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0},
350  {0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, -1, 1, 1, 0, 0, 0, 0},
351  {0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 1, 1, -1, -1},
352  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, 1, 1, -1}
353  };
354 
355  const double d[MFQ19]={19.0, 2394.0, 252.0, 10.0, 40.0, 10.0, 40.0, 10.0, 40.0, 36.0, 72.0, 12.0, 24.0, 4.0, 4.0, 4.0, 8.0, 8.0, 8.0}; //D=M*MT
356 
357 //Transformation from the velocity to momentum space
358  for(unsigned char k=0; k < MFQ19; k++) {
359  A[k]=0.0;
360  for(unsigned char n=0; n < MFQ19; n++)
361  A[k]+=f(n)*M[k][n]; //M
362  }
363 //Relaxation in the momentum space
364  for(unsigned char k=0; k < MFQ19; k++) {
365  if(s[k]!=0.0) {
366  MEQ=Meq_Incompr(Vector4, k);
367  A[k]=A[k]-s[k]*(A[k]-MEQ);
368  }
369  A[k]/=d[k];
370  }
371 //Transformation of results back to the velocity space
372  for(unsigned char k=0; k < MFQ19; k++) {
373  f(k)=0.0;
374  for(unsigned char n=0; n < MFQ19; n++)
375  f(k)+=A[n]*M[n][k]; //MT (M transposed)
376  }
377  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
378  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
379  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
380  }
381 }
382 
383 inline double MF::Solver_CPU::Collision::Meq_Qcompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
384  double x;
385 
386 // Equivalents to LBGK
387  const double alpha = 3.0;
388  const double beta = -11.0/2.0;
389  const double gamma = -0.5;
390 
391 //Optimized values of parameters according to D. D'Humieres et al 2002
392 //const double alpha=0.0;
393 //const double beta=-475.0/63.0;
394 //const double gamma=0.0;
395 
396  switch (k)
397  {
398  case 0: {x=RHO; break;}
399  case 1: {x=-11.0*RHO+19.0*RHO*(U*U+V*V+W*W);break;}
400  case 2: {x=RHO*alpha+beta*RHO*(U*U+V*V+W*W);break;}
401  case 3: {x=RHO*U;break;}
402  case 4: {x=-2.0/3.0*RHO*U;break;}
403  case 5: {x=RHO*V;break;}
404  case 6: {x=-2.0/3.0*RHO*V;break;}
405  case 7: {x=RHO*W;break;}
406  case 8: {x=-2.0/3.0*RHO*W;break;}
407  case 9: {x=RHO*(2.0*U*U-(V*V+W*W));break;}
408  case 10: {x=gamma*RHO*(2.0*U*U-(V*V+W*W));break;}
409  case 11: {x=RHO*(V*V-W*W);break;}
410  case 12: {x=gamma*RHO*(V*V-W*W);break;}
411  case 13: {x=RHO*(U*V);break;}
412  case 14: {x=RHO*(V*W);break;}
413  case 15: {x=RHO*(U*W);break;}
414  case 16: {x=0.0;break;}
415  case 17: {x=0.0;break;}
416  case 18: {x=0.0;break;}
417  default: {x=0.0;break;}
418  }
419  return (x);
420 }
421 
422 //MRT uncompressible He Luo 1997
423 inline double MF::Solver_CPU::Collision::Meq_Incompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
424  double x;
425  const double rho0_LB = m_CaseParameters_Ptr->ReferenceDensityLB_Rho0_LB;
426 
427 //Equivalents to LBGK
428  const double alpha = 3.0;
429  const double beta = -11.0/2.0;
430  const double gamma = -0.5;
431 
432 //Optimized values of parameters according to D. D'Humieres et al 2002
433 //const double alpha=0.0;
434 //const double beta=-475.0/63.0;
435 //const double gamma=0.0;
436 
437  switch (k)
438  {
439  case 0: {x=RHO; break;}
440  case 1: {x=-11.0*RHO+19.0*RHO*RHO*(U*U+V*V+W*W)/rho0_LB;break;}
441  case 2: {x=RHO*alpha+beta*RHO*RHO*(U*U+V*V+W*W)/rho0_LB;break;}
442  case 3: {x=RHO*U;break;}
443  case 4: {x=-2.0/3.0*RHO*U;break;}
444  case 5: {x=RHO*V;break;}
445  case 6: {x=-2.0/3.0*RHO*V;break;}
446  case 7: {x=RHO*W;break;}
447  case 8: {x=-2.0/3.0*RHO*W;break;}
448  case 9: {x=RHO*RHO*(2.0*U*U-(V*V+W*W))/rho0_LB;break;}
449  case 10: {x=gamma*RHO*RHO*(2.0*U*U-(V*V+W*W))/rho0_LB;break;}
450  case 11: {x=RHO*RHO*(V*V-W*W)/rho0_LB;break;}
451  case 12: {x=gamma*RHO*RHO*(V*V-W*W)/rho0_LB;break;}
452  case 13: {x=RHO*RHO*U*V/rho0_LB;break;}
453  case 14: {x=RHO*RHO*V*W/rho0_LB;break;}
454  case 15: {x=RHO*RHO*U*W/rho0_LB;break;}
455  case 16: {x=0.0;break;}
456  case 17: {x=0.0;break;}
457  case 18: {x=0.0;break;}
458  default: {x=0.0;break;}
459  }
460  return (x);
461 }
462 
463 //-----------------------------------------------------------------------------------------------------------------------------
464 // MRT 2------------------------------------------------------------------------------------------------------------------------
465 //-----------------------------------------------------------------------------------------------------------------------------
466 
467 void MF::Solver_CPU::Collision::Coll_MRT2_Qcompr(MF::Database::Node * pNode) {
468  double MEQ;
469  double A[MFQ19];
470  const double TAU = m_CaseParameters_Ptr->Tau;
471  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
472  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
473 
474  // From OpenLB. Use these relaxation times for higher stability
475  const double s[MFQ19] = {
476  /*s0*/ 0, // rho (conserved)
477  /*s1*/ 1, // bulk viscosity
478  /*s2*/ 1,
479  /*s3*/ 0, // rho*ux (conserved)
480  /*s4*/ 1,
481  /*s5*/ 0, // rho*uy (conserved)
482  /*s6*/ 1, // = s4
483  /*s7*/ 0, // rho*uz (conserved)
484  /*s8*/ 1, // = s4
485  /*s9*/ 1/TAU, //should be equal to s13, used to define nu
486  /*s10*/ 1,
487  /*s11*/ 1/TAU, // = s9,
488  /*s12*/ 1,
489  /*s13*/ 1/TAU, //should be equal to s9, used to define nu
490  /*s14*/ 1/TAU, // = s13,
491  /*s15*/ 1/TAU, // = s13,
492  /*s16*/ 1,
493  /*s17*/ 1, // = s16,
494  /*s18*/ 1 // = s16,
495  };
496 
497 //Matrix of coefficients for converting velocities to the momentum space and vice versa
498  const double M[MFQ19][MFQ19] = {
499  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
500  {-30, -11, -11, -11, -11, -11, -11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8},
501  {12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
502  {0, 1, 0, -1, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
503  {0, -4, 0, 4, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
504  {0, 0, 1, 0, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
505  {0, 0, -4, 0, 4, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
506  {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
507  {0, 0, 0, 0, 0, -4, 4, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
508  {0, 2, -1, 2, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
509  {0, -4, 2, -4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
510  {0, 0, 1, 0, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
511  {0, 0, -2, 0, -2, 2, 2, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
512  {0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0},
513  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1},
514  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0},
515  {0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, -1, 1, 1, 0, 0, 0, 0},
516  {0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 1, 1, -1, -1},
517  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, 1, 1, -1}
518  };
519 
520  const double d[MFQ19]={19.0, 2394.0, 252.0, 10.0, 40.0, 10.0, 40.0, 10.0, 40.0, 36.0, 72.0, 12.0, 24.0, 4.0, 4.0, 4.0, 8.0, 8.0, 8.0}; //D=M*MT
521 
522 //Transformation from the velocity to momentum space
523  for(unsigned char k=0; k < MFQ19; k++) {
524  A[k]=0.0;
525  for(unsigned char n=0; n < MFQ19; n++)
526  A[k]+=f(n)*M[k][n]; //M
527  }
528 //Relaxation in the momentum space
529  for(unsigned char k=0; k < MFQ19; k++) {
530  if(s[k]!=0.0) {
531  MEQ=Meq_Qcompr(Vector4, k);
532  A[k]=A[k]-s[k]*(A[k]-MEQ);
533  }
534  A[k]/=d[k];
535  }
536 //Transformation of results back to the velocity space
537  for(unsigned char k=0; k < MFQ19; k++) {
538  f(k)=0.0;
539  for(unsigned char n=0; n < MFQ19; n++)
540  f(k)+=A[n]*M[n][k]; //MT (M transposed)
541  }
542  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
543  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
544  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
545  }
546 }
547 
548 void MF::Solver_CPU::Collision::Coll_MRT2_Incompr(MF::Database::Node * pNode) {
549  double MEQ;
550  double A[MFQ19];
551  const double TAU = m_CaseParameters_Ptr->Tau;
552  const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
553  const MF::Database::Vec4<double> Vector4 = (*Thread_Ptr->m_pVRLBFunction)(pNode);
554 
555  // From OpenLB. Use these relaxation times for higher stability
556  const double s[MFQ19] = {
557  0, // rho (conserved)
558  1, // bulk viscosity
559  1,
560  0, // rho*ux (conserved)
561  1,
562  0, // rho*uy (conserved)
563  1, // = s4
564  0, // rho*uz (conserved)
565  1, // = s4
566  1/TAU, //should be equal to s13, used to define nu
567  1,
568  1/TAU, // = s9,
569  1,
570  1/TAU, //should be equal to s9, used to define nu
571  1/TAU, // = s13,
572  1/TAU, // = s13,
573  1,
574  1, // = s16,
575  1 // = s16,
576  };
577 
578 
579 //Matrix of coefficients for converting velocities to the momentum space and vice versa
580  const double M[MFQ19][MFQ19] = {
581  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
582  {-30, -11, -11, -11, -11, -11, -11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8},
583  {12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
584  {0, 1, 0, -1, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
585  {0, -4, 0, 4, 0, 0, 0, 1, -1, -1, 1, 1, 1, -1, -1, 0, 0, 0, 0},
586  {0, 0, 1, 0, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
587  {0, 0, -4, 0, 4, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, 1, -1, -1},
588  {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
589  {0, 0, 0, 0, 0, -4, 4, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1},
590  {0, 2, -1, 2, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
591  {0, -4, 2, -4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, -2, -2, -2, -2},
592  {0, 0, 1, 0, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
593  {0, 0, -2, 0, -2, 2, 2, 1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0},
594  {0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0},
595  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1},
596  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0},
597  {0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, -1, 1, 1, 0, 0, 0, 0},
598  {0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0, 1, 1, -1, -1},
599  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1, -1, 1, 1, -1}
600  };
601 
602  const double d[MFQ19]={19.0, 2394.0, 252.0, 10.0, 40.0, 10.0, 40.0, 10.0, 40.0, 36.0, 72.0, 12.0, 24.0, 4.0, 4.0, 4.0, 8.0, 8.0, 8.0}; //D=M*MT
603 
604 //Transformation from the velocity to momentum space
605  for(unsigned char k=0; k < MFQ19; k++) {
606  A[k]=0.0;
607  for(unsigned char n=0; n < MFQ19; n++)
608  A[k]+=f(n)*M[k][n]; //M
609  }
610 //Relaxation in the momentum space
611  for(unsigned char k=0; k < MFQ19; k++) {
612  if(s[k]!=0.0) {
613  MEQ=Meq_Incompr(Vector4, k);
614  A[k]=A[k]-s[k]*(A[k]-MEQ);
615  }
616  A[k]/=d[k];
617  }
618 //Transformation of results back to the velocity space
619  for(unsigned char k=0; k < MFQ19; k++) {
620  f(k)=0.0;
621  for(unsigned char n=0; n < MFQ19; n++)
622  f(k)+=A[n]*M[n][k]; //MT (M transposed)
623  }
624  // Swaps all opposite f(i) values. Required by swap-propagation procedure J. Latt 2007 tech. note. https://www.openlb.net/wp-content/uploads/2011/12/olb-tr1.pdf
625  for (auto LinkedDirection : MF::GU::LatticeParametersD3Q19::SwapDirections) {
626  std::swap(f(LinkedDirection[0]),f(LinkedDirection[1]));
627  }
628 }
629 
630 double MF::Solver_CPU::Collision::GetBGKFeq_Qcompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
631  return Feq_Qcompr(Vector4,k);
632 }
633 
634 double MF::Solver_CPU::Collision::GetBGKFeq_Incompr(const MF::Database::Vec4<double>& Vector4, const unsigned char& k) {
635  return Feq_Incompr(Vector4,k);
636 }
MF::Solver_CPU::Collision::m_CaseParameters_Ptr
static std::shared_ptr< MF::Solver_CPU::CaseParameters > m_CaseParameters_Ptr
Definition: Collision.h:28
MF::Solver_CPU::Collision::GetBGKFeq_Incompr
static double GetBGKFeq_Incompr(const MF::Database::Vec4< double > &Vector4, const unsigned char &k)
Gets FEQ from quasi incompressible BGK model.
Definition: Collision.cpp:634
V
#define V
Definition: Collision.cpp:15
MF::Database::Vec3::y
T y
Y direction value.
Definition: Vec3.h:20
MFQ19
#define MFQ19
Number of lattice directions D3Q19.
Definition: LatticeParametersD3Q19.h:14
MF::Database::Node::pMyThread_Ptr
std::shared_ptr< MF::Database::Thread > * pMyThread_Ptr
Pointer to Shared pointer to node MFThread parameters.
Definition: Node.h:27
RHO
#define RHO
Definition: Collision.cpp:17
MF::Solver_CPU::FlowTypeClass::Incompressible
MF::Solver_CPU::ModelTypeClass::FBGK
MF::Database::Node
Basic data structure for storing f(i) data for each computational grid node.
Definition: Node.h:26
MF::Solver_CPU::Collision::m_LatticeParameters_Ptr
static std::shared_ptr< MF::GU::LatticeParametersD3Q19 > m_LatticeParameters_Ptr
Definition: Collision.h:29
MF::Database::Vec3::z
T z
Z direction value.
Definition: Vec3.h:21
U
#define U
Definition: Collision.cpp:14
DirectionWeight_W
#define DirectionWeight_W(x)
Definition: Collision.cpp:21
MF::Database::Vec3::x
T x
X direction value.
Definition: Vec3.h:19
DirectionVectorComponent_Cx
#define DirectionVectorComponent_Cx(x)
Definition: Collision.cpp:18
f
#define f(x)
Definition: Collision.cpp:13
MF::Solver_CPU::Collision::GetBGKFeq_Qcompr
static double GetBGKFeq_Qcompr(const MF::Database::Vec4< double > &Vector4, const unsigned char &k)
Gets FEQ from quasi compressible BGK model.
Definition: Collision.cpp:630
MF::Solver_CPU::ModelTypeClass::MRT
MF::Solver_CPU::Collision::SetCollisionPointerToFunc
static void SetCollisionPointerToFunc(const std::shared_ptr< MF::Database::ThreadArray > &ThreadArray_Ptr)
Sets pointers to proper collision functions in MFThreads.
Definition: Collision.cpp:23
W
#define W
Definition: Collision.cpp:16
DirectionVectorComponent_Cz
#define DirectionVectorComponent_Cz(x)
Definition: Collision.cpp:20
DirectionVectorComponent_Cy
#define DirectionVectorComponent_Cy(x)
Definition: Collision.cpp:19
MF::Solver_CPU::ModelTypeClass::MRT2
MF::GU::LatticeParametersD3Q19::SwapDirections
static constexpr uint8_t SwapDirections[MFQ19_H][2]
Indexes of swap directions of f(i) in D3Q19.
Definition: LatticeParametersD3Q19.h:57