CFD Online Discussion Forums - WENO Code (1D Advection Equation)

Hi Everyone,

I have been trying to implement simple 5th order WENO scheme for a 1-D advection equation in a instantaneous spill situation ( my background is Groundwater Modeling). Therefore I assume that at t=0, C0=1 at cell 1 and for all other time its 0.

I have problems defining the boundary stencils and hence my code does not converge.. Could anyone help me to debug the code or send some suggestions!!

Carolyn

Fortran Code Below

!================================================= ================== ! Program to solve 1D advection equation using 5th order central WENO scheme ! dc/dt + udc/dx=0 C is the concentration and u is the !velocity ! Boundary Condition C(x,t0)=1 x=x0 and t=t0 ! Its a pule of conc C0 with velocity u !================================================= ==================

PROGRAM WENO1

IMPLICIT NONE

INTEGER I,J,Nmax,IT,ITN_max,IC_FLAG

PARAMETER (Nmax=1001,IC_FLAG=1)

REAL *8 U(Nmax),Phi(Nmax),DCDX(Nmax)

REAL *8 Lx,Dx,CFL,Dt,X0,U0,Tlast

REAL *8 Phi_1(Nmax),Phi_2(Nmax),Phi_3(Nmax)

REAL *8 PL(Nmax)

!=============Input variables====================================== ! OPEN(UNIT=1,FILE='Input.txt') ! READ (1,*) X0,Lx,Tlast,CFL,U0 ! CLOSE(UNIT=1) X0=0 Lx=10 Tlast=5 CFL=0.2 U0=0.8

!================================================= ================

!================================================= ================== ! Variable Initialization ! !================================================= ==================

Dx = (Lx-X0)/REAL(Nmax-1)

Dt = (Dx*CFL/U0)

!Dt=(Dx)**(5/3);

ITN_max = INT(TLast/Dt)

CALL INITIALIZE(Nmax,X0,Dx,U,U0,Phi)

OPEN(UNIT=2,FILE='Initial_Condition.dat',FORM='for matted')

DO I = 1, Nmax

WRITE(2,*)X0+(I-1)*Dx,Phi(I)

END DO

CLOSE(2)

!================================================= =========== ! Start of Time Loop ! Time Discretization done using 3rd order Runge-Kutta Method !================================================= ===========

WRITE(*,*)ITN_max

DO IT = 1, ITN_max

CALL WENO(Nmax,U,Phi,DCDX,Dx)

DO I = 1, Nmax

Phi_1(I) = Phi(I) + Dt*DCDX(I)

Phi_1(1)=0

END DO

CALL WENO(Nmax,U,Phi_1,DCDX,Dx)

DO I = 1, Nmax

Phi_2(I) = 3./4.*Phi(I) + 1./4.*Phi_1(I)+ 1./4.*Dt*DCDX(I)

Phi_2(1)=0

END DO

CALL WENO(Nmax,U,Phi_2,DCDX,Dx)

DO I = 1, Nmax

Phi_3(I) = 1./3.*Phi(I) + 2./3.*Phi_2(I)+ 2./3.*Dt*DCDX(I)

Phi_3(1)=0

END DO

!

DO I = 1, Nmax

Phi(I) = Phi_3(I)

END DO

END DO

!============= End of Time Loop============================================== ==========

OPEN(UNIT=3,FILE='SOLUTION.dat',FORM='formatted')

DO I = 1, Nmax

WRITE(3,*)X0+(I-1)*Dx,Phi(I)

END DO

CLOSE(3)

!STOP

END

!================================================= ================== ! Subroutine to Initialize the function ! !================================================= ==================

SUBROUTINE INITIALIZE(Nmax,X0,Dx,U,U0,Phi)

IMPLICIT NONE

INTEGER Nmax, I, IR

REAL *8 U(Nmax),Phi(Nmax),Dx,x,X0,PI, U0

x = X0

DO I = 1, Nmax

U(I) = U0

Phi(I)=0.0

x = x + Dx

END DO

Phi(1)=1

RETURN

END !=================End of subroutine Initialize====================================

