Demo entry 6655993

zv sd

   

Submitted by anonymous on Oct 28, 2017 at 16:48
Language: Fortran. Code size: 7.6 kB.

!*************************************************************************
!******   CODAGE DE L'EQUATION D'ADVECTION LINEAIRE SCALAIRE   ***********
!*************************************************************************
!*************************************************************************
!***** Comment editer votre fichier .f90 (Quels editeurs de texte ?) *****
!*****     - sous Windows : Crimson Editor...                        *****
!*****     - sous UNIX/Linux : nedit, emacs...                       *****
!*************************************************************************
program advection
  implicit none!ça annulle l'attribution standard des types aux variables

  !##### DECLARATION DES VARIABLES #####
  !***** Parametres *****       ! parameter : variables qui resteront constantes durant l'evolution du code
  integer, parameter :: nx1=101 ! integer : entier relatif ; nx1 : nombre de points du maillage en x

  !***** Generale *****
  integer :: i,j,pos,p

  !***** MESH (griglia di calcolo)  *****
  double precision :: delta_x1             			! double precision : reel / delta_x1 : pas d'espace
  double precision :: xg                   			! limite a gauche du domaine
  double precision :: xd                   			! limite a droite du domaine
  double precision :: long_x1              			! longueur du domaine
  double precision, dimension(1:nx1) :: x1 			! tableau de reels / discretisation du domaine
  double precision :: mu_n,integral,integral_34_1,sigma_n
  double precision :: integral_partial	
  integer :: k,l,m,n
  
  

  !***** DISCRETIZZAZIONE nel TEMPO -Discretisation en temps- *****
  integer :: nt                             ! nombre de pas de temps
  double precision :: delta_t               ! pas de temps

  !***** Champ de quantite advectee *****
  double precision, dimension(1:nx1) :: w1_nplus1,w1_n,w2_n  ! quantite advectee

  !***** Vitesse d'advection *****
  double precision :: adv,swip                               ! vitesse d'advection

  !***** Autres parametres *****
  double precision :: theta                             ! parametre theta
  double precision,parameter :: mu=0			! MEDIA
  double precision,parameter :: sigma=0.4		! SCARTO 
  
  !##### DEBUT DU CODE #####
  !***** Initialisation *****
  !----- DOMINIO DI CALCOLO  -----
  xg=-2.0D+00
  xd=2.0D+00

  long_x1=xd-xg
  delta_x1=long_x1/(nx1-1)

    do i=1,nx1
     x1(i)=(i-1)*delta_x1+xg
    enddo

  !----- Champ de quantite advectee - Etat initial du champ de quantite advectee
  !----- PARAMETRI DI CALCOLO -----
  adv=0.5 
  delta_t=0.05
  nt=10000
  theta=0.75   
  
  ! Calcul de la solution à l'état initial
  
  do i=1,nx1
     w1_n(i) = exp(-x1(i)*x1(i)/(2*(sigma**2)))
  enddo

  !Calcul de l'integral sous la gaussienne
  integral=0;
  do n=1,nx1-1
    integral=integral+(w1_n(n)+w1_n(n+1))*delta_x1/2
  enddo


  !Calcul de la moyenne partie de l'integral qui definie sigma
  integral_34_1=integral*0.341

  !C'est suffis calculer lintegral à l'etat initial parce que il reste constant
   
  

  !***** Calcul de la solution *****
  !----- Ouverture du fichier de sortie -----
  !+++++ Attention : - en FORTRAN les choix pour unit sont implicites jusque '09' +++++
  !+++++               par exemple, le '06' correspond a la sortie sur ecran      +++++
  !+++++             - a une unite correspond un et un seul fichier de sortie     +++++
  !+++++               jusqu'a ce que celui-ci ait ete ferme par l'instruction    +++++
  !+++++               'close'                                                    +++++
  
  open(unit=10,file='resultats_Nx101_theta_0_dt_005_nt=200.dat')
  
  !C'est un fichier dans lequel on sauvegarde les valeur finale de sigma et la valeur maximale de w
  open(unit=15,file='w1_max_sigma_Nx101_theta_1_dt_005_nt=200.txt')

  !----- En-tete de fichier pour Tecplot (a commenter pour Gnuplot) -----
  !write(10,*) 'TITLE="RESULTATS"'
  !write(10,*) 'VARIABLES="X","U"'
  !write(10,*) 'ZONE I=',nx1

  !+++++ Ecriture de la solution initiale dans le fichier de sortie +++++
  do i=1,nx1
     write(10,'(2(f8.4,1x))') x1(i),w1_n(i)!Il ecrit 2 variables (2), separés par un espace (1x) formatés avec 8 chiffres avant de la virgule et 4 après 
  enddo

  
  !+++++ Pour une sortie sur Gnuplot (a commenter pour Tecplot) +++++
  write(10,*)
  write(10,*)

 !----- Schema numerique -----
  do j=1,nt ! Avancement temporel (Boucle en temps)

     !+++++ Mise en place du schema de dicretisation en espace +++++
     do i=2,nx1-1 ! Attention au choix des bornes de la boucles : il ne faut pas deborder des dimensions du tableau
            w1_nplus1(i) = w1_n(i) - adv*(delta_t/(2*delta_x1))*  &
             (w1_n(i+1)-w1_n(i-1))+(theta/2)*(w1_n(i+1)-2*w1_n(i)+w1_n(i-1))
     enddo   
     
     ! Calcul de w1_nplus1(nx1) : on utilise la condition de périodicité
     
            w1_nplus1(nx1) = w1_n(nx1) - adv*(delta_t/(2*delta_x1))*  &
             (w1_n(2)-w1_n(nx1-1))+(theta/2)*(w1_n(2)-2*w1_n(nx1)+w1_n(nx1-1))
     
     ! Condition de périodicité
        w1_nplus1(1)=w1_nplus1(nx1)       
    
     ! Avancement de la solution 
	w1_n=w1_nplus1
   
  
  
     !+++++ Que faire des points que l'on n'a pas pu traiter = Conditions aux limites                 +++++
     !+++++ Dans la majorite des cas que vous aurez a etudier, on utilisera principalement            +++++
     !+++++ deux types de conditions limites :                                                        +++++
     !+++++ - les conditions de Dirichlet : w=0                                                       +++++
     !+++++ - les conditions de Von Neumann : dw/dn=0 , n etant la normale a la paroi                 +++++
     !+++++ Dans le projet qui vous est propose, on vous demande en plus des conditions periodiques : +++++
     !+++++ Astuce : recopier la sortie sur l'entree !                                                +++++

     !+++++ Ecriture des resultats dans fichier de sortie a chaque pas de temps +++++
     !+++++ Pour une sortie sur Tecplot (a commenter pour Gnuplot) +++++
     !write(10,*) 'ZONE I=',nx1
     !!!Pour une sortie sur Gnuplot (a commenter pour Tecplot) +++++
     write(10,*)
     write(10,*)

     !+++++ Ecriture +++++
     do i=1,nx1
        write(10,'(2(f8.4,1x))') x1(i),w1_n(i)
     enddo

                
  enddo

