Implement InitialCondition and store that in the systems (#159)

* Add `InitialCondition` * Use `InitialCondition` in `Semidiscretization` * Revert "Use `InitialCondition` in `Semidiscretization`" This reverts commit 479719b. * Store initial condition in systems * Adapt shapes * Remove `MergeShapes` * Start adapting example files * Fix undefined values * Only set current and not initial coordinates of solids in `restart_with` * Fix examples * Fix tests * Reformat code * Implement suggestions * Reformat code * Move smoothing kernel check before acceleration check The default acceleration depends on the smoothing kernel dimensionality, so passing a wrong smoothing kernel yielded "acceleration must be of length", which was confusing.
trixi-framework · May 26, 2023 · e19b806 · e19b806
1 parent e9e05c5
commit e19b806
Show file tree

Hide file tree

Showing 34 changed files with 523 additions and 586 deletions.
diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -36,39 +36,39 @@ state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0,
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
 
-setup = RectangularTank(particle_spacing, (water_width, water_height),
-                        (tank_width, tank_height), water_density,
-                        n_layers=boundary_layers, spacing_ratio=beta)
+tank = RectangularTank(particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
 
 # Move right boundary
 # Recompute the new water column width since the width has been rounded in `RectangularTank`.
-new_wall_position = (setup.n_particles_per_dimension[1] + 1) * particle_spacing
+new_wall_position = (tank.n_particles_per_dimension[1] + 1) * particle_spacing
 reset_faces = (false, true, false, false)
 positions = (0, new_wall_position, 0, 0)
 
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # ==========================================================================================
 # ==== Boundary models
 
-boundary_model = BoundaryModelDummyParticles(setup.boundary_densities,
-                                             setup.boundary_masses, state_equation,
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
                                              AdamiPressureExtrapolation(), smoothing_kernel,
                                              smoothing_length)
 
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,
-#                                              setup.boundary_masses)
+#                                              tank.boundary.mass)
 
 # ==========================================================================================
 # ==== Systems
 
-fluid_system = WeaklyCompressibleSPHSystem(setup, ContinuityDensity(), state_equation,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(setup.boundary_coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
@@ -100,7 +100,7 @@ sol = solve(ode, RDPK3SpFSAL49(),
 
 # Move right boundary
 positions = (0, tank_width, 0, 0)
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # Run full simulation
 tspan = (0.0, 5.7 / sqrt(9.81))

diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl
@@ -31,37 +31,37 @@ state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0,
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
 
-setup = RectangularTank(particle_spacing, (water_width, water_height, water_length),
-                        (tank_width, tank_height, tank_length), water_density,
-                        n_layers=boundary_layers, spacing_ratio=beta)
+tank = RectangularTank(particle_spacing, (water_width, water_height, water_length),
+                       (tank_width, tank_height, tank_length), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
 
 # Move right boundary
-new_wall_position = (setup.n_particles_per_dimension[1] + 1) * particle_spacing
+new_wall_position = (tank.n_particles_per_dimension[1] + 1) * particle_spacing
 reset_faces = (false, true, false, false, false, false)
 positions = (0, new_wall_position, 0, 0, 0, 0)
 
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # ==========================================================================================
 # ==== Boundary models
 
-boundary_model = BoundaryModelDummyParticles(setup.boundary_densities,
-                                             setup.boundary_masses, state_equation,
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
                                              AdamiPressureExtrapolation(), smoothing_kernel,
                                              smoothing_length)
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,
-#                                              setup.boundary_masses)
+#                                              tank.boundary.mass)
 
 # ==========================================================================================
 # ==== Systems
 
-fluid_system = WeaklyCompressibleSPHSystem(setup, ContinuityDensity(), state_equation,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity, 0.0))
 
-boundary_system = BoundarySPHSystem(setup.boundary_coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
@@ -91,7 +91,7 @@ sol = solve(ode, RDPK3SpFSAL49(),
 
 # Move right boundary
 positions = (0, tank_width, 0, 0)
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # Run full simulation
 tspan = (0.0, 5.7 / sqrt(9.81))

diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl
@@ -28,36 +28,36 @@ smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
 
-setup = RectangularTank(particle_spacing, (water_width, water_height),
-                        (tank_width, tank_height), water_density,
-                        n_layers=boundary_layers, spacing_ratio=beta)
+tank = RectangularTank(particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
 
 # Move water column
-for i in axes(setup.coordinates, 2)
-    setup.coordinates[:, i] .+= [0.5 * tank_width - 0.5 * water_width, 0.2]
+for i in axes(tank.fluid.coordinates, 2)
+    tank.fluid.coordinates[:, i] .+= [0.5 * tank_width - 0.5 * water_width, 0.2]
 end
 
 # ==========================================================================================
 # ==== Boundary models
 
-boundary_model = BoundaryModelDummyParticles(setup.boundary_densities,
-                                             setup.boundary_masses, state_equation,
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
                                              AdamiPressureExtrapolation(), smoothing_kernel,
                                              smoothing_length)
 
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,
-#                                              setup.boundary_masses)
+#                                              tank.boundary.mass)
 
 # ==========================================================================================
 # ==== Systems
 
-fluid_system = WeaklyCompressibleSPHSystem(setup, ContinuityDensity(), state_equation,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(setup.boundary_coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl
@@ -28,32 +28,32 @@ smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
 
-setup = RectangularTank(particle_spacing, (water_width, water_height),
-                        (tank_width, tank_height), water_density,
-                        n_layers=boundary_layers, spacing_ratio=beta)
+tank = RectangularTank(particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
 
 # ==========================================================================================
 # ==== Boundary models
 
-boundary_model = BoundaryModelDummyParticles(setup.boundary_densities,
-                                             setup.boundary_masses, state_equation,
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
                                              AdamiPressureExtrapolation(),
                                              smoothing_kernel,
                                              smoothing_length)
 
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,
-#                                              setup.boundary_masses)
+#                                              tank.boundary.mass)
 
 # ==========================================================================================
 # ==== Systems
 
-fluid_system = WeaklyCompressibleSPHSystem(setup, ContinuityDensity(), state_equation,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(setup.boundary_coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

diff --git a/examples/fsi/bending_beam_2d.jl b/examples/fsi/bending_beam_2d.jl
@@ -55,7 +55,7 @@ n_boundary_layers = 5
 # Since coordinates start in negative coordinate directions,
 # shift x and y coordinates outwards and downwards respectively.
 x_start_clamp_A = -(n_boundary_layers - 1) * solid_particle_spacing
-x_start_clamp_B = (beam.n_particles_per_dimension[1] + 1) * solid_particle_spacing
+x_start_clamp_B = (n_particles_per_dimension[1] + 1) * solid_particle_spacing
 
 n_particles_per_dimension_clamp = (n_boundary_layers, 10 * n_particles_y)
 n_fixed_particles = prod(n_particles_per_dimension_clamp) * 2
@@ -69,7 +69,7 @@ clamp_A = RectangularShape(solid_particle_spacing, n_particles_per_dimension_cla
 clamp_B = RectangularShape(solid_particle_spacing, n_particles_per_dimension_clamp,
                            (x_start_clamp_B, y_start_clamp), solid_density)
 
-solid = MergeShapes(beam, clamp_A, clamp_B)
+solid = InitialCondition(beam, clamp_A, clamp_B)
 
 # ==========================================================================================
 # ==== Fluid
@@ -100,7 +100,7 @@ fluid = RectangularShape(fluid_particle_spacing, n_particles_per_dimension,
 # ==========================================================================================
 # ==== Boundary models
 
-hydrodynamic_densites = water_density * ones(size(solid.densities))
+hydrodynamic_densites = water_density * ones(size(solid.density))
 hydrodynamic_masses = hydrodynamic_densites * solid_particle_spacing^2
 
 boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites,
@@ -116,20 +116,16 @@ boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites,
 # ==========================================================================================
 # ==== Systems
 
-solid_system = TotalLagrangianSPHSystem(solid.coordinates, solid.velocities,
-                                        solid.masses, solid.densities,
+solid_system = TotalLagrangianSPHSystem(solid,
                                         solid_smoothing_kernel, solid_smoothing_length,
                                         E, nu, boundary_model_solid,
                                         n_fixed_particles=n_fixed_particles,
                                         acceleration=(0.0, 0.0),
                                         penalty_force=PenaltyForceGanzenmueller(alpha=0.1))
 
-fluid_system = WeaklyCompressibleSPHSystem(fluid.coordinates,
-                                           zeros(size(fluid.coordinates)),
-                                           fluid.masses, fluid.densities,
+fluid_system = WeaklyCompressibleSPHSystem(fluid,
                                            ContinuityDensity(), state_equation,
-                                           fluid_smoothing_kernel,
-                                           fluid_smoothing_length,
+                                           fluid_smoothing_kernel, fluid_smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 

diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl
@@ -36,17 +36,17 @@ state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0,
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
 
-setup = RectangularTank(fluid_particle_spacing, (water_width, water_height),
-                        (tank_width, tank_height), water_density,
-                        n_layers=boundary_layers, spacing_ratio=beta)
+tank = RectangularTank(fluid_particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
 
 # Move right boundary
 # Recompute the new water column width since the width has been rounded in `RectangularTank`.
-new_wall_position = (setup.n_particles_per_dimension[1] + 1) * fluid_particle_spacing
+new_wall_position = (tank.n_particles_per_dimension[1] + 1) * fluid_particle_spacing
 reset_faces = (false, true, false, false)
 positions = (0, new_wall_position, 0, 0)
 
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # ==========================================================================================
 # ==== Solid
@@ -78,22 +78,22 @@ fixed_particles = RectangularShape(solid_particle_spacing,
                                    (n_particles_per_dimension[1], 1), (0.292, 0.0),
                                    solid_density)
 
-solid = MergeShapes(plate, fixed_particles)
+solid = InitialCondition(plate, fixed_particles)
 
 # ==========================================================================================
 # ==== Boundary models
 
-boundary_model = BoundaryModelDummyParticles(setup.boundary_densities,
-                                             setup.boundary_masses, state_equation,
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
                                              AdamiPressureExtrapolation(), smoothing_kernel,
                                              smoothing_length)
 
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,
-#                                              setup.boundary_masses)
+#                                              tank.boundary.mass)
 
 # For the FSI we need the hydrodynamic masses and densities in the solid boundary model
-hydrodynamic_densites = water_density * ones(size(solid.densities))
+hydrodynamic_densites = water_density * ones(size(solid.density))
 hydrodynamic_masses = hydrodynamic_densites * solid_particle_spacing^2
 
 solid_boundary_model = BoundaryModelDummyParticles(hydrodynamic_densites,
@@ -111,15 +111,14 @@ solid_boundary_model = BoundaryModelDummyParticles(hydrodynamic_densites,
 # ==========================================================================================
 # ==== Systems
 
-fluid_system = WeaklyCompressibleSPHSystem(setup, ContinuityDensity(), state_equation,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(setup.boundary_coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
 
-solid_system = TotalLagrangianSPHSystem(solid.coordinates, solid.velocities,
-                                        solid.masses, solid.densities,
+solid_system = TotalLagrangianSPHSystem(solid,
                                         solid_smoothing_kernel, solid_smoothing_length,
                                         E, nu,
                                         n_fixed_particles=n_particles_x,
@@ -156,7 +155,7 @@ sol = solve(ode, RDPK3SpFSAL49(),
 
 # Move right boundary
 positions = (0, tank_width, 0, 0)
-reset_wall!(setup, reset_faces, positions)
+reset_wall!(tank, reset_faces, positions)
 
 # Run full simulation
 tspan = (0.0, 1.0)