Use source terms instead of damping (#344)

* Use source terms instead of damping * Update docs of `Semidiscretization` * Update docs of systems * [skip ci] Add docs for source terms * Remove `neighborhood_search=GridNeighborhoodSearch` * Fix tests * Add check for positive density to `InitialCondition` * Add tests for source terms * Fix tests * Fix examples * Make n-body fast again * Implement suggestions * Add source terms after interactions * Implement suggestions * Reformat
trixi-framework · Jan 29, 2024 · e34b76c · e34b76c
1 parent 52740d7
commit e34b76c
Show file tree

Hide file tree

Showing 28 changed files with 221 additions and 124 deletions.
diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl
@@ -62,8 +62,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=10)

diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -65,8 +65,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=250)

diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl
@@ -56,8 +56,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=10)

diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl
@@ -56,8 +56,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/fluid/moving_wall_2d.jl b/examples/fluid/moving_wall_2d.jl
@@ -77,8 +77,7 @@ boundary_system_wall = BoundarySPHSystem(wall, boundary_model_wall,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_wall,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_wall)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl
@@ -52,7 +52,6 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 # ==========================================================================================
 # ==== Simulation
 semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch,
                           periodic_box_min_corner=[0.0, -0.25],
                           periodic_box_max_corner=[1.0, 0.75])
 ode = semidiscretize(semi, tspan)

diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl
@@ -63,8 +63,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)

diff --git a/examples/fluid/rectangular_tank_edac_2d.jl b/examples/fluid/rectangular_tank_edac_2d.jl
@@ -49,8 +49,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)

diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl
@@ -122,8 +122,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system, solid_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system, solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl
@@ -151,8 +151,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
 # ==========================================================================================
 # ==== Simulation
 semi = Semidiscretization(fluid_system, boundary_system_tank,
-                          boundary_system_gate, solid_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+                          boundary_system_gate, solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl
@@ -109,8 +109,7 @@ solid_system_2 = TotalLagrangianSPHSystem(sphere2,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=10)

diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl
@@ -64,8 +64,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, solid_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(fluid_system, solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl
@@ -68,15 +68,15 @@ masses = [
 # Offset sun momentum
 velocity[:, 1] = -velocity[:, 2:end] * masses[2:end] / SOLAR_MASS
 
-initial_condition = InitialCondition(; coordinates, velocity, density=0.0, mass=masses)
+initial_condition = InitialCondition(; coordinates, velocity, density=1.0, mass=masses)
 
 G = 1.0
 particle_system = NBodySystem(initial_condition, G)
 
 # ==========================================================================================
 # ==== Simulation
 
-semi = Semidiscretization(particle_system)
+semi = Semidiscretization(particle_system, neighborhood_search=nothing)
 
 # This is significantly faster than using OrdinaryDiffEq.
 function symplectic_euler!(velocity, coordinates, semi)

diff --git a/examples/n_body/n_body_solar_system.jl b/examples/n_body/n_body_solar_system.jl
@@ -23,15 +23,15 @@ masses = [
     1.99e30, 3.30e23, 4.87e24, 5.97e24, 6.42e23, 1.90e27, 5.68e26, 8.68e25, 1.02e26,
 ]
 
-initial_condition = InitialCondition(; coordinates, velocity, density=0.0, mass=masses)
+initial_condition = InitialCondition(; coordinates, velocity, density=1.0, mass=masses)
 
 G = 6.6743e-11
 particle_system = NBodySystem(initial_condition, G)
 
 # ==========================================================================================
 # ==== Simulation
 
-semi = Semidiscretization(particle_system)
+semi = Semidiscretization(particle_system, neighborhood_search=nothing)
 
 day = 24 * 3600.0
 year = 365day

diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl
@@ -18,10 +18,6 @@ TrixiParticles.timer_name(::NBodySystem) = "nbody"
 
 @inline Base.eltype(system::NBodySystem) = eltype(system.initial_condition.coordinates)
 
-@inline function TrixiParticles.add_acceleration!(dv, particle, system::NBodySystem)
-    return dv
-end
-
 function TrixiParticles.write_u0!(u0, system::NBodySystem)
     u0 .= system.initial_condition.coordinates
 
@@ -53,11 +49,8 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system,
                                   neighbor_system::NBodySystem)
     (; mass, G) = neighbor_system
 
-    system_coords = TrixiParticles.current_coordinates(u_particle_system,
-                                                       particle_system)
-
-    neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system,
-                                                         neighbor_system)
+    system_coords = TrixiParticles.current_coordinates(u_particle_system, particle_system)
+    neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, neighbor_system)
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff.
     TrixiParticles.for_particle_neighbor(particle_system, neighbor_system,

diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl
@@ -58,7 +58,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(solid_system, neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -50,11 +50,12 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
        PressureMirroring, PressureZeroing
 export BoundaryMovement
-export GridNeighborhoodSearch
+export GridNeighborhoodSearch, TrivialNeighborhoodSearch
 export examples_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape
 export VoxelSphere, RoundSphere, reset_wall!
+export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles

diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
@@ -126,6 +126,10 @@ struct InitialCondition{ELTYPE}
             densities = density_fun.(coordinates_svector)
         end
 
+        if any(densities .< eps())
+            throw(ArgumentError("density must be positive and larger than `eps()`"))
+        end
+
         if pressure isa AbstractVector
             if length(pressure) != n_particles
                 throw(ArgumentError("Expected: length(pressure) == size(coordinates, 2)\n" *

diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
@@ -1,27 +1,42 @@
 """
-    Semidiscretization(systems...; neighborhood_search=nothing, damping_coefficient=nothing)
+    Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch,
+                       periodic_box_min_corner=nothing, periodic_box_max_corner=nothing)
 
 The semidiscretization couples the passed systems to one simulation.
 
 The type of neighborhood search to be used in the simulation can be specified with
 the keyword argument `neighborhood_search`. A value of `nothing` means no neighborhood search.
 
+# Arguments
+- `systems`: Systems to be coupled in this semidiscretization
+
+# Keywords
+- `neighborhood_search`:    The type of neighborhood search to be used in the simulation.
+                            By default, the [`GridNeighborhoodSearch`](@ref) is used.
+                            Use [`TrivialNeighborhoodSearch`](@ref) to loop over all particles
+                            (no neighborhood search).
+- `periodic_box_min_corner`:    In order to use a periodic domain, pass the coordinates
+                                of the domain corner in negative coordinate directions.
+- `periodic_box_max_corner`:    In order to use a periodic domain, pass the coordinates
+                                of the domain corner in positive coordinate directions.
+
 # Examples
 ```julia
-semi = Semidiscretization(fluid_system, boundary_system; neighborhood_search=GridNeighborhoodSearch, damping_coefficient=nothing)
+semi = Semidiscretization(fluid_system, boundary_system)
+
+semi = Semidiscretization(fluid_system, boundary_system,
+                          neighborhood_search=TrivialNeighborhoodSearch)
 ```
 """
-struct Semidiscretization{S, RU, RV, NS, DC}
+struct Semidiscretization{S, RU, RV, NS}
     systems               :: S
     ranges_u              :: RU
     ranges_v              :: RV
     neighborhood_searches :: NS
-    damping_coefficient   :: DC
 
-    function Semidiscretization(systems...; neighborhood_search=nothing,
+    function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch,
                                 periodic_box_min_corner=nothing,
-                                periodic_box_max_corner=nothing,
-                                damping_coefficient=nothing)
+                                periodic_box_max_corner=nothing)
         sizes_u = [u_nvariables(system) * n_moving_particles(system)
                    for system in systems]
         ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i])
@@ -44,9 +59,8 @@ struct Semidiscretization{S, RU, RV, NS, DC}
                                for neighbor in systems)
                          for system in systems)
 
-        new{typeof(systems), typeof(ranges_u), typeof(ranges_v),
-            typeof(searches), typeof(damping_coefficient)}(systems, ranges_u, ranges_v,
-                                                           searches, damping_coefficient)
+        new{typeof(systems), typeof(ranges_u),
+            typeof(ranges_v), typeof(searches)}(systems, ranges_u, ranges_v, searches)
     end
 end
 
@@ -75,7 +89,6 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization)
         summary_line(io, "#systems", length(semi.systems))
         summary_line(io, "neighborhood search",
                      semi.neighborhood_searches |> eltype |> eltype |> nameof)
