diffusion_Dirichlet_Neumann.f90

Go to the documentation of this file.

program diffusion_dirichlet_neumann

  use choral_constants
  use choral

  implicit none

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

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

  !! penal = penalty coeff for the Dirichlet condition
  !! tol   = linear system inversion tolerance
  !!
  real(RP), parameter :: penal = 1e10_rp
  real(RP), parameter :: tol   = 1e-15_rp


  !!      SPACE DISCRETISATION
  !!
  !!   fe_v     = finite element method (volumic)
  !!   fe_s     = finite element method (surfacic)
  !!   qd_v     = quadrature method (volumic)  
  !!   qd_s     = quadrature method (surfacic)      
  !!   mesh_idx = index for the mesh
  !!
  integer         , parameter :: fe_v     = fe_p2_2d
  integer         , parameter :: fe_s     = fe_p2_1d
  integer         , parameter :: qd_v     = quad_gauss_trg_12
  integer         , parameter :: qd_s     = quad_gauss_edg_4
  character(LEN=1), parameter :: mesh_idx = "2"
  !!
  !! m        = mesh
  !! X_h      = finite element space
  !! qdm      = integration method
  !!
  type(mesh)     :: msh
  type(fespace)  :: X_h
  type(quadmesh) :: qdm
  !!
  !! msh_file 
  character(len=100), parameter :: mesh_file= &
       & trim(GMSH_DIR)//"square/square_"//mesh_idx//".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, aux


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

  character(LEN=300) :: shell_com

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


  !! !!!!!!!!!!!!!!!!!!!!!  INITIALISATION
  !!
  call choral_init(verb=verb)
  write(*,*)
  write(*,*)'diffusion_Dirichlet_Neumann: start'
  write(*,*)""
  write(*,*)"    EXAMPLE OF AN ELLIPTIC PROBLEM"
  write(*,*)""
  write(*,*)"           -\Delta(u) + u = f   on  \Omega"
  write(*,*)""
  write(*,*)"    WITH A MIXED DIRICHLET / NEUMANN BOUNDARY CONDITION "  
  write(*,*)""
  write(*,*)""

  write(*,*) ""
  write(*,*)"================================== MESH ASSEMBLING"
  msh = mesh(mesh_file, 'gmsh')
  
  call print(msh)

  write(*,*) ""
  write(*,*)"================================== FINITE ELEMENT&
       & SPACE ASSEMBLING"
  x_h =fespace(msh)
  call set(x_h, fe_v)
  call set(x_h, fe_s)
  call assemble(x_h)
  call print(x_h)

  write(*,*) ""
  write(*,*)"================================== INTEGRATION METHOD&
       & ON \Omega ASSEMBLING"
  qdm = quadmesh(msh)
  call set(qdm, qd_v)
  call assemble(qdm)
  call print(qdm)

  write(*,*) ""
  write(*,*)"================================== MASS AND STIFFNESS&
       & MATRIX ASSEMBLING"
  call diffusion_massmat(mass  , one_r3, x_h, qdm)
  call diffusion_stiffmat(stiff, emetric, x_h, qdm)


  write(*,*) ""
  write(*,*)"================================== LINEAR SYSTEM MATRIX"
  !!
  !! K = mass + stiff
  call add(k, stiff, 1._rp, mass, 1._rp)
  !!
  !! cleaning
  call clear(mass)
  call clear(stiff)
 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   1- COMPUTATION OF &
       & THE VOLUMIC SOURCE TERM F"
  !!
  !!  \int_\Omega f v_i dx
  call l2_product(rhs, f, x_h, qdm)

 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   2- COMPUTATION OF  &
       & THE SURFACIC SOURCE TERM G (NEUMANN CONDITION)"
  !!
  !!  \int_{Gamma, x=1} dx_u_1  v_i dx
  call diffusion_neumann_rhs(aux, dx_u_1, x_h, qd_s, x_ge_1) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, y=1} dy_u_1  v_i dx
  call diffusion_neumann_rhs(aux, dy_u_1, x_h, qd_s, y_ge_1) 
  rhs = rhs + aux

  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   3- PENALISATION OF  &
       & THE SYSTEM FOR THE DIRICHLET CONDITION"
  !!
  !!  Dirichlet condition on {Gamma, x=0} 
  call diffusion_dirichlet(k, rhs, u, x_h, penal, x_le_0) 
  !!
  !!  Dirichlet condition on {Gamma, y=0} 
  call diffusion_dirichlet(k, rhs, u, x_h, penal, y_le_0) 

  write(*,*) ""
  write(*,*)"================================== SOLVE LINEAR SYSTEM&
       & K U = RHS "
  !! 
  !! Solver setting
  kry = krylov(kry_cg, tol=tol, itmax=1000, verb=2)
  call print(kry)
  !!
  !! initial guess 
  call interp_scal_func(u_h, u, x_h)
  u_h = 0.0_rp
  !!
  !! linear system inversion
  call solve(u_h, kry, rhs, k)
  !!
  !! cleaning
  call clear(k)
  call freemem(rhs)

  write(*,*) ""
  write(*,*)"================================== POST-TREATMENT"
  !!
  !! graphical display
  !!
  !! output file def.
  shell_com = 'diff_Neumann_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), &
       & 'Diffusion 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)
  !!
  !! L2 and H1_0 numerical errors
  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

  call freemem(u_h)
  call freemem(aux)

contains 

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

    u = cos( x(1) ) * (1._rp + x(2)**2)

  end function u


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

  end function grad_u


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

    f = 2._rp * u(x) - 2._rp * cos( x(1) )


  end function f


  function x_le_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = -x(1)

  end function x_le_0


  function x_ge_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = x(1) - 1.0_rp

  end function x_ge_1

  function y_le_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = -x(2)

  end function y_le_0


  function y_ge_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = x(2) - 1.0_rp

  end function y_ge_1


  function dx_u_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = -sin( 1.0_rp ) * ( 1.0_rp + x(2)**2 )

  end function dx_u_1

  function dy_u_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = cos( x(1) ) * 2.0_rp

  end function dy_u_1

end program diffusion_dirichlet_neumann