140{
141 SlepcSession::getInstance(argc,argv);
142
143 cout << "\n";
144 cout << "=== EVP: Temporal spectra of the 2D Orr-Sommerfeld eqn =\n";
145 cout << "=== with a matrix problem assembled by hand using =\n";
146 cout << "=== the formulation of Tatsumi and Yoshimura (1990). =\n";
147 cout << "\n";
148
149 enum {v,w,F,G};
150
151
153
154 const double alpha ( 0.94 );
155 double Re ( 8200.0 );
157
158
159 const std::size_t NZ( 301 );
160 const std::size_t NY( 301 );
161
162
164
165 for ( std::size_t i = 0; i < NZ; ++i )
166 {
167 for ( std::size_t j = 0; j < NY; ++j )
168 {
169 double Z = base.coord(i,j).first;
170 double Y = base.coord(i,j).second;
171 double sum = 1.0 - Y*Y;
172 for ( std::size_t n = 0; n < 1000; ++n )
173 {
174 double C = (2*n+1)*M_PI/2.0;
175 double term = cos(C*Y)*(exp(C*(Z-A))+exp(-C*(Z+A)))/(1+exp(-2*C*A));
176 sum -= 4*std::pow(2./M_PI,3)*std::pow(-1.0,n)*term/std::pow(2*n+1,3);
177 }
178 base(i,j,0) = sum;
179 }
180 }
181 base.dump_gnu("./DATA/U.dat");
182
183
184
185 const std::size_t N( 4*NZ*NY );
186
187 cout << "There are " << N << " degrees of freedom in this problem.\n";
188
189
190 const double DZ = base.get_comp_step_sizes().first;
191 const double DY = base.get_comp_step_sizes().second;
192
193
194
195#ifdef DENSE
198#else
201#endif
202
203 std::size_t dof_per_i( 4*NY );
204
205 for ( std::size_t j = 0; j < NY; ++j )
206 {
207 const std::size_t i(0);
208
209 const std::size_t O( i*4*NY );
210 std::size_t k( j*4 );
211
212 double Z( base.coord(i,j).first );
213 double KZd( base.FnComp_Xd( Z ) );
214 double KZdd( base.FnComp_Xdd( Z ) );
215
216 a( O+k+0, O+k + v ) = 1.0;
217
218 a( O+k+1, O+k + w ) = 1.0;
219
220 a( O+k+2, O+k + dof_per_i + v ) = -5.0*(
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
221 + 4.0*(eye/(alpha*Re))*KZdd/(2.0*DZ);
222 a( O+k+2, O+k + 2*dof_per_i + v ) = 4.0*(
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
223 - 1.0*(eye/(alpha*Re))*KZdd/(2.0*DZ);
224 a( O+k+2, O+k + 3*dof_per_i + v ) = -1.0*(
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ);
225 a( O+k+2, O+k + F ) = -1.0;
226
227 a( O+k+3, O+k + dof_per_i + w ) = 10.0*(
eye/(
alpha*
Re))*KZd*KZd/(3*DZ*DZ);
228 a( O+k+3, O+k + 2*dof_per_i + w ) = -1.0*(
eye/(
alpha*
Re))*KZd*KZd/(3*DZ*DZ);
229 a( O+k+3, O+k + G ) = -1.0;
230 }
231
232 for ( std::size_t i = 1; i <= NZ - 2; ++i )
233 {
234 {
235 const std::size_t j(0);
236
237 std::size_t O( i*4*NY );
238 std::size_t k( j*4 );
239
240 double Y( base.coord(i,j).second );
241 double KYd( base.FnComp_Yd( Y ) );
242 double KYdd( base.FnComp_Ydd( Y ) );
243
244 a( O+k+0, O+k + v ) = 1.0;
245
246 a( O+k+1, O+k + w ) = 1.0;
247
248 a( O+k+2, O+k + 4 + v ) = 10.0*(
eye/(
alpha*
Re))*KYd*KYd/(3*DY*DY);
249 a( O+k+2, O+k + 8 + v ) = -1.0*(
eye/(
alpha*
Re))*KYd*KYd/(3*DY*DY);
250 a( O+k+2, O+k + F ) = -1.0;
251
252 a( O+k+3, O+k + 4 + w ) = -5.0*(
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
253 + 4.0*(eye/(alpha*Re))*KYdd/(2.0*DY);
254 a( O+k+3, O+k + 8 + w ) = 4.0*(
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
255 - 1.0*(eye/(alpha*Re))*KYdd/(2.0*DY);
256 a( O+k+3, O+k + 12 + w ) = -1.0*(
eye/(
alpha*
Re))*KYd*KYd/(DY*DY);
257 a( O+k+3, O+k + G ) = -1.0;
258 }
259 for ( std::size_t j = 1; j < NY-1; ++j )
260 {
261
262 double Z( base.coord(i,j).first );
263 double KZd( base.FnComp_Xd( Z ) );
264 double KZdd( base.FnComp_Xdd( Z ) );
265 double Y( base.coord(i,j).second );
266 double KYd( base.FnComp_Yd( Y ) );
267 double KYdd( base.FnComp_Ydd( Y ) );
268
269 double Uy = (base(i,j+1,0)-base(i,j-1,0))*KYd/(2*DY);
270 double Uz = (base(i+1,j,0)-base(i-1,j,0))*KZd/(2*DZ);
271 double Uyy = (base(i,j+1,0)-2*base(i,j,0)+base(i,j-1,0))*KYd*KYd/(DY*DY)
272 + (base(i,j+1,0)-base(i,j-1,0))*KYdd/(2*DY);
273 double Uzz = (base(i+1,j,0)-2*base(i,j,0)+base(i-1,j,0))*KZd*KZd/(DZ*DZ)
274 + (base(i+1,j,0)-base(i-1,j,0))*KZdd/(2*DZ);
275 double Uyz = (base(i+1,j+1,0)-base(i+1,j-1,0)
276 - base(i-1,j+1,0)+base(i-1,j-1,0))*KYd*KZd/(4*DY*DZ);
277
278 std::size_t O( i*4*NY );
279 std::size_t k( j*4 );
280
281
282
283
284 a( O+k+0, O+k + v ) = (
eye/(
alpha*
Re))*( - 2.0*KYd*KYd/(DY*DY)
285 - 2.0*KZd*KZd/(DZ*DZ)
287
288 a( O+k+0, O+k + dof_per_i + v ) = (
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
289 + (eye/(alpha*Re))*KZdd/(2.0*DZ);
290
291 a( O+k+0, O+k - dof_per_i + v ) = (
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
292 - (eye/(alpha*Re))*KZdd/(2.0*DZ);
293
294 a( O+k+0, O+k + 4 + v ) = (
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
295 + (eye/(alpha*Re))*KYdd/(2.0*DY);
296
297 a( O+k+0, O+k - 4 + v ) = (
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
298 - (eye/(alpha*Re))*KYdd/(2.0*DY);
299
300 a( O+k+0, O+k + F ) = -1.0;
301
302 b( O+k+0, O+k + v ) = 1.0;
303
304
305
306
307 a( O+k+1, O+k + w ) = (
eye/(
alpha*
Re))*( - 2.0*KYd*KYd/(DY*DY)
308 - 2.0*KZd*KZd/(DZ*DZ)
310
311 a( O+k+1, O+k + dof_per_i + w ) = (
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
312 + (eye/(alpha*Re))*KZdd/(2.0*DZ);
313
314 a( O+k+1, O+k - dof_per_i + w ) = (
eye/(
alpha*
Re))*KZd*KZd/(DZ*DZ)
315 - (eye/(alpha*Re))*KZdd/(2.0*DZ);
316
317 a( O+k+1, O+k + 4 + w ) = (
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
318 + (eye/(alpha*Re))*KYdd/(2.0*DY);
319
320 a( O+k+1, O+k - 4 + w ) = (
eye/(
alpha*
Re))*KYd*KYd/(DY*DY)
321 - (eye/(alpha*Re))*KYdd/(2.0*DY);
322
323 a( O+k+1, O+k + G ) = -1.0;
324
325 b( O+k+1, O+k + w ) = 1.0;
326
327
328
329
330 a( O+k+2, O+k + F ) = -2.0*KYd*KYd/(DY*DY) - alpha*alpha;
331
332 a( O+k+2, O+k + 4 + F ) = 1.0*KYd*KYd/(DY*DY) + KYdd/(2.0*DY);
333
334 a( O+k+2, O+k - 4 + F ) = 1.0*KYd*KYd/(DY*DY) - KYdd/(2.0*DY);
335
336
337 a( O+k+2, O+k + dof_per_i + 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
338
339 a( O+k+2, O+k + dof_per_i - 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
340
341 a( O+k+2, O+k - dof_per_i - 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
342
343 a( O+k+2, O+k - dof_per_i + 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
344
345
346 a( O+k+2, O+k + v ) = -2.0*Uyy;
347
348 a( O+k+2, O+k + 4 + v ) = -2.0*Uy*KYd/(2.0*DY);
349
350 a( O+k+2, O+k - 4 + v ) = +2.0*Uy*KYd/(2.0*DY);
351
352
353 a( O+k+2, O+k + w ) = -2.0*Uyz;
354
355 a( O+k+2, O+k + 4 + w ) = -2.0*Uz*KYd/(2.0*DY);
356
357 a( O+k+2, O+k - 4 + w ) = +2.0*Uz*KYd/(2.0*DY);
358
359
360
361
362
363 a( O+k+3, O+k + G ) = -2.0*KZd*KZd/(DZ*DZ) - alpha*alpha;
364
365 a( O+k+3, O+k + dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) + KZdd/(2.0*DZ);
366
367 a( O+k+3, O+k - dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) - KZdd/(2.0*DZ);
368
369
370 a( O+k+3, O+k + dof_per_i + 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
371
372 a( O+k+3, O+k + dof_per_i - 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
373
374 a( O+k+3, O+k - dof_per_i - 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
375
376 a( O+k+3, O+k - dof_per_i + 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
377
378
379 a( O+k+3, O+k + v ) = -2.0*Uyz;
380
381 a( O+k+3, O+k + dof_per_i + v ) = -2.0*Uy*KZd/(2.0*DZ);
382
383 a( O+k+3, O+k - dof_per_i + v ) = +2.0*Uy*KZd/(2.0*DZ);
384
385
386 a( O+k+3, O+k + w ) = -2.0*Uzz;
387
388 a( O+k+3, O+k + dof_per_i + w ) = -2.0*Uz*KZd/(2.0*DZ);
389
390 a( O+k+3, O+k - dof_per_i + w ) = +2.0*Uz*KZd/(2.0*DZ);
391 }
392 {
393
394
395 const std::size_t j(NY-1);
396 std::size_t O( i*4*NY );
397 const std::size_t k( j*4 );
398
399 double Y( base.coord(i,j).second );
400 double KYd( base.FnComp_Yd( Y ) );
401
402
403#ifdef V_EVEN_Y
404
405 a( O+k+0, O+k + v ) = 3.0*KYd/(2.0*DY);
406 a( O+k+0, O+k - 4 + v ) = -4.0*KYd/(2.0*DY);
407 a( O+k+0, O+k - 8 + v ) = 1.0*KYd/(2.0*DY);
408
409 a( O+k+2, O+k + F ) = 3.0*KYd/(2.0*DY);
410 a( O+k+2, O+k - 4 + F ) = -4.0*KYd/(2.0*DY);
411 a( O+k+2, O+k - 8 + F ) = 1.0*KYd/(2.0*DY);
412#elif defined(V_ODD_Y)
413
414 a( O+k+0, O+k + v ) = 1.0;
415
416 a( O+k+2, O+k + F ) = 1.0;
417#endif
418
419#ifdef W_EVEN_Y
420
421 a( O+k+1, O+k + w ) = 3.0*KYd/(2.0*DY);
422 a( O+k+1, O+k - 4 + w ) = -4.0*KYd/(2.0*DY);
423 a( O+k+1, O+k - 8 + w ) = 1.0*KYd/(2.0*DY);
424
425 a( O+k+3, O+k + G ) = 3.0*KYd/(2.0*DY);
426 a( O+k+3, O+k - 4 + G ) = -4.0*KYd/(2.0*DY);
427 a( O+k+3, O+k - 8 + G ) = 1.0*KYd/(2.0*DY);
428#elif defined(W_ODD_Y)
429
430 a( O+k+1, O+k + w ) = 1.0;
431
432 a( O+k+3, O+k + G ) = 1.0;
433#endif
434 }
435 }
436 {
437 for ( std::size_t j = 0; j < NY; ++j )
438 {
439
440 const std::size_t i(NZ-1);
441
442 const std::size_t O( i*4*NY );
443 std::size_t k( j*4 );
444
445 double Z( base.coord(i,j).first );
446 double KZd( base.FnComp_Xd( Z ) );
447
448#ifdef V_EVEN_Z
449
450 a( O+k+0, O+k + v ) = 3.0*KZd/(2.0*DZ);
451 a( O+k+0, O+k - dof_per_i + v ) = -4.0*KZd/(2.0*DZ);
452 a( O+k+0, O+k - 2*dof_per_i + v ) = 1.0*KZd/(2.0*DZ);
453
454 a( O+k+2, O+k + F ) = 3.0*KZd/(2.0*DZ);
455 a( O+k+2, O+k - dof_per_i + F ) = -4.0*KZd/(2.0*DZ);
456 a( O+k+2, O+k - 2*dof_per_i + F ) = 1.0*KZd/(2.0*DZ);
457#elif defined(V_ODD_Z)
458
459 a( O+k+0, O+k + v ) = 1.0;
460
461 a( O+k+2, O+k + F ) = 1.0;
462#endif
463
464#ifdef W_EVEN_Z
465
466 a( O+k+1, O+k + w ) = 3.0*KZd/(2.0*DZ);
467 a( O+k+1, O+k - dof_per_i + w ) = -4.0*KZd/(2.0*DZ);
468 a( O+k+1, O+k - 2*dof_per_i + w ) = 1.0*KZd/(2.0*DZ);
469
470 a( O+k+3, O+k + G ) = 3.0*KZd/(2.0*DZ);
471 a( O+k+3, O+k - dof_per_i + G ) = -4.0*KZd/(2.0*DZ);
472 a( O+k+3, O+k - 2*dof_per_i + G ) = 1.0*KZd/(2.0*DZ);
473#elif defined(W_ODD_Z)
474
475 a( O+k+1, O+k + w ) = 1.0;
476
477 a( O+k+3, O+k + G ) = 1.0;
478#endif
479 }
480 }
481
482
484#ifdef DENSE
486#else
487 SparseLinearEigenSystem<D_complex> system( &a, &b );
488 system.set_nev(1);
489 system.set_region(0.0, 1.0, -1.0, 1.0);
490 system.set_target(
D_complex( 0.244, 0.0 ) );
491 system.set_order("EPS_TARGET_IMAGINARY");
492#endif
493
494 system.eigensolve();
495
496 system.tag_eigenvalues_disc( + 1, 1.0 );
497 ev_c = system.get_tagged_eigenvalues();
498 TrackerFile spectrum(
"./DATA/spectrum_disc_2d.dat" );
499 spectrum.push_ptr( &ev_c, "evs" );
500 spectrum.update();
502
505
506 for ( std::size_t n = 0; n < ev_c.
size(); ++n )
507 {
508 for ( std::size_t i = 0; i < NZ; ++i )
509 {
510 for ( std::size_t j = 0; j < NY; ++j )
511 {
512 eigfn(i,j,v) = vec_mtx(n, i*4*NY + 4*j + v );
513 eigfn(i,j,w) = vec_mtx(n, i*4*NY + 4*j + w );
514 eigfn(i,j,F) = vec_mtx(n, i*4*NY + 4*j + F );
515 eigfn(i,j,G) = vec_mtx(n, i*4*NY + 4*j + G );
516 }
517 }
519 }
520
521 if ( std::abs( ev_c[0].
imag() ) < 1.e-4 )
522 {
523 cout << "\033[1;32;48m * PASSED \033[0m\n";
524 return 0;
525 }
526
527 cout << "\033[1;31;48m * FAILED \033[0m\n";
528 cout << " Final error = " << ev_c[0].imag() << "\n";
529 return 1;
530
531}
A linear Nth-order generalised eigensystem class.
A matrix class that constructs a DENSE matrix as a row major std::vector of DenseVectors.
An DenseVector class – a dense vector object.
void dump() const
Dump to std::cout.
std::size_t size() const
A pass-thru definition to get the size of the vector.
A matrix class that constructs a SPARSE matrix as a row major std::vector of SparseVectors.
double A(1.0)
initial hump amplitude
const D_complex eye(0., 1.0)
double Re
Globally define the Reynolds number and wavenumber.
DenseVector< double > imag(const DenseVector< D_complex > &X)
Return a double DENSE vector containing the imaginary part of a complex DENSE vector.
std::string stringify(const int &val)
Return an integer value as a string - useful for file naming.
std::complex< double > D_complex
A complex double precision number using std::complex.