Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@svchb
svchb committed Jul 18, 2024
1 parent c53246d commit 6bc5567
Show file tree
Hide file tree
Showing 7 changed files with 69 additions and 24 deletions.
14 changes: 7 additions & 7 deletions examples/fluid/falling_water_spheres_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,

# ==========================================================================================
# ==== Fluid
fluid_smoothing_length = 1.0 * fluid_particle_spacing
fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()

fluid_density_calculator = ContinuityDensity()
Expand All @@ -58,12 +58,12 @@ sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calc
surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
correction=AkinciFreeSurfaceCorrection(fluid_density))

# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator,
# state_equation, fluid_smoothing_kernel,
# fluid_smoothing_length,
# viscosity=viscosity,
# acceleration=(0.0, -gravity),
# correction=AkinciFreeSurfaceCorrection(fluid_density))
# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator,
# state_equation, fluid_smoothing_kernel,
# fluid_smoothing_length,
# viscosity=viscosity,
# acceleration=(0.0, -gravity),
# correction=AkinciFreeSurfaceCorrection(fluid_density))

# sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
# fluid_smoothing_length,
Expand Down
12 changes: 12 additions & 0 deletions src/general/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,3 +147,15 @@ end

# Only for systems requiring a mandatory callback
reset_callback_flag!(system) = system

function neighbor_number(::Val{D}, initial_particle_spacing, smoothing_length) where {D}
throw(ArgumentError("Unsupported dimension: $D"))
end

@inline function neighbor_number(::Val{2}, initial_particle_spacing, smoothing_length)
return pi * smoothing_length^2 / initial_particle_spacing^2
end

@inline function neighbor_number(::Val{3}, initial_particle_spacing, smoothing_length)
return 4.0 / 3.0 * pi * smoothing_length^3 / initial_particle_spacing^3
end
2 changes: 2 additions & 0 deletions src/schemes/boundary/boundary.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,3 +4,5 @@ include("open_boundary/boundary_zones.jl")
include("open_boundary/system.jl")
# Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`
# and the `TotalLagrangianSPHSystem` and are therefore included later.

@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel)
17 changes: 15 additions & 2 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,15 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
cache = (; create_cache_model(viscosity, n_particles, NDIMS)...,
create_cache_model(initial_density, density_calculator)...,
create_cache_model(correction, initial_density, NDIMS, n_particles)...,
(; colorfield=zeros(eltype(smoothing_length), n_particles))...)
(; colorfield_bnd=zeros(eltype(smoothing_length), n_particles),
colorfield=zeros(eltype(smoothing_length), n_particles),
neighbor_count=zeros(eltype(smoothing_length), n_particles))...,
(; neighbor_number=[0.0])...)

return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation,
density_calculator, smoothing_kernel,
smoothing_length, viscosity, correction, cache)
smoothing_length, viscosity, correction,
cache)
end

