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
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)
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.
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(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 i^{th }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 i^{th} 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)
cflv(1:4) : array of values related to Courant number:
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.
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:
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.
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
rp1 : user-supplied subroutine that implements the Riemann solver
The form of this subroutine is
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 i^{th} 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
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.
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)
info = 0
t = tstart
maxn = nv(1)
dt = dtv(1)
cflmax = 0.d0
dtmin = dt
dtmax = dt
nv(2) = 0
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
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
if (maxn.eq.0)
go to 900
do 100 n=1,maxn
told = t
if (told+dt .gt. tend)
dt = tend - toldif (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))
endif40 continue
dt2 = dt / 2.d0
thalf = t + dt2 midpoint in time for Strang splitting
t = told + dt time at end of step
extend data from grid to bordering
boundary cells:
call bc1(maxmx,meqn,mbc,mx,q,aux,told)
if (method(5).eq.1) then
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
endif
if (method(4) .eq. 1)
if (method(1) .eq. 1) then choose new time step if variable time step
endif
check to see if the Courant number was
too large:
if (cfl .le. cflv(1)) then
else reject this step
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
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)
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
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 i^{th} 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
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
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
endifdtdx(i) = dt / (dx*aux(i,mcapa))
else
dtdx(i) = dt/dx
endif
10 continue
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