Memosa-FVM  0.2
COMETESBGKDiscretizer.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 _COMETESBGKDISCRETIZER_H_
6 #define _COMETESBGKDISCRETIZER_H_
7 
8 #include "Mesh.h"
9 #include "Quadrature.h"
10 #include "DistFunctFields.h"
11 #include "CometMatrix.h"
12 #include "ArrowHeadMatrix.h"
13 #include "Array.h"
14 #include "Vector.h"
15 #include "StorageSite.h"
16 #include "CRConnectivity.h"
17 #include "GeomFields.h"
18 #include "NumType.h"
19 #include "COMETBC.h"
20 #include <math.h>
21 #include <map>
22 #include "MatrixJML.h"
23 #include "SquareMatrixESBGK.h"
24 #include "GradientModel.h"
25 #include "FluxLimiters.h"
26 #include "Limiters.h"
27 
28 #define SQR(x) ((x)*(x))
29 
30 template<class T>
31 class COMETESBGKDiscretizer
32 {
33 
34  public:
35  typedef typename NumTypeTraits<T>::T_Scalar T_Scalar;
36  typedef Vector<T_Scalar,3> VectorT3;
37  typedef Array<VectorT3> VectorT3Array;
38  typedef Array<T_Scalar> TArray;
39  typedef MatrixJML<T> TMatrix;
40  typedef CometMatrix<T> TComet;
41  typedef ArrowHeadMatrix<T,3> TArrow;
42  typedef SquareMatrixESBGK<T> TSquareESBGK;
43  typedef map<int,COMETBC<T>*> COMETBCMap;
44  typedef Array<int> IntArray;
45  typedef Array<bool> BoolArray;
46  typedef Vector<int,2> VecInt2;
47  typedef map<int,VecInt2> FaceToFg;
48  typedef Vector<T,5> VectorT5;
49  typedef Array<VectorT5> VectorT5Array;
50  typedef Vector<T,6> VectorT6;
51  typedef Array<VectorT6> VectorT6Array;
52  typedef Vector<T,10> VectorT10;
53  typedef Array<VectorT10> VectorT10Array;
54  typedef DistFunctFields<T> TDistFF;
55  typedef Quadrature<T> TQuad;
56  typedef Gradient<T> GradType;
57  typedef Array<GradType> GradArray;
58  typedef GradientModel<T> GradModelType;
59  typedef typename GradModelType::GradMatrixType GradMatrix;
60 
61  COMETESBGKDiscretizer(const Mesh& mesh, const GeomFields& geomfields, const StorageSite& solidFaces,
62  MacroFields& macroFields, TQuad& quadrature, TDistFF& dsfPtr,
63  TDistFF& dsfPtr1, TDistFF& dsfPtr2, TDistFF& dsfEqPtrES, TDistFF& dsfPtrRes, TDistFF& dsfPtrFAS,
64  const T dT, const int order, const bool transient,const T underRelaxation,
65  const T rho_init, const T T_init, const T MW, const int conOrder,
66  COMETBCMap& bcMap, map<int, vector<int> > faceReflectionArrayMap,
67  const IntArray& BCArray, const IntArray& BCfArray,const IntArray& ZCArray):
68  _mesh(mesh),
69  _geomFields(geomfields),
70  _cells(mesh.getCells()),
71  _faces(mesh.getFaces()),
72  _solidFaces(solidFaces),
73  _cellFaces(mesh.getCellFaces()),
74  _faceCells(mesh.getAllFaceCells()),
75  _areaMagField(_geomFields.areaMag),
76  _faceAreaMag((dynamic_cast<const TArray&>(_areaMagField[_faces]))),
77  _areaField(_geomFields.area),
78  _faceArea(dynamic_cast<const VectorT3Array&>(_areaField[_faces])),
79  _cellVolume(dynamic_cast<const TArray&>(_geomFields.volume[_cells])),
80  _macroFields(macroFields),
81  _quadrature(quadrature),
82  _dsfPtr(dsfPtr),
83  _dsfPtr1(dsfPtr1),
84  _dsfPtr2(dsfPtr2),
85  _dsfEqPtrES(dsfEqPtrES),
86  _dsfPtrRes(dsfPtrRes),
87  _dsfPtrFAS(dsfPtrFAS),
88  _dT(dT),
89  _order(order),
90  _transient(transient),
91  _underRelaxation(underRelaxation),
92  _rho_init(rho_init),
93  _T_init(T_init),
94  _MW(MW),
95  _conOrder(conOrder),
96  _bcMap(bcMap),
97  _faceReflectionArrayMap(faceReflectionArrayMap),
98  _BCArray(BCArray),
99  _BCfArray(BCfArray),
100  _ZCArray(ZCArray),
101  _aveResid(-1.),
102  _residChange(-1.),
103  _fgFinder(),
104  _numDir(_quadrature.getDirCount()),
105  _cx(dynamic_cast<const TArray&>(*_quadrature.cxPtr)),
106  _cy(dynamic_cast<const TArray&>(*_quadrature.cyPtr)),
107  _cz(dynamic_cast<const TArray&>(*_quadrature.czPtr)),
108  _wts(dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr)),
109  _density(dynamic_cast<TArray&>(_macroFields.density[_cells])),
110  _pressure(dynamic_cast<TArray&>(_macroFields.pressure[_cells])),
111  _velocity(dynamic_cast<VectorT3Array&>(_macroFields.velocity[_cells])),
112  _velocityResidual(dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[_cells])),
113  _velocityFASCorrection(0),
114  _temperature(dynamic_cast<TArray&>(_macroFields.temperature[_cells])),
115  _stress(dynamic_cast<VectorT6Array&>(_macroFields.Stress[_cells])),
116  _collisionFrequency(dynamic_cast<TArray&>(_macroFields.collisionFrequency[_cells])),
117  _coeffg(dynamic_cast<VectorT10Array&>(_macroFields.coeffg[_cells]))
118 
119  {
120  _fArrays = new TArray*[_numDir];
121  _fN1Arrays = new TArray*[_numDir];
122  _fN2Arrays = new TArray*[_numDir];
123  _fEqESArrays = new TArray*[_numDir];
124  _fResArrays = new TArray*[_numDir];
125  _fasArrays = new TArray*[_numDir];
126 
127  for(int direction=0;direction<_numDir;direction++)
128  {
129  Field& fnd = *_dsfPtr.dsf[direction];
130  Field& fN1nd = *_dsfPtr1.dsf[direction];
131  Field& fN2nd = *_dsfPtr2.dsf[direction];
132  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
133  Field& fndRes = *_dsfPtrRes.dsf[direction];
134  Field& fndFAS = *_dsfPtrFAS.dsf[direction];
135 
136  _fArrays[direction] = &dynamic_cast<TArray&>(fnd[_cells]);
137  _fEqESArrays[direction] = &dynamic_cast<TArray&>(fndEqES[_cells]);
138  _fResArrays[direction] = &dynamic_cast<TArray&>(fndRes[_cells]);
139  if (fN1nd.hasArray(_cells))
140  _fN1Arrays[direction] = &dynamic_cast<TArray&>(fN1nd[_cells]);
141  if (fN2nd.hasArray(_cells))
142  _fN2Arrays[direction] = &dynamic_cast<TArray&>(fN2nd[_cells]);
143  if (fndFAS.hasArray(_cells))
144  _fasArrays[direction] = &dynamic_cast<TArray&>(fndFAS[_cells]);
145  }
146 
147  if (_macroFields.velocityFASCorrection.hasArray(_cells))
148  _velocityFASCorrection = &dynamic_cast<VectorT3Array&>
149  (_macroFields.velocityFASCorrection[_cells]);
150 
151  }
152 
153  void COMETSolveFine(const int sweep, const int level)
154  {
155  const int cellcount=_cells.getSelfCount();
156  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
157  int start;
158 
159  if(sweep==1)
160  start=0;
161  if(sweep==-1)
162  start=cellcount-1;
163 
164  TArray Bvec(_numDir+3);
165  TArray Resid(_numDir+3);
166  TArrow AMat(_numDir+3);
167 
168  TArray fVal(_numDir);
169 
170  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
171 
172  for(int c=start;((c<cellcount)&&(c>-1));c+=sweep)
173  {
174  if (ibType[c] == Mesh::IBTYPE_FLUID){
175  for(int dir=0;dir<_numDir;dir++)
176  fVal[dir]=(*_fArrays[dir])[c];
177  if(_BCArray[c]==0)
178  {
179  Bvec.zero();
180  Resid.zero();
181  AMat.zero();
182 
183  if(_transient)
184  COMETUnsteady(c,&AMat,Bvec);
185 
186  COMETConvectionFine(c,AMat,Bvec,cellcount,gradMatrix);
187  COMETTest(c,&AMat,Bvec,fVal);
188  //COMETCollision(c,&AMat,Bvec);
189  //COMETMacro(c,&AMat,Bvec);
190 
191  if(level>0)
192  addFAS(c,Bvec);
193 
194  Resid=Bvec;
195 
196  AMat.Solve(Bvec);
197  Distribute(c,Bvec,Resid);
198  setBoundaryValFine(c,cellcount,gradMatrix);
199  //ComputeMacroparameters(c);
200  }
201  else if(_BCArray[c]==1)
202  {
203  Bvec.zero();
204  Resid.zero();
205  AMat.zero();
206 
207  if(_transient)
208  COMETUnsteady(c,&AMat,Bvec);
209 
210  COMETConvectionFine(c,AMat,Bvec,gradMatrix);
211  COMETTest(c,&AMat,Bvec,fVal);
212  //COMETCollision(c,&AMat,Bvec);
213  //COMETMacro(c,&AMat,Bvec);
214 
215  if(level>0)
216  addFAS(c,Bvec);
217 
218  Resid=Bvec;
219 
220  AMat.Solve(Bvec);
221  Distribute(c,Bvec,Resid);
222  setBoundaryValFine(c,cellcount,gradMatrix);
223  //ComputeMacroparameters(c);
224  }
225  else
226  throw CException("Unexpected value for boundary cell map.");
227  }
228  }
229  }
230 
231  void COMETSolve(const int sweep, const int level)
232  {
233  const int cellcount=_cells.getSelfCount();
234  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
235  int start;
236 
237  if(sweep==1)
238  start=0;
239  if(sweep==-1)
240  start=cellcount-1;
241 
242  TArray Bvec(_numDir+3);
243  TArray Resid(_numDir+3);
244  TArrow AMat(_numDir+3);
245 
246  TArray fVal(_numDir);
247 
248  for(int c=start;((c<cellcount)&&(c>-1));c+=sweep)
249  {
250  if (ibType[c] == Mesh::IBTYPE_FLUID){
251  for(int dir=0;dir<_numDir;dir++)
252  fVal[dir]=(*_fArrays[dir])[c];
253  if(_BCArray[c]==0)
254  {
255  Bvec.zero();
256  Resid.zero();
257  AMat.zero();
258 
259  if(_transient)
260  COMETUnsteady(c,&AMat,Bvec);
261 
262  COMETConvection(c,AMat,Bvec,cellcount);
263  COMETTest(c,&AMat,Bvec,fVal);
264  //COMETCollision(c,&AMat,Bvec);
265  //COMETMacro(c,&AMat,Bvec);
266 
267  if(level>0)
268  addFAS(c,Bvec);
269 
270  Resid=Bvec;
271 
272  AMat.Solve(Bvec);
273  Distribute(c,Bvec,Resid);
274  //ComputeMacroparameters(c);
275  }
276  else if(_BCArray[c]==1)
277  {
278  Bvec.zero();
279  Resid.zero();
280  AMat.zero();
281 
282  if(_transient)
283  COMETUnsteady(c,&AMat,Bvec);
284 
285  COMETConvection(c,AMat,Bvec);
286  COMETTest(c,&AMat,Bvec,fVal);
287  //COMETCollision(c,&AMat,Bvec);
288  //COMETMacro(c,&AMat,Bvec);
289 
290  if(level>0)
291  addFAS(c,Bvec);
292 
293  Resid=Bvec;
294 
295  AMat.Solve(Bvec);
296  Distribute(c,Bvec,Resid);
297  //ComputeMacroparameters(c);
298  }
299  else
300  throw CException("Unexpected value for boundary cell map.");
301  }
302  }
303  }
304 
305  template<class MatrixType>
306  void COMETUnsteady(const int cell, MatrixType Amat, TArray& BVec)
307  {
308  const T two(2.0);
309  const T onePointFive(1.5);
310  const T pointFive(0.5);
311 
312  int count = 1;
313 
314  for(int direction=0;direction<_numDir;direction++)
315  {
316  const T fbydT = _cellVolume[cell]/_dT; //pow(_nonDimLength,3);
317  Field& fnd = *_dsfPtr.dsf[direction];
318  const TArray& f = dynamic_cast<const TArray&>(fnd[_cells]);
319  Field& fN1nd = *_dsfPtr1.dsf[direction];
320  const TArray& fN1 = dynamic_cast<const TArray&>(fN1nd[_cells]);
321  Field& fN2nd = *_dsfPtr2.dsf[direction];
322  const TArray& fN2 = dynamic_cast<const TArray&>(fN2nd[_cells]);
323  if(_order>1)
324  {
325  Amat->getElement(count,count) -= fbydT*(onePointFive*f[cell]- two*fN1[cell]
326  + pointFive*fN2[cell]);
327  BVec[count-1] -= fbydT*onePointFive;
328  }
329  else
330  {
331  Amat->getElement(count,count) -= fbydT;
332  BVec[count-1] -= fbydT*(f[cell]- fN1[cell]);
333  }
334  count++;
335  }
336  }
337 
338  void COMETConvectionFine(const int cell, TArrow& Amat, TArray& BVec, const int cellcount, const GradMatrix& gMat)
339  {
340  const int neibcount=_cellFaces.getCount(cell);
341  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
342  const StorageSite& ibFaces = _mesh.getIBFaces();
343  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[_faces]);
344  GradArray Grads(_numDir);
345  Grads.zero();
346 
347  VectorT3 Gcoeff;
348  TArray limitCoeff1(_numDir);
349  TArray min1(_numDir);
350  TArray max1(_numDir);
351  TArray limitCoeff2(_numDir);
352  TArray min2(_numDir);
353  TArray max2(_numDir);
354  for(int dir=0;dir<_numDir;dir++)
355  {
356  limitCoeff1[dir]=T(1.e20);
357  const TArray& f = *_fArrays[dir];
358  min1[dir]=f[cell];
359  max1[dir]=f[cell];
360  }
361 
362  const VectorT3Array& faceCoords=
363  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_faces]);
364  const VectorT3Array& cellCoords=
365  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
366 
367  for(int j=0;j<neibcount;j++)
368  {
369  const int face=_cellFaces(cell,j);
370  int cell2=_faceCells(face,1);
371  if(cell2==cell)
372  cell2=_faceCells(face,0);
373 
374  Gcoeff=gMat.getCoeff(cell, cell2);
375 
376  for(int dir=0;dir<_numDir;dir++)
377  {
378  const TArray& f = *_fArrays[dir];
379  Grads[dir].accumulate(Gcoeff,f[cell2]-f[cell]);
380  if(min1[dir]>f[cell2])min1[dir]=f[cell2];
381  if(max1[dir]<f[cell2])max1[dir]=f[cell2];
382  }
383  }
384 
385  for(int j=0;j<neibcount;j++)
386  {
387  const int face=_cellFaces(cell,j);
388 
389  SuperbeeLimiter lf;
390  for(int dir=0;dir<_numDir;dir++)
391  {
392  const TArray& f = *_fArrays[dir];
393  GradType& grad=Grads[dir];
394 
395  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
396  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
397 
398  computeLimitCoeff(limitCoeff1[dir], f[cell], SOU, min1[dir], max1[dir], lf);
399  }
400  }
401 
402  for(int j=0;j<neibcount;j++)
403  {
404  const int face=_cellFaces(cell,j);
405  int cell2=_faceCells(face,1);
406  VectorT3 Af=_faceArea[face];
407  VectorT3 en = _faceArea[face]/_faceAreaMag[face];
408 
409  T flux;
410 
411  if(cell2==cell)
412  {
413  Af=Af*(-1.);
414  en=en*(-1.);
415  cell2=_faceCells(face,0);
416  }
417 
418  for(int dir=0;dir<_numDir;dir++)
419  {
420  limitCoeff2[dir]=T(1.e20);
421  const TArray& f = *_fArrays[dir];
422  min2[dir]=f[cell2];
423  max2[dir]=f[cell2];
424  }
425 
426  GradArray NeibGrads(_numDir);
427  NeibGrads.zero();
428 
429  const int neibcount1=_cellFaces.getCount(cell2);
430  for(int nj=0;nj<neibcount1;nj++)
431  {
432  const int f1=_cellFaces(cell2,nj);
433  int cell22=_faceCells(f1,1);
434  if(cell2==cell22)
435  cell22=_faceCells(f1,0);
436 
437  Gcoeff=gMat.getCoeff(cell2, cell22);
438 
439  for(int dir=0;dir<_numDir;dir++)
440  {
441  const TArray& f = *_fArrays[dir];
442  NeibGrads[dir].accumulate(Gcoeff,f[cell22]-f[cell2]);
443  if(min2[dir]>f[cell22])min2[dir]=f[cell22];
444  if(max2[dir]<f[cell22])max2[dir]=f[cell22];
445  }
446  }
447 
448 
449  for(int nj=0;nj<neibcount1;nj++)
450  {
451  const int f1=_cellFaces(cell2,nj);
452 
453  SuperbeeLimiter lf;
454  for(int dir=0;dir<_numDir;dir++)
455  {
456  const TArray& f = *_fArrays[dir];
457  GradType& neibGrad=NeibGrads[dir];
458 
459  VectorT3 fVec=faceCoords[f1]-cellCoords[cell2];
460  T SOU=(neibGrad[0]*fVec[0]+neibGrad[1]*fVec[1]+neibGrad[2]*fVec[2]);
461 
462  computeLimitCoeff(limitCoeff2[dir], f[cell2], SOU, min2[dir], max2[dir], lf);
463  }
464  }
465 
466  int count=1;
467  for(int dir=0;dir<_numDir;dir++)
468  {
469  const TArray& f = *_fArrays[dir];
470  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
471  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
472 
473  GradType& grad=Grads[count-1];
474  GradType& neibGrad=NeibGrads[count-1];
475 
476  if (_BCfArray[face]==7)
477  {
478  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[_solidFaces]);
479  const T uwall = v[1][0];
480  const T vwall = v[1][1];
481  const T wwall = v[1][2];
482  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
483  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
484  {
485  Amat.getElement(count,count)-=flux;
486  BVec[count-1]-=flux*f[cell];
487  }
488  else
489  {
490  const int ibFace = ibFaceIndex[face];
491  if (ibFace < 0)
492  throw CException("invalid ib face index");
493  Field& fnd = *_dsfPtr.dsf[dir];
494  const TArray& fIB = dynamic_cast<const TArray&>(fnd[ibFaces]);
495  BVec[count-1]-=flux*fIB[ibFace];
496  }
497  }
498  else if((_BCfArray[face]==0)||(_BCfArray[face]==-1))
499  {
500  if(c_dot_en>T_Scalar(0))
501  {
502  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
503  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
504  T r=gMat.computeR(grad,f,rVec,cell,cell2);
505  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
506  Amat.getElement(count,count)-=flux;
507  if(_conOrder==1)
508  BVec[count-1]-=flux*f[cell];
509  else
510  {
511  BVec[count-1]-=(flux*f[cell]+flux*SOU*limitCoeff1[dir]);
512  }
513  }
514  else
515  {
516  VectorT3 fVec=faceCoords[face]-cellCoords[cell2];
517  VectorT3 rVec=cellCoords[cell]-cellCoords[cell2];
518  T r=gMat.computeR(neibGrad,f,rVec,cell2,cell);
519  T SOU=(neibGrad[0]*fVec[0]+neibGrad[1]*fVec[1]+neibGrad[2]*fVec[2]);
520  if(_conOrder==1)
521  BVec[count-1]-=flux*f[cell2];
522  else
523  {
524  BVec[count-1]-=(flux*f[cell2]+flux*SOU*limitCoeff2[dir]);
525  }
526  }
527  }
528  else if(_BCfArray[face]==5)
529  {
530  if(c_dot_en>T_Scalar(0))
531  {
532  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
533  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
534  T r=gMat.computeR(grad,f,rVec,cell,cell2);
535  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
536  Amat.getElement(count,count)-=flux;
537  if(_conOrder==1)
538  BVec[count-1]-=flux*f[cell];
539  else
540  {
541  BVec[count-1]-=(flux*f[cell]+flux*SOU*limitCoeff1[dir]);
542  }
543  }
544  else
545  BVec[count-1]-=flux*f[cell2];
546  }
547  count++;
548  }
549  }
550  }
551 
552  void COMETConvection(const int cell, TArrow& Amat, TArray& BVec, const int cellcount)
553  {
554  const int neibcount=_cellFaces.getCount(cell);
555  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
556  const StorageSite& ibFaces = _mesh.getIBFaces();
557  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[_faces]);
558  for(int j=0;j<neibcount;j++)
559  {
560  const int face=_cellFaces(cell,j);
561  int cell2=_faceCells(face,1);
562  VectorT3 Af=_faceArea[face];
563  VectorT3 en = _faceArea[face]/_faceAreaMag[face];
564 
565  T flux;
566 
567  if(cell2==cell)
568  {
569  Af=Af*(-1.);
570  en=en*(-1.);
571  cell2=_faceCells(face,0);
572  }
573 
574  int count=1;
575  for(int dir=0;dir<_numDir;dir++)
576  {
577  const TArray& f = *_fArrays[dir];
578  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
579  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
580 
581  if (((ibType[cell] == Mesh::IBTYPE_FLUID) && (ibType[cell2] == Mesh::IBTYPE_BOUNDARY)) ||
582  ((ibType[cell2] == Mesh::IBTYPE_FLUID) && (ibType[cell] == Mesh::IBTYPE_BOUNDARY)))
583  {
584  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[_solidFaces]);
585  const T uwall = v[1][0];
586  const T vwall = v[1][1];
587  const T wwall = v[1][2];
588  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
589  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
590  {
591  Amat.getElement(count,count)-=flux;
592  BVec[count-1]-=flux*f[cell];
593  }
594  else
595  {
596  const int ibFace = ibFaceIndex[face];
597  if (ibFace < 0)
598  throw CException("invalid ib face index");
599  Field& fnd = *_dsfPtr.dsf[dir];
600  const TArray& fIB = dynamic_cast<const TArray&>(fnd[ibFaces]);
601  BVec[count-1]-=flux*fIB[ibFace];
602  }
603  }
604  else
605  {
606  if(c_dot_en>T_Scalar(0))
607  {
608  Amat.getElement(count,count)-=flux;
609  BVec[count-1]-=flux*f[cell];
610  }
611  else
612  BVec[count-1]-=flux*f[cell2];
613  }
614  count++;
615  }
616  }
617  }
618 
619  void COMETConvectionFine(const int cell, TArrow& Amat, TArray& BVec, const GradMatrix& gMat)
620  {
621 
622  const int neibcount=_cellFaces.getCount(cell);
623  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
624  const StorageSite& ibFaces = _mesh.getIBFaces();
625  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[_faces]);
626  const T one(1.0);
627  GradArray Grads(_numDir);
628  Grads.zero();
629 
630  VectorT3 Gcoeff;
631  TArray limitCoeff1(_numDir);
632  TArray min1(_numDir);
633  TArray max1(_numDir);
634  TArray limitCoeff2(_numDir);
635  TArray min2(_numDir);
636  TArray max2(_numDir);
637  for(int dir=0;dir<_numDir;dir++)
638  {
639  limitCoeff1[dir]=T(1.e20);
640  const TArray& f = *_fArrays[dir];
641  min1[dir]=f[cell];
642  max1[dir]=f[cell];
643  }
644 
645  const VectorT3Array& faceCoords=
646  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_faces]);
647  const VectorT3Array& cellCoords=
648  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
649 
650  for(int j=0;j<neibcount;j++)
651  {
652  const int face=_cellFaces(cell,j);
653  int cell2=_faceCells(face,1);
654  if(cell2==cell)
655  cell2=_faceCells(face,0);
656 
657  Gcoeff=gMat.getCoeff(cell, cell2);
658 
659  for(int dir=0;dir<_numDir;dir++)
660  {
661  const TArray& f = *_fArrays[dir];
662  Grads[dir].accumulate(Gcoeff,f[cell2]-f[cell]);
663  if(min1[dir]>f[cell2])min1[dir]=f[cell2];
664  if(max1[dir]<f[cell2])max1[dir]=f[cell2];
665  }
666  }
667 
668  for(int j=0;j<neibcount;j++)
669  {
670  const int face=_cellFaces(cell,j);
671 
672  SuperbeeLimiter lf;
673  for(int dir=0;dir<_numDir;dir++)
674  {
675  const TArray& f = *_fArrays[dir];
676  GradType& grad=Grads[dir];
677 
678  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
679  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
680 
681  computeLimitCoeff(limitCoeff1[dir], f[cell], SOU, min1[dir], max1[dir], lf);
682  }
683  }
684 
685  for(int j=0;j<neibcount;j++)
686  {
687  const int face=_cellFaces(cell,j);
688  int cell2=_faceCells(face,1);
689  VectorT3 Af=_faceArea[face];
690  VectorT3 en = _faceArea[face]/_faceAreaMag[face];
691 
692  T flux;
693 
694  if(cell2==cell)
695  {
696  Af=Af*(-1.);
697  en=en*(-1.);
698  cell2=_faceCells(face,0);
699  }
700 
701  for(int dir=0;dir<_numDir;dir++)
702  {
703  limitCoeff2[dir]=T(1.e20);
704  const TArray& f = *_fArrays[dir];
705  min2[dir]=f[cell2];
706  max2[dir]=f[cell2];
707  }
708 
709  GradArray NeibGrads(_numDir);
710  NeibGrads.zero();
711 
712  const int neibcount1=_cellFaces.getCount(cell2);
713  for(int nj=0;nj<neibcount1;nj++)
714  {
715  const int f1=_cellFaces(cell2,nj);
716  int cell22=_faceCells(f1,1);
717  if(cell2==cell22)
718  cell22=_faceCells(f1,0);
719 
720  Gcoeff=gMat.getCoeff(cell2, cell22);
721 
722  for(int dir=0;dir<_numDir;dir++)
723  {
724  const TArray& f = *_fArrays[dir];
725  NeibGrads[dir].accumulate(Gcoeff,f[cell22]-f[cell2]);
726  if(min2[dir]>f[cell22])min2[dir]=f[cell22];
727  if(max2[dir]<f[cell22])max2[dir]=f[cell22];
728  }
729  }
730 
731 
732  for(int nj=0;nj<neibcount1;nj++)
733  {
734  const int f1=_cellFaces(cell2,nj);
735 
736  SuperbeeLimiter lf;
737  for(int dir=0;dir<_numDir;dir++)
738  {
739  const TArray& f = *_fArrays[dir];
740  GradType& neibGrad=NeibGrads[dir];
741 
742  VectorT3 fVec=faceCoords[f1]-cellCoords[cell2];
743  T SOU=(neibGrad[0]*fVec[0]+neibGrad[1]*fVec[1]+neibGrad[2]*fVec[2]);
744 
745  computeLimitCoeff(limitCoeff2[dir], f[cell2], SOU, min2[dir], max2[dir], lf);
746  }
747  }
748 
749 
750 
751  if(_BCfArray[face]==2) //If the face in question is a reflecting wall
752  {
753  int Fgid=findFgId(face);
754  T uwall = (*(_bcMap[Fgid]))["specifiedXVelocity"];
755  T vwall = (*(_bcMap[Fgid]))["specifiedYVelocity"];
756  T wwall = (*(_bcMap[Fgid]))["specifiedZVelocity"];
757  T Twall = (*(_bcMap[Fgid]))["specifiedTemperature"];
758  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
759  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
760  const vector<int>& vecReflection=(*pos).second;
761  T alpha=(*(_bcMap[Fgid]))["accommodationCoefficient"];
762  T m1alpha = one-alpha;
763  const T pi(acos(-1.0));
764 
765  //first sweep - have to calculate wall number density
766  T Nmr(0.0);
767  T Dmr(0.0);
768  int count=1;
769  for(int dir1=0;dir1<_numDir;dir1++)
770  {
771  const TArray& f = *_fArrays[dir1];
772  const T fwall = 1.0/pow(pi*Twall,1.5)*exp(-(pow(_cx[dir1]-uwall,2.0)+pow(_cy[dir1]-vwall,2.0)+pow(_cz[dir1]-wwall,2.0))/Twall);
773  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
774  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
775  GradType& grad=Grads[count-1];
776  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
777  {
778  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
779  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
780  T r=gMat.computeR(grad,f,rVec,cell,cell2);
781  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
782  Amat.getElement(count,count)-=flux;
783  if(_conOrder==1)
784  {
785  BVec[count-1]-=flux*f[cell];
786  Nmr = Nmr + f[cell]*_wts[dir1]*(c_dot_en -wallV_dot_en);
787  }
788  else
789  {
790  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir1]*SOU);
791  Nmr = Nmr + (f[cell]+limitCoeff1[dir1]*SOU)*_wts[dir1]*(c_dot_en -wallV_dot_en);
792  }
793  }
794  else
795  { //have to move through all other directions
796  Dmr = Dmr - fwall*_wts[dir1]*(c_dot_en-wallV_dot_en);
797  }
798  count++;
799  }
800 
801  const T nwall = Nmr/Dmr; // wall number density for initializing Maxwellian
802 
803  //Second sweep
804  const T zero(0.0);
805  count=1;
806  for(int dir1=0;dir1<_numDir;dir1++)
807  {
808  const TArray& f = *_fArrays[dir1];
809  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
810  const T c1_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
811  if((c1_dot_en-wallV_dot_en)<T_Scalar(0))
812  {
813  const T coeff1 = 1.0/pow(pi*Twall,1.5)*exp(-(pow(_cx[dir1]-uwall,2.0)+pow(_cy[dir1]-vwall,2.0)+pow(_cz[dir1]-wwall,2.0))/Twall);
814  if(m1alpha!=zero)
815  {
816  const int direction_incident = vecReflection[dir1];
817  const TArray& dsfi = *_fArrays[direction_incident];
818  GradType& grad=Grads[direction_incident];
819  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
820  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
821  //Amat.getElement(count,direction_incident+1)-=flux*m1alpha;
822  if(_conOrder==1)
823  BVec[count-1]-=flux*m1alpha*dsfi[cell];
824  else
825  BVec[count-1]-=flux*m1alpha*(dsfi[cell]+limitCoeff1[direction_incident]*SOU);
826  }
827  BVec[count-1]-=flux*alpha*(nwall*coeff1);
828  }
829  count++;
830  }
831  }
832  /*
833  else if(_BCfArray[face]==3) //If the face in question is a inlet velocity face
834  {
835  int Fgid=findFgId(face);
836  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
837  const vector<int>& vecReflection=(*pos).second;
838  const T pi(acos(-1.0));
839 
840  //first sweep - have to calculate Dmr
841  const T uin = (*(_bcMap[Fgid]))["specifiedXVelocity"];
842  const T vin = (*(_bcMap[Fgid]))["specifiedYVelocity"];
843  const T win = (*(_bcMap[Fgid]))["specifiedZVelocity"];
844  const T Tin = (*(_bcMap[Fgid]))["specifiedTemperature"];
845  const T mdot = (*(_bcMap[Fgid]))["specifiedMassFlowRate"];
846  T Nmr(0.0);
847  T Dmr(0.0);
848  const T R=8314.0/_MW;
849  const T u_init=pow(2.0*R*_T_init,0.5);
850  int count=1;
851  for(int dir1=0;dir1<_numDir;dir1++)
852  {
853  Field& fnd = *_dsfPtr.dsf[dir1];
854  const TArray& f = dynamic_cast<const TArray&>(fnd[_cells]);
855  const T fwall = 1.0/pow(pi*Tin,1.5)*exp(-(pow(_cx[dir1]-uin,2.0)+pow(_cy[dir1]-vin,2.0)+pow(_cz[dir1]-win,2.0))/Tin);
856  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
857  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
858  if(c_dot_en>T_Scalar(0))
859  {
860  Amat(count,count)-=flux;
861  BVec[count-1]-=flux*f[cell];
862  }
863  else
864  { //have to move through all other directions
865  Dmr = Dmr + fwall*_wts[dir1]*c_dot_en;
866  }
867  count++;
868  }
869 
870  Nmr=mdot/(_rho_init*u_init);
871  const T nin = Nmr/Dmr; // wall number density for initializing Maxwellian
872 
873  //Second sweep
874  count=1;
875  for(int dir1=0;dir1<_numDir;dir1++)
876  {
877  Field& fnd = *_dsfPtr.dsf[dir1];
878  const TArray& f = dynamic_cast<const TArray&>(fnd[_cells]);
879  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
880  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
881  if(c_dot_en<T_Scalar(0))
882  {
883  const int direction_incident = vecReflection[dir1];
884  Field& fndi = *_dsfPtr.dsf[direction_incident];
885  const TArray& dsfi = dynamic_cast<const TArray&>(fndi[_cells]);
886  Amat(count,direction_incident+1)-=flux;
887  BVec[count-1]-=flux*(nin/pow(pi*Tin,1.5)*exp(-(pow(_cx[dir1]-uin,2.0)+pow(_cy[dir1]-vin,2.0)+pow(_cz[dir1]-win,2.0))/Tin)+dsfi[cell]);
888  }
889  count++;
890  }
891  }
892  */
893  else if(_BCfArray[face]==4) //if the face in question is zero derivative
894  {
895  int count=1;
896  for(int dir=0;dir<_numDir;dir++)
897  {
898  const TArray& f = *_fArrays[dir];
899  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
900  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
901  GradType& grad=Grads[count-1];
902 
903  if(c_dot_en>T_Scalar(0))
904  {
905  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
906  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
907  T r=gMat.computeR(grad,f,rVec,cell,cell2);
908  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
909  Amat.getElement(count,count)-=flux;
910  if(_conOrder==1)
911  BVec[count-1]-=flux*f[cell];
912  else
913  {
914  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir]*SOU);
915  }
916  }
917  else
918  {
919  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
920  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
921  T r=gMat.computeR(grad,f,rVec,cell,cell2);
922  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
923  Amat.getElement(count,count)-=flux;
924  if(_conOrder==1)
925  BVec[count-1]-=flux*f[cell];
926  else
927  {
928  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir]*SOU);
929  }
930  }
931  count++;
932  }
933  }
934  else if(_BCfArray[face]==5) //if the face in question is specified pressure
935  {
936  int count=1;
937  for(int dir=0;dir<_numDir;dir++)
938  {
939  const TArray& f = *_fArrays[dir];
940  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
941  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
942  GradType& grad=Grads[count-1];
943 
944  if(c_dot_en>T_Scalar(0))
945  {
946  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
947  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
948  T r=gMat.computeR(grad,f,rVec,cell,cell2);
949  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
950  Amat.getElement(count,count)-=flux;
951  if(_conOrder==1)
952  BVec[count-1]-=flux*f[cell];
953  else
954  {
955  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir]*SOU);
956  }
957  }
958  else
959  BVec[count-1]-=flux*f[cell2];
960  count++;
961  }
962  }
963  else if(_BCfArray[face]==6) //If the face in question is a symmetry wall
964  {
965  int Fgid=findFgId(face);
966  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
967  const vector<int>& vecReflection=(*pos).second;
968 
969  int count=1;
970  for(int dir1=0;dir1<_numDir;dir1++)
971  {
972  const TArray& f = *_fArrays[dir1];
973  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
974  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
975  GradType& grad=Grads[count-1];
976  if(c_dot_en>T_Scalar(0))
977  {
978  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
979  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
980  T r=gMat.computeR(grad,f,rVec,cell,cell2);
981  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
982  Amat.getElement(count,count)-=flux;
983  if(_conOrder==1)
984  BVec[count-1]-=flux*f[cell];
985  else
986  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir1]*SOU);
987  }
988  else
989  {
990  const int direction_incident = vecReflection[dir1];
991  const TArray& dsfi = *_fArrays[direction_incident];
992  BVec[count-1]-=flux*dsfi[cell];
993  }
994  count++;
995  }
996  }
997  else if(_BCfArray[face]==7) //if the face in question is an ibFace
998  {
999  int count=1;
1000  for(int dir=0;dir<_numDir;dir++)
1001  {
1002  const TArray& f = *_fArrays[dir];
1003  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1004  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1005 
1006  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[_solidFaces]);
1007  const T uwall = v[1][0];
1008  const T vwall = v[1][1];
1009  const T wwall = v[1][2];
1010  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
1011  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
1012  {
1013  Amat.getElement(count,count)-=flux;
1014  BVec[count-1]-=flux*f[cell];
1015  }
1016  else
1017  {
1018  const int ibFace = ibFaceIndex[face];
1019  if (ibFace < 0)
1020  throw CException("invalid ib face index");
1021  Field& fnd = *_dsfPtr.dsf[dir];
1022  const TArray& fIB = dynamic_cast<const TArray&>(fnd[ibFaces]);
1023  BVec[count-1]-=flux*fIB[ibFace];
1024  }
1025  count++;
1026  }
1027  }
1028  else if(_BCfArray[face]==0) //if the face in question is not reflecting
1029  {
1030  int count=1;
1031  for(int dir=0;dir<_numDir;dir++)
1032  {
1033  const TArray& f = *_fArrays[dir];
1034  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1035  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1036 
1037  GradType& grad=Grads[count-1];
1038  GradType& neibGrad=NeibGrads[count-1];
1039  if(c_dot_en>T_Scalar(0))
1040  {
1041  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
1042  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
1043  T r=gMat.computeR(grad,f,rVec,cell,cell2);
1044  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
1045  Amat.getElement(count,count)-=flux;
1046  if(_conOrder==1)
1047  BVec[count-1]-=flux*f[cell];
1048  else
1049  {
1050  BVec[count-1]-=flux*(f[cell]+limitCoeff1[dir]*SOU);
1051  }
1052  }
1053  else
1054  {
1055  VectorT3 fVec=faceCoords[face]-cellCoords[cell2];
1056  VectorT3 rVec=cellCoords[cell]-cellCoords[cell2];
1057  T r=gMat.computeR(neibGrad,f,rVec,cell2,cell);
1058  T SOU=(neibGrad[0]*fVec[0]+neibGrad[1]*fVec[1]+neibGrad[2]*fVec[2]);
1059  if(_conOrder==1)
1060  BVec[count-1]-=flux*f[cell2];
1061  else
1062  {
1063  BVec[count-1]-=flux*(f[cell2]+limitCoeff2[dir]*SOU);
1064  }
1065  }
1066  count++;
1067  }
1068  }
1069  else if(_BCfArray[face]==-1) //if the face in question is interface
1070  {
1071  int count=1;
1072  for(int dir=0;dir<_numDir;dir++)
1073  {
1074  const TArray& f = *_fArrays[dir];
1075  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1076  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1077 
1078  if(c_dot_en>T_Scalar(0))
1079  {
1080  Amat.getElement(count,count)-=flux;
1081  BVec[count-1]-=flux*f[cell];
1082  }
1083  else
1084  BVec[count-1]-=flux*f[cell2];
1085  count++;
1086  }
1087  }
1088  }
1089  }
1090 
1091  void COMETConvection(const int cell, TArrow& Amat, TArray& BVec)
1092  {
1093 
1094  const int neibcount=_cellFaces.getCount(cell);
1095  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
1096  const StorageSite& ibFaces = _mesh.getIBFaces();
1097  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[_faces]);
1098  const T one(1.0);
1099 
1100  for(int j=0;j<neibcount;j++)
1101  {
1102  const int face=_cellFaces(cell,j);
1103  int cell2=_faceCells(face,1);
1104  VectorT3 Af=_faceArea[face];
1105  VectorT3 en = _faceArea[face]/_faceAreaMag[face];
1106 
1107  T flux;
1108 
1109  if(cell2==cell)
1110  {
1111  Af=Af*(-1.);
1112  en=en*(-1.);
1113  cell2=_faceCells(face,0);
1114  }
1115 
1116  if(_BCfArray[face]==2) //If the face in question is a reflecting wall
1117  {
1118  int Fgid=findFgId(face);
1119  T uwall = (*(_bcMap[Fgid]))["specifiedXVelocity"];
1120  T vwall = (*(_bcMap[Fgid]))["specifiedYVelocity"];
1121  T wwall = (*(_bcMap[Fgid]))["specifiedZVelocity"];
1122  T Twall = (*(_bcMap[Fgid]))["specifiedTemperature"];
1123  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
1124  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
1125  const vector<int>& vecReflection=(*pos).second;
1126  T alpha=(*(_bcMap[Fgid]))["accommodationCoefficient"];
1127  T m1alpha = one-alpha;
1128  const T pi(acos(-1.0));
1129 
1130  //first sweep - have to calculate wall number density
1131  T Nmr(0.0);
1132  T Dmr(0.0);
1133  int count=1;
1134  for(int dir1=0;dir1<_numDir;dir1++)
1135  {
1136  const TArray& f = *_fArrays[dir1];
1137  const T fwall = 1.0/pow(pi*Twall,1.5)*exp(-(pow(_cx[dir1]-uwall,2.0)+pow(_cy[dir1]-vwall,2.0)+pow(_cz[dir1]-wwall,2.0))/Twall);
1138  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
1139  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
1140  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
1141  {
1142  Amat.getElement(count,count)-=flux;
1143  BVec[count-1]-=flux*f[cell];
1144  Nmr = Nmr + f[cell]*_wts[dir1]*(c_dot_en -wallV_dot_en);
1145  }
1146  else
1147  { //have to move through all other directions
1148  Dmr = Dmr - fwall*_wts[dir1]*(c_dot_en-wallV_dot_en);
1149  }
1150  count++;
1151  }
1152 
1153  const T nwall = Nmr/Dmr; // wall number density for initializing Maxwellian
1154 
1155  //Second sweep
1156  const T zero(0.0);
1157  count=1;
1158  for(int dir1=0;dir1<_numDir;dir1++)
1159  {
1160  const TArray& f = *_fArrays[dir1];
1161  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
1162  const T c1_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
1163  if((c1_dot_en-wallV_dot_en)<T_Scalar(0))
1164  {
1165  const T coeff1 = 1.0/pow(pi*Twall,1.5)*exp(-(pow(_cx[dir1]-uwall,2.0)+pow(_cy[dir1]-vwall,2.0)+pow(_cz[dir1]-wwall,2.0))/Twall);
1166  if(m1alpha!=zero)
1167  {
1168  const int direction_incident = vecReflection[dir1];
1169  const TArray& dsfi = *_fArrays[direction_incident];
1170  //Amat.getElement(count,direction_incident+1)-=flux*m1alpha;
1171  BVec[count-1]-=flux*m1alpha*dsfi[cell];
1172  }
1173  BVec[count-1]-=flux*alpha*(nwall*coeff1);
1174  }
1175  count++;
1176  }
1177  }
1178  /*
1179  else if(_BCfArray[face]==3) //If the face in question is a inlet velocity face
1180  {
1181  int Fgid=findFgId(face);
1182  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
1183  const vector<int>& vecReflection=(*pos).second;
1184  const T pi(acos(-1.0));
1185 
1186  //first sweep - have to calculate Dmr
1187  const T uin = (*(_bcMap[Fgid]))["specifiedXVelocity"];
1188  const T vin = (*(_bcMap[Fgid]))["specifiedYVelocity"];
1189  const T win = (*(_bcMap[Fgid]))["specifiedZVelocity"];
1190  const T Tin = (*(_bcMap[Fgid]))["specifiedTemperature"];
1191  const T mdot = (*(_bcMap[Fgid]))["specifiedMassFlowRate"];
1192  T Nmr(0.0);
1193  T Dmr(0.0);
1194  const T R=8314.0/_MW;
1195  const T u_init=pow(2.0*R*_T_init,0.5);
1196  int count=1;
1197  for(int dir1=0;dir1<_numDir;dir1++)
1198  {
1199  Field& fnd = *_dsfPtr.dsf[dir1];
1200  const TArray& f = dynamic_cast<const TArray&>(fnd[_cells]);
1201  const T fwall = 1.0/pow(pi*Tin,1.5)*exp(-(pow(_cx[dir1]-uin,2.0)+pow(_cy[dir1]-vin,2.0)+pow(_cz[dir1]-win,2.0))/Tin);
1202  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
1203  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
1204  if(c_dot_en>T_Scalar(0))
1205  {
1206  Amat(count,count)-=flux;
1207  BVec[count-1]-=flux*f[cell];
1208  }
1209  else
1210  { //have to move through all other directions
1211  Dmr = Dmr + fwall*_wts[dir1]*c_dot_en;
1212  }
1213  count++;
1214  }
1215 
1216  Nmr=mdot/(_rho_init*u_init);
1217  const T nin = Nmr/Dmr; // wall number density for initializing Maxwellian
1218 
1219  //Second sweep
1220  count=1;
1221  for(int dir1=0;dir1<_numDir;dir1++)
1222  {
1223  Field& fnd = *_dsfPtr.dsf[dir1];
1224  const TArray& f = dynamic_cast<const TArray&>(fnd[_cells]);
1225  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
1226  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
1227  if(c_dot_en<T_Scalar(0))
1228  {
1229  const int direction_incident = vecReflection[dir1];
1230  Field& fndi = *_dsfPtr.dsf[direction_incident];
1231  const TArray& dsfi = dynamic_cast<const TArray&>(fndi[_cells]);
1232  Amat(count,direction_incident+1)-=flux;
1233  BVec[count-1]-=flux*(nin/pow(pi*Tin,1.5)*exp(-(pow(_cx[dir1]-uin,2.0)+pow(_cy[dir1]-vin,2.0)+pow(_cz[dir1]-win,2.0))/Tin)+dsfi[cell]);
1234  }
1235  count++;
1236  }
1237  }
1238  */
1239  else if(_BCfArray[face]==4) //if the face in question is zero derivative
1240  {
1241  int count=1;
1242  for(int dir=0;dir<_numDir;dir++)
1243  {
1244  const TArray& f = *_fArrays[dir];
1245  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1246  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1247 
1248  if(c_dot_en>T_Scalar(0))
1249  {
1250  Amat.getElement(count,count)-=flux;
1251  BVec[count-1]-=flux*f[cell];
1252  }
1253  else
1254  {
1255  Amat.getElement(count,count)-=flux;
1256  BVec[count-1]-=flux*f[cell];
1257  }
1258  count++;
1259  }
1260  }
1261  else if(_BCfArray[face]==5) //if the face in question is specified pressure
1262  {
1263  int count=1;
1264  for(int dir=0;dir<_numDir;dir++)
1265  {
1266  const TArray& f = *_fArrays[dir];
1267  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1268  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1269 
1270  if(c_dot_en>T_Scalar(0))
1271  {
1272  Amat.getElement(count,count)-=flux;
1273  BVec[count-1]-=flux*f[cell];
1274  }
1275  else
1276  BVec[count-1]-=flux*f[cell2];
1277  count++;
1278  }
1279  }
1280  else if(_BCfArray[face]==6) //If the face in question is a symmetry wall
1281  {
1282  int Fgid=findFgId(face);
1283  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(Fgid);
1284  const vector<int>& vecReflection=(*pos).second;
1285 
1286  int count=1;
1287  for(int dir1=0;dir1<_numDir;dir1++)
1288  {
1289  const TArray& f = *_fArrays[dir1];
1290  flux=_cx[dir1]*Af[0]+_cy[dir1]*Af[1]+_cz[dir1]*Af[2];
1291  const T c_dot_en = _cx[dir1]*en[0]+_cy[dir1]*en[1]+_cz[dir1]*en[2];
1292  if(c_dot_en>T_Scalar(0))
1293  {
1294  Amat.getElement(count,count)-=flux;
1295  BVec[count-1]-=flux*f[cell];
1296  }
1297  else
1298  {
1299  const int direction_incident = vecReflection[dir1];
1300  const TArray& dsfi = *_fArrays[direction_incident];
1301  BVec[count-1]-=flux*dsfi[cell];
1302  }
1303  count++;
1304  }
1305  }
1306  else if(_BCfArray[face]==7) //if the face in question is an ibFace
1307  {
1308  int count=1;
1309  for(int dir=0;dir<_numDir;dir++)
1310  {
1311  const TArray& f = *_fArrays[dir];
1312  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1313  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1314 
1315  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[_solidFaces]);
1316  const T uwall = v[1][0];
1317  const T vwall = v[1][1];
1318  const T wwall = v[1][2];
1319  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
1320  if((c_dot_en-wallV_dot_en)>T_Scalar(0))
1321  {
1322  Amat.getElement(count,count)-=flux;
1323  BVec[count-1]-=flux*f[cell];
1324  }
1325  else
1326  {
1327  const int ibFace = ibFaceIndex[face];
1328  if (ibFace < 0)
1329  throw CException("invalid ib face index");
1330  Field& fnd = *_dsfPtr.dsf[dir];
1331  const TArray& fIB = dynamic_cast<const TArray&>(fnd[ibFaces]);
1332  BVec[count-1]-=flux*fIB[ibFace];
1333  }
1334  count++;
1335  }
1336  }
1337  else if(_BCfArray[face]==0) //if the face in question is not reflecting
1338  {
1339  int count=1;
1340  for(int dir=0;dir<_numDir;dir++)
1341  {
1342  const TArray& f = *_fArrays[dir];
1343  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1344  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1345 
1346  if(c_dot_en>T_Scalar(0))
1347  {
1348  Amat.getElement(count,count)-=flux;
1349  BVec[count-1]-=flux*f[cell];
1350  }
1351  else
1352  BVec[count-1]-=flux*f[cell2];
1353  count++;
1354  }
1355  }
1356  else if(_BCfArray[face]==-1) //if the face in question is interface
1357  {
1358  int count=1;
1359  for(int dir=0;dir<_numDir;dir++)
1360  {
1361  const TArray& f = *_fArrays[dir];
1362  flux=_cx[dir]*Af[0]+_cy[dir]*Af[1]+_cz[dir]*Af[2];
1363  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1364 
1365  if(c_dot_en>T_Scalar(0))
1366  {
1367  Amat.getElement(count,count)-=flux;
1368  BVec[count-1]-=flux*f[cell];
1369  }
1370  else
1371  BVec[count-1]-=flux*f[cell2];
1372  count++;
1373  }
1374  }
1375  }
1376  }
1377 
1378  template<class MatrixType>
1379  void COMETTest(const int cell, MatrixType Amat, TArray& BVec, TArray& fV)
1380  {
1381  const int order=_numDir;
1382 
1383  VectorT3Array& v = _velocity;
1384  TArray& density = _density;
1385 
1386  const T two(2.0);
1387 
1388  T coeff;
1389  int count = 1;
1390 
1391  for(int direction=0;direction<_numDir;direction++)
1392  {
1393  //const TArray& f = *_fArrays[direction];
1394  const TArray& fEqES = *_fEqESArrays[direction];
1395  coeff =_cellVolume[cell]*_collisionFrequency[cell];
1396 
1397  T C1=(_cx[direction]-v[cell][0]);
1398  T C2=(_cy[direction]-v[cell][1]);
1399  T C3=(_cz[direction]-v[cell][2]);
1400 
1401  Amat->getElement(count,order+1)+=coeff*fEqES[cell]*(two*_coeffg[cell][1]*C1-_coeffg[cell][2]);
1402  Amat->getElement(count,order+2)+=coeff*fEqES[cell]*(two*_coeffg[cell][3]*C2-_coeffg[cell][4]);
1403  Amat->getElement(count,order+3)+=coeff*fEqES[cell]*(two*_coeffg[cell][5]*C3-_coeffg[cell][6]);
1404  Amat->getElement(count,count)-=coeff;
1405 
1406  BVec[count-1]+=coeff*(fEqES[cell]-fV[direction]);
1407 
1408  Amat->getElement(_numDir+1,count)+=_wts[direction]*C1/density[cell];
1409  BVec[_numDir]+=_cx[direction]*_wts[direction]*fV[direction]/density[cell];
1410  Amat->getElement(_numDir+2,count)+=_wts[direction]*C2/density[cell];
1411  BVec[_numDir+1]+=_cy[direction]*_wts[direction]*fV[direction]/density[cell];
1412  Amat->getElement(_numDir+3,count)+=_wts[direction]*C3/density[cell];
1413  BVec[_numDir+2]+=_cz[direction]*_wts[direction]*fV[direction]/density[cell];
1414 
1415  count++;
1416  }
1417  Amat->getElement(_numDir+1,_numDir+1)-=1;
1418  BVec[_numDir]-=v[cell][0];
1419  Amat->getElement(_numDir+2,_numDir+2)-=1;
1420  BVec[_numDir+1]-=v[cell][1];
1421  Amat->getElement(_numDir+3,_numDir+3)-=1;
1422  BVec[_numDir+2]-=v[cell][2];
1423  }
1424 
1425  template<class MatrixType>
1426  void COMETCollision(const int cell, MatrixType Amat, TArray& BVec)
1427  {
1428  const int order=_numDir;
1429 
1430  VectorT3Array& v = _velocity;
1431 
1432  const T two(2.0);
1433 
1434  T coeff;
1435  int count = 1;
1436 
1437  for(int direction=0;direction<_numDir;direction++)
1438  {
1439  const TArray& f = *_fArrays[direction];
1440  coeff =_cellVolume[cell]*_collisionFrequency[cell];
1441 
1442  T C1=(_cx[direction]-v[cell][0]);
1443  T C2=(_cy[direction]-v[cell][1]);
1444  T C3=(_cz[direction]-v[cell][2]);
1445  T fGamma=_coeffg[cell][0]*exp(-_coeffg[cell][1]*SQR(C1)+_coeffg[cell][2]*C1
1446  -_coeffg[cell][3]*SQR(C2)+_coeffg[cell][4]*C2
1447  -_coeffg[cell][5]*SQR(C3)+_coeffg[cell][6]*C3
1448  +_coeffg[cell][7]*_cx[direction]*_cy[direction]
1449  +_coeffg[cell][8]*_cy[direction]*_cz[direction]
1450  +_coeffg[cell][9]*_cz[direction]*_cx[direction]);
1451 
1452  Amat->getElement(count,order+1)+=coeff*fGamma*(two*_coeffg[cell][1]*C1-_coeffg[cell][2]);
1453  Amat->getElement(count,order+2)+=coeff*fGamma*(two*_coeffg[cell][3]*C2-_coeffg[cell][4]);
1454  Amat->getElement(count,order+3)+=coeff*fGamma*(two*_coeffg[cell][5]*C3-_coeffg[cell][6]);
1455  Amat->getElement(count,count)-=coeff;
1456 
1457  BVec[count-1]+=coeff*fGamma;
1458  BVec[count-1]-=coeff*f[cell];
1459  count++;
1460  }
1461  }
1462 
1463  template<class MatrixType>
1464  void COMETMacro(const int cell, MatrixType Amat, TArray& BVec)
1465  {
1466 
1467  VectorT3Array& v = _velocity;
1468 
1469  T density(0.);
1470  for(int dir=0;dir<_numDir;dir++)
1471  {
1472  const TArray& f = *_fArrays[dir];
1473  density+=f[cell]*_wts[dir];
1474  }
1475 
1476  int count = 1;
1477 
1478  for(int dir=0;dir<_numDir;dir++)
1479  {
1480  const TArray& f = *_fArrays[dir];
1481  T C1=(_cx[dir]-v[cell][0]);
1482  T C2=(_cy[dir]-v[cell][1]);
1483  T C3=(_cz[dir]-v[cell][2]);
1484 
1485  Amat->getElement(_numDir+1,count)+=_wts[dir]*C1/density;
1486  BVec[_numDir]+=_cx[dir]*_wts[dir]*f[cell]/density;
1487  Amat->getElement(_numDir+2,count)+=_wts[dir]*C2/density;
1488  BVec[_numDir+1]+=_cy[dir]*_wts[dir]*f[cell]/density;
1489  Amat->getElement(_numDir+3,count)+=_wts[dir]*C3/density;
1490  BVec[_numDir+2]+=_cz[dir]*_wts[dir]*f[cell]/density;
1491  count++;
1492  }
1493  Amat->getElement(_numDir+1,_numDir+1)-=1;
1494  BVec[_numDir]-=v[cell][0];
1495  Amat->getElement(_numDir+2,_numDir+2)-=1;
1496  BVec[_numDir+1]-=v[cell][1];
1497  Amat->getElement(_numDir+3,_numDir+3)-=1;
1498  BVec[_numDir+2]-=v[cell][2];
1499  }
1500 
1501  void Distribute(const int cell, TArray& BVec, TArray& Rvec)
1502  {
1503  VectorT3Array& v = _velocity;
1504  VectorT3Array& vR = _velocityResidual;
1505 
1506  for(int direction=0;direction<_numDir;direction++)
1507  {
1508  TArray& f = *_fArrays[direction];
1509  TArray& fRes = *_fResArrays[direction];
1510  f[cell]-=_underRelaxation*BVec[direction];
1511  fRes[cell]=-Rvec[direction];
1512  }
1513  v[cell][0]-=_underRelaxation*BVec[_numDir];
1514  v[cell][1]-=_underRelaxation*BVec[_numDir+1];
1515  v[cell][2]-=_underRelaxation*BVec[_numDir+2];
1516  vR[cell][0]=-Rvec[_numDir];
1517  vR[cell][1]=-Rvec[_numDir+1];
1518  vR[cell][2]=-Rvec[_numDir+2];
1519  }
1520 
1521  void Distribute(const int cell, TArray& Rvec)
1522  {
1523  VectorT3Array& vR = _velocityResidual;
1524 
1525  for(int direction=0;direction<_numDir;direction++)
1526  {
1527  TArray& fRes = *_fResArrays[direction];
1528  fRes[cell]=-Rvec[direction];
1529  }
1530  vR[cell][0]=-Rvec[_numDir];
1531  vR[cell][1]=-Rvec[_numDir+1];
1532  vR[cell][2]=-Rvec[_numDir+2];
1533  }
1534 
1535  void ComputeMacroparameters(const int cell)
1536  {
1537 
1538  VectorT3Array& v = _velocity;
1539 
1540  const T zero(0.);
1541 
1542  _density[cell]=zero;
1543  _temperature[cell]=zero;
1544  _stress[cell][0]=0.0;_stress[cell][1]=0.0;_stress[cell][2]=0.0;
1545  _stress[cell][3]=0.0;_stress[cell][4]=0.0;_stress[cell][5]=0.0;
1546 
1547  for(int dir=0;dir<_numDir;dir++)
1548  {
1549  const TArray& f = *_fArrays[dir];
1550  _density[cell] = _density[cell]+_wts[dir]*f[cell];
1551  _temperature[cell]= _temperature[cell]+(SQR(_cx[dir])+SQR(_cy[dir])
1552  +SQR(_cz[dir]))*f[cell]*_wts[dir];
1553  }
1554 
1555  _temperature[cell]=_temperature[cell]-(SQR(v[cell][0])
1556  +SQR(v[cell][1])
1557  +SQR(v[cell][2]))*_density[cell];
1558  _temperature[cell]=_temperature[cell]/(1.5*_density[cell]);
1559  _pressure[cell]=_density[cell]*_temperature[cell];
1560 
1561  for(int dir=0;dir<_numDir;dir++)
1562  {
1563  const TArray& f = *_fArrays[dir];
1564  _stress[cell][0] +=SQR((_cx[dir]-v[cell][0]))*f[cell]*_wts[dir];
1565  _stress[cell][1] +=SQR((_cy[dir]-v[cell][1]))*f[cell]*_wts[dir];
1566  _stress[cell][2] +=SQR((_cz[dir]-v[cell][2]))*f[cell]*_wts[dir];
1567  _stress[cell][3] +=(_cx[dir]-v[cell][0])*(_cy[dir]-v[cell][1])*f[cell]*_wts[dir];
1568  _stress[cell][4] +=(_cy[dir]-v[cell][1])*(_cz[dir]-v[cell][2])*f[cell]*_wts[dir];
1569  _stress[cell][5] +=(_cz[dir]-v[cell][2])*(_cx[dir]-v[cell][0])*f[cell]*_wts[dir];
1570  }
1571  }
1572 
1573  void findResidFine(const bool plusFAS)
1574  {
1575  const int cellcount=_cells.getSelfCount();
1576  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
1577 
1578  TArray ResidSum(_numDir+3);
1579  TArray Bsum(_numDir+3);
1580  TArray temp(_numDir+3+1);
1581  T traceSum=0.;
1582  T ResidScalar=0.;
1583  ResidSum.zero();
1584  Bsum.zero();
1585  temp.zero();
1586 
1587  TArray Bvec(_numDir+3);
1588  TArray Resid(_numDir+3);
1589  TArrow AMat(_numDir+3);
1590 
1591  TArray fVal(_numDir);
1592 
1593  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
1594 
1595  for(int c=0;c<cellcount;c++)
1596  {
1597  if (ibType[c] == Mesh::IBTYPE_FLUID){
1598  for(int dir=0;dir<_numDir;dir++)
1599  fVal[dir]=(*_fArrays[dir])[c];
1600  if(_BCArray[c]==0)
1601  {
1602  Bvec.zero();
1603  Resid.zero();
1604  AMat.zero();
1605 
1606  if(_transient)
1607  COMETUnsteady(c,&AMat,Bvec);
1608 
1609  COMETConvectionFine(c,AMat,Bvec,cellcount,gradMatrix);
1610  COMETTest(c,&AMat,Bvec,fVal);
1611  //COMETCollision(c,&AMat,Bvec);
1612  //COMETMacro(c,&AMat,Bvec);
1613 
1614  if(plusFAS)
1615  addFAS(c,Bvec);
1616 
1617  traceSum+=AMat.getTraceAbs();
1618  Resid=Bvec;
1619  //Bvec.zero();
1620  Distribute(c,Resid);
1621 
1622  makeValueArray(c,Bvec);
1623 
1624  AMat.multiply(Resid,Bvec);
1625  Resid=Bvec;
1626 
1627  //ArrayAbs(Bvec);
1628  //ArrayAbs(Resid);
1629  ArrayAbs(Bvec,Resid);
1630  Bsum+=Bvec;
1631  ResidSum+=Resid;
1632  }
1633  else if(_BCArray[c]==1) //reflecting boundary
1634  {
1635  Bvec.zero();
1636  Resid.zero();
1637  AMat.zero();
1638 
1639  if(_transient)
1640  COMETUnsteady(c,&AMat,Bvec);
1641 
1642  COMETConvectionFine(c,AMat,Bvec,gradMatrix);
1643  //COMETConvection(c,AMat,Bvec);
1644  COMETTest(c,&AMat,Bvec,fVal);
1645  //COMETCollision(c,&AMat,Bvec);
1646  //COMETMacro(c,&AMat,Bvec);
1647 
1648  if(plusFAS)
1649  addFAS(c,Bvec);
1650 
1651  traceSum+=AMat.getTraceAbs();
1652  Resid=Bvec;
1653  //Bvec.zero();
1654  Distribute(c,Resid);
1655 
1656  makeValueArray(c,Bvec);
1657 
1658  AMat.multiply(Resid,Bvec);
1659  Resid=Bvec;
1660 
1661  //ArrayAbs(Bvec);
1662  //ArrayAbs(Resid);
1663  ArrayAbs(Bvec,Resid);
1664  Bsum+=Bvec;
1665  ResidSum+=Resid;
1666  }
1667  else
1668  throw CException("Unexpected value for boundary cell map.");
1669  }
1670  }
1671  for(int o=0;o<_numDir+3;o++)
1672  {
1673  temp[o]=ResidSum[o];
1674  }
1675  temp[_numDir+3]=traceSum;
1676 
1677  //MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, ResidSum.getData(), _numDir+3, MPI::DOUBLE, MPI::SUM);
1678  //MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &traceSum, 1, MPI::DOUBLE, MPI::SUM);
1679 #ifdef FVM_PARALLEL
1680  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, temp.getData(), _numDir+4, MPI::DOUBLE, MPI::SUM);
1681 #endif
1682  /*
1683  for(int o=0;o<_numDir+3;o++)
1684  {
1685  ResidScalar+=ResidSum[o];
1686  }
1687 
1688  ResidScalar/=traceSum;
1689  */
1690 
1691  for(int o=0;o<_numDir+3;o++)
1692  {
1693  ResidScalar+=temp[o];
1694  }
1695 
1696  ResidScalar/=temp[_numDir+3];
1697 
1698  if(_aveResid==-1)
1699  {_aveResid=ResidScalar;}
1700  else
1701  {
1702  _residChange=fabs(_aveResid-ResidScalar)/_aveResid;
1703  _aveResid=ResidScalar;
1704  }
1705  }
1706 
1707  void findResid(const bool plusFAS)
1708  {
1709  const int cellcount=_cells.getSelfCount();
1710  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
1711 
1712  TArray ResidSum(_numDir+3);
1713  TArray Bsum(_numDir+3);
1714  TArray temp(_numDir+3+1);
1715  T traceSum=0.;
1716  T ResidScalar=0.;
1717  ResidSum.zero();
1718  Bsum.zero();
1719  temp.zero();
1720 
1721  TArray Bvec(_numDir+3);
1722  TArray Resid(_numDir+3);
1723  TArrow AMat(_numDir+3);
1724 
1725  TArray fVal(_numDir);
1726 
1727  for(int c=0;c<cellcount;c++)
1728  {
1729  if (ibType[c] == Mesh::IBTYPE_FLUID){
1730  for(int dir=0;dir<_numDir;dir++)
1731  fVal[dir]=(*_fArrays[dir])[c];
1732  if(_BCArray[c]==0)
1733  {
1734  Bvec.zero();
1735  Resid.zero();
1736  AMat.zero();
1737 
1738  if(_transient)
1739  COMETUnsteady(c,&AMat,Bvec);
1740 
1741  COMETConvection(c,AMat,Bvec,cellcount);
1742  COMETTest(c,&AMat,Bvec,fVal);
1743  //COMETCollision(c,&AMat,Bvec);
1744  //COMETMacro(c,&AMat,Bvec);
1745 
1746  if(plusFAS)
1747  addFAS(c,Bvec);
1748 
1749  traceSum+=AMat.getTraceAbs();
1750  Resid=Bvec;
1751  //Bvec.zero();
1752  Distribute(c,Resid);
1753 
1754  makeValueArray(c,Bvec);
1755 
1756  AMat.multiply(Resid,Bvec);
1757  Resid=Bvec;
1758 
1759  //ArrayAbs(Bvec);
1760  //ArrayAbs(Resid);
1761  ArrayAbs(Bvec,Resid);
1762  Bsum+=Bvec;
1763  ResidSum+=Resid;
1764  }
1765  else if(_BCArray[c]==1) //reflecting boundary
1766  {
1767  Bvec.zero();
1768  Resid.zero();
1769  AMat.zero();
1770 
1771  if(_transient)
1772  COMETUnsteady(c,&AMat,Bvec);
1773 
1774  COMETConvection(c,AMat,Bvec);
1775  COMETTest(c,&AMat,Bvec,fVal);
1776  //COMETCollision(c,&AMat,Bvec);
1777  //COMETMacro(c,&AMat,Bvec);
1778 
1779  if(plusFAS)
1780  addFAS(c,Bvec);
1781 
1782  traceSum+=AMat.getTraceAbs();
1783  Resid=Bvec;
1784  //Bvec.zero();
1785  Distribute(c,Resid);
1786 
1787  makeValueArray(c,Bvec);
1788 
1789  AMat.multiply(Resid,Bvec);
1790  Resid=Bvec;
1791 
1792  //ArrayAbs(Bvec);
1793  //ArrayAbs(Resid);
1794  ArrayAbs(Bvec,Resid);
1795  Bsum+=Bvec;
1796  ResidSum+=Resid;
1797  }
1798  else
1799  throw CException("Unexpected value for boundary cell map.");
1800  }
1801  }
1802  for(int o=0;o<_numDir+3;o++)
1803  {
1804  temp[o]=ResidSum[o];
1805  }
1806  temp[_numDir+3]=traceSum;
1807 
1808  //MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, ResidSum.getData(), _numDir+3, MPI::DOUBLE, MPI::SUM);
1809  //MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &traceSum, 1, MPI::DOUBLE, MPI::SUM);
1810 #ifdef FVM_PARALLEL
1811  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, temp.getData(), _numDir+4, MPI::DOUBLE, MPI::SUM);
1812 #endif
1813  /*
1814  for(int o=0;o<_numDir+3;o++)
1815  {
1816  ResidScalar+=ResidSum[o];
1817  }
1818 
1819  ResidScalar/=traceSum;
1820  */
1821 
1822  for(int o=0;o<_numDir+3;o++)
1823  {
1824  ResidScalar+=temp[o];
1825  }
1826 
1827  ResidScalar/=temp[_numDir+3];
1828 
1829  if(_aveResid==-1)
1830  {_aveResid=ResidScalar;}
1831  else
1832  {
1833  _residChange=fabs(_aveResid-ResidScalar)/_aveResid;
1834  _aveResid=ResidScalar;
1835  }
1836  }
1837 
1838  T getResidChange() {return _residChange;}
1839  T getAveResid() {return _aveResid;}
1840 
1841  void setfgFinder()
1842  {
1843  foreach(const FaceGroupPtr fgPtr, _mesh.getBoundaryFaceGroups())
1844  {
1845  const FaceGroup& fg=*fgPtr;
1846  const int off=fg.site.getOffset();
1847  const int cnt=fg.site.getCount();
1848  const int id=fg.id;
1849  VecInt2 BegEnd;
1850 
1851  BegEnd[0]=off;
1852  BegEnd[1]=off+cnt;
1853  _fgFinder[id]=BegEnd;
1854  }
1855  }
1856 
1857  void addFAS(const int c, TArray& bVec)
1858  {
1859  int count=0;
1860  for(int dir=0;dir<_numDir;dir++)
1861  {
1862  TArray& fasArray=*_fasArrays[dir];
1863  bVec[count]-=fasArray[c];
1864  count+=1;
1865  }
1866  VectorT3Array& fasArray = *_velocityFASCorrection;
1867  bVec[_numDir]-=fasArray[c][0];
1868  bVec[_numDir+1]-=fasArray[c][1];
1869  bVec[_numDir+2]-=fasArray[c][2];
1870  }
1871 
1872  int findFgId(const int faceIndex)
1873  {
1874  FaceToFg::iterator id;
1875  for(id=_fgFinder.begin();id!=_fgFinder.end();id++)
1876  {
1877  if(id->second[0]<=faceIndex && id->second[1]>faceIndex)
1878  return id->first;
1879  }
1880  throw CException("Didn't find matching FaceGroup!");
1881  return -1;
1882  }
1883 
1884  void ArrayAbs(TArray& o)
1885  {
1886  int length=o.getLength();
1887  for(int i=0;i<length;i++)
1888  o[i]=fabs(o[i]);
1889  }
1890 
1891  void ArrayAbs(TArray& o1, TArray& o2)
1892  {
1893  int length=o1.getLength();
1894  for(int i=0;i<length;i++)
1895  {
1896  o1[i]=fabs(o1[i]);
1897  o2[i]=fabs(o2[i]);
1898  }
1899  }
1900 
1901  void makeValueArray(const int c, TArray& o)
1902  {
1903  for(int dir=0;dir<_numDir;dir++)
1904  {
1905  const TArray& f = *_fArrays[dir];
1906  o[dir]=f[c];
1907  }
1908  const VectorT3Array& v = _velocity;
1909  o[_numDir]=v[c][0];
1910  o[_numDir+1]=v[c][1];
1911  o[_numDir+2]=v[c][2];
1912  }
1913 
1914  void setBoundaryValFine(const int cell, const int cellcount, const GradMatrix& gMat)
1915  {
1916  const int neibcount=_cellFaces.getCount(cell);
1917  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[_cells]);
1918  const StorageSite& ibFaces = _mesh.getIBFaces();
1919  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[_faces]);
1920  GradArray Grads(_numDir);
1921  Grads.zero();
1922 
1923  VectorT3 Gcoeff;
1924  TArray limitCoeff1(_numDir);
1925  TArray min1(_numDir);
1926  TArray max1(_numDir);
1927  for(int dir=0;dir<_numDir;dir++)
1928  {
1929  limitCoeff1[dir]=T(1.e20);
1930  const TArray& f = *_fArrays[dir];
1931  min1[dir]=f[cell];
1932  max1[dir]=f[cell];
1933  }
1934 
1935  const VectorT3Array& faceCoords=
1936  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_faces]);
1937  const VectorT3Array& cellCoords=
1938  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells]);
1939 
1940  for(int j=0;j<neibcount;j++)
1941  {
1942  const int face=_cellFaces(cell,j);
1943  int cell2=_faceCells(face,1);
1944  if(cell2==cell)
1945  cell2=_faceCells(face,0);
1946 
1947  Gcoeff=gMat.getCoeff(cell, cell2);
1948 
1949  for(int dir=0;dir<_numDir;dir++)
1950  {
1951  const TArray& f = *_fArrays[dir];
1952  Grads[dir].accumulate(Gcoeff,f[cell2]-f[cell]);
1953  if(min1[dir]>f[cell2])min1[dir]=f[cell2];
1954  if(max1[dir]<f[cell2])max1[dir]=f[cell2];
1955  }
1956  }
1957 
1958  for(int j=0;j<neibcount;j++)
1959  {
1960  const int face=_cellFaces(cell,j);
1961 
1962  SuperbeeLimiter lf;
1963  for(int dir=0;dir<_numDir;dir++)
1964  {
1965  const TArray& f = *_fArrays[dir];
1966  GradType& grad=Grads[dir];
1967 
1968  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
1969  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
1970 
1971  computeLimitCoeff(limitCoeff1[dir], f[cell], SOU, min1[dir], max1[dir], lf);
1972  }
1973  }
1974 
1975  for(int j=0;j<neibcount;j++)
1976  {
1977  const int face=_cellFaces(cell,j);
1978  int cell2=_faceCells(face,1);
1979  VectorT3 Af=_faceArea[face];
1980  VectorT3 en = _faceArea[face]/_faceAreaMag[face];
1981 
1982  if(cell2==cell)
1983  {
1984  Af=Af*(-1.);
1985  en=en*(-1.);
1986  cell2=_faceCells(face,0);
1987  }
1988 
1989  if ((_BCfArray[face]==5)||(_BCfArray[face]==4)||(_BCfArray[face]==2))
1990  {
1991  for(int dir=0;dir<_numDir;dir++)
1992  {
1993  TArray& f = *_fArrays[dir];
1994  const T c_dot_en = _cx[dir]*en[0]+_cy[dir]*en[1]+_cz[dir]*en[2];
1995  GradType& grad=Grads[dir];
1996  if(c_dot_en>T_Scalar(0))
1997  {
1998  VectorT3 fVec=faceCoords[face]-cellCoords[cell];
1999  VectorT3 rVec=cellCoords[cell2]-cellCoords[cell];
2000  T r=gMat.computeR(grad,f,rVec,cell,cell2);
2001  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2]);
2002  f[cell2]=(f[cell]+limitCoeff1[dir]*SOU);
2003  }
2004  }
2005  }
2006  }
2007  }
2008 
2009  private:
2010  const Mesh& _mesh;
2011  const GeomFields& _geomFields;
2012  const StorageSite& _cells;
2013  const StorageSite& _faces;
2014  const StorageSite& _solidFaces;
2015  const CRConnectivity& _cellFaces;
2016  const CRConnectivity& _faceCells;
2017  const Field& _areaMagField;
2018  const TArray& _faceAreaMag;
2019  const Field& _areaField;
2020  const VectorT3Array& _faceArea;
2021  const TArray& _cellVolume;
2022  MacroFields& _macroFields;
2023  TQuad& _quadrature;
2024  TDistFF& _dsfPtr;
2025  TDistFF& _dsfPtr1;
2026  TDistFF& _dsfPtr2;
2027  TDistFF& _dsfEqPtrES;
2028  TDistFF& _dsfPtrRes;
2029  TDistFF& _dsfPtrFAS;
2030  const T _dT;
2031  const int _order;
2032  const bool _transient;
2033  const T _underRelaxation;
2034  const T _rho_init;
2035  const T _T_init;
2036  const T _MW;
2037  const int _conOrder;
2038  COMETBCMap& _bcMap;
2039  map<int, vector<int> > _faceReflectionArrayMap;
2040  const IntArray& _BCArray;
2041  const IntArray& _BCfArray;
2042  const IntArray& _ZCArray;
2043  T _aveResid;
2044  T _residChange;
2045  FaceToFg _fgFinder;
2046  const int _numDir;
2047  const TArray& _cx;
2048  const TArray& _cy;
2049  const TArray& _cz;
2050  const TArray& _wts;
2051  TArray& _density;
2052  TArray& _pressure;
2053  VectorT3Array& _velocity;
2054  VectorT3Array& _velocityResidual;
2055  VectorT3Array* _velocityFASCorrection;
2056  TArray& _temperature;
2057  VectorT6Array& _stress;
2058  TArray& _collisionFrequency;
2059  VectorT10Array& _coeffg;
2060 
2061  TArray** _fArrays;
2062  TArray** _fN1Arrays;
2063  TArray** _fN2Arrays;
2064  TArray** _fEqESArrays;
2065  TArray** _fResArrays;
2066  TArray** _fasArrays;
2067 
2068 };
2069 
2070 
2071 #endif
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
COMETESBGKDiscretizer::TDistFF
DistFunctFields< T > TDistFF
Definition: COMETESBGKDiscretizer.h:54
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
COMETESBGKDiscretizer::VectorT10Array
Array< VectorT10 > VectorT10Array
Definition: COMETESBGKDiscretizer.h:53
COMETESBGKDiscretizer::Distribute
void Distribute(const int cell, TArray &BVec, TArray &Rvec)
Definition: COMETESBGKDiscretizer.h:1501
COMETESBGKDiscretizer::COMETBCMap
map< int, COMETBC< T > * > COMETBCMap
Definition: COMETESBGKDiscretizer.h:43
COMETESBGKDiscretizer::_geomFields
const GeomFields & _geomFields
Definition: COMETESBGKDiscretizer.h:2011
COMETESBGKDiscretizer::COMETMacro
void COMETMacro(const int cell, MatrixType Amat, TArray &BVec)
Definition: COMETESBGKDiscretizer.h:1464
COMETESBGKDiscretizer::setBoundaryValFine
void setBoundaryValFine(const int cell, const int cellcount, const GradMatrix &gMat)
Definition: COMETESBGKDiscretizer.h:1914
COMETESBGKDiscretizer::VectorT3
Vector< T_Scalar, 3 > VectorT3
Definition: COMETESBGKDiscretizer.h:36
COMETESBGKDiscretizer::COMETSolveFine
void COMETSolveFine(const int sweep, const int level)
Definition: COMETESBGKDiscretizer.h:153
Field::hasArray
bool hasArray(const StorageSite &s) const
Definition: Field.cpp:37
FaceGroup
Definition: Mesh.h:28
COMETESBGKDiscretizer::ComputeMacroparameters
void ComputeMacroparameters(const int cell)
Definition: COMETESBGKDiscretizer.h:1535
Field
Definition: Field.h:14
COMETESBGKDiscretizer::TArrow
ArrowHeadMatrix< T, 3 > TArrow
Definition: COMETESBGKDiscretizer.h:41
COMETESBGKDiscretizer::findResid
void findResid(const bool plusFAS)
Definition: COMETESBGKDiscretizer.h:1707
GradientMatrix::computeR
X computeR(const Gradient< X > &g, const Array< X > &x, const Coord dist, int i, int j) const
Definition: GradientMatrix.h:97
Mesh
Definition: Mesh.h:49
GradientMatrix::getCoeff
Coord & getCoeff(const int i, const int j)
Definition: GradientMatrix.h:129
COMETESBGKDiscretizer::_cellFaces
const CRConnectivity & _cellFaces
Definition: COMETESBGKDiscretizer.h:2015
COMETESBGKDiscretizer::COMETConvection
void COMETConvection(const int cell, TArrow &Amat, TArray &BVec, const int cellcount)
Definition: COMETESBGKDiscretizer.h:552
COMETESBGKDiscretizer::COMETCollision
void COMETCollision(const int cell, MatrixType Amat, TArray &BVec)
Definition: COMETESBGKDiscretizer.h:1426
COMETESBGKDiscretizer::COMETSolve
void COMETSolve(const int sweep, const int level)
Definition: COMETESBGKDiscretizer.h:231
computeLimitCoeff
void computeLimitCoeff(T &lc, const T x, const T &dx, const T &min, const T &max, const LimitFunc &f)
Definition: FluxLimiters.h:72
COMETESBGKDiscretizer::TSquareESBGK
SquareMatrixESBGK< T > TSquareESBGK
Definition: COMETESBGKDiscretizer.h:42
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
COMETESBGKDiscretizer::COMETConvectionFine
void COMETConvectionFine(const int cell, TArrow &Amat, TArray &BVec, const GradMatrix &gMat)
Definition: COMETESBGKDiscretizer.h:619
FaceGroup::id
const int id
Definition: Mesh.h:41
COMETESBGKDiscretizer::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: COMETESBGKDiscretizer.h:35
COMETESBGKDiscretizer::_faceArea
const VectorT3Array & _faceArea
Definition: COMETESBGKDiscretizer.h:2020
ArrowHeadMatrix::multiply
void multiply(const XArray &x, XArray &b)
Definition: esbgkbase/ArrowHeadMatrix.h:108
Array::getData
virtual void * getData() const
Definition: Array.h:275
ArrowHeadMatrix::getElement
T_Scalar & getElement(const int i, const int j)
Definition: esbgkbase/ArrowHeadMatrix.h:50
SQR
#define SQR(x)
Definition: COMETESBGKDiscretizer.h:28
COMETESBGKDiscretizer::COMETESBGKDiscretizer
COMETESBGKDiscretizer(const Mesh &mesh, const GeomFields &geomfields, const StorageSite &solidFaces, MacroFields &macroFields, TQuad &quadrature, TDistFF &dsfPtr, TDistFF &dsfPtr1, TDistFF &dsfPtr2, TDistFF &dsfEqPtrES, TDistFF &dsfPtrRes, TDistFF &dsfPtrFAS, const T dT, const int order, const bool transient, const T underRelaxation, const T rho_init, const T T_init, const T MW, const int conOrder, COMETBCMap &bcMap, map< int, vector< int > > faceReflectionArrayMap, const IntArray &BCArray, const IntArray &BCfArray, const IntArray &ZCArray)
Definition: COMETESBGKDiscretizer.h:61
COMETESBGKDiscretizer::findResidFine
void findResidFine(const bool plusFAS)
Definition: COMETESBGKDiscretizer.h:1573
GeomFields::ibFaceIndex
Field ibFaceIndex
Definition: GeomFields.h:40
COMETESBGKDiscretizer::GradMatrix
GradModelType::GradMatrixType GradMatrix
Definition: COMETESBGKDiscretizer.h:59
COMETESBGKDiscretizer::_faceReflectionArrayMap
map< int, vector< int > > _faceReflectionArrayMap
Definition: COMETESBGKDiscretizer.h:2039
COMETESBGKDiscretizer::ArrayAbs
void ArrayAbs(TArray &o1, TArray &o2)
Definition: COMETESBGKDiscretizer.h:1891
GradientModel< T >::getGradientMatrix
static GradMatrixType & getGradientMatrix(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:454
COMETESBGKDiscretizer::VectorT5Array
Array< VectorT5 > VectorT5Array
Definition: COMETESBGKDiscretizer.h:49
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Array
Definition: Array.h:14
COMETESBGKDiscretizer::GradModelType
GradientModel< T > GradModelType
Definition: COMETESBGKDiscretizer.h:58
COMETESBGKDiscretizer::COMETConvectionFine
void COMETConvectionFine(const int cell, TArrow &Amat, TArray &BVec, const int cellcount, const GradMatrix &gMat)
Definition: COMETESBGKDiscretizer.h:338
MacroFields::velocity
Field velocity
Definition: MacroFields.h:15
StorageSite::getOffset
int getOffset() const
Definition: StorageSite.h:87
COMETESBGKDiscretizer::FaceToFg
map< int, VecInt2 > FaceToFg
Definition: COMETESBGKDiscretizer.h:47
COMETESBGKDiscretizer::addFAS
void addFAS(const int c, TArray &bVec)
Definition: COMETESBGKDiscretizer.h:1857
COMETESBGKDiscretizer::makeValueArray
void makeValueArray(const int c, TArray &o)
Definition: COMETESBGKDiscretizer.h:1901
COMETESBGKDiscretizer::COMETUnsteady
void COMETUnsteady(const int cell, MatrixType Amat, TArray &BVec)
Definition: COMETESBGKDiscretizer.h:306
COMETESBGKDiscretizer::GradArray
Array< GradType > GradArray
Definition: COMETESBGKDiscretizer.h:57
COMETESBGKDiscretizer::_solidFaces
const StorageSite & _solidFaces
Definition: COMETESBGKDiscretizer.h:2014
ArrowHeadMatrix::Solve
void Solve(XArray &bVec)
Definition: esbgkbase/ArrowHeadMatrix.h:71
COMETESBGKDiscretizer::_velocityFASCorrection
VectorT3Array * _velocityFASCorrection
Definition: COMETESBGKDiscretizer.h:2055
COMETESBGKDiscretizer::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: COMETESBGKDiscretizer.h:37
MacroFields::velocityFASCorrection
Field velocityFASCorrection
Definition: MacroFields.h:18
DistFunctFields::dsf
std::vector< Field * > dsf
Definition: DistFunctFields.h:36
COMETESBGKDiscretizer::_cells
const StorageSite & _cells
Definition: COMETESBGKDiscretizer.h:2012
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
COMETESBGKDiscretizer::Distribute
void Distribute(const int cell, TArray &Rvec)
Definition: COMETESBGKDiscretizer.h:1521
COMETESBGKDiscretizer::_faces
const StorageSite & _faces
Definition: COMETESBGKDiscretizer.h:2013
COMETESBGKDiscretizer::COMETTest
void COMETTest(const int cell, MatrixType Amat, TArray &BVec, TArray &fV)
Definition: COMETESBGKDiscretizer.h:1379
COMETESBGKDiscretizer::COMETConvection
void COMETConvection(const int cell, TArrow &Amat, TArray &BVec)
Definition: COMETESBGKDiscretizer.h:1091
COMETESBGKDiscretizer::_faceCells
const CRConnectivity & _faceCells
Definition: COMETESBGKDiscretizer.h:2016
COMETESBGKDiscretizer::findFgId
int findFgId(const int faceIndex)
Definition: COMETESBGKDiscretizer.h:1872
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
COMETESBGKDiscretizer::VectorT6Array
Array< VectorT6 > VectorT6Array
Definition: COMETESBGKDiscretizer.h:51
Array::getLength
int getLength() const
Definition: Array.h:87