Potential cooperation and integration between DiffEqOperators and SummationByPartsOperators

[xtalax]

Hi,

I'm working on DiffEqOperators.jl and am noticing that we have a lot of common functionality between these packages. Of note, we have general multidimensional handling, and more general boundary conditions. We don't however have Sumnation by parts derivative operators, or diffusion operators, but we would like to.

I was wondering if you could support with developing conversion routines between your DerivativeOperator and our DerivativeOperator, to allow for use of these operators with the rest of the growing ecosystem that uses DiffEqOperators, and to enable people to use them with multidimensional PDEs?

Also, I note that your operators are size (N, N), are they valid at [i, end] and [i, 1] without contributions from boundary conditions for example? This could help improve the stability of our neumann and robin boundary conditions.

[ranocha]

Hi @xtalax! I would be happy to see my package being used in DiffEqOperators.jl etc. and would be happy to help. Right now, this package does not provide multidimensional operators (you would need to convert them, e.g., via sparse and form the Kronecker products on your own).
Yes, the basic idea behing summation by parts (SBP) operators is to use adapted stencils near the boundaries to get an approximation of the first/second/... derivative without incorporating boundary data strongly (in contrast to the approach of DiffEqOperators.jl). However, you still need to inform the method about the boundary conditions, of course. This is done using weakly imposed boundary conditions, i.e., you add some additional terms to your discretization. I do not see any way to derive the form of stable boundary procedures programatically in general. Hence, I am not sure how "automatic" this conversion could be. If you just focus on a specific PDE, you can probably find some way of stable boundary conditions with SBP operators in the literature, e.g., for linear advection, the wave equation etc. You can find an example of this setup in the tutorials. Note that

mul!(du, D, u, -one(eltype(D)))

computes an approximation of the (negative) first derivative everywhere without using boundary data. Then,

du[begin] += (uL_func(t) - u[begin]) / left_boundary_weight(D)

adds an additional term to impose the inflow boundary condition weakly.

[xtalax]

As long as we can derive an expression Dx(u[begin]) => dot(D.lower_index_stencil, u[begin:length(D.lower_index_stencil)]) we can do neumann I think, diff eq operators rearranges this expression to get the value at the ghost node. MethodOfLines uses the same dot expression to create a replacement rule for the symbolic derivative in any boundary condition equation specified.

I'd like to convert to our derivative operator directly, then it should just work with all the multidimensional stuff, take a look at our derivative operator, you'll probably find it quite familiar.

I had a question about how the weights are stored, how would I get the weights for the interior, and for 1:num_boundary_points points away from u[begin]? Similarly for u[end]?

[ranocha]

I had a question about how the weights are stored, how would I get the weights for the interior, and for 1:num_boundary_points points away from u[begin]? Similarly for u[end]?

What do you mean by weights? The coefficients of the derivative operator?

[ranocha]

As long as we can derive an expression Dx(u[begin]) => dot(D.lower_index_stencil, u[begin:length(D.lower_index_stencil)]) we can do neumann I think

Not necessarily in a stable way. A homogeneous Neumann boundary condition for the wave equation can be written like

SummationByPartsOperators.jl/src/second_order_eqs/wave_eq.jl

Lines 48 to 49 in 21d463f

	if left_bc == Val(:HomogeneousNeumann)
	@inbounds ddu[1] += derivative_left(derivative, u, Val(1)) / left_boundary_weight(derivative)

The expression for the heat equation would basically be similar. However, this way of setting a homogeneous Neumann boundary condition does not result in a stable method for the BBM-BBM system, see Theorem 4.7 of https://dx.doi.org/10.4208/cicp.OA-2020-0119.