Program “burger.for”

Program description

The program was written in Fortran 90.

The scheme used is the Crank-Nicholson generalised scheme.

The program permits to realize two types of calculation, by means of two types of initialization.

The first initialization has rigid boundary conditions: the initial value of u is a sinusoidal curve from x=0 to 1. The boundary values are set to 0.

The second initialization has free boundary conditions. The initial value of u is a straight line from 1 to -1 while x=0 to 1. The boundary condition is free, so the fluxes in 0 and 1 are set to 0.

Input-Output

The input and output phasis are made by means of three files: input.txt, output.txt and burger.ini.

input

The input file contains all the necessary values to perform the calculation:

Variables	Use
Nu	Cinematic viscosity (m²/s)
delta_t	Time step (s)
n_iter	Number of time iterations
pas_sortie	Output time step (iterations)
n_points	Number of points in space (<121)
scheme	Choice of the scheme: 1: finite differences 2: finite elements 3: mass operator 4: 4pts-upwind
init	Type of initialization: 1: sinusoidal 2: straight line

output

The output file is divided into two sections: the first describes the parameters of the simulation and the second contains the results, in a matrix format (u in fonction of x and t).

format

The program needs to read the format of writing in a configuration file, burger.ini. To avoid problems at the writing, the maximum number of points in space is 121.

Solver

The program calculates the matrix corresponding to the scheme and storages it in a vector format.