!Je copie le vecteur w1_n dans le vecteur w2_n
  w2_n = w1_n
  
  !Bubble sort sur w2_n
  do k=1,nx1 
  	do l=k+1,nx1 
    	if (w2_n(k).gt.w2_n(l))then
        	swip=w2_n(k)
        	w2_n(k)=w2_n(l)
        	w2_n(l)=swip
        endif
    enddo
  enddo
  
  !Recherche de mu à l'etat final
  m=1;
  do while (w1_n(m).ne.w2_n(nx1))
    m=m+1
  end do
  pos=m;
  mu_n=x1(m)
  
  !Calcul de l'integral partial
  integral_partial=0
  p=pos
  do while (integral_partial.lt.integral_34_1)
    integral_partial=integral_partial+(w1_n(p)+w1_n(p+1))*delta_x1/2
    p=p+1
  end do
  sigma_n=x1(p)-x1(pos)

  !Ecriture de la valeur maximale de w et sigma à l'état final dans le fichier mu_sigma_finale.dat
  write(15,'(2(f8.4,1x))') w1_n(pos),sigma_n
  

  

  

  !----- Fermeture du fichier de sortie -----
  close(10)
  close(15)

end program advection
!*****************************************************************************
!***** Pour plus d'information sur le langage FORTRAN 90/95 vous pouvez  *****
!***** telecharger la documentation sur le site de l'IDRIS a l'adresse : *****
!*****   http://www.idris.fr/data/cours/lang/fortran/choix_doc.html      *****
!*****   Choisir l'edition "Fortran : notions de base" (1er niveau)      *****
!*****************************************************************************

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).