Memosa-FVM  0.2
ScatteringKernel.h
Go to the documentation of this file.
1 #ifndef _SCATTERINGKERNEL_H_
2 #define _SCATTERINGKERNEL_H_
3 
4 #include "Array.h"
5 #include "KSConnectivity.h"
6 #include <iostream>
7 #include <fstream>
8 #include <stdio.h>
9 #include <math.h>
10 
11 template<class T> class Kspace;
12 
13 template <class T>
14 class ScatteringKernel
15 {
16 
17  typedef Array<T> TArray;
18  typedef Array<int> IntArray;
19  typedef Array<bool> BArray;
20  typedef Vector<T,3> VectorT3;
21  typedef Array<VectorT3> VectorT3Array;
22  typedef kvol<T> Tkvol;
23  typedef pmode<T> Tmode;
24  typedef KSConnectivity<T> Tksconn;
25 
26  public:
27  ScatteringKernel(Kspace<T>& kspace) :
28  _kspace(kspace),
29  _type1Collisions(kspace.gettotmodes(),kspace.gettotmodes()),
30  _type2Collisions(kspace.gettotmodes(),kspace.gettotmodes()),
31  _maxPhi(0),
32  _maxDkl(0)
33  {_kspace.setScattKernel(*this);}
34 
35  void ReadType1(const string& NamePhonon2, const string& NamePhonon3, const T tol)
36  {
37  int dum(0);
38  const T hbar=6.582119e-16; // (eV s)
39  const T kb=8.617343e-5; // (eV/K)
40 
41  cout<<"Reading type I collisions..."<<endl;
42  /*
43  ifstream fp_in2;
44  ifstream fp_in3;
45  fp_in2.open(NamePhonon2,ios::in | ios::binary);
46  fp_in3.open(NamePhonon3,ios::in | ios::binary);
47  */
48 
49  FILE* fp_in2(NULL);
50  FILE* fp_in3(NULL);
51 
52  if(fp_in3)
53  fclose(fp_in3);
54  if(fp_in2)
55  fclose(fp_in2);
56 
57  int numAtt(10);
58  while(fp_in2==NULL && numAtt>0)
59  {
60  fp_in2=fopen(NamePhonon2.c_str(),"rb");
61  if(fp_in2==NULL)
62  cout<<"Error reading 2 "<<numAtt<<endl;
63  numAtt--;
64  }
65 
66  numAtt=10;
67  while(fp_in3==NULL && numAtt>0)
68  {
69  fp_in3=fopen(NamePhonon3.c_str(),"rb");
70  if(fp_in3==NULL)
71  cout<<"Error reading 3 "<<numAtt<<endl;
72  numAtt--;
73  }
74 
75  const int rowLen=_type1Collisions.getSelfSize();
76 
77  int row2(-1);
78  int nnz2(-1);
79  int row3(-1);
80  int nnz3(-1);
81  int index2(-1);
82  int index3(-1);
83  T dkl(-1);
84  T phi(-1);
85 
86  _type1Collisions.initSelfCount();
87  _type1Collisions.initOtherCount();
88 
89  TArray& w(_kspace.getFreqArray());
90 
91  cout<<"Counting I interactions..."<<endl;
92 
93  for(int i=0;i<rowLen;i++)
94  {
95  //if(i % 2000==0)
96  //cout<<"Row: "<<i<<endl;
97 
98  dum=fread(&row2,sizeof(int),1,fp_in2);
99  dum=fread(&nnz2,sizeof(int),1,fp_in2);
100  dum=fread(&row3,sizeof(int),1,fp_in3);
101  dum=fread(&nnz3,sizeof(int),1,fp_in3);
102  //const T n1=1./(exp(hbar*w[row2]/kb/300.)-1);
103 
104  for(int j=0;j<nnz2;j++)
105  {
106  dum=fread(&index2,sizeof(int),1,fp_in2);
107  dum=fread(&dkl,sizeof(T),1,fp_in2);
108  dum=fread(&index3,sizeof(int),1,fp_in3);
109  dum=fread(&phi,sizeof(T),1,fp_in3);
110 
111  //const T n2=1./(exp(hbar*w[index2]/kb/300.)-1);
112  //const T n3=1./(exp(hbar*w[index3]/kb/300.)-1);
113  //const T err=fabs((n1+1)*(n2+1)*n3-n1*n2*(n3+1))/n1;
114  const T err=fabs(w[row2]+w[index2]-w[index3])/w[row2];
115 
116  const int m1=row2%6;
117  const int m2=index2%6;
118  const int m3=index3%6;
119 
120  int zCnt(0);
121 
122 
123  if((m1==0)||(m1==3))
124  zCnt++;
125  if((m2==0)||(m2==3))
126  zCnt++;
127  if((m3==0)||(m3==3))
128  zCnt++;
129 
130 
131  //if((m1==0)||(m2==0)||(m3==0))
132  // zCnt=1;
133 
134  //if(w[row2]<.25e14)
135  // zCnt=1;
136 
137  if(phi<0)
138  {
139  cout<<row2<<" "<<index2<<endl;
140  throw CException("I: negative value for phi!.");
141  }
142 
143  if(dkl<0)
144  {
145  cout<<row2<<" "<<index2<<endl;
146  throw CException("I: negative value for dkl!.");
147  }
148 
149  //zCnt=0;
150  if(err<tol && zCnt!=1 && zCnt!=3 && dkl>0.)
151  {
152  _type1Collisions.addCountSelf(row2,1);
153  _type1Collisions.addCountOther(row3,1);
154  if(phi>_maxPhi)
155  _maxPhi=phi;
156  if(dkl>_maxDkl)
157  _maxDkl=dkl;
158  }
159 
160  }
161  }
162 
163  _type1Collisions.finishCountSelf();
164  _type1Collisions.finishCountOther();
165 
166  rewind(fp_in2);
167  rewind(fp_in3);
168 
169  cout<<"Entering I values..."<<endl;
170 
171  for(int i=0;i<rowLen;i++)
172  {
173  //if(i % 2000==0)
174  //cout<<"Row: "<<i<<endl;
175 
176  dum=fread(&row2,sizeof(int),1,fp_in2);
177  dum=fread(&nnz2,sizeof(int),1,fp_in2);
178  dum=fread(&row3,sizeof(int),1,fp_in3);
179  dum=fread(&nnz3,sizeof(int),1,fp_in3);
180  //const T n1=1./(exp(hbar*w[row2]/kb/300.)-1);
181 
182  for(int j=0;j<nnz2;j++)
183  {
184  dum=fread(&index2,sizeof(int),1,fp_in2);
185  dum=fread(&dkl,sizeof(T),1,fp_in2);
186  dum=fread(&index3,sizeof(int),1,fp_in3);
187  dum=fread(&phi,sizeof(T),1,fp_in3);
188 
189  //const T n2=1./(exp(hbar*w[index2]/kb/300.)-1);
190  //const T n3=1./(exp(hbar*w[index3]/kb/300.)-1);
191  //const T err=fabs((n1+1)*(n2+1)*n3-n1*n2*(n3+1))/n1;
192  const T err=fabs(w[row2]+w[index2]-w[index3])/w[row2];
193 
194  const int m1=row2%6;
195  const int m2=index2%6;
196  const int m3=index3%6;
197 
198  int zCnt(0);
199 
200  if((m1==0)||(m1==3))
201  zCnt++;
202  if((m2==0)||(m2==3))
203  zCnt++;
204  if((m3==0)||(m3==3))
205  zCnt++;
206 
207  //if((m1==0)||(m2==0)||(m3==0))
208  // zCnt=1;
209 
210  //if(w[row2]<.25e14)
211  // zCnt=1;
212 
213  //zCnt=0;
214  if(err<tol && zCnt!=1 && zCnt!=3 && dkl>0.)
215  {
216  _type1Collisions.addSelf(row2,index2,dkl);
217  _type1Collisions.addOther(row3,index3,phi);
218  }
219  }
220  }
221 
222  _type1Collisions.finishAddSelf();
223  _type1Collisions.finishAddOther();
224 
225  fclose(fp_in3);
226  fclose(fp_in2);
227 
228  cout<<"Type I complete: "<<_type1Collisions.getSelfNNZ()<<endl;
229 
230  }
231 
232  void ReadType2(const string& NamePhonon2, const string& NamePhonon3, const T tol)
233  {
234 
235  cout<<"Reading type II collisions..."<<endl;
236  const T hbar=6.582119e-16; // (eV s)
237  const T kb=8.617343e-5; // (eV/K)
238 
239  FILE* fp_in2(NULL);
240  FILE* fp_in3(NULL);
241 
242  if(fp_in3)
243  fclose(fp_in3);
244  if(fp_in2)
245  fclose(fp_in2);
246 
247  int numAtt(10);
248  while(fp_in2==NULL && numAtt>0)
249  {
250  fp_in2=fopen(NamePhonon2.c_str(),"rb");
251  if(fp_in2==NULL)
252  cout<<"Error reading 2 "<<numAtt<<endl;
253  numAtt--;
254  }
255 
256  numAtt=10;
257  while(fp_in3==NULL && numAtt>0)
258  {
259  fp_in3=fopen(NamePhonon3.c_str(),"rb");
260  if(fp_in3==NULL)
261  cout<<"Error reading 3 "<<numAtt<<endl;
262  numAtt--;
263  }
264 
265  TArray& w(_kspace.getFreqArray());
266  const int rowLen=_type2Collisions.getSelfSize();
267 
268  int row2(-1);
269  int nnz2(-1);
270  int row3(-1);
271  int nnz3(-1);
272  int index2(-1);
273  int index3(-1);
274  T dkl(-1);
275  T phi(-1);
276  int dum(0);
277 
278  _type2Collisions.initSelfCount();
279  _type2Collisions.initOtherCount();
280 
281  cout<<"Counting II interactions..."<<endl;
282 
283  for(int i=0;i<rowLen;i++)
284  {
285  //if(i % 2000==0)
286  //cout<<"Row: "<<i<<endl;
287 
288  dum=fread(&row2,sizeof(int),1,fp_in2);
289  dum=fread(&nnz2,sizeof(int),1,fp_in2);
290  dum=fread(&row3,sizeof(int),1,fp_in3);
291  dum=fread(&nnz3,sizeof(int),1,fp_in3);
292  //const T n1=1./(exp(hbar*w[row2]/kb/300.)-1);
293 
294  for(int j=0;j<nnz2;j++)
295  {
296  dum=fread(&index2,sizeof(int),1,fp_in2);
297  dum=fread(&dkl,sizeof(T),1,fp_in2);
298  dum=fread(&index3,sizeof(int),1,fp_in3);
299  dum=fread(&phi,sizeof(T),1,fp_in3);
300 
301  //const T n2=1./(exp(hbar*w[index2]/kb/300.)-1);
302  //const T n3=1./(exp(hbar*w[index3]/kb/300.)-1);
303  //const T err=fabs(n3*n2*(n1+1)-n1*(n2+1)*(n3+1))/n1;
304  const T err=fabs(w[row2]-w[index2]-w[index3])/w[index2];
305 
306  const int m1=row2%6;
307  const int m2=index2%6;
308  const int m3=index3%6;
309 
310  int zCnt(0);
311 
312  if((m1==0)||(m1==3))
313  zCnt++;
314  if((m2==0)||(m2==3))
315  zCnt++;
316  if((m3==0)||(m3==3))
317  zCnt++;
318 
319  if(phi<0)
320  {
321  cout<<row2<<" "<<index2<<endl;
322  throw CException("II: negative value for phi!.");
323  }
324 
325  if(dkl<0)
326  {
327  cout<<row2<<" "<<index2<<endl;
328  throw CException("II: negative value for dkl!.");
329  }
330 
331  //if((m1==0)||(m2==0)||(m3==0))
332  // zCnt=1;
333 
334  //if(w[row2]<.25e14)
335  // zCnt=1;
336 
337  //zCnt=0;
338  if(zCnt!=1 && zCnt!=3 && dkl>0. && err<tol)
339  {
340  _type2Collisions.addCountSelf(row2,1);
341  _type2Collisions.addCountOther(row3,1);
342  if(phi>_maxPhi)
343  _maxPhi=phi;
344  if(dkl>_maxDkl)
345  _maxDkl=dkl;
346  }
347  }
348  }
349 
350  _type2Collisions.finishCountSelf();
351  _type2Collisions.finishCountOther();
352 
353  rewind(fp_in2);
354  rewind(fp_in3);
355 
356  cout<<"Entering II values..."<<endl;
357 
358  for(int i=0;i<rowLen;i++)
359  {
360  //if(i % 2000==0)
361  //cout<<"Row: "<<i<<endl;
362 
363  dum=fread(&row2,sizeof(int),1,fp_in2);
364  dum=fread(&nnz2,sizeof(int),1,fp_in2);
365  dum=fread(&row3,sizeof(int),1,fp_in3);
366  dum=fread(&nnz3,sizeof(int),1,fp_in3);
367  //const T n1=1./(exp(hbar*w[row2]/kb/300.)-1);
368 
369  for(int j=0;j<nnz2;j++)
370  {
371  dum=fread(&index2,sizeof(int),1,fp_in2);
372  dum=fread(&dkl,sizeof(T),1,fp_in2);
373  dum=fread(&index3,sizeof(int),1,fp_in3);
374  dum=fread(&phi,sizeof(T),1,fp_in3);
375 
376  //const T n2=1./(exp(hbar*w[index2]/kb/300.)-1);
377  //const T n3=1./(exp(hbar*w[index3]/kb/300.)-1);
378  //const T err=fabs(n3*n2*(n1+1)-n1*(n2+1)*(n3+1))/n1;
379  const T err=fabs(w[row2]-w[index2]-w[index3])/w[index2];
380 
381  const int m1=row2%6;
382  const int m2=index2%6;
383  const int m3=index3%6;
384 
385  int zCnt(0);
386 
387  if((m1==0)||(m1==3))
388  zCnt++;
389  if((m2==0)||(m2==3))
390  zCnt++;
391  if((m3==0)||(m3==3))
392  zCnt++;
393 
394  //if((m1==0)||(m2==0)||(m3==0))
395  // zCnt=1;
396 
397  //if(w[row2]<.25e14)
398  // zCnt=1;
399 
400  //zCnt=0;
401  if(zCnt!=1 && zCnt!=3 && dkl>0. && err<tol)
402  {
403  _type2Collisions.addSelf(row2,index2,dkl);
404  _type2Collisions.addOther(row3,index3,phi);
405  }
406  }
407  }
408 
409  _type2Collisions.finishAddSelf();
410  _type2Collisions.finishAddOther();
411  fclose(fp_in3);
412  fclose(fp_in2);
413 
414  cout<<"Type II complete: "<<_type2Collisions.getSelfNNZ()<<endl;
415 
416  normalize();
417 
418  }
419 
420  void addFreqs()
421  {
422  //const T hbarJoule=1.054571726e-34;
423  TArray& w(_kspace.getFreqArray());
424  const int Rows=w.getLength();
425 
426  const IntArray& t1p2row=_type1Collisions.getSelfRow();
427  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
428  const IntArray& t1p2col=_type1Collisions.getSelfCol();
429  const IntArray& t1p3col=_type1Collisions.getOtherCol();
430  const IntArray& t2p2row=_type2Collisions.getSelfRow();
431  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
432  const IntArray& t2p2col=_type2Collisions.getSelfCol();
433  const IntArray& t2p3col=_type2Collisions.getOtherCol();
434 
435  TArray& t1phi=_type1Collisions.getNonConstOtherCoeffs();
436  TArray& t2phi=_type2Collisions.getNonConstOtherCoeffs();
437 
438  //type 1 collisions
439 
440  for(int i=0;i<Rows;i++)
441  {
442  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
443  {
444  const int j2=t1p2col[pos];
445  const int j3=t1p3col[pos];
446  T& phi=t1phi[pos];
447  phi=phi/w[i]/w[j2]/w[j3];
448  }
449  }
450 
451  //type 2 collisions
452  for(int i=0;i<Rows;i++)
453  {
454  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
455  {
456  const int j2=t2p2col[pos];
457  const int j3=t2p3col[pos];
458  T& phi=t2phi[pos];
459  //const T w3=w[i]-w[j2];
460  phi=phi/w[i]/w[j2]/w[j3];
461  }
462  }
463 
464  }
465 
466  void updateSourceTerm(const TArray& e, const TArray& w, TArray& S)
467  {
468 
469  const T hbarJoule=1.054571726e-34;
470  const T JouleToeV=6.24150974e18;
471  const T Acell=5.378395621705545e-20;
472 
473 
474  S.zero();
475  const int Rows=S.getLength();
476  const IntArray& t1p2row=_type1Collisions.getSelfRow();
477  const IntArray& t1p3row=_type1Collisions.getOtherRow();
478  const IntArray& t1p2col=_type1Collisions.getSelfCol();
479  const IntArray& t1p3col=_type1Collisions.getOtherCol();
480  const IntArray& t2p2row=_type2Collisions.getSelfRow();
481  const IntArray& t2p3row=_type2Collisions.getOtherRow();
482  const IntArray& t2p2col=_type2Collisions.getSelfCol();
483  const IntArray& t2p3col=_type2Collisions.getOtherCol();
484 
485  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
486  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
487  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
488  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
489 
490  //type 1 collisions
491  for(int i=0;i<Rows;i++)
492  {
493  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
494  {
495  const int j2=t1p2col[pos];
496  const int j3=t1p3col[pos];
497  const T dkl=t1dkl[pos];
498  const T phi=t1phi[pos];
499  const T n1=e[i]/w[i];
500  const T n2=e[j2]/w[j2];
501  const T n3=e[j3]/w[j3];
502  S[i]+=dkl*phi*(n1*(n3-n2)+n3*(n2+1));
503  }
504  }
505 
506  //type 2 collisions
507  for(int i=0;i<Rows;i++)
508  {
509  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
510  {
511  const int j2=t2p2col[pos];
512  const int j3=t2p3col[pos];
513  const T dkl=t2dkl[pos];
514  const T phi=t2phi[pos];
515  const T n1=e[i]/w[i];
516  const T n2=e[j2]/w[j2];
517  const T n3=e[j3]/w[j3];
518  S[i]+=dkl*phi*(n1*(n3-n2)+n3*(n2+1));
519  }
520  }
521 
522  const T preFac=Acell/6.283185307/hbarJoule*JouleToeV;
523  for(int i=0;i<Rows;i++)
524  S[i]*=preFac*w[i];
525 
526  }
527 
528  void updateSourceTermTest(const T Tl)
529  {
530 
531  const T hbarJoule=1.054571726e-34;
532  const T hbar=6.582119e-16;
533  //const T JouleToeV=6.24150974e18;
534  const T Acell=5.378395621705545e-20;
535  const T cv=1e-9*1e-9;
536  T temp=Tl;
537 
538  TArray S(_kspace.gettotmodes());
539  S.zero();
540  TArray e(_kspace.gettotmodes());
541  e.zero();
542  _kspace.getEquilibriumArray(e,Tl);
543  TArray& w(_kspace.getFreqArray());
544 
545  const int Rows=S.getLength();
546  const IntArray& t1p2row=_type1Collisions.getSelfRow();
547  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
548  const IntArray& t1p2col=_type1Collisions.getSelfCol();
549  const IntArray& t1p3col=_type1Collisions.getOtherCol();
550  const IntArray& t2p2row=_type2Collisions.getSelfRow();
551  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
552  const IntArray& t2p2col=_type2Collisions.getSelfCol();
553  const IntArray& t2p3col=_type2Collisions.getOtherCol();
554 
555  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
556  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
557  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
558  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
559 
560  //type 1 collisions
561  for(int i=0;i<Rows;i++)
562  {
563  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
564  {
565  const int j2=t1p2col[pos];
566  const int j3=t1p3col[pos];
567  const T dkl=t1dkl[pos];
568  const T phi=t1phi[pos]*Acell;
569  const T n1=e[i]/hbar/w[i];
570  const T n2=e[j2]/hbar/w[j2];
571  const T n3=e[j3]/hbar/w[j3];
572  S[i]+=dkl*phi*((n1+1)*(n2+1)*n3-n1*n2*(n3+1));
573  }
574  }
575 
576  //type 2 collisions
577  for(int i=0;i<Rows;i++)
578  {
579  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
580  {
581  const int j2=t2p2col[pos];
582  const int j3=t2p3col[pos];
583  const T dkl=t2dkl[pos];
584  const T phi=t2phi[pos]*Acell;
585  const T n1=e[i]/w[i]/hbar;
586  const T n2=e[j2]/w[j2]/hbar;
587  const T n3=e[j3]/w[j3]/hbar;
588  S[i]+=0.5*dkl*phi*(n3*n2*(n1+1)-n1*(n2+1)*(n3+1));
589  }
590  }
591 
592  const T preFac=1./2./3.141592653/pow(hbarJoule,2)*_maxDkl*_maxPhi*cv;
593  for(int i=0;i<Rows;i++)
594  S[i]*=preFac*w[i]*hbar;
595 
596  T defect(0);
597  for(int i=0;i<Rows;i++)
598  defect+=S[i];
599 
600  for(int i=0;i<Rows;i++)
601  e[i]+=S[i];
602 
603  int cnt(0);
604  T esum(0);
605  const int klen=_kspace.getlength();
606  for(int k=0;k<klen;k++)
607  {
608  Tkvol& kv=_kspace.getkvol(k);
609  const int modenum=kv.getmodenum();
610  T dk3=kv.getdk3();
611  for(int m=0;m<modenum;m++)
612  {
613  esum+=e[cnt]*dk3;
614  cnt++;
615  }
616  }
617 
618  _kspace.calcTemp(temp,esum);
619  cout<<"New Temp: "<<temp<<endl;
620 
621  }
622 
623  void updateSource(const int c, TArray& S)
624  {
625 
626  const T hbarJoule=1.054571726e-34;
627  const T hbar=6.582119e-16;
628  const T Acell=5.378395621705545e-20;
629  const T kb=8.617343e-5; // (eV/K)
630  T cOt(0);
631  S.zero();
632  const int Rows=S.getLength();
633  TArray e(Rows);
634  //TArray e0(Rows);
635  _kspace.geteCellVals(c,e);
636  //_kspace.gete0CellVals(c,e0);
637  const TArray& w(_kspace.getFreqArray());
638  const T Tlat(300.05);//Tl(_kspace.calcLatTemp(c));
639  const T Tl(_kspace.calcLatTemp(c));
640  TArray e0(Rows);
641  //_kspace.gete0CellVals(c,e0);
642  _kspace.getEquilibriumArray(e0,Tl);
643 
644  const IntArray& t1p2row=_type1Collisions.getSelfRow();
645  const IntArray& t1p3row=_type1Collisions.getOtherRow();
646  const IntArray& t1p2col=_type1Collisions.getSelfCol();
647  const IntArray& t1p3col=_type1Collisions.getOtherCol();
648  const IntArray& t2p2row=_type2Collisions.getSelfRow();
649  const IntArray& t2p3row=_type2Collisions.getOtherRow();
650  const IntArray& t2p2col=_type2Collisions.getSelfCol();
651  const IntArray& t2p3col=_type2Collisions.getOtherCol();
652 
653  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
654  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
655  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
656  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
657 
658  //type 1 collisions
659  for(int i=0;i<Rows;i++)
660  {
661  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
662  {
663  const int j2=t1p2col[pos];
664  const int j3=t1p3col[pos];
665  const T dkl=t1dkl[pos];
666  const T phi=t1phi[pos];
667  const T w3=w[i]+w[j2];
668  const T w1=w[j3]-w[j2];
669  const T n1=e[i]/w[i]/hbar;
670  const T n2=e[j2]/w[j2]/hbar;
671  const T n3=e[j3]/w[j3]/hbar;
672  const T n10=e0[i]/w[i]/hbar;
673  const T n20=e0[j2]/w[j2]/hbar;
674  const T n30=e0[j3]/w[j3]/hbar;
675  const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
676  const T dn1=n1-n10;
677  const T dn2=n2-n20;
678  const T dn3=(n3-n30);//*n012/n30;
679  T Large(n10);
680  T e032(0);
681  T frac(0);
682  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
683 
684  if(w1>0.)
685  {
686  e032=hbar*w1/(exp(hbar*w1/kb/Tl)-1.);
687  frac=e032/e0[i];
688  }
689 
690  if(n20>Large)
691  Large=n20;
692  if(n30>Large)
693  Large=n30;
694 
695  const T nsum=Large*(-dn1*dn2/Large+dn3/Large+dn1*dn3/Large
696  +dn2*dn3/Large-dn2*n10/Large+dn3*n10/Large-dn1*n20/Large
697  +dn3*n20/Large+dn1*n30/Large+dn2*n30/Large);
698  /*
699  if(j2>i)
700  {
701  S[i]+=dkl*phi*nsum;
702  S[j2]+=dkl*phi*nsum;
703  S[j3]-=0.5*dkl*phi*nsum;
704  //S[i]+=dkl*phi*((n1+1)*(n2+1)*n3-n1*n2*(n3+1));
705  }*/
706 
707  S[i]+=dkl*phi*nsum;
708  }
709  }
710 
711  T maxS(0);
712 
713  //type 2 collisions
714 
715  for(int i=0;i<Rows;i++)
716  {
717  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
718  {
719  const int j2=t2p2col[pos];
720  const int j3=t2p3col[pos];
721  const T dkl=t2dkl[pos];
722  const T phi=t2phi[pos];
723  const T w3=w[i]+w[j2];
724  const T w1=w[j3]-w[j2];
725  const T n1=e[i]/w[i]/hbar;
726  const T n2=e[j2]/w[j2]/hbar;
727  const T n3=e[j3]/w[j3]/hbar;
728  const T n10=e0[i]/w[i]/hbar;
729  const T n20=e0[j2]/w[j2]/hbar;
730  const T n30=e0[j3]/w[j3]/hbar;
731  const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
732  const T dn1=n1-n10;
733  const T dn2=n2-n20;
734  const T dn3=(n3-n30);//*n012/n30;
735  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
736  T e032(0);
737  T frac(0);
738 
739  if(w1>0.)
740  {
741  e032=hbar*w1/(exp(hbar*w1/kb/Tl)-1.);
742  frac=e032/e0[i];
743  }
744 
745  T Large(n10);
746 
747  if(n20>Large)
748  Large=n20;
749  if(n30>Large)
750  Large=n30;
751 
752  const T nsum=Large*(-dn1/Large-dn1*dn2/Large-dn1*dn3/Large
753  +dn2*dn3/Large-dn2*n10/Large-dn3*n10/Large
754  -dn1*n20/Large+dn3*n20/Large-dn1*n30/Large
755  +dn2*n30/Large);
756 
757  S[i]+=0.5*dkl*phi*nsum;
758  //S[i]+=0.5*dkl*phi*(n3*n2*(n1+1)-n1*(n2+1)*(n3+1));
759  }
760  }
761 
762  const T preFac=Acell/16./3.141592653*hbarJoule*_maxDkl*_maxPhi;
763  const int klen=_kspace.getlength();
764  for(int i=0;i<Rows;i++)
765  {
766  S[i]*=(preFac*w[i]*hbar);
767  if(fabs(S[i])>fabs(maxS))
768  maxS=fabs(S[i]);
769  }
770 
771  const T DK3=_kspace.getDK3();
772  T defect(0);
773  int cnt(0);
774  T eq(0);
775  const int cStart(_kspace.getGlobalIndex(c,0));
776  int ind(cStart);
777 
778  for(int k=0;k<klen;k++)
779  {
780  Tkvol& kv=_kspace.getkvol(k);
781  const int modenum=kv.getmodenum();
782  T dk3=kv.getdk3();
783  for(int m=0;m<modenum;m++)
784  {
785  defect+=((S[cnt]/maxS)*(dk3/DK3));
786  Tmode& mode=kv.getmode(m);
787  const T tau=mode.gettau();
788  cOt+=mode.calce0(Tl)*(dk3/DK3);
789  //const T e0=mode.calce0(Tl);
790  //eq+=(e0-e[ind])/tau*(dk3/DK3);
791  ind++;
792  cnt++;
793  }
794  }
795 
796  cOt*=DK3;
797  T f(maxS*defect/eq);
798  defect*=maxS*DK3;
799  //cout<<"Defect1: "<<defect<<endl;
800 
801  if(fabs(defect)>0.)
802  {
803  cnt=0;
804  T d(defect);
805  defect=0.;
806  ind=cStart;
807  for(int k=0;k<klen;k++)
808  {
809  Tkvol& kv=_kspace.getkvol(k);
810  const int modenum=kv.getmodenum();
811  T dk3=kv.getdk3();
812  for(int m=0;m<modenum;m++)
813  {
814  Tmode& mode=kv.getmode(m);
815  const T tau=mode.gettau();
816  const T cp=mode.calcde0dT(Tl);
817  const T e0(mode.calce0(Tl));
818  //S[cnt]+=f*(e0-e[ind])/tau;
819  S[cnt]=d*(S[cnt]/d-e0/cOt);
820  defect+=((S[cnt]/maxS)*(dk3/DK3));
821  ind++;
822  cnt++;
823  }
824  }
825  defect*=DK3*maxS;
826  //cout<<"Defect2: "<<defect<<endl;
827  }
828  else
829  {
830  S.zero();
831  }
832 
833  }
834 
835  void getTypeIsource(const int c, TArray& S, TArray& dS, const bool correct)
836  {
837 
838  const T hbarJoule=1.054571726e-34;
839  const T hbar=6.582119e-16;
840  const T Acell=5.378395621705545e-20;
841  const T kb=8.617343e-5; // (eV/K)
842  T cOt(0);
843  S.zero();
844  dS.zero();
845  const int Rows=S.getLength();
846  TArray Sl(Rows);
847  TArray dSl(Rows);
848  Sl.zero();
849  dSl.zero();
850  TArray e(Rows);
851  //TArray e0(Rows);
852  _kspace.geteCellVals(c,e);
853  //_kspace.gete0CellVals(c,e0);
854  const TArray& w(_kspace.getFreqArray());
855  //const T Tlat(300.0);//Tl(_kspace.calcLatTemp(c));
856  const T Tl(_kspace.calcLatTemp(c));
857  TArray e0(Rows);
858  //_kspace.gete0CellVals(c,e0);
859  _kspace.getEquilibriumArray(e0,Tl);
860  const T DK3=_kspace.getDK3();
861 
862  const IntArray& t1p2row=_type1Collisions.getSelfRow();
863  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
864  const IntArray& t1p2col=_type1Collisions.getSelfCol();
865  const IntArray& t1p3col=_type1Collisions.getOtherCol();
866  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
867  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
868 
869  //type 1 collisions
870 
871  for(int i=0;i<Rows;i++)
872  {
873  const int k1=floor(i/6);
874  const T dk31=_kspace.getkvol(k1).getdk3();
875  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
876  {
877  const int j2=t1p2col[pos];
878  const int k2=floor(j2/6);
879  const T dk32=_kspace.getkvol(k2).getdk3();
880  const int j3=t1p3col[pos];
881  const int k3=floor(j3/6);
882  const T dk33=_kspace.getkvol(k3).getdk3();
883  const T dkl=t1dkl[pos];
884  const T phi=t1phi[pos];
885  const T w3=w[i]+w[j2];
886  //const T w1=w[j3]-w[j2];
887  const T n1=e[i]/w[i]/hbar;
888  const T n2=e[j2]/w[j2]/hbar;
889  const T n3=e[j3]/w[j3]/hbar;
890  const T n10=e0[i]/w[i]/hbar;
891  const T n20=e0[j2]/w[j2]/hbar;
892  const T n30=e0[j3]/w[j3]/hbar;
893  //const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
894  const T dn1=n1-n10;
895  const T dn2=n2-n20;
896  const T dn3=(n3-n30);//*n012/n30;
897  T Large(n10);
898  //T e032(0);
899  //T frac(0);
900  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
901 
902  if(n20>Large)
903  Large=n20;
904  if(n30>Large)
905  Large=n30;
906 
907  const T nsum=(-dn1*dn2+dn3+dn1*dn3
908  +dn2*dn3-dn2*n10+dn3*n10-dn1*n20
909  +dn3*n20+dn1*n30+dn2*n30);
910 
911  //const T nsum1=(-dn1*n20+dn1*n30);
912  //const T nsum1=(n1+1)*(n20+1)*n30-n1*n20*(n30+1);
913 
914 
915  Sl[i]+=dkl*phi*nsum*w[i];
916  Sl[j2]+=dkl*phi*nsum*w[j2]*dk31/dk32;
917  Sl[j3]-=dkl*phi*nsum*w3*dk31/dk33;
918  //S[i]+=dkl*phi*((n1+1)*(n2+1)*n3-n1*n2*(n3+1));
919 
920  }
921  }
922 
923  T maxS(0);
924  const T preFac=Acell/16./3.141592653*hbarJoule*_maxPhi*hbar*_maxDkl;
925  const int klen=_kspace.getlength();
926  for(int i=0;i<Rows;i++)
927  {
928  Sl[i]*=preFac;//*w[i]*hbar);
929  dSl[i]*=preFac;
930  if(fabs(Sl[i])>fabs(maxS))
931  maxS=fabs(Sl[i]);
932  }
933 
934  long double defect(0);
935  int cnt(0);
936  T eq(0);
937  const int cStart(_kspace.getGlobalIndex(c,0));
938  int ind(cStart);
939 
940  for(int k=0;k<klen;k++)
941  {
942  Tkvol& kv=_kspace.getkvol(k);
943  const int modenum=kv.getmodenum();
944  T dk3=kv.getdk3();
945  for(int m=0;m<modenum;m++)
946  {
947  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
948  //if(fabs(Sl[cnt])*(dk3/DK3)>100)
949  // Sl[cnt]=0.;
950  Tmode& mode=kv.getmode(m);
951  //const T tau=mode.gettau();
952  cOt+=mode.calce0(Tl)*(dk3/DK3);
953  //const T e0=mode.calce0(Tl);
954  //eq+=(e0-e[ind])/tau*(dk3/DK3);
955  ind++;
956  cnt++;
957  }
958  }
959 
960  cOt*=DK3;
961  //T f(maxS*defect/eq);
962  defect*=maxS*DK3;
963  //cout<<"Defect1: "<<defect<<endl;
964  //cout<<"Temperature: "<<Tl<<endl;
965 
966 
967  if(correct)
968  {
969  cnt=0;
970  //T d(defect);
971  defect=0.;
972  ind=cStart;
973  for(int k=0;k<klen;k++)
974  {
975  Tkvol& kv=_kspace.getkvol(k);
976  const int modenum=kv.getmodenum();
977  T dk3=kv.getdk3();
978  for(int m=0;m<modenum;m++)
979  {
980  Tmode& mode=kv.getmode(m);
981  //const T tau=mode.gettau();
982  //const T cp=mode.calcde0dT(Tl);
983  //const T e0(mode.calce0(Tl));
984  //S[cnt]+=f*(e0-e[ind])/tau;
985  //Sl[cnt]=d*(Sl[cnt]/d-e0/cOt);
986  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
987  ind++;
988  cnt++;
989  }
990  }
991  defect*=DK3*maxS;
992  //cout<<"Defect2: "<<defect<<endl;
993  }
994 
995  for(int i=0;i<Rows;i++)
996  {
997  S[i]=Sl[i];
998  //dS[i]=dSl[i];
999  }
1000 
1001  }
1002 
1003  void updateSource2(const int c, TArray& S, TArray& dS)
1004  {
1005 
1006  const T hbarJoule=1.054571726e-34;
1007  const T hbar=6.582119e-16;
1008  const T Acell=5.378395621705545e-20;
1009  const T kb=8.617343e-5; // (eV/K)
1010  T cOt(0);
1011  S.zero();
1012  dS.zero();
1013  const int Rows=S.getLength();
1014  TArray Sl(Rows);
1015  TArray dSl(Rows);
1016  Sl.zero();
1017  dSl.zero();
1018  TArray e(Rows);
1019  //TArray e0(Rows);
1020  _kspace.geteCellVals(c,e);
1021  //_kspace.gete0CellVals(c,e0);
1022  const TArray& w(_kspace.getFreqArray());
1023  //const T Tlat(300.0);//Tl(_kspace.calcLatTemp(c));
1024  const T Tl(_kspace.calcLatTemp(c));
1025  TArray e0(Rows);
1026  //_kspace.gete0CellVals(c,e0);
1027  _kspace.getEquilibriumArray(e0,Tl);
1028  const T DK3=_kspace.getDK3();
1029 
1030  const IntArray& t1p2row=_type1Collisions.getSelfRow();
1031  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
1032  const IntArray& t1p2col=_type1Collisions.getSelfCol();
1033  const IntArray& t1p3col=_type1Collisions.getOtherCol();
1034  const IntArray& t2p2row=_type2Collisions.getSelfRow();
1035  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
1036  const IntArray& t2p2col=_type2Collisions.getSelfCol();
1037  const IntArray& t2p3col=_type2Collisions.getOtherCol();
1038 
1039  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
1040  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
1041  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
1042  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
1043 
1044  //type 1 collisions
1045 
1046  for(int i=0;i<Rows;i++)
1047  {
1048  //const int k1=floor(i/6);
1049  //const T dk31=_kspace.getkvol(k1).getdk3();
1050  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
1051  {
1052  const int j2=t1p2col[pos];
1053  //const int k2=floor(j2/6);
1054  //const T dk32=_kspace.getkvol(k2).getdk3();
1055  const int j3=t1p3col[pos];
1056  //const int k3=floor(j3/6);
1057  //const T dk33=_kspace.getkvol(k3).getdk3();
1058  const T dkl=t1dkl[pos];
1059  const T phi=t1phi[pos];
1060  //const T w3=w[i]+w[j2];
1061  //const T w1=w[j3]-w[j2];
1062  //const T n1=e[i]/w[i]/hbar;
1063  //const T n2=e[j2]/w[j2]/hbar;
1064  //const T n3=e[j3]/w[j3]/hbar;
1065  const T n10=1./(exp(hbar*w[i]/kb/300)-1.);
1066  const T n20=1./(exp(hbar*w[j2]/kb/300)-1.);
1067  const T n30=1./(exp(hbar*w[j3]/kb/300)-1.);
1068  //const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
1069  const T dn1=e[i];
1070  const T dn2=e[j2];
1071  const T dn3=e[j3];
1072  //T e032(0);
1073  //T frac(0);
1074  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
1075 
1076  const T nsum=(-dn1*dn2+dn3+dn1*dn3
1077  +dn2*dn3-dn2*n10+dn3*n10-dn1*n20
1078  +dn3*n20+dn1*n30+dn2*n30);
1079 
1080  //const T nsum1=(-dn1*n20+dn1*n30);
1081 
1082  const T ds1=-dn2+dn3-n20+n30;
1083  //const T ds2=-dn1+dn3-n10+n30;
1084  //const T ds3=1+dn1+dn2+n10+n20;
1085 
1086 
1087  Sl[i]+=dkl*phi*nsum;//*dk33/dk31;
1088  //Sl[j2]+=dkl*phi*nsum*dk31/dk32;
1089  //Sl[j3]-=dkl*phi*nsum*dk31/dk33;
1090  //Sl[j3]-=(dkl*phi*nsum*w[i]+dkl*phi*nsum*w[j2]);
1091 
1092  dSl[i]+=dkl*phi*ds1;
1093  //dSl[j2]+=dkl*phi*ds2*w[j2]*dk31/dk32;
1094  //dSl[j3]-=dkl*phi*ds3*w[j3]*dk31/dk33;
1095 
1096  }
1097  }
1098 
1099  T maxS(0);
1100 
1101  //type 2 collisions
1102 
1103  for(int i=0;i<Rows;i++)
1104  {
1105  //const int k1=floor(i/6);
1106  //const T dk31=_kspace.getkvol(k1).getdk3();
1107  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1108  {
1109  const int j2=t2p2col[pos];
1110  //const int k2=floor(j2/6);
1111  //const T dk32=_kspace.getkvol(k2).getdk3();
1112  const int j3=t2p3col[pos];
1113  //const int k3=floor(j3/6);
1114  //const T dk33=_kspace.getkvol(k3).getdk3();
1115  const T dkl=t2dkl[pos];
1116  const T phi=t2phi[pos];
1117  //const T w3=w[i]-w[j2];
1118  //const T w2=w[j3]-w[i];
1119  //const T w1=w[j3]-w[j2];
1120  //const T n1=e[i]/w[i]/hbar;
1121  //const T n2=e[j2]/w[j2]/hbar;
1122  //const T n3=e[j3]/w[j3]/hbar;
1123  const T n10=1./(exp(hbar*w[i]/kb/300)-1.);
1124  const T n20=1./(exp(hbar*w[j2]/kb/300)-1.);
1125  const T n30=1./(exp(hbar*w[j3]/kb/300)-1.);
1126  //const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
1127  //const T e012=hbar*w3/(exp(hbar*w3/kb/Tl)-1.);
1128  const T dn1=e[i];
1129  const T dn2=e[j2];
1130  const T dn3=e[j3];
1131  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
1132 
1133  const T nsum=(-dn1-dn1*dn2-dn1*dn3
1134  +dn2*dn3-dn2*n10-dn3*n10
1135  -dn1*n20+dn3*n20-dn1*n30
1136  +dn2*n30);
1137 
1138  const T ds1=-(1+dn2+dn3+n20+n30);
1139  //const T ds2=-dn1+dn3-n10+n30;
1140  //const T ds3=-dn1+dn2-n10+n20;
1141 
1142 
1143  //const T nsum1=(-dn1-dn1*n20-dn1*n30);
1144  //const T nsum2=(-dn2*n10+dn2*n30);
1145  //const T nsum3=(-dn3*n10+dn3*n20);
1146 
1147 
1148  Sl[i]+=0.5*dkl*phi*nsum;
1149  //Sl[i]+=0.5*dkl*phi*nsum*w[j2]+0.5*dkl*phi*nsum*w[j3];
1150  dSl[i]+=0.5*dkl*phi*ds1;
1151 
1152  //Sl[j2]-=0.5*dkl*phi*nsum*dk31/dk32;
1153  //dSl[j2]-=0.5*dkl*phi*ds2*w[j2]*dk31/dk32;
1154 
1155  //Sl[j3]-=0.5*dkl*phi*nsum*dk31/dk33;
1156  //dSl[j3]-=0.5*dkl*phi*ds3*w[j3]*dk31/dk33;
1157 
1158  }
1159  }
1160 
1161  const T preFac=Acell/16./3.141592653*hbarJoule*_maxPhi*_maxDkl;
1162  const int klen=_kspace.getlength();
1163  for(int i=0;i<Rows;i++)
1164  {
1165  Sl[i]*=preFac;//*w[i]*hbar);
1166  dSl[i]*=preFac;
1167  if(fabs(Sl[i])>fabs(maxS))
1168  maxS=fabs(Sl[i]);
1169  }
1170 
1171  long double defect(0);
1172  int cnt(0);
1173  T eq(0);
1174  const int cStart(_kspace.getGlobalIndex(c,0));
1175  int ind(cStart);
1176 
1177  for(int k=0;k<klen;k++)
1178  {
1179  Tkvol& kv=_kspace.getkvol(k);
1180  const int modenum=kv.getmodenum();
1181  T dk3=kv.getdk3();
1182  for(int m=0;m<modenum;m++)
1183  {
1184  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
1185  //if(fabs(Sl[cnt])*(dk3/DK3)>100)
1186  // Sl[cnt]=0.;
1187  Tmode& mode=kv.getmode(m);
1188  const T tau=mode.gettau();
1189  cOt+=mode.calce0(Tl)*(dk3/DK3);
1190  //const T e0=mode.calce0(Tl);
1191  eq+=(e0[ind]-e[ind])/tau*(dk3/DK3);
1192  ind++;
1193  cnt++;
1194  }
1195  }
1196 
1197  cOt*=DK3;
1198  //T f(maxS*defect/eq);
1199  //cout<<"Defect1: "<<defect<<endl;
1200  //cout<<"Temperature: "<<Tl<<endl;
1201 
1202  if(fabs(defect)>0.)
1203  {
1204  defect*=maxS*DK3;
1205  cnt=0;
1206  //T d(defect);
1207  defect=0.;
1208  ind=cStart;
1209  for(int k=0;k<klen;k++)
1210  {
1211  Tkvol& kv=_kspace.getkvol(k);
1212  const int modenum=kv.getmodenum();
1213  T dk3=kv.getdk3();
1214  for(int m=0;m<modenum;m++)
1215  {
1216  //Tmode& mode=kv.getmode(m);
1217  //const T tau=mode.gettau();
1218  //const T cp=mode.calcde0dT(Tl);
1219  //const T e0(mode.calce0(Tl));
1220  //S[cnt]-=f*(e0[ind]-e[ind])/tau;
1221  //Sl[cnt]=d*(Sl[cnt]/d-e0/cOt);
1222  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
1223  ind++;
1224  cnt++;
1225  }
1226  }
1227  //defect*=DK3*maxS;
1228  //cout<<"Defect2: "<<defect<<endl;
1229  //cout<<"f: "<<f<<endl;
1230  }
1231 
1232  for(int i=0;i<Rows;i++)
1233  {
1234  S[i]=Sl[i];
1235  //dS[i]=dSl[i];
1236  }
1237 
1238  }
1239 
1240  void getTypeIIsource(const int c, TArray& S, TArray& dS, const bool correct)
1241  {
1242 
1243  const T hbarJoule=1.054571726e-34;
1244  const T hbar=6.582119e-16;
1245  const T Acell=5.378395621705545e-20;
1246  const T kb=8.617343e-5; // (eV/K)
1247  T cOt(0);
1248  S.zero();
1249  dS.zero();
1250  const int Rows=S.getLength();
1251  TArray Sl(Rows);
1252  TArray dSl(Rows);
1253  Sl.zero();
1254  dSl.zero();
1255  TArray e(Rows);
1256  //TArray e0(Rows);
1257  _kspace.geteCellVals(c,e);
1258  //_kspace.gete0CellVals(c,e0);
1259  const TArray& w(_kspace.getFreqArray());
1260  //const T Tlat(300.0);//Tl(_kspace.calcLatTemp(c));
1261  const T Tl(_kspace.calcLatTemp(c));
1262  TArray e0(Rows);
1263  //_kspace.gete0CellVals(c,e0);
1264  _kspace.getEquilibriumArray(e0,Tl);
1265  const T DK3=_kspace.getDK3();
1266 
1267  const IntArray& t2p2row=_type2Collisions.getSelfRow();
1268  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
1269  const IntArray& t2p2col=_type2Collisions.getSelfCol();
1270  const IntArray& t2p3col=_type2Collisions.getOtherCol();
1271 
1272  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
1273  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
1274 
1275  T maxS(0);
1276 
1277  //type 2 collisions
1278 
1279  for(int i=0;i<Rows;i++)
1280  {
1281  const int k1=floor(i/6);
1282  const T dk31=_kspace.getkvol(k1).getdk3();
1283  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1284  {
1285  const int j2=t2p2col[pos];
1286  const int k2=floor(j2/6);
1287  const T dk32=_kspace.getkvol(k2).getdk3();
1288  const int j3=t2p3col[pos];
1289  const int k3=floor(j3/6);
1290  const T dk33=_kspace.getkvol(k3).getdk3();
1291  const T dkl=t2dkl[pos];
1292  const T phi=t2phi[pos];
1293  const T w3=w[i]-w[j2];
1294  //const T w2=w[j3]-w[i];
1295  //const T w1=w[j3]-w[j2];
1296  const T n1=e[i]/w[i]/hbar;
1297  const T n2=e[j2]/w[j2]/hbar;
1298  const T n3=e[j3]/w[j3]/hbar;
1299  const T n10=e0[i]/w[i]/hbar;
1300  const T n20=e0[j2]/w[j2]/hbar;
1301  const T n30=e0[j3]/w[j3]/hbar;
1302  //const T n012=1./(exp(hbar*w3/kb/Tl)-1.);
1303  //const T e012=hbar*w3/(exp(hbar*w3/kb/Tl)-1.);
1304  const T dn1=n1-n10;
1305  const T dn2=n2-n20;
1306  const T dn3=(n3-n30);//*n012/n30;
1307  //const T frac=w[j3]*(1-w[j2]/w[j3])/w[i];
1308 
1309  T Large(n10);
1310 
1311  if(n20>Large)
1312  Large=n20;
1313  if(n30>Large)
1314  Large=n30;
1315 
1316 
1317  const T nsum=(-dn1-dn1*dn2-dn1*dn3
1318  +dn2*dn3-dn2*n10-dn3*n10
1319  -dn1*n20+dn3*n20-dn1*n30
1320  +dn2*n30);
1321 
1322  const T ds1=-(1+dn2+dn3-n20-n30);
1323  const T ds2=-dn1+dn3-n10+n30;
1324  const T ds3=-dn1+dn2-n10+n20;
1325 
1326 
1327  //const T nsum1=(-dn1-dn1*n20-dn1*n30);
1328  //const T nsum1=(n30*n20*(n1+1)-n1*(n20+1)*(n30+1));
1329  //const T nsum2=Large*(-dn2*n10/Large+dn2*n30/Large);
1330  //const T nsum3=Large*(-dn3*n10/Large+dn3*n20/Large);
1331 
1332 
1333  Sl[i]+=0.5*dkl*phi*nsum*w[i];
1334  dSl[i]+=0.5*dkl*phi*ds1*w[i];
1335 
1336  //Sl[i]+=0.25*dkl*phi*nsum*w[j2]*dk31/dk32+0.25*dkl*phi*nsum*w3*dk31/dk33;
1337  //Sl[i]+=0.5*dkl*phi*nsum*w[j2]+0.5*dkl*phi*nsum*w3;
1338 
1339  Sl[j2]-=0.5*dkl*phi*nsum*w[j2]*dk31/dk32;
1340  dSl[j2]-=0.5*dkl*phi*ds2*w[j2]*dk31/dk32;
1341 
1342  Sl[j3]-=0.5*dkl*phi*nsum*w3*dk31/dk33;
1343  dSl[j3]-=0.5*dkl*phi*ds3*w3*dk31/dk33;
1344 
1345  //Sl[j2]-=0.25*dkl*phi*nsum*w[j2]*dk31/dk32;
1346  //Sl[j3]-=0.25*dkl*phi*nsum*w3*dk31/dk33;//*e012/e0[j3];
1347  //S[i]+=0.5*dkl*phi*(n3*n2*(n1+1)-n1*(n2+1)*(n3+1));
1348 
1349  //if(j2==0)
1350  //int yes(0);
1351  //if(j3==0)
1352  //int yes(0);
1353 
1354  }
1355  }
1356 
1357  const T preFac=Acell/16./3.141592653*hbarJoule*_maxPhi*hbar*_maxDkl;
1358  const int klen=_kspace.getlength();
1359  for(int i=0;i<Rows;i++)
1360  {
1361  Sl[i]*=preFac;//*w[i]*hbar);
1362  dSl[i]*=preFac;
1363  if(fabs(Sl[i])>fabs(maxS))
1364  maxS=fabs(Sl[i]);
1365  }
1366 
1367  long double defect(0);
1368  int cnt(0);
1369  T eq(0);
1370  const int cStart(_kspace.getGlobalIndex(c,0));
1371  int ind(cStart);
1372 
1373  for(int k=0;k<klen;k++)
1374  {
1375  Tkvol& kv=_kspace.getkvol(k);
1376  const int modenum=kv.getmodenum();
1377  T dk3=kv.getdk3();
1378  for(int m=0;m<modenum;m++)
1379  {
1380  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
1381  //if(fabs(Sl[cnt])*(dk3/DK3)>100)
1382  // Sl[cnt]=0.;
1383  Tmode& mode=kv.getmode(m);
1384  //const T tau=mode.gettau();
1385  cOt+=mode.calce0(Tl)*(dk3/DK3);
1386  //const T e0=mode.calce0(Tl);
1387  //eq+=(e0-e[ind])/tau*(dk3/DK3);
1388  ind++;
1389  cnt++;
1390  }
1391  }
1392 
1393  cOt*=DK3;
1394  //T f(maxS*defect/eq);
1395  defect*=maxS*DK3;
1396  //cout<<"Defect1: "<<defect<<endl;
1397  //cout<<"Temperature: "<<Tl<<endl;
1398 
1399  if(correct)
1400  {
1401  cnt=0;
1402  //T d(defect);
1403  defect=0.;
1404  ind=cStart;
1405  for(int k=0;k<klen;k++)
1406  {
1407  Tkvol& kv=_kspace.getkvol(k);
1408  const int modenum=kv.getmodenum();
1409  T dk3=kv.getdk3();
1410  for(int m=0;m<modenum;m++)
1411  {
1412  Tmode& mode=kv.getmode(m);
1413  //const T tau=mode.gettau();
1414  //const T cp=mode.calcde0dT(Tl);
1415  //const T e0(mode.calce0(Tl));
1416  //S[cnt]+=f*(e0-e[ind])/tau;
1417  //Sl[cnt]=d*(Sl[cnt]/d-e0/cOt);
1418  defect+=((Sl[cnt]/maxS)*(dk3/DK3));
1419  ind++;
1420  cnt++;
1421  }
1422  }
1423  defect*=DK3*maxS;
1424  //cout<<"Defect2: "<<defect<<endl;
1425  }
1426 
1427  for(int i=0;i<Rows;i++)
1428  {
1429  S[i]=Sl[i];
1430  dS[i]=dSl[i];
1431  }
1432 
1433  }
1434 
1435  ArrayBase* IterateToEquilibrium(const T Tl, const int totIts, const T tStep)
1436  {
1437 
1438  const T hbarJoule=1.054560652926899e-34;
1439  const T hbar=6.582119e-16;
1440  //const T JouleToeV=6.24150974e18;
1441  const T Acell=5.378395621705545e-20;
1442  const T kb=8.617343e-5; // (eV/K)
1443  //const T cv=1e-9*1e-9;
1444  T temp=Tl;
1445  //T Tlat(_kspace.calcLatTemp(0));
1446  T maxS(0);
1447 
1448  TArray S(_kspace.gettotmodes());
1449  S.zero();
1450  TArray e0(_kspace.gettotmodes());
1451  e0.zero();
1452  TArray* ep=new TArray(_kspace.gettotmodes());
1453  TArray& e=*ep;
1454  e.zero();
1455  //_kspace.getEquilibriumArray(e,Tl);
1456  //_kspace.geteCellVals(0,e);
1457  TArray& w(_kspace.getFreqArray());
1458  TArray t(_kspace.gettotmodes());
1459  t.zero();
1460  TArray gam1(_kspace.gettotmodes());
1461  gam1.zero();
1462  TArray gam2(_kspace.gettotmodes());
1463  gam2.zero();
1464  TArray gam(_kspace.gettotmodes());
1465  gam.zero();
1466 
1467  const int Rows=S.getLength();
1468  const IntArray& t1p2row=_type1Collisions.getSelfRow();
1469  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
1470  const IntArray& t1p2col=_type1Collisions.getSelfCol();
1471  const IntArray& t1p3col=_type1Collisions.getOtherCol();
1472  const IntArray& t2p2row=_type2Collisions.getSelfRow();
1473  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
1474  const IntArray& t2p2col=_type2Collisions.getSelfCol();
1475  const IntArray& t2p3col=_type2Collisions.getOtherCol();
1476 
1477  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
1478  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
1479  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
1480  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
1481  const T preFac=Acell/16.*hbarJoule*_maxPhi*_maxDkl/3.14159265358979323846264;
1482 
1483  for(int i=0;i<Rows;i++)
1484  e[i]=hbar*w[i]/(exp(hbar*w[i]/kb/Tl)-1);
1485 
1486  e[100]=hbar*w[100]/(exp(hbar*w[100]/kb/(Tl+1.))-1);
1487 
1488  for(int it=0;it<totIts;it++)
1489  {
1490  S.zero();
1491  T cOt(0);
1492  //type 1 collisions
1493 
1494  for(int i=0;i<Rows;i++)
1495  {
1496  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
1497  {
1498  const int j2=t1p2col[pos];
1499  const int j3=t1p3col[pos];
1500  const T dkl=t1dkl[pos];
1501  const T phi=t1phi[pos];
1502  //const T w3=w[i]+w[j2];
1503  const T n1=e[i]/hbar/w[i];
1504  const T n2=e[j2]/hbar/w[j2];
1505  const T n3=e[j3]/hbar/w[j3];
1506  const T n10=1./(exp(hbar*w[i]/kb/300.)-1);
1507  const T n20=1./(exp(hbar*w[j2]/kb/300.)-1);
1508  const T n30=1./(exp(hbar*w[j3]/kb/300.)-1);
1509  const T dn1=n1-n10;
1510  const T dn2=n2-n20;
1511  const T dn3=n3-n30;
1512 
1513  const T nsum=(-dn1*dn2+dn3+dn1*dn3
1514  +dn2*dn3-dn2*n10+dn3*n10-dn1*n20
1515  +dn3*n20+dn1*n30+dn2*n30);
1516 
1517  S[i]+=dkl*phi*nsum;
1518  t[i]+=dkl*phi*(n20-n30);
1519  gam1[i]+=dkl*phi*nsum;
1520  }
1521 
1522  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1523  {
1524  const int j2=t2p2col[pos];
1525  const int k2=floor(j2/6);
1526  //const T dk32=_kspace.getkvol(k2).getdk3();
1527  const int j3=t2p3col[pos];
1528  const int k3=floor(j3/6);
1529  //const T dk33=_kspace.getkvol(k3).getdk3();
1530  const T dkl=t2dkl[pos];
1531  const T phi=t2phi[pos];
1532  //const T w3=w[i]-w[j2];
1533  const T n1=e[i]/w[i]/hbar;
1534  const T n2=e[j2]/w[j2]/hbar;
1535  const T n3=e[j3]/w[j3]/hbar;
1536  const T n10=1./(exp(hbar*w[i]/kb/300.)-1);
1537  const T n20=1./(exp(hbar*w[j2]/kb/300.)-1);
1538  const T n30=1./(exp(hbar*w[j3]/kb/300.)-1);
1539  const T dn1=n1-n10;
1540  const T dn2=n2-n20;
1541  const T dn3=n3-n30;
1542 
1543  const T nsum=(-dn1-dn1*dn2-dn1*dn3
1544  +dn2*dn3-dn2*n10-dn3*n10
1545  -dn1*n20+dn3*n20-dn1*n30
1546  +dn2*n30);
1547 
1548  S[i]+=0.5*dkl*phi*nsum;
1549  //S[j2]-=0.5*dkl*phi*nsum*dk31/dk32;
1550  //S[j3]-=0.5*dkl*phi*nsum*dk31/dk33;
1551  t[i]+=0.5*dkl*phi*(n20+n30+1);
1552  gam2[i]+=0.5*dkl*phi*nsum;
1553  }
1554  //S[i]*=(preFac*w[i]*hbar);
1555  //e[i]+=S[i]*tStep;
1556 
1557  }
1558 
1559  //type 2 collisions
1560  /*
1561  for(int i=0;i<Rows;i++)
1562  {
1563  const int k1=floor(i/6);
1564  const T dk31=_kspace.getkvol(k1).getdk3();
1565  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1566  {
1567  const int j2=t2p2col[pos];
1568  const int k2=floor(j2/6);
1569  const T dk32=_kspace.getkvol(k2).getdk3();
1570  const int j3=t2p3col[pos];
1571  const int k3=floor(j3/6);
1572  const T dk33=_kspace.getkvol(k3).getdk3();
1573  const T dkl=t2dkl[pos];
1574  const T phi=t2phi[pos];
1575  const T w3=w[i]-w[j2];
1576  const T n1=e[i]/w[i]/hbar;
1577  const T n2=e[j2]/w[j2]/hbar;
1578  const T n3=e[j3]/w[j3]/hbar;
1579  const T n10=1./(exp(hbar*w[i]/kb/300.)-1);
1580  const T n20=1./(exp(hbar*w[j2]/kb/300.)-1);
1581  const T n30=1./(exp(hbar*w[j3]/kb/300.)-1);
1582  const T dn1=n1-n10;
1583  const T dn2=n2-n20;
1584  const T dn3=n3-n30;
1585 
1586  const T nsum=(-dn1-dn1*dn2-dn1*dn3
1587  +dn2*dn3-dn2*n10-dn3*n10
1588  -dn1*n20+dn3*n20-dn1*n30
1589  +dn2*n30);
1590 
1591  S[i]+=0.5*dkl*phi*nsum;
1592  //S[j2]-=0.5*dkl*phi*nsum*dk31/dk32;
1593  //S[j3]-=0.5*dkl*phi*nsum*dk31/dk33;
1594  t[i]+=0.5*dkl*phi*(n20+n30+1);
1595  gam2[i]+=0.5*dkl*phi*n10*(n20+1)*(n30+1);
1596  }
1597  //S[i]*=(preFac*w[i]*hbar);
1598  //e[i]+=S[i]*tStep;
1599  }
1600  */
1601 
1602  const int klen=_kspace.getlength();
1603  for(int i=0;i<Rows;i++)
1604  {
1605  S[i]*=(preFac*w[i]*hbar);
1606  t[i]*=preFac;
1607  t[i]=1./t[i];
1608  const T n10=1./(exp(hbar*w[i]/kb/300.)-1);
1609  gam[i]=n10*(n10+1)/preFac/(gam1[i]+gam2[i]);
1610  if(fabs(S[i])>fabs(maxS))
1611  maxS=fabs(S[i]);
1612 
1613  }
1614 
1615  const T DK3=_kspace.getDK3();
1616  T defect(0);
1617  int cnt(0);
1618  T maxRatio(0);
1619  int maxInd(-1);
1620  int minInd(-1);
1621  T minRatio(10);
1622  for(int k=0;k<klen;k++)
1623  {
1624  Tkvol& kv=_kspace.getkvol(k);
1625  const int modenum=kv.getmodenum();
1626  T dk3=kv.getdk3();
1627  for(int m=0;m<modenum;m++)
1628  {
1629  Tmode& mode=kv.getmode(m);
1630  const T tau=mode.gettau();
1631  //const T dif=tau-t[cnt];
1632  const T gamRatio=tau/gam[cnt];
1633  if(gamRatio>maxRatio)
1634  {
1635  maxRatio=gamRatio;
1636  maxInd=cnt;
1637  }
1638  if(gamRatio<minRatio)
1639  {
1640  minRatio=gamRatio;
1641  minInd=cnt;
1642  }
1643  defect+=(S[cnt]/maxS)*(dk3/DK3);
1644  cnt++;
1645  }
1646  }
1647 
1648  //cout<<"Max difference: "<<maxRatio<<" at: "<<maxInd<<endl;
1649  //cout<<"Min difference: "<<minRatio<<" at: "<<minInd<<endl;
1650 
1651 
1652  //T eq(0);
1653  if(it>-1)
1654  {
1655  cnt=0;
1656  for(int k=0;k<klen;k++)
1657  {
1658  Tkvol& kv=_kspace.getkvol(k);
1659  const int modenum=kv.getmodenum();
1660  T dk3=kv.getdk3();
1661  for(int m=0;m<modenum;m++)
1662  {
1663  Tmode& mode=kv.getmode(m);
1664  const T ez=mode.calce0(Tl);
1665  cOt+=ez*dk3/DK3;
1666  //eq+=(mode.calce0(temp)-e[cnt])*(dk3/tau/DK3);
1667  cnt++;
1668  }
1669  }
1670 
1671  defect=defect*maxS*DK3;
1672  cOt*=DK3;
1673  //cout<<"Defect1: "<<defect<<endl;
1674  cnt=0;
1675  //T d(defect);
1676  defect=0.;
1677  for(int k=0;k<klen;k++)
1678  {
1679  Tkvol& kv=_kspace.getkvol(k);
1680  const int modenum=kv.getmodenum();
1681  T dk3=kv.getdk3();
1682  for(int m=0;m<modenum;m++)
1683  {
1684  Tmode& mode=kv.getmode(m);
1685  //const T ez=mode.calce0(Tl);
1686  //S[cnt]=d*(S[cnt]/d-ez/cOt);
1687  defect+=(S[cnt]/maxS)*(dk3/DK3);
1688  cnt++;
1689  }
1690  }
1691  }
1692 
1693 
1694  for(int i=0;i<Rows;i++)
1695  e[i]+=S[i]*gam[i];
1696 
1697  cnt=0;
1698  T esum(0);
1699 
1700  for(int k=0;k<klen;k++)
1701  {
1702  Tkvol& kv=_kspace.getkvol(k);
1703  const int modenum=kv.getmodenum();
1704  T dk3=kv.getdk3();
1705  for(int m=0;m<modenum;m++)
1706  {
1707  esum+=e[cnt]*dk3;
1708  cnt++;
1709  }
1710  }
1711 
1712  _kspace.calcTemp(temp,esum);
1713  cout<<"New Temp: "<<temp<<endl;
1714  }
1715 
1716  _kspace.weightArray(e);
1717 
1718  return ep;
1719  }
1720 
1721  void correctDetailedBalance()
1722  {
1723  const IntArray& t1p2row=_type1Collisions.getSelfRow();
1724  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
1725  const IntArray& t1p2col=_type1Collisions.getSelfCol();
1726  const IntArray& t1p3col=_type1Collisions.getOtherCol();
1727 
1728  //const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
1729  TArray& t1phi=_type1Collisions.getNonConstOtherCoeffs();
1730  BArray Dups(t1phi.getLength());
1731  Dups=false;
1732 
1733  const int Rows(_kspace.gettotmodes());
1734  int duplicated(0);
1735  int removed(0);
1736 
1737  for(int i=0;i<Rows;i++)
1738  {
1739  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
1740  {
1741  const int j2=t1p2col[pos];
1742  //const int j3=t1p3col[pos];
1743  int p1(-1);
1744  //int p2(-1);
1745  //int p3(-1);
1746 
1747  //type1 for j2
1748  bool hasI1(false);
1749  if(true)
1750  {
1751  for(int pos2=t1p2row[j2];pos2<t1p2row[j2+1];pos2++)
1752  {
1753  const int jj2=t1p2col[pos2];
1754  if(i==jj2)
1755  {
1756  hasI1=true;
1757  p1=pos2;
1758  break;
1759  }
1760  }
1761  }
1762 
1763  if((hasI1))
1764  {
1765  Dups[pos]=true;
1766  Dups[p1]=true;
1767  t1phi[pos]*=0.5;
1768  duplicated++;
1769  }
1770  }
1771  }
1772 
1773  for(int i=0;i<Dups.getLength();i++)
1774  {
1775  if(!Dups[i])
1776  {
1777  //t1phi[i]=0;
1778  removed++;
1779  }
1780  }
1781 
1782  cout<<"Duplicated -- "<<duplicated<<endl;
1783  cout<<"Removed -- "<<removed<<endl;
1784 
1785  }
1786 
1787  ArrayBase* calculatePsi(const int totIts)
1788  {
1789 
1790  const int Rows=_kspace.gettotmodes();
1791  VectorT3Array* PsiP=new VectorT3Array(Rows);
1792  VectorT3Array& Psi=*PsiP;
1793  Psi.zero();
1794  VectorT3Array v(Rows);
1795  VectorT3Array sum(Rows);
1796  sum.zero();
1797  VectorT3Array Psi0(Rows);
1798  Psi0.zero();
1799  TArray Gam(Rows);
1800  Gam.zero();
1801  TArray n0(Rows);
1802  n0.zero();
1803  const T hbarJoule=1.054571726e-34;
1804  const T hbar=6.582119e-16;
1805  const T Acell=5.378395621705545e-20;
1806  const T kb=8.617343e-5; // (eV/K)
1807  const TArray& w(_kspace.getFreqArray());
1808  _kspace.getVelocities(v);
1809  const T Tl(300.0);
1810  //const T DK3=_kspace.getDK3();
1811 
1812  const IntArray& t1p2row=_type1Collisions.getSelfRow();
1813  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
1814  const IntArray& t1p2col=_type1Collisions.getSelfCol();
1815  const IntArray& t1p3col=_type1Collisions.getOtherCol();
1816  const IntArray& t2p2row=_type2Collisions.getSelfRow();
1817  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
1818  const IntArray& t2p2col=_type2Collisions.getSelfCol();
1819  const IntArray& t2p3col=_type2Collisions.getOtherCol();
1820 
1821  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
1822  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
1823  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
1824  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
1825 
1826  for(int i=0;i<Rows;i++)
1827  n0[i]=1./(exp(hbar*w[i]/kb/Tl)-1);
1828 
1829  //first calculate Gamma for each mode
1830 
1831  //type 1 collisions
1832 
1833  for(int i=0;i<Rows;i++)
1834  {
1835  const int k1=floor(i/6);
1836  //const T dk31=_kspace.getkvol(k1).getdk3();
1837  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
1838  {
1839  const int j2=t1p2col[pos];
1840  const int k2=floor(j2/6);
1841  //const T dk32=_kspace.getkvol(k2).getdk3();
1842  const int j3=t1p3col[pos];
1843  const int k3=floor(j3/6);
1844  //const T dk33=_kspace.getkvol(k3).getdk3();
1845  const T dkl=t1dkl[pos];
1846  const T phi=t1phi[pos];
1847 
1848  Gam[i]+=n0[i]*n0[j2]*(n0[j3]+1)*dkl*phi;
1849 
1850  }
1851  }
1852 
1853  //T maxS(0);
1854 
1855  //type 2 collisions
1856 
1857  for(int i=0;i<Rows;i++)
1858  {
1859  const int k1=floor(i/6);
1860  //const T dk31=_kspace.getkvol(k1).getdk3();
1861  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1862  {
1863  const int j2=t2p2col[pos];
1864  const int k2=floor(j2/6);
1865  //const T dk32=_kspace.getkvol(k2).getdk3();
1866  const int j3=t2p3col[pos];
1867  const int k3=floor(j3/6);
1868  //const T dk33=_kspace.getkvol(k3).getdk3();
1869  const T dkl=t2dkl[pos];
1870  const T phi=t2phi[pos];
1871 
1872  Gam[i]+=0.5*n0[i]*(n0[j2]+1)*(n0[j3]+1)*phi*dkl;
1873 
1874  }
1875  }
1876 
1877  const T preFac=Acell/16./3.141592653*hbarJoule*_maxPhi*_maxDkl;
1878  for(int i=0;i<Rows;i++)
1879  Gam[i]*=preFac;
1880 
1881  //Initialize psi at SMRT
1882 
1883  for(int i=0;i<Rows;i++)
1884  {
1885  Psi0[i]=hbarJoule*w[i]*n0[i]*(n0[i]+1)/Gam[i]/Tl*v[i];
1886  Psi[i]=hbarJoule*w[i]*n0[i]*(n0[i]+1)/Gam[i]/Tl*v[i];
1887  if(i==10)
1888  {
1889  Psi0[i]=hbarJoule*w[i]*n0[i]*(n0[i]+1)/Gam[i]/Tl*v[i]*1.1;
1890  Psi[i]=hbarJoule*w[i]*n0[i]*(n0[i]+1)/Gam[i]/Tl*v[i]*1.1;
1891  }
1892  }
1893 
1894  //Iterate to get Full Scattering answer
1895 
1896 
1897  for(int its=0;its<totIts;its++)
1898  {
1899  sum.zero();
1900  for(int i=0;i<Rows;i++)
1901  {
1902  const int k1=floor(i/6);
1903  //const T dk31=_kspace.getkvol(k1).getdk3();
1904  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
1905  {
1906  const int j2=t1p2col[pos];
1907  const int k2=floor(j2/6);
1908  //const T dk32=_kspace.getkvol(k2).getdk3();
1909  const int j3=t1p3col[pos];
1910  const int k3=floor(j3/6);
1911  //const T dk33=_kspace.getkvol(k3).getdk3();
1912  const T dkl=t1dkl[pos];
1913  const T phi=t1phi[pos];
1914 
1915  sum[i]+=n0[i]*n0[j2]*(n0[j3]+1)*(Psi[j3]-Psi[j2])*phi*dkl*preFac/Gam[i];
1916  }
1917  }
1918 
1919  //T maxS(0);
1920 
1921  //type 2 collisions
1922 
1923  for(int i=0;i<Rows;i++)
1924  {
1925  const int k1=floor(i/6);
1926  //const T dk31=_kspace.getkvol(k1).getdk3();
1927  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
1928  {
1929  const int j2=t2p2col[pos];
1930  const int k2=floor(j2/6);
1931  //const T dk32=_kspace.getkvol(k2).getdk3();
1932  const int j3=t2p3col[pos];
1933  const int k3=floor(j3/6);
1934  //const T dk33=_kspace.getkvol(k3).getdk3();
1935  const T dkl=t2dkl[pos];
1936  const T phi=t2phi[pos];
1937 
1938  sum[i]+=0.5*n0[i]*(n0[j2]+1)*(n0[j3]+1)*(Psi[j3]+Psi[j2])*phi*dkl*preFac/Gam[i];
1939 
1940  }
1941  Psi[i]=sum[i]+Psi0[i];
1942  }
1943  }
1944 
1945  return PsiP;
1946 
1947  }
1948 
1949  void scatterPhonons(const int c, const int totIts,
1950  TArray& C, TArray& B,
1951  const TArray& V, TArray& newE, const T cv)
1952  {
1953 
1954  const int Rows=_kspace.gettotmodes();
1955  TArray P(Rows);
1956  P.zero();
1957  TArray P0(Rows);
1958  P0.zero();
1959  TArray n0(Rows);
1960  n0.zero();
1961  TArray n(Rows);
1962  n.zero();
1963  TArray Q(Rows);
1964  Q.zero();
1965  T Tl(_kspace.calcLatTemp(c));
1966  //_kspace.getEquilibriumArray(n,Tl);
1967  _kspace.gete0CellVals(c,n0);
1968  _kspace.geteCellVals(c,n);
1969  const T hbarJoule=1.054571726e-34;
1970  //const T hbar=6.582119e-16;
1971  const T Acell=5.378395621705545e-20;
1972  //const T kb=8.617343e-5; // (eV/K)
1973  const T kbJoule=1.3806488e-23; // (J/K)
1974  //const TArray& w(_kspace.getFreqArray());
1975  //const T DK3=_kspace.getDK3();
1976 
1977  const IntArray& t1p2row=_type1Collisions.getSelfRow();
1978  //const IntArray& t1p3row=_type1Collisions.getOtherRow();
1979  const IntArray& t1p2col=_type1Collisions.getSelfCol();
1980  const IntArray& t1p3col=_type1Collisions.getOtherCol();
1981  const IntArray& t2p2row=_type2Collisions.getSelfRow();
1982  //const IntArray& t2p3row=_type2Collisions.getOtherRow();
1983  const IntArray& t2p2col=_type2Collisions.getSelfCol();
1984  const IntArray& t2p3col=_type2Collisions.getOtherCol();
1985 
1986  const TArray& t1dkl=_type1Collisions.getSelfCoeffs();
1987  const TArray& t1phi=_type1Collisions.getOtherCoeffs();
1988  const TArray& t2dkl=_type2Collisions.getSelfCoeffs();
1989  const TArray& t2phi=_type2Collisions.getOtherCoeffs();
1990 
1991  /*
1992  for(int i=0;i<Rows;i++)
1993  {
1994  n0[i]/=(hbar*w[i]); //Since units are in eV when it comes in
1995  C[i]/=(hbar*w[i]);
1996  B[i]/=(hbar*w[i]);
1997  n[i]/=(hbar*w[i]);
1998  }*/
1999 
2000  //first calculate Q for each mode
2001 
2002  //type 1 collisions
2003 
2004  cout<<"Temperature: "<<Tl<<endl;
2005 
2006  for(int i=0;i<Rows;i++)
2007  {
2008  //const int k1=floor(i/6);
2009  //const T dk31=_kspace.getkvol(k1).getdk3();
2010 
2011  //type I collisions
2012  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
2013  {
2014  const int j2=t1p2col[pos];
2015  //const int k2=floor(j2/6);
2016  //const T dk32=_kspace.getkvol(k2).getdk3();
2017  const int j3=t1p3col[pos];
2018  //const int k3=floor(j3/6);
2019  //const T dk33=_kspace.getkvol(k3).getdk3();
2020  const T dkl=t1dkl[pos];
2021  const T phi=t1phi[pos];
2022 
2023  Q[i]+=n0[i]*n0[j2]*(n0[j3]+1)*dkl*phi;
2024  }
2025 
2026  //type II collisions
2027  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
2028  {
2029  const int j2=t2p2col[pos];
2030  //const int k2=floor(j2/6);
2031  //const T dk32=_kspace.getkvol(k2).getdk3();
2032  const int j3=t2p3col[pos];
2033  //const int k3=floor(j3/6);
2034  //const T dk33=_kspace.getkvol(k3).getdk3();
2035  const T dkl=t2dkl[pos];
2036  const T phi=t2phi[pos];
2037 
2038  Q[i]+=0.5*n0[i]*(n0[j2]+1)*(n0[j3]+1)*phi*dkl;
2039  }
2040  }
2041 
2042 
2043  const T preFac=Acell/16./3.141592653*hbarJoule*_maxPhi*_maxDkl;
2044  for(int i=0;i<Rows;i++)
2045  Q[i]*=preFac;
2046 
2047 
2048  //Initialize P with constant terms
2049 
2050  for(int i=0;i<Rows;i++)
2051  {
2052  if(Q[i]>0)
2053  P0[i]=C[i]/Q[i]*kbJoule*Tl;
2054  else
2055  P0[i]=0.;
2056 
2057  P[i]=P0[i];
2058  }
2059 
2060  //Iterate to get Full Scattering answer
2061 
2062  for(int its=0;its<totIts;its++)
2063  {
2064 
2065  for(int i=0;i<Rows;i++)
2066  {
2067  T sum(0);
2068  for(int pos=t1p2row[i];pos<t1p2row[i+1];pos++)
2069  {
2070  const int j2=t1p2col[pos];
2071  const int j3=t1p3col[pos];
2072  const T dkl=t1dkl[pos];
2073  const T phi=t1phi[pos];
2074 
2075  if(Q[i]>0)
2076  sum+=n0[i]*n0[j2]*(n0[j3]+1)*(P[j3]-P[j2])*phi*dkl*preFac/Q[i];
2077  }
2078 
2079  for(int pos=t2p2row[i];pos<t2p2row[i+1];pos++)
2080  {
2081  const int j2=t2p2col[pos];
2082  const int j3=t2p3col[pos];
2083  const T dkl=t2dkl[pos];
2084  const T phi=t2phi[pos];
2085 
2086  if(Q[i]>0)
2087  sum+=0.5*n0[i]*(n0[j2]+1)*(n0[j3]+1)*(P[j3]+P[j2])*phi*dkl*preFac/Q[i];
2088 
2089  }
2090 
2091  P[i]=0.75*P[i]+0.25*(sum+P0[i]);
2092 
2093  }
2094 
2095  T maxChange(0);
2096  T minChange(10);
2097  for(int i=0;i<Rows;i++)
2098  {
2099  T factor(1+P[i]*(1+n0[i])/kbJoule/Tl);
2100 
2101  if(factor>maxChange)
2102  maxChange=factor;
2103  if(factor<minChange)
2104  minChange=factor;
2105  }
2106 
2107  cout<<"Max Change: "<<maxChange<<" Min Change: "<<minChange<<endl;
2108 
2109  }
2110 
2111  T maxChange(0);
2112  T minChange(10);
2113  for(int i=0;i<Rows;i++)
2114  {
2115  T factor(1+P[i]*(1+n0[i])/kbJoule/Tl);
2116  T rel(0.1);
2117 
2118  newE[i]=n0[i]*(rel*(factor-1)+1);
2119 
2120  if(factor>maxChange)
2121  maxChange=factor;
2122  if(factor<minChange)
2123  minChange=factor;
2124  }
2125 
2126  cout<<"Max Change: "<<maxChange<<" Min Change: "<<minChange<<endl;
2127 
2128  }
2129 
2130  private:
2131 
2132  void normalize()
2133  {
2134  _type1Collisions.multiplySelf(1./_maxDkl);
2135  _type1Collisions.multiplyOther(1./_maxPhi);
2136  _type2Collisions.multiplySelf(1./_maxDkl);
2137  _type2Collisions.multiplyOther(1./_maxPhi);
2138  }
2139 
2140  ScatteringKernel(const ScatteringKernel&);
2141  Kspace<T>& _kspace;
2142  KSConnectivity<T> _type1Collisions;
2143  KSConnectivity<T> _type2Collisions;
2144  T _maxPhi;
2145  T _maxDkl;
2146 
2147 };
2148 
2149 #endif
Array::zero
virtual void zero()
Definition: Array.h:281
ScatteringKernel::_type1Collisions
KSConnectivity< T > _type1Collisions
Definition: ScatteringKernel.h:2142
ScatteringKernel::_kspace
Kspace< T > & _kspace
Definition: ScatteringKernel.h:2141
Kspace
Definition: Kspace.h:28
ScatteringKernel::ScatteringKernel
ScatteringKernel(Kspace< T > &kspace)
Definition: ScatteringKernel.h:27
ScatteringKernel::VectorT3
Vector< T, 3 > VectorT3
Definition: ScatteringKernel.h:20
pmode::calce0
T calce0(T Tl)
Definition: pmode.h:88
kvol::getmode
Tmode & getmode(int n) const
Definition: kvol.h:44
kvol::getmodenum
int getmodenum()
Definition: kvol.h:43
ScatteringKernel::IterateToEquilibrium
ArrayBase * IterateToEquilibrium(const T Tl, const int totIts, const T tStep)
Definition: ScatteringKernel.h:1435
ScatteringKernel::Tksconn
KSConnectivity< T > Tksconn
Definition: ScatteringKernel.h:24
ScatteringKernel::ReadType2
void ReadType2(const string &NamePhonon2, const string &NamePhonon3, const T tol)
Definition: ScatteringKernel.h:232
ScatteringKernel::updateSource2
void updateSource2(const int c, TArray &S, TArray &dS)
Definition: ScatteringKernel.h:1003
ScatteringKernel::BArray
Array< bool > BArray
Definition: ScatteringKernel.h:19
ScatteringKernel::_type2Collisions
KSConnectivity< T > _type2Collisions
Definition: ScatteringKernel.h:2143
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Array
Definition: Array.h:14
pmode::calcde0dT
T calcde0dT(T Tl)
Definition: pmode.h:126
ScatteringKernel::ReadType1
void ReadType1(const string &NamePhonon2, const string &NamePhonon3, const T tol)
Definition: ScatteringKernel.h:35
ScatteringKernel::getTypeIIsource
void getTypeIIsource(const int c, TArray &S, TArray &dS, const bool correct)
Definition: ScatteringKernel.h:1240
ScatteringKernel::calculatePsi
ArrayBase * calculatePsi(const int totIts)
Definition: ScatteringKernel.h:1787
ScatteringKernel::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: ScatteringKernel.h:21
ScatteringKernel::updateSourceTermTest
void updateSourceTermTest(const T Tl)
Definition: ScatteringKernel.h:528
pmode
Definition: pmode.h:18
kvol
Definition: kvol.h:14
pmode::gettau
T gettau()
Definition: pmode.h:61
ScatteringKernel::updateSource
void updateSource(const int c, TArray &S)
Definition: ScatteringKernel.h:623
ScatteringKernel::IntArray
Array< int > IntArray
Definition: ScatteringKernel.h:18
ScatteringKernel::getTypeIsource
void getTypeIsource(const int c, TArray &S, TArray &dS, const bool correct)
Definition: ScatteringKernel.h:835
kvol::getdk3
T getdk3()
Definition: kvol.h:42
Array::getLength
int getLength() const
Definition: Array.h:87
ScatteringKernel::updateSourceTerm
void updateSourceTerm(const TArray &e, const TArray &w, TArray &S)
Definition: ScatteringKernel.h:466
ScatteringKernel::scatterPhonons
void scatterPhonons(const int c, const int totIts, TArray &C, TArray &B, const TArray &V, TArray &newE, const T cv)
Definition: ScatteringKernel.h:1949