Add dvdu_ode vector to custom_quantities (#879)

* add dvu_ode vector to custom quantities

* add for fluid and boundary

* fix tests

* fix tests

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* Combine PST and TVF into a unified framework (#884)

* Combine PST and TVF into a unified framework

* Require update callback for PST

* Fix WCSPH

* Update PST only in callback

* Fix EDAC

* Update docs

* Fix alle example files

* Fix tests

* Fix periodic channel

* Fix docs

* Update news

* Fix tests

* Fix example file

* Rename 'tlsph' to 'place_on_shell' (#814)

* rename

* format

* forgot some

* format

* naming

* forgot some more

* fix test

* incorporate review comments

* format

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

* Combine PST and TVF into a unified framework (#884)

* Combine PST and TVF into a unified framework

* Require update callback for PST

* Fix WCSPH

* Update PST only in callback

* Fix EDAC

* Update docs

* Fix alle example files

* Fix tests

* Fix periodic channel

* Fix docs

* Update news

* Fix tests

* Fix example file

* fix tests

* fix doctests

* implement suggestions

* fix tests

* fix tests

* implement suggestions

* fix

* fix

---------

Co-authored-by: LasNikas <niklas.nehe@web.de>
Co-authored-by: Sven Berger <berger.sven@gmail.com>
Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>

Author: 4 people

examples/postprocessing/interpolation_plane.jl

@@ -104,10 +104,10 @@ combined_plot = Plots.plot(plot1, plot2, plot3, plot_3d, layout=(2, 2),
                            size=(1000, 1500), margin=3mm)
 
 # If we want to save planes at regular intervals, we can use the postprocessing callback.
-# Note that the arguments `system, v_ode, u_ode, semi, t` are more powerful than the
+# Note that the arguments `system, dv_ode, du_ode, v_ode, u_ode, semi, t` are more powerful than the
 # documented arguments `system, data, t`, allowing us to use interpolation (which requires
 # a semidiscretization).
-function save_interpolated_plane(system, v_ode, u_ode, semi, t)
+function save_interpolated_plane(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     # Size of the patch to be interpolated
     interpolation_start = [0.0, 0.0]
     interpolation_end = [tank_size[1], tank_size[2]]

examples/postprocessing/postprocessing.jl

@@ -13,9 +13,9 @@ using CSV
 using DataFrames
 using JSON
 
-# Any custom function with the arguments `v, u, t, system` can be passed to the callback
+# Any custom function with the arguments `system, data, t` can be passed to the callback
 # to be called every 10th timestep. See example below:
-function hello(system, v_ode, u_ode, semi, t)
+function hello(system, data, t)
     # Will write "hello" and the current simulation time
     println("hello at ", t)
 

src/TrixiParticles.jl

@@ -23,7 +23,8 @@ using ReadVTK: ReadVTK
 using RecipesBase: RecipesBase, @series
 using Random: seed!
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
-                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!
+                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!,
+                 get_du
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt

src/callbacks/post_process.jl

@@ -229,6 +229,7 @@ end
 # `affect!`
 function (pp::PostprocessCallback)(integrator)
     @trixi_timeit timer() "apply postprocess cb" begin
+        dv_ode, du_ode = get_du(integrator).x
         vu_ode = integrator.u
         v_ode, u_ode = vu_ode.x
         semi = integrator.p
@@ -251,7 +252,7 @@ function (pp::PostprocessCallback)(integrator)
             system_index = system_indices(system, semi)
 
             for (key, f) in pp.func
-                result = custom_quantity(f, system, v_ode, u_ode, semi, t)
+                result = custom_quantity(f, system, dv_ode, du_ode, v_ode, u_ode, semi, t)
                 if result !== nothing
                     add_entry!(pp, string(key), t, result, filenames[system_index])
                     new_data = true

src/callbacks/solution_saving.jl

@@ -154,6 +154,7 @@ function (solution_callback::SolutionSavingCallback)(integrator)
     (; interval, output_directory, custom_quantities, write_meta_data, git_hash,
      verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback
 
+    dvdu_ode = get_du(integrator)
     vu_ode = integrator.u
     semi = integrator.p
     iter = get_iter(interval, integrator)
@@ -173,7 +174,7 @@ function (solution_callback::SolutionSavingCallback)(integrator)
         println("Writing solution to $output_directory at t = $(integrator.t)")
     end
 
-    @trixi_timeit timer() "save solution" trixi2vtk(vu_ode, semi, integrator.t;
+    @trixi_timeit timer() "save solution" trixi2vtk(dvdu_ode, vu_ode, semi, integrator.t;
                                                     iter, output_directory, prefix,
                                                     write_meta_data, git_hash=git_hash[],
                                                     max_coordinates, custom_quantities...)

src/callbacks/steady_state_reached.jl

@@ -131,11 +131,12 @@ end
 
     vu_ode = integrator.u
     v_ode, u_ode = vu_ode.x
+    dv_ode, du_ode = get_du(integrator).x
     semi = integrator.p
 
     # Calculate kinetic energy
     ekin = sum(semi.systems) do system
-        return kinetic_energy(system, v_ode, u_ode, semi, 0)
+        return kinetic_energy(system, dv_ode, du_ode, v_ode, u_ode, semi, 0)
     end
 
     if length(previous_ekin) == interval_size

src/general/custom_quantities.jl

@@ -3,7 +3,7 @@
 
 Returns the total kinetic energy of all particles in a system.
 """
-function kinetic_energy(system, v_ode, u_ode, semi, t)
+function kinetic_energy(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
 
     # TODO: `current_velocity` should only contain active particles
@@ -17,18 +17,20 @@ function kinetic_energy(system, v_ode, u_ode, semi, t)
     end
 end
 
-kinetic_energy(system::BoundarySystem, v_ode, u_ode, semi, t) = zero(eltype(system))
+function kinetic_energy(system::BoundarySystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
+    return zero(eltype(system))
+end
 
 """
     total_mass
 
 Returns the total mass of all particles in a system.
 """
-function total_mass(system, v_ode, u_ode, semi, t)
+function total_mass(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return sum(system.mass)
 end
 
-function total_mass(system::BoundarySystem, v_ode, u_ode, semi, t)
+function total_mass(system::BoundarySystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     # It does not make sense to return a mass for boundary systems.
     # The material density and therefore the physical mass of the boundary is not relevant
     # when simulating a solid, stationary wall. The boundary always behaves as if it had
@@ -47,12 +49,12 @@ end
 
 Returns the maximum pressure over all particles in a system.
 """
-function max_pressure(system::FluidSystem, v_ode, u_ode, semi, t)
+function max_pressure(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     return maximum(current_pressure(v, system))
 end
 
-function max_pressure(system, v_ode, u_ode, semi, t)
+function max_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
 
@@ -61,12 +63,12 @@ end
 
 Returns the minimum pressure over all particles in a system.
 """
-function min_pressure(system::FluidSystem, v_ode, u_ode, semi, t)
+function min_pressure(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     return minimum(current_pressure(v, system))
 end
 
-function min_pressure(system, v_ode, u_ode, semi, t)
+function min_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
 
@@ -75,13 +77,13 @@ end
 
 Returns the average pressure over all particles in a system.
 """
-function avg_pressure(system::FluidSystem, v_ode, u_ode, semi, t)
+function avg_pressure(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     sum_ = sum(current_pressure(v, system))
     return sum_ / nparticles(system)
 end
 
-function avg_pressure(system, v_ode, u_ode, semi, t)
+function avg_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
 
@@ -90,12 +92,12 @@ end
 
 Returns the maximum density over all particles in a system.
 """
-function max_density(system::FluidSystem, v_ode, u_ode, semi, t)
+function max_density(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     return maximum(current_density(v, system))
 end
 
-function max_density(system, v_ode, u_ode, semi, t)
+function max_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
 
@@ -104,12 +106,12 @@ end
 
 Returns the minimum density over all particles in a system.
 """
-function min_density(system::FluidSystem, v_ode, u_ode, semi, t)
+function min_density(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     return minimum(current_density(v, system))
 end
 
-function min_density(system, v_ode, u_ode, semi, t)
+function min_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
 
@@ -118,12 +120,12 @@ end
 
 Returns the average_density over all particles in a system.
 """
-function avg_density(system::FluidSystem, v_ode, u_ode, semi, t)
+function avg_density(system::FluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
     sum_ = sum(current_density(v, system))
     return sum_ / nparticles(system)
 end
 
-function avg_density(system, v_ode, u_ode, semi, t)
+function avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end

src/io/write_vtk.jl

@@ -48,15 +48,27 @@ trixi2vtk(sol.u[end], semi, 0.0, iter=1, my_custom_quantity=kinetic_energy)
 
 ```
 """
-function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="",
-                   write_meta_data=true, git_hash=compute_git_hash(),
+function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out",
+                   prefix="", write_meta_data=true, git_hash=compute_git_hash(),
+                   max_coordinates=Inf, custom_quantities...)
+
+    # The first argument is not necessary in most cases. Since it is usually not available to the user,
+    # this API wrapper makes it optional.
+    # Note that custom quantities using the fluid acceleration will not work and return NaN acceleration.
+    return trixi2vtk(fill!(similar(vu_ode), NaN), vu_ode, semi, t; iter, output_directory,
+                     prefix, write_meta_data, git_hash, max_coordinates,
+                     custom_quantities...)
+end
+
+function trixi2vtk(dvdu_ode, vu_ode, semi, t; iter=nothing, output_directory="out",
+                   prefix="", write_meta_data=true, git_hash=compute_git_hash(),
                    max_coordinates=Inf, custom_quantities...)
     (; systems) = semi
-    v_ode, u_ode = vu_ode.x
 
     # Update quantities that are stored in the systems. These quantities (e.g. pressure)
     # still have the values from the last stage of the previous step if not updated here.
     @trixi_timeit timer() "update systems" begin
+        v_ode, u_ode = vu_ode.x
         # Don't create sub-timers here to avoid cluttering the timer output
         @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
     end
@@ -67,24 +79,26 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix
         system_index = system_indices(system, semi)
         periodic_box = get_neighborhood_search(system, semi).periodic_box
 
-        trixi2vtk(system, v_ode, u_ode, semi, t, periodic_box;
+        trixi2vtk(system, dvdu_ode, vu_ode, semi, t, periodic_box;
                   system_name=filenames[system_index], output_directory, iter, prefix,
                   write_meta_data, git_hash, max_coordinates, custom_quantities...)
     end
 end
 
 # Convert data for a single TrixiParticle system to VTK format
-function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_directory="out",
-                   prefix="", iter=nothing, system_name=vtkname(system_),
-                   write_meta_data=true, max_coordinates=Inf, git_hash=compute_git_hash(),
-                   custom_quantities...)
+function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
+                   output_directory="out", prefix="", iter=nothing,
+                   system_name=vtkname(system_), write_meta_data=true, max_coordinates=Inf,
+                   git_hash=compute_git_hash(), custom_quantities...)
     mkpath(output_directory)
 
     # Skip empty systems
     if nparticles(system_) == 0
         return
     end
 
+    v_ode_, u_ode_ = vu_ode_.x
+
     # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
     v_ode, u_ode, system, semi = transfer2cpu(v_ode_, u_ode_, system_, semi_)
 
@@ -136,10 +150,16 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc
         end
 
         # Extract custom quantities for this system
-        for (key, quantity) in custom_quantities
-            value = custom_quantity(quantity, system, v_ode, u_ode, semi, t)
-            if value !== nothing
-                vtk[string(key)] = value
+        if !isempty(custom_quantities)
+            dv_ode, du_ode = dvdu_ode_.x
+            dv_ode, du_ode = transfer2cpu(dv_ode, du_ode)
+
+            for (key, quantity) in custom_quantities
+                value = custom_quantity(quantity, system, dv_ode, du_ode, v_ode, u_ode,
+                                        semi, t)
+                if value !== nothing
+                    vtk[string(key)] = value
+                end
             end
         end
 
@@ -150,33 +170,45 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc
 end
 
 function transfer2cpu(v_::AbstractGPUArray, u_, system_, semi_)
-    v = Adapt.adapt(Array, v_)
-    u = Adapt.adapt(Array, u_)
     semi = Adapt.adapt(Array, semi_)
     system_index = system_indices(system_, semi_)
     system = semi.systems[system_index]
 
+    v, u = transfer2cpu(v_, u_)
+
     return v, u, system, semi
 end
 
 function transfer2cpu(v_, u_, system_, semi_)
     return v_, u_, system_, semi_
 end
 
-function custom_quantity(quantity::AbstractArray, system, v_ode, u_ode, semi, t)
+function transfer2cpu(v_::AbstractGPUArray, u_)
+    v = Adapt.adapt(Array, v_)
+    u = Adapt.adapt(Array, u_)
+
+    return v, u
+end
+
+function transfer2cpu(v_, u_)
+    return v_, u_
+end
+
+function custom_quantity(quantity::AbstractArray, system, dv_ode, du_ode, v_ode, u_ode,
+                         semi, t)
     return quantity
 end
 
-function custom_quantity(quantity, system, v_ode, u_ode, semi, t)
+function custom_quantity(quantity, system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     # Check if `quantity` is a function of `system`, `v_ode`, `u_ode`, `semi` and `t`
     if !isempty(methods(quantity,
-                        (typeof(system), typeof(v_ode), typeof(u_ode),
-                         typeof(semi), typeof(t))))
-        return quantity(system, v_ode, u_ode, semi, t)
+                        (typeof(system), typeof(dv_ode), typeof(du_ode), typeof(v_ode),
+                         typeof(u_ode), typeof(semi), typeof(t))))
+        return quantity(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     end
 
     # Assume `quantity` is a function of `data`
-    data = system_data(system, v_ode, u_ode, semi)
+    data = system_data(system, dv_ode, du_ode, v_ode, u_ode, semi)
     return quantity(system, data, t)
 end
 

src/schemes/boundary/open_boundary/system.jl

@@ -484,7 +484,7 @@ end
 # When the neighbor is an open boundary system, just use the viscosity of the fluid `system` instead
 @inline viscosity_model(system, neighbor_system::OpenBoundarySPHSystem) = system.viscosity
 
-function system_data(system::OpenBoundarySPHSystem, v_ode, u_ode, semi)
+function system_data(system::OpenBoundarySPHSystem, dv_ode, du_ode, v_ode, u_ode, semi)
     v = wrap_v(v_ode, system, semi)
     u = wrap_u(u_ode, system, semi)
 

src/schemes/boundary/system.jl

@@ -453,23 +453,24 @@ function system_correction(system::BoundarySPHSystem{<:BoundaryModelDummyParticl
     return system.boundary_model.correction
 end
 
-function system_data(system::BoundarySPHSystem, v_ode, u_ode, semi)
+function system_data(system::BoundarySPHSystem, dv_ode, du_ode, v_ode, u_ode, semi)
+    dv = [current_acceleration(system, particle) for particle in eachparticle(system)]
     v = wrap_v(v_ode, system, semi)
     u = wrap_u(u_ode, system, semi)
 
     coordinates = current_coordinates(u, system)
-    velocity = current_velocity(v, system)
+    velocity = [current_velocity(v, system, particle) for particle in eachparticle(system)]
     density = current_density(v, system)
     pressure = current_pressure(v, system)
 
-    return (; coordinates, velocity, density, pressure)
+    return (; coordinates, velocity, density, pressure, acceleration=dv)
 end
 
 function available_data(::BoundarySPHSystem)
-    return (:coordinates, :velocity, :density, :pressure)
+    return (:coordinates, :velocity, :density, :pressure, :acceleration)
 end
 
-function system_data(system::BoundaryDEMSystem, v_ode, u_ode, semi)
+function system_data(system::BoundaryDEMSystem, dv_ode, du_ode, v_ode, u_ode, semi)
     (; coordinates, radius, normal_stiffness) = system
 
     return (; coordinates, radius, normal_stiffness)