adapt docs

Author: LasNikas

examples/fluid/play_ground.jl

@@ -0,0 +1,233 @@
+# Flow past a circular cylinder (vortex street), Tafuni et al. (2018).
+# Other literature using this validation:
+# Vacandio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998)
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+factor_d = 0.08 # Resolution in the paper is `0.01` (5M particles)
+
+const cylinder_diameter = 0.1
+particle_spacing = factor_d * cylinder_diameter
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 4
+
+# Make sure that the kernel support of fluid particles at an open boundary is always
+# fully sampled.
+# Note: Due to the dynamics at the inlets and outlets of open boundaries,
+# it is recommended to use `open_boundary_layers > boundary_layers`
+open_boundary_layers = 8
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 1.5)
+
+# Boundary geometry and initial fluid particle positions
+domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter)
+
+flow_direction = [1.0, 0.0]
+reynolds_number = 200
+const prescribed_velocity = 1.0
+const fluid_density = 1000.0
+sound_speed = 10 * prescribed_velocity
+
+boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+                 domain_size[2])
+
+pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
+                       n_layers=boundary_layers, velocity=[prescribed_velocity, 0.0],
+                       faces=(false, false, true, true))
+
+# Shift pipe walls in negative x-direction for the inflow
+pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
+
+n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1)
+
+cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2)
+cylinder = SphereShape(particle_spacing, cylinder_diameter / 2,
+                       cylinder_center, fluid_density, sphere_type=RoundSphere())
+
+fluid = setdiff(pipe.fluid, cylinder)
+
+# ==========================================================================================
+# ==== Fluid
+wcsph = true
+
+h_factor = 1.5
+smoothing_length = h_factor * particle_spacing
+smoothing_kernel = WendlandC2Kernel{2}()
+
+fluid_density_calculator = ContinuityDensity()
+
+kinematic_viscosity = prescribed_velocity * cylinder_diameter / reynolds_number
+
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7)
+viscosity = ViscosityAdami(nu=kinematic_viscosity)
+
+if wcsph
+    density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
+
+    fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                               state_equation, smoothing_kernel,
+                                               density_diffusion=density_diffusion,
+                                               smoothing_length, viscosity=viscosity,
+                                               pressure_acceleration=tensile_instability_control,
+                                               buffer_size=n_buffer_particles)
+else
+    state_equation = nothing
+    fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
+                                               sound_speed, viscosity=viscosity,
+                                               density_calculator=fluid_density_calculator,
+                                               buffer_size=n_buffer_particles)
+end
+
+# ==========================================================================================
+# ==== Open Boundary
+function velocity_function2d(pos, t)
+    return SVector(prescribed_velocity, 0.0)
+end
+
+open_boundary_model = BoundaryModelTafuni()
+
+boundary_type_in = InFlow()
+plane_in = ([0.0, 0.0], [0.0, domain_size[2]])
+inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers,
+                      density=fluid_density, particle_spacing,
+                      boundary_type=boundary_type_in)
+
+reference_velocity_in = velocity_function2d
+# At the inlet, neither pressure nor density are prescribed; instead,
+# these values are extrapolated from the fluid domain
+reference_pressure_in = nothing
+reference_density_in = nothing
+open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
+                                         boundary_model=open_boundary_model,
+                                         buffer_size=n_buffer_particles,
+                                         reference_density=reference_density_in,
+                                         reference_pressure=reference_pressure_in,
+                                         reference_velocity=reference_velocity_in)
+
+boundary_type_out = OutFlow()
+
+plane_out = ([pipe.fluid_size[1], 0.0], [pipe.fluid_size[1], domain_size[2]])
+outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction),
+                       open_boundary_layers, density=fluid_density, particle_spacing,
+                       boundary_type=boundary_type_out)
+
+# At the outlet, we allow the flow to exit freely without imposing any boundary conditions
+reference_velocity_out = nothing
+reference_pressure_out = nothing
+reference_density_out = nothing
+open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
+                                          boundary_model=open_boundary_model,
+                                          buffer_size=n_buffer_particles,
+                                          reference_density=reference_density_out,
+                                          reference_pressure=reference_pressure_out,
+                                          reference_velocity=reference_velocity_out)
+# ==========================================================================================
+# ==== Boundary
+boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             state_equation=state_equation,
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system_wall = BoundarySPHSystem(pipe.boundary, boundary_model)
+
+boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass,
+                                                      AdamiPressureExtrapolation(),
+                                                      state_equation=state_equation,
+                                                      viscosity=viscosity,
+                                                      smoothing_kernel, smoothing_length)
+
+boundary_system_cylinder = BoundarySPHSystem(cylinder, boundary_model_cylinder)
+
+# ==========================================================================================
+# ==== Postprocessing
+circle = SphereShape(particle_spacing, (cylinder_diameter + particle_spacing) / 2,
+                     cylinder_center, fluid_density, n_layers=1,
+                     sphere_type=RoundSphere())
+
+# Points for pressure interpolation, located at the wall interface
+const data_points = copy(circle.coordinates)
+const center = SVector(cylinder_center)
+
+calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing
+function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t)
+    force = zero(SVector{ndims(system), eltype(system)})
+
+    values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false,
+                                clip_negative_pressure=false)
+    pressure = Array(values.pressure)
+
+    for i in axes(data_points, 2)
+        point = TrixiParticles.current_coords(data_points, system, i)
+
+        # F = ∑ -p_i * A_i * n_i
+        force -= pressure[i] * particle_spacing .*
+                 TrixiParticles.normalize(point - center)
+    end
+
+    return 2 * force[2] / (fluid_density * prescribed_velocity^2 * cylinder_diameter)
+end
+
+calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing
+function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t)
+    force = zero(SVector{ndims(system), eltype(system)})
+
+    values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false,
+                                clip_negative_pressure=false)
+    pressure = Array(values.pressure)
+
+    for i in axes(data_points, 2)
+        point = TrixiParticles.current_coords(data_points, system, i)
+
+        # F = ∑ -p_i * A_i * n_i
+        force -= pressure[i] * particle_spacing .*
+                 TrixiParticles.normalize(point - center)
+    end
+
+    return 2 * force[1] / (fluid_density * prescribed_velocity^2 * cylinder_diameter)
+end
+
+# ==========================================================================================
+# ==== Simulation
+min_corner = minimum(pipe.boundary.coordinates .- 5 * particle_spacing, dims=2)
+max_corner = maximum(pipe.boundary.coordinates .+ 5 * particle_spacing, dims=2)
+cell_list = FullGridCellList(; min_corner, max_corner)
+
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list,
+                                                update_strategy=ParallelUpdate())
+
+semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out,
+                          boundary_system_wall, boundary_system_cylinder;
+                          neighborhood_search=neighborhood_search,
+                          parallelization_backend=PolyesterBackend())
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+output_directory = "out"
+
+saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory)
+
+pp_callback = PostprocessCallback(; dt=0.02,
+                                  f_l=calculate_lift_force, f_d=calculate_drag_force,
+                                  output_directory, filename="resulting_force",
+                                  write_csv=true, write_file_interval=10)
+
+extra_callback = nothing
+
+callbacks = CallbackSet(info_callback, UpdateCallback(), saving_callback,
+                        ParticleShiftingCallback(), # To obtain a near-uniform particle distribution in the wake
+                        pp_callback, extra_callback)
+
+sol = solve(ode, RDPK3SpFSAL35(),
+            abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

src/schemes/boundary/open_boundary/method_of_characteristics.jl

@@ -9,8 +9,7 @@ For more information about the method see [description below](@ref method_of_cha
 
 # Keywords
 - `extrapolate_reference_values=false`: If `true`, the reference values are extrapolated
-  from the fluid domain to the boundary particles. This is useful for open boundaries where
-  the reference values are not known a priori.
+  from the fluid domain to the boundary particles.
   **Note:** This feature is experimental and has not been fully validated yet.
   As of now, we are not aware of any published literature supporting its use.
 """

src/schemes/boundary/open_boundary/mirroring.jl

@@ -150,7 +150,10 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary,
             density[particle] = reference_value(reference_density, density[particle],
                                                 particle_coords, t)
         else
-            inverse_state_equation!(density, state_equation, pressure, particle)
+            f_d = L_inv * extrapolated_density_correction[]
+            df_d = f_d[two_to_end]
+
+            density[particle] = f_d[1] + dot(pos_diff, df_d)
         end
     end
 

src/schemes/boundary/open_boundary/system.jl

@@ -27,10 +27,10 @@ Open boundary system for in- and outflow particles.
                        or a scalar for a constant density over all particles.
 
 !!! note "Note"
-    When using the [`BoundaryModelTafuni()`](@ref), the reference values (`reference_velocity`,
-    `reference_pressure`, `reference_density`) can also be set to `nothing`
-    since this model allows for either assigning physical quantities a priori or extrapolating them
-    from the fluid domaim to the buffer zones (inflow and outflow) using ghost nodes.
+    The reference values (`reference_velocity`, `reference_pressure`, `reference_density`)
+    can also be set to `nothing`.
+    In this case, they will either be extrapolated from the fluid domain ([BoundaryModelTafuni](@ref BoundaryModelTafuni))
+    or evolved using the characteristic flow variables ([BoundaryModelLastiwka](@ref BoundaryModelLastiwka)).
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in future releases.
@@ -83,9 +83,6 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
                                reference_density=nothing)
     (; initial_condition) = boundary_zone
 
-    check_reference_values!(boundary_model, reference_density, reference_pressure,
-                            reference_velocity)
-
     buffer = SystemBuffer(nparticles(initial_condition), buffer_size)
 
     initial_condition = allocate_buffer(initial_condition, buffer)
@@ -261,6 +258,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_
     u = wrap_u(u_ode, system, semi)
     v = wrap_v(v_ode, system, semi)
 
+    @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
+
     # Update density, pressure and velocity based on the characteristic variables.
     # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971
     @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system,
@@ -270,8 +269,6 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_
                                                                                    u_ode,
                                                                                    semi, t)
 
-    @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
-
     return system
 end
 
@@ -443,25 +440,6 @@ end
 # instead of using the method of characteristics
 reference_value(value::Nothing, quantity, position, t) = quantity
 
-function check_reference_values!(boundary_model, reference_density, reference_pressure,
-                                 reference_velocity)
-    return boundary_model
-end
-
-function check_reference_values!(boundary_model::BoundaryModelLastiwka,
-                                 reference_density, reference_pressure, reference_velocity)
-
-                                 # TODO
-                                 return boundary_model
-    boundary_model.extrapolate_reference_values && return boundary_model
-
-    if any(isnothing.([reference_density, reference_pressure, reference_velocity]))
-        throw(ArgumentError("for `BoundaryModelLastiwka` all reference values must be specified"))
-    end
-
-    return boundary_model
-end
-
 # To account for boundary effects in the viscosity term of the RHS, use the viscosity model
 # of the neighboring particle systems.
 @inline function viscosity_model(system::OpenBoundarySPHSystem,