Add option to shift grid by half a point in MOLFiniteDifference for conservative method

[valentinsulzer]

To allow a MOLFiniteDifference method to be conservative, we just need to shift the grid by dx/2 so that x=0 and x=L are aligned with edges of the grid instead of cell centers. In this case the grid would be dx/2:dx:L-dx/2 instead of 0:dx:L.

For example this could be a keyword-argument in MOLFiniteDifference: grid_align="center" (default, current implementation), or "edge".

In terms of implementation, the only difference would be in the boundary conditions:

• Dt(u1) ~ Dx((u2-u1)/dx - a) where a is the flux at the left-hand boundary
• Dt(un) ~ Dx(b - (un-u(n-1))/dx) where b is the flux at the right-hand boundary.

If the boundary conditions are Neumann a and b are given directly. If the boundary conditions are Dirichlet a and b can be calculated by extrapolating to ghost nodes.
Therefore, in the general case

f(u(t,0),Dx(u(t,0)) ~ c

we can replace Dx(u(t,0)) with a and u(t,0) with u1-a*dx/2, which then gives a bc equation for the flux a (and similarly on the right-hand boundary).

(This is like a Finite Volume discretization, but on a regular grid the result is the same)

Will wait for discussions in this issue before trying to implement this

[ChrisRackauckas]

We will want arrays of grid_align per dependent variable for things like common discretizations of Navier-Stokes equations.

In terms of implementation, the only difference would be in the boundary conditions:

The points at which coefficient functions are evaluated also changes.

Other than that, yes agreed.

[valentinsulzer]

Update: Almost all of the existing code can be reused by doing

# Define the grid on which the dependent variables will be evaluated (see #378)
# "center" is recommended for Dirichlet BCs
# "edge" is recommended for Neumann BCs (spatial discretization is conservative)
if grid_align == "center"
       grid = space
       grid_indices = space_indices
elseif grid_align == "edge"
       # boundary conditions implementation assumes centered_order=2
       @assert discretization.centered_order==2
       # construct grid including ghost nodes beyond outer edges
       # e.g. space 0:dx:1 goes to grid -dx/2:dx:1+dx/2
       space_ext = map(s -> vcat(2s[1]-s[2],s,2s[end]-s[end-1]), space)
       grid = map(s -> (s[1:end-1]+s[2:end])/2, space_ext)
       # TODO: allow depvar-specific center/edge choice?
end

and using grid instead of space; then the bcs become equations for the ghost nodes u1 and uend: in the general case

f(u(t,0),Dx(u(t,0)) ~ c

we can replace Dx(u(t,0)) with (u1-u0)/dx and u(t,0) with (u1+u2)/2, which then gives a bc equation for the ghost node u1 (and similarly on the right-hand boundary).

i.e. if grid_align=="center", we have bc equations for u1 and un and interior equations for u2:un-1 (values at dx:dx:L-dx); if grid_align="edges", we have bc equations for u1 and un+1 (ghost nodes) and interior equations for u2:un (values of u at dx/2:dx:L-dx/2)