Setting discontinuous initial conditions

[jlchan]

As far as I can tell, if an initial condition is discontinuous across an element interface, it's ambiguous what value the node on the interface will be set to. It would be nice if there was a way to set discontinuous initial conditions in a non-ambiguous manner, or at least an elixir where compute_coefficients! is overloaded to demonstrate this process.

[sloede]

Hm. The whole point of the equations ecosystem is to be not aware of the underlying discretization. But it is a fair use case you're bringing up - any ideas how to implement this?

IIRC, @ranocha suggested that (for visualization purposes), we "shrink" the elements by something like 1e-12 to avoid ambiguities on element interfaces. I don't know if and how this can be make to work for curved elements, but if yes, we could rely on type traits - i.e., define has_discontinuities(initial_condition_with_jumps) = Val{true} - to active this feature for a particular IC.

[jlchan]

Out of curiosity, why is shrinking the elements is necessary for visualization purposes? We can handle discontinuities fine with the current recipes in main.

One approach I've thought of would be to provide an interface with element information (for example, cell nodal coordinates or even just the cell centroid). I'm not sure about the API for this yet though.

[ranocha]

My preferred solution would be to change the node locations at interfaces slightly (using nextfloat and prevfloat). This tiny change shouldn't change any smooth behavior but allows us to deal with discontinuous data in a well-defined way. On Cartesian meshes, the boundary node location of the left element should be moved to the left via prevfloat, the boundary node location of the right element should be moved to the right via nextfloat. On curved meshes, we should do something similar based on the normal direction at the interface.

I have used this approach successfully before, e.g., https://github.com/ranocha/SummationByPartsOperators.jl/blob/8b9854690d158d68a28fb4ba3df9b9a92206ef17/src/coupling.jl#L87-L98

[sloede]

Not sure if nextfloat and prevfloat are sensitive enough to deal with node locations that are obtained via some mapping. Any experience with that? Other than that, if @gregorgassner and @andrewwinters5000 have no objections, I can see this as a viable option.

A question I'd have is whether this would be default behavior (which might be confusing to user and violate the principle of least astonishment) or opt-in (which makes the code more complex and might lead to other confusion if e.g. visualization or other geometry-based functionality works differently). Slightly leaning towards the latter for now, with the option of eliminating choice if we see no complaints or downsides after a few months of use.

[jlchan]

I don't think it's common to map node locations when working with discontinuous initial conditions. The main conditions I'm thinking of are the nice 4-square discontinuous initial conditions found in Riemann problems.

@ranocha don't you need a sense of absolute direction to use this? How would it work for UnstructuredMesh2D or DGMulti-type meshes?

[jlchan]

Nevermind - read the comment on "normals". I could see that working (shifting the node opposite the direction of the outward normal). I wonder if this will be sensitive to roundoff for meshes with small element sizes, but it sounds OK otherwise.