diffusion_1D.f90

PROGRAM diffusion_1D
! *******************************************************************
!     Auteurs:
!     Date de creation:
!     Date de modification:
!
!     Ce programme permet de resoudre une equation de diffusion 1D:
!
!     dC/dt = K d2C/dx2
!
!     C(x,t): la concentration en kg/m3
!     K: la diffusivite en m2/s
!     t: le temps en s
!     x: la variable d'espace en m
!
! *******************************************************************
!
! ==========================
! Declaration des variables
! ==========================
!
  IMPLICIT NONE                      ! Toutes les variables doivent être declarees
!
  INTEGER :: i                       ! indice de la variable spatiale
  INTEGER :: n                       ! indice de la variable temporelle
  INTEGER :: NX                      ! nombre d'intervales suivant Ox
  INTEGER :: NT,NTU                  ! nombre d'iterations totales et demandées
!
  PARAMETER(NX = 100)                ! Fixe le nombre de points suivant Ox
  PARAMETER(NT = 2000)               ! Fixe le nombre d'iterations max
!
  REAl :: M0                         ! la masse initiale de polluant
  REAl :: c0                         ! la concentration initiale c(x,0)=c0
  REAL :: K                          ! la diffusivite
  REAL :: U0                          ! la speed of wind
  REAL :: tfin                       ! le temps max a atteindre
  REAL :: dt                         ! le pas de temps
  REAL :: dx                         ! le pas d'espace
  REAL :: L                          ! longueur du domaine de calcul
  REAL :: Mt                         ! Masse totale de colorant
  REAL :: per

!
  REAL, DIMENSION(0:NX, 0:NT)      :: ca   !  Solution analytique
  REAL, DIMENSION(0:NX)      :: x    ! la variable d'espace x
  REAL, DIMENSION(0:NT)      :: t    ! le temps
  REAL, DIMENSION(0:NX,0:NT) :: c    ! la concentration c(x,t)
!
  INTEGER           :: LENGTH        ! Longeur de la chaine de caracteres du nom des fichiers
  CHARACTER(LEN=4)  :: INT_TO_STRING ! Pour conversion entier -> chaine de caractères
  CHARACTER(LEN=35) :: NOMFICHANIM   ! Nom du fichier de resultats
  INTEGER :: LENGTH2
  CHARACTER(LEN=35) :: NOMFICHANIM2
!
! ===================================================
! Affectation des grandeurs physiques et géométriques
! ===================================================
!
  K = .001
  U0 = 0.5
  M0 = 1.
  L = 10.
  dx = L / NX
  c0 = M0 / dx
  per = ACOS(-1.) ! per, c'est pi!

!
! =================================
! Creation et affichage du maillage
! =================================
!
  WRITE(*,*) '--------------------------------'
  WRITE(*,*) 'Affichage des noeuds du maillage'
  WRITE(*,*) '--------------------------------'
  DO i=0,NX
     x(i) = i*dx
     WRITE(*,*) 'X(',i,')=', x(i)
  ENDDO
!
! ===================================
! Definition des paramètres temporels
! ===================================
!
!
! Definition du pas de temps
!
  WRITE(*,*) 'Donnez votre pas de temps, aussi connu sous le nom de dt :'
  READ(*,*) dt
!
! Calcul du nombre d'iterations a effectuer NTU
!
  WRITE(*,*) 'Quand est-ce qu''on s''arrete ?'
  READ(*,*) tfin
  NTU = tfin / dt
  IF (NTU.GT.NT) THEN
     WRITE(*,*) 'Désolé garcon, le temps max est: ', NT*dt
     NTU = NT
     WRITE(*,*) 'Le nombre d''iterations sera donc : ', NTU, ', désolé mec'
  ENDIF
!
! Definition du tableau temps
!
  DO n=0,NTU
     t(n) = n*dt
  ENDDO
!
! =====================================
! Mise à zéro du champ de concentration
! =====================================
!
  DO i=0,NX
     if (i<25) then
        c(i, 0) = MAX(0, i-15)
     else
        c(i, 0) = MAX(0, 35-i)
    end if
  ENDDO
!
! ==================
! Condition initiale
! ==================
!
  !c(NX/2, 0) = c0
!
! ================================
! Calcul du champ de concentration
! ================================
!
  WRITE(*,*) 'Calcul en cours, va chercher un café...'
!
! -----------------------------
! Debut des iterations en temps
! -----------------------------
  DO n=0,NTU
!
     WRITE(*,*) '-------------------------------'
     WRITE(*,*) 'Iteration n=',n,' avec t=',t(n)
     WRITE(*,*) '-------------------------------'
