Memosa-FVM  0.2
COMETDiscretizer.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 _COMETDISCRETIZER_H_
6 #define _COMETDISCRETIZER_H_
7 
8 #include "Mesh.h"
9 #include "Kspace.h"
10 #include "kvol.h"
11 #include "pmode.h"
12 #include "ArrowHeadMatrix.h"
13 #include "Array.h"
14 #include "Vector.h"
15 #include "StorageSite.h"
16 #include "CRConnectivity.h"
17 #include "GeomFields.h"
18 #include "NumType.h"
19 #include "COMETBC.h"
20 #include <math.h>
21 #include <map>
22 #include "MatrixJML.h"
23 #include "SquareMatrix.h"
24 #include "PhononMacro.h"
25 #include "SquareTensor.h"
26 #include "GradientModel.h"
27 #include "FluxLimiters.h"
28 #include <omp.h>
29 
30 template<class T>
31 class COMETDiscretizer
32 {
33 
34  public:
35  typedef typename NumTypeTraits<T>::T_Scalar T_Scalar;
36  typedef Kspace<T> Tkspace;
37  typedef kvol<T> Tkvol;
38  typedef pmode<T> Tmode;
39  typedef Vector<T_Scalar,3> VectorT3;
40  typedef Array<VectorT3> VectorT3Array;
41  typedef Array<T_Scalar> TArray;
42  typedef MatrixJML<T> TMatrix;
43  typedef ArrowHeadMatrix<T> TArrow;
44  typedef SquareMatrix<T> TSquare;
45  typedef map<int,COMETBC<T>*> COMETBCMap;
46  typedef COMETModelOptions<T> COpts;
47  typedef Array<int> IntArray;
48  typedef Array<bool> BoolArray;
49  typedef Vector<int,2> VecInt2;
50  typedef map<int,VecInt2> FaceToFg;
51  typedef typename Tmode::Refl_pair Refl_pair;
52  typedef SquareTensor<T,3> T3Tensor;
53  typedef KSConnectivity<T> TKConnectivity;
54  typedef TKConnectivity* TKCptr;
55  typedef vector<TKCptr> TKClist;
56  typedef map<int, TKClist> FgTKClistMap;
57  typedef Gradient<T> GradType;
58  typedef Array<GradType> GradArray;
59  typedef GradientModel<T> GradModelType;
60  typedef typename GradModelType::GradMatrixType GradMatrix;
61 
62  COMETDiscretizer(const Mesh& mesh, const GeomFields& geomfields,
63  PhononMacro& macro, Tkspace& kspace, COMETBCMap& bcMap,
64  const IntArray& BCArray, const IntArray& BCfArray, COpts& options,
65  const FgTKClistMap& FgToKsc):
66  _mesh(mesh),
67  _geomFields(geomfields),
68  _cells(mesh.getCells()),
69  _faces(mesh.getFaces()),
70  _cellFaces(mesh.getCellFaces()),
71  _faceCells(mesh.getAllFaceCells()),
72  _faceAreaMag((dynamic_cast<const TArray&>(_geomFields.areaMag[_faces]))),
73  _faceArea(dynamic_cast<const VectorT3Array&>(_geomFields.area[_faces])),
74  _cellVolume(dynamic_cast<const TArray&>(_geomFields.volume[_cells])),
75  _cellCoords(dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_cells])),
76  _macro(macro),
77  _kspace(kspace),
78  _bcMap(bcMap),
79  _BCArray(BCArray),
80  _BCfArray(BCfArray),
81  _aveResid(-1.),
82  _residChange(-1.),
83  _fgFinder(),
84  _options(options),
85  _FaceToKSC(FgToKsc),
86  _eArray(kspace.geteArray()),
87  _e0Array(kspace.gete0Array()),
88  _resArray(kspace.getResArray())
89  {}
90 
91  void COMETSolveFine(const int dir,const int level)
92  {
93  const int cellcount=_cells.getSelfCount();
94  int start;
95 
96  if(dir==1)
97  start=0;
98  if(dir==-1)
99  start=cellcount-1;
100  const int totalmodes=_kspace.gettotmodes();
101  TArray Bvec(totalmodes+1);
102  TArray Resid(totalmodes+1);
103  TArrow AMat(totalmodes+1);
104  const T newTol=_options.NewtonTol;
105  const int maxNew=_options.maxNewton;
106  const int minNew=_options.minNewton;
107 
108  const GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
109 
110  for(int c=start;((c<cellcount)&&(c>-1));c+=dir)
111  {
112 
113  if(_BCArray[c]==0) //no reflections at all--interior cell or temperature boundary
114  {
115  T dt=1;
116  int NewtonIters=0;
117 
118  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
119  {
120 
121  Bvec.zero();
122  Resid.zero();
123  AMat.zero();
124 
125  updateGhostFine(c, gradMatrix);
126  COMETConvectionFine(c,AMat,Bvec,gradMatrix);
127  COMETCollision(c,&AMat,Bvec);
128  COMETEquilibrium(c,&AMat,Bvec);
129 
130  if(_options.withNormal)
131  COMETShifted(c,&AMat,Bvec);
132 
133  if(level>0)
134  addFAS(c,Bvec);
135 
136  Resid=Bvec;
137  AMat.Solve(Bvec);
138 
139  Distribute(c,Bvec,Resid);
140  updatee0(c);
141  NewtonIters++;
142  dt=fabs(Bvec[totalmodes]);
143  }
144  }
145  else if(_BCArray[c]==1) //Implicit reflecting boundary only
146  {
147  T dt=1;
148  int NewtonIters=0;
149  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
150  {
151  TSquare AMatS(totalmodes+1);
152  Bvec.zero();
153  Resid.zero();
154  AMatS.zero();
155 
156  COMETConvection(c,AMatS,Bvec);
157  COMETCollision(c,&AMatS,Bvec);
158  COMETEquilibrium(c,&AMatS,Bvec);
159 
160  if(level>0)
161  addFAS(c,Bvec);
162 
163  Resid=Bvec;
164  AMatS.Solve(Bvec);
165 
166  Distribute(c,Bvec,Resid);
167  updatee0(c);
168  NewtonIters++;
169  dt=fabs(Bvec[totalmodes]);
170  }
171  }
172  else if(_BCArray[c]==2) //Explicit boundary only
173  {
174  T dt=1;
175  int NewtonIters=0;
176  updateGhostFine(c, gradMatrix);
177  //updateGhostCoarse(c);
178 
179  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
180  {
181 
182  Bvec.zero();
183  Resid.zero();
184  AMat.zero();
185 
186  COMETConvectionFine(c,AMat,Bvec,gradMatrix);
187  COMETCollision(c,&AMat,Bvec);
188  COMETEquilibrium(c,&AMat,Bvec);
189 
190  if(_options.withNormal)
191  COMETShifted(c,&AMat,Bvec);
192 
193  if(level>0)
194  addFAS(c,Bvec);
195 
196  Resid=Bvec;
197  AMat.Solve(Bvec);
198 
199  Distribute(c,Bvec,Resid);
200  dt=fabs(Bvec[totalmodes]);
201  updatee0(c);
202  updateGhostFine(c, gradMatrix);
203  if(!_FaceToKSC.empty())
204  correctInterface(c,Bvec);
205  NewtonIters++;
206  }
207  }
208  else if(_BCArray[c]==3) //Mix Implicit/Explicit reflecting boundary
209  {
210  T dt=1;
211  int NewtonIters=0;
212  updateGhostFine(c,gradMatrix);
213  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
214  {
215  TSquare AMatS(totalmodes+1);
216  Bvec.zero();
217  Resid.zero();
218  AMatS.zero();
219 
220  COMETConvection(c,AMatS,Bvec);
221  COMETCollision(c,&AMatS,Bvec);
222  COMETEquilibrium(c,&AMatS,Bvec);
223 
224  if(level>0)
225  addFAS(c,Bvec);
226 
227  Resid=Bvec;
228  AMatS.Solve(Bvec);
229 
230  Distribute(c,Bvec,Resid);
231  updatee0(c);
232  updateGhostFine(c, gradMatrix);
233  NewtonIters++;
234  dt=fabs(Bvec[totalmodes]);
235  }
236  }
237  else
238  throw CException("Unexpected value for boundary cell map.");
239 
240  }
241 
242  }
243 
244  void COMETSolveCoarse(const int dir,const int level)
245  {
246  const int cellcount=_cells.getSelfCount();
247  int start;
248 
249  if(dir==1)
250  start=0;
251  if(dir==-1)
252  start=cellcount-1;
253  const int totalmodes=_kspace.gettotmodes();
254  TArray Bvec(totalmodes+1);
255  TArray Resid(totalmodes+1);
256  TArrow AMat(totalmodes+1);
257  const T newTol=_options.NewtonTol;
258  const int maxNew=_options.maxNewton;
259  const int minNew=_options.minNewton;
260  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
261 
262  for(int c=start;((c<cellcount)&&(c>-1));c+=dir)
263  {
264 
265  if(_BCArray[c]==0) //no reflections at all--interior cell or temperature boundary
266  {
267  T dt=1;
268  T rscaled=10;
269  int NewtonIters=0;
270 
271  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
272  {
273 
274  Bvec.zero();
275  Resid.zero();
276  AMat.zero();
277 
278  COMETConvectionCoarse(c,AMat,Bvec);
279  COMETCollision(c,&AMat,Bvec);
280  COMETEquilibrium(c,&AMat,Bvec);
281  COMETSource(c,Bvec);
282 
283  if(_options.withNormal)
284  COMETShifted(c,&AMat,Bvec);
285 
286  if(level>0)
287  addFAS(c,Bvec);
288 
289  Resid=Bvec;
290  AMat.Solve(Bvec);
291  rscaled=scaledResid(Bvec,c);
292 
293  Distribute(c,Bvec,Resid);
294  updatee0(c);
295  if(!_FaceToKSC.empty())
296  correctInterface(c,Bvec);
297  NewtonIters++;
298  dt=0.5*fabs(Bvec[totalmodes]/Tl[c])+rscaled*0.5;
299  }
300 
301  }
302  else if(_BCArray[c]==1) //Implicit reflecting boundary only
303  {
304  T dt=1;
305  int NewtonIters=0;
306  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
307  {
308  TSquare AMatS(totalmodes+1);
309  Bvec.zero();
310  Resid.zero();
311  AMatS.zero();
312 
313  COMETConvection(c,AMatS,Bvec);
314  COMETCollision(c,&AMatS,Bvec);
315  COMETEquilibrium(c,&AMatS,Bvec);
316 
317  if(level>0)
318  addFAS(c,Bvec);
319 
320  Resid=Bvec;
321  AMatS.Solve(Bvec);
322 
323  Distribute(c,Bvec,Resid);
324  updatee0(c);
325  NewtonIters++;
326  dt=fabs(Bvec[totalmodes]);
327  }
328  }
329  else if(_BCArray[c]==2) //Explicit boundary only
330  {
331  T dt=1;
332  T rscaled=10;
333  int NewtonIters=0;
334  updateGhostCoarse(c);
335 
336  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
337  {
338 
339  Bvec.zero();
340  Resid.zero();
341  AMat.zero();
342 
343  COMETConvectionCoarse(c,AMat,Bvec);
344  COMETCollision(c,&AMat,Bvec);
345  COMETEquilibrium(c,&AMat,Bvec);
346  COMETSource(c,Bvec);
347 
348  if(_options.withNormal)
349  COMETShifted(c,&AMat,Bvec);
350 
351  if(level>0)
352  addFAS(c,Bvec);
353 
354  Resid=Bvec;
355  AMat.Solve(Bvec);
356  rscaled=scaledResid(Bvec,c);
357 
358  Distribute(c,Bvec,Resid);
359  dt=0.5*fabs(Bvec[totalmodes]/Tl[c])+0.5*rscaled;
360  updatee0(c);
361  updateGhostCoarse(c);
362  if(!_FaceToKSC.empty())
363  correctInterface(c,Bvec);
364  NewtonIters++;
365  }
366  }
367  else if(_BCArray[c]==3) //Mix Implicit/Explicit reflecting boundary
368  {
369  T dt=1;
370  int NewtonIters=0;
371  updateGhostCoarse(c);
372  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
373  {
374  TSquare AMatS(totalmodes+1);
375  Bvec.zero();
376  Resid.zero();
377  AMatS.zero();
378 
379  COMETConvection(c,AMatS,Bvec);
380  COMETCollision(c,&AMatS,Bvec);
381  COMETEquilibrium(c,&AMatS,Bvec);
382 
383  if(level>0)
384  addFAS(c,Bvec);
385 
386  Resid=Bvec;
387  AMatS.Solve(Bvec);
388 
389  Distribute(c,Bvec,Resid);
390  updatee0(c);
391  updateGhostCoarse(c);
392  NewtonIters++;
393  dt=fabs(Bvec[totalmodes]);
394  }
395  }
396  else
397  throw CException("Unexpected value for boundary cell map.");
398 
399  }
400 
401  }
402 
403  void COMETSolveFull(const int dir,const int level)
404  {
405  const int cellcount=_cells.getSelfCount();
406  int start;
407 
408  if(dir==1)
409  start=0;
410  if(dir==-1)
411  start=cellcount-1;
412  const int totalmodes=_kspace.gettotmodes();
413  TArray Bvec(totalmodes+1);
414  TArray Resid(totalmodes+1);
415  TArray s(totalmodes);
416  TArray ds(totalmodes);
417  TArrow AMat(totalmodes+1);
418  const T newTol=_options.NewtonTol;
419  const int maxNew=_options.maxNewton;
420  const int minNew=_options.minNewton;
421 
422  for(int c=start;((c<cellcount)&&(c>-1));c+=dir)
423  {
424 
425  if(_BCArray[c]==0) //no reflections at all--interior cell or temperature boundary
426  {
427  T dt=1;
428  int NewtonIters=0;
429 
430  for(int srcIts=0;srcIts<1;srcIts++)
431  {
432 
433  s.zero();
434  //_kspace.getSourceTerm(c,s,ds);
435 
436  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
437  {
438 
439  Bvec.zero();
440  Resid.zero();
441  AMat.zero();
442 
443  COMETConvectionCoarse(c,AMat,Bvec);
444  COMETCollision(c,&AMat,Bvec);
445  COMETEquilibrium(c,&AMat,Bvec);
446  //COMETFullScatt(c,s,Bvec);
447 
448  if(_options.withNormal)
449  COMETShifted(c,&AMat,Bvec);
450 
451  if(level>0)
452  addFAS(c,Bvec);
453 
454  Resid=Bvec;
455  AMat.Solve(Bvec);
456  //AMat.SolveDiag(Bvec);
457  //AMat.smoothGS(Bvec);
458 
459  Distribute(c,Bvec,Resid);
460  updatee0(c);
461  NewtonIters++;
462  dt=fabs(Bvec[totalmodes]);
463  }
464 
465  ScatterPhonons(c);
466 
467  }
468  }
469  else if(_BCArray[c]==1) //Implicit reflecting boundary only
470  {
471  T dt=1;
472  int NewtonIters=0;
473  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
474  {
475  TSquare AMatS(totalmodes+1);
476  Bvec.zero();
477  Resid.zero();
478  AMatS.zero();
479 
480  COMETConvection(c,AMatS,Bvec);
481  COMETCollision(c,&AMatS,Bvec);
482  COMETEquilibrium(c,&AMatS,Bvec);
483 
484  if(level>0)
485  addFAS(c,Bvec);
486 
487  Resid=Bvec;
488  AMatS.Solve(Bvec);
489 
490  Distribute(c,Bvec,Resid);
491  updatee0(c);
492  NewtonIters++;
493  dt=fabs(Bvec[totalmodes]);
494  }
495  }
496  else if(_BCArray[c]==2) //Explicit boundary only
497  {
498  T dt=1;
499  int NewtonIters=0;
500  updateGhostCoarse(c);
501 
502  for(int srcIts=0;srcIts<1;srcIts++)
503  {
504  s.zero();
505  //_kspace.getSourceTerm(c,s,ds);
506 
507  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
508  {
509 
510  Bvec.zero();
511  Resid.zero();
512  AMat.zero();
513 
514  COMETConvectionCoarse(c,AMat,Bvec);
515  COMETCollision(c,&AMat,Bvec);
516  COMETEquilibrium(c,&AMat,Bvec);
517  //COMETFullScatt(c,s,Bvec);
518 
519  if(_options.withNormal)
520  COMETShifted(c,&AMat,Bvec);
521 
522  if(level>0)
523  addFAS(c,Bvec);
524 
525  Resid=Bvec;
526  AMat.Solve(Bvec);
527  //AMat.SolveDiag(Bvec);
528  //AMat.smoothGS(Bvec);
529 
530  Distribute(c,Bvec,Resid);
531  dt=fabs(Bvec[totalmodes]);
532  updatee0(c);
533  updateGhostCoarse(c);
534  if(!_FaceToKSC.empty())
535  correctInterface(c,Bvec);
536  NewtonIters++;
537  }
538 
539  ScatterPhonons(c);
540 
541  }
542  }
543  else if(_BCArray[c]==3) //Mix Implicit/Explicit reflecting boundary
544  {
545  T dt=1;
546  int NewtonIters=0;
547  updateGhostCoarse(c);
548  while((dt>newTol && NewtonIters<maxNew) || NewtonIters<minNew)
549  {
550  TSquare AMatS(totalmodes+1);
551  Bvec.zero();
552  Resid.zero();
553  AMatS.zero();
554 
555  COMETConvection(c,AMatS,Bvec);
556  COMETCollision(c,&AMatS,Bvec);
557  COMETEquilibrium(c,&AMatS,Bvec);
558 
559  if(level>0)
560  addFAS(c,Bvec);
561 
562  Resid=Bvec;
563  AMatS.Solve(Bvec);
564 
565  Distribute(c,Bvec,Resid);
566  updatee0(c);
567  updateGhostCoarse(c);
568  NewtonIters++;
569  dt=fabs(Bvec[totalmodes]);
570  }
571  }
572  else
573  throw CException("Unexpected value for boundary cell map.");
574 
575  }
576 
577  }
578 
579  void COMETConvectionFine(const int cell0, TArrow& Amat,
580  TArray& BVec, const GradMatrix& gMat)
581  {
582 
583  const int neibcount=_cellFaces.getCount(cell0);
584  GradArray Grads(_kspace.gettotmodes());
585  TArray pointMin(_kspace.gettotmodes());
586  pointMin=-1;
587  TArray pointMax(_kspace.gettotmodes());
588  pointMax.zero();
589  TArray neibMin(_kspace.gettotmodes());
590  neibMin=-1;
591  TArray neibMax(_kspace.gettotmodes());
592  neibMax.zero();
593  TArray pointLim(_kspace.gettotmodes());
594  pointLim=100.;
595  TArray neibLim(_kspace.gettotmodes());
596  neibLim=100.;
597 
598  for(int k=0;k<_kspace.gettotmodes();k++)
599  Grads[k].zero();
600 
601  const VectorT3Array& faceCoords=
602  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_faces]);
603 
604  for(int j=0;j<neibcount;j++) //first loop to get grad and max/min vals
605  {
606  const int f=_cellFaces(cell0,j);
607  int cell1=_faceCells(f,1);
608  if(cell1==cell0)
609  cell1=_faceCells(f,0);
610  const VectorT3 Gcoeff=gMat.getCoeff(cell0, cell1);
611 
612  int c0ind=_kspace.getGlobalIndex(cell0,0);
613  int c1ind=_kspace.getGlobalIndex(cell1,0);
614 
615  for(int k=0;k<_kspace.gettotmodes();k++)
616  {
617  const T e1=_eArray[c1ind];
618  const T e0=_eArray[c0ind];
619  T& curMax=pointMax[k];
620  T& curMin=pointMin[k];
621  Grads[k].accumulate(Gcoeff, e1-e0);
622 
623  if(e1>curMax)
624  curMax=e1;
625 
626  if(curMin==-1)
627  curMin=e1;
628  else
629  {
630  if(e1<curMin)
631  curMin=e1;
632  }
633 
634  c0ind++;
635  c1ind++;
636  }
637  }
638 
639  for(int j=0;j<neibcount;j++) //second loop to calculate limiting coeff
640  {
641  const int f=_cellFaces(cell0,j);
642  int cell1=_faceCells(f,1);
643  if(cell1==cell0)
644  cell1=_faceCells(f,0);
645 
646  const VectorT3 dr0(faceCoords[f]-_cellCoords[cell0]);
647  int c0ind=_kspace.getGlobalIndex(cell0,0);
648  vanLeer lf;
649 
650  for(int k=0;k<_kspace.gettotmodes();k++)
651  {
652  const T minVal=pointMin[k];
653  const T maxVal=pointMax[k];
654  const T de0(Grads[k]*dr0);
655  T& cl=pointLim[k];
656  computeLimitCoeff2(cl, _eArray[c0ind], de0, minVal, maxVal, lf);
657  c0ind++;
658  }
659  }
660 
661  for(int j=0;j<neibcount;j++) //loop to construct point matrix
662  {
663  const int f=_cellFaces(cell0,j);
664  int cell1=_faceCells(f,1);
665  VectorT3 Af=_faceArea[f];
666  if(cell1==cell0)
667  {
668  cell1=_faceCells(f,0);
669  Af*=-1.;
670  }
671 
672  if(_BCfArray[f]==0) //interior face
673  {
674  const int klen=_kspace.getlength();
675 
676  GradArray NeibGrads(_kspace.gettotmodes());
677 
678  for(int k=0;k<_kspace.gettotmodes();k++)
679  NeibGrads[k].zero();
680 
681  neibMax.zero();
682  neibMin=-1.;
683  const int neibcount1=_cellFaces.getCount(cell1);
684  for(int nj=0;nj<neibcount1;nj++) //loop to make gradients and min/max
685  {
686  const int f1=_cellFaces(cell1,nj);
687  int cell11=_faceCells(f1,1);
688  if(cell1==cell11)
689  cell11=_faceCells(f1,0);
690 
691  const VectorT3 Gcoeff=gMat.getCoeff(cell1, cell11);
692 
693  int c1ind=_kspace.getGlobalIndex(cell1,0);
694  int c11ind=_kspace.getGlobalIndex(cell11,0);
695 
696  for(int k=0;k<_kspace.gettotmodes();k++)
697  {
698  const T e1=_eArray[c1ind];
699  const T e11=_eArray[c11ind];
700  T& curMax=neibMax[k];
701  T& curMin=neibMin[k];
702  NeibGrads[k].accumulate(Gcoeff,e11-e1);
703 
704  if(e11>curMax)
705  curMax=e11;
706 
707  if(curMin==-1)
708  curMin=e11;
709  else
710  {
711  if(e11<curMin)
712  curMin=e11;
713  }
714 
715  c11ind++;
716  c1ind++;
717  }
718  }
719 
720  neibLim=100.;
721  for(int nj=0;nj<neibcount1;nj++) //second loop to calculate limiting coeff
722  {
723  const int f1=_cellFaces(cell1,nj);
724  int cell11=_faceCells(f1,1);
725  if(cell1==cell11)
726  cell11=_faceCells(f1,0);
727 
728  const VectorT3 dr1(faceCoords[f1]-_cellCoords[cell1]);
729  int c1ind=_kspace.getGlobalIndex(cell1,0);
730  vanLeer lf;
731 
732  for(int k=0;k<_kspace.gettotmodes();k++)
733  {
734  const T minVal=neibMin[k];
735  const T maxVal=neibMax[k];
736  const T de1(NeibGrads[k]*dr1);
737  T& cl=neibLim[k];
738  computeLimitCoeff2(cl, _eArray[c1ind], de1, minVal, maxVal, lf);
739  if(_BCfArray[f]!=0)
740  NeibGrads[k].zero();
741  c1ind++;
742  }
743  }
744 
745 
746  T flux;
747 
748  for(int k=0;k<klen;k++)
749  {
750  Tkvol& kvol=_kspace.getkvol(k);
751  const int numModes=kvol.getmodenum();
752  for(int m=0;m<numModes;m++)
753  {
754  Tmode& mode=kvol.getmode(m);
755  const int count=mode.getIndex();
756  VectorT3 vg=mode.getv();
757  //T tau=mode.gettau();
758  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
759 
760  const int c0ind=_kspace.getGlobalIndex(cell0,count-1);
761  const int c1ind=_kspace.getGlobalIndex(cell1,count-1);
762 
763  const GradType& grad=Grads[count-1];
764  const GradType& neibGrad=NeibGrads[count-1];
765  //vanLeer vl;
766 
767  if(flux>T_Scalar(0))
768  {
769  const VectorT3 rVec=_cellCoords[cell1]-_cellCoords[cell0];
770  const VectorT3 fVec=faceCoords[f]-_cellCoords[cell0];
771  //const T r=gMat.computeR(grad,_eArray,rVec,c0ind,c1ind);
772 
773  T SOU=(fVec[0]*grad[0]+fVec[1]*grad[1]+
774  fVec[2]*grad[2])*pointLim[count-1];
775  Amat.getElement(count,count)-=flux;
776  BVec[count-1]-=flux*_eArray[c0ind]-flux*SOU;
777  }
778  else
779  {
780  const VectorT3 rVec=_cellCoords[cell0]-_cellCoords[cell1];
781  const VectorT3 fVec=faceCoords[f]-_cellCoords[cell1];
782  //const T r=gMat.computeR(grad,_eArray,rVec,c1ind,c0ind);
783 
784  T SOU=(fVec[0]*neibGrad[0]+fVec[1]*neibGrad[1]+
785  fVec[2]*neibGrad[2])*neibLim[count-1];
786  BVec[count-1]-=flux*_eArray[c1ind]-flux*SOU;
787  }
788 
789  }
790  }
791  }
792  else
793  {
794  const int klen=_kspace.getlength();
795  T flux(0);
796  for(int k=0;k<klen;k++)
797  {
798 
799  Tkvol& kvol=_kspace.getkvol(k);
800  const int numModes=kvol.getmodenum();
801  for(int m=0;m<numModes;m++)
802  {
803  Tmode& mode=kvol.getmode(m);
804  const int count=mode.getIndex();
805  const VectorT3 vg=mode.getv();
806  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
807 
808  if(flux>T_Scalar(0))
809  Amat.getElement(count,count)-=flux;
810 
811  BVec[count-1]-=flux*_eArray[_kspace.getGlobalIndex(cell1,count-1)];
812  }
813  }
814  }
815 
816  }
817  }
818 
819  void COMETConvectionCoarse(const int cell0, TArrow& Amat, TArray& BVec)
820  {
821  /* This is the COMET discretization for cells that do not
822  do not have a face which is a reflecting boundary. When
823  there is a face with a reflecting boundary, we can no longer
824  use an arrowhead matrix as the structure becomes unknown a priori.
825  */
826  const int neibcount=_cellFaces.getCount(cell0);
827 
828  for(int j=0;j<neibcount;j++)
829  {
830  const int f=_cellFaces(cell0,j);
831  int cell1=_faceCells(f,1);
832  VectorT3 Af=_faceArea[f];
833  if(cell1==cell0)
834  {
835  cell1=_faceCells(f,0);
836  Af*=-1.;
837  }
838  const int klen=_kspace.getlength();
839 
840  T flux;
841 
842  for(int k=0;k<klen;k++)
843  {
844 
845  Tkvol& kvol=_kspace.getkvol(k);
846  const int numModes=kvol.getmodenum();
847  for(int m=0;m<numModes;m++)
848  {
849  Tmode& mode=kvol.getmode(m);
850  const int count=mode.getIndex();
851  const VectorT3 vg=mode.getv();
852  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
853 
854  if(flux>T_Scalar(0))
855  {
856  Amat.getElement(count,count)+=flux;
857  BVec[count-1]-=flux*_eArray[_kspace.getGlobalIndex(cell0,count-1)];
858  }
859  else
860  BVec[count-1]-=flux*_eArray[_kspace.getGlobalIndex(cell1,count-1)];
861  }
862  }
863 
864  }
865  }
866 
867  void COMETConvection(const int cell, TSquare& Amat, TArray& BVec)
868  {
869  const int neibcount=_cellFaces.getCount(cell);
870 
871  for(int j=0;j<neibcount;j++)
872  {
873  const int f=_cellFaces(cell,j);
874  int cell2=_faceCells(f,1);
875  const int klen=_kspace.getlength();
876  VectorT3 Af=_faceArea[f];
877 
878  T flux;
879 
880  if(cell2==cell)
881  {
882  Af=Af*(-1.);
883  cell2=_faceCells(f,0);
884  }
885 
886  if(_BCfArray[f]==2) //If the face in question is an implicit reflecting face
887  {
888  int Fgid=findFgId(f);
889  T refl=(*(_bcMap[Fgid]))["specifiedReflection"];
890  T oneMinusRefl=1.-refl;
891 
892  //first sweep - have to make sumVdotA
893  T sumVdotA=0.;
894  for(int k=0;k<klen;k++)
895  {
896  Tkvol& kvol=_kspace.getkvol(k);
897  T dk3=kvol.getdk3();
898  const int numModes=kvol.getmodenum();
899  for(int m=0;m<numModes;m++)
900  {
901  Tmode& mode=kvol.getmode(m);
902  const int count=mode.getIndex();
903  VectorT3 vg=mode.getv();
904  Field& efield=mode.getfield();
905  TArray& eArray=dynamic_cast<TArray&>(efield[_cells]);
906  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
907  if(flux>T_Scalar(0))
908  {
909  Amat(count,count)-=flux;
910  BVec[count-1]-=flux*eArray[cell];
911  }
912  else
913  {
914  sumVdotA-=flux*dk3;
915  Refl_pair& rpairs=mode.getReflpair(Fgid);
916  const int kk=rpairs.second.second;
917  Tkvol& kkvol=_kspace.getkvol(kk);
918  const int ccount=kkvol.getmode(m).getIndex();
919  Field& eefield=kkvol.getmode(m).getfield();
920  TArray& eeArray=dynamic_cast<TArray&>(eefield[_cells]);
921  Amat(count,ccount)-=flux*refl;
922  BVec[count-1]-=eeArray[cell]*refl*flux;
923  }
924  }
925  }
926 
927  T inverseSumVdotA=1./sumVdotA;
928 
929  //Second sweep
930  for(int k=0;k<klen;k++)
931  {
932  Tkvol& kvol=_kspace.getkvol(k);
933  const int numModes=kvol.getmodenum();
934  for(int m=0;m<numModes;m++)
935  {
936  const int count=kvol.getmode(m).getIndex();
937  VectorT3 vg=kvol.getmode(m).getv();
938  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
939  if(flux<T_Scalar(0))
940  {
941  for(int kk=0;kk<klen;kk++)
942  {
943  Tkvol& kkvol=_kspace.getkvol(kk);
944  const int numMODES=kkvol.getmodenum();
945  T ddk3=kkvol.getdk3();
946  for(int mm=0;mm<numMODES;mm++)
947  {
948  VectorT3 vvg=kkvol.getmode(mm).getv();
949  T VdotA=vvg[0]*Af[0]+vvg[1]*Af[1]+vvg[2]*Af[2];
950  if(VdotA>T_Scalar(0))
951  {
952  const int ccount=kkvol.getmode(mm).getIndex();
953  Field& eefield=kkvol.getmode(mm).getfield();
954  TArray& eeArray=dynamic_cast<TArray&>(eefield[_cells]);
955  Amat(count,ccount)-=flux*VdotA*ddk3*oneMinusRefl*inverseSumVdotA;
956  BVec[count-1]-=flux*VdotA*ddk3
957  *eeArray[cell]*inverseSumVdotA*oneMinusRefl;
958  }
959  }
960  }
961  }
962  }
963  }
964 
965 
966  }
967  else //if the face in question is not implicit
968  {
969  for(int k=0;k<klen;k++)
970  {
971  Tkvol& kvol=_kspace.getkvol(k);
972  const int numModes=kvol.getmodenum();
973  for(int m=0;m<numModes;m++)
974  {
975  Tmode& mode=kvol.getmode(m);
976  const int count=mode.getIndex();
977  VectorT3 vg=mode.getv();
978  Field& efield=mode.getfield();
979  TArray& eArray=dynamic_cast<TArray&>(efield[_cells]);
980  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
981  if(flux>T_Scalar(0))
982  {
983  Amat(count,count)-=flux;
984  BVec[count-1]-=flux*eArray[cell];
985  }
986  else
987  BVec[count-1]-=flux*eArray[cell2];
988  }
989  }
990  }
991 
992  }
993 
994  }
995 
996  void COMETCollision(const int cell, TMatrix* Amat, TArray& BVec)
997  {
998  const int klen=_kspace.getlength();
999  const int totalmodes=_kspace.gettotmodes();
1000  const int order=totalmodes+1;
1001  TArray& Tlold=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1002  T coeff;
1003  int cellIndex=_kspace.getGlobalIndex(cell,0);
1004 
1005  for(int k=0;k<klen;k++)
1006  {
1007  Tkvol& kvol=_kspace.getkvol(k);
1008  const int numModes=kvol.getmodenum();
1009  for(int m=0;m<numModes;m++)
1010  {
1011  Tmode& mode=kvol.getmode(m);
1012  const int count=mode.getIndex();
1013  const T tau=_kspace.getTau(cellIndex);
1014  T de0dT=mode.calcde0dT(Tlold[cell]);
1015  coeff=_cellVolume[cell]/tau;
1016  Amat->getElement(count,order)-=coeff*de0dT;
1017  Amat->getElement(count,count)+=coeff;
1018  BVec[count-1]-=coeff*_eArray[cellIndex];
1019  BVec[count-1]+=coeff*_e0Array[cellIndex];
1020  cellIndex++;
1021  }
1022  }
1023 
1024  }
1025 
1026  void COMETEquilibrium(const int cell, TMatrix* Amat, TArray& BVec)
1027  {
1028  const int klen=_kspace.getlength();
1029  const int totalmodes=_kspace.gettotmodes();
1030  const int order=totalmodes+1;
1031  TArray& Tlold=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1032  const T tauTot=_kspace.getde0taudT(cell,Tlold[cell]);
1033  T coeff;
1034  const T DK3=_kspace.getDK3();
1035  int cellIndex=_kspace.getGlobalIndex(cell,0);
1036  const T hbar=6.582119e-16; // (eV s)
1037 
1038  for(int k=0;k<klen;k++)
1039  {
1040  Tkvol& kvol=_kspace.getkvol(k);
1041  const T dk3=kvol.getdk3();
1042  const int numModes=kvol.getmodenum();
1043  for(int m=0;m<numModes;m++)
1044  {
1045  Tmode& mode=kvol.getmode(m);
1046  const int count=mode.getIndex();
1047  const T tau=_kspace.getTau(cellIndex);
1048  //const T energy=mode.getomega()*hbar;
1049  coeff=(dk3/DK3)/tau;
1050  Amat->getElement(order,count)-=coeff;
1051  BVec[totalmodes]+=coeff*_eArray[cellIndex];
1052  BVec[totalmodes]-=coeff*_e0Array[cellIndex];
1053  cellIndex++;
1054  }
1055  }
1056  Amat->getElement(order,order)=tauTot;
1057  }
1058 
1059  void COMETSource(const int cell, TArray& BVec)
1060  {_kspace.addSource(cell, BVec,_cellVolume[cell]);}
1061 
1062  void COMETFullScatt(const int cell, TArray& s, TArray& BVec)
1063  {
1064  const int klen=_kspace.getlength();
1065  const int totalmodes=_kspace.gettotmodes();
1066  //const int order=totalmodes+1;
1067  //TArray& Tlold=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1068  //const T DK3=_kspace.getDK3();
1069  int cellIndex=_kspace.getGlobalIndex(cell,0);
1070  TArray ds(totalmodes);
1071  //_kspace.getSourceTerm(cell,s,ds);
1072 
1073  const T relFac(1);
1074  //const T impl(1e0);
1075 
1076  for(int k=0;k<klen;k++)
1077  {
1078  Tkvol& kvol=_kspace.getkvol(k);
1079  //const T dk3=kvol.getdk3();
1080  const int numModes=kvol.getmodenum();
1081  for(int m=0;m<numModes;m++)
1082  {
1083  Tmode& mode=kvol.getmode(m);
1084  const T tau=_kspace.getTau(cellIndex);
1085  const int count=mode.getIndex();
1086  //const T coeff=_cellVolume[cell]/tau;
1087 
1088  BVec[count-1]+=(s[count-1])*_cellVolume[cell];
1089  //-(_e0Array[cellIndex]-_eArray[cellIndex])/tau)*_cellVolume[cell];
1090 
1091  BVec[count-1]*=relFac;
1092  cellIndex++;
1093  }
1094  }
1095  }
1096 
1097  void ScatterPhonons(const int cell0)
1098  {
1099 
1100  TArray C(_kspace.gettotmodes());
1101  TArray B(_kspace.gettotmodes());
1102  TArray V(_kspace.gettotmodes());
1103  TArray newE(_kspace.gettotmodes());
1104  C.zero();
1105  B.zero();
1106  V.zero();
1107  newE.zero();
1108 
1109  const int neibcount=_cellFaces.getCount(cell0);
1110 
1111  for(int j=0;j<neibcount;j++)
1112  {
1113  const int f=_cellFaces(cell0,j);
1114  int cell1=_faceCells(f,1);
1115  VectorT3 Af=_faceArea[f];
1116  if(cell1==cell0)
1117  {
1118  cell1=_faceCells(f,0);
1119  Af*=-1.;
1120  }
1121  const int klen=_kspace.getlength();
1122 
1123  T flux;
1124  for(int k=0;k<klen;k++)
1125  {
1126 
1127  Tkvol& kvol=_kspace.getkvol(k);
1128  const int numModes=kvol.getmodenum();
1129  for(int m=0;m<numModes;m++)
1130  {
1131  Tmode& mode=kvol.getmode(m);
1132  const int count=mode.getIndex();
1133  const VectorT3 vg=mode.getv();
1134  flux=vg[0]*Af[0]+vg[1]*Af[1]+vg[2]*Af[2];
1135 
1136  if(flux>T_Scalar(0)) //outgoing
1137  {
1138  V[count-1]+=flux/_cellVolume[cell0];
1139  B[count-1]+=flux*_e0Array[_kspace.getGlobalIndex(cell0,count-1)]/_cellVolume[cell0];
1140  C[count-1]+=flux*_eArray[_kspace.getGlobalIndex(cell0,count-1)]/_cellVolume[cell0];
1141  }
1142  else //incoming
1143  C[count-1]+=flux*_eArray[_kspace.getGlobalIndex(cell1,count-1)]/_cellVolume[cell0];
1144 
1145  }
1146  }
1147  }
1148 
1149  _kspace.ScatterPhonons(cell0, 50, C, B, V, newE, _cellVolume[cell0]);
1150 
1151  }
1152 
1153  void COMETShifted(const int cell, TMatrix* Amat, TArray& BVec)
1154  { //adds to collision and equilibrium
1155  const int klen=_kspace.getlength();
1156  const int totalmodes=_kspace.gettotmodes();
1157  const int order=totalmodes+1;
1158  TArray& Tlold=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1159  const T tauTot=_kspace.getde0taudT(cell,Tlold[cell]);
1160  T coeff;
1161 
1162  for(int k=0;k<klen;k++)
1163  {
1164  Tkvol& kvol=_kspace.getkvol(k);
1165  const int numModes=kvol.getmodenum();
1166  for(int m=0;m<numModes;m++)
1167  {
1168  Tmode& mode=kvol.getmode(m);
1169  const int count=mode.getIndex();
1170  T tau=mode.gettauN();
1171  Field& efield=mode.getfield();
1172  Field& eShiftedfield=mode.geteShifted();
1173  TArray& eArray=dynamic_cast<TArray&>(efield[_cells]);
1174  TArray& eShifted=dynamic_cast<TArray&>(eShiftedfield[_cells]);
1175  coeff=_cellVolume[cell]/tau;
1176  Amat->getElement(count,count)-=coeff;
1177  Amat->getElement(order,count)+=coeff/tauTot;
1178  BVec[count-1]-=coeff*eArray[cell];
1179  BVec[count-1]+=coeff*eShifted[cell];
1180  BVec[totalmodes]+=coeff*eArray[cell]/tauTot;
1181  BVec[totalmodes]-=coeff*eShifted[cell]/tauTot;
1182  }
1183  }
1184 
1185  }
1186 
1187  void Distribute(const int cell, TArray& BVec, TArray& Rvec)
1188  {
1189  const int klen=_kspace.getlength();
1190  const int totalmodes=_kspace.gettotmodes();
1191  int cellIndex=_kspace.getGlobalIndex(cell,0);
1192 
1193  for(int k=0;k<klen;k++)
1194  {
1195  Tkvol& kvol=_kspace.getkvol(k);
1196  const int numModes=kvol.getmodenum();
1197  for(int m=0;m<numModes;m++)
1198  {
1199  Tmode& mode=kvol.getmode(m);
1200  const int count=mode.getIndex()-1;
1201  _eArray[cellIndex]+=BVec[count];
1202  _resArray[cellIndex]=-Rvec[count];
1203  cellIndex++;
1204  }
1205  }
1206 
1207  TArray& TlArray=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1208  TlArray[cell]+=BVec[totalmodes];
1209  TArray& deltaTArray=dynamic_cast<TArray&>(_macro.deltaT[_cells]);
1210  deltaTArray[cell]=BVec[totalmodes];
1211  TArray& TlResArray=dynamic_cast<TArray&>(_macro.TlResidual[_cells]);
1212  TlResArray[cell]=-Rvec[totalmodes];
1213  }
1214 
1215  void findResidFine()
1216  {
1217  const int cellcount=_cells.getSelfCount();
1218  const int totalmodes=_kspace.gettotmodes();
1219  TArray ResidSum(totalmodes+1);
1220  TArray Bsum(totalmodes+1);
1221  T ResidScalar=0.;
1222  T traceSum=0.;
1223  TArray Bvec(totalmodes+1);
1224  TArray Resid(totalmodes+1);
1225  TArray Dummy(totalmodes+1);
1226  TArrow AMat(totalmodes+1);
1227  ResidSum.zero();
1228  Bsum.zero();
1229  Dummy.zero();
1230 
1231  GradMatrix& gradMatrix=GradModelType::getGradientMatrix(_mesh,_geomFields);
1232 
1233  for(int c=0;c<cellcount;c++)
1234  {
1235  Bvec.zero();
1236  Resid.zero();
1237  Dummy.zero();
1238 
1239  if(_BCArray[c]==0 || _BCArray[c]==2) //Arrowhead
1240  {
1241 
1242  //updateGhostFine(c, gradMatrix);
1243  //updateGhostCoarse(c);
1244 
1245  AMat.zero();
1246 
1247  COMETConvectionFine(c,AMat,Bvec,gradMatrix);
1248  COMETCollision(c,&AMat,Bvec);
1249  COMETEquilibrium(c,&AMat,Bvec);
1250 
1251  traceSum+=AMat.getTraceAbs();
1252 
1253  Resid=Bvec;
1254  Bvec.zero();
1255  Distribute(c,Bvec,Resid);
1256  ArrayAbs(Resid);
1257  ResidSum+=Resid;
1258 
1259  AMat.multiply(Resid,Bvec);
1260  Resid=Bvec;
1261 
1262  makeValueArray(c,Bvec);
1263  AMat.multiply(Bvec,Dummy);
1264  Bvec=Dummy;
1265  ArrayAbs(Bvec);
1266  ArrayAbs(Resid);
1267  Bsum+=Bvec;
1268  }
1269  else if(_BCArray[c]==1 || _BCArray[c]==3) //General Dense
1270  {
1271  TSquare AMatS(totalmodes+1);
1272  AMatS.zero();
1273 
1274  COMETConvection(c,AMatS,Bvec);
1275  COMETCollision(c,&AMatS,Bvec);
1276  COMETEquilibrium(c,&AMatS,Bvec);
1277 
1278  traceSum+=AMatS.getTraceAbs();
1279  Resid=Bvec;
1280  Bvec.zero();
1281  Distribute(c,Bvec,Resid);
1282 
1283  AMatS.multiply(Resid,Bvec);
1284  Resid=Bvec;
1285 
1286  makeValueArray(c,Bvec);
1287  ArrayAbs(Resid);
1288  ArrayAbs(Bvec);
1289  Bsum+=Bvec;
1290  ResidSum+=Resid;
1291  }
1292  else
1293  throw CException("Unexpected value for boundary cell map.");
1294  }
1295 
1296  //traceSum=0; //added
1297  for(int o=0;o<totalmodes+1;o++)
1298  {
1299  ResidScalar+=ResidSum[o];
1300  //traceSum+=Bsum[o]; //added
1301  }
1302 
1303  ResidScalar/=traceSum;
1304 
1305  if(_aveResid==-1)
1306  {_aveResid=ResidScalar;}
1307  else
1308  {
1309  _residChange=fabs(_aveResid-ResidScalar)/_aveResid;
1310  _aveResid=ResidScalar;
1311  }
1312  }
1313 
1314  void findResidCoarse(const bool plusFAS)
1315  {
1316  const int cellcount=_cells.getSelfCount();
1317  const int totalmodes=_kspace.gettotmodes();
1318  TArray ResidSum(totalmodes+1);
1319  TArray Bsum(totalmodes+1);
1320  T ResidScalar=0.;
1321  T traceSum=0.;
1322  TArray Bvec(totalmodes+1);
1323  TArray Resid(totalmodes+1);
1324  TArray Dummy(totalmodes+1);
1325  TArrow AMat(totalmodes+1);
1326  ResidSum.zero();
1327  Bsum.zero();
1328  Dummy.zero();
1329 
1330  for(int c=0;c<cellcount;c++)
1331  {
1332  Bvec.zero();
1333  Resid.zero();
1334  Dummy.zero();
1335 
1336  if(_BCArray[c]==0 || _BCArray[c]==2) //Arrowhead
1337  {
1338  if(_BCArray[c]==2)
1339  updateGhostCoarse(c);
1340 
1341  AMat.zero();
1342 
1343  COMETConvectionCoarse(c,AMat,Bvec);
1344  COMETCollision(c,&AMat,Bvec);
1345  COMETEquilibrium(c,&AMat,Bvec);
1346  COMETSource(c,Bvec);
1347 
1348  if(plusFAS)
1349  addFAS(c,Bvec);
1350 
1351  traceSum+=AMat.getTraceAbs();
1352 
1353  Resid=Bvec;
1354  Bvec.zero();
1355  Distribute(c,Bvec,Resid);
1356  ArrayAbs(Resid);
1357  ResidSum+=Resid;
1358 
1359  AMat.multiply(Resid,Bvec);
1360  Resid=Bvec;
1361 
1362  makeValueArray(c,Bvec);
1363  AMat.multiply(Bvec,Dummy);
1364  Bvec=Dummy;
1365  ArrayAbs(Bvec);
1366  ArrayAbs(Resid);
1367  Bsum+=Bvec;
1368  }
1369  else if(_BCArray[c]==1 || _BCArray[c]==3) //General Dense
1370  {
1371  TSquare AMatS(totalmodes+1);
1372  AMatS.zero();
1373 
1374  COMETConvection(c,AMatS,Bvec);
1375  COMETCollision(c,&AMatS,Bvec);
1376  COMETEquilibrium(c,&AMatS,Bvec);
1377 
1378  if(plusFAS)
1379  addFAS(c,Bvec);
1380 
1381  traceSum+=AMatS.getTraceAbs();
1382  Resid=Bvec;
1383  Bvec.zero();
1384  Distribute(c,Bvec,Resid);
1385 
1386  AMatS.multiply(Resid,Bvec);
1387  Resid=Bvec;
1388 
1389  makeValueArray(c,Bvec);
1390  ArrayAbs(Resid);
1391  ArrayAbs(Bvec);
1392  Bsum+=Bvec;
1393  ResidSum+=Resid;
1394  }
1395  else
1396  throw CException("Unexpected value for boundary cell map.");
1397  }
1398 
1399  //traceSum=0; //added
1400  for(int o=0;o<totalmodes+1;o++)
1401  {
1402  ResidScalar+=ResidSum[o];
1403  //traceSum+=Bsum[o]; //added
1404  }
1405 
1406  ResidScalar/=traceSum;
1407 
1408  if(_aveResid==-1)
1409  {_aveResid=ResidScalar;}
1410  else
1411  {
1412  _residChange=fabs(_aveResid-ResidScalar)/_aveResid;
1413  _aveResid=ResidScalar;
1414  }
1415  }
1416 
1417  void findResidFull()
1418  {
1419  const int cellcount=_cells.getSelfCount();
1420  const int totalmodes=_kspace.gettotmodes();
1421  TArray ResidSum(totalmodes+1);
1422  TArray Bsum(totalmodes+1);
1423  T ResidScalar=0.;
1424  T traceSum=0.;
1425  TArray Bvec(totalmodes+1);
1426  TArray Resid(totalmodes+1);
1427  TArray Dummy(totalmodes+1);
1428  TArray s(totalmodes);
1429  TArray ds(totalmodes);
1430  TArrow AMat(totalmodes+1);
1431  ResidSum.zero();
1432  Bsum.zero();
1433  Dummy.zero();
1434 
1435  for(int c=0;c<cellcount;c++)
1436  {
1437  Bvec.zero();
1438  Resid.zero();
1439  Dummy.zero();
1440  s.zero();
1441  ds.zero();
1442 
1443  _kspace.getSourceTerm(c,s,ds);
1444 
1445  if(_BCArray[c]==0 || _BCArray[c]==2) //Arrowhead
1446  {
1447  if(_BCArray[c]==2)
1448  updateGhostCoarse(c);
1449 
1450  AMat.zero();
1451 
1452  COMETConvectionCoarse(c,AMat,Bvec);
1453  COMETCollision(c,&AMat,Bvec);
1454  // COMETEquilibrium(c,&AMat,Bvec);
1455  COMETFullScatt(c,s,Bvec);
1456 
1457  traceSum+=AMat.getTraceAbs();
1458 
1459  Resid=Bvec;
1460  Bvec.zero();
1461  Distribute(c,Bvec,Resid);
1462  ArrayAbs(Resid);
1463  ResidSum+=Resid;
1464 
1465  AMat.multiply(Resid,Bvec);
1466  Resid=Bvec;
1467 
1468  makeValueArray(c,Bvec);
1469  AMat.multiply(Bvec,Dummy);
1470  Bvec=Dummy;
1471  ArrayAbs(Bvec);
1472  ArrayAbs(Resid);
1473  Bsum+=Bvec;
1474  }
1475  else if(_BCArray[c]==1 || _BCArray[c]==3) //General Dense
1476  {
1477  TSquare AMatS(totalmodes+1);
1478  AMatS.zero();
1479 
1480  COMETConvection(c,AMatS,Bvec);
1481  COMETCollision(c,&AMatS,Bvec);
1482  COMETEquilibrium(c,&AMatS,Bvec);
1483 
1484  traceSum+=AMatS.getTraceAbs();
1485  Resid=Bvec;
1486  Bvec.zero();
1487  Distribute(c,Bvec,Resid);
1488 
1489  AMatS.multiply(Resid,Bvec);
1490  Resid=Bvec;
1491 
1492  makeValueArray(c,Bvec);
1493  ArrayAbs(Resid);
1494  ArrayAbs(Bvec);
1495  Bsum+=Bvec;
1496  ResidSum+=Resid;
1497  }
1498  else
1499  throw CException("Unexpected value for boundary cell map.");
1500  }
1501 
1502  //traceSum=0; //added
1503  for(int o=0;o<totalmodes+1;o++)
1504  {
1505  ResidScalar+=ResidSum[o];
1506  //traceSum+=Bsum[o]; //added
1507  }
1508 
1509  ResidScalar/=traceSum;
1510 
1511  if(_aveResid==-1)
1512  {_aveResid=ResidScalar;}
1513  else
1514  {
1515  _residChange=fabs(_aveResid-ResidScalar)/_aveResid;
1516  _aveResid=ResidScalar;
1517  }
1518  }
1519 
1520  TArray gatherResid(const int c)
1521  {
1522  const int klen=_kspace.getlength();
1523  const int totalmodes=_kspace.gettotmodes();
1524  const int order=totalmodes+1;
1525 
1526  TArray rArray(order);
1527 
1528  for(int k=0;k<klen;k++)
1529  {
1530  Tkvol& kvol=_kspace.getkvol(k);
1531  const int numModes=kvol.getmodenum();
1532  for(int m=0;m<numModes;m++)
1533  {
1534  Tmode& mode=kvol.getmode(m);
1535  const int count=mode.getIndex();
1536  Field& resfield=mode.getresid();
1537  TArray& resArray=dynamic_cast<TArray&>(resfield[_cells]);
1538  rArray[count]=resArray[c];
1539  }
1540  }
1541  rArray[order]=_macro.TlResidual[c];
1542  }
1543 
1544  T getResidChange() {return _residChange;}
1545  T getAveResid() {return _aveResid;}
1546 
1547  void setfgFinder()
1548  {
1549  foreach(const FaceGroupPtr fgPtr, _mesh.getBoundaryFaceGroups())
1550  {
1551  const FaceGroup& fg=*fgPtr;
1552  const int off=fg.site.getOffset();
1553  const int cnt=fg.site.getCount();
1554  const int id=fg.id;
1555  VecInt2 BegEnd;
1556 
1557  BegEnd[0]=off;
1558  BegEnd[1]=off+cnt;
1559  _fgFinder[id]=BegEnd;
1560  }
1561  }
1562 
1563  int findFgId(const int faceIndex)
1564  {
1565  FaceToFg::iterator id;
1566  for(id=_fgFinder.begin();id!=_fgFinder.end();id++)
1567  {
1568  if(id->second[0]<=faceIndex && id->second[1]>faceIndex)
1569  return id->first;
1570  }
1571  cout<<"Face index: "<<faceIndex<<endl;
1572  throw CException("Didn't find matching FaceGroup!");
1573  return -1;
1574  }
1575 
1576  void ArrayAbs(TArray& o)
1577  {
1578  int length=o.getLength();
1579  for(int i=0;i<length;i++)
1580  o[i]=fabs(o[i]);
1581  }
1582 
1583  void makeValueArray(const int c, TArray& o)
1584  {
1585  int klen=_kspace.getlength();
1586  const int totmodes=_kspace.gettotmodes();
1587  int cellIndex=_kspace.getGlobalIndex(c,0);
1588  for(int k=0;k<klen;k++)
1589  {
1590  Tkvol& kvol=_kspace.getkvol(k);
1591  const int numModes=kvol.getmodenum();
1592  for(int m=0;m<numModes;m++)
1593  {
1594  Tmode& mode=kvol.getmode(m);
1595  const int count=mode.getIndex()-1;
1596  o[count]=_eArray[cellIndex];
1597  cellIndex++;
1598  }
1599  }
1600  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1601  o[totmodes]=Tl[c];
1602  }
1603 
1604  void addFAS(const int c, TArray& bVec)
1605  {
1606  const int totmodes=_kspace.gettotmodes();
1607  _kspace.addFAS(c,bVec);
1608  TArray& fasArray=dynamic_cast<TArray&>(_macro.TlFASCorrection[_cells]);
1609  bVec[totmodes]-=fasArray[c];
1610  }
1611 
1612  void updatee0()
1613  {
1614  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1615 
1616  const T cellCount=_cells.getSelfCount();
1617  int klen=_kspace.getlength();
1618  for(int c=0;c<cellCount;c++)
1619  {
1620  int cellIndex=_kspace.getGlobalIndex(c,0);
1621  for(int k=0;k<klen;k++)
1622  {
1623  Tkvol& kvol=_kspace.getkvol(k);
1624  const int numModes=kvol.getmodenum();
1625  for(int m=0;m<numModes;m++)
1626  {
1627  Tmode& mode=kvol.getmode(m);
1628  _e0Array[cellIndex]=mode.calce0(Tl[c]);
1629  cellIndex++;
1630  }
1631  }
1632  }
1633  }
1634 
1635  void updatee0(const int c)
1636  {
1637  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1638 
1639  int klen=_kspace.getlength();
1640  int cellIndex=_kspace.getGlobalIndex(c,0);
1641  for(int k=0;k<klen;k++)
1642  {
1643  Tkvol& kvol=_kspace.getkvol(k);
1644  const int numModes=kvol.getmodenum();
1645  for(int m=0;m<numModes;m++)
1646  {
1647  Tmode& mode=kvol.getmode(m);
1648  _e0Array[cellIndex]=mode.calce0(Tl[c]);
1649  _kspace.updateTau(c,Tl[c]);
1650  cellIndex++;
1651  }
1652  }
1653  }
1654 
1655  void updateGhostFine(const int cell, const GradMatrix& gMat)
1656  {
1657  const int neibcount=_cellFaces.getCount(cell);
1658  const int klen=_kspace.getlength();
1659  //TArray SumEVdotA(neibcount);
1660  //TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1661  //T DK3=_kspace.getDK3();
1662 
1663  //SumEVdotA.zero();
1664 
1665  const VectorT3Array& faceCoords=
1666  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[_faces]);
1667 
1668  TArray pointMin(_kspace.gettotmodes());
1669  pointMin=-1;
1670  TArray pointMax(_kspace.gettotmodes());
1671  pointMax.zero();
1672  TArray pointLim(_kspace.gettotmodes());
1673  pointLim=100.;
1674 
1675  GradArray Grads(_kspace.gettotmodes());
1676  Grads.zero();
1677 
1678  VectorT3 Gcoeff;
1679 
1680  for(int j=0;j<neibcount;j++) //first loop to get grad and max/min vals
1681  {
1682  const int f=_cellFaces(cell,j);
1683  int cell1=_faceCells(f,1);
1684  if(cell1==cell)
1685  cell1=_faceCells(f,0);
1686  const VectorT3 Gcoeff=gMat.getCoeff(cell, cell1);
1687 
1688  int c0ind=_kspace.getGlobalIndex(cell,0);
1689  int c1ind=_kspace.getGlobalIndex(cell1,0);
1690 
1691  for(int k=0;k<_kspace.gettotmodes();k++)
1692  {
1693  const T e1=_eArray[c1ind];
1694  const T e0=_eArray[c0ind];
1695  T& curMax=pointMax[k];
1696  T& curMin=pointMin[k];
1697  Grads[k].accumulate(Gcoeff, e1-e0);
1698 
1699  if(e1>curMax)
1700  curMax=e1;
1701 
1702  if(curMin==-1)
1703  curMin=e1;
1704  else
1705  {
1706  if(e1<curMin)
1707  curMin=e1;
1708  }
1709 
1710  c0ind++;
1711  c1ind++;
1712  }
1713  }
1714 
1715  for(int j=0;j<neibcount;j++) //second loop to calculate limiting coeff
1716  {
1717  const int f=_cellFaces(cell,j);
1718  int cell1=_faceCells(f,1);
1719  if(cell1==cell)
1720  cell1=_faceCells(f,0);
1721 
1722  const VectorT3 dr0(faceCoords[f]-_cellCoords[cell]);
1723  int c0ind=_kspace.getGlobalIndex(cell,0);
1724  vanLeer lf;
1725 
1726  for(int k=0;k<_kspace.gettotmodes();k++)
1727  {
1728  const T minVal=pointMin[k];
1729  const T maxVal=pointMax[k];
1730  const T de0(Grads[k]*dr0);
1731  T& cl=pointLim[k];
1732  computeLimitCoeff2(cl, _eArray[c0ind], de0, minVal, maxVal, lf);
1733  c0ind++;
1734  }
1735  }
1736 
1737 
1738  for(int j=0;j<neibcount;j++)
1739  {
1740 
1741  const int f=_cellFaces(cell,j);
1742 
1743  if(_BCfArray[f]==1)
1744  {//temperature
1745  //int Fgid=findFgId(f);
1746  int cell2=_faceCells(f,1);
1747  VectorT3 Af=_faceArea[f];
1748 
1749  if(cell2==cell)
1750  {
1751  Af=Af*-1.;
1752  cell2=_faceCells(f,0);
1753  }
1754 
1755  int cellIndex=_kspace.getGlobalIndex(cell,0);
1756  int cell2Index=_kspace.getGlobalIndex(cell2,0);
1757  //vanLeer vl;
1758 
1759  for(int k=0;k<klen;k++)
1760  {
1761  Tkvol& kvol=_kspace.getkvol(k);
1762  const int numModes=kvol.getmodenum();
1763  //T dk3=kvol.getdk3();
1764  for(int m=0;m<numModes;m++)
1765  {
1766  Tmode& mode=kvol.getmode(m);
1767  VectorT3 vg=mode.getv();
1768  const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
1769  const int count=mode.getIndex();
1770 
1771  GradType& grad=Grads[count-1];
1772 
1773  if(VdotA>0)
1774  {
1775  const VectorT3 fVec=faceCoords[f]-_cellCoords[cell];
1776  const VectorT3 rVec=_cellCoords[cell2]-_cellCoords[cell];
1777  //const T r=gMat.computeR(grad,_eArray,rVec,cellIndex,cell2Index);
1778 
1779  T SOU=(fVec[0]*grad[0]+fVec[1]*grad[1]+
1780  fVec[2]*grad[2])*pointLim[count-1];
1781  _eArray[cell2Index]=_eArray[cellIndex]+SOU;
1782  }
1783  cellIndex++;
1784  cell2Index++;
1785  }
1786  }
1787 
1788  }
1789  else if(_BCfArray[f]==3)
1790  {//reflecting
1791  int Fgid=findFgId(f);
1792  const T refl=(*(_bcMap[Fgid]))["specifiedReflection"];
1793  int cell2=_faceCells(f,1);
1794  VectorT3 Af=_faceArea[f];
1795 
1796  if(cell2==cell)
1797  {
1798  Af=Af*-1.;
1799  cell2=_faceCells(f,0);
1800  }
1801 
1802  //const int cellstart=_kspace.getGlobalIndex(cell,0);
1803  const int cell2start=_kspace.getGlobalIndex(cell2,0);
1804  int cellIndex=_kspace.getGlobalIndex(cell,0);
1805  int cell2Index=_kspace.getGlobalIndex(cell2,0);
1806  vanLeer vl;
1807 
1808  //first loop to sum up the diffuse reflection
1809  for(int k=0;k<klen;k++)
1810  {
1811  Tkvol& kvol=_kspace.getkvol(k);
1812  const int numModes=kvol.getmodenum();
1813  //T dk3=kvol.getdk3();
1814  for(int m=0;m<numModes;m++)
1815  {
1816  Tmode& mode=kvol.getmode(m);
1817  VectorT3 vg=mode.getv();
1818  const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
1819  const int count=mode.getIndex();
1820 
1821  GradType& grad=Grads[count-1];
1822 
1823  if(VdotA>0)
1824  {
1825  VectorT3 rVec=_cellCoords[cell2]-_cellCoords[cell];
1826  VectorT3 fVec=faceCoords[f]-_cellCoords[cell];
1827  //T r=gMat.computeR(grad,_eArray,rVec,cellIndex,cell2Index);
1828  T SOU=(grad[0]*fVec[0]+grad[1]*fVec[1]+grad[2]*fVec[2])*pointLim[count-1];
1829  _eArray[cell2Index]=(_eArray[cellIndex]+SOU)*refl;
1830  }
1831  cellIndex++;
1832  cell2Index++;
1833  }
1834  }
1835 
1836  cellIndex=_kspace.getGlobalIndex(cell,0);
1837  cell2Index=_kspace.getGlobalIndex(cell2,0);
1838 
1839  //first loop to sum up the diffuse reflection
1840  for(int k=0;k<klen;k++)
1841  {
1842  Tkvol& kvol=_kspace.getkvol(k);
1843  const int numModes=kvol.getmodenum();
1844  //T dk3=kvol.getdk3();
1845  for(int m=0;m<numModes;m++)
1846  {
1847  Tmode& mode=kvol.getmode(m);
1848  VectorT3 vg=mode.getv();
1849  const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
1850  //const int count=mode.getIndex();
1851 
1852  if(VdotA<0)
1853  {
1854  Refl_pair& rpairs=mode.getReflpair(Fgid);
1855  const int kk=rpairs.second.second;
1856  Tkvol& kkvol=_kspace.getkvol(kk);
1857  Tmode& refMode=kkvol.getmode(m);
1858  const int rIndex=cell2start+refMode.getIndex()-1;
1859  _eArray[cell2Index]=_eArray[rIndex]*refl;
1860  }
1861  cellIndex++;
1862  cell2Index++;
1863  }
1864  }
1865 
1866  }
1867  }
1868 
1869  /*
1870  T wallTemp(Tl[cell]);
1871  for(int j=0;j<neibcount;j++)
1872  {
1873  const int f=_cellFaces(cell,j);
1874  if(_BCfArray[f]==3)
1875  {
1876  int cell2=_faceCells(f,1);
1877  VectorT3 Af=_faceArea[f];
1878  if(cell2==cell)
1879  {
1880  Af=Af*-1.;
1881  cell2=_faceCells(f,0);
1882  }
1883  T esum=SumEVdotA[j]*DK3;
1884  _kspace.calcTemp(wallTemp,esum,Af);
1885  SumEVdotA[j]=wallTemp;
1886  Tl[cell2]=wallTemp;
1887  }
1888  }
1889 
1890  //second loop to add in the diffuse reflection
1891  for(int k=0;k<klen;k++)
1892  {
1893  Tkvol& kvol=_kspace.getkvol(k);
1894  const int numModes=kvol.getmodenum();
1895  for(int m=0;m<numModes;m++)
1896  {
1897  Tmode& mode=kvol.getmode(m);
1898  VectorT3 vg=mode.getv();
1899  Field& eField=mode.getfield();
1900  TArray& eArray=dynamic_cast<TArray&>(eField[_cells]);
1901  for(int j=0;j<neibcount;j++)
1902  {
1903  const int f=_cellFaces(cell,j);
1904  if(_BCfArray[f]==3)
1905  {
1906  int Fgid=findFgId(f);
1907  const T refl=(*(_bcMap[Fgid]))["specifiedReflection"];
1908  const T oneMinusRefl=1.-refl;
1909  int cell2=_faceCells(f,1);
1910  VectorT3 Af=_faceArea[f];
1911 
1912  if(cell2==cell)
1913  {
1914  Af=Af*-1.;
1915  cell2=_faceCells(f,0);
1916  }
1917 
1918  const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
1919 
1920  if(VdotA<0)
1921  eArray[cell2]+=mode.calce0(SumEVdotA[j])*oneMinusRefl;
1922  else
1923  eArray[cell2]+=mode.calce0(SumEVdotA[j])*oneMinusRefl;
1924 
1925  }
1926  }
1927  }
1928  }
1929  */
1930  }
1931 
1932  void updateGhostCoarse(const int cell)
1933  {
1934  const int neibcount=_cellFaces.getCount(cell);
1935  const int klen=_kspace.getlength();
1936  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
1937  T DK3=_kspace.getDK3();
1938 
1939  for(int j=0;j<neibcount;j++)
1940  {
1941 
1942  const int f=_cellFaces(cell,j);
1943  if(_BCfArray[f]==3)
1944  {
1945  T SumEVdotA(0);
1946  int Fgid=findFgId(f);
1947  const T refl=(*(_bcMap[Fgid]))["specifiedReflection"];
1948  int cell2=_faceCells(f,1);
1949  VectorT3 Af=_faceArea[f];
1950 
1951  if(cell2==cell)
1952  {
1953  Af=Af*-1.;
1954  cell2=_faceCells(f,0);
1955  }
1956 
1957  const int cellstart=_kspace.getGlobalIndex(cell,0);
1958  int cellIndex=_kspace.getGlobalIndex(cell,0);
1959  int cell2Index=_kspace.getGlobalIndex(cell2,0);
1960 
1961  //first loop to sum up the diffuse reflection
1962  for(int k=0;k<klen;k++)
1963  {
1964  Tkvol& kvol=_kspace.getkvol(k);
1965  const int numModes=kvol.getmodenum();
1966  T dk3=kvol.getdk3();
1967  for(int m=0;m<numModes;m++)
1968  {
1969  Tmode& mode=kvol.getmode(m);
1970  VectorT3 vg=mode.getv();
1971  const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
1972 
1973  if(VdotA>0)
1974  {
1975  _eArray[cell2Index]=_eArray[cellIndex]*refl;
1976  SumEVdotA+=_eArray[cellIndex]*VdotA*(dk3/DK3);
1977  }
1978  else
1979  {
1980  if(refl==1.)
1981  {
1982  Refl_pair& rpairs=mode.getReflpair(Fgid);
1983  const int kk=rpairs.second.second;
1984  Tkvol& kkvol=_kspace.getkvol(kk);
1985  Tmode& refMode=kkvol.getmode(m);
1986  const int rIndex=cellstart+refMode.getIndex()-1;
1987  _eArray[cell2Index]=_eArray[rIndex]*refl;
1988  }
1989  else
1990  _eArray[cell2Index]=0;
1991  }
1992  cellIndex++;
1993  cell2Index++;
1994  }
1995  }
1996 
1997  //second loop to add in the diffuse reflection
1998  if(refl==0.)
1999  {
2000 
2001  T wallTemp=Tl[cell2];
2002  T esum=SumEVdotA*DK3;
2003  _kspace.calcTemp(wallTemp,esum,Af);
2004  SumEVdotA=wallTemp;
2005  Tl[cell2]=wallTemp;
2006 
2007  for(int k=0;k<klen;k++)
2008  {
2009  Tkvol& kvol=_kspace.getkvol(k);
2010  const int numModes=kvol.getmodenum();
2011  for(int m=0;m<numModes;m++)
2012  {
2013  Tmode& mode=kvol.getmode(m);
2014  VectorT3 vg=mode.getv();
2015  const int Index=mode.getIndex()-1;
2016  const int cell2Index=_kspace.getGlobalIndex(cell2,Index);
2017 
2018  const T oneMinusRefl=1.-refl;
2019  //const T VdotA=Af[0]*vg[0]+Af[1]*vg[1]+Af[2]*vg[2];
2020 
2021  _eArray[cell2Index]+=mode.calce0(Tl[cell2])*oneMinusRefl;
2022 
2023  }
2024  }
2025  }
2026 
2027  }
2028  }
2029 
2030  }
2031 
2032  void correctInterface(const int cell0, TArray& Bvec)
2033  {
2034  const int klen=_kspace.gettotmodes();
2035  const T DK3=_kspace.getDK3();
2036  TArray Correction(klen+1);
2037 
2038  const int neibcount=_cellFaces.getCount(cell0);
2039  for(int j=0;j<neibcount;j++)
2040  {
2041  const int f=_cellFaces(cell0,j);
2042  if(_BCfArray[f]==4)
2043  {
2044  int Fgid=findFgId(f);
2045  const FaceGroup& fg=_mesh.getFaceGroup(Fgid);
2046  const StorageSite& faces0=fg.site;
2047  const int offset=faces0.getOffset();
2048  int cell1=_faceCells(f,1);
2049  if(cell1==cell0)
2050  cell1=_faceCells(f,0);
2051  const TKConnectivity& TKC=*(_FaceToKSC.find(Fgid)->second)[f-offset];
2052  TKC.multiplySelf(Bvec,Correction,DK3);
2053  Bvec.zero();
2054  Distribute(cell1,Correction, Bvec);
2055  }
2056  }
2057 
2058  }
2059 
2060  void updateeShifted()
2061  {
2062  const int cellcount=_cells.getSelfCount();
2063  TArray& Tl=dynamic_cast<TArray&>(_macro.temperature[_cells]);
2064  VectorT3Array& lam=dynamic_cast<VectorT3Array&>(_macro.lam[_cells]);
2065 
2066  for(int c=0;c<cellcount;c++)
2067  {
2068  T3Tensor TensorSum;
2069  VectorT3 VectorSum;
2070  VectorT3 lambda;
2071  T TpreFactor=0.;
2072  T VpreFactor=0.;
2073 
2074  int klen=_kspace.getlength();
2075  for(int k=0;k<klen;k++)
2076  {
2077  Tkvol& kvol=_kspace.getkvol(k);
2078  const int numModes=kvol.getmodenum();
2079  for(int m=0;m<numModes;m++)
2080  {
2081  Tmode& mode=kvol.getmode(m);
2082  Field& eField=mode.getfield();
2083  TArray& eArray=dynamic_cast<TArray&>(eField[_cells]);
2084  TpreFactor+=mode.calcTensorPrefactor(Tl[c]);
2085  VpreFactor+=mode.calcVectorPrefactor(Tl[c],eArray[c]);
2086  }
2087  }
2088 
2089  TensorSum.zero();
2090  VectorSum.zero();
2091 
2092  for(int k=0;k<klen;k++)
2093  {
2094  Tkvol& kvol=_kspace.getkvol(k);
2095  T dk3=kvol.getdk3();
2096  VectorT3 Kvec=kvol.getkvec();
2097  T3Tensor TempTensor;
2098  outerProduct(Kvec,Kvec,TempTensor);
2099  TensorSum+=TempTensor*dk3*TpreFactor;
2100  VectorSum+=Kvec*dk3*VpreFactor;
2101  }
2102 
2103  lambda=inverse(TensorSum)*VectorSum;
2104  lam[c]=lambda;
2105 
2106  for(int k=0;k<klen;k++)
2107  {
2108  Tkvol& kvol=_kspace.getkvol(k);
2109  VectorT3 Kvec=kvol.getkvec();
2110  T shift=Kvec[0]*lambda[0]+Kvec[1]*lambda[1]+Kvec[2]*lambda[2];
2111  const int numModes=kvol.getmodenum();
2112  for(int m=0;m<numModes;m++)
2113  {
2114  Tmode& mode=kvol.getmode(m);
2115  Field& eShiftedField=mode.geteShifted();
2116  TArray& eShifted=dynamic_cast<TArray&>(eShiftedField[_cells]);
2117  eShifted[c]=mode.calcShifted(Tl[c],shift);
2118  }
2119  }
2120 
2121  }
2122  }
2123 
2124  void outerProduct(const VectorT3& v1, const VectorT3& v2, T3Tensor& out)
2125  {
2126  for(int i=0;i<3;i++)
2127  for(int j=0;j<3;j++)
2128  out(i,j)=v1[i]*v2[j];
2129  }
2130 
2131  T scaledResid(const TArray& de, const int c)
2132  {
2133  T scaled(0);
2134  int cellIndex(_kspace.getGlobalIndex(c,0));
2135  for(int i=0;i<de.getLength()-1;i++)
2136  {
2137  scaled+=fabs(de[i])/_eArray[cellIndex];
2138  cellIndex++;
2139  }
2140 
2141  return scaled/de.getLength();
2142  }
2143 
2144  private:
2145 
2146  const Mesh& _mesh;
2147  const GeomFields& _geomFields;
2148  const StorageSite& _cells;
2149  const StorageSite& _faces;
2150  const CRConnectivity& _cellFaces;
2151  const CRConnectivity& _faceCells;
2152  const TArray& _faceAreaMag;
2153  const VectorT3Array& _faceArea;
2154  const TArray& _cellVolume;
2155  const VectorT3Array& _cellCoords;
2156  PhononMacro& _macro;
2157  Tkspace& _kspace;
2158  COMETBCMap& _bcMap;
2159  const IntArray& _BCArray;
2160  const IntArray& _BCfArray;
2161  T _aveResid;
2162  T _residChange;
2163  FaceToFg _fgFinder;
2164  COpts _options;
2165  const FgTKClistMap& _FaceToKSC;
2166  TArray& _eArray;
2167  TArray& _e0Array;
2168  TArray& _resArray;
2169 
2170 };
2171 
2172 
2173 #endif
SquareMatrix::multiply
void multiply(const TArray &x, TArray &b)
Definition: SquareMatrix.h:194
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
pmode::Refl_pair
pair< Reflection, Reflection > Refl_pair
Definition: pmode.h:29
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
pmode::getReflpair
Refl_pair & getReflpair(int i)
Definition: pmode.h:71
COMETDiscretizer::_cellCoords
const VectorT3Array & _cellCoords
Definition: COMETDiscretizer.h:2155
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
COMETDiscretizer::makeValueArray
void makeValueArray(const int c, TArray &o)
Definition: COMETDiscretizer.h:1583
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
COMETDiscretizer::_macro
PhononMacro & _macro
Definition: COMETDiscretizer.h:2156
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
COMETDiscretizer::COMETConvection
void COMETConvection(const int cell, TSquare &Amat, TArray &BVec)
Definition: COMETDiscretizer.h:867
COMETDiscretizer::TArrow
ArrowHeadMatrix< T > TArrow
Definition: COMETDiscretizer.h:43
pmode::getv
Tvec getv()
Definition: pmode.h:59
Kspace
Definition: Kspace.h:28
COMETDiscretizer::_faceCells
const CRConnectivity & _faceCells
Definition: COMETDiscretizer.h:2151
COMETDiscretizer::COMETSolveFull
void COMETSolveFull(const int dir, const int level)
Definition: COMETDiscretizer.h:403
Kspace::updateTau
void updateTau(const int c, const T Tl)
Definition: Kspace.h:1491
COMETDiscretizer::VecInt2
Vector< int, 2 > VecInt2
Definition: COMETDiscretizer.h:49
COMETDiscretizer::_BCArray
const IntArray & _BCArray
Definition: COMETDiscretizer.h:2159
FaceGroup
Definition: Mesh.h:28
Kspace::getde0taudT
T getde0taudT(const int c, T Tl)
Definition: Kspace.h:646
pmode::calce0
T calce0(T Tl)
Definition: pmode.h:88
Field
Definition: Field.h:14
Kspace::getkvol
Tkvol & getkvol(int n) const
Definition: Kspace.h:390
COMETDiscretizer::COMETConvectionFine
void COMETConvectionFine(const int cell0, TArrow &Amat, TArray &BVec, const GradMatrix &gMat)
Definition: COMETDiscretizer.h:579
kvol::getmode
Tmode & getmode(int n) const
Definition: kvol.h:44
Kspace::getlength
int getlength() const
Definition: Kspace.h:391
PhononMacro::deltaT
Field deltaT
Definition: PhononMacro.h:22
COMETDiscretizer::ScatterPhonons
void ScatterPhonons(const int cell0)
Definition: COMETDiscretizer.h:1097
COMETDiscretizer::TKConnectivity
KSConnectivity< T > TKConnectivity
Definition: COMETDiscretizer.h:53
PhononMacro::temperature
Field temperature
Definition: PhononMacro.h:21
COMETDiscretizer::Distribute
void Distribute(const int cell, TArray &BVec, TArray &Rvec)
Definition: COMETDiscretizer.h:1187
COMETDiscretizer::_faces
const StorageSite & _faces
Definition: COMETDiscretizer.h:2149
COMETDiscretizer::T_Scalar
NumTypeTraits< T >::T_Scalar T_Scalar
Definition: COMETDiscretizer.h:35
Mesh::getFaceGroup
const FaceGroup & getFaceGroup(const int fgId) const
Definition: Mesh.cpp:1570
Mesh
Definition: Mesh.h:49
COMETDiscretizer::COMETSource
void COMETSource(const int cell, TArray &BVec)
Definition: COMETDiscretizer.h:1059
COMETDiscretizer::gatherResid
TArray gatherResid(const int c)
Definition: COMETDiscretizer.h:1520
kvol::getmodenum
int getmodenum()
Definition: kvol.h:43
Kspace::getDK3
T getDK3() const
Definition: Kspace.h:408
COMETDiscretizer::T3Tensor
SquareTensor< T, 3 > T3Tensor
Definition: COMETDiscretizer.h:52
GradientMatrix::getCoeff
Coord & getCoeff(const int i, const int j)
Definition: GradientMatrix.h:129
COMETDiscretizer::COMETFullScatt
void COMETFullScatt(const int cell, TArray &s, TArray &BVec)
Definition: COMETDiscretizer.h:1062
COMETDiscretizer::scaledResid
T scaledResid(const TArray &de, const int c)
Definition: COMETDiscretizer.h:2131
COMETDiscretizer::BoolArray
Array< bool > BoolArray
Definition: COMETDiscretizer.h:48
MatrixJML::getElement
virtual T & getElement(const int i, const int j)=0
COMETDiscretizer::IntArray
Array< int > IntArray
Definition: COMETDiscretizer.h:47
Kspace::calcTemp
void calcTemp(T &guess, const T e_sum, const Tvec An)
Definition: Kspace.h:443
COMETDiscretizer::COMETCollision
void COMETCollision(const int cell, TMatrix *Amat, TArray &BVec)
Definition: COMETDiscretizer.h:996
COMETDiscretizer::ArrayAbs
void ArrayAbs(TArray &o)
Definition: COMETDiscretizer.h:1576
COMETDiscretizer::COMETBCMap
map< int, COMETBC< T > * > COMETBCMap
Definition: COMETDiscretizer.h:45
COMETDiscretizer::GradArray
Array< GradType > GradArray
Definition: COMETDiscretizer.h:58
COMETDiscretizer::TKCptr
TKConnectivity * TKCptr
Definition: COMETDiscretizer.h:54
COMETDiscretizer::Refl_pair
Tmode::Refl_pair Refl_pair
Definition: COMETDiscretizer.h:51
FaceGroup::id
const int id
Definition: Mesh.h:41
COMETDiscretizer::TArray
Array< T_Scalar > TArray
Definition: COMETDiscretizer.h:41
Kspace::addFAS
void addFAS(const int c, TArray &Bvec)
Definition: Kspace.h:1273
pmode::gettauN
T gettauN()
Definition: pmode.h:62
ArrowHeadMatrix::multiply
void multiply(const XArray &x, XArray &b)
Definition: esbgkbase/ArrowHeadMatrix.h:108
pmode::getIndex
int getIndex()
Definition: pmode.h:73
PhononMacro::TlResidual
Field TlResidual
Definition: PhononMacro.h:24
Kspace::getGlobalIndex
int getGlobalIndex(const int cell, const int count)
Definition: Kspace.h:1216
ArrowHeadMatrix::getElement
T_Scalar & getElement(const int i, const int j)
Definition: esbgkbase/ArrowHeadMatrix.h:50
COMETDiscretizer::_BCfArray
const IntArray & _BCfArray
Definition: COMETDiscretizer.h:2160
COMETDiscretizer::outerProduct
void outerProduct(const VectorT3 &v1, const VectorT3 &v2, T3Tensor &out)
Definition: COMETDiscretizer.h:2124
COMETDiscretizer::TMatrix
MatrixJML< T > TMatrix
Definition: COMETDiscretizer.h:42
COMETDiscretizer::VectorT3
Vector< T_Scalar, 3 > VectorT3
Definition: COMETDiscretizer.h:39
COMETDiscretizer::_cellVolume
const TArray & _cellVolume
Definition: COMETDiscretizer.h:2154
COMETDiscretizer::addFAS
void addFAS(const int c, TArray &bVec)
Definition: COMETDiscretizer.h:1604
PhononMacro::TlFASCorrection
Field TlFASCorrection
Definition: PhononMacro.h:26
COMETDiscretizer::updatee0
void updatee0(const int c)
Definition: COMETDiscretizer.h:1635
COMETDiscretizer::COMETDiscretizer
COMETDiscretizer(const Mesh &mesh, const GeomFields &geomfields, PhononMacro &macro, Tkspace &kspace, COMETBCMap &bcMap, const IntArray &BCArray, const IntArray &BCfArray, COpts &options, const FgTKClistMap &FgToKsc)
Definition: COMETDiscretizer.h:62
COMETDiscretizer::COpts
COMETModelOptions< T > COpts
Definition: COMETDiscretizer.h:46
COMETDiscretizer::TSquare
SquareMatrix< T > TSquare
Definition: COMETDiscretizer.h:44
Kspace::getTau
const T getTau(const int index)
Definition: Kspace.h:1231
COMETDiscretizer::_cellFaces
const CRConnectivity & _cellFaces
Definition: COMETDiscretizer.h:2150
computeLimitCoeff2
void computeLimitCoeff2(T &lc, const T x, const T &dx, const T &min, const T &max, const LimitFunc &f)
Definition: FluxLimiters.h:85
COMETDiscretizer::correctInterface
void correctInterface(const int cell0, TArray &Bvec)
Definition: COMETDiscretizer.h:2032
COMETDiscretizer::_faceAreaMag
const TArray & _faceAreaMag
Definition: COMETDiscretizer.h:2152
GradientModel< T >::getGradientMatrix
static GradMatrixType & getGradientMatrix(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:454
COMETDiscretizer::COMETSolveFine
void COMETSolveFine(const int dir, const int level)
Definition: COMETDiscretizer.h:91
COMETDiscretizer::COMETEquilibrium
void COMETEquilibrium(const int cell, TMatrix *Amat, TArray &BVec)
Definition: COMETDiscretizer.h:1026
inverse
SquareMatrix< T, 2 > inverse(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:118
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Array
Definition: Array.h:14
pmode::getfield
Field & getfield()
Definition: pmode.h:74
StorageSite::getOffset
int getOffset() const
Definition: StorageSite.h:87
COMETDiscretizer::COMETConvectionCoarse
void COMETConvectionCoarse(const int cell0, TArrow &Amat, TArray &BVec)
Definition: COMETDiscretizer.h:819
COMETDiscretizer::GradModelType
GradientModel< T > GradModelType
Definition: COMETDiscretizer.h:59
COMETDiscretizer::_geomFields
const GeomFields & _geomFields
Definition: COMETDiscretizer.h:2147
Kspace::ScatterPhonons
void ScatterPhonons(const int c, const int totIts, TArray &C, TArray &B, const TArray &V, TArray &newE, const T cv)
Definition: Kspace.h:1476
SquareMatrix::Solve
void Solve(TArray &bVec)
Definition: SquareMatrix.h:34
COMETDiscretizer::COMETSolveCoarse
void COMETSolveCoarse(const int dir, const int level)
Definition: COMETDiscretizer.h:244
ArrowHeadMatrix::Solve
void Solve(XArray &bVec)
Definition: esbgkbase/ArrowHeadMatrix.h:71
pmode::calcde0dT
T calcde0dT(T Tl)
Definition: pmode.h:126
Kspace::gettotmodes
int gettotmodes()
Definition: Kspace.h:393
Kspace::addSource
void addSource(const int c, TArray &BVec, const T cv)
Definition: Kspace.h:1693
COMETDiscretizer::_faceArea
const VectorT3Array & _faceArea
Definition: COMETDiscretizer.h:2153
COMETDiscretizer::FgTKClistMap
map< int, TKClist > FgTKClistMap
Definition: COMETDiscretizer.h:56
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
pmode::getresid
Field & getresid()
Definition: pmode.h:78
COMETDiscretizer::updateGhostFine
void updateGhostFine(const int cell, const GradMatrix &gMat)
Definition: COMETDiscretizer.h:1655
COMETDiscretizer::GradType
Gradient< T > GradType
Definition: COMETDiscretizer.h:57
Vector::zero
void zero()
Definition: Vector.h:156
KSConnectivity::multiplySelf
void multiplySelf(const TArray &x, TArray &b, const T scale) const
Definition: KSConnectivity.h:80
kvol::getkvec
Tvec getkvec()
Definition: kvol.h:39
pmode
Definition: pmode.h:18
COMETDiscretizer::_cells
const StorageSite & _cells
Definition: COMETDiscretizer.h:2148
Kspace::getSourceTerm
void getSourceTerm(const int c, TArray &s, TArray &ds)
Definition: Kspace.h:1468
kvol
Definition: kvol.h:14
COMETDiscretizer::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: COMETDiscretizer.h:40
COMETDiscretizer::GradMatrix
GradModelType::GradMatrixType GradMatrix
Definition: COMETDiscretizer.h:60
COMETDiscretizer::TKClist
vector< TKCptr > TKClist
Definition: COMETDiscretizer.h:55
COMETDiscretizer::FaceToFg
map< int, VecInt2 > FaceToFg
Definition: COMETDiscretizer.h:50
COMETDiscretizer::updateGhostCoarse
void updateGhostCoarse(const int cell)
Definition: COMETDiscretizer.h:1932
COMETDiscretizer::COMETShifted
void COMETShifted(const int cell, TMatrix *Amat, TArray &BVec)
Definition: COMETDiscretizer.h:1153
COMETDiscretizer::_FaceToKSC
const FgTKClistMap & _FaceToKSC
Definition: COMETDiscretizer.h:2165
kvol::getdk3
T getdk3()
Definition: kvol.h:42
pmode::geteShifted
Field & geteShifted()
Definition: pmode.h:75
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
Array::getLength
int getLength() const
Definition: Array.h:87
COMETDiscretizer::findResidCoarse
void findResidCoarse(const bool plusFAS)
Definition: COMETDiscretizer.h:1314
COMETDiscretizer::findFgId
int findFgId(const int faceIndex)
Definition: COMETDiscretizer.h:1563
COMETDiscretizer::Tkspace
Kspace< T > Tkspace
Definition: COMETDiscretizer.h:36