SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
tests.cpp
Go to the documentation of this file.
1 #include "SPImain.hpp"
2 #include "tests.hpp"
3 #include <time.h>
4 #include <slepcsvd.h>
5 
6 void test_if_true(PetscBool test,std::string name){
7  if (test) { SPI::printf("\x1b[32m"+name+" test passed"+"\x1b[0m"); }
8  else{ std::cout<<"\x1b[31m"+name+" test failed"+"\x1b[0m"<<std::endl;}
9 }
10 
11 void test_if_close(PetscScalar value,PetscScalar golden, std::string name, PetscReal tol){
12  PetscReal valuer=PetscRealPart(value);
13  PetscReal goldenr=PetscRealPart(golden);
14  if ((goldenr-tol<=valuer) && (valuer<=goldenr+tol)) { SPI::printf("\x1b[32m"+name+" test passed"+"\x1b[0m"); }
15  else{ //std::cout<<"\x1b[31m"+name+" test failed"+"\x1b[0m"<<std::endl;
16  //std::cout<<" valuer="<<valuer<<std::endl;
17  //std::cout<<" goldenr="<<goldenr<<std::endl;
18  SPI::printf("\x1b[31m"+name+" test failed"+"\x1b[0m");
19  SPI::printf("\x1b[31m goldenr=%.20f \x1b[0m",goldenr);
20  SPI::printf("\x1b[31m valuer =%.20f \x1b[0m",valuer );
21  }
22 }
23 
24 int tests(){
25  PetscInt m=4, n=4;
26  PetscBool alltests=PETSC_FALSE;
27  // Vec tests
28  if(alltests){
29  SPI::printf("------------ Vec tests start-------------");
30  // initialize SPIVec and Init function
31  SPI::SPIVec X1(4,"X1"),X2(4,"X2"),X3("X3");
32 
33  // () operators
34  X1(0,1.+1.*PETSC_i);// PetscScalar
35  X1(1,1.); // const double
36  X1(2,1); // const int
37  X1.set(3,1.+1.*PETSC_i); // set function
38  X1(); // assemble
39  test_if_close(X1(0,PETSC_TRUE),1.,"SPIVec() 0");
40  test_if_close(X1(1,PETSC_TRUE),1.,"SPIVec() 1");
41  test_if_close(X1(2,PETSC_TRUE),1.,"SPIVec() 2");
42  test_if_close(X1(3,PETSC_TRUE),1.,"SPIVec() 3");
43 
44  SPI::SPIVec X5(X1,"X5_copy_X1"); // initialize with SPIVec
45  test_if_close(X5(2,PETSC_TRUE),1.,"SPIVec(SPIVec)");
46 
47  // equals operator
48  X2=X1;
49  test_if_close(X2(2,PETSC_TRUE),1.,"SPIVec=SPIVec");
50 
51  // +- operators
52  X2+=X1;
53  test_if_close(X2(2,PETSC_TRUE),2.,"SPIVec+=SPIVec");
54  X3 = X2+X1;
55  test_if_close(X3(2,PETSC_TRUE),3.,"SPIVec+SPIVec");
56  X3.axpy(1.,X1);
57  test_if_close(X3(2,PETSC_TRUE),4.,"SPIVec.axpy(PetscScalar,SPIVec)");
58  X3-=X1;
59  test_if_close(X3(2,PETSC_TRUE),3.,"SPIVec-=SPIVec");
60  SPI::SPIVec X4("X4=0");
61  X4 = X2-X1;
62  X4 -= X1;
63  test_if_close(X4(2,PETSC_TRUE),0.,"SPIVec-SPIVec");
64 
65  // * operators
66  SPI::SPIVec X6("X6=i*X1"), X7("X7"), X8("X8");
67  X6 = (X1*(0.+1.*PETSC_i));// PetscScalar
68  test_if_close(X6(3,PETSC_TRUE),-1.,"SPIVec*PetscScalar");
69  X7 = X1*7.; // double
70  test_if_close(X7(3,PETSC_TRUE),7.,"SPIVec*double");
71  X8 = X1*8; // int
72  test_if_close(X8(3,PETSC_TRUE),8.,"SPIVec*int");
73  X8*=0.+1.*PETSC_i;
74  test_if_close(X8(3,PETSC_TRUE),-8.,"SPIVec*=PetscScalar");
75  SPI::SPIVec X9("X9=X1*X2");
76  X9 = X1*X2;
77  test_if_close(X9(2,PETSC_TRUE),2.,"SPIVec*SPIVec");
78  }
79  if(alltests){
80  // Mat tests
81  SPI::printf("------------ Mat tests start ---------------");
82 
83  // test constructors
84  SPI::SPIMat D("D");
85  SPI::SPIMat I(SPI::eye(m),"I");
86  SPI::SPIMat B(m,"B");
87  SPI::SPIMat A2(m,n,"A2"), E(4*m,4*n,"E");
88 
89  // set operators
90  for (int i=0; i<m; i++){
91  A2(i,i,1.+PETSC_i);
92  }
93  A2(0,n-1,0.43);
94  A2();
95  test_if_close(PetscImaginaryPart(A2(2,2,PETSC_TRUE)),1.,"SPIMat(PetscInt,PetscInt,PETSC_TRUE)");
96  B(0,1,1);
97  B();
98  test_if_close(B(0,1,PETSC_TRUE),1.,"SPIMat(PetscInt,PetscInt,PetscInt)");
99  B=(3.4+PETSC_i*4.2)*A2+4.*A2;
100  test_if_close(B(1,1,PETSC_TRUE),3.2,"PetscScalar*SPIMat+double*SPIMat");
101  E(0,0,A2);
102  E(m,n,B);
103  E();
104  test_if_close(E(1,1,PETSC_TRUE),1., "SPIMat(PetscInt,PetscInt,SPIMat) 1");
105  test_if_close(E(4,7,PETSC_TRUE),3.182, "SPIMat(PetscInt,PetscInt,SPIMat) 2");
106  A2+=B;
107  test_if_close(A2(1,1,PETSC_TRUE),4.2, "SPIMat+=SPIMat");
108  SPI::SPIMat A2T("A2T");
109  A2.T(A2T);
110  test_if_close(A2T(3,0,PETSC_TRUE),3.612, "SPIMat.T(SPIMat)");
111  D = A2T*B;
112  test_if_close(D(0,3,PETSC_TRUE),-9.3912, "SPIMat*SPIMat");
113  D *= 4.;
114  test_if_close(D(0,3,PETSC_TRUE),4.*-9.3912, "SPIMat*=PetscScalar");
115  D.~SPIMat();
116  D = A2*I;
117  test_if_close(A2(0,3,PETSC_TRUE),3.612, "SPIMat*eye(SPIMat)");
118  B.T();
119  test_if_close(B(3,0,PETSC_TRUE),3.182, "SPIMat.T()");
120  D=I;
121  test_if_close(D(1,1,PETSC_TRUE),1., "SPIMat=SPIMat");
122  E(4,8,B);// insert
123  E();
124  test_if_close(E(7,8,PETSC_TRUE),3.182, "SPIMat(PetscInt,PetscInt,SPIMat) 3");
125  test_if_close(PetscImaginaryPart(E(6,10,PETSC_TRUE)),11.6, "SPIMat(PetscInt,PetscInt,SPIMat) 4");
126  E.H();
127  test_if_close(PetscImaginaryPart(E(10,6,PETSC_TRUE)),-11.6, "SPIMat.H()");
128  E.conj();
129  test_if_close(PetscImaginaryPart(E(10,6,PETSC_TRUE)),11.6, "SPIMat.conj()");
130  E.T();
131  test_if_close(PetscImaginaryPart(E(6,10,PETSC_TRUE)),11.6, "SPIMat.T()");
132  SPI::SPIMat F(E,"F");
133  test_if_close(PetscImaginaryPart(F(6,10,PETSC_TRUE)),11.6, "SPIMat(SPIMat)");
134  D.diag();
135  test_if_close(D.diag()(1,PETSC_TRUE),1., "diag(SPIMat)");
136  SPI::printf("------------ Mat tests end ---------------");
137  }
138  if(alltests){
139  SPI::printf("------------ A*x tests start ---------------");
140  SPI::SPIMat A(4,4,"A");
141  SPI::SPIVec x(4,"x"),b;
142 
143  for(int i=0; i<4; i++){
144  A(i,i,1.+i);
145  x(i,2.3+PETSC_i*1.);
146  }
147  A();
148  x();
149  b=A*x;
150  test_if_close(b(2,PETSC_TRUE),6.9,"SPIMat*SPIVec");
151  //A.print();
152  //x.print();
153  //b.print();
154  //Vec b2;
155  //VecCreate(PETSC_COMM_WORLD,&b2);
156  //VecSetSizes(b2,PETSC_DECIDE,4);
157  //VecSetType(b2,VECMPICUDA);
158  //MatMult(A.mat,x.vec,b2);
159  //VecView(b2,PETSC_VIEWER_STDOUT_WORLD);
160 
161  SPI::printf("------------ A*x tests end ---------------");
162  }
163 
164  if(alltests){
165  SPI::printf("------------ A*x=b tests start ---------------");
166  SPI::SPIMat A(4,4,"A");
167  SPI::SPIVec b(4,"x"),x;
168 
169  for(int i=0; i<4; i++){
170  A(i,i,1.+i);
171  b(i,2.);
172  }
173  A();
174  b();
175  x=b/A;
176  test_if_close(x(3,PETSC_TRUE),0.5,"SPIVec/SPIMat");
177  SPI::printf("------------ A*x=b tests end ---------------");
178  }
179  if(alltests){
180  SPI::printf("------------ Mat func tests start-------------");
181  SPI::SPIMat I(SPI::eye(4),"I-identity");
182  test_if_close(I(1,1,PETSC_TRUE),1.,"eye(PetscInt)");
183 
184  SPI::SPIVec two(3,"two");
185  for (int i=0; i<3; i++) two(i,2.);
186  two();
187 
188  SPI::SPIMat A("A");
189  A = SPI::diag(two);
190  test_if_close(A(2,2,PETSC_TRUE),2.,"diag(SPIVec)");
191 
192  SPI::SPIMat C("C");
193  C=SPI::kron(I,A);
194  test_if_close(C(6,6,PETSC_TRUE),2.,"kron(SPIMat,SPIMat) 1");
195 
196  SPI::SPIMat D(2,2,"D");
197  D(0,0,1.); D(0,1,2.);
198  D(1,0,3.); D(1,1,4.);
199  D();
200 
201  test_if_close(SPI::kron(D,I)(5,1,PETSC_TRUE),3.,"kron(SPIMat,SPIMat) 2");
202  SPI::printf("------------ Mat func tests end -------------");
203  }
204  if(alltests){
205  SPI::printf("------------ Mat eig tests start-------------");
206  SPI::SPIMat A(2,"A");
207  A(0,0,1.);
208  A(0,1,1.);
209  A(1,1,1.);
210  A();
211  SPI::SPIMat B(SPI::eye(2),"I-identity");
212  PetscScalar alpha;
213  SPI::SPIVec eig(2,"eig1");
214  SPI::SPIVec eig2(2,"eig2");
215  std::tie(alpha,eig,eig2) = SPI::eig(A,B,1.0+PETSC_i*0.5);
216  eig /= eig.max(); // normalize by max amplitude
217  test_if_close(alpha,1.,"eig(SPIMat,SPIMat,PetscScalar) 1",1.E-8);
218  PetscScalar alpha2;
219  std::tie(alpha2,eig2,eig) = SPI::eig(A,B,-1.0+PETSC_i*0.00005,1.E-19,10);
220  eig2 /= eig2.max(); // normalize by max amplitude
221  test_if_close(alpha2,1.,"eig(SPIMat,SPIMat,PetscScalar) 2",1.E-7);
222  // eig_right
223  std::tie(alpha,eig) = SPI::eig_right(A,B,1.0+PETSC_i*0.5);
224  test_if_close(alpha,1.,"eig_right(SPIMat,SPIMat,PetscScalar) 1",1.E-8);
225  std::tie(alpha2,eig) = SPI::eig_right(A,B,-1.0+PETSC_i*0.00005,1.E-19,10);
226  test_if_close(alpha2,1.,"eig_right(SPIMat,SPIMat,PetscScalar) 2",1.E-7);
227  // eig_init
228  std::tie(alpha,eig,eig2) = SPI::eig_init(A,B,1.0+PETSC_i*0.5,eig.conj(),eig);
229  test_if_close(alpha,1.,"eig_init(SPIMat,SPIMat,PetscScalar,SPIVec) 1",1.E-8);
230  std::tie(alpha2,eig2,eig) = SPI::eig_init(A,B,-1.0+PETSC_i*0.00005,eig.conj(),eig,1.E-19,10);
231  test_if_close(alpha2,1.,"eig_init(SPIMat,SPIMat,PetscScalar,SPIVec,PetscReal,PetscInt) 2",1.E-7);
232  // eig_init_right
233  std::tie(alpha,eig) = SPI::eig_init_right(A,B,1.0+PETSC_i*0.5,eig);
234  test_if_close(alpha,1.,"eig_init_right(SPIMat,SPIMat,PetscScalar,SPIVec) 1",1.E-8);
235  std::tie(alpha2,eig) = SPI::eig_init_right(A,B,-1.0+PETSC_i*0.00005,eig,1.E-19,10);
236  test_if_close(alpha2,1.,"eig_init_right(SPIMat,SPIMat,PetscScalar,SPIVec,PetscReal,PetscInt) 2",1.E-7);
237 
238  // TODO these cases don't work as it is a linear problem
239  //std::tie(alpha,eig) = SPI::polyeig({A,-B},1.0+0.5*PETSC_i);
240  //test_if_close(alpha,1.,"polyeig(std::vector<SPIMat>,PetscScalar) 1",1.E-8);
241  //std::tie(alpha,eig) = SPI::polyeig({A},1.0+0.5*PETSC_i);
242  //test_if_close(alpha,1.,"polyeig(std::vector<SPIMat>,PetscScalar) 2",1.E-8);
243 
244  SPI::SPIMat M(4,"M"),C(4,"C"),K(4,"K"); // (M*lambda^2 + C*lambda + K)*x = 0 from MatLab mathworks documentation
245  // M
246  M(0,0,3.);
247  M(1,1,1.);
248  M(2,2,3.);
249  M(3,3,1.);
250  // C
251  C(0,0,0.4); C(0,2,-0.3);
252 
253  C(2,0,-0.3); C(2,2,0.5); C(2,3,-0.2);
254  C(3,2,-0.2);C(3,3,0.2);
255  // K
256  K(0,0,-7.); K(0,1,2.); K(0,2,4);
257  K(1,0,2); K(1,1,-4); K(1,2,2);
258  K(2,0,4); K(2,1,2); K(2,2,-9); K(2,3,3);
259  K(3,2,3); K(3,3,-3);
260  // assemble
261  M();C();K();
262 
263 
264  PetscScalar alpha4;
265  SPI::SPIVec eig4(4),eig5(4);
266  std::tie(alpha4,eig4) = SPI::polyeig({K,C,M},-2.5+0.5*PETSC_i);
267  test_if_close(alpha4,-2.44985,"polyeig(std::vector<SPIMat>,PetscScalar) 1",1.E-5);
268  std::tie(alpha4,eig4) = SPI::polyeig({K,C,M},0.33+0.005*PETSC_i);
269  test_if_close(alpha4,0.3353,"polyeig(std::vector<SPIMat>,PetscScalar) 2",1.E-5);
270  // polyeig_init
271  //std::tie(alpha,eig) = SPI::polyeig_init({A,-B},1.0+0.5*PETSC_i,eig);
272  //test_if_close(alpha,1.,"polyeig_init(std::vector<SPIMat>,PetscScalar,SPIVec) 1",1.E-8);
273  //std::tie(alpha,eig) = SPI::polyeig_init({A},1.0+0.5*PETSC_i,eig);
274  //test_if_close(alpha,1.,"polyeig_init(std::vector<SPIMat>,PetscScalar,SPIVec) 2",1.E-8);
275  std::tie(alpha4,eig4) = SPI::polyeig_init({K,C,M},-2.5+0.5*PETSC_i,eig4);
276  test_if_close(alpha4,-2.44985,"polyeig_init(std::vector<SPIMat>,PetscScalar) 1",1.E-5);
277  std::tie(alpha4,eig4) = SPI::polyeig_init({K,C,M},0.33+0.005*PETSC_i,eig4);
278  test_if_close(alpha4,0.3353,"polyeig_init(std::vector<SPIMat>,PetscScalar) 2",1.E-5);
279 
280  SPI::printf("------------ Mat eig tests end -------------");
281  }
282 
283  if(alltests){// I/O using hdf5
284  SPI::printf("------------ I/O tests start -------------");
285  SPI::SPIVec A(2,"A_Vec");
286  A(0,1.);
287  A(1,2.);
288  A();
289  SPI::save(A,"saved_data.hdf5");
290  SPI::SPIVec B(2,"B_Vec");
291  B(0,1.+PETSC_i*0.5);
292  B(1,2.*PETSC_i);
293  B();
294  SPI::save(B,"saved_data.hdf5");
295  SPI::printf("------------ I/O tests end -------------");
296  }
297  if(alltests){
298  SPI::printf("------------ I/O tests2 start -------------");
299  SPI::SPIVec A_read(2,"A_Vec");
300  SPI::load(A_read,"saved_data.hdf5");
301  test_if_close(A_read(0,PETSC_TRUE),1.,"load(SPIVec,std::string) 1");
302  test_if_close(A_read(1,PETSC_TRUE),2.,"load(SPIVec,std::string) 2");
303  SPI::SPIVec B_read(2,"B_Vec");
304  SPI::load(B_read,"saved_data.hdf5");
305  test_if_close(PetscImaginaryPart(B_read(0,PETSC_TRUE)),0.5,"load(SPIVec,std::string) 3");
306  test_if_close(B_read(1,PETSC_TRUE),0.,"load(SPIVec,std::string) 4");
307  SPI::printf("------------ I/O tests2 end -------------");
308  }
309  if(alltests){// I/O using binary for Mat
310  SPI::printf("------------ I/O tests3 start -------------");
311  SPI::SPIMat A(2,2,"A");
312  A(0,0,1.);
313  A(1,0,2.);
314  A(1,1,3.+4.59*PETSC_i);
315  A();
316  SPI::save(A,"saved_data_mat.dat");
317  SPI::SPIMat B(2,2,"B");
318  B(0,0,3.);
319  B(1,0,4.);
320  B(1,1,5.+4.89*PETSC_i);
321  B();
322  SPI::save(B,"saved_data_mat.dat");
323  SPI::printf("------------ I/O tests3 end -------------");
324  }
325  if(alltests){
326  SPI::printf("------------ I/O tests4 start -------------");
327  SPI::SPIMat A_read(2,2,"A_Mat");
328  SPI::load(A_read,"saved_data_mat.dat");
329  test_if_close(A_read(0,0,PETSC_TRUE),1.,"load(SPIMat,std::string) 1");
330  test_if_close(A_read(1,0,PETSC_TRUE),2.,"load(SPIMat,std::string) 2");
331  test_if_close(A_read(1,1,PETSC_TRUE),3.,"load(SPIMat,std::string) 3");
332  test_if_close(PetscImaginaryPart(A_read(1,1,PETSC_TRUE)),4.59,"load(SPIMat,std::string) 4");
333  std::vector<SPI::SPIMat> AB(2,SPI::eye(2)*0.);
334  SPI::load(AB,"saved_data_mat.dat");
335  test_if_close(AB[0](0,0,PETSC_TRUE),1.,"load(std::vector<SPIMat>,std::string) 1");
336  test_if_close(AB[0](1,0,PETSC_TRUE),2.,"load(std::vector<SPIMat>,std::string) 2");
337  test_if_close(AB[0](1,1,PETSC_TRUE),3.,"load(std::vector<SPIMat>,std::string) 3");
338  test_if_close(PetscImaginaryPart(AB[0](1,1,PETSC_TRUE)),4.59,"load(std::vector<SPIMat>,std::string) 4");
339  test_if_close(AB[1](0,0,PETSC_TRUE),3.,"load(std::vector<SPIMat>,std::string) 5");
340  test_if_close(AB[1](1,0,PETSC_TRUE),4.,"load(std::vector<SPIMat>,std::string) 6");
341  test_if_close(AB[1](1,1,PETSC_TRUE),5.,"load(std::vector<SPIMat>,std::string) 7");
342  test_if_close(PetscImaginaryPart(AB[1](1,1,PETSC_TRUE)),4.89,"load(std::vector<SPIMat>,std::string) 8");
343 
344  //SPI::draw(AB[1]); // test failed TODO
345 
346  SPI::printf("------------ I/O tests4 end -------------");
347  }
348  if(alltests){
349  SPI::printf("------------ block test start -------------");
350  SPI::SPIMat A(2,"A");
351  A(0,0,2.);
352  A(0,1,3.);
353  A(1,1,4.);
354  A();
355  SPI::SPIMat B(SPI::eye(2),"I-identity");
356 
357  SPI::SPIMat block(SPI::block({{A,B},{A,B}}));
358  block();
359  test_if_close(block(0,0,PETSC_TRUE),2.,"block(std::vector<std::vector<SPIMat>>) 1");
360  test_if_close(block(0,2,PETSC_TRUE),1.,"block(std::vector<std::vector<SPIMat>>) 2");
361  test_if_close(block(3,1,PETSC_TRUE),4.,"block(std::vector<std::vector<SPIMat>>) 3");
362  test_if_close(block(3,3,PETSC_TRUE),1.,"block(std::vector<std::vector<SPIMat>>) 4");
363  SPI::printf("------------ block test end -------------");
364  }
365  if(alltests){
366  SPI::printf("------------ LST_temporal test start -------------");
367  // create grid and derivatives using chebyshev polynomials
368  PetscInt n=128;
369  SPI::SPIVec y(SPI::set_Cheby_y(n),"yCheby");
370  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
371 
372  // channel flow Plane Poiseuille solution
373  SPI::SPIMat U(SPI::diag(1.0-((grid.y)^2)),"U");
374  SPI::SPIMat Uy(SPI::diag(-2.*grid.y),"Uy");
375  SPI::SPIVec o(U.diag()*0.0,"o"); // zero vector
376  SPI::SPIbaseflow channel(U.diag(),o,o,o,Uy.diag(),o,o,o,o,o,o,o,o,o);
377 
378  // parameters for Orr-Sommerfield eq.
379  PetscScalar Re=2000.0;
380  PetscScalar alpha=1.0;
381  PetscScalar beta=0.0;
382 
383  // set parameters into param struct
384  SPI::SPIparams params("channel parameter");
385  params.Re = Re;
386  params.omega = 0.3121002078-0.0197986590*PETSC_i;
387  params.alpha = alpha;
388  params.beta = beta;
389 
390  // solve eigenvalue problem
391  SPI::SPIVec eigenfunction(grid.y.rows*4,"q");
392  PetscScalar eigenvalue;
393  eigenfunction.name = "eigenfunction";
394  std::tie(eigenvalue,eigenfunction) = SPI::LST_temporal(params,grid,channel);
395  //SPI::printfc("eigenvalue is: %.10f+%.10fi",eigenvalue);
396  //eigenfunction.print();
397  test_if_close(eigenvalue,0.3121002979-0.0197986590*PETSC_i,"LST_temporal 1",1e-9);
398  std::tie(eigenvalue,eigenfunction) = SPI::LST_temporal(params,grid,channel,eigenfunction);
399  //SPI::printfc(" eigenvalue is: %.10f+%.10fi",params.omega);
400  test_if_close(eigenvalue,0.3121002979-0.0197986590*PETSC_i,"LST_temporal 2",1e-9);
401  SPI::printf("------------ LST_temporal test end -------------");
402  }
403  if(alltests){
404  SPI::printf("------------ LST_spatial test start -------------");
405  PetscInt n=128;
406  SPI::SPIVec y(SPI::set_Cheby_y(n),"yCheby");
407  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
408  // channel flow Orr-Sommerfeld solution
409  SPI::SPIMat U(SPI::diag(1.0-((grid.y)^2)),"U");
410  SPI::SPIMat Uy(SPI::diag(-2.*grid.y),"Uy");
411  //SPI::SPIMat Uyy(SPI::diag(-2.*SPI::ones(grid.y.rows)),"Uyy");
412 
413  SPI::SPIparams params("channel parameter");
414  params.Re = 2000.0;
415  params.omega = 0.3;
416  params.alpha = 0.97875+0.044394*PETSC_i;
417  params.beta = 0.0;
418 
419  SPI::SPIVec eigenfunction(grid.y.rows*6,"q");
420  PetscScalar eigenvalue;
421  eigenfunction.name = "eigenfunction";
422 
423  SPI::SPIVec o(U.diag()*0.0,"o");
424  SPI::SPIbaseflow channel(U.diag(),o,o,o,Uy.diag(),o,o,o,o,o,o,o,o,o);
425 
426  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,channel);
427  test_if_close(eigenvalue,0.97875+0.044394*PETSC_i,"LST_spatial 1",1e-5);
428  test_if_close(params.alpha,0.97875+0.044394*PETSC_i,"LST_spatial 1",1e-5);
429  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,channel,eigenfunction); // with initial guess
430  test_if_close(eigenvalue,0.97875+0.044394*PETSC_i,"LST_spatial 1",1e-5);
431  test_if_close(params.alpha,0.97875+0.044394*PETSC_i,"LST_spatial 1",1e-5);
432  params.alpha = 0.34312+0.049677*PETSC_i;
433  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,channel,eigenfunction);
434  test_if_close(eigenvalue,0.34312+0.049677*PETSC_i,"LST_spatial 2",1e-5);
435  test_if_close(params.alpha,0.34312+0.049677*PETSC_i,"LST_spatial 2",1e-5);
436  params.alpha = 0.61+0.1*PETSC_i;
437  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,channel,eigenfunction);
438  test_if_close(eigenvalue,0.61167+0.140492*PETSC_i,"LST_spatial 3",1e-5);
439  test_if_close(params.alpha,0.61167+0.140492*PETSC_i,"LST_spatial 3",1e-5);
440  SPI::printf("------------ LST_spatial test end -------------");
441  }
442  if(alltests){
443  SPI::printf("------------ LSTNP_spatial test start -------------");
444  PetscInt n=64;
445  SPI::SPIVec y(SPI::set_Cheby_y(n),"yCheby");
446  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
447  // channel flow Orr-Sommerfeld solution
448  SPI::SPIVec U((1.0-(grid.y*grid.y)),"U");
449  SPI::SPIVec Uy((-2.*grid.y),"Uy");
450  //SPI::SPIMat Uyy(SPI::diag(-2.*SPI::ones(grid.y.rows)),"Uyy");
451 
452  SPI::SPIparams params("channel parameter");
453  params.Re = 2000.0;
454  params.omega = 0.3;
455  params.alpha = 0.97875+0.044394*PETSC_i;
456  params.beta = 0.0;
457 
458  SPI::SPIVec eigenfunction(grid.y.rows*16,"q");
459  SPI::SPIVec eigenfunction8(grid.y.rows*8,"q");
460  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
461  PetscScalar eigenvalue;
462  eigenfunction.name = "eigenfunction";
463 
464  SPI::SPIVec o(U*0.0,"o");
465  SPI::SPIbaseflow channel(U,o,o,o,Uy,o,o,o,o,o,o,o,o,o);
466 
467  PetscScalar cg;
468  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,channel);
469  SPI::printfc("cg = %g+%gi",cg);
470  test_if_close(eigenvalue,0.978748+0.04439397*PETSC_i,"LSTNP_spatial 1",1e-5);
471  test_if_close(params.alpha,0.978748+0.04439397*PETSC_i,"LSTNP_spatial 1",1e-5);
472  // LSTNP_spatial_right
473  std::tie(eigenvalue,eigenfunction8) = SPI::LSTNP_spatial_right(params,grid,channel);
474  test_if_close(params.alpha,0.978748+0.04439397*PETSC_i,"LSTNP_spatial_right 1",1e-5);
475  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
476  params.alpha = 0.978748+0.04439397*PETSC_i;
477  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,channel,eigenfunction.conj(),eigenfunction); // with initial guess
478  SPI::printfc("cg = %g+%gi",cg);
479  test_if_close(eigenvalue,0.978748+0.04439397*PETSC_i,"LSTNP_spatial 2",1e-5);
480  test_if_close(params.alpha,0.978748+0.04439397*PETSC_i,"LSTNP_spatial 2",1e-5);
481  std::tie(eigenvalue,eigenfunction8) = SPI::LSTNP_spatial_right(params,grid,channel,eigenfunction);
482  test_if_close(params.alpha,0.978748+0.04439397*PETSC_i,"LSTNP_spatial_right 2",1e-5);
483  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
484  params.alpha = 0.34305+0.0498376872*PETSC_i;
485  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,channel,leigenfunction,eigenfunction);
486  SPI::printfc("cg = %g+%gi",cg);
487  test_if_close(eigenvalue,0.34305+0.049837*PETSC_i,"LSTNP_spatial 3",1e-5);
488  std::tie(eigenvalue,eigenfunction8) = SPI::LSTNP_spatial_right(params,grid,channel,eigenfunction);
489  test_if_close(eigenvalue,0.34305+0.049837*PETSC_i,"LSTNP_spatial_right 3",1e-5);
490  //test_if_close(params.alpha,0.34305+0.049837*PETSC_i,"LSTNP_spatial 2",1e-5);
491  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
492  params.alpha = 0.6116672+0.140493*PETSC_i;
493  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,channel,leigenfunction,eigenfunction);
494  SPI::printfc("cg = %g+%gi",cg);
495  test_if_close(eigenvalue,0.6116672+0.140493*PETSC_i,"LSTNP_spatial 4",1e-5);
496  std::tie(eigenvalue,eigenfunction8) = SPI::LSTNP_spatial_right(params,grid,channel,eigenfunction);
497  test_if_close(eigenvalue,0.6116672+0.140493*PETSC_i,"LSTNP_spatial_right 4",1e-5);
498  //test_if_close(params.alpha,0.635797+0.08405*PETSC_i,"LSTNP_spatial 3",1e-5);
499  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
500  SPI::printf("------------ LSTNP_spatial test end -------------");
501  }
502  if(0){
503  SPI::printf("------------ Chebyshev derivatives start -----------");
504  PetscInt n=16;
505  //SPI::SPIVec y(SPI::set_Cheby_y(n),"yCheby");
506  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0,61.,n) ,"yCheby");
507  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
508  grid.Dyy.print();
509  SPI::printf("------------ Chebyshev derivatives end -----------");
510  }
511  if(alltests){
512  SPI::printf("------------ LSTNP Blasius boundary layer start -----------");
513  PetscInt n=168;
514  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,61.,n) ,"yCheby");
515  //SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
516  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
517  //SPI::printf("n = %d",n);
518  SPI::SPIVec y(SPI::set_FD_stretched_y(61.,n) ,"yFD");
519  SPI::SPIgrid1D grid(y,"grid",SPI::FD);
520  //grid.print();
521  //(grid.Dy*grid.y).print();
522 
523  SPI::SPIparams params("Blasius parameters");
524  params.Re = 400.0;
525  params.omega = 86.*params.Re/(1000000.);
526  params.nu = 1./params.Re;
527  params.x = params.Re;
528  params.x_start = params.x;
529  params.alpha = 0.094966+0.004564*PETSC_i;
530  //params.alpha = 0.106654+0.0018979*PETSC_i;
531  params.beta = 0.0;
532  //params.print();
533 
534  SPI::SPIVec eigenfunction(grid.y.rows*16,"q");
535  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
536  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
537  PetscScalar eigenvalue;
538  eigenfunction.name = "eigenfunction";
539 
540  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
541  if(0){ // set to parallel baseflow
542  SPI::SPIVec o(SPI::zeros(n),"zero");
543  bl_flow.Ux = o;
544  bl_flow.Uxy = o;
545  bl_flow.V = o;
546  bl_flow.Vy = o;
547  }
548  //bl_flow.print();
549 
550  PetscScalar cg;
551  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
552  SPI::printfc("cg = %g+%gi",cg);
553  //std::tie(eigenvalue,eig_vec) = SPI::LST_spatial(params,grid,bl_flow);
554  test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial 1",1e-5);
555  test_if_close(params.alpha,0.094966+0.004564*PETSC_i,"LSTNP_spatial 1",1e-5);
556  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
557  test_if_close(params.alpha,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 1",1e-5);
558  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
559 
560  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
561  SPI::printfc("cg = %g+%gi",cg);
562  test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial 2",1e-5);
563  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
564  test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 2",1e-5);
565  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
566 
567  params.alpha = 0.1067+0.001898*PETSC_i;
568  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
569  SPI::printfc("cg = %g+%gi",cg);
570  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 3",1e-5);
571  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
572  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 3",1e-5);
573  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
574  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
575  SPI::printfc("cg = %g+%gi",cg);
576  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 4",1e-5);
577  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
578  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
579  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 4",1e-5);
580 
581  SPI::printf("------------ LSTNP Blasius boundary layer end -----------");
582  }
583  if(0){
584  SPI::printf("------------ LSTNP Blasius boundary layer start -----------");
585  PetscInt n=168;
586  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,61.,n) ,"yCheby");
587  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
588  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
589  //SPI::printf("n = %d",n);
590  //SPI::SPIVec y(SPI::set_FD_stretched_y(61.,n) ,"yFD");
591  //SPI::SPIgrid1D grid(y,"grid",SPI::FD);
592  //grid.print();
593  //(grid.Dy*grid.y).print();
594 
595  SPI::SPIparams params("Blasius parameters");
596  params.Re = 400.0;
597  params.omega = 86.*params.Re/(1000000.);
598  params.nu = 1./params.Re;
599  params.x = params.Re;
600  params.x_start = params.x;
601  params.alpha = 0.094966+0.004564*PETSC_i;
602  //params.alpha = 0.106654+0.0018979*PETSC_i;
603  params.beta = 0.0;
604  //params.print();
605 
606  SPI::SPIVec eigenfunction(grid.y.rows*16,"q");
607  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
608  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
609  PetscScalar eigenvalue;
610  eigenfunction.name = "eigenfunction";
611 
612  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
613  if(0){ // set to parallel baseflow
614  SPI::SPIVec o(SPI::zeros(n),"zero");
615  bl_flow.Ux = o;
616  bl_flow.Uxy = o;
617  bl_flow.V = o;
618  bl_flow.Vy = o;
619  }
620  //bl_flow.print();
621 
622  PetscScalar cg;
623  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
624  SPI::printfc("cg = %g+%gi",cg);
625  //std::tie(eigenvalue,eig_vec) = SPI::LST_spatial(params,grid,bl_flow);
626  test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial 1",1e-5);
627  test_if_close(params.alpha,0.094966+0.004564*PETSC_i,"LSTNP_spatial 1",1e-5);
628  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
629  //test_if_close(params.alpha,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 1",1e-5);
630  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
631 
632  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
633  SPI::printfc("cg = %g+%gi",cg);
634  test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial 2",1e-5);
635  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
636  //test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 2",1e-5);
637  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
638 
639  params.alpha = 0.1067+0.001898*PETSC_i;
640  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
641  SPI::printfc("cg = %g+%gi",cg);
642  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 3",1e-5);
643  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
644  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 3",1e-5);
645  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
646  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
647  SPI::printfc("cg = %g+%gi",cg);
648  test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 4",1e-5);
649  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
650  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
651  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 4",1e-5);
652 
653  SPI::printf("------------ LSTNP Blasius boundary layer end -----------");
654  }
655  if(alltests){
656  SPI::printf("------------ LST Blasius boundary layer start -----------");
657  PetscInt n=200;
658  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
659  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
660  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
661  //SPI::printf("n = %d",n);
662  //SPI::SPIVec y(SPI::set_FD_stretched_y(21.,n) ,"yFD");
663  //SPI::SPIgrid1D grid(y,"grid",SPI::FD);
664  //grid.print();
665  //(grid.Dy*grid.y).print();
666 
667  SPI::SPIparams params("Blasius parameters");
668  params.Re = 1000.0/1.7208;
669  params.omega = 0.26/1.7208;
670  params.nu = 1./params.Re;
671  params.x = params.Re;
672  params.x_start = params.x;
673  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
674  params.beta = 0.0;
675 
676  SPI::SPIVec eigenfunction(grid.y.rows*8,"q");
677  PetscScalar eigenvalue;
678  eigenfunction.name = "eigenfunction";
679 
680  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
681  if(1){ // set to parallel baseflow
682  SPI::SPIVec o(SPI::zeros(n),"zero");
683  bl_flow.Ux = o;
684  bl_flow.Uxy = o;
685  bl_flow.V = o;
686  bl_flow.Vy = o;
687  }
688  //bl_flow.print();
689 
690  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
691  test_if_close(eigenvalue*1.7208,(0.74155+0.345132*PETSC_i),"LST_spatial 1",1e-5);
692  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LST_spatial 1",1e-5);
693  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
694 
695  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
696  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
697  test_if_close(eigenvalue*1.7208,(0.54213+0.083968*PETSC_i),"LST_spatial 2",1e-5);
698  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
699 
700  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
701  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
702  test_if_close(eigenvalue*1.7208,(0.29967+0.083968*PETSC_i),"LST_spatial 3",1e-5);
703  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
704 
705  SPI::printf("------------ LST Blasius boundary layer end -----------");
706  }
707  if(alltests){
708  SPI::printf("------------ LSTNP parallel Blasius boundary layer start -----------");
709  PetscInt n=1000; // needs 1700 if using finite difference with default stretching...., 1000 if using delta=1.01 stretching
710  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
711  //SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
712  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
713  //SPI::printf("n = %d",n);
714  SPI::SPIVec y(SPI::set_FD_stretched_y(21.,n,1.01) ,"yFD");
715  SPI::SPIgrid1D grid(y,"grid",SPI::FD);
716  //grid.print();
717  //(grid.Dy*grid.y).print();
718 
719  SPI::SPIparams params("Blasius parameters");
720  params.Re = 1000.0/1.7208;
721  params.omega = 0.26/1.7208;
722  params.nu = 1./params.Re;
723  params.x = params.Re;
724  params.x_start = params.x;
725  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
726  params.beta = 0.0;
727 
728  SPI::SPIVec eigenfunction(grid.y.rows*8,"q");
729  PetscScalar eigenvalue;
730  eigenfunction.name = "eigenfunction";
731 
732  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
733  if(1){ // set to parallel baseflow
734  SPI::SPIVec o(SPI::zeros(n),"zero");
735  bl_flow.Ux = o;
736  bl_flow.Uxy = o;
737  bl_flow.V = o;
738  bl_flow.Vy = o;
739  }
740  //bl_flow.print();
741 
742  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
743  test_if_close(eigenvalue*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right 1",1e-5);
744  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right 1",1e-5);
745  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
746 
747  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
748  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
749  test_if_close(eigenvalue*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial_right 2",1e-5);
750  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
751 
752  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
753  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
754  test_if_close(eigenvalue*1.7208,(0.29967+0.083968*PETSC_i),"LSTNP_spatial_right 3",1e-5);
755  SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue*1.7208);
756 
757  SPI::printf("------------ LSTNP parallel Blasius boundary layer end -----------");
758  }
759  if(alltests){
760  SPI::printf("------------ Fourier Collocated derivative operator start -----------");
761  PetscInt n=8; // must be even number for Fourier derivatives
762  SPI::SPIVec t(SPI::set_Fourier_t(2.0*M_PI,n) ,"tFT");
763  SPI::SPIMat Dt(SPI::set_D_Fourier(t),"Dt");
764  SPI::SPIVec y(sin(t));
765  SPI::SPIVec yp(cos(t));
766  test_if_close((Dt*y)(2,PETSC_TRUE),yp(2,PETSC_TRUE),"set_D_Fourier 1",1e-12);
767  test_if_close(SPI::L2((Dt*y)-yp),0.0,"set_D_Fourier 2",1e-12);
768  Dt.print();
769  SPI::printf("------------ Fourier Collocated derivative operator end -----------");
770  }
771  if(alltests){
772  SPI::printf("------------ BV Orthogonalize start -----------");
773  //SPI::SPIMat A(2,2,"A");
774  //A(0,0,1.);
775  //A(0,1,1.);
776  //A(1,0,2.);
777  //A(1,1,1.);
778  //A();
779  //A.print();
780  SPI::SPIVec A1(2,"A1"),A2(2,"A2");
781  A1(0,1.0+0.5*PETSC_i); A2(0,1.);
782  A1(1,2.0+0.5*PETSC_i); A2(1,1.);
783  A1();A2();
784  Vec As[]={A1.vec,A2.vec};
785  // orthogonalize
786  if(0){
787  PetscErrorCode ierr;
788  BV bv;
789  ierr = BVCreate(PETSC_COMM_WORLD,&bv); CHKERRXX(ierr);
790  PetscInt m=2;
791  BVSetSizesFromVec(bv,A1.vec,m);
792  ierr = BVSetFromOptions(bv);CHKERRXX(ierr);
793  ierr = BVInsertVecs(bv,0,&m,As,PETSC_TRUE);
794  SPI::SPIMat AorthH("AorthH");
795  ierr = BVCreateMat(bv,&AorthH.mat); CHKERRXX(ierr);
796  AorthH.rows=A1.rows;
797  AorthH.cols=m;
798  AorthH.print();
799  ierr = BVDestroy(&bv); CHKERRXX(ierr);
800  }
801  SPI::SPIMat Aorth(SPI::orthogonalize({A1,A2}));
802  test_if_close(Aorth(1,0,PETSC_TRUE),0.85280286542244166,"orthogonalize 1",1e-12);
803 
804  SPI::SPIMat Aorth2(2,2,"Aorth2");
805  Aorth2 = SPI::orthogonalize({A1,A2});
806  //Aorth.print();
807  test_if_close(Aorth2(1,0,PETSC_TRUE),0.85280286542244166,"orthogonalize 2",1e-12);
808  //ierr = MatSetType(AorthH.mat,MATMPIAIJ);CHKERRXX(ierr);
809  //AorthH.print();
810  //SPI::SPIMat Aorth(AorthH,"Aorth");
811  //AorthH.H();
812  //AorthH.print();
813  //(AorthH*A1).print();
814  //(AorthH*AorthH).print();
815 
816  SPI::printf("------------ BV Orthogonalize end -----------");
817  }
818  if(alltests){
819  SPI::printf("------------ Gram-Schmidt Orthogonalize start -----------");
820  PetscInt n=8;
821  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
822  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
823  SPI::SPIVec A1(grid.That*(1.0-(grid.y^2)),"A1"),A2(grid.That*(1.0+(grid.y^2)),"A2");
824  std::vector<SPI::SPIVec> A={A1,A2};
825  //SPI::SPIMat Aorth(SPI::orthogonalize(A,grid));
826  //test_if_close(Aorth(2,0,PETSC_TRUE),-0.484122918275927,"orthogonalize 1",1e-12);
827  //test_if_close(Aorth(2,1,PETSC_TRUE),1.082531754730548,"orthogonalize 2",1e-12);
828  std::vector<SPI::SPIVec> Aorth(SPI::orthogonalize(A,grid));
829  test_if_close(Aorth[0](2,PETSC_TRUE),-0.484122918275927,"orthogonalize 1",1e-12);
830  test_if_close(Aorth[1](2,PETSC_TRUE),1.082531754730548,"orthogonalize 2",1e-12);
831 
832  SPI::printf("------------ Gram-Schmidt Orthogonalize end -----------");
833  }
834  if(alltests){
835  SPI::printf("------------ Gram-Schmidt Orthogonalize start -----------");
836  PetscInt n=8;
837  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
838  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
839  SPI::SPIVec A1(grid.That*(1.0-(grid.y^2)),"A1"),A2(grid.That*(1.0+(grid.y^2)),"A2");
840  std::vector<SPI::SPIVec> A={A1,A2};
841  //SPI::SPIMat Aorth(SPI::orthogonalize(A,grid));
842  //test_if_close(Aorth(2,0,PETSC_TRUE),-0.484122918275927,"orthogonalize 1",1e-12);
843  //test_if_close(Aorth(2,1,PETSC_TRUE),1.082531754730548,"orthogonalize 2",1e-12);
844  std::vector<SPI::SPIVec> Aorth(SPI::orthogonalize(A,grid));
845  test_if_close(Aorth[0](2,PETSC_TRUE),-0.484122918275927,"orthogonalize 1",1e-12);
846  test_if_close(Aorth[1](2,PETSC_TRUE),1.082531754730548,"orthogonalize 2",1e-12);
847 
848  // now try multiple length projection and integration
849  SPI::SPIVec A1l(SPI::block({
850  {SPI::diag(A1),grid.O},
851  {grid.O,SPI::diag(A2)}
852  })()*SPI::ones(2*n),"A1l");
853  SPI::SPIVec A2l(SPI::block({
854  {SPI::diag(A2),grid.O},
855  {grid.O,SPI::diag(A1)}
856  })()*SPI::ones(2*n),"A2l");
857  //A1.print();
858  //A2.print();
859  std::vector<SPI::SPIVec> Al = {A1l,A2l};
860  std::vector<SPI::SPIVec> Alorth(SPI::orthogonalize(Al,grid));
861  test_if_close(Alorth[0](2,PETSC_TRUE),-0.228217732293819,"orthogonalize 3",1e-12);
862  test_if_close(Alorth[0](10,PETSC_TRUE),0.228217732293819,"orthogonalize 4",1e-12);
863  test_if_close(Alorth[1](0,PETSC_TRUE),0.714434508311760,"orthogonalize 5",1e-12);
864  test_if_close(Alorth[1](8,PETSC_TRUE),-0.306186217847897,"orthogonalize 6",1e-12);
865  //Alorth[0].print();
866  //Alorth[1].print();
867  SPI::SPIVec A1ll(SPI::block({
868  {SPI::diag(A1),grid.O,grid.O,grid.O},
869  {grid.O,SPI::diag(A2),grid.O,grid.O},
870  {grid.O,grid.O,SPI::diag(A1),grid.O},
871  {grid.O,grid.O,grid.O,SPI::diag(A2)},
872  })()*SPI::ones(4*n),"A1ll");
873  SPI::SPIVec A2ll(SPI::block({
874  {SPI::diag(A2),grid.O,grid.O,grid.O},
875  {grid.O,SPI::diag(A1),grid.O,grid.O},
876  {grid.O,grid.O,SPI::diag(A2),grid.O},
877  {grid.O,grid.O,grid.O,SPI::diag(A1)},
878  })()*SPI::ones(4*n),"A2ll");
879  std::vector<SPI::SPIVec> All = {A1ll,A2ll};
880  std::vector<SPI::SPIVec> Allorth(SPI::orthogonalize(All,grid));
881  test_if_close(Allorth[0](2,PETSC_TRUE),-0.161374306091976,"orthogonalize 7",1e-12);
882  test_if_close(Allorth[0](8,PETSC_TRUE),0.484122918275927,"orthogonalize 8",1e-12);
883  test_if_close(Allorth[0](10,PETSC_TRUE),0.161374306091976,"orthogonalize 9",1e-12);
884  test_if_close(Allorth[0](18,PETSC_TRUE),-0.161374306091976,"orthogonalize 10",1e-12);
885  test_if_close(Allorth[0](24,PETSC_TRUE),0.484122918275927,"orthogonalize 11",1e-12);
886  test_if_close(Allorth[0](26,PETSC_TRUE),0.161374306091976,"orthogonalize 12",1e-12);
887  test_if_close(Allorth[1](2,PETSC_TRUE),0.360843918243516,"orthogonalize 13",1e-12);
888  test_if_close(Allorth[1](8,PETSC_TRUE),-0.216506350946110,"orthogonalize 14",1e-12);
889  test_if_close(Allorth[1](10,PETSC_TRUE),-0.360843918243516,"orthogonalize 15",1e-12);
890  test_if_close(Allorth[1](18,PETSC_TRUE),0.360843918243516,"orthogonalize 16",1e-12);
891  test_if_close(Allorth[1](24,PETSC_TRUE),-0.216506350946110,"orthogonalize 17",1e-12);
892  test_if_close(Allorth[1](26,PETSC_TRUE),-0.360843918243516,"orthogonalize 18",1e-12);
893 
894  SPI::printf("------------ Gram-Schmidt Orthogonalize end -----------");
895  }
896  if(alltests){
897  SPI::printf("------------ Gram-Schmidt Orthogonalize 2D start -----------");
898  PetscInt n=8;
899  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
900  SPI::SPIVec t(SPI::set_Cheby_mapped_y(0.,3.,n) ,"tCheby");
901  SPI::SPIgrid2D grid(y,t,"grid",SPI::Chebyshev,SPI::Chebyshev);
902  SPI::SPIgrid2D gridFD(y,t,"grid",SPI::FD,SPI::FD);
903  SPI::SPIgrid2D gridUltraS(y,t,"grid",SPI::UltraS,SPI::FD);
904 
905  // test integration
906  SPI::SPIVec A1((1.0-(grid.y)),"A1"),A2(grid.grid1Dy.That*(1.0+(grid.y)),"A2");
907  SPI::SPIVec A1l(SPIVec1Dto2D(grid,A1));
908  SPI::SPIVec A2l(SPIVec1Dto2D(grid,A2));
909  //std::cout<<"intChebyshev = "<<SPI::integrate(A1l,grid)<<std::endl;
910  //std::cout<<"intFD = "<<SPI::integrate(A1l,gridFD)<<std::endl;
911  test_if_close(SPI::integrate(A1l,grid),6.0,"integrate SPIgrid2D Chebyshev 1",1e-12);
912  test_if_close(SPI::integrate(A1l,gridFD),6.0,"integrate SPIgrid2D FD 2",1e-12);
913  test_if_close(SPI::integrate(A2l,gridUltraS),6.0,"integrate SPIgrid2D UltraS 3",1e-12);
914  //std::cout<<"intUltraS = "<<SPI::integrate(A2l,gridUltraS)<<std::endl;
915  //std::cout<<"int_1D = "<<SPI::integrate(A1,grid.grid1Dy)<<std::endl;
916  test_if_close(SPI::integrate(A1,grid.grid1Dy),2.0,"integrate SPIgrid1D Chebyshev 4",1e-12);
917 
918  // test projection
919  SPI::SPIVec yp1(y+1.0);
920  SPI::SPIVec nym1(1.0-y);
921  A2l = SPIVec1Dto2D(grid,yp1);
922  SPI::SPIVec A1lp,A2lp;
923  A1lp = SPI::proj(A1l,A2l,grid);
924  //A2lp = SPI::proj(A1l,A2l,gridFD);
925  //A1lp.print();
926  //A2lp.print();
927  //SPI::proj(nym1,yp1,grid.grid1Dy).print();
928  test_if_close(A1lp(1,PETSC_TRUE),SPI::proj(nym1,yp1,grid.grid1Dy)(1,PETSC_TRUE),"proj SPIgrid2D 1",1e-12);
929  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
930  //std::cout<<"norm = "<<std::sqrt(SPI::integrate(SPI::abs(A1lp)^2,grid))<<std::endl;
931  //A1l = SPIVec1Dto2D(grid,nym1);
932  //A2l = SPIVec1Dto2D(grid,yp1);
933  std::vector<SPI::SPIVec> A12ls = {A1l,A2l};
934  std::vector<SPI::SPIVec> A12lso;
935  A12lso = SPI::orthogonalize(A12ls,grid);
936  //for(PetscInt i=0; i<A12lso.size(); ++i){
937  //A12lso[i].print();
938  //}
939  //for(PetscInt i=0; i<A12lso.size(); ++i){
940  //std::cout<<"norm = "<<std::sqrt(SPI::integrate(SPI::abs(A12lso[i])^2,grid))<<std::endl;
941  //}
942  test_if_close(std::sqrt(SPI::integrate(SPI::abs(A12lso[0])^2,grid)),1.0,"norm SPIgrid2D after projection 1",1e-12);
943  test_if_close(std::sqrt(SPI::integrate(SPI::abs(A12lso[1])^2,grid)),1.0,"norm SPIgrid2D after orthogonalize 1",1e-12);
944 
945 
946 
947 
948  /*
949  SPI::SPIVec A1(grid.That*(1.0-(grid.y^2)),"A1"),A2(grid.That*(1.0+(grid.y^2)),"A2");
950  std::vector<SPI::SPIVec> A={A1,A2};
951  std::vector<SPI::SPIVec> Aorth(SPI::orthogonalize(A,grid));
952  test_if_close(Aorth[0](2,PETSC_TRUE),-0.484122918275927,"orthogonalize 1",1e-12);
953  test_if_close(Aorth[1](2,PETSC_TRUE),1.082531754730548,"orthogonalize 2",1e-12);
954 
955  // now try multiple length projection and integration
956  SPI::SPIVec A1l(SPI::block({
957  {SPI::diag(A1),grid.O},
958  {grid.O,SPI::diag(A2)}
959  })()*SPI::ones(2*n),"A1l");
960  SPI::SPIVec A2l(SPI::block({
961  {SPI::diag(A2),grid.O},
962  {grid.O,SPI::diag(A1)}
963  })()*SPI::ones(2*n),"A2l");
964  std::vector<SPI::SPIVec> Al = {A1l,A2l};
965  std::vector<SPI::SPIVec> Alorth(SPI::orthogonalize(Al,grid));
966  test_if_close(Alorth[0](2,PETSC_TRUE),-0.228217732293819,"orthogonalize 3",1e-12);
967  test_if_close(Alorth[0](10,PETSC_TRUE),0.228217732293819,"orthogonalize 4",1e-12);
968  test_if_close(Alorth[1](0,PETSC_TRUE),0.714434508311760,"orthogonalize 5",1e-12);
969  test_if_close(Alorth[1](8,PETSC_TRUE),-0.306186217847897,"orthogonalize 6",1e-12);
970  SPI::SPIVec A1ll(SPI::block({
971  {SPI::diag(A1),grid.O,grid.O,grid.O},
972  {grid.O,SPI::diag(A2),grid.O,grid.O},
973  {grid.O,grid.O,SPI::diag(A1),grid.O},
974  {grid.O,grid.O,grid.O,SPI::diag(A2)},
975  })()*SPI::ones(4*n),"A1ll");
976  SPI::SPIVec A2ll(SPI::block({
977  {SPI::diag(A2),grid.O,grid.O,grid.O},
978  {grid.O,SPI::diag(A1),grid.O,grid.O},
979  {grid.O,grid.O,SPI::diag(A2),grid.O},
980  {grid.O,grid.O,grid.O,SPI::diag(A1)},
981  })()*SPI::ones(4*n),"A2ll");
982  std::vector<SPI::SPIVec> All = {A1ll,A2ll};
983  std::vector<SPI::SPIVec> Allorth(SPI::orthogonalize(All,grid));
984  test_if_close(Allorth[0](2,PETSC_TRUE),-0.161374306091976,"orthogonalize 7",1e-12);
985  test_if_close(Allorth[0](8,PETSC_TRUE),0.484122918275927,"orthogonalize 8",1e-12);
986  test_if_close(Allorth[0](10,PETSC_TRUE),0.161374306091976,"orthogonalize 9",1e-12);
987  test_if_close(Allorth[0](18,PETSC_TRUE),-0.161374306091976,"orthogonalize 10",1e-12);
988  test_if_close(Allorth[0](24,PETSC_TRUE),0.484122918275927,"orthogonalize 11",1e-12);
989  test_if_close(Allorth[0](26,PETSC_TRUE),0.161374306091976,"orthogonalize 12",1e-12);
990  test_if_close(Allorth[1](2,PETSC_TRUE),0.360843918243516,"orthogonalize 13",1e-12);
991  test_if_close(Allorth[1](8,PETSC_TRUE),-0.216506350946110,"orthogonalize 14",1e-12);
992  test_if_close(Allorth[1](10,PETSC_TRUE),-0.360843918243516,"orthogonalize 15",1e-12);
993  test_if_close(Allorth[1](18,PETSC_TRUE),0.360843918243516,"orthogonalize 16",1e-12);
994  test_if_close(Allorth[1](24,PETSC_TRUE),-0.216506350946110,"orthogonalize 17",1e-12);
995  test_if_close(Allorth[1](26,PETSC_TRUE),-0.360843918243516,"orthogonalize 18",1e-12);
996  */
997 
998  SPI::printf("------------ Gram-Schmidt Orthogonalize 2D end -----------");
999  }
1000  if(alltests){
1001  SPI::printf("------------ SPIMat.H(SPIVec) start -----------");
1002  SPI::SPIMat A(4,2,"A");
1003  A(0,0,4.); A(0,1,3.);
1004  A(1,0,3.+4.0*PETSC_i); A(1,1,3.);
1005  A(2,0,2.); A(2,1,3.);
1006  A(3,0,1.); A(3,1,3.);
1007  A();
1008  //A.print();
1009  SPI::SPIVec AHq(A.H(SPI::ones(4)));
1010  test_if_close(AHq(0,PETSC_TRUE),10.0,"SPIMat.H(SPIVec) 1",1e-12);
1011  test_if_close(PetscImaginaryPart(AHq(0,PETSC_TRUE)),-4.0,"SPIMat.H(SPIVec) 2",1e-12);
1012  test_if_close(AHq(1,PETSC_TRUE),12.0,"SPIMat.H(SPIVec) 3",1e-12);
1013  test_if_close(PetscImaginaryPart(AHq(1,PETSC_TRUE)),0.0,"SPIMat.H(SPIVec) 3",1e-12);
1014  SPI::printf("------------ SPIMat.H(SPIVec) end -----------");
1015  }
1016  if(alltests){
1017  SPI::printf("------------ long skinny mat * short vec SPIMat*SPIVec start -----------");
1018  PetscInt ny=3;
1019  SPI::SPIMat A(ny,2,"A");
1020  for(PetscInt i=0; i<ny; ++i){
1021  A(i,0,(PetscScalar)i);
1022  A(i,1,(PetscScalar)(2*i));
1023  }
1024  A();
1025  SPI::SPIVec x(SPI::arange(2));
1026  SPI::SPIVec y(A*x);
1027  test_if_close(y(1,PETSC_TRUE),2.0,"SPIMat*SPIVec 1",1e-12);
1028  SPI::printf("------------ long skinny mat * short vec SPIMat*SPIVec start -----------");
1029  }
1030  if(alltests){
1031  SPI::printf("------------ abs(SPIMat) start -----------");
1032  PetscInt ny=3;
1033  SPI::SPIMat A(ny,2,"A");
1034  for(PetscInt i=0; i<ny; ++i){
1035  A(i,0,-(PetscScalar)i+0.0*PETSC_i);
1036  A(i,1,-(PetscScalar)(2*i) - 2.0*PETSC_i);
1037  }
1038  A();
1039  //A.print();
1040  //SPI::abs(A).print();
1041  test_if_close(SPI::abs(A)(1,1,PETSC_TRUE),sqrt(8.0),"abs(SPIMat) 1",1e-12);
1042  test_if_close(SPI::abs(A)(2,1,PETSC_TRUE),sqrt(20.0),"abs(SPIMat) 2",1e-12);
1043  SPI::printf("------------ abs(SPIMat) end -----------");
1044  }
1045  if(alltests){
1046  SPI::printf("------------ inv(SPIMat) start -----------");
1047  SPI::SPIMat A(3,3,"A");
1048  A(0,0, 0.0); A(0,1,-3.0); A(0,2,-2.0);
1049  A(1,0, 1.0); A(1,1,-4.0); A(1,2,-2.0);
1050  A(2,0,-3.0); A(2,1, 4.0); A(2,2, 1.0);
1051  A();
1052  SPI::SPIVec yy(3);
1053  MatGetColumnVector(A.mat,yy.vec,1);
1054  SPI::SPIMat Ainv = inv(A);
1055  test_if_close(Ainv(2,1,PETSC_TRUE),9.0,"inv(SPIMat) 1",1e-12);
1056  SPI::printf("------------ inv(SPIMat) end -----------");
1057  }
1058  if(alltests){
1059  SPI::printf("------------ UltraSpherical ops start -----------");
1060  //SPI::SPIVec y(SPI::set_Cheby_y(5),"yCheby");
1061  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0,10,5),"yCheby");
1062  SPI::SPIMat S0(5,5,"S0"),D1(5,5,"D1");
1063  SPI::SPIMat S1(5,5,"S1"),D2(5,5,"D2");
1064  std::tie(S0,D1) = SPI::set_D_UltraS(y,1);
1065  std::tie(S1,D2) = SPI::set_D_UltraS(y,2);
1066  //S0.print();
1067  //(S1*D1).print();
1068  //(D2).print();
1069  SPI::SPIMat T(5,5,"T"),That(5,5,"That");
1070  std::tie(T,That) = SPI::set_T_That(y.rows);
1071  //T.print();
1072  //That.print();
1073  //(T*inv(S0)*D1*That).print();
1074  test_if_close((S1*D1)(1,2,PETSC_TRUE),0.2,"set_D_UltraS(SPIMat) 1",1e-12);
1075  test_if_close((S1*D1)(1,4,PETSC_TRUE),-0.2,"set_D_UltraS(SPIMat) 2",1e-12);
1076  test_if_close(D2(1,3,PETSC_TRUE),0.24,"set_D_UltraS(SPIMat) 3",1e-12);
1077  test_if_close((T*inv(S0)*D1*That)(1,2,PETSC_TRUE),0.282842712,"set_D_UltraS(SPIMat) 4",1e-7);
1078  test_if_close((T*inv(S0)*inv(S1)*D2*That)(1,3,PETSC_TRUE),-0.08,"set_D_UltraS(SPIMat) 4",1e-7);
1079  //(T*inv(S0)*inv(S1)*D2*That).print();
1080  //((T*inv(S0)*inv(S1)*D2*That)*y).print();
1081  //inv(S0).print();
1082  //inv(S1).print();
1083  //(inv(S0)*inv(S1)).print();
1084  test_if_close((inv(S0)*inv(S1))(2,4,PETSC_TRUE),16.0,"inv(S0)*inv(S1) 1",1e-7);
1085  SPI::printf("------------ UltraSpherical ops end -----------");
1086  }
1087  if(alltests){
1088  SPI::printf("------------ LST_temporal channel UltraS start -----------");
1089  PetscInt n=64;
1090  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1091  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1092  //grid.print();
1093  //exit(0);
1094  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1095  //SPI::printf("n = %d",n);
1096  //grid.print();
1097  //(grid.Dy*grid.y).print();
1098 
1099  SPI::SPIparams params("Channel parameters");
1100  params.Re = 2000.0;
1101  params.omega = 0.3121-0.01978*PETSC_i;
1102  params.nu = 1./params.Re;
1103  params.x = 1.0;
1104  params.x_start = params.x;
1105  params.alpha = 1.0;
1106  params.beta = 0.0;
1107 
1108  SPI::SPIVec eigenfunction(grid.y.rows*4,"q");
1109  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1110  PetscScalar eigenvalue;
1111  eigenfunction.name = "eigenfunction";
1112 
1113  //SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1114  SPI::SPIVec U((1.0-((grid.y)^2)),"U");
1115  SPI::SPIVec Uy((-2.*grid.y),"Uy");
1116  SPI::SPIVec o(U*0.0,"o"); // zero vector
1117  SPI::SPIbaseflow channel(U,o,o,o,Uy,o,o,o,o,o,o,o,o,o);
1118  if(0){ // set to parallel baseflow
1119  //SPI::SPIVec o(SPI::zeros(n),"zero");
1120  //bl_flow.Ux = o;
1121  //bl_flow.Uxy = o;
1122  //bl_flow.V = o;
1123  //bl_flow.Vy = o;
1124  }
1125  //bl_flow.print();
1126 
1127  PetscScalar cg;
1128  std::tie(eigenvalue,eigenfunction) = SPI::LST_temporal(params,grid,channel);
1129  SPI::printfc("eigenvalue = %g+%gi",eigenvalue);
1130  //std::tie(eigenvalue,eig_vec) = SPI::LST_spatial(params,grid,channel);
1131  test_if_close(eigenvalue,0.3121-0.01978*PETSC_i,"LST_temporal UltraS 1",1e-5);
1132  test_if_close(params.omega,0.3121+0.01978*PETSC_i,"LST_temporal UltraS 2",1e-5);
1133  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,channel);
1134  //test_if_close(params.alpha,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 1",1e-5);
1135  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
1136 
1137  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
1138  //SPI::printfc("cg = %g+%gi",cg);
1139  //test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial 2",1e-5);
1140  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
1141  //test_if_close(eigenvalue,0.094966+0.004564*PETSC_i,"LSTNP_spatial_right 2",1e-5);
1142  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
1143 
1144  //params.alpha = 0.1067+0.001898*PETSC_i;
1145  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1146  //SPI::printfc("cg = %g+%gi",cg);
1147  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 3",1e-5);
1148  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1149  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 3",1e-5);
1150  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
1151  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow,leigenfunction,eigenfunction); // with initial guess
1152  //SPI::printfc("cg = %g+%gi",cg);
1153  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial 4",1e-5);
1154  //SPI::printfc("eigenvalue is %.10f + %.10fi",eigenvalue);
1155  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow,eigenfunction); // with initial guess
1156  //test_if_close(eigenvalue,0.10665+0.0018979*PETSC_i,"LSTNP_spatial_right 4",1e-5);
1157 
1158  SPI::printf("------------ LST_temporal channel UltraS end -----------");
1159  }
1160  if(alltests){
1161  SPI::printf("------------ LST_spatial channel flow UltraS start -----------");
1162  PetscInt n=64;
1163  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1164  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1165  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1166  //grid.print();
1167  //exit(0);
1168  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1169  //SPI::printf("n = %d",n);
1170  //grid.print();
1171  //(grid.Dy*grid.y).print();
1172 
1173  SPI::SPIparams params("Blasius parameters");
1174  //params.Re = 400.0;
1175  //params.omega = 86.*params.Re/(1000000.);
1176  //params.Re = 1000.0/1.7208;
1177  //params.omega = 0.26/1.7208;
1178  params.Re = 2000.0;
1179  params.omega = 0.3;
1180  params.nu = 1./params.Re;
1181  params.x = params.Re;
1182  params.x_start = params.x;
1183  //params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1184  params.alpha = 0.97875+0.044394*PETSC_i;
1185  params.beta = 0.0;
1186 
1187  SPI::SPIVec eigenfunction(grid.y.rows*8,"q");
1188  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1189  SPI::SPIVec leigenfunction(grid.y.rows*8,"q");
1190  PetscScalar eigenvalue;
1191  eigenfunction.name = "eigenfunction";
1192 
1193  //SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1194  SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1195  if(0){ // set to parallel baseflow
1196  SPI::SPIVec o(SPI::zeros(n),"zero");
1197  bl_flow.Ux = o;
1198  bl_flow.Uxy = o;
1199  bl_flow.V = o;
1200  bl_flow.Vy = o;
1201  }
1202 
1203  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1204  test_if_close(params.alpha,(0.97875+0.044394*PETSC_i),"LST_spatial 1",1e-5);
1205 
1206  params.alpha = 0.34312+0.049677*PETSC_i;
1207  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1208  test_if_close(params.alpha,(0.34312+0.049677*PETSC_i),"LST_spatial 2",1e-5);
1209 
1210  params.alpha = 0.61167+0.140492*PETSC_i;
1211  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1212  test_if_close(params.alpha,(0.61167+0.140492*PETSC_i),"LST_spatial 3",1e-5);
1213 
1214  SPI::printf("------------ LST_spatial channel flow UltraS end -----------");
1215  }
1216  if(alltests){ // requires -mat_mumps_icntl_14 25 in command line call to work
1217  SPI::printf("------------ LST_spatial Blasius boundary layer UltraS start -----------");
1218  PetscInt n=200;
1219  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1220  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1221  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1222  //grid.print();
1223  //exit(0);
1224  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1225  //SPI::printf("n = %d",n);
1226  //grid.print();
1227  //(grid.Dy*grid.y).print();
1228 
1229  SPI::SPIparams params("Blasius parameters");
1230  //params.Re = 400.0;
1231  //params.omega = 86.*params.Re/(1000000.);
1232  params.Re = 1000.0/1.7208;
1233  params.omega = 0.26/1.7208;
1234  params.nu = 1./params.Re;
1235  params.x = params.Re;
1236  params.x_start = params.x;
1237  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1238  params.beta = 0.0;
1239  //params.print();
1240 
1241  SPI::SPIVec eigenfunction(grid.y.rows*8,"q");
1242  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1243  SPI::SPIVec leigenfunction(grid.y.rows*8,"q");
1244  PetscScalar eigenvalue;
1245  eigenfunction.name = "eigenfunction";
1246 
1247  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1248  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1249  if(1){ // set to parallel baseflow
1250  SPI::SPIVec o(SPI::zeros(n),"zero");
1251  bl_flow.Ux = o;
1252  bl_flow.Uxy = o;
1253  bl_flow.V = o;
1254  bl_flow.Vy = o;
1255  }
1256  //bl_flow.print();
1257 
1258  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1259  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1260  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LST_spatial 1",1e-5);
1261 
1262  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1263  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1264  test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LST_spatial 2",1e-5);
1265 
1266  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1267  std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1268  test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LST_spatial 3",1e-5);
1269 
1270  SPI::printf("------------ LST_spatial Blasius boundary layer UltraS end -----------");
1271  }
1272  if(0){
1273  // noticed block and eye commands were the slowest part of the LSTNP_spatial_right eigenvalue solver, so let's speed those up here
1274  SPI::printf("------------ eye and block timings start -----------");
1275  PetscInt n=900000;
1276  SPI::SPIMat I(SPI::eye(n*8),"I");
1277  SPI::SPIMat O(SPI::zeros(n*8,n*8),"O");
1278  std::cout<<"created I and O"<<std::endl;
1279  SPI::block({
1280  {I,O,O},
1281  {I,-I,I},
1282  {I,I,O}
1283  })();
1284  std::cout<<"created block"<<std::endl;
1285  SPI::block({
1286  {I,O,O},
1287  {O,O,O},
1288  {O,O,O}
1289  })();
1290  std::cout<<"created block P"<<std::endl;
1291  SPI::SPIMat B(3,3);
1292  B(0,0,1.0)();
1293  B(0,1,2.0)();
1294  SPI::kron(B,I);
1295  std::cout<<"created block one kron"<<std::endl;
1296  SPI::SPIMat B1(3,3);
1297  B1(0,0,1.0)();
1298  SPI::SPIMat B2(3,3);
1299  B2(0,1,2.0)();
1300  SPI::kron(B1,I) + SPI::kron(B2,I);
1301  std::cout<<"created block using 2 kron"<<std::endl;
1302 
1303  SPI::printf("------------ eye and block timings end -----------");
1304  }
1305  if(alltests){
1306  SPI::printf("------------ LSTNP_spatial_right Blasius boundary layer UltraS start -----------");
1307  PetscInt n=169;
1308  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1309  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1310  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1311  //grid.print();
1312  //exit(0);
1313  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1314  //SPI::printf("n = %d",n);
1315  //grid.print();
1316  //(grid.Dy*grid.y).print();
1317 
1318  SPI::SPIparams params("Blasius parameters");
1319  //params.Re = 400.0;
1320  //params.omega = 86.*params.Re/(1000000.);
1321  params.Re = 1000.0/1.7208;
1322  params.omega = 0.26/1.7208;
1323  params.nu = 1./params.Re;
1324  params.x = params.Re;
1325  params.x_start = params.x;
1326  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1327  params.beta = 0.0;
1328  //params.print();
1329 
1330  SPI::SPIVec eigenfunction(grid.y.rows*8,"q");
1331  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1332  //SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
1333  PetscScalar eigenvalue;
1334  eigenfunction.name = "eigenfunction";
1335 
1336  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1337  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1338  if(1){ // set to parallel baseflow
1339  SPI::SPIVec o(SPI::zeros(n),"zero");
1340  bl_flow.Ux = o;
1341  bl_flow.Uxy = o;
1342  bl_flow.V = o;
1343  bl_flow.Vy = o;
1344  }
1345  //bl_flow.print();
1346 
1347  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1348  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1349  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right 1",1e-5);
1350 
1351  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1352  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1353  test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial_right 2",1e-5);
1354 
1355  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1356  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1357  test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatial_right 3",1e-5);
1358 
1359  SPI::printf("------------ LSTNP_spatial_right Blasius boundary layer UltraS end -----------");
1360  }
1361  if(alltests){
1362  SPI::printf("------------ LST_spatial_cg Blasius boundary layer physical vs UltraS start -----------");
1363  PetscInt n=169;
1364  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1365  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1366  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1367  SPI::SPIgrid1D grid2(y,"grid",SPI::Chebyshev);
1368  //grid.print();
1369  //exit(0);
1370  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1371  //SPI::printf("n = %d",n);
1372  //grid.print();
1373  //(grid.Dy*grid.y).print();
1374 
1375  SPI::SPIparams params("Blasius parameters");
1376  //params.Re = 400.0;
1377  //params.omega = 86.*params.Re/(1000000.);
1378  params.Re = 1000.0/1.7208;
1379  params.omega = 0.26/1.7208;
1380  params.nu = 1./params.Re;
1381  params.x = params.Re;
1382  params.x_start = params.x;
1383  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1384  params.beta = 0.0;
1385  //params.print();
1386 
1387  SPI::SPIVec eigenfunction(grid.y.rows*8,"eigenfunction");
1388  SPI::SPIVec leigenfunction(grid.y.rows*8,"leigenfunction");
1389  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1390  //SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
1391  PetscScalar eigenvalue;
1392  eigenfunction.name = "eigenfunction";
1393 
1394  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1395  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1396  if(0){ // set to parallel baseflow
1397  SPI::SPIVec o(SPI::zeros(n),"zero");
1398  bl_flow.Ux = o;
1399  bl_flow.Uxy = o;
1400  bl_flow.V = o;
1401  bl_flow.Vy = o;
1402  }
1403  //bl_flow.print();
1404 
1405  // LST_spatial_cg UltraS
1406  PetscScalar cg;
1407  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LST_spatial_cg(params,grid,bl_flow);
1408  SPI::printfc("cg = %.10f + %.10fi",cg);
1409  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1410  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LST_spatial_cg UltraS 1",1e-5);
1411  SPI::SPIMat &O = grid.O;
1412  SPI::SPIMat T(SPI::block({
1413  {grid.T,O,O,O,O,O,O,O},
1414  {O,grid.T,O,O,O,O,O,O},
1415  {O,O,grid.T,O,O,O,O,O},
1416  {O,O,O,grid.T,O,O,O,O},
1417  {O,O,O,O,grid.T,O,O,O},
1418  {O,O,O,O,O,grid.T,O,O},
1419  {O,O,O,O,O,O,grid.T,O},
1420  {O,O,O,O,O,O,O,grid.T},
1421  })());
1422  SPI::SPIVec eigenfunctionout(T*eigenfunction,"UltraS");
1423  SPI::SPIVec leigenfunctionout(T*leigenfunction,"UltraSl");
1424  // normalize
1425  eigenfunctionout /= SPI::L2(eigenfunctionout);
1426  leigenfunctionout /= SPI::L2(leigenfunctionout);
1427  SPI::save(eigenfunctionout,"eigenfunction.h5");
1428  SPI::save(leigenfunctionout,"eigenfunction.h5");
1429  // LST_spatial physical
1430  SPI::SPIVec leigenfunction8(8*n,"Physicall");
1431  SPI::SPIVec eigenfunction8(8*n,"Physical");
1432  std::tie(eigenvalue,cg,leigenfunction8,eigenfunction8) = SPI::LST_spatial_cg(params,grid2,bl_flow);
1433  SPI::printfc("cg = %.10f + %.10fi",cg);
1434  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LST_spatial_cg physical 1",1e-5);
1435  SPI::save(eigenfunction8,"eigenfunction.h5");
1436  SPI::save(leigenfunction8,"eigenfunction.h5");
1437  SPI::save(grid.y,"eigenfunction.h5");
1438  // finite difference what is the group velocity of this mode?
1439  PetscScalar deltaOmega = 1e-2;
1440  params.omega += deltaOmega;
1441  PetscScalar deltaeigenvalue;
1442  std::tie(deltaeigenvalue,cg,leigenfunction8,eigenfunction8) = SPI::LST_spatial_cg(params,grid2,bl_flow);
1443  SPI::printfc("with delta_omega cg = %.10f + %.10fi",cg);
1444  cg = deltaOmega/(deltaeigenvalue-eigenvalue);
1445  SPI::printfc("with finite diff delta_omega cg = %.10f + %.10fi",cg);
1446  params.omega -= 2.0*deltaOmega;
1447  PetscScalar deltaeigenvaluem;
1448  std::tie(deltaeigenvaluem,cg,leigenfunction8,eigenfunction8) = SPI::LST_spatial_cg(params,grid,bl_flow);
1449  cg = (2.0*deltaOmega)/(deltaeigenvalue-deltaeigenvaluem);
1450  SPI::printfc("with finite diff delta_omega2 cg = %.10f + %.10fi",cg);
1451  // LSTNP_spatial physical
1452  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid2,bl_flow);
1453  //SPI::printfc("cg = %.10f + %.10fi",cg);
1454  //test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial physical 1",1e-4);
1455 
1456  // these following 2 cases don't work
1457  //params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1458  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1459  //SPI::printfc("cg = %.10f + %.10fi",cg);
1460  //test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial UltraS 2",1e-5);
1461  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LST_spatial_cg(params,grid2,bl_flow);
1462  //SPI::printfc("cg = %.10f + %.10fi",cg);
1463  //test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LST_spatial_cg 2",1e-5);
1464 
1465  //params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1466  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1467  //SPI::printfc("cg = %.10f + %.10fi",cg);
1468  //test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatial UltraS 3",1e-5);
1469  //std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LST_spatial_cg(params,grid,bl_flow);
1470  //SPI::printfc("cg = %.10f + %.10fi",cg);
1471  //test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LST_spatial_cg 3",1e-5);
1472 
1473  SPI::printf("------------ LST_spatial_cg Blasius boundary layer UltraS end -----------");
1474  }
1475  if(alltests){
1476  SPI::printf("------------ LSTNP_spatial Blasius boundary layer UltraS start -----------");
1477  PetscInt n=169;
1478  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1479  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1480  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1481  SPI::SPIgrid1D grid2(y,"grid",SPI::Chebyshev);
1482  //grid.print();
1483  //exit(0);
1484  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1485  //SPI::printf("n = %d",n);
1486  //grid.print();
1487  //(grid.Dy*grid.y).print();
1488 
1489  SPI::SPIparams params("Blasius parameters");
1490  //params.Re = 400.0;
1491  //params.omega = 86.*params.Re/(1000000.);
1492  params.Re = 1000.0/1.7208;
1493  params.omega = 0.26/1.7208;
1494  params.nu = 1./params.Re;
1495  params.x = params.Re;
1496  params.x_start = params.x;
1497  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1498  params.beta = 0.0;
1499  //params.print();
1500 
1501  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
1502  SPI::SPIVec leigenfunction(grid.y.rows*16,"leigenfunction");
1503  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1504  //SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
1505  PetscScalar eigenvalue;
1506  eigenfunction.name = "eigenfunction";
1507 
1508  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1509  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1510  if(1){ // set to parallel baseflow
1511  SPI::SPIVec o(SPI::zeros(n),"zero");
1512  bl_flow.Ux = o;
1513  bl_flow.Uxy = o;
1514  bl_flow.V = o;
1515  bl_flow.Vy = o;
1516  }
1517  //bl_flow.print();
1518 
1519  PetscScalar cg;
1520  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1521  SPI::printfc("cg = %.10f + %.10fi",cg);
1522  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1523  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial UltraS 1",1e-4);
1524  SPI::SPIMat &O = grid.O;
1525  SPI::SPIMat &t = grid.T;
1526  SPI::SPIMat T(SPI::block({
1527  {t,O,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1528  {O,t,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1529  {O,O,t,O,O,O,O,O,O,O,O,O,O,O,O,O},
1530  {O,O,O,t,O,O,O,O,O,O,O,O,O,O,O,O},
1531  {O,O,O,O,t,O,O,O,O,O,O,O,O,O,O,O},
1532  {O,O,O,O,O,t,O,O,O,O,O,O,O,O,O,O},
1533  {O,O,O,O,O,O,t,O,O,O,O,O,O,O,O,O},
1534  {O,O,O,O,O,O,O,t,O,O,O,O,O,O,O,O},
1535  {O,O,O,O,O,O,O,O,t,O,O,O,O,O,O,O},
1536  {O,O,O,O,O,O,O,O,O,t,O,O,O,O,O,O},
1537  {O,O,O,O,O,O,O,O,O,O,t,O,O,O,O,O},
1538  {O,O,O,O,O,O,O,O,O,O,O,t,O,O,O,O},
1539  {O,O,O,O,O,O,O,O,O,O,O,O,t,O,O,O},
1540  {O,O,O,O,O,O,O,O,O,O,O,O,O,t,O,O},
1541  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,t,O},
1542  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,O,t},
1543  })());
1544  SPI::SPIVec eigenfunctionout(T*eigenfunction,"UltraS");
1545  SPI::SPIVec leigenfunctionout(T*leigenfunction,"UltraSl");
1546  //SPI::save(eigenfunctionout,"eigenfunctionNP.h5");
1547  //SPI::save(leigenfunctionout,"eigenfunctionNP.h5");
1548  //SPI::save(grid.y,"eigenfunctionNP.h5");
1549  // LST_spatial physical
1550  SPI::SPIVec leigenfunction16(16*n);
1551  SPI::SPIVec eigenfunction16(16*n);
1552  std::tie(eigenvalue,cg,leigenfunction16,eigenfunction16) = SPI::LSTNP_spatial(params,grid2,bl_flow);
1553  SPI::printfc("cg = %.10f + %.10fi",cg);
1554  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial physical 1",1e-4);
1555  // LSTNP_spatial physical
1556  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid2,bl_flow);
1557  SPI::printfc("cg = %.10f + %.10fi",cg);
1558  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial physical 1",1e-4);
1559 
1560  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1561  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1562  SPI::printfc("cg = %.10f + %.10fi",cg);
1563  test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial UltraS 2",1e-5);
1564 
1565  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1566  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
1567  SPI::printfc("cg = %.10f + %.10fi",cg);
1568  test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatial UltraS 3",1e-5);
1569 
1570  SPI::printf("------------ LSTNP_spatial Blasius boundary layer UltraS end -----------");
1571  }
1572  if(alltests){
1573  SPI::printf("------------ LSTNP_spatial_right Blasius boundary layer UltraS start -----------");
1574  PetscInt n=169;
1575  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1576  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1577  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1578  SPI::SPIgrid1D grid2(y,"grid",SPI::Chebyshev);
1579  //grid.print();
1580  //exit(0);
1581  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1582  //SPI::printf("n = %d",n);
1583  //grid.print();
1584  //(grid.Dy*grid.y).print();
1585 
1586  SPI::SPIparams params("Blasius parameters");
1587  //params.Re = 400.0;
1588  //params.omega = 86.*params.Re/(1000000.);
1589  params.Re = 1000.0/1.7208;
1590  params.omega = 0.26/1.7208;
1591  params.nu = 1./params.Re;
1592  params.x = params.Re;
1593  params.x_start = params.x;
1594  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1595  params.beta = 0.0;
1596  //params.print();
1597 
1598  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
1599  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1600  //SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
1601  PetscScalar eigenvalue;
1602  eigenfunction.name = "eigenfunction";
1603 
1604  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1605  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1606  if(1){ // set to parallel baseflow
1607  SPI::SPIVec o(SPI::zeros(n),"zero");
1608  bl_flow.Ux = o;
1609  bl_flow.Uxy = o;
1610  bl_flow.V = o;
1611  bl_flow.Vy = o;
1612  }
1613  //bl_flow.print();
1614 
1615  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1616  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1617  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right UltraS 1",1e-4);
1618  SPI::SPIMat &O = grid.O;
1619  SPI::SPIMat &t = grid.T;
1620  SPI::SPIMat T(SPI::block({
1621  {t,O,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1622  {O,t,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1623  {O,O,t,O,O,O,O,O,O,O,O,O,O,O,O,O},
1624  {O,O,O,t,O,O,O,O,O,O,O,O,O,O,O,O},
1625  {O,O,O,O,t,O,O,O,O,O,O,O,O,O,O,O},
1626  {O,O,O,O,O,t,O,O,O,O,O,O,O,O,O,O},
1627  {O,O,O,O,O,O,t,O,O,O,O,O,O,O,O,O},
1628  {O,O,O,O,O,O,O,t,O,O,O,O,O,O,O,O},
1629  {O,O,O,O,O,O,O,O,t,O,O,O,O,O,O,O},
1630  {O,O,O,O,O,O,O,O,O,t,O,O,O,O,O,O},
1631  {O,O,O,O,O,O,O,O,O,O,t,O,O,O,O,O},
1632  {O,O,O,O,O,O,O,O,O,O,O,t,O,O,O,O},
1633  {O,O,O,O,O,O,O,O,O,O,O,O,t,O,O,O},
1634  {O,O,O,O,O,O,O,O,O,O,O,O,O,t,O,O},
1635  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,t,O},
1636  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,O,t},
1637  })());
1638  SPI::SPIVec eigenfunctionout(T*eigenfunction,"UltraS");
1639  //SPI::save(eigenfunctionout,"eigenfunctionNP.h5");
1640  //SPI::save(grid.y,"eigenfunctionNP.h5");
1641  // LST_spatial physical
1642  SPI::SPIVec leigenfunction16(16*n);
1643  SPI::SPIVec eigenfunction16(16*n);
1644  std::tie(eigenvalue,eigenfunction16) = SPI::LSTNP_spatial_right(params,grid2,bl_flow);
1645  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right physical 1",1e-4);
1646  // LSTNP_spatial physical
1647  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid2,bl_flow);
1648  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right physical 1",1e-4);
1649 
1650  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1651  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1652  test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial_right UltraS 2",1e-5);
1653 
1654  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1655  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1656  test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatial_right UltraS 3",1e-5);
1657 
1658  SPI::printf("------------ LSTNP_spatial_right Blasius boundary layer UltraS end -----------");
1659  }
1660  if(alltests){
1661  SPI::printf("------------ LSTNP_spatials_right Blasius boundary layer UltraS start -----------");
1662  PetscInt n=169;
1663  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1664  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1665  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1666  //SPI::SPIgrid1D grid2(y,"grid",SPI::Chebyshev);
1667  //grid.print();
1668  //exit(0);
1669  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1670  //SPI::printf("n = %d",n);
1671  //grid.print();
1672  //(grid.Dy*grid.y).print();
1673 
1674  SPI::SPIparams params("Blasius parameters");
1675  //params.Re = 400.0;
1676  //params.omega = 86.*params.Re/(1000000.);
1677  params.Re = 1000.0/1.7208;
1678  params.omega = 0.26/1.7208;
1679  params.nu = 1./params.Re;
1680  params.x = params.Re;
1681  params.x_start = params.x;
1682  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1683  params.beta = 0.0;
1684  //params.print();
1685 
1686  std::vector<SPI::SPIVec> eigenfunction(1); //grid.y.rows*16,"eigenfunction");
1687  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1688  eigenfunction[0] = SPI::ones(grid.y.rows*16);
1689  std::vector<PetscScalar> eigenvalue(1);
1690  std::vector<PetscScalar> eigenvalues;
1691  std::vector<SPI::SPIVec> eigenfunctions;
1692  eigenvalue[0] = params.alpha;
1693  //eigenfunction[0].name = "eigenfunction";
1694 
1695  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1696  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1697  if(1){ // set to parallel baseflow
1698  SPI::SPIVec o(SPI::zeros(n),"zero");
1699  bl_flow.Ux = o;
1700  bl_flow.Uxy = o;
1701  bl_flow.V = o;
1702  bl_flow.Vy = o;
1703  }
1704  //bl_flow.print();
1705 
1706  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatials_right(params,grid,bl_flow,eigenvalue,eigenfunction);
1707  eigenvalues.push_back(params.alpha);
1708  eigenfunctions.push_back(eigenfunction[0]);
1709  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1710  test_if_close(eigenvalue[0]*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatials_right UltraS 1",1e-4);
1711  //std::cout<<"eigenvalue = "<<eigenvalue[0]*1.7208<<std::endl;
1712  SPI::SPIMat &O = grid.O;
1713  SPI::SPIMat &t = grid.T;
1714  SPI::SPIMat T(SPI::block({
1715  {t,O,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1716  {O,t,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
1717  {O,O,t,O,O,O,O,O,O,O,O,O,O,O,O,O},
1718  {O,O,O,t,O,O,O,O,O,O,O,O,O,O,O,O},
1719  {O,O,O,O,t,O,O,O,O,O,O,O,O,O,O,O},
1720  {O,O,O,O,O,t,O,O,O,O,O,O,O,O,O,O},
1721  {O,O,O,O,O,O,t,O,O,O,O,O,O,O,O,O},
1722  {O,O,O,O,O,O,O,t,O,O,O,O,O,O,O,O},
1723  {O,O,O,O,O,O,O,O,t,O,O,O,O,O,O,O},
1724  {O,O,O,O,O,O,O,O,O,t,O,O,O,O,O,O},
1725  {O,O,O,O,O,O,O,O,O,O,t,O,O,O,O,O},
1726  {O,O,O,O,O,O,O,O,O,O,O,t,O,O,O,O},
1727  {O,O,O,O,O,O,O,O,O,O,O,O,t,O,O,O},
1728  {O,O,O,O,O,O,O,O,O,O,O,O,O,t,O,O},
1729  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,t,O},
1730  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,O,t},
1731  })());
1732  SPI::SPIVec eigenfunctionout(T*eigenfunction[0],"UltraS");
1733  SPI::save(eigenfunctionout,"eigenfunctionNP.h5");
1734  SPI::save(grid.y,"eigenfunctionNP.h5");
1735  // LST_spatial physical
1736  //std::vector<SPI::SPIVec> eigenfunction16;
1737  //std::tie(eigenvalue,eigenfunction16) = SPI::LSTNP_spatials_right(params,grid2,bl_flow);
1738  //test_if_close(eigenvalue[0]*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatials_right physical 1",1e-4);
1739  //std::cout<<"eigenvalue = "<<eigenvalue<<std::endl;
1740  // LSTNP_spatial physical
1741  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatials_right(params,grid2,bl_flow);
1742  //test_if_close(eigenvalue[0]*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatials_right physical 1",1e-4);
1743  //std::cout<<"eigenvalue = "<<eigenvalue<<std::endl;
1744 
1745  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1746  eigenvalue[0] = params.alpha;
1747  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatials_right(params,grid,bl_flow,eigenvalue,eigenfunction);
1748  eigenvalues.push_back(params.alpha);
1749  eigenfunctions.push_back(eigenfunction[0]);
1750  test_if_close(eigenvalue[0]*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatials_right UltraS 2",1e-5);
1751  //std::cout<<"eigenvalue = "<<eigenvalue[0]*1.7208<<std::endl;
1752 
1753  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1754  eigenvalue[0] = params.alpha;
1755  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatials_right(params,grid,bl_flow,eigenvalue,eigenfunction);
1756  eigenvalues.push_back(params.alpha);
1757  eigenfunctions.push_back(eigenfunction[0]);
1758  test_if_close(eigenvalue[0]*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatials_right UltraS 3",1e-4);
1759  //std::cout<<"eigenvalue = "<<eigenvalue[0]*1.7208<<std::endl;
1760  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatials_right(params,grid,bl_flow,eigenvalue,eigenfunction);
1761  //test_if_close(eigenvalue[0]*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatials_right UltraS 3",1e-5);
1762  //std::cout<<"eigenvalue = "<<eigenvalue[0]*1.7208<<std::endl;
1763 
1764  std::vector<PetscScalar> eigenvalues2;
1765  std::vector<SPI::SPIVec> eigenfunctions2;
1766  std::tie(eigenvalues2,eigenfunctions2) = SPI::LSTNP_spatials_right(params,grid,bl_flow,eigenvalues,eigenfunctions);
1767  test_if_close(eigenvalues2[0]*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatials_right UltraS 1",1e-4);
1768  test_if_close(eigenvalues2[1]*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatials_right UltraS 2",1e-5);
1769  test_if_close(eigenvalues2[2]*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatials_right UltraS 3",1e-5);
1770 
1771  SPI::printf("------------ LSTNP_spatials_right Blasius boundary layer UltraS end -----------");
1772  }
1773  if(alltests){
1774  SPI::printf("------------ UltraS int start -----------");
1775  PetscInt n=11;
1776  SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1777  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1778  SPI::SPIVec a(grid.That*(y^2));
1779  SPI::SPIVec a2(y^2);
1780  //std::cout<<"integrate = "<<integrate_coeffs(a)<<std::endl;
1781  test_if_close(integrate_coeffs(a),2.0/3.0,"integrate_coeffs 1",1e-12);
1782  test_if_close(integrate_coeffs(a),2.0/3.0,"integrate_coeffs 1",1e-12);
1783  a = (grid.That*((y^2)+(4.0*y)+1.0));
1784  test_if_close(integrate_coeffs(a),8.0/3.0,"integrate_coeffs 2",1e-12);
1785  a = (grid.That*((y^2)+(4.0*y)));
1786  test_if_close(integrate_coeffs(a),2.0/3.0,"integrate_coeffs 3",1e-12);
1787  a = (grid.That*(4.0*(y^2)+(y)+10.0));
1788  test_if_close(integrate_coeffs(a),68.0/3.0,"integrate_coeffs 4",1e-12);
1789 
1790  SPI::SPIVec y2(SPI::set_Cheby_mapped_y(0.,10.,n) ,"yCheby");
1791  SPI::SPIgrid1D grid2(y2,"grid",SPI::UltraS);
1792  a = grid2.That*(4.0*(y2^2)+(y2)+10.0);
1793  test_if_close(integrate_coeffs(a,grid2.y),4450.0/3.0,"integrate_coeffs 5",1e-12);
1794 
1795  SPI::SPIVec y3(SPI::set_Cheby_mapped_y(-10.,10.,n) ,"yCheby");
1796  SPI::SPIgrid1D grid3(y3,"grid",SPI::UltraS);
1797  a = grid3.That*(4.0*(y3^2)+(y3)+10.0);
1798  test_if_close(integrate_coeffs(a,grid3.y),8600.0/3.0,"integrate_coeffs 6",1e-12);
1799 
1800  SPI::SPIVec y4(SPI::set_Cheby_mapped_y(-1.,10.,n) ,"yCheby");
1801  SPI::SPIgrid1D grid4(y4,"grid",SPI::UltraS);
1802  a = grid4.That*(4.0*(y4^2)+(y4)+10.0);
1803  test_if_close(integrate_coeffs(a,grid4.y),8965.0/6.0,"integrate_coeffs 7",1e-12);
1804 
1805  SPI::SPIVec y5(SPI::set_Cheby_mapped_y(-1.,10.,n) ,"yCheby");
1806  SPI::SPIgrid1D grid5(y5,"grid",SPI::UltraS);
1807  SPI::SPIgrid1D grid5_2(y5,"grid",SPI::Chebyshev);
1808  a = grid5.That*(4.0*(y5^2)+(y5)+10.0);
1809  a2 = (4.0*(y5^2)+(y5)+10.0);
1810  test_if_close(integrate(a,grid5),8965.0/6.0,"integrate 1 UltraS",1e-12);
1811  test_if_close(integrate(a2,grid5_2),8965.0/6.0,"integrate 1 Chebyshev",1e-12);
1812 
1813  SPI::SPIVec y6(SPI::set_Cheby_mapped_y(-1.,10.,201) ,"yCheby");
1814  SPI::SPIgrid1D grid6(y6,"grid",SPI::Chebyshev);
1815  a2 = (4.0*(y6^2)+(y6)+10.0);
1816  test_if_close(integrate(a2,grid6),8965.0/6.0,"integrate 2",1e-1);
1817 
1818  SPI::SPIVec y7(SPI::set_Cheby_mapped_y(0.,21.,41) ,"yCheby");
1819  SPI::SPIgrid1D grid7(y7,"grid",SPI::UltraS);
1820  SPI::SPIVec a7(41*4,"a3");
1821  //SPI::SPIVec tmp(grid7.That*(1.0-(y7^2)));
1822  SPI::SPIVec tmp(grid7.That*(1.0-(y7^2)+y7*y7*y7+y7*PETSC_i));
1823  for(PetscInt i=0; i<4; ++i){
1824  for(PetscInt j=0; j<41; ++j){
1825  a7(i*41+j,tmp(j,PETSC_TRUE));
1826  }
1827  }
1828  a7();
1829  //test_if_close(integrate(a7,grid7),-3066.0*4.0,"integrate 3",1e-11);
1830  test_if_close(integrate(a7,grid7),182217.0,"integrate 3",1e-11);
1831  test_if_close(PetscImaginaryPart(integrate(a7,grid7)),441.0*2.0,"integrate 3 imaginary",1e-11);
1832  SPI::printf("------------ UltraS int end -----------");
1833  }
1834  if(alltests){ // and timing of LSTNP_spatials_right vs LSTNP_spatial_right
1835  SPI::printf("------------ LSTNP_spatials_right Blasius boundary layer UltraS start -----------");
1836  PetscInt n=119;
1837  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,21.,n) ,"yCheby");
1838  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
1839  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
1840  //grid.print();
1841  //exit(0);
1842  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
1843  //SPI::printf("n = %d",n);
1844  //grid.print();
1845  //(grid.Dy*grid.y).print();
1846 
1847  SPI::SPIparams params("Blasius parameters");
1848  //params.Re = 400.0;
1849  //params.omega = 86.*params.Re/(1000000.);
1850  params.Re = 1000.0/1.7208;
1851  params.omega = 0.26/1.7208;
1852  params.nu = 1./params.Re;
1853  params.x = params.Re;
1854  params.x_start = params.x;
1855  params.alpha = (0.74155+0.345132*PETSC_i)/1.7208;
1856  params.beta = 0.0;
1857  //params.print();
1858 
1859  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
1860  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
1861  //SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
1862  PetscScalar eigenvalue;
1863  eigenfunction.name = "eigenfunction";
1864 
1865  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
1866  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
1867  if(1){ // set to parallel baseflow
1868  SPI::SPIVec o(SPI::zeros(n),"zero");
1869  bl_flow.Ux = o;
1870  bl_flow.Uxy = o;
1871  bl_flow.V = o;
1872  bl_flow.Vy = o;
1873  }
1874  //bl_flow.print();
1875  // timing
1876  time_t timer1,timer2,timer3;
1877  double seconds;
1878  time(&timer1); // get current time
1879 
1880  std::vector<PetscScalar> alphas_guess(3);
1881  std::vector<SPI::SPIVec> eigenfunction_guess(3);
1882  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1883  alphas_guess[0] = eigenvalue;
1884  eigenfunction_guess[0] = eigenfunction;
1885  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
1886  test_if_close(params.alpha*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatial_right 1",1e-4);
1887 
1888  params.alpha = (0.54213+0.083968*PETSC_i)/1.7208;
1889  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1890  alphas_guess[1] = eigenvalue;
1891  eigenfunction_guess[1] = eigenfunction;
1892  test_if_close(params.alpha*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatial_right 2",1e-5);
1893 
1894  params.alpha = (0.29967+0.230773*PETSC_i)/1.7208;
1895  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
1896  alphas_guess[2] = eigenvalue;
1897  eigenfunction_guess[2] = eigenfunction;
1898  test_if_close(params.alpha*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatial_right 3",1e-3);
1899  time(&timer2); // get current time
1900 
1901  std::vector<PetscScalar> alpha_spatials_right;
1902  std::vector<SPI::SPIVec> eig_vecs_spatials_right;
1903  std::tie(alpha_spatials_right,eig_vecs_spatials_right) = SPI::LSTNP_spatials_right(params,grid,bl_flow,alphas_guess,eigenfunction_guess);
1904  test_if_close(alpha_spatials_right[0]*1.7208,(0.74155+0.345132*PETSC_i),"LSTNP_spatials_right 1",1e-4);
1905  test_if_close(alpha_spatials_right[1]*1.7208,(0.54213+0.083968*PETSC_i),"LSTNP_spatials_right 2",1e-5);
1906  test_if_close(alpha_spatials_right[2]*1.7208,(0.29967+0.230773*PETSC_i),"LSTNP_spatials_right 3",1e-3);
1907  //std::cout<<"alpha_spatials_right[0] = "<<alpha_spatials_right[0]<<std::endl;
1908  //std::cout<<"alpha_spatials_right[1] = "<<alpha_spatials_right[1]<<std::endl;
1909  //std::cout<<"alpha_spatials_right[2] = "<<alpha_spatials_right[2]<<std::endl;
1910  //std::cout<<"alphas_guess[0] = "<<alphas_guess[0]<<std::endl;
1911  //std::cout<<"alphas_guess[1] = "<<alphas_guess[1]<<std::endl;
1912  //std::cout<<"alphas_guess[2] = "<<alphas_guess[2]<<std::endl;
1913  time(&timer3); // get current time
1914  seconds = difftime(timer2,timer1);
1915  SPI::printf("%.10f seconds for LSTNP_spatial_right 3 solves on Blasius boundary layer",seconds);
1916  seconds = difftime(timer3,timer2);
1917  SPI::printf("%.10f seconds for LSTNP_spatials_right 3 solves on Blasius boundary layer",seconds);
1918 
1919  SPI::printf("------------ LSTNP_spatials_right Blasius boundary layer UltraS end -----------");
1920  }
1921  if(alltests){
1922  SPI::printf("------------ A*x=b grid UltraS start -----------");
1923  PetscInt n=20;
1924  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,1.,n) ,"yCheby");
1925  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev);
1926  //SPI::SPIgrid1D grid(y,"grid",SPI::FD);
1927  SPI::SPIgrid1D grid2(y,"grid",SPI::UltraS);
1928  SPI::SPIMat A;
1929  A = grid.Dyy + grid.Dy + grid.I;
1930  SPI::SPIVec b(SPI::zeros(n));
1931  std::vector<PetscInt> rowBCs = {0*n,n-1};
1932  A.eye_rows(rowBCs);
1933  b(rowBCs[0],0.0);
1934  b(rowBCs[1],1.0);
1935  b();
1936  SPI::SPIVec x;
1937  x = SPI::solve(A,b);
1938  //x.print();
1939 
1940  // exact solution
1941  SPI::SPIVec xexact(SPI::zeros(n));
1942  PetscScalar i=PETSC_i;
1943  PetscScalar a=std::exp(0.5)/(2.0*i*std::sin(std::sqrt(3.0)/2.0));
1944  //std::cout<<"a = "<<a<<std::endl;
1945  xexact = SPI::exp(-0.5*grid.y)*((a*SPI::exp((i*std::sqrt(3.0)/2.0)*grid.y)) - (a*SPI::exp((-i*std::sqrt(3.0)/2.0)*grid.y )));
1946  //(x-xexact).print();
1947  std::cout<<"error physical = "<<SPI::L2(x-xexact)<<std::endl;
1948 
1949  // UltraS solution
1950  SPI::SPIMat A2(n,n);
1951  //MatMPIAIJSetPreallocation(A.mat,20,NULL,20,NULL);
1952  A2 = grid2.Dyy + grid2.Dy + grid2.I;
1953  //MatMPIAIJSetPreallocation(A.mat,20,NULL,20,NULL);
1954  b = SPI::zeros(n);
1955  rowBCs = {n-2,n-1};
1956  A2.zero_rows(rowBCs);
1957  MatSetOption(A2.mat,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1958  for(PetscInt j=0; j<n; ++j){
1959  A2(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
1960  A2(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the wall
1961  }
1962  b(rowBCs[0],0.0);
1963  b(rowBCs[1],1.0);
1964  b();
1965  A2();
1966  x = SPI::solve(A2,b);
1967  std::cout<<"error UltraS = "<<SPI::L2((grid.T*x)-xexact)<<std::endl;
1968  xexact = (SPI::exp(0.5-grid.y/2.0)/std::sin(std::sqrt(3.0)/2.0))*SPI::sin(std::sqrt(3.0)/2.0 * grid.y);
1969  std::cout<<"error UltraS = "<<SPI::L2((grid.T*x)-xexact)<<std::endl;
1970  //std::cout<<"erf(1) = "<<std::erf(1.0)<<std::endl;
1971  //std::cout<<"erf(1) = "<<std::erf((double)1.0)<<std::endl;
1972  //std::cout<<"erf(1) = "<<std::erf(std::complex<double>(1.0))<<std::endl;
1973  //SPI::erf(SPI::ones(n)).print();
1974  //
1975  // solve u'' + x*u' + u = (4*y + 17)*e^(4*x) with u(0) = 1 and u(1) = e^4
1976  // solution is u=e^(4*y)
1977  A = grid.Dyy + SPI::diag(y)*grid.Dy + grid.I;
1978  b = (4.0*y + 17.0) * SPI::exp(4.0*y);
1979  rowBCs = {0*n,n-1};
1980  A.eye_rows(rowBCs);
1981  b(rowBCs[0],1.0);
1982  b(rowBCs[1],std::exp(4.0));
1983  b();
1984  x = SPI::solve(A,b);
1985 
1986  // exact solution
1987  xexact = SPI::exp(4.0*y);
1988 
1989  std::cout<<"error non-const = "<<SPI::L2(x - xexact)<<std::endl;
1990 
1991  // solve non-const using UltraS
1992  A = grid2.Dyy + grid2.S1S0That*(SPI::diag(y)*(grid2.T*grid2.S0invS1inv*grid2.Dy)) + grid2.I;
1993  b = grid2.S1S0That*(((4.0*y) + 17.0) * SPI::exp(4.0*y));
1994  rowBCs = {n-2,n-1};
1995  A.zero_rows(rowBCs);
1996  MatSetOption(A.mat,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1997  for(PetscInt j=0; j<n; ++j){
1998  A(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
1999  A(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the wall
2000  }
2001  A();
2002  b(rowBCs[0],1.0);
2003  b(rowBCs[1],std::exp(4.0));
2004  b();
2005  x = SPI::solve(A,b);
2006 
2007  std::cout<<"error non-const UltraS = "<<SPI::L2(grid.T*x - xexact)<<std::endl;
2008 
2009  SPI::printf("------------ A*x=b grid UltraS end -----------");
2010  }
2011  if(alltests){ // and timing of LSTNP_spatials_right vs LSTNP_spatial_right
2012  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS start -----------");
2013  PetscInt n=169;
2014  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,61.,n) ,"yCheby");
2015  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
2016  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
2017  //grid.print();
2018  //exit(0);
2019  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
2020  //SPI::printf("n = %d",n);
2021  //grid.print();
2022  //(grid.Dy*grid.y).print();
2023 
2024  SPI::SPIparams params("Blasius parameters");
2025  params.Re = 400.0;
2026  params.omega = 86.*params.Re/(1000000.);
2027  params.nu = 1./params.Re;
2028  params.x = params.Re;
2029  params.x_start = params.x;
2030  params.alpha = (0.094966+0.004564*PETSC_i);
2031  params.beta = 0.0;
2032  //params.print();
2033 
2034  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
2035  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
2036  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
2037  PetscScalar eigenvalue,cg;
2038  eigenfunction.name = "eigenfunction";
2039 
2040  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
2041  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
2042  if(0){ // set to parallel baseflow
2043  SPI::SPIVec o(SPI::zeros(n),"zero");
2044  bl_flow.Ux = o;
2045  bl_flow.Uxy = o;
2046  bl_flow.V = o;
2047  bl_flow.Vy = o;
2048  }
2049  //bl_flow.print();
2050  // timing
2051  time_t timer1,timer2;
2052  double seconds;
2053 
2054  time(&timer1); // get current time
2055  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
2056  time(&timer2); // get current time
2057  seconds = difftime(timer2,timer1);
2058  SPI::printf("%.10f seconds for LSTNP_spatial_right 1 solves on Blasius boundary layer",seconds);
2059  test_if_close(params.alpha,(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatial_right 1",1e-8);
2060 
2061  std::vector<PetscScalar> alphas_guess(2);
2062  std::vector<SPI::SPIVec> eigenfunction_guess(2);
2063  time(&timer1); // get current time
2064  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
2065  time(&timer2); // get current time
2066  seconds = difftime(timer2,timer1);
2067  SPI::printf("%.10f seconds for LSTNP_spatial_right 1 solves on Blasius boundary layer",seconds);
2068  alphas_guess[0] = eigenvalue;
2069  eigenfunction_guess[0] = eigenfunction;
2070  //std::tie(eigenvalue,eigenfunction) = SPI::LST_spatial(params,grid,bl_flow);
2071  test_if_close(params.alpha,(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatial_right 1",1e-8);
2072  //SPI::printf("LSTNP_spatial_right 1 non-Parallel Blasius boundary layer eigenvalue = %.10f + %.10fi",eigenvalue);
2073 
2074  time(&timer1); // get current time
2075  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
2076  time(&timer2); // get current time
2077  seconds = difftime(timer2,timer1);
2078  SPI::printf("%.10f seconds for LSTNP_spatial 2 solves on Blasius boundary layer",seconds);
2079  test_if_close(params.alpha,(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatial 2",1e-8);
2080 
2081  params.alpha = (0.10665444+0.00189793*PETSC_i);
2082  time(&timer1); // get current time
2083  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
2084  time(&timer2); // get current time
2085  seconds = difftime(timer2,timer1);
2086  SPI::printf("%.10f seconds for LSTNP_spatial_right 3 solves on Blasius boundary layer",seconds);
2087  alphas_guess[1] = eigenvalue;
2088  eigenfunction_guess[1] = eigenfunction;
2089  test_if_close(params.alpha,(0.106654447306241+0.001897936897103*PETSC_i),"LSTNP_spatial_right 3",1e-8);
2090  //SPI::printf("LSTNP_spatial_right 2 non-Parallel Blasius boundary layer eigenvalue = %.10f + %.10fi",eigenvalue);
2091 
2092  time(&timer1); // get current time
2093  std::tie(eigenvalue,cg,leigenfunction,eigenfunction) = SPI::LSTNP_spatial(params,grid,bl_flow);
2094  time(&timer2); // get current time
2095  seconds = difftime(timer2,timer1);
2096  SPI::printf("%.10f seconds for LSTNP_spatial 4 solves on Blasius boundary layer",seconds);
2097  test_if_close(params.alpha,(0.106654447306241+0.001897936897103*PETSC_i),"LSTNP_spatial 4",1e-8);
2098 
2099  std::vector<PetscScalar> eigenvalues(2);
2100  std::vector<SPI::SPIVec> eigenfunctions(2);
2101  time(&timer1); // get current time
2102  std::tie(eigenvalues,eigenfunctions) = SPI::LSTNP_spatials_right(params,grid,bl_flow,alphas_guess,eigenfunction_guess);
2103  time(&timer2); // get current time
2104  seconds = difftime(timer2,timer1);
2105  SPI::printf("%.10f seconds for LSTNP_spatials_right 5,6 solves on Blasius boundary layer",seconds);
2106  test_if_close(eigenvalues[0],(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatials_right 5",1e-8);
2107  test_if_close(eigenvalues[1],(0.106654447306241+0.001897936897103*PETSC_i),"LSTNP_spatials_right 6",1e-8);
2108 
2109  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS end -----------");
2110  }
2111  if(alltests){ // and timing of LSTNP_spatials_right vs LSTNP_spatial_right
2112  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS start -----------");
2113  PetscInt n=169;
2114  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,61.,n) ,"yCheby");
2115  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
2116  SPI::SPIgrid1D grid(y,"grid",SPI::UltraS);
2117  //grid.print();
2118  //exit(0);
2119  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
2120  //SPI::printf("n = %d",n);
2121  //grid.print();
2122  //(grid.Dy*grid.y).print();
2123 
2124  SPI::SPIparams params("Blasius parameters");
2125  params.Re = 400.0;
2126  params.omega = 86.*params.Re/(1000000.);
2127  params.nu = 1./params.Re;
2128  params.x = params.Re;
2129  params.x_start = params.x;
2130  params.alpha = (0.094966+0.004564*PETSC_i);
2131  params.beta = 0.0;
2132  //params.print();
2133 
2134  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
2135  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
2136  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
2137  PetscScalar eigenvalue,cg;
2138  eigenfunction.name = "eigenfunction";
2139 
2140  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
2141  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
2142  if(0){ // set to parallel baseflow
2143  SPI::SPIVec o(SPI::zeros(n),"zero");
2144  bl_flow.Ux = o;
2145  bl_flow.Uxy = o;
2146  bl_flow.V = o;
2147  bl_flow.Vy = o;
2148  }
2149  //bl_flow.print();
2150  // timing
2151  time_t timer1,timer2;
2152  double seconds;
2153 
2154  time(&timer1); // get current time
2155  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right2(params,grid,bl_flow);
2156  time(&timer2); // get current time
2157  seconds = difftime(timer2,timer1);
2158  SPI::printf("%.10f seconds for LSTNP_spatial_right 1 solves on Blasius boundary layer",seconds);
2159  test_if_close(params.alpha,(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatial_right 1",1e-8);
2160  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS end -----------");
2161  }
2162  if(0){ // and timing of LSTNP_spatials_right vs LSTNP_spatial_right
2163  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS start -----------");
2164  PetscInt n=169;
2165  SPI::SPIVec y(SPI::set_Cheby_mapped_y(0.,61.,n) ,"yCheby");
2166  //SPI::SPIVec y(SPI::set_Cheby_mapped_y(-1.,1.,n) ,"yCheby");
2167  //SPI::SPIgrid1D grid0(y,"grid",SPI::UltraS); // good now
2168  //SPI::SPIgrid1D grid2(y,"grid",SPI::Chebyshev); // good now
2169  //SPI::SPIgrid1D grid3(y,"grid",SPI::FD); // good now
2170  //SPI::SPIVec y2;
2171  //2.0*y;
2172  //2.0+y;
2173  //2.0/y;
2174  //y/2.0;
2175  //y2 = y;
2176  // test grid making routines
2177  //SPI::SPIMat Dy_1(SPI::set_D(y,1)); // this is good now
2178  //SPI::SPIMat Dy_2(SPI::set_D_Chebyshev(y,1,PETSC_TRUE)); // this is good now
2179  //SPI::SPIMat Dyy_1(SPI::set_D(y,2)); // this is good now
2180  //SPI::SPIMat Dyy_2(SPI::set_D_Chebyshev(y,2,PETSC_TRUE)); // this is good now
2181  //SPI::SPIVec xi = (SPI::ones(n));// good
2182  //SPI::SPIMat I(SPI::eye(n),"I");// good
2183  //SPI::SPIVec s(SPI::arange(4));// good
2184  //SPI::SPIMat D(n);// good
2185  //D = SPI::zeros(n,n);
2186  //SPI::SPIVec Coeffs(SPI::get_D_Coeffs(s,1)); // this is good now
2187  //Coeffs.~SPIVec();
2188  //Coeffs = SPI::arange(4);
2189  //Coeffs.~SPIVec();
2190  //Coeffs = SPI::arange(4);
2191  //D*=4.0;
2192  //SPI::map_D(D,y,1,4); // good now
2193  //SPI::SPIMat Dy(SPI::map_D(D,y,1,4)); // good now
2194  //SPI::SPIMat D1(SPI::eye(n)); // good
2195  //SPI::diag(1./(D*y)); // good
2196  //SPI::diag(1./(D*y))*D; // good now;
2197  //D*D; // good now
2198  //y^3;
2199  //y^y;
2200  //D*y;
2201  //D+D;
2202  //SPI::SPIMat Dy1(SPI::set_D(y,1,4,PETSC_TRUE));
2203  //-1.*(D*y);
2204  //SPI::diag(y^2);
2205  //s *= y; // good
2206  //SPI::factorial(3);
2207  //SPI::SPIMat A(3);
2208  //A(0,0,1.0); A(0,1,1.0); A(0,2,1.0);
2209  //A(1,0,1.0); A(1,1,2.0); A(1,2,3.0);
2210  //A(2,0,1.0); A(2,1,4.0); A(2,2,9.0);
2211  //SPI::SPIVec b(3);
2212  //b(1,1.0);
2213  //A();b();
2214  //SPI::solve(A,b); // good now
2215  //b/A; // good now
2216  //SPI::block({
2217  //{A,A},
2218  //{A,A}
2219  //})();
2220  //SPI::set_T_That(y.rows); // good now
2221  //SPI::meshgrid(y,y); // good now
2222  //D%D; // good now
2223  //SPI::SPIMat A2(4,3); // good
2224  //A2();
2225  //SPI::SPIMat A3(A2); // good
2226  //SPI::SPIMat Y1,Y2;
2227  //std::tie(Y1,Y2) = SPI::meshgrid(y,y); // good
2228  //cos(Y1%acos(Y2));
2229  //SPI::SPIMat Y3;
2230  //Y1.T(Y3);
2231  //SPI::SPIgrid1D grid(y,"grid",SPI::UltraS); // works
2232  SPI::SPIgrid1D grid(y,"grid",SPI::Chebyshev); // works
2233  //grid.print();
2234  //exit(0);
2235  //SPI::SPIVec y(SPI::linspace(0.,61.,n) ,"yFD");
2236  //SPI::printf("n = %d",n);
2237  //grid.print();
2238  //(grid.Dy*grid.y).print();
2239 
2240  SPI::SPIparams params("Blasius parameters");
2241  params.Re = 400.0;
2242  params.omega = 86.*params.Re/(1000000.);
2243  params.nu = 1./params.Re;
2244  params.x = params.Re;
2245  params.x_start = params.x;
2246  params.alpha = (0.094966+0.004564*PETSC_i);
2247  params.beta = 0.0;
2248  //params.print();
2249 
2250  SPI::SPIVec eigenfunction(grid.y.rows*16,"eigenfunction");
2251  SPI::SPIVec eigenfunction2(grid.y.rows*8,"eigenfunction");
2252  //SPI::SPIVec eig_vec(grid.y.rows*8,"q");
2253  SPI::SPIVec leigenfunction(grid.y.rows*16,"q");
2254  PetscScalar eigenvalue,eigenvalue2,cg;
2255  eigenfunction.name = "eigenfunction";
2256 
2257  SPI::SPIbaseflow bl_flow = SPI::blasius(params,grid);
2258  //SPI::SPIbaseflow bl_flow = SPI::channel(params,grid);
2259  if(0){ // set to parallel baseflow
2260  SPI::SPIVec o(SPI::zeros(n),"zero");
2261  bl_flow.Ux = o;
2262  bl_flow.Uxy = o;
2263  bl_flow.V = o;
2264  bl_flow.Vy = o;
2265  }
2266  //bl_flow.print();
2267  // timing
2268  time_t timer1,timer2;
2269  double seconds;
2270 
2271  time(&timer1); // get current time
2272  params.print();
2273  //grid.print();
2274  //(grid.Dy*grid.y).print();
2275  //(grid.Dyy*grid.y).print();
2276  //bl_flow.print();
2277  //grid.T.print();
2278  //grid.That.print();
2279  //std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right2(params,grid,bl_flow); // good
2280  std::tie(eigenvalue,eigenfunction) = SPI::LSTNP_spatial_right(params,grid,bl_flow);
2281  std::tie(eigenvalue2,eigenfunction2) = SPI::LST_spatial(params,grid,bl_flow);
2282  time(&timer2); // get current time
2283  seconds = difftime(timer2,timer1);
2284  SPI::printf("%.10f seconds for LSTNP_spatial_right 1 solves on Blasius boundary layer",seconds);
2285  test_if_close(eigenvalue,(0.094966355495876+0.004564261943353*PETSC_i),"LSTNP_spatial_right 1",1e-8);
2286  SPI::printf("------------ LSTNP_spatials_right non-Parallel Blasius boundary layer UltraS end -----------");
2287 
2288  SPI::SPIMat M(4),C(4),K(4);
2289  M(0,0,3.0); M(1,1,1.0); M(2,2,3.0); M(3,3,1.0); M();
2290  C(0,0,0.4); C(0,2,-0.3);
2291  C(2,0,-0.3); C(2,2,0.5); C(2,3,-0.2);
2292  C(3,2,-0.2); C(3,3,0.2);
2293  C();
2294  K(0,0,-7.0); K(0,1, 2.0); K(0,2, 4.0);
2295  K(1,0, 2.0); K(1,1,-4.0); K(1,2, 2.0);
2296  K(2,0, 4.0); K(2,1, 2.0); K(2,2,-9.0); K(2,3, 3.0);
2297  K(3,2, 3.0); K(3,3,-3.0);
2298  K();
2299  SPI::SPIVec eigenfunction_4(4);
2300  std::tie(eigenvalue,eigenfunction_4) = SPI::polyeig({K,C,M},-2.4498); // good
2301  SPI::printfc("eigenvalue = %.10f + %.10f",eigenvalue);
2302  K.zero_rows({1,2});
2303  M.eye_row(1);
2304  K = M*K;
2305  //eigenfunction_4.print();
2306  //(0.1828*eigenfunction_4/eigenfunction_4(0,PETSC_TRUE)).print();
2307  //(eigenvalue*eigenvalue*(M*eigenfunction_4) + eigenvalue*(C*eigenfunction_4) + (K*eigenfunction_4)).print();
2308  }
2309  if(alltests){
2310  SPI::printf("------------ SVD start -----------");
2311  SPI::SPIVec x(4,"x");
2312  x(0,0.0); x(1,1.0); x(2,2.0); x(3,3.0);
2313  x();
2314  SPI::SPIVec y(4,"y");
2315  y(0,-1.0); y(1,0.2); y(2,0.9); y(3,2.1);
2316  y();
2317  //y.print();
2318  SPI::SPIMat A(4,2,"A");
2319  A(0,0,0); A(0,1,1.0);
2320  A(1,0,1); A(1,1,1.0);
2321  A(2,0,2); A(2,1,1.0);
2322  A(3,0,3); A(3,1,1.0);
2323  A();
2324 
2325  std::vector<PetscReal> sigma;
2326  std::vector<SPI::SPIVec> u;
2327  std::vector<SPI::SPIVec> v;
2328 
2329  //std::cout<<"made it here before svd"<<std::endl;
2330  //A.print();
2331  std::tie(sigma,u,v) = SPI::svd(A);
2332  //SPI::svd(A);
2333  //std::cout<<"made it here after svd"<<std::endl;
2334  SPI::SPIVec xi(SPI::zeros(2));
2335  for(PetscInt j=0; j<2; ++j){
2336  xi += ((y.dot(u[j]))/sigma[j]) * v[j];
2337  //std::cout<<"sigma["<<j<<"] = "<<sigma[j]<<std::endl;
2338  //u[j].print();
2339  //v[j].print();
2340  }
2341  //xi.print();
2342  test_if_close(xi(0,PETSC_TRUE),1.0,"SVD 1",1e-14);
2343  test_if_close(xi(1,PETSC_TRUE),-0.95,"SVD 2",1e-14);
2344 
2345  //(A*v[0] - sigma[0]*u[0]).print();
2346  //(A*v[1] - sigma[1]*u[1]).print();
2347  //(A*xi).print();
2348 
2349  SPI::printf("------------ SVD end -----------");
2350  }
2351  if(alltests){
2352  SPI::printf("------------ least squares start -----------");
2353  SPI::SPIVec x;
2354  //SPI::SPIVec x(4,"x");
2355  //x(0,0.0); x(1,1.0); x(2,2.0); x(3,3.0);
2356  //x();
2357  SPI::SPIVec y(4,"y");
2358  y(0,-1.0); y(1,0.2); y(2,0.9); y(3,2.1);
2359  y();
2360  SPI::SPIMat A(4,2,"A");
2361  A(0,0,0); A(0,1,1.0);
2362  A(1,0,1); A(1,1,1.0);
2363  A(2,0,2); A(2,1,1.0);
2364  A(3,0,3); A(3,1,1.0);
2365  A();
2366 
2367  x = lstsq(A,y);
2368  test_if_close(x(0,PETSC_TRUE),1.0,"lstsq 1",1e-14);
2369  test_if_close(x(1,PETSC_TRUE),-0.95,"lstsq 2",1e-14);
2370  //x.print();
2371 
2372  SPI::printf("------------ least squares end -----------");
2373  }
2374  if(alltests){
2375  SPI::printf("------------ set_col start -----------");
2376  SPI::SPIVec x1(4,"x1");
2377  x1(0,0.0); x1(1,1.0); x1(2,2.0); x1(3,3.0);
2378  x1();
2379  //x1.print();
2380  SPI::SPIVec x2(4,"x2");
2381  x2(0,1.0); x2(1,1.0); x2(2,1.0); x2(3,1.0);
2382  x2();
2383  //x2.print();
2384  SPI::SPIMat A(4,2,"A");
2385  A.set_col(0,x1);
2386  A.set_col(1,x2);
2387  A();
2388  //A.print();
2389  test_if_close(A(2,0,PETSC_TRUE),2.0,"set_col 1",1e-14);
2390  test_if_close(A(2,1,PETSC_TRUE),1.0,"set_col 2",1e-14);
2391  SPI::printf("------------ set_col end -----------");
2392  }
2393  if(alltests){
2394  SPI::printf("------------ SPIMat(std::vector<SPIVec>) and lstsq(std::vector<SPIVec>,SPIVec) start -----------");
2395  SPI::SPIVec y(4,"y");
2396  y(0,-1.0); y(1,0.2); y(2,0.9); y(3,2.1);
2397  y();
2398  SPI::SPIVec x1(4,"x1");
2399  x1(0,0.0); x1(1,1.0); x1(2,2.0); x1(3,3.0);
2400  x1();
2401  SPI::SPIVec x2(4,"x2");
2402  x2(0,1.0); x2(1,1.0); x2(2,1.0); x2(3,1.0);
2403  x2();
2404  std::vector<SPI::SPIVec> x = {x1,x2};
2405  SPI::SPIMat A(x);
2406  //A.print();
2407  test_if_close(A(2,0,PETSC_TRUE),2.0,"set_col 1",1e-14);
2408  test_if_close(A(2,1,PETSC_TRUE),1.0,"set_col 2",1e-14);
2409 
2410  SPI::SPIVec xi(SPI::lstsq(x,y));
2411  test_if_close(xi(0,PETSC_TRUE),1.0,"lstsq(std::vector<SPIVec>,SPIVec) 1",1e-14);
2412  test_if_close(xi(1,PETSC_TRUE),-0.95,"lstsq(std::vector<SPIVec>,SPIVec) 2",1e-14);
2413  SPI::printf("------------ SPIMat(std::vector<SPIVec>) and lstsq(std::vector<SPIVec>,SPIVec) end -----------");
2414  }
2415  if(alltests){
2416  SPI::printf("------------ SPIMat save start -----------");
2417  SPI::SPIMat A(4,2,"A");
2418  A(0,0,0.2+3.1*PETSC_i); A(0,1,1.0+4.0*PETSC_i);
2419  A(1,0,1.0+4.2*PETSC_i); A(1,1,2.0+3.0*PETSC_i);
2420  A(2,0,2.0+3.3*PETSC_i); A(2,1,3.0+2.0*PETSC_i);
2421  A(3,0,3.0+4.4*PETSC_i); A(3,1,4.0+1.0*PETSC_i);
2422  A();
2423  A.print();
2424  SPI::save(A,"A.dat");
2425  //SPI::load(A,"A.dat");
2426  //A.print();
2427  SPI::printf("------------ SPIMat save end -----------");
2428  }
2429  if(alltests){
2430  SPI::printf("------------ SPIgrid2D.avgt start -----------");
2431  PetscInt ny=400;
2432  PetscInt nt=4;
2433  SPI::SPIVec y(SPI::set_FD_stretched_y(61.,ny,1.01) ,"yFD");
2434  SPI::SPIVec t(SPI::set_Fourier_t((2.0*M_PI)/0.0344,nt) ,"t");
2435  SPI::SPIgrid2D grid(y,t,"grid",SPI::FD,SPI::Fourier);
2436  //(grid.avgt*SPI::ones(400*4)).print();
2437  test_if_close(grid.avgt(1,1,PETSC_TRUE),0.25,"SPIgrid2D.avgt 1",1e-14);
2438  //t.print();
2439  SPI::printf("------------ SPIgrid2D.avgt end -----------");
2440  }
2441  if(alltests){
2442  SPI::printf("------------ SPIgrid interp1D_Mat start -----------");
2443  PetscInt n1=8,n2=11;
2444  SPI::SPIVec y1(SPI::linspace(0,4,n1),"y1");
2445  SPI::SPIVec y2(SPI::linspace(0,2,n2),"y2");
2446  SPI::SPIgrid1D grid1(y1,"grid1",SPI::FD);
2447  SPI::SPIgrid1D grid2(y2,"grid2",SPI::FD);
2448  //grid1.print();
2449  //(grid1.Dyy*(y1*y1)).print();
2450  SPI::SPIMat interp(SPI::interp1D_Mat(grid1,grid2));
2451  std::cout<<"L2(y2-interp*y1) = "<<SPI::L2(grid2.y,interp*grid1.y)<<std::endl;
2452  std::cout<<"L2(sin(y2)-interp*sin(y1)) = "<<SPI::L2(SPI::sin(grid2.y),interp*SPI::sin(grid1.y))<<std::endl;
2453  interp.print();
2454 
2455  interp = SPI::interp1D_Mat(grid1.y,grid2.y);
2456  interp.print();
2457  SPI::printf("------------ SPIgrid interp1D_Mat end -----------");
2458  }
2459  if(alltests){
2460  SPI::printf("------------ SPIgrid set_D_periodic start -----------");
2461  PetscInt nt=64;
2462  SPI::SPIVec t(SPI::set_Fourier_t(4,nt),"t");
2463  SPI::SPIMat Dt(SPI::set_D_periodic(t,1),"Dt");
2464  (0.5*PETSC_PI*SPI::cos(0.5*PETSC_PI*t)).print();
2465  (Dt*(SPI::sin(0.5*PETSC_PI*t))).print();
2466  SPI::printf("------------ SPIgrid set_D_periodic end -----------");
2467  }
2468  if(1){
2469  SPI::printf("------------ SPIgrid dft and dft_dftinv_Ihalf_Ihalfn start -----------");
2470  PetscInt nt=8;
2471  SPI::SPIMat FT(SPI::dft(nt),"FT");
2472  SPI::SPIMat FTinv(nt,nt,"FTinv");
2473  SPI::SPIMat Ihalf(nt,nt,"Ihalf");
2474  SPI::SPIMat Ihalfn(nt,nt,"Ihalfn");
2475  test_if_close(FT(5,7,PETSC_TRUE),-0.0883883-0.0883883*PETSC_i,"dft 1",1e-7);
2476  test_if_close(FT(6,7,PETSC_TRUE),0.-0.125*PETSC_i,"dft 2",1e-7);
2477  std::tie(FT,FTinv,Ihalf,Ihalfn) = SPI::dft_dftinv_Ihalf_Ihalfn(nt);
2478  test_if_close(FT(5,7,PETSC_TRUE),-0.0883883-0.0883883*PETSC_i,"dft_dftinv_Ihalf_Ihalfn 1",1e-7);
2479  test_if_close(FT(6,7,PETSC_TRUE),0.-0.125*PETSC_i,"dft_dftinv_Ihalf_Ihalfn 2",1e-7);
2480  test_if_close(FTinv(5,7,PETSC_TRUE),-0.707106781+0.0707106781*PETSC_i,"dft_dftinv_Ihalf_Ihalfn 3",1e-7);
2481  test_if_close(FTinv(6,7,PETSC_TRUE),0.+1.*PETSC_i,"dft_dftinv_Ihalf_Ihalfn 4",1e-7);
2482  test_if_close(Ihalf(6,6,PETSC_TRUE),0.,"dft_dftinv_Ihalf_Ihalfn 5",1e-7);
2483  test_if_close(Ihalf(1,1,PETSC_TRUE),1.,"dft_dftinv_Ihalf_Ihalfn 6",1e-7);
2484  test_if_close(Ihalfn(6,6,PETSC_TRUE),1.,"dft_dftinv_Ihalf_Ihalfn 7",1e-7);
2485  test_if_close(Ihalfn(1,1,PETSC_TRUE),0.,"dft_dftinv_Ihalf_Ihalfn 8",1e-7);
2486  SPI::printf("------------ SPIgrid dft and dft_dftinv_Ihalf_Ihalfn end -----------");
2487  }
2488 
2489 
2490  return 0;
2491 }
SPI::sin
SPIMat sin(const SPIMat &A)
take the sin of each element in a matrix
Definition: SPIMat.cpp:2025
SPI::polyeig
std::tuple< PetscScalar, SPIVec > polyeig(const std::vector< SPIMat > &As, const PetscScalar target, const PetscReal tol=-1., const PetscInt max_iter=-1)
solve general polynomial eigenvalue problem of (A0 + A1*alpha + A2*alpha^2 + ...)*x = 0 and return a ...
Definition: SPIMat.cpp:1686
SPI::polyeig_init
std::tuple< PetscScalar, SPIVec > polyeig_init(const std::vector< SPIMat > &As, const PetscScalar target, const SPIVec &qr, const PetscReal tol=-1., const PetscInt max_iter=-1)
solve general polynomial eigenvalue problem of (A0 + A1*alpha + A2*alpha^2 + ...)*x = 0 and return a ...
Definition: SPIMat.cpp:1754
SPI::SPIMat::print
PetscInt print()
print mat to screen using PETSC_VIEWER_STDOUT_WORLD
Definition: SPIMat.cpp:498
SPI::SPIbaseflow::Uxy
SPIVec Uxy
Definition: SPIbaseflow.hpp:44
test_if_true
void test_if_true(PetscBool test, std::string name)
Definition: tests.cpp:6
SPI::SPIparams::Re
PetscScalar Re
Reynolds number.
Definition: SPIparams.hpp:13
SPI::SPIMat::rows
PetscInt rows
number of rows in mat
Definition: SPIMat.hpp:18
SPI::SPIbaseflow::V
SPIVec V
Definition: SPIbaseflow.hpp:40
SPI::SPIgrid1D::Dyy
SPIMat Dyy
2nd derivative operator with respect to y
Definition: SPIgrid.hpp:52
SPI::SPIMat::zero_rows
SPIMat & zero_rows(std::vector< PetscInt > rows)
set rows to zero
Definition: SPIMat.cpp:475
SPI::integrate
PetscScalar integrate(const SPIVec &a, SPIgrid1D &grid)
integrate a vector of chebyshev Coefficients on a physical grid
Definition: SPIgrid.cpp:873
SPI::SPIMat::mat
Mat mat
petsc Mat data
Definition: SPIMat.hpp:28
SPI::eig_right
std::tuple< PetscScalar, SPIVec > eig_right(const SPIMat &A, const SPIMat &B, const PetscScalar target, const PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1072
SPI::SPIparams::alpha
PetscScalar alpha
alpha streamwise wavenumber
Definition: SPIparams.hpp:15
SPI::lstsq
SPIVec lstsq(const SPIMat &A, SPIVec &y)
solve least squares problem of A*x=y for skinny A matrix using SVD
Definition: SPIMat.cpp:905
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
SPImain.hpp
SPI::SPIparams::omega
PetscScalar omega
omega, temporal frequency (rad/s)
Definition: SPIparams.hpp:16
SPI::LST_spatial_cg
std::tuple< PetscScalar, PetscScalar, SPIVec, SPIVec > LST_spatial_cg(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow)
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:349
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::FD
@ FD
finite difference grid
Definition: SPIgrid.hpp:26
SPI::SPIVec
Definition: SPIVec.hpp:13
SPI::printf
PetscInt printf(std::string msg,...)
Definition: SPIprint.cpp:5
SPI::Fourier
@ Fourier
Fourier transform collocated grid.
Definition: SPIgrid.hpp:28
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::SPIbaseflow::Vy
SPIVec Vy
Definition: SPIbaseflow.hpp:47
SPI::LST_temporal
std::tuple< PetscScalar, SPIVec > LST_temporal(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, SPIVec q=SPIVec())
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:5
SPI::SPIgrid2D::y
SPIVec y
grid for wall-normal dimension
Definition: SPIgrid.hpp:83
SPI::ones
SPIVec ones(const PetscInt rows)
Definition: SPIVec.cpp:605
SPI::save
PetscInt save(const SPIMat &A, const std::string filename)
save matrix to filename in binary format (see Petsc documentation for format ) Format is (from Petsc ...
Definition: SPIMat.cpp:1924
SPI::kron
SPIMat kron(const SPIMat &A, const SPIMat &B)
set kronecker inner product of two matrices
Definition: SPIMat.cpp:780
SPI::L2
PetscReal L2(SPIVec x1, const SPIVec x2, NormType type=NORM_2)
calculate the norm of the difference between and vectors.
Definition: SPIVec.cpp:859
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::LSTNP_spatial_right2
std::tuple< PetscScalar, SPIVec > LSTNP_spatial_right2(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, SPIVec qr=SPIVec())
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:1276
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
tests
int tests()
Definition: tests.cpp:24
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::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::SPIMat::diag
SPIVec diag()
get diagonal of matrix
Definition: SPIMat.cpp:445
SPI::SPIgrid1D::That
SPIMat That
Chebyshev operator taking it from physical space to Chebyshev coefficients.
Definition: SPIgrid.hpp:59
SPI::SPIparams::print
virtual PetscInt print()
print all variables in SPIparams
Definition: SPIparams.cpp:23
SPI::SPIparams::x
PetscScalar x
current streamwise position
Definition: SPIparams.hpp:18
SPI::SPIVec::name
std::string name
Vec name (important for hdf5 i/o)
Definition: SPIVec.hpp:26
SPI::blasius
SPIbaseflow blasius(SPIparams &params, SPIgrid1D &grid)
Definition: SPIbaseflow.cpp:90
SPI::eye
SPIMat eye(const PetscInt n)
create, form, and return identity matrix of size n
Definition: SPIMat.cpp:697
SPI::inv
SPIMat inv(const SPIMat &A)
Solve linear system A*Ainv=B using MatMatSolve.
Definition: SPIMat.cpp:603
SPI::LST_spatial
std::tuple< PetscScalar, SPIVec > LST_spatial(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, SPIVec q=SPIVec())
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:147
SPI::orthogonalize
std::vector< SPIVec > orthogonalize(std::vector< SPIVec > &x, SPIgrid1D &grid)
Definition: SPIgrid.cpp:1323
SPI::SPIgrid1D::S0invS1inv
SPIMat S0invS1inv
[in] inverse of S0^-1 * S1^-1
Definition: SPIgrid.hpp:56
SPI::exp
SPIVec exp(const SPIVec &A)
take the exp of each element in a vector
Definition: SPIVec.cpp:712
SPI::SPIparams::nu
PetscScalar nu
kinematic viscosity (typically 1/Re)
Definition: SPIparams.hpp:20
SPI::SPIbaseflow::Ux
SPIVec Ux
Definition: SPIbaseflow.hpp:41
SPI::SPIMat::~SPIMat
~SPIMat()
destructor to delete memory
Definition: SPIMat.cpp:508
SPI::diag
SPIMat diag(const SPIVec &diag, const PetscInt k=0)
set diagonal of matrix
Definition: SPIMat.cpp:728
SPI::printfc
PetscInt printfc(std::string msg, const PetscScalar val)
Definition: SPIprint.cpp:29
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::LSTNP_spatials_right
std::tuple< std::vector< PetscScalar >, std::vector< SPIVec > > LSTNP_spatials_right(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, std::vector< PetscScalar > &alphas, std::vector< SPIVec > &qrs)
solve the local stability theory problem for the linearized Navier-Stokes equations using non-paralle...
Definition: SPILST.cpp:1482
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::channel
SPIbaseflow channel(SPIparams &params, SPIgrid1D &grid)
Definition: SPIbaseflow.cpp:253
SPI::integrate_coeffs
PetscScalar integrate_coeffs(const SPIVec &a)
integrate a vector of chebyshev Coefficients on a grid
Definition: SPIVec.cpp:813
SPI::LSTNP_spatial
std::tuple< PetscScalar, PetscScalar, SPIVec, SPIVec > LSTNP_spatial(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, SPIVec ql=SPIVec(), SPIVec qr=SPIVec())
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:715
SPI::SPIMat::eye_row
SPIMat & eye_row(const PetscInt row)
set a row to zero and set 1 in diagonal entry
Definition: SPIMat.cpp:458
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::SPIVec::axpy
SPIVec & axpy(const PetscScalar a, const SPIVec &X)
Definition: SPIVec.cpp:184
SPI::solve
SPIVec solve(const SPIMat &A, const SPIVec &b)
Solve linear system, Ax=b using solve(A,b) notation.
Definition: SPIMat.cpp:596
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::SPIparams
Definition: SPIparams.hpp:8
SPI::SPIgrid2D
Definition: SPIgrid.hpp:72
SPI::eig_init
std::tuple< PetscScalar, SPIVec, SPIVec > eig_init(const SPIMat &A, const SPIMat &B, const PetscScalar target, const SPIVec &ql, const SPIVec &qr, PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1217
SPI::sqrt
SPIVec sqrt(const SPIVec &A)
take the atanh of each element in a vector
Definition: SPIVec.cpp:741
test_if_close
void test_if_close(PetscScalar value, PetscScalar golden, std::string name, PetscReal tol)
Definition: tests.cpp:11
SPI::SPIgrid1D::S1S0That
SPIMat S1S0That
UltraSpherical helper matrix S1*S0*That for baseflow.
Definition: SPIgrid.hpp:55
SPI::load
PetscInt load(SPIMat &A, const std::string filename)
load matrix from filename from binary format (works with save(SPIMat,std::string) function
Definition: SPIMat.cpp:1966
SPI::SPIMat
Definition: SPIMat.hpp:17
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::SPIMat::conj
SPIMat & conj()
elemenwise conjugate current matrix
Definition: SPIMat.cpp:435
SPI::SPIparams::beta
PetscScalar beta
beta spanwise wavenumber
Definition: SPIparams.hpp:14
tests.hpp
SPI::SPIVec::conj
SPIVec & conj()
Definition: SPIVec.cpp:419
SPI::SPIVec::rows
PetscInt rows
number of rows in vec
Definition: SPIVec.hpp:14
SPI::LSTNP_spatial_right
std::tuple< PetscScalar, SPIVec > LSTNP_spatial_right(SPIparams &params, SPIgrid1D &grid, SPIbaseflow &baseflow, SPIVec qr=SPIVec())
solve the local stability theory problem for the linearized Navier-Stokes equations using parallel ba...
Definition: SPILST.cpp:1046
SPI::proj
SPIVec proj(SPIVec &u, SPIVec &v, SPIgrid1D &grid)
Definition: SPIgrid.cpp:1220
SPI::SPIparams::x_start
PetscScalar x_start
streamwise starting location
Definition: SPIparams.hpp:17
SPI::SPIVec::dot
PetscScalar dot(SPIVec y)
take the inner product of two vectors
Definition: SPIVec.cpp:459
SPI::SPIgrid2D::grid1Dy
SPIgrid1D grid1Dy
Definition: SPIgrid.hpp:105
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::eig_init_right
std::tuple< PetscScalar, SPIVec > eig_init_right(const SPIMat &A, const SPIMat &B, const PetscScalar target, const SPIVec &qr, PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1368
SPI::eig
std::tuple< PetscScalar, SPIVec, SPIVec > eig(const SPIMat &A, const SPIMat &B, const PetscScalar target, const PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:928
SPI::cos
SPIMat cos(const SPIMat &A)
take the cos of each element in a matrix
Definition: SPIMat.cpp:2027
SPI::SPIgrid1D::I
SPIMat I
identity matrix same size as derivative operators
Definition: SPIgrid.hpp:65
SPI::SPIVec::vec
Vec vec
petsc Vec data
Definition: SPIVec.hpp:21
SPI::SPIbaseflow
Definition: SPIbaseflow.hpp:17
SPI::SPIMat::set_col
SPIMat & set_col(const PetscInt col, const SPIVec &v)
Definition: SPIMat.cpp:86
SPI::SPIMat::T
PetscInt T(SPIMat &A)
A = Transpose(*this.mat) operation with initialization of A.
Definition: SPIMat.cpp:392
SPI::SPIVec::max
PetscScalar max()
Definition: SPIVec.cpp:425
SPI::svd
std::tuple< std::vector< PetscReal >, std::vector< SPIVec >, std::vector< SPIVec > > svd(const SPIMat &A)
solve SVD problem of A=U*E*V^H for skinny A matrix
Definition: SPIMat.cpp:851
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::SPIMat::eye_rows
SPIMat & eye_rows(std::vector< PetscInt > rows)
set rows to zero and set main diagonal to 1
Definition: SPIMat.cpp:482
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