CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

WENO Code (1D Advection Equation)

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   February 19, 2007, 17:09
Default WENO Code (1D Advection Equation)
  #1
Carolyn
Guest
 
Posts: n/a
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========================== ! ! !================================================= ==================
  Reply With Quote

Old   February 20, 2007, 01:54
Default Re: WENO Code (1D Advection Equation)
  #2
rt
Guest
 
Posts: n/a
valid code, with eno, venleer, roe methods (not weno, weno is similar with eno), is avalible, if u need give ur email
  Reply With Quote

Old   February 20, 2007, 02:03
Default Re: WENO Code (1D Advection Equation)
  #3
Carolyn
Guest
 
Posts: n/a
My email is carolyn.msu@gmail.com

Thanks

Carolyn
  Reply With Quote

Old   February 20, 2007, 05:06
Default Re: WENO Code (1D Advection Equation)
  #4
idee
Guest
 
Posts: n/a
I need it please forward to me:

research2phase@yahoo.com

Your help very much appreciated.

  Reply With Quote

Old   February 22, 2007, 22:49
Default Re: WENO Code (1D Advection Equation)
  #5
yang, x.h.
Guest
 
Posts: n/a
I need it please forward to me:

yangxh@ynao.ac.cn

Thank you.

  Reply With Quote

Old   March 10, 2007, 04:45
Default Re: WENO Code (1D Advection Equation)
  #6
dusky.he
Guest
 
Posts: n/a
I also need it , could you forward to me? qunwuhe@gmail.com

Thanks
  Reply With Quote

Old   March 11, 2007, 14:21
Default Re: WENO Code (1D Advection Equation)
  #7
Shiv
Guest
 
Posts: n/a
Could you please forward the code to me too. I am using ENO and it will help me in debuggin my code. My email id is shiv4u@gmail.com Thank you in advance

  Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
advection equation venkataramana OpenFOAM 1 September 4, 2012 03:12
2D linear advection equation mcaro Main CFD Forum 2 February 11, 2011 13:30
orlandi code for poisson equation with FFT ( paid tutorial) HaKu Main CFD Forum 0 June 29, 2009 15:40
F77 Code for Energy Equation Carlos Main CFD Forum 0 April 2, 2002 14:57
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 15:56


All times are GMT -4. The time now is 01:10.