Inconsistency in RiemannSolver interface

Right now we define and use it as follows

• build a rsolver = ApproxRiemannSolver(flux, speed, equation) object
• somewhere in the code call it with broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
• optional: BoundaryConditions are built from an ApproxRiemannSolver instance and are latter called with broadcast_boundaries!(bdryconds, rsolver, cache, mesh)

Leaving out the BoundaryConditions case for now: The ApproxRiemannSolver tries to mimic a @with_signature implementation, hence, we split out the lax_friedrich_flux definition into a separate method. However, it seems unintuitive to not define which flux method to use when construction an ApproxRiemmannSolver. This also adds code duplication and potentially bugs, because we (should) use the same solver in all instances, or shouldn't we?

Regarding BoundaryConditions: For this to work we should not really have to worry about the actual flux method and instead just provide an instance of ApproxRiemannSolver.

It is not urgent, but we should keep it in mind for the next refactoring.