You need to include the following as a data file, "indata"
102
0.2
1.0
0.5
nx (maximum 102)
Tfinal
riemann data ul
ur
c burgers equation by godunov method
program burger1
c SOLVE du/dt +(1/2) du^2/dx = 0
c BY 1st ORDER GODUNOV
implicit double precision(a-h,o-z)
dimension u(102)
common/grid/x(102)
common/Igrid/nx
common/steps/dt,dx,dtdx
common/rdata/ul,ur
open(unit=1,file="indata")
read(1,*)nx
read(1,*)Tfinal
read(1,*)ul
read(1,*)ur
close(1)
c ONE CELL AT EACH END USED FOR RIEMANN DATA
dx=1.0d0/(nx-2)
do 01 i=1,nx
x(i)=(i-1-0.5)*dx
01 continue
c INITIAL DATA
t=0.0d0
do 10 i=1,51
u(i)=ul
10 continue
do 12 i=52,nx
u(i)=ur
12 continue
do 999 niter=1,100
c CFL ESTIMATE
sound=1.0
do i=1,nx
sound=dmax1(sound,u(i))
enddo
dt=0.90d0*dx/sound
tnew=t+dt
if (tnew .gt. Tfinal) then
dt=Tfinal-t
tnew=t+dt
endif
dtdx=dt/dx
call step(u)
t=tnew
if (t .ge. Tfinal) then
go to 1001
endif
999 continue
c WRITE DATA
1001 write(6,888)t
888 format('','FINAL TIME = ',f12.7)
open(unit=1,file='solution')
do 1100 i=1,nx
write(1,1066) x(i), u(i)
1100 continue
close(1)
1066 format('',f10.6,1x,f12.7)
stop
end
subroutine step(u)
implicit double precision(a-h,o-z)
dimension up(102),u(102),ulhalf(102),urhalf(102),
+ us(102)
common/grid/x(102)
common/Igrid/nx
common/steps/dt,dx,dtdx
common/rdata/ul,ur
c STATEMENT FUNCTIONS
sgn(w)=dsign(1.0d0,w)
f(z)=0.5*z*z
dfdu(w,z)=(f(z)-f(w))/(z-w)
c LEFT AND RIGHT PREDICTORS
c OUR NOTATION FOR LEFT & RIGHT STATES:
c
c --UR_{i-1}-+-UL_i-----UR_i-+-UL_{i+1}--
c
c ---+---------------+---------------+---
c i-1 i i+1
do i=1,nx-1
ulhalf(i)=u(i)
urhalf(i)=u(i+1)
enddo
c AT RIGHT BOUNDARY, RIEMANN DATA
ulhalf(nx)=ur
urhalf(nx)=ur
c STATIONARY STATE AT i+1/2
do 200 i=1,nx-1
c case 1 Trivial
if (urhalf(i) .eq. ulhalf(i)) then
us(i)=urhalf(i)
endif
c case 2 Shock
if (ulhalf(i) .gt. urhalf(i)) then
speed=0.5d0*(ulhalf(i)+urhalf(i))
if (speed .gt. 0.d0) us(i)=ulhalf(i)
if (speed .lt. 0.d0) us(i)=urhalf(i)
if (speed .eq. 0.d0) us(i)=0.0
endif
c case 3 Rarefaction
c suffers dog leg and entropy violation
if (ulhalf(i) .lt. urhalf(i)) then
tail=ulhalf(i)
head=urhalf(i)
if (tail .gt. 0.d0) then
us(i)=ulhalf(i)
else if (head .lt. 0.d0) then
us(i)=urhalf(i)
else
us(i)=0.0
endif
endif
200 continue
C AT THE BOUNDARY
us(nx)=ur
c CONSERVATIVE DIFFERENCE
do 250 i=2,nx
up(i)=u(i)-dtdx*(f(us(i))-f(us(i-1)))
250 continue
c RIEMANN DATA
up(1)=ul
do 300 i=1,nx
u(i)=up(i)
300 continue
return
end
/* For the c hackers, from netlib's f2c */
/* -- translated by f2c (version 19970219).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
doublereal x[102];
} grid_;
#define grid_1 grid_
struct {
integer nx;
} igrid_;
#define igrid_1 igrid_
struct {
doublereal dt, dx, dtdx;
} steps_;
#define steps_1 steps_
struct {
doublereal ul, ur;
} rdata_;
#define rdata_1 rdata_
/* Table of constant values */
static integer c__3 = 3;
static integer c__1 = 1;
static integer c__5 = 5;
/* burgers equation by godunov method */
/* Main program */ MAIN__()
{
/* Format strings */
static char fmt_888[] = "(\002\002,\002FINAL TIME = \002,f12.7)";
static char fmt_1066[] = "(\002\002,f10.6,1x,f12.7)";
/* System generated locals */
integer i__1;
doublereal d__1, d__2;
olist o__1;
cllist cl__1;
/* Builtin functions */
integer f_open(), s_rsle(), do_lio(), e_rsle(), f_clos(), s_wsfe(),
do_fio(), e_wsfe();
/* Subroutine */ int s_stop();
/* Local variables */
extern /* Subroutine */ int step_();
static doublereal tnew;
static integer i__;
static doublereal t, u[102];
static integer niter;
static doublereal sound, tfinal;
/* Fortran I/O blocks */
static cilist io___1 = { 0, 1, 0, 0, 0 };
static cilist io___2 = { 0, 1, 0, 0, 0 };
static cilist io___4 = { 0, 1, 0, 0, 0 };
static cilist io___5 = { 0, 1, 0, 0, 0 };
static cilist io___12 = { 0, 6, 0, fmt_888, 0 };
static cilist io___13 = { 0, 1, 0, fmt_1066, 0 };
/* SOLVE du/dt +(1/2) du^2/dx = 0 */
/* BY 1st ORDER GODUNOV */
o__1.oerr = 0;
o__1.ounit = 1;
o__1.ofnmlen = 6;
o__1.ofnm = "indata";
o__1.orl = 0;
o__1.osta = 0;
o__1.oacc = 0;
o__1.ofm = 0;
o__1.oblnk = 0;
f_open(&o__1);
s_rsle(&io___1);
do_lio(&c__3, &c__1, (char *)&igrid_1.nx, (ftnlen)sizeof(integer));
e_rsle();
s_rsle(&io___2);
do_lio(&c__5, &c__1, (char *)&tfinal, (ftnlen)sizeof(doublereal));
e_rsle();
s_rsle(&io___4);
do_lio(&c__5, &c__1, (char *)&rdata_1.ul, (ftnlen)sizeof(doublereal));
e_rsle();
s_rsle(&io___5);
do_lio(&c__5, &c__1, (char *)&rdata_1.ur, (ftnlen)sizeof(doublereal));
e_rsle();
cl__1.cerr = 0;
cl__1.cunit = 1;
cl__1.csta = 0;
f_clos(&cl__1);
/* ONE CELL AT EACH END USED FOR RIEMANN DATA */
steps_1.dx = 1. / (igrid_1.nx - 2);
i__1 = igrid_1.nx;
for (i__ = 1; i__ <= i__1; ++i__) {
grid_1.x[i__ - 1] = (i__ - 1 - (float).5) * steps_1.dx;
/* L1: */
}
/* INITIAL DATA */
t = 0.;
for (i__ = 1; i__ <= 51; ++i__) {
u[i__ - 1] = rdata_1.ul;
/* L10: */
}
i__1 = igrid_1.nx;
for (i__ = 52; i__ <= i__1; ++i__) {
u[i__ - 1] = rdata_1.ur;
/* L12: */
}
for (niter = 1; niter <= 100; ++niter) {
/* CFL ESTIMATE */
sound = (float)1.;
i__1 = igrid_1.nx;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MAX */
d__1 = sound, d__2 = u[i__ - 1];
sound = max(d__1,d__2);
}
steps_1.dt = steps_1.dx * .9 / sound;
tnew = t + steps_1.dt;
if (tnew > tfinal) {
steps_1.dt = tfinal - t;
tnew = t + steps_1.dt;
}
steps_1.dtdx = steps_1.dt / steps_1.dx;
step_(u);
t = tnew;
if (t >= tfinal) {
goto L1001;
}
/* L999: */
}
/* WRITE DATA */
L1001:
s_wsfe(&io___12);
do_fio(&c__1, (char *)&t, (ftnlen)sizeof(doublereal));
e_wsfe();
o__1.oerr = 0;
o__1.ounit = 1;
o__1.ofnmlen = 8;
o__1.ofnm = "solution";
o__1.orl = 0;
o__1.osta = 0;
o__1.oacc = 0;
o__1.ofm = 0;
o__1.oblnk = 0;
f_open(&o__1);
i__1 = igrid_1.nx;
for (i__ = 1; i__ <= i__1; ++i__) {
s_wsfe(&io___13);
do_fio(&c__1, (char *)&grid_1.x[i__ - 1], (ftnlen)sizeof(doublereal));
do_fio(&c__1, (char *)&u[i__ - 1], (ftnlen)sizeof(doublereal));
e_wsfe();
/* L1100: */
}
cl__1.cerr = 0;
cl__1.cunit = 1;
cl__1.csta = 0;
f_clos(&cl__1);
s_stop("", 0L);
} /* MAIN__ */
/* Subroutine */ int step_(u)
doublereal *u;
{
/* System generated locals */
integer i__1;
/* Local variables */
static doublereal head, tail;
static integer i__;
static doublereal speed, up[102], us[102], ulhalf[102], urhalf[102];
/* STATEMENT FUNCTIONS */
/* LEFT AND RIGHT PREDICTORS */
/* OUR NOTATION FOR LEFT & RIGHT STATES: */
/* --UR_{i-1}-+-UL_i-----UR_i-+-UL_{i+1}-- */
/* ---+---------------+---------------+--- */
/* i-1 i i+1 */
/* Parameter adjustments */
--u;
/* Function Body */
i__1 = igrid_1.nx - 1;
for (i__ = 1; i__ <= i__1; ++i__) {
ulhalf[i__ - 1] = u[i__];
urhalf[i__ - 1] = u[i__ + 1];
}
/* AT RIGHT BOUNDARY, RIEMANN DATA */
ulhalf[igrid_1.nx - 1] = rdata_1.ur;
urhalf[igrid_1.nx - 1] = rdata_1.ur;
/* STATIONARY STATE AT i+1/2 */
i__1 = igrid_1.nx - 1;
for (i__ = 1; i__ <= i__1; ++i__) {
/* case 1 Trivial */
if (urhalf[i__ - 1] == ulhalf[i__ - 1]) {
us[i__ - 1] = urhalf[i__ - 1];
}
/* case 2 Shock */
if (ulhalf[i__ - 1] > urhalf[i__ - 1]) {
speed = (ulhalf[i__ - 1] + urhalf[i__ - 1]) * .5;
if (speed > 0.) {
us[i__ - 1] = ulhalf[i__ - 1];
}
if (speed < 0.) {
us[i__ - 1] = urhalf[i__ - 1];
}
if (speed == 0.) {
us[i__ - 1] = (float)0.;
}
}
/* case 3 Rarefaction */
/* suffers dog leg and entropy violation */
if (ulhalf[i__ - 1] < urhalf[i__ - 1]) {
tail = ulhalf[i__ - 1];
head = urhalf[i__ - 1];
if (tail > 0.) {
us[i__ - 1] = ulhalf[i__ - 1];
} else if (head < 0.) {
us[i__ - 1] = urhalf[i__ - 1];
} else {
us[i__ - 1] = (float)0.;
}
}
/* L200: */
}
/* AT THE BOUNDARY */
us[igrid_1.nx - 1] = rdata_1.ur;
/* CONSERVATIVE DIFFERENCE */
i__1 = igrid_1.nx;
for (i__ = 2; i__ <= i__1; ++i__) {
up[i__ - 1] = u[i__] - steps_1.dtdx * (us[i__ - 1] * (float).5 * us[
i__ - 1] - us[i__ - 2] * (float).5 * us[i__ - 2]);
/* L250: */
}
/* RIEMANN DATA */
up[0] = rdata_1.ul;
i__1 = igrid_1.nx;
for (i__ = 1; i__ <= i__1; ++i__) {
u[i__] = up[i__ - 1];
/* L300: */
}
return 0;
} /* step_ */