Notebook

Wave Tank Inflow Boundary Conditions¶

In [26]:

%pylab inline

Populating the interactive namespace from numpy and matplotlib

The next cell defines a wave that is the sum of two plane waves, with normals at angles $\theta$ and $-\theta$ to the $x$ -axis. So if the two amplitudes A[0] and A[1] are equal, the $y$ -component of velocity $v$ is identically zero along side walls placed at any integer multiples of width.

In [27]:

import clawpack.visclaw.colormaps as C

theta = pi/8.  
h0 = 10.
grav = 9.81
c = sqrt(grav*h0)
width = 10.
L =  2*width*sin(theta) 
print "Wavelength in direction of wave propagation L = %g" % L
print "Linearized propagation speed (if A small relative to h0): c = %g" % c

def wave(x,y,t,A):
    eta0 = A[0]*sin(2*pi*(x*cos(theta) + y*sin(theta) - c*t)/L)
    eta1 = A[1]*sin(2*pi*(x*cos(theta) - y*sin(theta) - c*t)/L)
    h = h0 + eta0 + eta1
    u = sqrt(grav/h0) * (eta0 + eta1) * cos(theta)
    v = sqrt(grav/h0) * (eta0 - eta1) * sin(theta)
    return h,u,v

Wavelength in direction of wave propagation L = 7.65367
Linearized propagation speed (if A small relative to h0): c = 9.90454

In [28]:

x = linspace(-20,20,201)
y = linspace(-width,width,101)
X,Y = meshgrid(x,y)

def plot_wave(t,A):
    h,u,v = wave(X,Y,t,A)
    pcolor(X,Y,h,cmap=C.blue_white_red)
    contour(X,Y,v,colors='k')
    title('t = %8.3f, pcolor=h, contours=v' % t)
    axis('scaled')

Here are the two waves separately:

In [29]:

figure(figsize=(10,6))
subplot(121)
plot_wave(0,[1,0])
subplot(122)
plot_wave(0,[0,1])

And the linear combination with equal amplitudes, at two times:

In [30]:

plot_wave(0,[1,1])
figure()
plot_wave(0.5,[1,1])

To implement in-coming wave in the bc2 subroutine:¶

Let $h_1, \mu_1$ be the depth and momentum ( $\mu=hu$ ) in the first interior cell at time $t_n$ . The topography value $B_1$ is extrapolated to ghost cells, so $B_0 = B_1$ and the Riemann problem at this interface has no jump in topography.

Let $\eta(y,t)$ be the desired surface elevation at the left boundary (if there were no outgoing waves), i.e. eta0 + eta1 in the function above. What we use is the variation in elevation in this right-going 2-wave, which is independent of outgoing waves.

We want to specify $h_0, \mu_0$ in the ghost cell at each $y=y_j$ in such a way that the incoming wave has a jump in surface elevation equal to $\Delta \eta = \eta(y_j,t_{n+1}) - \eta(y_j,t_n)$ .

We don't care what the outgoing 1-wave is in the resulting Riemann solution since the ghost cell values are over-written every time step, so we can choose $h_0, \mu_0$ so that the Riemann solution consists solely of a 2-wave with $h_1 - h_0 = \Delta \eta$ . This implies we should choose $h_0 = h_1 - \Delta \eta$ .

Since the Roe solver is being used, the 2-wave is proportional to the eigenvector $[1, \bar c]^T$ , where $\bar c = \sqrt{g(h_0 + h_1)/2}$ , which can be computed since $h_0$ is known. Hence $\Delta\mu = \bar c \Delta h = \bar c \Delta\eta$ . This gives $\mu_0 = \mu_1 - \Delta\mu = \mu_1 - \bar c\Delta\eta$ .

The Fortran fragment implementing this is below. (It is in Clawpack 4.x form following the code in the repository, for 5.x the indices of q need to be re-ordered.)

This assumes the function eta(y,t) is provided.

    do j = 1-mbc, my+mbc
        yj = ylower + (j-0.5d0) * dy
        d_eta = eta(yj,time+dt) - eta(yj,time)
        h1 = q(1,j,1)
        h0 = h1 - d_eta
        cbar = sqrt(grav*0.5d0*(h0+h1))
        mu0 = q(1,j,2) - cbar * d_eta
        do ibc=1,mbc
            q(1-ibc,j,1) = h0
            q(1-ibc,j,2) = mu0
            enddo
        enddo

Please check this for errors since I haven't implemented it, but I've used similar things for other problems lately and the general idea should be right.

In [30]: