Merge branch 'main' into Validation_dam

.github/workflows/ci.yml

@@ -12,6 +12,7 @@ on:
       - 'NEWS.md'
       - 'README.md'
       - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
       - 'docs/**'
   pull_request:
       paths-ignore:
@@ -22,6 +23,7 @@ on:
       - 'NEWS.md'
       - 'README.md'
       - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
       - 'docs/**'
   workflow_dispatch:
 

docs/src/callbacks.md

@@ -4,3 +4,13 @@
 Modules = [TrixiParticles]
 Pages = map(file -> joinpath("callbacks", file), readdir(joinpath("..", "src", "callbacks")))
 ```
+
+# [Custom Quantities](@id custom_quantities)
+
+The following pre-defined custom quantities can be used with the
+[`SolutionSavingCallback`](@ref) and [`PostprocessCallback`](@ref).
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = ["general/custom_quantities.jl"]
+```

src/callbacks/post_process.jl

@@ -12,10 +12,15 @@ anonymous as the function name will be used as part of the name of the value.
 The callback can be triggered either by a fixed number of time steps (`interval`) or by
 a fixed interval of simulation time (`dt`).
 
-
 # Keywords
-- `funcs...`: Functions to be executed at specified intervals during the simulation.
-              The functions must have the arguments `(v, u, t, system)`.
+- `funcs...`:   Functions to be executed at specified intervals during the simulation.
+                Each function must have the arguments `(v, u, t, system)`,
+                and will be called for every system, where `v` and `u` are the
+                wrapped solution arrays for the corresponding system and `t` is
+                the current simulation time. Note that working with these `v`
+                and `u` arrays requires undocumented internal functions of
+                TrixiParticles. See [Custom Quantities](@ref custom_quantities)
+                for a list of pre-defined functions that can be used here.
 - `interval=0`: Specifies the number of time steps between each invocation of the callback.
                 If set to `0`, the callback will not be triggered based on time steps.
                 Either `interval` or `dt` must be set to something larger than 0.
@@ -34,15 +39,11 @@ a fixed interval of simulation time (`dt`).
 
 # Examples
 ```jldoctest; output = false
-function example_function(v, u, t, system)
-    println("test_func ", t)
-end
-
 # Create a callback that is triggered every 100 time steps
-postprocess_callback = PostprocessCallback(interval=100, example_quantity=example_function)
+postprocess_callback = PostprocessCallback(interval=100, example_quantity=kinetic_energy)
 
 # Create a callback that is triggered every 0.1 simulation time units
-postprocess_callback = PostprocessCallback(dt=0.1, example_quantity=example_function)
+postprocess_callback = PostprocessCallback(dt=0.1, example_quantity=kinetic_energy)
 
 # output
 ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
@@ -373,56 +374,3 @@ function add_entry!(pp, entry_key, t, value, system_name)
     # Add the new entry to the list
     push!(entries, value)
 end
-
-function kinetic_energy(v, u, t, system)
-    return sum(each_moving_particle(system)) do particle
-        velocity = current_velocity(v, system, particle)
-        return 0.5 * system.mass[particle] * dot(velocity, velocity)
-    end
-end
-
-function total_mass(v, u, t, system)
-    return sum(each_moving_particle(system)) do particle
-        return system.mass[particle]
-    end
-end
-
-function max_pressure(v, u, t, system)
-    return maximum(particle -> particle_pressure(v, system, particle),
-                   each_moving_particle(system))
-end
-
-function min_pressure(v, u, t, system)
-    return minimum(particle -> particle_pressure(v, system, particle),
-                   each_moving_particle(system))
-end
-
-function avg_pressure(v, u, t, system)
-    if n_moving_particles(system) == 0
-        return 0.0
-    end
-
-    sum_ = sum(particle -> particle_pressure(v, system, particle),
-               each_moving_particle(system))
-    return sum_ / n_moving_particles(system)
-end
-
-function max_density(v, u, t, system)
-    return maximum(particle -> particle_density(v, system, particle),
-                   each_moving_particle(system))
-end
-
-function min_density(v, u, t, system)
-    return minimum(particle -> particle_density(v, system, particle),
-                   each_moving_particle(system))
-end
-
-function avg_density(v, u, t, system)
-    if n_moving_particles(system) == 0
-        return 0.0
-    end
-
-    sum_ = sum(particle -> particle_density(v, system, particle),
-               each_moving_particle(system))
-    return sum_ / n_moving_particles(system)
-end

src/callbacks/solution_saving.jl

@@ -4,6 +4,7 @@
                            output_directory="out", append_timestamp=false, max_coordinates=2^15,
                            custom_quantities...)
 
+
 Callback to save the current numerical solution in VTK format in regular intervals.
 Either pass `interval` to save every `interval` time steps,
 or pass `dt` to save in intervals of `dt` in terms of integration time by adding
@@ -28,33 +29,34 @@ To ignore a custom quantity for a specific system, return `nothing`.
 - `custom_quantities...`:       Additional user-defined quantities.
 - `write_meta_data`:            Write meta data.
 - `verbose=false`:              Print to standard IO when a file is written.
-- `max_coordinates=2^15`        The coordinates of particles will be clipped if their absolute values exceed this threshold.
+- `max_coordinates=2^15`:       The coordinates of particles will be clipped if their
+                                absolute values exceed this threshold.
+- `custom_quantities...`:   Additional custom quantities to include in the VTK output.
+                            Each custom quantity must be a function of `(v, u, t, system)`,
+                            which will be called for every system, where `v` and `u` are the
+                            wrapped solution arrays for the corresponding system and `t` is
+                            the current simulation time. Note that working with these `v`
+                            and `u` arrays requires undocumented internal functions of
+                            TrixiParticles. See [Custom Quantities](@ref custom_quantities)
+                            for a list of pre-defined custom quantities that can be used here.
 
 # Examples
 ```jldoctest; output = false, filter = [r"output directory:.*", r"\s+│"]
