EVPHarmonicComplex_lapack.cpp File Reference

Solves the harmonic equation. More...

#include <Utility.h>
#include <EVP_bundle.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 ()

Detailed Description

Solves the harmonic equation.

as an eigenvalue problem for over a path in the complex plane with homogeneous boundary conditions for , returning any eigenvalue(s) with absolute value less than 10. The resulting eigenvalues should approach as with corrections where and is the number of nodal points in the FD representation. The complex path is parametrised by a real parameter "x" and a 2nd order central difference representation is used. The matrix problem is solved for all eigenvalues within a specified distance of the origin of the comlplex plane by calling the LAPACK generalised eigenvalue routine.

Definition in file EVPHarmonicComplex_lapack.cpp.

Function Documentation

main()

int main

(

)

Definition at line 46 of file EVPHarmonicComplex_lapack.cpp.

{
  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
    DenseMatrix<D_complex> a( N, N, 0.0 );
    DenseMatrix<D_complex> b( N, N, 0.0 );
    // 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
    DenseLinearEigenSystem<D_complex> system( &a, &b );
    system.set_calc_eigenvectors( true );
 
    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::DenseLinearEigenSystem< _Type >::eigensolve(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvalues(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvectors(), CppNoddy::Timer::get_time(), CppNoddy::LinearEigenSystem_base::set_calc_eigenvectors(), CppNoddy::Timer::start(), CppNoddy::Timer::stop(), CppNoddy::DenseLinearEigenSystem< _Type >::tag_eigenvalues_disc(), and CppNoddy::Utility::uniform_node_vector().