eigen_Ventcel_2D.f90

Go to the documentation of this file.

program eigen_ventcel_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.


  !! eig    = eigen-problem solver
  !!
  !! 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   = 5
  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_1  = volumic  finite element method in Omega
  !!   fe_2  = surfacic finite element method on Gamma 
  !!
  !!   qd_1  = volumic  quadrature method in Omega
  !!   qd_2  = surfacic quadrature method on Gamma  
  !!
  integer, parameter :: fe_1 = fe_p2_2d
  integer, parameter :: qd_1 = quad_gauss_trg_12
  integer, parameter :: fe_2 = fe_p2_1d
  integer, parameter :: qd_2 = quad_gauss_edg_4
  !!
  !! msh      = mesh
  !! X_h      = finite element space
  !! qdm_1    = integration method on Omega
  !! qdm_2    = integration method on Gamma
  !!
  type(mesh)     :: msh
  type(fespace)  :: X_h
  type(quadmesh) :: qdm_1, qdm_2
  !!
  !! msh_file 
  character(len=100), parameter :: &
       & mesh_file= trim(GMSH_DIR)//"disk/disk2_3.msh"

  
  !!  LINEAR SYSTEM
  !!
  !!  S1  = stiffness matrix
  !!  S2  = stiffness matrix on the boundary
  !!  M2  = mass matrix on the boundary
  !!  K   = ( S1 + S2 )- sigma * M2, 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)     :: S1, S2, M2, 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=1)
  write(*,*)
  write(*,*)'eigen_Ventcel_2D: start'
  write(*,*)""
  write(*,*)"    RESOLUTION OF THE EIGEN-PROBLEM"
  write(*,*)""
  write(*,*)"     Delta u = 0   on  Omega"
  write(*,*)""
  write(*,*)"     Delta_B u + \partial_n u + lambda u = 0 on Gamma"  
  write(*,*)""
  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_1)   !! set the volumic  elements with fe_1
  call set(x_h, fe_2)   !! set the boundary elements with fe_2
  call assemble(x_h)
  if (verb>1) call print(x_h)

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

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE&
       & INTEGRATION METHOD ON Gamma"
  qdm_2 = quadmesh(msh)
  call set(qdm_2, qd_2) 
  call assemble(qdm_2)
  if (verb>1) call print(qdm_2)

  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE&
       & MATRICES S1, S2, M2"
  !!
  !! M2 = mass matrix on Gamma
  !! 'one_R3' is the density function equal to 1 on \Omega
  call diffusion_massmat(m2, one_r3, x_h, qdm_2)
  !!
  !! S1 = stiffness matrix on Omega
  !! S2 = stiffness matrix on Gamma
  !! 'EMetric' is the Euclidian metric on \Omega
  call diffusion_stiffmat(s1, emetric, x_h, qdm_1)
  call diffusion_stiffmat(s2, emetric, x_h, qdm_2)

  write(*,*) ""
  write(*,*)"================================== COMPUTE THE SYSTEM &
       &MATRIX K = S1 + S2 - sigma * M2"
  !!
  call copy(k, s1)
  call add_submatrix(k, s2)
  call add_submatrix(k, m2, -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(*,*)"================================== 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), SP)
  end do

  write(*,*) ""
  write(*,*) "Numerical errors |lambda - lambda_h| "
  !!
  !! Exact eigen values
  !!
  lambda(1) = 0._rp
  lambda(2) = 2.0_rp
  lambda(3) = 2.0_rp
  lambda(4) = 6.0_rp
  lambda(5) = 6.0_rp
  lambda(6) = 12.0_rp
  !!
  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_ventcel.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, m2, x)
  end subroutine pmv_m

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

    call solve(x, kry, y, k, pc)

    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_ventcel_2d