-# Save every 100 time steps.
+# Save every 100 time steps
 saving_callback = SolutionSavingCallback(interval=100)
 
-# Save in intervals of 0.1 in terms of simulation time.
+# Save in intervals of 0.1 in terms of simulation time
 saving_callback = SolutionSavingCallback(dt=0.1)
 
-# Additionally store the norm of the particle velocity for fluid systems as "v_mag".
-using LinearAlgebra
-function v_mag(v, u, t, system)
-    # Ignore for other systems.
-    return nothing
-end
-function v_mag(v, u, t, system::WeaklyCompressibleSPHSystem)
-    return [norm(v[1:ndims(system), i]) for i in axes(v, 2)]
-end
-saving_callback = SolutionSavingCallback(dt=0.1, v_mag=v_mag)
+# Additionally store the kinetic energy of each system as "my_custom_quantity"
+saving_callback = SolutionSavingCallback(dt=0.1, my_custom_quantity=kinetic_energy)
 
 # output
 ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
 │ SolutionSavingCallback                                                                           │
 │ ══════════════════════                                                                           │
 │ dt: ……………………………………………………………………… 0.1                                                              │
-│ custom quantities: ……………………………… [:v_mag => v_mag]                                                │
+│ custom quantities: ……………………………… [:my_custom_quantity => TrixiParticles.kinetic_energy]           │
 │ save initial solution: …………………… yes                                                              │
 │ save final solution: ………………………… yes                                                              │
 │ output directory: ………………………………… *path ignored with filter regex above*                           │

src/general/custom_quantities.jl

