Numerical solution of the 2D linear partial differential equation of first order

Question

It is an actual geophysical problem where we study the liquid flow. We measure at each 2D grid point and each time interval three components of the liquid velocity and we want to compute how the liquid level changes over time.

If I understand the problem correctly (which is beyond this question, unless you specifically know the answer), we can assume that the problem is time-independent, that is we find the liquid level as a function of the spatial coordinates $x,y$ at a specific time $t$.

This problem can be rewritten as an equation for each time step as:

$$ (V_{x}\frac{\partial h}{\partial x}+h\frac{\partial V_{x}}{\partial x}) + (V_{y}\frac{\partial h}{\partial y}+h\frac{\partial V_{y}}{\partial y}) = - V_{z} $$

Here the known quantities are the 3D velocities at each 2D grid point $V_{x},V_{y},V_{z}$.

We can also numerically compute from the known quantities $\frac{\partial V_{x}}{\partial x}, \frac{\partial V_{y}}{\partial y}$.

Our objective is to compute $h=h(x,y)$ at each grid point.

If it was a 1D problem (i.e. only $x$ dependent) then the problem can be solved using the standard numerical techniques derived from the Euler’s method. But how to approach this problem in the 2D case (i.e. $x,y$ dependent)?

To mathematically simplify the problem we can rewrite the equation above like this:

$$ \frac{\partial h}{\partial x}a(x)+hb(x) + \frac{\partial h}{\partial y}c(y)+hd(y) = - e(x,y) $$ where $a,b,c,d,e$ are numerically known and $h=h(x,y)$ is the objective.

Answer 1

Given what seems to be the physics of the situation (shallow water incompressible flow with $h$ being the depth) surely we have
$$ V_z(x,y,t)\equiv \frac{\partial h(x,y,t)}{\partial t}. $$ So, given the $V(x,y,t)$'s and $h(x,y,t)$ at at time $t$, you have all the ingredients to find $h(x,y, t+\Delta t)$ by Euler's method.

Added Edit: I see now, after the OP's comments, what is being asked. It is not the time evolution that we want, but to solve have the equation $$ ({\bf v}\cdot \nabla) h(x,y) + h(x,y) (\nabla \cdot {\bf v})= {\rm known} $$ Here ${\bf v}= (v_x, v_y)$ is known, and we have that $h=0$ on some boundary. If we want to find $h$ at the point $(x,y)$ we need to start at $x,y$ and compute (numerically) the curve $x(s), y(s)$ such that $x(0)=x$,$y(0)=y$ and $$ \frac{d x}{ds}= -v_x(x,y)\\ \frac{d y}{ds}= -v_y(x,y). $$ and follow it until at $s=s_0$ it reaches the boundary where $h=0$. These curves--- backtracting the fluid flow lines--- are the characteristics. The we solve the first order equation $$ -\frac{d h}{ds}+ (\nabla\cdot {\bf v})h(s) = v_z(x(s),y(s)) $$ with initial data $h(s_0)=0$. Here the diveregence $(\nabla\cdot {\bf v})$ is known at each $s$ because we can compute $\nabla\cdot {\bf v}$ at each point $x(s) y(s)$ from the given velocity-field data. Then $h(x,y)= h(x(0),y(0))$ solves the problem because the derivative along the curve is $$ \frac{d h}{ds} = -({\bf v}\cdot \nabla) h. $$

Now that we know that the is an answer to the problem we can find $h(x,y)$ much more directly by turning it into a matrix problem. For an $N$-by-$N$ grid Treat the knowns and unknowns as arrays $v_i(n,m)$ ($i=x,y,z$), $h(n,m)$ and the arrays as vectors with $N^2$ entries. Then the problem is linear $$ W h= V_z $$ where $W$ is an $N^2$-by$N^2$ matrix arising from discretizing $$ \nabla\cdot (h {\bf v})= V_z. $$ by nearest neighbour differences. Now invert the matrix to write $h=W^{-1} V_z$. As the $W$ matrix is sparse I would recommend using conjugate-gradient algorithm