# Simple numerical model of tide with an elliptic scheme

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

# Now we can compute components associated with zonal and meridional speed
#    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)]
for i in range(grid_size)])
Y_pos=np.array([[dy[i,0]/2. + dy[i,:j].sum()
for j in range(grid_size)]
for i in range(grid_size)])

# 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.

• 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 Aug 26 '18 at 15:54
• Hello ! Thanks for the answer. I'm trying that right away ! – Koekoe94 Aug 28 '18 at 22:10
• 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 Aug 28 '18 at 23:04