@@ -0,0 +1,96 @@
+"""
+    kinetic_energy
+
+Returns the total kinetic energy of all particles in a system.
+"""
+function kinetic_energy(v, u, t, system)
+    if n_moving_particles(system) == 0
+        return 0.0
+    end
+
+    return sum(each_moving_particle(system)) do particle
+        velocity = current_velocity(v, system, particle)
+        return 0.5 * system.mass[particle] * dot(velocity, velocity)
+    end
+end
+
+"""
+    total_mass
+
+Returns the total mass of all particles in a system.
+"""
+function total_mass(v, u, t, system)
+    return sum(each_moving_particle(system)) do particle
+        return system.mass[particle]
+    end
+end
+
+"""
+    max_pressure
+
+Returns the maximum pressure over all particles in a system.
+"""
+function max_pressure(v, u, t, system)
+    return maximum(particle -> particle_pressure(v, system, particle),
+                   each_moving_particle(system))
+end
+
+"""
+    min_pressure
+
+Returns the minimum pressure over all particles in a system.
+"""
+function min_pressure(v, u, t, system)
+    return minimum(particle -> particle_pressure(v, system, particle),
+                   each_moving_particle(system))
+end
+
+"""
+    avg_pressure
+
+Returns the average pressure over all particles in a system.
+"""
+function avg_pressure(v, u, t, system)
+    if n_moving_particles(system) == 0
+        return 0.0
+    end
+
+    sum_ = sum(particle -> particle_pressure(v, system, particle),
+               each_moving_particle(system))
+    return sum_ / n_moving_particles(system)
+end
+
+"""
+    max_density
+
+Returns the maximum density over all particles in a system.
+"""
+function max_density(v, u, t, system)
+    return maximum(particle -> particle_density(v, system, particle),
+                   each_moving_particle(system))
+end
+
+"""
+    min_density
+
+Returns the minimum density over all particles in a system.
+"""
+function min_density(v, u, t, system)
+    return minimum(particle -> particle_density(v, system, particle),
+                   each_moving_particle(system))
+end
+
+"""
+    avg_density
+
+Returns the average_density over all particles in a system.
+"""
+function avg_density(v, u, t, system)
+    if n_moving_particles(system) == 0
+        return 0.0
+    end
+
+    sum_ = sum(particle -> particle_density(v, system, particle),
+               each_moving_particle(system))
+    return sum_ / n_moving_particles(system)
+end

src/general/general.jl

@@ -25,3 +25,4 @@ include("initial_condition.jl")
 include("system.jl")
 include("interpolation.jl")
 include("file_system.jl")
+include("custom_quantities.jl")

src/visualization/write2vtk.jl

@@ -1,32 +1,42 @@
 """
     trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="",
-              write_meta_data=true, custom_quantities...)
+              write_meta_data=true, max_coordinates=Inf, custom_quantities...)
 
 Convert Trixi simulation data to VTK format.
 
 # Arguments
-- `vu_ode`: Solution of the TrixiParticles ODE system at one time step. This expects an `ArrayPartition` as returned in the examples as `sol`.
+- `vu_ode`: Solution of the TrixiParticles ODE system at one time step.
+            This expects an `ArrayPartition` as returned in the examples as `sol.u[end]`.
 - `semi`:   Semidiscretization of the TrixiParticles simulation.
 - `t`:      Current time of the simulation.
 
 # Keywords
-- `iter`:                 Iteration number when multiple iterations are to be stored in separate files.
-- `output_directory`:     Output directory path. Defaults to `"out"`.
-- `prefix`:               Prefix for output files. Defaults to an empty string.
-- `write_meta_data`:      Write meta data.
-- `custom_quantities...`: Additional custom quantities to include in the VTK output. TODO.
-- `max_coordinates=Inf`   The coordinates of particles will be clipped if their absolute values exceed this threshold.
-
+- `iter=nothing`:           Iteration number when multiple iterations are to be stored in
+                            separate files. This number is just appended to the filename.
+- `output_directory="out"`: Output directory path.
+- `prefix=""`:              Prefix for output files.
+- `write_meta_data=true`:   Write meta data.
+- `max_coordinates=Inf`     The coordinates of particles will be clipped if their absolute
+                            values exceed this threshold.
+- `custom_quantities...`:   Additional custom quantities to include in the VTK output.
+                            Each custom quantity must be a function of `(v, u, t, system)`,
+                            which will be called for every system, where `v` and `u` are the
+                            wrapped solution arrays for the corresponding system and `t` is
+                            the current simulation time. Note that working with these `v`
+                            and `u` arrays requires undocumented internal functions of
+                            TrixiParticles. See [Custom Quantities](@ref custom_quantities)
+                            for a list of pre-defined custom quantities that can be used here.
 
 # Example
 ```jldoctest; output = false, setup = :(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), tspan=(0.0, 0.01), callbacks=nothing))
 trixi2vtk(sol.u[end], semi, 0.0, iter=1, output_directory="output", prefix="solution")
 
+# Additionally store the kinetic energy of each system as "my_custom_quantity"
+trixi2vtk(sol.u[end], semi, 0.0, iter=1, my_custom_quantity=kinetic_energy)
+
 # output
 
 ```
-
-TODO: example for custom_quantities
 """
 function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="",
                    write_meta_data=true, max_coordinates=Inf, custom_quantities...)