EVPOrrSommerfeldSparse_slepcz.cpp File Reference

Solves the following linear eigenvalue problem for values that satisfy : More...

#include <Generic_bundle.h>
#include <SparseLinearEigenSystem.h>
#include <SlepcSession.h>

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Detailed Description

Solves the following linear eigenvalue problem for values that satisfy :

subject to where , and . The matrix problem is constructed manually in this case, using second-order finite differences. A sparse eigenvalue routine is employed, via SLEPc. These values approximately correspond to the first neutral temporal mode in plane Poiseuille flow, therefore the test to be satisfied is that an eigenvalue exists with .

Definition in file EVPOrrSommerfeldSparse_slepcz.cpp.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 24 of file EVPOrrSommerfeldSparse_slepcz.cpp.

{
  SlepcSession::getInstance(argc,argv);
 
  cout << "\n";
  cout << "=== EVP: Temporal spectra of the Orr-Sommerfeld eqn =\n";
  cout << "===  with a matrix problem assembled by hand and \n";
  cout << "===  eigenproblem solved with SLEPc sparse solver.\n";
  cout << "\n";
 
  // discretise with these many nodal points
  const std::size_t nodes( 1001 );
  // we'll solve as TWO second order problems
  const std::size_t N( 2 * nodes );
  // domain boundaries
  const double left = -1.0;
  const double right = 1.0;
  // spatial step for a uniform mesh
  const double d = ( right - left ) / ( nodes - 1 );
 
  // matrices for the EVP, initialised with zeroes
  SparseMatrix<D_complex> a( N, N );
  SparseMatrix<D_complex> b( N, N );
 
  // streamwise wavenumber and Reynolds number
  const double alpha ( 1.02 );
  const double Re ( 5772.2 );
  const D_complex I( 0.0, 1.0 );
 
  // boundary conditions at the left boundary
  a( 0, 0 ) = 1.0;           // phi( left ) = 0
  a( 1, 0 ) = -1.5 / d;      // phi'( left ) = 0
  a( 1, 2 ) = 2.0 / d;
  a( 1, 4 ) = -0.5 / d;
  // fill the interior nodes
  for ( std::size_t i = 1; i <= nodes - 2; ++i )
  {
    // position in the channel
    const double y = left + i * d;
    // Poiseuille flow profile
    const double U = ( 1.0 - y * y );
    const double Udd = -2.0;
 
    // the first quation at the i'th nodal point
    std::size_t row = 2 * i;
    a( row, row ) = -2.0 / ( d * d ) - alpha * alpha;
    a( row, row - 2 ) = 1.0 / ( d * d );
    a( row, row + 2 ) = 1.0 / ( d * d );
    a( row, row + 1 ) = -1.0;
 
    row += 1;
    // the second equation at the i'th nodal point
    a( row, row ) = -2.0 / ( d * d ) - alpha * alpha - I * alpha * Re * U;
    a( row, row - 2 ) = 1.0 / ( d * d );
    a( row, row + 2 ) = 1.0 / ( d * d );
    a( row, row - 1 ) = I * alpha * Re * Udd;
 
    b( row, row ) = - I * alpha * Re;
  }
  // boundary conditions at right boundary
  a( N - 2, N - 2 ) = 1.5 / d;
  a( N - 2, N - 4 ) = -2.0 / d;
  a( N - 2, N - 6 ) = 0.5 / d; // psi'( right ) = 0
  a( N - 1, N - 2 ) = 1.0;     // psi( right ) = 0
 
  /* Note: current KSP configuration uses an LU preconditioner and MUMPS
    which doesn't mind zeros on the diagonal of a. */
 
 
  // a vector for storing the eigenvalues
  DenseVector<D_complex> lambdas;
  SparseLinearEigenSystem<D_complex> system( &a, &b );
  system.set_target(D_complex(0.2,0.0));
  system.set_nev(4);
  system.set_order( "EPS_TARGET_MAGNITUDE" );
 
  try
  {
    system.eigensolve();
  }
  catch (const std::runtime_error &error )
  {
    cout << " \033[1;31;48m  * FAILED THROUGH EXCEPTION BEING RAISED \033[0m\n";
    return 1;
  }
  // tag any eigenvalues with imaginary part > -0.1
  system.set_shift( D_complex( 0.0, -0.1 ) );
  system.tag_eigenvalues_upper( + 1 );
  lambdas = system.get_tagged_eigenvalues();
  lambdas.dump();
  double min_growth_rate( lambdas[ 0 ].imag() );
  // make sure we have a near neutral mode
  const double tol = 1.e-3;
 
  std::string dirname("./DATA");
  mkdir( dirname.c_str(), S_IRWXU );
 
  TrackerFile spectrum( "./DATA/spectrum.dat" );
  spectrum.push_ptr( &lambdas, "evs" );
  spectrum.update();
 
  // eigenfunctions
  DenseMatrix<D_complex> eigenvectors;
  eigenvectors = system.get_tagged_eigenvectors();
  DenseVector<double> ynodes = Utility::uniform_node_vector(left,right,nodes);
  OneD_Node_Mesh<D_complex> mesh( ynodes, 4 );
  // non-zero eigenvalue list => dump index 0 mode
  if ( lambdas.nelts() > 0 ) {    
    int index(0);
    for ( unsigned i = 0; i < nodes; ++i ) {
      mesh(i,0) = eigenvectors(index,2*i+0); //phi_i
      mesh(i,1) = eigenvectors(index,2*i+1); //psi_i
    }
    for ( unsigned i = 1; i < nodes-1; ++i ) {
      mesh(i,2) = (mesh(i+1,0)-mesh(i-1,0))/(2*d); //u_i=phi'
      mesh(i,3) = -I*alpha*mesh(i,0);              //v_i=-i*alpha*phi
    }
  }
  mesh.normalise(0);
  mesh.dump_gnu("./DATA/eigenmode.dat");
  
  if ( std::abs( min_growth_rate ) < tol )
  {
    cout << "\033[1;32;48m  * PASSED \033[0m\n";
    return 0;
  }
 
  cout << "\033[1;31;48m  * FAILED \033[0m\n";
  cout << "    Final error = " << min_growth_rate << "\n";
  return 1;
  
}

References CppNoddy::DenseVector< _Type >::dump(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::dump_gnu(), CppNoddy::DenseVector< _Type >::nelts(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::normalise(), CppNoddy::TrackerFile::push_ptr(), U, CppNoddy::Utility::uniform_node_vector(), and CppNoddy::TrackerFile::update().