-        summary_line(io, "damping coefficient", semi.damping_coefficient)
         summary_line(io, "total #particles", sum(nparticles.(semi.systems)))
         summary_footer(io)
     end
@@ -310,11 +323,10 @@ function kick!(dv_ode, v_ode, u_ode, semi, t)
         @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
                                                                               semi, t)
 
-        @trixi_timeit timer() "gravity and damping" gravity_and_damping!(dv_ode, v_ode,
-                                                                         semi)
-
         @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode,
                                                                        semi)
+
+        @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi)
     end
 
     return dv_ode
@@ -374,25 +386,28 @@ function update_nhs(u_ode, semi)
     end
 end
 
-function gravity_and_damping!(dv_ode, v_ode, semi)
-    (; damping_coefficient) = semi
-
-    # Set velocity and add acceleration for each system
+function add_source_terms!(dv_ode, v_ode, u_ode, semi)
     foreach_system(semi) do system
         dv = wrap_v(dv_ode, system, semi)
         v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
 
         @threaded for particle in each_moving_particle(system)
-            # This can be dispatched per system
+            # Dispatch by system type to exclude boundary systems
             add_acceleration!(dv, particle, system)
-            add_damping_force!(dv, damping_coefficient, v, particle, system)
+            add_source_terms_inner!(dv, v, u, particle, system, source_terms(system))
         end
     end
 
     return dv_ode
 end
 
-@inline function add_acceleration!(dv, particle, system)
+@inline source_terms(system) = nothing
+@inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms
+
+@inline add_acceleration!(dv, particle, system) = dv
+
+@inline function add_acceleration!(dv, particle, system::Union{FluidSystem, SolidSystem})
     (; acceleration) = system
 
     for i in 1:ndims(system)
@@ -402,20 +417,56 @@ end
     return dv
 end
 
-@inline add_acceleration!(dv, particle, system::BoundarySPHSystem) = dv
+@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_)
+    coords = current_coords(u, system, particle)
+    velocity = current_velocity(v, system, particle)
+    density = particle_density(v, system, particle)
+    pressure = particle_pressure(v, system, particle)
 
-@inline function add_damping_force!(dv, damping_coefficient, v, particle,
-                                    system::FluidSystem)
-    for i in 1:ndims(system)
-        dv[i, particle] -= damping_coefficient * v[i, particle]
+    source = source_terms_(coords, velocity, density, pressure)
+
+    # Loop over `eachindex(source)`, so that users could also pass source terms for
+    # the density when using `ContinuityDensity`.
+    for i in eachindex(source)
+        dv[i, particle] += source[i]
     end
 
     return dv
 end
 
-# Currently no damping for non-fluid systems
-@inline add_damping_force!(dv, damping_coefficient, v, particle, system) = dv
-@inline add_damping_force!(dv, ::Nothing, v, particle, system::FluidSystem) = dv
+@inline add_source_terms_inner!(dv, v, u, particle, system, source_terms_::Nothing) = dv
+
+@doc raw"""
+    SourceTermDamping(; damping_coefficient)
+
+A source term to be used when a damping step is required before running a full simulation.
+The term ``-c \cdot v_a`` is added to the acceleration ``\frac{\mathrm{d}v_a}{\mathrm{d}t}``
+of particle ``a``, where ``c`` is the damping coefficient and ``v_a`` is the velocity of
+particle ``a``.
+
+# Keywords
+- `damping_coefficient`:    The coefficient ``d`` above. A higher coefficient means more
+                            damping. A coefficient of `1e-4` is a good starting point for
+                            damping a fluid at rest.
+
+# Examples
+```julia
+source_terms = SourceTermDamping(; damping_coefficient=1e-4)
+```
+"""
+struct SourceTermDamping{ELTYPE}
+    damping_coefficient::ELTYPE
+
+    function SourceTermDamping(; damping_coefficient)
+        return new{typeof(damping_coefficient)}(damping_coefficient)
+    end
+end
+
+@inline function (source_term::SourceTermDamping)(coords, velocity, density, pressure)
+    (; damping_coefficient) = source_term
+
+    return -damping_coefficient * velocity
+end
 
 function system_interaction!(dv_ode, v_ode, u_ode, semi)
     # Call `interact!` for each pair of systems