Diffusion operator doesn't handle discontinuous diffusivity

[majosm]

Making an issue to organize the discussion about this.

[majosm]

Since the path forward for this isn't 100% obvious yet (and @dshtey2 was wondering how to proceed in the meantime), I started thinking about workarounds. I was wondering if a reasonable workaround would be to smooth out alpha a little, such that the nodal values on a shared face are the same on both sides?

I hacked this into my 1D test case (heat equation with Dirichlet BCs, u=0 on the left, u=1 on the right, interface in the middle, 8 elements across, order 2), and it gave results that matched up pretty well with the expected steady state solution (apart from being smooth at the interface instead of having a kink):

I'm not sure how to generalize the setup code for this to higher dimensions, though.

[inducer]

My suggestion is to use one of the schemes from this paper by Ayuso and Marini, for example incomplete interior penalty Galerkin:

@majosm You only need the diffusive part, though for CNS I think we ultimately will want to use the entire advection-diffusion part. (cc @MTCam)

This scheme needs a little bit of rewriting to express the jump [[v]] on the test function, since grudge only lets you express the "local" test function.

In a sense these are overkill, because they guarantee invertibility of the linear system (which we don't need). We could just go with what the Hesthaven/Warbuton book does, which strikes me as fairly ad-hoc. (nailed together out of "gradient" and "divergence" operators) This needs a decision.

@lukeolson Do you agree?

cc @thomasgibson

[inducer]

My postdoc @xywei dug up this paper that claims to have an LDG-based scheme that claims to tolerate discontinuous coefficients.

[majosm]

My postdoc @xywei dug up this paper that claims to have an LDG-based scheme that claims to tolerate discontinuous coefficients.

So I'm still in the middle of trying to implement this, but I noticed something interesting while doing so. I redefined q to be alpha*grad(u) instead of sqrt(alpha)*grad(u), and all of a sudden it gets the right answer:

(Still using internal alpha values only for the fluxes.)

Guessing this has something to do with the fact that q is no longer discontinuous?

[majosm]

Here are some results from applying the diffusion operator to the exact steady-state solution using three different definitions of q.

q = sqrt(alpha)*grad(u) (the current implementation)

q = alpha*grad(u) (the working version from the previous post)

q = grad(u) (the Bassi/Rebay definition)

It's interesting that all three of them compute q correctly, they just subsequently choke when computing u_t. Looking a little more closely still, it appears that it's the flux terms for the u_t equation that are causing trouble:

[majosm]

The above implementations were using interior alpha value when computing the fluxes. If I use fully central fluxes (i.e., (alpha_int*q_int + alpha_ext*q_ext)/2), then the first two don't work, but Bassi/Rebay does:

It seems like this just comes down to ensuring that the fluxes don't have discontinuities?

[majosm]

Here is a table summarizing all of the experiments I tried:

Based on this and our discussion this morning it looks like Bassi/Rebay with full central fluxes is the way to go.

I'm still not sure about the "averaging discontinuous quantities" hypothesis. I suspect it's probably something more subtle than that. One reason to doubt it is I tried initializing my heat equation example with discontinuous u, and it still works (with both 'correct' implementations), even though they average u in order to compute the flux for q. 🤷 Could be because u becomes continuous instantly after that? Dunno.

Edit: Fixed error in table.

[inducer]

I'm still not sure about the "averaging discontinuous quantities" hypothesis.

If you're interested in pursuing this, could precisely state this hypothesis? (I'm still not sure I fully understand what you mean.)

Otherwise, I agree, BR1 seems like a reasonable choice.

[majosm]

If you're interested in pursuing this,

🤷 Happy enough to implement this and move on to other stuff that needs attention, I think. 🙂

could precisely state this hypothesis? (I'm still not sure I fully understand what you mean.)

I'm mostly just looking at these experiments empirically and trying to make associations. By "averaging discontinuous quantity" I meant whether the numerical flux computation is averaging something that is expected to have a jump across the interface (in the analytical sense, not the DG sense). For example, in the second q definition with full central flux, the flux for the q equation is (alpha*u)_avg*normal. At the interface, alpha is discontinuous and u is continuous, so alpha*u will have a jump. Whereas in the third definition you just have u_avg, so no jump. And for the other equation in the third definition + fully central version you have (alpha*q)_avg, where alpha and q are both discontinuous, but their product is the physical flux (up to a negative sign), so no jump there either.

One other thing I also just noticed that seems weird about the first two definitions of q is that if alpha has a jump, their fluxes (sqrt(alpha)*u and alpha*u, respectively) won't be continuous across the interface, even analytically. So in that sense they're kind of non-conservative? Does that have any implications when trying to apply DG to them?

[inducer]

My postdoc @xywei thought about this for a bit and figured out that it is possible to permit discontinuous fluxes even when "sqrt-splitting" the alpha. Here's his write-up: https://notes.tiker.net/s/kAZuOs3gj#

[majosm]

https://notes.tiker.net/s/kAZuOs3gj#

I looked through this, and I have a couple of questions: link