EVPHarmonicComplexSparse_slepcz.cpp File Reference

#include <Utility.h>
#include <EVP_bundle.h>
#include <SparseLinearEigenSystem.h>
#include <SlepcSession.h>

Go to the source code of this file.

Namespaces
namespace	CppNoddy
	A collection of OO numerical routines aimed at simple (typical) applied problems in continuum mechanics.
namespace	CppNoddy::Example

Functions
const D_complex	CppNoddy::Example::eye (0., 1.0)
const double	CppNoddy::Example::delta (0.5)
D_complex	CppNoddy::Example::z (double x)
D_complex	CppNoddy::Example::zx (double x)
D_complex	CppNoddy::Example::zxx (double x)
int	main (int argc, char *argv[])

Function Documentation

main()

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

Definition at line 49 of file EVPHarmonicComplexSparse_slepcz.cpp.

{
 SlepcSession::getInstance(argc,argv);
  
  cout << "\n";
  cout << "=== EVP: Harmonic equation solved using LAPACK  =====\n";
  cout << "===  with a manually assembled matrix problem.\n";
  cout << "===  The problem is solved along a path in the complex\n";
  cout << "===  plane, deformed away from the real axis.\n";
  cout << "\n";
 
  cout.precision( 12 );
  cout << " Number of nodal points : Leading eigenvalue error : Total CPU time taken (ms) \n";
  bool failed = false;
  size_t N = 4;
  // a vector for the eigenvalues
  DenseVector<D_complex> lambdas;
  DenseMatrix<D_complex> eigenvectors;
  
  Timer timer;
  for ( int i = 2; i < 11; ++i )
  {
    N = ( size_t ) ( std::pow( 2., i ) );
    const double delta = 1. / ( N - 1 );
    const double delta2 = delta * delta;
    // matrix problem
    SparseMatrix<D_complex> a( N, N);
    SparseMatrix<D_complex> b( N, N);
    // boundary conditions at f(0) = 0
    a( 0, 0 ) = 1.0;
    a( 0, 1 ) = 0.0;
    for ( unsigned j = 1; j < N-1; ++j ) {
      // Finite difference representation of f''(x)
      double x = j*delta;
      a( j, j-1 ) =  (1.0/pow(Example::zx(x),2.0))/delta2 + (Example::zxx(x)/pow(Example::zx(x),3.0))/(2*delta);
      a( j, j)    = -(2.0/pow(Example::zx(x),2.0))/delta2;
      a( j, j+1 ) =  (1.0/pow(Example::zx(x),2.0))/delta2 - (Example::zxx(x)/pow(Example::zx(x),3.0))/(2*delta);
 
      // not a generalised problem - but we'll apply that routine anyway
      b( j, j )   =  -1.0;
    }
    // boundary conditions at f(1) = 0
    a( N - 1, N - 1 ) = 1.0;
    a( N - 1, N - 2 ) = 0.0;
    // a vector for the eigenvectors - although we won't use them
    SparseLinearEigenSystem<D_complex> system( &a, &b );
    system.set_calc_eigenvectors( true );
    system.set_nev(4);
    system.set_order( "EPS_TARGET_MAGNITUDE" ); //default
    
    timer.start();
    try
    {
      system.eigensolve();
    } catch (const std::runtime_error &error ) {
      cout << " \033[1;31;48m  * FAILED THROUGH EXCEPTION BEING RAISED \033[0m\n";
      return 1;
    }
    
    system.tag_eigenvalues_disc( + 1, 10. );
    lambdas = system.get_tagged_eigenvalues();
    
    eigenvectors = system.get_tagged_eigenvectors();
 
    cout << "    " << N << " : " << lambdas[ 0 ].real() - M_PI * M_PI << " +i " << lambdas[0].imag() 
         << " : " << timer.get_time() << "\n";
    timer.stop();
  }
 
  // real parameterisation of complex path
  DenseVector<double> xNodes( Utility::uniform_node_vector( 0.0, 1.0, N ) );
  // complex path
  DenseVector<D_complex> zNodes( N, D_complex(0.0,0.0) );
  // mesh of complex data along a complex path
  OneD_Node_Mesh<D_complex, D_complex> mesh( zNodes, 1 );
  int index(0); // first and only eigenvalue returned
  for ( unsigned j = 0; j < N; ++j ) {
    // store the nodes in the complex plane
    mesh.coord(j) = Example::z(xNodes[j]);
    mesh(j,0) = eigenvectors(index,j+0); //f_j
  }
  // write the complex solution on the complex path
  mesh.dump_gnu("/DATA/complexMesh.dat"); 
  const double tol = 1.e-4;
  if ( abs( lambdas[ 0 ].real() - M_PI * M_PI ) > tol )
    failed = true;
 
  if ( failed )
  {
    cout << "\033[1;31;48m  * FAILED \033[0m\n";
    cout << "    Final error = " << abs( lambdas[ 0 ].real() - M_PI * M_PI ) << "\n";
    return 1;
  }
 
  cout << "\033[1;32;48m  * PASSED \033[0m\n";
  return 0;
}

References CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::coord(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::dump_gnu(), CppNoddy::Timer::get_time(), CppNoddy::Timer::start(), CppNoddy::Timer::stop(), and CppNoddy::Utility::uniform_node_vector().