Correct RHS by subtracting RHS of a steady state solution

[efaulhaber]

Some simulations (like the baroclinic instability, for which I specifically implemented this feature in #979) rely on a steady state being preserved for a long time, while the solution is only slightly diverging from the steady state.

If these simulations have a source term that causes well-balanced errors, it may be necessary to use a relatively high resolution just to preserve the steady state. This problem can be mitigated by evaluating and saving the RHS of the exact steady-state solution at the beginning of the simulation, and subtracting it in every RHS evaluation. By construction, this steady state will then be preserved exactly (to machine rounding errors, of course).

The baroclinic instability, for example, produces entirely unusable results for a resolution of 8x8x4 cells per cube face with a polydeg of 3. Only with a polydeg of 4 or higher, the results are slowly becoming usable.
This feature makes even the polydeg 3 simulation produce usable (although badly resolved) results, and eliminates most of the grid imprinting in the polydeg 4 simulation.
Here is a comparison of the polydeg 3 simulation with (right) and without (left) this feature.

Here is a picture of how it's supposed to look like (16x16x8 cells per cube face, polydeg 5). Note that this is probably not at the same simulation time (I had these pictures lying around). This picture is just to show that the right picture looks surprisingly good.

Do we want to have this feature in Trixi?
#979 contains a proof-of-concept implementation of this feature.

[ranocha]

I'm not sure whether we really need this feature since it's basically a couple of lines like

u_steady_state = compute_coefficients(steady_state_function, 0.0, semi)
let du_steady_state = similar(u_steady_state)
  Trixi.rhs!(du_steady_state, u_steady_state, semi, 0.0)
  global function corrected_rhs!(du, u, semi, t)
    Trixi.rhs!(du, u, semi, t)
    du .-= du_steady_state
  end
end
u0 = compute_coefficients(0.0, semi)
ode = ODEProblem(new_rhs!, u0, semi, tspan))

(assuming that no AMR is used).

From my point of view, this appears to be a very specialized feature that quite easy to implement in an elixir - if somebody knows what they want to have, it's clear what to do.

[efaulhaber]

Wow, this is amazing! I haven't even remotely thought about doing something like this. That's way easier and cleaner than the stuff I did in my PR. Right now, I don't need it to work with AMR. Maybe we can come up with a solution that allows for AMR later.

Why exactly do we need the let block here?

[ranocha]

Why exactly do we need the let block here?

We don't need it to make this work, just to make it fast since we would be using a non-constant global variable (du_steady_state) without the let block.

[efaulhaber]

I see, thanks! This one's as fast as my complicated implementation in #979:

u_steady_state = compute_coefficients(steady_state_baroclinic_instability, 0.0, semi)
let du_steady_state = similar(u_steady_state)
  Trixi.rhs!(du_steady_state, u_steady_state, semi, 0.0)

  global function corrected_rhs!(du, u, semi, t)
    Trixi.rhs!(du, u, semi, t)
    Trixi.@trixi_timeit Trixi.timer() "rhs correction" begin
      @batch for i in eachindex(du)
        du[i] -= du_steady_state[i]
      end
    end
  end
end

u0 = compute_coefficients(0.0, semi)
ode = ODEProblem(corrected_rhs!, u0, tspan, semi)