 subroutine copyq1(maxmx,meqn,mbc,mx,q1,q2)

copy the contents of q1 into q2 c

implicit double precision (a-h,o-z)
dimension q1(1-mbc:maxmx+mbc, meqn)
dimension q2(1-mbc:maxmx+mbc, meqn)

do 10 i = 1-mbc, mx+mbc

do 10 m=1,meqn
q2(i,m) = q1(i,m)

10 continue

return

end subroutine claw1(maxmx,meqn,mwaves,mbc,mx, q,aux,dx,tstart,tend,dtv,cflv,nv,method,mthlim, work,mwork,info,bc1,rp1,src1)

Introduction

Solves a hyperbolic system of conservation laws in one space dimension of the general form
capa * q_t + A q_x = psi
The "capacity function" capa(x) and source term psi are optional (see below). For a more complete description see the documentation in the directory claw/doc, especially note1.ps and note2.ps. Sample driver programs and user-supplied subroutines are available  see the the directory claw/clawpack/1d/example for one example, and codes in claw/applications for more extensive examples. The user must supply the following subroutines: cbc1, rp1 subroutines specifying the boundary conditions and Riemann solver. These are described in greater detail below. In addition, if the equation contains source terms psi, then the user  must provide: src1 subroutine that solves capa * q_t = psi   over a single time step. These routines must be declared EXTERNAL in the main program. For description of the calling sequences, see below.

The user must supply the following subroutines: bc1, rp1 subroutines specifying the boundary conditions and Riemann solver. These are described in greater detail below. In addition, if the equation contains source terms psi, then the user must provide: src1 subroutine that solves capa * q_t = psi over a single time step. These routines must be declared EXTERNAL in the main program. For description of the calling sequences, see below.

Description of parameters

maxmx : maximum number of interior grid points in x, and is used in declarations of the array q.

meqn    : number of equations in the system of conservation laws.

mwaves : number of waves that result from the solution of each Riemann problem. Often mwaves = meqn but for some problems these may be different.

mbc      :  number of "ghost cells" that must be added on to each side of the domain to handle boundary conditions. The cells actually in the physical domain are labelled from 1 to mx in x. The arrays are dimensioned actually indexed from 1-mbc to mx+mbc. For the methods currently implemented, mbc = 2 should be used. If the user implements another method that has a larger stencil and hence requires more ghost cells, a larger value of mbc could be used. q is extended from the physical domain to the ghost cells by the user-supplied routine bc1.

mx     :   number of grid cells in the x-direction, in the physical domain. In addition there are mbc grid cells along each edge of the grid that are used for boundary conditions. Must have mx .le. maxmx q(1-mbc:maxmx+mbc, meqn)
On input: initial data at time tstart.
On output: final solution at time tend. q(i,m) = value of mth component in the i'th cell.
Values within the physical domain are in q(i,m) for i = 1,2,...,mx mbc extra cells on each end are needed for boundary conditions as specified in the routine bc1.

method(1:7)   :   array of values specifying the numerical method to use

method(1) = 0 if fixed size time steps are to be taken. In this case, dt = dtv(1) in all steps.
= 1 if variable time steps are to be used. In this case, dt = dtv(1) in the first step and thereafter the value cflv(2) is used to choose the next time step based on the maximum wave speed seen in the previous step. Note that since this value comes from the previous step, the Courant number will not in general be exactly equal to the desired value If the actual Courant number in the next step is greater than 1, then this step is redone with a smaller dt.

method(2) = 1 if Godunov's method is to be used, with no 2nd order corrections.
= 2 if second order correction terms are to be added, with a flux limiter as specified by mthlim.
= 3 if "third order" correction terms are to be added, based on my paper my paper "A high-resolution conservative algorithm for advection in incompressible flow". This is currently recommended only for problems with smooth solutions, using no limiter (mthlim = 0)

method(3) is not used in one-dimension.

