0'th order Polynomials / First Order Finite Volume

[DanielDoehring]

Following the discussion in the Trixi Student's slack channel & the advantages of a TVD spatial discretization, I think it's a good idea to the implement constant polynomials.
As mentioned in the Slack discussion, this is not really possible using Gauss-Lobatto nodes but rather Gauss-Legendre.
That said, doing this properly requires IMO the implementation of this new basis and in principle a completely new approach, shifting from nodal to modal.
An advantage of the modal approach is also the possibility to realize relatively easily limiting for linear polynomials.

I already tried hacking this into the existing code

Changes in Trixi.jl/src/solvers/dgsem/basis_lobatto_legendre.jl

# From FLUXO (but really from blue book by Kopriva)
function gauss_lobatto_nodes_weights(n_nodes::Integer)
  # From Kopriva's book
  n_iterations = 10
  tolerance = 1e-15

  # Initialize output
  nodes = zeros(n_nodes)
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1

+  if N > 0
    # Calculate values at boundary
    nodes[1] = -1.0
    nodes[end] = 1.0
    weights[1] = 2 / (N * (N + 1))
    weights[end] = weights[1]

    # Calculate interior values
    if N > 1
      cont1 = pi/N
      cont2 = 3/(8 * N * pi)

      # Use symmetry -> only left side is computed
      for i in 1:(div(N + 1, 2) - 1)
        # Calculate node
        # Initial guess for Newton method
        nodes[i+1] = -cos(cont1*(i+0.25) - cont2/(i+0.25))

        # Newton iteration to find root of Legendre polynomial (= integration node)
        for k in 0:n_iterations
          q, qder, _ = calc_q_and_l(N, nodes[i+1])
          dx = -q/qder
          nodes[i+1] += dx
          if abs(dx) < tolerance * abs(nodes[i+1])
            break
          end
        end

        # Calculate weight
        _, _, L = calc_q_and_l(N, nodes[i+1])
        weights[i+1] = weights[1] / L^2

        # Set nodes and weights according to symmetry properties
        nodes[N+1-i] = -nodes[i+1]
        weights[N+1-i] = weights[i+1]
      end
    end

    # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
    if n_nodes % 2 == 1
      _, _, L = calc_q_and_l(N, 0)
      nodes[div(N, 2) + 1] = 0.0
      weights[div(N, 2) + 1] = weights[1] / L^2
    end
+  else # N = 0: Use Gauss-LEGENDRE (NOT Lobatto; Fallback case)
+    # https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature
+    #nodes[1] = 0.0
+    #weights[1] = 2.0
+
+    # Some dummy values (apparently not used)
+    nodes[1] = 42.0
+    weights[1] = 42.0
+  end

  return nodes, weights
end

but it is a bit difficult to actually judge whether that is valid.
Reason for this is that the convergence for the elixir_euler_convergence_pure_fv.jl testcase is roughly order 1 in $L^2$ norm, but not in $L^\infty$:

Small convergence study

This is done by changing the number of gridcells

 16 Cells (initial_refinement_level=4)

 L2 error:       3.50685142e-02   4.63749047e-02   7.40729755e-02
 Linf error:     4.92996392e-02   6.47267216e-02   1.04992214e-01


 32 Cells (initial_refinement_level=5)

 L2 error:       2.41263401e-02   2.85610126e-02   3.45726622e-02
 Linf error:     3.41023005e-02   4.05880510e-02   5.00591691e-02


 64 Cells (initial_refinement_level=6)

 L2 error:       1.44475748e-02   1.65376964e-02   1.63696274e-02
 Linf error:     2.06680979e-02   2.34825289e-02   2.38935423e-02


 128 Cells (initial_refinement_level=7)

 L2 error:       7.97491581e-03   9.09979034e-03   8.32683878e-03
 Linf error:     1.14997467e-02   1.30339731e-02   1.20321406e-02

What is promising, however, is that the oscillations present in the $N\geq 1$ case are gone:

$N=1$:

$N=0$:

For constant polynomials the computation of the Volume Term can/should be omitted.

[DanielDoehring]

What also requires a closer look is if the initial condition gets for this hack actually properly integrated, i.e.,

$$U_i(t_0) = \frac{1}{\Delta x} \int_{x_{i-\frac12}}^{x_{i+\frac12}} u_0(x) \mathrm d x$$

or if just the value at cell center is assigned:

$$ U_i(t_0) = u_0(x_i). $$

[DanielDoehring]

Addressed in #1489