To solve the system AX=B, it uses the subroutine “GAUSS”, a subroutine designed to solve the systems with a band matrix. This subroutine was furnished by Mrs. Maubourguet from the IMFT laboratory (Toulouse) (http://www.imft.fr).

The listing and the arguments of the subroutine are detailed below:

c*****************************************************
c MODE D'UTILISATION : CALL GAUSS (A,B,LV,MX) *
c*****************************************************
C Donnees : A Coefficients de la matrice bande, *
* variable a une dimension dont les *
* valeurs numeriques doivent etre *
* fournies conformement au schema *
* ci-dessous. *
* Les points exterieurs a la bande ne *
* sont pas pris en compte pdt le calcul *
* A la sortie de GAUSS, A est detruite *
* B Termes du 2eme membre, variable a un *
* indice. A la sortie de GAUSS, la *
* solution se trouve dans B. *
* LV Ordre de la matrice A. *
* LB Largeur de la bande = 2 MX + 1; *
* MX est appelee "demi largeur de bande" *
* *
******************************************************

SUBROUTINE GAUSS(A,B,LV,MX)
DIMENSION A(1),B(1)
LMOI=LV-MX
ICODE=1
MX2=2*MX+1
MX1=MX+1
M22=MX2-1
L=MX1
DO K=1,LV
DIVB=1./A(L)
LMX=L+MX
IF(K .GT. LMOI) THEN
LMX=LMX-ICODE
ICODE=ICODE+1
ENDIF
DO J=L,LMX
A(J)=A(J)*DIVB
ENDDO
B(K)=B(K)*DIVB
IF(K.NE.LV) THEN
MI22=0
DO II=ICODE,MX
MI22=MI22+M22
LMI22=L+MI22
IF (A(LMI22).NE.0) THEN
COEFF=A(LMI22)
DO J=L,LMX
JMI22=J+MI22
A(JMI22)=A(JMI22)-A(J)*COEFF
ENDDO
III=II-ICODE+1
KIII=K+III
B(KIII)=B(KIII)-B(K)*COEFF
ENDIF
ENDDO
L=L+MX2
ENDIF
ENDDO
ICODE=0
I=LV
L=MX1+1+(LV-1)*MX2
DO WHILE (I.GT.1)
SIGMA=0.
L=L-MX2
LMX=L+MX-1
IF(I.GT.LMOI) THEN
LMX=L+ICODE
ICODE=ICODE+1
ENDIF
II=I
DO J=L,LMX
SIGMA=SIGMA+A(J)*B(II)
II=II+1
ENDDO
I=I-1
B(I)=B(I)-SIGMA
ENDDO
RETURN
END

If you want to download the source code of the subroutine, click here.

Program listing

The listing of the program is entirely commented and uses the different files mentionned above to run. It is compiled in Fortran 90 and designed to run on a Windows environment, but is also compliant with other OS, with an appropriate compiler.

It manages the different types of schemes, initializations and boundary conditions detailed in the description of the scheme.

You can download the program with the appropriate files and configurations by clicking on the button below. The download file is under the ZIP format (size 8Ko). The program is not compiled, so you will need a Fortran90 compiler to use it.

Download program

The listing of the program is detailed below:

c Program burger.for: resolution de l'equation de Burgers 1D
IMPLICIT NONE

c***********************************************************
c declaration des variables *
c***********************************************************
!variables de boucle
INTEGER i, j, t
!declaration des constantes
REAL PI

!pas de temps et d'espace, visc., param. num.
INTEGER n, k !n_iter, n_espace
REAL deltax, deltat, s ,r
REAL q, delta !param num schema
REAL nu !visc. cin.

!choix du schema et de l'initialisation
INTEGER klim !indce de changement de sgn
INTEGER scheme !var de schema
INTEGER init !var d'initialisation
CHARACTER*10 name_scheme
CHARACTER*10 name_init

!variables de sortie des resultats
INTEGER pas_sortie, compteur_time, cur_time
CHARACTER*2000 frmt

!definition des matrices et vecteurs
REAL mat,d,u0
ALLOCATABLE mat(:),d(:),u0(:)

c definition des constantes ********************************
PI = 3.14159265359

c***********************************************************
c initialisation *
c***********************************************************
c**** lecture des donnees **********************************
!donnees
OPEN (1,FILE='input.txt')
READ (1,*) nu
READ (1,*) deltat
READ (1,*) n
READ (1,*) pas_sortie
READ (1,*) k
READ (1,*) scheme
READ (1,*) init
CLOSE(1)
!format d'ecriture
OPEN (2,FILE='burger.ini')
READ (2,*) frmt
CLOSE (2)

c---- calculs: pas d'espace et param------------------------
deltax=1/FLOAT(k-1)
s=nu*deltat/(deltax*deltax)
r = deltat/deltax

c---- allocation de la memoire------------------------------
ALLOCATE (mat(5*k),u0(k),d(k))

c**** initialisation des coefficients **********************
!crank nicholson fin. diff.
IF (scheme.EQ.1) THEN
delta=0
q=0
name_scheme=' fd.'
!crank nicholson fin. elmt.
ELSE IF(scheme.EQ.2) THEN
delta=1.0/6.0
q=0
name_scheme=' fem.'
!crank nicholson mass operator
ELSE IF(scheme.EQ.3) THEN
delta=0.12
q=0
name_scheme=' massop'
!crank nicholson 4pt. upwind
ELSE IF(scheme.EQ.4) THEN
delta=0
q=0.5
name_scheme=' 4pt-up'
END IF

c**** solution initiale ************************************
c---- calcul du point zero----------------------------------
klim=(k-1)/2+1

c---- initialisation sinusoidale ---------------------------
IF (init.EQ.1) THEN
name_init=' sinus'
!conditions initiales (pts: 1 à k)
DO j=1,k
u0(j)=SIN(2*PI*(j-1)*deltax)
END DO
!conditions aux limites (rigides)
u0(1)=0
u0(k)=0
!point zero
u0(klim)=0

c---- initialisation droite --------------------------------
ELSEIF (init.EQ.2) THEN
name_init=' droite'
!conditions initiales (pts: 1 à k)
DO j=1,k
u0(j)=-2*(j-1)*deltax + 1
END DO
!conditions aux limites (libres)
u0(1)=1
u0(k)=-1
!point zero
u0(klim)=0
END IF

c**** preparation et ecriture dans fichier resultats *******
compteur_time = 1
OPEN (10,FILE='output.txt')
WRITE (10,20) 'visc. :',nu
WRITE (10,20) 't_step:',deltat
WRITE (10,30) 'n_iter:',n
WRITE (10,30) 'p_sort:',pas_sortie
WRITE (10,30) 'n_pts :',k
WRITE (10,40) 'schema:',name_scheme
WRITE (10,40) 'init :',name_init
WRITE (10,*) ' '
WRITE (10,frmt) 0.0,((i-1)*deltax,i=1,k)
WRITE (10,frmt) 0.0,(u0(j),j=1,k)

c***********************************************************
c noyau du programme *
c***********************************************************
c**** boucle en temps **************************************
DO t=1,n

c**** systeme AX=B *****************************************

c-----initialisation de la matrice A------------------------
!CL amont-----------------------------------------
!sinusoidale
IF (init.EQ.1) THEN
j=1
mat(5*j-4)=0
mat(5*j-3)=0
mat(5*j-2)=1-2*delta+s+(0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)
mat(5*j) =0
j=2
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s-(0.25+0.5*q)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)
mat(5*j) =0

