3
$\begingroup$

I have recently been given the task of developping a coastal model for tides. I have a few requirements, it must not be time stepping in order to be as fast as possible and in finite differences because something as simple as possible is asked.

For a first start, I tried implementing linearized shallow water equations and supposed a time dependancy in the form of a complex exponential. This gave me the following system :

(1) $-i \omega \hat u -f \hat v +g \partial_x \hat h= -r \hat u$

(2) $-i \omega \hat v +f \hat u +g \partial_y\hat h= -r \hat v$

(3) $-i \omega \hat h + \partial_x H\hat u + \partial_y H\hat v = 0$

where my unknows are $\hat u$,$\hat v$,$\hat h$, and only two parameters possibly showing spatial variation : $r$ and $H$. I tried to discretize it with finite differences but I quickly realize that I needed to be very careful with my boundary conditions. Because the 3 variables that I wanted to determined where closely linked to one another and that actually $\hat u$ and $\hat v$ could easily be determined by the derivatives of $\hat h$.

So I modified my system in order to solve only this variable and came with the following set of equations :

(4) $ \hat u = \frac{-g(r-i \omega)\partial_x \hat h -fg\partial_y \hat h}{f²+(r-i\omega)²}$

(5) $ \hat v = \frac{-g(r-i \omega)\partial_y \hat h +fg\partial_x \hat h}{f²+(r-i\omega)²}$

(6) $ i \omega \hat h = \partial_x(\frac{H}{f²+(r-i\omega)²}(-g(r-i \omega)\partial_x \hat h -fg\partial_y \hat h)) + \partial_y(\frac{H}{f²+(r-i\omega)²}(-g(r-i \omega)\partial_y \hat h +fg\partial_x \hat h)) $

Then I tried to approximate the derivates with centered first order scheme. I tried to reproduce a Kelvin coastal wave on a reduced domain (10*10 cells). However, I am at best able not to crash it and most of the time have very strong oscillations in the center.

Here is the code I use :

#******Packages import*#
import numpy as np
import matplotlib.pyplot as plt

#******solving function**#

def solver_water_elevation(eta,depth,dx,dy,r,f=2e-4,sigma=2e-5,g=9.81):

    """Solver of linearized elliptic shallow water equation for complex coastal
    problems. We use a jacobi iterator in a first place. Then we will refine it.
    Variables here are mostly complex-valued.
    INPUTS:
    eta : Overelevation (m)
    depth : depth (m)
    dx,dy : Size of grid cell
    r : Linear friction coefficient (m.s-1)
    f : Coriolis frequency (s-1)
    sigma : Frequency of the considered wave (s-1)
    g : Gravitationnal acceleration (m.s-2)
    OUTPUTS :
    eta : See above"""

    # First we determine the local friction coefficient in m-1
    friction=r/depth
    #Then we will proceed step by step
    # First we compute gradients of eta
    etadx=(eta[1:-1,2]-eta[1:-1,:-2])/(2.*dx[1:-1,1:-1])
    etady=(eta[2:,1:-1]-eta[:-2,1:-1])/(2.*dy[1:-1,1:-1])

    # Now we can compute components associated with zonal and meridional speed
    u=-g*(friction[1:-1,1:-1]-1j*sigma)*etadx-f*g*etady
    v=-g*(friction[1:-1,1:-1]-1j*sigma)*etady+f*g*etadx
#    print('v',v)

    # Now we multiply them so that they correspond to the right terms
    multiplicator=depth[1:-1,1:-1]/(f**2+(friction[1:-1,1:-1]-1j*sigma)**2)
    print('multiplicator',multiplicator)
    #Now we compute the gradient of those terms
    dxu=(multiplicator[1:-1,2:]*u[1:-1,2:]-multiplicator[1:-1,:-2]*u[1:-1,:-2])\
        /(2.*dx[2:-2,2:-2])
    dyv=(multiplicator[2:,1:-1]*v[2:,1:-1]-multiplicator[:-2,1:-1]*v[:-2,1:-1])\
        /(2.*dy[2:-2,2:-2])
    # And finally, determine eta
    eta[2:-2,2:-2]=-1j/sigma*(dxu+dyv)

    return eta