@doc raw"""
Expand Down Expand Up @@ -507,3 +511,12 @@ end
@inline function correction_matrix(system::BoundarySystem, particle)
extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle)
end

function initialize_boundary_model!(boundary_model::BoundaryModelDummyParticles,
initial_particle_spacing)
(; smoothing_kernel, smoothing_length) = boundary_model
boundary_model.cache.neighbor_number[1] = neighbor_number(Val(ndims(boundary_model)),
initial_particle_spacing,
compact_support(smoothing_kernel,
smoothing_length))
end
17 changes: 12 additions & 5 deletions src/schemes/boundary/system.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,8 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM,
ismoving, adhesion_coefficient, cache, buffer)
ELTYPE = eltype(coordinates)

initialize_boundary_model!(boundary_model, initial_condition.particle_spacing)

new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition),
typeof(coordinates), typeof(movement), typeof(ismoving),
typeof(cache)}(initial_condition, coordinates, boundary_model, movement,
Expand Down Expand Up @@ -397,10 +399,15 @@ function initialize!(system::BoundarySPHSystem, neighborhood_search)
(; smoothing_kernel, smoothing_length) = system.boundary_model

foreach_point_neighbor(system, system,
system_coords, system_coords,
neighborhood_search, points=eachparticle(system)) do particle, neighbor, pos_diff, distance
system.boundary_model.cache.colorfield[particle] += kernel(smoothing_kernel,
distance,
smoothing_length)
system_coords, system_coords,
neighborhood_search,
points=eachparticle(system)) do particle, neighbor, pos_diff,
distance
system.boundary_model.cache.colorfield_bnd[particle] += kernel(smoothing_kernel,
distance,
smoothing_length)
system.boundary_model.cache.neighbor_count[particle] += 1
end
end

initialize_boundary_model!(boundary_model, initial_particle_spacing) = boundary_model
20 changes: 14 additions & 6 deletions src/schemes/fluid/surface_normal_sph.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,15 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system,
# First we need to calculate the smoothed colorfield values
# TODO: move colorfield to extra step
# TODO: this is only correct for a single fluid
colorfield = zeros(eltype(neighbor_system), nparticles(neighbor_system))
# colorfield = zeros(eltype(neighbor_system), nparticles(neighbor_system))
# colorfield_bnd = zeros(eltype(neighbor_system), nparticles(neighbor_system))

foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
colorfield[neighbor] += kernel(smoothing_kernel, distance, smoothing_length)
neighbor_system.boundary_model.cache.colorfield[neighbor] += kernel(smoothing_kernel,
distance,
smoothing_length)
end

# foreach_point_neighbor(neighbor_system, neighbor_system,
Expand All @@ -93,14 +95,20 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system,
# end

# Since we don't want to calculate the unused boundary weight sum we normalize against the maximum value
colorfield = colorfield ./ (2.0 * maximum(colorfield))
# colorfield = colorfield ./ (2.0 * maximum(colorfield))

# println(neighbor_system.boundary_model.cache.neighbor_number)

colorfield_bnd = neighbor_system.boundary_model.cache.colorfield_bnd
colorfield = neighbor_system.boundary_model.cache.colorfield

# println(colorfield)

# foreach_point_neighbor(system, neighbor_system,
# system_coords, neighbor_system_coords,
# nhs) do particle, neighbor, pos_diff, distance
# colorfield[neighbor] = colorfield[neighbor] /
# (colorfield[neighbor] +
# neighbor_system.boundary_model.cache.colorfield[neighbor])
# neighbor_system.boundary_model.cache.colorfield[neighbor] = colorfield[neighbor] /
# (colorfield[neighbor] + colorfield_bnd[neighbor])
# end

# colorfield = colorfield / (colorfield + neighbor_system.boundary_model.cache.colorfield)
Expand Down
11 changes: 7 additions & 4 deletions src/visualization/write2vtk.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,16 +358,19 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system;
write_meta_data=true)
if write_meta_data
vtk["boundary_model"] = "BoundaryModelDummyParticles"
vtk["smoothing_kernel"] = type2string(system.boundary_model.smoothing_kernel)
vtk["smoothing_length"] = system.boundary_model.smoothing_length
vtk["density_calculator"] = type2string(system.boundary_model.density_calculator)
vtk["state_equation"] = type2string(system.boundary_model.state_equation)
vtk["smoothing_kernel"] = type2string(model.smoothing_kernel)
vtk["smoothing_length"] = model.smoothing_length
vtk["density_calculator"] = type2string(model.density_calculator)
vtk["state_equation"] = type2string(model.state_equation)
vtk["viscosity_model"] = type2string(model.viscosity)
end

vtk["hydrodynamic_density"] = [particle_density(v, system, particle)
for particle in eachparticle(system)]
vtk["pressure"] = model.pressure
vtk["colorfield_bnd"] = model.cache.colorfield_bnd
vtk["colorfield"] = model.cache.colorfield
vtk["neighbor_count"] = model.cache.neighbor_count

if model.viscosity isa ViscosityAdami
vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :)
Expand Down

0 comments on commit 6bc5567

Please sign in to comment.