SciML / DiffEqOperators.jl

Linear operators for discretizations of differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/DiffEqOperators/stable/
Other
283 stars 74 forks source link

Feature request: Support Summation-By-Parts operators #506

Open eschnett opened 2 years ago

eschnett commented 2 years ago

Summation By Parts is the discrete equivalent of integrating by parts. This is in many cases a very convenient property to have in a discretization operator. Technically, one needs to choose boundary locations, discretization operators, and the dot product (norm) in a consistent manner.

On a practical level, one needs to choose stencil near boundaries in a particular manner. The straightforward choice to use slightly off-centred stencils with the same order of accuracy does not work.

This paper lists and examines such operators for up to 8th order with various properties. The original idea is much older than this paper, of course.

Different from the current setup, one typically applies the difference stencil to all grid points, i.e. there are no ghost points. Boundary conditions are usually applied weakly, i.e. by modifying the RHS of the PDE that is solved. Is this setup (no ghost points) possible?

xtalax commented 2 years ago

This may prove useful to ongoing work in MethodOfLines.jl as well as here. Do you have a link to the pseudocode as I can't seem to find it on arxiv?

eschnett commented 2 years ago

The source code is distributed as part of the paper source. This link downloads a tarball that contains the LaTeX source, as well as the coefficients in the coeffs directory. I realize that there is no Fortran (or pseudo) code included that would easily tell how to use the coefficients. I will write some.

eschnett commented 2 years ago

Here is some Julia code with comments. Please let me know if you have further questions. I wrote the code for documentation (I did not test it).

ranocha commented 2 years ago

@eschnett If you are interested in SBP operators in Julia, you could look at my package SummationByPartsOperators.jl. I have coded up many coefficients there (but not yours, which I should get fixed...). This package is mostly for teaching and simple research projects, not optimized for high performance with multiple dimensions (yet).