# 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.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post mohamadaliv CFX 11 October 19, 2011 04:30 Rich Main CFD Forum 0 December 16, 2009 15:01 Tony Limjuco Main CFD Forum 4 April 2, 2005 00:41 sajar Main CFD Forum 9 March 4, 2004 05:55 Lans Main CFD Forum 13 October 27, 1998 11:20

All times are GMT -4. The time now is 12:51.