This demo example shows ...



In [25]:

    
%pylab inline









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.



In [26]:

    
import numpy as np
import time
import matplotlib.pyplot as plt
from IPython.display import clear_output
from mpl_toolkits.mplot3d.axes3d import Axes3D

First, define some constants we need



In [27]:

    
# Number of grid cells in each direction of the grid
m = 21
n = 21

# Grid spacing
dx = 10/float(m)
dy = 10/float(n)

# The grid itself
Q = np.zeros((3, n, m))
Q[0, :, :] = 1.0
Q[0, n/2, m/2] = 3.0

# Gravitational constant
g = 9.80665

Vi ønsker å løse grundtvannslikningen eksplisitt

Grundtvannslikningene kan skrives

\begin{aligned} \left[ \begin{array}{c} h\\\ hu\\\ hv \end{array} \right]_t + \left[ \begin{array}{c} hu\\\ hu^2 + \tfrac{1}{2}gh^2\\\ huv \end{array} \right]_x + \left[ \begin{array}{c} hv\\\ huv\\\ hv^2 + \tfrac{1}{2}gh^2 \end{array} \right]_y = \left[ \begin{array}{c} 0\\\ 0\\\ 0 \end{array} \right]. \end{aligned}

Her er $h$ vannhøyden, $hu$ fluxen langs x-aksen, og $hv$ fluxen langs y-aksen. Vi kan skrive det mer kompakt på vektor form som

$Q_t + F(Q)_x + G(Q)_y = 0$

Og diskretisere ved å beregne fluxen inn og ut av hver celle:

$$ Q_t = - \frac{ F\_{i+1/2}(Q) - F\_{i-1/2}(Q) }{ \Delta x } - \frac{ G\_{j+1/2}(Q) - G\_{j-1/2}(Q) }{ \Delta y } $$

Vi kan approksimere $F(Q)$ ved integrasjonspunktet mellom to celler med hjelp av central-upwind fluks funksjonen

\begin{aligned} F_{j+1/2}(Q^k) = \frac{ a^+\_{j+1/2} F(Q^{k-}\_{j+1/2}) - a^-\_{j+1/2} F(Q^{k+}\_{j+1/2}) }{ a^+\_{j+1/2} - a^-\_{j+1/2} } + \frac{ a^+\_{j+1/2} a^-\_{j+1/2} }{ a^+\_{j+1/2} - a^-\_{j+1/2} } \left( Q^{k+}\_{j+1/2}-Q^{k-}\_{j+1/2} \right) \end{aligned}

Her er + brukt for å vise at verdien tas fra høyre side av interfacet, og - for å vise at den tas fra venstre side av interfacet, og $a$ er største/minste egenverdi for integrasjonspunktet (tilsvarer største bølgehastighet)



In [28]:

    
#Central upwind flux function
def CU(Fl, Ql , Fr, Qr, ap, am):
    return (ap*Fr - am*Fl) / (ap-am) + (ap*am)/(ap-am)*(Qr-Ql)

#Find the discrete F-flux at an integration point
def F(Ql, Qr):
    #Find the maximum / minimum local wave speeds
    cl = sqrt(g*Ql[0])
    cr = sqrt(g*Qr[0])
    ul = Ql[1]/Ql[0]
    ur = Qr[1]/Qr[0]
    am = min(min(ul-cl, ur-cr), 0.0)
    ap = max(max(ul+cl, ur+cr), 0.0)
    
    #Define the local helper function Fh
    def Fh(Q):
        return np.array([
            Q[1],
            Q[1]*Q[1]/Q[0] + 0.5*g*Q[0]*Q[0],
            Q[1]*Q[2]/Q[0]
            ])
    
    Fl = Fh(Ql)
    Fr = Fh(Qr)
    
    return CU(Fl, Ql, Fr, Qr, ap, am)

#Find the discrete G-flux at an integration point
def G(Ql, Qr):
    #Find the maximum / minimum local wave speeds
    cl = sqrt(g*Ql[0])
    cr = sqrt(g*Qr[0])
    vl = Ql[2]/Ql[0]
    vr = Qr[2]/Qr[0]
    am = min(min(vl-cl, vr-cr), 0.0)
    ap = max(max(vl+cl, vr+cr), 0.0)
    
    #Define the local helper function Gh
    def Gh(Q):
        return np.array([
            Q[2],
            Q[1]*Q[2]/Q[0],
            Q[2]*Q[2]/Q[0] + 0.5*g*Q[0]*Q[0]
            ])
    
    Gl = Gh(Ql)
    Gr = Gh(Qr)
    
    return CU(Gl, Ql, Gr, Qr, ap, am)

Vi beregner fluxen inn og ut av hver celle som følger

$$ Q^{k+1} = Q^k + \Delta t \left(\frac{F_{j-1/2}(Q^k) - F_{j+1/2}(Q^k)}{\Delta x} + \frac{G_{j-1/2}(Q^k) - G_{j+1/2}(Q^k)}{\Delta y} \right) = Q^k + \Delta t \cdot R(Q^k) $$



In [29]:

    
fig, ax = plt.subplots()

dt = 0.05


#Perform timestepping
for i in range(0, 100):
    f_flux = np.zeros((3, n, m-1))
    g_flux = np.zeros((3, n-1, m))
    
    # Loop over all intefaces and create the flux
    for k in range(n):
        for l in range(m-1):
            f_flux[:, k, l] = F(Q[:, k, l], Q[:, k, l+1])
    for k in range(n-1):
        for l in range(m):
            g_flux[:, k, l] = G(Q[:, k, l], Q[:, k+1, l])
    
    # Loop over all internal cells and sum fluxes
    for k in range(1, n-1):
        for l in range(1, m-1):
            Q[:, k, l] = Q[:, k, l] + dt * ( (f_flux[:, k, l-1] - f_flux[:, k, l])/dx + (g_flux[:, k-1, l] - g_flux[:, k, l])/dy )
            
    
    #Boundary conditions (wall)
    Q[:, :, 0] = Q[:, :, 1]
    Q[:, :, -1] = Q[:, :, -2]
    Q[1, :, 0] = -Q[1, :, 0]
    Q[1, :, -1] = -Q[1, :, -1]
    
    Q[:, 0, :] = Q[:, 1, :]
    Q[:, -1, :] = Q[:, -2, :]
    Q[2, 0, :] = -Q[2, 0, :]
    Q[2, -1, :] = -Q[2, -1, :]
    
    if (i % 10 == 0):
        #clear_output()
        ax.cla()
        x = np.arange(0, m)
        y = np.arange(0, n)
        x, y = np.meshgrid(x, y)
        ax = fig.gca(projection='3d')
        ax.plot_surface(x, y, Q[0, :, :], rstride=1, cstride=1, antialiased=True, cmap=cm.coolwarm);
        #time.sleep(0.2)
        display(fig)
plt.close()

$\Delta t$ er begrenset med følgende CFL-betingelse

$$ \Delta t \lt \frac{1}{4} \text{min}\left\\{\Delta x / |a_x|, \Delta y / |a_y| \right\\} $$

hvor $a_{x}$ og $a_y$ er største / minste egenverdi i domenet for $F$ og $G$ respektivt.

Vi trenger i tillegg randkrav. I vårt tilfelle benytter vi vegg randkrav, hvor vi benytter skygge-celler:

$Q_{-1,-1} = Q_{0,0}, \quad Q_{m+1, n+1} = Q_{m, n}$



In [29]: