diffusion_Neumann.f90

Go to the documentation of this file.

program diffusion_neumann

  use choral_constants
  use choral

  implicit none

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

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

  !! Diffusivity coefficients
  !!
  real(RP), parameter :: L_1 = 1.0_rp
  real(RP), parameter :: L_2 = 1.0_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
  !!
  !!  precc_type = preconditioning type
  !!
  integer, parameter :: prec_type = pc_0
  !!
  !!  kry   = krylov method def.
  !!  mass  = mass matrix
  !!  stiff = stiffness matrix
  !!  K     = linear system matrix ( = mass + stiff )
  !!  pc    = preconditionner 
  !!  rhs   = right hand side
  !!  u_h   = numerical solution
  !!
  type(krylov) :: kry
  type(csr)    :: mass, stiff, K
  type(precond)   :: pc
  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=1)
  write(*,*)
  write(*,*)'diffusion_Neumann: start'
  write(*,*)""
  write(*,*)"    EXAMPLE OF AN ELLIPTIC PROBLEM"
  write(*,*)""
  write(*,*)"           -div(B(x) grad u) + rho(x)u = f   on  \Omega"
  write(*,*)""
  write(*,*)"    WITH A NEUMANN BOUNDARY CONDITION "  
  write(*,*)""
  write(*,*)"            B(x)grad u . n = g   on   \partial\Omega"
  write(*,*)""
  write(*,*)"    AND WITH AN ANISOTROPIC DIFFUSIVITY TENSOR B(x)."
  write(*,*)""
  write(*,*)"    The problem data 'B', 'u', 'f' and 'g' "
  write(*,*)"    have been derived in"
  write(*,*)"    choral/maple/example_diffusion_Neumann.mw"
  write(*,*)""
  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  , rho, x_h, qdm)
  call diffusion_stiffmat(stiff, bl_form_b, x_h, qdm)


  write(*,*) ""
  write(*,*)"================================== LINEAR SYSTEM MATRIX&
       & AND PRECONDITIONING"
  !!
  !! K = mass + stiff
  call add(k, stiff, 1._rp, mass, 1._rp)
  !!
  !! cleaning
  call clear(mass)
  call clear(stiff)
  !!
  !! preconditioning
  pc = precond(k, prec_type)

 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   1- COMPUTATION OF 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(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   2- COMPUTATION OF G"
  !!
  !!  \int_{Gamma, x=0} g  u_i dx
  !! The function 'g' here is defined below after the 'contains' keyword
  !! 'g_x_0' is a function that defines the boundary subdomain
  !! {x\in \Gamma, x=0}
  !!
  call diffusion_neumann_rhs(aux, g_x_0, x_h, qd_s, x_le_0) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, x=1} g  u_i dx
  call diffusion_neumann_rhs(aux, g_x_1, x_h, qd_s, x_ge_1) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, y=0} g  u_i dx
  call diffusion_neumann_rhs(aux, g_y_0, x_h, qd_s, y_le_0) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, y=1} g  u_i dx
  call diffusion_neumann_rhs(aux, g_y_1, x_h, qd_s, y_ge_1) 
  rhs = rhs + aux

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

  write(*,*) ""
  write(*,*)"================================== POST-TREATMENT"
  if (verb>1) then
     !!
     !! 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)

  end if

  !!
  !! 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

  call freemem(u_h)
  call freemem(aux)

contains 

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

    u = sin(pi*x(1) ) * sin(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 * cos(pi*x(1)) * sin(pi*x(2))
    grad_u(2) = pi * sin(pi*x(1)) * cos(pi*x(2))
    grad_u(3) = 0._rp

  end function grad_u


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

    real(RP) :: cx, sx, cy, sy, aux
    real(RP) :: t, ct, st

    t  = theta(x)
    ct = cos(t)
    st = sin(t)

    cx = cos(pi*x(1))
    sx = sin(pi*x(1))
    cy = cos(pi*x(2))
    sy = sin(pi*x(2))

    aux =-4.0_rp *ct*st
    f   = aux*cx*cy

    aux = 2.0_rp*ct*st - 2.0_rp*ct**2 + 1.0_rp
    f   = f + aux*cx*sy

    aux =-2.0_rp*ct*st - 2.0_rp*ct**2 + 1.0_rp
    f = f + aux*sx*cy

    f = f * (l_1-l_2)

    aux = 2.0_rp*(l_1+l_2)
    f = f + aux*sx*sy

    f = f * pi**2 * 0.5_rp + rho(x)*u(x)

  end function f

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

    rho = sum(x*x)
    rho = log( 1.0_rp + rho)

  end function rho

  !! Bilinear form associated with the diffuion tensor B(x)
  !! \f$ (x, \xi_1, \xi_2) \mapsto B(x)\xi_1 \cdot \xi_2 \f$
  !!   
  function bl_form_b(x, w1, w2) result(res)
    real(RP)                           :: res
    real(RP), dimension(3), intent(in) :: x, w1, w2

    real(RP), dimension(2,2) :: B_x
    real(RP), dimension(2  ) :: B_w1

    b_x = b(x)
    call matvecprod(b_w1, b_x, w1(1:2))
    
    res = b_w1(1)*w2(1) + b_w1(2)*w2(2)

  end function bl_form_b


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

    theta = 0.5_rp * pi * ( x(1) + x(2) )

  end function theta

  function b(x) 
    real(RP), dimension(2,2)           :: B
    real(RP), dimension(3), intent(in) :: x

    real(RP) :: t, c, s

    t = theta(x)
    c = cos(t)
    s = sin(t)

    t =  l_1 * c**2 + l_2 * s**2
    b(1,1) = t
    b(2,2) = l_1 + l_2  - t

    t = (l_1 - l_2) * c*s
    b(1,2) = t
    b(2,1) = t

  end function b
  


  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 g_x_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    real(RP) :: t, c2, s2

    t  = theta(x)
    c2 = cos(t)**2
    s2 = 1.0_rp - c2

    r =  - pi * sin(pi*x(2)) * (l_1*c2 + l_2*s2)

  end function g_x_0


  function g_x_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = g_x_0(x)

  end function g_x_1


  function g_y_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    real(RP) :: t, c2, s2

    t  = theta(x)
    c2 = cos(t)**2
    s2 = 1.0_rp - c2

    r = - pi * sin(pi*x(1)) * (l_1*s2 + l_2*c2)

  end function g_y_0


  function g_y_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = g_y_0(x)

  end function g_y_1


end program diffusion_neumann