Closed boundary conditions in finite difference method for diffusive-advective equation

Question

I am implementing a finite difference method in solving the diffusive-advective equation: $$ u_t + v \cdot u_x = D\cdot u_{xx} $$ (v, D are constants). Planning to use the operator splitting method (see below), I am now focusing on the advective part.

I chose staggered leapfrog method (from Numerical Recipes book) $$ u^{n+1}_j = u^{n-1}_j - \frac{v\cdot\Delta t }{\Delta x}(u^n_{j+1} - u^n_{j-1}) \;\;\;\; $$ because it should be Courant-stable in case of conserved flux; this turns out to be the case, since it works perfectly using periodic conditions. The problems show when trying to change the boundaries. My question, in brief, is: how to explicitly implement these conditions in such a method? Below, I am going to provide more details on my failings.

No-flux conditions

Let us start from a more general case: as it is pointed out here for diffusion-advection equation, Robin conditions are required to achieve closed boundaries $$ v\cdot u - Du_x = 0 $$ I tried to discretize them: $$ v\cdot u_j - \frac{D}{\Delta x} (c_j - c_{j-1}) = 0 $$ However, when D=0 the only condition remaining is that $$ v\cdot u = 0 $$ thus it is to me unclear how to interpret this: imposing u=0 would result in the (failing) strategy for absorbing conditions (see below), while v is constant (and I don't understand where I could impose it equals to zero).

Absorbing conditions

Leapfrog method is meant for flux-conserving situations. In fact, imposing

u(j=0,n) = 0 
u(jmax,n) = 0

creates artificial instabilities. The leapfrog updating has got three parts: one for the previous (two steps before) value in the point, one increasing this value proportionally to the uphill point, the other reducing proportionally to the downhill one. This is also what one could expect from a physical point of view: some quantity enters, some goes out. The reason why the simple imposition of null boundaries fails is now evident: if the downhill point is null, there is no reducing term, so the point before the last explodes; the point before it has now a big reducing term, so goes to zero; the point before then has a small reducing term and so on, resulting in an alternance of increasing and decreasing terms (see figure 1).

Operator splitting

As explained in the Numerical Recipes book (cited above), an equation like $$ u_t = \mathcal{L} u $$ where a decomposition in linear operators $$ \mathcal{L}=\sum^m_{i}\mathcal{L}_i $$ is possible can be resolved applying, for each step, the sequence $$ \mathcal{U}_1(\mathcal{U}_2(\dots\mathcal{U}_m(u^n)\dots)) = u^{n+1} $$ provided that each operator U solves for a term in the sum and is stable.

---

Even if I know alternatives are possible, which contemplate diffusion-advection function in general, before digging into them I am willing to solve this doubt of mine.

Answer 1

If you consider only the pure advective part and let, say, $v>0$, then you can obtain "no flux condition" at the right point only by setting (redefining) $v=0$ there. Of course, the solution can then develop a "boundary layer" there and it can become discontinuous at the right point as the solution will preserve there the value given by initial condition, i.e. $u_{jmax}^{n+1} = u^0_{jmax}$.

A natural boundary condition for constant speed $v>0$ is that you let the solution "flow out" through the right point of interval that can be viewed as outflow boundary. It can be obtained by applying your leapfrog method for $u_{jmax}^{n+1}$ with the auxiliary value $u_{jmax+1}^n= u_{jmax}^n$ (a constant extrapolation) or $u_{jmax+1}^n= 2 u_{jmax}^n - u_{jmax-1}^n$ (a linear extrapolation). One can see such treatment as "do nothing boundary conditions" that is appropriate for outflow boundary and pure advection. It is like you are solving the pure advection (the differential equation) also at the right point.

Just a remark - if you use a proper upwind method instead of leapfrog method, you even need no auxiliary values outside of your domain to compute $u_{jmax}^{n+1}$.

For the inflow boundary, i.e. the left point in the case $v>0$, you must prescribe the values $u_0^{n+1}$ directly by "Dirichlet boundary conditions" and not solving the differential equation.

For the case of advection-diffusion equation you can in theory define the no-flux boundary conditions as you write with $D \neq 0$, but it means for $v>0$ that you prescribe an "inward diffusion", consequently a corresponding gradient, that will cancel the advection flux, so the total flux will be zero. Such situation is hard to find in practice. More typical is to prescribe something only for the diffusion flux and to use the previous ``do nothing'' approach (or an extrapolation) for the advective flux.