Implement BR1 scheme

[sloede]

TODO

• Support for AMR
• Support for boundary conditions
• Support for source terms
• Time step calculation for parabolic semi
• Time step calculation for hyperbolic-parabolic semi
• Add LDG (local DG) variant
• Add "warning: experimental" to all exported symbols to allow for breaking changes later

Maybe TODO (or in a subsequent PR)

• Improve support for boundary conditions (i.e., with less copy&paste and fewer hardcoded magic numbers)
• Adequate in-source documentation (mainly docstrings)
• High-level usage information in the docs
• Make SemidiscretizationParabolicAuxVars dimension agnostic
• Implement 3D
• Implement 1D
• Investigate (and improve?) performance
• Implement Navier-Stokes
• Allow gradient calculation with different sets of variables (e.g., primitive or entropy, instead of conservative)

Open questions/issues

• Should we keep the current overall structure? We have a lot of code reuse and feature orthogonality, but pay for it by having to solve ndims + 1 additional full systems (one gradient equation for each dimension + the parabolic system)
• Should we at least make the orientation a part of the gradient equations type?

[ranocha]

Thanks a lot for prototyping this, @sloede!

[ranocha]

As far as I know, the common choice is to use the weak form volume integral for the parabolic terms, even if something more complex is used for the hyperbolic terms. Would it be better to use create_solver_auxvars here?

[andrewwinters5000]

There is no hard set rule for this, it is largely a matter of convention (or data convenience) which form is used. When we use LGL nodes the weak and strong form volume integrals are equivalent so we can use whichever we like. It is just one of those decisions we need to make and then be consistent with it moving forward :)

[ranocha]

What I meant was that we shouldn't use something like shock-capturing volume integrals here

[sloede]

The question for me is: is it the "common" or the "only sensible" choice? For example, we could make it the default to use the weak formulation by add a keyword argument volume_integral_parabolic=VolumeIntegralWeakFormulation() that allows the user to override it.

So if there are is no difference and no proofs depend on using the strong form vs the weak form, I'll hardcode the weak form integral here as well.

Would it be better to use create_solver_auxvars here?

To me, the "auxvars solver" describes the concept of the solver I need to create when computing parabolic problems by getting the gradients as "auxiliary variables" from a separate solver. Here, the shoe doesn't really fit, so for now I'd leave it as is (also given the fact that BR1 is already very solver specific).

[sloede]

What I meant was that we shouldn't use something like shock-capturing volume integrals here

This has convinced me. I made the weak form the default in 665e99f and will allow users to override this with volume_integral_parabolic.

[ranocha]

We will need to have the possibility to use different numerical fluxes for the first and second step (computation of auxiliary variables and doing stuff with them) for methods different from BR1.

[sloede]

Absolutely. That's one of the reasons I wanted the parabolic terms to be handled by a separate solver. I originally planned to postpone this to a second PR, but I think we can also put it in here. Maybe we can discuss this next Monday?

[andrewwinters5000]

There is no hard set rule for this, it is largely a matter of convention (or data convenience) which form is used. When we use LGL nodes the weak and strong form volume integrals are equivalent so we can use whichever we like. It is just one of those decisions we need to make and then be consistent with it moving forward :)

[ranocha]

Suggested change

	solver_parabolic = DGSEM(solver.basis, flux_central, volume_integral_parabolic, solver.mortar)
	solver_parabolic = DGSEM(solver.basis, flux_central, volume_integral_parabolic)

We should probably not use Jesse's EC/ES mortars for parabolic terms...

[ranocha]

There might be three (slightly orthogonal) approaches to solve the current ugliness:

• Reuse less infrastructure by not using the gradient equations as an additional hyperbolic system and re-write the necessary parts in the parabolic semidiscretization.
• Allow passing general parameters (or nothing) to the semidiscretization and from there to the initial conditions, boundary conditions etc. This will be a breaking change unless we try to be very clever to keep backward compatibility. But it would be a great feature to have anyways. Then, we could just pass the heat equation as parameter to the BCs of the gradient equations.
• Switch to a primal formulation instead of a flux formulation and implement that part from scratch.

[sloede]

Thanks for these suggestions! Here are my initial thoughts (just spitballing)...

Reuse less infrastructure by not using the gradient equations as an additional hyperbolic system and re-write the necessary parts in the parabolic semidiscretization.

That would solve the current ugliness of BCs implemented in the gradient equations. But we'd still need to find a way to define/implement BCs for the gradients. How and where would we do that then?

Allow passing general parameters (or nothing) to the semidiscretization and from there to the initial conditions, boundary conditions etc. This will be a breaking change unless we try to be very clever to keep backward compatibility. But it would be a great feature to have anyways. Then, we could just pass the heat equation as parameter to the BCs of the gradient equations.

This sounds like (potentially) a lot of boilerplate and requires even more complicated passing around of "unknown" data (unknown in the sense that a user reading the implementation does not know what is happening). Also, I would hate it if everyone was required to allow for this "extra" argument (that defaults to nothing) - how would it be possible to allow for two function signatures to be used as ICs/BCs? And where would we define what kind of data can be passed to these functions?

Switch to a primal formulation instead of a flux formulation and implement that part from scratch.

That I don't understand. What does it mean to use a primal formulation and what would it look like?

[ranocha]

This sounds like (potentially) a lot of boilerplate and requires even more complicated passing around of "unknown" data (unknown in the sense that a user reading the implementation does not know what is happening). Also, I would hate it if everyone was required to allow for this "extra" argument (that defaults to nothing) - how would it be possible to allow for two function signatures to be used as ICs/BCs? And where would we define what kind of data can be passed to these functions?

I would pass parameters, x, t (or whatever abbreviation of parameters we like) instead of x, t to the ICs/BCs/sources. Being able to pass additional parameters (without type restriction) is exactly what we use from OrdinaryDiffEq, where we pass our semidiscretization as parameters. This will remove the necessity to hard-code parameters inside functions or create callable structs or closures.

That I don't understand. What does it mean to use a primal formulation and what would it look like?

It would look like the FD methods of https://doi.org/10.1007/s10915-009-9301-5. I'm not aware of an article using this approach with AMR, though, but I'm sure it's possible to figure that out.

[jlchan]

Closing; superceded by #1136 and #1156