Following on from the earlier question I am trying to derive a finite-difference scheme for the advection equation which is conservative. It was suggested that for advection equation with variable velocity,
$$ \frac{\partial u}{\partial x} + \frac{\partial}{\partial x}\left( \boldsymbol{v}(x)u\right) = 0 $$ that the chain rule should not be applied if conservation is required because the flux of $u$ is the fundamental quantity.
The product $\boldsymbol{v}(x)⋅∇u$ is not conservative transport if $\boldsymbol{v}(x)$ is not divergence-free (i.e., constant in 1D). You should stick with the conservative form and resist the urge to apply the chain rule when evaluating conservation properties.
Question
Is the following a reasonable approach for discretization?
- I have kept the flux together (i.e the product $\boldsymbol{v}_ju_j$) while applying the discretization.
- However to write as a matrix equation where $u_j$ is the unknown they eventually must be separated.
- I have a slight concern that this is incorrect because the only difference between the constant velocity case and the varying velocity case is that that $\boldsymbol{v}$ becomes a vector.
- If we applied the chain rule we could get exactly the same equation but with the additional $u \frac{\partial\boldsymbol{v(x)}}{\partial x}$ term which is treated as a source term.
For example, if we apply Crank-Nicolson scheme to the advection equation,
$$ \frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} + \left[ \frac{1-\beta}{2\Delta x} \left( \boldsymbol{v}_{j+1}^{n} u_{j+1}^{n} - \boldsymbol{v}_{j-1}^{n} u_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \boldsymbol{v}_{j+1}^{n+1} u_{j+1}^{n+1} - \boldsymbol{v}_{j-1}^{n+1} u_{j-1}^{n+1} \right) \right] = 0 $$
Using the following substitution, $$ \boldsymbol{r}_j^n = \frac{\boldsymbol{v}_j^n}{2}\frac{\Delta t}{\Delta x} $$ and move the $n+1$ terms to the l.h.s gives,
$$ u_{j}^{n+1} + \beta\boldsymbol{r}_{j+1}^{n+1} u_{j+1}^{n+1} - \beta\boldsymbol{r}_{j-1}^{n+1} u_{j-1}^{n+1} = u_{j}^{n} - (1- \beta)\boldsymbol{r}_{j+1}^{n} u_{j+1}^{n} + (1- \beta)\boldsymbol{r}_{j-1}^{n} u_{j-1}^{n} $$
We can write this as the matrix equation, $$ \begin{pmatrix} 1 & \beta r_2^{n+1} & & 0 \\ -\beta r_1^{n+1} & 1 & \beta r_3^{n+1} & \\ & \ddots & \ddots & \ddots \\ & -\beta r_{J-2}^{n+1} & 1 & \beta r_J^{n+1} \\ 0 & & -\beta r_{J-1}^{n+1} & 1 \\ \end{pmatrix} \begin{pmatrix} u_{1}^{n+1} \\ u_{2}^{n+1} \\ \vdots \\ u_{J-1}^{n+1} \\ u_{J}^{n+1} \end{pmatrix} = \begin{pmatrix} 1 & -(1 - \beta)r_2^n & & 0 \\ (1 - \beta) r_1^n & 1 & -(1 - \beta)r_3^n & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta) r_{J-2}^n & 1 & -(1 - \beta)r_J^n \\ 0 & & (1 - \beta) r_{J-1}^n & 1 \\ \end{pmatrix} \begin{pmatrix} u_{1}^{n} \\ u_{2}^{n} \\ \vdots \\ u_{J-1}^{n} \\ u_{J}^{n} \end{pmatrix} $$
Update I Corrected a sign error suggested by @Geoff-Oxberry.