Merge branch 'edac-tvf' into permeable-BC-lastiwka

examples/fluid/taylor_green_vortex_2d.jl

@@ -9,10 +9,12 @@ using TrixiParticles
 using OrdinaryDiffEq
 
 # ==========================================================================================
-# ==== Fluid
-
+# ==== Resolution
 particle_spacing = 0.02
 
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 5.0)
 reynolds_number = 100.0
 
 box_length = 1.0
@@ -23,13 +25,6 @@ U = 1.0 # m/s
 
 b = -8pi^2 / reynolds_number
 
-fluid_density = 1.0
-sound_speed = 10U
-
-nu = U * box_length / reynolds_number
-
-background_pressure = sound_speed^2 * fluid_density
-
 p(pos) = -U^2 * exp(2 * b * 0) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4
 p(pos, t) = -U^2 * exp(2 * b * t) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4
 
@@ -38,6 +33,15 @@ v_y(pos) = U * exp(b * 0) * sin(2pi * pos[1]) * cos(2pi * pos[2])
 v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2])
 v_y(pos, t) = U * exp(b * t) * sin(2pi * pos[1]) * cos(2pi * pos[2])
 
+# ==========================================================================================
+# ==== Fluid
+fluid_density = 1.0
+sound_speed = 10U
+
+nu = U * box_length / reynolds_number
+
+background_pressure = sound_speed^2 * fluid_density
+
 smoothing_length = 1.0 * particle_spacing
 smoothing_kernel = SchoenbergQuinticSplineKernel{2}()
 
@@ -46,21 +50,24 @@ fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.
                          init_velocity=(0.0, 0.0))
 
 # Add small random displacement to the particles to avoid stagnant streamlines.
-seed!(42);
-fluid.coordinates .+= rand((-particle_spacing / 5):1e-5:(particle_spacing / 5),
-                           size(fluid.coordinates))
+#seed!(42);
+#fluid.coordinates .+= rand((-particle_spacing / 5):1e-5:(particle_spacing / 5),
+#                           size(fluid.coordinates))
 
 fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
                                            sound_speed, initial_pressure_function=p,
                                            initial_velocity_function=(v_x, v_y),
                                            transport_velocity=TransportVelocityAdami(background_pressure),
-                                           viscosity=ViscosityAdami(nu))
+                                           viscosity=ViscosityAdami(;nu))
+
+# ==========================================================================================
+# ==== Simulation
+
 semi = Semidiscretization(fluid_system,
                           neighborhood_search=GridNeighborhoodSearch,
                           periodic_box_min_corner=[0.0, 0.0],
                           periodic_box_max_corner=[box_length, box_length])
 
-tspan = (0.0, 5.0)
 ode = semidiscretize(semi, tspan)
 
 dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu))

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -149,6 +149,3 @@ end
 end
 
 @inline average_pressure(system, particle) = 0.0
-
-@inline viscosity_function(system) = system.viscosity
-@inline viscosity_function(system::BoundarySPHSystem) = system.boundary_model.viscosity

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -201,17 +201,14 @@ end
 function update_quantities!(system::EntropicallyDampedSPHSystem, v, u,
                             v_ode, u_ode, semi, t)
     summation_density!(system, semi, u, u_ode, system.density)
-    update_average_pressure!(system, system.transport_velocity, system_index,
-                             v_ode, u_ode, semi)
+    update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi)
 end
 
 function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi)
     return system
 end
 
-function update_average_pressure!(system, ::TransportVelocityAdami, system_index,
-                                  v_ode, u_ode, semi)
-    (; systems, neighborhood_searches) = semi
+function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode, semi)
     (; cache) = system
     (; pressure_average, neighbor_counter) = cache
 
@@ -221,15 +218,14 @@ function update_average_pressure!(system, ::TransportVelocityAdami, system_index
     u = wrap_u(u_ode, system, semi)
 
     # Use all other systems for the average pressure
-    @trixi_timeit timer() "compute average pressure" foreach_enumerate(systems) do (neighbor_system_index,
-                                                                                    neighbor_system)
-        u_neighbor_system = wrap_u(u_ode, neighbor_neighbor_system, semi)
-        v_neighbor_system = wrap_v(v_ode, neighbor_neighbor_system, semi)
+    @trixi_timeit timer() "compute average pressure" foreach_system(semi) do neighbor_system
+        u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
+        v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)
 
         system_coords = current_coordinates(u, system)
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
-        neighborhood_search = neighborhood_searches[system_index][neighbor_system_index]
+        neighborhood_search = neighborhood_searches(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff.
         for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords,

src/schemes/fluid/transport_velocity.jl

@@ -5,25 +5,27 @@ struct TransportVelocityAdami{ELTYPE}
     end
 end
 
-@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi)
-    return system
-end
+function update_transport_velocity!(u, ode, semi, t)
+    v_ode, u_ode = ode.u.x
+
+    foreach_system(semi) do system
+        update_transport_velocity!(system, v_ode, semi)
+    end
 
-@inline function update_transport_velocity!(system::FluidSystem, system_index,
-                                            v_ode, u_ode, semi)
-    update_transport_velocity!(system, system_index, v_ode, u_ode, semi,
-                               system.transport_velocity)
+    return ode
 end
 
-@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi,
-                                            ::Nothing)
-    return system
+@inline update_transport_velocity!(system, v_ode, semi) = system
+
+@inline function update_transport_velocity!(system::FluidSystem, v_ode, semi)
+    update_transport_velocity!(system, v_ode, semi, system.transport_velocity)
 end
 
-@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi,
-                                            ::TransportVelocityAdami)
+@inline update_transport_velocity!(system, v_ode, semi, ::Nothing) = system
+
+@inline function update_transport_velocity!(system, v_ode, semi, ::TransportVelocityAdami)
     (; cache) = system
-    v = wrap_v(v_ode, system_index, system, semi)
+    v = wrap_v(v_ode, system, semi)
 
     for particle in each_moving_particle(system)
         for i in 1:ndims(system)