# Advice on numerical solution for 2D hyperbolic PDE with zero flux boundary conditions

I would like to numerically solve a hyperbolic PDE of the form

$\frac{\partial\theta_t}{\partial t}(x,y)+\frac{\partial\left[\theta_t \gamma_t^x\right]}{\partial x}(x,y)+\frac{\partial\left[\theta_t \gamma_t^y\right]}{\partial y}(x,y)=0,$

which is very similar to the 2D advection equation, except that the partial derivatives w.r.t $x$ and $y$ are products of my dependent variable $\theta_t(x,y)$ with another function $\gamma_t^x(x,y)$ or $\gamma_t^y(x,y)$.

I am a biologist by training and new to this type of thing, so any advice (at all!) would be much appreciated. I have written a straightforward explicit upwind differencing scheme in Matlab to numerically solve this (difference equation below). My solution almost works (I'm replicating a result from a paper), but the density is not conserved, and much of it seems 'attracted' to the boundaries. I suspect I have implemented my zero-flux boundary conditions incorrectly, and this is what I would like advice on.

$\theta_{t+\Delta t}(x,y) \approx \theta_{t}(x,y) -\Delta t \left(\frac{\theta_{t}(x+\Delta x,y)\gamma_t^x(x+\Delta x,y)-\theta_{t}(x,y)\gamma_t^x(x,y)}{\Delta x} + \frac{\theta_{t}(x,y+\Delta y)\gamma_t^y(x,y+\Delta y)-\theta_{t}(x,y)\gamma_t^y(x,y)}{\Delta y} \right).$

In particular, because I have a product of two functions that are all dependent on $t$, $x$ and $y$, it seems harder to solve for a zero-flux boundary condition because I have extra unknowns at my ghost points.

Thanks!

Clarification: Here's how I applied the finite difference product rule (from wikipedia: $\Delta(fg)=f\Delta g + g\Delta f + \Delta g \Delta f$), please correct me if I've used it incorrectly. To keep it short I've just focussed on $x$ and ignored $y$ and dropped the denominator (since everything is over $\Delta x$):

$\frac{\partial\left[\theta_t \gamma_t^x\right]}{\partial x}(x)\approx\theta_t(x)\left[\gamma_t^x(x+\Delta x)-\gamma_t^x(x)\right]+\gamma_t^x(x)\left[\theta_t(x+\Delta x)-\theta_t(x)\right]+\left[\theta_t(x+\Delta x)-\theta_t(x)\right]\left[\gamma_t^x(x+\Delta x)-\gamma_t^x(x)\right].$

After expanding brackets I get

$=\theta_t(x)\gamma_t^x(x+\Delta x)-\theta_t(x)\gamma_t^x(x)+\theta_t(x+\Delta x)\gamma_t^x(x)-\theta_t(x)\gamma_t^x(x) -\theta_t(x)\gamma_t^x(x+\Delta x) -\theta_t(x+\Delta x)\gamma_t^x(x)+\theta_t(x+\Delta x)\gamma_t^x(x+\Delta x)+\theta_t(x)\gamma_t^x(x),$

And after cancelling I get

$=\theta_t(x+\Delta x)\gamma_t^x(x+\Delta x)-\theta_t(x)\gamma_t^x(x).$

# About notations in numerical methods

There are a few mistakes in the way you write your equations. These are only details, but for someone who is used to it, it can a bit disturbing.

• $\theta_t$, is either written as : $\theta(t)$ if written in an analytical equation or $\theta^n := \theta(t^n)$ if written in a time-discretised equation
• $q(x,y)$ is fine when working in an analytical equation, but it should be written $q_{i,j}$ if working in a space-discretised equation. $q$ stands for whichever quantity that may apply (that notation is used by a part of the CFD community, mathematicians would use $u$ instead).

First, let me clarify your surprise about the form of the equation you are solving.

The advection equation can be written in different forms. The form you seem to be most familiar with is called the advective form :

$\partial_t q + (\underline{u} \cdot \underline{\nabla}) q = 0$

Which would give, if written in the partial derivative formulation :

$\partial_t q + u_x \partial_x q + u_y \partial_y q = 0$

After space-discretising :

$\partial_t q_{i,j}(t) + u_{i,j} [D_x q]_{i,j} + v_{i,j} [D_y q]_{i,j} = 0$

With :

• $q$ : the "working variable" i.e. the variable you want to integrate
• $\underline{u} := u \underline{e}_x + v \underline{e}_y$ a speed function
• $D_x$ and $D_y$ your differentiation schemes. For instance : $D_x q := \frac{q_{i+1} - q_{i}}{\Delta x}$

## Conservative form

The form of the advection equation you have is called the conservative form or the divergent form :

$\partial_t q + \underline{\nabla} \cdot \underline{f}(q) = 0$

Using partial derivatives :

$\partial_t q + \partial_x f(q) + \partial_y f(q) = 0$

After space-discretising :

$\partial_t q_{i,j}(t) = - \frac{1}{\Delta x} \left [ f(q_{i+\frac{1}{2},j}) - f(q_{i-\frac{1}{2},j}) \right ] - \frac{1}{\Delta y} \left [ f(q_{i,j+\frac{1}{2}}) - f(q_{i,j-\frac{1}{2}}) \right ]$

Where :

• $f(q)$ is a numerical flux. In your case, $f(\theta^n) := \theta^n \gamma^n$. Note that the numerical evaluation of $f(q)$ requires to know values of $q$ located between nodes e.g. $q_{i+\frac{1}{2},j}$

I am a bit confused that you seem to indicate that $\gamma_t^x \neq \gamma_t^y$... Could you give me more information on this, please ?

The conservative and advective forms correspond to exactly the same analytical equation i.e. an advection equation.

## Why different forms ?

There are many different forms of the advection equation. All forms correspond to the exact same analytical equation. However, each form may not conserve the quantities as well as others. In particular, the conservative form is very good at conserving the working variable ($\theta$ in your case). Let me get more into details :

Let $I^n$ be defined as $I^n := \sum_{i=0}^{N_i} \sum_{j=0}^{N_j} \theta_{i,j}^n$ with:

• $N_i$, $N_j$ the total number of columns resp. lines in your computational domain
• $n$ the number of the time step
• $i$, $j$ the number of the column resp. line of the current node in the sum

If you use the conservative form, then you will have $I^n = I^0 \forall n$. For example, $I^{1000} = I^0$, which basically means that you have "as much" $q$ at iteration $n=1000$ than when you started your simulation (i.e. $n=0$).

This is perhaps the behaviour you are expecting, but in finite differences strange things happen. In particular, quantity may vanish or appear. Indeed, if you were using the advective form then you would most likely have $I^n \neq I^0$.

If you want to read a discussion on the conservative properties of an equation, you can refer to : Morinishi 1995, "Conservative properties of finite difference schemes for incompressible flow".

Note that we are speaking finite differences here. Even when using the conservative form of the advection equation, we discretise in the style of finite differences. It would also be possible to do finite volumes, but we are not doing that here.

You seem to want :

• an Euler Explicit time integration scheme
• a conservative form of the equation
• no-flux BC

## Euler Explicit time integration scheme

$q^{n+1} = q^{n} + rhs \, \Delta t$

## A conservative form of the equation

$\partial_t \theta_{i,j}(t) = - \frac{1}{\Delta x} \left [ f(\theta_{i+\frac{1}{2},j}) - f(\theta_{i-\frac{1}{2},j}) \right ] - \frac{1}{\Delta y} \left [ f(\theta_{i,j+\frac{1}{2}}) - f(\theta_{i,j-\frac{1}{2}}) \right ]$

With: $f(\theta^n) := \theta^n \gamma^n$

If second order accuracy is good enough for you, then :

$\theta_{i+\frac{1}{2},j} = \frac{\theta_{i+1,j} + \theta_{i,j}} { 2 } + O(\Delta x^2)$

## No-flux BC

The best way to apply it is to use a layer of mirror nodes, as you did. That is, we need to set $q_{-1,j} := q_{0,j}$ on the left border and do the same on all four borders. I work in C, hence the bottom left node of my "physical" domain is $(0,0)$ and the top-right node is $(N_i-1,N_j-1)$. In other words, $q_{0,j}$ is a "physical" node and $q_{-1,j}$ is a ghost node (or mirror node in our case).

# Scheme chosen

The scheme that was chosen by OP @MichaelAndrewBentley is 2nd order Lax-Wendroff, which is of the form :

$u^{n+1} = b_{-1} u_{i-1}^n + b_0 u_i^n + b_1 u_{i+1}^n + O(\Delta t^2) + O(\Delta x^2)$.

With : $b_{-1} := \frac{1}{2} c (1+c)$, $b_0 := 1-c^2$, $b_1 := -\frac{1}{2} c (1-c)$, $c := \frac{a \Delta t}{\Delta x}$

More details can be found in LeVeque's 2002 book "Finite Volume Methods for Hyperbolic Problems"

• Thank you very much @G.Lorieul for such a detailed explanation, it was all very helpful, including the reference, which I just read. The notation guidance is particularly appreciated - I will switch to the CFD notation when I write up. – Michael Andrew Bentley Aug 21 '14 at 18:58
• I will try to clarify the question you asked. The functions $\gamma^x(x,y)$ and $\gamma^y(x,y)$ are different flux functions for my dependent variable $\theta(x,y,t)$ in the $x$ and $y$ direction respectively. They are functions of space only (not time as it looks like I originally indicated). If I combine my notation with your notation for the flux vector I would have $\underline{f}(\theta(x,y,t))=(\gamma^x(x,y),\gamma^y(x,y))$ would that be correct? In divergence form I would then have $\partial_t \theta(x,y,t)+\underline{\nabla}\cdot\underline{f}(\theta(x,y,t))=0$ – Michael Andrew Bentley Aug 21 '14 at 18:58
• Finaly, as it happens I got my code working yesterday. I used a 2nd order Lax-Wendroff finite volume method described in Leveque (2002) 'Finite volume methods for hyperbolic problems'. The form is as you specified in your answer though. I have used operator splitting to deal with the dimensions and ghost cells to define the boundary conditions. Sound good? Any comments? Thanks again! – Michael Andrew Bentley Aug 21 '14 at 18:59
• @MichaelAndrewBentley : yes, this does sound good. BTW, LeVeque's book is a great reference. Just make sure you do not get confused between finite differencing of an equation in conservative form and finite volume methods. The numerical equations of finite volumes/differences are very much alike, however they are not the same. That being said, using Lax-Wendroff for your problem is fine. Just know that it will blow up if a discontinuity appears in your solution ; more details in Leveque's book. BTW, I would very much appreciate if you gave the ref of the article you are trying to replicate. – Gael Lorieul Aug 22 '14 at 7:20
• Also, matlab can be ok for prototyping, but it is not the right tool to be using if you intend to run heavy simulations. That is, matlab runs much slower as other languages. In particular, people tend to prefer either C/C++ or fortran90 (which is a completely different language as compared to fortran77). You can get an idea of the performance of some of the most-popular tools in the table here : julialang.org Be careful that benchmarking different languages is not an easy task, and hence if two languages perform more or less the same, do not try to guess which is the best of both. – Gael Lorieul Aug 22 '14 at 7:23

If we expand the derivatives, the new equation reads $$\frac{\partial\theta(x,y)}{\partial t} + \theta(x,y)\frac{\partial \gamma^x(x,y)}{\partial x} + \gamma^x \frac{\partial \theta(x,y)}{\partial x} + \theta(x,y)\frac{\partial \gamma^y(x,y)}{\partial y} + \gamma^y \frac{\partial \theta(x,y)}{\partial y} = 0 \enspace .$$ If we know use discretization scheme that you proposed, we get $$\theta_{t+\Delta t}(x,y) = \theta_{t}(x,y) - \Delta t\left[\frac{\theta_t(x,y)\gamma_t^x(x+\Delta x,y) + \gamma_t^x(x,y)\theta_t(x+\Delta x, y) - 2\theta_t(x,y)\gamma_t^x(x,y)}{\Delta x} + \frac{\theta_t(x,y)\gamma_t^y(x, y + \Delta y) + \gamma_t^y(x,y)\theta_t(x, y+\Delta y) - 2\theta_t(x,y)\gamma_t^y(x,y)}{\Delta y}\right] \enspace .$$ I think that you are not taking into account the product rule of finite differences.

You also mentioned that you might be applying the boundary conditions in a wrong way. How are you applying the boundary conditions?

• thanks for your help. I did actually use the FD product rule and found that some terms cancelled leaving the equation I wrote above. For my boundary conditions I have tried creating a single layer of ghost points around the grid containing the same values as those in the 'real' locations next to them. I also tried implementing if statements to specify no loss of mass if I am on the border of my domain (just a rectangle). I've been trying to solve for the values of my ghost cells, but I can't because I have products with unknown gamma(x+deltax) in them. – Michael Andrew Bentley Jul 17 '14 at 19:34
• Well, if you used the FD product rule I cannot see why our expressions are note the same. – nicoguaro Jul 18 '14 at 2:35
• I've added a clarification to my question to show how I derived it. Thanks. – Michael Andrew Bentley Jul 18 '14 at 8:26