method(4) = 0 to suppress printing
= 1 to print dt and Courant number every time step

method(5) = 0 if there is no source term psi. In this case the subroutine src2 is never called so a dummy parameter can be given.
= 1 if there is a source term. In this case the subroutine src2 must be provided.

method(6) = 0 if there is no capacity function capa.
= mcapa > 0 if there is a capacity function. In this case aux(i,mcapa) is the capacity of the ith cell and you must also specify method(7) .ge. mcapa and set aux.

method(7) = 0 if there is no aux array used.
= maux > 0 if there are maux auxiliary variables.

The recommended choice of methods for most problems is method(1) = 1, method(2) = 2.

aux(1-mbc:maxmx+mbc, maux)   :   array of auxiliary variables that are used in specifying the problem.
If method(7) = 0 then there are no auxiliary variables and aux can be a dummy variable.
If method(7) = maux > 0 then there are maux auxiliary variables and aux must be dimensioned as above.
Capacity functions are one particular form of auxiliary variable. These arise in some applications, e.g. variable coefficients in advection or acoustics problems. See Clawpack Note # 5 for examples.
If method(6) = 0 then there is no capacity function.
If method(6) = mcapa > 0 then there is a capacity function and capa(i), the "capacity" of the ith cell, is assumed to be stored in aux(i,mcapa). In this case we require method(7).ge.mcapa.

dx : grid spacing in x. (for a computation in ax <= x <= bx, set dx = (bx-ax)/mx.)

tstart : initial time.

tend desired final time (on input) ; actual time reached (on output).

dtv(1:5)   :   array of values related to the time step: (Note: method(1)=1 indicates variable size time steps)

dtv(1) = value of dt to be used in all steps if method(1) = 0
value of dt to use in first step if method(1) = 1
dtv(2) = unused if method(1) = 0.
= maximum dt allowed if method(1) = 1.
dtv(3) = smallest dt used (on output)
dtv(4) = largest dt used (on output)
dtv(5) = dt used in last step (on output)

cflv(1:4)  : array of values related to Courant number:

cflv(1) = maximum Courant number to be allowed. With variable time steps the step is repeated if the Courant number is larger than this value. With fixed time steps the routine aborts. Usually cflv(1)=1.0 should work.
cflv(2) = unused if method(1) = 0.
= desired Courant number if method(1) = 1. Should be somewhat less than cflv(1), e.g. 0.9
cflv(3) = largest Courant number observed (on output).
cflv(4) = Courant number in last step (on output).

nv(1:2)  :   array of values related to the number of time steps:

nv(1) = unused if method(1) = 0
= maximum number of time steps allowed if method(1) = 1
nv(2) = number of time steps taken (on output).

mthlim(1:mwaves)   :   array of values specifying the flux limiter to be used in each wave family mw. Often the same value will be used for each value of mw, but in some cases it may be desirable to use different limiters. For example, for the Euler equations the superbee limiter might be used for the contact discontinuity (mw=2) while another limiter is used for the nonlinear waves. Several limiters are built in and others can be added by modifying the subroutine philim.

mthlim(mw) = 0 for no limiter
= 1 for minmod
= 2 for superbee
= 3 for van Leer
= 4 for monotonized centered

work(mwork)   :   double precision work array of length at least mwork

mwork      :     length of work array. Must be at least (maxmx + 2*mbc) * (2 + 3*meqn + mwaves + meqn*mwaves) If mwork is too small then the program returns with info = 4 and prints the necessary value of mwork to unit 6.

info   :   output value yielding error information:

= 0 if normal return.
= 1 if mx.gt.maxmx or mbc.lt.2
= 2 if method(1)=0 and dt doesn't divide (tend - tstart).
= 3 if method(1)=1 and cflv(2) > cflv(1).
= 4 if mwork is too small.
= 11 if the code attempted to take too many time steps, n > nv(1). This could only happen if method(1) = 1 (variable time steps).
= 12 if the method(1)=0 and the Courant number is greater than 1 in some time step.

