PROGRAM ADVEKTION USE Dislin Implicit none !_______________________________________________________________________________ #include'fftw3.f' !Befehlzum Kompilieren: ifort Advektion.f90 -fpp -lfftw3 -WARN -ldislin !_______________________________________________________________________________ INTEGER k,i PARAMETER N=2048, n_iter=80000, screen_every=100 DOUBLE PRECISION:: pi,deltat,c DOUBLE COMPLEX:: imi INTEGER*8 plan_hin,plan_her, plan_dummy DOUBLE PRECISION, DIMENSION(N):: in, dummy_feld DOUBLE PRECISION, DIMENSION(N/2+1):: kvec DOUBLE COMPLEX, Dimension(N/2+1):: out,K1,K2,K3,K4,us, dummy_c_feld !_______________________________________________________________________________ !Programm zur Berechnung der Advektionsgleichung c=0.01 imi=(0.0d0,1.0d0) pi=DACOS(-1.0d0) Deltat=0.01 !_______________________________________________________________________________ !Initialisierung des Feldes open(10,File='Exponentialfunktion.dat') DO k=1,N in(k)=exp(-2.0*pi*(2.0*pi*(k-1)/real(N)-pi/2.0)**2) write(10,*) 2.0*pi/REAL(n)*(k-1),in(k) END DO close(10) ! Wellenzahlen initialisieren DO k=1,N/2+1 kvec(k)=k-1 END DO !_______________________________________________________________________________ !_______________________________________________________________________________ ! Dislin initialisieren CALL METAFL ('XWIN') ! Dislin Initialisierung CALL PAGE(2500,2500) CALL SCRMOD('REVERS') CALL SCLMOD('FULL') CALL DISINI() !_______________________________________________________________________________ !_______________________________________________________________________________ ! Plaene initialisieren CALL dfftw_plan_dft_r2c_1d(plan_hin,N,in,out,FFTW_ESTIMATE) CALL dfftw_plan_dft_c2r_1d(plan_her,N,out,in,FFTW_ESTIMATE) CALL dfftw_plan_dft_c2r_1d(plan_dummy,N,dummy_c_feld,dummy_feld,FFTW_ESTIMATE) !_______________________________________________________________________________ ! in den Fourierraum CALL dfftw_execute_dft_r2c(plan_hin, in, out) out=out/DBLE(N)**0.5d0 !_______________________________________________________________________________ ! Runge-Kutta-Verfahren für Zeitableitung do k=1,N/2+1 us(k)=out(k) end do DO i=1,n_iter ! dislin, eine dummy Kopie vom Feld IF(MOD(i,screen_every).EQ.0) THEN dummy_c_feld=us CALL dfftw_execute_dft_c2r(plan_dummy,dummy_c_feld,dummy_feld) dummy_feld=dummy_feld/DBLE(N)**0.5d0 CALL display(dummy_feld,n,i,c,deltat) !Pause END IF CALL RechteHand(us,k1,N,c,kvec) CALL RechteHand(us+deltat/2.0*k1,k2,N,c,kvec) CALL RechteHand(us+deltat/2.0*k2,k3,N,c,kvec) CALL RechteHand(us+deltat*k3,k4,N,c,kvec) us=us+(k1/6.0+k2/3.0+k3/3.0+k4/6.0)*Deltat END DO !_______________________________________________________________________________ !_______________________________________________________________________________ ! in den Ortsraum CALL dfftw_execute_dft_c2r(plan_her,us,in) in=in/DBLE(N)**0.5d0 open(14,File='Ergebnis.dat') DO k=1,N write(14,*) 2.0*pi/REAL(n)*(k-1),in(k) END DO close(14) !_______________________________________________________________________________ !Zerstörung der Plaene CALL dfftw_destroy_plan(plan_hin) CALL dfftw_destroy_plan(plan_her) !_______________________________________________________________________________ CALL DISFIN() END PROGRAM ADVEKTION !******************************************************************************* !******************************************************************************* SUBROUTINE RechteHand(in,out,n,c,kvec) !_______________________________________________________________________________ !Rechte-Hand-Seite INTEGER :: n,k DOUBLE COMPLEX,DIMENSION(N/2+1):: in,out DOUBLE PRECISION,DIMENSION(N/2+1):: kvec DOUBLE PRECISION:: c DOUBLE COMPLEX:: imi imi=DCMPLX(0.0d0,1.0d0) ! Ableitung berechnen DO k=1,N/2+1 out(k)=-c*imi*kvec(k)*in(k) END DO !_______________________________________________________________________________ END SUBROUTINE RechteHand !*********************************************************************************** !*********************************************************************************** SUBROUTINE display(feld_in,n,zeit,c,deltat) USE dislin IMPLICIT NONE INTEGER :: i,n,zeit,IC DOUBLE PRECISION, DIMENSION(n) :: feld_in,xwert REAL, DIMENSION(n) :: x_achse, y_achse,y2_achse,y3_achse REAL :: pi, mini, maxi DOUBLE PRECISION:: c,deltat pi=ACOS(-1.0) DO i=1,n x_achse(i)=2.0*pi/REAL(n)*(i-1) xwert(i)=2.0*pi*(i-1)/REAL(N)-(c*zeit*deltat) DO WHILE (xwert(i)<0) xwert(i)=xwert(i)+2*pi END DO y_achse(i)=REAL(feld_in(i)) y2_achse(i)=REAL(EXP(-2.0*pi*(xwert(i)-pi/2.0)**2)) y3_achse(i)=REAL(EXP(-2.0*pi*(2.0*pi*(i-1)/REAL(N)-pi/2.0)**2)) END DO CALL AXSPOS(450,1800) CALL AXSLEN(1200,1200) CALL AXSSCL('LIN','XY') CALL NAME('x-axis','X') CALL NAME('y-axis','Y') CALL LABDIG(-1,'X') CALL TICKS(10,'XY') CALL TITLIN('Demonstration',1) CALL TITLIN('Advektionsgleichung',2) mini=MINVAL(y_achse) maxi=MAXVAL(y_achse+0.1) CALL GRAF(-0.1, 2.0*pi, 0.0, 1.0 , mini, maxi ,mini, (maxi-mini)/10.0) CALL SETRGB(0.7,0.7,0.7) !CALL GRID(1,1) CALL COLOR('FORE') CALL TITLE() IC=INTRGB(0.95,0.95,0.95) CALL AXSBGD(IC) !CALL SETRGB(0.7,0.7,0.7) CALL GRID(1,1) CALL COLOR('Green') CALL CURVE(x_achse,y_achse,n) CALL COLOR('Blue') CALL CURVE(x_achse,y2_achse,n) CALL COLOR('Red') CALL CURVE(x_achse,y3_achse,n) CALL ENDGRF() END SUBROUTINE !*********************************************************************************** !***********************************************************************************