!
! Calcul de la solution au pas de temps n+1
!
    DO i=1,NX-1
        c(i, n+1) = U0*dt/(dx*2) * ( c(i-1, n) - c(i+1, n)) + c(i, n)
        ca(i, n+1) = c0 * (i*dx - U0*n*dt)
    ENDDO

!
! Conditions aux limites
!
    c(0, n+1) = 0

! ---------------------------
! Fin des iterations en temps
! ---------------------------
!
  ENDDO

  ! Vérification de la masse totale
    Mt = 0
    DO i=1,NX-1
        Mt = Mt + c(i, NTU)
    ENDDO
    Mt = Mt + .5 * (c(0, NTU) + c(NX, NTU))
    WRITE(*,*) 'M0 = ', M0, ', Mt = ', Mt*dx
!
  WRITE(*,*) 'Victoire ! Le calcul est termine avec succes \o/'
!
!
! ==========================================================================
! Ecriture du champ de concentration dans des fichiers pour chaque iteration
! ==========================================================================
!
  WRITE(*,*) 'Ecriture des resultats dans des fichiers ...'
!
  DO n=0,NTU
!
! Definition du nom du fichier
!
!
     IF (n.LT.10)                    WRITE (int_to_string, '(I1)') n
     IF ((n.LT.100).AND.(n.GE.10))   WRITE (int_to_string, '(I2)') n
     IF ((n.LT.1000).AND.(n.GE.100)) WRITE (int_to_string, '(I3)') n
     IF ((n.LT.10000).AND.(n.GE.1000)) WRITE (int_to_string, '(I4)') n
!
     NOMFICHANIM='diffusion_1D_'//int_to_string
     LENGTH=INDEX(NOMFICHANIM,' ')-1
     IF (LENGTH.LT.0) LENGTH=LEN(NOMFICHANIM)
     NOMFICHANIM=NOMFICHANIM(1:LENGTH)//'.vtk'


     NOMFICHANIM2='diffusion_analytique_'//int_to_string
     LENGTH2=INDEX(NOMFICHANIM2,' ')-1
     IF (LENGTH2.LT.0) LENGTH2=LEN(NOMFICHANIM2)
     NOMFICHANIM2=NOMFICHANIM2(1:LENGTH2)//'.vtk'

!
! Ouverture du fichier
!
     OPEN(                      &
          UNIT   = 9,           &
          FILE   = NOMFICHANIM, &
          STATUS = 'UNKNOWN',   &
          FORM   = 'FORMATTED'  &
          )

    OPEN(                      &
          UNIT   = 10,           &
          FILE   = NOMFICHANIM2, &
          STATUS = 'UNKNOWN',   &
          FORM   = 'FORMATTED'  &
          )
!
! Ecriture dans le fichier
!
     WRITE(9,'(''# vtk DataFile Version 2.0'')')
     WRITE(9,'(''2D scalar data'')')
     WRITE(9,'(''           '')')
!
     WRITE(9,'(''ASCII'')')
!
     WRITE(9,'(''DATASET STRUCTURED_GRID'')')
     WRITE(9,'(''DIMENSIONS'',3i7)') NX+1,1,1
     WRITE(9,'(''POINTS '',i7,'' float'')') (NX +1)*1*1
     DO i=0,NX
        WRITE(9,*) x(i),0,0
     ENDDO
!
     WRITE(9,'(''POINT_DATA '',i7)') (NX+1)*1*1
!
     WRITE(9,'(''SCALARS concentration float'')')
     WRITE(9,'(''LOOKUP_TABLE default'')')
     DO i=0,NX
        WRITE(9,*) c(i,n)
     ENDDO


     ! Solution analytique
     WRITE(10,'(''# vtk DataFile Version 2.0'')')
     WRITE(10,'(''2D scalar data'')')
     WRITE(10,'(''           '')')
!
     WRITE(10,'(''ASCII'')')
!
     WRITE(10,'(''DATASET STRUCTURED_GRID'')')
     WRITE(10,'(''DIMENSIONS'',3i7)') NX+1,1,1
     WRITE(10,'(''POINTS '',i7,'' float'')') (NX +1)*1*1
     DO i=0,NX
        WRITE(10,*) x(i),0,0
     ENDDO
!
     WRITE(10,'(''POINT_DATA '',i7)') (NX+1)*1*1
!
     WRITE(10,'(''SCALARS concentration float'')')
     WRITE(10,'(''LOOKUP_TABLE default'')')
     DO i=0,NX
        WRITE(10,*) ca(i, n)
     ENDDO
!
  ENDDO
!
!     Fermeture du fichier
!
  CLOSE(9)
  CLOSE(10)
!
!     FIN DU PROGRAMME
!----------------------------------------------------------------------
END PROGRAM diffusion_1D