Ode reoder

docs/src/general/semidiscretization.md

@@ -2,5 +2,6 @@
 
 ```@autodocs
 Modules = [TrixiParticles]
-Pages = [joinpath("general", "semidiscretization.jl")]
+Pages = [joinpath("general", "semidiscretization.jl"),
+         joinpath("general", "source_terms.jl")]
 ```

src/callbacks/split_integration.jl

@@ -214,46 +214,14 @@ end
 
 # Update the systems before calling `interact!` to compute forces.
 function update_systems_split!(semi, v_ode, u_ode, t)
-    # First update step before updating the NHS
-    # (for example for writing the current coordinates in the solid system)
-    foreach_system(semi) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-
-        update_positions!(system, v, u, v_ode, u_ode, semi, t)
-    end
-
-    # Second update step.
-    # This is used to calculate density and pressure of the fluid systems
-    # before updating the boundary systems,
-    # since the fluid pressure is needed by the Adami interpolation.
-    foreach_system(semi) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-
-        update_quantities!(system, v, u, v_ode, u_ode, semi, t)
-    end
-
-    # Perform correction and pressure calculation
-    foreach_system(semi) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-
-        update_pressure!(system, v, u, v_ode, u_ode, semi, t)
-    end
-
     # No `update_boundary_interpolation!` for performance reasons, or we will lose
     # a lot of the speedup that we can gain with split integration.
     # We assume that the TLSPH particles move so little during the substeps
     # that the extrapolated pressure/density values can be treated as constant.
-
-    # Final update step for all remaining systems
-    foreach_system(semi) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-
-        update_final!(system, v, u, v_ode, u_ode, semi, t)
-    end
+    update_systems!(v_ode, u_ode, semi, t;
+                    update_nhs=false,
+                    update_boundary_interpolation=false,
+                    update_inter_system=false)
 end
 
 function system_interaction_split!(dv_ode_split, v_ode, u_ode, semi,

src/general/abstract_system.jl

@@ -16,6 +16,39 @@ abstract type AbstractBoundarySystem{NDIMS} <: AbstractSystem{NDIMS} end
 timer_name(::AbstractBoundarySystem) = "boundary"
 vtkname(system::AbstractBoundarySystem) = "boundary"
 
+# Helpers for iterating system collections without allocations.
+@inline function system_indices(system, systems_or_semi)
+    # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict`
+    # is ~30x slower because `hash(::System)` is very slow.
+    systems = hasfield(typeof(systems_or_semi), :systems) ? systems_or_semi.systems :
+              systems_or_semi
+    index = findfirst(==(system), systems)
+
+    if isnothing(index)
+        throw(ArgumentError("system is not in the semidiscretization"))
+    end
+
+    return index
+end
+
+# This is just for readability to loop over all systems without allocations.
+@inline function foreach_system(f, systems_or_semi)
+    systems = hasfield(typeof(systems_or_semi), :systems) ? systems_or_semi.systems :
+              systems_or_semi
+    return foreach_noalloc(f, systems)
+end
+
+# This is just for readability to loop over all systems with indices without allocations.
+@inline function foreach_system_indexed(f, systems_or_semi)
+    systems = hasfield(typeof(systems_or_semi), :systems) ? systems_or_semi.systems :
+              systems_or_semi
+    indices = ntuple(identity, Val(length(systems)))
+
+    return foreach_noalloc(indices, systems) do (index, system)
+        f(index, system)
+    end
+end
+
 # Number of integrated variables in the first component of the ODE system (coordinates)
 @inline u_nvariables(system) = ndims(system)
 

src/general/general.jl

@@ -8,4 +8,5 @@ include("initial_condition.jl")
 include("buffer.jl")
 include("interpolation.jl")
 include("custom_quantities.jl")
+include("source_terms.jl")
 include("time_integration.jl")