eigen_Laplacian_2D.f90

Go to the documentation of this file.

program eigen_laplacian_2d

  use choral
  use choral_constants

  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!      VARAIBLE DEFINITION : BEGIN
  !!
  !!  

  !! GENERAL SETTINGS  
  !!
  !! verb        = verbosity level
  !! WITH_MUMPS  = use MUMPS direct solver for linear system inversions
  !!
  integer, parameter :: VERB        = 2
  logical, parameter :: WITH_MUMPS  = .false.

  !! Arnoldi algorithm parameters
  !! nev    = number of computed eigenvalues
  !! size   = problem size (number of DOF)
  !! sigma  = shift factor
  !! tol    = tolerance for resolution   
  !! itMax  = maximal number of iterations
  !!
  !! lambda = exact eigenvalues (analytic)
  !!
  type(eigen) :: eig
  !!
  integer , parameter :: nev   = 6
  integer , parameter :: itMax = 100
  real(RP), parameter :: tol   = 1e-7_rp
  real(RP), parameter :: sigma =-1.0_rp 
  integer             :: size
  !!
  integer , parameter            :: nev_exact = 6
  real(RP), dimension(nev_exact) :: lambda


  !!      SPACE DISCRETISATION
  !!
  !!   fe     = finite element method 
  !!   quad   = quadrature method 
  !!
  integer, parameter :: fe   = fe_p2_2d
  integer, parameter :: quad = quad_gauss_trg_12
  !!
  !! msh      = mesh
  !! X_h      = finite element space
  !! qdm      = integration method
  !!
  type(mesh)     :: msh
  type(fespace)  :: X_h
  type(quadmesh) :: qdm
  !!
  !! msh_file = mesh file
  character(len=100), parameter :: &
       & mesh_file= trim(GMSH_DIR)//"square/square_3.msh"

  
  !!       LINEAR SYSTEM
  !!
  !!  M  = mass matrix
  !!  stiff = stiffness matrix
  !!  K     = matrix S - sigma*M, sigma the shift factor
  !!
  !!  lu  = mumps decomposition for K   ( if WITH_MUMPS )
  !!  kry = krylov method def.          ( else )
  !!  pc  = preconditioning for the CG krylov
  !!
  !!  Kinv = pointer to the system inversion K x = y
  !!
  type(csr)     :: M, S, K
  type(mumps)   :: LU
  type(krylov)  :: kry
  type(precond) :: pc
  !!
  procedure(rntornxl), pointer :: Kinv => null()

  !!   OTHER VARS
  !!
  real(RP)  :: err
  integer   :: ii
  character(LEN=300) :: shell_com

  !!
  !!   VARAIBLE DEFINITION : END
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  !! !!!!!!!!!!!!!!!!!!!!!  INITIALISATION
  !!
  call choral_init(verb=verb)
  write(*,*)
  write(*,*)'eigen_Laplacian_2D: start'
  write(*,*)""
  write(*,*)"    RESOLUTION OF THE EIGEN-PROBLEM"
  write(*,*)""
  write(*,*)"           -Delta u = lambda u   on  \Omega"
  write(*,*)""
  write(*,*)"    WITH A HOMOGENEOUS NEUMANN BOUNDARY CONDITION "  
  write(*,*)""
  write(*,*)"            grad u . n  = 0   on   \partial\Omega"
  write(*,*)""

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE MESH "
  msh = mesh(mesh_file, 'gmsh')
  if (verb>1) call print(msh)

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE&
       & FINITE ELEMENT SPACE"
  x_h = fespace(msh)
  call set(x_h, fe)
  call assemble(x_h)
  if (verb>1) call print(x_h)

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE&
       & INTEGRATION METHOD ON \Omega"
  qdm = quadmesh(msh)
  call set(qdm, quad)
  call assemble(qdm)
  if (verb>1) call print(qdm)

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE&
       & MASS AND STIFFNESS MATRICES"
  !!
  !! 'one_R3' is the density function equal to 1 on \Omega
  call diffusion_massmat(m  , one_r3, x_h, qdm)
  !!
  !! 'EMetric' is the Euclidian metric on \Omega
  call diffusion_stiffmat(s, emetric, x_h, qdm)

  write(*,*) ""
  write(*,*)"================================== COMPUTE THE SYSTEM &
       &MATRIX K = S - sigma * M"
  !!
  call add(k, s, 1._rp, m, -sigma)

  write(*,*) ""
  write(*,*)"================================== LINEAR SOLVER K*x = y"
  !!
  if (with_mumps) then
     write(*,*) " K = LDLT,  decomposition with MUMPS"
     lu = mumps(k, mumps_ldlt_sdp, verb = 0)

     kinv => kinv_mumps

  else
     write(*,*) " Krylov CG settings"
     kry = krylov(kry_cg, tol=tol*1.0e-2_rp, itmax=1000, verb=0)

     pc = precond(k, pc_jacobi)

     kinv => kinv_pcg
     
  end if


  write(*,*) ""
  write(*,*)"================================== EIGEN-SOLVER DEF."
  !!
  !! create the eigen-solver 'eig'
  size = x_h%nbDof
  eig  = eigen(size = size, nev = nev, sym=.true., std=.false.)
  !!
  !! set the eigen-solver 'eig' 
  call set(eig, itmax=itmax, tol=tol, verb=2)
  call set(eig, which = eig_which_sm)
  call set(eig, mode = eig_shift_invert, shift = sigma)
  if (verb>1) call print(eig)
  write(*,*) ""
 

  write(*,*) ""
  write(*,*)"================================== NUMERICAL RESOLUTION"
  !!
  !! Numerical resolution
  call solve(eig, b=pmv_m, op=kinv)
  !!
  call residual()

  write(*,*) ""
  write(*,*)"================================== POST-PROCESSING"
  print*
  write(*,*) "Computed eigenvalues"
  do ii=1, nev
     write(*,*) " i =", int(ii,1), ", lambda = ", &
          & real(eig%lambda(ii)/PI**2, SP), "* Pi**2"
  end do

  write(*,*) ""
  write(*,*) "Numerical errors |lambda - lambda_h| "
  !!
  !! Exact eigen values
  !!
  lambda(1) = 0._rp
  lambda(2) = pi**2
  lambda(3) = pi**2
  lambda(4) = 2.0_rp * pi**2
  lambda(5) = 4.0_rp * pi**2
  lambda(6) = 5.0_rp * pi**2
  !!
  do ii=1, min(nev, nev_exact)
     err = abs( eig%lambda(ii) - lambda(ii)) 
     write(*,*) " i =", int(ii,1), ", error = ", &
          & real(err, sp)
  end do

  if (verb>1) then
     !!
     !! graphical display
     !!
     !! output file def.
     shell_com = 'eig_Laplace.msh'
     !!
     !! first write  the 'finite element mesh'
     call write(x_h, trim(shell_com), 'gmsh')
     !!
     !! second save the eigenvectors u_h
     !! each eigenvector is a finite element function
     do ii=1, nev
        call gmsh_addview(x_h, eig%v(:,ii), &
             & trim(shell_com), &
             & 'Eigen-functions', real(ii, RP), ii)
     end do
     !!
     !! shell command to visualise with gmsh
     shell_com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'&
          & //' '//trim(shell_com)
     call system(shell_com)

  end if