!droite
ELSEIF (init.EQ.2) THEN
j=1
mat(5*j-4)=0
mat(5*j-3)=0
mat(5*j-2)=1-delta+0.5*s+(q/6-0.25)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)
mat(5*j) =0
j=2
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s-(0.25+q/3)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)
mat(5*j) =0
END IF

!initialisation sol. positive---------------------
DO j=3, klim-1
mat(5*j-4)=(q/6.0)*r*u0(j-2)
mat(5*j-3)=delta-0.5*s-(0.25+0.5*q)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)
mat(5*j) =0
END DO

!initialisation point zero------------------------
j=klim
mat(5*j-4)=(q/6.0)*r*u0(j-2)
mat(5*j-3)=0.5*(delta-0.5*s-(0.25+0.5*q)*r*u0(j-1)+
& delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1))
mat(5*j-2)=0.5*(1-2*delta+s+(0.5*q)*r*u0(j)+
& 1-2*delta+s+(-0.5*q)*r*u0(j))
mat(5*j-1)=0.5*(delta-0.5*s+(0.25-q/6.0)*r*u0(j+1)+
& delta-0.5*s+(0.25+0.5*q)*r*u0(j+1))
mat(5*j-0)=-(q/6.0)*r*u0(j+2)

!initialisation sol. negative---------------------
DO j=klim+1, k-2
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(-0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25+0.5*q)*r*u0(j+1)
mat(5*j-0)=-(q/6.0)*r*u0(j+2)
END DO

!CL aval------------------------------------------
!sinusoidale
IF (init.EQ.1) THEN
j=k-1
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(-0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25+0.5*q)*r*u0(j+1)
mat(5*j-0)=0
j=k
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(-0.5*q)*r*u0(j)
mat(5*j-1)=0
mat(5*j-0)=0

!droite
ELSEIF (init.EQ.2) THEN
j=k-1
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1)
mat(5*j-2)=1-2*delta+s+(-0.5*q)*r*u0(j)
mat(5*j-1)=delta-0.5*s+(0.25+q/3)*r*u0(j+1)
mat(5*j-0)=0
j=k
mat(5*j-4)=0
mat(5*j-3)=delta-0.5*s+(-0.25+q/6.0)*r*u0(j-1)
mat(5*j-2)=1-delta+0.5*s+(0.25-q/6)*r*u0(j)
mat(5*j-1)=0
mat(5*j-0)=0

END IF

c-----initialisation du vecteur B---------------------------
!sinus
IF (init.EQ.1) THEN
d(1)=(delta+0.5*s)*u0(2)+(1-2*delta-s)*u0(1)
DO j=2,k-1
d(j)=(delta+0.5*s)*(u0(j-1)+u0(j+1))
& +(1-2*delta-s)*u0(j)
END DO
d(k)=(delta+0.5*s)*u0(k-1)+(1-2*delta-s)*u0(k)

!droite
ELSEIF (init.EQ.2) THEN
d(1)=(delta+0.5*s)*u0(2)+(1-delta-0.5*s)*u0(1)
DO j=2,k-1
d(j)=(delta+0.5*s)*(u0(j-1)+u0(j+1))
& +(1-2*delta-s)*u0(j)
END DO
d(k)=(delta+0.5*s)*u0(k-1)+(1-delta-0.5*s)*u0(k)

END IF

c---- resolution du systeme---------------------------------
CALL GAUSS(mat,d,k,2) !voir subroutine gauss

c**** postprocessing ***************************************
c---- sauvegarde de la solution-----------------------------
!affectation
DO j=1,k
u0(j)=d(j)
END DO

c---- integration des conditions limites
! sinusoidale
IF (init.EQ.1) THEN
u0(1)=0
u0(k)=0
! droite
ELSEIF (init.EQ.2) THEN
u0(1)=1
u0(k)=-1
END IF

c---- ecriture de la solution-------------------------------
cur_time= t / (compteur_time * pas_sortie)
IF (cur_time.EQ.1) THEN
WRITE(10,frmt) t*deltat,(u0(j),j=1,k)
compteur_time = compteur_time + 1
END IF

c**** fin de la boucle en temps ****************************
END DO

c***********************************************************
c terminaison *
c***********************************************************
c---- fermeture du fichier ---------------------------------
CLOSE (10)

c---- liberation de la memoire
DEALLOCATE (mat,u0,d)

STOP

c---- definition des formats d'ecritue----------------------
10 FORMAT(F10.8)
20 FORMAT(A7,F7.4)
30 FORMAT(A7,I7)
40 FORMAT(A7,A7)
END