!================================================= ================================ ! Subroutine to solve 1st Spatial Derivative using 5th-Order Central WENO Scheme !================================================= ================================

SUBROUTINE WENO(Nmax,U,CW,DCDX,Dx)

IMPLICIT NONE

INTEGER Nmax, I, IL, IR, ILL, IRR

REAL *8 U(Nmax),CW(Nmax),DCDX(Nmax),Dx

REAL *8 ROESPD(NMAX)

REAL *8 EFLX_1(NMAX),EFLX_2(NMAX),EFLX_3(NMAX)

REAL *8 WFLX(NMAX)

REAL *8 EPS,d1,d2,d3,alpha1,alpha2,alpha3

REAL *8 beta1,beta2,beta3,w1,w2,w3,SUM_alpha

REAL *8 a1(8),a2(8)

EPS = 1.0E-6

DO I = 1, Nmax

IL = I -1

ILL = I-2

IR = I + 1

IRR = I + 2

IF (IR.GT.Nmax) IR = IR-Nmax+1

IF (IRR.GT.Nmax) IRR = IRR-Nmax+1

IF (IL.LT.1) IL=IL+Nmax-1

IF (ILL.LT.1) ILL=ILL+Nmax-1

EFLX_1(I) = (1.0/3.0)*CW(I) + (5.0/6.0)*CW(IR)- (1./6.)*CW(IRR)

EFLX_2(I) = (-1.0/6.0)*CW(IL) + (5.0/6.0)*CW(I) +(1.0/3.0)*CW(IR)

EFLX_3(I) = (1.0/3.0)*CW(ILL) - (7.0/6.0)*CW(IL) + (11.0/6.0)*CW(I)

beta1 = (13.0/12.0)*(CW(I)-2.0*CW(IR)+ CW(IRR))**2.0 +(1.0/4.0)*(3.0*CW(I)&

-4.0*CW(IR)+CW(IRR))**2.0

beta2 = (13.0/12.0)*(CW(IL)-2.0*CW(I)+CW(IR))**2.0 +(1.0/4.0)*(CW(IL)-CW(IR))**2

beta3 = (13.0/12.0)*(CW(ILL)-2.0*CW(IL)+ CW(I))**2.0 +(1.0/4.0)*(CW(ILL)&

-4.0*CW(IL)+ 3.0*CW(I))**2.0

d1 = 3.0/10.0

d2 = 3.0/5.0

d3 = 1.0/10.0

alpha1 = d1/(EPS+beta1)**2.0

alpha2 = d2/(EPS+beta2)**2.0

alpha3 = d3/(EPS+beta3)**2.0

SUM_alpha = alpha1 + alpha2 + alpha3

w1 = alpha1/SUM_alpha

w2 = alpha2/SUM_alpha

w3 = alpha3/SUM_alpha

WFLX(I) = w1*EFLX_1(I) + w2*EFLX_2(I) + w3*EFLX_3(I)

END DO

DO I = 1, Nmax

IL = I - 1

IR = I + 1

IF (IR.GT.Nmax) IR=IR-Nmax+1

IF (IRR.GT.Nmax) IRR=IRR-Nmax+1

IF (IL.LT.1) IL=IL+Nmax-1

IF (ILL.LT.1)ILL=ILL+Nmax-1

DCDX(I) = -U(I)*(WFLX(I)-WFLX(IL))/DX

DCDX(Nmax)=0

END DO

!================================================= ============ ! Boundary Stencil Definition ! Since boundary stencil (1,2,3 and Nmax-1,Nmax) not !possible to define using 5 point stencil !================================================= ============

DO I=2,3

IRR=I+2

IR=I+1

DCDX(I)=-U(I)*(CW(I+1)-CW(I))/(DX)

END DO

DO I=Nmax-2,Nmax-1

IL=I-1

ILL=I-2

DCDX(I)=-U(I)*(CW(I)-CW(I-1))/(DX)

END DO

DCDX(Nmax)=0

DCDX(1)= -U(1)*(CW(2)-CW(1))/(DX)

RETURN

END

!=================End of Subroutine WENO========================== ! ! !================================================= ==================