! ! Copyright Gloria Faccanoni, janvier 2016 ! ! ! Resolution du ! probleme de Cauchy ! d_t u(x,t) + c d_x u(x,t) = 0 ! u(x,t=0)=g(x) ! pour t>0 et x in [0;L] avec conditions au bord periodiques ! ! Creer une directory data dans laquelle seront sauvegarde les sorties ! ! --- Pour compiler && executer ! gfortran transport.f90 -o transport.o && ./transport.o ! ! --- Pour voir les ``films'' avec gnuplot ! cd data ! gnuplot plot_centree.gnu ! gnuplot plot_decentreegauche.gnu ! gnuplot plot_decentreedroite.gnu ! gnuplot plot_upwind.gnu ! gnuplot plot_laxfriedrichs.gnu ! gnuplot plot_laxwendroff.gnu ! gnuplot plot_beamwarming.gnu ! gnuplot plot_fromm.gnu ! gnuplot plot_antidiffusif.gnu ! program transport implicit none double precision, parameter :: pi=2.0*acos(0.0) integer, parameter :: Nx = 500 ! nombre de points double precision, parameter :: L = 8.0 ! domaine double precision, parameter :: c = 1. ! vitesse du transport double precision, parameter :: Tmax = 24.0 ! borne en temps integer, parameter :: Smax = 40 ! nombre de sauvegardes double precision, parameter :: cfl = 0.99 ! coefficient CFL integer :: Nt, j, sauv double precision :: dx, dt, t double precision, dimension(Nx) :: x double precision, dimension(Nx) :: exacte double precision, dimension(Nx) :: centree double precision, dimension(Nx) :: decentreegauche double precision, dimension(Nx) :: decentreedroite double precision, dimension(Nx) :: upwind double precision, dimension(Nx) :: laxfriedrichs double precision, dimension(Nx) :: laxwendroff double precision, dimension(Nx) :: beamwarming double precision, dimension(Nx) :: fromm double precision, dimension(Nx) :: antidiffusif Nt = 0 ! numero de la sauvegarde dx = L/Nx ! pas d'espace dt = cfl*abs(c)*dx ! pas de temps t = 0.0 ! instant x = (/(j*dx, j=1,Nx)/) ! position exacte = (/(0., j=1,Nx)/) ! solution exacte centree = (/(0., j=1,Nx)/) ! solution avec schema centre decentreegauche = (/(0., j=1,Nx)/) ! solution avec schema decentre a gauche decentreedroite = (/(0., j=1,Nx)/) ! solution avec schema decentre a droite upwind = (/(0., j=1,Nx)/) ! solution avec schema upwind laxfriedrichs = (/(0., j=1,Nx)/) ! solution avec schema de Lax-Friedrichs laxwendroff = (/(0., j=1,Nx)/) ! solution avec schema de Lax-Wendroff beamwarming = (/(0., j=1,Nx)/) ! solution avec schema de Beam-Warming fromm = (/(0., j=1,Nx)/) ! solution avec schema de Fromm antidiffusif = (/(0., j=1,Nx)/) ! solution avec schema AntiDiffusif ! Initialisations call Calcul_Exacte(x,t,exacte) centree = exacte decentreegauche = exacte decentreedroite = exacte upwind = exacte laxfriedrichs = exacte laxwendroff = exacte beamwarming = exacte fromm = exacte antidiffusif = exacte ! Sauvegarde de la CI sauv = 0 call sauvegarde(sauv,x,exacte,centree, decentreegauche,decentreedroite,upwind,laxfriedrichs, laxwendroff, beamwarming,& fromm, antidiffusif) ! ********************************************************************** ! MARCHE EN TEMPS ! ********************************************************************** do while (tfloor((t-dt)*Smax/Tmax) ) then sauv=sauv+1 call sauvegarde(sauv,x,exacte,centree, decentreegauche,decentreedroite,upwind,laxfriedrichs, laxwendroff, beamwarming,& fromm, antidiffusif) end if end do call scriptgnuplot(sauv) ! ********************************************************************** ! SOUS-PROGRAMMES ! ********************************************************************** contains subroutine sauvegarde(sauv,x,exacte,centree, decentreegauche,decentreedroite,& upwind,laxfriedrichs, laxwendroff, beamwarming, fromm, antidiffusif) integer :: j integer, intent(in) :: sauv double precision, dimension(Nx), intent(in) :: x, exacte,centree, decentreegauche,decentreedroite,upwind,& laxfriedrichs, laxwendroff, beamwarming, fromm, antidiffusif character(len=30) :: file_name ! nome du fichier a sauvegarder write(file_name,'("./data/transport",i3.3,".dat")') sauv open(unit=11,file=file_name) do j=1,Nx write(11,*) x(j), & exacte(j), & centree(j), & decentreegauche(j), & decentreedroite(j), & upwind(j), & laxfriedrichs(j), & laxwendroff(j), & beamwarming(j), & fromm(j), & antidiffusif(j) end do close(11) print "(A,I3,A)", " ===== Sauvegarde num.", sauv, " =====" end subroutine subroutine Calcul_Exacte(x,t,exacte) double precision, dimension(Nx), intent(in) :: x double precision, intent(in) :: t double precision, dimension(Nx), intent(out) :: exacte double precision, dimension(Nx) :: xi xi = x-c*t - floor((x-c*t)/L)*L where (xiL*0.833) exacte = 0.0 elsewhere exacte = 1.0 end where end subroutine Calcul_Exacte subroutine Calcul_Centree(centree) double precision, dimension(Nx), intent(inout) :: centree integer :: j double precision, dimension(Nx) :: centree_P, centree_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques centree_P(1:Nx-1) = centree(2:Nx) centree_P(Nx) = centree(1) centree_M(2:Nx) = centree(1:Nx-1) centree_M(1) = centree(Nx) ! schema centree = centree-alpha*0.5*(centree_P-centree_M) end subroutine Calcul_Centree subroutine Calcul_Decentree_Gauche(decentreegauche) double precision, dimension(Nx), intent(inout) :: decentreegauche integer :: j double precision, dimension(Nx) :: decentreegauche_P, decentreegauche_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques decentreegauche_M(2:Nx) = decentreegauche(1:Nx-1) decentreegauche_M(1) = decentreegauche(Nx) ! schema decentreegauche = decentreegauche-alpha*(decentreegauche-decentreegauche_M) end subroutine Calcul_Decentree_Gauche subroutine Calcul_Decentree_Droite(decentreedroite) double precision, dimension(Nx), intent(inout) :: decentreedroite integer :: j double precision, dimension(Nx) :: decentreedroite_P, decentreedroite_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques decentreedroite_P(1:Nx-1) = decentreedroite(2:Nx) decentreedroite_P(Nx) = decentreedroite(1) ! schema decentreedroite = decentreedroite-alpha*(decentreedroite_P-decentreedroite) end subroutine Calcul_Decentree_Droite subroutine Calcul_Upwind(upwind) double precision, dimension(Nx), intent(inout) :: upwind double precision, dimension(Nx) :: upwind_P, upwind_M ! conditions periodiques upwind_P(1:Nx-1) = upwind(2:Nx) upwind_P(Nx) = upwind(1) upwind_M(2:Nx) = upwind(1:Nx-1) upwind_M(1) = upwind(Nx) ! schema upwind = upwind-(dt/dx)*(c+abs(c))*0.5*(upwind-upwind_M)-(dt/dx)*(c-abs(c))*0.5*(upwind_P-upwind) end subroutine Calcul_Upwind subroutine Calcul_Lax_Friedrichs(laxfriedrichs) double precision, dimension(Nx), intent(inout) :: laxfriedrichs double precision, dimension(Nx) :: laxfriedrichs_P, laxfriedrichs_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques laxfriedrichs_P(1:Nx-1) = laxfriedrichs(2:Nx) laxfriedrichs_P(Nx) = laxfriedrichs(1) laxfriedrichs_M(2:Nx) = laxfriedrichs(1:Nx-1) laxfriedrichs_M(1) = laxfriedrichs(Nx) ! schema laxfriedrichs = laxfriedrichs_P + laxfriedrichs_M - alpha*(laxfriedrichs_P-laxfriedrichs_M) laxfriedrichs = 0.5*laxfriedrichs end subroutine Calcul_Lax_Friedrichs subroutine Calcul_Lax_Wendroff(laxwendroff) double precision, dimension(Nx), intent(inout) :: laxwendroff double precision, dimension(Nx) :: laxwendroff_P, laxwendroff_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques laxwendroff_P(1:Nx-1) = laxwendroff(2:Nx) laxwendroff_P(Nx) = laxwendroff(1) laxwendroff_M(2:Nx) = laxwendroff(1:Nx-1) laxwendroff_M(1) = laxwendroff(Nx) ! schema laxwendroff = laxwendroff & - 0.5*alpha*(laxwendroff_P - laxwendroff_M) & + 0.5*alpha**2*(laxwendroff_P-2*laxwendroff+laxwendroff_M) end subroutine Calcul_Lax_Wendroff subroutine Calcul_Beam_Warming(beamwarming) double precision, dimension(Nx), intent(inout) :: beamwarming double precision, dimension(Nx) :: beamwarming_P, beamwarming_PP, beamwarming_M, beamwarming_MM, g_P, g_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques beamwarming_P(1:Nx-1) = beamwarming(2:Nx) beamwarming_P(Nx) = beamwarming(1) beamwarming_PP(1:Nx-1) = beamwarming_P(2:Nx) beamwarming_PP(Nx) = beamwarming_P(1) beamwarming_M(2:Nx) = beamwarming(1:Nx-1) beamwarming_M(1) = beamwarming(Nx) beamwarming_MM(2:Nx) = beamwarming_M(1:Nx-1) beamwarming_MM(1) = beamwarming_M(Nx) ! schema if (c>0) then beamwarming=beamwarming-dt/dx*(beamwarming-beamwarming_M+0.5*(1.-alpha)*(beamwarming-2.*beamwarming_M+beamwarming_MM)) else ! a implementer beamwarming = 0.0*beamwarming end if end subroutine Calcul_Beam_Warming subroutine Calcul_Fromm(fromm) double precision, dimension(Nx), intent(inout) :: fromm double precision, dimension(Nx) :: fromm_P, fromm_PP, fromm_M, fromm_MM double precision :: alpha alpha = c*dt/dx ! conditions periodiques fromm_P(1:Nx-1) = fromm(2:Nx) fromm_P(Nx) = fromm(1) fromm_PP(1:Nx-1) = fromm_P(2:Nx) fromm_PP(Nx) = fromm_P(1) fromm_M(2:Nx) = fromm(1:Nx-1) fromm_M(1) = fromm(Nx) fromm_MM(2:Nx) = fromm_M(1:Nx-1) fromm_MM(1) = fromm_M(Nx) ! schema if (c>0) then fromm = alpha*(alpha-1.)*0.25*fromm_MM & + alpha*(5.-alpha)*0.25*fromm_M & + (1.-alpha)*(alpha+4.)*0.25*fromm & + alpha*(alpha-1.)*0.25*fromm_P else ! a implementer fromm = 0.0*fromm end if end subroutine Calcul_Fromm subroutine Calcul_Anti_Diffusif(antidiffusif) double precision, dimension(Nx), intent(inout) :: antidiffusif double precision, dimension(Nx) :: antidiffusif_P, antidiffusif_PP, antidiffusif_M, antidiffusif_MM double precision, dimension(Nx) :: A_P, A_M, B_P, B_M, g_P, g_M double precision :: alpha alpha = c*dt/dx ! conditions periodiques antidiffusif_P(1:Nx-1) = antidiffusif(2:Nx) antidiffusif_P(Nx) = antidiffusif(1) !antidiffusif_PP(1:Nx-1) = antidiffusif_P(2:Nx) !antidiffusif_PP(Nx) = antidiffusif_P(1) antidiffusif_M(2:Nx) = antidiffusif(1:Nx-1) antidiffusif_M(1) = antidiffusif(Nx) antidiffusif_MM(2:Nx) = antidiffusif_M(1:Nx-1) antidiffusif_MM(1) = antidiffusif_M(Nx) ! schema if (c>0) then A_P = max(antidiffusif_M ,antidiffusif )+(antidiffusif -max(antidiffusif_M ,antidiffusif ))/alpha A_M = max(antidiffusif_MM,antidiffusif_M)+(antidiffusif_M-max(antidiffusif_MM,antidiffusif_M))/alpha B_P = min(antidiffusif_M ,antidiffusif )+(antidiffusif -min(antidiffusif_M ,antidiffusif ))/alpha B_M = min(antidiffusif_MM,antidiffusif_M)+(antidiffusif_M-min(antidiffusif_MM,antidiffusif_M))/alpha where(antidiffusif_P<=A_P) g_P=A_P elsewhere(antidiffusif_P>B_P) g_P=B_P elsewhere g_P=antidiffusif_P end where where(antidiffusif<=A_M) g_M=A_M elsewhere(antidiffusif>B_M) g_M=B_M elsewhere g_M=antidiffusif end where antidiffusif = antidiffusif - alpha*(g_P-g_M) else ! a implementer antidiffusif = 0.0*antidiffusif end if end subroutine Calcul_Anti_Diffusif !******************************************************************* ! Ecriture des fichier "plot_XXX.gnu" de commandes gnuplot subroutine scriptgnuplot(sauv) integer, intent(in) :: sauv integer :: i, j, h character(len=30) :: file_name open(unit=13,file="./data/plot_centree.gnu") open(unit=14,file="./data/plot_decentreegauche.gnu") open(unit=15,file="./data/plot_decentreedroite.gnu") open(unit=16,file="./data/plot_upwind.gnu") open(unit=17,file="./data/plot_laxfriedrichs.gnu") open(unit=18,file="./data/plot_laxwendroff.gnu") open(unit=19,file="./data/plot_beamwarming.gnu") open(unit=20,file="./data/plot_fromm.gnu") open(unit=21,file="./data/plot_antidiffusif.gnu") do i=13,21 write(i,'("set yrange[-0.25:1.25];")') h=i-10 write(i,'(a58,i2,a6)') "plot 'transport000.dat' u 1:2 w l, 'transport000.dat' u 1:",h," w lp;" write(i,'("pause -1;")') do j=1,sauv h=i-10 write(i,'(a15,i3.3,a27,i3.3,a10,i2,a6)') "plot 'transport",j,".dat' u 1:2 w l, 'transport",j,".dat' u 1:",h," w lp;" write(i,'("pause 0.25;")') end do write(i,'("pause -1;")') close(i) end do open(unit=30,file="./data/plot_comparaison_finale.gnu") write(i,'("set yrange[-0.25:1.25];")') write(30,'(a18,i3.3,a6)') "fichier='transport",sauv,".dat';" write(30,*) "plot fichier u 1: 2 w l ti 'Exacte'\" !write(30,*) " , fichier u 1: 3 w lp ti 'Centre'\" !write(30,*) " , fichier u 1: 4 w lp ti 'Decentre a gauche'\" !write(30,*) " , fichier u 1: 5 w lp ti 'Decentre a droite'\" write(30,*) " , fichier u 1: 6 w lp ti 'Upwind'\" write(30,*) " , fichier u 1: 7 w lp ti 'Lax-Friedrichs'\" write(30,*) " , fichier u 1: 8 w lp ti 'Lax-Wendroff'\" write(30,*) " , fichier u 1: 9 w lp ti 'Beam-Warmimg'\" write(30,*) " , fichier u 1:10 w lp ti 'Fromm'\" write(30,*) " , fichier u 1:11 w lp ti 'Antidiffusif'\" write(30,*) " ;" write(30,'("pause -1;")') write(30,*) "set term png" write(30,*) "set output 'comparaison_finale.png'" write(30,*) "rep" close(30) end subroutine scriptgnuplot end program transport