Initialize rectangular tank with hydrostatic pressure

src/setups/rectangular_shape.jl

@@ -32,8 +32,8 @@ rectangular = RectangularShape(particle_spacing, (5, 4, 7), (1.0, 2.0, 3.0), 100
 """
 function RectangularShape(particle_spacing, n_particles_per_dimension,
                           min_coordinates, density; pressure=0.0, tlsph=false,
-                          init_velocity=ntuple(_ -> 0.0,
-                                               length(n_particles_per_dimension)),
+                          init_velocity=ntuple(_ -> 0.0, length(n_particles_per_dimension)),
+                          acceleration=nothing, state_equation=nothing,
                           loop_order=:x_first)
     if particle_spacing < eps()
         throw(ArgumentError("`particle_spacing` needs to be positive and larger than $(eps())"))
@@ -58,8 +58,24 @@ function RectangularShape(particle_spacing, n_particles_per_dimension,
                                            loop_order=loop_order)
     velocities = init_velocity .* ones(ELTYPE, size(coordinates))
 
-    densities = density * ones(ELTYPE, n_particles)
-    masses = density * particle_spacing^NDIMS * ones(ELTYPE, n_particles)
+    if acceleration === nothing && state_equation === nothing
+        densities = density * ones(ELTYPE, n_particles)
+    elseif (acceleration isa AbstractVector || acceleration isa Tuple) &&
+           state_equation !== nothing
+        # Initialize hydrostatic pressure
+        pressure = Vector{ELTYPE}(undef, n_particles)
+        initialize_pressure!(pressure, particle_spacing, SVector(acceleration),
+                             state_equation, n_particles_per_dimension, loop_order)
+
+        # Get density from inverse state equation
+        densities = inverse_state_equation.(Ref(state_equation), pressure)
+    else
+        throw(ArgumentError("`acceleration` and `state_equation` must either be " *
+                            "used both or both must be `nothing`"))
+    end
+
+    particle_volume = particle_spacing^NDIMS
+    masses = particle_volume * densities
 
     return InitialCondition(coordinates, velocities, masses, densities, pressure=pressure,
                             particle_spacing=particle_spacing)
@@ -157,3 +173,81 @@ end
     coordinates[2, particle] = min_coordinates[2] + (y - 0.5) * particle_spacing
     coordinates[3, particle] = min_coordinates[3] + (z - 0.5) * particle_spacing
 end
+
+# 2D
+function initialize_pressure!(pressure, particle_spacing, acceleration, state_equation,
+                              n_particles_per_dimension::NTuple{2}, loop_order)
+    n_particles_x, n_particles_y = n_particles_per_dimension
+
+    @inline density(pressure_) = inverse_state_equation(state_equation, pressure_)
+
+    if loop_order !== :x_first
+        throw(ArgumentError("hydrostatic pressure calculation is only supported with loop" *
+                            "order `:x_first`"))
+    end
+
+    # Start with the highest index and loop backwards in order to start at the fluid surface
+    particle = prod(n_particles_per_dimension)
+
+    # The hydrostatic pressure is given by the ODE `dp/dr = rho * g`, where `r` is the
+    # distance to the surface (water depth), `rho` is the density and `g` is the
+    # acceleration.
+    # For high water columns (and especially for higher compressibility), this yields
+    # better results than just assuming `rho` to be constant.
+    #
+    # To solve this ODE, we use the explicit Euler method in each dimension
+    # independently, matching the particle indexing.
+    # This allows diagonal accelerations as well (albeit only on negative coordinate
+    # directions).
+    pressure_x = 0.0
+    pressure_x -= 0.5particle_spacing * acceleration[1] * density(pressure_x)
+    for x in n_particles_x:-1:1
+        # For the integration in y-direction, start at `pressure_x`
+        pressure_y = pressure_x
+        pressure_y -= 0.5particle_spacing * acceleration[2] * density(pressure_y)
+        for y in n_particles_y:-1:1
+            pressure[particle] = pressure_y
+            particle -= 1
+
+            # Explicit Euler step
+            pressure_y -= particle_spacing * acceleration[2] * density(pressure_y)
+        end
+
+        # Explicit Euler step
+        pressure_x -= particle_spacing * acceleration[1] * density(pressure_x)
+    end
+end
+
+# 3D
+function initialize_pressure!(pressure, particle_spacing, acceleration, state_equation,
+                              n_particles_per_dimension::NTuple{3}, loop_order)
+    n_particles_x, n_particles_y, n_particles_z = n_particles_per_dimension
+
+    @inline density(pressure_) = inverse_state_equation(state_equation, pressure_)
+
+    if loop_order !== :x_first
+        throw(ArgumentError("hydrostatic pressure calculation is only supported with loop" *
+                            "order `:x_first`"))
+    end
+
+    # Start with the highest index and loop backwards in order to start at the fluid surface
+    particle = prod(n_particles_per_dimension)
+
+    pressure_x = 0.0
+    for x in n_particles_x:-1:1
+        pressure_y = pressure_x
+        for y in n_particles_y:-1:1
+            pressure_z = pressure_y
+            for z in n_particles_z:-1:1
+                pressure[particle] = pressure_z
+                particle -= 1
+
+                pressure_z -= particle_spacing * acceleration[3] * density(pressure_z)
+            end
+
+            pressure_y -= particle_spacing * acceleration[2] * density(pressure_y)
+        end
+
+        pressure_x -= particle_spacing * acceleration[1] * density(pressure_x)
+    end
+end

src/setups/rectangular_tank.jl

@@ -58,12 +58,12 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real}
     n_layers                  :: Int
     n_particles_per_dimension :: NTuple{NDIMS, Int}
 
-    function RectangularTank(particle_spacing, fluid_size, tank_size,
-                             fluid_density; pressure=0.0,
-                             n_layers=1, spacing_ratio=1.0,
+    function RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density;
+                             pressure=0.0, n_layers=1, spacing_ratio=1.0,
                              init_velocity=zeros(length(fluid_size)),
                              boundary_density=fluid_density,
-                             faces=Tuple(trues(2 * length(fluid_size))))
+                             faces=Tuple(trues(2 * length(fluid_size))),
+                             acceleration=nothing, state_equation=nothing)
         NDIMS = length(fluid_size)
         ELTYPE = eltype(particle_spacing)
         fluid_size_ = Tuple(ELTYPE.(fluid_size))
@@ -114,7 +114,8 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real}
 
         fluid = RectangularShape(particle_spacing, n_particles_per_dim, zeros(NDIMS),
                                  fluid_density, init_velocity=init_velocity,
-                                 pressure=pressure)
+                                 pressure=pressure,
+                                 acceleration=acceleration, state_equation=state_equation)
 
         boundary = InitialCondition(boundary_coordinates, boundary_velocities,
                                     boundary_masses, boundary_densities,