Résolution numérique de l'équation de convection bidimentionelle par différences finis

Ce code permet de résoudre numériquement par différences finis l'équation de convection bidimensionnelle. Le code implémente un schéma décentré et facile à le convertir pour qu'il soit volume finis.

!-----------------------------------------------------------------------------!
!  FD Method for solving convection equation in 2D
!     du/dt = vx(du/dx )+ vy(du/dy) + q(x,y)
!     Upwind scheme + 1D Euler
!     Drichlet b.c.
!-----------------------------------------------------------------------------!
! Written by Elmiloud Chaabelasri
!            University Mohammed First
!            Faculty of sciences, Oujda, Morocco
!            chaabelasri@gmail.com
! 
! Last updated: Jul. 06, 2019
!-----------------------------------------------------------------------------!

      program Conv2d
      implicit none
      integer::i,j,k,nx,ny,nt,ns,nf
      real*8 ::dx,dy,dt,x0,xL,y0,yL,t,Tmax,alpha
      real*8 ::vx,vy
      real*8,allocatable ::u(:,:),x(:),y(:),fluxx(:,:),fluxy(:,:)
      
      !Domain
      x0 =-1.0d0 !left
      xL = 1.0d0 !right
      
      y0 =-1.0d0 !bottom
      yL = 1.0d0 !up
      
      !number of points
      nx = 50
      ny = 50
      
      !Velocity fields
      vx=1.d0
      vy=1.d0
      
      !grid spacing (spatial)
      dx = (xL-x0)/dfloat(nx)
      dy = (yL-y0)/dfloat(ny)
      
      !spatial coordinates 
      allocate(x(0:nx))
      do i=0,nx
      x(i) = x0 + dfloat(i)*dx
      end do
      
      allocate(y(0:ny))
      do j=0,ny
      y(j) = y0 + dfloat(j)*dy
      end do
      
      
      !maximum time desired
      Tmax = 1.d0
      
      !diffusion constant:
      alpha = 1.0d0
      
      !time step
      dt = 0.005d0
      
      !number of points in time
      nt = nint(Tmax/dt)
      
      !number of snapshots to plot
      ns = nint(dfloat(nt)/10)
      
      !frequency for plotting
      nf = nint(dfloat(nt)/dfloat(ns))
      
      !u: convected variable 
      allocate(u(0:nx,0:ny))
      allocate(fluxx(0:nx,0:ny))
      allocate(fluxy(0:nx,0:ny))
      
      !initial condition
      t = 0.0d0
      do j=0,ny
        do i=0,nx
         u(i,j) = 0.0d0
         if( x(i).gt.-0.9d0 .and. x(i).lt.-0.60 .and.
     &       y(j).gt.-0.9d0 .and. y(j).lt.-0.6d0 )then
             u(i,j) = 1.d0
          endif
        end do
      end do
      
      !Plot initial condition
      open(19,file='Initial.plt')
      write(19,*) 'variables ="x","y","u"'
      write(19,100)'zone f=point i=',nx+1,',j=',ny+1,',t="time',t,'"'
      do j=0,ny
      do i=0,nx
      write(19,*) x(i),y(j),u(i,j)
      end do
      end do
      close(19)
      
!      open(18,file='EvolutionShow.plt')
!      write(18,*) 'variables ="x","y","u"'
!      write(18,100)'zone f=point i=',nx+1,',j=',ny+1,',t="time',t,'"'
!      do j=0,ny
!      do i=0,nx
!      write(18,*) x(i),y(j),u(i,j)
!      end do
!      end do

      !time integration
      do k=1,nt

        ! Calcul flux at interfaces
        call CalculFlux(nx,ny,dx,dy,dt,t,alpha,x,y,u,fluxx,fluxy)
      
        !Time marching
       do j=1,ny-1
 do i=1,nx-1
          u(i,j) =  u(i,j) - (vx*dt/dx)*(fluxx(i,j) - fluxx(i-1,j) )
     &                       - (vy*dt/dy)*(fluxy(i,j) - fluxy(i,j-1) )
 end do
       end do
       
       
       !Boundary conditions
       !Bottom  and top
 do i=1,nx
          u(i,0) = u(i,1)
          u(i,ny) = u(i,nx-1)
 end do

       !right and left
 do j=1,ny
          u(0,j) = u(1,j)
          u(nx-1,j) = u(nx-1,j)
 end do

        !update t
        t = t+dt
      
          !plot field
!        if (mod(k,nf).eq.0) then
!        write(18,100)'zone f=point i=',nx+1,',j=',ny+1,',t="time',t,'"'
!        do j=0,ny
!   do i=0,nx
!          write(18,*) x(i),y(j),u(i,j)
!   end do
!        end do
!        end if

       print*,k,t,maxval(u)
      end do

      close(18)
      
      !Plot initial condition
      open(18,file='Finale.plt')
      write(18,*) 'variables ="x","y","u"'
      write(18,100)'zone f=point i=',nx+1,',j=',ny+1,',t="time',t,'"'
      do j=0,ny
      do i=0,nx
      write(18,*) x(i),y(j),u(i,j)
      end do
      end do


100    format(a16,i8,a4,i8,a10,f10.4,a3)
      end

!-----------------------------------------------------------------------------!
! Calcul flux with DF
!-----------------------------------------------------------------------------!
      subroutine CalculFlux(nx,ny,dx,dy,dt,t,alpha,x,y,u,fluxx,fluxy)
      implicit none
      integer::nx,ny,i,j
      real*8 ::dx,dy,dt,t,alpha
      real*8 ::u(0:nx,0:ny),x(0:nx),y(0:ny)
      real*8 ::fluxx(0:nx,0:ny),fluxy(0:nx,0:ny)

      
      do j=0,ny
        do i=0,nx
         fluxx(i,j) = u(i,j)
         fluxy(i,j) = u(i,j)
        end do
      end do
      

      
      end