SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
SPILST.cpp
Go to the documentation of this file.
1 #include "SPILST.hpp"
2 
3 namespace SPI{
5  std::tuple<PetscScalar, SPIVec> LST_temporal(
6  SPIparams &params,
7  SPIgrid1D &grid,
8  SPIbaseflow &baseflow,
9  SPIVec q
10  ){
11  PetscInt n = grid.ny;
12  PetscScalar Re = params.Re;
13  PetscScalar alpha = params.alpha;
14  PetscScalar omega = params.omega;
15  PetscScalar beta = params.beta;
16  PetscScalar i=PETSC_i;
17  PetscScalar k=alpha*alpha+beta*beta;
18  PetscScalar k2=k*k;
19  const SPIMat &O = grid.O;
20  const SPIMat &I = grid.I;
21  SPIMat U = SPI::diag(baseflow.U);
22  SPIMat Uy = SPI::diag(baseflow.Uy);
23  if(grid.ytype==UltraS){// modify baseflow to be in C^(2) coefficient space
24  //SPI::SPIMat S1S0That = grid.S1*grid.S0*grid.That;
25  U = grid.S1S0That*U*grid.T;
26  Uy = grid.S1S0That*Uy*grid.T;
27  }
28  SPI::SPIMat d((i*alpha*U)+((k2/Re)*I)-(1./Re)*grid.Dyy,"d");
29  SPI::SPIMat L=-i*SPI::block({
30  {d, Uy, O, i*alpha*I}, // u-mom
31  {O, d, O, grid.Dy }, // v-mom
32  {O, O, d, i*beta*I }, // w-mom
33  {i*alpha*I, grid.Dy, i*beta*I, O } // continuity
34  })();//,"L");
35  SPI::SPIMat M=SPI::block({
36  {I, O, O, O}, // u-mom
37  {O, I, O, O}, // v-mom
38  {O, O, I, O}, // w-mom
39  {O, O, O, O} // continuity
40  })();//,"M");
41  // set BCs
42  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
43  //PetscInt rowBCs[] = {0,n-1,2*n,3*n-1,3*n,4*n-1}; // u,v,w at wall and freestream
44  L();
45  M();
46  //for(PetscInt rowi : rowBCs){
48  //L.eye_row(rowi);
49  //M.eye_row(rowi);
50  //L(rowi,rowi,1.0);
51  //M(rowi,rowi,60.0);
52  //M();
53  //L();
54  //for(PetscInt j=0; j<4*n; ++j){
55  //L(rowi,j,0.0,INSERT_VALUES);
56  //M(rowi,j,0.0,INSERT_VALUES);
57  //}
58  //}
59  //for(PetscInt rowi : rowBCs){
60  //L(rowi,rowi,1.0);
61  //M(rowi,rowi,60.0);
62  //}
63  if(grid.ytype==UltraS){
64  std::vector<PetscInt> rowBCs = {n-2,n-1,2*n-2,2*n-1,3*n-2,3*n-1}; // u,v,w at wall and freestream and dv'/dy = 0 at wall // before Permutation matrix reordering
65  //std::vector<PetscInt> rowBCs = {0,1,n,n+1,2*n,2*n+1}; // u,v,w at wall and freestream
66  L.zero_rows(rowBCs);
67  M.zero_rows(rowBCs);
68  for(PetscInt j=0; j<n; ++j){
69  L(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
70  L(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the freestream
71  L(rowBCs[2],n+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
72  L(rowBCs[3],n+j,grid.T(n-1,j,PETSC_TRUE)); // v at the freestream
73  L(rowBCs[4],2*n+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
74  L(rowBCs[5],2*n+j,grid.T(n-1,j,PETSC_TRUE));// w at the freestream
75  //M(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
76  //M(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the freestream
77  //M(rowBCs[2],n+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
78  //M(rowBCs[3],n+j,grid.T(n-1,j,PETSC_TRUE)); // v at the freestream
79  //M(rowBCs[4],2*n+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
80  //M(rowBCs[5],2*n+j,grid.T(n-1,j,PETSC_TRUE));// w at the freestream
81  //L(rowBCs[6],n+j,(grid.T*inv(grid.S0)*inv(grid.S1)*grid.Dy)(0,j,PETSC_TRUE));// v at the freestream
82  }
83  L();
84  M();
85  // reorder matrix rows to reduce LU factorization due to extra numerical pivoting
86  // otherwise you need to add -mat_mumps_icntl_14 25 to the simulation command line call to increase factorization space
87  //SPIMat o2xnm2 = zeros(2,n-2);
88  //SPIMat onm2x2 = zeros(n-2,2);
89  //SPIMat o2x2 = zeros(2,2);
90  //SPIMat Ptmp = block({
91  //{zeros(2,n-2),zeros(2,2)},
92  //{eye(n-2),zeros(n-2,2)}
93  //});
94  //Ptmp(0,n-2,1.0);
95  //Ptmp(1,n-1,1.0);
96  //Ptmp();
97  SPIMat P = block({
98  {grid.P,O,O,O},
99  {O,grid.P,O,O},
100  {O,O,grid.P,O},
101  //{O,O,O,eye(n)},
102  {O,O,O,eye(n)},
103  })();
104  //P.print();
105  L = P*L;
106  M = P*M;
107 
108  //L.print();
109 
110  //save(L,"L_from_LST_temporal_UltraS.dat");
111  //save(M,"M_from_LST_temporal_UltraS.dat");
112  }
113  else{
114  std::vector<PetscInt> rowBCs = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
115  L.eye_rows(rowBCs);
116  M.eye_rows(rowBCs);
117  }
118  L();
119  M();
120  //L.print();
121  //M.print();
122  //L();
123  //M();
124  SPI::SPIVec eig_vec(grid.y.rows*4,"q");
125  SPI::SPIVec eigl_vec(grid.y.rows*4,"q");
126  //PetscScalar omega;
127  // std::tie(eigenvalue,eigenfunction) = SPI::eig(L,M,0.3-0.0001*i); // doesn't work because M is singular
128  if(q.flag_init){
129  std::tie(omega,eigl_vec,eig_vec) = SPI::eig_init(M,L,1./(params.omega),q.conj(),q);
130  omega = 1./omega; // invert
131  }else{
132  //std::cout<<"here right before eig"<<std::endl;
133  //std::tie(omega,eigl_vec,eig_vec) = SPI::eig(M,L,1./params.omega);
134  //std::tie(omega,eig_vec) = SPI::eig_right(M,L,1./params.omega);
135  //std::cout<<"here right after eig"<<std::endl;
136  //std::tie(omega,eigl_vec,eig_vec) = SPI::eig(L,M,params.omega);
137  std::tie(omega,eig_vec) = SPI::eig_right(L,M,params.omega);
138  //std::tie(omega,eig_vec) = SPI::eig_right(M,L,1./params.omega);
139  //omega = 1./omega; // invert
140  //std::cout<<"here right after eig 2"<<std::endl;
141  }
142  params.omega = omega;
143  //SPI::printfc("ω is: %g+%gi",omega);
144  return std::make_tuple(omega,eig_vec);
145  }
147  std::tuple<PetscScalar, SPIVec> LST_spatial(
148  SPIparams &params,
149  SPIgrid1D &grid,
150  SPIbaseflow &baseflow,
151  SPIVec q
152  ){
153  PetscInt n = grid.ny;
154  PetscScalar Re = params.Re;
155  PetscScalar alpha = params.alpha;
156  PetscScalar omega = params.omega;
157  PetscScalar beta = params.beta;
158  PetscScalar i=PETSC_i;
159  SPIMat &O = grid.O;
160  SPIMat &I = grid.I;
161  SPIMat U = SPI::diag(baseflow.U);
162  SPIMat Uy = SPI::diag(baseflow.Uy);
163  if(grid.ytype==SPI::UltraS){
164  SPIMat &S1S0That = grid.S1S0That;
165  SPIMat &T = grid.T;
166  U = S1S0That*U*T;
167  Uy = S1S0That*Uy*T;
168  }
169  SPIMat &Dy = grid.Dy;
170  SPIMat &Dyy = grid.Dyy;
171  SPIMat d(Dyy + i*Re*omega*I - beta*beta*I,"d");
172  //SPI::SPIMat d("d");
173  //d = Dyy + i*Re*omega*I - beta*beta*I;
174  if(1){
175  SPI::SPIMat L0(SPI::block({
176  {d, -Re*Uy, O, O},
177  {O, d, O, -Re*Dy},
178  {O, O, d, -i*Re*beta*I},
179  {O, Dy, i*beta*I, O}
180  })(),"L0");
181  SPI::SPIMat L1(SPI::block({
182  {-i*Re*U, O, O, -i*Re*I},
183  {O, -i*Re*U, O, O},
184  {O, O, -i*Re*U, O},
185  {i*I, O, O, O}
186  })(),"L1");
187  SPI::SPIMat L2(SPI::block({
188  {-I, O, O, O},
189  {O, -I, O, O},
190  {O, O, -I, O},
191  {O, O, O, O}
192  })(),"L2");
193  // set BCs
194  if(grid.ytype==SPI::UltraS){
195  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
196  std::vector<PetscInt> rowBCs = {n-2,n-1,2*n-2,2*n-1,3*n-2,3*n-1}; // u,v,w at wall and freestream
197  L0.zero_rows(rowBCs);
198  L1.zero_rows(rowBCs);
199  L2.zero_rows(rowBCs);
200  for(PetscInt j=0; j<n; ++j){
201  L0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
202  L0(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the freestream
203  L0(rowBCs[2],n+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
204  L0(rowBCs[3],n+j,grid.T(n-1,j,PETSC_TRUE)); // v at the freestream
205  L0(rowBCs[4],2*n+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
206  L0(rowBCs[5],2*n+j,grid.T(n-1,j,PETSC_TRUE)); // w at the freestream
207  }
208  // assemble
209  L0();
210  L1();
211  L2();
212  SPIMat P = block({
213  {grid.P,O,O,O},
214  {O,grid.P,O,O},
215  {O,O,grid.P,O},
216  {O,O,O,eye(n)},
217  })();
218  // reorder rows for UltraSpherical
219  L0 = P*L0;
220  L1 = P*L1;
221  L2 = P*L2;
222  //std::cout<<"Re-ordered rows for UltraS"<<std::endl;
223  //save(U,"U_from_LST_spatial_UltraS.dat");
224  //save(Uy,"Uy_from_LST_spatial_UltraS.dat");
225  //save(baseflow.U,"U_and_Uy_from_LST_spatial_UltraS.h5");
226  //save(baseflow.Uy,"U_and_Uy_from_LST_spatial_UltraS.h5");
227  //baseflow.Uy.print();
228  //save(L0,"L0_from_LST_spatial_UltraS.dat");
229  //save(L1,"L1_from_LST_spatial_UltraS.dat");
230  //save(L2,"L2_from_LST_spatial_UltraS.dat");
231  }else{
232  PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
233  for(PetscInt rowi : rowBCs){
234  //SPI::printf(std::to_string(rowi));
235  L0.eye_row(rowi);
236  L1.zero_row(rowi);
237  L2.zero_row(rowi);
238  //L0(rowi,rowi,1.0);
239  //L2(rowi,rowi,60.0);
240  }
241  }
242  // assemble
243  L0();
244  L1();
245  L2();
246  //L0.print();
247  //L1.print();
248  //L2.print();
249  SPI::SPIVec eig_vec(grid.y.rows*4,"q");
250  //SPI::SPIVec eigl_vec(grid.y.rows*4,"q");
251  //PetscScalar omega;
252  // std::tie(eigenvalue,eigenfunction) = SPI::eig(L,M,0.3-0.0001*i); // doesn't work because M is singular
253  if(q.flag_init){
254  std::tie(alpha,eig_vec) = SPI::polyeig_init({L0,L1,L2},alpha,q);
255  }else{
256  std::tie(alpha,eig_vec) = SPI::polyeig({L0,L1,L2},alpha);
257  }
258  //alpha = alpha; // invert
259  params.alpha = alpha;
260  //SPI::printfc("α is: %g+%gi",alpha);
261  return std::make_tuple(alpha,eig_vec);
262  }
263  else if(0){
264  SPI::SPIMat L(SPI::block({
265  {O, O, O, O, I, O, O, O},
266  {O, O, O, O, O, I, O, O},
267  {O, O, O, O, O, O, I, O},
268  {O, O, O, O, O, O, O, I},
269  {d, -Re*Uy, O, O, -i*Re*U, O, O, -i*Re*I},
270  {O, d, O, -Re*Dy, O, -i*Re*U, O, O},
271  {O, O, d, -i*Re*beta*I,O, O, -i*Re*U, O},
272  {O, Dy, i*beta*I, O, i*I, O, O, O},
273  }),"L");
274  SPI::SPIMat M(SPI::block({
275  {I, O, O, O, O, O, O, O},
276  {O, I, O, O, O, O, O, O},
277  {O, O, I, O, O, O, O, O},
278  {O, O, O, I, O, O, O, O},
279  {O, O, O, O, I, O, O, O},
280  {O, O, O, O, O, I, O, O},
281  {O, O, O, O, O, O, I, O},
282  {O, O, O, O, O, O, O, O},
283  }),"M");
284  // set BCs
285  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
286  PetscInt rowBCs[] = {4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
287  for(PetscInt rowi : rowBCs){
288  //SPI::printf(std::to_string(rowi));
289  L.zero_row(rowi);
290  M.zero_row(rowi);
291  L(rowi,rowi,1.0);
292  M(rowi,rowi,60.0*i);
293  }
294  SPI::SPIVec eig_vec(grid.y.rows*8,"q");
295  SPI::SPIVec eigl_vec(grid.y.rows*8,"q");
296  if(q.flag_init){
297  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig_init(L,M,alpha,q.conj(),q);
298  }else{
299  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig(M,L,1./alpha);
300  }
301  alpha = alpha; // invert
302  params.alpha = alpha;
303  //SPI::printfc("α is: %g+%gi",alpha);
304  return std::make_tuple(alpha,eig_vec);
305  }
306  else{
307  d = 1./Re*(Dyy - beta*beta*I) + i*omega*I;
308  SPI::SPIMat L(SPI::block({
309  //u v w p vx wx
310  {O, i*Dy, -beta*I, O, O, O}, // continuity
311  {O, O, O, O, I, O}, // v-sub
312  {O, O, O, O, O, I}, // w-sub
313  {-i*d, i*(Uy-U*Dy),beta*U, O, -1./Re*Dy, -i*beta/Re*I}, // u-mom
314  {O, i*Re*d, O, -i*Re*Dy, Re*U, O}, // v-mom
315  {O, O, i*Re*d, beta*Re*I, O, Re*U}, // w-mom
316  }),"L");
317  SPI::SPIMat M(SPI::block({
318  {I, O, O, O, O, O},
319  {O, I, O, O, O, O},
320  {O, O, I, O, O, O},
321  {O, O, O, I, O, O},
322  {O, O, O, O, I, O},
323  {O, O, O, O, O, I}
324  }),"M");
325  // set BCs
326  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
327  PetscInt rowBCs[] = {0*n,1*n-1,1*n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
328  for(PetscInt rowi : rowBCs){
329  //SPI::printf(std::to_string(rowi));
330  L.zero_row(rowi);
331  M.zero_row(rowi);
332  L(rowi,rowi,1.0);
333  M(rowi,rowi,60.0*i);
334  }
335  SPI::SPIVec eig_vec(grid.y.rows*8,"q");
336  SPI::SPIVec eigl_vec(grid.y.rows*8,"q");
337  if(q.flag_init){
338  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig_init(L,M,alpha,q.conj(),q,1e-16,20000);
339  }else{
340  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig(L,M,alpha,1e-12,20000);
341  }
342  //alpha = 1./alpha; // invert
343  params.alpha = alpha;
344  //SPI::printfc("α is: %g+%gi",alpha);
345  return std::make_tuple(alpha,eig_vec);
346  }
347  }
349  std::tuple<PetscScalar, PetscScalar, SPIVec, SPIVec> LST_spatial_cg(
350  SPIparams &params,
351  SPIgrid1D &grid,
352  SPIbaseflow &baseflow
353  ){
354  PetscInt n = grid.ny;
355  PetscScalar Re = params.Re;
356  PetscScalar alpha = params.alpha;
357  PetscScalar omega = params.omega;
358  PetscScalar beta = params.beta;
359  PetscScalar i=PETSC_i;
360  SPIMat &O = grid.O;
361  SPIMat &I = grid.I;
362  SPIMat U = SPI::diag(baseflow.U);
363  SPIMat Uy = SPI::diag(baseflow.Uy);
364  if(grid.ytype==SPI::UltraS){
365  SPIMat &S1S0That = grid.S1S0That;
366  SPIMat &T = grid.T;
367  U = S1S0That*U*T;
368  Uy = S1S0That*Uy*T;
369  }
370  SPIMat &Dy = grid.Dy;
371  SPIMat &Dyy = grid.Dyy;
372  SPIMat d(Dyy + i*Re*omega*I - beta*beta*I,"d");
373  //SPI::SPIMat d("d");
374  //d = Dyy + i*Re*omega*I - beta*beta*I;
375  if(1){
376  SPI::SPIMat L0(SPI::block({
377  {d, -Re*Uy, O, O},
378  {O, d, O, -Re*Dy},
379  {O, O, d, -i*Re*beta*I},
380  {O, Dy, i*beta*I, O}
381  })(),"L0");
382  SPI::SPIMat L1(SPI::block({
383  {-i*Re*U, O, O, -i*Re*I},
384  {O, -i*Re*U, O, O},
385  {O, O, -i*Re*U, O},
386  {i*I, O, O, O}
387  })(),"L1");
388  SPI::SPIMat L2(SPI::block({
389  {-I, O, O, O},
390  {O, -I, O, O},
391  {O, O, -I, O},
392  {O, O, O, O}
393  })(),"L2");
394  SPI::SPIMat L2_noBCs(SPI::block({
395  {-I, O, O, O},
396  {O, -I, O, O},
397  {O, O, -I, O},
398  {O, O, O, O}
399  })(),"L2");
400  SPI::SPIMat L2_noBCs_physical(SPI::block({
401  {-eye(n), O, O, O},
402  {O, -eye(n), O, O},
403  {O, O, -eye(n), O},
404  {O, O, O, O}
405  })(),"L2");
406  //std::vector<PetscInt> rowBCs_physical = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
407  //L2_noBCs_physical.zero_rows(rowBCs_physical);
408  // set BCs
409  if(grid.ytype==SPI::UltraS){
410  //SPIMat S1S0T(grid.S1*grid.S0*grid.T,"S1*S0*T");
411  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
412  std::vector<PetscInt> rowBCs = {n-2,n-1,2*n-2,2*n-1,3*n-2,3*n-1}; // u,v,w at wall and freestream
413  L0.zero_rows(rowBCs);
414  L1.zero_rows(rowBCs);
415  L2.zero_rows(rowBCs);
416  for(PetscInt j=0; j<n; ++j){
417  L0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
418  L0(rowBCs[1],j,grid.T(n-1,j,PETSC_TRUE)); // u at the freestream
419  L0(rowBCs[2],n+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
420  L0(rowBCs[3],n+j,grid.T(n-1,j,PETSC_TRUE)); // v at the freestream
421  L0(rowBCs[4],2*n+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
422  L0(rowBCs[5],2*n+j,grid.T(n-1,j,PETSC_TRUE)); // w at the freestream
423  }
424  // assemble
425  L0();
426  L1();
427  L2();
428  SPIMat P = block({
429  {grid.P,O,O,O},
430  {O,grid.P,O,O},
431  {O,O,grid.P,O},
432  {O,O,O,eye(n)},
433  })();
434  // reorder rows for UltraSpherical
435  L0 = P*L0;
436  L1 = P*L1;
437  L2 = P*L2;
438  //std::cout<<"Re-ordered rows for UltraS"<<std::endl;
439  //save(U,"U_from_LST_spatial_UltraS.dat");
440  //save(Uy,"Uy_from_LST_spatial_UltraS.dat");
441  //save(baseflow.U,"U_and_Uy_from_LST_spatial_UltraS.h5");
442  //save(baseflow.Uy,"U_and_Uy_from_LST_spatial_UltraS.h5");
443  //baseflow.Uy.print();
444  //save(L0,"L0_from_LST_spatial_UltraS.dat");
445  //save(L1,"L1_from_LST_spatial_UltraS.dat");
446  //save(L2,"L2_from_LST_spatial_UltraS.dat");
447  }else{
448  PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
449  for(PetscInt rowi : rowBCs){
450  //SPI::printf(std::to_string(rowi));
451  L0.eye_row(rowi);
452  L1.zero_row(rowi);
453  L2.zero_row(rowi);
454  //L0(rowi,rowi,1.0);
455  //L2(rowi,rowi,60.0);
456  }
457  }
458  // assemble
459  L0();
460  L1();
461  L2();
462  //L0.print();
463  //L1.print();
464  //L2.print();
465  SPI::SPIVec eig_vec(n*8,"q");
466  SPI::SPIVec eigl_vec(n*8,"ql");
467  SPIMat O4(zeros(4*n,4*n));
468  //SPIMat I4(eye(4*n));
469  SPIMat I4(block({
470  {I,O,O,O},
471  {O,I,O,O},
472  {O,O,I,O},
473  {O,O,O,I}}
474  )(),"I4");
475  SPIMat L(block({
476  {O4, I4},
477  {L0, L1}})(),"L");
478  SPIMat M(block({
479  {I4, O4},
480  {O4, -L2}})(),"M");
481  //SPIMat M_noBCs(block({
482  //{I4, O4},
483  //{O4, -L2_noBCs}})(),"M noBCs");
484  //SPI::SPIVec eigl_vec(grid.y.rows*4,"q");
485  //PetscScalar omega;
486  // std::tie(eigenvalue,eigenfunction) = SPI::eig(L,M,0.3-0.0001*i); // doesn't work because M is singular
487  std::tie(alpha,eigl_vec,eig_vec) = eig(L,M,alpha);
488  //alpha = alpha; // invert
489  params.alpha = alpha;
490  // get cg
491  SPIMat dLdomega4(block({
492  {i*Re*I,O, O, O},
493  {O, i*Re*I, O, O},
494  {O, O, i*Re*I, O},
495  {O, O, O, O},
496  })(),"dLdomega4");
497  // set BCs for physical values
498  //dLdomega4_physical.zero_rows(rowBCs_physical);
499  SPIMat dLdomega(block({
500  {O4, O4},
501  {dLdomega4, O4}})(),"dLdomega");
502  //SPIVec eigl_vecconj(eigl_vec.conj());
503  //eigl_vec.conj();
504  PetscScalar cg;
505  if(grid.ytype==UltraS){
506  SPIMat M_noBCs_physical(block({
507  {eye(4*n), O4},
508  {O4, -L2_noBCs_physical}})(),"M noBCs physical");
509  SPIMat dLdomega4_physical(block({
510  {i*Re*eye(n), O, O, O},
511  {O, i*Re*eye(n),O, O},
512  {O, O, i*Re*eye(n),O},
513  {O, O, O, O},
514  })(),"dLdomega4");
515  SPIMat dLdomega_physical(block({
516  {O4, O4},
517  {dLdomega4_physical, O4}})(),"dLdomega");
518  //cg = ((M*eig_vec).dot(eigl_vec)) / ((dLdomega*eig_vec).dot(eigl_vec));
519  //SPIVec eigl_vecconj(eigl_vec.conj());
520  //eigl_vec.conj();
521  //SPIMat S0invS1inv(block({
522  //{grid.S0invS1inv,O,O,O,O,O,O,O},
523  //{O,grid.S0invS1inv,O,O,O,O,O,O},
524  //{O,O,grid.S0invS1inv,O,O,O,O,O},
525  //{O,O,O,grid.S0invS1inv,O,O,O,O},
526  //{O,O,O,O,grid.S0invS1inv,O,O,O},
527  //{O,O,O,O,O,grid.S0invS1inv,O,O},
528  //{O,O,O,O,O,O,grid.S0invS1inv,O},
529  //{O,O,O,O,O,O,O,grid.S0invS1inv},
530  //})());
531  SPIMat T(block({
532  {grid.T,O,O,O,O,O,O,O},
533  {O,grid.T,O,O,O,O,O,O},
534  {O,O,grid.T,O,O,O,O,O},
535  {O,O,O,grid.T,O,O,O,O},
536  {O,O,O,O,grid.T,O,O,O},
537  {O,O,O,O,O,grid.T,O,O},
538  {O,O,O,O,O,O,grid.T,O},
539  {O,O,O,O,O,O,O,grid.T},
540  })());
541  //SPIMat That(block({
542  //{grid.That,O,O,O,O,O,O,O},
543  //{O,grid.That,O,O,O,O,O,O},
544  //{O,O,grid.That,O,O,O,O,O},
545  //{O,O,O,grid.That,O,O,O,O},
546  //{O,O,O,O,grid.That,O,O,O},
547  //{O,O,O,O,O,grid.That,O,O},
548  //{O,O,O,O,O,O,grid.That,O},
549  //{O,O,O,O,O,O,O,grid.That},
550  //})());
551  //cg = integrate(((eig_vec))*conj(S0invS1inv*M*eigl_vec),grid) / integrate((S0invS1inv*(dLdomega*eig_vec))*conj(eigl_vec),grid);
552  //SPI::printfc("in LST_spatial_cg UltraS cg = %.10f + %.10fi",cg);
553  //cg = integrate(That*((T*((eig_vec)))*conj(T*(S0invS1inv*(M*eigl_vec)))),grid) / integrate(That*((T*(S0invS1inv*(dLdomega*eig_vec)))*conj(T*eigl_vec)),grid);
554  //SPI::printfc("in LST_spatial_cg UltraS cg = %.10f + %.10fi",cg);
555  //SPIgrid1D grid2(grid.y,"grid",SPI::Chebyshev);
556  //cg = integrate(((T*(S0invS1inv*(M*eig_vec)))*conj(T*eigl_vec)),grid2) / integrate(((T*(S0invS1inv*(dLdomega*eig_vec)))*conj(T*eigl_vec)),grid2);
557  //SPI::printfc("in LST_spatial_cg UltraSp cg = %.10f + %.10fi",cg);
558  //cg = ((T*((eig_vec))).dot(T*(S0invS1inv*(M*eigl_vec)))) / ((T*(S0invS1inv*(dLdomega*eig_vec))).dot(T*eigl_vec));
559  //SPI::printfc("in LST_spatial_cg UltraS to physical dot cg = %.10f + %.10fi",cg);
560  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",eig_vec.dot(eigl_vec));
561  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",eig_vec.dot(eigl_vec));
562  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",eig_vec.dot(eig_vec));
563  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",(T*eig_vec).dot(T*eig_vec));
564  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",(T*eig_vec).dot(T*eigl_vec));
565  //SPI::printfc("in LST_spatial_cg vecdot cg = %.10f + %.10fi",(T*eig_vec).dot(T*eigl_vec));
566  //cg = integrate((S0invS1inv*(M_noBCs*eig_vec))*conj(eigl_vec),grid) / integrate((S0invS1inv*(dLdomega*eig_vec))*conj(eigl_vec),grid);
567  //SPI::printfc("in LST_spatial_cg noBCs cg = %.10f + %.10fi",cg);
568  cg = (M_noBCs_physical*(T*eig_vec)).dot(T*eigl_vec) / ((dLdomega_physical*(T*eig_vec))).dot(T*eigl_vec);
569  SPI::printfc("in LST_spatial_cg physical noBCs cg = %.10f + %.10fi",cg);
570  //cg = (T*eig_vec).dot(M_noBCs_physical*(T*eigl_vec)) / ((dLdomega_physical*(T*eig_vec))).dot(T*eigl_vec);
571  //SPI::printfc("in LST_spatial_cg physical noBCs cg = %.10f + %.10fi",cg);
572  //SPIVec eig_vec2(eig_vec);
573  //SPIVec eigl_vec2(eigl_vec);
574  //eig_vec2 /= (T*eig_vec).dot(M_noBCs_physical*(T*eigl_vec));
575  //cg = 1. / ((dLdomega_physical*(T*eig_vec2))).dot(T*eigl_vec);
576  //SPI::printfc("in LST_spatial_cg physical noBCs normed cg = %.10f + %.10fi",cg);
577  //SPIVec eig_vec_p(T*eig_vec);
578  //SPIVec eigl_vec_p(T*eigl_vec);
579  //eig_vec_p /= SPI::L2(eig_vec_p);
580  //eigl_vec_p /= SPI::L2(eigl_vec_p);
581  //cg = (((eig_vec_p)).dot(M_noBCs_physical*eigl_vec_p)) / ((dLdomega_physical*(eig_vec_p)).dot(eigl_vec_p));
582  //SPI::printfc("in LST_spatial_cg physical noBCs normed2 cg = %.10f + %.10fi",cg);
583  //SPIVec eig_vec_2H(conj(eig_vec));
584  //SPIVec eig_vec_2(eig_vec);
585  //SPIVec eigl_vec_2H(conj(eigl_vec));
586  //SPIVec eigl_vec_2(eigl_vec);
587  //std::cout<<"in LST_spatial_cg UltraS normed2 to norm value = "<<(integrate(eig_vec_2*conj(eig_vec_2),grid))<<std::endl;
588  //std::cout<<"in LST_spatial_cg UltraSp normed2 to norm value = "<<(integrate(That*((T*eig_vec_2)*conj(T*eig_vec_2)),grid))<<std::endl;
589  //std::cout<<"in LST_spatial_cg UltraS normed2 to sqrt(norm value) = "<<sqrt(integrate(eig_vec_2*conj(eig_vec_2),grid))<<std::endl;
590  //PetscScalar normval = integrate(eig_vec_2*conj(eig_vec_2),grid);
591  //if(PetscRealPart(normval)<=0.){
592  //eig_vec_2 /= sqrt(normval);
593  //}else{
594  //eig_vec_2 /= sqrt(normval);
595  //}
596  //eig_vec_2H = conj(eig_vec_2);
597  //eig_vec_2.conj();
598  //eigl_vec_2 = ((1./integrate(eigl_vec_2*conj(eigl_vec_2),grid))*eigl_vec_2);
599  //eigl_vec_2H = conj(eigl_vec_2);
600  //eigl_vec_2.conj();
601  //SPI::printfc("in LST_spatial_cg UltraS normed2 norm = %.10f + %.10fi",integrate(eig_vec_2*conj(eig_vec_2),grid));
602  //SPI::printfc("in LST_spatial_cg UltraS normed2 norm = %.10f + %.10fi",integrate(eigl_vec_2*conj(eigl_vec_2),grid));
603  //cg = integrate(((eig_vec_2))*((M*eigl_vec_2).conj()),grid) / integrate((dLdomega*(eig_vec_2))*conj(eigl_vec_2),grid);
604  //SPI::printfc("in LST_spatial_cg UltraS normed2 cg = %.10f + %.10fi",cg);
605  //M.H();
606  //eig_vec2 /= integrate(eig_vec2*(S0invS1inv*(M*conj(eigl_vec))),grid);
607  //SPI::printfc("in LST_spatial_cg physical noBCs normed normval = %.10f + %.10fi",integrate(eig_vec2*(S0invS1inv*(M*conj(eigl_vec))),grid));
608  //cg = 1. / integrate((S0invS1inv*(dLdomega*(eig_vec2)))*conj(eigl_vec),grid);
609  //M.H();
610  //SPI::printfc("in LST_spatial_cg physical noBCs normed cg = %.10f + %.10fi",cg);
611  //grid.ytype=SPI::Chebyshev;
612  //cg = integrate((M_noBCs_physical*(T*eig_vec))*conj(T*eigl_vec),grid) / integrate((dLdomega_physical*(T*eig_vec))*conj(T*eigl_vec),grid);
613  //SPI::printfc("in LST_spatial_cg physical integrate cg = %.10f + %.10fi",cg);
614  //grid.ytype=SPI::UltraS;
615  }else{
616  //cg = ((M*eig_vec).dot(eigl_vec)) / ((dLdomega*eig_vec).dot(eigl_vec));
617  //SPIVec eigl_vecconj(eigl_vec.conj());
618  //eigl_vec.conj();
619  cg = integrate((M*eig_vec)*conj(eigl_vec),grid) / integrate((dLdomega*eig_vec)*conj(eigl_vec),grid);
620  SPI::printfc("in LST_spatial_cg physical integrate cg = %.10f + %.10fi",cg);
621  cg = (((eig_vec)).dot(M*eigl_vec)) / ((dLdomega*(eig_vec)).dot(eigl_vec));
622  SPI::printfc("in LST_spatial_cg physical dot cg = %.10f + %.10fi",cg);
623  }
624  //SPI::printfc("α is: %g+%gi",alpha);
625  return std::make_tuple(alpha,cg,eigl_vec,eig_vec);
626  }
627  else if(0){
628  SPI::SPIMat L(SPI::block({
629  {O, O, O, O, I, O, O, O},
630  {O, O, O, O, O, I, O, O},
631  {O, O, O, O, O, O, I, O},
632  {O, O, O, O, O, O, O, I},
633  {d, -Re*Uy, O, O, -i*Re*U, O, O, -i*Re*I},
634  {O, d, O, -Re*Dy, O, -i*Re*U, O, O},
635  {O, O, d, -i*Re*beta*I,O, O, -i*Re*U, O},
636  {O, Dy, i*beta*I, O, i*I, O, O, O},
637  }),"L");
638  SPI::SPIMat M(SPI::block({
639  {I, O, O, O, O, O, O, O},
640  {O, I, O, O, O, O, O, O},
641  {O, O, I, O, O, O, O, O},
642  {O, O, O, I, O, O, O, O},
643  {O, O, O, O, I, O, O, O},
644  {O, O, O, O, O, I, O, O},
645  {O, O, O, O, O, O, I, O},
646  {O, O, O, O, O, O, O, O},
647  }),"M");
648  // set BCs
649  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
650  PetscInt rowBCs[] = {4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
651  for(PetscInt rowi : rowBCs){
652  //SPI::printf(std::to_string(rowi));
653  L.zero_row(rowi);
654  M.zero_row(rowi);
655  L(rowi,rowi,1.0);
656  M(rowi,rowi,60.0*i);
657  }
658  SPI::SPIVec eig_vec(grid.y.rows*8,"q");
659  SPI::SPIVec eigl_vec(grid.y.rows*8,"q");
660  //if(q.flag_init){
661  //std::tie(alpha,eigl_vec,eig_vec) = SPI::eig_init(L,M,alpha,q.conj(),q);
662  //}else{
663  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig(M,L,1./alpha);
664  //}
665  alpha = alpha; // invert
666  params.alpha = alpha;
667  //SPI::printfc("α is: %g+%gi",alpha);
668  //return std::make_tuple(alpha,eig_vec);
669  return std::make_tuple(alpha,alpha,eig_vec,eig_vec);
670  }
671  else{
672  d = 1./Re*(Dyy - beta*beta*I) + i*omega*I;
673  SPI::SPIMat L(SPI::block({
674  //u v w p vx wx
675  {O, i*Dy, -beta*I, O, O, O}, // continuity
676  {O, O, O, O, I, O}, // v-sub
677  {O, O, O, O, O, I}, // w-sub
678  {-i*d, i*(Uy-U*Dy),beta*U, O, -1./Re*Dy, -i*beta/Re*I}, // u-mom
679  {O, i*Re*d, O, -i*Re*Dy, Re*U, O}, // v-mom
680  {O, O, i*Re*d, beta*Re*I, O, Re*U}, // w-mom
681  }),"L");
682  SPI::SPIMat M(SPI::block({
683  {I, O, O, O, O, O},
684  {O, I, O, O, O, O},
685  {O, O, I, O, O, O},
686  {O, O, O, I, O, O},
687  {O, O, O, O, I, O},
688  {O, O, O, O, O, I}
689  }),"M");
690  // set BCs
691  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
692  PetscInt rowBCs[] = {0*n,1*n-1,1*n,2*n-1,2*n,3*n-1}; // u,v,w at wall and freestream
693  for(PetscInt rowi : rowBCs){
694  //SPI::printf(std::to_string(rowi));
695  L.zero_row(rowi);
696  M.zero_row(rowi);
697  L(rowi,rowi,1.0);
698  M(rowi,rowi,60.0*i);
699  }
700  SPI::SPIVec eig_vec(grid.y.rows*8,"q");
701  SPI::SPIVec eigl_vec(grid.y.rows*8,"q");
702  //if(q.flag_init){
703  //std::tie(alpha,eigl_vec,eig_vec) = SPI::eig_init(L,M,alpha,q.conj(),q,1e-16,20000);
704  //}else{
705  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig(L,M,alpha,1e-12,20000);
706  //}
707  //alpha = 1./alpha; // invert
708  params.alpha = alpha;
709  //SPI::printfc("α is: %g+%gi",alpha);
710  //return std::make_tuple(alpha,eig_vec);
711  return std::make_tuple(alpha,alpha,eig_vec,eig_vec);
712  }
713  }
715  std::tuple<PetscScalar, PetscScalar, SPIVec, SPIVec> LSTNP_spatial(
716  SPIparams &params,
717  SPIgrid1D &grid,
718  SPIbaseflow &baseflow,
719  SPIVec ql,
720  SPIVec qr
721  ){
722  PetscInt ny = grid.ny;
723  PetscScalar Re = params.Re;
724  PetscScalar Reinv = 1./Re;
725  PetscScalar alpha = params.alpha;
726  PetscScalar omega = params.omega;
727  PetscScalar beta = params.beta;
728  PetscScalar i=PETSC_i;
729  SPIMat &O = grid.O;
730  SPIMat &I = grid.I;
731  SPIMat U = SPI::diag(baseflow.U);
732  SPIMat V = SPI::diag(baseflow.V);
733  SPIMat Ux = SPI::diag(baseflow.Ux);
734  SPIMat Uy = SPI::diag(baseflow.Uy);
735  SPIMat Uxy = SPI::diag(baseflow.Uxy);
736  SPIMat Vy = SPI::diag(baseflow.Vy);
737  SPIMat W = SPI::diag(baseflow.W);
738  SPIMat Wx = SPI::diag(baseflow.Wx);
739  SPIMat Wy = SPI::diag(baseflow.Wy);
740  SPIMat Wxy = SPI::diag(baseflow.Wxy);
741  SPIMat P = SPI::diag(baseflow.P);
742  if(grid.ytype==UltraS){ // map baseflow to C^(2) coefficient space using S0, S1, T, That operators
743  //SPIMat S1S0That = grid.S1*grid.S0*grid.That;
744  SPIMat &S1S0That = grid.S1S0That;
745  //SPIMat &S1 = grid.S1;
746  SPIMat &T = grid.T;
747  //SPIMat &That = grid.That;
748  U = S1S0That*U*T;
749  V = S1S0That*V*T;
750  Ux = S1S0That*Ux*T;
751  Uy = S1S0That*Uy*T;
752  Uxy = S1S0That*Uxy*T;
753  Vy = S1S0That*Vy*T;
754  W = S1S0That*W*T;
755  Wx = S1S0That*Wx*T;
756  Wy = S1S0That*Wy*T;
757  Wxy = S1S0That*Wxy*T;
758  P = S1S0That*P*T;
759  //std::cout<<"Altered the baseflow"<<std::endl;
760  }
761  SPIMat &Dy = grid.Dy;
762  SPIMat &Dyy = grid.Dyy;
763  SPIMat VDy(V*Dy,"V*Dy");
764  if(grid.ytype==UltraS){
765  //SPIMat Vp = SPI::diag(baseflow.V); // physical to physical baseflow
766  //SPIMat &S0invS1inv = grid.S0invS1inv;
767  //SPIMat Dyp = grid.T*((S0invS1inv*Dy)*grid.That); // from physical->Chebyshev->C^(2)->Chebyshev->physical
768  //VDy = grid.S1S0That*((grid.T*(S0invS1inv*V))*(grid.T*(S0invS1inv*Dy)));
769  //VDy = grid.S1S0That*((Vp*Dyp)*grid.T); // from Chebyshev to C^(2)
770  //VDy = O;
771  //VDy = V*S0invS1inv*Dy;
772  VDy = grid.S1S0That*(diag(baseflow.V) * (grid.T*grid.S0invS1inv*Dy*grid.That)) * grid.T;
773  }
774  SPIMat d(VDy - Reinv*Dyy + i*beta*W + (-i*omega + beta*beta*Reinv)*I,"d");
775  if(1){
776  SPIMat ReinvI = Reinv*I;
777  SPIMat ABC2(SPI::block({
778  {ReinvI, O, O, O},
779  {O, ReinvI, O, O},
780  {O, O, ReinvI, O},
781  {O, O, O, O}
782  })(),"ABC2");
783  SPIMat O4(zeros(4*ny,4*ny));
784  SPIMat D2(O4,"D2");
785  SPIMat dA2(O4,"dA2");
786  SPIMat iU = i*U;
787  SPIMat iI = i*I;
788  SPIMat ABC1(SPI::block({
789  {iU, O, O, iI},
790  {O, iU, O, O },
791  {O, O, iU, O },
792  {iI, O, O, O }
793  })(),"ABC1");
794  SPIMat D1(O4,"D1");
795  SPIMat iUx = i*Ux;
796  SPIMat dA1(SPI::block({
797  {iUx, O, O, O },
798  {O, iUx, O, O },
799  {O, O, iUx, O },
800  {O, O, O, O }
801  })(),"dA1");
802  SPIMat ibetaI = i*beta*I;
803  SPIMat ABC0(SPI::block({
804  {d+Ux, Uy, O, O },
805  {O, d+Vy, O, Dy },
806  {Wx, Wy, d, ibetaI},
807  {O, Dy, ibetaI, O }
808  })(),"ABC0");
809  SPIMat D0(SPI::block({
810  {U, O, O, I },
811  {O, U, O, O },
812  {O, O, U, O },
813  {I, O, O, O }
814  })(),"D0");
815  SPIMat ibetaWx = i*beta*Wx;
816  SPIMat dA0(SPI::block({
817  {ibetaWx, Uxy, O, O },
818  {O, ibetaWx, O, O },
819  {O, Wxy, ibetaWx, O },
820  {O, O, O, O }
821  })(),"dA0");
822  // set BCs
823  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
824  //PetscInt rowBCs[] = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
825  //for(PetscInt rowi : rowBCs){
826  // //SPI::printf(std::to_string(rowi));
827  // ABC2.zero_row(rowi);
828  // ABC1.zero_row(rowi);
829  // dA1.zero_row(rowi);
830  // ABC0.zero_row(rowi);
831  // D0.zero_row(rowi);
832  // dA0.zero_row(rowi);
833 
834  // ABC2(rowi,rowi,1.0);
835  // ABC1(rowi,rowi,1.0);
836  // dA1(rowi,rowi,1.0);
837  // ABC0(rowi,rowi,1.0);
838  // D0(rowi,rowi,1.0);
839  // dA0(rowi,rowi,1.0);
840  //}
841  // alternative way to set the BCs
842  if(grid.ytype==UltraS){
843  std::vector<PetscInt> rowBCs = {1*ny-2,1*ny-1,2*ny-2,2*ny-1,3*ny-2,3*ny-1}; // u,v,w at wall and freestream
844  ABC2.zero_rows(rowBCs);
845  ABC1.zero_rows(rowBCs);
846  dA1.zero_rows(rowBCs);
847  ABC0.zero_rows(rowBCs);
848  D0.zero_rows(rowBCs);
849  dA0.zero_rows(rowBCs);
850  for(PetscInt j=0; j<ny; ++j){
851  ABC0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
852  ABC0(rowBCs[1],j,grid.T(ny-1,j,PETSC_TRUE)); // u at the freestream
853  ABC0(rowBCs[2],ny+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
854  ABC0(rowBCs[3],ny+j,grid.T(ny-1,j,PETSC_TRUE)); // v at the freestream
855  ABC0(rowBCs[4],2*ny+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
856  ABC0(rowBCs[5],2*ny+j,grid.T(ny-1,j,PETSC_TRUE)); // w at the freestream
857  }
858  ABC2();
859  ABC1();
860  dA1();
861  ABC0();
862  D0();
863  dA0();
864  //std::cout<<"set BCs"<<std::endl;
865  //ABC0.print();
866  //grid.T.print();
867  }
868  else{
869  std::vector<PetscInt> rowBCs = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
870  ABC2.zero_rows(rowBCs); // is eye_rows in python? but isn't is supposed to be zero?
871  ABC1.zero_rows(rowBCs); // is eye_rows in python... but isn't it supposed to be zero?
872  dA1.zero_rows(rowBCs); // is eye_rows in python... but isn't it supposed to be zero?
873  ABC0.eye_rows(rowBCs);
874  D0.zero_rows(rowBCs); // is eye_rows in python... but isn't it supposed to be zero?
875  dA0.zero_rows(rowBCs); // is eye_rows in python... but isn't it supposed to be zero?
876  }
877 
878  // inflate for derivatives in streamwise direction
879  SPIMat L0(block({
880  {ABC0, D0 },
881  {dA0, ABC0}})(),"L0");
882  SPIMat L1(block({
883  {ABC1, D1 },
884  {dA1, ABC1}})(),"L1");
885  SPIMat L2(block({
886  {ABC2, D2 },
887  {dA2, ABC2}})(),"L2");
888  // inflate due to polynomial eigenvalue problem
889  SPIMat O8(zeros(8*ny,8*ny));
890  //SPIMat I8(eye(8*ny));
891  SPIMat I8(block({
892  {I,O,O,O,O,O,O,O},
893  {O,I,O,O,O,O,O,O},
894  {O,O,I,O,O,O,O,O},
895  {O,O,O,I,O,O,O,O},
896  {O,O,O,O,I,O,O,O},
897  {O,O,O,O,O,I,O,O},
898  {O,O,O,O,O,O,I,O},
899  {O,O,O,O,O,O,O,I},
900  })(),"I8");
901  SPIMat L(block({
902  {O8, I8},
903  {L0, L1}})(),"L");
904  SPIMat M(block({
905  {I8, O8},
906  {O8, -L2}})(),"M");
907  //L.print();
908  //M.print();
909  //SPIMat L(block({
910  //{O8, I8},
911  //{-L0, -L1}}),"L");
912  //SPIMat M(block({
913  //{I8, O8},
914  //{O8, L2}}),"M");
915  //SPIMat L(block({
916  //{-L0, O8},
917  //{O8, I8}}),"L");
918  //SPIMat M(block({
919  //{L1, L2},
920  //{I8, O8}}),"M");
921  //PetscScalar a=10.5,b=20.5;
922  //SPIMat L(block({
923  //{-b*L0, a*I8},
924  //{-a*L0, -a*L1+b*I8}}),"L");
925  //SPIMat M(block({
926  //{a*I8+b*L1, b*L2},
927  //{b*I8, a*L2}}),"M");
928 
929 
930  SPI::SPIVec eig_vec(grid.y.rows*16,"q");
931  SPI::SPIVec eigl_vec(grid.y.rows*16,"q");
932  if(ql.flag_init){
933  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig_init(L,M,alpha,ql,qr);
934  }else{
935  std::tie(alpha,eigl_vec,eig_vec) = SPI::eig(L,M,alpha);
936  //std::tie(alpha,eig_vec) = SPI::eig_right(L,M,alpha);
937  }
938  //SPIMat L(block(
939  //alpha = alpha; // invert
940  params.alpha = alpha;
941  //SPI::printfc("α is: %g+%gi",alpha);
942  SPIMat dLdomega4(block({
943  {-i*I, O, O, O},
944  {O, -i*I, O, O},
945  {O, O, -i*I, O},
946  {O, O, O, O},
947  })(),"dLdomega 4nx4n");
948  // inflate for non-parallel
949  SPIMat dLdomega8(block({
950  {dLdomega4, O4},
951  {O4, dLdomega4}
952  })(),"dLdomega 8nx8n");
953  // inflate for polynomial eigenvalue problem
954  SPIMat dLdomega(block({
955  {O8, O8},
956  {dLdomega8, O8}
957  })(),"dLdomega 16nx16n");
958  PetscScalar cg;
959  if(grid.ytype==UltraS){
960  // TODO this group velocity is incorrect for the test case... need to fix!
961  //SPIMat &s = grid.S0invS1inv;
962  SPIMat &t = grid.T;
963  //SPIMat S(block({
964  // {s,O,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
965  // {O,s,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
966  // {O,O,s,O,O,O,O,O,O,O,O,O,O,O,O,O},
967  // {O,O,O,s,O,O,O,O,O,O,O,O,O,O,O,O},
968  // {O,O,O,O,s,O,O,O,O,O,O,O,O,O,O,O},
969  // {O,O,O,O,O,s,O,O,O,O,O,O,O,O,O,O},
970  // {O,O,O,O,O,O,s,O,O,O,O,O,O,O,O,O},
971  // {O,O,O,O,O,O,O,s,O,O,O,O,O,O,O,O},
972  // {O,O,O,O,O,O,O,O,s,O,O,O,O,O,O,O},
973  // {O,O,O,O,O,O,O,O,O,s,O,O,O,O,O,O},
974  // {O,O,O,O,O,O,O,O,O,O,s,O,O,O,O,O},
975  // {O,O,O,O,O,O,O,O,O,O,O,s,O,O,O,O},
976  // {O,O,O,O,O,O,O,O,O,O,O,O,s,O,O,O},
977  // {O,O,O,O,O,O,O,O,O,O,O,O,O,s,O,O},
978  // {O,O,O,O,O,O,O,O,O,O,O,O,O,O,s,O},
979  // {O,O,O,O,O,O,O,O,O,O,O,O,O,O,O,s}
980  // })(),"S0invS1inv");
981  SPIMat T(block({
982  {t,O,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
983  {O,t,O,O,O,O,O,O,O,O,O,O,O,O,O,O},
984  {O,O,t,O,O,O,O,O,O,O,O,O,O,O,O,O},
985  {O,O,O,t,O,O,O,O,O,O,O,O,O,O,O,O},
986  {O,O,O,O,t,O,O,O,O,O,O,O,O,O,O,O},
987  {O,O,O,O,O,t,O,O,O,O,O,O,O,O,O,O},
988  {O,O,O,O,O,O,t,O,O,O,O,O,O,O,O,O},
989  {O,O,O,O,O,O,O,t,O,O,O,O,O,O,O,O},
990  {O,O,O,O,O,O,O,O,t,O,O,O,O,O,O,O},
991  {O,O,O,O,O,O,O,O,O,t,O,O,O,O,O,O},
992  {O,O,O,O,O,O,O,O,O,O,t,O,O,O,O,O},
993  {O,O,O,O,O,O,O,O,O,O,O,t,O,O,O,O},
994  {O,O,O,O,O,O,O,O,O,O,O,O,t,O,O,O},
995  {O,O,O,O,O,O,O,O,O,O,O,O,O,t,O,O},
996  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,t,O},
997  {O,O,O,O,O,O,O,O,O,O,O,O,O,O,O,t}
998  })(),"T");
999  SPIMat dLdomega4_physical(block({
1000  {-i*eye(ny), O, O, O},
1001  {O, -i*eye(ny), O, O},
1002  {O, O, -i*eye(ny), O},
1003  {O, O, O, O},
1004  })(),"dLdomega 4nx4n");
1005  // inflate for non-parallel
1006  SPIMat dLdomega8_physical(block({
1007  {dLdomega4_physical,O4},
1008  {O4, dLdomega4_physical}
1009  })(),"dLdomega 8nx8n");
1010  // inflate for polynomial eigenvalue problem
1011  SPIMat dLdomega_physical(block({
1012  {O8, O8},
1013  {dLdomega8_physical,O8}
1014  })(),"dLdomega 16nx16n");
1015  SPIMat ABC2_physical(SPI::block({
1016  {Reinv*eye(ny), O, O, O},
1017  {O, Reinv*eye(ny), O, O},
1018  {O, O, Reinv*eye(ny), O},
1019  {O, O, O, O}
1020  })(),"ABC2");
1021  SPIMat L2_physical(block({
1022  {ABC2_physical, D2 },
1023  {dA2, ABC2_physical}})(),"L2");
1024  SPIMat M_physical(block({
1025  {eye(8*ny), O8},
1026  {O8, -L2_physical}})(),"M");
1027  //SPIVec eigl_vecconj(eigl_vec.conj());
1028  //eigl_vec.conj();
1029  //SPIVec tmp((S16*(M*eig_vec))*eigl_vecconj);
1030  //PetscScalar numerator = integrate((S*(M*eig_vec))*(eigl_vecconj),grid);
1031  //PetscScalar denominator = integrate((S*(dLdomega*eig_vec))*(eigl_vecconj),grid);
1032  //eigl_vec.conj();
1033  //PetscScalar numerator = (T*S*(M*eig_vec)).dot(T*eigl_vec);
1034  //PetscScalar denominator = (T*S*(dLdomega*eig_vec)).dot(T*eigl_vec);
1035  PetscScalar numerator = (((M_physical*(T*eig_vec))).dot(T*eigl_vec));
1036  PetscScalar denominator = (((dLdomega_physical*(T*eig_vec))).dot(T*eigl_vec));
1037  cg = numerator/denominator;
1038  }
1039  else{
1040  cg = ((M*eig_vec).dot(eigl_vec)) / ((dLdomega*eig_vec).dot(eigl_vec));
1041  }
1042  return std::make_tuple(alpha,cg,eigl_vec,eig_vec);
1043  }
1044  }
1046  std::tuple<PetscScalar, SPIVec> LSTNP_spatial_right(
1047  SPIparams &params,
1048  SPIgrid1D &grid,
1049  SPIbaseflow &baseflow,
1050  SPIVec qr
1051  ){
1052  PetscInt ny = grid.ny;
1053  PetscScalar Re = params.Re;
1054  PetscScalar Reinv = 1.0/Re;
1055  PetscScalar alpha = params.alpha;
1056  PetscScalar omega = params.omega;
1057  PetscScalar beta = params.beta;
1058  PetscScalar i=PETSC_i;
1059  SPIMat &O = grid.O;
1060  SPIMat &I = grid.I;
1061  SPIMat U = SPI::diag(baseflow.U);
1062  SPIMat V = SPI::diag(baseflow.V);
1063  SPIMat Ux = SPI::diag(baseflow.Ux);
1064  SPIMat Uy = SPI::diag(baseflow.Uy);
1065  SPIMat Uxy = SPI::diag(baseflow.Uxy);
1066  SPIMat Vy = SPI::diag(baseflow.Vy);
1067  SPIMat W = SPI::diag(baseflow.W);
1068  SPIMat Wx = SPI::diag(baseflow.Wx);
1069  SPIMat Wy = SPI::diag(baseflow.Wy);
1070  SPIMat Wxy = SPI::diag(baseflow.Wxy);
1071  SPIMat P = SPI::diag(baseflow.P);
1072  if(grid.ytype==SPI::UltraS){
1073  SPIMat &S1S0That = grid.S1S0That;
1074  SPIMat &T = grid.T;
1075  U = S1S0That*U*T;
1076  V = S1S0That*V*T;
1077  Ux = S1S0That*Ux*T;
1078  Uy = S1S0That*Uy*T;
1079  Uxy = S1S0That*Uxy*T;
1080  Vy = S1S0That*Vy*T;
1081  W = S1S0That*W*T;
1082  Wx = S1S0That*Wx*T;
1083  Wy = S1S0That*Wy*T;
1084  Wxy = S1S0That*Wxy*T;
1085  P = S1S0That*P*T;
1086  }
1087  SPIMat &Dy = grid.Dy;
1088  SPIMat &Dyy = grid.Dyy;
1089  SPIMat VDy(V*Dy,"V*Dy");
1090  if(grid.ytype==UltraS){
1091  //SPIMat &S0invS1inv = grid.S0invS1inv;
1092  //VDy = V*S0invS1inv*Dy;
1093  //VDy = grid.S1S0That*((grid.T*(S0invS1inv*V))*(grid.T*(S0invS1inv*Dy)));
1094  VDy = grid.S1S0That*(diag(baseflow.V) * (grid.T*grid.S0invS1inv*Dy*grid.That)) * grid.T;
1095  //VDy = grid.O;
1096  }
1097  SPIMat d(VDy - Reinv*Dyy + i*beta*W + (-i*omega + beta*beta*Reinv)*I,"d");
1098  if(1){
1099  SPIMat ReinvI(Reinv*I);
1100  SPIMat ABC2(SPI::block({
1101  {ReinvI, O, O, O},
1102  {O, ReinvI, O, O},
1103  {O, O, ReinvI, O},
1104  {O, O, O, O}
1105  })(),"ABC2");
1106  SPIMat O4(zeros(4*ny,4*ny));
1107  SPIMat D2(O4,"D2");
1108  SPIMat dA2(O4,"dA2");
1109  SPIMat iU(i*U);
1110  SPIMat iI(i*I);
1111  SPIMat ABC1(SPI::block({
1112  {iU, O, O, iI },
1113  {O, iU, O, O },
1114  {O, O, iU, O },
1115  {iI, O, O, O }
1116  })(),"ABC1");
1117  SPIMat D1(O4,"D1");
1118  SPIMat iUx(i*Ux);
1119  SPIMat dA1(SPI::block({
1120  {iUx, O, O, O },
1121  {O, iUx, O, O },
1122  {O, O, iUx, O },
1123  {O, O, O, O }
1124  })(),"dA1");
1125  SPIMat ibetaI(i*beta*I);
1126  SPIMat ABC0(SPI::block({
1127  {d+Ux, Uy, O, O },
1128  {O, d+Vy, O, Dy },
1129  {Wx, Wy, d, ibetaI},
1130  {O, Dy, ibetaI,O }
1131  })(),"ABC0");
1132  SPIMat D0(SPI::block({
1133  {U, O, O, I },
1134  {O, U, O, O },
1135  {O, O, U, O },
1136  {I, O, O, O }
1137  })(),"D0");
1138  SPIMat ibetaWx(i*beta*Wx);
1139  SPIMat dA0(SPI::block({
1140  {ibetaWx, Uxy, O, O },
1141  {O, ibetaWx, O, O },
1142  {O, Wxy, ibetaWx, O },
1143  {O, O, O, O }
1144  })(),"dA0");
1145  // set BCs
1146  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
1147  //PetscInt rowBCs[] = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
1148  //for(PetscInt rowi : rowBCs){
1149  // //SPI::printf(std::to_string(rowi));
1150  // ABC2.zero_row(rowi);
1151  // ABC1.zero_row(rowi);
1152  // dA1.zero_row(rowi);
1153  // ABC0.zero_row(rowi);
1154  // D0.zero_row(rowi);
1155  // dA0.zero_row(rowi);
1156 
1157  // ABC2(rowi,rowi,1.0);
1158  // ABC1(rowi,rowi,1.0);
1159  // dA1(rowi,rowi,1.0);
1160  // ABC0(rowi,rowi,1.0);
1161  // D0(rowi,rowi,1.0);
1162  // dA0(rowi,rowi,1.0);
1163  //}
1164  // alternative way to set the BCs
1165  if(grid.ytype==SPI::UltraS){
1166  std::vector<PetscInt> rowBCs = {1*ny-2,1*ny-1,2*ny-2,2*ny-1,3*ny-2,3*ny-1}; // u,v,w at wall and freestream
1167  ABC2.zero_rows(rowBCs);
1168  ABC1.zero_rows(rowBCs);
1169  dA1.zero_rows(rowBCs);
1170  ABC0.zero_rows(rowBCs);
1171  D0.zero_rows(rowBCs);
1172  dA0.zero_rows(rowBCs);
1173  for(PetscInt j=0; j<ny; ++j){
1174  ABC0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
1175  ABC0(rowBCs[1],j,grid.T(ny-1,j,PETSC_TRUE)); // u at the freestream
1176  ABC0(rowBCs[2],ny+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
1177  ABC0(rowBCs[3],ny+j,grid.T(ny-1,j,PETSC_TRUE)); // v at the freestream
1178  ABC0(rowBCs[4],2*ny+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
1179  ABC0(rowBCs[5],2*ny+j,grid.T(ny-1,j,PETSC_TRUE)); // w at the freestream
1180  }
1181  // assemble
1182  ABC0();
1183  // SPIMat P = block({
1184  // {grid.P,O,O,O},
1185  // {O,grid.P,O,O},
1186  // {O,O,grid.P,O},
1187  // {O,O,O,eye(ny)},
1188  // })();
1189  // // reorder rows for UltraSpherical
1190  // ABC2 = P*ABC2;
1191  // ABC1 = P*ABC1;
1192  // dA1 = P*dA1;
1193  // ABC0 = P*ABC0;
1194  // D0 = P*D0;
1195  // dA0 = P*dA0;
1196  }else{
1197  std::vector<PetscInt> rowBCs = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
1198  ABC2.zero_rows(rowBCs);// is eye_rows in python...
1199  ABC1.zero_rows(rowBCs);// is eye_rows in python...
1200  dA1.zero_rows(rowBCs);// is eye_rows in python...
1201  ABC0.eye_rows(rowBCs);
1202  D0.zero_rows(rowBCs);// is eye_rows in python...
1203  dA0.zero_rows(rowBCs);// is eye_rows in python...
1204  }
1205 
1206  // inflate for derivatives in streamwise direction
1207  SPIMat L0(block({
1208  {ABC0, D0 },
1209  {dA0, ABC0}})(),"L0");
1210  SPIMat L1(block({
1211  {ABC1, D1 },
1212  {dA1, ABC1}})(),"L1");
1213  SPIMat L2(block({
1214  {ABC2, D2 },
1215  {dA2, ABC2}})(),"L2");
1216  // inflate due to polynomial eigenvalue problem
1217  //if(grid.ytype==SPI::UltraS){
1218  //SPI::SPIVec eig_vec_tmp(grid.y.rows*8,"q");
1219  //std::tie(alpha,eig_vec_tmp) = SPI::polyeig({L0,L1,L2},alpha);
1220  //return std::make_tuple(alpha,eig_vec_tmp);
1221  //}
1222  SPIMat O8(zeros(8*ny,8*ny));
1223  //SPIMat I8(eye(8*ny));
1224  SPIMat I8(block({
1225  {I,O,O,O,O,O,O,O},
1226  {O,I,O,O,O,O,O,O},
1227  {O,O,I,O,O,O,O,O},
1228  {O,O,O,I,O,O,O,O},
1229  {O,O,O,O,I,O,O,O},
1230  {O,O,O,O,O,I,O,O},
1231  {O,O,O,O,O,O,I,O},
1232  {O,O,O,O,O,O,O,I},
1233  })(),"I8");
1234  SPIMat L(block({
1235  {O8, I8},
1236  {L0, L1}})(),"L");
1237  SPIMat M(block({
1238  {I8, O8},
1239  {O8, -L2}})(),"M");
1240  //SPIMat L(block({
1241  //{O8, I8},
1242  //{-L0, -L1}}),"L");
1243  //SPIMat M(block({
1244  //{I8, O8},
1245  //{O8, L2}}),"M");
1246  //SPIMat L(block({
1247  //{-L0, O8},
1248  //{O8, I8}}),"L");
1249  //SPIMat M(block({
1250  //{L1, L2},
1251  //{I8, O8}}),"M");
1252  //PetscScalar a=10.5,b=20.5;
1253  //SPIMat L(block({
1254  //{-b*L0, a*I8},
1255  //{-a*L0, -a*L1+b*I8}}),"L");
1256  //SPIMat M(block({
1257  //{a*I8+b*L1, b*L2},
1258  //{b*I8, a*L2}}),"M");
1259 
1260 
1261  SPI::SPIVec eig_vec(grid.y.rows*16,"q");
1262  if(qr.flag_init){
1263  std::tie(alpha,eig_vec) = SPI::eig_init_right(L,M,alpha,qr,1e-30);
1264  }else{
1265  std::tie(alpha,eig_vec) = SPI::eig_right(L,M,alpha,1e-30);
1266  }
1267  eig_vec.rows=4*ny;
1268  //SPIMat L(block(
1269  //alpha = alpha; // invert
1270  params.alpha = alpha;
1271  //SPI::printfc("α is: %g+%gi",alpha);
1272  return std::make_tuple(alpha,eig_vec);
1273  }
1274  }
1276  std::tuple<PetscScalar, SPIVec> LSTNP_spatial_right2(
1277  SPIparams &params,
1278  SPIgrid1D &grid,
1279  SPIbaseflow &baseflow,
1280  SPIVec qr
1281  ){
1282  PetscInt ny = grid.ny;
1283  PetscScalar Re = params.Re;
1284  PetscScalar Reinv = 1.0/Re;
1285  PetscScalar alpha = params.alpha;
1286  PetscScalar omega = params.omega;
1287  PetscScalar beta = params.beta;
1288  PetscScalar i=PETSC_i;
1289  SPIMat &O = grid.O;
1290  SPIMat &I = grid.I;
1291  SPIMat U = SPI::diag(baseflow.U);
1292  SPIMat V = SPI::diag(baseflow.V);
1293  SPIMat Ux = SPI::diag(baseflow.Ux);
1294  SPIMat Uy = SPI::diag(baseflow.Uy);
1295  SPIMat Uxy = SPI::diag(baseflow.Uxy);
1296  SPIMat Vy = SPI::diag(baseflow.Vy);
1297  SPIMat W = SPI::diag(baseflow.W);
1298  SPIMat Wx = SPI::diag(baseflow.Wx);
1299  SPIMat Wy = SPI::diag(baseflow.Wy);
1300  SPIMat Wxy = SPI::diag(baseflow.Wxy);
1301  SPIMat P = SPI::diag(baseflow.P);
1302  if(grid.ytype==SPI::UltraS){
1303  SPIMat &S1S0That = grid.S1S0That;
1304  SPIMat &T = grid.T;
1305  U = S1S0That*U*T;
1306  V = S1S0That*V*T;
1307  Ux = S1S0That*Ux*T;
1308  Uy = S1S0That*Uy*T;
1309  Uxy = S1S0That*Uxy*T;
1310  Vy = S1S0That*Vy*T;
1311  W = S1S0That*W*T;
1312  Wx = S1S0That*Wx*T;
1313  Wy = S1S0That*Wy*T;
1314  Wxy = S1S0That*Wxy*T;
1315  P = S1S0That*P*T;
1316  }
1317  SPIMat &Dy = grid.Dy;
1318  SPIMat &Dyy = grid.Dyy;
1319  SPIMat VDy(V*Dy,"V*Dy");
1320  if(grid.ytype==UltraS){
1321  VDy = grid.S1S0That*(diag(baseflow.V) * (grid.T*grid.S0invS1inv*Dy*grid.That)) * grid.T;
1322  }
1323  SPIMat d(VDy - Reinv*Dyy + i*beta*W + (-i*omega + beta*beta*Reinv)*I,"d");
1324  if(1){
1325  SPIMat ReinvI(Reinv*I);
1326  SPIMat ABC2(SPI::block({
1327  {ReinvI, O, O, O},
1328  {O, ReinvI, O, O},
1329  {O, O, ReinvI, O},
1330  {O, O, O, O}
1331  })(),"ABC2");
1332  SPIMat O4(zeros(4*ny,4*ny));
1333  SPIMat D2(O4,"D2");
1334  SPIMat dA2(O4,"dA2");
1335  SPIMat iU(i*U);
1336  SPIMat iI(i*I);
1337  SPIMat ABC1(SPI::block({
1338  {iU, O, O, iI },
1339  {O, iU, O, O },
1340  {O, O, iU, O },
1341  {iI, O, O, O }
1342  })(),"ABC1");
1343  SPIMat D1(O4,"D1");
1344  SPIMat iUx(i*Ux);
1345  SPIMat dA1(SPI::block({
1346  {iUx, O, O, O },
1347  {O, iUx, O, O },
1348  {O, O, iUx, O },
1349  {O, O, O, O }
1350  })(),"dA1");
1351  SPIMat ibetaI(i*beta*I);
1352  SPIMat ABC0(SPI::block({
1353  {d+Ux, Uy, O, O },
1354  {O, d+Vy, O, Dy },
1355  {Wx, Wy, d, ibetaI},
1356  {O, Dy, ibetaI,O }
1357  })(),"ABC0");
1358  SPIMat D0(SPI::block({
1359  {U, O, O, I },
1360  {O, U, O, O },
1361  {O, O, U, O },
1362  {I, O, O, O }
1363  })(),"D0");
1364  SPIMat ibetaWx(i*beta*Wx);
1365  SPIMat dA0(SPI::block({
1366  {ibetaWx, Uxy, O, O },
1367  {O, ibetaWx, O, O },
1368  {O, Wxy, ibetaWx, O },
1369  {O, O, O, O }
1370  })(),"dA0");
1371  // alternative way to set the BCs
1372  if(grid.ytype==SPI::UltraS){
1373  std::vector<PetscInt> rowBCs = {1*ny-2,1*ny-1,2*ny-2,2*ny-1,3*ny-2,3*ny-1}; // u,v,w at wall and freestream
1374  ABC2.zero_rows(rowBCs);
1375  ABC1.zero_rows(rowBCs);
1376  dA1.zero_rows(rowBCs);
1377  ABC0.zero_rows(rowBCs);
1378  D0.zero_rows(rowBCs);
1379  dA0.zero_rows(rowBCs);
1380  for(PetscInt j=0; j<ny; ++j){
1381  ABC0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
1382  ABC0(rowBCs[1],j,grid.T(ny-1,j,PETSC_TRUE)); // u at the freestream
1383  ABC0(rowBCs[2],ny+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
1384  ABC0(rowBCs[3],ny+j,grid.T(ny-1,j,PETSC_TRUE)); // v at the freestream
1385  ABC0(rowBCs[4],2*ny+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
1386  ABC0(rowBCs[5],2*ny+j,grid.T(ny-1,j,PETSC_TRUE)); // w at the freestream
1387  }
1388  // assemble
1389  ABC0();
1390  SPIMat P = block({
1391  {grid.P,O,O,O},
1392  {O,grid.P,O,O},
1393  {O,O,grid.P,O},
1394  {O,O,O,eye(ny)},
1395  })();
1396  // reorder rows for UltraSpherical
1397  ABC2 = P*ABC2;
1398  ABC1 = P*ABC1;
1399  dA1 = P*dA1;
1400  ABC0 = P*ABC0;
1401  D0 = P*D0;
1402  dA0 = P*dA0;
1403  }else{
1404  std::vector<PetscInt> rowBCs = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
1405  ABC2.zero_rows(rowBCs);// is eye_rows in python...
1406  ABC1.zero_rows(rowBCs);// is eye_rows in python...
1407  dA1.zero_rows(rowBCs);// is eye_rows in python...
1408  ABC0.eye_rows(rowBCs);
1409  D0.zero_rows(rowBCs);// is eye_rows in python...
1410  dA0.zero_rows(rowBCs);// is eye_rows in python...
1411  }
1412 
1413  // inflate for derivatives in streamwise direction
1414  SPIMat L0(block({
1415  {ABC0, D0 },
1416  {dA0, ABC0}})(),"L0");
1417  SPIMat L1(block({
1418  {ABC1, D1 },
1419  {dA1, ABC1}})(),"L1");
1420  SPIMat L2(block({
1421  {ABC2, D2 },
1422  {dA2, ABC2}})(),"L2");
1423  // inflate due to polynomial eigenvalue problem
1424  //if(grid.ytype==SPI::UltraS){
1425  //SPI::SPIVec eig_vec_tmp(grid.y.rows*8,"q");
1426  //std::tie(alpha,eig_vec_tmp) = SPI::polyeig({L0,L1,L2},alpha);
1427  //return std::make_tuple(alpha,eig_vec_tmp);
1428  //}
1429  SPIMat O8(zeros(8*ny,8*ny));
1430  //SPIMat I8(eye(8*ny));
1431  SPIMat I8(block({
1432  {I,O,O,O,O,O,O,O},
1433  {O,I,O,O,O,O,O,O},
1434  {O,O,I,O,O,O,O,O},
1435  {O,O,O,I,O,O,O,O},
1436  {O,O,O,O,I,O,O,O},
1437  {O,O,O,O,O,I,O,O},
1438  {O,O,O,O,O,O,I,O},
1439  {O,O,O,O,O,O,O,I},
1440  })(),"I8");
1441  SPIMat L(block({
1442  {O8, I8},
1443  {L0, L1}})(),"L");
1444  SPIMat M(block({
1445  {I8, O8},
1446  {O8, -L2}})(),"M");
1447  //SPIMat L(block({
1448  //{O8, I8},
1449  //{-L0, -L1}})(),"L");
1450  //SPIMat M(block({
1451  //{I8, O8},
1452  //{O8, L2}})(),"M");
1453  //SPIMat L(block({
1454  //{-L0, O8},
1455  //{O8, I8}})(),"L");
1456  //SPIMat M(block({
1457  //{L1, L2},
1458  //{I8, O8}})(),"M");
1459  //PetscScalar a=10.5,b=20.5;
1460  //SPIMat L(block({
1461  //{-b*L0, a*I8},
1462  //{-a*L0, -a*L1+b*I8}})(),"L");
1463  //SPIMat M(block({
1464  //{a*I8+b*L1, b*L2},
1465  //{b*I8, a*L2}})(),"M");
1466 
1467 
1468  SPI::SPIVec eig_vec(grid.y.rows*16,"q");
1469  if(qr.flag_init){
1470  std::tie(alpha,eig_vec) = SPI::eig_init_right(L,M,alpha,qr);
1471  }else{
1472  std::tie(alpha,eig_vec) = SPI::eig_right(L,M,alpha);
1473  }
1474  //SPIMat L(block(
1475  //alpha = alpha; // invert
1476  params.alpha = alpha;
1477  //SPI::printfc("α is: %g+%gi",alpha);
1478  return std::make_tuple(alpha,eig_vec);
1479  }
1480  }
1482  std::tuple<std::vector<PetscScalar>, std::vector<SPIVec>> LSTNP_spatials_right(
1483  SPIparams &params,
1484  SPIgrid1D &grid,
1485  SPIbaseflow &baseflow,
1486  std::vector<PetscScalar> &alphas,
1487  std::vector<SPIVec> &qrs
1488  ){
1489  PetscInt ny = grid.ny;
1490  PetscScalar Re = params.Re;
1491  PetscScalar Reinv = 1./Re;
1492  //PetscScalar alpha = alphas[0];
1493  std::vector<PetscScalar> _alphas;
1494  std::vector<SPIVec> eig_vecs;
1495  PetscScalar omega = params.omega;
1496  PetscScalar beta = params.beta;
1497  PetscScalar i=PETSC_i;
1498  SPIMat &O = grid.O;
1499  SPIMat &I = grid.I;
1500  SPIMat U = SPI::diag(baseflow.U);
1501  SPIMat V = SPI::diag(baseflow.V);
1502  SPIMat Ux = SPI::diag(baseflow.Ux);
1503  SPIMat Uy = SPI::diag(baseflow.Uy);
1504  SPIMat Uxy = SPI::diag(baseflow.Uxy);
1505  SPIMat Vy = SPI::diag(baseflow.Vy);
1506  SPIMat W = SPI::diag(baseflow.W);
1507  SPIMat Wx = SPI::diag(baseflow.Wx);
1508  SPIMat Wy = SPI::diag(baseflow.Wy);
1509  SPIMat Wxy = SPI::diag(baseflow.Wxy);
1510  SPIMat P = SPI::diag(baseflow.P);
1511  if(grid.ytype==SPI::UltraS){
1512  SPIMat &S1S0That = grid.S1S0That;
1513  SPIMat &T = grid.T;
1514  U = S1S0That*U*T;
1515  V = S1S0That*V*T;
1516  Ux = S1S0That*Ux*T;
1517  Uy = S1S0That*Uy*T;
1518  Uxy = S1S0That*Uxy*T;
1519  Vy = S1S0That*Vy*T;
1520  W = S1S0That*W*T;
1521  Wx = S1S0That*Wx*T;
1522  Wy = S1S0That*Wy*T;
1523  Wxy = S1S0That*Wxy*T;
1524  P = S1S0That*P*T;
1525  }
1526  SPIMat &Dy = grid.Dy;
1527  SPIMat &Dyy = grid.Dyy;
1528  SPIMat VDy(V*Dy,"V*Dy");
1529  if(grid.ytype==UltraS){
1530  //SPIMat &S0invS1inv = grid.S0invS1inv;
1531  //VDy = V*S0invS1inv*Dy;
1532  //VDy = grid.S1S0That*((grid.T*(S0invS1inv*V))*(grid.T*(S0invS1inv*Dy)));
1533  //VDy = grid.O;
1534  VDy = grid.S1S0That*(diag(baseflow.V) * (grid.T*grid.S0invS1inv*Dy*grid.That)) * grid.T;
1535  }
1536  SPIMat d(VDy - Reinv*Dyy + i*beta*W + (-i*omega + beta*beta*Reinv)*I,"d");
1537  if(1){
1538  SPIMat ABC2(SPI::block({
1539  {Reinv*I, O, O, O},
1540  {O, Reinv*I,O, O},
1541  {O, O, Reinv*I,O},
1542  {O, O, O, O}
1543  })(),"ABC2");
1544  SPIMat O4(zeros(4*ny,4*ny));
1545  SPIMat D2(O4,"D2");
1546  SPIMat dA2(O4,"dA2");
1547  SPIMat ABC1(SPI::block({
1548  {i*U, O, O, i*I},
1549  {O, i*U, O, O },
1550  {O, O, i*U, O },
1551  {i*I, O, O, O }
1552  })(),"ABC1");
1553  SPIMat D1(O4,"D1");
1554  SPIMat dA1(SPI::block({
1555  {i*Ux, O, O, O },
1556  {O, i*Ux, O, O },
1557  {O, O, i*Ux, O },
1558  {O, O, O, O }
1559  })(),"dA1");
1560  SPIMat ABC0(SPI::block({
1561  {d+Ux, Uy, O, O },
1562  {O, d+Vy, O, Dy },
1563  {Wx, Wy, d, i*beta*I},
1564  {O, Dy, i*beta*I,O }
1565  })(),"ABC0");
1566  SPIMat D0(SPI::block({
1567  {U, O, O, I },
1568  {O, U, O, O },
1569  {O, O, U, O },
1570  {I, O, O, O }
1571  })(),"D0");
1572  SPIMat dA0(SPI::block({
1573  {i*beta*Wx, Uxy, O, O },
1574  {O, i*beta*Wx, O, O },
1575  {O, Wxy, i*beta*Wx, O },
1576  {O, O, O, O }
1577  })(),"dA0");
1578  // set BCs
1579  //PetscInt rowBCs[] = {0,n-1,n,2*n-1,2*n,3*n-1,4*n,5*n-1,5*n,6*n-1,6*n,7*n-1}; // u,v,w at wall and freestream
1580  //PetscInt rowBCs[] = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
1581  //for(PetscInt rowi : rowBCs){
1582  // //SPI::printf(std::to_string(rowi));
1583  // ABC2.zero_row(rowi);
1584  // ABC1.zero_row(rowi);
1585  // dA1.zero_row(rowi);
1586  // ABC0.zero_row(rowi);
1587  // D0.zero_row(rowi);
1588  // dA0.zero_row(rowi);
1589 
1590  // ABC2(rowi,rowi,1.0);
1591  // ABC1(rowi,rowi,1.0);
1592  // dA1(rowi,rowi,1.0);
1593  // ABC0(rowi,rowi,1.0);
1594  // D0(rowi,rowi,1.0);
1595  // dA0(rowi,rowi,1.0);
1596  //}
1597  // alternative way to set the BCs
1598  if(grid.ytype==SPI::UltraS){
1599  std::vector<PetscInt> rowBCs = {1*ny-2,1*ny-1,2*ny-2,2*ny-1,3*ny-2,3*ny-1}; // u,v,w at wall and freestream
1600  ABC2.zero_rows(rowBCs);
1601  ABC1.zero_rows(rowBCs);
1602  dA1.zero_rows(rowBCs);
1603  ABC0.zero_rows(rowBCs);
1604  D0.zero_rows(rowBCs);
1605  dA0.zero_rows(rowBCs);
1606  for(PetscInt j=0; j<ny; ++j){
1607  ABC0(rowBCs[0],j,grid.T(0,j,PETSC_TRUE)); // u at the wall
1608  ABC0(rowBCs[1],j,grid.T(ny-1,j,PETSC_TRUE)); // u at the freestream
1609  ABC0(rowBCs[2],ny+j,grid.T(0,j,PETSC_TRUE)); // v at the wall
1610  ABC0(rowBCs[3],ny+j,grid.T(ny-1,j,PETSC_TRUE)); // v at the freestream
1611  ABC0(rowBCs[4],2*ny+j,grid.T(0,j,PETSC_TRUE)); // w at the wall
1612  ABC0(rowBCs[5],2*ny+j,grid.T(ny-1,j,PETSC_TRUE)); // w at the freestream
1613  }
1614  // assemble
1615  ABC0();
1616  // SPIMat P = block({
1617  // {grid.P,O,O,O},
1618  // {O,grid.P,O,O},
1619  // {O,O,grid.P,O},
1620  // {O,O,O,eye(ny)},
1621  // })();
1622  // // reorder rows for UltraSpherical
1623  // ABC2 = P*ABC2;
1624  // ABC1 = P*ABC1;
1625  // dA1 = P*dA1;
1626  // ABC0 = P*ABC0;
1627  // D0 = P*D0;
1628  // dA0 = P*dA0;
1629  }else{
1630  std::vector<PetscInt> rowBCs = {0*ny,1*ny-1,1*ny,2*ny-1,2*ny,3*ny-1}; // u,v,w at wall and freestream
1631  ABC2.zero_rows(rowBCs);// is eye_rows in python...
1632  ABC1.zero_rows(rowBCs);// is eye_rows in python...
1633  dA1.zero_rows(rowBCs);// is eye_rows in python...
1634  ABC0.eye_rows(rowBCs);
1635  D0.zero_rows(rowBCs);// is eye_rows in python...
1636  dA0.zero_rows(rowBCs);// is eye_rows in python...
1637  }
1638 
1639  // inflate for derivatives in streamwise direction
1640  SPIMat L0(block({
1641  {ABC0, D0 },
1642  {dA0, ABC0}})(),"L0");
1643  SPIMat L1(block({
1644  {ABC1, D1 },
1645  {dA1, ABC1}})(),"L1");
1646  SPIMat L2(block({
1647  {ABC2, D2 },
1648  {dA2, ABC2}})(),"L2");
1649  // inflate due to polynomial eigenvalue problem
1650  //if(grid.ytype==SPI::UltraS){
1651  //SPI::SPIVec eig_vec_tmp(grid.y.rows*8,"q");
1652  //std::tie(alpha,eig_vec_tmp) = SPI::polyeig({L0,L1,L2},alpha);
1653  //return std::make_tuple(alpha,eig_vec_tmp);
1654  //}
1655  SPIMat O8(zeros(8*ny,8*ny));
1656  //SPIMat I8(eye(8*ny));
1657  SPIMat I8(block({
1658  {I,O,O,O,O,O,O,O},
1659  {O,I,O,O,O,O,O,O},
1660  {O,O,I,O,O,O,O,O},
1661  {O,O,O,I,O,O,O,O},
1662  {O,O,O,O,I,O,O,O},
1663  {O,O,O,O,O,I,O,O},
1664  {O,O,O,O,O,O,I,O},
1665  {O,O,O,O,O,O,O,I},
1666  })(),"I8");
1667  SPIMat L(block({
1668  {O8, I8},
1669  {L0, L1}})(),"L");
1670  SPIMat M(block({
1671  {I8, O8},
1672  {O8, -L2}})(),"M");
1673  //SPIMat L(block({
1674  //{O8, I8},
1675  //{-L0, -L1}}),"L");
1676  //SPIMat M(block({
1677  //{I8, O8},
1678  //{O8, L2}}),"M");
1679  //SPIMat L(block({
1680  //{-L0, O8},
1681  //{O8, I8}}),"L");
1682  //SPIMat M(block({
1683  //{L1, L2},
1684  //{I8, O8}}),"M");
1685  //PetscScalar a=10.5,b=20.5;
1686  //SPIMat L(block({
1687  //{-b*L0, a*I8},
1688  //{-a*L0, -a*L1+b*I8}}),"L");
1689  //SPIMat M(block({
1690  //{a*I8+b*L1, b*L2},
1691  //{b*I8, a*L2}}),"M");
1692 
1693 
1694  //SPI::SPIVec eig_vec(grid.y.rows*16,"q");
1695  //if(qr.flag_init){
1696  std::tie(_alphas,eig_vecs) = SPI::eig_init_rights(L,M,alphas,qrs);
1697  //}else{
1698  //std::tie(alphas[0],eig_vecs[0]) = SPI::eig_right(L,M,alpha);
1699  //}
1700  //SPIMat L(block(
1701  //alpha = alpha; // invert
1702  //params.alpha = alpha;
1703  //SPI::printfc("α is: %g+%gi",alpha);
1704  return std::make_tuple(_alphas,eig_vecs);
1705  }
1706  }
1707 
1708 }
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::SPIbaseflow::Uxy
SPIVec Uxy
Definition: SPIbaseflow.hpp:44
SPI::SPIparams::Re
PetscScalar Re
Reynolds number.
Definition: SPIparams.hpp:13
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::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::SPIgrid1D
Class to contain various grid parameters.
Definition: SPIgrid.hpp:35
SPI::SPIgrid1D::T
SPIMat T
Chebyshev operator taking it from Chebyshev coefficients to physical space.
Definition: SPIgrid.hpp:58
SPI::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::SPIVec
Definition: SPIVec.hpp:13
SPI::SPIbaseflow::Vy
SPIVec Vy
Definition: SPIbaseflow.hpp:47
SPI::SPIVec::flag_init
PetscBool flag_init
flag if it has been initialized
Definition: SPIVec.hpp:25
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
Definition: SPIbaseflow.hpp:16
SPI::conj
SPIVec conj(const SPIVec &A)
Definition: SPIVec.cpp:622
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::SPIbaseflow::Wy
SPIVec Wy
Definition: SPIbaseflow.hpp:50
SPI::SPIbaseflow::P
SPIVec P
Definition: SPIbaseflow.hpp:52
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::SPIbaseflow::Uy
SPIVec Uy
Definition: SPIbaseflow.hpp:43
SPI::SPIbaseflow::Wx
SPIVec Wx
Definition: SPIbaseflow.hpp:49
SPI::SPIgrid1D::That
SPIMat That
Chebyshev operator taking it from physical space to Chebyshev coefficients.
Definition: SPIgrid.hpp:59
SPI::eye
SPIMat eye(const PetscInt n)
create, form, and return identity matrix of size n
Definition: SPIMat.cpp:697
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::SPIgrid1D::S0invS1inv
SPIMat S0invS1inv
[in] inverse of S0^-1 * S1^-1
Definition: SPIgrid.hpp:56
SPI::dot
PetscScalar dot(SPIVec x, SPIVec y)
take the inner product of the two vectors (i.e. y^H x) where ^H is the complex conjugate transpose
Definition: SPIVec.cpp:788
SPI::SPIbaseflow::Ux
SPIVec Ux
Definition: SPIbaseflow.hpp:41
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::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::SPIbaseflow::Wxy
SPIVec Wxy
Definition: SPIbaseflow.hpp:51
SPI::SPIgrid1D::ny
PetscInt ny
Definition: SPIgrid.hpp:39
SPI::SPIgrid1D::O
SPIMat O
zero matrix same size as derivative operators
Definition: SPIgrid.hpp:64
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::eig_init_rights
std::tuple< std::vector< PetscScalar >, std::vector< SPIVec > > eig_init_rights(const SPIMat &A, const SPIMat &B, const std::vector< PetscScalar > targets, const std::vector< SPIVec > &qrs, 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:1521
SPI::SPIbaseflow::W
SPIVec W
Definition: SPIbaseflow.hpp:48
SPI::SPIbaseflow::U
SPIVec U
Definition: SPIbaseflow.hpp:39
SPI::SPIgrid1D::P
SPIMat P
row permutation matrix for UltraSpherical operators to shift rows from bottom to top to reduce LU fac...
Definition: SPIgrid.hpp:57
SPI::SPIparams
Definition: SPIparams.hpp:8
SPI::SPIgrid1D::ytype
gridtype ytype
type of grid
Definition: SPIgrid.hpp:48
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::SPIgrid1D::S1S0That
SPIMat S1S0That
UltraSpherical helper matrix S1*S0*That for baseflow.
Definition: SPIgrid.hpp:55
SPILST.hpp
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::SPIparams::beta
PetscScalar beta
beta spanwise wavenumber
Definition: SPIparams.hpp:14
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::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::SPIMat::zero_row
SPIMat & zero_row(const PetscInt row)
set a row to zero
Definition: SPIMat.cpp:451
SPI::SPIgrid1D::I
SPIMat I
identity matrix same size as derivative operators
Definition: SPIgrid.hpp:65
SPI::SPIbaseflow
Definition: SPIbaseflow.hpp:17
SPI::zeros
SPIMat zeros(const PetscInt m, const PetscInt n)
create, form, and return zeros matrix of size mxn
Definition: SPIMat.cpp:708
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
Generated by   doxygen 1.8.17