Fix rectangular shape input handling

src/TrixiParticles.jl

@@ -21,7 +21,7 @@ using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using ReadVTK: ReadVTK
 using RecipesBase: RecipesBase, @series
-using Random: seed!
+using Random: MersenneTwister
 using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
                  get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!,
                  add_tstop!

src/setups/rectangular_shape.jl

@@ -92,7 +92,7 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord
         throw(ArgumentError("`min_coordinates` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
-    if density !== nothing && any(density .< eps())
+    if density !== nothing && !(density isa Function) && any(density .< eps())
         throw(ArgumentError("`density` needs to be positive and larger than $(eps())"))
     end
 
@@ -105,15 +105,19 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord
                                            place_on_shell, loop_order)
 
     if !isnothing(coordinates_perturbation)
-        seed!(1)
         amplitude = coordinates_perturbation * particle_spacing
-        coordinates .+= rand((-amplitude):(particle_spacing * 1e-3):(amplitude),
+        coordinates .+= rand(MersenneTwister(1),
+                             (-amplitude):(particle_spacing * 1e-3):(amplitude),
                              NDIMS, n_particles)
     end
 
     # Allow zero acceleration with state equation, but interpret `nothing` acceleration
     # with state equation as a likely mistake.
     if acceleration isa AbstractVector || acceleration isa Tuple
+        if length(acceleration) != NDIMS
+            throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
+        end
+
         if pressure != 0.0
             throw(ArgumentError("`pressure` cannot be used together with `acceleration` " *
                                 "and `state_equation` (hydrostatic pressure gradient)"))
@@ -233,6 +237,11 @@ function initialize_pressure!(pressure, particle_spacing, acceleration, density_
     # Dimension in which the acceleration is acting
     accel_dim = findfirst(a -> abs(a) > eps(), acceleration)
 
+    if accel_dim === nothing
+        fill!(pressure, zero(eltype(pressure)))
+        return pressure
+    end
+
     # Compute 1D pressure gradient with explicit Euler method
     factor = particle_spacing * abs(acceleration[accel_dim])
 

test/setups/rectangular_shape.jl

@@ -48,6 +48,33 @@
                 @test shape.coordinates == expected_coords[i]
             end
         end
+
+        @testset "Function Density" begin
+            shape = RectangularShape(0.1, (2, 1), (0.0, 0.0),
+                                     density=coords -> 1000.0 + coords[1])
+
+            @test shape.density ≈ [1000.05, 1000.15]
+        end
+
+        @testset "Coordinates Perturbation Does Not Reset Random State" begin
+            Random.seed!(42)
+            first_random_number = rand()
+            next_random_number = rand()
+
+            Random.seed!(42)
+            @test rand() == first_random_number
+
+            RectangularShape(0.1, (2, 2), (0.0, 0.0), density=1.0,
+                             coordinates_perturbation=0.1)
+
+            @test rand() == next_random_number
+        end
+
+        @testset "Errors" begin
+            @test_throws ArgumentError RectangularShape(0.1, (2, 2), (0.0, 0.0),
+                                                       density=1000.0,
+                                                       acceleration=(0.0, -9.81, 0.0))
+        end
     end
 
     # Only show all of these nested testsets in case of errors
@@ -123,6 +150,14 @@
                 @test shape.pressure ≈ 4.71 * 1000.0 * vec(reverse(pressure'))
             end
         end
+
+        @testset "Zero Acceleration" begin
+            shape = RectangularShape(particle_spacing, (2, 5), (0.0, 0.0),
+                                     density=1000.0, acceleration=(0.0, 0.0))
+
+            @test shape.pressure == zeros(10)
+            @test shape.density == 1000 * ones(10)
+        end
     end
 
     # Use `@trixi_testset` to isolate the mock functions in a separate namespace
@@ -186,6 +221,16 @@
                                                          shape.pressure)
             @test shape.mass == particle_spacing^2 * shape.density
         end
+
+        @testset "Zero Acceleration" begin
+            shape = RectangularShape(particle_spacing, (2, 5), (0.0, 0.0);
+                                     acceleration=(0.0, 0.0), state_equation)
+
+            @test shape.pressure == zeros(10)
+            @test shape.density ==
+                  TrixiParticles.inverse_state_equation.(Ref(state_equation),
+                                                         shape.pressure)
+        end
     end
 end