Fixing the conventions

[ChrisRackauckas]

The issues that @francispoulin pointed out are pretty gruesome. While things "work", what was discovered is that the operators we build are changing conventions left and right. This is because "the standard" operators tend to act this way: people naturally discretize differently for different BCs. @shivin9 I think was warning me about this, but I didn't realize until way to late...

To fix this issue, I triaged with @dextorious and @jagot and we tentatively have a plan:

1. Always assume the discretization on on-grid unless by default.
2. Always assume the boundaries are the first and last points by default.

Here's the damage. Note that Neumann is unchanged by these choices.

• Currently Dirichlet0 and Dirichlet assume the boundaries are off the grid. We can instead take the values of the boundary to be u[1] and u[end]. Using this definition, the operator needs no hidden constant to operate correctly since all of the information is in the vector. We define for higher order methods via constant extrapolation of the end point value. Should we keep Dirichlet0 and document it, given its importance as the stencil matrix?
• Periodic BCs would then have to operate a little differently. The current operator can be renamed to CenteredPeriodic since it's the center-grid discretization for periodic BCs. Then, we need to modify the stencil so at the boundaries it skips over the repeat. For example, at the left end it should be u[end-2],u[1],u[2] instead of grabbing from the other end. It should check if u[1] != u[end] and warn accordingly.
• Document that Neumann works for center-grid discretization as well.
• It's clear how Robin BCs are defined for this, but what we're doing probably doesn't follow it.
• Rename be BC argument to derivative since now in both Dirichlet and Robin the BC value comes from the end point of the array itself.

Why these choices?

1. Having the end points in the grid makes for easier plotting.
2. It makes Dirichlet with time dependence a lot simpler to implement and understand.
3. The total memory cost is trivial.

This needs to get comments before action.

[dextorious]

I'm in favor of all the proposed changes because the CPU/memory costs are trivial and I think it's important to have a consistent, easy to understand convention for users. Once we have that, if there's interest, we can work towards exposing more customizability, but it's important to get the basics right and to make the user-facing API as simple as reasonably possible.

It's very important that we get feedback now, especially from users who oppose these changes. Ultimately, a preference for on-grid vs off-grid boundaries is a somewhat subjective tradeoff and there is no globally optimal choice. We need to make a decision one way or the other for consistency, but there are two ways this can go so if you're a user of this package and prefer a different set of conventions, please speak up now.

[ChrisRackauckas]

Here's a visualization of the different choices and the effect on the operator:

Dirichlet0: Boundaries off the array (special case). On-grid
A = [2 1 0 0 0
     1 2 1 0 0
     0 1 2 1 0
     0 0 1 2 1
     0 0 0 1 2]
Intuition: The "1" off the grid is multiplied by the BC (0) so it's safely ignored.

Dirichlet: Boundaries off the array. On-grid
A = [2 1 0 0 0
     1 2 1 0 0
     0 1 2 1 0
     0 0 1 2 1
     0 0 0 1 2]
+ B which is zeros(u) with B[1] and B[end] the BCs

Intuition: The BC is off the grid so we need to add it to the end points for the "1" that was off the grid.

Dirichet: Boundaries on the array. On-grid
A = [1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1]
Intuition: The BC is in the array, so `A*u` with just the stencil works.

Dirichet: Boundaries off the array. Center-grid (constant extrapolation)
Doesn't work.
A = [1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1]
Seems like it would, but it doesn't actually have the BC in the calculation so
it would be unstable.

Neumann0: Boundaries on the array. On-grid
A = [2 2 0 0 0 0 0
     1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1
     0 0 0 0 0 2 2]
Intuition: Reflect the "1" the falls off over the boundary to get "2 2"

Neumann0: Boundaries off the array. Center-grid
A = [2 2 0 0 0 0 0
     1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1
     0 0 0 0 0 2 2]
Intuition: Reflect the "1" the falls off over the boundary to get "2 2"

Periodic: Boundaries on the array. Center-grid
A = [2 1 0 0 0 0 1
     1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1
     1 0 0 0 0 1 2]
Intuition: The "1" wraps around to the other side

Periodic: Boundaries on the array. On-grid
A = [2 1 0 0 0 1 0
     1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1
     0 1 0 0 0 1 2]
Intuition: The "1" wraps around to the other side, but u[1]=u[end] so you need
to skip it.

Periodic: Boundaries on the left side of the array. On-grid
A = [2 1 0 0 0 0 1
     1 2 1 0 0 0 0
     0 1 2 1 0 0 0
     0 0 1 2 1 0 0
     0 0 0 1 2 1 0
     0 0 0 0 1 2 1
     1 0 0 0 0 1 2]
Intuition: The "1" wraps around to the other side.

Some points to note:

1. For Neumann, on-grid and center grid is always the same. Neumann derivatives always need "hidden values" for the derivatives.
2. For Dirichlet, center grid doesn't work.
3. For Dirichlet with non-zero BC with BCs off the array, the operation is not a matrix multiplication.
4. The memory difference having the boundaries on and off is just 2 for the arrays. Since the actual stencils we are using are lazy, there's no cost. So here the difference shown is 5x5 vs 7x7, but that's pseudo because we don't make the matrix!
5. Dirichlet with values on grid is no longer tridiagonal (but still banded)
6. Periodic with bcs on grid is a little weird looking but all of them are not banded, so it doesn't really hurt the sparsity pattern.

As you can see from the visualization, Dirichlet with values off is the "standard" version people use, Neumann under either interpretation is "standard", periodic with a center discretization is "standard".

[ChrisRackauckas]

Closing this because it's incorrect. Solving for points on the boundary will make the derivative numerically unstable, so we instead want to do SciML/DifferentialEquations.jl#260 (comment)