what's wrong about my code for 2d burgers equation
 

Thank you for your read firstly!

To simulate the the 2 dimensional Burgers equation with initial condition:

u_t + (u*u/2.0)_x + (u*u/2.0)_y =0. The initial condition is:u(0,x,y)=-1.0, 0.5, -0.2, 0.8 at the four quadrant respectively.

I use 5order weno + 3order Runge-Kutta + Lax-Friedrich flux splitting.

But through running in Fortran6.5 and plotting in tecplot10, the figure shows that the four piece of initial value don't work properly, which means that the figure just reflects the distribution of initial value in 3D.

Who could give my some suggestion or advance! If you have time can you send your code to my email box: morxio@hotmail.com ?

Thank you in advance!

my code:

Re: what's wrong about my code for 2d burgers equa
 

C C C======================================== C=============FUNCTION f(x)============== C========================================

function f(x)

implicit double precision (a-h,o-z)

double precision :: f

f=x*x/2.0

end C======================================== C=============FUNCTION fp(x)============= C========================================

function fp(x)

implicit double precision (a-h,o-z)

double precision :: fp

fp=0.5*(f(x)+x)

end C======================================== C=============FUNCTION fm(x)============= C========================================

function fm(x)

implicit double precision (a-h,o-z)

double precision :: fm

fm=0.5*(f(x)-x)

end C======================================== C=============FUNCTION g(x)============== C========================================

function g(x)

implicit double precision (a-h,o-z)

double precision :: g

g=x*x/2.0

end C======================================== C=============FUNCTION gp(x)============= C========================================

function gp(x)

implicit double precision (a-h,o-z)

double precision :: gp

gp=0.5*(g(x)+x)

end C======================================== C=============FUNCTION gm(x)============= C========================================

function gm(x)

implicit double precision (a-h,o-z)

double precision :: gm

gm=0.5*(g(x)-x)

end C------------------------------------------------------------------------- C----------------------------------SUB1----------------------------------- C-------------------------------------------------------------------------

subroutine sub1(u,ftemp,gtemp,flux,glux,up1_2p,up1_2m,tu,

& dx,dy,nx,ighost)

implicit double precision (a-h,o-z)

double precision u(-ighost:ighost+nx,-ighost:ighost+nx)

double precision up1_2m(1:nx),up1_2p(1:nx),tu(-ighost:ighost+nx)

double precision flux(-ighost:ighost+nx)

double precision glux(-ighost:ighost+nx)

double precision ftemp(-ighost:ighost+nx,-ighost:ighost+nx)

double precision gtemp(-ighost:ighost+nx,-ighost:ighost+nx)

double precision fluxm(1:nx)

double precision fluxp(1:nx)

double precision gluxm(1:nx)

double precision gluxp(1:nx)

double precision,external :: f,g,gm,fm,gp,fp

do j=0,nx

do i=-ighost,ighost+nx

tu(i)=fm(u(i,j))

enddo

call reconstruction(tu,nx,up1_2m,up1_2p,ighost)

do i=1,nx

fluxm(i)=up1_2p(i)

enddo

enddo

do j=0,nx

do i=-ighost,ighost+nx

tu(i)=fp(u(i,j))

enddo

call reconstruction(tu,nx,up1_2m,up1_2p,ighost)

do i=1,nx

fluxp(i)=up1_2m(i)

enddo

enddo

do i=1,nx

flux(i)=fluxm(i)+fluxp(i)

ftemp(i,j)=flux(i)-flux(i-1)

enddo C===========================================

do i=0,nx

do j=-ighost,ighost+nx

tu(j)=fm(u(i,j))

enddo

call reconstruction(tu,nx,up1_2m,up1_2p,ighost)

do j=1,nx

gluxm(j)=up1_2p(j)

enddo

enddo

do i=0,nx

do j=-ighost,ighost+nx

tu(j)=fp(u(i,j))

enddo

call reconstruction(tu,nx,up1_2m,up1_2p,ighost)

do j=1,nx

gluxp(j)=up1_2m(j)

enddo

enddo

do j=1,nx

glux(j)=gluxm(j)+gluxp(j)

