CppNoddy  0.92
Loading...
Searching...
No Matches
EVP2DOrrSommerfeld_slepcz_superludist.cpp
Go to the documentation of this file.
1/// \file EVP2DOrrSommerfeld_slepcz.cpp
2/// \ingroup Tests
3/// \ingroup EVP
4/// Solves the linear eigenvalue problem for values \f$ c \f$ formulated
5/// in the paper of Tatsumi and Yoshimura (J Fluid Mechanics, 1990).
6/// The macros defined below determine the symmetry of the sought eigenmode
7/// again using the classification introduced by Tatsumi and Yoshimura.
8/// The test is for the existence of a near (within 1.e-4) neutral wave
9/// for a duct flow of aspect ratio 6.
10
11
12#include <cassert>
13
14#include <Utility.h>
15#include <SparseMatrix.h>
16#include <SparseLinearEigenSystem.h>
17#include <TwoD_Mapped_Node_Mesh.h>
18#include <TwoD_Node_Mesh.h>
19#include <SlepcSession.h>
20#include <TrackerFile.h>
21
22//#define DENSE
23#define UNIFORM
24#define TYPEI
25
26#if defined(TYPEI)
27 #define V_EVEN_Y
28 #define V_EVEN_Z
29#endif
30#if defined(TYPEII)
31 #define V_EVEN_Y
32 #define V_ODD_Z
33#endif
34#if defined(TYPEIII)
35 #define V_ODD_Y
36 #define V_EVEN_Z
37#endif
38#if defined(TYPEIV)
39 #define V_ODD_Y
40 #define V_ODD_Z
41#endif
42
43#if defined(V_ODD_Y)
44 #define W_EVEN_Y
45#elif defined(V_EVEN_Y)
46 #define W_ODD_Y
47#endif
48
49#if defined(V_EVEN_Z)
50 #define W_ODD_Z
51#elif defined(V_ODD_Z)
52 #define W_EVEN_Z
53#endif
54
55using namespace CppNoddy;
56using namespace std;
57
58namespace Example
59{
60
61 // instanitate a mapped mesh
62 template <typename _Type>
63 class My_Mesh : public TwoD_Mapped_Node_Mesh<_Type>
64 {
65 public:
66 My_Mesh( double left, double right, double bottom, double top, std::size_t nz, std::size_t ny, std::size_t nv ) : TwoD_Mapped_Node_Mesh<_Type>( left, right, bottom, top, nz, ny, nv )
67 {
68 // mesh stretching parameters
69 BX = 10.0;
70 BY = 10.0;
71 CX = 0.2;
72 CY = 0.8;
73 this -> init_mapping();
74 }
75
76#if defined(UNIFORM)
77 // the computational x coord as a function of physical x
78 double FnComp_X( const double& x ) const
79 {
80 return x;
81 }
82 double FnComp_Xd( const double& x ) const
83 {
84 return 1.0;
85 }
86 double FnComp_Xdd( const double& x ) const
87 {
88 return 0.0;
89 }
90 // the computational y coord as a function of physical y
91 double FnComp_Y( const double& y ) const
92 {
93 return y;
94 }
95 double FnComp_Yd( const double& y ) const
96 {
97 return 1.0;
98 }
99 double FnComp_Ydd( const double& y ) const
100 {
101 return 0.0;
102 }
103#else
104 // the computational x coord as a function of physical x
105 double FnComp_X( const double& x ) const
106 {
107 return BX + x - BX*exp(-(x - this -> m_left)/CX );
108 }
109 double FnComp_Xd( const double& x ) const
110 {
111 return 1 + BX*exp(-(x - this -> m_left)/CX )/CX;
112 }
113 double FnComp_Xdd( const double& x ) const
114 {
115 return - BX*exp(-(x - this -> m_left)/CX )/(CX*CX);
116 }
117 // the computational y coord as a function of physical y
118 double FnComp_Y( const double& y ) const
119 {
120 return BY + y - BY*exp(-(y - this -> m_bottom)/CY );
121 }
122 double FnComp_Yd( const double& y ) const
123 {
124 return 1 + BY*exp(-(y - this -> m_bottom)/CY )/CY;
125 }
126 double FnComp_Ydd( const double& y ) const
127 {
128 return - BY*exp(-(y - this -> m_bottom)/CY )/(CY*CY);
129 }
130#endif
131
132 private:
133 double BX,CX,BY,CY;
134 };
135
136}
137
138
139int main( int argc, char *argv[] )
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 // duct aspect ratio
152 const double A(6.0);
153 // streamwise wavenumber and Reynolds number
154 const double alpha ( 0.94 );
155 double Re ( 8200.0 );
156 const D_complex eye( 0.0, 1.0 );
157
158 // discretise with these many nodal points
159 const std::size_t NZ( 301 );
160 const std::size_t NY( 301 );
161 //
162 //
163 Example::My_Mesh<double> base( -A, 0.0, -1.0, 0.0, NZ, NY, 1 );
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 // we'll solve as FOUR equations
185 const std::size_t N( 4*NZ*NY );
186
187 cout << "There are " << N << " degrees of freedom in this problem.\n";
188
189 // spatial step for the uniform computational-domain mesh
190 const double DZ = base.get_comp_step_sizes().first;
191 const double DY = base.get_comp_step_sizes().second;
192 //std::cout << " DZ = " << DZ << " DY = " << DY << "\n";
193
194 // matrices for the EVP, initialised with zeroes
195#ifdef DENSE
196 DenseMatrix<D_complex> a( N, N, 0.0 );
197 DenseMatrix<D_complex> b( N, N, 0.0 );
198#else
199 SparseMatrix<D_complex> a( N, N );
200 SparseMatrix<D_complex> b( N, N );
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 // left boundary conditions for the duct
209 const std::size_t O( i*4*NY );
210 std::size_t k( j*4 );
211 // coordinate mapping
212 double Z( base.coord(i,j).first );
213 double KZd( base.FnComp_Xd( Z ) );
214 double KZdd( base.FnComp_Xdd( Z ) );
215 // v = 0
216 a( O+k+0, O+k + v ) = 1.0;
217 // w = 0
218 a( O+k+1, O+k + w ) = 1.0;
219 // definition of F with w=0 for all y and 2nd order 1-sided diff.
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 // definition of G -- building in w_z = 0 with high order scheme
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 // bottom boundary conditions for the duct
237 std::size_t O( i*4*NY );
238 std::size_t k( j*4 );
239 // coordinate mapping
240 double Y( base.coord(i,j).second );
241 double KYd( base.FnComp_Yd( Y ) );
242 double KYdd( base.FnComp_Ydd( Y ) );
243 // v=0 on y=0
244 a( O+k+0, O+k + v ) = 1.0;
245 // w=0 on y=0
246 a( O+k+1, O+k + w ) = 1.0;
247 // definition of F, building in v_y=0 with high order scheme
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 // definition of G with w=0 for all y and 2nd order 1-sided diff.
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 // coordinate mapping
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 // interior equations within the duct
278 std::size_t O( i*4*NY );
279 std::size_t k( j*4 );
280 //
281 // L{v}=F
282 //
283 // v_i,j
284 a( O+k+0, O+k + v ) = (eye/(alpha*Re))*( - 2.0*KYd*KYd/(DY*DY)
285 - 2.0*KZd*KZd/(DZ*DZ)
286 - alpha*alpha ) + base(i,j,0);
287 // v_i+1,j
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 // v_i-1,j
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 // v_i,j+1
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 // v_i,j-1
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 // F_i,j
300 a( O+k+0, O+k + F ) = -1.0;
301 // ev -- wave speed c
302 b( O+k+0, O+k + v ) = 1.0;
303 //
304 // L{w}=G
305 //
306 // w_i,j
307 a( O+k+1, O+k + w ) = (eye/(alpha*Re))*( - 2.0*KYd*KYd/(DY*DY)
308 - 2.0*KZd*KZd/(DZ*DZ)
309 - alpha*alpha ) + base(i,j,0);
310 // w_i+1,j
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 // w_i-1,j
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 // w_i,j+1
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 // w_i,j-1
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 // G_i,j
323 a( O+k+1, O+k + G ) = -1.0;
324 // ev -- wave speed c
325 b( O+k+1, O+k + w ) = 1.0;
326 //
327 // E(y,z){v} = O(y,z){w}
328 //
329 // F_i,j
330 a( O+k+2, O+k + F ) = -2.0*KYd*KYd/(DY*DY) - alpha*alpha;
331 // F_i,j+1
332 a( O+k+2, O+k + 4 + F ) = 1.0*KYd*KYd/(DY*DY) + KYdd/(2.0*DY);
333 // F_i,j-1
334 a( O+k+2, O+k - 4 + F ) = 1.0*KYd*KYd/(DY*DY) - KYdd/(2.0*DY);
335 //
336 // G_i+1,j+1
337 a( O+k+2, O+k + dof_per_i + 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
338 // G_i+1,j-1
339 a( O+k+2, O+k + dof_per_i - 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
340 // G_i-1,j-1
341 a( O+k+2, O+k - dof_per_i - 4 + G ) = 0.25*KYd*KZd/(DY*DZ);
342 // G_i-1,j+1
343 a( O+k+2, O+k - dof_per_i + 4 + G ) = -0.25*KYd*KZd/(DY*DZ);
344 //
345 // v_i,j
346 a( O+k+2, O+k + v ) = -2.0*Uyy;
347 // v_i,j+1
348 a( O+k+2, O+k + 4 + v ) = -2.0*Uy*KYd/(2.0*DY);
349 // v_i,j-1
350 a( O+k+2, O+k - 4 + v ) = +2.0*Uy*KYd/(2.0*DY);
351 //
352 // w_i,j
353 a( O+k+2, O+k + w ) = -2.0*Uyz;
354 // w_i,j+1
355 a( O+k+2, O+k + 4 + w ) = -2.0*Uz*KYd/(2.0*DY);
356 // w_i,j-1
357 a( O+k+2, O+k - 4 + w ) = +2.0*Uz*KYd/(2.0*DY);
358 //
359 //
360 // E(z,y){w} = O(z,y){v}
361 //
362 // G_i,j
363 a( O+k+3, O+k + G ) = -2.0*KZd*KZd/(DZ*DZ) - alpha*alpha;
364 // G_i+1,j
365 a( O+k+3, O+k + dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) + KZdd/(2.0*DZ);
366 // G_i-1,j
367 a( O+k+3, O+k - dof_per_i + G ) = 1.0*KZd*KZd/(DZ*DZ) - KZdd/(2.0*DZ);
368 //
369 // F_i+1,j+1
370 a( O+k+3, O+k + dof_per_i + 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
371 // F_i+1,j-1
372 a( O+k+3, O+k + dof_per_i - 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
373 // F_i-1,j-1
374 a( O+k+3, O+k - dof_per_i - 4 + F ) = 0.25*KYd*KZd/(DY*DZ);
375 // F_i-1,j+1
376 a( O+k+3, O+k - dof_per_i + 4 + F ) = -0.25*KYd*KZd/(DY*DZ);
377 //
378 // v_i,j
379 a( O+k+3, O+k + v ) = -2.0*Uyz;
380 // v_i+1,j
381 a( O+k+3, O+k + dof_per_i + v ) = -2.0*Uy*KZd/(2.0*DZ);
382 // v_i-1,j
383 a( O+k+3, O+k - dof_per_i + v ) = +2.0*Uy*KZd/(2.0*DZ);
384 //
385 // w_i,j
386 a( O+k+3, O+k + w ) = -2.0*Uzz;
387 // w_i+1,j
388 a( O+k+3, O+k + dof_per_i + w ) = -2.0*Uz*KZd/(2.0*DZ);
389 // w_i-1,j
390 a( O+k+3, O+k - dof_per_i + w ) = +2.0*Uz*KZd/(2.0*DZ);
391 }
392 {
393 // top boundary
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 // coordinate mapping
399 double Y( base.coord(i,j).second );
400 double KYd( base.FnComp_Yd( Y ) );
401 //double KYdd( base.FnComp_Ydd( Y ) );
402 //
403#ifdef V_EVEN_Y
404 // v_y = 0
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 // v_yyy = 0 => F_y = 0
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 // v=0
414 a( O+k+0, O+k + v ) = 1.0;
415 // v_yy =0 => F=0
416 a( O+k+2, O+k + F ) = 1.0;
417#endif
418
419#ifdef W_EVEN_Y
420 // w_y = 0
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 // w_yyy = 0 => G_y = 0
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 // w = 0
430 a( O+k+1, O+k + w ) = 1.0;
431 // G = 0
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 // right hand boundary
440 const std::size_t i(NZ-1);
441 // left boundary conditions for the duct
442 const std::size_t O( i*4*NY );
443 std::size_t k( j*4 );
444 // coordinate mapping
445 double Z( base.coord(i,j).first );
446 double KZd( base.FnComp_Xd( Z ) );
447 //double KZdd( base.FnComp_Xdd( Z ) );
448#ifdef V_EVEN_Z
449 // v_z = 0
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 // v_zzz = 0 => F_z = 0
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 // v=0
459 a( O+k+0, O+k + v ) = 1.0;
460 // v_zz =0 => F=0
461 a( O+k+2, O+k + F ) = 1.0;
462#endif
463
464#ifdef W_EVEN_Z
465 // w_z = 0
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 // w_zzz = 0 => G_z = 0
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 // w = 0
475 a( O+k+1, O+k + w ) = 1.0;
476 // w_zz = 0 => G = 0
477 a( O+k+3, O+k + G ) = 1.0;
478#endif
479 }
480 }
481
482 // a vector for the eigenvalues
483 DenseVector<D_complex> ev_c;
484#ifdef DENSE
485 DenseLinearEigenSystem<D_complex> system( &a, &b );
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();
501 ev_c.dump();
502
503 DenseMatrix<D_complex> vec_mtx = system.get_tagged_eigenvectors();
504 Example::My_Mesh<D_complex> eigfn( -A, 0.0, -1.0, 0.0, NZ, NY, 4 );
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 }
518 eigfn.dump_gnu("./DATA/eigfn_"+Utility::stringify(n,4)+".dat");
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}
main
int main()
Definition: ArcCircle.cpp:39
SparseLinearEigenSystem.h
Specification of the sparse linear eigensystem class.
SparseMatrix.h
A matrix class that constructs a SPARSE matrix as an STL Vector of SparseVectors, inheriting from Mat...
TrackerFile.h
A class that can be passed pointers to scalar/vector/mesh objects, then their contents are written to...
TwoD_Mapped_Node_Mesh.h
A specification for a two dimensional mapped mesh object.
TwoD_Node_Mesh.h
A specification for a two dimensional mesh object.
Utility.h
A spec for a collection of utility functions.
CppNoddy::DenseLinearEigenSystem
A linear Nth-order generalised eigensystem class.
Definition: DenseLinearEigenSystem.h:19
CppNoddy::DenseLinearEigenSystem::tag_eigenvalues_disc
void tag_eigenvalues_disc(const int &val, const double &radius)
Tag those eigenvalues that are within a disc centred at a point in the complex plane.
Definition: DenseLinearEigenSystem.cpp:386
CppNoddy::DenseLinearEigenSystem::get_tagged_eigenvalues
DenseVector< D_complex > get_tagged_eigenvalues() const
Get the the tagged eigenvalues.
Definition: DenseLinearEigenSystem.cpp:337
CppNoddy::DenseLinearEigenSystem::eigensolve
void eigensolve()
Solve the matrix linear eigensystem.
Definition: DenseLinearEigenSystem.cpp:28
CppNoddy::DenseLinearEigenSystem::get_tagged_eigenvectors
DenseMatrix< D_complex > get_tagged_eigenvectors() const
Get the the tagged eigenvectors.
Definition: DenseLinearEigenSystem.cpp:359
CppNoddy::DenseMatrix
A matrix class that constructs a DENSE matrix as a row major std::vector of DenseVectors.
Definition: DenseMatrix.h:25
CppNoddy::DenseVector
An DenseVector class – a dense vector object.
Definition: DenseVector.h:34
CppNoddy::DenseVector::dump
void dump() const
Dump to std::cout.
Definition: DenseVector.cpp:64
CppNoddy::DenseVector::size
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:330
CppNoddy::SparseMatrix
A matrix class that constructs a SPARSE matrix as a row major std::vector of SparseVectors.
Definition: SparseMatrix.h:31
CppNoddy::TrackerFile::push_ptr
void push_ptr(double *scalar, std::string desc="")
Definition: TrackerFile.cpp:27
CppNoddy::TwoD_Mapped_Node_Mesh
A two dimensional (mapped) mesh utility object.
Definition: TwoD_Mapped_Node_Mesh.h:21
CppNoddy::TwoD_Mapped_Node_Mesh::m_left
double m_left
Definition: TwoD_Mapped_Node_Mesh.h:301
CppNoddy::TwoD_Mapped_Node_Mesh::get_comp_step_sizes
std::pair< double, double > get_comp_step_sizes() const
Definition: TwoD_Mapped_Node_Mesh.cpp:155
CppNoddy::TwoD_Mapped_Node_Mesh::init_mapping
void init_mapping()
Sometimes a useful mapping is painful to invert analytically.
Definition: TwoD_Mapped_Node_Mesh.cpp:75
CppNoddy::TwoD_Mapped_Node_Mesh::m_bottom
double m_bottom
Definition: TwoD_Mapped_Node_Mesh.h:301
CppNoddy::TwoD_Mapped_Node_Mesh::dump_gnu
void dump_gnu(std::string filename) const
A simple method for dumping data to a file for gnuplot surface plotting.
CppNoddy::TwoD_Mapped_Node_Mesh::coord
std::pair< double, double > coord(const std::size_t nodex, const std::size_t nodey) const
Access the physical nodal position - as a pair.
Definition: TwoD_Mapped_Node_Mesh.h:371
Example::My_Mesh
Definition: EVP2DOrrSommerfeld_slepcz_mumps.cpp:64
Example::My_Mesh::FnComp_Ydd
double FnComp_Ydd(const double &y) const
Mapping function that provides the second derivative of the computational Y coordinate as a function ...
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:99
Example::My_Mesh::My_Mesh
My_Mesh(double left, double right, double bottom, double top, std::size_t nz, std::size_t ny, std::size_t nv)
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:66
Example::My_Mesh::FnComp_Xd
double FnComp_Xd(const double &x) const
Mapping function that provides the first derivative of the computational m_X coordinate as a function...
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:82
Example::My_Mesh::FnComp_X
double FnComp_X(const double &x) const
Mapping function that provides a computational X coordinate from a physical coordinate.
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:78
Example::My_Mesh::FnComp_Xdd
double FnComp_Xdd(const double &x) const
Mapping function that provides the second derivative of the computational X coordinate as a function ...
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:86
Example::My_Mesh::FnComp_Yd
double FnComp_Yd(const double &y) const
Mapping function that provides the first derivative of the computational Y coordinate as a function o...
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:95
Example::My_Mesh::FnComp_Y
double FnComp_Y(const double &y) const
Mapping function that provides the computational Y coordinate as a function of the physical coordinat...
Definition: EVP2DOrrSommerfeld_slepcz_superludist.cpp:91
CppNoddy::Utility::stringify
std::string stringify(const int &val)
Return an integer value as a string - useful for file naming.
Definition: Utility.cpp:193
CppNoddy
A collection of OO numerical routines aimed at simple (typical) applied problems in continuum mechani...
Definition: ArcLength_base.h:13
CppNoddy::D_complex
std::complex< double > D_complex
A complex double precision number using std::complex.
Definition: Types.h:98
Generated on Tue Jan 1 1980 00:00:00 for CppNoddy by   doxygen 1.9.5

© 2012

R.E. Hewitt