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

Question

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.

Answer 1

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.

Answer 2

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.