Note: if info.ne.0, then tend is reset to the value of t actually reached and q contains the value of the solution at this time.

User-supplied subroutines

bc1   :   subroutine that specifies the boundary conditions. This subroutine should extend the values of q from cells 1:mx to the mbc ghost cells along each edge of the domain.

The form of this subroutine is

subroutine bc1(maxmx,meqn,mbc,mx,q,aux,t)
implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)

rp1    :   user-supplied subroutine that implements the Riemann solver

The form of this subroutine is

subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
implicit double precision (a-h,o-z)
dimension ql(1-mbc:maxmx+mbc, meqn)
dimension qr(1-mbc:maxmx+mbc, meqn)
dimension auxl(1-mbc:maxmx+mbc, *)
dimension auxr(1-mbc:maxmx+mbc, *)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)

On input, ql contains the state vector at the left edge of each cell , qr contains the state vector at the right edge of each cell , auxl contains auxiliary values at the left edge of each cell and auxr contains auxiliary values at the right edge of each cell.

Note that the ith Riemann problem has left state qr(i-1,:) and right state ql(i,:) In the standard clawpack routines, this Riemann solver is called with ql=qr=q along this slice. More flexibility is allowed in case the user wishes to implement another solution method that requires left and rate states at each interface. c If method(7)=maux > 0 then the auxiliary variables along this slice are passed in using auxl and auxr. Again, in the standard routines auxl=auxr=aux in the call to rp1.

On output, wave(i,m,mw) is the m'th component of the jump across wave number mw in the ith Riemann problem. s(i,mw) is the wave speed of wave number mw in the ith Riemann problem. amdq(i,m) = m'th component of A^- Delta q, apdq(i,m) = m'th component of A^+ Delta q, the decomposition of the flux difference f(qr(i-1)) - f(ql(i)) into leftgoing and rightgoing parts respectively.

It is assumed that each wave consists of a jump discontinuity propagating at a single speed, as results, for example, from a Roe approximate Riemann solver. An entropy fix can be included into the specification of amdq and apdq.

src1   :   subroutine for the source terms that solves the equation capa * q_t = psi over time dt.

If method(5)=0 then the equation does not contain a source term and this routine is never called. A dummy argument can be used with many compilers, or provide a dummy subroutine that does nothing (such a subroutine can be found in clawpack/1d/misc/src1xx.f)

The form of this subroutine is

subroutine src1(maxmx,meqn,mbc,mx,q,aux,t,dt)
implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)

If method(7)=0 or the auxiliary variables are not needed in this solver, then the latter dimension statement can be omitted, but aux should still appear in the argument list.

On input, q(i,m) contains the data for solving the source term equation. On output, q(i,m) should have been replaced by the solution to the source term equation after a step of length dt.

Beginning of claw1 code

implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)
dimension work(mwork)
dimension mthlim(mwaves),method(7),dtv(5),cflv(4),nv(2)

• Initialisation

info = 0
t = tstart
maxn = nv(1)
dt = dtv(1)
cflmax = 0.d0
dtmin = dt
dtmax = dt
nv(2) = 0

• Check for errors in data

if (mx .gt. maxmx) then

info = 1
go to 900

endif

if (method(1) .eq. 0) then fixed size time steps. Compute the number of steps:

maxn = (tend - tstart + 1d-10) / dt

if (dabs(maxn*dt - (tend-tstart)) .gt. 1d-8) then
dt doesn't divide time interval integer number of times
info = 2
go to 900
endif

endif

if (method(1).eq.1 .and. cflv(1).gt.1.d0) then

info = 3
go to 900

endif

• Partition work array into pieces for passing into step1:

