<p>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.</p>

<p>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 :</p>

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

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

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

<p>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$.</p>

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

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

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

<p>(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)) $</p>

<p>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. </p>

<p>Here is the code I use :</p>

<pre><code>#******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&gt;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')
</code></pre>

<p>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.</p>

<p>Thanks for your help, and sorry if it seems really basic.</p>

CC BY-SA 4.0
2018-08-26T15:54:10.617Z
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.
AtmosphericPrisonEscape
CC BY-SA 4.0
2018-08-28T22:10:06.817Z
Hello ! Thanks for the answer. I'm trying that right away !
Koekoe94
CC BY-SA 4.0
2018-08-28T23:04:56.757Z
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.
Koekoe94
<models><tides>
Simple numerical model of tide with an elliptic scheme
Koekoe94
# Simple numerical model of tide with an elliptic scheme

(Asked by: username_0 on 2018-08-24)

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.

**Tags:** models, tides

> **Comments:**
> **username_1**: 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.  
> **username_0**: Hello ! Thanks for the answer. I'm trying that right away !  
> **username_0**: 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.  

---

## 0 answers

