PROGRAM Burgers USE Dislin Implicit none !_______________________________________________________________________________ #include'fftw3.f' !Befehlzum Kompilieren: ifort Burgers.f90 -fpp -lfftw3 -WARN -ldislin !_______________________________________________________________________________ INTEGER k,i,g PARAMETER N=2048, n_iter=80000, screen_every=200 INTEGER*8 plan_hin,plan_her, plan_dummy DOUBLE PRECISION :: pi,deltat,c,amplitude DOUBLE PRECISION, DIMENSION(N) :: in, dummy_feld,sinusse DOUBLE PRECISION, DIMENSION(N/2+1):: kvec DOUBLE COMPLEX :: imi DOUBLE COMPLEX, Dimension(N/2+1) :: out,K1,K2,K3,K4,us, dummy_c_feld !_______________________________________________________________________________ !Programm zur Berechnung der Advektionsgleichung c=0.01 !Viskosität imi=(0.0d0,1.0d0) pi=DACOS(-1.0d0) Deltat=0.0001 !_______________________________________________________________________________ !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) DO g=1,12 amplitude=(-1.0)**g!*(real(g))**(0.1) sinusse(k)=amplitude*sin(g*2.0*pi*(k-1)/real(N)+3.0*real(g)) in(k)=in(k)+sinusse(k) END DO in(k)=in(k)+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,us) us=us/DBLE(N)**0.5d0 !_______________________________________________________________________________ ! Runge-Kutta-Verfahren für Zeitableitung 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) END IF CALL RechteHand(us,k1,N,c,kvec,plan_hin,plan_her) CALL RechteHand(us+deltat/2.0*k1,k2,N,c,kvec,plan_hin,plan_her) CALL RechteHand(us+deltat/2.0*k2,k3,N,c,kvec,plan_hin,plan_her) CALL RechteHand(us+deltat*k3,k4,N,c,kvec,plan_hin,plan_her) 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 Burgers !******************************************************************************* !******************************************************************************* SUBROUTINE RechteHand(in,out,n,c,kvec,plan_hin,plan_her) #include'fftw3.f' !_______________________________________________________________________________ !Rechte-Hand-Seite INTEGER :: N,k,i INTEGER*8 plan_hin,plan_her DOUBLE COMPLEX,DIMENSION(N/2+1):: in,out,dusus,duscomp,reserve DOUBLE PRECISION,DIMENSION(N/2+1):: kvec DOUBLE PRECISION,DIMENSION(N):: dususreal,usreal,dusreal DOUBLE PRECISION:: c DOUBLE COMPLEX:: imi imi=DCMPLX(0.0d0,1.0d0) !___________________________________________________________ ! 1. Ableitung berechnen DO k=1,N/2+1 duscomp(k)=imi*kvec(k)*in(k) reserve(k)=in(k) END DO !___________________________________________________________ !Übergang in Realraum zur Berechnung der Nichtlinearität CALL dfftw_execute_dft_c2r(plan_her,reserve,usreal) usreal=usreal/DBLE(N)**0.5d0 CALL dfftw_execute_dft_c2r(plan_her,duscomp,dusreal) dusreal=dusreal/DBLE(N)**0.5d0 ! ! DO i=1,N dususreal(i)=usreal(i)*dusreal(i) END DO ! ! !____________________________________________________________ ! !Übergang zurück in den Fourier-Raum CALL dfftw_execute_dft_r2c(plan_hin,dususreal,dusus) dusus=dusus/DBLE(N)**0.5d0 DO k=1,N/2+1 out(k)=-c*(kvec(k))**2*in(k)-dusus(k) End DO END SUBROUTINE RechteHand !*********************************************************************************** !*********************************************************************************** SUBROUTINE display(feld_in,N) USE dislin IMPLICIT NONE INTEGER :: i,N DOUBLE PRECISION, DIMENSION(N) :: feld_in REAL, DIMENSION(n) :: x_achse, y_achse REAL :: pi, mini, maxi pi=ACOS(-1.0) DO i=1,N x_achse(i)=2.0*pi/REAL(n)*(i-1) y_achse(i)=REAL(feld_in(i)) END DO CALL ERASE() 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('Waermeleitungsgleichung',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() CALL GRID(1,1) CALL COLOR('Red') CALL CURVE(x_achse,y_achse,n) CALL ENDGRF() END SUBROUTINE !*********************************************************************************** !***********************************************************************************