def simple_convergence_error(eta,eta_old):
    """Compare the value of 2 array of eta and returns the maximum              relative error.

    INPUTS :
    eta : New array 
    eta_old :old array
    OUTPUTS :
    error : Maximal relative difference between the two arrays"""

    return np.max(abs((eta-eta_old)))

def sol_Kelvin_wave(depth,grid_size,dx,dy,f=2e-4,
                    sigma=2.31484e-5,g=9.81,h_0=1.+0.5j) :
    """ Here we compute the analytical solution of Kelvin waves
    for the conditions given in the BC.

    INPUTS:
    depth : Static water depth, in meters
    grid_size : size of the grid (in number of cells)
    dx,dy : Size of a grid cell in x and y directions, in meters
    f : Coriolis frequency, in s-1
    sigma : Frequency of the given tide in s-1
    g : Gravity acceleration, in m.s-2
    h_0 : Initial overelevation, in m (Complex valued)
    OUTPUTs:
    eta : Overelevation of water"""

    # First we set some useful values
    c=np.sqrt(g*depth)
    omega=sigma*2.*np.pi # We get in rad.s-1
    k=omega/c # The wavevector

    # We must determine the distance to boundaries
    X_pos=np.array([[dx[0,j]/2. + dx[:i,j].sum() 
        for j in range(grid_size[1])] 
        for i in range(grid_size[0])])
    Y_pos=np.array([[dy[i,0]/2. + dy[i,:j].sum()
        for j in range(grid_size[1])] 
        for i in range(grid_size[0])])

    # Computing the value
    eta=h_0*np.exp(1.j*k[:,:]*X_pos[:,:])\
        *np.exp(-f*Y_pos[:,:]/c[:,:])

    return eta

#******Defining variables*****#
nx,ny=10,10
eta=np.zeros((nx,ny),dtype=complex)
dx,dy=10.*np.ones((nx,ny)),10.*np.ones((nx,ny))
sigma=2e-4 # frequency
f=2e-4 # Coriolis frequency
r=0.#2.5e-3 # Friction coefficient
depth=np.ones((nx,ny)) # Constant depth

#***I try to initialize it with a analytic solution***#
eta_test=sol_Kelvin_wave(depth,(nx,ny),dx,dy,f,sigma)
# And now I set the boundary conditions
eta[:2,:]=eta_test[:2,:] # Bottom
eta[-2:,:]=eta_test[-2:,:] # Upper side
eta[:,:2]=eta_test[:,:2] # Left side
eta[:,-2:]=eta_test[:,-2:] #Right side
#*****Calling the solver**********#

error=1. # Initialzing the error
tol=1e-8 # Tolerance to error
nbr_iter=0

 while error>tol : # We use a simple iterative solver
    # Counting the iterations
    nbr_iter+=1
    eta_old=eta.copy()# Saving the old value
    #Calling the solver
    eta=solver_water_elevation(eta,depth,dx,dy,r,sigma=sigma*2.*np.pi,f=f)
    #Checking the error
    error=simple_convergence_error(eta,eta_old)
    print('Computed in ',nbr_iter,' iterations')

Is there a clear error in my code ? Or am I just trying to solve my problem with unadapted tools ? It seems very basic but still it is clearly running out of my hands. And I found only little litterature about elliptic equations for tide computing. So I guessed it was not very efficient but I was totally unable to find a good explanation of why it was so inefficient.

Thanks for your help, and sorry if it seems really basic.

$\endgroup$
  • $\begingroup$ One thing to look for is your approximation of the derivatives: In hyperbolic problems the symmetric first order derivative is numerically unstable. Oscillations occur at discontinuities and are amplified. This might be also the case for your elliptic scheme. Try simple upwinding. $\endgroup$ – AtmosphericPrisonEscape Aug 26 '18 at 15:54
  • $\begingroup$ Hello ! Thanks for the answer. I'm trying that right away ! $\endgroup$ – Koekoe94 Aug 28 '18 at 22:10
  • $\begingroup$ So I tried to do a not centered approximation, but it does not seem to work neither. I guess I did a mistake in my equation and I am trying to solve a not physical system. $\endgroup$ – Koekoe94 Aug 28 '18 at 23:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.