i0f = 1
i0wave = i0f + (maxmx + 2*mbc) * meqn
i0s = i0wave + (maxmx + 2*mbc) * meqn * mwaves
i0dtdx = i0s + (maxmx + 2*mbc) * mwaves
i0qwork = i0dtdx + (maxmx + 2*mbc)
i0amdq = i0qwork + (maxmx + 2*mbc) * meqn
i0apdq = i0amdq + (maxmx + 2*mbc) * meqn
i0dtdx = i0apdq + (maxmx + 2*mbc) * meqn
i0end = i0dtdx + (maxmx + 2*mbc) - 1

if (mwork .lt. i0end) then

write(6,*) 'mwork must be increased to ',i0end
info = 4
go to 900

endif

• Main loop

if (maxn.eq.0)

go to 900

do 100 n=1,maxn

told = t

if (told+dt .gt. tend)
dt = tend - told

if (method(1).eq.1) then
save old q in case we need to retake step with smaller dt:
call copyq1(maxmx,meqn,mbc,mx,q,work(i0qwork))
endif

40 continue
dt2 = dt / 2.d0
thalf = t + dt2 midpoint in time for Strang splitting
t = told + dt time at end of step

• main steps in algorithm, using Strang splitting for source term:

extend data from grid to bordering boundary cells:
call bc1(maxmx,meqn,mbc,mx,q,aux,told)
if (method(5).eq.1) then

with source term: use Strang splitting
call src1(maxmx,meqn,mbc,mx,q,aux,told,dt2)

endif

take a step on the homogeneous conservation law:
call step1(maxmx,meqn,mwaves,mbc,mx,q,aux,dx,dt, & method,mthlim,cfl,work(i0f),work(i0wave), & work(i0s),work(i0amdq),work(i0apdq),work(i0dtdx), & rp1)
if (method(5).eq.1) then

source terms over second half time step:
call src1(maxmx,meqn,mbc,mx,q,aux,thalf,dt2)

endif

if (method(4) .eq. 1)

write(6,601) n,cfl,dt,t
601 format('CLAW1... Step',i4, & ' Courant number =',f6.3,' dt =',d12.4, & ' t =',d12.4)

if (method(1) .eq. 1) then choose new time step if variable time step

if (cfl .gt. 0.d0) then
dt = dmin1(dtv(2), dt * cflv(2)/cfl)
dtmin = dmin1(dt,dtmin)
dtmax = dmax1(dt,dtmax)
else
dt = dtv(2)
endif

endif

check to see if the Courant number was too large:
if (cfl .le. cflv(1)) then

accept this step
cflmax = dmax1(cfl,cflmax)

else reject this step

t = told
call copyq1(maxmx,meqn,mbc,mx,work(i0qwork),q)

if (method(4) .eq. 1) then
write(6,602)
602 format('CLAW1 rejecting step... ', & 'Courant number too large')
endif

if (method(1).eq.1) then if variable dt, go back and take a smaller step
go to 40
else if fixed dt, give up and return
cflmax = dmax1(cfl,cflmax)
go to 900
endif

endif

see if we are done:
nv(2) = nv(2) + 1
if (t .ge. tend)
go to 900

100 continue
900 continue

• return information

if (method(1).eq.1 .and. t.lt.tend .and. nv(2) .eq. maxn) then
too many timesteps
info = 11
endif

if (method(1).eq.0 .and. cflmax .gt. cflv(1)) then
Courant number too large with fixed dt
info = 12
endif

tend = t
cflv(3) = cflmax
cflv(4) = cfl
dtv(3) = dtmin
dtv(4) = dtmax
dtv(5) = dt

return
end subroutine step1(maxmx,meqn,mwaves,mbc,mx,q,aux,dx,dt, method,mthlim,cfl,f,wave,s,amdq,apdq,dtdx,rp1)

Introduction

Take one time step, updating q.

method(1) = 1 ==> Godunov method
method(1) = 2 ==> Slope limiter
method mthlim(p) controls what limiter is used in the pth family

Parameters

