Spectral Integral Suite in C++
Ex_09.cpp
Go to the documentation of this file.
1 
70 #define SIS_USE_LAPACK
71 #include <fstream>
72 #include <iostream>
73 #include <sis.hpp>
74 
75 using namespace std;
76 
77 int main() {
78  using namespace sis;
79  int bre;
80  // Number of Chebyshev polynomials
81  N = 91;
82  sis_setup();
83 
84  // Number of Eigenvalues to compute:
85  int num_vals = 40;
86 
87  valarray<double> y, Uy, U, Uyy(N + 1);
88  // Set in cheb-point
89  setChebPts(y);
90 
91  // Velocity and derivative for PP flow
92  U = 1.0 - pow(y, 2.0);
93  Uy = -2.0 * y;
94  Uyy = -2.0;
95  double Re = 2000;
96  double kx = 1.0;
97  double kz = 0.0;
98  double k2 = kx * kx + kz * kz;
99  double k4 = k2 * k2;
100  complex<double> ii(0.0, 1.0);
101  Linop<double> Delta, Delta2, Dy;
102  LinopMat<complex<double> > Lmat(2, 2), Mmat(2, 2);
103 
104  Delta.n = 2;
105  Delta.set();
106  Delta.coef << 1.0, 0.0, -k2;
107 
108  Delta2.n = 4;
109  Delta2.set();
110  Delta2.coef << 1.0, 0.0, -2 * k2, 0.0, k4;
111 
112  Dy.n = 1;
113  Dy.set();
114  Dy.coef << 1.0, 0.0;
115 
116  // Linop<complex<double> > tempOp;
117  // Eigen::ArrayXcd tempArray;
118  // tempOp = -ii * kx * U + Delta;
119 
120  // For first type:
121  Lmat << (-ii * kx * U * Delta) + (ii * kx * Uyy) + (Delta2 / Re), 0.0, //
122  (-ii * kz * Uy), (-ii * kx * U) + (Delta / Re);
123 
124  Mmat << Delta, 0.0, //
125  0.0, 1.0;
126  BcMat<std::complex<double> > bc(6, 2);
127  bc.L << 1.0, 0.0, //
128  0.0, 1.0, //
129  Dy, 0.0, //
130  1.0, 0.0, //
131  0.0, 1.0, //
132  Dy, 0.0;
133  bc.eval << -1.0, -1.0, //
134  -1.0, -1.0, //
135  -1.0, -1.0, //
136  1.0, 1.0, //
137  1.0, 1.0, //
138  1.0, 1.0;
139 
140 
141  GeneralizedEigenSolver<complex<double> > eigs;
142  eigs.compute(Lmat, Mmat, num_vals, bc);
143  eigs.keepConverged();
144  eigs.sortByLargestReal();
145  num_vals = eigs.eigenvalues.size();
146  cout << "Eigenvalues in method 1: \n" << eigs.eigenvalues << "\n";
147  ofstream outf;
148  outf.open("data/Ex_15_1.txt");
149  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> out_to_file(num_vals,
150  2);
151  out_to_file.col(0) = eigs.eigenvalues.real();
152  out_to_file.col(1) = eigs.eigenvalues.imag();
153  outf << out_to_file;
154  outf.close();
155  std::cout << "Done 1 :::::::" << '\n';
156  // For second type:
157  Lmat.resize(4, 4);
158  Mmat.resize(4, 4);
159  /*
160  Old Mmat:
161  Mmat << 1.0, 0.0, 0.0, 0.0, //
162  0.0, 1.0, 0.0, 0.0, //
163  0.0, 0.0, 1.0, 0.0, //
164  0.0, 0.0, 0.0, 0.0;
165 */
166  Mmat << 1.0, 0.0, 0.0, 0.0, //
167  0.0, 1.0, 0.0, 0.0, //
168  0.0, 0.0, 1.0, 0.0, //
169  0.0, 0.0, 0.0, 0.0 * Delta;
170 
171  //
172  Lmat << (-ii * kx * U) + (Delta / Re), -Uy, 0.0, -ii * kx, //
173  0.0, (-ii * kx * U) + (Delta / Re), 0.0, -Dy, //
174  0.0, 0.0, (-ii * kx * U) + (Delta / Re), -ii * kz, //
175  ii * kx, Dy, ii * kz, 0.0;
176  // ii*kx, 2.0*ii * kx * Uy + Dy, ii*kz, Delta; //Mj's suggestion
177  /*
178  Old bcs:
179  Lmat.BcVec[0] = bc_dir;
180  Lmat.BcVec[1] = bc_dir;
181  Lmat.BcVec[2] = bc_dir;
182  Lmat.BcVec[3] = bc_p;
183  */
184  BcMat<std::complex<double> > Lbc(8, 4);
185  Lbc.L << 1.0, 0.0, 0.0, 0.0, //
186  1.0, 0.0, 0.0, 0.0, //
187  0.0, 1.0, 0.0, 0.0, //
188  0.0, 1.0, 0.0, 0.0, //
189  0.0, Dy, 0.0, 0.0, //
190  0.0, Dy, 0.0, 0.0, //
191  0.0, 0.0, 1.0, 0.0, //
192  0.0, 0.0, 1.0, 0.0;
193  Lbc.eval << 1.0, 0.0, 0.0, 0.0, //
194  -1.0, 0.0, 0.0, 0.0, //
195  0.0, 1.0, 0.0, 0.0, //
196  0.0, -1.0, 0.0, 0.0, //
197  0.0, 1.0, 0.0, 0.0, //
198  0.0, -1.0, 0.0, 0.0, //
199  0.0, 0.0, 1.0, 0.0, //
200  0.0, 0.0, -1.0, 0.0;
201  // std::cout << "done in " << __LINE__ << '\n';
202  num_vals = 4 * (N + 1);
203  eigs.compute(Lmat, Mmat, num_vals, Lbc);
204  eigs.keepConverged();
205  eigs.sortByLargestReal();
206  cout << "Eigenvalues in method 2: \n" << eigs.eigenvalues << "\n" << flush;
207  num_vals = eigs.eigenvalues.size();
208  outf.open("data/Ex_15_2.txt");
209  out_to_file.resize(num_vals, 2);
210  out_to_file.col(0) = eigs.eigenvalues.real();
211  out_to_file.col(1) = eigs.eigenvalues.imag();
212  outf << out_to_file;
213  outf.close();
214 return 0;
215  // Lastly, solve with conventional boundary conditions.
216  Mmat.resize(4, 4);
217  Mmat << 1.0, 0.0, 0.0, 0.0, //
218  0.0, 1.0, 0.0, 0.0, //
219  0.0, 0.0, 1.0, 0.0, //
220  0.0, 0.0, 0.0, 0.0;
221 
222  Lbc.resize(7, 4);
223  Lbc.L << 1.0, 0.0, 0.0, 0.0, //
224  1.0, 0.0, 0.0, 0.0, //
225  0.0, 1.0, 0.0, 0.0, //
226  0.0, 1.0, 0.0, 0.0, // // //
227  0.0, 0.0, 1.0, 0.0, //
228  0.0, 0.0, 1.0, 0.0, //
229  0.0, 0.0, 0.0, 1.0;
230 
231  Lbc.eval << 1.0, 0.0, 0.0, 0.0, //
232  -1.0, 0.0, 0.0, 0.0, //
233  0.0, 1.0, 0.0, 0.0, //
234  0.0, -1.0, 0.0, 0.0, //
235  0.0, 0.0, 1.0, 0.0, //
236  0.0, 0.0, -1.0, 0.0, //
237  0.0, 0.0, 0.0, -1.0;
238  num_vals = 4 * (N + 1);
239  ind = 2;
240  eigs.compute(Lmat, Mmat, num_vals, Lbc);
241  eigs.keepConverged();
242  eigs.sortByLargestReal();
243  cout << "Eigenvalues in method 3: \n" << eigs.eigenvalues << "\n" << flush;
244  num_vals = eigs.eigenvalues.size();
245  outf.open("data/Ex_15_3.txt");
246  out_to_file.resize(num_vals, 2);
247  out_to_file.col(0) = eigs.eigenvalues.real();
248  out_to_file.col(1) = eigs.eigenvalues.imag();
249  outf << out_to_file;
250  outf.close();
251  //return 0;
252 }
sis::sis_setup
void sis_setup()
Definition: sis.hpp:954
sis::BcMat
BcMat will hold general Boundary conditions as LinopMats at evealuation points, as given by operator ...
Definition: sis.hpp:529
sis::BcMat::eval
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > eval
Definition: sis.hpp:8832
sis::Linop
Linop This class creates a Linear operator to solve TPBVPs.
Definition: sis.hpp:2560
std
Definition: sis.hpp:260
ii
complex< double > ii(0.0, 1.0)
main
int main()
Definition: Ex_09.cpp:77
sis::GeneralizedEigenSolver::eigenvalues
Eigen::Matrix< std::complex< T >, Eigen::Dynamic, 1 > eigenvalues
Definition: sis.hpp:3257
sis::LinopMat
This class represents a block matrix operator. It is a matrix of operators.
Definition: sis.hpp:528
sis
Definition: sis.hpp:461
sis::Linop::n
int n
The order of the Linear differential operator.
Definition: sis.hpp:2566
sis::Linop::set
void set()
Sets the the size of coefficient array based on the order of differential equation. Suppose order is 3, then the size of array needed is 4.
Definition: sis.hpp:2733
sis::Linop::coef
Eigen::Matrix< T, Eigen::Dynamic, 1 > coef
Stores the coefficients in the differential equation.
Definition: sis.hpp:2569
sis::setChebPts
void setChebPts(std::valarray< T > &in)
This function sets points to evaluate a function so that a DCT will give represent the same function ...
Definition: sis.hpp:920
std::pow
valarray< complex< T > > pow(valarray< complex< T > > base, T power)
Definition: sis.hpp:453
sis::BcMat::L
LinopMat< T > L
Definition: sis.hpp:8831
sis::GeneralizedEigenSolver::compute
void compute(Linop< T > L, Linop< T > M, int num_vals)
Call this with an input Linear operator to solve for eigenvalues and vectors. The number of Eigen val...
Definition: sis.hpp:3267
sis::N
int N
Specifies number of Chebyshev polynomials, default N = 31.
Definition: sis.hpp:472
sis::BcMat::resize
void resize(int m_, int n_)
Definition: sis.hpp:8848
ind
int ind
Definition: sis.hpp:256
sis::GeneralizedEigenSolver
This class will solve the generalized eigenvalue problem for two linear operators. One of them can be singular.
Definition: sis.hpp:3253
sis::y
std::valarray< SIS_TYPE > y(N+1)
sis::LinopMat::resize
void resize(int r_, int c_)
This clears all contents in the LinopMat, and then creates a fresh LinopMat of size r_ x c_...
Definition: sis.hpp:8214