Let me give a part of your answer, I would need some more indications from your side to answer you fully. So please read, and write some comments so that I can complete my answer.
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).
About the conservative and advective forms of an equation
First, let me clarify your surprise about the form of the equation you are solving.
Advective form
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.
About your particular problem
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"
Edit : added the scheme used by @MichaelAndrewBentley, see comments below.