gtemp(i,j)=glux(j)-glux(j-1)

enddo

end C------------------------------------------------------------------------- C-------------------------PERIOD VALUE------------------------------------ C-------------------------------------------------------------------------

subroutine period(u,nx,ighost)

implicit double precision (a-h,o-z)

double precision u(-ighost:nx+ighost,-ighost:nx+ighost)

do j=-ighost,nx+ighost

do i=0,-ighost,-1

u(i,j)=u(i+nx,j)

enddo

do i=nx+1,nx+ighost

u(i,j)=u(i-nx,j)

enddo

enddo

do i=-ighost,nx+ighost

do j=0,-ighost,-1

u(i,j)=u(i,j+nx)

enddo

do j=nx+1,nx+ighost

u(i,j)=u(i,j-nx)

enddo

enddo

end C------------------------------------------------------------------------- C-------------------------Reconstruction---------------------------------- C-------------------------------------------------------------------------

subroutine reconstruction(u,nx,up1_2m,up1_2p,ighost)

implicit double precision (a-h,o-z)

double precision u(-ighost:nx+ighost)

double precision up1_2m(1:nx),up1_2p(1:nx)

eps=1.0D-10

do i=1,nx

d1=1.0/10

d2=6.0/10

d3=3.0/10

u1=1.0*u(i-2)/3.0-7.0*u(i-1)/6.0+11.0*u(i)/6.0

u2=-1.0*u(i-1)/6.0+5.0*u(i)/6.0+1.0*u(i+1)/3.0

u3=1.0*u(i)/3.0+5.0*u(i+1)/6.0-1.0*u(i+2)/6.0

t1=u(i-2)-2.0*u(i-1)+u(i)

t2=u(i-2)-4.0*u(i-1)+3.0*u(i)

b1=13.0*t1*t1/3.0+t2*t2

t1=u(i-1)-2.0*u(i)+u(i+1)

t2=u(i+1)-u(i-1)

b2=13.0*t1*t1/3.0+t2*t2

t1=u(i)-2.0*u(i+1)+u(i+2)

t2=3.0*u(i)-4.0*u(i+1)+u(i+2)

b3=13.0*t1*t1/3.0+t2*t2

a1=d1/(eps+b1)/(eps+b1)

a2=d2/(eps+b2)/(eps+b2)

a3=d3/(eps+b3)/(eps+b3)

aa=a1+a2+a3

w1=a1/aa

w2=a2/aa

w3=a3/aa

up1_2p(i)=w1*u1+w2*u2+w3*u3

d1=3.0/10.0

d2=6.0/10.0

d3=1.0/10.0

u_1=-1.0*u(i-1)/6.0+5.0*u(i)/6.0+1.0*u(i+1)/3.0

u_2=1.0*u(i)/3.0+5.0*u(i+1)/6.0-1.0*u(i+2)/6.0

u_3=11.0*u(i+1)/6.0-7.0*u(i+2)/6.0+1.0*u(i+3)/3.0

t1=u(i-1)-2.0*u(i)+u(i+1)

t2=u(i-1)-4.0*u(i)+3.0*u(i+1)

b1=13.0*t1*t1/3.0+t2*t2

t1=u(i)-2.0*u(i+1)+u(i+2)

t2=u(i+2)-u(i)

b2=13.0*t1*t1/3.0+t2*t2

t1=u(i+1)-2.0*u(i+2)+u(i+3)

t2=3.0*u(i+1)-4.0*u(i+2)+u(i+3)

b3=13.0*t1*t1/3.0+t2*t2

a_1=d3/(eps+b1)/(eps+b1)

a_2=d2/(eps+b2)/(eps+b2)

a_3=d1/(eps+b3)/(eps+b3)

a_a=a_1+a_2+a_3

w_1=a_1/a_a

w_2=a_2/a_a

w_3=a_3/a_a

up1_2m(i)=w_1*u_1+w_2*u_2+w_3*u_3

enddo

end C------------------------------------------------------------------------- C-----------------------------MAIN---------------------------------------- C-------------------------------------------------------------------------

program main

implicit double precision (a-h,o-z)

parameter (nx=80,ighost=6)

