Fix higher order upwind differences and separate centered and upwind order specification in MOLFiniteDifference

[archermarx]

This PR addresses #308 and in doing so separates the specification of difference order for upwind and central differences in the MOLFiniteDifference constructor, allowing one to do this:

dx = 0.01
order = 1
discretization_1 = MOLFiniteDifference(dx, order)     # Sets upwind difference and central difference order to 1
discretization_2 = MOLFiniteDifference(dx)   # Sets both difference orders to 2 by default
discretization_3 = MOLFiniteDifference(dx; upwind_order = 2, centered_order = 4)  # Specify different orders for each difference type

This fixes the issue in #308 where specifying higher order upwind differences doesn't do anything:

Before (1D Linear advection)

After

4th order pure upwinding leads to massive spurious oscillations, the correction of which requires the use of a biased upwind operator or a flux/slope limiter of some sort, which is work for another time

[ChrisRackauckas]

I think it's natural to default to 2nd order central and 1st order upwinding?

[archermarx]

Yep you're right. I've been swamped with work but will try to push this along soon.

[archermarx]

I changed how the default order is specified so that if you type MOLFiniteDifference(dxs; order = 2), for example, both orders get set to 2, but if you just type MOLFiniteDifference(dxs) you get a default 1st order upwind and 2nd order central diff

[ChrisRackauckas]

@emmanuellujan can you review?

[emmanuellujan]

I like the idea. Here are some comments: ideally you should be able to specify the numerical scheme of each derivative of each of your equations. Particularly, you should be able to decide if you want to use upwind or central difference, and the order of each scheme for each derivative. Currently, this is not implemented. The reason why this prototype code was using, as default, first-order upwind in first-order derivatives is due to the stability of the convection equation, which is a very commonly used equation. (A more detailed description can be found here [1]. If u Δt / Δx = 1, where u represents the speed at which information propagates, we get the exact solution, and we avoid having numerical viscosity. Also, If you implement the code here [2], you can change it a little bit to compare upwind vs central difference, and see the impact of changing dt, dx, and u). A similar reasoning can be used to show the convenience of discretizing diffusion equations, i.e. second-order derivatives, using central differences. In conclusion, I think it is reasonable to set as default, first-order upwind for odd-order derivatives, and central difference for even-order derivatives. Currently I am working with the MOL_discretization.jl, I will add this contribution to my code.
 
[1] http://bender.astro.sunysb.edu/classes/numerical_methods/lectures/advection.pdf
[2] https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/01_Step_1.ipynb

[archermarx]

Yeah i totally understand the reasoning, and we can see in the second-order upwind derivative we no longer have a total variation-preserving situation in the presence of shocks, so a first-order default is highly reasonable. Still, in problems where we do not expect to feature strong shocks, a second-order upwind is sometimes useful (biased upwind operators may also be useful). Ideally, we'd have some form of slope/flux limiter so that we can mix and match numerical scheme and limiter to get the desired behavior for the problem at hand. Is there any plan to implement slope limiters in the future so we can get TVD, 2nd-order behavior with advection problems? I'm open to trying this out myself if there isn't

[ChrisRackauckas]

Slope limiters are already implemented. That's a separate issue since it's in the time stepping method, so it's in:

https://diffeq.sciml.ai/stable/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)

[ChrisRackauckas]

Sounds like this is a good improvement, so I'm merging. Thanks @archermarx!

[archermarx]

Great, happy to help! And super glad to hear we have slope limiters. I'm gonna try migrating my existing code over to the PDEsystem framework soon. I consciously aped a lot of DifferentialEquations.jl design patterns when designing it so it'll be great to be working with the real thing.