Memosa-FVM  0.2
COMETBoundaryConditions.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #ifndef _COMETBOUNDARYCONDITIONS_H_
6 #define _COMETBOUNDARYCONDITIONS_H_
7 #include <stdio.h>
8 #include "Mesh.h"
9 #include <vector>
10 #include <cmath>
11 #include "NumType.h"
12 #include "Array.h"
13 #include "Vector.h"
14 #include "Field.h"
15 #include "StorageSite.h"
16 #include "CRConnectivity.h"
17 #include "MultiFieldMatrix.h"
18 #include "CRMatrix.h"
19 #include "FluxJacobianMatrix.h"
20 #include "DiagonalMatrix.h"
21 #include "GradientModel.h"
22 #include "FluxLimiters.h"
23 #include "Limiters.h"
24 
25 #include "StressTensor.h"
26 template<class X,class Diag,class OffDiag>
27 class COMETBoundaryConditions
28 {
29  public :
30 
31  typedef typename NumTypeTraits<X>::T_Scalar T_Scalar;
32 
33  typedef Array<int> IntArray;
34 
35  typedef Array<T_Scalar> TArray;
36  typedef Vector<T_Scalar,3> VectorT3;
37  typedef Array<VectorT3> VectorT3Array;
38 
39  typedef Array<X> XArray;
40  typedef Vector<X,3> VectorX3; //new for esbgk
41  typedef Array<VectorX3> VectorX3Array;
42  typedef StressTensor<X> StressTensorX6;
43  typedef Array<StressTensorX6> StressTensorArray;
44  typedef Gradient<T_Scalar> GradType;
45  typedef Array<GradType> GradArray;
46  typedef GradientModel<T_Scalar> GradModelType;
47  typedef typename GradModelType::GradMatrixType GradMatrix;
48 
49  COMETBoundaryConditions(const StorageSite& faces,
50  const Mesh& mesh,
51  const GeomFields& geomFields,
52  const Quadrature<X>& quadrature,
53  MacroFields& macroFields,
54  DistFunctFields<X>& dsfPtr):
55 
56  _faces(faces),
57  _allFaces(mesh.getFaces()),
58  _mesh(mesh),
59  _geomFields(geomFields),
60  _cells(mesh.getCells()),
61  _cellFaces(mesh.getCellFaces()),
62  _ibType(dynamic_cast<const IntArray&>(geomFields.ibType[_cells])),
63  _quadrature(quadrature),
64  _macroFields(macroFields),
65  _dsfPtr(dsfPtr),
66  _faceCells(mesh.getFaceCells(_faces)),
67  _allFaceCells(mesh.getAllFaceCells()),
68  _cellCells(mesh.getCellCells()),
69  _areaMagField(geomFields.areaMag),
70  _faceAreaMag(dynamic_cast<const TArray&>(_areaMagField[_faces])),
71  _areaField(geomFields.area),
72  _faceArea(dynamic_cast<const VectorT3Array&>(_areaField[_faces]))
73 
74  {}
75  COMETModelOptions<X>& getOptions() {return _options;}
76 
77  void applyPressureInletBC(int f,const X& inletTemperature,const X& inletPressure,GradMatrix& gMat) const
78  {
79 
80  int surface(-1);
81  const double pi=_options.pi;
82  const int c0 = _faceCells(f,0);
83  const int c1 = _faceCells(f,1);
84 
85  if (_ibType[c0] != Mesh::IBTYPE_FLUID)
86  return;
87 
88  GradArray Grads(_quadrature.getDirCount());
89  Grads.zero();
90 
91  VectorT3 Gcoeff;
92  TArray limitCoeff1(_quadrature.getDirCount());
93  TArray min1(_quadrature.getDirCount());
94  TArray max1(_quadrature.getDirCount());
95 
96  for(int dir=0;dir<_quadrature.getDirCount();dir++)
97  {
98  limitCoeff1[dir]=T_Scalar(1.e20);
99  Field& fnd = *_dsfPtr.dsf[dir];
100  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
101  min1[dir]=dsf[c0];
102  max1[dir]=dsf[c0];
103  }
104 
105  const VectorT3Array& faceCoords=
106  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_allFaces]);
107  const VectorT3Array& cellCoords=
108  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
109 
110 
111  const int neibcount=_cellFaces.getCount(c0);
112  for(int j=0;j<neibcount;j++)
113  {
114 
115  const int face=_cellFaces(c0,j);
116 
117  int cell2=_allFaceCells(face,1);
118  if(cell2==c0)
119  cell2=_allFaceCells(face,0);
120 
121  if(cell2==c1)surface=face;
122  Gcoeff=gMat.getCoeff(c0, cell2);
123  for(int dir=0;dir<_quadrature.getDirCount();dir++)
124  {
125  Field& fnd = *_dsfPtr.dsf[dir];
126  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
127  Grads[dir].accumulate(Gcoeff,dsf[cell2]-dsf[c0]);
128  if(min1[dir]>dsf[cell2])min1[dir]=dsf[cell2];
129  if(max1[dir]<dsf[cell2])max1[dir]=dsf[cell2];
130  }
131  }
132 
133  for(int j=0;j<neibcount;j++)
134  {
135  const int face=_cellFaces(c0,j);
136 
137  SuperbeeLimiter lf;
138  for(int dir=0;dir<_quadrature.getDirCount();dir++)
139  {
140  Field& fnd = *_dsfPtr.dsf[dir];
141  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
142  GradType& grad=Grads[dir];
143  VectorT3 fVec=faceCoords[face]-cellCoords[c0];
144  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
145  computeLimitCoeff(limitCoeff1[dir], dsf[c0], SOU, min1[dir], max1[dir], lf);
146  }
147  }
148 
149  const int numDirections = _quadrature.getDirCount();
150  const XArray& cx = dynamic_cast<const XArray&>(*_quadrature.cxPtr);
151  const XArray& cy = dynamic_cast<const XArray&>(*_quadrature.cyPtr);
152  const XArray& cz = dynamic_cast<const XArray&>(*_quadrature.czPtr);
153  const XArray& wts = dynamic_cast<const XArray&>(*_quadrature.dcxyzPtr);
154  XArray& density = dynamic_cast<XArray&>(_macroFields.density[_cells]);
155  XArray& pressure = dynamic_cast<XArray&>(_macroFields.pressure[_cells]);
156  XArray& temperature = dynamic_cast<XArray&>(_macroFields.temperature[_cells]);
157  //XArray& heatFlux = dynamic_cast<XArray&>(_macroFields.heatFlux[_cells]);
158  VectorX3Array& v = dynamic_cast<VectorX3Array&>(_macroFields.velocity[_cells]);
159 
160 
161  X uwallnew = 0.0;
162 
163  const X Tin = inletTemperature;
164  const X Pin = inletPressure;
165  const X nin = Pin/Tin; // wall number density for initializing Maxwellian
166 
167  //store boundary values
168  temperature[c1]=inletTemperature;
169  pressure[c1]=inletPressure;
170  density[c1]=nin;
171  v[c1][0]=0.0;
172  v[c1][1]=0.0;
173  v[c1][2]=0.0;
174 
175  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
176 
177 
178  //update velocity
179  for (int j=0; j<numDirections; j++)
180  {
181  Field& fnd = *_dsfPtr.dsf[j];
182  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
183  dsf[c1] = nin/pow(pi*Tin,1.5)*exp(-(pow(cx[j]-v[c1][0],2.0)+pow(cy[j]-v[c1][1],2.0)+pow(cz[j]-v[c1][2],2.0))/Tin);
184 
185  const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
186 
187  if (abs(en[0]) == T_Scalar(1.0))
188  {
189  if( c_dot_en > T_Scalar(0.0))
190  {
191  uwallnew=uwallnew+cx[j]*dsf[c0]*wts[j];
192  }
193  else {uwallnew=uwallnew+cx[j]*dsf[c1]*wts[j];}
194  }
195  else if (abs(en[1]) == T_Scalar(1.0))
196  {
197  if( c_dot_en > T_Scalar(0.0))
198  {
199  uwallnew=uwallnew+cy[j]*dsf[c0]*wts[j];
200  }
201  else {uwallnew=uwallnew+cy[j]*dsf[c1]*wts[j];}
202  }
203  else if (abs(en[2]) == T_Scalar(1.0))
204  {
205  if( c_dot_en > T_Scalar(0.0))
206  {
207  uwallnew=uwallnew+cz[j]*dsf[c0]*wts[j];
208  }
209  else {uwallnew=uwallnew+cz[j]*dsf[c1]*wts[j];}
210  }
211  }
212  /*
213  if (abs(en[0]) == T_Scalar(1.0))
214  v[c1][0]=uwallnew/nin;
215  else if (abs(en[1]) == T_Scalar(1.0)) //incoming {vwallnew=vwallnew+cy[j]*dsf[c1]*wts[j];}
216  v[c1][1]=uwallnew/nin;
217  else if(abs(en[2]) == T_Scalar(1.0))
218  v[c1][2]=uwallnew/nin;//uwall=uwall+relax*(uwallnew-uwall);
219  */
220 
221 
222  //update f
223  for (int j=0; j<numDirections; j++)
224  {
225  Field& fnd = *_dsfPtr.dsf[j];
226  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
227  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
228  const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
229  if (c_dot_en< T_Scalar(0.0))
230  {
231  dsf[c1] = nin/pow(pi*Tin,1.5)*exp(-(pow(cx[j]-v[c1][0],2.0)+pow(cy[j]-v[c1][1],2.0)+pow(cz[j]-v[c1][2],2.0))/Tin);
232  // dsf[c0]=dsf[c1];// change value of fluid boundary-cell
233  }
234  else
235  {
236  GradType& grad=Grads[j];
237  VectorT3 fVec=faceCoords[surface]-cellCoords[c0];
238  VectorT3 rVec=cellCoords[c1]-cellCoords[c0];
239  T_Scalar r=gMat.computeR(grad,dsf,rVec,c0,c1);
240  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
241  //dsf[c1]=dsf[c0];
242  //dsf[c1]=dsf[c0]+superbee(r)*SOU;
243  dsf[c1]=dsf[c0]+limitCoeff1[j]*SOU;
244  }
245  }
246 
247  }
248 
249  void applyPressureInletBC(const FloatValEvaluator<X>& bTemperature,const FloatValEvaluator<X>& bPresssure) const
250  {
251  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
252  for (int i=0; i<_faces.getCount();i++)
253  applyPressureInletBC(i,bTemperature[i],bPresssure[i],gradMatrix);
254 
255  }
256 
257  //outlet Boundary Condition
258  void applyPressureOutletBC(int f,const X& outletTemperature,const X& outletPressure) const
259  {
260 
261  const double pi=_options.pi;
262  const int c0 = _faceCells(f,0);
263  const int c1 = _faceCells(f,1);
264 
265  if (_ibType[c0] != Mesh::IBTYPE_FLUID)
266  return;
267 
268  const int numDirections = _quadrature.getDirCount();
269  const XArray& cx = dynamic_cast<const XArray&>(*_quadrature.cxPtr);
270  const XArray& cy = dynamic_cast<const XArray&>(*_quadrature.cyPtr);
271  const XArray& cz = dynamic_cast<const XArray&>(*_quadrature.czPtr);
272  XArray& density = dynamic_cast<XArray&>(_macroFields.density[_cells]);
273  VectorX3Array& v = dynamic_cast<VectorX3Array&>(_macroFields.velocity[_cells]);
274  XArray& pressure = dynamic_cast<XArray&>(_macroFields.pressure[_cells]);
275  XArray& temperature = dynamic_cast<XArray&>(_macroFields.temperature[_cells]);
276 
277  // X Tout = outletTemperature;
278  const X Pout = outletPressure;
279  X Pcell= pressure[c0];
280  X gamma =_options.SpHeatRatio;
281  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
282  X nout =density[c0];
283 
284  if (Pcell > Pout){
285  const X avel2 =gamma*Pcell/density[c0]; //velcotiy of sound square
286  nout =density[c0]-(Pcell-Pout)/avel2;
287  X ubulk=pow(pow(v[c0][0],2.0)+pow(v[c0][1],2.0)+pow(v[c0][2],2.0),0.5);
288  X ucoeff=1.0/ubulk*(ubulk+(Pcell-Pout)/(pow(2.0*avel2,0.5)*density[c0]));
289  if(abs(en[0])==T_Scalar(1.0))
290  v[c1][0] = v[c0][0]*ucoeff;
291  //v[c1][0]=v[c0][0]+(pressure[c0]-Pout)/(pow(2.0*avel2,0.5)*density[c0]); //exit velocity
292  else if(abs(en[1])==T_Scalar(1.0))
293  v[c1][1] = v[c0][1]*ucoeff;
294  else if(abs(en[2])==T_Scalar(1.0))
295  v[c1][2] = v[c0][2]*ucoeff;
296 
297  density[c1]=nout; pressure[c1]=Pout; //update macroPr
298  temperature[c1]=Pout/nout;
299  }
300  else{
301 
302  density[c1]=density[c0];
303  pressure[c1]=pressure[c0];
304  temperature[c1]=temperature[c0];
305 
306  }
307  for (int j=0; j<numDirections; j++)
308  {
309  Field& fnd = *_dsfPtr.dsf[j];
310  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
311 
312  const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
313  if (c_dot_en < T_Scalar(0.0))
314  {
315  dsf[c1] = nout/pow(pi*temperature[c1],1.5)*exp(-(pow(cx[j]-v[c1][0],2.0)+pow(cy[j]-v[c1][1],2.0)+pow(cz[j]-v[c1][2],2.0))/temperature[c1]);
316  // dsf[c0]=dsf[c1];// change value of fluid boundary-cell
317  }
318  else{dsf[c1]=dsf[c0];}
319 
320  }
321 
322 
323 }
324 
325  void applyPressureOutletBC(const X& bTemperature, const X& bPressure) const
326  {
327  for (int i=0; i<_faces.getCount();i++)
328  applyPressureOutletWallBC(i,bTemperature,bPressure);
329  }
330 
331 
332  void applyPressureOutletBC(const FloatValEvaluator<X>& bTemperature,const FloatValEvaluator<X>& bPresssure) const
333  {
334  for (int i=0; i<_faces.getCount();i++)
335  applyPressureOutletBC(i,bTemperature[i],bPresssure[i]);
336 
337  }
338 
339  // Real wall with momentum accomodation coefficient
340  void applyRealWallBC(int f,const VectorX3& WallVelocity,const X& WallTemperature,const X& accommodationCoefficient,const vector<int>& vecReflection,GradMatrix& gMat) const
341  {
342  int surface(-1);
343  const double pi=_options.pi;
344  //const double epsilon=_options.epsilon_ES;
345 
346  const int c0 = _faceCells(f,0);
347  const int c1 = _faceCells(f,1);
348 
349  if (_ibType[c0] != Mesh::IBTYPE_FLUID)
350  return;
351 
352  GradArray Grads(_quadrature.getDirCount());
353  Grads.zero();
354 
355  VectorT3 Gcoeff;
356  TArray limitCoeff1(_quadrature.getDirCount());
357  TArray min1(_quadrature.getDirCount());
358  TArray max1(_quadrature.getDirCount());
359 
360  for(int dir=0;dir<_quadrature.getDirCount();dir++)
361  {
362  limitCoeff1[dir]=T_Scalar(1.e20);
363  Field& fnd = *_dsfPtr.dsf[dir];
364  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
365  min1[dir]=dsf[c0];
366  max1[dir]=dsf[c0];
367  }
368 
369  const VectorT3Array& faceCoords=
370  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_allFaces]);
371  const VectorT3Array& cellCoords=
372  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
373 
374  const int neibcount=_cellFaces.getCount(c0);
375  for(int j=0;j<neibcount;j++)
376  {
377 
378  const int face=_cellFaces(c0,j);
379 
380  int cell2=_allFaceCells(face,1);
381  if(cell2==c0)
382  cell2=_allFaceCells(face,0);
383 
384  if(cell2==c1)surface=face;
385  Gcoeff=gMat.getCoeff(c0, cell2);
386  for(int dir=0;dir<_quadrature.getDirCount();dir++)
387  {
388  Field& fnd = *_dsfPtr.dsf[dir];
389  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
390  Grads[dir].accumulate(Gcoeff,dsf[cell2]-dsf[c0]);
391  if(min1[dir]>dsf[cell2])min1[dir]=dsf[cell2];
392  if(max1[dir]<dsf[cell2])max1[dir]=dsf[cell2];
393  }
394  }
395 
396  for(int j=0;j<neibcount;j++)
397  {
398  const int face=_cellFaces(c0,j);
399 
400  SuperbeeLimiter lf;
401  for(int dir=0;dir<_quadrature.getDirCount();dir++)
402  {
403  Field& fnd = *_dsfPtr.dsf[dir];
404  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
405  GradType& grad=Grads[dir];
406  VectorT3 fVec=faceCoords[face]-cellCoords[c0];
407  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
408  computeLimitCoeff(limitCoeff1[dir], dsf[c0], SOU, min1[dir], max1[dir], lf);
409  }
410  }
411 
412  const int numDirections = _quadrature.getDirCount();
413  const XArray& cx = dynamic_cast<const XArray&>(*_quadrature.cxPtr);
414  const XArray& cy = dynamic_cast<const XArray&>(*_quadrature.cyPtr);
415  const XArray& cz = dynamic_cast<const XArray&>(*_quadrature.czPtr);
416  const XArray& wts= dynamic_cast<const XArray&>(*_quadrature.dcxyzPtr);
417  VectorX3Array& v = dynamic_cast<VectorX3Array&>(_macroFields.velocity[_cells]);
418  XArray& density = dynamic_cast<XArray&>(_macroFields.density[_cells]);
419  XArray& temperature = dynamic_cast<XArray&>(_macroFields.temperature[_cells]);
420  const X uwall = WallVelocity[0];
421  const X vwall = WallVelocity[1];
422  const X wwall = WallVelocity[2];
423  const X Twall = WallTemperature;
424  const X alpha = accommodationCoefficient;
425  X m1alpha=1.0-alpha;
426  v[c1][0]=uwall;
427  v[c1][1]=vwall;
428  v[c1][2]=wwall;
429 
430  temperature[c1]=Twall;
431  X Nmr(0.0) ;
432  X Dmr(0.0) ;
433  X incomFlux(0.0);
434  for (int j=0; j<numDirections; j++)
435  {
436  Field& fnd = *_dsfPtr.dsf[j];
437  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
438  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
439  const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
440  const X wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
441  const X fwall = 1.0/pow(pi*Twall,1.5)*exp(-(pow(cx[j]-uwall,2.0)+pow(cy[j]-vwall,2.0)+pow(cz[j]-wwall,2.0))/Twall);
442  if (c_dot_en -wallV_dot_en < T_Scalar(0.0)) //incoming
443  {
444  Dmr = Dmr - fwall*wts[j]*(c_dot_en-wallV_dot_en);
445  }
446  else
447  {
448  GradType& grad=Grads[j];
449  VectorT3 fVec=faceCoords[surface]-cellCoords[c0];
450  VectorT3 rVec=cellCoords[c1]-cellCoords[c0];
451  T_Scalar r=gMat.computeR(grad,dsf,rVec,c0,c1);
452  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
453  Nmr = Nmr + (dsf[c0]+limitCoeff1[j]*SOU)*wts[j]*(c_dot_en-wallV_dot_en);
454  }
455  }
456 
457  const X nwall = Nmr/Dmr; // wall number density for initializing Maxwellian
458  density[c1]=nwall;
459 
460  for (int j=0; j<numDirections; j++)
461  {
462  Field& fnd = *_dsfPtr.dsf[j];
463  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
464  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
465  const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
466  const X wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
467 
468  if (c_dot_en-wallV_dot_en < T_Scalar(0.0))
469  {
470  const int direction_incident = vecReflection[j];
471  Field& fndi = *_dsfPtr.dsf[direction_incident];
472  const XArray& dsfi = dynamic_cast<const XArray&>(fndi[_cells]);
473 
474  dsf[c1] = alpha*nwall/pow(pi*Twall,1.5)*exp(-(pow(cx[j]-uwall,2.0)+pow(cy[j]-vwall,2.0)+pow(cz[j]-wwall,2.0))/Twall)+m1alpha*dsfi[c0];
475  }
476  else
477  {
478  GradType& grad=Grads[j];
479  VectorT3 fVec=faceCoords[surface]-cellCoords[c0];
480  VectorT3 rVec=cellCoords[c1]-cellCoords[c0];
481  T_Scalar r=gMat.computeR(grad,dsf,rVec,c0,c1);
482  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
483  dsf[c1]=dsf[c0]+limitCoeff1[j]*SOU;
484  }
485  }
486  }
487 
488  void applyRealWallBC(const VectorX3& bVelocity,const X& bTemperature,const X& accomCoeff,const vector<int>& vecReflection) const
489  {
490  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
491  for (int i=0; i<_faces.getCount();i++)
492  applyRealWallBC(i,bVelocity,bTemperature,accomCoeff,vecReflection,gradMatrix);
493  }
494 
495 
496  void applyRealWallBC(const FloatValEvaluator<VectorX3>& bVelocity,const FloatValEvaluator<X>& bTemperature,const FloatValEvaluator<X>& accomCoeff,const vector<int>& vecReflection) const
497  {
498  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
499  for (int i=0; i<_faces.getCount();i++)
500  applyRealWallBC(i,bVelocity[i],bTemperature[i],accomCoeff[i],vecReflection,gradMatrix);
501 
502  }
503 
504  void applyZeroGradientBC(int f, GradMatrix& gMat) const
505  {
506  const int c0 = _faceCells(f,0);
507  const int c1 = _faceCells(f,1);
508 
509  if (_ibType[c0] != Mesh::IBTYPE_FLUID)
510  return;
511 
512  int surface(-1);
513 
514  GradArray Grads(_quadrature.getDirCount());
515  Grads.zero();
516 
517  VectorT3 Gcoeff;
518  TArray limitCoeff1(_quadrature.getDirCount());
519  TArray min1(_quadrature.getDirCount());
520  TArray max1(_quadrature.getDirCount());
521 
522  for(int dir=0;dir<_quadrature.getDirCount();dir++)
523  {
524  limitCoeff1[dir]=T_Scalar(1.e20);
525  Field& fnd = *_dsfPtr.dsf[dir];
526  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
527  min1[dir]=dsf[c0];
528  max1[dir]=dsf[c0];
529  }
530 
531  const VectorT3Array& faceCoords=
532  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_allFaces]);
533  const VectorT3Array& cellCoords=
534  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
535 
536  const int neibcount=_cellFaces.getCount(c0);
537  for(int j=0;j<neibcount;j++)
538  {
539 
540  const int face=_cellFaces(c0,j);
541 
542  int cell2=_allFaceCells(face,1);
543  if(cell2==c0)
544  cell2=_allFaceCells(face,0);
545 
546  //const int cell2=_cellCells(c0,j);
547  if(cell2==c1)surface=face;
548  Gcoeff=gMat.getCoeff(c0, cell2);
549  for(int dir=0;dir<_quadrature.getDirCount();dir++)
550  {
551  Field& fnd = *_dsfPtr.dsf[dir];
552  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
553  Grads[dir].accumulate(Gcoeff,dsf[cell2]-dsf[c0]);
554  if(min1[dir]>dsf[cell2])min1[dir]=dsf[cell2];
555  if(max1[dir]<dsf[cell2])max1[dir]=dsf[cell2];
556  }
557  }
558 
559  for(int j=0;j<neibcount;j++)
560  {
561  const int face=_cellFaces(c0,j);
562 
563  SuperbeeLimiter lf;
564  for(int dir=0;dir<_quadrature.getDirCount();dir++)
565  {
566  Field& fnd = *_dsfPtr.dsf[dir];
567  TArray& dsf = dynamic_cast< TArray&>(fnd[_cells]);
568  GradType& grad=Grads[dir];
569 
570  VectorT3 fVec=faceCoords[face]-cellCoords[c0];
571  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
572 
573  computeLimitCoeff(limitCoeff1[dir], dsf[c0], SOU, min1[dir], max1[dir], lf);
574  }
575  }
576 
577  const int numDirections = _quadrature.getDirCount();
578 
579  for (int j=0; j<numDirections; j++)
580  {
581  Field& fnd = *_dsfPtr.dsf[j];
582  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
583  GradType& grad=Grads[j];
584  VectorT3 fVec=faceCoords[surface]-cellCoords[c0];
585  VectorT3 rVec=cellCoords[c1]-cellCoords[c0];
586  T_Scalar r=gMat.computeR(grad,dsf,rVec,c0,c1);
587  T_Scalar SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
588  dsf[c1]=dsf[c0]+limitCoeff1[j]*SOU;
589  }
590  }
591  void applyZeroGradientBC() const
592  {
593  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
594  for (int i=0; i<_faces.getCount();i++)
595  applyZeroGradientBC(i,gradMatrix);
596  }
597 
598  void applyNSInterfaceBC( )const
599 //const FloatValEvaluator<X>& bTemperature,const FloatValEvaluator<X>& bPressure,const FloatValEvaluator<VectorX3>& bVelocity, const FloatValEvaluator<VectorX3>& bStress) const
600  {
601 
602  for (int i=0; i<_faces.getCount();i++)
603  applyNSInterfaceBC(i); //,bTemperature[i],bPressure[i],bVelocity[i],bStress[i]);
604 
605 
606 
607  }
608 
609 
610 
611  void applyNSInterfaceBC(int f)const //X& inletTemperature,const X& inletPressure, const VectorX3& inletVelocity, const VectorX3&tau) const
612  //StressTensor<X>&tau) const
613  {
614 
615  const double pi=_options.pi;
616  const int c0 = _faceCells(f,0);
617  const int c1 = _faceCells(f,1);
618 
619  if (_ibType[c0] != Mesh::IBTYPE_FLUID)
620  return;
621 
622 
623  const int numDirections = _quadrature.getDirCount();
624  const XArray& cx = dynamic_cast<const XArray&>(*_quadrature.cxPtr);
625  const XArray& cy = dynamic_cast<const XArray&>(*_quadrature.cyPtr);
626  const XArray& cz = dynamic_cast<const XArray&>(*_quadrature.czPtr);
627 
628  //XArray& density = dynamic_cast<XArray&>(_macroFields.density[_cells]);
629  //XArray& pressure = dynamic_cast<XArray&>(_macroFields.pressure[_cells]);
630  XArray& temperature = dynamic_cast<XArray&>(_macroFields.temperature[_cells]);
631  //VectorX3Array& v = dynamic_cast<VectorX3Array&>(_macroFields.velocity[_cells]);
632 
633  XArray& viscosity = dynamic_cast<XArray&>(_macroFields.viscosity[_cells]);
634 
635 
636  VectorX3Array& IntVel = dynamic_cast<VectorX3Array&>(_macroFields.InterfaceVelocity[_faces]);
637  XArray& IntPress = dynamic_cast<XArray&>(_macroFields.InterfacePressure[_faces]);
638  //XArray& IntDens = dynamic_cast<XArray&>(_macroFields.InterfaceDensity[_faces]);
639  StressTensorArray& tau = dynamic_cast<StressTensorArray&>(_macroFields.InterfaceStress[_faces]);
640 
641  const X Tin = temperature[c0];
642  const X Pin = IntPress[f];
643  const X nin = Pin/Tin; // wall number density for initializing Maxwellian
644 
645 
646  const X u_in=IntVel[f][0];
647  const X v_in=IntVel[f][1];
648  const X w_in=IntVel[f][2];
649 
650  const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
651 
652  // Extra for CE expression
653  const X rho_init=_options["rho_init"];
654  const X T_init= _options["T_init"];
655  const X R=8314.0/_options["molecularWeight"];
656  const X nondim_length=_options["nonDimLt"];
657  const X mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
658 
659  //update f to Maxwellian
660  for (int j=0; j<numDirections; j++)
661  {
662  Field& fnd = *_dsfPtr.dsf[j];
663  XArray& dsf = dynamic_cast< XArray&>(fnd[_cells]);
664  //const VectorT3 en = _faceArea[f]/_faceAreaMag[f];
665  //const X c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
666  //if (c_dot_en< T_Scalar(0.0))
667  {
668  //cout << "nin " <<nin <<endl;
669  dsf[c1] = nin/pow(pi*Tin,1.5)*exp(-(pow(cx[j]-u_in,2.0)+pow(cy[j]-v_in,2.0)+pow(cz[j]-w_in,2.0))/Tin)
670  *(1-viscosity[c0]/mu0/(nin*pow(Tin,2.0))*(pow(cx[j]-u_in,2.0)*tau[f][0]+pow(cy[j]-v_in,2.0)*tau[f][1]+pow(cz[j]-w_in,2.0)*tau[f][2])
671  +2.0*(cx[j]-u_in)*(cy[j]-v_in)*tau[f][3]+2.0*(cy[j]-v_in)*(cz[j]-w_in)*tau[f][4]+2.0*(cz[j]-w_in)*(cx[j]-u_in)*tau[f][5]); //update f to ChapmannEnskog
672  }
673  //else{dsf[c1]=dsf[c0];}//outgoing
674 
675  }
676  //update f to ChapmannEnskog
677 
678 
679 //dsf[c1]=dsf[c1]*(1-2*cx[j]*cy[j]*slopeL*viscosity[c]/mu0) //velocity tangential only
680 //C2=pow(cx[j]-u_in,2.0)+pow(cy[j]-v_in,2.0)+pow(cz[j]-w_in,2.0)
681 //dsf[c1]=dsf[c1]*(1+viscosity[c1]/(mu0*Pr*nin*pow(Tin,2.0))*((cx[j]-u_in)*dtbdx+(cy[j]-v_in)*dtbdy+(cz[j]-w_in)*dtbdz)*(C2-2.5))
682 //temperature difference too
683  }
684 
685 
686 
687 
688  protected:
689 
690  const StorageSite& _faces;
691  const StorageSite& _allFaces;
692  const Mesh& _mesh;
693  const GeomFields& _geomFields;
694  const StorageSite& _cells;
695  const CRConnectivity& _cellFaces;
696  const Array<int>& _ibType;
697  const Quadrature<X>& _quadrature;
698  MacroFields& _macroFields;
699  const DistFunctFields<X>& _dsfPtr;
700  const CRConnectivity& _faceCells;
701  const CRConnectivity& _allFaceCells;
702  const CRConnectivity& _cellCells;
703  const Field& _areaMagField;
704  const TArray& _faceAreaMag;
705  const Field& _areaField;
706  const VectorT3Array& _faceArea;
707  COMETModelOptions<X> _options;
708 };
709 
710 #endif
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
COMETBoundaryConditions::applyRealWallBC
void applyRealWallBC(const FloatValEvaluator< VectorX3 > &bVelocity, const FloatValEvaluator< X > &bTemperature, const FloatValEvaluator< X > &accomCoeff, const vector< int > &vecReflection) const
Definition: COMETBoundaryConditions.h:496
COMETBoundaryConditions::_faceArea
const VectorT3Array & _faceArea
Definition: COMETBoundaryConditions.h:706
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
COMETBoundaryConditions::applyRealWallBC
void applyRealWallBC(const VectorX3 &bVelocity, const X &bTemperature, const X &accomCoeff, const vector< int > &vecReflection) const
Definition: COMETBoundaryConditions.h:488
COMETBoundaryConditions::StressTensorArray
Array< StressTensorX6 > StressTensorArray
Definition: COMETBoundaryConditions.h:43
Field
Definition: Field.h:14
COMETBoundaryConditions::_faceCells
const CRConnectivity & _faceCells
Definition: COMETBoundaryConditions.h:700
COMETBoundaryConditions::_options
COMETModelOptions< X > _options
Definition: COMETBoundaryConditions.h:707
COMETBoundaryConditions::_dsfPtr
const DistFunctFields< X > & _dsfPtr
Definition: COMETBoundaryConditions.h:699
GradientMatrix::computeR
X computeR(const Gradient< X > &g, const Array< X > &x, const Coord dist, int i, int j) const
Definition: GradientMatrix.h:97
COMETBoundaryConditions::applyPressureInletBC
void applyPressureInletBC(const FloatValEvaluator< X > &bTemperature, const FloatValEvaluator< X > &bPresssure) const
Definition: COMETBoundaryConditions.h:249
Quadrature::cyPtr
TArray * cyPtr
Definition: Quadrature.h:42
Mesh
Definition: Mesh.h:49
COMETBoundaryConditions::getOptions
COMETModelOptions< X > & getOptions()
Definition: COMETBoundaryConditions.h:75
COMETBoundaryConditions::applyZeroGradientBC
void applyZeroGradientBC(int f, GradMatrix &gMat) const
Definition: COMETBoundaryConditions.h:504
Quadrature::getDirCount
int getDirCount() const
Definition: Quadrature.h:56
GradientMatrix::getCoeff
Coord & getCoeff(const int i, const int j)
Definition: GradientMatrix.h:129
MacroFields::InterfaceVelocity
Field InterfaceVelocity
Definition: MacroFields.h:39
COMETBoundaryConditions::applyPressureOutletBC
void applyPressureOutletBC(const X &bTemperature, const X &bPressure) const
Definition: COMETBoundaryConditions.h:325
Quadrature::dcxyzPtr
TArray * dcxyzPtr
Definition: Quadrature.h:52
MacroFields::InterfacePressure
Field InterfacePressure
Definition: MacroFields.h:40
computeLimitCoeff
void computeLimitCoeff(T &lc, const T x, const T &dx, const T &min, const T &max, const LimitFunc &f)
Definition: FluxLimiters.h:72
Quadrature::czPtr
TArray * czPtr
Definition: Quadrature.h:47
MacroFields::viscosity
Field viscosity
Definition: MacroFields.h:20
COMETBoundaryConditions::GradType
Gradient< T_Scalar > GradType
Definition: COMETBoundaryConditions.h:44
MacroFields::temperature
Field temperature
Definition: MacroFields.h:22
COMETBoundaryConditions::_cellFaces
const CRConnectivity & _cellFaces
Definition: COMETBoundaryConditions.h:695
COMETBoundaryConditions::applyPressureOutletBC
void applyPressureOutletBC(int f, const X &outletTemperature, const X &outletPressure) const
Definition: COMETBoundaryConditions.h:258
COMETBoundaryConditions::applyPressureOutletBC
void applyPressureOutletBC(const FloatValEvaluator< X > &bTemperature, const FloatValEvaluator< X > &bPresssure) const
Definition: COMETBoundaryConditions.h:332
COMETBoundaryConditions::_cellCells
const CRConnectivity & _cellCells
Definition: COMETBoundaryConditions.h:702
Quadrature::cxPtr
TArray * cxPtr
Definition: Quadrature.h:37
COMETBoundaryConditions::_allFaceCells
const CRConnectivity & _allFaceCells
Definition: COMETBoundaryConditions.h:701
COMETBoundaryConditions::applyRealWallBC
void applyRealWallBC(int f, const VectorX3 &WallVelocity, const X &WallTemperature, const X &accommodationCoefficient, const vector< int > &vecReflection, GradMatrix &gMat) const
Definition: COMETBoundaryConditions.h:340
COMETBoundaryConditions::T_Scalar
NumTypeTraits< X >::T_Scalar T_Scalar
Definition: COMETBoundaryConditions.h:31
COMETBoundaryConditions::applyNSInterfaceBC
void applyNSInterfaceBC(int f) const
Definition: COMETBoundaryConditions.h:611
GradientModel::getGradientMatrix
static GradMatrixType & getGradientMatrix(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:454
COMETBoundaryConditions::VectorT3
Vector< T_Scalar, 3 > VectorT3
Definition: COMETBoundaryConditions.h:36
MacroFields::velocity
Field velocity
Definition: MacroFields.h:15
COMETBoundaryConditions::_quadrature
const Quadrature< X > & _quadrature
Definition: COMETBoundaryConditions.h:697
MacroFields::pressure
Field pressure
Definition: MacroFields.h:19
MacroFields::InterfaceStress
Field InterfaceStress
Definition: MacroFields.h:41
DistFunctFields::dsf
std::vector< Field * > dsf
Definition: DistFunctFields.h:36
COMETBoundaryConditions::applyPressureInletBC
void applyPressureInletBC(int f, const X &inletTemperature, const X &inletPressure, GradMatrix &gMat) const
Definition: COMETBoundaryConditions.h:77
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
COMETBoundaryConditions::GradMatrix
GradModelType::GradMatrixType GradMatrix
Definition: COMETBoundaryConditions.h:47
COMETBoundaryConditions::COMETBoundaryConditions
COMETBoundaryConditions(const StorageSite &faces, const Mesh &mesh, const GeomFields &geomFields, const Quadrature< X > &quadrature, MacroFields &macroFields, DistFunctFields< X > &dsfPtr)
Definition: COMETBoundaryConditions.h:49
COMETBoundaryConditions::GradModelType
GradientModel< T_Scalar > GradModelType
Definition: COMETBoundaryConditions.h:46
MacroFields::density
Field density
Definition: MacroFields.h:21