# Solving numerically a linearized system of elliptic (?) Navier-Stokes equation (Shallow Water Derived)

For my PhD Thesis, my advisor asked me to build a solver inspired from the article "Optimal Control Theory Applied to an Objective Analysis of a Tidal Current Mapping by HF Radar, J-L Devenon, 1989". In this article, in the numerical fluid-solving part, the author linearize the shallow water 2-D equations by removing the advection terms and using a linear friction coefficient. Then he uses a Fourier transform on the equation to finally end up with the system:

\begin{aligned} -i\sigma\vec{U}+\vec{f}\times\vec{U}+g\vec{\nabla}h&=-r\vec{U}\\\\ i\sigma h+\vec{\nabla}\cdot\left(H\vec{U}\right)&=0 \end{aligned}

Here $i$ is the basis for the imaginary number, $\sigma$ the frequency of the given Fourier component, $\vec{f}$ the Coriolis force vector, $g$ the gravitational acceleration, $h$ the local sea-level elevation compared to "mean local level", $r$ the linear friction coefficient, $H$ the local depth and $\vec{U}$ the local velocity.

Then, this system is simplified in a singled vector equation of the form :

$$F(\vec{V})=\lambda\vec{V}+\vec{f}\times\vec{V}+\vec{\nabla}\left(\vec{\nabla}\cdot\left(\gamma\vec{V}\right)\right)=\vec{0}$$

Where:

$$\lambda = rH-i\sigma,\quad \gamma=\frac{-igH}{\sigma}$$

And $\vec{V}$ is the transport, given as the product of the speed (supposed constant on the entire water-column) and the local height of the water-column( supposing the fluctuations $h$ are small compared to the overall depth).

He then solves the Dirichlet problems with specified boundary conditions. For this, he uses an iterative relaxation method (from Smith, G.D, Numerical solution of partial differential equations: Finite differences methods, 1978). And use only centered discretisation schemes.

I then wrote my own code for solving it. My first attempts where on a rectangular bassin with fixed boundary conditions. However It did not seem to work with any condition that differed from 0. I also tried to figure out analytical solutions but the only (very obvious) solution I found was the null solutions everywhere. Here is an example of the code I used:

############################################
# First test of implementing Devenon's     #
# Equation system and solving them in a clo#
# -sed basin                               #
############################################

#*****Packages import**********************#
import numpy as np

#****Functions*****************************#

def  relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,omega):
"""Function use to solve 2D linearized equation
of Navier-Stokes for a rigid lid. From Devenon 1989
. We're trying an overrelaxation technique with omega
as the relaxation parameter.

INPUT:
-----
U,V: Zonal and meridional speed, size nx*ny, also the output,complex valued
nx: Number of grid cell in x direction
ny: Number of grid cell in y direction
dx: Size of a grid cell in x direction
dy: Size of a grid cell in y direction
b: Forcing source function, of size nx x ny
gamma: Local complex phase speed
lamda : Local complex wave vector
omega: Relaxation factor, must belong to [1:2]

OUTPUT:
-------

# First we compute the non-homogeneous terms

#U_b=f*V[1:-1,1:-1]+ (gamma[:-2,:-2]*V[:-2,:-2] +
#       gamma[2:,2:]*V[2:,2:]-gamma[2:,:-2]*V[2:,:-2]-
#       gamma[:-2,2:]*V[:-2,2:])/(4*dy*dx)

#V_b= -f*U[1:-1,1:-1]+ (gamma[:-2,:-2]*U[:-2,:-2] +
#         gamma[2:,2:]*U[2:,2:]- gamma[2:,:-2]*U[2:,:-2]-
#         gamma[:-2,2:]*U[:-2,2:])/(4*dx*dy)

U_b=V[1:-1,1:-1]*0.
V_b=V[1:-1,1:-1]*0.
# Now we can use relaxation technique with an external term
U[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[2:,1:-1]*U[2:,1:-1]+
2*(1-omega)*gamma[1:-1,1:-1]*U[1:-1,1:-1]+(omega-1)*gamma[:-2,1:-1]*U[:-2,1:-1]) \
/dx**2  +1/lamda[1:-1,1:-1]*U_b

V[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[1:-1,2:]*V[1:-1,2:]+
2*(1-omega)*(gamma[1:-1,1:-1]*V[1:-1,1:-1])+(omega-1)*gamma[1:-1,:-2]*V[1:-1,:-2]) \
/dy**2 +1/lamda[1:-1,1:-1]*V_b

return U,V

#*****Parameters and variables*************#
nx = 40
ny = 40 # Number of grid points in each direction
nt  = 100 # Number of iteration of the relaxation method
xmin = 0.
xmax = 12000.
ymin = 0.
ymax = 12000. # Dimension (in meters) of the grid

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
# I understood that each speed would be complex-valued
# However, I wondered if U and V could only be the real and imaginary
# Part of the same variable in this case
U,V  = np.zeros((nx,ny),dtype=complex),np.zeros((nx,ny),dtype=complex)
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)
X, Y = np.meshgrid(x, y)

H=10 # Setting constant depth for testing
sigma=2.31484e-5 # Approximately semidiurnal frequency,  for tide (in hertz)
g=9.81 # Gravity acceleration
K=2.5e-2 #Friction coefficient
f=2e-4
lamda=np.ones((nx,ny),dtype=complex)*(K-sigma*1j)# (Friction and wavenumber)
gamma=np.ones((nx,ny),dtype=complex)*-1j*g*H/sigma #(Gravity and wavenumber)

#Initial conditions
U[0,:]=1.
U[-1,:]=1.
U[:,0]=1.
U[:,-1]=1.

V[0,:]=1.
V[-1,:]=1.
V[:,0]=1.
V[:,-1]=1.

#********Solving****************************#

for i in range(nt):
U,V=relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,2.)


For any conditions, I used except for the constant ones I found whether a very stable behaviour (the system stayed exactly in its initial position) or a very unstable one (I have instabilities growing from the sides that propagate, very quickly giving me aberrant values). This behaviour seems to be driven by the parameters lambda, gamma as well as the size of the cells. However no combination I have tried yet helped me finding a solution. This leads me to believe that I did a very bad job at coding this.

I tried searching the site ( and the web in general for answers) and found many links for nice EDP and even Navier-Stokes solvers. Though I did not get anything that I could link to my problem, may be due to a poor comprehension of it. This being said, I have another constraint that my advisor would like this to become an educational and simple problem. I, therefore, have to understand the internal functioning of my solver.

Then I have a few questions for you :

1) The first is obvious: Why doesn't it work? Did I do a bad job at coding? What should I correct to have a correct working algorithm?

2) The second is less obvious: Can I write my vector V in a complex manner, using one component of the transport as the real part and the second as the imaginary part? This seems obvious to my advisor but it does not look obvious to me at all.

3) The paper refers to it as an "elliptic" problem, close to a Helmholtz one. However, I just see two ODE's with a coupling term. Can I talk about an elliptic problem in this case?

• This would seem to be something you should have longer discussions with your adviser about! – Wolfgang Bangerth Apr 7 '18 at 4:02