PDEONE Example B: Shallow Water Flow
Supposed one wanted to solve the following system of partial
differential equations on the interval [-400,+400] (units are in
cm).
where u and
denote the horizontal velocity and depth of the fluid, H=H(x)
describes the surface over which the flow occurs, and g=980 is the
vertical acceleration due to gravity. These equations we
recognized as the shallow water equations. We desire the solution
of the PDE above with boundary conditions given by
With initial conditions and bottom topography given by
This corresponds to the second example in the original reference
to PDEONE (given above). This is quite simple to solve with the
pdeone gateway in octave. First the user must specify an octave
function "F.m" which evaluates the right hand side of the PDE
expressed in the following form
The incoming arguments to the octave function F.m are t,x,u,ux,
and duxx. Here t is the scalar time, x the scalar grid point, u a
vector of the NPDE unknowns, ux a vector of the x derivatives of
these unknowns, and duxx is a matrix expressing the
following derivatives
With the background we can write our F.m function for evaluating the right hand
side of our PDEs. In addition to F.m, PDEONE requires the
information on the diffusion coefficients D. This is provided in
the form of a function
D.m.
The routine "F.m" calls another octave function d_bottom_H.m to implement the derivative of
the bottom topography.
PDEONE requires the boundary conditions of the form
To implement this in octave we must specify three boundary
functions with input arguments t,x, and u (with the same meanings as
above):
BNDRY_ALPHA.m,
BNDRY_BETA.m, and
BNDRY_GAMMA.m.
Finally the user must specify the initial conditions. In this
example these are specified with the function shallow_init.m, where the bottom is described
in the octave function bottom_H.m. The two inputs to the initial
condition function are t and x. Here t is a scalar input of time
but in this case x is a vector of grid points where the initial
conditions are desired. Note, that to work properly this function
must be vectorized meaning that it can return a matrix
(of dimension [NPDE,Number of Grid Points]) of initial
conditions.
Once this code is written a person can call "pdeone.m". Here we
present an example pdeone_Script.m driver file. Note: Since
this problem is in Cartesian the constant c=0. Since this is the
default, there is no need to specify anything additional.
Finally, we can reap the benefits of our work by viewing the
solution to our problem:
Note the downstream shock that develops. This last plot is a
replica of the plot given in the Sincovec and Madsen paper and is
itself a duplicate of Figure 14b found in
Nonlinear shallow fluid flow over an isolated ridge
by D. D. Houghton and A. Kashara
Comm. Pure Appl. Math. Vol. 21 (1968) 1-23
Due to the ease of use of the octave interface one can quickly
recognize the stationary shock locations, refine the grid around
these locations and view the results.
John Weatherwax
Last modified: Sat Jan 7 10:58:30 EST 2006