amdq(1-mbc:maxmx+mbc, meqn)   :   left-going flux-differences
apdq(1-mbc:maxmx+mbc, meqn)    :   right-going flux-differences e.g.
amdq(i,m)     :  m'th component of A^- \Delta q from ith Riemann problem (between cells i-1 and i).

wave(1-mbc:maxmx+mbc, meqn, mwaves)     :    waves from solution of Riemann problems
wave(i,m,mw)      :     mth component of jump in q across wave in family mw in Riemann problem between states i-1 and i.

s(1-mbc:maxmx+mbc, mwaves) :  wave speeds, s(i,mw)
:  speed of wave in family mw in Riemann problem between states i-1 and i.

f(1-mbc:maxmx+mbc, meqn)      :   correction fluxes for second order method f(i,m)
:   mth component of flux at left edge of ith cell

Beginning of the code

implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)
dimension f(1-mbc:maxmx+mbc, meqn)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension dtdx(1-mbc:maxmx+mbc)
dimension method(7),mthlim(mwaves)
logical limit

• check if any limiters are used:

limit = .false.

do 5 mw=1,mwaves

if (mthlim(mw) .gt. 0) limit = .true.

5 continue

mcapa = method(6)

do 10 i=1-mbc,mx+mbc

if (mcapa.gt.0) then

if (aux(i,mcapa) .le. 0.d0) then
write(6,*) 'Error -- capa must be positive'
stop
endif

dtdx(i) = dt / (dx*aux(i,mcapa))

else

dtdx(i) = dt/dx

endif

10 continue

• Solve Riemann problem at each interface

call rp1(maxmx,meqn,mwaves,mbc,mx,q,q,aux,aux,wave,s,amdq,apdq)

Modify q for Godunov update:
Note this may not correspond to a conservative flux-differencing for equations not in conservation form. It is conservative if amdq + apdq = f(q(i)) - f(q(i-1)).

do 40 i=1,mx+1

do 40 m=1,meqn
q(i,m) = q(i,m) - dtdx(i)*apdq(i,m)
q(i-1,m) = q(i-1,m) - dtdx(i-1)*amdq(i,m)

40 continue

Compute maximum wave speed:
cfl = 0.d0
do 50 i=1,mx+1
cfl = dmax1(cfl, dtdx(i)*dabs(s(i,1)), dtdx(i)*dabs(s(i,mwaves)))
50 continue

if (method(2) .eq. 1)

go to 900

Compute correction fluxes for second order q_{xx} terms:
do 100 m = 1, meqn
do 100 i = 1-mbc, mx+mbc
f(i,m) = 0.d0
100 continue

Apply limiter to waves: if (limit) call limiter(maxmx,meqn,mwaves,mbc,mx,wave,s,mthlim)
do 120 i=1,mx+1
do 120 m=1,meqn
do 110 mw=1,mwaves
dtdxave = 0.5d0 * (dtdx(i-1) + dtdx(i))
f(i,m) = f(i,m) + 0.5d0 * dabs(s(i,mw)) * (1.d0 - dabs(s(i,mw))*dtdxave) * wave(i,m,mw)

Third order corrections:
(still experimental... works well for smooth solutions with no limiters but not well with limiters so far)

if (method(2).lt.3)

go to 110
if (s(i,mw) .gt. 0.d0) then
dq2 = wave(i,m,mw) - wave(i-1,m,mw)
else
dq2 = wave(i+1,m,mw) - wave(i,m,mw)
endif
f(i,m) = f(i,m) - s(i,mw)/6.d0 * (1.d0 - (s(i,mw)*dtdxave)**2) * dq2

110 continue
120 continue
140 continue

Update q by differencing correction fluxes
(Note: Godunov update has already been performed above)

do 150 m=1,meqn
do 150 i=1,mx
q(i,m) = q(i,m) - dtdx(i) * (f(i+1,m) - f(i,m))
150 continue
900 continue

return
end