Save additional node variables in SaveSolutionCallback (#2298)

* Save additional node variables in `SaveSolutionCallback`

* dispatch functions correctly, document new dict

* use vorticity computation

* continue

* working version

* parallel (to be tested!)

* clean up

* mpi tests

* 3d fix

* 3d

* access non-exported func

* comment

* remove doubled method

* simplify

* use singular

* need to figure out which elixirs are MPI-parallel, but single I/O

* add example for coverage

* Update examples/tree_2d_dgsem/elixir_advection_timeintegration_adaptive.jl

* not needed, mpi tests windows tiem out

* send data

* fmt

* add ":limiting_coefficient" as "extra_node_variables"

* Subcell + CNS

* fmt

* fix path

* test vals

* Update src/callbacks_step/save_solution.jl

* comments

* Update src/solvers/dg.jl

Co-authored-by: Benjamin Bolm <[email protected]>

* docu + news entry

* fix merge errors

* fmt

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

---------

Co-authored-by: Benjamin Bolm <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>

Author: 3 people

NEWS.md

@@ -9,8 +9,18 @@ for human readability.
 
 #### Added
 
+- Arbitrary solution-dependent quantities can now be saved in the `SaveSolutionCallback` (and thus later on visualized) ([#2298]).
+
 #### Changed
 
+- When using the `VolumeIntegralSubcellLimiting` with the `SubcellLimiterIDP` the
+  `:limiting_coefficient` must be explicitly provided to the `SaveSolutionCallback` via
+  ```julia
+  save_sol_cb = SaveSolutionCallback(interval = 42,
+                                     extra_node_variables = (:limiting_coefficient,))
+  ```
+  i.e., is no longer automatically saved ([#2298]).
+
 #### Deprecated
 
 #### Removed

examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -80,7 +80,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
 save_solution = SaveSolutionCallback(interval = 300,
                                      save_initial_solution = true,
-                                     save_final_solution = true)
+                                     save_final_solution = true,
+                                     extra_node_variables = (:limiting_coefficient,))
 
 stepsize_callback = StepsizeCallback(cfl = 0.5)
 

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl

@@ -144,7 +144,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = 1000,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
 
 stepsize_callback = StepsizeCallback(cfl = 0.8)
 

examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl

@@ -102,10 +102,52 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+# Add `:vorticity` to `extra_node_variables` tuple ...
+extra_node_variables = (:vorticity,)
+
+# ... and specify the function `get_node_variable` for this symbol, 
+# with first argument matching the symbol (turned into a type via `Val`) for dispatching.
+# Note that for parabolic(-extended) equations, `equations_parabolic` and `cache_parabolic`
+# must be declared as the last two arguments of the function to match the expected signature.
+function Trixi.get_node_variable(::Val{:vorticity}, u, mesh, equations, dg, cache,
+                                 equations_parabolic, cache_parabolic)
+    n_nodes = nnodes(dg)
+    n_elements = nelements(dg, cache)
+    # By definition, the variable must be provided at every node of every element!
+    # Otherwise, the `SaveSolutionCallback` will crash.
+    vorticity_array = zeros(eltype(cache.elements),
+                            n_nodes, n_nodes, # equivalent: `ntuple(_ -> n_nodes, ndims(mesh))...,`
+                            n_elements)
+
+    @unpack viscous_container = cache_parabolic
+    @unpack gradients = viscous_container
+    gradients_x, gradients_y = gradients
+
+    # We can accelerate the computation by thread-parallelizing the loop over elements
+    # by using the `@threaded` macro.
+    Trixi.@threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+
+            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg,
+                                        i, j, element)
+            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg,
+                                        i, j, element)
+
+            vorticity_nodal = vorticity(u_node, (gradients_1, gradients_2),
+                                        equations_parabolic)
+            vorticity_array[i, j, element] = vorticity_nodal
+        end
+    end
+
+    return vorticity_array
+end
+
 save_solution = SaveSolutionCallback(dt = 1.0,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = extra_node_variables) # Supply the additional `extra_node_variables` here
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,

examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl

@@ -132,9 +132,8 @@ function Trixi.analyze(::Val{:energy_potential}, du, u_euler, t,
                                               u_gravity) do u, i, j, element,
                                                             equations_euler, dg,
                                                             equations_gravity, u_gravity
-        u_euler_local = Trixi.get_node_vars(u_euler, equations_euler, dg, i, j, element)
-        u_gravity_local = Trixi.get_node_vars(u_gravity, equations_gravity, dg, i, j,
-                                              element)
+        u_euler_local = get_node_vars(u_euler, equations_euler, dg, i, j, element)
+        u_gravity_local = get_node_vars(u_gravity, equations_gravity, dg, i, j, element)
         # OBS! subtraction is specific to Jeans instability test where rho0 = 1.5e7
         # For formula of potential energy see
         # "Galactic Dynamics" by Binney and Tremaine, 2nd ed., equation (2.18)

examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -95,7 +95,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
 
 stepsize_callback = StepsizeCallback(cfl = 0.6)
 

examples/t8code_3d_dgsem/elixir_euler_ec.jl

@@ -68,9 +68,37 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+# Add `:thermodynamic_entropy` to `extra_node_variables` tuple ...
+extra_node_variables = (:thermodynamic_entropy,)
+
+# ... and specify the function `get_node_variable` for this symbol, 
+# with first argument matching the symbol (turned into a type via `Val`) for dispatching.
+function Trixi.get_node_variable(::Val{:thermodynamic_entropy}, u, mesh, equations,
+                                 dg, cache)
+    n_nodes = nnodes(dg)
+    n_elements = nelements(dg, cache)
+    # By definition, the variable must be provided at every node of every element!
+    # Otherwise, the `SaveSolutionCallback` will crash.
+    entropy_array = zeros(eltype(cache.elements),
+                          ntuple(_ -> n_nodes, ndims(mesh))..., # equivalent: `n_nodes, n_nodes, n_nodes`
+                          n_elements)
+
+    # We can accelerate the computation by thread-parallelizing the loop over elements
+    # by using the `@threaded` macro.
+    Trixi.@threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+            entropy_array[i, j, k, element] = Trixi.entropy_thermodynamic(u_node, equations)
+        end
+    end
+
+    return entropy_array
+end
 save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
-                                     save_final_solution = true)
+                                     save_final_solution = true,
+                                     extra_node_variables = extra_node_variables) # Supply the additional `extra_node_variables` here
 
 stepsize_callback = StepsizeCallback(cfl = 1.0)
 

examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl

@@ -69,10 +69,40 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+# Add `:mach` to `extra_node_variables` tuple ...
+extra_node_variables = (:mach,)
+
+# ... and specify the function `get_node_variable` for this symbol, 
+# with first argument matching the symbol (turned into a type via `Val`) for dispatching.
+function Trixi.get_node_variable(::Val{:mach}, u, mesh, equations, dg, cache)
+    n_nodes = nnodes(dg)
+    n_elements = nelements(dg, cache)
+    # By definition, the variable must be provided at every node of every element!
+    # Otherwise, the `SaveSolutionCallback` will crash.
+    mach_array = zeros(eltype(cache.elements),
+                       n_nodes, n_nodes, # equivalent: `ntuple(_ -> n_nodes, ndims(mesh))...,`
+                       n_elements)
+
+    # We can accelerate the computation by thread-parallelizing the loop over elements
+    # by using the `@threaded` macro.
+    Trixi.@threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            rho, v1, v2, p = prim2cons(u_node, equations)
+            c = sqrt(equations.gamma * p / rho) # speed of sound
+            v_magnitude = sqrt(v1^2 + v2^2)
+
+            mach_array[i, j, element] = v_magnitude / c
+        end
+    end
+
+    return mach_array
+end
 save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = extra_node_variables) # Supply the additional `extra_node_variables` here
 
 amr_indicator = IndicatorHennemannGassner(semi,
                                           alpha_max = 0.5,

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl

@@ -77,7 +77,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
 
 stepsize_callback = StepsizeCallback(cfl = 0.3)
 

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl

@@ -66,7 +66,8 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     solution_variables = cons2prim)
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
 
 save_restart = SaveRestartCallback(interval = 1000,
                                    save_final_restart = true)