SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
SPIgrid.cpp
Go to the documentation of this file.
1 #include "SPIgrid.hpp"
2 
3 namespace SPI{
4 
6  PetscInt factorial(
7  PetscInt n
8  ) {
9  PetscInt value = 1;
10  for(int i=1; i<=n; ++i) value *= i;
11  return value;
12  }
13 
15  SPIVec get_D_Coeffs(
16  SPIVec &s,
17  PetscInt d
18  ){
19  // get_D_coeffs function with inputs s,d
20  PetscInt N = s.rows;
21  SPIVec spow(N,"spow");
22  spow = ones(N);
23  SPIMat A(N,"A");
24  for(PetscInt i=0; i<N; i++){
25  //std::cout<<"i = "<<i<<std::endl;
26  if(i>=1){
27  //spow = (s^i); // doesn't work.... :(
28  spow *= s;
29  }
30  //SPI::draw(spow);
31  //spow.print();
32  for(PetscInt j=0; j<N; j++){
33  //std::cout<<"i,j = "<<i<<","<<j<<std::endl;
34  A(i,j,spow(j,PETSC_TRUE));
35  }
36  }
37  A();
38  //A.print();
39  SPIVec b(zeros(N),"b");
40  b(d,factorial(d));
41  b();
42  //b.print();
43  SPIVec x(N,"x");
44  x = (b/A);
45  //x();
46  //x.print();
47 
48  return x;
49  }
50 
52  SPIMat map_D(
53  SPIMat D,
54  SPIVec y,
55  PetscInt d,
56  PetscInt order
57  ){
58  // map_D
59  if(d==1){
60  //SPI::SPIVec dydxi(D*y);
61  //SPI::SPIMat dxidy(SPI::diag(1./(D*y)));
62  //dxidy.print();
63  SPIMat Dy(diag(1./(D*y))*D);
64  //(Dy*y).print();
65  return Dy;
66  }
67  else if(d==2){
68  SPIMat D1(set_D(y,1,order,PETSC_TRUE),"D1");
69  SPIVec dxidy(1./(D1*y),"dxidy");
70  SPIVec d2xidy2(-1.*(D*y)*(dxidy^3),"d2xidy2");
71  SPIMat Dy(diag(dxidy*dxidy)*D + (diag(d2xidy2)*D1),"Dy");
72  //(Dy*(y)).print();
73  return Dy;
74  }
75  else{
76  SPI::printf("--------------Something wrong with map_D call -------------");
77  exit(0);
78  }
79  }
80 
82  SPIMat set_D(
83  SPIVec &y,
84  PetscInt d,
85  PetscInt order,
86  PetscBool uniform
87  ){
88  //PetscInt d=1;
89  // set_D
90  //PetscInt order=4;
91  //PetscInt n=21;
92  PetscInt n=y.size();
93  SPIVec xi = (linspace(0.,1.,n));
94  PetscScalar h = xi(1,PETSC_TRUE)-xi(0,PETSC_TRUE);
95  SPIVec one(ones(n),"ones");
96  SPIMat I(eye(n),"I");
97  PetscInt N = order+d;
98  if(N>n) {
99  PetscErrorCode ierr=1;
100  CHKERRXX(ierr);
101  }
102  PetscInt Nm1 = N-1;
103  if ((d%2) != 0) Nm1 += 1; // increase for odd derivative
104  SPIVec s(arange(Nm1)-(Nm1-1)/2,"s"); // set stencil
105  PetscInt smax = s(s.rows-1,PETSC_TRUE).real();
106 
107  SPIVec Coeffs(get_D_Coeffs(s,d),"Coeffs");
108 
109  SPIMat D(n,"D");
110  D();
111  for(PetscInt i=0; i<s.rows; i++){
112  //diag_to_add.~SPIMat(); // destroy to free memory
113  PetscInt k=(PetscInt)s(i,PETSC_TRUE).real();
114  PetscInt nmk=n-std::abs(k);
115  D += diag(Coeffs(i,PETSC_TRUE)*ones(nmk),k);
116 
117  //SPI::diag(Coeffs(0,PETSC_TRUE)*SPI::one(10),-1).print()
118  }
119  //D.print();
120  // BCs
121  D.ierr = MatSetOption(D.mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);CHKERRXX(D.ierr);
122  for(PetscInt i=0; i<smax; i++){
123  // for ith row
124  s.~SPIVec(); // deallocate
125  if((d%2)!=0){// odd derivative
126  s = (arange(Nm1-1)-i); // stencil for shifted diff of order-1
127  }
128  else{
129  s = arange(Nm1)-i;// stencil for shifted diff of order-1
130  }
131  Coeffs.~SPIVec();
132  Coeffs = get_D_Coeffs(s,d);
133  D.zero_row(i);
134  for(PetscInt j=0; j<s.rows; j++){
135  PetscInt sj=s(j,PETSC_TRUE).real();
136  //SPI::printf("setting D(%d,%d)=%g",i,sj+i,Coeffs(j,PETSC_TRUE));
137  D(i,sj+i,Coeffs(j,PETSC_TRUE));
138  }
139  D();
140  // for -ith-1 row
141  s.~SPIVec(); // deallocate
142  if((d%2)!=0){// odd derivative
143  s = -(arange(Nm1-1)-i); // stencil for shifted diff of order-1
144  }
145  else{
146  s = -(arange(Nm1)-i);// stencil for shifted diff of order-1
147  }
148  Coeffs.~SPIVec();
149  Coeffs = get_D_Coeffs(s,d);
150  D.zero_row(D.rows-1-i);
151  for(PetscInt j=0; j<s.rows; j++){
152  PetscInt sj=s(j,PETSC_TRUE).real();
153  //SPI::printf("setting D(%d,%d)=%g",i,sj+i,Coeffs(j,PETSC_TRUE));
154  D(D.rows-1-i,D.cols-1+sj-i,Coeffs(j,PETSC_TRUE));
155  }
156  D();
157  }
158  D.ierr = MatSetOption(D.mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);CHKERRXX(D.ierr);
159  D();
160  D*=(1./std::pow(h,d));
161  D();
162  if(uniform==PETSC_FALSE){
163  SPIMat Dy(map_D(D,y,d,order));
164  return Dy;
165  }
166  else{
167  return D;
168  }
169  }
170 
172  SPIMat set_D_periodic(
173  SPIVec &y,
174  PetscInt d,
175  PetscInt order
176  ){
177  PetscInt n=y.size();
178  PetscScalar h = y(1,PETSC_TRUE)-y(0,PETSC_TRUE);
179  SPIVec one(ones(n),"ones");
180  SPIMat I(eye(n),"I");
181  PetscInt N = order+d;
182  if(N>n) {
183  PetscErrorCode ierr=1;
184  CHKERRXX(ierr);
185  }
186  PetscInt Nm1 = N-1;
187  if ((d%2) != 0) Nm1 += 1; // increase for odd derivative
188  SPIVec s(arange(Nm1)-(Nm1-1)/2,"s"); // set stencil
189  PetscInt smax = s(s.rows-1,PETSC_TRUE).real();
190 
191  SPIVec Coeffs(get_D_Coeffs(s,d),"Coeffs");
192 
193  SPIMat D(n,"D");
194  D();
195  for(PetscInt i=0; i<s.rows; i++){
196  //diag_to_add.~SPIMat(); // destroy to free memory
197  PetscInt k=(PetscInt)s(i,PETSC_TRUE).real();
198  PetscInt nmk=n-std::abs(k);
199  D += diag(Coeffs(i,PETSC_TRUE)*ones(nmk),k);
200  if(k>0){
201  D += diag(Coeffs(i,PETSC_TRUE)*ones(k),-nmk);
202  }
203  else if(k<0){
204  D += diag(Coeffs(i,PETSC_TRUE)*ones(-k),k+n);
205  }
206 
207  }
208  D();
209  D*=(1./std::pow(h,d));
210  D();
211  return D;
212  }
214  SPIVec set_FD_stretched_y(
215  PetscScalar y_max,
216  PetscInt ny,
217  PetscScalar delta
218  ){
219  //SPI::SPIVec y((y_max*(1.+(SPI::tanh(delta*((linspace(0.,y_max,ny)/y_max) - 1.))/tanh(delta)))));
220  //return y;
221  return y_max*(1.+(SPI::tanh(delta*((linspace(0.,y_max,ny)/y_max) - 1.))/tanh(delta)));
222  }
223 
225  SPIMat set_D_Chebyshev(
226  SPIVec &x,
227  PetscInt d,
228  PetscBool need_map
229  ){
230  if(need_map){
231  SPIMat D(map_D_Chebyshev(x,d));
232  return D;
233  }
234  else{
235  PetscInt N=x.rows-1;
236  //PetscInt order=-1;
237  SPIMat D(zeros(N+1,N+1));
238  SPIVec c(ones(N+1));
239  c(0,2.);
240  c(N,2.);
241  for(PetscInt j=0; j<N+1; j++){
242  PetscScalar cj=c(j,PETSC_TRUE);
243  PetscScalar xj=x(j,PETSC_TRUE);
244  for(PetscInt k=0; k<N+1; k++){
245  PetscScalar ck=c(k,PETSC_TRUE);
246  PetscScalar xk=x(k,PETSC_TRUE);
247  if(j!=k){
248  D(j,k,cj*pow(-1.,(PetscScalar)((double)(j+k)+0.0*PETSC_i)) / (ck*(xj-xk)));
249  }
250  else if((j==k) && ((j!=0) && (j!=N))){
251  D(j,k,-xj/(2.*(1.-(pow(xj,2)))));
252  }
253  else if((j==k) && (j==0)){
254  D(j,k,-1.0*(2.*pow((PetscScalar)((double)N+0.0*PETSC_i),2.)+1.)/6.); // make sure xj goes from -1 to 1
255  }
256  else if((j==k) && (j==N)){
257  D(j,k,1.0*(2.*pow((PetscScalar)((double)N+0.0*PETSC_i),2.)+1.)/6.);
258  }
259  else{
260  std::cout<<"you messed up in set_D_Chebyshev"<<std::endl;
261  exit(0);
262  }
263  }
264 
265  }
266  D();
267  if(d==2){
268  SPIMat D2(D*D,"Dyy");
269  return D2;
270  }else{
271  return D;
272  }
273  }
274  }
275 
276 
278  SPIMat map_D_Chebyshev(
279  SPIVec &x,
280  PetscInt d
281  ){
282  PetscInt N=x.rows-1;
283  SPIVec xi(set_Cheby_y(N+1));
284  if(d==1){
285  SPIMat D1(set_D_Chebyshev(xi,1,PETSC_FALSE));
286  SPIMat D(diag(1./(D1*x))*D1,"Dy");
287  return D;
288  }
289  else if(d==2){
290  SPIMat D1(set_D_Chebyshev(xi,1,PETSC_FALSE));
291  SPIMat D2(D1*D1);
292  SPIVec dxdxi(D1*x);
293  SPIVec dxidx(1./dxdxi);
294  SPIVec d2xidx2(-1.*(D2*x)*(dxidx^3));
295  SPIMat D((diag(dxidx^2))*D2 + (diag(d2xidx2)*D1),"Dyy");
296  return D;
297  }
298  else{
299  std::cout<<"order of derivative is not implemented in map_D_Chebyshev"<<std::endl;
300  exit(0);
301  }
302  }
303 
305  SPIVec set_Cheby_stretched_y(
306  PetscScalar y_max,
307  PetscInt ny,
308  PetscScalar yi
309  ){
310  SPIVec xi(set_Cheby_y(ny)); // default grid [-1,1]
311  PetscScalar a=yi*y_max/(y_max-2.*yi);
312  PetscScalar b=1.+2.*a/y_max;
313  SPIVec y(a*(1.+xi)/(b-xi),"y");
314  return y;
315  }
316 
318  SPIVec set_Cheby_mapped_y(
319  PetscScalar a,
320  PetscScalar b,
321  PetscInt ny
322  ){
323  SPIVec xi(set_Cheby_y(ny)); // default grid [-1,1]
324  SPIVec y(b*(xi+1.0)/2.0+a*(1.0-xi)/2.0,"y"); // grid from [a,b]
325  return y;
326  }
327 
329  SPIVec set_Cheby_y(
330  PetscInt ny
331  ){
332  //PetscScalar pi = 2.*std::acos(0.0);
333  SPIVec y(cos(PETSC_PI*arange((PetscScalar)ny-1.,-1.,-1.)/((PetscScalar)ny-1.)),"y");
334  return y;
335  }
336 
338  std::tuple<SPIMat,SPIMat> set_D_UltraS(
339  SPIVec &x,
340  PetscInt d
341  ){
342  PetscInt n=x.rows;
343  PetscScalar a=x(0,PETSC_TRUE);
344  PetscScalar b=x(n-1,PETSC_TRUE);
345  SPIMat dxidx(diag(ones(n)*(2.0/(b-a))),"dxi/dx"); // assuming it is made from set_Cheby_mapped_y (identity matrix if created from set_Cheby_y)
346  //SPIVec xi(cos(PETSC_PI*arange((PetscScalar)n-1.,-1.,-1.)/((PetscScalar)n-1.)),"y");
347  if(d==1){
348  SPIMat S0(diag(ones(n)*0.5) + diag(ones(n-2)*-0.5,2),"S0");
349  S0(0,0,1.0);
350  S0();
351  //std::cout<<"in set_D_UltraS"<<std::endl;
352  //S0.print();
353  SPIMat D1((diag(arange(1,n),1))*dxidx,"D1");
354  return std::make_tuple(S0,D1);
355  }
356  else if(d==2){
357  SPIMat S1(diag(1.0/(arange(n)+1)) + diag(-1.0/(arange(2,n)+1),2),"S1");
358  SPIMat D2((2.0*diag(arange(2,n),2))*(dxidx*dxidx),"D2");
359  return std::make_tuple(S1,D2);
360  }
361  else{
362  SPI::printf("Warning: d>2 not implemented in set_D_UltraS function");
363  exit(1);
364  }
365  }
367  std::tuple<SPIMat,SPIMat> set_T_That(
368  PetscInt n
369  ){
370  PetscInt N=n-1;
371  SPIVec xi(cos(PETSC_PI*arange((PetscScalar)n-1.,-1.,-1.)/((PetscScalar)n-1.)),"y");
372  SPIMat X(n,n);
373  SPIMat Ns(n,n); // X and Ns meshgrid for T and That creation
374  SPIVec ns(arange(n));
375  std::tie(X,Ns) = meshgrid(xi,ns);
376  SPIMat T(cos(Ns%acos(X)),"T");
377  SPIMat That; // initialized in next line
378  T.T(That); // initialize That and set it to be T.T()
379  SPIMat Thatcopy(That);
380  for(PetscInt i=0; i<n; ++i) That(i,0,Thatcopy(i,0,PETSC_TRUE)/2.0);
381  for(PetscInt i=0; i<n; ++i) That(i,n-1,Thatcopy(i,n-1,PETSC_TRUE)/2.0);
382  That();
383  Thatcopy=That;
384  for(PetscInt i=0; i<n; ++i) That(0,i,Thatcopy(0,i,PETSC_TRUE)/2.0);
385  for(PetscInt i=0; i<n; ++i) That(n-1,i,Thatcopy(n-1,i,PETSC_TRUE)/2.0);
386  That();
387  That *= (2.0/((PetscScalar)N));
388  return std::make_tuple(T,That);
389  }
390 
392  SPIVec set_Fourier_t(
393  PetscScalar T,
394  PetscInt nt
395  ){
396  PetscScalar step = T/(PetscScalar)nt;
397  SPIVec t(nt);
398  PetscScalar value=0.;
399  for(PetscInt i=0; i<nt; ++i){
400  t(i,value);
401  value += step;
402  }
403  t();
404  return t;
405  }
407  SPIMat set_D_Fourier(
408  SPIVec t,
409  PetscInt d
410  ){
411  //PetscScalar pi = PETSC_PI;
412  PetscScalar dt = t(1,PETSC_TRUE)-t(0,PETSC_TRUE); // uniform grid spacing
413  PetscInt npts=t.rows;
414  if(npts%2==0){ // only works with even number of grid points
415  PetscScalar nptss = (PetscScalar)npts;
416  PetscScalar T = t(npts-1,PETSC_TRUE) + dt;
417  SPIMat D;
418  SPIMat N,J;
419  SPIVec n(arange(0,npts));
420  SPIVec j(arange(0,npts));
421  std::tie(N,J) = meshgrid(n,j);
422  SPIMat NmJ(N-J,"NmJ");
423  D = ((PETSC_PI/T)*((-1.0)^(NmJ)))/tan((PETSC_PI/nptss)*(NmJ));
424  //std::cout<<"PETSC_PI/T = "<<PETSC_PI/T<<std::endl;
425  //((-1.0)^(NmJ)).print();
426  //((PETSC_PI/T)*((-1.0)^(NmJ))).print();
427  //tan((PETSC_PI/nptss)*(NmJ)).print();
428  //((((PETSC_PI/T)*((-1.0)^(NmJ))).real())/tan((PETSC_PI/nptss)*(NmJ)).real()).print();
429  //((((PETSC_PI/T)*((-1.0)^(NmJ))).real())*(1.0/(tan((PETSC_PI/nptss)*(NmJ)).real()))).print();
430  //std::cout<<"made it here"<<std::endl;
431  //(1.0/(tan((PETSC_PI/nptss)*(NmJ)).real())).print();
432  //std::cout<<"made it here2"<<std::endl;
433 
434  //D.print();
435  D.ierr = MatDiagonalSet(D.mat,zeros(npts).vec,INSERT_VALUES);CHKERRXX(D.ierr);
436  D();
437  D.real();
438  if(d==2) D = D*D; // then return Dyy
439  //D.print();
440  return D;
441  }
442  else{
443  exit(1);
444  }
445  }
446 
448  SPIgrid1D::SPIgrid1D(){}
450  SPIgrid1D::SPIgrid1D(
451  SPIVec &y,
452  std::string name,
453  gridtype _ytype
454  ){
455  // set grid name and type
456  this->name = name;
457  this->ytype = _ytype;
458  // set grid
459  this->set_grid(y);
460  // set respective derivatives
461  this->set_derivatives();
462  // set respective operators
463  this->set_operators();
464  }
465 
466 
468  void SPIgrid1D::print(){
469  SPI::printf("---------------- "+this->name+" start --------------------------");
470  if(this->flag_set_grid) {
471  //PetscInt ny2 = this->ny;
472  //SPI::printf("ny = %D",ny2);
473  this->y.print();
474  }
475  if(this->flag_set_derivatives){
476  this->Dy.print();
477  this->Dyy.print();
478  if(this->ytype==UltraS){
479  this->S0.print();
480  this->S1.print();
481  this->T.print();
482  this->That.print();
483  }
484  }
485  if(this->flag_set_operators){
486  this->O.print();
487  this->I.print();
488  }
489  SPI::printf("---------------- "+this->name+" done ---------------------------");
490  }
491 
493  void SPIgrid1D::set_grid(
494  SPIVec &y
495  ){
496  this->y=y;
497  this->y.name=std::string("y");
498  this->ny = y.rows;
499  this->flag_set_grid=PETSC_TRUE;
500  }
501 
503  void SPIgrid1D::set_derivatives(
504  PetscInt order
505  ){
506  if(this->ytype==FD){
507  this->Dy=set_D(this->y,1); // default of fourth order nonuniform grid
508  this->Dy.name=std::string("Dy");
509  this->Dyy=set_D(this->y,2); // default of fourth order nonuniform grid
510  this->Dyy.name=std::string("Dyy");
511  this->Dy.real(); // just take only real part
512  this->Dyy.real(); // just take only real part
513  this->flag_set_derivatives=PETSC_TRUE;
514  }
515  else if(this->ytype==FDperiodic){
516  this->Dy=set_D_periodic(this->y,1); // default of fourth order nonuniform grid
517  this->Dy.name=std::string("Dy");
518  this->Dyy=set_D_periodic(this->y,2); // default of fourth order nonuniform grid
519  this->Dyy.name=std::string("Dyy");
520  this->Dy.real(); // just take only real part
521  this->Dyy.real(); // just take only real part
522  this->flag_set_derivatives=PETSC_TRUE;
523  }
524  else if(this->ytype==Chebyshev){
525  std::tie(T,That) = set_T_That(y.rows); // get Chebyshev operators to take back and forth from physical space
526  this->Dy=set_D_Chebyshev(this->y,1,PETSC_TRUE); // default Chebyshev operator on non-uniform grid
527  this->Dy.name=std::string("Dy");
528  this->Dyy=set_D_Chebyshev(this->y,2,PETSC_TRUE); // default Chebyshev operator on non-uniform grid
529  this->Dyy.name=std::string("Dyy");
530  this->Dy.real(); // just take only real part
531  this->Dyy.real(); // just take only real part
532  this->flag_set_derivatives=PETSC_TRUE;
533  }
534  else if(this->ytype==Fourier){
535  this->Dy=set_D_Fourier(this->y,1); // default Fourier operator on uniform grid
536  this->Dy.name=std::string("Dy");
537  this->Dyy=set_D_Fourier(this->y,2); // default Chebyshev operator on uniform grid
538  this->Dyy.name=std::string("Dyy");
539  this->Dy.real(); // just take only real part
540  this->Dyy.real(); // just take only real part
541  this->flag_set_derivatives=PETSC_TRUE;
542  }
543  else if(this->ytype==UltraS){
544  std::tie(T,That) = set_T_That(y.rows); // get Chebyshev operators to take back and forth from physical space
545  std::tie(S0,Dy) = set_D_UltraS(this->y,1); // default UltraSpherical operators on non-uniform grid
546  std::tie(S1,Dyy) = set_D_UltraS(this->y,2); // default UltraSpherical operators on non-uniform grid
547  this->Dy.real(); // just take only real part
548  this->Dyy.real(); // just take only real part
549  this->Dy = S1*Dy; // make it output C^(2) coefficient space
550  this->S0.real(); // just take only real part
551  this->S1.real(); // just take only real part
552  this->T.real(); // just take only real part
553  this->That.real(); // just take only real part
554  this->S1S0That = this->S1*this->S0*this->That;
555  this->S0invS1inv = inv(this->S1*this->S0);
556 
557  this->Dy.name=std::string("Dy");
558  this->S0.name=std::string("S0");
559  this->S1.name=std::string("S1");
560  this->Dyy.name=std::string("Dyy");
561  this->T.name=std::string("T");
562  this->That.name=std::string("That");
563  this->S1S0That.name=std::string("S1*S0*That");
564  this->S0invS1inv.name=std::string("inv(S0)*inv(S1)");
565  //this->PS1S0That.name=std::string("PS1S0That");
566  this->flag_set_derivatives=PETSC_TRUE;
567  }
568  else{
569  SPI::printf("Warning this type of ytype grid is not implemented in SPIgrid1D.set_derivatives");
570  }
571  }
572 
574  void SPIgrid1D::set_operators(){
575  PetscInt m = Dy.rows;
576  PetscInt n = Dy.cols;
577  this->O = zeros(m,n);; // default Chebyshev operator on non-uniform grid
578  this->O.name = "zero";
579  if(this->ytype==UltraS){
580  this->P = block({
581  {zeros(2,ny-2),zeros(2,2)},
582  {eye(ny-2),zeros(ny-2,2)}
583  })();
584  P(0,ny-2,1.0);
585  P(1,ny-1,1.0);
586  P();
587  this->I = S1*S0; // get it up to second order coefficients (C^(2))
588  // permute such that the bottom two rows are now the top two rows (good for LU factorization and pivoting)
589  //PS1S0That = P*S1*S0*That;
590  //Dy = P*Dy;
591  //Dyy = P*Dyy;
592  //S0 = P*S0;
593  //S1 = P*S1;
594  }else{
595  this->I = eye(m); // default Chebyshev operator on non-uniform grid
596  }
597  this->I.name="eye";
598  this->flag_set_operators=PETSC_TRUE;
599  }
600 
602  SPIgrid1D::~SPIgrid1D(){
603  // destroy grid variables and reset flags
604  if (this->flag_set_grid){
605  this->y.~SPIVec();
606  this->flag_set_grid=PETSC_FALSE;
607  }
608  // destroy derivative operators and reset flags
609  if (this->flag_set_derivatives){
610  this->Dy.~SPIMat();
611  this->Dyy.~SPIMat();
612  this->S0.~SPIMat();
613  this->S1.~SPIMat();
614  this->T.~SPIMat();
615  this->That.~SPIMat();
616  this->S1S0That.~SPIMat();
617  this->S0invS1inv.~SPIMat();
618  this->flag_set_derivatives=PETSC_FALSE;
619  }
620  // destroy operators and reset flags
621  if (this->flag_set_operators){
622  this->O.~SPIMat();
623  this->I.~SPIMat();
624  this->P.~SPIMat();
625  this->FTinv.~SPIMat();
626  this->FT.~SPIMat();
627  this->Ihalf.~SPIMat();
628  this->Ihalfn.~SPIMat();
629  this->flag_set_operators=PETSC_FALSE;
630  }
631  }
633  SPIgrid2D::SPIgrid2D(
634  SPIVec &y,
635  SPIVec &t,
636  std::string name,
637  gridtype y_gridtype,
638  gridtype t_gridtype
639  ){
640  // set grid name and type
641  this->name = name;
642  this->ytype = y_gridtype;
643  this->ttype = t_gridtype;
644  // set grid
645  this->set_grid(y,t);
646  // set respective derivatives
647  this->set_derivatives();
648  // set respective operators
649  this->set_operators();
650  }
651 
652 
654  void SPIgrid2D::print(){
655  SPI::printf("---------------- "+this->name+" start --------------------------");
656  if(this->flag_set_grid) {
657  //PetscInt ny2 = this->ny;
658  //SPI::printf("ny = %D",ny2);
659  this->y.print();
660  this->t.print();
661  }
662  if(this->flag_set_derivatives){
663  this->Dy.print();
664  this->Dyy.print();
665  this->Dt.print();
666  }
667  if(this->flag_set_operators){
668  this->O.print();
669  this->I.print();
670  this->avgt.print();
671  }
672  SPI::printf("---------------- "+this->name+" done ---------------------------");
673  }
674 
676  void SPIgrid2D::set_grid(
677  SPIVec &y,
678  SPIVec &t
679  ){
680  this->y=y;
681  this->y.name=std::string("y");
682  this->ny = y.rows;
683  this->t=t;
684  this->t.name=std::string("t");
685  this->nt = t.rows;
686  this->grid1Dy.set_grid(y);
687  this->grid1Dt.set_grid(t);
688  this->flag_set_grid=PETSC_TRUE;
689  }
690 
692  void SPIgrid2D::set_derivatives(
693  PetscInt order
694  ){
695  //std::cout<<" grid3D set derivatives 1"<<std::endl;
696  //PetscInt ny=this->ny;
697  //PetscInt nt=this->nt;
698  // set 2D operators in wall-normal
699  grid1Dy.ytype=ytype;
700  grid1Dy.set_derivatives(order);
701  grid1Dt.ytype=ttype;
702  grid1Dt.set_derivatives(order);
703  // if(this->ytype==FD){
704  // this->Dy2D=set_D(this->y,1,order); // default of fourth order nonuniform grid
705  // this->Dy2D.name=std::string("Dy2D");
706  // this->Dyy2D=set_D(this->y,2,order); // default of fourth order nonuniform grid
707  // this->Dyy2D.name=std::string("Dyy2D");
708  // this->Dy2D.real(); // just take only real part
709  // this->Dyy2D.real(); // just take only real part
710  // this->flag_set_derivatives=PETSC_TRUE;
711  // }
712  // else if(this->ytype==FT){
713  // this->Dy2D=set_D_Fourier(this->y,1); // default Fourier operator on uniform grid
714  // this->Dy2D.name=std::string("Dy2D");
715  // this->Dyy2D=set_D_Fourier(this->y,2); // default Chebyshev operator on uniform grid
716  // this->Dyy2D.name=std::string("Dyy2D");
717  // this->Dy2D.real(); // just take only real part
718  // this->Dyy2D.real(); // just take only real part
719  // this->flag_set_derivatives=PETSC_TRUE;
720  // }
721  // else if(this->ytype==Chebyshev){
722  // this->Dy2D=set_D_Chebyshev(this->y,1,PETSC_TRUE); // default Chebyshev operator on non-uniform grid
723  // this->Dy2D.name=std::string("Dy2D");
724  // this->Dyy2D=set_D_Chebyshev(this->y,2,PETSC_TRUE); // default Chebyshev operator on non-uniform grid
725  // this->Dyy2D.name=std::string("Dyy2D");
726  // this->Dy2D.real(); // just take only real part
727  // this->Dyy2D.real(); // just take only real part
728  // this->flag_set_derivatives=PETSC_TRUE;
729  // }
730  // set 2D operators in time
731  // if(this->ttype==FD){
732  // this->Dt2D=set_D(this->t,1); // default of fourth order nonuniform grid
733  // this->Dt2D.name=std::string("Dt2D");
734  // this->Dt2D.real(); // just take only real part
735  // this->flag_set_derivatives=PETSC_TRUE;
736  // }
737  // else if(this->ttype==FT){
738  // this->Dt2D=set_D_Fourier(this->t,1); // default Fourier operator on uniform grid
739  // this->Dt2D.name=std::string("Dt2D");
740  // this->Dt2D.real(); // just take only real part
741  // this->flag_set_derivatives=PETSC_TRUE;
742  // }
743  // else if(this->ttype==Chebyshev){
744  // this->Dt2D=set_D_Chebyshev(this->y,1,PETSC_TRUE); // default Chebyshev operator on non-uniform grid
745  // this->Dt2D.name=std::string("Dt2D");
746  // this->Dt2D.real(); // just take only real part
747  // this->flag_set_derivatives=PETSC_TRUE;
748  // }
749  // set 3D operators
750  this->Dy = kron(eye(nt),this->grid1Dy.Dy);
751  this->Dy.name = "Dy";
752  this->Dyy = kron(eye(nt),this->grid1Dy.Dyy);
753  this->Dyy.name = "Dyy";
754  this->Dt = kron(this->grid1Dt.Dy,eye(ny));
755  this->Dt.name = "Dt";
756  // set 3D operators
757  if((ytype==Chebyshev) || (ytype==UltraS)){
758  this->T = kron(eye(nt),this->grid1Dy.T);
759  this->T.name = "T";
760  this->That = kron(eye(nt),this->grid1Dy.That);
761  this->That.name = "That";
762  }
763  if(ytype==UltraS){
764  this->S1S0That = kron(eye(nt),this->grid1Dy.S1S0That);
765  this->S1S0That.name = "S1S0That";
766  this->S0invS1inv = kron(eye(nt),this->grid1Dy.S0invS1inv);
767  this->S0invS1inv.name = "S0invS1inv";
768  }
769  }
770 
772  void SPIgrid2D::set_operators(){
773  //PetscInt ny = Dy.rows;
774  //PetscInt n = Dy.cols;
775  PetscInt n=(this->ny)*(this->nt);
776  this->O=zeros(n,n);; // zero matrix of size ny*nt x ny*nt
777  this->O.name="zero";
778  this->I=eye(n);
779  this->I.name="eye"; // identity matrix of size ny*nt x ny*nt
780  // set average in time operator
781  this->avgt=zeros(n,n);
782  PetscScalar val = 1.0/(PetscScalar)(this->nt);
783  //this->avgt /= (double)(this->nt);
784  for(PetscInt ti=0; ti<nt; ++ti){
785  for(PetscInt yi=0; yi<ny; ++yi){
786  PetscInt ii = ti*ny + yi;
787  for(PetscInt jj=yi; jj<n; jj+=this->ny){
788  this->avgt(ii,jj,val);
789  }
790  }
791  }
792  this->avgt();
793  this->grid1Dy.set_operators(); // in case they are needed
794  this->grid1Dt.set_operators(); // in case they are needed
795  // set operators for Fourier transform
796  std::tie(this->grid1Dt.FT,this->grid1Dt.FTinv,this->grid1Dt.Ihalf,this->grid1Dt.Ihalfn) = dft_dftinv_Ihalf_Ihalfn(this->nt);
797  this->FT = kron(this->grid1Dt.FT,eye(ny));
798  this->FTinv = kron(this->grid1Dt.FTinv,eye(ny));
799  this->Ihalf = kron(this->grid1Dt.Ihalf,eye(ny));
800  this->Ihalfn = kron(this->grid1Dt.Ihalfn,eye(ny));
801  this->FT4 = kron(eye(4),this->FT);
802  this->FTinv4 = kron(eye(4),this->FTinv);
803  this->Ihalf4 = kron(eye(4),this->Ihalf);
804  this->Ihalfn4 = kron(eye(4),this->Ihalfn);
805  this->flag_set_operators=PETSC_TRUE;
806  }
807 
809  SPIgrid2D::~SPIgrid2D(){
810  // destroy grid variables and reset flags
811  if (this->flag_set_grid){
812  this->y.~SPIVec();
813  this->t.~SPIVec();
814  this->flag_set_grid=PETSC_FALSE;
815  }
816  // destroy derivative operators and reset flags
817  if (this->flag_set_derivatives){
818  //this->Dy2D.~SPIMat();
819  //this->Dyy2D.~SPIMat();
820  //this->Dt2D.~SPIMat();
821  this->Dy.~SPIMat();
822  this->Dyy.~SPIMat();
823  this->Dt.~SPIMat();
824  this->T.~SPIMat();
825  this->That.~SPIMat();
826  this->S1S0That.~SPIMat();
827  this->S0invS1inv.~SPIMat();
828  this->flag_set_derivatives=PETSC_FALSE;
829  }
830  // destroy operators and reset flags
831  if (this->flag_set_operators){
832  this->O.~SPIMat();
833  this->I.~SPIMat();
834  this->avgt.~SPIMat();
835  this->FT.~SPIMat();
836  this->FTinv.~SPIMat();
837  this->Ihalf.~SPIMat();
838  this->Ihalfn.~SPIMat();
839  this->FT4.~SPIMat();
840  this->FTinv4.~SPIMat();
841  this->Ihalf4.~SPIMat();
842  this->Ihalfn4.~SPIMat();
843  this->flag_set_operators=PETSC_FALSE;
844  }
845  grid1Dy.~SPIgrid1D();
846  grid1Dt.~SPIgrid1D();
847  }
848 
850  SPIVec SPIVec1Dto2D(
851  SPIgrid2D &grid2D,
852  SPIVec &u
853  ){
854  PetscInt nu=u.rows; // should be equal to ny
855  PetscInt ny=grid2D.ny;
856  if(nu!=ny) exit(1);
857  PetscInt nt=grid2D.nt;
858  SPIVec U(ny*nt,u.name);
859  std::vector<PetscScalar> utmp(ny);
860  for(PetscInt j=0; j<ny; ++j){
861  utmp[j] = u(j,PETSC_TRUE);
862  }
863  for(PetscInt i=0; i<nt; ++i){
864  for(PetscInt j=0; j<ny; ++j){
865  //std::cout<<" i,j,j+i*ny = "<<i<<","<<j<<","<<j+i*ny<<std::endl;
866  U(j+i*ny,utmp[j]);
867  }
868  }
869  U();
870  return U;
871  }
873  PetscScalar integrate(
874  const SPIVec &a,
875  SPIgrid1D &grid
876  ){
877  if(grid.ytype==UltraS){ // if using UltraSpherical, then it is in Chebyshev coefficients
878  PetscInt ny=grid.y.rows;
879  PetscInt nyx=a.rows;
880  PetscInt ni=nyx/ny;
881  PetscScalar val=0.0;
882  SPIVec atmp(ny,"atmp");
883  for(PetscInt i=0; i<ni; ++i){
884  for(PetscInt j=0; j<ny; ++j){
885  atmp(j,a(ny*i + j,PETSC_TRUE));
886  }
887  atmp();
888  val += integrate_coeffs(atmp,grid.y);
889  }
890  return val;
891  }
892  else if(grid.ytype==Chebyshev){ // if using UltraSpherical, then it is in Chebyshev coefficients
893  PetscInt ny=grid.y.rows;
894  PetscInt nyx=a.rows;
895  PetscInt ni=nyx/ny;
896  PetscScalar val=0.0;
897  SPIVec atmp(ny,"atmp");
898  SPIVec atmp2(ny,"atmp");
899  for(PetscInt i=0; i<ni; ++i){
900  for(PetscInt j=0; j<ny; ++j){
901  atmp(j,a(ny*i + j,PETSC_TRUE));
902  }
903  atmp();
904  //((grid.That)*atmp).print();
905  atmp2 = ((grid.That)*atmp);
906  val += integrate_coeffs(atmp2,grid.y);
907  }
908  return val;
909  }
910  else{ // otherwise they are physical values, let's integrate using trapezoidal rule
911  PetscInt ny=grid.y.rows;
912  PetscInt nyx=a.rows;
913  PetscInt ni=nyx/ny;
914  PetscScalar val=0.0;
915  SPIVec atmp(ny,"atmp");
916  for(PetscInt i=0; i<ni; ++i){
917  for(PetscInt j=0; j<ny; ++j){
918  atmp(j,a(ny*i + j,PETSC_TRUE));
919  }
920  atmp();
921  val += trapz(atmp,grid.y);
922  }
923  return val;
924  }
925  }
927  PetscScalar integrate(
928  const SPIVec &a,
929  SPIgrid2D &grid
930  ){
931  if((grid.ytype==UltraS) && (grid.ttype==FD)){ // if using UltraSpherical, then it is in Chebyshev coefficients
932  PetscInt ny=grid.ny;
933  PetscInt nt=grid.nt;
934  PetscInt nytx=a.rows;
935  PetscInt ni=nytx/(ny*nt);
936  PetscScalar val=0.0;
937  PetscScalar val2=0.0;
938  SPIVec atmp(ny,"atmp");
939  //SPIVec atmp2(ny,"atmp2");
940  SPIVec atmpt(nt,"atmpt");
941  //SPIVec atmpt2(nt,"atmpt2");
942  for(PetscInt i=0; i<ni; ++i){
943  for(PetscInt j=0; j<nt; ++j){
944  for(PetscInt k=0; k<ny; ++k){
945  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
946  }
947  atmp();
948  //((grid.That)*atmp).print();
949  //atmp2 = ((grid.grid1Dy.That)*atmp);
950  val2 = integrate_coeffs(atmp,grid.y);
951  //val += val2;
952  atmpt(j,val2);
953  }
954  atmpt();
955  val2 = trapz(atmpt,grid.t);
956  val += val2;
957  }
958  return val;
959  }
960  else if((grid.ytype==UltraS) && (grid.ttype==Chebyshev)){ // if using UltraSpherical, then it is in Chebyshev coefficients
961  PetscInt ny=grid.ny;
962  PetscInt nt=grid.nt;
963  PetscInt nytx=a.rows;
964  PetscInt ni=nytx/(ny*nt);
965  PetscScalar val=0.0;
966  PetscScalar val2=0.0;
967  SPIVec atmp(ny,"atmp");
968  //SPIVec atmp2(ny,"atmp2");
969  SPIVec atmpt(nt,"atmpt");
970  SPIVec atmpt2(nt,"atmpt2");
971  for(PetscInt i=0; i<ni; ++i){
972  for(PetscInt j=0; j<nt; ++j){
973  for(PetscInt k=0; k<ny; ++k){
974  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
975  }
976  atmp();
977  //((grid.That)*atmp).print();
978  //atmp2 = ((grid.grid1Dy.That)*atmp);
979  val2 = integrate_coeffs(atmp,grid.y);
980  //val += val2;
981  atmpt(j,val2);
982  }
983  atmpt();
984  atmpt2 = ((grid.grid1Dy.That)*atmpt);
985  val2 = integrate_coeffs(atmpt,grid.t);
986  val += val2;
987  }
988  return val;
989  }
990  else if((grid.ytype==Chebyshev) && (grid.ttype==Chebyshev)){ // if using UltraSpherical, then it is in Chebyshev coefficients
991  PetscInt ny=grid.ny;
992  PetscInt nt=grid.nt;
993  PetscInt nytx=a.rows;
994  PetscInt ni=nytx/(ny*nt);
995  PetscScalar val=0.0;
996  PetscScalar val2=0.0;
997  SPIVec atmp(ny,"atmp");
998  SPIVec atmp2(ny,"atmp2");
999  SPIVec atmpt(nt,"atmpt");
1000  SPIVec atmpt2(nt,"atmpt2");
1001  for(PetscInt i=0; i<ni; ++i){
1002  for(PetscInt j=0; j<nt; ++j){
1003  for(PetscInt k=0; k<ny; ++k){
1004  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1005  }
1006  atmp();
1007  //((grid.That)*atmp).print();
1008  atmp2 = ((grid.grid1Dy.That)*atmp);
1009  val2 = integrate_coeffs(atmp2,grid.y);
1010  //val += val2;
1011  atmpt(j,val2);
1012  }
1013  atmpt();
1014  atmpt2 = ((grid.grid1Dt.That)*atmpt);
1015  val2 = integrate_coeffs(atmpt2,grid.t);
1016  val += val2;
1017  }
1018  //val2 = integrate_coeffs(atmp2,grid.y);
1019  return val;
1020  }
1021  else if((grid.ytype==UltraS) && (grid.ttype==Fourier)){
1022  PetscInt ny=grid.ny;
1023  PetscInt nt=grid.nt;
1024  PetscInt nytx=a.rows;
1025  PetscInt ni=nytx/(ny*nt);
1026  PetscScalar val=0.0;
1027  PetscScalar val2=0.0;
1028  SPIVec atmp(ny,"atmp");
1029  //SPIVec atmp2(ny,"atmp2");
1030  //SPIVec atmpt(nt,"atmpt");
1031  //SPIVec atmpt2(nt,"atmpt2");
1032  SPIVec atmpt(nt+1,"atmpt");
1033  SPIVec t(nt+1,"t");
1034  for(PetscInt i=0; i<nt; ++i){
1035  t(i,grid.t(i,PETSC_TRUE));
1036  }
1037  PetscScalar dt = grid.t(1,PETSC_TRUE) - grid.t(0,PETSC_TRUE);
1038  t(nt,grid.t(nt-1,PETSC_TRUE)+dt);
1039  t();
1040  for(PetscInt i=0; i<ni; ++i){
1041  for(PetscInt j=0; j<nt; ++j){
1042  for(PetscInt k=0; k<ny; ++k){
1043  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1044  }
1045  atmp();
1046  //((grid.That)*atmp).print();
1047  //atmp2 = ((grid.grid1Dy.That)*atmp);
1048  val2 = integrate_coeffs(atmp,grid.y);
1049  //val += val2;
1050  atmpt(j,val2);
1051  if(j==0){ atmpt(nt,val2); }
1052  }
1053  atmpt();
1054  //atmpt2 = ((grid.grid1Dt.That)*atmpt);
1055  //val2 = integrate_coeffs(atmpt2,grid.t);
1056  val2 = trapz(atmpt,t);
1057  val += val2;
1058  }
1059  std::cout<<"integrate val = "<<val<<std::endl;
1060  return val;
1061  }
1062  else if((grid.ytype==Chebyshev) && (grid.ttype==Fourier)){
1063  PetscInt ny=grid.ny;
1064  PetscInt nt=grid.nt;
1065  PetscInt nytx=a.rows;
1066  PetscInt ni=nytx/(ny*nt);
1067  PetscScalar val=0.0;
1068  PetscScalar val2=0.0;
1069  SPIVec atmp(ny,"atmp");
1070  SPIVec atmp2(ny,"atmp2");
1071  //SPIVec atmpt(nt,"atmpt");
1072  //SPIVec atmpt2(nt,"atmpt2");
1073  SPIVec atmpt(nt+1,"atmpt");
1074  SPIVec t(nt+1,"t");
1075  for(PetscInt i=0; i<nt; ++i){
1076  t(i,grid.t(i,PETSC_TRUE));
1077  }
1078  PetscScalar dt = grid.t(1,PETSC_TRUE) - grid.t(0,PETSC_TRUE);
1079  t(nt,grid.t(nt-1,PETSC_TRUE)+dt);
1080  t();
1081  for(PetscInt i=0; i<ni; ++i){
1082  for(PetscInt j=0; j<nt; ++j){
1083  for(PetscInt k=0; k<ny; ++k){
1084  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1085  }
1086  atmp();
1087  //((grid.That)*atmp).print();
1088  atmp2 = ((grid.grid1Dy.That)*atmp);
1089  val2 = integrate_coeffs(atmp2,grid.y);
1090  //val += val2;
1091  atmpt(j,val2);
1092  if(j==0){ atmpt(nt,val2); }
1093  }
1094  atmpt();
1095  //atmpt2 = ((grid.grid1Dt.That)*atmpt);
1096  //val2 = integrate_coeffs(atmpt2,grid.t);
1097  val2 = trapz(atmpt,t);
1098  val += val2;
1099  }
1100  std::cout<<"integrate val = "<<val<<std::endl;
1101  return val;
1102  }
1103  else if((grid.ytype==FD) && (grid.ttype==FD)){ // otherwise they are physical values, let's integrate using trapezoidal rule
1104  PetscInt ny=grid.ny;
1105  PetscInt nt=grid.nt;
1106  PetscInt nytx=a.rows;
1107  PetscInt ni=nytx/(ny*nt);
1108  PetscScalar val=0.0;
1109  PetscScalar val2=0.0;
1110  SPIVec atmp(ny,"atmp");
1111  //SPIVec atmp2(ny,"atmp2");
1112  SPIVec atmpt(nt,"atmpt");
1113  //SPIVec atmpt2(nt,"atmpt2");
1114  for(PetscInt i=0; i<ni; ++i){
1115  for(PetscInt j=0; j<nt; ++j){
1116  for(PetscInt k=0; k<ny; ++k){
1117  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1118  }
1119  atmp();
1120  //((grid.That)*atmp).print();
1121  //atmp2 = ((grid.grid1Dy.That)*atmp);
1122  val2 = trapz(atmp,grid.y);
1123  //val += val2;
1124  atmpt(j,val2);
1125  }
1126  atmpt();
1127  //atmpt2 = ((grid.grid1Dt.That)*atmpt);
1128  val2 = trapz(atmpt,grid.t);
1129  val += val2;
1130  }
1131  //val2 = integrate_coeffs(atmp2,grid.y);
1132  return val;
1133  }
1134  else if((grid.ytype==FD) && (grid.ttype==FDperiodic)){ // otherwise they are physical values, let's integrate using trapezoidal rule assuming periodic endpoint
1135  //std::cout<<" integrating FD and FDperiodic"<<std::endl;
1136  PetscInt ny=grid.ny;
1137  PetscInt nt=grid.nt;
1138  PetscInt nytx=a.rows;
1139  PetscInt ni=nytx/(ny*nt);
1140  PetscScalar val=0.0;
1141  PetscScalar val2=0.0;
1142  SPIVec atmp(ny,"atmp");
1143  //SPIVec atmp2(ny,"atmp2");
1144  SPIVec atmpt(nt+1,"atmpt");
1145  SPIVec t(nt+1,"t");
1146  for(PetscInt i=0; i<nt; ++i){
1147  t(i,grid.t(i,PETSC_TRUE));
1148  }
1149  PetscScalar dt = grid.t(1,PETSC_TRUE) - grid.t(0,PETSC_TRUE);
1150  t(nt,grid.t(nt-1,PETSC_TRUE)+dt);
1151  t();
1152  //SPIVec atmpt2(nt,"atmpt2");
1153  for(PetscInt i=0; i<ni; ++i){
1154  for(PetscInt j=0; j<nt; ++j){
1155  for(PetscInt k=0; k<ny; ++k){
1156  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1157  }
1158  atmp();
1159  //((grid.That)*atmp).print();
1160  //atmp2 = ((grid.grid1Dy.That)*atmp);
1161  val2 = trapz(atmp,grid.y);
1162  //val += val2;
1163  atmpt(j,val2);
1164  if(j==0){ atmpt(nt,val2); }
1165  }
1166  atmpt();
1167  //atmpt2 = ((grid.grid1Dt.That)*atmpt);
1168  val2 = trapz(atmpt,t);
1169  val += val2;
1170  }
1171  //val2 = integrate_coeffs(atmp2,grid.y);
1172  return val;
1173  }
1174  else if((grid.ytype==FD) && (grid.ttype==Fourier)){ // otherwise they are physical values, let's integrate using trapezoidal rule assuming periodic endpoint
1175  PetscInt ny=grid.ny;
1176  PetscInt nt=grid.nt;
1177  PetscInt nytx=a.rows;
1178  PetscInt ni=nytx/(ny*nt);
1179  PetscScalar val=0.0;
1180  PetscScalar val2=0.0;
1181  SPIVec atmp(ny,"atmp");
1182  //SPIVec atmp2(ny,"atmp2");
1183  SPIVec atmpt(nt+1,"atmpt");
1184  SPIVec t(nt+1,"t");
1185  for(PetscInt i=0; i<nt; ++i){
1186  t(i,grid.t(i,PETSC_TRUE));
1187  }
1188  PetscScalar dt = grid.t(1,PETSC_TRUE) - grid.t(0,PETSC_TRUE);
1189  t(nt,grid.t(nt-1,PETSC_TRUE)+dt);
1190  t();
1191  //SPIVec atmpt2(nt,"atmpt2");
1192  for(PetscInt i=0; i<ni; ++i){
1193  for(PetscInt j=0; j<nt; ++j){
1194  for(PetscInt k=0; k<ny; ++k){
1195  atmp(k,a(ny*nt*i + ny*j + k,PETSC_TRUE));
1196  }
1197  atmp();
1198  //((grid.That)*atmp).print();
1199  //atmp2 = ((grid.grid1Dy.That)*atmp);
1200  val2 = trapz(atmp,grid.y);
1201  //val += val2;
1202  atmpt(j,val2);
1203  if(j==0){ atmpt(nt,val2); }
1204  }
1205  atmpt();
1206  //atmpt2 = ((grid.grid1Dt.That)*atmpt);
1207  val2 = trapz(atmpt,t);
1208  val += val2;
1209  }
1210  //val2 = integrate_coeffs(atmp2,grid.y);
1211  return val;
1212  }
1213  else{
1214  std::cout<<"integrate mixed types not available"<<std::endl;
1215  exit(1);
1216  return 0;
1217  }
1218  }
1219  /* \brief project using inner product for Gram-Schmidt process */
1220  SPIVec proj(
1221  SPIVec &u,
1222  SPIVec &v,
1223  SPIgrid1D &grid
1224  ){
1225  if(grid.ytype==UltraS){
1226  PetscInt n = u.rows/grid.y.rows;
1227  SPIVec utmp, vtmp;
1228  SPIMat T, That;
1229  if(n==1){
1230  T = grid.T;
1231  That = grid.That;
1232  utmp=T*u;
1233  vtmp=T*v;
1234  }else if(n==2){
1235  T = block({
1236  {grid.T,grid.O},
1237  {grid.O,grid.T},
1238  })();
1239  That = block({
1240  {grid.That,grid.O},
1241  {grid.O,grid.That},
1242  })();
1243  utmp=T*u;
1244  vtmp=T*v;
1245  }else if(n==4){
1246  T = block({
1247  {grid.T,grid.O,grid.O,grid.O},
1248  {grid.O,grid.T,grid.O,grid.O},
1249  {grid.O,grid.O,grid.T,grid.O},
1250  {grid.O,grid.O,grid.O,grid.T},
1251  })();
1252  That = block({
1253  {grid.That,grid.O,grid.O,grid.O},
1254  {grid.O,grid.That,grid.O,grid.O},
1255  {grid.O,grid.O,grid.That,grid.O},
1256  {grid.O,grid.O,grid.O,grid.That},
1257  })();
1258  utmp=T*u;
1259  vtmp=T*v;
1260  }
1261  return (integrate(That*(conj(utmp)*vtmp),grid)/integrate(That*(conj(utmp)*utmp),grid)) * u;
1262  }else{
1263  return (integrate(conj(u)*v,grid)/integrate(conj(u)*u,grid)) * u;
1264  }
1265  }
1266  /* \brief project using inner product for Gram-Schmidt process */
1267  SPIVec proj(
1268  SPIVec &u,
1269  SPIVec &v,
1270  SPIgrid2D &grid
1271  ){
1272  if(grid.ytype==UltraS){
1273  std::cout<<"in proj 1"<<std::endl;
1274  PetscInt n = u.rows/(grid.y.rows*grid.t.rows);
1275  SPIVec utmp, vtmp;
1276  SPIMat T, That;
1277  if(n==1){
1278  std::cout<<"in proj 2"<<std::endl;
1279  T = grid.T;
1280  That = grid.That;
1281  utmp=T*u;
1282  vtmp=T*v;
1283  }else if(n==2){
1284  std::cout<<"in proj 3"<<std::endl;
1285  T = block({
1286  {grid.T,grid.O},
1287  {grid.O,grid.T},
1288  })();
1289  That = block({
1290  {grid.That,grid.O},
1291  {grid.O,grid.That},
1292  })();
1293  utmp=T*u;
1294  vtmp=T*v;
1295  }else if(n==4){
1296  std::cout<<"in proj 4"<<std::endl;
1297  T = block({
1298  {grid.T,grid.O,grid.O,grid.O},
1299  {grid.O,grid.T,grid.O,grid.O},
1300  {grid.O,grid.O,grid.T,grid.O},
1301  {grid.O,grid.O,grid.O,grid.T},
1302  })();
1303  std::cout<<"in proj 5"<<std::endl;
1304  That = block({
1305  {grid.That,grid.O,grid.O,grid.O},
1306  {grid.O,grid.That,grid.O,grid.O},
1307  {grid.O,grid.O,grid.That,grid.O},
1308  {grid.O,grid.O,grid.O,grid.That},
1309  })();
1310  std::cout<<"in proj 6"<<std::endl;
1311  utmp=T*u;
1312  std::cout<<"in proj 7"<<std::endl;
1313  vtmp=T*v;
1314  std::cout<<"in proj 8"<<std::endl;
1315  }
1316  std::cout<<"in proj 9"<<std::endl;
1317  return (integrate(That*(conj(utmp)*vtmp),grid)/integrate(That*(conj(utmp)*utmp),grid)) * u;
1318  }else{
1319  return (integrate(conj(u)*v,grid)/integrate(conj(u)*u,grid)) * u;
1320  }
1321  }
1322  /* \brief orthogonalize a basis dense matrix from an array of vec using Gram-Schmidt */
1323  std::vector<SPIVec> orthogonalize(
1324  std::vector<SPIVec> &x,
1325  SPIgrid1D &grid
1326  ){
1327  //PetscInt m=x[0].rows; // number of rows
1328  PetscInt n=x.size(); // number of columns
1329  //SPIMat Q(m,n,"Q");
1330  std::vector<SPIVec> qi(n);
1331  // copy x[i]
1332  for(PetscInt i=0; i<n; ++i){
1333  qi[i] = x[i];
1334  }
1335  // project
1336  for(PetscInt i=0; i<n; ++i){
1337  if(i>0){
1338  for(PetscInt j=1; j<=i; ++j){
1339  //qi[i] -= proj(qi[j-1],x[i],grid); // if using Classical Gram-Schmidt orthogonalization
1340  qi[i] -= proj(qi[j-1],qi[i],grid); // if using Modified Gram-Schmidt orthogonalization
1341  //std::cout<<" projecting i,j = "<<i<<","<<j<<std::endl;
1342  }
1343  }
1344  }
1345  // normalize
1346  if(grid.ytype==UltraS){
1347  SPIVec qtmp;
1348  PetscInt ni = qi[0].rows/(grid.y.rows);
1349  if(ni==1){
1350  for(PetscInt i=0; i<n; ++i){
1351  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1352  qtmp = SPI::abs(grid.T*qi[i]);
1353  qtmp = grid.That*(qtmp*qtmp);
1354  qi[i] /= sqrt(integrate(qtmp,grid));
1355  }
1356  }
1357  else if(ni==2){
1358  SPIMat T(block({
1359  {grid.T,grid.O},
1360  {grid.O,grid.T},
1361  })(),"T");
1362  SPIMat That(block({
1363  {grid.That,grid.O},
1364  {grid.O,grid.That},
1365  })(),"That");
1366  for(PetscInt i=0; i<n; ++i){
1367  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1368  qtmp = SPI::abs(T*qi[i]);
1369  qtmp = That*(qtmp*qtmp);
1370  qi[i] /= sqrt(integrate(qtmp,grid));
1371  }
1372  }
1373  else if(ni==4){
1374  SPIMat T(block({
1375  {grid.T,grid.O,grid.O,grid.O},
1376  {grid.O,grid.T,grid.O,grid.O},
1377  {grid.O,grid.O,grid.T,grid.O},
1378  {grid.O,grid.O,grid.O,grid.T},
1379  })(),"T");
1380  SPIMat That(block({
1381  {grid.That,grid.O,grid.O,grid.O},
1382  {grid.O,grid.That,grid.O,grid.O},
1383  {grid.O,grid.O,grid.That,grid.O},
1384  {grid.O,grid.O,grid.O,grid.That},
1385  })(),"That");
1386  for(PetscInt i=0; i<n; ++i){
1387  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1388  qtmp = SPI::abs(T*qi[i]);
1389  qtmp = That*(qtmp*qtmp);
1390  qi[i] /= sqrt(integrate(qtmp,grid));
1391  }
1392  }
1393  else{
1394  SPI::printf("orthogonalize not implemented yet");
1395  exit(1);
1396  }
1397  }
1398  else{
1399  for(PetscInt i=0; i<n; ++i){
1400  qi[i] /= sqrt(integrate(SPI::abs(qi[i])^2,grid));
1401  }
1402  }
1403  return qi;
1404  }
1405  /* \brief orthogonalize a basis dense matrix from an array of vec using SLEPc BV */
1406  std::vector<SPIVec> orthogonalize(
1407  std::vector<SPIVec> &x,
1408  SPIgrid2D &grid
1409  ){
1410  //std::cout<<"in orthogonalize 1"<<std::endl;
1411  //PetscInt m=x[0].rows; // number of rows
1412  PetscInt n=x.size(); // number of columns
1413  //SPIMat Q(m,n,"Q");
1414  std::vector<SPIVec> qi(n);
1415  // copy x[i]
1416  for(PetscInt i=0; i<n; ++i){
1417  qi[i] = x[i];
1418  }
1419  //std::cout<<"in orthogonalize 2"<<std::endl;
1420  // project
1421  for(PetscInt i=0; i<n; ++i){
1422  if(i>0){
1423  for(PetscInt j=1; j<=i; ++j){
1424  //qi[i] -= proj(qi[j-1],x[i],grid);
1425  //qi[i] -= proj(qi[j-1],x[i],grid); // if using Classical Gram-Schmidt orthogonalization
1426  qi[i] -= proj(qi[j-1],qi[i],grid); // if using Modified Gram-Schmidt orthogonalization
1427  //std::cout<<" projecting i,j = "<<i<<","<<j<<std::endl;
1428  }
1429  }
1430  }
1431  //std::cout<<"in orthogonalize 3"<<std::endl;
1432  // normalize
1433  if(grid.ytype==UltraS){
1434  SPIVec qtmp;
1435  PetscInt ni = qi[0].rows/(grid.y.rows*grid.t.rows);
1436  if(ni==1){
1437  for(PetscInt i=0; i<n; ++i){
1438  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1439  qtmp = SPI::abs(grid.T*qi[i]);
1440  qtmp = grid.That*(qtmp*qtmp);
1441  qi[i] /= sqrt(integrate(qtmp,grid));
1442  }
1443  }
1444  else if(ni==2){
1445  SPIMat T(block({
1446  {grid.T,grid.O},
1447  {grid.O,grid.T},
1448  })(),"T");
1449  SPIMat That(block({
1450  {grid.That,grid.O},
1451  {grid.O,grid.That},
1452  })(),"That");
1453  for(PetscInt i=0; i<n; ++i){
1454  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1455  qtmp = SPI::abs(T*qi[i]);
1456  qtmp = That*(qtmp*qtmp);
1457  qi[i] /= sqrt(integrate(qtmp,grid));
1458  }
1459  }
1460  else if(ni==4){
1461  //std::cout<<"in orthogonalize 4"<<std::endl;
1462  SPIMat T(block({
1463  {grid.T,grid.O,grid.O,grid.O},
1464  {grid.O,grid.T,grid.O,grid.O},
1465  {grid.O,grid.O,grid.T,grid.O},
1466  {grid.O,grid.O,grid.O,grid.T},
1467  })(),"T");
1468  //std::cout<<"in orthogonalize 5"<<std::endl;
1469  SPIMat That(block({
1470  {grid.That,grid.O,grid.O,grid.O},
1471  {grid.O,grid.That,grid.O,grid.O},
1472  {grid.O,grid.O,grid.That,grid.O},
1473  {grid.O,grid.O,grid.O,grid.That},
1474  })(),"That");
1475  //std::cout<<"in orthogonalize 6"<<std::endl;
1476  for(PetscInt i=0; i<n; ++i){
1477  //std::cout<<" norm(q) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid));
1478  qtmp = SPI::abs(T*qi[i]);
1479  qtmp = That*(qtmp*qtmp);
1480  qi[i] /= sqrt(integrate(qtmp,grid));
1481  }
1482  //std::cout<<"in orthogonalize 7"<<std::endl;
1483  }
1484  else{
1485  SPI::printf("orthogonalize not implemented yet");
1486  exit(1);
1487  }
1488  }
1489  else{
1490  for(PetscInt i=0; i<n; ++i){
1491  //std::cout<<"sqrt(integrate(SPI::abs(qi[i])^2,grid)) = "<<sqrt(integrate(SPI::abs(qi[i])^2,grid))<<std::endl;
1492  //std::cout<<"sqrt(integrate(SPI::abs(qi[i])*SPI::abs(qi[i]),grid)) = "<<sqrt(integrate(SPI::abs(qi[i])*SPI::abs(qi[i]),grid))<<std::endl;
1493  //std::cout<<"abs(qi[i])(1),qi[i](1) = "<<SPI::abs(qi[i])(1,PETSC_TRUE)<<", "<<qi[i](1,PETSC_TRUE)<<std::endl;
1494  //std::cout<<"(abs(qi[i])*abs(qi[i]))(1) = "<<(SPI::abs(qi[i])*SPI::abs(qi[i]))(1,PETSC_TRUE)<<std::endl;
1495  //std::cout<<"sum(abs(qi[i])*abs(qi[i])) = "<<SPI::sum(SPI::abs(qi[i])*SPI::abs(qi[i]))<<std::endl;
1496  //std::cout<<"integrate(abs(qi[i])*abs(qi[i])) = "<<SPI::integrate(SPI::abs(qi[i])*SPI::abs(qi[i]),grid)<<std::endl;
1497  //(SPI::abs(qi[i])*SPI::abs(qi[i])).print();
1498  qi[i] /= sqrt(integrate(SPI::abs(qi[i])^2,grid));
1499  }
1500  }
1501  //std::cout<<"in orthogonalize 8"<<std::endl;
1502  //qi[0].print();
1503  return qi;
1504  }
1505  /* \brief create matrix to interpolate from grid1 to grid2 \returns out matrix such that u2(grid2.y) = out*u1(grid1.y) */
1506  SPIMat interp1D_Mat(
1507  SPIgrid1D &grid1,
1508  SPIgrid1D &grid2
1509  ){
1510  PetscInt n1=grid1.ny;
1511  PetscInt n2=grid2.ny;
1512  SPIMat I_n2Xn1(n2,n1,"I");
1513  for(PetscInt i=0; i<n2; i++){
1514  PetscScalar y2i = grid2.y(i,PETSC_TRUE);
1515  PetscBool flag=PETSC_TRUE;
1516  for(PetscInt j=0; j<n1; j++){
1517  PetscScalar y1j = grid1.y(j,PETSC_TRUE);
1518  if(y1j.real()>y2i.real()){
1519  if(flag){
1520  I_n2Xn1(i,j-1,1.0);
1521  flag = PETSC_FALSE;
1522  }
1523  }
1524  else if(y1j==y2i){
1525  if(flag){
1526  I_n2Xn1(i,j,1.0);
1527  flag = PETSC_FALSE;
1528  }
1529  }
1530  }
1531  }
1532  I_n2Xn1();
1533  SPIMat Deltay(diag(grid2.y - (I_n2Xn1*grid1.y)),"Deltay");
1534  SPIMat interp_n1_to_n2(I_n2Xn1 + (Deltay*(I_n2Xn1*grid1.Dy)) + 0.5*(Deltay*Deltay*(I_n2Xn1*grid1.Dyy)),"interp");
1535  return interp_n1_to_n2;
1536  }
1537  /* \brief create matrix to interpolate from grid1 to grid2 \returns out matrix such that u2(grid2.y) = out*u1(grid1.y) */
1538  SPIMat interp1D_Mat(
1539  SPIVec &y1,
1540  SPIVec &y2
1541  ){
1542  PetscInt n1=y1.rows;
1543  PetscInt n2=y2.rows;
1544  SPIMat Dy1(set_D(y1,1),"Dy");
1545  SPIMat Dyy1(set_D(y1,2),"Dy");
1546  SPIMat I_n2Xn1(n2,n1,"I"); // get nearest neighbor below the value (always positive dy)
1547  // nearest neighbor below the value (TODO improve by finding nearest neighbor)
1548  for(PetscInt i=0; i<n2; i++){
1549  PetscScalar y2i = y2(i,PETSC_TRUE);
1550  PetscBool flag=PETSC_TRUE;
1551  for(PetscInt j=0; j<n1; j++){
1552  PetscScalar y1j = y1(j,PETSC_TRUE);
1553  if(y1j.real()>y2i.real()){
1554  if(flag){
1555  I_n2Xn1(i,j-1,1.0);
1556  flag = PETSC_FALSE;
1557  }
1558  }
1559  else if(y1j==y2i){
1560  if(flag){
1561  I_n2Xn1(i,j,1.0);
1562  flag = PETSC_FALSE;
1563  }
1564  }
1565  }
1566  }
1567  I_n2Xn1();
1568  SPIMat Deltay(diag(y2 - (I_n2Xn1*y1)),"Deltay");
1569  SPIMat interp_n1_to_n2(I_n2Xn1 + (Deltay*(I_n2Xn1*Dy1)) + 0.5*(Deltay*Deltay*(I_n2Xn1*Dyy1)),"interp");
1570  return interp_n1_to_n2;
1571  }
1572  /* \brief create a discrete fourier transform transformation operator \returns DFT matrix*/
1573  SPIMat dft(
1574  PetscInt nt
1575  ){
1576  PetscScalar Nt = (PetscScalar)nt;
1577  PetscScalar e = (PetscExpScalar(-2.0*PETSC_PI*PETSC_i/Nt));
1578  SPIMat W(nt,nt,"W");
1579  for(PetscInt i=0; i<nt; ++i){
1580  for(PetscInt j=0; j<nt; ++j){
1581  W(i,j,PetscPowScalar(e,i*j));
1582  }
1583  }
1584  W(); // assemble
1585  W /= Nt;
1586  return W;
1587  }
1588  /* \brief create a discrete fourier transform transformation operator and it's associated inverse along with positive and negative identity like operators\returns DFT matrix, inv(DFT), and I_half, and I_halfn (Ihalf + Ihalfn = eye(nt))*/
1589  std::tuple<SPIMat,SPIMat,SPIMat,SPIMat> dft_dftinv_Ihalf_Ihalfn(
1590  PetscInt nt
1591  ){
1592  PetscScalar Nt = (PetscScalar)nt;
1593  SPIMat FT(dft(nt),"FT");
1594  SPIMat FTinv(nt,nt,"inv(FT)");
1595  FT.H(FTinv);
1596  FTinv *= Nt;
1597  SPIMat Ihalf(nt,nt,"Ihalf");
1598  SPIMat Ihalfn(nt,nt,"Ihalf");
1599  PetscInt Nthalf = nt/2;
1600  if(nt%2 == 0) // nt is even
1601  Nthalf = nt/2;
1602  else // odd
1603  Nthalf = nt/2 + 1;
1604  for(PetscInt i=0; i<nt; ++i){
1605  if(i==0){
1606  Ihalf(i,i,0.5);
1607  Ihalfn(i,i,0.5);
1608  }
1609  else if(i < Nthalf){
1610  Ihalf(i,i,1.0);
1611  }
1612  else if(i >= Nthalf){
1613  Ihalfn(i,i,1.0);
1614  }
1615  }
1616  // assemble identity operators for splitting positive and negative wavenumbers
1617  Ihalf();
1618  Ihalfn();
1619 
1620  return std::make_tuple(FT,FTinv,Ihalf,Ihalfn);
1621  }
1622 
1623 }
SPI::SPIgrid2D::FTinv
SPIMat FTinv
inverse Fourier Transform from wavenumber space to physical space for time
Definition: SPIgrid.hpp:90
SPI::SPIgrid2D::S0invS1inv
SPIMat S0invS1inv
inverse of S0^-1 * S1^-1
Definition: SPIgrid.hpp:88
SPI::FDperiodic
@ FDperiodic
finite difference grid on periodic grid
Definition: SPIgrid.hpp:27
SPI::get_D_Coeffs
SPIVec get_D_Coeffs(SPIVec &s, PetscInt d)
get the coefficients of the given stencil.
Definition: SPIgrid.cpp:15
SPI::SPIgrid2D::Ihalf
SPIMat Ihalf
positive wavenumbers from FT
Definition: SPIgrid.hpp:91
SPI::SPIMat::print
PetscInt print()
print mat to screen using PETSC_VIEWER_STDOUT_WORLD
Definition: SPIMat.cpp:498
SPI::SPIgrid1D::FT
SPIMat FT
Fourier Transform operator.
Definition: SPIgrid.hpp:60
SPI::SPIgrid2D::FT
SPIMat FT
Fourier Transform operator from physical space to wavenumber space for time.
Definition: SPIgrid.hpp:89
SPI::SPIMat::rows
PetscInt rows
number of rows in mat
Definition: SPIMat.hpp:18
SPI::SPIgrid1D::Ihalfn
SPIMat Ihalfn
negative wavenumbers from FT
Definition: SPIgrid.hpp:63
SPI::SPIgrid1D::Dyy
SPIMat Dyy
2nd derivative operator with respect to y
Definition: SPIgrid.hpp:52
SPI::integrate
PetscScalar integrate(const SPIVec &a, SPIgrid1D &grid)
integrate a vector of chebyshev Coefficients on a physical grid
Definition: SPIgrid.cpp:873
SPI::SPIgrid2D::flag_set_operators
PetscBool flag_set_operators
flag if set_operators has been executed
Definition: SPIgrid.hpp:117
SPI::SPIgrid1D::flag_set_grid
PetscBool flag_set_grid
flag if set_grid has been executed
Definition: SPIgrid.hpp:67
SPI::SPIMat::mat
Mat mat
petsc Mat data
Definition: SPIMat.hpp:28
SPI::SPIgrid1D
Class to contain various grid parameters.
Definition: SPIgrid.hpp:35
SPI::SPIgrid1D::T
SPIMat T
Chebyshev operator taking it from Chebyshev coefficients to physical space.
Definition: SPIgrid.hpp:58
SPI::SPIVec::real
SPIVec & real()
take the real part of the vector
Definition: SPIVec.cpp:449
SPI::SPIgrid2D::~SPIgrid2D
~SPIgrid2D()
destructor of saved SPIVec and SPIMat
Definition: SPIgrid.cpp:809
SPI::SPIgrid1D::SPIgrid1D
SPIgrid1D()
constructor with at no arguments
Definition: SPIgrid.cpp:448
SPI::SPIgrid1D::y
SPIVec y
grid
Definition: SPIgrid.hpp:47
SPI::SPIgrid1D::Dy
SPIMat Dy
1st derivative operator with respect to y
Definition: SPIgrid.hpp:51
SPI::Chebyshev
@ Chebyshev
Chebyshev collocated grid.
Definition: SPIgrid.hpp:29
SPI::arange
SPIVec arange(const PetscScalar begin, const PetscScalar end, const PetscScalar stepsize=1)
return a range of number of equally spaced points between begin and end by step size step
Definition: SPIVec.cpp:667
SPI::SPIMat::name
std::string name
Matrix name.
Definition: SPIMat.hpp:33
SPI::FD
@ FD
finite difference grid
Definition: SPIgrid.hpp:26
SPI::SPIVec
Definition: SPIVec.hpp:13
SPI::SPIgrid2D::flag_set_derivatives
PetscBool flag_set_derivatives
flag if set_derivatives has been executed
Definition: SPIgrid.hpp:116
SPI::printf
PetscInt printf(std::string msg,...)
Definition: SPIprint.cpp:5
SPI::Fourier
@ Fourier
Fourier transform collocated grid.
Definition: SPIgrid.hpp:28
SPI::SPIgrid2D::O
SPIMat O
zero matrix same size as derivative operators of size ny*nt x ny*nt
Definition: SPIgrid.hpp:111
SPI::set_D_periodic
SPIMat set_D_periodic(SPIVec &y, PetscInt d, PetscInt order=4)
set the derivative operator for the proper periodic grid assuming uniform discretization.
Definition: SPIgrid.cpp:172
SPI
Definition: SPIbaseflow.hpp:16
SPI::SPIgrid2D::y
SPIVec y
grid for wall-normal dimension
Definition: SPIgrid.hpp:83
SPI::conj
SPIVec conj(const SPIVec &A)
Definition: SPIVec.cpp:622
SPI::SPIgrid2D::print
void print()
saves grid to internal grid
Definition: SPIgrid.cpp:654
SPI::set_D
SPIMat set_D(SPIVec &y, PetscInt d, PetscInt order=4, PetscBool uniform=PETSC_FALSE)
set the derivative operator for the proper y grid if uniform=false. Uses map_D function
Definition: SPIgrid.cpp:82
SPI::ones
SPIVec ones(const PetscInt rows)
Definition: SPIVec.cpp:605
SPI::SPIMat::ierr
PetscErrorCode ierr
ierr for various routines and operators
Definition: SPIMat.hpp:29
SPI::kron
SPIMat kron(const SPIMat &A, const SPIMat &B)
set kronecker inner product of two matrices
Definition: SPIMat.cpp:780
SPI::SPIgrid1D::set_derivatives
void set_derivatives(PetscInt order=4)
sets derivatives Dy and Dyy using saved grid
Definition: SPIgrid.cpp:503
SPI::SPIgrid2D::ny
PetscInt ny
Definition: SPIgrid.hpp:75
SPI::dft_dftinv_Ihalf_Ihalfn
std::tuple< SPIMat, SPIMat, SPIMat, SPIMat > dft_dftinv_Ihalf_Ihalfn(PetscInt nt)
Definition: SPIgrid.cpp:1589
SPI::SPIMat::H
PetscInt H(SPIMat &A)
A = Hermitian Transpose(*this.mat) operation with initialization of A (tranpose and complex conjugate...
Definition: SPIMat.cpp:415
SPI::SPIVec::print
PetscInt print()
Definition: SPIVec.cpp:470
SPI::UltraS
@ UltraS
UltraSpherical grid and derivatives.
Definition: SPIgrid.hpp:30
SPI::set_Cheby_mapped_y
SPIVec set_Cheby_mapped_y(PetscScalar a, PetscScalar b, PetscInt ny)
create a stretched Chebyshev grid from [a,b] using default Chebfun like mapping
Definition: SPIgrid.cpp:318
SPI::linspace
SPIVec linspace(const PetscScalar begin, const PetscScalar end, const PetscInt rows)
return linspace of number of rows equally spaced points between begin and end
Definition: SPIVec.cpp:645
SPI::SPIgrid2D::ttype
gridtype ttype
type of grid for time dimension
Definition: SPIgrid.hpp:98
SPI::SPIgrid2D::set_derivatives
void set_derivatives(PetscInt order=4)
sets derivatives Dy and Dyy using saved grid
Definition: SPIgrid.cpp:692
SPI::set_D_UltraS
std::tuple< SPIMat, SPIMat > set_D_UltraS(SPIVec &x, PetscInt d=1)
set a UltraSpherical operator acting with respect to the collocated grid (keeps everything in UltraSp...
Definition: SPIgrid.cpp:338
SPI::SPIgrid1D::That
SPIMat That
Chebyshev operator taking it from physical space to Chebyshev coefficients.
Definition: SPIgrid.hpp:59
SPI::SPIgrid1D::set_grid
void set_grid(SPIVec &y)
saves grid to internal grid
Definition: SPIgrid.cpp:493
SPI::SPIgrid2D::Ihalfn
SPIMat Ihalfn
negative wavenumbers from FT
Definition: SPIgrid.hpp:92
SPI::SPIgrid2D::I
SPIMat I
identity matrix same size as derivative operators of size ny*nt x ny*nt
Definition: SPIgrid.hpp:112
SPI::SPIgrid2D::Dyy
SPIMat Dyy
2nd derivative operator with respect to y of size ny*nt x ny*nt
Definition: SPIgrid.hpp:108
SPI::SPIVec::name
std::string name
Vec name (important for hdf5 i/o)
Definition: SPIVec.hpp:26
SPI::SPIgrid2D::Dt
SPIMat Dt
1st derivative operator with respect to t of size ny*nt x ny*nt
Definition: SPIgrid.hpp:109
SPI::map_D_Chebyshev
SPIMat map_D_Chebyshev(SPIVec &x, PetscInt d=1)
map a Chebyshev collocated operator acting with respect to the stretched collocated grid
Definition: SPIgrid.cpp:278
SPI::eye
SPIMat eye(const PetscInt n)
create, form, and return identity matrix of size n
Definition: SPIMat.cpp:697
SPI::SPIgrid1D::~SPIgrid1D
~SPIgrid1D()
destructor of saved SPIVec and SPIMat
Definition: SPIgrid.cpp:602
SPI::inv
SPIMat inv(const SPIMat &A)
Solve linear system A*Ainv=B using MatMatSolve.
Definition: SPIMat.cpp:603
SPI::orthogonalize
std::vector< SPIVec > orthogonalize(std::vector< SPIVec > &x, SPIgrid1D &grid)
Definition: SPIgrid.cpp:1323
SPI::SPIgrid2D::FT4
SPIMat FT4
Fourier Transform operator from physical space to wavenumber space for time. For 4 variable state vec...
Definition: SPIgrid.hpp:93
SPI::SPIgrid1D::S0invS1inv
SPIMat S0invS1inv
[in] inverse of S0^-1 * S1^-1
Definition: SPIgrid.hpp:56
SPI::SPIgrid2D::grid1Dt
SPIgrid1D grid1Dt
Definition: SPIgrid.hpp:105
SPI::SPIMat::~SPIMat
~SPIMat()
destructor to delete memory
Definition: SPIMat.cpp:508
SPI::tan
SPIMat tan(const SPIMat &A)
take the tan of each element in a matrix
Definition: SPIMat.cpp:2031
SPI::diag
SPIMat diag(const SPIVec &diag, const PetscInt k=0)
set diagonal of matrix
Definition: SPIMat.cpp:728
SPI::dft
SPIMat dft(PetscInt nt)
Definition: SPIgrid.cpp:1573
SPI::set_Fourier_t
SPIVec set_Fourier_t(PetscScalar T, PetscInt ny)
create a Fourier grid from [0,T]
Definition: SPIgrid.cpp:392
SPI::SPIVec1Dto2D
SPIVec SPIVec1Dto2D(SPIgrid2D &grid2D, SPIVec &u)
expand a 1D vector to a 2D vector copying data along time dimension
Definition: SPIgrid.cpp:850
SPI::SPIgrid1D::Ihalf
SPIMat Ihalf
positive wavenumbers from FT
Definition: SPIgrid.hpp:62
SPI::SPIgrid2D::set_operators
void set_operators()
sets zero and identity operators for grid
Definition: SPIgrid.cpp:772
SPI::SPIgrid2D::FTinv4
SPIMat FTinv4
inverse Fourier Transform from wavenumber space to physical space for time. For 4 variable state vect...
Definition: SPIgrid.hpp:94
SPI::trapz
PetscScalar trapz(const SPIVec y)
trapezoidal integration of y with unity coordinate spacing,
Definition: SPIVec.cpp:898
SPI::factorial
PetscInt factorial(PetscInt n)
compute the factorial of n. This is needed for get_D_Coeffs function
Definition: SPIgrid.cpp:6
SPI::SPIgrid1D::ny
PetscInt ny
Definition: SPIgrid.hpp:39
SPI::SPIgrid2D::SPIgrid2D
SPIgrid2D(SPIVec &y, SPIVec &t, std::string name="SPIgrid2D", gridtype y_gridtype=FD, gridtype t_gridtype=Fourier)
constructor with no arguments (set default values)
Definition: SPIgrid.cpp:633
SPI::SPIVec::size
PetscInt size()
Definition: SPIVec.cpp:96
SPI::SPIgrid1D::O
SPIMat O
zero matrix same size as derivative operators
Definition: SPIgrid.hpp:64
SPI::abs
SPIMat abs(const SPIMat &A)
take the abs of each element in a matrix
Definition: SPIMat.cpp:2045
SPI::integrate_coeffs
PetscScalar integrate_coeffs(const SPIVec &a)
integrate a vector of chebyshev Coefficients on a grid
Definition: SPIVec.cpp:813
SPI::SPIgrid2D::Ihalfn4
SPIMat Ihalfn4
negative wavenumbers from FT. For 4 variable state vector
Definition: SPIgrid.hpp:96
SPI::set_Cheby_stretched_y
SPIVec set_Cheby_stretched_y(PetscScalar y_max, PetscInt ny, PetscScalar yi=10.)
create a stretched Chebyshev grid from [0,y_max] with yi being the midpoint location for the number o...
Definition: SPIgrid.cpp:305
SPI::SPIgrid2D::Dy
SPIMat Dy
1st derivative operator with respect to y of size ny*nt x ny*nt
Definition: SPIgrid.hpp:107
SPI::map_D
SPIMat map_D(SPIMat D, SPIVec y, PetscInt d, PetscInt order=4)
map the derivative operator to the proper y grid
Definition: SPIgrid.cpp:52
SPI::SPIgrid1D::P
SPIMat P
row permutation matrix for UltraSpherical operators to shift rows from bottom to top to reduce LU fac...
Definition: SPIgrid.hpp:57
SPI::set_FD_stretched_y
SPIVec set_FD_stretched_y(PetscScalar y_max, PetscInt ny, PetscScalar delta=2.0001)
set stretched grid from [0,y_max] using tanh stretching for use with finite difference operators
Definition: SPIgrid.cpp:214
SPI::SPIgrid2D::avgt
SPIMat avgt
average in time operator
Definition: SPIgrid.hpp:113
SPI::set_Cheby_y
SPIVec set_Cheby_y(PetscInt ny)
create a Chebyshev grid from [-1,1]
Definition: SPIgrid.cpp:329
SPI::SPIgrid2D::That
SPIMat That
transform from physical to chebyshev
Definition: SPIgrid.hpp:86
SPI::SPIgrid2D
Definition: SPIgrid.hpp:72
SPI::SPIgrid1D::ytype
gridtype ytype
type of grid
Definition: SPIgrid.hpp:48
SPI::SPIgrid2D::name
std::string name
name of grid
Definition: SPIgrid.hpp:81
SPI::SPIgrid2D::Ihalf4
SPIMat Ihalf4
positive wavenumbers from FT. For 4 variable state vector
Definition: SPIgrid.hpp:95
SPI::sqrt
SPIVec sqrt(const SPIVec &A)
take the atanh of each element in a vector
Definition: SPIVec.cpp:741
SPI::SPIVec::~SPIVec
~SPIVec()
Definition: SPIVec.cpp:483
SPI::SPIMat::real
SPIMat & real()
take the real part of the vector
Definition: SPIMat.cpp:440
SPI::SPIgrid1D::S1S0That
SPIMat S1S0That
UltraSpherical helper matrix S1*S0*That for baseflow.
Definition: SPIgrid.hpp:55
SPI::SPIgrid1D::flag_set_operators
PetscBool flag_set_operators
flag if set_operators has been executed
Definition: SPIgrid.hpp:69
SPI::pow
SPIVec pow(const SPIVec &A, SPIVec &B)
take the pow of each element in the vectors
Definition: SPIVec.cpp:775
SPI::SPIgrid1D::print
void print()
saves grid to internal grid
Definition: SPIgrid.cpp:468
SPI::gridtype
gridtype
enumeration of grid types
Definition: SPIgrid.hpp:25
SPI::SPIMat
Definition: SPIMat.hpp:17
SPI::SPIgrid1D::S0
SPIMat S0
UltraSpherical helper matrix S_0 takes chebyshev coefficients and outputs C^(1) coefficients.
Definition: SPIgrid.hpp:53
SPI::SPIgrid1D::flag_set_derivatives
PetscBool flag_set_derivatives
flag if set_derivatives has been executed
Definition: SPIgrid.hpp:68
SPI::block
SPIMat block(const Block2D< SPIMat > Blocks)
set block matrices using an input array of size rows*cols. Fills rows first
Definition: SPIMat.cpp:1857
SPI::tanh
SPIVec tanh(const SPIVec &A)
take the tanh of each element in a vector
Definition: SPIVec.cpp:727
SPI::acos
SPIMat acos(const SPIMat &A)
take the arccos of each element in a matrix
Definition: SPIMat.cpp:2029
SPI::SPIgrid2D::t
SPIVec t
grid for time dimension
Definition: SPIgrid.hpp:84
SPI::SPIVec::rows
PetscInt rows
number of rows in vec
Definition: SPIVec.hpp:14
SPI::proj
SPIVec proj(SPIVec &u, SPIVec &v, SPIgrid1D &grid)
Definition: SPIgrid.cpp:1220
SPI::SPIgrid2D::set_grid
void set_grid(SPIVec &y, SPIVec &t)
saves grid to internal grid
Definition: SPIgrid.cpp:676
SPI::SPIgrid1D::FTinv
SPIMat FTinv
inverse Fourier Transform
Definition: SPIgrid.hpp:61
SPI::SPIgrid2D::grid1Dy
SPIgrid1D grid1Dy
Definition: SPIgrid.hpp:105
SPI::SPIgrid1D::set_operators
void set_operators()
sets zero and identity operators for grid
Definition: SPIgrid.cpp:574
SPI::SPIgrid2D::ytype
gridtype ytype
type of grid for wall-normal dimension
Definition: SPIgrid.hpp:97
SPI::set_D_Fourier
SPIMat set_D_Fourier(SPIVec t, PetscInt d=1)
create a Fourier derivative operator acting on grid t
Definition: SPIgrid.cpp:407
SPI::SPIgrid2D::nt
PetscInt nt
Definition: SPIgrid.hpp:75
SPI::SPIgrid2D::flag_set_grid
PetscBool flag_set_grid
flag if set_grid has been executed
Definition: SPIgrid.hpp:115
SPI::SPIgrid1D::name
std::string name
name of grid
Definition: SPIgrid.hpp:45
SPI::meshgrid
std::tuple< SPIMat, SPIMat > meshgrid(SPIVec &x, SPIVec &y)
Definition: SPIMat.cpp:1894
SPI::SPIgrid1D::S1
SPIMat S1
UltraSpherical helper matrices S_1 takes C^(1) coefficients and outputs C^(2) coefficients.
Definition: SPIgrid.hpp:54
SPI::SPIgrid2D::S1S0That
SPIMat S1S0That
UltraSpherical helper matrix S1*S0*That for baseflow.
Definition: SPIgrid.hpp:87
SPI::SPIgrid2D::T
SPIMat T
transform from chebyshev to physical
Definition: SPIgrid.hpp:85
SPI::cos
SPIMat cos(const SPIMat &A)
take the cos of each element in a matrix
Definition: SPIMat.cpp:2027
SPI::SPIMat::zero_row
SPIMat & zero_row(const PetscInt row)
set a row to zero
Definition: SPIMat.cpp:451
SPI::SPIgrid1D::I
SPIMat I
identity matrix same size as derivative operators
Definition: SPIgrid.hpp:65
SPI::SPIMat::T
PetscInt T(SPIMat &A)
A = Transpose(*this.mat) operation with initialization of A.
Definition: SPIMat.cpp:392
SPIgrid.hpp
SPI::zeros
SPIMat zeros(const PetscInt m, const PetscInt n)
create, form, and return zeros matrix of size mxn
Definition: SPIMat.cpp:708
SPI::interp1D_Mat
SPIMat interp1D_Mat(SPIgrid1D &grid1, SPIgrid1D &grid2)
Definition: SPIgrid.cpp:1506
SPI::SPIMat::cols
PetscInt cols
number of columns in mat
Definition: SPIMat.hpp:19
SPI::set_D_Chebyshev
SPIMat set_D_Chebyshev(SPIVec &x, PetscInt d=1, PetscBool need_map=PETSC_FALSE)
set a Chebyshev collocated operator acting with respect to the collocated grid
Definition: SPIgrid.cpp:225
SPI::set_T_That
std::tuple< SPIMat, SPIMat > set_T_That(PetscInt n)
set a T and That operators acting with respect to the Chebyshev collocated grid, That: physical -> Ch...
Definition: SPIgrid.cpp:367
Generated by   doxygen 1.8.17