Memosa-FVM  0.2
KineticModel.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #ifndef _KINETICMODEL_H_
6 #define _KINETICMODEL_H_
7 
8 #ifdef FVM_PARALLEL
9 #include <mpi.h>
10 #endif
11 
12 #include <stdio.h>
13 #include <map>
14 #include <cmath>
15 #include <vector>
16 #include "Model.h"
17 #include "Array.h"
18 #include "Vector.h"
19 
20 #include "Mesh.h"
21 
22 #include "Quadrature.h"
23 #include "DistFunctFields.h"
24 
25 #include "MacroFields.h"
26 #include "FlowFields.h"
27 
28 #include "KineticBC.h"
29 #include "TimeDerivativeDiscretization_Kmodel.h"
30 #include "CollisionTermDiscretization.h"
31 #include "ConvectionDiscretization_Kmodel.h"
32 #include "KineticBoundaryConditions.h"
33 #include "GenericKineticBCS.h"
34 #include "GenericKineticIBDiscretization.h"
35 
36 #include "Linearizer.h"
37 #include "CRConnectivity.h"
38 #include "LinearSystem.h"
39 #include "MultiFieldMatrix.h"
40 #include "CRMatrix.h"
41 #include "CRMatrixTranspose.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 KineticModel : 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 Array2D<T> TArray2D;
58  typedef Vector<T,3> VectorT3;
59  typedef Array<VectorT3> VectorT3Array;
60  typedef StressTensor<T> StressTensorT6;
61  typedef Array<StressTensorT6> StressTensorArray;
62  typedef std::vector<Field*> stdVectorField;
63  typedef DistFunctFields<T> TDistFF;
64 
65  typedef Vector<T,5> VectorT5;
66  typedef Array<VectorT5> VectorT5Array;
67  typedef Vector<T,6> VectorT6;
68  typedef Array<VectorT6> VectorT6Array;
69  typedef Vector<T,10> VectorT10;
70  typedef Array<VectorT10> VectorT10Array;
71  //typedef SquareMatrix<T,5> SqMatrix5;
72 
73  typedef std::map<int,KineticBC<T>*> KineticBCMap;
74  typedef std::map<int,KineticVC<T>*> KineticVCMap;
79  //MacroFields& macroFields;
80 
81  KineticModel(const MeshList& meshes, const GeomFields& geomFields, MacroFields& macroFields, const Quadrature<T>& quad):
82 
83  Model(meshes),
84  _geomFields(geomFields),
85  _quadrature(quad),
86  _macroFields(macroFields),
87  _dsfPtr(_meshes,_quadrature,"dsf_"),
88  _dsfPtr1(_meshes,_quadrature,"dsf1_"),
89  _dsfPtr2(_meshes,_quadrature,"dsf2_"),
90  _dsfEqPtr(_meshes,_quadrature,"dsfEq_"),
91  _dsfEqPtrES(_meshes,_quadrature,"dsfEqES_"),
92  _initialKmodelNorm(),
93  _niters(0)
94  {
95 
96  const int numMeshes = _meshes.size();
97  for (int n=0; n<numMeshes; n++)
98  {
99  const Mesh& mesh = *_meshes[n];
100 
101  KineticVC<T> *vc(new KineticVC<T>());
102  vc->vcType = "flow";
103  _vcMap[mesh.getID()] = vc;
104  }
105  init();
106  SetBoundaryConditions();
107 
108  //weightedMaxwellian(1.0,0.00,0.00,1.0,1.0);
109  InitializeMacroparameters();
110  initializeMaxwellian();
111 
112  ComputeMacroparameters(); //calculate density,velocity,temperature
113 
114  ComputeCollisionfrequency(); //calculate viscosity, collisionFrequency
115  initializeMaxwellianEq(); //equilibrium distribution
116 
117  //EquilibriumDistributionBGK();
118  //callBoundaryConditions();
119 
120  }
121 
122  void init()
123 
124  {
125  const int numMeshes = _meshes.size();
126  for (int n=0; n<numMeshes; n++)
127  {
128  const Mesh& mesh = *_meshes[n];
129 
130  const KineticVC<T>& vc = *_vcMap[mesh.getID()];
131 
132  const StorageSite& cells = mesh.getCells();
133 
134  const int nCells = cells.getCountLevel1();
135  shared_ptr<VectorT3Array> vCell(new VectorT3Array(nCells));
136 
137  VectorT3 initialVelocity;
138  initialVelocity[0] = _options["initialXVelocity"];
139  initialVelocity[1] = _options["initialYVelocity"];
140  initialVelocity[2] = _options["initialZVelocity"];
141  *vCell = initialVelocity;
142  _macroFields.velocity.addArray(cells,vCell);
143 
144 
145  shared_ptr<TArray> pCell(new TArray(nCells));
146  *pCell = _options["operatingPressure"];
147  _macroFields.pressure.addArray(cells,pCell);
148 
149 
150  shared_ptr<TArray> rhoCell(new TArray(nCells));
151  *rhoCell = vc["density"];
152  _macroFields.density.addArray(cells,rhoCell);
153 
154  shared_ptr<TArray> muCell(new TArray(nCells));
155  *muCell = vc["viscosity"];
156  _macroFields.viscosity.addArray(cells,muCell);
157 
158  shared_ptr<TArray> tempCell(new TArray(cells.getCountLevel1()));
159  *tempCell = _options["operatingTemperature"];
160  _macroFields.temperature.addArray(cells,tempCell);
161 
162  shared_ptr<TArray> collFreqCell(new TArray(cells.getCountLevel1()));
163  *collFreqCell = vc["viscosity"];
164  _macroFields.collisionFrequency.addArray(cells,collFreqCell);
165 
166 
167  //coeffs for perturbed BGK distribution function
168  shared_ptr<VectorT5Array> coeffCell(new VectorT5Array(cells.getCountLevel1()));
169  VectorT5 initialCoeff;
170  initialCoeff[0] = 1.0;
171  initialCoeff[1] = 1.0;
172  initialCoeff[2] = 0.0;
173  initialCoeff[3] = 0.0;
174  initialCoeff[4] = 0.0;
175  *coeffCell = initialCoeff;
176  _macroFields.coeff.addArray(cells,coeffCell);
177 
178  //coeffs for perturbed BGK distribution function
179  shared_ptr<VectorT10Array> coeffgCell(new VectorT10Array(cells.getCountLevel1()));
180  VectorT10 initialCoeffg;
181  initialCoeffg[0] = 1.0;
182  initialCoeffg[1] = 1.0;
183  initialCoeffg[2] = 0.0;
184  initialCoeffg[3] = 1.0;
185  initialCoeffg[4] = 0.0;
186  initialCoeffg[5] = 1.0;
187  initialCoeffg[6] = 0.0;
188  initialCoeffg[7] = 0.0;
189  initialCoeffg[8] = 0.0;
190  initialCoeffg[9] = 0.0;
191  *coeffgCell = initialCoeffg;
192  _macroFields.coeffg.addArray(cells,coeffgCell);
193 
194  // used for ESBGK equilibrium distribution function
195  shared_ptr<TArray> tempxxCell(new TArray(cells.getCountLevel1()));
196  *tempxxCell = _options["operatingTemperature"]/3;
197  _macroFields.Txx.addArray(cells,tempxxCell);
198 
199  shared_ptr<TArray> tempyyCell(new TArray(cells.getCountLevel1()));
200  *tempyyCell = _options["operatingTemperature"]/3;
201  _macroFields.Tyy.addArray(cells,tempyyCell);
202 
203  shared_ptr<TArray> tempzzCell(new TArray(cells.getCountLevel1()));
204  *tempzzCell = _options["operatingTemperature"]/3;
205  _macroFields.Tzz.addArray(cells,tempzzCell);
206 
207  shared_ptr<TArray> tempxyCell(new TArray(cells.getCountLevel1()));
208  *tempxyCell = 0.0;
209  _macroFields.Txy.addArray(cells,tempxyCell);
210 
211  shared_ptr<TArray> tempyzCell(new TArray(cells.getCountLevel1()));
212  *tempyzCell = 0.0;
213  _macroFields.Tyz.addArray(cells,tempyzCell);
214 
215  shared_ptr<TArray> tempzxCell(new TArray(cells.getCountLevel1()));
216  *tempzxCell = 0.0;
217  _macroFields.Tzx.addArray(cells,tempzxCell);
218 
219  //Entropy and Entropy Generation Rate for switching
220  shared_ptr<TArray> EntropyCell(new TArray(cells.getCountLevel1()));
221  *EntropyCell = 0.0;
222  _macroFields.Entropy.addArray(cells,EntropyCell);
223 
224  shared_ptr<TArray> EntropyGenRateCell(new TArray(cells.getCountLevel1()));
225  *EntropyGenRateCell = 0.0;
226  _macroFields.EntropyGenRate.addArray(cells,EntropyGenRateCell);
227 
228  shared_ptr<TArray> EntropyGenRateColl(new TArray(cells.getCountLevel1()));
229  *EntropyGenRateColl = 0.0;
230  _macroFields.EntropyGenRate_Collisional.addArray(cells,EntropyGenRateColl);
231 
232 
233  //Pxx,Pyy,Pzz,Pxy,Pxz,Pyz
234  shared_ptr<VectorT6Array> stressCell(new VectorT6Array(nCells));
235  VectorT6 initialstress;
236  initialstress[0] = 1.0;
237  initialstress[1] = 1.0;
238  initialstress[2] = 1.0;
239  initialstress[3] = 0.0;
240  initialstress[4] = 0.0;
241  initialstress[5] = 0.0;
242  *stressCell = initialstress;
243  _macroFields.Stress.addArray(cells,stressCell);
244 
245  //Knq=M300+M120+M102 for Couette with uy
246  shared_ptr<TArray> KnqCell(new TArray(cells.getCountLevel1()));
247  *KnqCell = 0.0;
248  _macroFields.Knq.addArray(cells,KnqCell);
249 
250 
251  //higher order moments of distribution function
252  /*
253  shared_ptr<VectorT3Array> M300Cell(new VectorT3Array(cells.getCount()));
254  VectorT3 initialM300;
255  initialM300[0] = 0.0;
256  initialM300[1] = 0.0;
257  initialM300[2] = 0.0;
258  *M300Cell = initialM300;
259  _macroFields.M300.addArray(cells,M300Cell);
260 
261  shared_ptr<VectorT3Array> M030Cell(new VectorT3Array(cells.getCount()));
262  VectorT3 initialM030;
263  initialM030[0] = 0.0;
264  initialM030[1] = 0.0;
265  initialM030[2] = 0.0;
266  *M300Cell = initialM030;
267  _macroFields.M030.addArray(cells,M030Cell);
268  */
269 
270  const int numDirections = _quadrature.getDirCount();
271  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
272  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
273  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
274  //FILE * pFile;
275  //pFile=fopen("ref_incMEMOSA.txt","w");
276  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups()){
277  const FaceGroup& fg = *fgPtr;
278 
279  if((fg.groupType == "symmetry")||(fg.groupType == "realwall")||(fg.groupType == "velocity-inlet")){
280 
281  const StorageSite& faces = fg.site;
282 
283  const Field& areaMagField = _geomFields.areaMag;
284  const TArray& faceAreaMag = dynamic_cast<const TArray &>(areaMagField[faces]);
285  const Field& areaField = _geomFields.area;
286  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
287 
288  const VectorT3 en = faceArea[0]/faceAreaMag[0];
289  vector<int> tempVec(numDirections);
290 
291  for (int j=0; j<numDirections; j++){
292  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
293  const T cx_incident = cx[j] - 2.0*c_dot_en*en[0];
294  const T cy_incident = cy[j] - 2.0*c_dot_en*en[1];
295  const T cz_incident = cz[j] - 2.0*c_dot_en*en[2];
296  int direction_incident=0;
297  T Rdotprod=1e54;
298  T dotprod=0.0;
299  for (int js=0; js<numDirections; js++){
300  dotprod=pow(cx_incident-cx[js],2)+pow(cy_incident-cy[js],2)+pow(cz_incident-cz[js],2);
301  if (dotprod< Rdotprod){
302  Rdotprod =dotprod;
303  direction_incident=js;}
304  }
305  tempVec[j] = direction_incident;
306  //fprintf(pFile,"%d %d %d \n",fg.id, j,direction_incident);
307 
308  }
309  const int fgid=fg.id;
310  _faceReflectionArrayMap[fgid] = tempVec; //add to map
311 
312  }
313  }
314  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups()){
315  const FaceGroup& fg = *fgPtr;
316  if(fg.groupType == "NSinterface"){
317  const StorageSite& Intfaces = fg.site;
318 
319  shared_ptr<VectorT3Array> InterfaceVelFace(new VectorT3Array(Intfaces.getCount()));
320  InterfaceVelFace ->zero();
321  _macroFields.InterfaceVelocity.addArray(Intfaces,InterfaceVelFace);
322 
323  shared_ptr<StressTensorArray> InterfaceStressFace(new StressTensorArray(Intfaces.getCount()));
324  InterfaceStressFace ->zero();
325  _macroFields.InterfaceStress.addArray(Intfaces,InterfaceStressFace);
326 
327  shared_ptr<TArray> InterfacePressFace(new TArray(Intfaces.getCount()));
328  *InterfacePressFace = _options["operatingPressure"];
329  _macroFields.InterfacePressure.addArray(Intfaces,InterfacePressFace);
330 
331  shared_ptr<TArray> InterfaceDensityFace(new TArray(Intfaces.getCount()));
332  *InterfaceDensityFace =vc["density"];
333  _macroFields.InterfaceDensity.addArray(Intfaces,InterfaceDensityFace);
334 
335 
336  }
337 
338  //fclose(pFile);
339 
340 
341  } //end of loop through meshes
342  _niters =0;
343  _initialKmodelNorm = MFRPtr();
344  //_initialKmodelvNorm = MFRPtr();
345 
346  }
347  }
348 
349  void InitializeMacroparameters()
350  { const int numMeshes =_meshes.size();
351  for (int n=0; n<numMeshes; n++)
352  {
353  const Mesh& mesh = *_meshes[n];
354  const StorageSite& cells = mesh.getCells();
355  const int nCells = cells.getCountLevel1();
356  // const double pi(3.14159);
357  const double pi=_options.pi;
358  TArray& Entropy = dynamic_cast<TArray&>(_macroFields.Entropy[cells]);
359  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
360  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
361 
362  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
363  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[cells]);
364 
365  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
366  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
367  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
368  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
369  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
370  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
371  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
372  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
373  //if ( MPI::COMM_WORLD.Get_rank() == 0 ) {cout << "ncells="<<nCells<<endl;}
374 
375  TArray& Knq = dynamic_cast<TArray&>(_macroFields.Knq[cells]);
376  //initialize density,velocity
377  for(int c=0; c<nCells;c++)
378  {
379  density[c] =1.0;
380  v[c][0]=0.0;
381  v[c][1]=0.0;
382  v[c][2]=0.0;
383  temperature[c]=1.0;
384  pressure[c]=temperature[c]*density[c];
385 
386  //BGK
387  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
388  coeff[c][1]=1/temperature[c];
389  coeff[c][2]=0.0;coeff[c][3]=0.0;coeff[c][4]=0.0;
390 
391  Entropy[c]=0.0;
392  if(_options.fgamma ==2){
393  //ESBGK
394  coeffg[c][0]=coeff[c][0];
395  coeffg[c][1]=coeff[c][1];
396  coeffg[c][2]=coeff[c][2];
397  coeffg[c][3]=coeff[c][1];
398  coeffg[c][4]=coeff[c][3];
399  coeffg[c][5]=coeff[c][1];
400  coeffg[c][6]=coeff[c][4];
401  coeffg[c][7]=0.0;
402  coeffg[c][8]=0.0;
403  coeffg[c][9]=0.0;
404  }
405 
406  Txx[c]=0.5;
407  Tyy[c]=0.5;
408  Tzz[c]=0.5;
409  Txy[c]=0.0;
410  Tyz[c]=0.0;
411  Tzx[c]=0.0;
412 
413  Knq[c]=0.0;
414  }
415  }
416  }
417  void InitializeFgammaCoefficients()
418  {
419  const int numMeshes =_meshes.size();
420  for (int n=0; n<numMeshes; n++)
421  {
422  const Mesh& mesh = *_meshes[n];
423  const StorageSite& cells = mesh.getCells();
424  const int nCells = cells.getCountLevel1();
425  // const double pi(3.14159);
426  const double pi=_options.pi;
427 
428  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
429  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
430  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
431  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
432  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
433  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
434  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
435  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
436  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
437  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
438 
439  for(int c=0; c<nCells;c++)
440  {
441 
442  //BGK
443  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
444  coeff[c][1]=1/temperature[c];
445  coeff[c][2]=0.0;coeff[c][3]=0.0;coeff[c][4]=0.0;
446 
447  if(_options.fgamma ==2){
448  //ESBGK
449  coeffg[c][0]=coeff[c][0];
450  coeffg[c][1]=coeff[c][1];
451  coeffg[c][2]=coeff[c][2];
452  coeffg[c][3]=coeff[c][1];
453  coeffg[c][4]=coeff[c][3];
454  coeffg[c][5]=coeff[c][1];
455  coeffg[c][6]=coeff[c][4];
456  coeffg[c][7]=0.0;
457  coeffg[c][8]=0.0;
458  coeffg[c][9]=0.0;
459 
460  Txx[c]=0.5*temperature[c];
461  Tyy[c]=0.5*temperature[c];
462  Tzz[c]=0.5*temperature[c];
463  Txy[c]=0.0;
464  Tyz[c]=0.0;
465  Tzx[c]=0.0;
466  }
467  }
468  }
469  }
470 
471  void ComputeMacroparameters()
472  {
473  //FILE * pFile;
474  //pFile = fopen("distfun_mf.txt","w");
475  const int numMeshes = _meshes.size();
476  for (int n=0; n<numMeshes; n++)
477  {
478 
479  const Mesh& mesh = *_meshes[n];
480  const StorageSite& cells = mesh.getCells();
481  const int nCells = cells.getCountLevel1(); //
482 
483 
484  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
485  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
486  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
487  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[cells]);
488  const int N123 = _quadrature.getDirCount();
489 
490  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
491  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
492  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
493  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
494 
495  VectorT6Array& stress = dynamic_cast<VectorT6Array&>(_macroFields.Stress[cells]);
496 
497  //initialize density,velocity,temperature to zero
498  for(int c=0; c<nCells;c++)
499  {
500  density[c]=0.0;
501  v[c][0]=0.0;
502  v[c][1]=0.0;
503  v[c][2]=0.0;
504  temperature[c]=0.0;
505  stress[c][0]=0.0;stress[c][1]=0.0;stress[c][2]=0.0;
506  stress[c][3]=0.0;stress[c][4]=0.0;stress[c][5]=0.0;
507 
508  }
509 
510 
511  for(int j=0;j<N123;j++){
512 
513  Field& fnd = *_dsfPtr.dsf[j];
514  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
515 
516  //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]);
517 
518  for(int c=0; c<nCells;c++){
519  density[c] = density[c]+wts[j]*f[c];
520  v[c][0]= v[c][0]+(cx[j]*f[c])*wts[j];
521  v[c][1]= v[c][1]+(cy[j]*f[c])*wts[j];
522  v[c][2]= v[c][2]+(cz[j]*f[c])*wts[j];
523  temperature[c]= temperature[c]+(pow(cx[j],2.0)+pow(cy[j],2.0)
524  +pow(cz[j],2.0))*f[c]*wts[j];
525 
526  }
527 
528  }
529 
530 
531 
532  for(int c=0; c<nCells;c++){
533  v[c][0]=v[c][0]/density[c];
534  v[c][1]=v[c][1]/density[c];
535  v[c][2]=v[c][2]/density[c];
536  temperature[c]=temperature[c]-(pow(v[c][0],2.0)
537  +pow(v[c][1],2.0)
538  +pow(v[c][2],2.0))*density[c];
539  temperature[c]=temperature[c]/(1.5*density[c]);
540  pressure[c]=density[c]*temperature[c];
541  }
542 
543  //Find Pxx,Pyy,Pzz,Pxy,Pyz,Pzx, etc in field
544 
545  for(int j=0;j<N123;j++){
546  Field& fnd = *_dsfPtr.dsf[j];
547  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
548  for(int c=0; c<nCells;c++){
549  stress[c][0] +=pow((cx[j]-v[c][0]),2.0)*f[c]*wts[j];
550  stress[c][1] +=pow((cy[j]-v[c][1]),2.0)*f[c]*wts[j];
551  stress[c][2] +=pow((cz[j]-v[c][2]),2.0)*f[c]*wts[j];
552  stress[c][3] +=(cx[j]-v[c][0])*(cy[j]-v[c][1])*f[c]*wts[j];
553  stress[c][4] +=(cy[j]-v[c][1])*(cz[j]-v[c][2])*f[c]*wts[j];
554  stress[c][5] +=(cz[j]-v[c][2])*(cx[j]-v[c][0])*f[c]*wts[j];
555 
556  }}
557 
558 
559  }// end of loop over nmeshes
560  //fclose(pFile);
561  }
562 
563 
564 
565  void ComputeMacroparametersESBGK()
566  {
567 
568  const int numMeshes = _meshes.size();
569  for (int n=0; n<numMeshes; n++)
570  {
571  const Mesh& mesh = *_meshes[n];
572  const StorageSite& cells = mesh.getCells();
573  const int nCells = cells.getCountLevel1();
574 
575  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
576  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
577 
578  TArray& Txx = dynamic_cast<TArray&>(_macroFields.Txx[cells]);
579  TArray& Tyy = dynamic_cast<TArray&>(_macroFields.Tyy[cells]);
580  TArray& Tzz = dynamic_cast<TArray&>(_macroFields.Tzz[cells]);
581  TArray& Txy = dynamic_cast<TArray&>(_macroFields.Txy[cells]);
582  TArray& Tyz = dynamic_cast<TArray&>(_macroFields.Tyz[cells]);
583  TArray& Tzx = dynamic_cast<TArray&>(_macroFields.Tzx[cells]);
584 
585  const int N123 = _quadrature.getDirCount();
586 
587  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
588  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
589  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
590  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
591 
592  const double Pr=_options.Prandtl;
593  //cout <<"Prandlt" <<Pr<<endl;
594 
595  //initialize density,velocity,temperature to zero
596  for(int c=0; c<nCells;c++)
597  {
598  Txx[c]=0.0;
599  Tyy[c]=0.0;
600  Tzz[c]=0.0;
601  Txy[c]=0.0;
602  Tyz[c]=0.0;
603  Tzx[c]=0.0;
604  }
605  for(int j=0;j<N123;j++){
606 
607  Field& fnd = *_dsfPtr.dsf[j];
608  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
609  Field& fndEq = *_dsfEqPtr.dsf[j];
610  const TArray& fgam = dynamic_cast<const TArray&>(fndEq[cells]);
611  for(int c=0; c<nCells;c++){
612  Txx[c]=Txx[c]+pow(cx[j]-v[c][0],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
613  Tyy[c]=Tyy[c]+pow(cy[j]-v[c][1],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j] ;
614  Tzz[c]=Tzz[c]+pow(cz[j]-v[c][2],2)*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
615  //Txy[c]=Txy[c]+(cx[j])*(cy[j])*((1-1/Pr)*f[c]+1/Pr*fgam[c])*wts[j];
616  //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];
617 
618 
619  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];
620  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];
621  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];
622  }
623  }
624 
625  for(int c=0; c<nCells;c++){
626  Txx[c]=Txx[c]/density[c];
627  Tyy[c]=Tyy[c]/density[c];
628  Tzz[c]=Tzz[c]/density[c];
629  Txy[c]=Txy[c]/density[c];
630  Tyz[c]=Tyz[c]/density[c];
631  Tzx[c]=Tzx[c]/density[c];
632  }
633 
634  }
635  }
636 
637 
638 
639  /*
640  * Collision frequency
641  *
642  *
643  */
644 
645  void ComputeCollisionfrequency() {
646  const int numMeshes = _meshes.size();
647  for (int n=0; n<numMeshes; n++)
648  {
649 
650  const T rho_init=_options["rho_init"];
651  const T T_init= _options["T_init"];
652  const T mu_w= _options["mu_w"];
653  const T Tmuref= _options["Tmuref"];
654  const T muref= _options["muref"];
655  const T R=8314.0/_options["molecularWeight"];
656  const T nondim_length=_options["nonDimLt"];
657 
658  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
659 
660 
661  const Mesh& mesh = *_meshes[n];
662  const StorageSite& cells = mesh.getCells();
663  const int nCells = cells.getCountLevel1();
664 
665  TArray& density = dynamic_cast<TArray&>(_macroFields.density[cells]);
666  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[cells]);
667  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[cells]);
668 
669  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
670  for(int c=0; c<nCells;c++)
671  {
672  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
673  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
674  }
675 
676  if(_options.fgamma==2){
677  for(int c=0; c<nCells;c++)
678  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
679  }
680 
681  }
682  }
683 
684  void MomentHierarchy() {
685  const int numMeshes = _meshes.size();
686  for (int n=0; n<numMeshes; n++)
687  {
688  const int Knq_dir=_options.Knq_direction;
689  const Mesh& mesh = *_meshes[n];
690  const StorageSite& cells = mesh.getCells();
691  const int nCells = cells.getCountLevel1();
692  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
693  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
694  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
695  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
696  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
697  TArray& Knq = dynamic_cast<TArray&>(_macroFields.Knq[cells]);
698  const int num_directions = _quadrature.getDirCount();
699  if (Knq_dir ==0){
700  for(int j=0;j<num_directions;j++){
701  Field& fnd = *_dsfPtr.dsf[j];
702  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
703  for(int c=0; c<nCells;c++){
704  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));
705  }
706  }}
707  else if(Knq_dir ==1){
708  for(int j=0;j<num_directions;j++){
709  Field& fnd = *_dsfPtr.dsf[j];
710  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
711  for(int c=0; c<nCells;c++){
712  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));
713  }
714  }}
715 
716  else if(Knq_dir ==2){
717  for(int j=0;j<num_directions;j++){
718  Field& fnd = *_dsfPtr.dsf[j];
719  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
720  for(int c=0; c<nCells;c++){
721  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));
722  }
723  }}
724  }
725  }
726  void EntropyGeneration() {
727  const int numMeshes = _meshes.size();
728  for (int n=0; n<numMeshes; n++)
729  {
730  const T rho_init=_options["rho_init"];
731  const T T_init= _options["T_init"];
732  const T molwt=_options["molecularWeight"]*1E-26/6.023;
733  const T R=8314.0/_options["molecularWeight"];
734  const T u_init=pow(2.0*R*T_init,0.5);
735  const T Planck=_options.Planck;
736  const T h3bm4u3=pow(Planck,3)/ pow(molwt,4)*rho_init/pow(u_init,3);
737  //cout << "h3bm4u3 " << h3bm4u3 <<endl;
738  //cout <<" u_init "<<u_init<<" rho_init "<<rho_init<<endl;
739  const Mesh& mesh = *_meshes[n];
740  const StorageSite& cells = mesh.getCells();
741  const int nCells = cells.getCountLevel1();
742  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
743  TArray& Entropy = dynamic_cast<TArray&>(_macroFields.Entropy[cells]);
744  TArray& EntropyGenRate_Collisional = dynamic_cast<TArray&>(_macroFields.EntropyGenRate_Collisional[cells]);
745  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
746  for(int c=0; c<nCells;c++){
747  Entropy[c]=0.0;EntropyGenRate_Collisional[c]=0.0;
748  }
749  const int num_directions = _quadrature.getDirCount();
750  if (_options.fgamma ==2){
751  for(int j=0;j<num_directions;j++){
752  Field& fnd = *_dsfPtr.dsf[j];
753  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
754  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
755  const TArray& fgam = dynamic_cast<const TArray&>(feqES[cells]);
756  for(int c=0; c<nCells;c++){
757  Entropy[c]=Entropy[c]+f[c]*wts[j]*(1-log(h3bm4u3*f[c]));
758  EntropyGenRate_Collisional[c]+= (f[c]-fgam[c])*collisionFrequency[c]*log(h3bm4u3*f[c])*wts[j];
759  }
760  }
761  }
762  else{
763  for(int j=0;j<num_directions;j++){
764  Field& fnd = *_dsfPtr.dsf[j];
765  Field& feq = *_dsfEqPtr.dsf[j];
766  const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
767  const TArray& fgam = dynamic_cast<const TArray&>(feq[cells]);
768  for(int c=0; c<nCells;c++){
769  Entropy[c]=Entropy[c]+f[c]*wts[j]*(1-log(h3bm4u3*f[c]));
770  EntropyGenRate_Collisional[c]+=(f[c]-fgam[c])*collisionFrequency[c]*(log(h3bm4u3*f[c]))*wts[j];
771  }
772  }
773  }
774 
775 
776  }
777  }
778 
779  void initializeMaxwellianEq()
780  {
781  const int numMeshes = _meshes.size();
782  for (int n=0; n<numMeshes; n++)
783  {
784  const Mesh& mesh = *_meshes[n];
785  const StorageSite& cells = mesh.getCells();
786  const int nCells = cells.getCountLevel1();
787  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
788  const VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
789  const VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
790 
791  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
792  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
793  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
794  const int numFields= _quadrature.getDirCount();
795 
796  for(int j=0;j< numFields;j++){
797  Field& fndEq = *_dsfEqPtr.dsf[j];
798  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
799  for(int c=0; c<nCells;c++){
800  fEq[c]=coeff[c][0]*exp(-coeff[c][1]*(pow(cx[j]-v[c][0],2)+pow(cy[j]-v[c][1],2)
801  +pow(cz[j]-v[c][2],2))+coeff[c][2]*(cx[j]-v[c][0])
802  +coeff[c][3]*(cy[j]-v[c][1])+coeff[c][4]*(cz[j]-v[c][2]));
803  }
804  }
805 
806  if(_options.fgamma==2){
807 
808  for(int j=0;j< numFields;j++){
809  Field& fndEqES = *_dsfEqPtrES.dsf[j];
810  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
811  for(int c=0; c<nCells;c++){
812  T Cc1=(cx[j]-v[c][0]);
813  T Cc2=(cy[j]-v[c][1]);
814  T Cc3=(cz[j]-v[c][2]);
815  fEqES[c]=coeffg[c][0]*exp(-coeffg[c][1]*pow(Cc1,2)+coeffg[c][2]*Cc1
816  -coeffg[c][3]*pow(Cc2,2)+coeffg[c][4]*Cc2
817  -coeffg[c][5]*pow(Cc3,2)+coeffg[c][6]*Cc3
818  +coeffg[c][7]*cx[j]*cy[j]+coeffg[c][8]*cy[j]*cz[j]
819  +coeffg[c][9]*cz[j]*cx[j]);
820  }
821  }
822  }
823 
824 
825 
826  }
827  }
828 
829  void NewtonsMethodBGK(const int ktrial)
830  {
831  const int numMeshes = _meshes.size();
832  for (int n=0; n<numMeshes; n++)
833  {
834  //cout << " NewtonsMethod" <<endl;
835  const T tolx=_options["ToleranceX"];
836  const T tolf=_options["ToleranceF"];
837  const int sizeC=5;
838  const Mesh& mesh = *_meshes[n];
839  const StorageSite& cells = mesh.getCells();
840  const int nCells = cells.getCountLevel1();
841 
842  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
843  const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
844  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
845 
846  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
847 
848  for(int c=0; c<nCells;c++){
849 
850  for (int trial=0;trial<ktrial;trial ++){
851  SquareMatrix<T,sizeC> fjac(0);
852  SquareMatrix<T,sizeC> fjacinv(0);
853  VectorT5 fvec;
854 
855 
856  fvec[0]=density[c];
857  fvec[1]=density[c]*v[c][0];
858  fvec[2]=density[c]*v[c][1];
859  fvec[3]=density[c]*v[c][2];
860  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));
861 
862 
863  setJacobianBGK(fjac,fvec,coeff[c],v[c],c);
864 
865  //solve using GE or inverse
866  T errf=0.;
867  for (int row=0;row<sizeC;row++){errf+=fabs(fvec[row]);}
868 
869  if(errf <= tolf)
870  break;
871  VectorT5 pvec;
872  for (int row=0;row<sizeC;row++){pvec[row]=-fvec[row];}//rhs
873 
874 
875  //solve Ax=b for x
876  //p=GE_elim(fjac,p,3);
877  VectorT5 xvec;
878  fjacinv=inverseGauss(fjac,sizeC);
879 
880 
881 
882  for (int row=0;row<sizeC;row++){
883  xvec[row]=0.0;
884  for (int col=0;col<sizeC;col++){
885  xvec[row]+=fjacinv(row,col)*pvec[col];}
886  }
887  //check for convergence, update
888 
889  T errx=0.;
890  for (int row=0;row<sizeC;row++){
891  errx +=fabs(xvec[row]);
892  coeff[c][row]+= xvec[row];
893  }
894 
895 
896  if(errx <= tolx)
897  break;
898 
899  }
900 
901 
902  }
903  }
904 
905  }
906 
907  void EquilibriumDistributionBGK()
908  {
909  const int ktrial=_options.NewtonsMethod_ktrial;
910  const int numMeshes = _meshes.size();
911  for (int n=0; n<numMeshes; n++)
912  {
913  const Mesh& mesh = *_meshes[n];
914  const StorageSite& cells = mesh.getCells();
915  const int nCells = cells.getCountLevel1();
916 
917  //const double pi=_options.pi;
918  //const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
919  //const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
920 
921  //initialize coeff
922  VectorT5Array& coeff = dynamic_cast<VectorT5Array&>(_macroFields.coeff[cells]);
923 
924  /*
925  for(int c=0; c<nCells;c++){
926  coeff[c][0]=density[c]/pow((pi*temperature[c]),1.5);
927  coeff[c][1]=1/temperature[c];
928  coeff[c][2]=0.0;
929  coeff[c][3]=0.0;
930  coeff[c][4]=0.0;
931  }
932  */
933  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
934  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
935  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
936  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
937  const int numFields= _quadrature.getDirCount();
938 
939  //call Newtons Method
940 
941  NewtonsMethodBGK(ktrial);
942 
943  //calculate perturbed maxwellian for BGK
944  for(int j=0;j< numFields;j++){
945  Field& fndEq = *_dsfEqPtr.dsf[j];
946  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
947  for(int c=0; c<nCells;c++){
948  fEq[c]=coeff[c][0]*exp(-coeff[c][1]*(pow(cx[j]-v[c][0],2)+pow(cy[j]-v[c][1],2)
949  +pow(cz[j]-v[c][2],2))+coeff[c][2]*(cx[j]-v[c][0])
950  +coeff[c][3]*(cy[j]-v[c][1])+coeff[c][4]*(cz[j]-v[c][2]));
951  }
952 
953  }
954  }
955  }
956 
957 
958  void setJacobianBGK(SquareMatrix<T,5>& fjac, VectorT5& fvec, const VectorT5& xn,const VectorT3& v,const int c)
959  {
960 
961 
962  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
963  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
964  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
965  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
966 
967  const TArray2D& malphaBGK = dynamic_cast<const TArray2D&>(*_quadrature.malphaBGKPtr);
968  const int numFields= _quadrature.getDirCount();
969  VectorT5 mexp;
970 
971  for(int j=0;j< numFields;j++){
972  T Cconst=pow(cx[j]-v[0],2.0)+pow(cy[j]-v[1],2.0)+pow(cz[j]-v[2],2.0);
973  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];
974 
975  for (int row=0;row<5;row++){
976  fvec[row]+= -Econst*malphaBGK(j,row); //smm
977 
978  //fvec[row]=tvec[row]+fvec[row]; //mma
979  }
980 
981  mexp[0]=-Econst/xn[0];
982  mexp[1]=Econst*Cconst;
983  mexp[2]=-Econst*(cx[j]-v[0]);
984  mexp[3]=-Econst*(cy[j]-v[1]);
985  mexp[4]=-Econst*(cz[j]-v[2]);
986  for (int row=0;row<5;row++){
987  for (int col=0;col<5;col++){
988  fjac(row,col)+=malphaBGK(j,row)*mexp[col]; //new
989  }
990  }
991 
992 
993  }
994 
995 
996 
997  }
998 
999  void NewtonsMethodESBGK(const int ktrial)
1000  {
1001  // cout<< "Inside Newtons Method" <<endl;
1002  const int numMeshes = _meshes.size();
1003  for (int n=0; n<numMeshes; n++)
1004  {
1005  const T tolx=_options["ToleranceX"];
1006  const T tolf=_options["ToleranceF"];
1007  const int sizeC=10;
1008  const Mesh& mesh = *_meshes[n];
1009  const StorageSite& cells = mesh.getCells();
1010  const int nCells = cells.getCountLevel1();
1011 
1012  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1013 
1014  const TArray& Txx = dynamic_cast<const TArray&>(_macroFields.Txx[cells]);
1015  const TArray& Tyy = dynamic_cast<const TArray&>(_macroFields.Tyy[cells]);
1016  const TArray& Tzz = dynamic_cast<const TArray&>(_macroFields.Tzz[cells]);
1017  const TArray& Txy = dynamic_cast<const TArray&>(_macroFields.Txy[cells]);
1018  const TArray& Tyz = dynamic_cast<const TArray&>(_macroFields.Tyz[cells]);
1019  const TArray& Tzx = dynamic_cast<const TArray&>(_macroFields.Tzx[cells]);
1020 
1021  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1022 
1023  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1024 
1025  for(int c=0; c<nCells;c++){
1026  // {cout <<"NM:ES" <<c <<endl;}
1027  for (int trial=0;trial<ktrial;trial ++){
1028  SquareMatrix<T,sizeC> fjac(0);
1029  SquareMatrix<T,sizeC> fjacinv(0);
1030  Vector<T,sizeC> fvec;
1031  // if (c==_options.printCellNumber){cout <<"trial" <<trial <<endl;}
1032 
1033  fvec[0]=density[c];
1034  fvec[1]=density[c]*v[c][0];
1035  fvec[2]=density[c]*v[c][1];
1036  fvec[3]=density[c]*v[c][2];
1037  fvec[4]=density[c]*(pow(v[c][0],2)+Txx[c]);
1038  fvec[5]=density[c]*(pow(v[c][1],2)+Tyy[c]);
1039  fvec[6]=density[c]*(pow(v[c][2],2)+Tzz[c]);
1040  fvec[7]=density[c]*(v[c][0]*v[c][1]+Txy[c]);
1041  fvec[8]=density[c]*(v[c][1]*v[c][2]+Tyz[c]);
1042  fvec[9]=density[c]*(v[c][2]*v[c][0]+Tzx[c]);
1043 
1044  //calculate Jacobian
1045  setJacobianESBGK(fjac,fvec,coeffg[c],v[c],c);
1046 
1047 
1048  //solve using GaussElimination
1049  T errf=0.; //and Jacobian matrix in fjac.
1050  for (int row=0;row<sizeC;row++){errf+=fabs(fvec[row]);}
1051  if(errf <= tolf)
1052  break;
1053  Vector<T,sizeC> pvec;
1054  for (int row=0;row<sizeC;row++){pvec[row]=-fvec[row];}
1055 
1056 
1057  //solve Ax=b for x
1058  //p=GE_elim(fjac,p,3);
1059  Vector<T,sizeC> xvec;
1060  fjacinv=inverseGauss(fjac,sizeC);
1061  for (int row=0;row<sizeC;row++){
1062  xvec[row]=0.0;
1063  for (int col=0;col<sizeC;col++){
1064  xvec[row]+=fjacinv(row,col)*pvec[col];
1065 
1066  }
1067  }
1068  /*
1069  if (c==_options.printCellNumber){
1070  cout << " cg0 "<<coeffg[c][0]<<" cg1 "<<coeffg[c][1]<<" cg2 "<<coeffg[c][2] << endl;
1071  cout <<" cg3 " <<coeffg[c][3]<< " cg4 "<<coeffg[c][4]<<" cg5 "<<coeffg[c][5]<<" cg6 "<<coeffg[c][6] << endl;
1072  cout <<" cg7 " <<coeffg[c][7]<< " cg8 "<<coeffg[c][8]<<" cg9 "<<coeffg[c][9]<<endl;
1073 
1074 
1075  //cout << " fvec-ESBGK " << fvec[4] <<fvec[5]<<fvec[6]<<fvec[7]<<fvec[8]<<fvec[9] <<endl;
1076  FILE * pFile;
1077  pFile = fopen("fvecfjac.dat","wa");
1078  //fprintf(pFile,"%s %d \n","trial",trial);
1079  for (int mat_col=0;mat_col<sizeC;mat_col++){fprintf(pFile,"%12.4E",fvec[mat_col]);}
1080  fprintf(pFile,"\n");
1081  for (int mat_row=0;mat_row<sizeC;mat_row++){
1082  for (int mat_col=0;mat_col<sizeC;mat_col++){
1083  fprintf(pFile,"%12.4E",fjac(mat_row,mat_col));}
1084  fprintf(pFile,"\n");}
1085  // fprintf(pFile,"done \n");
1086  //inverse
1087  for (int mat_row=0;mat_row<sizeC;mat_row++){
1088  for (int mat_col=0;mat_col<sizeC;mat_col++){
1089  fprintf(pFile,"%12.4E",fjacinv(mat_row,mat_col));}
1090  fprintf(pFile,"\n");}
1091  //solution
1092  for (int mat_col=0;mat_col<sizeC;mat_col++){
1093  fprintf(pFile,"%12.4E",pvec[mat_col]);}
1094  }
1095  */
1096 
1097  //check for convergence, update
1098  T errx=0.;//%Check root convergence.
1099  for (int row=0;row<sizeC;row++){
1100  errx +=fabs(xvec[row]);
1101  coeffg[c][row]+= xvec[row];
1102  }
1103  //if (c==_options.printCellNumber){cout <<"errx "<<errx<<endl;}
1104  if(errx <= tolx)
1105  break;
1106 
1107  }
1108 
1109  }
1110  }
1111  }
1112  void EquilibriumDistributionESBGK()
1113  {
1114  ComputeMacroparametersESBGK();
1115  const int ktrial=_options.NewtonsMethod_ktrial;
1116  const int numMeshes = _meshes.size();
1117  for (int n=0; n<numMeshes; n++)
1118  {
1119  const Mesh& mesh = *_meshes[n];
1120  const StorageSite& cells = mesh.getCells();
1121  const int nCells = cells.getCountLevel1();
1122 
1123 
1124  //const VectorT5Array& coeff = dynamic_cast<const VectorT5Array&>(_macroFields.coeff[cells]);
1125  //initialize coeffg
1126  VectorT10Array& coeffg = dynamic_cast<VectorT10Array&>(_macroFields.coeffg[cells]);
1127  /*
1128  for(int c=0; c<nCells;c++){
1129  coeffg[c][0]=coeff[c][0];
1130  coeffg[c][1]=coeff[c][1];
1131  coeffg[c][2]=coeff[c][2];
1132  coeffg[c][3]=coeff[c][1];
1133  coeffg[c][4]=coeff[c][3];
1134  coeffg[c][5]=coeff[c][1];
1135  coeffg[c][6]=coeff[c][4];
1136  coeffg[c][7]=0.0;
1137  coeffg[c][8]=0.0;
1138  coeffg[c][9]=0.0;
1139  }
1140  */
1141  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1142  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1143  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1144  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1145  const int numFields= _quadrature.getDirCount();
1146 
1147  NewtonsMethodESBGK(ktrial);
1148 
1149  for(int j=0;j< numFields;j++){
1150  Field& fndEqES = *_dsfEqPtrES.dsf[j];
1151  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
1152  for(int c=0; c<nCells;c++){
1153  T Cc1=(cx[j]-v[c][0]);
1154  T Cc2=(cy[j]-v[c][1]);
1155  T Cc3=(cz[j]-v[c][2]);
1156  fEqES[c]=coeffg[c][0]*exp(-coeffg[c][1]*pow(Cc1,2)+coeffg[c][2]*Cc1
1157  -coeffg[c][3]*pow(Cc2,2)+coeffg[c][4]*Cc2
1158  -coeffg[c][5]*pow(Cc3,2)+coeffg[c][6]*Cc3
1159  +coeffg[c][7]*cx[j]*cy[j]+coeffg[c][8]*cy[j]*cz[j]
1160  +coeffg[c][9]*cz[j]*cx[j]);
1161  }
1162 
1163  }
1164  }
1165  }
1166 
1167  void setJacobianESBGK(SquareMatrix<T,10>& fjac, VectorT10& fvec, const VectorT10& xn,const VectorT3& v,const int c)
1168  {
1169  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1170  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1171  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1172  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
1173 
1174  const TArray2D& malphaESBGK = dynamic_cast<const TArray2D&>(*_quadrature.malphaESBGKPtr);
1175  const int numFields= _quadrature.getDirCount();
1176  VectorT10 mexp;
1177 
1178  for(int j=0;j< numFields;j++){
1179  T Cc1=cx[j]-v[0];
1180  T Cc2=cy[j]-v[1];
1181  T Cc3=cz[j]-v[2];
1182  T Econst=xn[0]*exp(-xn[1]*pow(Cc1,2)+xn[2]*Cc1-xn[3]*pow(Cc2,2)+ xn[4]*Cc2
1183  -xn[5]*pow(Cc3,2)+xn[6]*Cc3
1184  +xn[7]*cx[j]*cy[j]+xn[8]*cy[j]*cz[j]+xn[9]*cz[j]*cx[j])*wts[j];
1185 
1186  for (int row=0;row<10;row++){
1187  fvec[row]+= -Econst*malphaESBGK(j,row); //smm
1188  }
1189 
1190  mexp[0]=-Econst/xn[0];
1191  mexp[1]=Econst*pow(Cc1,2);
1192  mexp[2]=-Econst*Cc1;
1193  mexp[3]=Econst*pow(Cc2,2);
1194  mexp[4]=-Econst*Cc2;
1195  mexp[5]=Econst*pow(Cc3,2);
1196  mexp[6]=-Econst*Cc3;
1197 
1198  mexp[7]=-Econst*cx[j]*cy[j];
1199  mexp[8]=-Econst*cy[j]*cz[j];
1200  mexp[9]=-Econst*cz[j]*cx[j];
1201 
1202  for (int row=0;row<10;row++){
1203  for (int col=0;col<10;col++){
1204  fjac(row,col)+=malphaESBGK(j,row)*mexp[col]; //new
1205  }
1206  }
1207 
1208 
1209 
1210 
1211  }
1212 
1213 
1214  }
1215 
1216  void initializeMaxwellian()
1217  {
1218  const int numMeshes = _meshes.size();
1219  for (int n=0; n<numMeshes; n++)
1220  {
1221  const Mesh& mesh = *_meshes[n];
1222  const StorageSite& cells = mesh.getCells();
1223  const int nCells = cells.getCountLevel1();
1224 
1225  const double pi=_options.pi;
1226  const TArray& density = dynamic_cast<const TArray&>(_macroFields.density[cells]);
1227  const TArray& temperature = dynamic_cast<const TArray&>(_macroFields.temperature[cells]);
1228  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
1229 
1230  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1231  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1232  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1233  const int numFields= _quadrature.getDirCount();
1234 
1235  for(int j=0;j< numFields;j++){
1236  Field& fnd = *_dsfPtr.dsf[j];
1237  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
1238  for(int c=0; c<nCells;c++){
1239  f[c]=density[c]/pow((pi*temperature[c]),1.5)*
1240  exp(-(pow((cx[j]-v[c][0]),2.0)+pow((cy[j]-v[c][1]),2.0)+
1241  pow((cz[j]-v[c][2]),2.0))/temperature[c]);
1242 
1243 
1244  }
1245 
1246  if (_options.transient)
1247  //updateTime();
1248 
1249  {
1250  Field& fnd1 = *_dsfPtr1.dsf[j];
1251  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
1252  for (int c=0;c<nCells;c++)
1253  f1[c] = f[c];
1254  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
1255  if (_options.timeDiscretizationOrder > 1)
1256  {
1257  Field& fnd2 = *_dsfPtr2.dsf[j];
1258  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
1259  for (int c=0;c<nCells;c++)
1260  f2[c] = f[c];
1261  }
1262  }
1263 
1264 
1265  }
1266  }
1267  }
1268 
1269  void weightedMaxwellian(double weight1,double uvel1,double vvel1,double wvel1,double uvel2,double vvel2,double wvel2,double temp1,double temp2)
1270  {
1271  const int numMeshes = _meshes.size();
1272  for (int n=0; n<numMeshes; n++)
1273  {
1274  const Mesh& mesh = *_meshes[n];
1275  const StorageSite& cells = mesh.getCells();
1276  const int nCells = cells.getCountLevel1();
1277  //double pi(acos(-1.0));
1278  const double pi=_options.pi;
1279  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1280  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1281  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1282  const int numFields= _quadrature.getDirCount();
1283 
1284  for(int j=0;j< numFields;j++){
1285  Field& fnd = *_dsfPtr.dsf[j];
1286  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
1287  for(int c=0; c<nCells;c++){
1288  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)
1289  +(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);
1290  }
1291  if (_options.transient)
1292  {
1293  Field& fnd1 = *_dsfPtr1.dsf[j];
1294  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
1295  for (int c=0;c<nCells;c++)
1296  f1[c] = f[c];
1297  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
1298  if (_options.timeDiscretizationOrder > 1)
1299  {
1300  Field& fnd2 = *_dsfPtr2.dsf[j];
1301  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
1302  for (int c=0;c<nCells;c++)
1303  f2[c] = f[c];
1304  }
1305  }
1306  }
1307  }
1308  }
1309  void weightedMaxwellian(double weight1,double uvel1,double uvel2,double temp1,double temp2)
1310  {
1311  const double vvel1=0.0;
1312  const double wvel1=0.0;
1313  const double vvel2=0.0;
1314  const double wvel2=0.0;
1315  const int numMeshes = _meshes.size();
1316  for (int n=0; n<numMeshes; n++)
1317  {
1318  const Mesh& mesh = *_meshes[n];
1319  const StorageSite& cells = mesh.getCells();
1320  const int nCells = cells.getCountLevel1();
1321  //double pi(acos(-1.0));
1322  const double pi=_options.pi;
1323  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1324  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1325  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1326  const int numFields= _quadrature.getDirCount();
1327 
1328  for(int j=0;j< numFields;j++){
1329  Field& fnd = *_dsfPtr.dsf[j];
1330  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
1331  for(int c=0; c<nCells;c++){
1332  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)
1333  +(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);
1334  }
1335  if (_options.transient)
1336  {
1337  Field& fnd1 = *_dsfPtr1.dsf[j];
1338  TArray& f1 = dynamic_cast< TArray&>(fnd1[cells]);
1339  for (int c=0;c<nCells;c++)
1340  f1[c] = f[c];
1341  //cout << "discretization order " << _options.timeDiscretizationOrder << endl ;
1342  if (_options.timeDiscretizationOrder > 1)
1343  {
1344  Field& fnd2 = *_dsfPtr2.dsf[j];
1345  TArray& f2 = dynamic_cast< TArray&>(fnd2[cells]);
1346  for (int c=0;c<nCells;c++)
1347  f2[c] = f[c];
1348  }
1349  }
1350  }
1351  }
1352  }
1353  KineticBCMap& getBCMap() {return _bcMap;}
1354  KineticVCMap& getVCMap() {return _vcMap;}
1355 
1356  KineticModelOptions<T>& getOptions() {return _options;}
1357  const map<int, vector<int> >& getFaceReflectionArrayMap() const { return _faceReflectionArrayMap;}
1358 
1359 
1360 
1361  // const vector<int>& vecReflection = _faceReflectionArrayMap[faceID]
1362 map<string,shared_ptr<ArrayBase> >&
1363  getPersistenceData()
1364  {
1365  _persistenceData.clear();
1366 
1367  Array<int>* niterArray = new Array<int>(1);
1368  (*niterArray)[0] = _niters;
1369  _persistenceData["niters"]=shared_ptr<ArrayBase>(niterArray);
1370 
1371  if (_initialKmodelNorm)
1372  {
1373  // _persistenceData["initialKmodelNorm"] =_initialKmodelNorm->getArrayPtr(_macroFields.pressure);
1374  const Field& dsfField = *_dsfPtr.dsf[0];
1375  _persistenceData["initialKmodelNorm"] =_initialKmodelNorm->getArrayPtr(dsfField);
1376 
1377  }
1378  else
1379  {
1380  Array<T>* xArray = new Array<T>(1);
1381  xArray->zero();
1382  _persistenceData["initialKmodelNorm"]=shared_ptr<ArrayBase>(xArray);
1383 
1384  }
1385  return _persistenceData;
1386  }
1387 
1388  void restart()
1389  {
1390  if (_persistenceData.find("niters") != _persistenceData.end())
1391  {
1392  shared_ptr<ArrayBase> rp = _persistenceData["niters"];
1393  ArrayBase& r = *rp;
1394  Array<int>& niterArray = dynamic_cast<Array<int>& >(r);
1395  _niters = niterArray[0];
1396  }
1397 
1398  if (_persistenceData.find("initialKmodelNorm") != _persistenceData.end())
1399  {
1400  shared_ptr<ArrayBase> r = _persistenceData["initialKmodelNorm"];
1401  _initialKmodelNorm = MFRPtr(new MultiFieldReduction());
1402  Field& dsfField = *_dsfPtr.dsf[0];
1403  _initialKmodelNorm->addArray(dsfField,r);
1404  //_initialKmodelNorm->addArray(_dsfPtr,r);
1405  }
1406  }
1407 
1408  void SetBoundaryConditions()
1409  {
1410  const int numMeshes = _meshes.size();
1411  for (int n=0; n<numMeshes; n++)
1412  {
1413  const Mesh& mesh = *_meshes[n];
1414 
1415  KineticVC<T> *vc(new KineticVC<T>());
1416  vc->vcType = "flow";
1417  _vcMap[mesh.getID()] = vc;
1418  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1419  {
1420  const FaceGroup& fg = *fgPtr;
1421  if (_bcMap.find(fg.id) == _bcMap.end())
1422  {
1423  KineticBC<T> *bc(new KineticBC<T>());
1424 
1425  _bcMap[fg.id] = bc;
1426  if((fg.groupType == "wall"))
1427  {
1428  bc->bcType = "WallBC";
1429  }
1430  else if((fg.groupType == "realwall"))
1431  {
1432  bc->bcType = "RealWallBC";
1433  }
1434  else if (fg.groupType == "velocity-inlet")
1435  {
1436  bc->bcType = "VelocityInletBC";
1437  }
1438  else if (fg.groupType == "pressure-inlet")
1439  {
1440  bc->bcType = "PressureInletBC";
1441  }
1442  else if (fg.groupType == "pressure-outlet")
1443  {
1444  bc->bcType = "PressureOutletBC";
1445  }
1446  else if ((fg.groupType == "symmetry"))
1447  {
1448  bc->bcType = "SymmetryBC";
1449  }
1450 
1451  else if((fg.groupType =="zero-gradient "))
1452  {
1453  bc->bcType = "ZeroGradBC";
1454  }
1455  else
1456  throw CException("KineticModel: unknown face group type "
1457  + fg.groupType);
1458  }
1459  }
1460  /*
1461  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1462  {
1463  const FaceGroup& fg = *fgPtr;
1464  if (_bcMap.find(fg.id) == _bcMap.end())
1465  { KineticBC<T> *bc(new KineticBC<T>());
1466 
1467  _bcMap[fg.id] = bc;
1468 
1469  if ((fg.groupType == "NSinterface"))
1470  {
1471  bc->bcType = "NSInterfaceBC";
1472  }
1473  }
1474  }
1475  */
1476  }
1477  }
1478 
1479 
1480  void initKineticModelLinearization(LinearSystem& ls, int direction)
1481  {
1482  const int numMeshes = _meshes.size();
1483  for (int n=0; n<numMeshes; n++)
1484  {
1485  const Mesh& mesh = *_meshes[n];
1486  const StorageSite& cells = mesh.getCells();
1487 
1488  Field& fnd = *_dsfPtr.dsf[direction];
1489  //const TArray& f = dynamic_cast<const TArray&>(fnd[cells]);
1490 
1491  MultiField::ArrayIndex vIndex(&fnd,&cells); //dsf is in 1 direction
1492 
1493  ls.getX().addArray(vIndex,fnd.getArrayPtr(cells));
1494 
1495  const CRConnectivity& cellCells = mesh.getCellCells();
1496 
1497  shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells)); //diagonal off diagonal, phi DiagTensorT3,T,VectorT3 for mometum in fvmbase
1498 
1499  ls.getMatrix().addMatrix(vIndex,vIndex,m);
1500  }
1501  //mesh.getBoudaryfaces and interfacefaces?
1502  }
1503 
1504 
1505  void linearizeKineticModel(LinearSystem& ls, int direction)
1506  {
1507  // _velocityGradientModel.compute();
1508  DiscrList discretizations;
1509 
1510 
1511 
1512  Field& fnd = *_dsfPtr.dsf[direction];
1513  //Field& feq = *_dsfEqPtr.dsf[direction];
1514  //if(_options.ESBGK_fgamma){feq = *_dsfEqPtrES.dsf[direction];}
1515  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1516  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1517  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1518 
1519 
1520  if(_options.fgamma==2){
1521  Field& feqES = *_dsfEqPtrES.dsf[direction];
1522  shared_ptr<Discretization>
1523  sdEQ(new CollisionTermDiscretization<T,T,T>
1524  (_meshes, _geomFields,
1525  fnd,feqES,
1526  _macroFields.collisionFrequency));
1527  discretizations.push_back(sdEQ);}
1528  else{
1529  Field& feq = *_dsfEqPtr.dsf[direction];
1530  shared_ptr<Discretization>
1531  sd(new CollisionTermDiscretization<T,T,T>
1532  (_meshes, _geomFields,
1533  fnd,feq,
1534  _macroFields.collisionFrequency));
1535  discretizations.push_back(sd);
1536  }
1537 
1538 
1539  shared_ptr<Discretization>
1540  cd(new ConvectionDiscretization_Kmodel<T,T,T>
1541  (_meshes,
1542  _geomFields,
1543  fnd,
1544  cx[direction],
1545  cy[direction],
1546  cz[direction],
1547  _options.CentralDifference
1548  //_options["nonDimLt"],
1549  //_options["nonDimLx"],_options["nonDimLy"],_options["nonDimLz"],
1550  ));
1551  discretizations.push_back(cd);
1552 
1553  if (_options.transient)
1554  {
1555  // const int direction(0);
1556  Field& fnd = *_dsfPtr.dsf[direction];
1557  Field& fnd1 = *_dsfPtr1.dsf[direction];
1558  Field& fnd2 = *_dsfPtr2.dsf[direction];
1559 
1560 
1561  shared_ptr<Discretization>
1562  td(new TimeDerivativeDiscretization_Kmodel<T,T,T>
1563  (_meshes,_geomFields,
1564  fnd,fnd1,fnd2,
1565  _options["timeStep"],
1566  _options["nonDimLt"],
1567  _options.timeDiscretizationOrder));
1568 
1569  discretizations.push_back(td);
1570 
1571  }
1572  if(_options.ibm_enable){
1573  shared_ptr<Discretization>
1574  ibm(new GenericKineticIBDiscretization<T,T,T>
1575  (_meshes,_geomFields,fnd,
1576  cx[direction],
1577  cy[direction],
1578  cz[direction],
1579  _macroFields));
1580  discretizations.push_back(ibm);
1581  }
1582 
1583  Linearizer linearizer;
1584 
1585  linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
1586  ls.getX(), ls.getB());
1587 
1588  // boundary conditions
1589  const double epsilon=_options.epsilon_ES;
1590  const int numMeshes = _meshes.size();
1591  for (int n=0; n<numMeshes; n++)
1592  {
1593  const Mesh& mesh = *_meshes[n];
1594  const StorageSite& cells = mesh.getCells();
1595 
1596  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1597  {
1598  const FaceGroup& fg = *fgPtr;
1599  const StorageSite& faces = fg.site;
1600  const int nFaces = faces.getCount();
1601 
1602  const KineticBC<T>& bc = *_bcMap[fg.id];
1603  Field& fnd = *_dsfPtr.dsf[direction]; //field in a direction
1604  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
1605  BaseGenericKineticBCS<T,T,T> gkbc(faces, mesh, _geomFields,
1606  fnd,
1607  ls.getMatrix(),
1608  ls.getX(),
1609  ls.getB());
1610 
1611 
1612  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1613  const Field& areaMagField = _geomFields.areaMag;
1614  const TArray& faceAreaMag = dynamic_cast<const TArray &>(areaMagField[faces]);
1615  const Field& areaField = _geomFields.area;
1616  const VectorT3Array& faceArea=dynamic_cast<const VectorT3Array&>(areaField[faces]);
1617 
1618  FloatValEvaluator<VectorT3>
1619  bVelocity(bc.getVal("specifiedXVelocity"),
1620  bc.getVal("specifiedYVelocity"),
1621  bc.getVal("specifiedZVelocity"),
1622  faces);
1623 
1624  if (( bc.bcType == "ZeroGradBC"))
1625  {
1626  for(int f=0; f< nFaces; f++)
1627  {const int c1= faceCells(f,1);// boundary cell
1628  T bvalue =dsf[c1];
1629  gkbc.applyDirichletBC(f,bvalue);
1630  // gkbc.applyExtrapolationBC(f);
1631  }
1632  }
1633  else if (( bc.bcType == "VelocityInletBC")){
1634  for(int f=0; f< nFaces; f++)
1635  {
1636  const VectorT3 en = faceArea[f]/faceAreaMag[f];
1637  const T c_dot_en = cx[direction]*en[0]+cy[direction]*en[1]+cz[direction]*en[2];
1638  if(c_dot_en < T_Scalar(epsilon))
1639  //incoming direction - dirchlet bc
1640  { const int c1= faceCells(f,1);
1641  T bvalue =dsf[c1];
1642  gkbc.applyDirichletBC(f,bvalue);
1643  }
1644  else{
1645  //outgoing direction - extrapolation bc
1646  gkbc.applyExtrapolationBC(f);
1647  }
1648  }
1649  }
1650  //if ((bc.bcType == "WallBC")||(bc.bcType=="PressureInletBC")|| (bc.bcType=="PressureOutletBC")|| (bc.bcType=="VelocityInletBC"))
1651 
1652  else{
1653  for(int f=0; f< nFaces; f++)
1654  {
1655  const VectorT3 en = faceArea[f]/faceAreaMag[f];
1656  const T c_dot_en = cx[direction]*en[0]+cy[direction]*en[1]+cz[direction]*en[2];
1657  const VectorT3 WallVelocity = bVelocity[f];
1658  const T uwall = WallVelocity[0]; const T vwall = WallVelocity[1];
1659  const T wwall = WallVelocity[2]; const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
1660  if(c_dot_en -wallV_dot_en < T_Scalar(epsilon))
1661  //incoming direction - dirchlet bc
1662  { const int c1= faceCells(f,1);// boundary cell
1663  T bvalue =dsf[c1];
1664  gkbc.applyDirichletBC(f,bvalue);
1665  }
1666  else{
1667  //outgoing direction - extrapolation bc
1668  gkbc.applyExtrapolationBC(f);
1669  }
1670 
1671  }
1672  }
1673 
1674 
1675 
1676  }
1677 
1678  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1679  {
1680  const FaceGroup& fg = *fgPtr;
1681  const StorageSite& faces = fg.site;
1682  const int nFaces = faces.getCount();
1683  //const KineticBC<T>& bc = *_bcMap[fg.id];
1684  //Field& fnd = *_dsfPtr.dsf[direction]; //field in a direction
1685  //TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
1686  BaseGenericKineticBCS<T,T,T> gkbc(faces, mesh, _geomFields,
1687  fnd,
1688  ls.getMatrix(),
1689  ls.getX(),
1690  ls.getB());
1691  for(int f=0; f< nFaces; f++)
1692  {gkbc.applyInterfaceBC(f);}//do nothign
1693 
1694  }
1695 
1696 
1697  }
1698  }
1699 
1700 
1701  void computeIBFaceDsf(const StorageSite& solidFaces,const int method,const int RelaxDistribution=0)
1702  {
1703  typedef CRMatrixTranspose<T,T,T> IMatrix;
1704  typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
1705  if (method==1){
1706  const int numMeshes = _meshes.size();
1707  const int numFields= _quadrature.getDirCount();
1708  for (int direction = 0; direction < numFields; direction++)
1709  {
1710  Field& fnd = *_dsfPtr.dsf[direction];
1711  const TArray& pV =
1712  dynamic_cast<const TArray&>(fnd[solidFaces]);
1713  #ifdef FVM_PARALLEL
1714  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,pV.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
1715  #endif
1716 
1717  for (int n=0; n<numMeshes; n++)
1718  {
1719  const Mesh& mesh = *_meshes[n];
1720  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
1721 
1722  const StorageSite& cells = mesh.getCells();
1723  const StorageSite& ibFaces = mesh.getIBFaces();
1724 
1725  GeomFields::SSPair key1(&ibFaces,&cells);
1726  const IMatrix& mIC =
1727  dynamic_cast<const IMatrix&>
1728  (*_geomFields._interpolationMatrices[key1]);
1729 
1730  IMatrix mICV(mIC);
1731 
1732 
1733  GeomFields::SSPair key2(&ibFaces,&solidFaces);
1734  const IMatrix& mIP =
1735  dynamic_cast<const IMatrix&>
1736  (*_geomFields._interpolationMatrices[key2]);
1737 
1738  IMatrix mIPV(mIP);
1739 
1740 
1741  shared_ptr<TArray> ibV(new TArray(ibFaces.getCount()));
1742 
1743  const TArray& cV =
1744  dynamic_cast<const TArray&>(fnd[cells]);
1745 
1746  ibV->zero();
1747 
1748  mICV.multiplyAndAdd(*ibV,cV);
1749  mIPV.multiplyAndAdd(*ibV,pV);
1750 
1751 #if 0
1752  ofstream debugFile;
1753  stringstream ss(stringstream::in | stringstream::out);
1754  ss << MPI::COMM_WORLD.Get_rank();
1755  string fname1 = "IBVelocity_proc" + ss.str() + ".dat";
1756  debugFile.open(fname1.c_str());
1757 
1758  //debug use
1759  const Array<int>& ibFaceList = mesh.getIBFaceList();
1760  const StorageSite& faces = mesh.getFaces();
1761  const VectorT3Array& faceCentroid =
1762  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
1763  const double angV = 1.0;
1764  VectorT3 center;
1765  center[0]=0.;
1766  center[1]=0.;
1767  center[2]=0.;
1768 
1769  for(int f=0; f<ibFaces.getCount();f++){
1770  int fID = ibFaceList[f];
1771  debugFile << "f=" << f << setw(10) << " fID = " << fID << " faceCentroid = " << faceCentroid[fID] << " ibV = " << (*ibV)[f] << endl;
1772  }
1773 
1774  debugFile.close();
1775 #endif
1776 
1777  fnd.addArray(ibFaces,ibV);
1778  }
1779  }
1780  }
1781  for (int n=0; n<numMeshes; n++)
1782  {
1783  const Mesh& mesh = *_meshes[n];
1784  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
1785 
1786  const StorageSite& cells = mesh.getCells();
1787  const StorageSite& ibFaces = mesh.getIBFaces();
1788 
1789  GeomFields::SSPair key1(&ibFaces,&cells);
1790  const IMatrix& mIC =
1791  dynamic_cast<const IMatrix&>
1792  (*_geomFields._interpolationMatrices[key1]);
1793 
1794  IMatrixV3 mICV3(mIC);
1795 
1796  GeomFields::SSPair key2(&ibFaces,&solidFaces);
1797  const IMatrix& mIP =
1798  dynamic_cast<const IMatrix&>
1799  (*_geomFields._interpolationMatrices[key2]);
1800 
1801  IMatrixV3 mIPV3(mIP);
1802 
1803  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
1804 
1805 
1806  const VectorT3Array& cVel =
1807  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
1808  const VectorT3Array& sVel =
1809  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
1810 
1811  ibVvel->zero();
1812 
1813  //velocity interpolation (cells+solidfaces)
1814  mICV3.multiplyAndAdd(*ibVvel,cVel);
1815  mIPV3.multiplyAndAdd(*ibVvel,sVel);
1816  _macroFields.velocity.addArray(ibFaces,ibVvel);
1817 
1818 
1819  }
1820  }
1821  }
1822 
1823  if (method==2){
1824  const int numMeshes = _meshes.size();
1825  const int nSolidFaces = solidFaces.getCount();
1826 
1827  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
1828  *muSolid =0;
1829  _macroFields.viscosity.addArray(solidFaces,muSolid);
1830 
1831  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
1832  *nueSolid =0;
1833  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
1834 
1835  const T rho_init=_options["rho_init"];
1836  const T T_init= _options["T_init"];
1837  const T mu_w= _options["mu_w"];
1838  const T Tmuref= _options["Tmuref"];
1839  const T muref= _options["muref"];
1840  const T R=8314.0/_options["molecularWeight"];
1841  const T nondim_length=_options["nonDimLt"];
1842 
1843  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
1844 
1845  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
1846  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
1847  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
1848  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
1849 
1850  for(int c=0; c<nSolidFaces;c++)
1851  {
1852  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
1853  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
1854  }
1855 
1856  if(_options.fgamma==2){
1857  for(int c=0; c<nSolidFaces;c++)
1858  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
1859  }
1860 
1861 
1862 
1863 
1864  for (int n=0; n<numMeshes; n++)
1865  {
1866  const Mesh& mesh = *_meshes[n];
1867  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
1868 
1869  const StorageSite& cells = mesh.getCells();
1870  const StorageSite& ibFaces = mesh.getIBFaces();
1871 
1872  GeomFields::SSPair key1(&ibFaces,&cells);
1873  const IMatrix& mIC =
1874  dynamic_cast<const IMatrix&>
1875  (*_geomFields._interpolationMatrices[key1]);
1876 
1877  IMatrix mICV(mIC);
1878  IMatrixV3 mICV3(mIC);
1879 
1880  GeomFields::SSPair key2(&ibFaces,&solidFaces);
1881  const IMatrix& mIP =
1882  dynamic_cast<const IMatrix&>
1883  (*_geomFields._interpolationMatrices[key2]);
1884 
1885  IMatrix mIPV(mIP);
1886  IMatrixV3 mIPV3(mIP);
1887 
1888  shared_ptr<TArray> ibVtemp(new TArray(ibFaces.getCount()));
1889  shared_ptr<TArray> ibVnue(new TArray(ibFaces.getCount()));
1890  shared_ptr<TArray> ibVdensity(new TArray(ibFaces.getCount()));
1891  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
1892 
1893  const TArray& cTemp =
1894  dynamic_cast<TArray&>(_macroFields.temperature[cells]);
1895  const VectorT3Array& cVel =
1896  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
1897  const TArray& cDensity =
1898  dynamic_cast<TArray&>(_macroFields.density[cells]);
1899  const TArray& sDensity =
1900  dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
1901  const TArray& cNue =
1902  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
1903  const TArray& sNue =
1904  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
1905  const TArray& sTemp =
1906  dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
1907  const VectorT3Array& sVel =
1908  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
1909 
1910  ibVnue->zero();
1911  ibVtemp->zero();
1912  ibVvel->zero();
1913  ibVdensity->zero();
1914 
1915  //nue interpolation (cells)
1916  mICV.multiplyAndAdd(*ibVnue,cNue);
1917  mIPV.multiplyAndAdd(*ibVnue,sNue);
1918  _macroFields.collisionFrequency.addArray(ibFaces,ibVnue);
1919  //temperature interpolation (cells+solidfaces)
1920  mICV.multiplyAndAdd(*ibVtemp,cTemp);
1921  mIPV.multiplyAndAdd(*ibVtemp,sTemp);
1922  _macroFields.temperature.addArray(ibFaces,ibVtemp);
1923  //density interpolation (cells+solidfaces)
1924  mICV.multiplyAndAdd(*ibVdensity,cDensity);
1925  mIPV.multiplyAndAdd(*ibVdensity,sDensity);
1926  _macroFields.density.addArray(ibFaces,ibVdensity);
1927  //velocity interpolation (cells+solidfaces)
1928  mICV3.multiplyAndAdd(*ibVvel,cVel);
1929  mIPV3.multiplyAndAdd(*ibVvel,sVel);
1930  _macroFields.velocity.addArray(ibFaces,ibVvel);
1931 
1932 
1933  }
1934  }
1935 
1936  const int f_out = 3;
1937  if (f_out ==1){
1938  //Step 2 Find fgamma using macroparameters
1939  const int numMeshes = _meshes.size();
1940  for (int n=0; n<numMeshes; n++)
1941  {
1942  const Mesh& mesh = *_meshes[n];
1943  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
1944  const int numDirections = _quadrature.getDirCount();
1945  const StorageSite& ibFaces = mesh.getIBFaces();
1946  const int nibFaces=ibFaces.getCount();
1947  const double pi=_options.pi;
1948  const TArray& ibTemp =
1949  dynamic_cast<TArray&>(_macroFields.temperature[ibFaces]);
1950  const VectorT3Array& ibVel =
1951  dynamic_cast<VectorT3Array&>(_macroFields.velocity[ibFaces]);
1952  const TArray& ibDensity =
1953  dynamic_cast<TArray&>(_macroFields.density[ibFaces]);
1954 
1955  for (int j=0; j<numDirections; j++)
1956  {
1957  shared_ptr<TArray> ibFndPtrEqES(new TArray(nibFaces));
1958  TArray& ibFndEqES= *ibFndPtrEqES;
1959  ibFndPtrEqES->zero();
1960  Field& fndEqES = *_dsfEqPtrES.dsf[j];
1961 
1962  for (int i=0; i<nibFaces; i++)
1963  {
1964  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
1965  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
1966  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
1967  const T ibu = ibVel[i][0];
1968  const T ibv = ibVel[i][1];
1969  const T ibw = ibVel[i][2];
1970  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]);
1971  }
1972  fndEqES.addArray(ibFaces,ibFndPtrEqES);
1973  }
1974  }
1975  }
1976  }
1977  else if(f_out==2)
1978  {
1979  //Step 2 Find fgamma using interpolation (only ESBGK for now)
1980  for (int n=0; n<numMeshes; n++)
1981  {
1982  const int numFields= _quadrature.getDirCount();
1983  for (int direction = 0; direction < numFields; direction++)
1984  {
1985  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
1986  const Mesh& mesh = *_meshes[n];
1987  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
1988 
1989  const StorageSite& cells = mesh.getCells();
1990  const StorageSite& ibFaces = mesh.getIBFaces();
1991 
1992  GeomFields::SSPair key1(&ibFaces,&cells);
1993  const IMatrix& mIC =
1994  dynamic_cast<const IMatrix&>
1995  (*_geomFields._interpolationMatrices[key1]);
1996 
1997  IMatrix mICV(mIC);
1998 
1999  GeomFields::SSPair key2(&ibFaces,&solidFaces);
2000  const IMatrix& mIP =
2001  dynamic_cast<const IMatrix&>
2002  (*_geomFields._interpolationMatrices[key2]);
2003 
2004  IMatrix mIPV(mIP);
2005 
2006  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
2007 
2008  const TArray& cf =
2009  dynamic_cast<const TArray&>(fndEqES[cells]);
2010 
2011  ibVf->zero();
2012 
2013  //distribution function interpolation (cells)
2014  mICV.multiplyAndAdd(*ibVf,cf);
2015  fndEqES.addArray(ibFaces,ibVf);
2016 
2017  }
2018  }
2019  }
2020  }
2021  //Step3: Relax Distribution function from ibfaces to solid face
2022  for (int n=0; n<numMeshes; n++)
2023  {
2024  const int numDirections = _quadrature.getDirCount();
2025  for (int j=0; j<numDirections; j++)
2026  {
2027  Field& fnd = *_dsfPtr.dsf[j];
2028  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
2029 
2030 //#ifdef FVM_PARALLEL
2031 // MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsf.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
2032 //#endif
2033  const Mesh& mesh = *_meshes[n];
2034  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2035  const StorageSite& ibFaces = mesh.getIBFaces();
2036  TArray& dsfIB = dynamic_cast< TArray&>(fnd[ibFaces]);
2037  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2038  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[ibFaces]);
2039  const StorageSite& faces = mesh.getFaces();
2040  const StorageSite& cells = mesh.getCells();
2041  const CRConnectivity& faceCells = mesh.getAllFaceCells();
2042  const CRConnectivity& ibFacesTosolidFaces
2043  = mesh.getConnectivity(ibFaces,solidFaces);
2044  const IntArray& ibFaceIndices = mesh.getIBFaceList();
2045  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2046  const IntArray& sFCRow = ibFacesTosolidFaces.getRow();
2047  const IntArray& sFCCol = ibFacesTosolidFaces.getCol();
2048  const int nibFaces = ibFaces.getCount();
2049  const int nFaces = faces.getCount();
2050  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2051  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2052  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2053  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2054 
2055  for(int f=0; f<nibFaces; f++)
2056  {
2057  dsfIB[f]=0.0;
2058  double distIBSolidInvSum(0.0);
2059  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2060  {
2061  const int c = sFCCol[nc];
2062  const int faceIB= ibFaceIndices[f];
2063  const VectorT3Array& solidFaceCentroid =
2064  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2065  const VectorT3Array& faceCentroid =
2066  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2067 
2068  double distIBSolid (0.0);
2069  // based on distance - will be thought
2070  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
2071  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
2072  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
2073  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
2074  }
2075  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2076  {
2077  const int c = sFCCol[nc];
2078  const VectorT3Array& solidFaceCentroid =
2079  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2080  const VectorT3Array& faceCentroid =
2081  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2082  const TArray& nue =
2083  dynamic_cast<TArray&>(_macroFields.collisionFrequency[ibFaces]);
2084  const int faceIB= ibFaceIndices[f];
2085  // const T coeff = iCoeffs[nc];
2086  double time_to_wall (0.0);
2087  double distIBSolid (0.0);
2088  const T uwall = v[c][0];
2089  const T vwall = v[c][1];
2090  const T wwall = v[c][2];
2091  // based on distance - will be thought
2092  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
2093  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
2094  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
2095  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])));
2096  if(time_to_wall<0)
2097  time_to_wall = 0;
2098 
2099  dsfIB[f] += (dsfEqES[f]-(dsfEqES[f]-dsf[c])*exp(-time_to_wall*nue[f]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
2100 
2101  }
2102 
2103  }
2104  }
2105  }
2106  }
2107  }
2108  if (method==3){
2109  const int numMeshes = _meshes.size();
2110  for (int n=0; n<numMeshes; n++)
2111  {
2112  const Mesh& mesh = *_meshes[n];
2113  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2114 
2115  const StorageSite& cells = mesh.getCells();
2116  const StorageSite& ibFaces = mesh.getIBFaces();
2117 
2118  GeomFields::SSPair key1(&ibFaces,&cells);
2119  const IMatrix& mIC =
2120  dynamic_cast<const IMatrix&>
2121  (*_geomFields._interpolationMatrices[key1]);
2122 
2123  IMatrixV3 mICV3(mIC);
2124 
2125  GeomFields::SSPair key2(&ibFaces,&solidFaces);
2126  const IMatrix& mIP =
2127  dynamic_cast<const IMatrix&>
2128  (*_geomFields._interpolationMatrices[key2]);
2129 
2130  IMatrixV3 mIPV3(mIP);
2131 
2132  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
2133 
2134 
2135  const VectorT3Array& cVel =
2136  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
2137  const VectorT3Array& sVel =
2138  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2139 
2140  ibVvel->zero();
2141 
2142  //velocity interpolation (cells+solidfaces)
2143  mICV3.multiplyAndAdd(*ibVvel,cVel);
2144  mIPV3.multiplyAndAdd(*ibVvel,sVel);
2145  _macroFields.velocity.addArray(ibFaces,ibVvel);
2146 
2147 
2148  }
2149  }
2150  const int nSolidFaces = solidFaces.getCount();
2151 
2152  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
2153  *muSolid =0;
2154  _macroFields.viscosity.addArray(solidFaces,muSolid);
2155 
2156  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
2157  *nueSolid =0;
2158  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
2159 
2160  const T rho_init=_options["rho_init"];
2161  const T T_init= _options["T_init"];
2162  const T mu_w= _options["mu_w"];
2163  const T Tmuref= _options["Tmuref"];
2164  const T muref= _options["muref"];
2165  const T R=8314.0/_options["molecularWeight"];
2166  const T nondim_length=_options["nonDimLt"];
2167 
2168  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
2169 
2170  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
2171  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
2172  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
2173  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
2174 
2175  for(int c=0; c<nSolidFaces;c++)
2176  {
2177  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
2178  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
2179  }
2180 
2181  if(_options.fgamma==2){
2182  for(int c=0; c<nSolidFaces;c++)
2183  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
2184  }
2185 
2186  //Step 2 Find fgamma using interpolation (only ESBGK for now)
2187  const int numFields= _quadrature.getDirCount();
2188  for (int direction = 0; direction < numFields; direction++)
2189  {
2190  shared_ptr<TArray> ibVf(new TArray(solidFaces.getCount()));
2191  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
2192  ibVf->zero();
2193  for (int n=0; n<numMeshes; n++)
2194  {
2195  const Mesh& mesh = *_meshes[n];
2196  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2197 
2198  const StorageSite& cells = mesh.getCells();
2199  const StorageSite& ibFaces = mesh.getIBFaces();
2200 
2201  GeomFields::SSPair key1(&solidFaces,&cells);
2202  const IMatrix& mIC =
2203  dynamic_cast<const IMatrix&>
2204  (*_geomFields._interpolationMatrices[key1]);
2205 
2206  IMatrix mICV(mIC);
2207 
2208  const TArray& cf =
2209  dynamic_cast<const TArray&>(fndEqES[cells]);
2210 
2211  ibVf->zero();
2212 
2213  //distribution function interpolation (cells)
2214  mICV.multiplyAndAdd(*ibVf,cf);
2215  }
2216  }
2217  fndEqES.addArray(solidFaces,ibVf);
2218  }
2219  for (int n=0; n<numMeshes; n++)
2220  {
2221  const int numDirections = _quadrature.getDirCount();
2222  for (int j=0; j<numDirections; j++)
2223  {
2224  Field& fnd = *_dsfPtr.dsf[j];
2225  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
2226  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2227  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[solidFaces]);
2228 #ifdef FVM_PARALLEL
2229  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsf.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
2230  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,dsfEqES.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
2231 #endif
2232 
2233  const Mesh& mesh = *_meshes[n];
2234  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2235  const StorageSite& ibFaces = mesh.getIBFaces();
2236  const StorageSite& faces = mesh.getFaces();
2237  const StorageSite& cells = mesh.getCells();
2238  const CRConnectivity& faceCells = mesh.getAllFaceCells();
2239  const CRConnectivity& ibFacesTosolidFaces
2240  = mesh.getConnectivity(ibFaces,solidFaces);
2241  const IntArray& ibFaceIndices = mesh.getIBFaceList();
2242  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2243  const IntArray& sFCRow = ibFacesTosolidFaces.getRow();
2244  const IntArray& sFCCol = ibFacesTosolidFaces.getCol();
2245  const int nibFaces = ibFaces.getCount();
2246  const int nFaces = faces.getCount();
2247  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2248  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2249  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2250  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2251  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
2252  ibVf->zero();
2253  TArray& ibVfA= *ibVf;
2254  for(int f=0; f<nibFaces; f++)
2255  {
2256  double distIBSolidInvSum(0.0);
2257  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2258  {
2259  const int c = sFCCol[nc];
2260  const int faceIB= ibFaceIndices[f];
2261  const VectorT3Array& solidFaceCentroid =
2262  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2263  const VectorT3Array& faceCentroid =
2264  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2265 
2266  double distIBSolid (0.0);
2267  // based on distance - will be thought
2268  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
2269  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
2270  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
2271  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
2272  }
2273  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2274  {
2275  const int c = sFCCol[nc];
2276  const VectorT3Array& solidFaceCentroid =
2277  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2278  const VectorT3Array& faceCentroid =
2279  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2280  const TArray& nue =
2281  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
2282  const TArray& nueC =
2283  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
2284  const int faceIB= ibFaceIndices[f];
2285  const T uwall = v[c][0];
2286  const T vwall = v[c][1];
2287  const T wwall = v[c][2];
2288  // const T coeff = iCoeffs[nc];
2289  double time_to_wall (0.0);
2290  double distIBSolid (0.0);
2291  // based on distance - will be thought
2292  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[c][0]),2)+
2293  pow((faceCentroid[faceIB][1]-solidFaceCentroid[c][1]),2)+
2294  pow((faceCentroid[faceIB][2]-solidFaceCentroid[c][2]),2));
2295  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])));
2296  if(time_to_wall<0)
2297  time_to_wall = 0;
2298  ibVfA[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
2299  }
2300 
2301  }
2302 
2303  fnd.addArray(ibFaces,ibVf);
2304  }
2305  }
2306  }
2307  }
2308  }
2309 
2310  void computeSolidFacePressure(const StorageSite& solidFaces)
2311  {
2312  typedef CRMatrixTranspose<T,T,T> IMatrix;
2313  const int numMeshes = _meshes.size();
2314  for (int n=0; n<numMeshes; n++)
2315  {
2316  shared_ptr<TArray> ibP(new TArray(solidFaces.getCount()));
2317  ibP->zero();
2318  const Mesh& mesh = *_meshes[n];
2319  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2320 
2321  const StorageSite& cells = mesh.getCells();
2322 
2323  GeomFields::SSPair key1(&solidFaces,&cells);
2324  const IMatrix& mIC =
2325  dynamic_cast<const IMatrix&>
2326  (*_geomFields._interpolationMatrices[key1]);
2327 
2328  IMatrix mICV(mIC);
2329 
2330  const TArray& cP =
2331  dynamic_cast<TArray&>(_macroFields.pressure[cells]);
2332 
2333  ibP->zero();
2334 
2335 
2336  //nue interpolation (cells)
2337  mICV.multiplyAndAdd(*ibP,cP);
2338  }
2339  _macroFields.pressure.addArray(solidFaces,ibP);
2340  }
2341 #ifdef FVM_PARALLEL
2342  TArray& pressure = dynamic_cast<TArray&>(_macroFields.pressure[solidFaces]);
2343  MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE,pressure.getData(),solidFaces.getCount() , MPI::DOUBLE, MPI::SUM);
2344 #endif
2345  }
2346 
2347  void computeSolidFaceDsf(const StorageSite& solidFaces,const int method,const int RelaxDistribution=0)
2348  {
2349  typedef CRMatrixTranspose<T,T,T> IMatrix;
2350  typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
2351  const int numFields= _quadrature.getDirCount();
2352  if (method==1){
2353  const int numMeshes = _meshes.size();
2354  for (int direction = 0; direction < numFields; direction++) {
2355  Field& fnd = *_dsfPtr.dsf[direction];
2356  shared_ptr<TArray> ibV(new TArray(solidFaces.getCount()));
2357  ibV->zero();
2358  for (int n=0; n<numMeshes; n++)
2359  {
2360  const Mesh& mesh = *_meshes[n];
2361  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2362 
2363  const StorageSite& cells = mesh.getCells();
2364 
2365 
2366  GeomFields::SSPair key1(&solidFaces,&cells);
2367  const IMatrix& mIC =
2368  dynamic_cast<const IMatrix&>
2369  (*_geomFields._interpolationMatrices[key1]);
2370 
2371  IMatrix mICV(mIC);
2372  const TArray& cV =
2373  dynamic_cast<const TArray&>(fnd[cells]);
2374 
2375  ibV->zero();
2376 
2377  mICV.multiplyAndAdd(*ibV,cV);
2378 #if 0
2379  ofstream debugFile;
2380  stringstream ss(stringstream::in | stringstream::out);
2381  ss << MPI::COMM_WORLD.Get_rank();
2382  string fname1 = "IBVelocity_proc" + ss.str() + ".dat";
2383  debugFile.open(fname1.c_str());
2384 
2385  //debug use
2386  const Array<int>& ibFaceList = mesh.getIBFaceList();
2387  const StorageSite& faces = mesh.getFaces();
2388  const VectorT3Array& faceCentroid =
2389  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2390  const double angV = 1.0;
2391  VectorT3 center;
2392  center[0]=0.;
2393  center[1]=0.;
2394  center[2]=0.;
2395 
2396  for(int f=0; f<ibFaces.getCount();f++){
2397  int fID = ibFaceList[f];
2398  debugFile << "f=" << f << setw(10) << " fID = " << fID << " faceCentroid = " << faceCentroid[fID] << " ibV = " << (*ibV)[f] << endl;
2399  }
2400 
2401  debugFile.close();
2402 #endif
2403 
2404  }
2405 
2406  }
2407  fnd.addArray(solidFaces,ibV);
2408  }
2409  }
2410  if (method==2){
2411  // Step0: Compute Interpolation Matrices from (only) Cells to IBFaces
2412  const int numMeshes = _meshes.size();
2413  for (int n=0; n<numMeshes; n++)
2414  {
2415  const Mesh& mesh = *_meshes[n];
2416 
2417  const int numFields= _quadrature.getDirCount();
2418  for (int direction = 0; direction < numFields; direction++)
2419  {
2420  Field& fnd = *_dsfPtr.dsf[direction];
2421  Field& fndEqES = *_dsfEqPtrES.dsf[direction];
2422  const Mesh& mesh = *_meshes[n];
2423  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2424 
2425  const StorageSite& cells = mesh.getCells();
2426  const StorageSite& faces = mesh.getFaces();
2427  const StorageSite& ibFaces = mesh.getIBFaces();
2428 
2429  GeomFields::SSPair key1(&faces,&cells);
2430  const IMatrix& mIC =
2431  dynamic_cast<const IMatrix&>
2432  (*_geomFields._interpolationMatrices[key1]);
2433 
2434  IMatrix mICV(mIC);
2435 
2436  const TArray& cf =
2437  dynamic_cast<const TArray&>(fnd[cells]);
2438  const TArray& cfEq =
2439  dynamic_cast<const TArray&>(fndEqES[cells]);
2440 
2441  shared_ptr<TArray> ibVf(new TArray(ibFaces.getCount()));
2442  ibVf->zero();
2443 
2444  if (_options.fgamma==2){
2445  shared_ptr<TArray> ibVfEq(new TArray(ibFaces.getCount()));
2446  ibVfEq->zero();
2447  mICV.multiplyAndAdd(*ibVfEq,cfEq);
2448  fndEqES.addArray(ibFaces,ibVfEq);
2449  }
2450  mICV.multiplyAndAdd(*ibVf,cf);
2451  fnd.addArray(ibFaces,ibVf);
2452 
2453  }
2454  }
2455  }
2456  const int nSolidFaces = solidFaces.getCount();
2457 
2458  shared_ptr<TArray> muSolid(new TArray(nSolidFaces));
2459  *muSolid =0;
2460  _macroFields.viscosity.addArray(solidFaces,muSolid);
2461 
2462  shared_ptr<TArray> nueSolid(new TArray(nSolidFaces));
2463  *nueSolid =0;
2464  _macroFields.collisionFrequency.addArray(solidFaces,nueSolid);
2465 
2466  const T rho_init=_options["rho_init"];
2467  const T T_init= _options["T_init"];
2468  const T mu_w= _options["mu_w"];
2469  const T Tmuref= _options["Tmuref"];
2470  const T muref= _options["muref"];
2471  const T R=8314.0/_options["molecularWeight"];
2472  const T nondim_length=_options["nonDimLt"];
2473 
2474  const T mu0=rho_init*R* T_init*nondim_length/pow(2*R* T_init,0.5);
2475 
2476  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
2477  TArray& viscosity = dynamic_cast<TArray&>(_macroFields.viscosity[solidFaces]);
2478  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
2479  TArray& collisionFrequency = dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
2480 
2481  for(int c=0; c<nSolidFaces;c++)
2482  {
2483  viscosity[c]= muref*pow(temperature[c]*T_init/ Tmuref,mu_w); // viscosity power law
2484  collisionFrequency[c]=density[c]*temperature[c]/viscosity[c]*mu0;
2485  }
2486 
2487  if(_options.fgamma==2){
2488  for(int c=0; c<nSolidFaces;c++)
2489  collisionFrequency[c]=_options.Prandtl*collisionFrequency[c];
2490  }
2491 
2492  //Step 1 Interpolate Macroparameters and f to IBface
2493  for (int n=0; n<numMeshes; n++)
2494  {
2495  const Mesh& mesh = *_meshes[n];
2496  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2497 
2498  const StorageSite& cells = mesh.getCells();
2499  const StorageSite& ibFaces = mesh.getIBFaces();
2500 
2501  GeomFields::SSPair key1(&ibFaces,&cells);
2502  const IMatrix& mIC =
2503  dynamic_cast<const IMatrix&>
2504  (*_geomFields._interpolationMatrices[key1]);
2505 
2506  IMatrix mICV(mIC);
2507  IMatrixV3 mICV3(mIC);
2508 
2509  GeomFields::SSPair key2(&ibFaces,&solidFaces);
2510  const IMatrix& mIP =
2511  dynamic_cast<const IMatrix&>
2512  (*_geomFields._interpolationMatrices[key2]);
2513 
2514  IMatrix mIPV(mIP);
2515  IMatrixV3 mIPV3(mIP);
2516 
2517  shared_ptr<TArray> ibVtemp(new TArray(ibFaces.getCount()));
2518  shared_ptr<TArray> ibVnue(new TArray(ibFaces.getCount()));
2519  shared_ptr<TArray> ibVdensity(new TArray(ibFaces.getCount()));
2520  shared_ptr<VectorT3Array> ibVvel(new VectorT3Array(ibFaces.getCount()));
2521 
2522  const TArray& cTemp =
2523  dynamic_cast<TArray&>(_macroFields.temperature[cells]);
2524  const VectorT3Array& cVel =
2525  dynamic_cast<VectorT3Array&>(_macroFields.velocity[cells]);
2526  const TArray& cDensity =
2527  dynamic_cast<TArray&>(_macroFields.density[cells]);
2528  const TArray& sDensity =
2529  dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
2530  const TArray& cNue =
2531  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
2532  const TArray& sNue =
2533  dynamic_cast<TArray&>(_macroFields.collisionFrequency[solidFaces]);
2534  const TArray& sTemp =
2535  dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
2536  const VectorT3Array& sVel =
2537  dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2538 
2539  ibVnue->zero();
2540  ibVtemp->zero();
2541  ibVvel->zero();
2542  ibVdensity->zero();
2543 
2544  //nue interpolation (cells)
2545  mICV.multiplyAndAdd(*ibVnue,cNue);
2546  mIPV.multiplyAndAdd(*ibVnue,sNue);
2547  _macroFields.collisionFrequency.addArray(ibFaces,ibVnue);
2548  //temperature interpolation (cells+solidfaces)
2549  mICV.multiplyAndAdd(*ibVtemp,cTemp);
2550  mIPV.multiplyAndAdd(*ibVtemp,sTemp);
2551  _macroFields.temperature.addArray(ibFaces,ibVtemp);
2552  //density interpolation (cells+solidfaces)
2553  mICV.multiplyAndAdd(*ibVdensity,cDensity);
2554  mIPV.multiplyAndAdd(*ibVdensity,sDensity);
2555  _macroFields.density.addArray(ibFaces,ibVdensity);
2556  //velocity interpolation (cells+solidfaces)
2557  mICV3.multiplyAndAdd(*ibVvel,cVel);
2558  mIPV3.multiplyAndAdd(*ibVvel,sVel);
2559  _macroFields.velocity.addArray(ibFaces,ibVvel);
2560 
2561 
2562  }
2563  }
2564 
2565  if (_options.fgamma==1){
2566  //Step 2 Find fgamma using macroparameters
2567  const int numMeshes = _meshes.size();
2568  for (int n=0; n<numMeshes; n++)
2569  {
2570  const Mesh& mesh = *_meshes[n];
2571  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2572  const int numDirections = _quadrature.getDirCount();
2573  const StorageSite& ibFaces = mesh.getIBFaces();
2574  const int nibFaces=ibFaces.getCount();
2575  const double pi=_options.pi;
2576  const TArray& ibTemp =
2577  dynamic_cast<TArray&>(_macroFields.temperature[ibFaces]);
2578  const VectorT3Array& ibVel =
2579  dynamic_cast<VectorT3Array&>(_macroFields.velocity[ibFaces]);
2580  const TArray& ibDensity =
2581  dynamic_cast<TArray&>(_macroFields.density[ibFaces]);
2582 
2583  for (int j=0; j<numDirections; j++)
2584  {
2585  shared_ptr<TArray> ibFndPtrEqES(new TArray(nibFaces));
2586  TArray& ibFndEqES= *ibFndPtrEqES;
2587 
2588  ibFndPtrEqES->zero();
2589 
2590  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2591 
2592  for (int i=0; i<nibFaces; i++)
2593  {
2594  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2595  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2596  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2597  const T ibu = ibVel[i][0];
2598  const T ibv = ibVel[i][1];
2599  const T ibw = ibVel[i][2];
2600  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]);
2601  }
2602  fndEqES.addArray(ibFaces,ibFndPtrEqES);
2603  }
2604  }
2605  }
2606  }
2607 
2608  //Step3: Relax Distribution function from ibfaces to solid face
2609  const int numDirections = _quadrature.getDirCount();
2610  for (int j=0; j<numDirections; j++)
2611  {
2612  const int nSolidFaces = solidFaces.getCount();
2613  shared_ptr<TArray> solidFndPtr(new TArray(nSolidFaces));
2614  solidFndPtr->zero();
2615  TArray& solidFnd= *solidFndPtr;
2616  Field& fnd = *_dsfPtr.dsf[j];
2617  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2618  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2619  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2620  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2621  for (int n=0; n<numMeshes; n++)
2622  {
2623  const Mesh& mesh = *_meshes[n];
2624  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2625  const StorageSite& ibFaces = mesh.getIBFaces();
2626  const CRConnectivity& solidFacesToibFaces
2627  = mesh.getConnectivity(solidFaces,ibFaces);
2628  const IntArray& ibFaceIndices = mesh.getIBFaceList();
2629  const IntArray& sFCRow = solidFacesToibFaces.getRow();
2630  const IntArray& sFCCol = solidFacesToibFaces.getCol();
2631  TArray& dsf = dynamic_cast< TArray&>(fnd[ibFaces]);
2632  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[ibFaces]);
2633  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2634 
2635  for(int f=0; f<nSolidFaces; f++)
2636  {
2637  double distIBSolidInvSum(0.0);
2638  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2639  {
2640  const StorageSite& faces = mesh.getFaces();
2641  const int c = sFCCol[nc];
2642  const int faceIB= ibFaceIndices[c];
2643  const VectorT3Array& solidFaceCentroid =
2644  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2645  const VectorT3Array& faceCentroid =
2646  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2647 
2648  double distIBSolid (0.0);
2649  // based on distance - will be thought
2650  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[f][0]),2)+
2651  pow((faceCentroid[faceIB][1]-solidFaceCentroid[f][1]),2)+
2652  pow((faceCentroid[faceIB][2]-solidFaceCentroid[f][2]),2));
2653  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
2654  }
2655  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2656  {
2657  const int c = sFCCol[nc];
2658  const StorageSite& faces = mesh.getFaces();
2659  const VectorT3Array& solidFaceCentroid =
2660  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2661  const VectorT3Array& faceCentroid =
2662  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
2663  const int faceIB= ibFaceIndices[c];
2664  T time_to_wall (0.0);
2665  T distIBSolid (0.0);
2666  distIBSolid = sqrt(pow((faceCentroid[faceIB][0]-solidFaceCentroid[f][0]),2)+
2667  pow((faceCentroid[faceIB][1]-solidFaceCentroid[f][1]),2)+
2668  pow((faceCentroid[faceIB][2]-solidFaceCentroid[f][2]),2));
2669  // based on distance - will be thought
2670  const T uwall = v[f][0];
2671  const T vwall = v[f][1];
2672  const T wwall = v[f][2];
2673  const TArray& nue =
2674  dynamic_cast<TArray&>(_macroFields.collisionFrequency[ibFaces]);
2675  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])));
2676  if(time_to_wall<0)
2677  time_to_wall = 0;
2678 
2679  solidFnd[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
2680  }
2681 
2682  }
2683  }
2684  }
2685  fnd.addArray(solidFaces,solidFndPtr);
2686  }
2687  }
2688  if (method==3){
2689  const int numDirections = _quadrature.getDirCount();
2690  for (int j=0; j<numDirections; j++)
2691  {
2692  const int nSolidFaces = solidFaces.getCount();
2693  shared_ptr<TArray> solidFndPtr(new TArray(nSolidFaces));
2694  solidFndPtr->zero();
2695  TArray& solidFnd= *solidFndPtr;
2696  Field& fnd = *_dsfPtr.dsf[j];
2697  Field& fndEqES = *_dsfEqPtrES.dsf[j];
2698  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2699  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2700  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2701  const int numMeshes = _meshes.size();
2702  for (int n=0; n<numMeshes; n++)
2703  {
2704  const Mesh& mesh = *_meshes[n];
2705  if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
2706  const StorageSite& cells = mesh.getCells();
2707  const StorageSite& ibFaces = mesh.getIBFaces();
2708  const CRConnectivity& solidFacesToCells
2709  = mesh.getConnectivity(solidFaces,cells);
2710  const IntArray& sFCRow = solidFacesToCells.getRow();
2711  const IntArray& sFCCol = solidFacesToCells.getCol();
2712  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2713  TArray& dsfEqES = dynamic_cast< TArray&>(fndEqES[cells]);
2714  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
2715 
2716  for(int f=0; f<nSolidFaces; f++)
2717  {
2718  double distIBSolidInvSum(0.0);
2719  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2720  {
2721  const StorageSite& faces = mesh.getFaces();
2722  const int c = sFCCol[nc];
2723  const VectorT3Array& solidFaceCentroid =
2724  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2725  const VectorT3Array& faceCentroid =
2726  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[cells]);
2727 
2728  double distIBSolid (0.0);
2729  // based on distance - will be thought
2730  distIBSolid = sqrt(pow((faceCentroid[c][0]-solidFaceCentroid[f][0]),2)+
2731  pow((faceCentroid[c][1]-solidFaceCentroid[f][1]),2)+
2732  pow((faceCentroid[c][2]-solidFaceCentroid[f][2]),2));
2733  distIBSolidInvSum += 1/pow(distIBSolid,RelaxDistribution);
2734  }
2735  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
2736  {
2737  const int c = sFCCol[nc];
2738  const StorageSite& faces = mesh.getFaces();
2739  const VectorT3Array& solidFaceCentroid =
2740  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[solidFaces]);
2741  const VectorT3Array& faceCentroid =
2742  dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[cells]);
2743  T time_to_wall (0.0);
2744  T distIBSolid (0.0);
2745  const T uwall = v[f][0];
2746  const T vwall = v[f][1];
2747  const T wwall = v[f][2];
2748  distIBSolid = sqrt(pow((faceCentroid[c][0]-solidFaceCentroid[f][0]),2)+
2749  pow((faceCentroid[c][1]-solidFaceCentroid[f][1]),2)+
2750  pow((faceCentroid[c][2]-solidFaceCentroid[f][2]),2));
2751  // based on distance - will be thought
2752 
2753  const TArray& nue =
2754  dynamic_cast<TArray&>(_macroFields.collisionFrequency[cells]);
2755  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])));
2756  if(time_to_wall<0)
2757  time_to_wall = 0;
2758 
2759  solidFnd[f] += (dsfEqES[c]-(dsfEqES[c]-dsf[c])*exp(-time_to_wall*nue[c]))/(pow(distIBSolid,RelaxDistribution)*distIBSolidInvSum);
2760  }
2761 
2762  }
2763  }
2764  }
2765  fnd.addArray(solidFaces,solidFndPtr);
2766 
2767  }
2768  }
2769  }
2770 
2771 
2772 
2773 
2774  void correctMassDeficit()
2775  {
2776 
2777  const int numMeshes = _meshes.size();
2778  T netFlux(0.0);
2779 
2780  for (int n=0; n<numMeshes; n++)
2781  {
2782  const Mesh& mesh = *_meshes[n];
2783  const StorageSite& ibFaces = mesh.getIBFaces();
2784  const StorageSite& cells = mesh.getCells();
2785  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2786  const StorageSite& faces = mesh.getFaces();
2787  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
2788  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
2789  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
2790  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
2791  const int numDirections = _quadrature.getDirCount();
2792  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[faces]);
2793  const VectorT3Array& faceArea =
2794  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
2795  const TArray& faceAreaMag =
2796  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
2797  const CRConnectivity& faceCells = mesh.getAllFaceCells();
2798  const int nibFaces = ibFaces.getCount();
2799 
2800  for(int f=0; f<nibFaces; f++)
2801  {
2802  const int c0 = faceCells(f,0);
2803  const int c1 = faceCells(f,1);
2804  if (((ibType[c0] == Mesh::IBTYPE_FLUID) && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
2805  ((ibType[c1] == Mesh::IBTYPE_FLUID) && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
2806  {
2807  const int ibFace = ibFaceIndex[f];
2808  if (ibFace < 0)
2809  throw CException("invalid ib face index");
2810  if (ibType[c0] == Mesh::IBTYPE_FLUID)
2811  {
2812  const VectorT3 en = faceArea[f]/faceAreaMag[f];
2813  for (int j=0; j<numDirections; j++)
2814  {
2815  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
2816  Field& fnd = *_dsfPtr.dsf[j];
2817  TArray& dsf = dynamic_cast<TArray&>(fnd[ibFaces]);
2818 
2819  netFlux -= dsf[f]*c_dot_en*wts[j]/abs(c_dot_en);
2820  }
2821 
2822  }
2823  else
2824  {
2825  const VectorT3 en = faceArea[f]/faceAreaMag[f];
2826  for (int j=0; j<numDirections; j++)
2827  {
2828  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
2829  Field& fnd = *_dsfPtr.dsf[j];
2830  TArray& dsf = dynamic_cast<TArray&>(fnd[ibFaces]);
2831 
2832  netFlux += dsf[f]*c_dot_en*wts[j]/abs(c_dot_en);
2833  }
2834  }
2835  }
2836  }
2837  }
2838 
2839 #ifdef FVM_PARALLEL
2840  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &netFlux, 1, MPI::DOUBLE, MPI::SUM);
2841 #endif
2842 
2843  T volumeSum(0.);
2844 
2845  for (int n=0; n<numMeshes; n++)
2846  {
2847  const Mesh& mesh = *_meshes[n];
2848  const StorageSite& cells = mesh.getCells();
2849  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2850  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
2851  for(int c=0; c<cells.getSelfCount(); c++)
2852  if (ibType[c] == Mesh::IBTYPE_FLUID)
2853  volumeSum += cellVolume[c];
2854  }
2855 #ifdef FVM_PARALLEL
2856  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &volumeSum, 1, MPI::DOUBLE, MPI::SUM);
2857 #endif
2858 
2859  netFlux /= volumeSum;
2860  for (int n=0; n<numMeshes; n++)
2861  {
2862  const Mesh& mesh = *_meshes[n];
2863  const StorageSite& cells = mesh.getCells();
2864  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2865  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
2866  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
2867 
2868  for(int c=0; c<cells.getSelfCount(); c++)
2869  {
2870  if (ibType[c] == Mesh::IBTYPE_FLUID){
2871  const int numDirections = _quadrature.getDirCount();
2872  T cellMass(0.0);
2873  for (int j=0; j<numDirections; j++)
2874  {
2875  Field& fnd = *_dsfPtr.dsf[j];
2876  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2877  cellMass += wts[j]*dsf[c];
2878  }
2879 
2880  for (int j=0; j<numDirections; j++)
2881  {
2882  Field& fnd = *_dsfPtr.dsf[j];
2883  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2884  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
2885  TArray& fgam = dynamic_cast< TArray&>(feqES[cells]);
2886  fgam[c] = fgam[c]*(1+netFlux*cellVolume[c]/cellMass);
2887  dsf[c] = dsf[c]*(1+netFlux*cellVolume[c]/cellMass);
2888  }
2889  }
2890 
2891  }
2892  }
2893  }
2894 
2895  void correctMassDeficit2(double n1,double n2)
2896  {
2897 
2898  const int numMeshes = _meshes.size();
2899  T netFlux(0.0);
2900 
2901  netFlux=n2-n1;
2902 
2903  T volumeSum(0.);
2904 
2905  for (int n=0; n<numMeshes; n++)
2906  {
2907  const Mesh& mesh = *_meshes[n];
2908  const StorageSite& cells = mesh.getCells();
2909  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2910  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
2911  for(int c=0; c<cells.getSelfCount(); c++)
2912  if (ibType[c] == Mesh::IBTYPE_FLUID)
2913  volumeSum += cellVolume[c];
2914  }
2915 #ifdef FVM_PARALLEL
2916  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &volumeSum, 1, MPI::DOUBLE, MPI::SUM);
2917 #endif
2918 
2919  netFlux /= volumeSum;
2920  for (int n=0; n<numMeshes; n++)
2921  {
2922  const Mesh& mesh = *_meshes[n];
2923  const StorageSite& cells = mesh.getCells();
2924  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2925  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
2926  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
2927 
2928  for(int c=0; c<cells.getSelfCount(); c++)
2929  {
2930  if (ibType[c] == Mesh::IBTYPE_FLUID){
2931  const int numDirections = _quadrature.getDirCount();
2932  T cellMass(0.0);
2933  for (int j=0; j<numDirections; j++)
2934  {
2935  Field& fnd = *_dsfPtr.dsf[j];
2936  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2937  cellMass += wts[j]*dsf[c];
2938  }
2939 
2940  for (int j=0; j<numDirections; j++)
2941  {
2942  Field& fnd = *_dsfPtr.dsf[j];
2943  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2944  Field& feqES = *_dsfEqPtrES.dsf[j]; //for fgamma_2
2945  TArray& fgam = dynamic_cast< TArray&>(feqES[cells]);
2946  fgam[c] = fgam[c]*(1+netFlux*cellVolume[c]/cellMass);
2947  dsf[c] = dsf[c]*(1+netFlux*cellVolume[c]/cellMass);
2948  }
2949  }
2950 
2951  }
2952  }
2953  }
2954  const double ConservationofMassCheck()
2955  {
2956  const int numMeshes = _meshes.size();
2957  T ndens_tot(0.0) ;
2958 
2959  for (int n=0; n<numMeshes; n++)
2960  {
2961  const Mesh& mesh = *_meshes[n];
2962  const StorageSite& cells = mesh.getCells();
2963  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
2964  const StorageSite& faces = mesh.getFaces();
2965  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
2966  const int numDirections = _quadrature.getDirCount();
2967  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[faces]);
2968  const VectorT3Array& faceArea =
2969  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
2970  const TArray& faceAreaMag =
2971  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
2972  const CRConnectivity& faceCells = mesh.getAllFaceCells();
2973  const int nFaces = faces.getCount();
2974 
2975  for(int c=0; c<cells.getCountLevel1(); c++)
2976  {
2977  if (ibType[c] == Mesh::IBTYPE_FLUID)
2978  {
2979  for (int j=0; j<numDirections; j++)
2980  {
2981  Field& fnd = *_dsfPtr.dsf[j];
2982  TArray& dsf = dynamic_cast< TArray&>(fnd[cells]);
2983  ndens_tot += wts[j]*dsf[c];
2984  }
2985  }
2986  }
2987  }
2988  cout << "Hello, I have" << ndens_tot << "number density";
2989  return ndens_tot;
2990  }
2991 
2992  void ConservationofMFSolid(const StorageSite& solidFaces, const int output =0, bool perUnitArea=0) const
2993  {
2994  FILE * pFile;
2995 
2996  shared_ptr<VectorT6Array> Stress(new VectorT6Array(solidFaces.getCount()));
2997  Stress->zero();
2998  _macroFields.Stress.addArray(solidFaces,Stress);
2999 
3000  shared_ptr<VectorT3Array> Force(new VectorT3Array(solidFaces.getCount()));
3001  Force->zero();
3002  _macroFields.force.addArray(solidFaces,Force);
3003 
3004  shared_ptr<TArray> TemperatureIB(new TArray(solidFaces.getCount()));
3005  TemperatureIB->zero();
3006  _macroFields.temperatureIB.addArray(solidFaces,TemperatureIB);
3007 
3008  const double pi=_options.pi;
3009  const double epsilon=_options.epsilon_ES;
3010  const int nSolidFaces = solidFaces.getCount();
3011  for (int i=0; i<nSolidFaces; i++)
3012  {
3013  const int numDirections = _quadrature.getDirCount();
3014  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3015  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3016  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3017  const VectorT3Array& solidFaceCentroid =
3018  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[solidFaces]);
3019  const VectorT3Array& solidFaceArea =
3020  dynamic_cast<const VectorT3Array&>(_geomFields.area[solidFaces]);
3021  const TArray& solidFaceAreaMag =
3022  dynamic_cast<const TArray&>(_geomFields.areaMag[solidFaces]);
3023  const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
3024  VectorT3Array& v = dynamic_cast<VectorT3Array&>(_macroFields.velocity[solidFaces]);
3025  TArray& density = dynamic_cast<TArray&>(_macroFields.density[solidFaces]);
3026  TArray& temperature = dynamic_cast<TArray&>(_macroFields.temperature[solidFaces]);
3027  TArray& tempIB = dynamic_cast<TArray&>(_macroFields.temperatureIB[solidFaces]);
3028  VectorT3Array& force = dynamic_cast<VectorT3Array&>(_macroFields.force[solidFaces]);
3029  VectorT6Array& stress = dynamic_cast<VectorT6Array&>(_macroFields.Stress[solidFaces]);
3030  const T Lx=_options["nonDimLx"];
3031  const T Ly=_options["nonDimLy"];
3032  const T Lz=_options["nonDimLz"];
3033 
3034  const T uwall = v[i][0];
3035  const T vwall = v[i][1];
3036  const T wwall = v[i][2];
3037  const T Twall = temperature[i];
3038 
3039  T Nmr(0.0) ;
3040  T Dmr(0.0) ;
3041  T incomFlux(0.0);
3042  T TempWall(0.0);
3043  T mWall(0.0);
3044 
3045  for (int j=0; j<numDirections; j++)
3046  {
3047  Field& fnd = *_dsfPtr.dsf[j];
3048  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
3049  const VectorT3 en = solidFaceArea[i]/solidFaceAreaMag[i];
3050  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
3051  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
3052  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);
3053 
3054  if (c_dot_en-wallV_dot_en > 0) //outgoing
3055  {
3056  Dmr = Dmr - fwall*wts[j]*(c_dot_en-wallV_dot_en);
3057  incomFlux=incomFlux-dsf[i]*wts[j]*(c_dot_en-wallV_dot_en);
3058  }
3059  else
3060  {
3061  Nmr = Nmr + dsf[i]*wts[j]*(c_dot_en-wallV_dot_en);
3062  if(output==1){
3063  TempWall = TempWall + dsf[i]*(pow(cx[j]-uwall,2.0)+pow(cy[j]-vwall,2.0)+pow(cz[j]-wwall,2.0))*wts[j]*0.5;
3064  density[i] = density[i] + dsf[i]*wts[j]*0.5;
3065  stress[i][0] +=pow((cx[j]-uwall),2.0)*dsf[i]*wts[j]*0.5;
3066  stress[i][1] +=pow((cy[j]-vwall),2.0)*dsf[i]*wts[j]*0.5;
3067  stress[i][2] +=pow((cz[j]-wwall),2.0)*dsf[i]*wts[j]*0.5;
3068  stress[i][3] +=(cx[j]-uwall)*(cy[j]-vwall)*dsf[i]*wts[j]*0.5;
3069  stress[i][4] +=(cy[j]-vwall)*(cz[j]-wwall)*dsf[i]*wts[j]*0.5;
3070  stress[i][5] +=(cx[j]-uwall)*(cz[j]-wwall)*dsf[i]*wts[j]*0.5;
3071  }
3072  }
3073  }
3074  const T nwall = Nmr/Dmr; // incoming wall number density for initializing Maxwellian
3075 
3076  for (int j=0; j<numDirections; j++)
3077  {
3078  Field& fnd = *_dsfPtr.dsf[j];
3079  TArray& dsf = dynamic_cast< TArray&>(fnd[solidFaces]);
3080  const VectorT3 en = solidFaceArea[i]/solidFaceAreaMag[i];
3081  const T c_dot_en = cx[j]*en[0]+cy[j]*en[1]+cz[j]*en[2];
3082  const T wallV_dot_en = uwall*en[0]+vwall*en[1]+wwall*en[2];
3083  if (c_dot_en-wallV_dot_en > 0)
3084  {
3085  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);
3086  if(output==1){
3087  TempWall = TempWall + dsf[i]*(pow(cx[j]-uwall,2.0)+pow(cy[j]-vwall,2.0)+pow(cz[j]-wwall,2.0))*wts[j]*0.5;
3088  density[i] = density[i] + dsf[i]*wts[j]*0.5;
3089  stress[i][0] +=pow((cx[j]-uwall),2.0)*dsf[i]*wts[j]*0.5;
3090  stress[i][1] +=pow((cy[j]-vwall),2.0)*dsf[i]*wts[j]*0.5;
3091  stress[i][2] +=pow((cz[j]-wwall),2.0)*dsf[i]*wts[j]*0.5;
3092  stress[i][3] +=(cx[j]-uwall)*(cy[j]-vwall)*dsf[i]*wts[j]*0.5;
3093  stress[i][4] +=(cy[j]-vwall)*(cz[j]-wwall)*dsf[i]*wts[j]*0.5;
3094  stress[i][5] +=(cx[j]-uwall)*(cz[j]-wwall)*dsf[i]*wts[j]*0.5;
3095  }
3096  }
3097  else
3098  dsf[i]=dsf[i];
3099  }
3100  if(output==1){
3101  const VectorT3& Af = solidFaceArea[i];
3102  force[i][0] = Af[0]*Ly*Lz*stress[i][0] + Af[1]*Lz*Lx*stress[i][3] + Af[2]*Lx*Ly*stress[i][5];
3103  force[i][1] = Af[0]*Ly*Lz*stress[i][3] + Af[1]*Lz*Lx*stress[i][1] + Af[2]*Lx*Ly*stress[i][4];
3104  force[i][2] = Af[0]*Ly*Lz*stress[i][5] + Af[1]*Lz*Lx*stress[i][4] + Af[2]*Ly*Ly*stress[i][2];
3105 
3106  pFile = fopen("WallTemperature.dat","a");
3107  fprintf(pFile,"%E %E %E %E\n",solidFaceCentroid[i][0],solidFaceCentroid[i][1],solidFaceCentroid[i][2],TempWall);
3108  fclose(pFile);
3109  }
3110  }
3111  }
3112 
3113 void MacroparameterIBCell(const StorageSite& solidFaces) const
3114  {
3115 // typedef CRMatrixTranspose<T,T,T> IMatrix;
3116 // typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
3117  const int numMeshes = _meshes.size();
3118 // for (int n=0; n<numMeshes; n++)
3119 // {
3120 // const Mesh& mesh = *_meshes[n];
3121 // if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
3122 
3123 // const StorageSite& cells = mesh.getCells();
3124 // const StorageSite& ibFaces = mesh.getIBFaces();
3125 
3126 // GeomFields::SSPair key1(&ibFaces,&cells);
3127 // const IMatrix& mIC =
3128 // dynamic_cast<const IMatrix&>
3129 // (*_geomFields._interpolationMatrices[key1]);
3130 
3131 // IMatrix mICV(mIC);
3132 // IMatrixV3 mICV3(mIC);
3133 
3134 // GeomFields::SSPair key2(&ibFaces,&solidFaces);
3135 // const IMatrix& mIP =
3136 // dynamic_cast<const IMatrix&>
3137 // (*_geomFields._interpolationMatrices[key2]);
3138 
3139 // IMatrix mIPV(mIP);
3140 // IMatrixV3 mIPV3(mIP);
3141 
3142 // shared_ptr<TArray> ibVtemp(new TArray(ibFaces.getCount()));
3143 
3144 // const TArray& cTemp =
3145 // dynamic_cast<TArray&>(_macroFields.temperature[cells]);
3146 // const TArray& sTemp =
3147 // dynamic_cast<TArray&>(_macroFields.temperatureIB[solidFaces]);
3148 
3149 // ibVtemp->zero();
3150 
3151  //temperature interpolation (cells+solidfaces)
3152 // mICV.multiplyAndAdd(*ibVtemp,cTemp);
3153 // mIPV.multiplyAndAdd(*ibVtemp,sTemp);
3154 // _macroFields.temperature.addArray(ibFaces,ibVtemp);
3155 // }
3156 // }
3157  for (int n=0;n<numMeshes;n++)
3158  {
3159  const Mesh& mesh = *_meshes[n];
3160  const CRConnectivity& faceCells = mesh.getAllFaceCells();
3161  const StorageSite& ibFaces = mesh.getIBFaces();
3162  const StorageSite& cells = mesh.getCells();
3163  const int nCells = cells.getCountLevel1();
3164  const StorageSite& faces = mesh.getFaces();
3165 // const TArray& TempIB =
3166 // dynamic_cast<TArray&>(_macroFields.temperature[ibFaces]);
3167  TArray& TempB =
3168  dynamic_cast<TArray&>(_macroFields.pressure[cells]);
3169  const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
3170  const IntArray& ibFaceIndex = dynamic_cast<const IntArray&>(_geomFields.ibFaceIndex[faces]);
3171 
3172  TArray xB(nCells);
3173  TArray wB(nCells);
3174 
3175  xB.zero();
3176  wB.zero();
3177 
3178  const int nFaces = faces.getCount();
3179  for(int f=0; f<nFaces; f++)
3180  {
3181  const int c0 = faceCells(f,0);
3182  const int c1 = faceCells(f,1);
3183 
3184  if (((ibType[c0] == Mesh::IBTYPE_FLUID) && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
3185  ((ibType[c1] == Mesh::IBTYPE_FLUID) && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
3186  {
3187  // this is an iBFace, determine which cell is interior and which boundary
3188 
3189  const int ibFace = ibFaceIndex[f];
3190  if (ibFace < 0)
3191  throw CException("invalid ib face index");
3192 // const T xFace = TempIB[ibFace];
3193 
3194  if (ibType[c0] == Mesh::IBTYPE_FLUID)
3195  {
3196  xB[c1] = xB[c0];
3197  wB[c1]++;
3198  }
3199 
3200  else
3201  {
3202  xB[c0] = xB[c1];
3203  wB[c0]++;
3204  }
3205  }
3206  else if ((ibType[c0] == Mesh::IBTYPE_FLUID) &&
3207  (ibType[c1] == Mesh::IBTYPE_FLUID))
3208  {
3209  // leave as is
3210  }
3211  else
3212  {
3213  // leave as is
3214  }
3215  }
3216 
3217  // set the phi for boundary cells as average of the ib face values
3218  for(int c=0; c<nCells; c++)
3219  {
3220  if (wB[c] > 0)
3221  cout << "wb value" << wB[c];
3222  TempB[c] = xB[c] / T_Scalar(wB[c]);
3223 
3224  }
3225 
3226  }
3227  }
3228 
3229  void updateTime()
3230  {
3231  const int numMeshes = _meshes.size();
3232  for (int n=0;n<numMeshes;n++)
3233  {
3234  const Mesh& mesh = *_meshes[n];
3235 
3236  const StorageSite& cells = mesh.getCells();
3237  const int numFields= _quadrature.getDirCount();
3238  //callBoundaryConditions(); //new
3239  for (int direction = 0; direction < numFields; direction++)
3240  {
3241  Field& fnd = *_dsfPtr.dsf[direction];
3242  Field& fndN1 = *_dsfPtr1.dsf[direction];
3243  TArray& f = dynamic_cast<TArray&>(fnd[cells]);
3244  TArray& fN1 = dynamic_cast<TArray&>(fndN1[cells]);
3245  if (_options.timeDiscretizationOrder > 1)
3246  {
3247  Field& fndN2 = *_dsfPtr2.dsf[direction];
3248  TArray& fN2 = dynamic_cast<TArray&>(fndN2[cells]);
3249  fN2 = fN1;
3250  }
3251  fN1 = f;
3252  }
3253 #ifdef FVM_PARALLEL
3254  if ( MPI::COMM_WORLD.Get_rank() == 0 )
3255  {cout << "updated time" <<endl;}
3256 
3257 #endif
3258 #ifndef FVM_PARALLEL
3259  cout << "updated time" <<endl;
3260 #endif
3261  //ComputeMacroparameters(); //update macroparameters
3262  //ComputeCollisionfrequency();
3263  //if (_options.fgamma==0){initializeMaxwellianEq();}
3264  //else{ EquilibriumDistributionBGK();}
3265  //if (_options.fgamma==2){EquilibriumDistributionESBGK();}
3266 
3267  }
3268  }
3269 
3270 
3271  void callBoundaryConditions()
3272 
3273  {
3274  const int numMeshes = _meshes.size();
3275  for (int n=0; n<numMeshes; n++)
3276  {
3277  const Mesh& mesh = *_meshes[n];
3278 
3279  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3280  {
3281  const FaceGroup& fg = *fgPtr;
3282  const StorageSite& faces = fg.site;
3283  //const int nFaces = faces.getCount();
3284 
3285  const KineticBC<T>& bc = *_bcMap[fg.id];
3286 
3287  KineticBoundaryConditions<T,T,T> kbc(faces, mesh,_geomFields,_quadrature,_macroFields,_dsfPtr);
3288 
3289  FloatValEvaluator<VectorT3> bVelocity(bc.getVal("specifiedXVelocity"),
3290  bc.getVal("specifiedYVelocity"),
3291  bc.getVal("specifiedZVelocity"),
3292  faces);
3293  // FloatValEvaluator<StressTensor<T>>
3294  FloatValEvaluator<VectorT3> bStress(bc.getVal("specifiedTauxx"),bc.getVal("specifiedTauyy"),
3295  bc.getVal("specifiedTauzz"),faces);
3296  //bc.getVal("specifiedTauxy"),
3297  // bc.getVal("specifiedTauyz"),bc.getVal("specifiedTauzx"),faces);
3298 
3299  FloatValEvaluator<T> mdot(bc.getVal("specifiedMassFlowRate"),faces);
3300  FloatValEvaluator<T> bTemperature(bc.getVal("specifiedTemperature"),faces);
3301  FloatValEvaluator<T> bPressure(bc.getVal("specifiedPressure"),faces);
3302  FloatValEvaluator<T> accomCoeff(bc.getVal("accommodationCoefficient"),faces);
3303  if (bc.bcType == "WallBC")
3304  {
3305  kbc.applyDiffuseWallBC(bVelocity,bTemperature);
3306  }
3307  else if (bc.bcType == "RealWallBC")
3308  {
3309  //kbc.applyRealWallBC(bVelocity,bTemperature,accomCoeff);
3310  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(fg.id);
3311  const vector<int>& vecReflection=(*pos).second;
3312  kbc.applyRealWallBC(bVelocity,bTemperature,accomCoeff,vecReflection);
3313  }
3314  else if(bc.bcType=="SymmetryBC")
3315  {
3316  //kbc.applySpecularWallBC(); //old boundary works only for cartesian-type quadrature
3317  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(fg.id);
3318  const vector<int>& vecReflection=(*pos).second;
3319  kbc.applySpecularWallBC(vecReflection);
3320  }
3321  else if(bc.bcType=="ZeroGradBC")
3322  {
3323  kbc.applyZeroGradientBC();
3324  }
3325  else if(bc.bcType=="PressureInletBC")
3326  {
3327  kbc.applyPressureInletBC(bTemperature,bPressure);
3328  }
3329 
3330  else if(bc.bcType=="VelocityInletBC")
3331  {
3332  map<int, vector<int> >::iterator pos = _faceReflectionArrayMap.find(fg.id);
3333  const vector<int>& vecReflection=(*pos).second;
3334  kbc.applyInletBC(bVelocity,bTemperature,mdot,vecReflection);
3335  }
3336  else if(bc.bcType=="PressureOutletBC")
3337  {
3338  kbc.applyPressureOutletBC(bTemperature,bPressure);
3339  }
3340 
3341  }
3342  foreach(const FaceGroupPtr igPtr, mesh.getInterfaceGroups())
3343  {
3344 
3345  const FaceGroup& ig = *igPtr;
3346  const StorageSite& faces = ig.site;
3347  //const int nFaces = faces.getCount();
3348 
3349  //const KineticBC<T>& bc = *_bcMap[fg.id];
3350 
3351  KineticBoundaryConditions<T,T,T> kbc(faces, mesh,_geomFields,_quadrature,_macroFields,_dsfPtr);
3352  if(ig.groupType=="NSinterface")
3353  {
3354  kbc.applyNSInterfaceBC();//bTemperature,bPressure,bVelocity,bStress);
3355  }
3356  }
3357  }//end of loop through meshes
3358  }
3359 
3360  bool advance(const int niter,const int updated=0)
3361  {
3362  for(int n=0; n<niter; n++)
3363  {
3364  const int N123 =_quadrature.getDirCount();
3365  MFRPtr rNorm;
3366  //MFRPtr vNorm;
3367  //const TArray& cx= dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3368  //const TArray& wts= dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
3369 
3370  //callBoundaryConditions();
3371 
3372  /*
3373  _macroFields.velocity.syncLocal();
3374  _macroFields.temperature.syncLocal();
3375  _macroFields.density .syncLocal();
3376  */
3377  for(int direction=0; direction<N123;direction++)
3378  {
3379  LinearSystem ls;
3380  initKineticModelLinearization(ls, direction);
3381  ls.initAssembly();
3382  linearizeKineticModel(ls,direction);
3383  ls.initSolve();
3384 
3385  //const T* kNorm=_options.getKineticLinearSolver().solve(ls);
3386  MFRPtr kNorm(_options.getKineticLinearSolver().solve(ls));
3387 
3388 
3389  if (!rNorm)
3390  rNorm = kNorm;
3391  else
3392  {
3393  // find the array for the 0the direction residual and
3394  // add the current residual to it
3395  Field& fn0 = *_dsfPtr.dsf[0];
3396  Field& fnd = *_dsfPtr.dsf[direction];
3397 
3398  ArrayBase& rArray0 = (*rNorm)[fn0];
3399  ArrayBase& rArrayd = (*kNorm)[fnd];//*wts[direction];
3400  rArray0 += rArrayd;//*wts[direction];
3401 
3402  // ArrayBase& vArray0 = (*vNorm)[fn0];
3403  // vArray0 += rArrayd;//*cx[direction]*wts[direction];
3404  }
3405 
3406  ls.postSolve();
3407  ls.updateSolution();
3408 
3409  _options.getKineticLinearSolver().cleanup();
3410 
3411  }
3412 
3413 
3414 
3415  if (!_initialKmodelNorm||updated==1) _initialKmodelNorm = rNorm;
3416  //if (!_initialKmodelvNorm) _initialKmodelvNorm = vNorm;
3417  if (_niters < 5||updated==2)
3418  {
3419  _initialKmodelNorm->setMax(*rNorm);
3420  // _initialKmodelvNorm->setMax(*vNorm);
3421 
3422  }
3423 
3424  MFRPtr normRatio((*rNorm)/(*_initialKmodelNorm));
3425  // MFRPtr vnormRatio((*vNorm)/(*_initialKmodelvNorm));
3426 #ifdef FVM_PARALLEL
3427  if ( MPI::COMM_WORLD.Get_rank() == 0 ){
3428  if (_options.printNormalizedResiduals){
3429  cout << _niters << ": " << *normRatio << endl;}
3430  else{
3431  cout << _niters << ": " << *rNorm <<endl; }}
3432 #endif
3433 #ifndef FVM_PARALLEL
3434  if (_options.printNormalizedResiduals){
3435  cout << _niters << ": " << *normRatio << endl;}
3436  else{
3437  cout << _niters << ": " << *rNorm <<endl; }
3438 #endif
3439  _niters++;
3440  //break here
3441 
3442  callBoundaryConditions();
3443  //cout << "called boundary"<<endl;
3444  ComputeMacroparameters(); //update macroparameters
3445  ComputeCollisionfrequency();
3446 
3447  //update equilibrium distribution function 0-maxwellian, 1-BGK,2-ESBGK
3448  if (_options.fgamma==0){initializeMaxwellianEq();}
3449  else{ EquilibriumDistributionBGK();}
3450  //cout << "called BGk" <<endl;
3451  if (_options.fgamma==2){EquilibriumDistributionESBGK();}
3452 
3453 
3454  if ((*rNorm < _options.absoluteTolerance)||(*normRatio < _options.relativeTolerance )){
3455  //&& ((*vNorm < _options.absoluteTolerance)||(*vnormRatio < _options.relativeTolerance )))
3456  return true;
3457  }
3458 
3459  return false;
3460 
3461 
3462  }
3463 
3464 
3465  }
3466 
3467 
3468 
3469 
3470  void OutputDsfBLOCK(const char* filename)
3471  {
3472  FILE * pFile;
3473  pFile = fopen(filename,"w");
3474  int N1=_quadrature.getNVCount();
3475  int N2=_quadrature.getNthetaCount();
3476  int N3=_quadrature.getNphiCount();
3477  fprintf(pFile,"%s \n", "VARIABLES= cx, cy, cz, f,fEq,fES");
3478  fprintf(pFile, "%s %i %s %i %s %i \n","ZONE I=", N3,",J=",N2,",K=",N1);
3479  fprintf(pFile, "%s \n","F=BLOCK, VARLOCATION=(NODAL,NODAL,NODAL,NODAL,NODAL,NODAL)");
3480  const int numMeshes = _meshes.size();
3481  const int cellno=_options.printCellNumber;
3482  for (int n=0; n<numMeshes; n++){
3483  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3484  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3485  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3486  const int numFields= _quadrature.getDirCount();
3487  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cx[j]);}
3488  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cy[j]);}
3489  for(int j=0;j< numFields;j++){fprintf(pFile,"%E \n",cz[j]);}
3490 
3491  const Mesh& mesh = *_meshes[n];
3492  const StorageSite& cells = mesh.getCells();
3493  for(int j=0;j< numFields;j++){
3494  Field& fnd = *_dsfPtr.dsf[j];
3495  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
3496  fprintf(pFile,"%E\n",f[cellno]);
3497  }
3498  for(int j=0;j< numFields;j++){
3499  Field& fEqnd = *_dsfEqPtr.dsf[j];
3500  TArray& fEq = dynamic_cast< TArray&>(fEqnd[cells]);
3501  fprintf(pFile,"%E\n",fEq[cellno]);
3502  }
3503  if(_options.fgamma==2){
3504  for(int j=0;j< numFields;j++){
3505 
3506  Field& fndEqES = *_dsfEqPtrES.dsf[j];
3507  TArray& fEqES = dynamic_cast< TArray&>(fndEqES[cells]);
3508  fprintf(pFile,"%E\n",fEqES[cellno]);
3509  }}
3510  else{
3511  for(int j=0;j< numFields;j++){
3512 
3513  Field& fndEq = *_dsfEqPtr.dsf[j];
3514  TArray& fEq = dynamic_cast< TArray&>(fndEq[cells]);
3515  fprintf(pFile,"%E\n",fEq[cellno]);
3516  }}
3517 
3518  }
3519  fclose(pFile);
3520  }
3521 
3522 
3523  void computeSurfaceForce(const StorageSite& solidFaces, bool perUnitArea, bool IBM=0)
3524  {
3525  typedef CRMatrixTranspose<T,T,T> IMatrix;
3526 
3527  const int nSolidFaces = solidFaces.getCount();
3528 
3529  boost::shared_ptr<VectorT3Array>
3530  forcePtr( new VectorT3Array(nSolidFaces));
3531  VectorT3Array& force = *forcePtr;
3532 
3533  force.zero();
3534  _macroFields.force.addArray(solidFaces,forcePtr);
3535 
3536  const VectorT3Array& solidFaceArea =
3537  dynamic_cast<const VectorT3Array&>(_geomFields.area[solidFaces]);
3538 
3539  const TArray& solidFaceAreaMag =
3540  dynamic_cast<const TArray&>(_geomFields.areaMag[solidFaces]);
3541 
3542  const int numMeshes = _meshes.size();
3543  for (int n=0; n<numMeshes; n++)
3544  {
3545  const Mesh& mesh = *_meshes[n];
3546  const StorageSite& cells = mesh.getCells();
3547 
3548  const VectorT3Array& v = dynamic_cast<const VectorT3Array&>(_macroFields.velocity[cells]);
3549  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3550  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3551  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3552  const TArray& wts = dynamic_cast<const TArray&>(*_quadrature.dcxyzPtr);
3553 
3554  const CRConnectivity& solidFacesToCells
3555  = mesh.getConnectivity(solidFaces,cells);
3556  const IntArray& sFCRow = solidFacesToCells.getRow();
3557  const IntArray& sFCCol = solidFacesToCells.getCol();
3558 
3559  const T Lx=_options["nonDimLx"];
3560  const T Ly=_options["nonDimLy"];
3561  const T Lz=_options["nonDimLz"];
3562 
3563 
3564  const int N123= _quadrature.getDirCount();
3565 
3566  const int selfCount = cells.getSelfCount();
3567  for(int f=0; f<nSolidFaces; f++){
3568 
3569  StressTensor<T> stress = NumTypeTraits<StressTensor<T> >::getZero();
3570 
3571  if (IBM){
3572  GeomFields::SSPair key1(&solidFaces,&cells);
3573  const IMatrix& mIC =
3574  dynamic_cast<const IMatrix&>
3575  (*_geomFields._interpolationMatrices[key1]);
3576  const Array<T>& iCoeffs = mIC.getCoeff();
3577  for(int j=0;j<N123;j++){
3578  Field& fnd = *_dsfPtr.dsf[j];
3579  const TArray& f_dsf = dynamic_cast<const TArray&>(fnd[cells]);
3580  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
3581  {
3582 
3583  const int c = sFCCol[nc];
3584  const T coeff = iCoeffs[nc];
3585  stress[0] -=coeff*pow((cx[j]-v[c][0]),2.0)*f_dsf[c]*wts[j];
3586  stress[1] -=coeff*pow((cy[j]-v[c][1]),2.0)*f_dsf[c]*wts[j];
3587  stress[2] -=coeff*pow((cz[j]-v[c][2]),2.0)*f_dsf[c]*wts[j];
3588  stress[3] -=coeff*(cx[j]-v[c][0])*(cy[j]-v[c][1])*f_dsf[c]*wts[j];
3589  stress[4] -=coeff*(cy[j]-v[c][1])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
3590  stress[5] -=coeff*(cx[j]-v[c][0])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
3591  }
3592  }
3593  }
3594  else
3595  {
3596  for(int j=0;j<N123;j++){
3597  Field& fnd = *_dsfPtr.dsf[j];
3598  const TArray& f_dsf = dynamic_cast<const TArray&>(fnd[cells]);
3599  for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
3600  {
3601 
3602  const int c = sFCCol[nc];
3603  if ( c < selfCount ){
3604  stress[0] +=pow((cx[j]-v[c][0]),2.0)*f_dsf[c]*wts[j];
3605  stress[1] +=pow((cy[j]-v[c][1]),2.0)*f_dsf[c]*wts[j];
3606  stress[2] +=pow((cz[j]-v[c][2]),2.0)*f_dsf[c]*wts[j];
3607  stress[3] +=(cx[j]-v[c][0])*(cy[j]-v[c][1])*f_dsf[c]*wts[j];
3608  stress[4] +=(cy[j]-v[c][1])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
3609  stress[5] +=(cx[j]-v[c][0])*(cz[j]-v[c][2])*f_dsf[c]*wts[j];
3610  }
3611  }
3612  }
3613 
3614  }
3615 
3616 
3617  const VectorT3& Af = solidFaceArea[f];
3618  force[f][0] = Af[0]*Ly*Lz*stress[0] + Af[1]*Lz*Lx*stress[3] + Af[2]*Lx*Ly*stress[5];
3619  force[f][1] = Af[0]*Ly*Lz*stress[3] + Af[1]*Lz*Lx*stress[1] + Af[2]*Lx*Ly*stress[4];
3620  force[f][2] = Af[0]*Ly*Lz*stress[5] + Af[1]*Lz*Lx*stress[4] + Af[2]*Ly*Ly*stress[2];
3621  if (perUnitArea){
3622  force[f] /= solidFaceAreaMag[f];}
3623  }
3624  }
3625  }
3626 
3627  void OutputDsfPOINT() //, const char* filename)
3628  {
3629  FILE * pFile;
3630  pFile = fopen("cxyz0.plt","w");
3631  int N1=_quadrature.getNVCount();
3632  int N2=_quadrature.getNthetaCount();
3633  int N3=_quadrature.getNphiCount();
3634  fprintf(pFile,"%s \n", "VARIABLES= cx, cy, cz, f,");
3635  fprintf(pFile, "%s %i %s %i %s %i \n","ZONE I=", N3,",J=",N2,"K=",N1);
3636  fprintf(pFile,"%s\n","F=POINT");
3637  const int numMeshes = _meshes.size();
3638  const int cellno=_options.printCellNumber;
3639  for (int n=0; n<numMeshes; n++)
3640  {
3641  const TArray& cx = dynamic_cast<const TArray&>(*_quadrature.cxPtr);
3642  const TArray& cy = dynamic_cast<const TArray&>(*_quadrature.cyPtr);
3643  const TArray& cz = dynamic_cast<const TArray&>(*_quadrature.czPtr);
3644  const int numFields= _quadrature.getDirCount();
3645 
3646  const Mesh& mesh = *_meshes[n];
3647  const StorageSite& cells = mesh.getCells();
3648  for(int j=0;j< numFields;j++)
3649  {
3650  Field& fnd = *_dsfPtr.dsf[j];
3651  TArray& f = dynamic_cast< TArray&>(fnd[cells]);
3652  Field& fEqnd = *_dsfEqPtr.dsf[j];
3653  TArray& fEq = dynamic_cast< TArray&>(fEqnd[cells]);
3654  fprintf(pFile,"%E %E %E %E %E\n",cx[j],cy[j],cz[j],f[cellno],fEq[cellno]);
3655  }
3656  }
3657  }
3658 
3659  const DistFunctFields<T>& getdsf() const { return _dsfPtr;}
3660  const DistFunctFields<T>& getdsf1() const { return _dsfPtr1;}
3661  const DistFunctFields<T>& getdsf2() const { return _dsfPtr2;}
3662  const DistFunctFields<T>& getdsfEq() const { return _dsfEqPtr;}
3663  const DistFunctFields<T>& getdsfEqES() const { return _dsfEqPtrES;}
3664 
3665 
3666  private:
3667  //shared_ptr<Impl> _impl;
3668 
3669  const GeomFields& _geomFields;
3670  const Quadrature<T>& _quadrature;
3671 
3672  MacroFields& _macroFields;
3673  DistFunctFields<T> _dsfPtr;
3674  DistFunctFields<T> _dsfPtr1;
3675  DistFunctFields<T> _dsfPtr2;
3676  DistFunctFields<T> _dsfEqPtr;
3677  DistFunctFields<T> _dsfEqPtrES;
3678 
3679  KineticBCMap _bcMap;
3680  KineticVCMap _vcMap;
3681 
3682  KineticModelOptions<T> _options;
3683 
3684 
3685  MFRPtr _initialKmodelNorm;
3686  //MFRPtr _initialKmodelvNorm;
3687  int _niters;
3688  map<int, vector<int> > _faceReflectionArrayMap;
3689  map<string,shared_ptr<ArrayBase> > _persistenceData;
3690 };
3691 
3692 #endif
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
Array::zero
virtual void zero()
Definition: Array.h:281
KineticModel::getdsf1
const DistFunctFields< T > & getdsf1() const
Definition: KineticModel.h:3660
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
LinearSystem::initAssembly
void initAssembly()
Definition: LinearSystem.cpp:32
KineticModel::getdsfEq
const DistFunctFields< T > & getdsfEq() const
Definition: KineticModel.h:3662
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
KineticBoundaryConditions::applyPressureOutletBC
void applyPressureOutletBC(int f, const X &outletTemperature, const X &outletPressure) const
Definition: KineticBoundaryConditions.h:611
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
KineticModel::getVCMap
KineticVCMap & getVCMap()
Definition: KineticModel.h:1354
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
KineticModel::MomentHierarchy
void MomentHierarchy()
Definition: KineticModel.h:684
KineticModel::_geomFields
const GeomFields & _geomFields
Definition: KineticModel.h:3669
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
KineticModel::updateTime
void updateTime()
Definition: KineticModel.h:3229
KineticModel::_faceReflectionArrayMap
map< int, vector< int > > _faceReflectionArrayMap
Definition: KineticModel.h:3688
KineticModel::ConservationofMFSolid
void ConservationofMFSolid(const StorageSite &solidFaces, const int output=0, bool perUnitArea=0) const
Definition: KineticModel.h:2992
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
KineticBoundaryConditions::applyInletBC
void applyInletBC(int f, const VectorX3 &InletVelocity, const X &InletTemperature, const X &InletMdot, const vector< int > &vecReflection) const
Definition: KineticBoundaryConditions.h:522
KineticModel::KineticBCMap
std::map< int, KineticBC< T > * > KineticBCMap
Definition: KineticModel.h:73
KineticModel::_macroFields
MacroFields & _macroFields
Definition: KineticModel.h:3672
KineticModel::linearizeKineticModel
void linearizeKineticModel(LinearSystem &ls, int direction)
Definition: KineticModel.h:1505
BaseGenericKineticBCS::applyInterfaceBC
void applyInterfaceBC(const int f) const
Definition: GenericKineticBCS.h:175
KineticModel::initializeMaxwellianEq
void initializeMaxwellianEq()
Definition: KineticModel.h:779
FaceGroup
Definition: Mesh.h:28
KineticModel::getPersistenceData
map< string, shared_ptr< ArrayBase > > & getPersistenceData()
Definition: KineticModel.h:1363
KineticModel::correctMassDeficit2
void correctMassDeficit2(double n1, double n2)
Definition: KineticModel.h:2895
KineticModel::ComputeMacroparameters
void ComputeMacroparameters()
Definition: KineticModel.h:471
KineticModel::VectorT10
Vector< T, 10 > VectorT10
Definition: KineticModel.h:69
Field
Definition: Field.h:14
KineticModel::_bcMap
KineticBCMap _bcMap
Definition: KineticModel.h:3679
MacroFields::Stress
Field Stress
Definition: MacroFields.h:37
Mesh::getIBFaceList
const Array< int > & getIBFaceList() const
Definition: Mesh.cpp:686
KineticBoundaryConditions::applyDiffuseWallBC
void applyDiffuseWallBC(int f, const VectorX3 &WallVelocity, const X &WallTemperature) const
Definition: KineticBoundaryConditions.h:65
MacroFields::EntropyGenRate_Collisional
Field EntropyGenRate_Collisional
Definition: MacroFields.h:35
KineticModel::getBCMap
KineticBCMap & getBCMap()
Definition: KineticModel.h:1353
Mesh
Definition: Mesh.h:49
KineticModel::_persistenceData
map< string, shared_ptr< ArrayBase > > _persistenceData
Definition: KineticModel.h:3689
KineticModel::correctMassDeficit
void correctMassDeficit()
Definition: KineticModel.h:2774
MacroFields::InterfaceVelocity
Field InterfaceVelocity
Definition: MacroFields.h:39
KineticModel::_quadrature
const Quadrature< T > & _quadrature
Definition: KineticModel.h:3670
MacroFields::EntropyGenRate
Field EntropyGenRate
Definition: MacroFields.h:34
KineticBoundaryConditions::applyPressureInletBC
void applyPressureInletBC(int f, const X &inletTemperature, const X &inletPressure) const
Definition: KineticBoundaryConditions.h:408
KineticModel::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: KineticModel.h:59
MacroFields::Entropy
Field Entropy
Definition: MacroFields.h:33
MacroFields::temperatureIB
Field temperatureIB
Definition: MacroFields.h:23
TimeDerivativeDiscretization_Kmodel
Definition: TimeDerivativeDiscretization_Kmodel.h:21
KineticBoundaryConditions::applyRealWallBC
void applyRealWallBC(int f, const VectorX3 &WallVelocity, const X &WallTemperature, const X &accommodationCoefficient, const vector< int > &vecReflection) const
Definition: KineticBoundaryConditions.h:158
KineticModel::getOptions
KineticModelOptions< T > & getOptions()
Definition: KineticModel.h:1356
MacroFields::InterfacePressure
Field InterfacePressure
Definition: MacroFields.h:40
KineticModel::getdsf
const DistFunctFields< T > & getdsf() const
Definition: KineticModel.h:3659
KineticModel::MacroparameterIBCell
void MacroparameterIBCell(const StorageSite &solidFaces) const
Definition: KineticModel.h:3113
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
FloatVarDict::getVal
FloatVal< T > getVal(const string varName) const
Definition: FloatVarDict.h:85
KineticModel::_dsfPtr2
DistFunctFields< T > _dsfPtr2
Definition: KineticModel.h:3675
MacroFields::viscosity
Field viscosity
Definition: MacroFields.h:20
KineticModel::getdsf2
const DistFunctFields< T > & getdsf2() const
Definition: KineticModel.h:3661
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
KineticModel::getFaceReflectionArrayMap
const map< int, vector< int > > & getFaceReflectionArrayMap() const
Definition: KineticModel.h:1357
FaceGroup::id
const int id
Definition: Mesh.h:41
LinearSystem::updateSolution
void updateSolution()
Definition: LinearSystem.cpp:265
KineticModel::OutputDsfBLOCK
void OutputDsfBLOCK(const char *filename)
Definition: KineticModel.h:3470
KineticModel::stdVectorField
std::vector< Field * > stdVectorField
Definition: KineticModel.h:62
KineticModel::_dsfPtr1
DistFunctFields< T > _dsfPtr1
Definition: KineticModel.h:3674
KineticBoundaryConditions::applySpecularWallBC
void applySpecularWallBC(int f, const vector< int > &vecReflection) const
Definition: KineticBoundaryConditions.h:287
MacroFields::force
Field force
Definition: MacroFields.h:36
KineticModel::EquilibriumDistributionBGK
void EquilibriumDistributionBGK()
Definition: KineticModel.h:907
Array::getData
virtual void * getData() const
Definition: Array.h:275
KineticModel::VectorT6
Vector< T, 6 > VectorT6
Definition: KineticModel.h:67
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
inverseGauss
SquareMatrix< T, N > inverseGauss(SquareMatrix< T, N > &a, int size)
Definition: MatrixOperation.h:201
LinearSystem::initSolve
void initSolve()
Definition: LinearSystem.cpp:40
KineticModel::_vcMap
KineticVCMap _vcMap
Definition: KineticModel.h:3680
KineticBC::bcType
string bcType
Definition: KineticBC.h:32
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
MacroFields::coeff
Field coeff
Definition: MacroFields.h:31
KineticModel::VectorT5Array
Array< VectorT5 > VectorT5Array
Definition: KineticModel.h:66
KineticModel::_dsfPtr
DistFunctFields< T > _dsfPtr
Definition: KineticModel.h:3673
DiscrList
vector< shared_ptr< Discretization > > DiscrList
Definition: Discretization.h:31
KineticModel::computeSolidFaceDsf
void computeSolidFaceDsf(const StorageSite &solidFaces, const int method, const int RelaxDistribution=0)
Definition: KineticModel.h:2347
KineticModel::TDistFF
DistFunctFields< T > TDistFF
Definition: KineticModel.h:63
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCellCells
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
Model
Definition: Model.h:13
GeomFields::ibFaceIndex
Field ibFaceIndex
Definition: GeomFields.h:40
KineticModel::KineticModel
KineticModel(const MeshList &meshes, const GeomFields &geomFields, MacroFields &macroFields, const Quadrature< T > &quad)
Definition: KineticModel.h:81
KineticModel::_dsfEqPtrES
DistFunctFields< T > _dsfEqPtrES
Definition: KineticModel.h:3677
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
ConvectionDiscretization_Kmodel.h
GeomFields::volume
Field volume
Definition: GeomFields.h:26
KineticModel::EntropyGeneration
void EntropyGeneration()
Definition: KineticModel.h:726
KineticModel::NewtonsMethodBGK
void NewtonsMethodBGK(const int ktrial)
Definition: KineticModel.h:829
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
KineticModel::InitializeMacroparameters
void InitializeMacroparameters()
Definition: KineticModel.h:349
KineticModel::getdsfEqES
const DistFunctFields< T > & getdsfEqES() const
Definition: KineticModel.h:3663
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
KineticModel::_initialKmodelNorm
MFRPtr _initialKmodelNorm
Definition: KineticModel.h:3685
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
KineticModel::ComputeMacroparametersESBGK
void ComputeMacroparametersESBGK()
Definition: KineticModel.h:565
KineticModel::VectorT5
Vector< T, 5 > VectorT5
Definition: KineticModel.h:65
KineticModel::setJacobianBGK
void setJacobianBGK(SquareMatrix< T, 5 > &fjac, VectorT5 &fvec, const VectorT5 &xn, const VectorT3 &v, const int c)
Definition: KineticModel.h:958
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
KineticModel::VectorT6Array
Array< VectorT6 > VectorT6Array
Definition: KineticModel.h:68
MacroFields::velocity
Field velocity
Definition: MacroFields.h:15
MultiFieldMatrix::addMatrix
void addMatrix(const Index &rowI, const Index &colI, shared_ptr< Matrix > m)
Definition: MultiFieldMatrix.h:137
KineticModel::OutputDsfPOINT
void OutputDsfPOINT()
Definition: KineticModel.h:3627
KineticModel::advance
bool advance(const int niter, const int updated=0)
Definition: KineticModel.h:3360
MacroFields::pressure
Field pressure
Definition: MacroFields.h:19
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
KineticModel::IntArray
Array< int > IntArray
Definition: KineticModel.h:55
MacroFields::InterfaceStress
Field InterfaceStress
Definition: MacroFields.h:41
KineticModel::ConservationofMassCheck
const double ConservationofMassCheck()
Definition: KineticModel.h:2954
KineticModel::setJacobianESBGK
void setJacobianESBGK(SquareMatrix< T, 10 > &fjac, VectorT10 &fvec, const VectorT10 &xn, const VectorT3 &v, const int c)
Definition: KineticModel.h:1167
KineticModel::computeIBFaceDsf
void computeIBFaceDsf(const StorageSite &solidFaces, const int method, const int RelaxDistribution=0)
Definition: KineticModel.h:1701
KineticModel::VectorT3
Vector< T, 3 > VectorT3
Definition: KineticModel.h:58
KineticModel::computeSurfaceForce
void computeSurfaceForce(const StorageSite &solidFaces, bool perUnitArea, bool IBM=0)
Definition: KineticModel.h:3523
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
LinearSystem::getB
MultiField & getB()
Definition: LinearSystem.h:33
KineticModel::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: KineticModel.h:54
KineticModel::SetBoundaryConditions
void SetBoundaryConditions()
Definition: KineticModel.h:1408
KineticModel::TArray2D
Array2D< T > TArray2D
Definition: KineticModel.h:57
KineticModel::NewtonsMethodESBGK
void NewtonsMethodESBGK(const int ktrial)
Definition: KineticModel.h:999
KineticModel::StressTensorArray
Array< StressTensorT6 > StressTensorArray
Definition: KineticModel.h:61
MFRPtr
shared_ptr< MultiFieldReduction > MFRPtr
Definition: MultiFieldReduction.h:56
GeomFields::area
Field area
Definition: GeomFields.h:23
MultiField::addArray
void addArray(const ArrayIndex &aIndex, shared_ptr< ArrayBase > a)
Definition: MultiField.cpp:270
KineticModel::_dsfEqPtr
DistFunctFields< T > _dsfEqPtr
Definition: KineticModel.h:3676
ConvectionDiscretization_Kmodel
Definition: ConvectionDiscretization_Kmodel.h:18
Mesh::getInterfaceGroups
const FaceGroupList & getInterfaceGroups() const
Definition: Mesh.h:190
LinearSystem::getX
MultiField & getX()
Definition: LinearSystem.h:32
KineticModel::initKineticModelLinearization
void initKineticModelLinearization(LinearSystem &ls, int direction)
Definition: KineticModel.h:1480
KineticModel::weightedMaxwellian
void weightedMaxwellian(double weight1, double uvel1, double uvel2, double temp1, double temp2)
Definition: KineticModel.h:1309
KineticModel::StressTensorT6
StressTensor< T > StressTensorT6
Definition: KineticModel.h:60
KineticModel::initializeMaxwellian
void initializeMaxwellian()
Definition: KineticModel.h:1216
MacroFields::collisionFrequency
Field collisionFrequency
Definition: MacroFields.h:24
KineticModel::VectorT10Array
Array< VectorT10 > VectorT10Array
Definition: KineticModel.h:70
Field::getArrayPtr
shared_ptr< ArrayBase > getArrayPtr(const StorageSite &)
Definition: Field.cpp:63
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
KineticModel::ComputeCollisionfrequency
void ComputeCollisionfrequency()
Definition: KineticModel.h:645
Mesh::getID
int getID() const
Definition: Mesh.h:106
KineticModel::TArray
Array< T > TArray
Definition: KineticModel.h:56
TimeDerivativeDiscretization_Kmodel.h
KineticModel::KineticVCMap
std::map< int, KineticVC< T > * > KineticVCMap
Definition: KineticModel.h:74
Linearizer::linearize
virtual void linearize(DiscrList &discretizations, const MeshList &meshes, MultiFieldMatrix &matrix, MultiField &x, MultiField &b)
Definition: Linearizer.cpp:17
MacroFields::density
Field density
Definition: MacroFields.h:21
KineticModel::callBoundaryConditions
void callBoundaryConditions()
Definition: KineticModel.h:3271
KineticModel::_options
KineticModelOptions< T > _options
Definition: KineticModel.h:3682
KineticModel::EquilibriumDistributionESBGK
void EquilibriumDistributionESBGK()
Definition: KineticModel.h:1112
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
LinearSystem::getMatrix
MultiFieldMatrix & getMatrix()
Definition: LinearSystem.h:37
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
KineticModel::weightedMaxwellian
void weightedMaxwellian(double weight1, double uvel1, double vvel1, double wvel1, double uvel2, double vvel2, double wvel2, double temp1, double temp2)
Definition: KineticModel.h:1269
KineticModel::InitializeFgammaCoefficients
void InitializeFgammaCoefficients()
Definition: KineticModel.h:417
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
MacroFields::InterfaceDensity
Field InterfaceDensity
Definition: MacroFields.h:42
KineticVC::vcType
string vcType
Definition: KineticBC.h:43
KineticModel::computeSolidFacePressure
void computeSolidFacePressure(const StorageSite &solidFaces)
Definition: KineticModel.h:2310