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

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

 April 26, 2007, 20:03 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:

 April 26, 2007, 20:04 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

 April 26, 2007, 23:50 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?

 April 27, 2007, 10:38 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.