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

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

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

Reply
 
LinkBack Thread Tools Display Modes
Old   April 26, 2007, 20:03
Default what's wrong about my code for 2d burgers equation
  #1
morxio
Guest
 
Posts: n/a
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:

  Reply With Quote

Old   April 26, 2007, 20:04
Default Re: what's wrong about my code for 2d burgers equa
  #2
morxio
Guest
 
Posts: n/a
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

  Reply With Quote

Old   April 26, 2007, 23:50
Default Re: what's wrong about my code for 2d burgers equa
  #3
Harish
Guest
 
Posts: n/a
Did you test it for 1-d burgers equation with an initial discontinuity ? Did you validate your code for continuous problems?

  Reply With Quote

Old   April 27, 2007, 10:38
Default Re: what's wrong about my code for 2d burgers equa
  #4
morxio
Guest
 
Posts: n/a
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.
  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
CEL code for simulating the equation of motion of a vibrating rigid body mohamadaliv CFX 11 October 19, 2011 04:30
Viscosity and the Energy Equation Rich Main CFD Forum 0 December 16, 2009 15:01
solution to Burger's Equation using fortran Tony Limjuco Main CFD Forum 4 April 2, 2005 00:41
exact solution of burger's equation sajar Main CFD Forum 9 March 4, 2004 05:55
What kind of Cmmercial CFD code you feel well? Lans Main CFD Forum 13 October 27, 1998 11:20


All times are GMT -4. The time now is 13:11.