contains 

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  subroutine pmv_m(y, x)
    real(RP), dimension(:), intent(out) :: y
    real(RP), dimension(:), intent(in)  :: x
    
    call matvecprod(y, m, x)

  end subroutine pmv_m

  !! Solve K x = y
  !!
  subroutine kinv_pcg(y, ierr, x)
    real(RP), dimension(:), intent(inout) :: y
    logical               , intent(out)   :: ierr
    real(RP), dimension(:), intent(in)    :: x

    call solve(y, kry, x, k)

    ierr = kry%ierr

  end subroutine kinv_pcg

  !! Solve LU*x = y, LU = LDLT MUMPS decomposition 
  !!                      of K = (S1 + S2 - sigma*M2) 
  !!
  subroutine kinv_mumps(x, ierr, y)
    real(RP), dimension(:), intent(inout) :: x
    logical               , intent(out)   :: ierr
    real(RP), dimension(:), intent(in)    :: y

    call solve(x, lu, y)
    ierr = .false.
  end subroutine kinv_mumps

  !! Compute and display the residuals on the
  !! computed eigenvalues and eigenvectors
  !!
  subroutine residual()
    real(RP), dimension(size) :: u, v, w
    real(RP) :: res, nu
    logical  :: ierr

    write(*,*)
    write(*,*) "Residual | Ku - lambda Mu |_l2 / | u |_l2"
    !!
    do ii=1, nev

       u = eig%v(:,ii)

       call matvecprod(v, k , u)
       call pmv_m(w, u)

       w = v - eig%lambda(ii)*w
       res = sqrt( sum(w*w) / sum(u*u) )

       print*, "  i = ", int(ii,1), " : res = ", real(res, sp)

    end do

    write(*,*)
    write(*,*) "Residual | \nu u -  K^{-1} M u |_l2 / | u |_l2"
    !!
    do ii=1, nev
       u = eig%v(:,ii)

       call pmv_m(v , u)
       call kinv(w, ierr, v)


       nu = 1.0_rp / ( eig%lambda(ii) - sigma )
       w = w - nu*u
       res = sqrt( sum(w*w) / sum(u*u) )

       print*, "  i = ", int(ii,1), " : res = ", real(res, sp)

    end do

  end subroutine residual


end program eigen_laplacian_2d