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)$.

Since I can calculate the $\gamma_t^x(x,y)$ and $\gamma_t^y(x,y)$ functions and their derivatives exactly, I don't need to approximate these and can just use the product rule to expand the terms with these in. My PDE is therefore of the form

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

I need to use a numeric scheme to approximate the derivatives with $\theta_t(x,y)$ in them. Numeric schemes which work for solving a 2D advection equation (where the gammas are absent and the partial derivatives are simply multiplied by a constant) do not seem to work for me. I have tried several (upwind, downwind, adaptive, etc.), but can't find one which yields conservation of mass for $\theta_{t}$.

Any advice would be greatly appreciated.

• How is $\gamma_t$ unknown ? – Piyush Grover Jul 18 '14 at 17:49
• Thanks for replying. Your question made me realise that my gammas aren't unknown and I can calculate their derivatives exactly - see updated question for details. What I meant was that previously I was using the finite difference product rule - which meant that sometimes both gamma and theta were outside of the grid (i.e. both undefined) and I couldn't solve for my ghost points because I had too many unknowns. Unfortunately, having changed my code to account for this I still don't get conservation of mass for this system and I don't know why! – Michael Andrew Bentley Jul 19 '14 at 23:14
• Gamma needs to satisfy the no-penetration condition (i.e. $\hat n.\gamma_t=0$) at the boundary for mass conservation of $\theta_t$. Is that something you have verified ? – Piyush Grover Jul 20 '14 at 1:43
• Thanks for helping @PiyushGrover. Sorry for being slow, I'm learning the maths here as I go! Do you mean, for example, that $\nabla\gamma_{t}^{x} \cdot \bf{\hat{j}}$ $=0$ along the x-axis (i.e. for all $y=0$)? Which implies $\frac{\partial \gamma_{t}^{x}}{\partial y}=0$ right? I interpret this as there should be no change in $\gamma_{t}^{x}$ the y-direction (i.e. no flow over the x-axis when $y=0$), which sounds reasonable. If this is what you meant, then the test fails, since I calculated $\frac{\partial \gamma_{t}^{x}}{\partial y}\neq0$ when $y=0$. – Michael Andrew Bentley Jul 20 '14 at 17:44
• It depends on the domain.See my answer below. – Piyush Grover Jul 21 '14 at 2:48

The reason you don't get conservation is that you've used the product rule before discretizing, so conservation would require an exact cancellation of truncation errors in the different product terms (which generally won't happen).

Instead, you should directly discretize the conservative form of the equation. Using the usual notation for finite differences and setting $$(\gamma^x,\gamma^y)=(u,v)$$, this would be something like:

$$\theta^{n+1}_{i,j} = \theta^{n-1}_{i,j} - \frac{\Delta t}{\Delta x}(\theta^n_{i+1,j}u^n_{i+1,j} - \theta^n_{i-1,j}u^n_{i-1,j}) - \frac{\Delta t}{\Delta y}(\theta^n_{i,j+1}v^n_{i-1,j+1} - \theta^n_{i,j}u^n_{i-1,j}).$$

It's easy to check that if you sum $$\theta^{n+1}$$ over the whole grid, all of the fluxes cancel out except for those at the boundaries. Here I've used centered differences in time and space. This will be fine if your initial data is smooth and well-resolved, and if second-order accuracy (in $$\Delta t$$ and $$\Delta x$$) is sufficient. For more complicated or demanding situations, you might refer to (for instance) Randall LeVeque's Finite Volume Methods for Hyperbolic Problems and the Clawpack software.

Note: In the literature I'm familiar with, the conservation law you've written is referred to as the 2D advection equation. The non-conservative version is sometimes given the same name, or may be referred to as the color equation.

Too long for a comment, so answering here: Call $\gamma_t=[\gamma^x_t\:\:\gamma^y_t]'$

Integrating the original equation over the domain, say $\Omega$, we get

$\dfrac{d||\theta_t||}{dt}=-\int_{\Omega} \nabla.(\gamma_t \theta_t) d\Omega$

Now using Divergence theorem on R.H.S, this gives

R.H.S $=-\int_{\partial\Omega} (\hat n.\gamma_t)\theta_t dS$, where $\hat n$ is the normal to

boundary, and $\partial\Omega$ is the boundary.

Hence, for mass conservation, i.e. L.H.S to be zero, we need the dot product on R.H.S to be identically zero.