double precision pu(-ighost:nx+ighost,-ighost:nx+ighost)

double precision su(-ighost:nx+ighost,-ighost:nx+ighost)

double precision u1(-ighost:nx+ighost,-ighost:nx+ighost)

double precision u2(-ighost:nx+ighost,-ighost:nx+ighost)

double precision tu(-ighost:nx+ighost)

double precision up1_2m(0:nx),up1_2p(0:nx)

double precision ftemp(-ighost:ighost+nx,-ighost:ighost+nx)

double precision gtemp(-ighost:ighost+nx,-ighost:ighost+nx)

double precision x(-ighost:nx+ighost),y(-ighost:nx+ighost)

double precision u0(-ighost:nx+ighost,-ighost:nx+ighost)

dx=0.025

dy=0.025

dt=0.000625

open(unit=1,file='d:\data.plt',status='unknown')

write(*,*) 'This is a program for two dimensional Burgers equation

&with four initial value.'

write(*,*) 'Please input the maxium time T='

read(*,*) supt

do i=-ighost,nx+ighost

x(i)=(i-1)*dx+dx/2.0

y(i)=(i-1)*dy+dy/2.0

enddo

do i=-ighost,ighost+nx

do j=-ighost, ighost+nx

if (x(i).gt.1.0) then

if (y(j).gt.1.0) then

u0(i,j)=-1.0

else

u0(i,j)=0.8

endif

else

if (y(j).le.1.0) then

u0(i,j)=0.5

else

u0(i,j)=-0.2

endif

endif

enddo

enddo

do i=-ighost,ighost+nx

do j=-ighost, ighost+nx

pu(i,j)=u0(i,j)

enddo

enddo

C-=======================================

t=0

do 100 while(t.lt.supt)

call sub1(pu,ftemp,gtemp,flux,glux,up1_2p,up1_2m,tu,

& dx,dy,nx,ighost)

do i=1,nx

do j=1,nx

u1(i,j)=pu(i,j)-dt*(ftemp(i,j)/dx+gtemp(i,j)/dy)

enddo

enddo

call period(u1,nx,ighost) C-------------------------------------------

call sub1(u1,ftemp,gtemp,flux,glux,up1_2p,up1_2m,tu,

& dx,dy,nx,ighost)

do i=1,nx

do j=1,nx

u2(i,j)=3.0/4.0*pu(i,j)+1.0/4.0*u1(i,j)-

& 1.0/4.0*dt*(ftemp(i,j)/dx+gtemp(i,j)/dy)

enddo

enddo

call period(u2,nx,ighost) C-------------------------------------------

call sub1(u2,ftemp,gtemp,flux,glux,up1_2p,up1_2m,tu,

& dx,dy,nx,ighost)

do i=1,nx

do j=1,nx

t1=1.0/3

t2=2.0/3

su(i,j)=t1*pu(i,j)+t2*u2(i,j)

& -t2*dt*(ftemp(i,j)/dx+gtemp(i,j)/dy)

enddo

enddo

call period(su,nx,ighost) C-------------------------------------------

do i=-ighost,nx+ighost

do j=-ighost,nx+ighost

pu(i,j)=su(i,j)

enddo

enddo

t=t+dt

if(t+dt.gt.supt) then

dt=supt-t+1

endif

write(*,*) t 100 continue

write(1,*) 'TITLE = "EXAMPLE: 3D GEOMETRIES" '

write(1,*) 'VARIABLES = "X", "Y", "Z"'

write(1,*) 'ZONE T="Floor", I= ',nx,'J= ',nx,'F=POINT'

do i=1,nx

do j=1,nx

write(1,101) x(i),y(j),su(i,j)

enddo

enddo 101 format(1x,e20.10,e20.10,e20.10)

end

Re: what's wrong about my code for 2d burgers equa
 

Did you test it for 1-d burgers equation with an initial discontinuity ? Did you validate your code for continuous problems?

Re: what's wrong about my code for 2d burgers equa
 

It is impossible for the 1d conservation law equation with 4 pieces of initial value in four quadrant respectively.

And the Weno scheme with Runge-kutta and LF flux splittig works wonderful.