poisson_2d.f90

Go to the documentation of this file.

program poisson_2d

  use choral
  use choral_constants

  implicit none

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

  !! verb    = verbosity level
  !!
  integer, parameter :: verb = 2

  !!      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
  !!
  !!  kry   = krylov method def.
  !!  mass  = mass matrix
  !!  stiff = stiffness matrix
  !!  K     = linear system matrix ( = mass + stiff )
  !!  rhs   = right hand side
  !!  u_h   = numerical solution
  !!
  type(krylov) :: kry
  type(csr)    :: mass, stiff, K
  real(RP), dimension(:), allocatable :: rhs, u_h


  !!        NUMERICAL ERRORS
  !!
  real(RP) :: err_L2, err_H1

  character(LEN=300) :: shell_com

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


  !! !!!!!!!!!!!!!!!!!!!!!  INITIALISATION
  !!
  call choral_init(verb=verb)
  write(*,*)
  write(*,*)'poisson_2d: start'
  write(*,*)""
  write(*,*)"    RESOLUTION OF THE POISSON PROBLEM"
  write(*,*)""
  write(*,*)"           -Delta u + u = f   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(mass  , one_r3, x_h, qdm)
  !!
  !! 'EMetric' is the Euclidian metric on \Omega
  call diffusion_stiffmat(stiff, emetric, x_h, qdm)
  !! 
  !! K = mass + stiff
  call add(k, stiff, 1._rp, mass, 1._rp)
 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE 'F'"
  !!
  !! \int_\Omega f u_i dx
  !! The function 'f' here is defined below after the 'contains' keyword
  !!
  call l2_product(rhs, f, x_h, qdm)

  write(*,*) ""
  write(*,*)"================================== SOLVE THE LINEAR SYSTEM&
       & K U = F "
  !! 
  !! Solver setting
  kry = krylov(kry_cg, tol=1e-6_rp, itmax=1000, verb=verb)
  if (verb>1) call print(kry)
  !!
  !! allocate memory for the solution
  call allocmem(u_h, x_h%nbDof)
  u_h = 0.0_rp
  !!
  !! linear system inversion
  call solve(u_h, kry, rhs, k)

  write(*,*) ""
  write(*,*)"================================== POST-TREATMENT"
  !!
  !! L2 and H1_0 numerical errors
  !! The exact solution 'u' and its gradient 'grad_y'
  !! are  defined below after the 'contains' keyword
  !!
  write(*,*) ""
  write(*,*) "Numerical errors"
  err_l2 = l2_dist(u, u_h, x_h, qdm)
  err_h1 = l2_dist_grad(grad_u, u_h, x_h, qdm)
  write(*,*)"  Error | u - u_h|_L2            =", err_l2
  write(*,*)"  Error |\nabla (u - u_h)|_L2    =", err_h1

  if (verb>=1) then
     !!
     !! graphical display
     !!
     write(*,*) ""
     write(*,*) "Graphical display"
     !! output file def.
     shell_com = 'poisson_2d_u.msh'
     !!
     !! first write  the 'finite element mesh'
     call write(x_h, trim(shell_com), 'gmsh')
     !!
     !! second save the finite element solution u_h
     call gmsh_addview(x_h, u_h, &
          & trim(shell_com), &
          & 'Poisson problem: numerical solution', 0.0_rp, 1)
     !!
     !! 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

  call freemem(u_h)
  call freemem(rhs)

contains 

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  function f(x) 
    real(RP)                           :: f
    real(RP), dimension(3), intent(in) :: x

    f = (2._rp * pi**2 + 1._rp) * cos(pi*x(1)) * cos(pi*x(2))

  end function f


  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  function u(x) 
    real(RP)                           :: u
    real(RP), dimension(3), intent(in) :: x

    u = cos(pi*x(1)) * cos(pi*x(2))

  end function u

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  function grad_u(x) 
    real(RP), dimension(3)             :: grad_u
    real(RP), dimension(3), intent(in) :: x
    
    grad_u(1) =-pi * sin(pi*x(1)) * cos(pi*x(2))
    grad_u(2) =-pi * cos(pi*x(1)) * sin(pi*x(2))
    grad_u(3) = 0._rp

  end function grad_u


end program poisson_2d