Program “burger.for”

Contents | Listing | Download

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.

UP

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 (m2/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.

UP

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-MXICODE=1MX2=2*MX+1MX1=MX+1M22=MX2-1L=MX1DO K=1,LVDIVB=1./A(L)LMX=L+MXIF(K .GT. LMOI) THENLMX=LMX-ICODEICODE=ICODE+1ENDIFDO J=L,LMXA(J)=A(J)*DIVBENDDOB(K)=B(K)*DIVBIF(K.NE.LV) THENMI22=0DO II=ICODE,MXMI22=MI22+M22LMI22=L+MI22IF (A(LMI22).NE.0) THENCOEFF=A(LMI22)DO J=L,LMXJMI22=J+MI22A(JMI22)=A(JMI22)-A(J)*COEFFENDDOIII=II-ICODE+1KIII=K+IIIB(KIII)=B(KIII)-B(K)*COEFFENDIFENDDOL=L+MX2ENDIFENDDOICODE=0I=LVL=MX1+1+(LV-1)*MX2DO WHILE (I.GT.1)SIGMA=0.L=L-MX2LMX=L+MX-1IF(I.GT.LMOI) THENLMX=L+ICODEICODE=ICODE+1ENDIFII=IDO J=L,LMXSIGMA=SIGMA+A(J)*B(II)II=II+1ENDDOI=I-1B(I)=B(I)-SIGMAENDDORETURNEND

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

UP

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 1DIMPLICIT NONE c***********************************************************c declaration des variables *c***********************************************************!variables de boucleINTEGER i, j, t!declaration des constantesREAL PI !pas de temps et d'espace, visc., param. num.INTEGER n, k !n_iter, n_espaceREAL deltax, deltat, s ,rREAL q, delta !param num schemaREAL nu !visc. cin.!choix du schema et de l'initialisationINTEGER klim !indce de changement de sgnINTEGER scheme !var de schemaINTEGER init !var d'initialisationCHARACTER*10 name_schemeCHARACTER*10 name_init!variables de sortie des resultatsINTEGER pas_sortie, compteur_time, cur_timeCHARACTER*2000 frmt!definition des matrices et vecteursREAL mat,d,u0ALLOCATABLE mat(:),d(:),u0(:)c definition des constantes ********************************PI = 3.14159265359c***********************************************************c initialisation *c***********************************************************c**** lecture des donnees **********************************!donneesOPEN (1,FILE='input.txt')READ (1,*) nuREAD (1,*) deltatREAD (1,*) nREAD (1,*) pas_sortieREAD (1,*) kREAD (1,*) schemeREAD (1,*) initCLOSE(1)!format d'ecritureOPEN (2,FILE='burger.ini')READ (2,*) frmtCLOSE (2)c---- calculs: pas d'espace et param------------------------deltax=1/FLOAT(k-1)s=nu*deltat/(deltax*deltax)r = deltat/deltaxc---- allocation de la memoire------------------------------ALLOCATE (mat(5*k),u0(k),d(k))c**** initialisation des coefficients **********************!crank nicholson fin. diff.IF (scheme.EQ.1) THENdelta=0q=0name_scheme=' fd.'!crank nicholson fin. elmt.ELSE IF(scheme.EQ.2) THENdelta=1.0/6.0q=0name_scheme=' fem.'!crank nicholson mass operatorELSE IF(scheme.EQ.3) THENdelta=0.12q=0name_scheme=' massop'!crank nicholson 4pt. upwindELSE IF(scheme.EQ.4) THENdelta=0q=0.5name_scheme=' 4pt-up'END IFc**** solution initiale ************************************c---- calcul du point zero----------------------------------klim=(k-1)/2+1c---- initialisation sinusoidale ---------------------------IF (init.EQ.1) THENname_init=' sinus'!conditions initiales (pts: 1 à k)DO j=1,ku0(j)=SIN(2*PI*(j-1)*deltax)END DO!conditions aux limites (rigides)u0(1)=0u0(k)=0!point zerou0(klim)=0c---- initialisation droite --------------------------------ELSEIF (init.EQ.2) THENname_init=' droite'!conditions initiales (pts: 1 à k)DO j=1,ku0(j)=-2*(j-1)*deltax + 1END DO!conditions aux limites (libres)u0(1)=1u0(k)=-1!point zerou0(klim)=0END IFc**** preparation et ecriture dans fichier resultats *******compteur_time = 1OPEN (10,FILE='output.txt')WRITE (10,20) 'visc. :',nuWRITE (10,20) 't_step:',deltatWRITE (10,30) 'n_iter:',nWRITE (10,30) 'p_sort:',pas_sortieWRITE (10,30) 'n_pts :',kWRITE (10,40) 'schema:',name_schemeWRITE (10,40) 'init :',name_initWRITE (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,nc**** systeme AX=B *****************************************c-----initialisation de la matrice A------------------------!CL amont-----------------------------------------!sinusoidaleIF (init.EQ.1) THENj=1mat(5*j-4)=0mat(5*j-3)=0mat(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) =0j=2mat(5*j-4)=0mat(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!droiteELSEIF (init.EQ.2) THENj=1mat(5*j-4)=0mat(5*j-3)=0mat(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) =0j=2mat(5*j-4)=0mat(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) =0END IF!initialisation sol. positive---------------------DO j=3, klim-1mat(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) =0END DO!initialisation point zero------------------------j=klimmat(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-2mat(5*j-4)=0mat(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------------------------------------------!sinusoidaleIF (init.EQ.1) THENj=k-1mat(5*j-4)=0mat(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)=0j=kmat(5*j-4)=0mat(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)=0mat(5*j-0)=0!droiteELSEIF (init.EQ.2) THENj=k-1mat(5*j-4)=0mat(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)=0j=kmat(5*j-4)=0mat(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)=0mat(5*j-0)=0END IFc-----initialisation du vecteur B---------------------------!sinusIF (init.EQ.1) THENd(1)=(delta+0.5*s)*u0(2)+(1-2*delta-s)*u0(1)DO j=2,k-1d(j)=(delta+0.5*s)*(u0(j-1)+u0(j+1))& +(1-2*delta-s)*u0(j)END DOd(k)=(delta+0.5*s)*u0(k-1)+(1-2*delta-s)*u0(k)!droiteELSEIF (init.EQ.2) THENd(1)=(delta+0.5*s)*u0(2)+(1-delta-0.5*s)*u0(1)DO j=2,k-1d(j)=(delta+0.5*s)*(u0(j-1)+u0(j+1))& +(1-2*delta-s)*u0(j)END DOd(k)=(delta+0.5*s)*u0(k-1)+(1-delta-0.5*s)*u0(k)END IFc---- resolution du systeme---------------------------------CALL GAUSS(mat,d,k,2) !voir subroutine gaussc**** postprocessing ***************************************c---- sauvegarde de la solution-----------------------------!affectationDO j=1,ku0(j)=d(j)END DOc---- integration des conditions limites! sinusoidaleIF (init.EQ.1) THENu0(1)=0u0(k)=0! droiteELSEIF (init.EQ.2) THENu0(1)=1u0(k)=-1END IFc---- ecriture de la solution-------------------------------cur_time= t / (compteur_time * pas_sortie)IF (cur_time.EQ.1) THENWRITE(10,frmt) t*deltat,(u0(j),j=1,k)compteur_time = compteur_time + 1END IFc**** fin de la boucle en temps ****************************END DOc***********************************************************c terminaison *c***********************************************************c---- fermeture du fichier ---------------------------------CLOSE (10)c---- liberation de la memoireDEALLOCATE (mat,u0,d)STOPc---- definition des formats d'ecritue----------------------10 FORMAT(F10.8)20 FORMAT(A7,F7.4)30 FORMAT(A7,I7)40 FORMAT(A7,A7)END

UP