Memosa-FVM  0.2
esbgkbase/COMETModel.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 
6 #ifndef _COMETMODEL_H_
7 #define _COMETMODEL_H_
8 
9 #include <stdio.h>
10 #include <map>
11 #include <cmath>
12 #include <vector>
13 
14 #ifdef FVM_PARALLEL
15  #include <mpi.h>
16 #endif
17 
18 #include "Model.h"
19 
20 #include "Array.h"
21 #include "Vector.h"
22 
23 #include "Mesh.h"
24 
25 #include "Quadrature.h"
26 #include "DistFunctFields.h"
27 
28 #include "MacroFields.h"
29 #include "FlowFields.h"
30 
31 #include "COMETBC.h"
32 #include "COMETBoundaryConditions.h"
33 #include "COMETESBGKDiscretizer.h"
34 
35 #include "Linearizer.h"
36 #include "CRConnectivity.h"
37 #include "LinearSystem.h"
38 #include "Matrix.h"
39 #include "MultiField.h"
40 #include "MultiFieldMatrix.h"
41 #include "CRMatrix.h"
42 #include "FluxJacobianMatrix.h"
43 #include "DiagonalMatrix.h"
44 
45 #include "MatrixOperation.h"
46 #include "NumType.h"
47 
48 #include "StressTensor.h"
49 
50 template<class T>
51 class COMETModel : public Model
52 {
53  public:
54  typedef typename NumTypeTraits<T>:: T_Scalar T_Scalar;
55  typedef Array<int> IntArray;
56  typedef Array<T> TArray;
57  typedef shared_ptr<TArray> TArrptr;
58  typedef Array<bool> BArray;
59  typedef Array2D<T> TArray2D;
60  typedef Vector<T,3> VectorT3;
61  typedef Array<VectorT3> VectorT3Array;
62  typedef shared_ptr<VectorT3Array> VT3Ptr;
63  typedef StressTensor<T> StressTensorT6;
64  typedef Array<StressTensorT6> StressTensorArray;
65  typedef std::vector<Field*> stdVectorField;
66  typedef DistFunctFields<T> TDistFF;
67 
68  typedef Vector<T,5> VectorT5;
69  typedef Array<VectorT5> VectorT5Array;
70  typedef Vector<T,6> VectorT6;
71  typedef Array<VectorT6> VectorT6Array;
72  typedef Vector<T,10> VectorT10;
73  typedef Array<VectorT10> VectorT10Array;
74 
75  typedef shared_ptr<MeshList> MshLstPtr;
76  typedef shared_ptr<Mesh> MeshPtr;
77  typedef shared_ptr<GeomFields> GeoFldsPtr;
78  typedef shared_ptr<StorageSite> SSPtr;
79  typedef shared_ptr<CRConnectivity> CRPtr;
80 
81  typedef Array<int> BCfaceArray;
82  typedef shared_ptr<BCfaceArray> BfacePtr;
83  typedef vector<BfacePtr> BCfaceList;
84  typedef Array<int> BCcellArray;
85  typedef shared_ptr<BCcellArray> BCellPtr;
86  typedef vector<BCellPtr> BCcellList;
87 
88  typedef Quadrature<T> TQuad;
89 
90  typedef std::map<int,COMETBC<T>*> COMETBCMap;
91  typedef std::map<int,COMETVC<T>*> COMETVCMap;
92  typedef COMETModel<T> TCOMET;
93  typedef shared_ptr<TCOMET> TCOMETPtr;
94 
95  typedef MultiField::ArrayIndex Index;
96  typedef pair<Index,Index> EntryIndex;
97  typedef pair<const StorageSite*, const StorageSite*> SSPair;
98 
99  typedef map<EntryIndex,shared_ptr<Matrix> > MatrixMap;
100  typedef map<Index,int> MatrixSizeMap;
101  typedef map<const Mesh*,int> SizeMap;
102  typedef map<const StorageSite*,StorageSite*> SiteMap;
103 
104  typedef map<SSPair,shared_ptr<Array<int> > > MatrixMappersMap;
105 
106  typedef map<Index,shared_ptr<StorageSite> > StorageSiteMap;
107  typedef map<const StorageSite*,shared_ptr<StorageSite> > GhostStorageSiteMap;
108 
113  //MacroFields& macroFields;
114 
115  COMETModel(const MeshList& meshes, const int level, GeomFields& geomFields, MacroFields& macroFields, Quadrature<T>& quad, const int ibm=0,
116  GeomFields* finestGeomFields=NULL, const MeshList* finestMeshes=NULL, MacroFields* finestMacroFields=NULL):
117 
118  Model(meshes),
119  _level(level),
120  _geomFields(geomFields),
121  _quadrature(quad),
122  _macroFields(macroFields),
123  _dsfPtr(_meshes,_quadrature,"dsf_"),
124  _dsfPtr1(_meshes,_quadrature,"dsf1_"),
125  _dsfPtr2(_meshes,_quadrature,"dsf2_"),
126  _dsfEqPtr(_meshes,_quadrature,"dsfEq_"),
127  _dsfEqPtrES(_meshes,_quadrature,"dsfEqES_"),
128  _dsfPtr0(_meshes,_quadrature,"dsf0_"),
129  _dsfPtrInj(_meshes,_quadrature,"dsfInj_"),
130  _dsfPtrRes(_meshes,_quadrature,"dsfRes_"),
131  _dsfPtrFAS(_meshes,_quadrature,"dsfFAS_"),
132  _coarseGeomFields("coarse"),
133  _initialKmodelNorm(),
134  _niters(0),
135  _residual(0.0),
136  _initialResidual(0.0),
137  _ibm(ibm),
138  _finestGeomFields(finestGeomFields ? *finestGeomFields : _geomFields),
139  _finestMeshes(finestMeshes ? *finestMeshes : _meshes),
140  _finestMacroFields(finestMacroFields ? *finestMacroFields : _macroFields)
141  {
142 
143  const int numMeshes = _meshes.size();
144  for (int n=0; n<numMeshes; n++)
145  {
146  const Mesh& mesh = *_meshes[n];
147  const StorageSite& faces=mesh.getFaces();
148  const StorageSite& cells=mesh.getCells();
149  const int faceCount=faces.getCount();
150  const int cellCount=cells.getSelfCount();
151 
152  BfacePtr BFptr(new BCfaceArray(faceCount));
153  BFptr->zero();
154  _BFaces.push_back(BFptr);
155 
156  BCellPtr BCptr(new BCcellArray(cellCount));
157  _BCells.push_back(BCptr);
158  BCptr->zero();
159 
160  BCellPtr ZCptr(new BCcellArray(cellCount));
161  _ZCells.push_back(ZCptr);
162  ZCptr->zero();
163  if(_level==0)
164  {
165  COMETVC<T> *vc(new COMETVC<T>());
166  vc->vcType = "flow";
167  _vcMap[mesh.getID()] = vc;
168  }
169  }
170 
171  if(_level==0)
172  {
173  SetBoundaryConditions();
174  init();
175  InitializeMacroparameters();
176  initializeMaxwellian();
177  initializeFineMaxwellian();
178  ComputeMacroparameters(); //calculate density,velocity,temperature
179  ComputeFineMacroparameters();
180  ComputeCollisionfrequency(); //calculate viscosity, collisionFrequency
181  initializeMaxwellianEq(); //equilibrium distribution
182  }
183  }
184 
185  void init()
186 
187  {
188  const int numMeshes = _meshes.size();
189  for (int n=0; n<numMeshes; n++)
190  {
191  const Mesh& mesh = *_meshes[n];
192  const Mesh& fMesh = *_finestMeshes[n];
193 
194  const COMETVC<T>& vc = *_vcMap[mesh.getID()];
195 
196  const StorageSite& cells = mesh.getCells();
197  const StorageSite& fCells = fMesh.getCells();
198 
199  const int nCells = cells.getCountLevel1();
200  shared_ptr<VectorT3Array> vCell(new VectorT3Array(nCells));
201 
202  VectorT3 initialVelocity;
203  initialVelocity[0] = _options["initialXVelocity"];
204  initialVelocity[1] = _options["initialYVelocity"];
205  initialVelocity[2] = _options["initialZVelocity"];
206  *vCell = initialVelocity;
207  _macroFields.velocity.addArray(cells,vCell);
208 
209  shared_ptr<IntArray> fineToCoarseCell(new IntArray(nCells));
210  *fineToCoarseCell = -1;
211  _geomFields.fineToCoarse.addArray(cells,fineToCoarseCell);
212 
213  if((_ibm==1)&&(_level==0))
214  {
215  _geomFields.ibType.syncLocal();
216  shared_ptr<Array<Vector<int,25> > >finestToCoarseCell(new Array<Vector<int,25> >(fCells.getCountLevel1()));
217  Vector<int,25> initialIndex;
218  for(int k=0;k<25;k++)
219  initialIndex[k]=-1;
220  *finestToCoarseCell = initialIndex;
221  _finestGeomFields.finestToCoarse.addArray(fCells,finestToCoarseCell);
222  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
223  Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
224  for(int c=0;c<nCells;c++)
225  FinestToCoarse[c][_level]=c;
226  }
227 
228  shared_ptr<TArray> pCell(new TArray(nCells));
229  *pCell = _options["operatingPressure"];
230  _macroFields.pressure.addArray(cells,pCell);
231 
232  shared_ptr<TArray> rhoCell(new TArray(nCells));
233  *rhoCell = 0.;
234  //*rhoCell = vc["density"];
235  _macroFields.density.addArray(cells,rhoCell);
236 
237  shared_ptr<TArray> muCell(new TArray(nCells));
238  *muCell = 0.;
239  //*muCell = vc["viscosity"];
240  _macroFields.viscosity.addArray(cells,muCell);
241 
242  shared_ptr<TArray> tempCell(new TArray(nCells));
243  *tempCell = _options["operatingTemperature"];
244  _macroFields.temperature.addArray(cells,tempCell);
245 
246  shared_ptr<TArray> collFreqCell(new TArray(nCells));
247  *collFreqCell = 0.;
248  //*collFreqCell = vc["viscosity"];
249  _macroFields.collisionFrequency.addArray(cells,collFreqCell);
250 
251  //coeffs for perturbed BGK distribution function
252  shared_ptr<VectorT5Array> coeffCell(new VectorT5Array(nCells));
253  VectorT5 initialCoeff;
254  initialCoeff[0] = 1.0;
255  initialCoeff[1] = 1.0;
256  initialCoeff[2] = 0.0;
257  initialCoeff[3] = 0.0;
258  initialCoeff[4] = 0.0;
259  *coeffCell = initialCoeff;
260  _macroFields.coeff.addArray(cells,coeffCell);
261 
262  //coeffs for perturbed BGK distribution function
263  shared_ptr<VectorT10Array> coeffgCell(new VectorT10Array(nCells));
264  VectorT10 initialCoeffg;
265  initialCoeffg[0] = 1.0;
266  initialCoeffg[1] = 1.0;
267  initialCoeffg[2] = 0.0;
268  initialCoeffg[3] = 1.0;
269  initialCoeffg[4] = 0.0;
270  initialCoeffg[5] = 1.0;
271  initialCoeffg[6] = 0.0;
272  initialCoeffg[7] = 0.0;
273  initialCoeffg[8] = 0.0;
274  initialCoeffg[9] = 0.0;
275  *coeffgCell = initialCoeffg;
276  _macroFields.coeffg.addArray(cells,coeffgCell);
277 
278  // used for ESBGK equilibrium distribution function
279  shared_ptr<TArray> tempxxCell(new TArray(cells.getCountLevel1()));
280  *tempxxCell = _options["operatingTemperature"]/3;
281  _macroFields.Txx.addArray(cells,tempxxCell);
282 
283  shared_ptr<TArray> tempyyCell(new TArray(cells.getCountLevel1()));
284  *tempyyCell = _options["operatingTemperature"]/3;
285  _macroFields.Tyy.addArray(cells,tempyyCell);
286 
287  shared_ptr<TArray> tempzzCell(new TArray(cells.getCountLevel1()));
288  *tempzzCell = _options["operatingTemperature"]/3;
289  _macroFields.Tzz.addArray(cells,tempzzCell);
290 
291  shared_ptr<TArray> tempxyCell(new TArray(cells.getCountLevel1()));
292  *tempxyCell = 0.0;
293  _macroFields.Txy.addArray(cells,tempxyCell);
294 
295  shared_ptr<TArray> tempyzCell(new TArray(cells.getCountLevel1()));
296  *tempyzCell = 0.0;
297  _macroFields.Tyz.addArray(cells,tempyzCell);
298 
299  shared_ptr<TArray> tempzxCell(new TArray(cells.getCountLevel1()));
300  *tempzxCell = 0.0;
301  _macroFields.Tzx.addArray(cells,tempzxCell);
302 
303  //Entropy and Entropy Generation Rate for switching
304  shared_ptr<TArray> EntropyCell(new TArray(cells.getCountLevel1()));
305  *EntropyCell = 0.0;
306  _macroFields.Entropy.addArray(cells,EntropyCell);
307 
308  shared_ptr<TArray> EntropyGenRateCell(new TArray(cells.getCountLevel1()));
309  *EntropyGenRateCell = 0.0;
310  _macroFields.EntropyGenRate.addArray(cells,EntropyGenRateCell);
311 
312  shared_ptr<TArray> EntropyGenRateColl(new TArray(cells.getCountLevel1()));
313  *EntropyGenRateColl = 0.0;
314  _macroFields.EntropyGenRate_Collisional.addArray(cells,EntropyGenRateColl);
315 
316  //Pxx,Pyy,Pzz,Pxy,Pxz,Pyz
317  shared_ptr<VectorT6Array> stressCell(new VectorT6Array(nCells));
318  VectorT6 initialstress;
319  initialstress[0] = 1.0;
320  initialstress[1] = 1.0;
321  initialstress[2] = 1.0;
322  initialstress[3] = 0.0;
323  initialstress[4] = 0.0;
324  initialstress[5] = 0.0;
325  *stressCell = initialstress;
326  _macroFields.Stress.addArray(cells,stressCell);
327 
328  //Knq=M300+M120+M102 for Couette with uy
329  shared_ptr<TArray> KnqCell(new TArray(cells.getCountLevel1()));
330  *KnqCell = 0.0;
331  _macroFields.Knq.addArray(cells,KnqCell);
332 
333 
334  //higher order moments of distribution function
335  /*
336  shared_ptr<VectorT3Array> M300Cell(new VectorT3Array(cells.getCount()));
337  VectorT3 initialM300;
338  initialM300[0] = 0.0;
339  initialM300[1] = 0.0;
340  initialM300[2] = 0.0;
341  *M300Cell = initialM300;
342  _macroFields.M300.addArray(cells,M300Cell);
343 
344  shared_ptr<VectorT3Array> M030Cell(new VectorT3Array(cells.getCount()));
345  VectorT3 initialM030;
346  initialM030[0] = 0.0;
347  initialM030[1] = 0.0;
348  initialM030[2] = 0.0;
349  *M300Cell = initialM030;
350  _macroFields.M030.addArray(cells,M030Cell);
351  */
352  //if(MPI::COMM_WORLD.Get_rank()==0)
353  //cout<<"array for fields created"<<endl;
354  const int numDirections = _quadrature.getDirCount();
355  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
356  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
357  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
358  //FILE * pFile;
359  //pFile=fopen("ref_incMEMOSA.txt","w");
360  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups()){
361  const FaceGroup& fg = *fgPtr;
362 
363  if((_bcMap[fg.id]->bcType == "SymmetryBC")||(_bcMap[fg.id]->bcType == "RealWallBC")||(_bcMap[fg.id]->bcType == "VelocityInletBC")){
364 
365  const StorageSite& faces = fg.site;
366 
367  const Field& areaMagField = _geomFields.areaMag;
368  const TArray& faceAreaMag = dynamic_cast<const TArray &>(areaMagField[faces]);
369  const Field& areaField = _geomFields.area;
370  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
371 
372  const VectorT3 en = faceArea[0]/faceAreaMag[0];
373  vector<int> tempVec(numDirections);
374 
375  for (int j=0; j<numDirections; j++){
376  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
377  const T cx_incident = cx[j] - 2.0*c_dot_en*en[0];
378  const T cy_incident = cy[j] - 2.0*c_dot_en*en[1];
379  const T cz_incident = cz[j] - 2.0*c_dot_en*en[2];
380  int direction_incident=0;
381  T Rdotprod=1e54;
382  T dotprod=0.0;
383  for (int js=0; js<numDirections; js++){
384  dotprod=pow(cx_incident-cx[js],2)+pow(cy_incident-cy[js],2)+pow(cz_incident-cz[js],2);
385  if (dotprod< Rdotprod){
386  Rdotprod =dotprod;
387  direction_incident=js;}
388  }
389  tempVec[j] = direction_incident;
390  //fprintf(pFile,"%d %d %d \n",fg.id, j,direction_incident);
391 
392  }
393  const int fgid=fg.id;
394  _faceReflectionArrayMap[fgid] = tempVec; //add to map
395 
396  }
397  }
398 
399  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups()){
400  const FaceGroup& fg = *fgPtr;
401  if(fg.groupType == "NSinterface"){
402  const StorageSite& Intfaces = fg.site;
403 
404  shared_ptr<VectorT3Array> InterfaceVelFace(new VectorT3Array(Intfaces.getCount()));
405  InterfaceVelFace ->zero();
406  _macroFields.InterfaceVelocity.addArray(Intfaces,InterfaceVelFace);
407 
408  shared_ptr<StressTensorArray> InterfaceStressFace(new StressTensorArray(Intfaces.getCount()));
409  InterfaceStressFace ->zero();
410  _macroFields.InterfaceStress.addArray(Intfaces,InterfaceStressFace);
411 
412  shared_ptr<TArray> InterfacePressFace(new TArray(Intfaces.getCount()));
413  *InterfacePressFace = _options["operatingPressure"];
414  _macroFields.InterfacePressure.addArray(Intfaces,InterfacePressFace);
415 
416  shared_ptr<TArray> InterfaceDensityFace(new TArray(Intfaces.getCount()));
417  *InterfaceDensityFace =vc["density"];
418  _macroFields.InterfaceDensity.addArray(Intfaces,InterfaceDensityFace);
419 
420 
421  }
422 
423  //fclose(pFile);
424 
425 
426  } //end of loop through meshes
427 
428  BCcellArray& BCArray=*(_BCells[n]);
429  BCfaceArray& BCfArray=*(_BFaces[n]);
430  BCcellArray& ZCArray=*(_ZCells[n]);
431  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
432  {
433  const FaceGroup& fg = *fgPtr;
434  if((_bcMap[fg.id]->bcType == "WallBC")||(_bcMap[fg.id]->bcType == "RealWallBC"))
435  {
436  const StorageSite& faces = fg.site;
437  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
438  const int faceCount=faces.getCount();
439  const int offSet=faces.getOffset();
440 
441  for(int i=offSet;i<offSet+faceCount;i++)
442  BCfArray[i]=2;
443 
444  for(int i=0;i<faceCount;i++)
445  {
446  int cell1=BfaceCells(i,0);
447  BCArray[cell1]=1;
448  }
449  }
450  else if(_bcMap[fg.id]->bcType == "VelocityInletBC")
451  {
452  const StorageSite& faces = fg.site;
453  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
454  const int faceCount=faces.getCount();
455  const int offSet=faces.getOffset();
456 
457  for(int i=offSet;i<offSet+faceCount;i++)
458  BCfArray[i]=3;
459 
460  for(int i=0;i<faceCount;i++)
461  {
462  int cell1=BfaceCells(i,0);
463  BCArray[cell1]=1;
464  }
465  }
466  else if(_bcMap[fg.id]->bcType == "ZeroGradBC")
467  {
468  const StorageSite& faces = fg.site;
469  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
470  const int faceCount=faces.getCount();
471  const int offSet=faces.getOffset();
472 
473  for(int i=offSet;i<offSet+faceCount;i++)
474  BCfArray[i]=4;
475 
476  for(int i=0;i<faceCount;i++)
477  {
478  int cell1=BfaceCells(i,0);
479  ZCArray[cell1]=1;
480  }
481  }
482  else if((_bcMap[fg.id]->bcType == "PressureInletBC")||(_bcMap[fg.id]->bcType == "PressureOutletBC"))
483  {
484  const StorageSite& faces = fg.site;
485  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
486  const int faceCount=faces.getCount();
487  const int offSet=faces.getOffset();
488 
489  for(int i=offSet;i<offSet+faceCount;i++)
490  BCfArray[i]=5;
491  }
492  else if(_bcMap[fg.id]->bcType == "SymmetryBC")
493  {
494  const StorageSite& faces = fg.site;
495  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
496  const int faceCount=faces.getCount();
497  const int offSet=faces.getOffset();
498 
499  for(int i=offSet;i<offSet+faceCount;i++)
500  BCfArray[i]=6;
501 
502  for(int i=0;i<faceCount;i++)
503  {
504  int cell1=BfaceCells(i,0);
505  BCArray[cell1]=1;
506  }
507  }
508  else
509  {
510  const StorageSite& faces = fg.site;
511  const int faceCount=faces.getCount();
512  const int offSet=faces.getOffset();
513 
514  for(int i=offSet;i<offSet+faceCount;i++)
515  BCfArray[i]=0;
516  }
517  }
518  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
519  {
520  const FaceGroup& fg = *fgPtr;
521  const StorageSite& faces = fg.site;
522  const CRConnectivity& BfaceCells=mesh.getFaceCells(faces);
523  const int faceCount=faces.getCount();
524  const int offSet=faces.getOffset();
525 
526  for(int i=offSet;i<offSet+faceCount;i++)
527  BCfArray[i]=-1;
528  /*
529  if(MPI::COMM_WORLD.Get_rank()==1)
530  {
531  int fC = (mesh.getFaces()).getCount();
532  cout<<"level,rank,facecount,iID,offSet,ISize = "<<_level<<" "<<MPI::COMM_WORLD.Get_rank()<<" "<<fC<<" "<<fg.id<<" "<<offSet<<" "<<(offSet+faceCount)<<endl;
533  }
534  */
535  }
536 
537  const StorageSite& faces = mesh.getFaces();
538  const int faceCount = faces.getCount();
539  const CRConnectivity& faceCells=mesh.getFaceCells(faces);
540  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
541  for(int i=0;i<faceCount;i++)
542  {
543  const int c0 = faceCells(i,0);
544  const int c1 = faceCells(i,1);
545  if (((ibType[c0] == Mesh::IBTYPE_FLUID) && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
546  ((ibType[c1] == Mesh::IBTYPE_FLUID) && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
547  {
548  BCfArray[i]=7;
549  }
550  }
551 
552  int count = 0;
553  for(int c=0;c<cells.getSelfCount();c++)
554  {
555  if(ibType[c] != Mesh::IBTYPE_FLUID)
556  count++;
557  }
558 #ifdef FVM_PARALLEL
559  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &count, 1, MPI::INT, MPI::SUM);
560  if((MPI::COMM_WORLD.Get_rank()==0)&&(_level==0))
561  cout<<"number of non-fluid cells in mesh at level "<<_level<<" = "<<count<<endl;
562 #endif
563 
564 #ifndef FVM_PARALLEL
565  if(_level==0)
566  cout<<"number of non-fluid cells in mesh at level "<<_level<<" = "<<count<<endl;
567 #endif
568  _niters =0;
569  _initialKmodelNorm = MFRPtr();
570  //_initialKmodelvNorm = MFRPtr();
571 
572  }
573  }
574 
575  void MakeCoarseModel(TCOMET* finerModel)
576  {
577 
578  if(_options.AgglomerationMethod=="FaceArea")
579  {
580  int maxLevs=finerModel->getOptions().maxLevels;
581  int thisLevel=(finerModel->getLevel())+1;
582 
583  if(thisLevel<maxLevs) //assumes # of levels will always work for the mesh
584  {
585  MeshList* newMeshesPtr=new MeshList;
586  TQuad* newQuadPtr=new TQuad();
587  MacroFields* newMacroPtr=new MacroFields("coarse");
588 
589  newQuadPtr->CopyQuad(finerModel->getQuadrature());
590 
591  MakeCoarseMesh1(finerModel->getMeshList(),
592  finerModel->getGeomFields(),
593  *newMeshesPtr);
594 
595  _geomFields.fineToCoarse.syncLocal();
596  const Mesh& mesh = *_meshes[0];
597  const StorageSite& cells = mesh.getCells();
598  const int nCells = cells.getCount();
599  Field& FineToCoarseField=(finerModel->getGeomFields()).fineToCoarse;
600  const IntArray& coarseIndex=dynamic_cast<const IntArray&>(FineToCoarseField[cells]);
601 
602  /*
603  if(MPI::COMM_WORLD.Get_rank()==1)
604  for(int c=0;c<nCells;c++)
605  cout<<" after sync, rank, level, cell no and finetocoarse = "<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<c<<" "<<coarseIndex[c]<<endl;
606  */
607 
608  syncGhostCoarsening(finerModel->getMeshList(),
609  finerModel->getGeomFields(),
610  *newMeshesPtr);
611  const int numMeshes =_meshes.size();
612  for (int n=0; n<numMeshes; n++)
613  {
614  const Mesh& mesh = *_meshes[n];
615  const StorageSite& fineSite = mesh.getCells();
616  StorageSite& coarseSite = *_siteMap[&fineSite];
617  const StorageSite::ScatterMap& fineScatterMap = fineSite.getScatterMap();
618  const StorageSite::ScatterMap& fineScatterMapLevel1 = fineSite.getScatterMapLevel1();
619  StorageSite::ScatterMap& coarseScatterMap = coarseSite.getScatterMap();
620 
621  foreach(const StorageSite::ScatterMap::value_type& pos, fineScatterMap)
622  {
623  const StorageSite& fineOSite = *pos.first;
624 
625 #ifdef FVM_PARALLEL
626  // the ghost site will not have its corresponding coarse
627  // site created yet so we create it here
628  if (_siteMap.find(&fineOSite) == _siteMap.end())
629  {
630  shared_ptr<StorageSite> ghostSite
631  (new StorageSite(-1));
632  ghostSite->setGatherProcID ( fineOSite.getGatherProcID() );
633  ghostSite->setScatterProcID( fineOSite.getScatterProcID() );
634  ghostSite->setTag( fineOSite.getTag() );
635  StorageSite& coarseOSite = *ghostSite;
636  _siteMap[&fineOSite]=&coarseOSite;
637  _sharedSiteMap[&fineOSite]=ghostSite;
638  }
639 #endif
640 
641  StorageSite& coarseOSite = *_siteMap[&fineOSite];
642 
643  SSPair sskey(&fineSite,&fineOSite);
644  coarseScatterMap[&coarseOSite] = _coarseScatterMaps[sskey];
645  }
646 
647  foreach(const StorageSite::ScatterMap::value_type& pos, fineScatterMapLevel1)
648  {
649  const StorageSite& fineOSite = *pos.first;
650  SSPair sskey(&fineSite,&fineOSite);
651  if (_coarseScatterMaps.find(sskey) != _coarseScatterMaps.end())
652  {
653 
654 #ifdef FVM_PARALLEL
655  // the ghost site will not have its corresponding coarse
656  // site created yet so we create it here
657  if (_siteMap.find(&fineOSite) == _siteMap.end())
658  {
659  shared_ptr<StorageSite> ghostSite
660  (new StorageSite(-1));
661  ghostSite->setGatherProcID ( fineOSite.getGatherProcID() );
662  ghostSite->setScatterProcID( fineOSite.getScatterProcID() );
663  ghostSite->setTag( fineOSite.getTag() );
664  StorageSite& coarseOSite = *ghostSite;
665  _siteMap[&fineOSite]=&coarseOSite;
666  _sharedSiteMap[&fineOSite]=ghostSite;
667  }
668 #endif
669 
670  StorageSite& coarseOSite = *_siteMap[&fineOSite];
671 
672  coarseScatterMap[&coarseOSite] = _coarseScatterMaps[sskey];
673  }
674  }
675 
676  const StorageSite::GatherMap& fineGatherMap = fineSite.getGatherMap();
677  const StorageSite::GatherMap& fineGatherMapLevel1 = fineSite.getGatherMapLevel1();
678  StorageSite::GatherMap& coarseGatherMap = coarseSite.getGatherMap();
679  foreach(const StorageSite::GatherMap::value_type& pos, fineGatherMap)
680  {
681  const StorageSite& fineOSite = *pos.first;
682  StorageSite& coarseOSite = *_siteMap[&fineOSite];
683  SSPair sskey(&fineSite,&fineOSite);
684 
685  coarseGatherMap[&coarseOSite] = _coarseGatherMaps[sskey];
686  }
687 
688  foreach(const StorageSite::GatherMap::value_type& pos, fineGatherMapLevel1)
689  {
690  const StorageSite& fineOSite = *pos.first;
691  SSPair sskey(&fineSite,&fineOSite);
692  if (_coarseGatherMaps.find(sskey) != _coarseGatherMaps.end())
693  {
694  foreach(SiteMap::value_type tempPos, _siteMap)
695  {
696  const StorageSite& tempOSite = *tempPos.first;
697  if(fineOSite.getTag()==tempOSite.getTag())
698  {
699  //StorageSite& coarseOSite = *_siteMap[&fineOSite];
700  StorageSite& coarseOSite = *_siteMap[&tempOSite];
701  coarseGatherMap[&coarseOSite] = _coarseGatherMaps[sskey];
702  }
703  }
704  }
705  }
706 
707  }
708 
709  int newCount= MakeCoarseMesh2(finerModel->getMeshList(),
710  finerModel->getGeomFields(),_coarseGeomFields,
711  *newMeshesPtr);
712 
713  TCOMET* newModelPtr=new COMETModel(*newMeshesPtr,thisLevel,
714  _coarseGeomFields,
715  *newMacroPtr,*newQuadPtr);
716 #ifdef FVM_PARALLEL
717  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &newCount, 1, MPI::INT, MPI::SUM);
718  if(MPI::COMM_WORLD.Get_rank()==0)
719  cout<<"Number of cells in level "<<thisLevel<<" is "<<newCount<<endl;
720 #endif
721 
722 #ifndef FVM_PARALLEL
723  cout<<"Number of cells in level "<<thisLevel<<" is "<<newCount<<endl;
724 #endif
725 
726  newModelPtr->setFinerLevel(finerModel);
727  finerModel->setCoarserLevel(newModelPtr);
728  newModelPtr->getOptions()=finerModel->getOptions();
729  newModelPtr->getBCMap()=finerModel->getBCMap();
730  newModelPtr->getVCMap()=finerModel->getVCMap();
731 
732  newModelPtr->init();
733  newModelPtr->InitializeMacroparameters();
734  newModelPtr->initializeMaxwellian();
735  newModelPtr->initializeCoarseMaxwellian();
736  newModelPtr->ComputeMacroparameters();
737  newModelPtr->ComputeCoarseMacroparameters();
738  newModelPtr->ComputeCollisionfrequency();
739  newModelPtr->initializeMaxwellianEq();
740 
741  if(newCount>_options.minCells)
742  newModelPtr->MakeCoarseModel(newModelPtr);
743  else
744  _options.maxLevels=newModelPtr->getLevel();
745  }
746  }
747  else if(_options.AgglomerationMethod=="AMG")
748  throw CException("Have not implemented AMG agglomeration method.");
749  else
750  throw CException("Unknown agglomeration method.");
751  }
752 
753  void MakeIBCoarseModel(TCOMET* finerModel, const StorageSite& solidFaces)
754  {
755 
756  if(_options.AgglomerationMethod=="FaceArea")
757  {
758  int maxLevs=finerModel->getOptions().maxLevels;
759  int thisLevel=(finerModel->getLevel())+1;
760 
761  if(thisLevel<maxLevs) //assumes # of levels will always work for the mesh
762  {
763  MeshList* newMeshesPtr=new MeshList;
764  TQuad* newQuadPtr=new TQuad();
765  MacroFields* newMacroPtr=new MacroFields("coarse");
766 
767  newQuadPtr->CopyQuad(finerModel->getQuadrature());
768 
769  MakeCoarseMesh1(finerModel->getMeshList(),
770  finerModel->getGeomFields(),
771  *newMeshesPtr);
772 
773  _geomFields.fineToCoarse.syncLocal();
774  const Mesh& mesh = *_meshes[0];
775  const StorageSite& cells = mesh.getCells();
776  const int nCells = cells.getCount();
777  Field& FineToCoarseField=(finerModel->getGeomFields()).fineToCoarse;
778  const IntArray& coarseIndex=dynamic_cast<const IntArray&>(FineToCoarseField[cells]);
779 
780  /*
781  if(MPI::COMM_WORLD.Get_rank()==1)
782  for(int c=0;c<nCells;c++)
783  cout<<" after sync, rank, level, cell no and finetocoarse = "<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<c<<" "<<coarseIndex[c]<<endl;
784  */
785 
786  syncGhostCoarsening(finerModel->getMeshList(),
787  finerModel->getGeomFields(),
788  *newMeshesPtr);
789  const int numMeshes =_meshes.size();
790  for (int n=0; n<numMeshes; n++)
791  {
792  const Mesh& mesh = *_meshes[n];
793  const StorageSite& fineSite = mesh.getCells();
794  StorageSite& coarseSite = *_siteMap[&fineSite];
795  const StorageSite::ScatterMap& fineScatterMap = fineSite.getScatterMap();
796  const StorageSite::ScatterMap& fineScatterMapLevel1 = fineSite.getScatterMapLevel1();
797  StorageSite::ScatterMap& coarseScatterMap = coarseSite.getScatterMap();
798 
799  foreach(const StorageSite::ScatterMap::value_type& pos, fineScatterMap)
800  {
801  const StorageSite& fineOSite = *pos.first;
802 
803 #ifdef FVM_PARALLEL
804  // the ghost site will not have its corresponding coarse
805  // site created yet so we create it here
806  if (_siteMap.find(&fineOSite) == _siteMap.end())
807  {
808  shared_ptr<StorageSite> ghostSite
809  (new StorageSite(-1));
810  ghostSite->setGatherProcID ( fineOSite.getGatherProcID() );
811  ghostSite->setScatterProcID( fineOSite.getScatterProcID() );
812  ghostSite->setTag( fineOSite.getTag() );
813  StorageSite& coarseOSite = *ghostSite;
814  _siteMap[&fineOSite]=&coarseOSite;
815  _sharedSiteMap[&fineOSite]=ghostSite;
816  }
817 #endif
818 
819  StorageSite& coarseOSite = *_siteMap[&fineOSite];
820 
821  SSPair sskey(&fineSite,&fineOSite);
822  coarseScatterMap[&coarseOSite] = _coarseScatterMaps[sskey];
823  }
824 
825  foreach(const StorageSite::ScatterMap::value_type& pos, fineScatterMapLevel1)
826  {
827  const StorageSite& fineOSite = *pos.first;
828  SSPair sskey(&fineSite,&fineOSite);
829  if (_coarseScatterMaps.find(sskey) != _coarseScatterMaps.end())
830  {
831 
832 #ifdef FVM_PARALLEL
833  // the ghost site will not have its corresponding coarse
834  // site created yet so we create it here
835  if (_siteMap.find(&fineOSite) == _siteMap.end())
836  {
837  shared_ptr<StorageSite> ghostSite
838  (new StorageSite(-1));
839  ghostSite->setGatherProcID ( fineOSite.getGatherProcID() );
840  ghostSite->setScatterProcID( fineOSite.getScatterProcID() );
841  ghostSite->setTag( fineOSite.getTag() );
842  StorageSite& coarseOSite = *ghostSite;
843  _siteMap[&fineOSite]=&coarseOSite;
844  _sharedSiteMap[&fineOSite]=ghostSite;
845  }
846 #endif
847 
848  StorageSite& coarseOSite = *_siteMap[&fineOSite];
849 
850  coarseScatterMap[&coarseOSite] = _coarseScatterMaps[sskey];
851  }
852  }
853 
854  const StorageSite::GatherMap& fineGatherMap = fineSite.getGatherMap();
855  const StorageSite::GatherMap& fineGatherMapLevel1 = fineSite.getGatherMapLevel1();
856  StorageSite::GatherMap& coarseGatherMap = coarseSite.getGatherMap();
857  foreach(const StorageSite::GatherMap::value_type& pos, fineGatherMap)
858  {
859  const StorageSite& fineOSite = *pos.first;
860  StorageSite& coarseOSite = *_siteMap[&fineOSite];
861  SSPair sskey(&fineSite,&fineOSite);
862 
863  coarseGatherMap[&coarseOSite] = _coarseGatherMaps[sskey];
864  }
865  foreach(const StorageSite::GatherMap::value_type& pos, fineGatherMapLevel1)
866  {
867  const StorageSite& fineOSite = *pos.first;
868  SSPair sskey(&fineSite,&fineOSite);
869  if (_coarseGatherMaps.find(sskey) != _coarseGatherMaps.end())
870  {
871  foreach(SiteMap::value_type tempPos, _siteMap)
872  {
873  const StorageSite& tempOSite = *tempPos.first;
874  if(fineOSite.getTag()==tempOSite.getTag())
875  {
876  //StorageSite& coarseOSite = *_siteMap[&fineOSite];
877  StorageSite& coarseOSite = *_siteMap[&tempOSite];
878  coarseGatherMap[&coarseOSite] = _coarseGatherMaps[sskey];
879  }
880  }
881  }
882  }
883 
884  }
885 
886  int newCount= MakeCoarseMesh2(finerModel->getMeshList(),
887  finerModel->getGeomFields(),_coarseGeomFields,
888  *newMeshesPtr);
889 
890  _coarseGeomFields.ibType.syncLocal();
891 
892  TCOMET* newModelPtr=new COMETModel(*newMeshesPtr,thisLevel,
893  _coarseGeomFields,
894  *newMacroPtr,*newQuadPtr,1,&_finestGeomFields,&_finestMeshes,&_finestMacroFields);
895 #ifdef FVM_PARALLEL
896  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &newCount, 1, MPI::INT, MPI::SUM);
897  if(MPI::COMM_WORLD.Get_rank()==0)
898  cout<<"Number of cells in level "<<thisLevel<<" is "<<newCount<<endl;
899 #endif
900 
901 #ifndef FVM_PARALLEL
902  cout<<"Number of cells in level "<<thisLevel<<" is "<<newCount<<endl;
903 #endif
904 
905  newModelPtr->setFinerLevel(finerModel);
906  finerModel->setCoarserLevel(newModelPtr);
907  newModelPtr->getOptions()=finerModel->getOptions();
908  newModelPtr->getBCMap()=finerModel->getBCMap();
909  newModelPtr->getVCMap()=finerModel->getVCMap();
910 
911  for (int n=0; n<numMeshes; n++)
912  {
913  const Mesh& mesh = *_meshes[n];
914  const StorageSite& fineIBFaces = mesh.getIBFaces();
915  if(fineIBFaces.getCount()>0)
916  {
917  StorageSite& coarseIBFaces = *_siteMap[&fineIBFaces];
918  for(int dir=0;dir<_quadrature.getDirCount();dir++)
919  {
920  Field& fnd = *_dsfPtr.dsf[dir];
921  const TArray& fIB = dynamic_cast<const TArray&>(fnd[fineIBFaces]);
922  shared_ptr<TArray> cIBV(new TArray(coarseIBFaces.getCount()));
923  cIBV->zero();
924  DistFunctFields<T>& coarserdsf = _coarserLevel->getdsf();
925  Field& cfnd = *coarserdsf.dsf[dir];
926  cfnd.addArray(coarseIBFaces,cIBV);
927  TArray& cIB = dynamic_cast<TArray&>(cfnd[coarseIBFaces]);
928  for(int i=0;i<coarseIBFaces.getCount();i++)
929  cIB[i]=fIB[i];
930  }
931  }
932 
933  shared_ptr<VectorT3Array> coarseSolidVel(new VectorT3Array(solidFaces.getCount()));
934  const VectorT3Array& fineSolidVel = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[solidFaces]);
935  *coarseSolidVel = fineSolidVel;
936  newMacroPtr->velocity.addArray(solidFaces,coarseSolidVel);
937 
938  shared_ptr<TArray> coarseSolidDensity(new TArray(solidFaces.getCount()));
939  const TArray& fineSolidDensity = dynamic_cast<const TArray&>(_macroFields.density[solidFaces]);
940  *coarseSolidDensity = fineSolidDensity;
941  newMacroPtr->density.addArray(solidFaces,coarseSolidDensity);
942 
943  shared_ptr<TArray> coarseSolidTemperature(new TArray(solidFaces.getCount()));
944  const TArray& fineSolidTemperature = dynamic_cast<const TArray&>(_macroFields.temperature[solidFaces]);
945  *coarseSolidTemperature = fineSolidTemperature;
946  newMacroPtr->temperature.addArray(solidFaces,coarseSolidTemperature);
947  }
948 
949  newModelPtr->init();
950  newModelPtr->InitializeMacroparameters();
951  newModelPtr->initializeMaxwellian();
952  newModelPtr->initializeCoarseMaxwellian();
953  newModelPtr->ComputeMacroparameters();
954  newModelPtr->ComputeCoarseMacroparameters();
955  newModelPtr->ComputeCollisionfrequency();
956  newModelPtr->initializeMaxwellianEq();
957 
958  if(newCount>_options.minCells)
959  newModelPtr->MakeIBCoarseModel(newModelPtr,solidFaces);
960  else
961  _options.maxLevels=newModelPtr->getLevel();
962  }
963  }
964  else if(_options.AgglomerationMethod=="AMG")
965  throw CException("Have not implemented AMG agglomeration method.");
966  else
967  throw CException("Unknown agglomeration method.");
968  }
969 
970  void MakeCoarseIndex(const StorageSite& solidFaces)
971  {
972 
973  const int maxLevs=_options.maxLevels;
974  int thisLevel=_level+1;
975 
976  const Mesh& mesh = *_meshes[0];
977  const StorageSite& cells = mesh.getCells();
978  const int nCells = cells.getCount();
979  const Mesh& fMesh = *_finestMeshes[0];
980  const StorageSite& fCells = fMesh.getCells();
981  const int nFCells = fCells.getCount();
982  Field& FineToCoarseField=_geomFields.fineToCoarse;
983  const IntArray& FineToCoarse=dynamic_cast<const IntArray&>(FineToCoarseField[cells]);
984  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
985  Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
986 
987  if(_level==0)
988  MakeParallel();
989  else
990  {
991  for(int dir=0;dir<_quadrature.getDirCount();dir++)
992  {
993  Field& fnd = *_dsfPtr.dsf[dir];
994  shared_ptr<TArray> cSV(new TArray(solidFaces.getCount()));
995  cSV->zero();
996  fnd.addArray(solidFaces,cSV);
997  }
998  }
999 
1000  if(thisLevel<maxLevs) //assumes # of levels will always work for the mesh
1001  {
1002  for(int c=0;c<nFCells;c++)
1003  FinestToCoarse[c][_level+1]=FineToCoarse[FinestToCoarse[c][_level]];
1004 
1005  _coarserLevel->MakeCoarseIndex(solidFaces);
1006  }
1007  }
1008 
1009  void MakeParallel()
1010  {
1011  Field::syncLocalVectorFields( _dsfPtr.dsf );
1012 #if 0
1013  for(int dir=0;dir<_quadrature.getDirCount();dir++)
1014  {
1015  Field& fnd = *_dsfPtr.dsf[dir];
1016  fnd.syncLocal();
1017  }
1018 #endif
1019  }
1020 
1021  void InitializeMacroparameters()
1022  { const int numMeshes =_meshes.size();
1023  for (int n=0; n<numMeshes; n++)
1024  {
1025  const Mesh& mesh = *_meshes[n];
1026  const StorageSite& cells = mesh.getCells();
1027  const int nCells = cells.getCount();
1028  // const double pi(3.14159);
1029  const double pi=_options.pi;
1030  TArray& Entropy = dynamic_cast<TArray&>(_macroFields.Entropy[cells]);
1031  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
1032  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
1033 
1034  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1035  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[cells]);
1036 
1037  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
1038  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1039  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
1040  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
1041  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
1042  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
1043  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
1044  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
1045  //if ( MPI::COMM_WORLD.Get_rank() == 0 ) {cout << "ncells="<<nCells<<endl;}
1046 
1047  TArray& Knq = dynamic_cast<TArray&>(_macroFields.Knq[cells]);
1048  //initialize density,velocity
1049  for(int c=0; c<nCells;c++)
1050  {
1051  density[c] =1.0;
1052  v[c][0]=0.0;
1053  v[c][1]=0.0;
1054  v[c][2]=0.0;
1055  temperature[c]=1.0;
1056  pressure[c]=temperature[c]*density[c];
1057 
1058  //BGK
1059  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
1060  coeff[c][1]=1/temperature[c];
1061  coeff[c][2]=0.0;coeff[c][3]=0.0;coeff[c][4]=0.0;
1062 
1063  Entropy[c]=0.0;
1064  if(_options.fgamma ==2){
1065  //ESBGK
1066  coeffg[c][0]=coeff[c][0];
1067  coeffg[c][1]=coeff[c][1];
1068  coeffg[c][2]=coeff[c][2];
1069  coeffg[c][3]=coeff[c][1];
1070  coeffg[c][4]=coeff[c][3];
1071  coeffg[c][5]=coeff[c][1];
1072  coeffg[c][6]=coeff[c][4];
1073  coeffg[c][7]=0.0;
1074  coeffg[c][8]=0.0;
1075  coeffg[c][9]=0.0;
1076  }
1077 
1078  Txx[c]=0.5;
1079  Tyy[c]=0.5;
1080  Tzz[c]=0.5;
1081  Txy[c]=0.0;
1082  Tyz[c]=0.0;
1083  Tzx[c]=0.0;
1084 
1085  Knq[c]=0.0;
1086  }
1087  }
1088  }
1089  void InitializeFgammaCoefficients()
1090  {
1091  const int numMeshes =_meshes.size();
1092  for (int n=0; n<numMeshes; n++)
1093  {
1094  const Mesh& mesh = *_meshes[n];
1095  const StorageSite& cells = mesh.getCells();
1096  const int nCells = cells.getCount();
1097  // const double pi(3.14159);
1098  const double pi=_options.pi;
1099 
1100  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
1101  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1102  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
1103  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1104  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
1105  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
1106  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
1107  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
1108  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
1109  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
1110 
1111  for(int c=0; c<nCells;c++)
1112  {
1113 
1114  //BGK
1115  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
1116  coeff[c][1]=1/temperature[c];
1117  coeff[c][2]=0.0;coeff[c][3]=0.0;coeff[c][4]=0.0;
1118 
1119  if(_options.fgamma ==2){
1120  //ESBGK
1121  coeffg[c][0]=coeff[c][0];
1122  coeffg[c][1]=coeff[c][1];
1123  coeffg[c][2]=coeff[c][2];
1124  coeffg[c][3]=coeff[c][1];
1125  coeffg[c][4]=coeff[c][3];
1126  coeffg[c][5]=coeff[c][1];
1127  coeffg[c][6]=coeff[c][4];
1128  coeffg[c][7]=0.0;
1129  coeffg[c][8]=0.0;
1130  coeffg[c][9]=0.0;
1131 
1132  Txx[c]=0.5*temperature[c];
1133  Tyy[c]=0.5*temperature[c];
1134  Tzz[c]=0.5*temperature[c];
1135  Txy[c]=0.0;
1136  Tyz[c]=0.0;
1137  Tzx[c]=0.0;
1138  }
1139  }
1140  }
1141  }
1142 
1143  void ComputeMacroparameters()
1144  {
1145  //FILE * pFile;
1146  //pFile = fopen("distfun_mf.txt","w");
1147  const int numMeshes = _meshes.size();
1148  for (int n=0; n<numMeshes; n++)
1149  {
1150 
1151  const Mesh& mesh = *_meshes[n];
1152  const StorageSite& cells = mesh.getCells();
1153  const int nCells = cells.getCount(); //
1154 
1155 
1156  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
1157  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1158  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
1159  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[cells]);
1160  const int N123 = _quadrature.getDirCount();
1161 
1162  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1163  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1164  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1165  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1166 
1167  VectorT6Array& stress = dynamic_cast<VectorT6Array&>(_macroFields.Stress[cells]);
1168 
1169  //initialize density,velocity,temperature to zero
1170  for(int c=0; c<nCells;c++)
1171  {
1172  density[c]=0.0;
1173  v[c][0]=0.0;
1174  v[c][1]=0.0;
1175  v[c][2]=0.0;
1176  temperature[c]=0.0;
1177  stress[c][0]=0.0;stress[c][1]=0.0;stress[c][2]=0.0;
1178  stress[c][3]=0.0;stress[c][4]=0.0;stress[c][5]=0.0;
1179 
1180  }
1181 
1182 
1183  for(int j=0;j<N123;j++){
1184 
1185  Field& fnd = *_dsfPtr.dsf[j];
1186  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1187 
1188  //fprintf(pFile,"%d %12.6f %E %E %E %E \n",j,dcxyz[j],cx[j],cy[j],f[80],density[80]+dcxyz[j]*f[80]);
1189 
1190  for(int c=0; c<nCells;c++){
1191  density[c] = density[c]+wts[j]*f[c];
1192  v[c][0]= v[c][0]+(cx[j]*f[c])*wts[j];
1193  v[c][1]= v[c][1]+(cy[j]*f[c])*wts[j];
1194  v[c][2]= v[c][2]+(cz[j]*f[c])*wts[j];
1195  temperature[c]= temperature[c]+(pow(cx[j],2.0)+pow(cy[j],2.0)
1196  +pow(cz[j],2.0))*f[c]*wts[j];
1197 
1198  }
1199 
1200  }
1201 
1202 
1203 
1204  for(int c=0; c<nCells;c++){
1205  v[c][0]=v[c][0]/density[c];
1206  v[c][1]=v[c][1]/density[c];
1207  v[c][2]=v[c][2]/density[c];
1208  temperature[c]=temperature[c]-(pow(v[c][0],2.0)
1209  +pow(v[c][1],2.0)
1210  +pow(v[c][2],2.0))*density[c];
1211  temperature[c]=temperature[c]/(1.5*density[c]);
1212  pressure[c]=density[c]*temperature[c];
1213  }
1214 
1215  //Find Pxx,Pyy,Pzz,Pxy,Pyz,Pzx, etc in field
1216 
1217  for(int j=0;j<N123;j++){
1218  Field& fnd = *_dsfPtr.dsf[j];
1219  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1220  for(int c=0; c<nCells;c++){
1221  stress[c][0] +=pow((cx[j]-v[c][0]),2.0)*f[c]*wts[j];
1222  stress[c][1] +=pow((cy[j]-v[c][1]),2.0)*f[c]*wts[j];
1223  stress[c][2] +=pow((cz[j]-v[c][2]),2.0)*f[c]*wts[j];
1224  stress[c][3] +=(cx[j]-v[c][0])*(cy[j]-v[c][1])*f[c]*wts[j];
1225  stress[c][4] +=(cy[j]-v[c][1])*(cz[j]-v[c][2])*f[c]*wts[j];
1226  stress[c][5] +=(cz[j]-v[c][2])*(cx[j]-v[c][0])*f[c]*wts[j];
1227 
1228  }}
1229 
1230 
1231  }// end of loop over nmeshes
1232  //fclose(pFile);
1233  }
1234 
1235  void ComputeCOMETMacroparameters()
1236  {
1237  //FILE * pFile;
1238  //pFile = fopen("distfun_mf.txt","w");
1239  const int numMeshes = _meshes.size();
1240  for (int n=0; n<numMeshes; n++)
1241  {
1242 
1243  const Mesh& mesh = *_meshes[n];
1244  const StorageSite& cells = mesh.getCells();
1245  const int nCells = cells.getCount(); //
1246 
1247 
1248  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
1249  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1250  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1251  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[cells]);
1252  const int N123 = _quadrature.getDirCount();
1253 
1254  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1255  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1256  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1257  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1258 
1259  VectorT6Array& stress = dynamic_cast<VectorT6Array&>(_macroFields.Stress[cells]);
1260 
1261  //initialize density,velocity,temperature to zero
1262  for(int c=0; c<nCells;c++)
1263  {
1264  density[c]=0.0;
1265  temperature[c]=0.0;
1266  stress[c][0]=0.0;stress[c][1]=0.0;stress[c][2]=0.0;
1267  stress[c][3]=0.0;stress[c][4]=0.0;stress[c][5]=0.0;
1268 
1269  }
1270 
1271 
1272  for(int j=0;j<N123;j++){
1273 
1274  Field& fnd = *_dsfPtr.dsf[j];
1275  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1276 
1277  //fprintf(pFile,"%d %12.6f %E %E %E %E \n",j,dcxyz[j],cx[j],cy[j],f[80],density[80]+dcxyz[j]*f[80]);
1278 
1279  for(int c=0; c<nCells;c++){
1280  density[c] = density[c]+wts[j]*f[c];
1281  temperature[c]= temperature[c]+(pow(cx[j],2.0)+pow(cy[j],2.0)
1282  +pow(cz[j],2.0))*f[c]*wts[j];
1283 
1284  }
1285 
1286  }
1287 
1288 
1289 
1290  for(int c=0; c<nCells;c++){
1291  temperature[c]=temperature[c]-(pow(v[c][0],2.0)
1292  +pow(v[c][1],2.0)
1293  +pow(v[c][2],2.0))*density[c];
1294  temperature[c]=temperature[c]/(1.5*density[c]);
1295  pressure[c]=density[c]*temperature[c];
1296  }
1297 
1298  //Find Pxx,Pyy,Pzz,Pxy,Pyz,Pzx, etc in field
1299 
1300  for(int j=0;j<N123;j++){
1301  Field& fnd = *_dsfPtr.dsf[j];
1302  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1303  for(int c=0; c<nCells;c++){
1304  stress[c][0] +=pow((cx[j]-v[c][0]),2.0)*f[c]*wts[j];
1305  stress[c][1] +=pow((cy[j]-v[c][1]),2.0)*f[c]*wts[j];
1306  stress[c][2] +=pow((cz[j]-v[c][2]),2.0)*f[c]*wts[j];
1307  stress[c][3] +=(cx[j]-v[c][0])*(cy[j]-v[c][1])*f[c]*wts[j];
1308  stress[c][4] +=(cy[j]-v[c][1])*(cz[j]-v[c][2])*f[c]*wts[j];
1309  stress[c][5] +=(cz[j]-v[c][2])*(cx[j]-v[c][0])*f[c]*wts[j];
1310 
1311  }}
1312 
1313 
1314  }// end of loop over nmeshes
1315  //fclose(pFile);
1316  }
1317 
1318  void ComputeFineMacroparameters()
1319  {
1320  const int numMeshes = _meshes.size();
1321  const T zero(0.0);
1322  for (int n=0; n<numMeshes; n++)
1323  {
1324 
1325  const Mesh& mesh = *_meshes[n];
1326  const StorageSite& cells = mesh.getCells();
1327  const int nCells = cells.getCount(); //
1328 
1329  VectorT3 zeroVelocity;
1330  zeroVelocity[0] = zero;
1331  zeroVelocity[1] = zero;
1332  zeroVelocity[2] = zero;
1333 
1334  shared_ptr<VectorT3Array> vRCell(new VectorT3Array(nCells));
1335  *vRCell = zeroVelocity;
1336  _macroFields.velocityResidual.addArray(cells,vRCell);
1337 
1338  /*
1339  VectorT3Array& vR = dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[cells]);
1340 
1341  for(int c=0; c<nCells;c++)
1342  {
1343  vR[c][0]=0.0;
1344  vR[c][1]=0.0;
1345  vR[c][2]=0.0;
1346  }
1347  */
1348 
1349  }// end of loop over nmeshes
1350  //fclose(pFile);
1351  }
1352 
1353  void ComputeCoarseMacroparameters()
1354  {
1355  const int numMeshes = _meshes.size();
1356  const T zero(0.0);
1357  for (int n=0; n<numMeshes; n++)
1358  {
1359 
1360  const Mesh& mesh = *_meshes[n];
1361  const StorageSite& cells = mesh.getCells();
1362  const int nCells = cells.getCount(); //
1363 
1364  VectorT3 zeroVelocity;
1365  zeroVelocity[0] = zero;
1366  zeroVelocity[1] = zero;
1367  zeroVelocity[2] = zero;
1368 
1369  shared_ptr<VectorT3Array> vRCell(new VectorT3Array(nCells));
1370  *vRCell = zeroVelocity;
1371  _macroFields.velocityResidual.addArray(cells,vRCell);
1372 
1373  shared_ptr<VectorT3Array> vICell(new VectorT3Array(nCells));
1374  *vICell = zeroVelocity;
1375  _macroFields.velocityInjected.addArray(cells,vICell);
1376 
1377  shared_ptr<VectorT3Array> vFCell(new VectorT3Array(nCells));
1378  *vFCell = zeroVelocity;
1379  _macroFields.velocityFASCorrection.addArray(cells,vFCell);
1380 
1381  /*
1382  VectorT3Array& vR = dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[cells]);
1383  VectorT3Array& vI = dynamic_cast<VectorT3Array&>(_macroFields.velocityInjected[cells]);
1384  VectorT3Array& vF = dynamic_cast<VectorT3Array&>(_macroFields.velocityFASCorrection[cells]);
1385 
1386  for(int c=0; c<nCells;c++)
1387  {
1388  vR[c][0]=0.0;
1389  vR[c][1]=0.0;
1390  vR[c][2]=0.0;
1391 
1392  vI[c][0]=0.0;
1393  vI[c][1]=0.0;
1394  vI[c][2]=0.0;
1395 
1396  vF[c][0]=0.0;
1397  vF[c][1]=0.0;
1398  vF[c][2]=0.0;
1399  }
1400  */
1401 
1402  }// end of loop over nmeshes
1403  //fclose(pFile);
1404  }
1405 
1406  void ComputeMacroparametersESBGK()
1407  {
1408 
1409  const int numMeshes = _meshes.size();
1410  for (int n=0; n<numMeshes; n++)
1411  {
1412  const Mesh& mesh = *_meshes[n];
1413  const StorageSite& cells = mesh.getCells();
1414  const int nCells = cells.getCount();
1415 
1416  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1417  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1418 
1419  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
1420  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
1421  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
1422  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
1423  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
1424  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
1425 
1426  const int N123 = _quadrature.getDirCount();
1427 
1428  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1429  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1430  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1431  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1432 
1433  const double Pr=_options.Prandtl;
1434  //cout <<"Prandlt" <<Pr<<endl;
1435 
1436  //initialize density,velocity,temperature to zero
1437  for(int c=0; c<nCells;c++)
1438  {
1439  Txx[c]=0.0;
1440  Tyy[c]=0.0;
1441  Tzz[c]=0.0;
1442  Txy[c]=0.0;
1443  Tyz[c]=0.0;
1444  Tzx[c]=0.0;
1445  }
1446  for(int j=0;j<N123;j++){
1447 
1448  Field& fnd = *_dsfPtr.dsf[j];
1449  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1450  Field& fndEq = *_dsfEqPtr.dsf[j];
1451  const TArray& fgam = dynamic_cast<const TArray&>(fndEq[cells]);
1452  for(int c=0; c<nCells;c++){
1453  Txx[c]=Txx[c]+pow(cx[j]-v[c][0],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1454  Tyy[c]=Tyy[c]+pow(cy[j]-v[c][1],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j] ;
1455  Tzz[c]=Tzz[c]+pow(cz[j]-v[c][2],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1456  //Txy[c]=Txy[c]+(cx[j])*(cy[j])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1457  //Tyz[c]=Tyz[c]+(cy[j])*(cz[j])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j]; //Tzx[c]=Tzx[c]+(cz[j])*(cx[j])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1458 
1459 
1460  Txy[c]=Txy[c]+(cx[j]-v[c][0])*(cy[j]-v[c][1])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1461  Tyz[c]=Tyz[c]+(cy[j]-v[c][1])*(cz[j]-v[c][2])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1462  Tzx[c]=Tzx[c]+(cz[j]-v[c][2])*(cx[j]-v[c][0])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
1463  }
1464  }
1465 
1466  for(int c=0; c<nCells;c++){
1467  Txx[c]=Txx[c]/density[c];
1468  Tyy[c]=Tyy[c]/density[c];
1469  Tzz[c]=Tzz[c]/density[c];
1470  Txy[c]=Txy[c]/density[c];
1471  Tyz[c]=Tyz[c]/density[c];
1472  Tzx[c]=Tzx[c]/density[c];
1473  }
1474 
1475  }
1476  }
1477 
1478 
1479 
1480  /*
1481  * Collision frequency
1482  *
1483  *
1484  */
1485 
1486  void ComputeCollisionfrequency() {
1487  const int numMeshes = _meshes.size();
1488  for (int n=0; n<numMeshes; n++)
1489  {
1490 
1491  const T rho_init=_options["rho_init"];
1492  const T T_init= _options["T_init"];
1493  const T mu_w= _options["mu_w"];
1494  const T Tmuref= _options["Tmuref"];
1495  const T muref= _options["muref"];
1496  const T R=8314.0/_options["molecularWeight"];
1497  const T nondim_length=_options["nonDimLt"];
1498 
1499  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
1500 
1501 
1502  const Mesh& mesh = *_meshes[n];
1503  const StorageSite& cells = mesh.getCells();
1504  const int nCells = cells.getCount();
1505 
1506  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
1507  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[cells]);
1508  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1509 
1510  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
1511  for(int c=0; c<nCells;c++)
1512  {
1513  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
1514  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
1515  }
1516 
1517  if(_options.fgamma==2){
1518  for(int c=0; c<nCells;c++)
1519  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
1520  }
1521 
1522  }
1523  }
1524 
1525  void MomentHierarchy() {
1526  const int numMeshes = _meshes.size();
1527  for (int n=0; n<numMeshes; n++)
1528  {
1529  const int Knq_dir=_options.Knq_direction;
1530  const Mesh& mesh = *_meshes[n];
1531  const StorageSite& cells = mesh.getCells();
1532  const int nCells = cells.getCount();
1533  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1534  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1535  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1536  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1537  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
1538  TArray& Knq = dynamic_cast<TArray&>(_macroFields.Knq[cells]);
1539  const int num_directions = _quadrature.getDirCount();
1540  if (Knq_dir ==0){
1541  for(int j=0;j<num_directions;j++){
1542  Field& fnd = *_dsfPtr.dsf[j];
1543  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1544  for(int c=0; c<nCells;c++){
1545  Knq[c]=Knq[c]+0.5*f[c]*wts[j]*(pow(cx[j]-v[c][0],3.0)+(cx[j]-v[c][0])*pow(cy[j]-v[c][1],2.0)+(cx[j]-v[c][0])*pow(cz[j]-v[c][2],2.0));
1546  }
1547  }}
1548  else if(Knq_dir ==1){
1549  for(int j=0;j<num_directions;j++){
1550  Field& fnd = *_dsfPtr.dsf[j];
1551  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1552  for(int c=0; c<nCells;c++){
1553  Knq[c]=Knq[c]+0.5*f[c]*wts[j]*(pow(cy[j]-v[c][1],3.0)+(cy[j]-v[c][1])*pow(cx[j]-v[c][0],2.0)+(cy[j]-v[c][1])*pow(cz[j]-v[c][2],2.0));
1554  }
1555  }}
1556 
1557  else if(Knq_dir ==2){
1558  for(int j=0;j<num_directions;j++){
1559  Field& fnd = *_dsfPtr.dsf[j];
1560  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1561  for(int c=0; c<nCells;c++){
1562  Knq[c]=Knq[c]+0.5*f[c]*wts[j]*(pow(cz[j]-v[c][2],3.0)+(cz[j]-v[c][2])*pow(cx[j]-v[c][0],2.0)+(cz[j]-v[c][2])*pow(cy[j]-v[c][1],2.0));
1563  }
1564  }}
1565  }
1566  }
1567  void EntropyGeneration() {
1568  const int numMeshes = _meshes.size();
1569  for (int n=0; n<numMeshes; n++)
1570  {
1571  const T rho_init=_options["rho_init"];
1572  const T T_init= _options["T_init"];
1573  const T molwt=_options["molecularWeight"]*1E-26/6.023;
1574  const T R=8314.0/_options["molecularWeight"];
1575  const T u_init=pow(2.0*R*T_init,0.5);
1576  const T Planck=_options.Planck;
1577  const T h3bm4u3=pow(Planck,3)/ pow(molwt,4)*rho_init/pow(u_init,3);
1578  //cout << "h3bm4u3 " << h3bm4u3 <<endl;
1579  //cout <<" u_init "<<u_init<<" rho_init "<<rho_init<<endl;
1580  const Mesh& mesh = *_meshes[n];
1581  const StorageSite& cells = mesh.getCells();
1582  const int nCells = cells.getCount();
1583  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1584  TArray& Entropy = dynamic_cast<TArray&>(_macroFields.Entropy[cells]);
1585  TArray& EntropyGenRate_Collisional = dynamic_cast<TArray&>(_macroFields.EntropyGenRate_Collisional[cells]);
1586  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
1587  for(int c=0; c<nCells;c++){
1588  Entropy[c]=0.0;EntropyGenRate_Collisional[c]=0.0;
1589  }
1590  const int num_directions = _quadrature.getDirCount();
1591  if (_options.fgamma ==2){
1592  for(int j=0;j<num_directions;j++){
1593  Field& fnd = *_dsfPtr.dsf[j];
1594  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
1595  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1596  const TArray& fgam = dynamic_cast<const TArray&>(feqES[cells]);
1597  for(int c=0; c<nCells;c++){
1598  Entropy[c]=Entropy[c]+f[c]*wts[j]*(1-log(h3bm4u3*f[c]));
1599  EntropyGenRate_Collisional[c]+= (f[c]-fgam[c])*collisionFrequency[c]*log(h3bm4u3*f[c])*wts[j];
1600  }
1601  }
1602  }
1603  else{
1604  for(int j=0;j<num_directions;j++){
1605  Field& fnd = *_dsfPtr.dsf[j];
1606  Field& feq = *_dsfEqPtr.dsf[j];
1607  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1608  const TArray& fgam = dynamic_cast<const TArray&>(feq[cells]);
1609  for(int c=0; c<nCells;c++){
1610  Entropy[c]=Entropy[c]+f[c]*wts[j]*(1-log(h3bm4u3*f[c]));
1611  EntropyGenRate_Collisional[c]+=(f[c]-fgam[c])*collisionFrequency[c]*(log(h3bm4u3*f[c]))*wts[j];
1612  }
1613  }
1614  }
1615 
1616 
1617  }
1618  }
1619 
1620  void initializeMaxwellianEq()
1621  {
1622  const int numMeshes = _meshes.size();
1623  for (int n=0; n<numMeshes; n++)
1624  {
1625  const Mesh& mesh = *_meshes[n];
1626  const StorageSite& cells = mesh.getCells();
1627  const int nCells = cells.getCount();
1628  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1629  const VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
1630  const VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1631 
1632  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1633  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1634  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1635  const int numFields= _quadrature.getDirCount();
1636 
1637  for(int j=0;j< numFields;j++){
1638  Field& fndEq = *_dsfEqPtr.dsf[j];
1639  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
1640  for(int c=0; c<nCells;c++){
1641  fEq[c]=coeff[c][0]*exp(-coeff[c][1]*(pow(cx[j]-v[c][0],2)+pow(cy[j]-v[c][1],2)
1642  +pow(cz[j]-v[c][2],2))+coeff[c][2]*(cx[j]-v[c][0])
1643  +coeff[c][3]*(cy[j]-v[c][1])+coeff[c][4]*(cz[j]-v[c][2]));
1644  }
1645  }
1646 
1647  if(_options.fgamma==2){
1648 
1649  for(int j=0;j< numFields;j++){
1650  Field& fndEqES = *_dsfEqPtrES.dsf[j];
1651  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
1652  for(int c=0; c<nCells;c++){
1653  T Cc1=(cx[j]-v[c][0]);
1654  T Cc2=(cy[j]-v[c][1]);
1655  T Cc3=(cz[j]-v[c][2]);
1656  fEqES[c]=coeffg[c][0]*exp(-coeffg[c][1]*pow(Cc1,2)+coeffg[c][2]*Cc1
1657  -coeffg[c][3]*pow(Cc2,2)+coeffg[c][4]*Cc2
1658  -coeffg[c][5]*pow(Cc3,2)+coeffg[c][6]*Cc3
1659  +coeffg[c][7]*cx[j]*cy[j]+coeffg[c][8]*cy[j]*cz[j]
1660  +coeffg[c][9]*cz[j]*cx[j]);
1661  }
1662  }
1663  }
1664 
1665 
1666 
1667  }
1668  }
1669 
1670  void NewtonsMethodBGK(const int ktrial)
1671  {
1672  const int numMeshes = _meshes.size();
1673  for (int n=0; n<numMeshes; n++)
1674  {
1675  //cout << " NewtonsMethod" <<endl;
1676  const T tolx=_options["ToleranceX"];
1677  const T tolf=_options["ToleranceF"];
1678  const int sizeC=5;
1679  const Mesh& mesh = *_meshes[n];
1680  const StorageSite& cells = mesh.getCells();
1681  const int nCells = cells.getCount();
1682 
1683  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1684  const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
1685  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1686 
1687  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
1688 
1689  for(int c=0; c<nCells;c++){
1690 
1691  for (int trial=0;trial<ktrial;trial ++){
1692  SquareMatrix<T,sizeC> fjac(0);
1693  SquareMatrix<T,sizeC> fjacinv(0);
1694  VectorT5 fvec;
1695 
1696 
1697  fvec[0]=density[c];
1698  fvec[1]=density[c]*v[c][0];
1699  fvec[2]=density[c]*v[c][1];
1700  fvec[3]=density[c]*v[c][2];
1701  fvec[4]=1.5*density[c]*temperature[c]+density[c]*(pow(v[c][0],2)+pow(v[c][1],2)+pow(v[c][2],2.0));
1702 
1703 
1704  setJacobianBGK(fjac,fvec,coeff[c],v[c],c);
1705 
1706  //solve using GE or inverse
1707  T errf=0.;
1708  for (int row=0;row<sizeC;row++){errf+=fabs(fvec[row]);}
1709 
1710  if(errf <= tolf)
1711  break;
1712  VectorT5 pvec;
1713  for (int row=0;row<sizeC;row++){pvec[row]=-fvec[row];}//rhs
1714 
1715 
1716  //solve Ax=b for x
1717  //p=GE_elim(fjac,p,3);
1718  VectorT5 xvec;
1719  fjacinv=inverseGauss(fjac,sizeC);
1720 
1721 
1722 
1723  for (int row=0;row<sizeC;row++){
1724  xvec[row]=0.0;
1725  for (int col=0;col<sizeC;col++){
1726  xvec[row]+=fjacinv(row,col)*pvec[col];}
1727  }
1728  //check for convergence, update
1729 
1730  T errx=0.;
1731  for (int row=0;row<sizeC;row++){
1732  errx +=fabs(xvec[row]);
1733  coeff[c][row]+= xvec[row];
1734  }
1735 
1736 
1737  if(errx <= tolx)
1738  break;
1739 
1740  }
1741 
1742 
1743  }
1744  }
1745 
1746  }
1747 
1748  void EquilibriumDistributionBGK()
1749  {
1750  const int ktrial=_options.NewtonsMethod_ktrial;
1751  const int numMeshes = _meshes.size();
1752  for (int n=0; n<numMeshes; n++)
1753  {
1754  const Mesh& mesh = *_meshes[n];
1755  const StorageSite& cells = mesh.getCells();
1756  const int nCells = cells.getCount();
1757 
1758  //const double pi=_options.pi;
1759  //const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1760  //const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
1761 
1762  //initialize coeff
1763  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
1764 
1765  /*
1766  for(int c=0; c<nCells;c++){
1767  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
1768  coeff[c][1]=1/temperature[c];
1769  coeff[c][2]=0.0;
1770  coeff[c][3]=0.0;
1771  coeff[c][4]=0.0;
1772  }
1773  */
1774  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1775  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1776  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1777  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1778  const int numFields= _quadrature.getDirCount();
1779 
1780  //call Newtons Method
1781 
1782  NewtonsMethodBGK(ktrial);
1783 
1784  //calculate perturbed maxwellian for BGK
1785  for(int j=0;j< numFields;j++){
1786  Field& fndEq = *_dsfEqPtr.dsf[j];
1787  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
1788  for(int c=0; c<nCells;c++){
1789  fEq[c]=coeff[c][0]*exp(-coeff[c][1]*(pow(cx[j]-v[c][0],2)+pow(cy[j]-v[c][1],2)
1790  +pow(cz[j]-v[c][2],2))+coeff[c][2]*(cx[j]-v[c][0])
1791  +coeff[c][3]*(cy[j]-v[c][1])+coeff[c][4]*(cz[j]-v[c][2]));
1792  }
1793 
1794  }
1795  }
1796  }
1797 
1798 
1799  void setJacobianBGK(SquareMatrix<T,5>& fjac, VectorT5& fvec, const VectorT5& xn,const VectorT3& v,const int c)
1800  {
1801 
1802 
1803  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1804  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1805  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1806  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1807 
1808  const TArray2D& malphaBGK = dynamic_cast<const TArray2D&>(*_quadrature.malphaBGKPtr);
1809  const int numFields= _quadrature.getDirCount();
1810  VectorT5 mexp;
1811 
1812  for(int j=0;j< numFields;j++){
1813  T Cconst=pow(cx[j]-v[0],2.0)+pow(cy[j]-v[1],2.0)+pow(cz[j]-v[2],2.0);
1814  T Econst=xn[0]*exp(-xn[1]*Cconst+xn[2]*(cx[j]-v[0])+xn[3]*(cy[j]-v[1])+xn[4]*(cz[j]-v[2]))*wts[j];
1815 
1816  for (int row=0;row<5;row++){
1817  fvec[row]+= -Econst*malphaBGK(j,row); //smm
1818 
1819  //fvec[row]=tvec[row]+fvec[row]; //mma
1820  }
1821 
1822  mexp[0]=-Econst/xn[0];
1823  mexp[1]=Econst*Cconst;
1824  mexp[2]=-Econst*(cx[j]-v[0]);
1825  mexp[3]=-Econst*(cy[j]-v[1]);
1826  mexp[4]=-Econst*(cz[j]-v[2]);
1827  for (int row=0;row<5;row++){
1828  for (int col=0;col<5;col++){
1829  fjac(row,col)+=malphaBGK(j,row)*mexp[col]; //new
1830  }
1831  }
1832 
1833 
1834  }
1835 
1836 
1837 
1838  }
1839 
1840  void NewtonsMethodESBGK(const int ktrial)
1841  {
1842  // cout<< "Inside Newtons Method" <<endl;
1843  const int numMeshes = _meshes.size();
1844  for (int n=0; n<numMeshes; n++)
1845  {
1846  const T tolx=_options["ToleranceX"];
1847  const T tolf=_options["ToleranceF"];
1848  const int sizeC=10;
1849  const Mesh& mesh = *_meshes[n];
1850  const StorageSite& cells = mesh.getCells();
1851  const int nCells = cells.getCount();
1852 
1853  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1854 
1855  const TArray& Txx = dynamic_cast<const TArray&>(_macroFields.Txx[cells]);
1856  const TArray& Tyy = dynamic_cast<const TArray&>(_macroFields.Tyy[cells]);
1857  const TArray& Tzz = dynamic_cast<const TArray&>(_macroFields.Tzz[cells]);
1858  const TArray& Txy = dynamic_cast<const TArray&>(_macroFields.Txy[cells]);
1859  const TArray& Tyz = dynamic_cast<const TArray&>(_macroFields.Tyz[cells]);
1860  const TArray& Tzx = dynamic_cast<const TArray&>(_macroFields.Tzx[cells]);
1861 
1862  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1863 
1864  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1865 
1866  for(int c=0; c<nCells;c++){
1867  // {cout <<"NM:ES" <<c <<endl;}
1868  for (int trial=0;trial<ktrial;trial ++){
1869  SquareMatrix<T,sizeC> fjac(0);
1870  SquareMatrix<T,sizeC> fjacinv(0);
1871  Vector<T,sizeC> fvec;
1872  // if (c==_options.printCellNumber){cout <<"trial" <<trial <<endl;}
1873 
1874  fvec[0]=density[c];
1875  fvec[1]=density[c]*v[c][0];
1876  fvec[2]=density[c]*v[c][1];
1877  fvec[3]=density[c]*v[c][2];
1878  fvec[4]=density[c]*(pow(v[c][0],2)+Txx[c]);
1879  fvec[5]=density[c]*(pow(v[c][1],2)+Tyy[c]);
1880  fvec[6]=density[c]*(pow(v[c][2],2)+Tzz[c]);
1881  fvec[7]=density[c]*(v[c][0]*v[c][1]+Txy[c]);
1882  fvec[8]=density[c]*(v[c][1]*v[c][2]+Tyz[c]);
1883  fvec[9]=density[c]*(v[c][2]*v[c][0]+Tzx[c]);
1884 
1885  //calculate Jacobian
1886  setJacobianESBGK(fjac,fvec,coeffg[c],v[c],c);
1887 
1888 
1889  //solve using GaussElimination
1890  T errf=0.; //and Jacobian matrix in fjac.
1891  for (int row=0;row<sizeC;row++){errf+=fabs(fvec[row]);}
1892  if(errf <= tolf)
1893  break;
1894  Vector<T,sizeC> pvec;
1895  for (int row=0;row<sizeC;row++){pvec[row]=-fvec[row];}
1896 
1897 
1898  //solve Ax=b for x
1899  //p=GE_elim(fjac,p,3);
1900  Vector<T,sizeC> xvec;
1901  fjacinv=inverseGauss(fjac,sizeC);
1902  for (int row=0;row<sizeC;row++){
1903  xvec[row]=0.0;
1904  for (int col=0;col<sizeC;col++){
1905  xvec[row]+=fjacinv(row,col)*pvec[col];
1906 
1907  }
1908  }
1909  /*
1910  if (c==_options.printCellNumber){
1911  cout << " cg0 "<<coeffg[c][0]<<" cg1 "<<coeffg[c][1]<<" cg2 "<<coeffg[c][2] << endl;
1912  cout <<" cg3 " <<coeffg[c][3]<< " cg4 "<<coeffg[c][4]<<" cg5 "<<coeffg[c][5]<<" cg6 "<<coeffg[c][6] << endl;
1913  cout <<" cg7 " <<coeffg[c][7]<< " cg8 "<<coeffg[c][8]<<" cg9 "<<coeffg[c][9]<<endl;
1914 
1915 
1916  //cout << " fvec-ESBGK " << fvec[4] <<fvec[5]<<fvec[6]<<fvec[7]<<fvec[8]<<fvec[9] <<endl;
1917  FILE * pFile;
1918  pFile = fopen("fvecfjac.dat","wa");
1919  //fprintf(pFile,"%s %d \n","trial",trial);
1920  for (int mat_col=0;mat_col<sizeC;mat_col++){fprintf(pFile,"%12.4E",fvec[mat_col]);}
1921  fprintf(pFile,"\n");
1922  for (int mat_row=0;mat_row<sizeC;mat_row++){
1923  for (int mat_col=0;mat_col<sizeC;mat_col++){
1924  fprintf(pFile,"%12.4E",fjac(mat_row,mat_col));}
1925  fprintf(pFile,"\n");}
1926  // fprintf(pFile,"done \n");
1927  //inverse
1928  for (int mat_row=0;mat_row<sizeC;mat_row++){
1929  for (int mat_col=0;mat_col<sizeC;mat_col++){
1930  fprintf(pFile,"%12.4E",fjacinv(mat_row,mat_col));}
1931  fprintf(pFile,"\n");}
1932  //solution
1933  for (int mat_col=0;mat_col<sizeC;mat_col++){
1934  fprintf(pFile,"%12.4E",pvec[mat_col]);}
1935  }
1936  */
1937 
1938  //check for convergence, update
1939  T errx=0.;//%Check root convergence.
1940  for (int row=0;row<sizeC;row++){
1941  errx +=fabs(xvec[row]);
1942  coeffg[c][row]+= xvec[row];
1943  }
1944  //if (c==_options.printCellNumber){cout <<"errx "<<errx<<endl;}
1945  if(errx <= tolx)
1946  break;
1947 
1948  }
1949 
1950  }
1951  }
1952  }
1953  void EquilibriumDistributionESBGK()
1954  {
1955  ComputeMacroparametersESBGK();
1956  const int ktrial=_options.NewtonsMethod_ktrial;
1957  const int numMeshes = _meshes.size();
1958  for (int n=0; n<numMeshes; n++)
1959  {
1960  const Mesh& mesh = *_meshes[n];
1961  const StorageSite& cells = mesh.getCells();
1962  const int nCells = cells.getCount();
1963 
1964 
1965  //const VectorT5Array& coeff = dynamic_cast<const VectorT5Array&>(_macroFields.coeff[cells]);
1966  //initialize coeffg
1967  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1968  /*
1969  for(int c=0; c<nCells;c++){
1970  coeffg[c][0]=coeff[c][0];
1971  coeffg[c][1]=coeff[c][1];
1972  coeffg[c][2]=coeff[c][2];
1973  coeffg[c][3]=coeff[c][1];
1974  coeffg[c][4]=coeff[c][3];
1975  coeffg[c][5]=coeff[c][1];
1976  coeffg[c][6]=coeff[c][4];
1977  coeffg[c][7]=0.0;
1978  coeffg[c][8]=0.0;
1979  coeffg[c][9]=0.0;
1980  }
1981  */
1982  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1983  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1984  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1985  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1986  const int numFields= _quadrature.getDirCount();
1987 
1988  NewtonsMethodESBGK(ktrial);
1989 
1990  for(int j=0;j< numFields;j++){
1991  Field& fndEqES = *_dsfEqPtrES.dsf[j];
1992  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
1993  for(int c=0; c<nCells;c++){
1994  T Cc1=(cx[j]-v[c][0]);
1995  T Cc2=(cy[j]-v[c][1]);
1996  T Cc3=(cz[j]-v[c][2]);
1997  fEqES[c]=coeffg[c][0]*exp(-coeffg[c][1]*pow(Cc1,2)+coeffg[c][2]*Cc1
1998  -coeffg[c][3]*pow(Cc2,2)+coeffg[c][4]*Cc2
1999  -coeffg[c][5]*pow(Cc3,2)+coeffg[c][6]*Cc3
2000  +coeffg[c][7]*cx[j]*cy[j]+coeffg[c][8]*cy[j]*cz[j]
2001  +coeffg[c][9]*cz[j]*cx[j]);
2002  }
2003 
2004  }
2005  }
2006  }
2007 
2008  void setJacobianESBGK(SquareMatrix<T,10>& fjac, VectorT10& fvec, const VectorT10& xn,const VectorT3& v,const int c)
2009  {
2010  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2011  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2012  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2013  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
2014 
2015  const TArray2D& malphaESBGK = dynamic_cast<const TArray2D&>(*_quadrature.malphaESBGKPtr);
2016  const int numFields= _quadrature.getDirCount();
2017  VectorT10 mexp;
2018 
2019  for(int j=0;j< numFields;j++){
2020  T Cc1=cx[j]-v[0];
2021  T Cc2=cy[j]-v[1];
2022  T Cc3=cz[j]-v[2];
2023  T Econst=xn[0]*exp(-xn[1]*pow(Cc1,2)+xn[2]*Cc1-xn[3]*pow(Cc2,2)+ xn[4]*Cc2
2024  -xn[5]*pow(Cc3,2)+xn[6]*Cc3
2025  +xn[7]*cx[j]*cy[j]+xn[8]*cy[j]*cz[j]+xn[9]*cz[j]*cx[j])*wts[j];
2026 
2027  for (int row=0;row<10;row++){
2028  fvec[row]+= -Econst*malphaESBGK(j,row); //smm
2029  }
2030 
2031  mexp[0]=-Econst/xn[0];
2032  mexp[1]=Econst*pow(Cc1,2);
2033  mexp[2]=-Econst*Cc1;
2034  mexp[3]=Econst*pow(Cc2,2);
2035  mexp[4]=-Econst*Cc2;
2036  mexp[5]=Econst*pow(Cc3,2);
2037  mexp[6]=-Econst*Cc3;
2038 
2039  mexp[7]=-Econst*cx[j]*cy[j];
2040  mexp[8]=-Econst*cy[j]*cz[j];
2041  mexp[9]=-Econst*cz[j]*cx[j];
2042 
2043  for (int row=0;row<10;row++){
2044  for (int col=0;col<10;col++){
2045  fjac(row,col)+=malphaESBGK(j,row)*mexp[col]; //new
2046  }
2047  }
2048 
2049 
2050 
2051 
2052  }
2053 
2054 
2055  }
2056 
2057  void initializeMaxwellian()
2058  {
2059  const int numMeshes = _meshes.size();
2060  for (int n=0; n<numMeshes; n++)
2061  {
2062  const Mesh& mesh = *_meshes[n];
2063  const StorageSite& cells = mesh.getCells();
2064  const int nCells = cells.getCount();
2065 
2066  const double pi=_options.pi;
2067  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
2068  const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
2069  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
2070 
2071  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2072  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2073  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2074  const int numFields= _quadrature.getDirCount();
2075 
2076  for(int j=0;j< numFields;j++){
2077  Field& fnd = *_dsfPtr.dsf[j];
2078  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
2079  for(int c=0; c<nCells;c++){
2080  f[c]=density[c]/pow((pi*temperature[c]),1.5)*
2081  exp(-(pow((cx[j]-v[c][0]),2.0)+pow((cy[j]-v[c][1]),2.0)+
2082  pow((cz[j]-v[c][2]),2.0))/temperature[c]);
2083 
2084 
2085  }
2086 
2087  if (_options.transient)
2088  //updateTime();
2089 
2090  {
2091  Field& fnd1 = *_dsfPtr1.dsf[j];
2092  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
2093  for (int c=0;c<nCells;c++)
2094  f1[c] = f[c];
2095  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
2096  if (_options.timeDiscretizationOrder > 1)
2097  {
2098  Field& fnd2 = *_dsfPtr2.dsf[j];
2099  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
2100  for (int c=0;c<nCells;c++)
2101  f2[c] = f[c];
2102  }
2103  }
2104 
2105 
2106  }
2107  }
2108  }
2109 
2110  void initializeFineMaxwellian()
2111  {
2112  const int numMeshes = _meshes.size();
2113  const T zero(0.);
2114  for (int n=0; n<numMeshes; n++)
2115  {
2116  const Mesh& mesh = *_meshes[n];
2117  const StorageSite& cells = mesh.getCells();
2118  const int nCells = cells.getCount();
2119 
2120  const double pi=_options.pi;
2121  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
2122  const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
2123  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
2124 
2125  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2126  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2127  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2128  const int numFields= _quadrature.getDirCount();
2129 
2130  for(int j=0;j< numFields;j++){
2131  Field& fnd0 = *_dsfPtr0.dsf[j];
2132  Field& fndRes = *_dsfPtrRes.dsf[j];
2133  TArray& f0 = dynamic_cast< TArray&>(fnd0[cells]);
2134  TArray& fRes = dynamic_cast< TArray&>(fndRes[cells]);
2135  for(int c=0; c<nCells;c++){
2136  f0[c]=density[c]/pow((pi*temperature[c]),1.5)*
2137  exp(-(pow((cx[j]-v[c][0]),2.0)+pow((cy[j]-v[c][1]),2.0)+
2138  pow((cz[j]-v[c][2]),2.0))/temperature[c]);
2139  fRes[c]=zero;
2140 
2141  }
2142  }
2143  }
2144  }
2145 
2146  void initializeCoarseMaxwellian()
2147  {
2148  const int numMeshes = _meshes.size();
2149  const T zero(0.);
2150  for (int n=0; n<numMeshes; n++)
2151  {
2152  const Mesh& mesh = *_meshes[n];
2153  const StorageSite& cells = mesh.getCells();
2154  const int nCells = cells.getCount();
2155  const int numFields= _quadrature.getDirCount();
2156 
2157  for(int j=0;j< numFields;j++){
2158  Field& fnd0 = *_dsfPtr0.dsf[j];
2159  Field& fndFAS = *_dsfPtrFAS.dsf[j];
2160  Field& fndRes = *_dsfPtrRes.dsf[j];
2161  TArray& f0 = dynamic_cast< TArray&>(fnd0[cells]);
2162  TArray& fFAS = dynamic_cast< TArray&>(fndFAS[cells]);
2163  TArray& fRes = dynamic_cast< TArray&>(fndRes[cells]);
2164  for(int c=0; c<nCells;c++){
2165  f0[c]=zero;
2166  fFAS[c]=zero;
2167  fRes[c]=zero;
2168 
2169  }
2170  }
2171  }
2172  }
2173 
2174  void weightedMaxwellian(double weight1,double uvel1,double vvel1,double wvel1,double uvel2,double vvel2,double wvel2,double temp1,double temp2)
2175  {
2176  const int numMeshes = _meshes.size();
2177  for (int n=0; n<numMeshes; n++)
2178  {
2179  const Mesh& mesh = *_meshes[n];
2180  const StorageSite& cells = mesh.getCells();
2181  const int nCells = cells.getCount();
2182  //double pi(acos(-1.0));
2183  const double pi=_options.pi;
2184  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2185  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2186  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2187  const int numFields= _quadrature.getDirCount();
2188 
2189  for(int j=0;j< numFields;j++){
2190  Field& fnd = *_dsfPtr.dsf[j];
2191  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
2192  for(int c=0; c<nCells;c++){
2193  f[c]=weight1*1.0/pow((pi*1.0),1.5)*exp(-(pow((cx[j]-uvel1),2.0)+pow((cy[j]-vvel1),2.0)+pow((cz[j]-wvel1),2.0))/temp1)
2194  +(1-weight1)*1.0/pow((pi*1.0),1.5)*exp(-(pow((cx[j]-uvel2),2.0)+pow((cy[j]-vvel2),2.0)+pow((cz[j]-wvel2),2.0))/temp2);
2195  }
2196  if (_options.transient)
2197  {
2198  Field& fnd1 = *_dsfPtr1.dsf[j];
2199  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
2200  for (int c=0;c<nCells;c++)
2201  f1[c] = f[c];
2202  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
2203  if (_options.timeDiscretizationOrder > 1)
2204  {
2205  Field& fnd2 = *_dsfPtr2.dsf[j];
2206  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
2207  for (int c=0;c<nCells;c++)
2208  f2[c] = f[c];
2209  }
2210  }
2211  }
2212  }
2213  }
2214  void weightedMaxwellian(double weight1,double uvel1,double uvel2,double temp1,double temp2)
2215  {
2216  const double vvel1=0.0;
2217  const double wvel1=0.0;
2218  const double vvel2=0.0;
2219  const double wvel2=0.0;
2220  const int numMeshes = _meshes.size();
2221  for (int n=0; n<numMeshes; n++)
2222  {
2223  const Mesh& mesh = *_meshes[n];
2224  const StorageSite& cells = mesh.getCells();
2225  const int nCells = cells.getCount();
2226  //double pi(acos(-1.0));
2227  const double pi=_options.pi;
2228  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2229  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2230  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2231  const int numFields= _quadrature.getDirCount();
2232 
2233  for(int j=0;j< numFields;j++){
2234  Field& fnd = *_dsfPtr.dsf[j];
2235  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
2236  for(int c=0; c<nCells;c++){
2237  f[c]=weight1*1.0/pow((pi*1.0),1.5)*exp(-(pow((cx[j]-uvel1),2.0)+pow((cy[j]-vvel1),2.0)+pow((cz[j]-wvel1),2.0))/temp1)
2238  +(1-weight1)*1.0/pow((pi*1.0),1.5)*exp(-(pow((cx[j]-uvel2),2.0)+pow((cy[j]-vvel2),2.0)+pow((cz[j]-wvel2),2.0))/temp2);
2239  }
2240  if (_options.transient)
2241  {
2242  Field& fnd1 = *_dsfPtr1.dsf[j];
2243  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
2244  for (int c=0;c<nCells;c++)
2245  f1[c] = f[c];
2246  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
2247  if (_options.timeDiscretizationOrder > 1)
2248  {
2249  Field& fnd2 = *_dsfPtr2.dsf[j];
2250  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
2251  for (int c=0;c<nCells;c++)
2252  f2[c] = f[c];
2253  }
2254  }
2255  }
2256  }
2257  }
2258  COMETBCMap& getBCMap() {return _bcMap;}
2259  COMETVCMap& getVCMap() {return _vcMap;}
2260 
2261  COMETModelOptions<T>& getOptions() {return _options;}
2262  const map<int, vector<int> >& getFaceReflectionArrayMap() const { return _faceReflectionArrayMap;}
2263 
2264 
2265 
2266  // const vector<int>& vecReflection = _faceReflectionArrayMap[faceID]
2267 map<string,shared_ptr<ArrayBase> >&
2268  getPersistenceData()
2269  {
2270  _persistenceData.clear();
2271 
2272  Array<int>* niterArray = new Array<int>(1);
2273  (*niterArray)[0] = _niters;
2274  _persistenceData["niters"]=shared_ptr<ArrayBase>(niterArray);
2275 
2276  if (_initialKmodelNorm)
2277  {
2278  // _persistenceData["initialKmodelNorm"] =_initialKmodelNorm->getArrayPtr(_macroFields.pressure);
2279  const Field& dsfField = *_dsfPtr.dsf[0];
2280  _persistenceData["initialKmodelNorm"] =_initialKmodelNorm->getArrayPtr(dsfField);
2281 
2282  }
2283  else
2284  {
2285  Array<T>* xArray = new Array<T>(1);
2286  xArray->zero();
2287  _persistenceData["initialKmodelNorm"]=shared_ptr<ArrayBase>(xArray);
2288 
2289  }
2290  return _persistenceData;
2291  }
2292 
2293  void restart()
2294  {
2295  if (_persistenceData.find("niters") != _persistenceData.end())
2296  {
2297  shared_ptr<ArrayBase> rp = _persistenceData["niters"];
2298  ArrayBase& r = *rp;
2299  Array<int>& niterArray = dynamic_cast<Array<int>& >(r);
2300  _niters = niterArray[0];
2301  }
2302 
2303  if (_persistenceData.find("initialKmodelNorm") != _persistenceData.end())
2304  {
2305  shared_ptr<ArrayBase> r = _persistenceData["initialKmodelNorm"];
2306  _initialKmodelNorm = MFRPtr(new MultiFieldReduction());
2307  Field& dsfField = *_dsfPtr.dsf[0];
2308  _initialKmodelNorm->addArray(dsfField,r);
2309  //_initialKmodelNorm->addArray(_dsfPtr,r);
2310  }
2311  }
2312 
2313  void SetBoundaryConditions()
2314  {
2315  const int numMeshes = _meshes.size();
2316  for (int n=0; n<numMeshes; n++)
2317  {
2318  const Mesh& mesh = *_meshes[n];
2319 
2320  COMETVC<T> *vc(new COMETVC<T>());
2321  vc->vcType = "flow";
2322  _vcMap[mesh.getID()] = vc;
2323  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2324  {
2325  const FaceGroup& fg = *fgPtr;
2326  if (_bcMap.find(fg.id) == _bcMap.end())
2327  {
2328  COMETBC<T> *bc(new COMETBC<T>());
2329 
2330  _bcMap[fg.id] = bc;
2331  if((fg.groupType == "wall"))
2332  {
2333  bc->bcType = "WallBC";
2334  }
2335  else if((fg.groupType == "realwall"))
2336  {
2337  bc->bcType = "RealWallBC";
2338  }
2339  else if (fg.groupType == "velocity-inlet")
2340  {
2341  bc->bcType = "VelocityInletBC";
2342  }
2343  else if (fg.groupType == "pressure-inlet")
2344  {
2345  bc->bcType = "PressureInletBC";
2346  }
2347  else if (fg.groupType == "pressure-outlet")
2348  {
2349  bc->bcType = "PressureOutletBC";
2350  }
2351  else if ((fg.groupType == "symmetry"))
2352  {
2353  bc->bcType = "SymmetryBC";
2354  }
2355 
2356  else if((fg.groupType =="zero-gradient "))
2357  {
2358  bc->bcType = "ZeroGradBC";
2359  }
2360  else
2361  throw CException("COMETModel: unknown face group type "
2362  + fg.groupType);
2363  }
2364  }
2365  /*
2366  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2367  {
2368  const FaceGroup& fg = *fgPtr;
2369  if (_bcMap.find(fg.id) == _bcMap.end())
2370  { COMETBC<T> *bc(new COMETBC<T>());
2371 
2372  _bcMap[fg.id] = bc;
2373 
2374  if ((fg.groupType == "NSinterface"))
2375  {
2376  bc->bcType = "NSInterfaceBC";
2377  }
2378  }
2379  }
2380  */
2381  }
2382  }
2383 
2384  void updateTime()
2385  {
2386  const int numMeshes = _meshes.size();
2387  for (int n=0;n<numMeshes;n++)
2388  {
2389  const Mesh& mesh = *_meshes[n];
2390 
2391  const StorageSite& cells = mesh.getCells();
2392  const int numFields= _quadrature.getDirCount();
2393  //callBoundaryConditions(); //new
2394  for (int direction = 0; direction < numFields; direction++)
2395  {
2396  Field& fnd = *_dsfPtr.dsf[direction];
2397  Field& fndN1 = *_dsfPtr1.dsf[direction];
2398  TArray& f = dynamic_cast<TArray&>(fnd[cells]);
2399  TArray& fN1 = dynamic_cast<TArray&>(fndN1[cells]);
2400  if (_options.timeDiscretizationOrder > 1)
2401  {
2402  Field& fndN2 = *_dsfPtr2.dsf[direction];
2403  TArray& fN2 = dynamic_cast<TArray&>(fndN2[cells]);
2404  fN2 = fN1;
2405  }
2406  fN1 = f;
2407  }
2408 #ifdef FVM_PARALLEL
2409  if ( MPI::COMM_WORLD.Get_rank() == 0 )
2410  {cout << "updated time" <<endl;}
2411 
2412 #endif
2413 #ifndef FVM_PARALLEL
2414  cout << "updated time" <<endl;
2415 #endif
2416  //ComputeMacroparameters(); //update macroparameters
2417  //ComputeCollisionfrequency();
2418  //if (_options.fgamma==0){initializeMaxwellianEq();}
2419  //else{ EquilibriumDistributionBGK();}
2420  //if (_options.fgamma==2){EquilibriumDistributionESBGK();}
2421 
2422  }
2423  }
2424 
2425 
2426  void callCOMETBoundaryConditions()
2427  {
2428  const int numMeshes = _meshes.size();
2429  for (int n=0; n<numMeshes; n++)
2430  {
2431  const Mesh& mesh = *_meshes[n];
2432 
2433  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2434  {
2435  const FaceGroup& fg = *fgPtr;
2436  const StorageSite& faces = fg.site;
2437  //const int nFaces = faces.getCount();
2438  const COMETBC<T>& bc = *_bcMap[fg.id];
2439  COMETBoundaryConditions<T,T,T> cbc(faces, mesh,_geomFields,_quadrature,_macroFields,_dsfPtr);
2440  FloatValEvaluator<VectorT3> bVelocity(bc.getVal("specifiedXVelocity"),
2441  bc.getVal("specifiedYVelocity"),
2442  bc.getVal("specifiedZVelocity"),
2443  faces);
2444  FloatValEvaluator<T> accomCoeff(bc.getVal("accommodationCoefficient"),faces);
2445  FloatValEvaluator<T> bTemperature(bc.getVal("specifiedTemperature"),faces);
2446  FloatValEvaluator<T> bPressure(bc.getVal("specifiedPressure"),faces);
2447  if(bc.bcType=="PressureInletBC")
2448  {
2449  cbc.applyPressureInletBC(bTemperature,bPressure);
2450  }
2451  else if(bc.bcType=="PressureOutletBC")
2452  {
2453  cbc.applyPressureOutletBC(bTemperature,bPressure);
2454  }
2455  else if (bc.bcType == "RealWallBC")
2456  {
2457  //kbc.applyRealWallBC(bVelocity,bTemperature,accomCoeff);
2458  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(fg.id);
2459  const vector<int>& vecReflection=(*pos).second;
2460  cbc.applyRealWallBC(bVelocity,bTemperature,accomCoeff,vecReflection);
2461  }
2462  /*
2463  else if(bc.bcType=="SymmetryBC")
2464  {
2465  //kbc.applySpecularWallBC(); //old boundary works only for cartesian-type quadrature
2466  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(fg.id);
2467  const vector<int>& vecReflection=(*pos).second;
2468  cbc.applySpecularWallBC(vecReflection);
2469  }
2470  */
2471  else if(bc.bcType=="ZeroGradBC")
2472  {
2473  cbc.applyZeroGradientBC();
2474  }
2475  }
2476  foreach(const FaceGroupPtr igPtr, mesh.getInterfaceGroups())
2477  {
2478 
2479  const FaceGroup& ig = *igPtr;
2480  const StorageSite& faces = ig.site;
2481  //const int nFaces = faces.getCount();
2482 
2483 
2484  COMETBoundaryConditions<T,T,T> cbc(faces, mesh,_geomFields,_quadrature,_macroFields,_dsfPtr);
2485  if(ig.groupType=="NSinterface")
2486  {
2487  cbc.applyNSInterfaceBC();//bTemperature,bPressure,bVelocity,bStress);
2488  }
2489  }
2490  }//end of loop through meshes
2491  }
2492 
2493  void OutputDsfBLOCK(const char* filename)
2494  {
2495  FILE * pFile;
2496  pFile = fopen(filename,"w");
2497  int N1=_quadrature.getNVCount();
2498  int N2=_quadrature.getNthetaCount();
2499  int N3=_quadrature.getNphiCount();
2500  fprintf(pFile,"%s \n", "VARIABLES= cx, cy, cz, f,fEq,fES");
2501  fprintf(pFile, "%s %i %s %i %s %i \n","ZONE I=", N3,",J=",N2,",K=",N1);
2502  fprintf(pFile, "%s \n","F=BLOCK, VARLOCATION=(NODAL,NODAL,NODAL,NODAL,NODAL,NODAL)");
2503  const int numMeshes = _meshes.size();
2504  const int cellno=_options.printCellNumber;
2505  for (int n=0; n<numMeshes; n++){
2506  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2507  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2508  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2509  const int numFields= _quadrature.getDirCount();
2510  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cx[j]);}
2511  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cy[j]);}
2512  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cz[j]);}
2513 
2514  const Mesh& mesh = *_meshes[n];
2515  const StorageSite& cells = mesh.getCells();
2516  for(int j=0;j< numFields;j++){
2517  Field& fnd = *_dsfPtr.dsf[j];
2518  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
2519  fprintf(pFile,"%E\n",f[cellno]);
2520  }
2521  for(int j=0;j< numFields;j++){
2522  Field& fEqnd = *_dsfEqPtr.dsf[j];
2523  TArray& fEq = dynamic_cast< TArray&>(fEqnd[cells]);
2524  fprintf(pFile,"%E\n",fEq[cellno]);
2525  }
2526  if(_options.fgamma==2){
2527  for(int j=0;j< numFields;j++){
2528 
2529  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2530  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
2531  fprintf(pFile,"%E\n",fEqES[cellno]);
2532  }}
2533  else{
2534  for(int j=0;j< numFields;j++){
2535 
2536  Field& fndEq = *_dsfEqPtr.dsf[j];
2537  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
2538  fprintf(pFile,"%E\n",fEq[cellno]);
2539  }}
2540 
2541  }
2542  fclose(pFile);
2543  }
2544 
2545  void MakeCoarseMesh1(const MeshList& inMeshes, GeomFields& inGeomFields,
2546  MeshList& outMeshes)
2547  {
2548 
2549  int smallestMesh=-1;
2550  const int numMeshes=inMeshes.size();
2551 
2552  for(int n=0;n<numMeshes;n++)
2553  {
2554  const Mesh& mesh=*inMeshes[n];
2555  const int dim=mesh.getDimension();
2556  Mesh* newMeshPtr=new Mesh(dim);
2557 
2558  outMeshes.push_back(newMeshPtr);
2559 
2560  const StorageSite& inCells=mesh.getCells();
2561  StorageSite& outCells=newMeshPtr->getCells();
2562  StorageSite& outFaces=newMeshPtr->getFaces();
2563  const StorageSite& inFaces=mesh.getFaces();
2564  const int inCellCount=inCells.getSelfCount();
2565  const int inCellTotal=inCells.getCount();
2566  const int inFaceCount=inFaces.getCount();
2567  const int inGhost=inCellTotal-inCellCount;
2568  int coarseCount=0;
2569  _siteMap[&inCells]=&outCells;
2570  Field& FineToCoarseField=inGeomFields.fineToCoarse;
2571  IntArray& FineToCoarse=dynamic_cast<IntArray&>(FineToCoarseField[inCells]);
2572 
2573  const CRConnectivity& inCellinFaces=mesh.getCellFaces();
2574  const CRConnectivity& inFaceinCells=mesh.getFaceCells(inFaces);
2575  Field& areaMagField=inGeomFields.areaMag;
2576  Field& areaField=inGeomFields.area;
2577  const TArray& areaMagArray=dynamic_cast<const TArray&>(areaMagField[inFaces]);
2578  const VectorT3Array& areaArray=dynamic_cast<const VectorT3Array&>(areaField[inFaces]);
2579  const BCfaceArray& inBCfArray=*(_BFaces[n]);
2580 
2581  const IntArray& ibType = dynamic_cast<const IntArray&>(inGeomFields.ibType[inCells]);
2582 
2583  //first sweep to make initial pairing
2584  int pairWith;
2585  Array<bool> marker(inFaceCount);
2586  Array<bool> marked(inCellCount);
2587  const T zero(0.);
2588  for(int c=0;c<inCellCount;c++)
2589  {
2590  if((FineToCoarse[c]<0)&&(ibType[c] == Mesh::IBTYPE_FLUID)) //dont bother if im already paired
2591  {
2592  //loop through all neighbors to find pairing
2593  const int neibCount=inCellinFaces.getCount(c);
2594  pairWith=-1;
2595  T maxArea=0.;
2596  int c2;
2597  for(int i=0; i<inFaceCount; i++)
2598  {
2599  marker[i] = false;
2600  }
2601  for(int i=0; i<inCellCount; i++)
2602  {
2603  marked[i] = false;
2604  }
2605  for(int neib=0;neib<neibCount;neib++)
2606  {
2607  const int f=inCellinFaces(c,neib);
2608 
2609  if(inBCfArray[f]==0) //not a boundary face
2610  {
2611  if(c==inFaceinCells(f,1))
2612  c2=inFaceinCells(f,0);
2613  else
2614  c2=inFaceinCells(f,1);
2615 
2616  VectorT3 tempArea;
2617  tempArea[0]=zero;
2618  tempArea[1]=zero;
2619  tempArea[2]=zero;
2620  if((FineToCoarse[c2]==-1)&&(!marked[c2]))
2621  {
2622  marker[f]=true;
2623  marked[c2]=true;
2624  for(int face=0;face<inFaceCount;face++)
2625  {
2626  if((c==inFaceinCells(face,0))&&(c2==inFaceinCells(face,1)))
2627  tempArea+=areaArray[face];
2628  else if((c2==inFaceinCells(face,0))&&(c==inFaceinCells(face,1)))
2629  tempArea-=areaArray[face];
2630  }
2631  }
2632 
2633  if((FineToCoarse[c2]==-1)&&(marker[f]))
2634  if(mag(tempArea)>maxArea)
2635  {
2636  pairWith=c2;
2637  maxArea=mag(tempArea);
2638  }
2639 
2640  /*
2641  if(FineToCoarse[c2]==-1)
2642  if(areaMagArray[f]>maxArea)
2643  {
2644  pairWith=c2;
2645  maxArea=areaMagArray[f];
2646  }
2647  */
2648  }
2649  }
2650 
2651  if(pairWith!=-1)
2652  {
2653  FineToCoarse[c]=coarseCount;
2654  FineToCoarse[pairWith]=coarseCount;
2655  coarseCount++;
2656  }
2657  }
2658  }
2659 
2660  //second sweep to group stragglers, or group with self
2661  for(int c=0;c<inCellCount;c++)
2662  {
2663  if((FineToCoarse[c]==-1)&&(ibType[c] == Mesh::IBTYPE_FLUID))
2664  {
2665  const int neibCount=inCellinFaces.getCount(c);
2666  T maxArea=0.;
2667  int c2,c2perm;
2668  pairWith=-2;
2669  for(int i=0; i<inFaceCount; i++)
2670  {
2671  marker[i] = false;
2672  }
2673  for(int i=0; i<inCellCount; i++)
2674  {
2675  marked[i] = false;
2676  }
2677 
2678  for(int neib=0;neib<neibCount;neib++)
2679  {
2680  const int f=inCellinFaces(c,neib);
2681 
2682  if(inBCfArray[f]==0) //not a boundary face
2683  {
2684  if(c==inFaceinCells(f,1))
2685  c2=inFaceinCells(f,0);
2686  else
2687  c2=inFaceinCells(f,1);
2688 
2689  VectorT3 tempArea;
2690  tempArea[0]=zero;
2691  tempArea[1]=zero;
2692  tempArea[2]=zero;
2693  if(!marked[c2])
2694  {
2695  marker[f]=true;
2696  marked[c2]=true;
2697  for(int face=0;face<inFaceCount;face++)
2698  {
2699  if((c==inFaceinCells(face,0))&&(c2==inFaceinCells(face,1)))
2700  tempArea+=areaArray[face];
2701  else if((c2==inFaceinCells(face,0))&&(c==inFaceinCells(face,1)))
2702  tempArea-=areaArray[face];
2703  }
2704  }
2705 
2706  if(marker[f])
2707  if(mag(tempArea)>maxArea)
2708  {
2709  pairWith=FineToCoarse[c2]; //coarse level cell
2710  c2perm=c2; //fine level cell
2711  maxArea=mag(tempArea);
2712  }
2713 
2714  /*
2715  if(areaMagArray[f]>maxArea)
2716  {
2717  pairWith=FineToCoarse[c2]; //coarse level cell
2718  c2perm=c2; //fine level cell
2719  maxArea=areaMagArray[f];
2720  }
2721  */
2722  }
2723  }
2724 
2725  if(pairWith==-2)
2726  {
2727  FineToCoarse[c]=coarseCount;
2728  coarseCount++;
2729  }
2730  else
2731  {
2732  if(FineToCoarse[c2perm]==-1)
2733  {
2734  FineToCoarse[c]=coarseCount;
2735  FineToCoarse[c2perm]=coarseCount;
2736  coarseCount++;
2737  }
2738  else
2739  FineToCoarse[c]=pairWith;
2740  }
2741  }
2742  }
2743 
2744  for(int c=0;c<inCellCount;c++)
2745  {
2746  if(ibType[c] != Mesh::IBTYPE_FLUID)
2747  {
2748  FineToCoarse[c]=coarseCount;
2749  coarseCount++;
2750  }
2751  }
2752 
2753  int coarseGhost=coarseCount;
2754  _coarseSizes[&mesh]=coarseCount;
2755 
2756  /*
2757  if(MPI::COMM_WORLD.Get_rank()==1)
2758  for(int c=0;c<inCells.getCount();c++)
2759  cout<<" in makecoarsemesh1 before boundary, rank,level,cell,index = "<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<c<<" "<<FineToCoarse[c]<<endl;
2760  */
2761 
2762  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2763  {
2764  const FaceGroup& fg = *fgPtr;
2765  const StorageSite& faces = fg.site;
2766  const int nFaces = faces.getCount();
2767 
2768  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
2769 
2770  for(int f=0; f< nFaces; f++)
2771  {
2772  const int c1= faceCells(f,1);// boundary cell
2773  FineToCoarse[c1]=coarseGhost;
2774  coarseGhost++;
2775  }
2776  }
2777 
2778  /*
2779  if(MPI::COMM_WORLD.Get_rank()==1)
2780  for(int c=0;c<inCells.getCount();c++)
2781  cout<<" in makecoarsemesh1, rank, level,cell,index = "<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<c<<" "<<FineToCoarse[c]<<endl;
2782  */
2783  }
2784  }
2785 
2786  int MakeCoarseMesh2(const MeshList& inMeshes, GeomFields& inGeomFields,
2787  GeomFields& coarseGeomFields,MeshList& outMeshes)
2788  {
2789 
2790  int smallestMesh=-1;
2791  const int numMeshes=inMeshes.size();
2792  for(int n=0;n<numMeshes;n++)
2793  {
2794  const Mesh& mesh=*inMeshes[n];
2795  const int dim=mesh.getDimension();
2796  //Mesh* newMeshPtr=new Mesh(dim);
2797 
2798  //outMeshes.push_back(newMeshPtr);
2799  Mesh& newMeshPtr=*outMeshes[n];
2800 
2801  const StorageSite& inCells=mesh.getCells();
2802  StorageSite& outCells=newMeshPtr.getCells();
2803  StorageSite& outFaces=newMeshPtr.getFaces();
2804  StorageSite& outIBFaces=newMeshPtr.getIBFaces();
2805  const IntArray& ibType = dynamic_cast<const IntArray&>(inGeomFields.ibType[inCells]);
2806  const StorageSite& inFaces=mesh.getFaces();
2807  const StorageSite& inIBFaces=mesh.getIBFaces();
2808  const int inCellCount=inCells.getSelfCount();
2809  const int inCellTotal=inCells.getCount();
2810  const int inFaceCount=inFaces.getCount();
2811  const int inGhost=inCellTotal-inCellCount;
2812  int outGhost = 0;
2813  int coarseCount=0;
2814  Field& FineToCoarseField=inGeomFields.fineToCoarse;
2815  const IntArray& FineToCoarse=dynamic_cast<const IntArray&>(FineToCoarseField[inCells]);
2816  const CRConnectivity& inCellinFaces=mesh.getCellFaces();
2817  const CRConnectivity& inFaceinCells=mesh.getFaceCells(inFaces);
2818  Field& areaMagField=inGeomFields.areaMag;
2819  const TArray& areaMagArray=dynamic_cast<const TArray&>(areaMagField[inFaces]);
2820  const BCfaceArray& inBCfArray=*(_BFaces[n]);
2821 
2822  /*** creating coarse level cells ***/
2823  coarseCount = _coarseSizes[&mesh];
2824  int interfaceCellsLevel0 = _coarseGhostSizes[&mesh];
2825 
2826  int boundaryCell=0;
2827  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2828  {
2829  const FaceGroup& fg = *fgPtr;
2830  const StorageSite& faces = fg.site;
2831  const int nFaces = faces.getCount();
2832 
2833  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
2834 
2835  for(int f=0; f< nFaces; f++)
2836  {
2837  const int c1= faceCells(f,1);// boundary cell
2838  if(boundaryCell<FineToCoarse[c1])
2839  boundaryCell=FineToCoarse[c1];
2840  }
2841  }
2842  boundaryCell++;
2843  if(boundaryCell!=1)
2844  boundaryCell-=coarseCount;
2845  else
2846  boundaryCell=0;
2847 
2848  for(int c=0; c< inCells.getCountLevel1(); c++)
2849  {
2850  if(outGhost<FineToCoarse[c])
2851  outGhost=FineToCoarse[c];
2852  }
2853  outGhost++;
2854  outGhost-=coarseCount;
2855 
2856  int interfaceCells = outGhost - boundaryCell;
2857 
2858  /*
2859  if(MPI::COMM_WORLD.Get_rank()==1)
2860  cout<<"rank,level,inCellinternal, incelltotal,outCellInternal, outCellExternal, outInterface ="<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<inCellCount<<" and "<<inCellTotal<<" and "<<coarseCount<<" and "<<outGhost<<" "<<interfaceCells<<endl;
2861  */
2862 
2863  /*
2864  if(MPI::COMM_WORLD.Get_rank()==1)
2865  for(int c=0;c<inCellTotal;c++)
2866  cout<<" in makecoarsemesh2, rank,level, cell,index = "<<MPI::COMM_WORLD.Get_rank()<<" "<<_level<<" "<<c<<" and "<<FineToCoarse[c]<<endl;
2867  */
2868 
2869  outCells.setCount(coarseCount,boundaryCell+interfaceCellsLevel0);
2870  outCells.setCountLevel1(coarseCount+outGhost);
2871 
2872  /*** created coarse level cells ***/
2873 
2874  //make the coarse cell to fine cell connectivity.
2875  CRPtr CoarseToFineCells=CRPtr(new CRConnectivity(outCells,inCells));
2876  CoarseToFineCells->initCount();
2877 
2878  for(int c=0;c<inCells.getCountLevel1();c++)
2879  CoarseToFineCells->addCount(FineToCoarse[c],1);
2880 
2881  CoarseToFineCells->finishCount();
2882 
2883  for(int c=0;c<inCells.getCountLevel1();c++)
2884  CoarseToFineCells->add(FineToCoarse[c],c);
2885 
2886  CoarseToFineCells->finishAdd();
2887 
2888  /*** connectivity between itself (cells) and its finer mesh cells ***/
2889  newMeshPtr.setConnectivity(outCells,inCells,CoarseToFineCells);
2890 
2891  CRPtr FineFacesCoarseCells=CRPtr(new CRConnectivity(inFaces,outCells));
2892  FineFacesCoarseCells->initCount();
2893 
2894  //count surviving faces
2895  int survivingFaces=0;
2896  int coarse0, coarse1;
2897  for(int f=0;f<inFaceCount;f++)
2898  {
2899  coarse0=FineToCoarse[inFaceinCells(f,0)];
2900  coarse1=FineToCoarse[inFaceinCells(f,1)];
2901  if(coarse0!=coarse1)
2902  {
2903  survivingFaces++;
2904  FineFacesCoarseCells->addCount(f,2);
2905  }
2906  }
2907 
2908  FineFacesCoarseCells->finishCount();
2909 
2910  //make non-zero's
2911  int fc0,fc1,cc0,cc1;
2912  for(int f=0;f<inFaceCount;f++)
2913  {
2914  fc0=inFaceinCells(f,0);
2915  fc1=inFaceinCells(f,1);
2916  cc0=FineToCoarse[fc0];
2917  cc1=FineToCoarse[fc1];
2918  if(cc0!=cc1)
2919  {
2920  FineFacesCoarseCells->add(f,cc0);
2921  FineFacesCoarseCells->add(f,cc1);
2922  }
2923  }
2924 
2925  FineFacesCoarseCells->finishAdd();
2926 
2927  CRPtr CoarseCellsFineFaces=FineFacesCoarseCells->getTranspose();
2928  CRPtr CellCellCoarse=CoarseCellsFineFaces->multiply(*FineFacesCoarseCells,true);
2929 
2930  /*** coarse level cellcell connectivity created ***/
2931 
2932  /*
2933  int counter=0;
2934  BArray counted(outCells.getCount());
2935  counted=false;
2936  for(int c=0;c<outCells.getCount();c++)
2937  {
2938  counted[c]=true;
2939  const int neibs=CellCellCoarse->getCount(c);
2940  for(int n=0;n<neibs;n++)
2941  {
2942  const int c1=(*CellCellCoarse)(c,n);
2943  if(!counted[c1])
2944  counter++;
2945  }
2946  }
2947 
2948  //outFaces.setCount(counter);
2949  */
2950 
2951  int countFaces=0;
2952  int cCell0, cCell1;
2953  for(int f=0;f<inFaceCount;f++)
2954  {
2955  cCell0=FineToCoarse[inFaceinCells(f,0)];
2956  cCell1=FineToCoarse[inFaceinCells(f,1)];
2957  if(cCell0!=cCell1)
2958  {
2959  countFaces++;
2960  }
2961  }
2962 
2963  outFaces.setCount(countFaces);
2964 
2965  int inFaceGhost = 0;
2966  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2967  inFaceGhost+=(*fgPtr).site.getCount();
2968  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2969  inFaceGhost+=(*fgPtr).site.getCount();
2970  const int inInteriorFaces = inFaceCount - inFaceGhost;
2971 
2972  const int del = inFaceCount - outFaces.getCount();
2973 
2974  //const int interiorCount=outFaces.getCount()-inGhost;
2975  const int interiorCount=inInteriorFaces-del;
2976  //if(MPI::COMM_WORLD.Get_rank()==1)
2977  //cout<<"level,outfaces, outghost, totalfaces = "<<_level<<" "<<interiorCount<<" "<<inGhost<<" "<<countFaces<<endl;
2978  const StorageSite& interiorFaces=newMeshPtr.createInteriorFaceGroup(interiorCount);
2979 
2980  int inOffset=interiorCount;
2981  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2982  {
2983  const FaceGroup& fg=*fgPtr;
2984  const int size=fg.site.getCount();
2985  newMeshPtr.createBoundaryFaceGroup(size,inOffset,fg.id,fg.groupType);
2986  inOffset+=size;
2987  }
2988 
2989  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2990  {
2991  const FaceGroup& fg=*fgPtr;
2992  const int size=fg.site.getCount();
2993  newMeshPtr.createInterfaceGroup(size,inOffset,fg.id);
2994  inOffset+=size;
2995  }
2996 
2997  CRPtr CoarseFaceCoarseCell=CRPtr(new CRConnectivity(outFaces,outCells));
2998  CoarseFaceCoarseCell->initCount();
2999 
3000  survivingFaces=0;
3001  for(int f=0;f<inFaceCount;f++)
3002  {
3003  coarse0=FineToCoarse[inFaceinCells(f,0)];
3004  coarse1=FineToCoarse[inFaceinCells(f,1)];
3005  if(coarse0!=coarse1)
3006  {
3007  CoarseFaceCoarseCell->addCount(survivingFaces,2);
3008  survivingFaces++;
3009  }
3010  }
3011 
3012  CoarseFaceCoarseCell->finishCount();
3013 
3014  //make non-zero's
3015  survivingFaces=0;
3016  for(int f=0;f<inFaceCount;f++)
3017  {
3018  fc0=inFaceinCells(f,0);
3019  fc1=inFaceinCells(f,1);
3020  cc0=FineToCoarse[fc0];
3021  cc1=FineToCoarse[fc1];
3022  if(cc0!=cc1)
3023  {
3024  CoarseFaceCoarseCell->add(survivingFaces,cc0);
3025  CoarseFaceCoarseCell->add(survivingFaces,cc1);
3026  survivingFaces++;
3027  }
3028  }
3029 
3030  CoarseFaceCoarseCell->finishAdd();
3031 
3032  CRPtr CoarseCellCoarseFace=CoarseFaceCoarseCell->getTranspose();
3033 
3034  newMeshPtr.setConnectivity(outCells,outFaces,CoarseCellCoarseFace);
3035  newMeshPtr.setConnectivity(outFaces,outCells,CoarseFaceCoarseCell);
3036 
3037  Field& coarseIbTypeField=coarseGeomFields.ibType;
3038  shared_ptr<IntArray> ibTypePtr(new IntArray(outCells.getCountLevel1()));
3039  *ibTypePtr = Mesh::IBTYPE_FLUID;
3040  coarseIbTypeField.addArray(outCells,ibTypePtr);
3041 
3042  IntArray& coarseIBType = dynamic_cast<IntArray&>(coarseGeomFields.ibType[outCells]);
3043 
3044  for(int c=0;c<inCellCount;c++)
3045  if(ibType[c] != Mesh::IBTYPE_FLUID)
3046  coarseIBType[FineToCoarse[c]]=ibType[c];
3047 
3048  outIBFaces.setCount(inIBFaces.getCount());
3049  _siteMap[&inIBFaces]=&outIBFaces;
3050 
3051  shared_ptr<IntArray> ibFaceIndexPtr(new IntArray(outFaces.getCount()));
3052  *ibFaceIndexPtr = -1;
3053  coarseGeomFields.ibFaceIndex.addArray(outFaces,ibFaceIndexPtr);
3054 
3055  const IntArray& fineIBFaceIndex=dynamic_cast<const IntArray&>(inGeomFields.ibFaceIndex[inFaces]);
3056  IntArray& coarseIBFaceIndex=dynamic_cast<IntArray&>(coarseGeomFields.ibFaceIndex[outFaces]);
3057 
3058  CRPtr CoarseFacesFineFaces=CRPtr(new CRConnectivity(outFaces,inFaces));
3059  CoarseFacesFineFaces->initCount();
3060 
3061  survivingFaces=0;
3062  for(int f=0;f<inFaceCount;f++)
3063  {
3064  int fc0=inFaceinCells(f,0);
3065  int fc1=inFaceinCells(f,1);
3066  const int cc0=FineToCoarse[fc0];
3067  const int cc1=FineToCoarse[fc1];
3068 
3069  if(cc1!=cc0)
3070  {
3071  CoarseFacesFineFaces->addCount(survivingFaces,1);
3072  survivingFaces++;
3073  }
3074  }
3075  //cout<<"hello 3 rank "<<MPI::COMM_WORLD.Get_rank()<<endl;
3076  CoarseFacesFineFaces->finishCount();
3077 
3078  survivingFaces=0;
3079  for(int f=0;f<inFaceCount;f++)
3080  {
3081  int fc0=inFaceinCells(f,0);
3082  int fc1=inFaceinCells(f,1);
3083  const int cc0=FineToCoarse[fc0];
3084  const int cc1=FineToCoarse[fc1];
3085  if(cc1!=cc0)
3086  {
3087  CoarseFacesFineFaces->add(survivingFaces,f);
3088  coarseIBFaceIndex[survivingFaces]=fineIBFaceIndex[f];
3089  survivingFaces++;
3090  }
3091  }
3092  //cout<<"hello 4 rank "<<MPI::COMM_WORLD.Get_rank()<<endl;
3093  CoarseFacesFineFaces->finishAdd();
3094 
3095  /*
3096  if(MPI::COMM_WORLD.Get_rank()==1)
3097  {
3098  for(int f=0;f<outFaces.getCount();f++)
3099  {
3100  cout<<"level,rank,face,neighbor = "<<_level<<" "<<MPI::COMM_WORLD.Get_rank()<<" "<<f<<" "<<(*CoarseFaceCoarseCell)(f,0)<<" "<<(*CoarseFaceCoarseCell)(f,1)<<endl;
3101  }
3102  }
3103  */
3104 
3105  //cout<<"hello 5 rank "<<MPI::COMM_WORLD.Get_rank()<<endl;
3106 
3107  //now make the geom fields
3108  const int outCellsCount=outCells.getSelfCount();
3109  TArrptr outCellVolumePtr=TArrptr(new TArray(outCellsCount));
3110  TArray& outCV=*outCellVolumePtr;
3111  outCV=0.;
3112 
3113  Field& VolumeField=inGeomFields.volume;
3114  Field& coarseVolumeField=coarseGeomFields.volume;
3115  const TArray& inCV=dynamic_cast<const TArray&>(VolumeField[inCells]);
3116 
3117  for(int c=0;c<outCellsCount;c++)
3118  {
3119  const int fineCount=CoarseToFineCells->getCount(c);
3120  for(int i=0;i<fineCount;i++)
3121  {
3122  int fc=(*CoarseToFineCells)(c,i);
3123  outCV[c]+=inCV[fc];
3124  }
3125  }
3126 
3127  coarseVolumeField.addArray(outCells,outCellVolumePtr);
3128 
3129  //cout<<"hello 6"<<endl;
3130 
3131  const int outFacesCount=outFaces.getCount();
3132  VT3Ptr outFaceAreaPtr=VT3Ptr(new VectorT3Array(outFacesCount));
3133  VectorT3Array& outFA=*outFaceAreaPtr;
3134  TArrptr outFaceAreaMagPtr=TArrptr(new TArray(outFacesCount));
3135  TArray& outFAMag=*outFaceAreaMagPtr;
3136 
3137  Field& FaceAreaField=inGeomFields.area;
3138  Field& coarseFaceAreaField=coarseGeomFields.area;
3139  Field& coarseAreaMagField=coarseGeomFields.areaMag;
3140  const VectorT3Array& inFA=
3141  dynamic_cast<const VectorT3Array&>(FaceAreaField[inFaces]);
3142 
3143  VectorT3 myZero;
3144  myZero[0]=0.;
3145  myZero[1]=0.;
3146  myZero[2]=0.;
3147 
3148  outFA=myZero;
3149  outFAMag=0.;
3150  for(int f=0;f<outFacesCount;f++)
3151  {
3152  const int fineCount=CoarseFacesFineFaces->getCount(f);
3153  const int cCell0=(*CoarseFaceCoarseCell)(f,0);
3154  for(int i=0;i<fineCount;i++)
3155  {
3156  const int fFace=(*CoarseFacesFineFaces)(f,i);
3157  const int fCell0=inFaceinCells(fFace,0);
3158  const int CCell0=FineToCoarse[fCell0];
3159 
3160  //must make sure the area vector is pointing
3161  //from c0 to c1
3162  if(CCell0==cCell0)
3163  outFA[f]+=inFA[fFace];
3164  else
3165  outFA[f]-=inFA[fFace];
3166  outFAMag[f]+=areaMagArray[fFace];
3167  }
3168  }
3169 
3170  coarseFaceAreaField.addArray(outFaces,outFaceAreaPtr);
3171  coarseAreaMagField.addArray(outFaces,outFaceAreaMagPtr);
3172 
3173  //cout<<"hello 7"<<endl;
3174  /*
3175  Field& ibTypeField=inGeomFields.ibType;
3176  Field& coarseIbTypeField=coarseGeomFields.ibType;
3177  shared_ptr<IntArray> ibTypePtr(new IntArray(outCells.getSelfCount()));
3178  *ibTypePtr = Mesh::IBTYPE_FLUID;
3179  coarseIbTypeField.addArray(outCells,ibTypePtr);
3180  */
3181  if(smallestMesh<0)
3182  smallestMesh=outCells.getSelfCount();
3183  else
3184  {
3185  if(outCells.getSelfCount()<smallestMesh)
3186  smallestMesh=outCells.getSelfCount();
3187  }
3188  //cout<<"hello 8"<<endl;
3189  /*
3190  //This is for checking purposes only
3191  cout<<"Coarse Faces to Fine Faces"<<endl;
3192  for(int f=0;f<outFaces.getCount();f++)
3193  {
3194  const int neibs=CoarseFacesFineFaces->getCount(f);
3195  for(int n=0;n<neibs;n++)
3196  cout<<f<<" "<<(*CoarseFacesFineFaces)(f,n)<<endl;
3197  cout<<endl;
3198  }
3199  cout<<"Coarse Cells to Coarse Faces"<<endl;
3200  for(int c=0;c<outCells.getCount();c++)
3201  {
3202  const int neibs=CoarseCellCoarseFace->getCount(c);
3203  for(int n=0;n<neibs;n++)
3204  cout<<c<<" "<<(*CoarseCellCoarseFace)(c,n)<<endl;
3205  cout<<endl;
3206  }
3207  */
3208  }
3209  //cout<<"hello 6"<<endl;
3210  return smallestMesh;
3211  }
3212 
3213 
3214  void syncGhostCoarsening(const MeshList& inMeshes, GeomFields& inGeomFields,
3215  MeshList& outMeshes)
3216  {
3217  //const int xLen = coarseIndexField.getLength();
3218 
3219  //#pragma omp parallel for
3220  const int numMeshes=inMeshes.size();
3221  for(int n=0;n<numMeshes;n++)
3222  {
3223  const Mesh& mesh=*inMeshes[n];
3224  const int dim=mesh.getDimension();
3225  Mesh& newMeshPtr=*outMeshes[n];
3226 
3227  const StorageSite& inCells=mesh.getCells();
3228  const StorageSite& site=mesh.getCells();
3229  //const StorageSite& site=newMeshPtr.getCells();
3230  const int inCellCount=inCells.getSelfCount();
3231  const int inCellTotal=inCells.getCount();
3232 
3233  Field& FineToCoarseField=inGeomFields.fineToCoarse;
3234  IntArray& coarseIndex=dynamic_cast<IntArray&>(FineToCoarseField[inCells]);
3235  IntArray tempIndex(inCells.getCountLevel1());
3236  for(int c=0;c<inCells.getCountLevel1();c++)
3237  tempIndex[c]=coarseIndex[c];
3238 
3239  int coarseGhostSize=0;
3240  int tempGhostSize=0;
3241  //const int coarseSize = _coarseSizes.find(rowIndex)->second;
3242  //const int coarseSize = _coarseSizes[&mesh];
3243  int coarseSize = -1;
3244  coarseSize = _coarseSizes[&mesh];
3245 
3246  int boundaryCell=0;
3247  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3248  {
3249  const FaceGroup& fg = *fgPtr;
3250  const StorageSite& faces = fg.site;
3251  const int nFaces = faces.getCount();
3252 
3253  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
3254 
3255  for(int f=0; f< nFaces; f++)
3256  {
3257  const int c1= faceCells(f,1);// boundary cell
3258  //if(MPI::COMM_WORLD.Get_rank()==1)
3259  //cout<<"fgid, boundarycell, boundary coarse index = "<<fg.id<<" "<<c1<<" "<<coarseIndex[c1]<<endl;
3260  if(boundaryCell<coarseIndex[c1])
3261  boundaryCell=coarseIndex[c1];
3262  }
3263  }
3264  boundaryCell++;
3265  if(boundaryCell!=1)
3266  coarseSize = boundaryCell;
3267 
3268  //const StorageSite& site = *rowIndex.second;
3269 
3270  const StorageSite::GatherMap& gatherMap = site.getGatherMap();
3271  const StorageSite::GatherMap& gatherMapLevel1 = site.getGatherMapLevel1();
3272 
3273  // collect all the toIndices arrays for each storage site from
3274  // both gatherMap and gatherMapLevel1
3275 
3276  typedef map<const StorageSite*, vector<const Array<int>* > > IndicesMap;
3277  IndicesMap toIndicesMap;
3278  IndicesMap tempIndicesMap;
3279 
3280  foreach(const StorageSite::GatherMap::value_type pos, gatherMap)
3281  {
3282  const StorageSite& oSite = *pos.first;
3283  const Array<int>& tempIndices = *pos.second;
3284  const Array<int>& toIndices = *pos.second;
3285 
3286  tempIndicesMap[&oSite].push_back(&tempIndices);
3287  toIndicesMap[&oSite].push_back(&toIndices);
3288  }
3289 
3290  foreach(const StorageSite::GatherMap::value_type pos, gatherMapLevel1)
3291  {
3292  const StorageSite& oSite = *pos.first;
3293  const Array<int>& toIndices = *pos.second;
3294 
3295  int found=0;
3296  foreach(const StorageSite::GatherMap::value_type posLevel0, gatherMap)
3297  {
3298  const StorageSite& oSiteLevel0 = *posLevel0.first;
3299  if(oSite.getTag()==oSiteLevel0.getTag())
3300  {
3301  toIndicesMap[&oSiteLevel0].push_back(&toIndices);
3302  found=1;
3303  }
3304  }
3305 
3306  if(found==0)
3307  toIndicesMap[&oSite].push_back(&toIndices);
3308  }
3309 
3310  foreach(IndicesMap::value_type pos, tempIndicesMap)
3311  {
3312  const StorageSite& oSite = *pos.first;
3313  const vector<const Array<int>* > tempIndicesArrays = pos.second;
3314 
3315 
3316  map<int,int> otherToMyMapping;
3317  UnorderedSet gatherSet;
3318 
3319  foreach(const Array<int>* tempIndicesPtr, tempIndicesArrays)
3320  {
3321  const Array<int>& tempIndices = *tempIndicesPtr;
3322  const int nGhostRows = tempIndices.getLength();
3323  for(int ng=0; ng<nGhostRows; ng++)
3324  {
3325  const int fineIndex = tempIndices[ng];
3326  const int coarseOtherIndex = tempIndex[fineIndex];
3327 
3328  if (coarseOtherIndex < 0)
3329  continue;
3330 
3331 
3332  if (otherToMyMapping.find(coarseOtherIndex) !=
3333  otherToMyMapping.end())
3334  {
3335 
3336  tempIndex[fineIndex] = otherToMyMapping[coarseOtherIndex];
3337  }
3338  else
3339  {
3340  tempIndex[fineIndex] = tempGhostSize+coarseSize;
3341  otherToMyMapping[coarseOtherIndex] = tempIndex[fineIndex];
3342  gatherSet.insert( tempIndex[fineIndex] );
3343  tempGhostSize++;
3344  }
3345  }
3346  }
3347  }
3348 
3349  foreach(IndicesMap::value_type pos, toIndicesMap)
3350  {
3351  const StorageSite& oSite = *pos.first;
3352  const vector<const Array<int>* > toIndicesArrays = pos.second;
3353 
3354 
3355  map<int,int> otherToMyMapping;
3356  UnorderedSet gatherSet;
3357 
3358  foreach(const Array<int>* toIndicesPtr, toIndicesArrays)
3359  {
3360  const Array<int>& toIndices = *toIndicesPtr;
3361  const int nGhostRows = toIndices.getLength();
3362  for(int ng=0; ng<nGhostRows; ng++)
3363  {
3364  const int fineIndex = toIndices[ng];
3365  const int coarseOtherIndex = coarseIndex[fineIndex];
3366 
3367  if (coarseOtherIndex < 0)
3368  continue;
3369 
3370  if (otherToMyMapping.find(coarseOtherIndex) !=
3371  otherToMyMapping.end())
3372  {
3373  coarseIndex[fineIndex] = otherToMyMapping[coarseOtherIndex];
3374  }
3375  else
3376  {
3377  coarseIndex[fineIndex] = coarseGhostSize+coarseSize;
3378  otherToMyMapping[coarseOtherIndex] = coarseIndex[fineIndex];
3379  gatherSet.insert( coarseIndex[fineIndex] );
3380  coarseGhostSize++;
3381  }
3382  //if(MPI::COMM_WORLD.Get_rank()==1)
3383  //cout<<"level,fineIndex, coarseIndex = "<<_level<<" "<<fineIndex<<" "<<coarseIndex[fineIndex]<<endl;
3384  }
3385  //if(MPI::COMM_WORLD.Get_rank()==1)
3386  //cout<<endl<<endl;
3387  }
3388 
3389  const int coarseMappersSize = otherToMyMapping.size();
3390 
3391  shared_ptr<Array<int> > coarseToIndices(new Array<int>(coarseMappersSize));
3392 
3393  for(int n = 0; n < gatherSet.size(); n++)
3394  {
3395  (*coarseToIndices)[n] = gatherSet.getData().at(n);
3396  }
3397 
3398  SSPair sskey(&site,&oSite);
3399  _coarseGatherMaps [sskey] = coarseToIndices;
3400  }
3401 
3402  const StorageSite::ScatterMap& scatterMap = site.getScatterMap();
3403  const StorageSite::ScatterMap& scatterMapLevel1 = site.getScatterMapLevel1();
3404 
3405  IndicesMap fromIndicesMap;
3406 
3407  foreach(const StorageSite::GatherMap::value_type pos, scatterMap)
3408  {
3409  const StorageSite& oSite = *pos.first;
3410  const Array<int>& fromIndices = *pos.second;
3411 
3412  fromIndicesMap[&oSite].push_back(&fromIndices);
3413  }
3414 
3415  foreach(const StorageSite::GatherMap::value_type pos, scatterMapLevel1)
3416  {
3417  const StorageSite& oSite = *pos.first;
3418  const Array<int>& fromIndices = *pos.second;
3419 
3420  int found=0;
3421  foreach(const StorageSite::ScatterMap::value_type posLevel0, scatterMap)
3422  {
3423  const StorageSite& oSiteLevel0 = *posLevel0.first;
3424  if(oSite.getTag()==oSiteLevel0.getTag())
3425  {
3426  fromIndicesMap[&oSiteLevel0].push_back(&fromIndices);
3427  found=1;
3428  }
3429  }
3430 
3431  if(found==0)
3432  fromIndicesMap[&oSite].push_back(&fromIndices);
3433  }
3434 
3435  foreach(IndicesMap::value_type pos, fromIndicesMap)
3436  {
3437  const StorageSite& oSite = *pos.first;
3438  const vector<const Array<int>* > fromIndicesArrays = pos.second;
3439 
3440  UnorderedSet scatterSet;
3441 
3442  foreach(const Array<int>* fromIndicesPtr, fromIndicesArrays)
3443  {
3444  const Array<int>& fromIndices = *fromIndicesPtr;
3445  const int nGhostRows = fromIndices.getLength();
3446  for(int ng=0; ng<nGhostRows; ng++)
3447  {
3448  const int fineIndex = fromIndices[ng];
3449  const int coarseOtherIndex = coarseIndex[fineIndex];
3450  if (coarseOtherIndex >= 0)
3451  scatterSet.insert( coarseOtherIndex );
3452  }
3453 
3454  }
3455 
3456  const int coarseMappersSize = scatterSet.size();
3457 
3458  shared_ptr<Array<int> > coarseFromIndices(new Array<int>(coarseMappersSize));
3459 
3460  for(int n = 0; n < scatterSet.size(); n++ ) {
3461  (*coarseFromIndices)[n] = scatterSet.getData().at(n);
3462  }
3463 
3464  SSPair sskey(&site,&oSite);
3465  _coarseScatterMaps[sskey] = coarseFromIndices;
3466 
3467  }
3468  //_coarseGhostSizes[rowIndex]=coarseGhostSize;
3469  _coarseGhostSizes[&mesh]=coarseGhostSize;
3470  }
3471  }
3472 
3473  void computeIBFaceDsf(const StorageSite& solidFaces,const int method,const int RelaxDistribution=0)
3474  {
3475  typedef CRMatrixTranspose<T,T,T> IMatrix;
3476  typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
3477  if (method==1){
3478  const int numMeshes = _meshes.size();
3479  const int numFields= _quadrature.getDirCount();
3480  for (int direction = 0; direction < numFields; direction++)
3481  {
3482  Field& fnd = *_dsfPtr.dsf[direction];
3483  const TArray& pV =
3484  dynamic_cast<const TArray&>(fnd[solidFaces]);
3485  #ifdef FVM_PARALLEL
3486  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,pV.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
3487  #endif
3488 
3489  for (int n=0; n<numMeshes; n++)
3490  {
3491  const Mesh& mesh = *_meshes[n];
3492  const Mesh& fMesh = *_finestMeshes[n];
3493  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3494 
3495  const StorageSite& cells = mesh.getCells();
3496  const StorageSite& fCells = fMesh.getCells();
3497  const StorageSite& ibFaces = mesh.getIBFaces();
3498  const StorageSite& fIbFaces = fMesh.getIBFaces();
3499 
3500  GeomFields::SSPair key1(&fIbFaces,&fCells);
3501  const IMatrix& mIC =
3502  dynamic_cast<const IMatrix&>
3503  (*_finestGeomFields._interpolationMatrices[key1]);
3504 
3505  IMatrix mICV(mIC);
3506 
3507 
3508  GeomFields::SSPair key2(&fIbFaces,&solidFaces);
3509  const IMatrix& mIP =
3510  dynamic_cast<const IMatrix&>
3511  (*_finestGeomFields._interpolationMatrices[key2]);
3512 
3513  IMatrix mIPV(mIP);
3514 
3515 
3516  shared_ptr<TArray> ibV(new TArray(ibFaces.getCount()));
3517 
3518  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
3519  const Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<const Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
3520 
3521  const TArray& cV =
3522  dynamic_cast<const TArray&>(fnd[cells]);
3523  TArray cFV(fCells.getCountLevel1());
3524  for(int c=0;c<fCells.getCountLevel1();c++)
3525  cFV[c]=cV[FinestToCoarse[c][_level]];
3526 
3527  ibV->zero();
3528 
3529  mICV.multiplyAndAdd(*ibV,cFV);
3530  mIPV.multiplyAndAdd(*ibV,pV);
3531 
3532 #if 0
3533  ofstream debugFile;
3534  stringstream ss(stringstream::in | stringstream::out);
3535  ss << MPI::COMM_WORLD.Get_rank();
3536  string fname1 = "IBVelocity_proc" + ss.str() + ".dat";
3537  debugFile.open(fname1.c_str());
3538 
3539  //debug use
3540  const Array<int>& ibFaceList = mesh.getIBFaceList();
3541  const StorageSite& faces = mesh.getFaces();
3542  const VectorT3Array& faceCentroid =
3543  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
3544  const double angV = 1.0;
3545  VectorT3 center;
3546  center[0]=0.;
3547  center[1]=0.;
3548  center[2]=0.;
3549 
3550  for(int f=0; f<ibFaces.getCount();f++){
3551  int fID = ibFaceList[f];
3552  debugFile << "f=" << f << setw(10) << " fID = " << fID << " faceCentroid = " << faceCentroid[fID] << " ibV = " << (*ibV)[f] << endl;
3553  }
3554 
3555  debugFile.close();
3556 #endif
3557 
3558  fnd.addArray(ibFaces,ibV);
3559  }
3560  }
3561  }
3562  /*
3563  for (int n=0; n<numMeshes; n++)
3564  {
3565  const Mesh& mesh = *_meshes[n];
3566  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3567 
3568  const StorageSite& cells = mesh.getCells();
3569  const StorageSite& cells = mesh.getCells();
3570  const StorageSite& ibFaces = mesh.getIBFaces();
3571  const StorageSite& fIbFaces = fMesh.getIBFaces();
3572 
3573  GeomFields::SSPair key1(&fIbFaces,&fCells);
3574  const IMatrix& mIC =
3575  dynamic_cast<const IMatrix&>
3576  (*_finestGeomFields._interpolationMatrices[key1]);
3577 
3578  IMatrixV3 mICV3(mIC);
3579 
3580  GeomFields::SSPair key2(&fIbFaces,&solidFaces);
3581  const IMatrix& mIP =
3582  dynamic_cast<const IMatrix&>
3583  (*_finestGeomFields._interpolationMatrices[key2]);
3584 
3585  IMatrixV3 mIPV3(mIP);
3586 
3587  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
3588 
3589 
3590  const VectorT3Array& cVel =
3591  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
3592  const VectorT3Array& sVel =
3593  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
3594 
3595  ibVvel->zero();
3596 
3597  //velocity interpolation (cells+solidfaces)
3598  mICV3.multiplyAndAdd(*ibVvel,cVel);
3599  mIPV3.multiplyAndAdd(*ibVvel,sVel);
3600  _macroFields.velocity.addArray(ibFaces,ibVvel);
3601 
3602 
3603  }
3604  }
3605  */
3606  }
3607 
3608  if (method==2){
3609  const int numMeshes = _meshes.size();
3610  const int nSolidFaces = solidFaces.getCount();
3611 
3612  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
3613  *muSolid =0;
3614  _macroFields.viscosity.addArray(solidFaces,muSolid);
3615 
3616  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
3617  *nueSolid =0;
3618  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
3619 
3620  const T rho_init=_options["rho_init"];
3621  const T T_init= _options["T_init"];
3622  const T mu_w= _options["mu_w"];
3623  const T Tmuref= _options["Tmuref"];
3624  const T muref= _options["muref"];
3625  const T R=8314.0/_options["molecularWeight"];
3626  const T nondim_length=_options["nonDimLt"];
3627 
3628  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
3629 
3630  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
3631  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
3632  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
3633  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
3634 
3635  for(int c=0; c<nSolidFaces;c++)
3636  {
3637  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
3638  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
3639  }
3640 
3641  if(_options.fgamma==2){
3642  for(int c=0; c<nSolidFaces;c++)
3643  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
3644  }
3645 
3646 
3647 
3648 
3649  for (int n=0; n<numMeshes; n++)
3650  {
3651  const Mesh& mesh = *_meshes[n];
3652  const Mesh& fMesh = *_finestMeshes[n];
3653  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3654 
3655  const StorageSite& cells = mesh.getCells();
3656  const StorageSite& fCells = fMesh.getCells();
3657  const StorageSite& ibFaces = mesh.getIBFaces();
3658  const StorageSite& fIbFaces = fMesh.getIBFaces();
3659 
3660  GeomFields::SSPair key1(&fIbFaces,&fCells);
3661  const IMatrix& mIC =
3662  dynamic_cast<const IMatrix&>
3663  (*_finestGeomFields._interpolationMatrices[key1]);
3664 
3665  IMatrix mICV(mIC);
3666  IMatrixV3 mICV3(mIC);
3667 
3668  GeomFields::SSPair key2(&fIbFaces,&solidFaces);
3669  const IMatrix& mIP =
3670  dynamic_cast<const IMatrix&>
3671  (*_finestGeomFields._interpolationMatrices[key2]);
3672 
3673  IMatrix mIPV(mIP);
3674  IMatrixV3 mIPV3(mIP);
3675 
3676  shared_ptr<TArray> ibVtemp(new TArray(ibFaces.getCount()));
3677  shared_ptr<TArray> ibVnue(new TArray(ibFaces.getCount()));
3678  shared_ptr<TArray> ibVdensity(new TArray(ibFaces.getCount()));
3679  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
3680 
3681  const TArray& cTemp =
3682  dynamic_cast<TArray&>(_macroFields.temperature[cells]);
3683  const VectorT3Array& cVel =
3684  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
3685  const TArray& cDensity =
3686  dynamic_cast<TArray&>(_macroFields.density[cells]);
3687  const TArray& sDensity =
3688  dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
3689  const TArray& cNue =
3690  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
3691  const TArray& sNue =
3692  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
3693  const TArray& sTemp =
3694  dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
3695  const VectorT3Array& sVel =
3696  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
3697 
3698  ibVnue->zero();
3699  ibVtemp->zero();
3700  ibVvel->zero();
3701  ibVdensity->zero();
3702 
3703 
3704  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
3705  const Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<const Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
3706 
3707  TArray cFTemp(fCells.getCount());
3708  VectorT3Array cFVel(fCells.getCount());
3709  TArray cFDensity(fCells.getCount());
3710  TArray cFNue(fCells.getCount());
3711 
3712  for(int c=0;c<fCells.getCount();c++)
3713  {
3714  cFTemp[c]=cTemp[FinestToCoarse[c][_level]];
3715  cFVel[c]=cVel[FinestToCoarse[c][_level]];
3716  cFDensity[c]=cDensity[FinestToCoarse[c][_level]];
3717  cFNue[c]=cNue[FinestToCoarse[c][_level]];
3718  }
3719 
3720  //nue interpolation (cells)
3721  mICV.multiplyAndAdd(*ibVnue,cFNue);
3722  mIPV.multiplyAndAdd(*ibVnue,sNue);
3723  _macroFields.collisionFrequency.addArray(ibFaces,ibVnue);
3724  //temperature interpolation (cells+solidfaces)
3725  mICV.multiplyAndAdd(*ibVtemp,cFTemp);
3726  mIPV.multiplyAndAdd(*ibVtemp,sTemp);
3727  _macroFields.temperature.addArray(ibFaces,ibVtemp);
3728  //density interpolation (cells+solidfaces)
3729  mICV.multiplyAndAdd(*ibVdensity,cFDensity);
3730  mIPV.multiplyAndAdd(*ibVdensity,sDensity);
3731  _macroFields.density.addArray(ibFaces,ibVdensity);
3732  //velocity interpolation (cells+solidfaces)
3733  mICV3.multiplyAndAdd(*ibVvel,cFVel);
3734  mIPV3.multiplyAndAdd(*ibVvel,sVel);
3735  _macroFields.velocity.addArray(ibFaces,ibVvel);
3736 
3737 
3738  }
3739  }
3740 
3741  const int f_out = 3;
3742  if (f_out ==1){
3743  //Step 2 Find fgamma using macroparameters
3744  const int numMeshes = _meshes.size();
3745  for (int n=0; n<numMeshes; n++)
3746  {
3747  const Mesh& mesh = *_meshes[n];
3748  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3749  const int numDirections = _quadrature.getDirCount();
3750  const StorageSite& ibFaces = mesh.getIBFaces();
3751  const int nibFaces=ibFaces.getCount();
3752  const double pi=_options.pi;
3753  const TArray& ibTemp =
3754  dynamic_cast<TArray&>(_macroFields.temperature[ibFaces]);
3755  const VectorT3Array& ibVel =
3756  dynamic_cast<VectorT3Array&>(_macroFields.velocity[ibFaces]);
3757  const TArray& ibDensity =
3758  dynamic_cast<TArray&>(_macroFields.density[ibFaces]);
3759 
3760  for (int j=0; j<numDirections; j++)
3761  {
3762  shared_ptr<TArray> ibFndPtrEqES(new TArray(nibFaces));
3763  TArray& ibFndEqES= *ibFndPtrEqES;
3764  ibFndPtrEqES->zero();
3765  Field& fndEqES = *_dsfEqPtrES.dsf[j];
3766 
3767  for (int i=0; i<nibFaces; i++)
3768  {
3769  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3770  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3771  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3772  const T ibu = ibVel[i][0];
3773  const T ibv = ibVel[i][1];
3774  const T ibw = ibVel[i][2];
3775  ibFndEqES[i]=ibDensity[i]/pow(pi*ibTemp[i],1.5)*exp(-(pow(cx[j]-ibu,2.0)+pow(cy[j]-ibv,2.0)+pow(cz[j]-ibw,2.0))/ibTemp[i]);
3776  }
3777  fndEqES.addArray(ibFaces,ibFndPtrEqES);
3778  }
3779  }
3780  }
3781  }
3782  else if(f_out==2)
3783  {
3784  //Step 2 Find fgamma using interpolation (only ESBGK for now)
3785  for (int n=0; n<numMeshes; n++)
3786  {
3787  const int numFields= _quadrature.getDirCount();
3788  for (int direction = 0; direction < numFields; direction++)
3789  {
3790  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
3791  const Mesh& mesh = *_meshes[n];
3792  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3793 
3794  const StorageSite& cells = mesh.getCells();
3795  const StorageSite& ibFaces = mesh.getIBFaces();
3796 
3797  GeomFields::SSPair key1(&ibFaces,&cells);
3798  const IMatrix& mIC =
3799  dynamic_cast<const IMatrix&>
3800  (*_geomFields._interpolationMatrices[key1]);
3801 
3802  IMatrix mICV(mIC);
3803 
3804  GeomFields::SSPair key2(&ibFaces,&solidFaces);
3805  const IMatrix& mIP =
3806  dynamic_cast<const IMatrix&>
3807  (*_geomFields._interpolationMatrices[key2]);
3808 
3809  IMatrix mIPV(mIP);
3810 
3811  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
3812 
3813  const TArray& cf =
3814  dynamic_cast<const TArray&>(fndEqES[cells]);
3815 
3816  ibVf->zero();
3817 
3818  //distribution function interpolation (cells)
3819  mICV.multiplyAndAdd(*ibVf,cf);
3820  fndEqES.addArray(ibFaces,ibVf);
3821 
3822  }
3823  }
3824  }
3825  }
3826  //Step3: Relax Distribution function from ibfaces to solid face
3827  for (int n=0; n<numMeshes; n++)
3828  {
3829  const int numDirections = _quadrature.getDirCount();
3830  for (int j=0; j<numDirections; j++)
3831  {
3832  Field& fnd = *_dsfPtr.dsf[j];
3833  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
3834 
3835 #ifdef FVM_PARALLEL
3836  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsf.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
3837 #endif
3838  const Mesh& mesh = *_meshes[n];
3839  const Mesh& fMesh = *_finestMeshes[n];
3840  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3841  const StorageSite& ibFaces = mesh.getIBFaces();
3842  const StorageSite& fIbFaces = fMesh.getIBFaces();
3843  TArray& dsfIB = dynamic_cast< TArray&>(fnd[ibFaces]);
3844  Field& fndEqES = *_dsfEqPtrES.dsf[j];
3845  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[ibFaces]);
3846  const StorageSite& faces = fMesh.getFaces();
3847  const StorageSite& cells = mesh.getCells();
3848  const CRConnectivity& faceCells = mesh.getAllFaceCells();
3849  const CRConnectivity& ibFacesTosolidFaces
3850  = fMesh.getConnectivity(fIbFaces,solidFaces);
3851  const IntArray& ibFaceIndices = fMesh.getIBFaceList();
3852  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
3853  const IntArray& sFCRow = ibFacesTosolidFaces.getRow();
3854  const IntArray& sFCCol = ibFacesTosolidFaces.getCol();
3855  const int nibFaces = ibFaces.getCount();
3856  const int nFaces = faces.getCount();
3857  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3858  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3859  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3860  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
3861 
3862  for(int f=0; f<nibFaces; f++)
3863  {
3864  dsfIB[f]=0.0;
3865  double distIBSolidInvSum(0.0);
3866  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
3867  {
3868  const int c = sFCCol[nc];
3869  const int faceIB= ibFaceIndices[f];
3870  const VectorT3Array& solidFaceCentroid =
3871  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[solidFaces]);
3872  const VectorT3Array& faceCentroid =
3873  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[faces]);
3874 
3875  double distIBSolid (0.0);
3876  // based on distance - will be thought
3877  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
3878  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
3879  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
3880  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
3881  }
3882  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
3883  {
3884  const int c = sFCCol[nc];
3885  const VectorT3Array& solidFaceCentroid =
3886  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[solidFaces]);
3887  const VectorT3Array& faceCentroid =
3888  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[faces]);
3889  const TArray& nue =
3890  dynamic_cast<TArray&>(_macroFields.collisionFrequency[ibFaces]);
3891  const int faceIB= ibFaceIndices[f];
3892  // const T coeff = iCoeffs[nc];
3893  double time_to_wall (0.0);
3894  double distIBSolid (0.0);
3895  const T uwall = v[c][0];
3896  const T vwall = v[c][1];
3897  const T wwall = v[c][2];
3898  // based on distance - will be thought
3899  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
3900  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
3901  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
3902  time_to_wall = -1*(pow(distIBSolid,2)/((cx[j]-uwall)*(faceCentroid[faceIB][0]-solidFaceCentroid[c][0])+(cy[j]-vwall)*(faceCentroid[faceIB][1]-solidFaceCentroid[c][1])+(cz[j]-wwall)*(faceCentroid[faceIB][2]-solidFaceCentroid[c][2])));
3903  if(time_to_wall<0)
3904  time_to_wall = 0;
3905 
3906  dsfIB[f] += (dsfEqES[f]-(dsfEqES[f]-dsf[c])*exp(-time_to_wall*nue[f]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
3907 
3908  }
3909 
3910  }
3911  }
3912  }
3913  }
3914  }
3915  if (method==3){
3916  const int numMeshes = _meshes.size();
3917  for (int n=0; n<numMeshes; n++)
3918  {
3919  const Mesh& mesh = *_meshes[n];
3920  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3921 
3922  const StorageSite& cells = mesh.getCells();
3923  const StorageSite& ibFaces = mesh.getIBFaces();
3924 
3925  GeomFields::SSPair key1(&ibFaces,&cells);
3926  const IMatrix& mIC =
3927  dynamic_cast<const IMatrix&>
3928  (*_geomFields._interpolationMatrices[key1]);
3929 
3930  IMatrixV3 mICV3(mIC);
3931 
3932  GeomFields::SSPair key2(&ibFaces,&solidFaces);
3933  const IMatrix& mIP =
3934  dynamic_cast<const IMatrix&>
3935  (*_geomFields._interpolationMatrices[key2]);
3936 
3937  IMatrixV3 mIPV3(mIP);
3938 
3939  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
3940 
3941 
3942  const VectorT3Array& cVel =
3943  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
3944  const VectorT3Array& sVel =
3945  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
3946 
3947  ibVvel->zero();
3948 
3949  //velocity interpolation (cells+solidfaces)
3950  mICV3.multiplyAndAdd(*ibVvel,cVel);
3951  mIPV3.multiplyAndAdd(*ibVvel,sVel);
3952  _macroFields.velocity.addArray(ibFaces,ibVvel);
3953 
3954 
3955  }
3956  }
3957  const int nSolidFaces = solidFaces.getCount();
3958 
3959  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
3960  *muSolid =0;
3961  _macroFields.viscosity.addArray(solidFaces,muSolid);
3962 
3963  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
3964  *nueSolid =0;
3965  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
3966 
3967  const T rho_init=_options["rho_init"];
3968  const T T_init= _options["T_init"];
3969  const T mu_w= _options["mu_w"];
3970  const T Tmuref= _options["Tmuref"];
3971  const T muref= _options["muref"];
3972  const T R=8314.0/_options["molecularWeight"];
3973  const T nondim_length=_options["nonDimLt"];
3974 
3975  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
3976 
3977  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
3978  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
3979  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
3980  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
3981 
3982  for(int c=0; c<nSolidFaces;c++)
3983  {
3984  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
3985  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
3986  }
3987 
3988  if(_options.fgamma==2){
3989  for(int c=0; c<nSolidFaces;c++)
3990  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
3991  }
3992 
3993  //Step 2 Find fgamma using interpolation (only ESBGK for now)
3994  const int numFields= _quadrature.getDirCount();
3995  for (int direction = 0; direction < numFields; direction++)
3996  {
3997  shared_ptr<TArray> ibVf(new TArray(solidFaces.getCount()));
3998  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
3999  ibVf->zero();
4000  for (int n=0; n<numMeshes; n++)
4001  {
4002  const Mesh& mesh = *_meshes[n];
4003  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4004 
4005  const StorageSite& cells = mesh.getCells();
4006  const StorageSite& ibFaces = mesh.getIBFaces();
4007 
4008  GeomFields::SSPair key1(&solidFaces,&cells);
4009  const IMatrix& mIC =
4010  dynamic_cast<const IMatrix&>
4011  (*_geomFields._interpolationMatrices[key1]);
4012 
4013  IMatrix mICV(mIC);
4014 
4015  const TArray& cf =
4016  dynamic_cast<const TArray&>(fndEqES[cells]);
4017 
4018  ibVf->zero();
4019 
4020  //distribution function interpolation (cells)
4021  mICV.multiplyAndAdd(*ibVf,cf);
4022  }
4023  }
4024  fndEqES.addArray(solidFaces,ibVf);
4025  }
4026  for (int n=0; n<numMeshes; n++)
4027  {
4028  const int numDirections = _quadrature.getDirCount();
4029  for (int j=0; j<numDirections; j++)
4030  {
4031  Field& fnd = *_dsfPtr.dsf[j];
4032  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
4033  Field& fndEqES = *_dsfEqPtrES.dsf[j];
4034  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[solidFaces]);
4035 #ifdef FVM_PARALLEL
4036  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsf.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
4037  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsfEqES.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
4038 #endif
4039 
4040  const Mesh& mesh = *_meshes[n];
4041  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4042  const StorageSite& ibFaces = mesh.getIBFaces();
4043  const StorageSite& faces = mesh.getFaces();
4044  const StorageSite& cells = mesh.getCells();
4045  const CRConnectivity& faceCells = mesh.getAllFaceCells();
4046  const CRConnectivity& ibFacesTosolidFaces
4047  = mesh.getConnectivity(ibFaces,solidFaces);
4048  const IntArray& ibFaceIndices = mesh.getIBFaceList();
4049  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4050  const IntArray& sFCRow = ibFacesTosolidFaces.getRow();
4051  const IntArray& sFCCol = ibFacesTosolidFaces.getCol();
4052  const int nibFaces = ibFaces.getCount();
4053  const int nFaces = faces.getCount();
4054  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4055  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4056  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4057  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
4058  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
4059  ibVf->zero();
4060  TArray& ibVfA= *ibVf;
4061  for(int f=0; f<nibFaces; f++)
4062  {
4063  double distIBSolidInvSum(0.0);
4064  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4065  {
4066  const int c = sFCCol[nc];
4067  const int faceIB= ibFaceIndices[f];
4068  const VectorT3Array& solidFaceCentroid =
4069  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
4070  const VectorT3Array& faceCentroid =
4071  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
4072 
4073  double distIBSolid (0.0);
4074  // based on distance - will be thought
4075  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
4076  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
4077  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
4078  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
4079  }
4080  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4081  {
4082  const int c = sFCCol[nc];
4083  const VectorT3Array& solidFaceCentroid =
4084  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
4085  const VectorT3Array& faceCentroid =
4086  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
4087  const TArray& nue =
4088  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
4089  const TArray& nueC =
4090  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
4091  const int faceIB= ibFaceIndices[f];
4092  const T uwall = v[c][0];
4093  const T vwall = v[c][1];
4094  const T wwall = v[c][2];
4095  // const T coeff = iCoeffs[nc];
4096  double time_to_wall (0.0);
4097  double distIBSolid (0.0);
4098  // based on distance - will be thought
4099  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
4100  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
4101  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
4102  time_to_wall = -1*(pow(distIBSolid,2)/((cx[j]-uwall)*(faceCentroid[faceIB][0]-solidFaceCentroid[c][0])+(cy[j]-vwall)*(faceCentroid[faceIB][1]-solidFaceCentroid[c][1])+(cz[j]-wwall)*(faceCentroid[faceIB][2]-solidFaceCentroid[c][2])));
4103  if(time_to_wall<0)
4104  time_to_wall = 0;
4105  ibVfA[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
4106  }
4107 
4108  }
4109 
4110  fnd.addArray(ibFaces,ibVf);
4111  }
4112  }
4113  }
4114  }
4115  }
4116 
4117  void computeSolidFacePressure(const StorageSite& solidFaces)
4118  {
4119  typedef CRMatrixTranspose<T,T,T> IMatrix;
4120  const int numMeshes = _meshes.size();
4121  for (int n=0; n<numMeshes; n++)
4122  {
4123  shared_ptr<TArray> ibP(new TArray(solidFaces.getCount()));
4124  ibP->zero();
4125  const Mesh& mesh = *_meshes[n];
4126  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4127 
4128  const StorageSite& cells = mesh.getCells();
4129 
4130  GeomFields::SSPair key1(&solidFaces,&cells);
4131  const IMatrix& mIC =
4132  dynamic_cast<const IMatrix&>
4133  (*_geomFields._interpolationMatrices[key1]);
4134 
4135  IMatrix mICV(mIC);
4136 
4137  const TArray& cP =
4138  dynamic_cast<TArray&>(_macroFields.pressure[cells]);
4139 
4140  ibP->zero();
4141 
4142 
4143  //nue interpolation (cells)
4144  mICV.multiplyAndAdd(*ibP,cP);
4145  }
4146  _macroFields.pressure.addArray(solidFaces,ibP);
4147  }
4148 #ifdef FVM_PARALLEL
4149  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[solidFaces]);
4150  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,pressure.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
4151 #endif
4152  }
4153 
4154  void computeSolidFaceDsf(const StorageSite& solidFaces,const int method,const int RelaxDistribution=0)
4155  {
4156  typedef CRMatrixTranspose<T,T,T> IMatrix;
4157  typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
4158  const int numFields= _quadrature.getDirCount();
4159  if (method==1){
4160  const int numMeshes = _meshes.size();
4161  for (int direction = 0; direction < numFields; direction++) {
4162  Field& fnd = *_dsfPtr.dsf[direction];
4163  shared_ptr<TArray> ibV(new TArray(solidFaces.getCount()));
4164  ibV->zero();
4165  for (int n=0; n<numMeshes; n++)
4166  {
4167  const Mesh& mesh = *_meshes[n];
4168  const Mesh& fMesh = *_finestMeshes[n];
4169  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4170 
4171  const StorageSite& cells = mesh.getCells();
4172  const StorageSite& fCells = fMesh.getCells();
4173 
4174  GeomFields::SSPair key1(&solidFaces,&fCells);
4175  const IMatrix& mIC =
4176  dynamic_cast<const IMatrix&>
4177  (*_finestGeomFields._interpolationMatrices[key1]);
4178 
4179  IMatrix mICV(mIC);
4180 
4181  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
4182  const Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<const Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
4183 
4184  const TArray& cV =
4185  dynamic_cast<const TArray&>(fnd[cells]);
4186  TArray cFV(fCells.getCountLevel1());
4187  for(int c=0;c<fCells.getCountLevel1();c++)
4188  cFV[c]=cV[FinestToCoarse[c][_level]];
4189 
4190  ibV->zero();
4191 
4192  mICV.multiplyAndAdd(*ibV,cFV);
4193 #if 0
4194  ofstream debugFile;
4195  stringstream ss(stringstream::in | stringstream::out);
4196  ss << MPI::COMM_WORLD.Get_rank();
4197  string fname1 = "IBVelocity_proc" + ss.str() + ".dat";
4198  debugFile.open(fname1.c_str());
4199 
4200  //debug use
4201  const Array<int>& ibFaceList = mesh.getIBFaceList();
4202  const StorageSite& faces = mesh.getFaces();
4203  const VectorT3Array& faceCentroid =
4204  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
4205  const double angV = 1.0;
4206  VectorT3 center;
4207  center[0]=0.;
4208  center[1]=0.;
4209  center[2]=0.;
4210 
4211  for(int f=0; f<ibFaces.getCount();f++){
4212  int fID = ibFaceList[f];
4213  debugFile << "f=" << f << setw(10) << " fID = " << fID << " faceCentroid = " << faceCentroid[fID] << " ibV = " << (*ibV)[f] << endl;
4214  }
4215 
4216  debugFile.close();
4217 #endif
4218 
4219  }
4220 
4221  }
4222  fnd.addArray(solidFaces,ibV);
4223  }
4224  }
4225  if (method==2){
4226 
4227  // Step0: Compute Interpolation Matrices from (only) Cells to IBFaces
4228  const int numMeshes = _meshes.size();
4229  for (int n=0; n<numMeshes; n++)
4230  {
4231  const Mesh& mesh = *_meshes[n];
4232  const Mesh& fMesh = *_finestMeshes[n];
4233 
4234  const int numFields= _quadrature.getDirCount();
4235  for (int direction = 0; direction < numFields; direction++)
4236  {
4237  Field& fnd = *_dsfPtr.dsf[direction];
4238  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
4239  const Mesh& mesh = *_meshes[n];
4240  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4241 
4242  const StorageSite& cells = mesh.getCells();
4243  const StorageSite& fCells = fMesh.getCells();
4244  const StorageSite& faces = mesh.getFaces();
4245  const StorageSite& fFaces = fMesh.getFaces();
4246  const StorageSite& ibFaces = mesh.getIBFaces();
4247  const StorageSite& fIbFaces = fMesh.getIBFaces();
4248 
4249  //GeomFields::SSPair key1(&faces,&cells);
4250  GeomFields::SSPair key1(&fFaces,&fCells);
4251  //GeomFields::SSPair key1(&fIbFaces,&fCells);
4252  const IMatrix& mIC =
4253  dynamic_cast<const IMatrix&>
4254  (*_finestGeomFields._interpolationMatrices[key1]);
4255 
4256  IMatrix mICV(mIC);
4257 
4258  const TArray& cf =
4259  dynamic_cast<const TArray&>(fnd[cells]);
4260  const TArray& cfEq =
4261  dynamic_cast<const TArray&>(fndEqES[cells]);
4262 
4263  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
4264  const Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<const Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
4265 
4266  TArray cFf(fCells.getCount());
4267  TArray cFfEq(fCells.getCount());
4268 
4269  for(int c=0;c<fCells.getCount();c++)
4270  {
4271  cFf[c]=cf[FinestToCoarse[c][_level]];
4272  cFfEq[c]=cfEq[FinestToCoarse[c][_level]];
4273  }
4274 
4275  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
4276  ibVf->zero();
4277  if (_options.fgamma==2){
4278  shared_ptr<TArray> ibVfEq(new TArray(ibFaces.getCount()));
4279  ibVfEq->zero();
4280  mICV.multiplyAndAdd(*ibVfEq,cFfEq);
4281  fndEqES.addArray(ibFaces,ibVfEq);
4282  }
4283  mICV.multiplyAndAdd(*ibVf,cFf);
4284  fnd.addArray(ibFaces,ibVf);
4285 
4286  }
4287  }
4288  }
4289  const int nSolidFaces = solidFaces.getCount();
4290 
4291  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
4292  *muSolid =0;
4293  _macroFields.viscosity.addArray(solidFaces,muSolid);
4294 
4295  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
4296  *nueSolid =0;
4297  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
4298 
4299  const T rho_init=_options["rho_init"];
4300  const T T_init= _options["T_init"];
4301  const T mu_w= _options["mu_w"];
4302  const T Tmuref= _options["Tmuref"];
4303  const T muref= _options["muref"];
4304  const T R=8314.0/_options["molecularWeight"];
4305  const T nondim_length=_options["nonDimLt"];
4306 
4307  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
4308 
4309  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
4310  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
4311  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
4312  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
4313 
4314  for(int c=0; c<nSolidFaces;c++)
4315  {
4316  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
4317  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
4318  }
4319 
4320  if(_options.fgamma==2){
4321  for(int c=0; c<nSolidFaces;c++)
4322  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
4323  }
4324 
4325  //Step 1 Interpolate Macroparameters and f to IBface
4326  for (int n=0; n<numMeshes; n++)
4327  {
4328  const Mesh& mesh = *_meshes[n];
4329  const Mesh& fMesh = *_finestMeshes[n];
4330  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4331 
4332  const StorageSite& cells = mesh.getCells();
4333  const StorageSite& fCells = fMesh.getCells();
4334  const StorageSite& ibFaces = mesh.getIBFaces();
4335  const StorageSite& fIbFaces = fMesh.getIBFaces();
4336 
4337  GeomFields::SSPair key1(&fIbFaces,&fCells);
4338  const IMatrix& mIC =
4339  dynamic_cast<const IMatrix&>
4340  (*_finestGeomFields._interpolationMatrices[key1]);
4341 
4342  IMatrix mICV(mIC);
4343  IMatrixV3 mICV3(mIC);
4344 
4345  GeomFields::SSPair key2(&fIbFaces,&solidFaces);
4346  const IMatrix& mIP =
4347  dynamic_cast<const IMatrix&>
4348  (*_finestGeomFields._interpolationMatrices[key2]);
4349 
4350  IMatrix mIPV(mIP);
4351  IMatrixV3 mIPV3(mIP);
4352 
4353  shared_ptr<TArray> ibVtemp(new TArray(ibFaces.getCount()));
4354  shared_ptr<TArray> ibVnue(new TArray(ibFaces.getCount()));
4355  shared_ptr<TArray> ibVdensity(new TArray(ibFaces.getCount()));
4356  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
4357 
4358  const TArray& cTemp =
4359  dynamic_cast<TArray&>(_macroFields.temperature[cells]);
4360  const VectorT3Array& cVel =
4361  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
4362  const TArray& cDensity =
4363  dynamic_cast<TArray&>(_macroFields.density[cells]);
4364  const TArray& sDensity =
4365  dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
4366  const TArray& cNue =
4367  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
4368  const TArray& sNue =
4369  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
4370  const TArray& sTemp =
4371  dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
4372  const VectorT3Array& sVel =
4373  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
4374 
4375  ibVnue->zero();
4376  ibVtemp->zero();
4377  ibVvel->zero();
4378  ibVdensity->zero();
4379 
4380 
4381  Field& FinestToCoarseField=_finestGeomFields.finestToCoarse;
4382  const Array<Vector<int,25> >& FinestToCoarse=dynamic_cast<const Array<Vector<int,25> >&>(FinestToCoarseField[fCells]);
4383 
4384  TArray cFTemp(fCells.getCount());
4385  VectorT3Array cFVel(fCells.getCount());
4386  TArray cFDensity(fCells.getCount());
4387  TArray cFNue(fCells.getCount());
4388 
4389  for(int c=0;c<fCells.getCount();c++)
4390  {
4391  cFTemp[c]=cTemp[FinestToCoarse[c][_level]];
4392  cFVel[c]=cVel[FinestToCoarse[c][_level]];
4393  cFDensity[c]=cDensity[FinestToCoarse[c][_level]];
4394  cFNue[c]=cNue[FinestToCoarse[c][_level]];
4395  }
4396 
4397  //nue interpolation (cells)
4398  mICV.multiplyAndAdd(*ibVnue,cFNue);
4399  mIPV.multiplyAndAdd(*ibVnue,sNue);
4400  _macroFields.collisionFrequency.addArray(ibFaces,ibVnue);
4401  //temperature interpolation (cells+solidfaces)
4402  mICV.multiplyAndAdd(*ibVtemp,cFTemp);
4403  mIPV.multiplyAndAdd(*ibVtemp,sTemp);
4404  _macroFields.temperature.addArray(ibFaces,ibVtemp);
4405  //density interpolation (cells+solidfaces)
4406  mICV.multiplyAndAdd(*ibVdensity,cFDensity);
4407  mIPV.multiplyAndAdd(*ibVdensity,sDensity);
4408  _macroFields.density.addArray(ibFaces,ibVdensity);
4409  //velocity interpolation (cells+solidfaces)
4410  mICV3.multiplyAndAdd(*ibVvel,cFVel);
4411  mIPV3.multiplyAndAdd(*ibVvel,sVel);
4412  _macroFields.velocity.addArray(ibFaces,ibVvel);
4413 
4414 
4415  }
4416  }
4417  if (_options.fgamma==1){
4418  //Step 2 Find fgamma using macroparameters
4419  const int numMeshes = _meshes.size();
4420  for (int n=0; n<numMeshes; n++)
4421  {
4422  const Mesh& mesh = *_meshes[n];
4423  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4424  const int numDirections = _quadrature.getDirCount();
4425  const StorageSite& ibFaces = mesh.getIBFaces();
4426  const int nibFaces=ibFaces.getCount();
4427  const double pi=_options.pi;
4428  const TArray& ibTemp =
4429  dynamic_cast<TArray&>(_macroFields.temperature[ibFaces]);
4430  const VectorT3Array& ibVel =
4431  dynamic_cast<VectorT3Array&>(_macroFields.velocity[ibFaces]);
4432  const TArray& ibDensity =
4433  dynamic_cast<TArray&>(_macroFields.density[ibFaces]);
4434 
4435  for (int j=0; j<numDirections; j++)
4436  {
4437  shared_ptr<TArray> ibFndPtrEqES(new TArray(nibFaces));
4438  TArray& ibFndEqES= *ibFndPtrEqES;
4439 
4440  ibFndPtrEqES->zero();
4441 
4442  Field& fndEqES = *_dsfEqPtrES.dsf[j];
4443 
4444  for (int i=0; i<nibFaces; i++)
4445  {
4446  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4447  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4448  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4449  const T ibu = ibVel[i][0];
4450  const T ibv = ibVel[i][1];
4451  const T ibw = ibVel[i][2];
4452  ibFndEqES[i]=ibDensity[i]/pow(pi*ibTemp[i],1.5)*exp(-(pow(cx[j]-ibu,2.0)+pow(cy[j]-ibv,2.0)+pow(cz[j]-ibw,2.0))/ibTemp[i]);
4453  }
4454  fndEqES.addArray(ibFaces,ibFndPtrEqES);
4455  }
4456  }
4457  }
4458  }
4459 
4460  //Step3: Relax Distribution function from ibfaces to solid face
4461  const int numDirections = _quadrature.getDirCount();
4462  for (int j=0; j<numDirections; j++)
4463  {
4464  const int nSolidFaces = solidFaces.getCount();
4465  shared_ptr<TArray> solidFndPtr(new TArray(nSolidFaces));
4466  solidFndPtr->zero();
4467  TArray& solidFnd= *solidFndPtr;
4468  Field& fnd = *_dsfPtr.dsf[j];
4469  Field& fndEqES = *_dsfEqPtrES.dsf[j];
4470  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4471  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4472  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4473  for (int n=0; n<numMeshes; n++)
4474  {
4475  const Mesh& mesh = *_meshes[n];
4476  const Mesh& fMesh = *_finestMeshes[n];
4477  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4478  const StorageSite& ibFaces = mesh.getIBFaces();
4479  const StorageSite& fIbFaces = fMesh.getIBFaces();
4480  const CRConnectivity& solidFacesToibFaces
4481  = fMesh.getConnectivity(solidFaces,fIbFaces);
4482  const IntArray& ibFaceIndices = fMesh.getIBFaceList();
4483  const IntArray& sFCRow = solidFacesToibFaces.getRow();
4484  const IntArray& sFCCol = solidFacesToibFaces.getCol();
4485  TArray& dsf = dynamic_cast< TArray&>(fnd[ibFaces]);
4486  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[ibFaces]);
4487  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
4488 
4489  for(int f=0; f<nSolidFaces; f++)
4490  {
4491  double distIBSolidInvSum(0.0);
4492  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4493  {
4494  const StorageSite& faces = fMesh.getFaces();
4495  const int c = sFCCol[nc];
4496  const int faceIB= ibFaceIndices[c];
4497  const VectorT3Array& solidFaceCentroid =
4498  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[solidFaces]);
4499  const VectorT3Array& faceCentroid =
4500  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[faces]);
4501 
4502  double distIBSolid (0.0);
4503  // based on distance - will be thought
4504  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[f][0]),2)+
4505  pow((faceCentroid[faceIB][1]-solidFaceCentroid[f][1]),2)+
4506  pow((faceCentroid[faceIB][2]-solidFaceCentroid[f][2]),2));
4507  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
4508  }
4509  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4510  {
4511  const int c = sFCCol[nc];
4512  const StorageSite& faces = fMesh.getFaces();
4513  const VectorT3Array& solidFaceCentroid =
4514  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[solidFaces]);
4515  const VectorT3Array& faceCentroid =
4516  dynamic_cast<const VectorT3Array&> (_finestGeomFields.coordinate[faces]);
4517  const int faceIB= ibFaceIndices[c];
4518  T time_to_wall (0.0);
4519  T distIBSolid (0.0);
4520  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[f][0]),2)+
4521  pow((faceCentroid[faceIB][1]-solidFaceCentroid[f][1]),2)+
4522  pow((faceCentroid[faceIB][2]-solidFaceCentroid[f][2]),2));
4523  // based on distance - will be thought
4524  const T uwall = v[f][0];
4525  const T vwall = v[f][1];
4526  const T wwall = v[f][2];
4527  const TArray& nue =
4528  dynamic_cast<TArray&>(_macroFields.collisionFrequency[ibFaces]);
4529  time_to_wall = (pow(distIBSolid,2)/((cx[j]-uwall)*(faceCentroid[faceIB][0]-solidFaceCentroid[f][0])+(cy[j]-vwall)*(faceCentroid[faceIB][1]-solidFaceCentroid[f][1])+(cz[j]-wwall)*(faceCentroid[faceIB][2]-solidFaceCentroid[f][2])));
4530  if(time_to_wall<0)
4531  time_to_wall = 0;
4532 
4533  solidFnd[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
4534  }
4535 
4536  }
4537  }
4538  }
4539  fnd.addArray(solidFaces,solidFndPtr);
4540  }
4541  }
4542  if (method==3){
4543  const int numDirections = _quadrature.getDirCount();
4544  for (int j=0; j<numDirections; j++)
4545  {
4546  const int nSolidFaces = solidFaces.getCount();
4547  shared_ptr<TArray> solidFndPtr(new TArray(nSolidFaces));
4548  solidFndPtr->zero();
4549  TArray& solidFnd= *solidFndPtr;
4550  Field& fnd = *_dsfPtr.dsf[j];
4551  Field& fndEqES = *_dsfEqPtrES.dsf[j];
4552  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4553  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4554  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4555  const int numMeshes = _meshes.size();
4556  for (int n=0; n<numMeshes; n++)
4557  {
4558  const Mesh& mesh = *_meshes[n];
4559  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
4560  const StorageSite& cells = mesh.getCells();
4561  const StorageSite& ibFaces = mesh.getIBFaces();
4562  const CRConnectivity& solidFacesToCells
4563  = mesh.getConnectivity(solidFaces,cells);
4564  const IntArray& sFCRow = solidFacesToCells.getRow();
4565  const IntArray& sFCCol = solidFacesToCells.getCol();
4566  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4567  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[cells]);
4568  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
4569 
4570  for(int f=0; f<nSolidFaces; f++)
4571  {
4572  double distIBSolidInvSum(0.0);
4573  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4574  {
4575  const StorageSite& faces = mesh.getFaces();
4576  const int c = sFCCol[nc];
4577  const VectorT3Array& solidFaceCentroid =
4578  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
4579  const VectorT3Array& faceCentroid =
4580  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[cells]);
4581 
4582  double distIBSolid (0.0);
4583  // based on distance - will be thought
4584  distIBSolid = sqrt(pow((faceCentroid[c][0]-solidFaceCentroid[f][0]),2)+
4585  pow((faceCentroid[c][1]-solidFaceCentroid[f][1]),2)+
4586  pow((faceCentroid[c][2]-solidFaceCentroid[f][2]),2));
4587  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
4588  }
4589  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
4590  {
4591  const int c = sFCCol[nc];
4592  const StorageSite& faces = mesh.getFaces();
4593  const VectorT3Array& solidFaceCentroid =
4594  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
4595  const VectorT3Array& faceCentroid =
4596  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[cells]);
4597  T time_to_wall (0.0);
4598  T distIBSolid (0.0);
4599  const T uwall = v[f][0];
4600  const T vwall = v[f][1];
4601  const T wwall = v[f][2];
4602  distIBSolid = sqrt(pow((faceCentroid[c][0]-solidFaceCentroid[f][0]),2)+
4603  pow((faceCentroid[c][1]-solidFaceCentroid[f][1]),2)+
4604  pow((faceCentroid[c][2]-solidFaceCentroid[f][2]),2));
4605  // based on distance - will be thought
4606 
4607  const TArray& nue =
4608  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
4609  time_to_wall = (pow(distIBSolid,2)/((cx[j]-uwall)*(faceCentroid[c][0]-solidFaceCentroid[f][0])+(cy[j]-vwall)*(faceCentroid[c][1]-solidFaceCentroid[f][1])+(cz[j]-wwall)*(faceCentroid[c][2]-solidFaceCentroid[f][2])));
4610  if(time_to_wall<0)
4611  time_to_wall = 0;
4612 
4613  solidFnd[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
4614  }
4615 
4616  }
4617  }
4618  }
4619  fnd.addArray(solidFaces,solidFndPtr);
4620 
4621  }
4622  }
4623  }
4624 
4625 
4626 
4627 
4628  void correctMassDeficit()
4629  {
4630 
4631  const int numMeshes = _meshes.size();
4632  T netFlux(0.0);
4633 
4634  for (int n=0; n<numMeshes; n++)
4635  {
4636  const Mesh& mesh = *_meshes[n];
4637  const StorageSite& ibFaces = mesh.getIBFaces();
4638  const StorageSite& cells = mesh.getCells();
4639  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4640  const StorageSite& faces = mesh.getFaces();
4641  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4642  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4643  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4644  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
4645  const int numDirections = _quadrature.getDirCount();
4646  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[faces]);
4647  const VectorT3Array& faceArea =
4648  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
4649  const TArray& faceAreaMag =
4650  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
4651  const CRConnectivity& faceCells = mesh.getAllFaceCells();
4652  const int nibFaces = ibFaces.getCount();
4653 
4654  for(int f=0; f<nibFaces; f++)
4655  {
4656  const int c0 = faceCells(f,0);
4657  const int c1 = faceCells(f,1);
4658  if (((ibType[c0] == Mesh::IBTYPE_FLUID) && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
4659  ((ibType[c1] == Mesh::IBTYPE_FLUID) && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
4660  {
4661  const int ibFace = ibFaceIndex[f];
4662  if (ibFace < 0)
4663  throw CException("invalid ib face index");
4664  if (ibType[c0] == Mesh::IBTYPE_FLUID)
4665  {
4666  const VectorT3 en = faceArea[f]/faceAreaMag[f];
4667  for (int j=0; j<numDirections; j++)
4668  {
4669  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
4670  Field& fnd = *_dsfPtr.dsf[j];
4671  TArray& dsf = dynamic_cast<TArray&>(fnd[ibFaces]);
4672 
4673  netFlux -= dsf[f]*c_dot_en*wts[j]/abs(c_dot_en);
4674  }
4675 
4676  }
4677  else
4678  {
4679  const VectorT3 en = faceArea[f]/faceAreaMag[f];
4680  for (int j=0; j<numDirections; j++)
4681  {
4682  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
4683  Field& fnd = *_dsfPtr.dsf[j];
4684  TArray& dsf = dynamic_cast<TArray&>(fnd[ibFaces]);
4685 
4686  netFlux += dsf[f]*c_dot_en*wts[j]/abs(c_dot_en);
4687  }
4688  }
4689  }
4690  }
4691  }
4692 
4693 #ifdef FVM_PARALLEL
4694  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &netFlux, 1, MPI::DOUBLE, MPI::SUM);
4695 #endif
4696 
4697  T volumeSum(0.);
4698 
4699  for (int n=0; n<numMeshes; n++)
4700  {
4701  const Mesh& mesh = *_meshes[n];
4702  const StorageSite& cells = mesh.getCells();
4703  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4704  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
4705  for(int c=0; c<cells.getSelfCount(); c++)
4706  if (ibType[c] == Mesh::IBTYPE_FLUID)
4707  volumeSum += cellVolume[c];
4708  }
4709 #ifdef FVM_PARALLEL
4710  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &volumeSum, 1, MPI::DOUBLE, MPI::SUM);
4711 #endif
4712 
4713  netFlux /= volumeSum;
4714  for (int n=0; n<numMeshes; n++)
4715  {
4716  const Mesh& mesh = *_meshes[n];
4717  const StorageSite& cells = mesh.getCells();
4718  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4719  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
4720  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
4721 
4722  for(int c=0; c<cells.getSelfCount(); c++)
4723  {
4724  if (ibType[c] == Mesh::IBTYPE_FLUID){
4725  const int numDirections = _quadrature.getDirCount();
4726  T cellMass(0.0);
4727  for (int j=0; j<numDirections; j++)
4728  {
4729  Field& fnd = *_dsfPtr.dsf[j];
4730  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4731  cellMass += wts[j]*dsf[c];
4732  }
4733 
4734  for (int j=0; j<numDirections; j++)
4735  {
4736  Field& fnd = *_dsfPtr.dsf[j];
4737  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4738  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
4739  TArray& fgam = dynamic_cast< TArray&>(feqES[cells]);
4740  fgam[c] = fgam[c]*(1+netFlux*cellVolume[c]/cellMass);
4741  dsf[c] = dsf[c]*(1+netFlux*cellVolume[c]/cellMass);
4742  }
4743  }
4744 
4745  }
4746  }
4747  }
4748 
4749  void correctMassDeficit2(double n1,double n2)
4750  {
4751 
4752  const int numMeshes = _meshes.size();
4753  T netFlux(0.0);
4754 
4755  netFlux=n2-n1;
4756 
4757  T volumeSum(0.);
4758 
4759  for (int n=0; n<numMeshes; n++)
4760  {
4761  const Mesh& mesh = *_meshes[n];
4762  const StorageSite& cells = mesh.getCells();
4763  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4764  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
4765  for(int c=0; c<cells.getSelfCount(); c++)
4766  if (ibType[c] == Mesh::IBTYPE_FLUID)
4767  volumeSum += cellVolume[c];
4768  }
4769 #ifdef FVM_PARALLEL
4770  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &volumeSum, 1, MPI::DOUBLE, MPI::SUM);
4771 #endif
4772 
4773  netFlux /= volumeSum;
4774  for (int n=0; n<numMeshes; n++)
4775  {
4776  const Mesh& mesh = *_meshes[n];
4777  const StorageSite& cells = mesh.getCells();
4778  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4779  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
4780  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
4781 
4782  for(int c=0; c<cells.getSelfCount(); c++)
4783  {
4784  if (ibType[c] == Mesh::IBTYPE_FLUID){
4785  const int numDirections = _quadrature.getDirCount();
4786  T cellMass(0.0);
4787  for (int j=0; j<numDirections; j++)
4788  {
4789  Field& fnd = *_dsfPtr.dsf[j];
4790  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4791  cellMass += wts[j]*dsf[c];
4792  }
4793 
4794  for (int j=0; j<numDirections; j++)
4795  {
4796  Field& fnd = *_dsfPtr.dsf[j];
4797  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4798  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
4799  TArray& fgam = dynamic_cast< TArray&>(feqES[cells]);
4800  fgam[c] = fgam[c]*(1+netFlux*cellVolume[c]/cellMass);
4801  dsf[c] = dsf[c]*(1+netFlux*cellVolume[c]/cellMass);
4802  }
4803  }
4804 
4805  }
4806  }
4807  }
4808  const double ConservationofMassCheck()
4809  {
4810  const int numMeshes = _meshes.size();
4811  T ndens_tot(0.0) ;
4812 
4813  for (int n=0; n<numMeshes; n++)
4814  {
4815  const Mesh& mesh = *_meshes[n];
4816  const StorageSite& cells = mesh.getCells();
4817  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
4818  const StorageSite& faces = mesh.getFaces();
4819  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
4820  const int numDirections = _quadrature.getDirCount();
4821  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[faces]);
4822  const VectorT3Array& faceArea =
4823  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
4824  const TArray& faceAreaMag =
4825  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
4826  const CRConnectivity& faceCells = mesh.getAllFaceCells();
4827  const int nFaces = faces.getCount();
4828 
4829  for(int c=0; c<cells.getCountLevel1(); c++)
4830  {
4831  if (ibType[c] == Mesh::IBTYPE_FLUID)
4832  {
4833  for (int j=0; j<numDirections; j++)
4834  {
4835  Field& fnd = *_dsfPtr.dsf[j];
4836  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
4837  ndens_tot += wts[j]*dsf[c];
4838  }
4839  }
4840  }
4841  }
4842  cout << "Hello, I have" << ndens_tot << "number density";
4843  return ndens_tot;
4844  }
4845 
4846  void ConservationofMFSolid(const StorageSite& solidFaces) const
4847  {
4848  const double pi=_options.pi;
4849  const double epsilon=_options.epsilon_ES;
4850  const int nSolidFaces = solidFaces.getCount();
4851  for (int i=0; i<nSolidFaces; i++)
4852  {
4853  const int numDirections = _quadrature.getDirCount();
4854  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
4855  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
4856  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
4857  const VectorT3Array& solidFaceCentroid =
4858  dynamic_cast<const VectorT3Array&>(_finestGeomFields.coordinate[solidFaces]);
4859  const VectorT3Array& solidFaceArea =
4860  dynamic_cast<const VectorT3Array&>(_finestGeomFields.area[solidFaces]);
4861  const TArray& solidFaceAreaMag =
4862  dynamic_cast<const TArray&>(_finestGeomFields.areaMag[solidFaces]);
4863  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
4864  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_finestMacroFields.velocity[solidFaces]);
4865  TArray& density = dynamic_cast<TArray&>(_finestMacroFields.density[solidFaces]);
4866  TArray& temperature = dynamic_cast<TArray&>(_finestMacroFields.temperature[solidFaces]);
4867  const T uwall = v[i][0];
4868  const T vwall = v[i][1];
4869  const T wwall = v[i][2];
4870  const T Twall = temperature[i];
4871 
4872  T Nmr(0.0) ;
4873  T Dmr(0.0) ;
4874  T incomFlux(0.0);
4875  for (int j=0; j<numDirections; j++)
4876  {
4877  Field& fnd = *_dsfPtr.dsf[j];
4878  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
4879  const VectorT3 en = solidFaceArea[i]/solidFaceAreaMag[i];
4880  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
4881  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
4882  const T 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);
4883 
4884  if (c_dot_en-wallV_dot_en > 0) //incoming
4885  {
4886  Dmr = Dmr - fwall*wts[j]*(c_dot_en-wallV_dot_en);
4887  incomFlux=incomFlux-dsf[i]*wts[j]*(c_dot_en-wallV_dot_en);
4888  }
4889  else
4890  {
4891  Nmr = Nmr + dsf[i]*wts[j]*(c_dot_en-wallV_dot_en);
4892  }
4893  }
4894  const T nwall = Nmr/Dmr; // wall number density for initializing Maxwellian
4895 
4896  density[i]=nwall;
4897 
4898  for (int j=0; j<numDirections; j++)
4899  {
4900  Field& fnd = *_dsfPtr.dsf[j];
4901  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
4902  const VectorT3 en = solidFaceArea[i]/solidFaceAreaMag[i];
4903  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
4904  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
4905  if (c_dot_en-wallV_dot_en > 0)
4906  {
4907  dsf[i] = 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);
4908  }
4909  else
4910  dsf[i]=dsf[i];
4911  }
4912  }
4913  }
4914 
4915  void doSweeps(const int sweeps, const int num)
4916  {
4917  for(int sweepNo=0;sweepNo<sweeps;sweepNo++)
4918  smooth(num);
4919  }
4920 
4921  void doSweeps(const int sweeps, const int num, const StorageSite& solidFaces)
4922  {
4923  for(int sweepNo=0;sweepNo<sweeps;sweepNo++)
4924  smooth(num,solidFaces);
4925  }
4926 
4927  void smooth(const int num,const StorageSite& solidFaces)
4928  {
4929  const int numDir=_quadrature.getDirCount();
4930  const int numMeshes=_meshes.size();
4931  for(int msh=0;msh<numMeshes;msh++)
4932  {
4933  const Mesh& mesh=*_meshes[msh];
4934  const BCcellArray& BCArray=*(_BCells[msh]);
4935  const BCfaceArray& BCfArray=*(_BFaces[msh]);
4936  const BCcellArray& ZCArray=*(_ZCells[msh]);
4937  COMETESBGKDiscretizer<T> CDisc(mesh,_geomFields,solidFaces,_macroFields,_quadrature,
4938  _dsfPtr,_dsfPtr1,_dsfPtr2,_dsfEqPtrES,_dsfPtrRes,_dsfPtrFAS,
4939  _options["timeStep"],_options.timeDiscretizationOrder,
4940  _options.transient,_options.underRelaxation,_options["rho_init"],
4941  _options["T_init"],_options["molecularWeight"],_options.conOrder,
4942  _bcMap,_faceReflectionArrayMap,BCArray,BCfArray,ZCArray);
4943 
4944  CDisc.setfgFinder();
4945 
4946  MakeParallel();
4947 
4948  CDisc.COMETSolve(1,_level); //forward
4949 
4950  MakeParallel();
4951 
4952  //callCOMETBoundaryConditions();
4953  ComputeCOMETMacroparameters();
4954  ComputeCollisionfrequency();
4955  //update equilibrium distribution function 0-maxwellian, 1-BGK,2-ESBGK
4956  if (_options.fgamma==0){initializeMaxwellianEq();}
4957  else{ EquilibriumDistributionBGK();}
4958  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
4959  computeSolidFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
4960  ConservationofMFSolid(solidFaces);
4961  computeIBFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
4962 
4963  CDisc.COMETSolve(-1,_level); //reverse
4964  if((num==1)||(num==0&&_level==0))
4965  {
4966  MakeParallel();
4967  //callCOMETBoundaryConditions();
4968  ComputeCOMETMacroparameters();
4969  ComputeCollisionfrequency();
4970  //update equilibrium distribution function 0-maxwellian, 1-BGK,2-ESBGK
4971  if (_options.fgamma==0){initializeMaxwellianEq();}
4972  else{ EquilibriumDistributionBGK();}
4973  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
4974  computeSolidFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
4975  ConservationofMFSolid(solidFaces);
4976  computeIBFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
4977  }
4978  }
4979  }
4980 
4981  void smooth(const int num)
4982  {
4983  const int numDir=_quadrature.getDirCount();
4984  const int numMeshes=_meshes.size();
4985  for(int msh=0;msh<numMeshes;msh++)
4986  {
4987  const Mesh& mesh=*_meshes[msh];
4988  const BCcellArray& BCArray=*(_BCells[msh]);
4989  const BCfaceArray& BCfArray=*(_BFaces[msh]);
4990  const BCcellArray& ZCArray=*(_ZCells[msh]);
4991  shared_ptr<StorageSite> solidFaces(new StorageSite(-1));
4992  COMETESBGKDiscretizer<T> CDisc(mesh,_geomFields,*solidFaces,_macroFields,_quadrature,
4993  _dsfPtr,_dsfPtr1,_dsfPtr2,_dsfEqPtrES,_dsfPtrRes,_dsfPtrFAS,
4994  _options["timeStep"],_options.timeDiscretizationOrder,
4995  _options.transient,_options.underRelaxation,_options["rho_init"],
4996  _options["T_init"],_options["molecularWeight"],_options.conOrder,
4997  _bcMap,_faceReflectionArrayMap,BCArray,BCfArray,ZCArray);
4998 
4999  CDisc.setfgFinder();
5000 
5001  if(_level==0)callCOMETBoundaryConditions();
5002  MakeParallel();
5003 
5004  if(_level==0)
5005  CDisc.COMETSolveFine(1,_level); //forward
5006  else
5007  CDisc.COMETSolve(1,_level); //forward
5008  //callCOMETBoundaryConditions();
5009  ComputeCOMETMacroparameters();
5010  ComputeCollisionfrequency();
5011  //update equilibrium distribution function 0-maxwellian, 1-BGK,2-ESBGK
5012  if (_options.fgamma==0){initializeMaxwellianEq();}
5013  else{ EquilibriumDistributionBGK();}
5014  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
5015 
5016  if(_level==0)callCOMETBoundaryConditions();
5017  MakeParallel();
5018 
5019 
5020  if(_level==0)
5021  CDisc.COMETSolveFine(-1,_level); //reverse
5022  else
5023  CDisc.COMETSolve(-1,_level); //forward
5024  if((num==1)||(num==0&&_level==0))
5025  {
5026  //callCOMETBoundaryConditions();
5027  ComputeCOMETMacroparameters();
5028  ComputeCollisionfrequency();
5029  //update equilibrium distribution function 0-maxwellian, 1-BGK,2-ESBGK
5030  if (_options.fgamma==0){initializeMaxwellianEq();}
5031  else{ EquilibriumDistributionBGK();}
5032  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
5033  }
5034 
5035  }
5036  }
5037 
5038  T updateResid(const bool addFAS)
5039  {
5040  const int numMeshes=_meshes.size();
5041  T lowResid=-1.;
5042  T currentResid;
5043  for(int msh=0;msh<numMeshes;msh++)
5044  {
5045  const Mesh& mesh=*_meshes[msh];
5046  const BCcellArray& BCArray=*(_BCells[msh]);
5047  const BCfaceArray& BCfArray=*(_BFaces[msh]);
5048  const BCcellArray& ZCArray=*(_ZCells[msh]);
5049  shared_ptr<StorageSite> solidFaces(new StorageSite(-1));
5050  COMETESBGKDiscretizer<T> CDisc(mesh,_geomFields,*solidFaces,_macroFields,_quadrature,
5051  _dsfPtr,_dsfPtr1,_dsfPtr2,_dsfEqPtrES,_dsfPtrRes,_dsfPtrFAS,
5052  _options["timeStep"],_options.timeDiscretizationOrder,
5053  _options.transient,_options.underRelaxation,_options["rho_init"],
5054  _options["T_init"],_options["molecularWeight"], _options.conOrder,
5055  _bcMap,_faceReflectionArrayMap,BCArray,BCfArray,ZCArray);
5056 
5057  CDisc.setfgFinder();
5058  const int numDir=_quadrature.getDirCount();
5059 
5060  if(_level==0)callCOMETBoundaryConditions();
5061  MakeParallel();
5062 
5063  if(_level==0)
5064  CDisc.findResidFine(addFAS);
5065  else
5066  CDisc.findResid(addFAS);
5067  currentResid=CDisc.getAveResid();
5068 
5069  if(lowResid<0)
5070  lowResid=currentResid;
5071  else
5072  if(currentResid<lowResid)
5073  lowResid=currentResid;
5074  }
5075  return lowResid;
5076  }
5077 
5078  T updateResid(const bool addFAS,const StorageSite& solidFaces)
5079  {
5080  const int numMeshes=_meshes.size();
5081  T lowResid=-1.;
5082  T currentResid;
5083  for(int msh=0;msh<numMeshes;msh++)
5084  {
5085  const Mesh& mesh=*_meshes[msh];
5086  const BCcellArray& BCArray=*(_BCells[msh]);
5087  const BCfaceArray& BCfArray=*(_BFaces[msh]);
5088  const BCcellArray& ZCArray=*(_ZCells[msh]);
5089  COMETESBGKDiscretizer<T> CDisc(mesh,_geomFields,solidFaces,_macroFields,_quadrature,
5090  _dsfPtr,_dsfPtr1,_dsfPtr2,_dsfEqPtrES,_dsfPtrRes,_dsfPtrFAS,
5091  _options["timeStep"],_options.timeDiscretizationOrder,
5092  _options.transient,_options.underRelaxation,_options["rho_init"],
5093  _options["T_init"],_options["molecularWeight"],_options.conOrder,
5094  _bcMap,_faceReflectionArrayMap,BCArray,BCfArray,ZCArray);
5095 
5096  CDisc.setfgFinder();
5097  const int numDir=_quadrature.getDirCount();
5098 
5099  MakeParallel();
5100 
5101  CDisc.findResid(addFAS);
5102  currentResid=CDisc.getAveResid();
5103 
5104  if(lowResid<0)
5105  lowResid=currentResid;
5106  else
5107  if(currentResid<lowResid)
5108  lowResid=currentResid;
5109  }
5110  return lowResid;
5111  }
5112 
5113  void cycle()
5114  {
5115  if(_level+1<_options.maxLevels)
5116  doSweeps(_options.preSweeps,1);
5117  else
5118  doSweeps(_options.preCoarsestSweeps,1);
5119  if(_level+1<_options.maxLevels)
5120  {
5121  if(_level==0)
5122  updateResid(false);
5123  else
5124  updateResid(true);
5125 
5126  injectResid();
5127  _coarserLevel->ComputeCOMETMacroparameters();
5128  _coarserLevel->ComputeCollisionfrequency();
5129  if (_options.fgamma==0){_coarserLevel->initializeMaxwellianEq();}
5130  else{_coarserLevel->EquilibriumDistributionBGK();}
5131  if (_options.fgamma==2){_coarserLevel->EquilibriumDistributionESBGK();}
5132 
5133  _coarserLevel->makeFAS();
5134  _coarserLevel->cycle();
5135  correctSolution();
5136 
5137  ComputeCOMETMacroparameters();
5138  ComputeCollisionfrequency();
5139  if (_options.fgamma==0){initializeMaxwellianEq();}
5140  else{EquilibriumDistributionBGK();}
5141  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
5142  }
5143  if(_level+1<_options.maxLevels)
5144  doSweeps(_options.postSweeps,0);
5145  else
5146  {
5147  if(_options.postCoarsestSweeps==1)
5148  doSweeps(_options.postCoarsestSweeps,0);
5149  else if(_options.postCoarsestSweeps>1)
5150  doSweeps(_options.postCoarsestSweeps,1);
5151  }
5152  }
5153 
5154 
5155  void cycle(const StorageSite& solidFaces)
5156  {
5157  if(_level+1<_options.maxLevels)
5158  doSweeps(_options.preSweeps,1,solidFaces);
5159  else
5160  doSweeps(_options.preCoarsestSweeps,1,solidFaces);
5161  if(_level+1<_options.maxLevels)
5162  {
5163  if(_level==0)
5164  updateResid(false,solidFaces);
5165  else
5166  updateResid(true,solidFaces);
5167 
5168  injectResid();
5169  _coarserLevel->ComputeCOMETMacroparameters();
5170  _coarserLevel->ComputeCollisionfrequency();
5171  if (_options.fgamma==0){_coarserLevel->initializeMaxwellianEq();}
5172  else{_coarserLevel->EquilibriumDistributionBGK();}
5173  if (_options.fgamma==2){_coarserLevel->EquilibriumDistributionESBGK();}
5174  _coarserLevel->MakeParallel();
5175  _coarserLevel->computeSolidFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
5176  _coarserLevel->ConservationofMFSolid(solidFaces);
5177  _coarserLevel->computeIBFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
5178 
5179  _coarserLevel->makeFAS(solidFaces);
5180  _coarserLevel->cycle(solidFaces);
5181  correctSolution();
5182 
5183  ComputeCOMETMacroparameters();
5184  ComputeCollisionfrequency();
5185  if (_options.fgamma==0){initializeMaxwellianEq();}
5186  else{EquilibriumDistributionBGK();}
5187  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
5188  MakeParallel();
5189  computeSolidFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
5190  ConservationofMFSolid(solidFaces);
5191  computeIBFaceDsf(solidFaces,_options.method,_options.relaxDistribution);
5192  }
5193  if(_level+1<_options.maxLevels)
5194  doSweeps(_options.postSweeps,0,solidFaces);
5195  else
5196  {
5197  if(_options.postCoarsestSweeps==1)
5198  doSweeps(_options.postCoarsestSweeps,0,solidFaces);
5199  else if(_options.postCoarsestSweeps>1)
5200  doSweeps(_options.postCoarsestSweeps,1,solidFaces);
5201  }
5202  }
5203 
5204  void injectResid()
5205  {
5206  const int numMeshes = _meshes.size();
5207 
5208  for (int n=0; n<numMeshes; n++)
5209  {
5210  const Mesh& finerMesh=*_meshes[n];
5211  const Mesh& coarserMesh=*(_coarserLevel->getMeshList())[n];
5212  MacroFields& coarserMacro=_coarserLevel->getMacro();
5213  GeomFields& coarserGeomFields=_coarserLevel->getGeomFields();
5214  DistFunctFields<T>& coarserdsf = _coarserLevel->getdsf();
5215  DistFunctFields<T>& coarserdsf0 = _coarserLevel->getdsf0();
5216  DistFunctFields<T>& coarserdsfFAS = _coarserLevel->getdsfFAS();
5217  const StorageSite& finerCells=finerMesh.getCells();
5218  const StorageSite& coarserCells=coarserMesh.getCells();
5219  const CRConnectivity& CoarserToFiner=coarserMesh.getConnectivity(coarserCells,finerCells);
5220  //const TArray& coarserVol=dynamic_cast<TArray&>(_geomFields.volume[coarserCells]);
5221  const TArray& coarserVol=dynamic_cast<TArray&>(coarserGeomFields.volume[coarserCells]);
5222  const TArray& finerVol=dynamic_cast<TArray&>(_geomFields.volume[finerCells]);
5223 
5224  const int cellCount=coarserCells.getSelfCount();
5225  const int numDir=_quadrature.getDirCount();
5226  for(int dir=0;dir<numDir;dir++)
5227  {
5228  Field& fnd = *_dsfPtr.dsf[dir];
5229  Field& fndRes = *_dsfPtrRes.dsf[dir];
5230  Field& cfnd = *coarserdsf.dsf[dir];
5231  Field& cfndInj = *coarserdsf0.dsf[dir];
5232  Field& cfndFAS = *coarserdsfFAS.dsf[dir];
5233  TArray& coarserVar=dynamic_cast<TArray&>(cfnd[coarserCells]);
5234  TArray& coarserInj=dynamic_cast<TArray&>(cfndInj[coarserCells]);
5235  TArray& coarserFAS=dynamic_cast<TArray&>(cfndFAS[coarserCells]);
5236  TArray& finerVar=dynamic_cast<TArray&>(fnd[finerCells]);
5237  TArray& finerRes=dynamic_cast<TArray&>(fndRes[finerCells]);
5238 
5239  for(int c=0;c<cellCount;c++)
5240  {
5241  const int fineCount=CoarserToFiner.getCount(c);
5242  coarserVar[c]=0.;
5243  coarserFAS[c]=0.;
5244 
5245  for(int fc=0;fc<fineCount;fc++)
5246  {
5247  const int cell=CoarserToFiner(c,fc);
5248  coarserVar[c]+=finerVar[cell]*finerVol[cell];
5249  coarserFAS[c]+=finerRes[cell];
5250  }
5251  coarserVar[c]/=coarserVol[c];
5252  coarserInj[c]=coarserVar[c];
5253  }
5254  }
5255 
5256  VectorT3Array& coarserVar=dynamic_cast<VectorT3Array&>(coarserMacro.velocity[coarserCells]);
5257  VectorT3Array& coarserInj=dynamic_cast<VectorT3Array&>(coarserMacro.velocityInjected[coarserCells]);
5258  VectorT3Array& coarserFAS=dynamic_cast<VectorT3Array&>(coarserMacro.velocityFASCorrection[coarserCells]);
5259  VectorT3Array& finerVar=dynamic_cast<VectorT3Array&>(_macroFields.velocity[finerCells]);
5260  VectorT3Array& finerRes=dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[finerCells]);
5261 
5262  for(int c=0;c<cellCount;c++)
5263  {
5264  const int fineCount=CoarserToFiner.getCount(c);
5265  coarserVar[c]=0.;
5266  coarserFAS[c]=0.;
5267 
5268  for(int fc=0;fc<fineCount;fc++)
5269  {
5270  const int cell=CoarserToFiner(c,fc);
5271  coarserVar[c]+=finerVar[cell]*finerVol[cell];
5272  coarserFAS[c]+=finerRes[cell];
5273  }
5274  coarserVar[c]/=coarserVol[c];
5275  coarserInj[c]=coarserVar[c];
5276  }
5277  }
5278  }
5279 
5280  void makeFAS()
5281  {
5282  updateResid(false);
5283 
5284  const int numMeshes = _meshes.size();
5285  for (int n=0; n<numMeshes; n++)
5286  {
5287  const Mesh& mesh=*_meshes[n];
5288  const StorageSite& cells=mesh.getCells();
5289 
5290  const int numDir = _quadrature.getDirCount();
5291  for(int dir=0;dir<numDir;dir++)
5292  {
5293  Field& fndRes = *_dsfPtrRes.dsf[dir];
5294  Field& fndFAS = *_dsfPtrFAS.dsf[dir];
5295  TArray& fRes = dynamic_cast<TArray&>(fndRes[cells]);
5296  TArray& fFAS = dynamic_cast<TArray&>(fndFAS[cells]);
5297 
5298  fFAS-=fRes;
5299  }
5300 
5301  VectorT3Array& vR = dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[cells]);
5302  VectorT3Array& vF = dynamic_cast<VectorT3Array&>(_macroFields.velocityFASCorrection[cells]);
5303 
5304  vF-=vR;
5305  }
5306  }
5307 
5308  void makeFAS(const StorageSite& solidFaces)
5309  {
5310  updateResid(false,solidFaces);
5311 
5312  const int numMeshes = _meshes.size();
5313  for (int n=0; n<numMeshes; n++)
5314  {
5315  const Mesh& mesh=*_meshes[n];
5316  const StorageSite& cells=mesh.getCells();
5317 
5318  const int numDir = _quadrature.getDirCount();
5319  for(int dir=0;dir<numDir;dir++)
5320  {
5321  Field& fndRes = *_dsfPtrRes.dsf[dir];
5322  Field& fndFAS = *_dsfPtrFAS.dsf[dir];
5323  TArray& fRes = dynamic_cast<TArray&>(fndRes[cells]);
5324  TArray& fFAS = dynamic_cast<TArray&>(fndFAS[cells]);
5325 
5326  fFAS-=fRes;
5327  }
5328 
5329  VectorT3Array& vR = dynamic_cast<VectorT3Array&>(_macroFields.velocityResidual[cells]);
5330  VectorT3Array& vF = dynamic_cast<VectorT3Array&>(_macroFields.velocityFASCorrection[cells]);
5331 
5332  vF-=vR;
5333  }
5334  }
5335 
5336  void correctSolution()
5337  {
5338  const int numMeshes = _meshes.size();
5339 
5340  for (int n=0; n<numMeshes; n++)
5341  {
5342  const Mesh& finerMesh=*_meshes[n];
5343  const Mesh& coarserMesh=*(_coarserLevel->getMeshList())[n];
5344  MacroFields& coarserMacro=_coarserLevel->getMacro();
5345  DistFunctFields<T>& coarserdsf = _coarserLevel->getdsf();
5346  DistFunctFields<T>& coarserdsf0 = _coarserLevel->getdsf0();
5347  const StorageSite& finerCells=finerMesh.getCells();
5348  const StorageSite& coarserCells=coarserMesh.getCells();
5349  const CRConnectivity& CoarserToFiner=coarserMesh.getConnectivity(coarserCells,finerCells);
5350 
5351  const int cellCount=coarserCells.getSelfCount();
5352 
5353  const int numDir=_quadrature.getDirCount();
5354  for(int dir=0;dir<numDir;dir++)
5355  {
5356  Field& fnd = *_dsfPtr.dsf[dir];
5357  Field& cfnd = *coarserdsf.dsf[dir];
5358  Field& cfndInj = *coarserdsf0.dsf[dir];
5359  TArray& finerArray = dynamic_cast<TArray&>(fnd[finerCells]);
5360  TArray& coarserArray = dynamic_cast<TArray&>(cfnd[coarserCells]);
5361  TArray& injArray = dynamic_cast<TArray&>(cfndInj[coarserCells]);
5362 
5363  for(int c=0;c<cellCount;c++)
5364  {
5365  const int fineCount=CoarserToFiner.getCount(c);
5366  const T correction=coarserArray[c]-injArray[c];
5367 
5368  for(int fc=0;fc<fineCount;fc++)
5369  finerArray[CoarserToFiner(c,fc)]+=correction;
5370  }
5371  }
5372 
5373  VectorT3Array& coarserArray=dynamic_cast<VectorT3Array&>(coarserMacro.velocity[coarserCells]);
5374  VectorT3Array& injArray=dynamic_cast<VectorT3Array&>(coarserMacro.velocityInjected[coarserCells]);
5375  VectorT3Array& finerArray=dynamic_cast<VectorT3Array&>(_macroFields.velocity[finerCells]);
5376 
5377  for(int c=0;c<cellCount;c++)
5378  {
5379  const int fineCount=CoarserToFiner.getCount(c);
5380  const VectorT3 correction=coarserArray[c]-injArray[c];
5381 
5382  for(int fc=0;fc<fineCount;fc++)
5383  finerArray[CoarserToFiner(c,fc)]+=correction;
5384  }
5385  }
5386  }
5387 
5388  void advance(const int iters)
5389  {
5390  callCOMETBoundaryConditions();
5391  _residual=updateResid(false);
5392  _initialResidual=_residual;
5393  T residualRatio(1.0);
5394 
5395 #ifdef FVM_PARALLEL
5396  if ( MPI::COMM_WORLD.Get_rank() == 0 )
5397  cout<<"Initial Residual:"<<_initialResidual<<" ResidualRatio: "<<residualRatio<<endl;
5398 #endif
5399 
5400 
5401 #ifndef FVM_PARALLEL
5402  cout << "Initial Residual:"<<_initialResidual<<" ResidualRatio: "<<residualRatio<<endl;
5403 #endif
5404 
5405 
5406 
5407  int niters=0;
5408  const T absTol=_options.absoluteTolerance;
5409  const T relTol=_options.relativeTolerance;
5410  const int show=_options.showResidual;
5411 
5412  while((niters<iters) && ((_residual>absTol)&&(residualRatio>relTol)))
5413  {
5414  cycle();
5415  niters++;
5416  _residual=updateResid(false);
5417  if(niters==1)
5418  _initialResidual=_residual;
5419  residualRatio=_residual/_initialResidual;
5420 #ifdef FVM_PARALLEL
5421  if((niters%show==0)&&(MPI::COMM_WORLD.Get_rank()==0))
5422  cout<<"Iteration:"<<niters<<" Residual:"<<_residual<<" ResidualRatio: "<<residualRatio<<endl;
5423 #endif
5424 
5425 #ifndef FVM_PARALLEL
5426  if(niters%show==0)
5427  cout<<"Iteration:"<<niters<<" Residual:"<<_residual<<" ResidualRatio: "<<residualRatio<<endl;
5428 #endif
5429  }
5430  callCOMETBoundaryConditions();
5431  //cout<<endl<<"Total Iterations:"<<niters<<" Residual:"<<_residual<<endl;
5432  }
5433 
5434  T advance(const int iters,const StorageSite& solidFaces)
5435  {
5436  callCOMETBoundaryConditions();
5437  _residual=updateResid(false,solidFaces);
5438  _initialResidual=_residual;
5439  T residualRatio(1.0);
5440 
5441 #ifdef FVM_PARALLEL
5442  if ( MPI::COMM_WORLD.Get_rank() == 0 )
5443  cout<<"Initial Residual:"<<_initialResidual<<" ResidualRatio: "<<residualRatio<<endl;
5444 #endif
5445 
5446 
5447 #ifndef FVM_PARALLEL
5448  cout << "Initial Residual:"<<_initialResidual<<" ResidualRatio: "<<residualRatio<<endl;
5449 #endif
5450 
5451 
5452 
5453  int niters=0;
5454  const T absTol=_options.absoluteTolerance;
5455  const T relTol=_options.relativeTolerance;
5456  const int show=_options.showResidual;
5457 
5458  while((niters<iters) && ((_residual>absTol)&&(residualRatio>relTol)))
5459  {
5460  cycle(solidFaces);
5461  niters++;
5462  _residual=updateResid(false,solidFaces);
5463  if(niters==1)
5464  _initialResidual=_residual;
5465  residualRatio=_residual/_initialResidual;
5466 #ifdef FVM_PARALLEL
5467  if((niters%show==0)&&(MPI::COMM_WORLD.Get_rank()==0))
5468  cout<<"Iteration:"<<niters<<" Residual:"<<_residual<<" ResidualRatio: "<<residualRatio<<endl;
5469 #endif
5470 
5471 #ifndef FVM_PARALLEL
5472  if(niters%show==0)
5473  cout<<"Iteration:"<<niters<<" Residual:"<<_residual<<" ResidualRatio: "<<residualRatio<<endl;
5474 #endif
5475  }
5476  callCOMETBoundaryConditions();
5477  return _residual;
5478  //cout<<endl<<"Total Iterations:"<<niters<<" Residual:"<<_residual<<endl;
5479  }
5480 
5481  void computeSurfaceForce(const StorageSite& solidFaces, bool perUnitArea, bool IBM=0)
5482  {
5483  typedef CRMatrixTranspose<T,T,T> IMatrix;
5484 
5485  const int nSolidFaces = solidFaces.getCount();
5486 
5487  boost::shared_ptr<VectorT3Array>
5488  forcePtr( new VectorT3Array(nSolidFaces));
5489  VectorT3Array& force = *forcePtr;
5490 
5491  force.zero();
5492  _macroFields.force.addArray(solidFaces,forcePtr);
5493 
5494  const VectorT3Array& solidFaceArea =
5495  dynamic_cast<const VectorT3Array&>(_geomFields.area[solidFaces]);
5496 
5497  const TArray& solidFaceAreaMag =
5498  dynamic_cast<const TArray&>(_geomFields.areaMag[solidFaces]);
5499 
5500  const int numMeshes = _meshes.size();
5501  for (int n=0; n<numMeshes; n++)
5502  {
5503  const Mesh& mesh = *_meshes[n];
5504  const StorageSite& cells = mesh.getCells();
5505 
5506  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
5507  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
5508  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
5509  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
5510  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
5511 
5512  const CRConnectivity& solidFacesToCells
5513  = mesh.getConnectivity(solidFaces,cells);
5514  const IntArray& sFCRow = solidFacesToCells.getRow();
5515  const IntArray& sFCCol = solidFacesToCells.getCol();
5516 
5517  const T Lx=_options["nonDimLx"];
5518  const T Ly=_options["nonDimLy"];
5519  const T Lz=_options["nonDimLz"];
5520 
5521 
5522  const int N123= _quadrature.getDirCount();
5523 
5524  const int selfCount = cells.getSelfCount();
5525  for(int f=0; f<nSolidFaces; f++){
5526 
5527  StressTensor<T> stress = NumTypeTraits<StressTensor<T> >::getZero();
5528 
5529  if (IBM){
5530  const VectorT3Array& vs = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[solidFaces]);
5531  GeomFields::SSPair key1(&solidFaces,&cells);
5532  const IMatrix& mIC =
5533  dynamic_cast<const IMatrix&>
5534  (*_geomFields._interpolationMatrices[key1]);
5535  const Array<T>& iCoeffs = mIC.getCoeff();
5536  for(int j=0;j<N123;j++){
5537  Field& fnd = *_dsfPtr.dsf[j];
5538  const TArray& f_dsf = dynamic_cast<const TArray&>(fnd[cells]);
5539  const TArray& fs_dsf = dynamic_cast<const TArray&>(fnd[solidFaces]);
5540  stress[0] -=pow((cx[j]-vs[f][0]),2.0)*fs_dsf[f]*wts[j];
5541  stress[1] -=pow((cy[j]-vs[f][1]),2.0)*fs_dsf[f]*wts[j];
5542  stress[2] -=pow((cz[j]-vs[f][2]),2.0)*fs_dsf[f]*wts[j];
5543  stress[3] -=(cx[j]-vs[f][0])*(cy[j]-vs[f][1])*fs_dsf[f]*wts[j];
5544  stress[4] -=(cy[j]-vs[f][1])*(cz[j]-vs[f][2])*fs_dsf[f]*wts[j];
5545  stress[5] -=(cx[j]-vs[f][0])*(cz[j]-vs[f][2])*fs_dsf[f]*wts[j];
5546  /*
5547  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
5548  {
5549 
5550  const int c = sFCCol[nc];
5551  const T coeff = iCoeffs[nc];
5552  stress[0] -=coeff*pow((cx[j]-v[c][0]),2.0)*f_dsf[c]*wts[j];
5553  stress[1] -=coeff*pow((cy[j]-v[c][1]),2.0)*f_dsf[c]*wts[j];
5554  stress[2] -=coeff*pow((cz[j]-v[c][2]),2.0)*f_dsf[c]*wts[j];
5555  stress[3] -=coeff*(cx[j]-v[c][0])*(cy[j]-v[c][1])*f_dsf[c]*wts[j];
5556  stress[4] -=coeff*(cy[j]-v[c][1])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
5557  stress[5] -=coeff*(cx[j]-v[c][0])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
5558  }
5559  */
5560  }
5561  }
5562  else
5563  {
5564  for(int j=0;j<N123;j++){
5565  Field& fnd = *_dsfPtr.dsf[j];
5566  const TArray& f_dsf = dynamic_cast<const TArray&>(fnd[cells]);
5567  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
5568  {
5569 
5570  const int c = sFCCol[nc];
5571  if ( c < selfCount ){
5572  stress[0] +=pow((cx[j]-v[c][0]),2.0)*f_dsf[c]*wts[j];
5573  stress[1] +=pow((cy[j]-v[c][1]),2.0)*f_dsf[c]*wts[j];
5574  stress[2] +=pow((cz[j]-v[c][2]),2.0)*f_dsf[c]*wts[j];
5575  stress[3] +=(cx[j]-v[c][0])*(cy[j]-v[c][1])*f_dsf[c]*wts[j];
5576  stress[4] +=(cy[j]-v[c][1])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
5577  stress[5] +=(cx[j]-v[c][0])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
5578  }
5579  }
5580  }
5581 
5582  }
5583 
5584 
5585  const VectorT3& Af = solidFaceArea[f];
5586  force[f][0] = Af[0]*Ly*Lz*stress[0] + Af[1]*Lz*Lx*stress[3] + Af[2]*Lx*Ly*stress[5];
5587  force[f][1] = Af[0]*Ly*Lz*stress[3] + Af[1]*Lz*Lx*stress[1] + Af[2]*Lx*Ly*stress[4];
5588  force[f][2] = Af[0]*Ly*Lz*stress[5] + Af[1]*Lz*Lx*stress[4] + Af[2]*Ly*Ly*stress[2];
5589  if (perUnitArea){
5590  force[f] /= solidFaceAreaMag[f];}
5591  }
5592  }
5593  }
5594 
5595  void OutputDsfPOINT() //, const char* filename)
5596  {
5597  FILE * pFile;
5598  pFile = fopen("cxyz0.plt","w");
5599  int N1=_quadrature.getNVCount();
5600  int N2=_quadrature.getNthetaCount();
5601  int N3=_quadrature.getNphiCount();
5602  fprintf(pFile,"%s \n", "VARIABLES= cx, cy, cz, f,");
5603  fprintf(pFile, "%s %i %s %i %s %i \n","ZONE I=", N3,",J=",N2,"K=",N1);
5604  fprintf(pFile,"%s\n","F=POINT");
5605  const int numMeshes = _meshes.size();
5606  const int cellno=_options.printCellNumber;
5607  for (int n=0; n<numMeshes; n++)
5608  {
5609  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
5610  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
5611  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
5612  const int numFields= _quadrature.getDirCount();
5613 
5614  const Mesh& mesh = *_meshes[n];
5615  const StorageSite& cells = mesh.getCells();
5616  for(int j=0;j< numFields;j++)
5617  {
5618  Field& fnd = *_dsfPtr.dsf[j];
5619  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
5620  Field& fEqnd = *_dsfEqPtr.dsf[j];
5621  TArray& fEq = dynamic_cast< TArray&>(fEqnd[cells]);
5622  fprintf(pFile,"%E %E %E %E %E\n",cx[j],cy[j],cz[j],f[cellno],fEq[cellno]);
5623  }
5624  }
5625  }
5626 
5627  DistFunctFields<T>& getdsf() { return _dsfPtr;}
5628  const DistFunctFields<T>& getdsf1() const { return _dsfPtr1;}
5629  const DistFunctFields<T>& getdsf2() const { return _dsfPtr2;}
5630  const DistFunctFields<T>& getdsfEq() const { return _dsfEqPtr;}
5631  const DistFunctFields<T>& getdsfEqES() const { return _dsfEqPtrES;}
5632  DistFunctFields<T>& getdsf0() { return _dsfPtr0;}
5633  DistFunctFields<T>& getdsfInj() { return _dsfPtrInj;}
5634  DistFunctFields<T>& getdsfRes() { return _dsfPtrRes;}
5635  DistFunctFields<T>& getdsfFAS() { return _dsfPtrFAS;}
5636  void setBCMap(COMETBCMap* bcMap) {_bcMap=*bcMap;}
5637  void setCoarserLevel(TCOMET* cl) {_coarserLevel=cl;}
5638  void setFinerLevel(TCOMET* fl) {_finerLevel=fl;}
5639  int getLevel() {return _level;}
5640  const MeshList& getMeshList() {return _meshes;}
5641  GeomFields& getGeomFields() {return _geomFields;}
5642  TQuad& getQuadrature() {return _quadrature;}
5643  MacroFields& getMacro() {return _macroFields;}
5644  T getResidual() {return _residual;}
5645 
5646 
5647 
5648  private:
5649  //shared_ptr<Impl> _impl;
5650 
5651  const int _level;
5652  GeomFields& _geomFields;
5653  GeomFields _coarseGeomFields;
5654  GeomFields& _finestGeomFields;
5655  const MeshList& _finestMeshes;
5656  Quadrature<T>& _quadrature;
5657  const int _ibm;
5658 
5659  MacroFields& _macroFields;
5660  MacroFields& _finestMacroFields;
5661  TCOMET* _finestLevel;
5662  TCOMET* _coarserLevel;
5663  TCOMET* _finerLevel;
5664  DistFunctFields<T> _dsfPtr;
5665  DistFunctFields<T> _dsfPtr1;
5666  DistFunctFields<T> _dsfPtr2;
5667  DistFunctFields<T> _dsfEqPtr;
5668  DistFunctFields<T> _dsfEqPtrES;
5669  DistFunctFields<T> _dsfPtr0;
5670  DistFunctFields<T> _dsfPtrInj;
5671  DistFunctFields<T> _dsfPtrRes;
5672  DistFunctFields<T> _dsfPtrFAS;
5673 
5674  COMETBCMap _bcMap;
5675  COMETVCMap _vcMap;
5676 
5677  COMETModelOptions<T> _options;
5678  T _residual;
5679  T _initialResidual;
5680 
5681 
5682  MFRPtr _initialKmodelNorm;
5683  BCcellList _BCells;
5684  BCfaceList _BFaces;
5685  BCcellList _ZCells;
5686  int _niters;
5687  map<int, vector<int> > _faceReflectionArrayMap;
5688  map<string,shared_ptr<ArrayBase> > _persistenceData;
5689 
5690  //MatrixSizeMap _coarseSizes;
5691  //MatrixSizeMap _coarseGhostSizes;
5692  SizeMap _coarseSizes;
5693  SizeMap _coarseGhostSizes;
5694  SiteMap _siteMap;
5695  GhostStorageSiteMap _sharedSiteMap;
5696  MatrixMappersMap _coarseScatterMaps;
5697  MatrixMappersMap _coarseGatherMaps;
5698 };
5699 
5700 #endif
COMETModel::_finestMeshes
const MeshList & _finestMeshes
Definition: esbgkbase/COMETModel.h:5655
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
epsilon
#define epsilon
Definition: Mesh.cpp:17
CRConnectivity::getCol
const Array< int > & getCol() const
Definition: CRConnectivity.h:183
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
COMETModel::getPersistenceData
map< string, shared_ptr< ArrayBase > > & getPersistenceData()
Definition: esbgkbase/COMETModel.h:2268
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
COMETModel::setJacobianBGK
void setJacobianBGK(SquareMatrix< T, 5 > &fjac, VectorT5 &fvec, const VectorT5 &xn, const VectorT3 &v, const int c)
Definition: esbgkbase/COMETModel.h:1799
COMETModel::SizeMap
map< const Mesh *, int > SizeMap
Definition: esbgkbase/COMETModel.h:101
COMETModel::BCcellArray
Array< int > BCcellArray
Definition: esbgkbase/COMETModel.h:84
COMETModel::_dsfPtrInj
DistFunctFields< T > _dsfPtrInj
Definition: esbgkbase/COMETModel.h:5670
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
COMETModel::computeIBFaceDsf
void computeIBFaceDsf(const StorageSite &solidFaces, const int method, const int RelaxDistribution=0)
Definition: esbgkbase/COMETModel.h:3473
COMETModel::_geomFields
GeomFields & _geomFields
Definition: esbgkbase/COMETModel.h:5652
COMETModel::updateResid
T updateResid(const bool addFAS)
Definition: esbgkbase/COMETModel.h:5038
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
COMETModel::EntryIndex
pair< Index, Index > EntryIndex
Definition: esbgkbase/COMETModel.h:96
COMETModel::OutputDsfBLOCK
void OutputDsfBLOCK(const char *filename)
Definition: esbgkbase/COMETModel.h:2493
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
Quadrature::CopyQuad
void CopyQuad(TQuad &copyFromQuad)
Definition: Quadrature.h:442
quad
static Cell< Quad > quad
Definition: Cell.cpp:210
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
MacroFields::coeffg
Field coeffg
Definition: MacroFields.h:32
COMETModel::getdsfInj
DistFunctFields< T > & getdsfInj()
Definition: esbgkbase/COMETModel.h:5633
COMETModel::TCOMETPtr
shared_ptr< TCOMET > TCOMETPtr
Definition: esbgkbase/COMETModel.h:93
COMETESBGKDiscretizer::COMETSolveFine
void COMETSolveFine(const int sweep, const int level)
Definition: COMETESBGKDiscretizer.h:153
COMETModel::ConservationofMassCheck
const double ConservationofMassCheck()
Definition: esbgkbase/COMETModel.h:4808
COMETModel::smooth
void smooth(const int num)
Definition: esbgkbase/COMETModel.h:4981
COMETModel::BCellPtr
shared_ptr< BCcellArray > BCellPtr
Definition: esbgkbase/COMETModel.h:85
COMETModel::COMETModel
COMETModel(const MeshList &meshes, const int level, GeomFields &geomFields, MacroFields &macroFields, Quadrature< T > &quad, const int ibm=0, GeomFields *finestGeomFields=NULL, const MeshList *finestMeshes=NULL, MacroFields *finestMacroFields=NULL)
Definition: esbgkbase/COMETModel.h:115
COMETModel::COMETBCMap
std::map< int, COMETBC< T > * > COMETBCMap
Definition: esbgkbase/COMETModel.h:90
COMETModel::initializeCoarseMaxwellian
void initializeCoarseMaxwellian()
Definition: esbgkbase/COMETModel.h:2146
COMETVC::vcType
string vcType
Definition: esbgkbase/COMETBC.h:43
COMETModel::stdVectorField
std::vector< Field * > stdVectorField
Definition: esbgkbase/COMETModel.h:65
FaceGroup
Definition: Mesh.h:28
COMETModel::_dsfPtr2
DistFunctFields< T > _dsfPtr2
Definition: esbgkbase/COMETModel.h:5666
COMETModel::makeFAS
void makeFAS(const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:5308
COMETModel::TArrptr
shared_ptr< TArray > TArrptr
Definition: esbgkbase/COMETModel.h:57
COMETModel::weightedMaxwellian
void weightedMaxwellian(double weight1, double uvel1, double vvel1, double wvel1, double uvel2, double vvel2, double wvel2, double temp1, double temp2)
Definition: esbgkbase/COMETModel.h:2174
COMETModel::_quadrature
Quadrature< T > & _quadrature
Definition: esbgkbase/COMETModel.h:5656
COMETModel::setJacobianESBGK
void setJacobianESBGK(SquareMatrix< T, 10 > &fjac, VectorT10 &fvec, const VectorT10 &xn, const VectorT3 &v, const int c)
Definition: esbgkbase/COMETModel.h:2008
COMETModel::TArray
Array< T > TArray
Definition: esbgkbase/COMETModel.h:56
COMETModel::MshLstPtr
shared_ptr< MeshList > MshLstPtr
Definition: esbgkbase/COMETModel.h:75
COMETModel::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: esbgkbase/COMETModel.h:54
COMETModel::getGeomFields
GeomFields & getGeomFields()
Definition: esbgkbase/COMETModel.h:5641
COMETModel::_coarseGatherMaps
MatrixMappersMap _coarseGatherMaps
Definition: esbgkbase/COMETModel.h:5697
Field
Definition: Field.h:14
COMETModel::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: esbgkbase/COMETModel.h:97
Mesh::createInteriorFaceGroup
const StorageSite & createInteriorFaceGroup(const int size)
Definition: Mesh.cpp:259
COMETModel::getdsfRes
DistFunctFields< T > & getdsfRes()
Definition: esbgkbase/COMETModel.h:5634
MacroFields::velocityInjected
Field velocityInjected
Definition: MacroFields.h:17
StorageSite::getGatherMapLevel1
const GatherMap & getGatherMapLevel1() const
Definition: StorageSite.h:74
StorageSite::getScatterMapLevel1
const ScatterMap & getScatterMapLevel1() const
Definition: StorageSite.h:73
MacroFields::Stress
Field Stress
Definition: MacroFields.h:37
COMETESBGKDiscretizer::findResid
void findResid(const bool plusFAS)
Definition: COMETESBGKDiscretizer.h:1707
Mesh::getIBFaceList
const Array< int > & getIBFaceList() const
Definition: Mesh.cpp:686
MacroFields::EntropyGenRate_Collisional
Field EntropyGenRate_Collisional
Definition: MacroFields.h:35
COMETModel::computeSolidFaceDsf
void computeSolidFaceDsf(const StorageSite &solidFaces, const int method, const int RelaxDistribution=0)
Definition: esbgkbase/COMETModel.h:4154
Mesh
Definition: Mesh.h:49
COMETBoundaryConditions::applyZeroGradientBC
void applyZeroGradientBC(int f, GradMatrix &gMat) const
Definition: COMETBoundaryConditions.h:504
COMETModel::BCcellList
vector< BCellPtr > BCcellList
Definition: esbgkbase/COMETModel.h:86
COMETModel::GhostStorageSiteMap
map< const StorageSite *, shared_ptr< StorageSite > > GhostStorageSiteMap
Definition: esbgkbase/COMETModel.h:107
mag
T mag(const Vector< T, 3 > &a)
Definition: Vector.h:260
MacroFields::InterfaceVelocity
Field InterfaceVelocity
Definition: MacroFields.h:39
COMETModel::_options
COMETModelOptions< T > _options
Definition: esbgkbase/COMETModel.h:5677
StorageSite::setCount
void setCount(const int selfCount, const int nGhost=0)
Definition: StorageSite.h:42
UnorderedSet::insert
void insert(int x)
Definition: MultiFieldMatrix.h:25
COMETModel::setBCMap
void setBCMap(COMETBCMap *bcMap)
Definition: esbgkbase/COMETModel.h:5636
Mesh::setConnectivity
void setConnectivity(const StorageSite &rowSite, const StorageSite &colSite, shared_ptr< CRConnectivity > conn)
Definition: Mesh.cpp:352
MacroFields::EntropyGenRate
Field EntropyGenRate
Definition: MacroFields.h:34
MacroFields::Entropy
Field Entropy
Definition: MacroFields.h:33
COMETESBGKDiscretizer::COMETSolve
void COMETSolve(const int sweep, const int level)
Definition: COMETESBGKDiscretizer.h:231
COMETModel::initializeFineMaxwellian
void initializeFineMaxwellian()
Definition: esbgkbase/COMETModel.h:2110
COMETModel::TDistFF
DistFunctFields< T > TDistFF
Definition: esbgkbase/COMETModel.h:66
COMETModel::VectorT6
Vector< T, 6 > VectorT6
Definition: esbgkbase/COMETModel.h:70
MacroFields::InterfacePressure
Field InterfacePressure
Definition: MacroFields.h:40
MacroFields::velocityResidual
Field velocityResidual
Definition: MacroFields.h:16
Mesh::createBoundaryFaceGroup
const StorageSite & createBoundaryFaceGroup(const int size, const int offset, const int id, const string &boundaryType)
Definition: Mesh.cpp:278
COMETModel::ComputeCollisionfrequency
void ComputeCollisionfrequency()
Definition: esbgkbase/COMETModel.h:1486
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
COMETModel::SiteMap
map< const StorageSite *, StorageSite * > SiteMap
Definition: esbgkbase/COMETModel.h:102
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
COMETModel::BfacePtr
shared_ptr< BCfaceArray > BfacePtr
Definition: esbgkbase/COMETModel.h:82
FloatVarDict::getVal
FloatVal< T > getVal(const string varName) const
Definition: FloatVarDict.h:85
COMETModel::setFinerLevel
void setFinerLevel(TCOMET *fl)
Definition: esbgkbase/COMETModel.h:5638
COMETModel::TCOMET
COMETModel< T > TCOMET
Definition: esbgkbase/COMETModel.h:92
COMETModel::getdsfFAS
DistFunctFields< T > & getdsfFAS()
Definition: esbgkbase/COMETModel.h:5635
COMETModel::doSweeps
void doSweeps(const int sweeps, const int num, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:4921
COMETModel::ComputeFineMacroparameters
void ComputeFineMacroparameters()
Definition: esbgkbase/COMETModel.h:1318
COMETModel::GeoFldsPtr
shared_ptr< GeomFields > GeoFldsPtr
Definition: esbgkbase/COMETModel.h:77
MacroFields::viscosity
Field viscosity
Definition: MacroFields.h:20
COMETModel::_dsfPtrFAS
DistFunctFields< T > _dsfPtrFAS
Definition: esbgkbase/COMETModel.h:5672
GeomFields::fineToCoarse
Field fineToCoarse
Definition: GeomFields.h:41
COMETModel::correctMassDeficit2
void correctMassDeficit2(double n1, double n2)
Definition: esbgkbase/COMETModel.h:4749
sqrt
Tangent sqrt(const Tangent &a)
Definition: Tangent.h:317
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
MacroFields::temperature
Field temperature
Definition: MacroFields.h:22
COMETModel::initializeMaxwellianEq
void initializeMaxwellianEq()
Definition: esbgkbase/COMETModel.h:1620
COMETModel::syncGhostCoarsening
void syncGhostCoarsening(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes)
Definition: esbgkbase/COMETModel.h:3214
Mesh::getCellFaces
const CRConnectivity & getCellFaces() const
Definition: Mesh.cpp:454
COMETModel::ComputeCOMETMacroparameters
void ComputeCOMETMacroparameters()
Definition: esbgkbase/COMETModel.h:1235
FaceGroup::id
const int id
Definition: Mesh.h:41
COMETModel::_faceReflectionArrayMap
map< int, vector< int > > _faceReflectionArrayMap
Definition: esbgkbase/COMETModel.h:5687
COMETModel::MatrixMappersMap
map< SSPair, shared_ptr< Array< int > > > MatrixMappersMap
Definition: esbgkbase/COMETModel.h:104
COMETBC::bcType
string bcType
Definition: esbgkbase/COMETBC.h:32
COMETModel::getdsf2
const DistFunctFields< T > & getdsf2() const
Definition: esbgkbase/COMETModel.h:5629
MacroFields::force
Field force
Definition: MacroFields.h:36
COMETModel::getOptions
COMETModelOptions< T > & getOptions()
Definition: esbgkbase/COMETModel.h:2261
COMETModel::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: esbgkbase/COMETModel.h:61
COMETBoundaryConditions::applyPressureOutletBC
void applyPressureOutletBC(int f, const X &outletTemperature, const X &outletPressure) const
Definition: COMETBoundaryConditions.h:258
COMETModel::getFaceReflectionArrayMap
const map< int, vector< int > > & getFaceReflectionArrayMap() const
Definition: esbgkbase/COMETModel.h:2262
Mesh::createInterfaceGroup
const StorageSite & createInterfaceGroup(const int size, const int offset, const int id)
Definition: Mesh.cpp:268
Array::getData
virtual void * getData() const
Definition: Array.h:275
StorageSite::getTag
int getTag() const
Definition: StorageSite.h:84
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
COMETModel::BCfaceList
vector< BfacePtr > BCfaceList
Definition: esbgkbase/COMETModel.h:83
StorageSite::getScatterProcID
int getScatterProcID() const
Definition: StorageSite.h:82
COMETModel::_finestGeomFields
GeomFields & _finestGeomFields
Definition: esbgkbase/COMETModel.h:5654
StorageSite::getGatherProcID
int getGatherProcID() const
Definition: StorageSite.h:83
inverseGauss
SquareMatrix< T, N > inverseGauss(SquareMatrix< T, N > &a, int size)
Definition: MatrixOperation.h:201
COMETModel::_sharedSiteMap
GhostStorageSiteMap _sharedSiteMap
Definition: esbgkbase/COMETModel.h:5695
COMETModel::TArray2D
Array2D< T > TArray2D
Definition: esbgkbase/COMETModel.h:59
COMETModel::SSPtr
shared_ptr< StorageSite > SSPtr
Definition: esbgkbase/COMETModel.h:78
COMETModel::getdsf
DistFunctFields< T > & getdsf()
Definition: esbgkbase/COMETModel.h:5627
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
MacroFields::coeff
Field coeff
Definition: MacroFields.h:31
COMETModel::NewtonsMethodBGK
void NewtonsMethodBGK(const int ktrial)
Definition: esbgkbase/COMETModel.h:1670
COMETModel::NewtonsMethodESBGK
void NewtonsMethodESBGK(const int ktrial)
Definition: esbgkbase/COMETModel.h:1840
COMETModel::StorageSiteMap
map< Index, shared_ptr< StorageSite > > StorageSiteMap
Definition: esbgkbase/COMETModel.h:106
COMETESBGKDiscretizer::findResidFine
void findResidFine(const bool plusFAS)
Definition: COMETESBGKDiscretizer.h:1573
COMETModel::updateResid
T updateResid(const bool addFAS, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:5078
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
COMETModel::getVCMap
COMETVCMap & getVCMap()
Definition: esbgkbase/COMETModel.h:2259
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
COMETModel::computeSolidFacePressure
void computeSolidFacePressure(const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:4117
COMETModel::_macroFields
MacroFields & _macroFields
Definition: esbgkbase/COMETModel.h:5659
COMETModel::Index
MultiField::ArrayIndex Index
Definition: esbgkbase/COMETModel.h:95
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
COMETModel::StressTensorArray
Array< StressTensorT6 > StressTensorArray
Definition: esbgkbase/COMETModel.h:64
COMETModel::MakeCoarseMesh1
void MakeCoarseMesh1(const MeshList &inMeshes, GeomFields &inGeomFields, MeshList &outMeshes)
Definition: esbgkbase/COMETModel.h:2545
Model
Definition: Model.h:13
UnorderedSet::size
int size() const
Definition: MultiFieldMatrix.h:34
GeomFields::ibFaceIndex
Field ibFaceIndex
Definition: GeomFields.h:40
COMETModel::getMacro
MacroFields & getMacro()
Definition: esbgkbase/COMETModel.h:5643
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::volume
Field volume
Definition: GeomFields.h:26
COMETModel::ComputeMacroparameters
void ComputeMacroparameters()
Definition: esbgkbase/COMETModel.h:1143
COMETModel::_persistenceData
map< string, shared_ptr< ArrayBase > > _persistenceData
Definition: esbgkbase/COMETModel.h:5688
COMETModel::computeSurfaceForce
void computeSurfaceForce(const StorageSite &solidFaces, bool perUnitArea, bool IBM=0)
Definition: esbgkbase/COMETModel.h:5481
COMETModel::MakeIBCoarseModel
void MakeIBCoarseModel(TCOMET *finerModel, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:753
COMETModel::EquilibriumDistributionESBGK
void EquilibriumDistributionESBGK()
Definition: esbgkbase/COMETModel.h:1953
COMETModel::ComputeMacroparametersESBGK
void ComputeMacroparametersESBGK()
Definition: esbgkbase/COMETModel.h:1406
COMETModel::VectorT6Array
Array< VectorT6 > VectorT6Array
Definition: esbgkbase/COMETModel.h:71
COMETModel::VectorT3
Vector< T, 3 > VectorT3
Definition: esbgkbase/COMETModel.h:60
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
COMETModel::cycle
void cycle(const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:5155
COMETModel::CRPtr
shared_ptr< CRConnectivity > CRPtr
Definition: esbgkbase/COMETModel.h:79
COMETModel::setCoarserLevel
void setCoarserLevel(TCOMET *cl)
Definition: esbgkbase/COMETModel.h:5637
StorageSite::getScatterMap
const ScatterMap & getScatterMap() const
Definition: StorageSite.h:58
COMETModel::getMeshList
const MeshList & getMeshList()
Definition: esbgkbase/COMETModel.h:5640
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
COMETModel::callCOMETBoundaryConditions
void callCOMETBoundaryConditions()
Definition: esbgkbase/COMETModel.h:2426
COMETModel::VectorT5
Vector< T, 5 > VectorT5
Definition: esbgkbase/COMETModel.h:68
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
COMETModel::_dsfEqPtrES
DistFunctFields< T > _dsfEqPtrES
Definition: esbgkbase/COMETModel.h:5668
COMETModel::BArray
Array< bool > BArray
Definition: esbgkbase/COMETModel.h:58
StorageSite::ScatterMap
map< const StorageSite *, shared_ptr< Array< int > > > ScatterMap
Definition: StorageSite.h:23
UnorderedSet::getData
const vector< int > & getData() const
Definition: MultiFieldMatrix.h:33
COMETModel::_finestMacroFields
MacroFields & _finestMacroFields
Definition: esbgkbase/COMETModel.h:5660
COMETModel::VectorT10Array
Array< VectorT10 > VectorT10Array
Definition: esbgkbase/COMETModel.h:73
COMETModel::VT3Ptr
shared_ptr< VectorT3Array > VT3Ptr
Definition: esbgkbase/COMETModel.h:62
COMETModel::StressTensorT6
StressTensor< T > StressTensorT6
Definition: esbgkbase/COMETModel.h:63
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
MacroFields::velocity
Field velocity
Definition: MacroFields.h:15
COMETModel::SetBoundaryConditions
void SetBoundaryConditions()
Definition: esbgkbase/COMETModel.h:2313
StorageSite::getOffset
int getOffset() const
Definition: StorageSite.h:87
COMETModel::advance
T advance(const int iters, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:5434
MacroFields::pressure
Field pressure
Definition: MacroFields.h:19
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getGatherMap
const GatherMap & getGatherMap() const
Definition: StorageSite.h:59
MacroFields::InterfaceStress
Field InterfaceStress
Definition: MacroFields.h:41
COMETModel::MatrixMap
map< EntryIndex, shared_ptr< Matrix > > MatrixMap
Definition: esbgkbase/COMETModel.h:99
GeomFields::finestToCoarse
Field finestToCoarse
Definition: GeomFields.h:42
COMETModel::_dsfPtr0
DistFunctFields< T > _dsfPtr0
Definition: esbgkbase/COMETModel.h:5669
MacroFields::velocityFASCorrection
Field velocityFASCorrection
Definition: MacroFields.h:18
DistFunctFields::dsf
std::vector< Field * > dsf
Definition: DistFunctFields.h:36
COMETModel::MakeCoarseMesh2
int MakeCoarseMesh2(const MeshList &inMeshes, GeomFields &inGeomFields, GeomFields &coarseGeomFields, MeshList &outMeshes)
Definition: esbgkbase/COMETModel.h:2786
COMETBoundaryConditions::applyPressureInletBC
void applyPressureInletBC(int f, const X &inletTemperature, const X &inletPressure, GradMatrix &gMat) const
Definition: COMETBoundaryConditions.h:77
Field::syncLocalVectorFields
static void syncLocalVectorFields(std::vector< Field * > &dsf)
Definition: Field.cpp:702
COMETModel::_dsfPtrRes
DistFunctFields< T > _dsfPtrRes
Definition: esbgkbase/COMETModel.h:5671
COMETModel::getQuadrature
TQuad & getQuadrature()
Definition: esbgkbase/COMETModel.h:5642
COMETModel::_dsfPtr
DistFunctFields< T > _dsfPtr
Definition: esbgkbase/COMETModel.h:5664
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
COMETModel::_dsfPtr1
DistFunctFields< T > _dsfPtr1
Definition: esbgkbase/COMETModel.h:5665
COMETModel::VectorT10
Vector< T, 10 > VectorT10
Definition: esbgkbase/COMETModel.h:72
COMETModel::getdsfEq
const DistFunctFields< T > & getdsfEq() const
Definition: esbgkbase/COMETModel.h:5630
MFRPtr
shared_ptr< MultiFieldReduction > MFRPtr
Definition: MultiFieldReduction.h:56
COMETModel::IntArray
Array< int > IntArray
Definition: esbgkbase/COMETModel.h:55
GeomFields::area
Field area
Definition: GeomFields.h:23
COMETModel::MatrixSizeMap
map< Index, int > MatrixSizeMap
Definition: esbgkbase/COMETModel.h:100
StorageSite::GatherMap
map< const StorageSite *, shared_ptr< Array< int > > > GatherMap
Definition: StorageSite.h:24
Mesh::getInterfaceGroups
const FaceGroupList & getInterfaceGroups() const
Definition: Mesh.h:190
COMETModel::_coarseGeomFields
GeomFields _coarseGeomFields
Definition: esbgkbase/COMETModel.h:5653
MacroFields::collisionFrequency
Field collisionFrequency
Definition: MacroFields.h:24
COMETModel::getBCMap
COMETBCMap & getBCMap()
Definition: esbgkbase/COMETModel.h:2258
COMETModel::_dsfEqPtr
DistFunctFields< T > _dsfEqPtr
Definition: esbgkbase/COMETModel.h:5667
COMETModel::smooth
void smooth(const int num, const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:4927
COMETModel::InitializeMacroparameters
void InitializeMacroparameters()
Definition: esbgkbase/COMETModel.h:1021
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
COMETModel::_coarseScatterMaps
MatrixMappersMap _coarseScatterMaps
Definition: esbgkbase/COMETModel.h:5696
StorageSite::setCountLevel1
void setCountLevel1(const int countLevel1)
Definition: StorageSite.h:69
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
COMETModel::getdsfEqES
const DistFunctFields< T > & getdsfEqES() const
Definition: esbgkbase/COMETModel.h:5631
COMETModel::getdsf0
DistFunctFields< T > & getdsf0()
Definition: esbgkbase/COMETModel.h:5632
COMETModel::EquilibriumDistributionBGK
void EquilibriumDistributionBGK()
Definition: esbgkbase/COMETModel.h:1748
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
COMETModel::COMETVCMap
std::map< int, COMETVC< T > * > COMETVCMap
Definition: esbgkbase/COMETModel.h:91
COMETModel::InitializeFgammaCoefficients
void InitializeFgammaCoefficients()
Definition: esbgkbase/COMETModel.h:1089
Mesh::getID
int getID() const
Definition: Mesh.h:106
COMETModel::MeshPtr
shared_ptr< Mesh > MeshPtr
Definition: esbgkbase/COMETModel.h:76
COMETModel::ComputeCoarseMacroparameters
void ComputeCoarseMacroparameters()
Definition: esbgkbase/COMETModel.h:1353
COMETModel::advance
void advance(const int iters)
Definition: esbgkbase/COMETModel.h:5388
COMETModel::MakeCoarseIndex
void MakeCoarseIndex(const StorageSite &solidFaces)
Definition: esbgkbase/COMETModel.h:970
MacroFields::density
Field density
Definition: MacroFields.h:21
COMETModel::MakeCoarseModel
void MakeCoarseModel(TCOMET *finerModel)
Definition: esbgkbase/COMETModel.h:575
COMETModel::BCfaceArray
Array< int > BCfaceArray
Definition: esbgkbase/COMETModel.h:81
COMETModel::doSweeps
void doSweeps(const int sweeps, const int num)
Definition: esbgkbase/COMETModel.h:4915
COMETModel::getdsf1
const DistFunctFields< T > & getdsf1() const
Definition: esbgkbase/COMETModel.h:5628
COMETModel::TQuad
Quadrature< T > TQuad
Definition: esbgkbase/COMETModel.h:88
COMETModel::ConservationofMFSolid
void ConservationofMFSolid(const StorageSite &solidFaces) const
Definition: esbgkbase/COMETModel.h:4846
COMETModel::weightedMaxwellian
void weightedMaxwellian(double weight1, double uvel1, double uvel2, double temp1, double temp2)
Definition: esbgkbase/COMETModel.h:2214
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
COMETModel::VectorT5Array
Array< VectorT5 > VectorT5Array
Definition: esbgkbase/COMETModel.h:69
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
MacroFields::InterfaceDensity
Field InterfaceDensity
Definition: MacroFields.h:42
Array::getLength
int getLength() const
Definition: Array.h:87