Evolve: align stepper logic between all ERK implementations...

Evolve: align stepper logic between all ERK implementations (!189)

This assumes that the first row in the `a` matrix of the Butcher tableau for ERK methods has only zeros, which seems to be the case for the all the methods have atm.

---

Problem explanation:

So far the LSERK4 and the remaining ERK steppers worked slightly different.
The former used for the first substep the initial data array `u0`, and `u!`  for the remaining steps. That shortcut avoids a copy and appears natural for the LSERK4 algorithm.
OTOH, in the ERK stepper we first copy `u0` into `u!`, then accumulate some substep `k`s into `u!`, and then call the RHS with `u!`. (the reason we always copied `u0` over was that we did not account for the above assumption on `a` before)

The problem was now that in the GRHD evolutions we sometimes alter the state vector in the RHS computation.
This influenced the results to the extent that we could not run the (newly added) reftest only with the LSERK4 method.
What was happening was that when using the ERK method we never altered `u0` directly, only `u!`, whereas `u0` was altered in case the GRHD RHS decided to do so on the first step.
Subsequent substeps of the ERK method depend upon `u0`, and so this broke things.

src/evolve.jl

@@ -225,33 +225,24 @@ Stepper for Low Storage five-stage fourth-order Explicit Runge Kutta scheme.
 Taken from Hesthaven, Warburton (2007), section 3.4.
 """
 function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::COEFFS_LSERK4, callback_substep)
-
   @unpack nstages, a, b, c = coeffs
-
   kim1, ki = tmp_k
-  N = length(u!)
 
   for i = 1:nstages
-
     ti = t + c[i] * dt
     F(ki, i == 1 ? u0 : u!, ti)
-
-    @. ki *= dt
+    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
+    !success && return false
     if i == 1
-      @. u! = u0 + b[i] * ki
+      @. u! = u0 + b[i] * dt * ki
     else
       ai = a[i]
       @. ki += ai * kim1
-      @. u! += b[i] * ki
+      @. u! += b[i] * dt * ki
     end
-
-    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
-    !success && return false
-
     if i < nstages
       @. kim1 = ki
     end
-
   end # for
 
   return true
@@ -266,23 +257,23 @@ Stepper for Explicit Runge Kutta scheme.
 function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::TimeStepAlgorithm, callback_substep)
   @unpack nstages, a, b, c = coeffs
 
-  for n = 1:nstages
-    tn = t + c[n] * dt
+  for i = 1:nstages
+    ti = t + c[i] * dt
+    # assumes that a[1,:] are all zeros
+    F(tmp_k[i], i == 1 ? u0 : u!, ti)
+    i == nstages && break
+    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
+    !success && return false
     u! .= u0
-    for i = 1:n-1
-      u! .+= a[n, i] * tmp_k[i]
+    for j = 1:i
+      @. u! += a[i+1, j] * dt * tmp_k[j]
     end
-    F(tmp_k[n], u!, tn)
-    tmp_k[n] .*= dt
-    success = callback_substep(u!, tn, n)
-    !success && return false
   end
 
   u! .= u0
-  for n = 1:nstages
-    u! .+= b[n] * tmp_k[n]
+  for i = 1:nstages
+    @. u! += b[i] * dt * tmp_k[i]
   end
-
   success = callback_substep(u!, t+dt, nstages)
 
   return true

test/IntegrationTests/refs/grhd_bondi_infall_avmda/grhd_bondi_infall_avmda.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 0.02143872868864732
+polytrope_gamma = "$(4/3)"
+
+[GRHD]
+bc = "from_id"
+id = "bondi_accretion_infall"
+variables0d_analyze_error = ["ρ", "D"]
+variables1d_analyze_error = ["ρ", "D"]
+analyze_error_reference_solution = "bondi_accretion"
+c2p_enforce_causal_atm = false
+
+[Mesh]
+range = [ 1.8, 20.0 ]
+n = 7
+k = 16
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d = [ "D", "Sr", "τ" ]
+aligned_ts = "$(collect(range(1.0,10.0,step=1.0)))"
+
+[Evolution]
+cfl = 0.2
+tend = 10.0
+method = "rk4"
+
+[HRSC]
+method = "av"
+av_method = "mda"
+

test/IntegrationTests/refs/testconfig.toml

@@ -126,6 +126,9 @@ reltol1d = 1e-6
 variables0d  = [ "D_err_norm" ]
 variables1d  = [ "D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]
 
+[grhd_bondi_infall_avmda]
+variables1d  = [ "D", "Sr", "τ", ]
+
 [fv_grhd_bondi_accretion]
 variables1d  = [ "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]