From cc713b49abab7906c5d41f542d9fee63c6da7c98 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Jun 2023 14:06:21 +0200 Subject: [PATCH 01/49] working --- .../fluid/weakly_compressible_sph/rhs.jl | 25 ++++++++++++------- .../fluid/weakly_compressible_sph/system.jl | 6 ++++- 2 files changed, 21 insertions(+), 10 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 341ce4a79..879c09265 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -4,12 +4,15 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::WeaklyCompressibleSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem) @unpack density_calculator, state_equation, viscosity, smoothing_length, - correction = particle_system + correction, dv_pressure, dv_viscosity = particle_system @unpack sound_speed = state_equation system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + set_zero!(dv_pressure) + set_zero!(dv_viscosity) + # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, @@ -32,23 +35,27 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = particle_system.mass[particle] m_b = neighbor_system.mass[neighbor] - dv_pressure = pressure_correction * (-m_b * + tmp_a = extract_svector(dv_pressure, particle_system, particle) + dv_pressure[:, particle] = tmp_a + pressure_correction * (-m_b * (particle_system.pressure[particle] / rho_a^2 + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) - dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, + tmp_b = extract_svector(dv_viscosity, particle_system, particle) + dv_viscosity[:, particle] = tmp_b + viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b) + continuity_equation!(dv, density_calculator, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + particle_system, neighbor_system) + end + + for particle in each_moving_particle(particle_system) for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity[i] + dv[i, particle] += dv_pressure[i, particle] + dv_viscosity[i, particle] end - - continuity_equation!(dv, density_calculator, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - particle_system, neighbor_system) end return dv diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 8b09c3fb6..abaf7f338 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -20,6 +20,8 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] pressure :: Array{ELTYPE, 1} # [particle] + dv_pressure :: Array{ELTYPE, 2} # [dimension, particle] + dv_viscosity :: Array{ELTYPE, 2} # [dimension, particle] density_calculator :: DC state_equation :: SE smoothing_kernel :: K @@ -42,6 +44,8 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} mass = copy(initial_condition.mass) pressure = Vector{ELTYPE}(undef, n_particles) + dv_pressure = Array{ELTYPE, 2}(undef, NDIMS, n_particles) + dv_viscosity = Array{ELTYPE, 2}(undef, NDIMS, n_particles) if ndims(smoothing_kernel) != NDIMS throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) @@ -59,7 +63,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} return new{NDIMS, ELTYPE, typeof(density_calculator), typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), typeof(cache) - }(initial_condition, mass, pressure, density_calculator, state_equation, + }(initial_condition, mass, pressure, dv_pressure, dv_viscosity, density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity, acceleration_, correction, cache) end From 923c71ccb4daf67540dc2560688ed1639f9804c4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Jun 2023 16:16:52 +0200 Subject: [PATCH 02/49] combine --- .../fluid/weakly_compressible_sph/rhs.jl | 79 ++++++------------- .../weakly_compressible_sph/viscosity.jl | 8 +- 2 files changed, 24 insertions(+), 63 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 879c09265..25e4c6229 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -2,11 +2,12 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, - neighbor_system::WeaklyCompressibleSPHSystem) - @unpack density_calculator, state_equation, viscosity, smoothing_length, + neighbor_system) + @unpack density_calculator, state_equation, smoothing_length, correction, dv_pressure, dv_viscosity = particle_system @unpack sound_speed = state_equation + viscosity = neighbor_system isa WeaklyCompressibleSPHSystem ? particle_system.viscosity : neighbor_system.boundary_model.viscosity system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -29,22 +30,19 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, rho_mean) - # Pressure forces grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) - m_a = particle_system.mass[particle] - m_b = neighbor_system.mass[neighbor] + m_a = hydrodynamic_mass(particle_system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) tmp_a = extract_svector(dv_pressure, particle_system, particle) - dv_pressure[:, particle] = tmp_a + pressure_correction * (-m_b * - (particle_system.pressure[particle] / rho_a^2 + - neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) + dv_pressure[:, particle] = tmp_a + calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) tmp_b = extract_svector(dv_viscosity, particle_system, particle) dv_viscosity[:, particle] = tmp_b + viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b) + sound_speed, m_a, m_b, rho_mean) continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, @@ -61,6 +59,21 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) + return (-m_b * + (particle_system.pressure[particle] / rho_a^2 + + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * pressure_correction +end + +@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) + @unpack boundary_model = neighbor_system + + return boundary_particle_impact(particle, neighbor, boundary_model, + v_particle_system, v_neighbor_system, + particle_system, neighbor_system, + pos_diff, distance, m_b) +end + @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -83,51 +96,3 @@ end particle_system, neighbor_system) return dv end - -# Fluid-boundary and fluid-solid interaction -function interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::WeaklyCompressibleSPHSystem, - neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}) - @unpack density_calculator, state_equation, smoothing_length = particle_system - @unpack sound_speed = state_equation - @unpack boundary_model = neighbor_system - @unpack viscosity = boundary_model - - system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - # Loop over all pairs of particles and neighbors within the kernel cutoff. - for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0. - distance < sqrt(eps()) && return - - # In fluid-solid interaction, use the "hydrodynamic mass" of the solid particles - # corresponding to the rest density of the fluid and not the material density. - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - dv_viscosity = viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, particle, - neighbor, pos_diff, distance, sound_speed, m_a, m_b) - - # Boundary forces - dv_boundary = boundary_particle_impact(particle, neighbor, boundary_model, - v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - pos_diff, distance, m_b) - - for i in 1:ndims(particle_system) - dv[i, particle] += dv_boundary[i] + dv_viscosity[i] - end - - continuity_equation!(dv, density_calculator, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - particle_system, neighbor_system) - end - - return dv -end diff --git a/src/schemes/fluid/weakly_compressible_sph/viscosity.jl b/src/schemes/fluid/weakly_compressible_sph/viscosity.jl index aa82ee983..b49aed102 100644 --- a/src/schemes/fluid/weakly_compressible_sph/viscosity.jl +++ b/src/schemes/fluid/weakly_compressible_sph/viscosity.jl @@ -2,7 +2,7 @@ struct NoViscosity end @inline function (::NoViscosity)(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b) + sound_speed, m_a, m_b, rho_mean) return SVector(ntuple(_ -> 0.0, Val(ndims(particle_system)))) end @@ -62,17 +62,13 @@ end v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b) + distance, sound_speed, m_a, m_b, rho_mean) @unpack smoothing_length = particle_system v_a = current_velocity(v_particle_system, particle_system, particle) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 - pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, smoothing_length) if pi_ab < eps() From 3de3008de0a01bc1575c1e4072cf40e6d0d1b2f1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Jun 2023 16:41:38 +0200 Subject: [PATCH 03/49] format --- src/schemes/boundary/system.jl | 4 ++ .../fluid/weakly_compressible_sph/rhs.jl | 49 +++++++++++++------ .../fluid/weakly_compressible_sph/system.jl | 7 ++- .../weakly_compressible_sph/viscosity.jl | 3 +- .../solid/total_lagrangian_sph/system.jl | 4 ++ 5 files changed, 49 insertions(+), 18 deletions(-) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 78eab465d..e9af50622 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -211,3 +211,7 @@ function restart_with!(system, model, ::ContinuityDensity, v, u) return system end + +function system_viscosity(system::BoundarySPHSystem) + return system.boundary_model.viscosity +end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 25e4c6229..eac036e98 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -7,7 +7,7 @@ function interact!(dv, v_particle_system, u_particle_system, correction, dv_pressure, dv_viscosity = particle_system @unpack sound_speed = state_equation - viscosity = neighbor_system isa WeaklyCompressibleSPHSystem ? particle_system.viscosity : neighbor_system.boundary_model.viscosity + viscosity = system_viscosity(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -36,18 +36,25 @@ function interact!(dv, v_particle_system, u_particle_system, m_b = hydrodynamic_mass(neighbor_system, neighbor) tmp_a = extract_svector(dv_pressure, particle_system, particle) - dv_pressure[:, particle] = tmp_a + calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) + dv_pressure[:, particle] = tmp_a + + calc_pressure_eq(pressure_correction, m_b, particle, + particle_system, v_particle_system, + neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, + pos_diff, distance, grad_kernel) tmp_b = extract_svector(dv_viscosity, particle_system, particle) - dv_viscosity[:, particle] = tmp_b + viscosity_correction * viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_mean) + dv_viscosity[:, particle] = tmp_b + + viscosity_correction * + viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) continuity_equation!(dv, density_calculator, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - particle_system, neighbor_system) + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + particle_system, neighbor_system) end for particle in each_moving_particle(particle_system) @@ -59,19 +66,29 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end -@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) +@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, + v_particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel) return (-m_b * - (particle_system.pressure[particle] / rho_a^2 + - neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * pressure_correction + (particle_system.pressure[particle] / rho_a^2 + + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * + pressure_correction end -@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) +@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, + v_particle_system, neighbor, + neighbor_system::Union{BoundarySPHSystem, + TotalLagrangianSPHSystem}, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel) @unpack boundary_model = neighbor_system return boundary_particle_impact(particle, neighbor, boundary_model, - v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - pos_diff, distance, m_b) + v_particle_system, v_neighbor_system, + particle_system, neighbor_system, + pos_diff, distance, m_b) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index abaf7f338..77e49508a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -63,7 +63,8 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} return new{NDIMS, ELTYPE, typeof(density_calculator), typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), typeof(cache) - }(initial_condition, mass, pressure, dv_pressure, dv_viscosity, density_calculator, state_equation, + }(initial_condition, mass, pressure, dv_pressure, dv_viscosity, + density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity, acceleration_, correction, cache) end @@ -279,3 +280,7 @@ function restart_with!(system, ::ContinuityDensity, v, u) return system end + +function system_viscosity(system::WeaklyCompressibleSPHSystem) + return system.viscosity +end diff --git a/src/schemes/fluid/weakly_compressible_sph/viscosity.jl b/src/schemes/fluid/weakly_compressible_sph/viscosity.jl index b49aed102..65b51f697 100644 --- a/src/schemes/fluid/weakly_compressible_sph/viscosity.jl +++ b/src/schemes/fluid/weakly_compressible_sph/viscosity.jl @@ -62,7 +62,8 @@ end v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, rho_mean) + distance, sound_speed, m_a, m_b, + rho_mean) @unpack smoothing_length = particle_system v_a = current_velocity(v_particle_system, particle_system, particle) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index c0d3b57b6..b177b16f5 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -457,3 +457,7 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u) # This is dispatched in the boundary system.jl file restart_with!(system, system.boundary_model, v, u) end + +function system_viscosity(system::TotalLagrangianSPHSystem) + return system.boundary_model.viscosity +end From da9e34af3c69f18e61cb062a263a9ad4cce3b331 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Jun 2023 18:11:24 +0200 Subject: [PATCH 04/49] write values --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 8 ++++---- src/visualization/write2vtk.jl | 2 ++ 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index eac036e98..ff40d0587 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -35,16 +35,16 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - tmp_a = extract_svector(dv_pressure, particle_system, particle) - dv_pressure[:, particle] = tmp_a + + tmp_pressure = extract_svector(dv_pressure, particle_system, particle) + dv_pressure[:, particle] = tmp_pressure + calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) - tmp_b = extract_svector(dv_viscosity, particle_system, particle) - dv_viscosity[:, particle] = tmp_b + + tmp_viscosity = extract_svector(dv_viscosity, particle_system, particle) + dv_viscosity[:, particle] = tmp_viscosity + viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 37ec6ca37..ba3186338 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -74,6 +74,8 @@ function write2vtk!(vtk, v, u, t, system::WeaklyCompressibleSPHSystem) vtk["density"] = [particle_density(v, system, particle) for particle in eachparticle(system)] vtk["pressure"] = system.pressure + vtk["dv_pressure"] = view(system.dv_pressure, 1:ndims(system), :) + vtk["dv_viscosity"] = view(system.dv_viscosity, 1:ndims(system), :) return vtk end From 3c0d1b2ae466a9d0511665271529ea2790d80150 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 20 Jun 2023 18:50:44 +0200 Subject: [PATCH 05/49] fix --- src/schemes/solid/total_lagrangian_sph/rhs.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index f74164666..a3bcfcaa9 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -84,11 +84,15 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = particle_density(v_particle_system, particle_system, particle) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + # use `m_a` to get the same viscosity as for the fluid-solid direction. dv_viscosity = viscosity(neighbor_system, particle_system, v_neighbor_system, v_particle_system, neighbor, particle, - pos_diff, distance, sound_speed, m_b, m_a) + pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces dv_boundary = boundary_particle_impact(neighbor, particle, boundary_model, From 05fae66104de9eb0f42e5517706bf90720758688 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 Jun 2023 14:35:58 +0200 Subject: [PATCH 06/49] add consistent form of pressure eq --- .../fluid/weakly_compressible_sph/rhs.jl | 25 +++++++++++++++---- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index ff40d0587..961aa0f61 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -41,7 +41,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel) + pos_diff, distance, grad_kernel, density_calculator) tmp_viscosity = extract_svector(dv_viscosity, particle_system, particle) dv_viscosity[:, particle] = tmp_viscosity + @@ -66,15 +66,30 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +# As shown in "Variational and momentum preservation aspects of Smooth +# Particle Hydrodynamic formulations" by Boney and Lok, 1999 for a consistent formulation this form has to be used with ContinuityDensity @inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel) + grad_kernel, ::ContinuityDensity) return (-m_b * - (particle_system.pressure[particle] / rho_a^2 + - neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * - pressure_correction + (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) / (rho_a * rho_b) * grad_kernel) * + pressure_correction +end + +# As shown in "Variational and momentum preservation aspects of Smooth +# Particle Hydrodynamic formulations" by Boney and Lok, 1999 for a consistent formulation this form has to be used with SummationDensity +@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, + v_particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, ::SummationDensity) + +return (-m_b * + (particle_system.pressure[particle] / rho_a^2 + + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * + pressure_correction end @inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, From b4c3f35748d2ce0269ddb7de0a7fa4fd63a56308 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 Jun 2023 14:42:40 +0200 Subject: [PATCH 07/49] typo --- .../fluid/weakly_compressible_sph/rhs.jl | 31 ++++++++++--------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 961aa0f61..4ec224e85 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -41,7 +41,8 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel, density_calculator) + pos_diff, distance, grad_kernel, + density_calculator) tmp_viscosity = extract_svector(dv_viscosity, particle_system, particle) dv_viscosity[:, particle] = tmp_viscosity + @@ -67,29 +68,31 @@ function interact!(dv, v_particle_system, u_particle_system, end # As shown in "Variational and momentum preservation aspects of Smooth -# Particle Hydrodynamic formulations" by Boney and Lok, 1999 for a consistent formulation this form has to be used with ContinuityDensity +# Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation +# this form has to be used with ContinuityDensity @inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel, ::ContinuityDensity) return (-m_b * - (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) / (rho_a * rho_b) * grad_kernel) * - pressure_correction + (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) / + (rho_a * rho_b) * grad_kernel) * + pressure_correction end # As shown in "Variational and momentum preservation aspects of Smooth -# Particle Hydrodynamic formulations" by Boney and Lok, 1999 for a consistent formulation this form has to be used with SummationDensity +# Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation +# this form has to be used with SummationDensity @inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, - neighbor_system::WeaklyCompressibleSPHSystem, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel, ::SummationDensity) - -return (-m_b * - (particle_system.pressure[particle] / rho_a^2 + - neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * - pressure_correction + v_particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, ::SummationDensity) + return (-m_b * + (particle_system.pressure[particle] / rho_a^2 + + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * + pressure_correction end @inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, From 53d8b2d4094cdb3d522a06b4d34ee0c383e8fabf Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 28 Jun 2023 01:04:50 +0200 Subject: [PATCH 08/49] optimize --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index e40dc0ad4..99c9e3a3d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -35,17 +35,13 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - tmp_pressure = extract_svector(dv_pressure, particle_system, particle) - dv_pressure[:, particle] = tmp_pressure + - calc_pressure_eq(pressure_correction, m_b, particle, + @views dv_pressure[:, particle] += calc_pressure_eq(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) - tmp_viscosity = extract_svector(dv_viscosity, particle_system, particle) - dv_viscosity[:, particle] = tmp_viscosity + - viscosity_correction * + @views dv_viscosity[:, particle] += viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, From c2f00c8dd3aaf76c38ed0db48c4ab7ec8832723a Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 28 Jun 2023 17:05:21 +0200 Subject: [PATCH 09/49] remove global dv arrays --- .../fluid/weakly_compressible_sph/rhs.jl | 107 ++++++++++-------- .../fluid/weakly_compressible_sph/system.jl | 6 +- src/visualization/write2vtk.jl | 2 - 3 files changed, 61 insertions(+), 54 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 99c9e3a3d..c11a5005a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -4,58 +4,71 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) @unpack density_calculator, state_equation, smoothing_length, - correction, dv_pressure, dv_viscosity = particle_system + correction = particle_system @unpack sound_speed = state_equation viscosity = system_viscosity(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - - set_zero!(dv_pressure) - set_zero!(dv_viscosity) - - # Loop over all pairs of particles and neighbors within the kernel cutoff. - for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - # Only consider particles with a distance > 0. - distance < sqrt(eps()) && return - - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 - - # Determine correction values - viscosity_correction, pressure_correction = free_surface_correction(correction, - particle_system, - rho_mean) - - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - @views dv_pressure[:, particle] += calc_pressure_eq(pressure_correction, m_b, particle, - particle_system, v_particle_system, - neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel) - - @views dv_viscosity[:, particle] += viscosity_correction * - viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_mean) - - continuity_equation!(dv, density_calculator, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - particle_system, neighbor_system, grad_kernel) - end + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + @threaded for particle in each_moving_particle(particle_system) + particle_coords = extract_svector(system_coords, particle_system, particle) + + #Note: don't move outside of @threaded block causes a bug (@2/2023) + NDIMS = ndims(particle_system) + + dv_viscosity = zero(SVector{NDIMS, eltype(particle_system)}) + dv_pressure = zero(SVector{NDIMS, eltype(particle_system)}) + #dv_surface_tension = zero(SVector{NDIMS, eltype(particle_system)}) + + for neighbor in eachneighbor(particle_coords, neighborhood_search) + neighbor_coords = extract_svector(neighbor_system_coords, neighbor_system, + neighbor) + + pos_diff = particle_coords - neighbor_coords + distance2 = dot(pos_diff, pos_diff) + + + if eps() < distance2 <= compact_support(particle_system, neighbor_system)^2 + distance = sqrt(distance2) + + rho_a = particle_density(v_particle_system, particle_system, particle) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + + + # Determine correction values + viscosity_correction, pressure_correction = free_surface_correction(correction, + particle_system, + rho_mean) + + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + + m_a = hydrodynamic_mass(particle_system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + dv_pressure += calc_pressure_eq(pressure_correction, m_b, particle, + particle_system, v_particle_system, + neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, + pos_diff, distance, grad_kernel) + + + dv_viscosity += viscosity_correction * + viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) + + continuity_equation!(dv, density_calculator, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + particle_system, neighbor_system, grad_kernel) + end + end - for particle in each_moving_particle(particle_system) - for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i, particle] + dv_viscosity[i, particle] + for i in 1:NDIMS + dv[i, particle] += dv_pressure[i] + dv_viscosity[i] end end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index d033d0dde..a036f8465 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -20,8 +20,6 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] pressure :: Array{ELTYPE, 1} # [particle] - dv_pressure :: Array{ELTYPE, 2} # [dimension, particle] - dv_viscosity :: Array{ELTYPE, 2} # [dimension, particle] density_calculator :: DC state_equation :: SE smoothing_kernel :: K @@ -44,8 +42,6 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} mass = copy(initial_condition.mass) pressure = Vector{ELTYPE}(undef, n_particles) - dv_pressure = Array{ELTYPE, 2}(undef, NDIMS, n_particles) - dv_viscosity = Array{ELTYPE, 2}(undef, NDIMS, n_particles) if ndims(smoothing_kernel) != NDIMS throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) @@ -70,7 +66,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, COR, C} return new{NDIMS, ELTYPE, typeof(density_calculator), typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), typeof(cache) - }(initial_condition, mass, pressure, dv_pressure, dv_viscosity, + }(initial_condition, mass, pressure, density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity, acceleration_, correction, cache) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 87515fe5e..12843ef66 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -74,8 +74,6 @@ function write2vtk!(vtk, v, u, t, system::WeaklyCompressibleSPHSystem) vtk["density"] = [particle_density(v, system, particle) for particle in eachparticle(system)] vtk["pressure"] = system.pressure - vtk["dv_pressure"] = view(system.dv_pressure, 1:ndims(system), :) - vtk["dv_viscosity"] = view(system.dv_viscosity, 1:ndims(system), :) return vtk end From 0c1792ba0740feea8450542ca2970329e4970369 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 29 Jun 2023 11:27:46 +0200 Subject: [PATCH 10/49] add example for debugging --- .../fluid/weakly_compressible_sph/rhs.jl | 7 ++ src/visualization/write2vtk.jl | 83 +++++++++++++++++- test.vtu | Bin 0 -> 36342 bytes 3 files changed, 89 insertions(+), 1 deletion(-) create mode 100644 test.vtu diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index c11a5005a..56c89d110 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -11,6 +11,9 @@ function interact!(dv, v_particle_system, u_particle_system, system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + # Example + # debug_array = zeros(eltype(particle_system), ndims(particle_system), nparticles(particle_system)) + @threaded for particle in each_moving_particle(particle_system) particle_coords = extract_svector(system_coords, particle_system, particle) @@ -69,8 +72,12 @@ function interact!(dv, v_particle_system, u_particle_system, for i in 1:NDIMS dv[i, particle] += dv_pressure[i] + dv_viscosity[i] + # example + # debug_array[i, particle] = dv_pressure[i] end end + # example + # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, custom_value=Dict("debug"=>debug_array), prefix="debug") return dv end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 12843ef66..d0a5d9e75 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -1,3 +1,25 @@ +""" + trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", custom_quantities...) + +Converts Trixi simulation data to VTK format. + +# Arguments +- `vu_ode::ODESolution`: Solution of the Trixi ODE system. +- `semi::SemiDiscretization`: Semi-discretization of the TrixiParticles simulation. +- `t::Float64`: Current time of the simulation. +- `iter::Union{Nothing, Int}` (optional): Iteration number. Defaults to `nothing`. +- `output_directory::AbstractString` (optional): Output directory path. Defaults to `"out"`. +- `prefix::AbstractString` (optional): Prefix for output files. Defaults to an empty string. +- `custom_quantities...`: Additional custom quantities to include in the VTK output. + +# Details +This function converts TrixiParticles simulation data to VTK format. +It iterates over each system in the semi-discretization and calls the overloaded `trixi2vtk` function for each system. + +# Example +```julia +trixi2vtk(vu, semi, 0.0, iter=1, output_directory="output", prefix="solution", velocity=compute_velocity) +""" function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", custom_quantities...) @unpack systems = semi @@ -19,8 +41,33 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix end end +""" + trixi2vtk(v, u, t, system; output_directory="out", prefix="", iter=nothing, + system_name=vtkname(system), custom_quantities...) + +Converts a single Trixi system data to VTK format. + +# Arguments +- `v::AbstractVector`: Vector of solution variables. +- `u::AbstractVector`: Vector of unknowns. +- `t::Float64`: Current time of the simulation. +- `system::TrixiSystem`: Trixi system to convert to VTK. +- `output_directory::AbstractString` (optional): Output directory path. Defaults to "out". +- `prefix::AbstractString` (optional): Prefix for output files. Defaults to an empty string. +- `iter::Union{Nothing, Int}` (optional): Iteration number. Defaults to nothing. +- `system_name::AbstractString` (optional): Name of the system for the output file. Defaults to the VTK name of the system. +- `custom_quantities...`: Additional custom quantities to include in the VTK output. + +# Details +This function converts a single Trixi system's data to VTK format. It creates the necessary VTK file, +writes the solution variables and unknowns to the file, and includes additional custom quantities if provided. + +# Example +trixi2vtk(v, u, 0.0, fluid_system, output_directory="output", prefix="solution", velocity=compute_velocity) + +""" function trixi2vtk(v, u, t, system; output_directory="out", prefix="", iter=nothing, - system_name=vtkname(system), + system_name=vtkname(system), custom_value=nothing, custom_quantities...) mkpath(output_directory) @@ -41,6 +88,16 @@ function trixi2vtk(v, u, t, system; output_directory="out", prefix="", iter=noth # Store particle index vtk["index"] = eachparticle(system) + if custom_value !== nothing + for (key, value) in custom_value + if axes(value, 1) == ndims(system) + vtk[string(key)] = view(value, 1:ndims(system), :) + else + vtk[string(key)] = value + end + end + end + # Extract custom quantities for this system for (key, func) in custom_quantities value = func(v, u, t, system) @@ -51,6 +108,30 @@ function trixi2vtk(v, u, t, system; output_directory="out", prefix="", iter=noth end end +""" + trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates") + +Converts coordinate data to VTK format. + +# Arguments +- `coordinates::AbstractMatrix`: Matrix of coordinate data. +- `output_directory::AbstractString` (optional): Output directory path. Defaults to `"out"`. +- `prefix::AbstractString` (optional): Prefix for the output file. Defaults to an empty string. +- `filename::AbstractString` (optional): Name of the output file. Defaults to `"coordinates"`. + +# Details +This function converts coordinate data to VTK format. +It creates a VTK file with the specified filename and saves the coordinate data as points in the VTK file. +Each coordinate is treated as a vertex cell in the VTK file. + +# Returns +- `file::AbstractString`: Path to the generated VTK file. + +# Example +```julia +coordinates = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] +vtk_file = trixi2vtk(coordinates, output_directory="output", prefix="data", filename="coords") +""" function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates") mkpath(output_directory) file = prefix === "" ? joinpath(output_directory, filename) : diff --git a/test.vtu b/test.vtu new file mode 100644 index 0000000000000000000000000000000000000000..038b61c2d0b4751ea6b27ea01c87d023d1ce4653 GIT binary patch literal 36342 zcmeHv30RZI*1tP-gIX7$AfQyJ6$ubo4M|%ov=yP&s!)+VP%g+85h9xw6%~*w$|3>U z`l;=eg1}WlSw#p-NUF9h6#)q(dO^0x5|-pUGf61e zzcb;ID4_>^wGLB*Xx{#Qmg^1m4c2Q>{oMUMy#2f^*9RZ++-$aI^;{;?&sm{>bG9& z3#u!OI!oDnLYePx4^X=oxQl|8}$68c_HkSksK zHLZO8UGV~KT@S*yme5aQ3hVRtqhfKy^;(|(K?hwAS*~{t2%!3TP(dmFo}M%+j2NsZ z$fS@6Go_O7%4BbUXAM#nNmk(Q@8?H#KjeMb`w;#2*J1c)>oDA6Ks5e!E&ocIWcjqe zAjrtT%-rlx1?_?FWd@wFV$o^%OJYn-ji-qr5oRb;vcG6g|3tm7Z1G{LufO{rD%a@G zlxt*aW;kszzoXn&WO-2iXn!c`PxsNt)O?H4pQ`s2QP{e`SpUjani&}XKb#qFIQgM} zNsQr^X=0>i^6Tco05WsQPLiy=slFb+zMQUWsilz!GiXSJ-*XIDDulONwz zYy95vc-0#&iyfA%Sz}dIZNIh2C9P+3!Sof05QNLEeA1mx;8#k(Slkyrl z;o@hv#XXv$E=^HvRKInUVB~hGP+(oI+Yu$&8a1)6bc7$(;7eI8-X|VRrikyTdg(-t zhm}qc#8*|L^hNvAD3j5pBm2cA&3*4r^i7B|Cmy#=^g3@5IGnK)9eG)wDPnqNY#JyX z@6Hte#ERN1>WdWIjk3HfZc%Tvn$Q$mJP%8!S)ncWnfjZ5>kGS#w zmj5g$&tm>EqtivZ#G>c%oht*6Xy48}^GMB{ZNq4-9^iMpwCW96XY1b?o&;SsZ3A>C z8lkJvXaxk^CqU~ElWtagZm2odnet>yF!7pcUIy=_RfAB=*56|?0ceLSffO|JJidqP z0A0`rbnluMpW{bQ8i+~^50?1br1hyi;*D4RL-_A}+%?G?IpZyeaQ`K|EybkG^lJS` z{9k??w3-aNYI@!d7T*&24k*Hm2H#y(OjY&!uefKvyW<6J~tXBeeBssTYl$R<$4 z+QXoT1}{S7qI`o0Seg00ks3ENcT(%pagu|{pEFR}GQGAur5 z{CoIvTo$NHJxU+c6{(Dlkk>}&lJ=o))Mn@g9s~kX>kmL1+;kdrl36@D^l0a9Q_x4Z zH)sy*4Kycy8@kRcI5&+Sf%eE5mC)_m4&9P8&@{P8iQBCv;Eh!k6_yzKXZy3#{;iKc zD>*S!OdCyk^+WM=T7-M6W7}VT)c60X)bxn9@B>)*-!!4b;5 z0$qV4Xv1856S~ZspwGe3+ei_vx}b=N;f_W-)16O1A6eJHj)KWB%|Uz_Xd?{a+&Jr^ zBNVy@x~0mf>j#<^x=`YBV)ZX|e;doY^X%I~{X1V||2rS|M$ji|4;+;bY-!xCTp#>k zJMF-szxsFr98%;n`8cGK%sk)TYkFPrpl(0T_{KV3WNlvI)F5}}NZPND7fjD2^j*`N zHiG#4`uOQkM?@l(Mxi#h*g2St>`tnmydZn5RjXvq9C?M{mB(5zv}1BIwDB9)dE}mt zTfrndhl`Z%{;jtn`_lBmNV?|XpIc%M7>~H_tB+{fLz6t#pzkfN(jCkg& z1OIC7Z)F8XK2odEMSDnD4d%s21>@jyz-_ejy$_Zr)PpXA|1s+3K$nyqk1Xvj(B>&_ z1)J&4S1ygbs-X(sC$)MH10L?j#SD3R?(9&qt`ZhVT*XI#w#tcMPPjgnIYx9}> zyQnLpM?!v|_*IH?V7-~+zig3;#q2J#DROeunaWR0II(_GC&b@=yI%0^e){FnQTl@) zigv+4oE)S-I4CH3wihYH>Cj^m;m*MfO7{)tkUMA$$9nCt-1m zE}4H3X+R=ZI)6IFCjV8M6p0YDG6}nMyO&mjr@;gs+ksuEd+C5n7^R^je)n}a{p20c z|?}3nni_2qJ5p&?%u%3E139!zt9ZqRhNx`1);j z&2>nJ3<~s2`>yrCiQ@tXq)N9ea1XQoi;@}&YO}1o24lw2Z5TU znR))LQ_Nph4RPX&&V{c>!DBb51XZgk33VtA)C*EUVYe01>}P*$D*+hw&<7fs*buH?ug5D0$+_G6(5rj#W{L8BxE@|2uJ zaE}or0I1Dvv8Zbb;3E8nW>vSpdfJfbS^Rx(^T*Vfo|EyH!F?8qs&-**MHY&xavh?$ zZ-Z6mlma48++)#`fn}Il0u}FmM@;x0&EF3jGqmWduCYKrZ8sC2|?v!%g5avLCyH%fPP% z_=SdY>!<`H;YBnXAMbKQa_wPt6;)1wWN zPm*Sj!ja2oGv)OChe~Ygk2vc}LJF-oXLo*U<3|uuocg-ol7tk;zV4s0MuN|S9W{FU$<0MtWg(nM-5{|xE^QU-``U*8>Ovqz9D!np+DU%B4WNt zFL|~dY0CNKdZej5MMoL0K@!hm<^&7c+E*rM=18SBTW4z@@orEqu7

2=jzo1F9cC zs@;@VI^VeAmg#xg0k<&i+x$EChpo<51fPBpub$cc#$&ZY=j|bfF<|sX*(XiD*5Zu$brxf>jo=QS%mL%$TwyB>@RSYlMpwij47yN^DAtoxPeag|Aw2yg^bzwzm*ADX5r*#5G4aW z!3CeCC_j~jq61T+0~$NES|rGZKzqISlKOFg zG^1t;Po*Fxb&7%sGo@}-{rFuSOc6PGzbJb81dsI85E)C!gaIrUNh2HewP)A7SAdWZ zhypvs1?Tz21GghBzZP$nlUgFrK#D5dI1Kf-iB%A2UC{uSNU4Mj)d;mh1Ju$beLT)} zLp01fMBiP8+X~!zlG%2-=rlQa@%{sXZxrtugjQuiAka%xpbM--n3b%S{IV)1ZBisY zC!Y+}EaZ6%oQlF7luIpFH1Ym2U3(uU=%8`5QACnmX;ieQ`rtKFxc7qPzlb@kUwwJ%$G9U6wa8*3eVpH;*r~yt^m#-$P^&F#6YkI zumj*Bzzzr~jbb!4If!r+UMb>2@r%gXK>OuOnbZzLn_&g5gYiNCt?7Xio~0-kgRYXT9z|v%{6grhubLB z(j5`aU?`&ZSet^{2;71vrNiQID&`2w@8%RiJMKnCLa-y{;+V2`W`%nWwl+t4I5Nq?GLv=LXBbwb5FxwLYPa$ zucPpbzyZr6VcH4FXvjeU4S7V2BHCc0W+^-&XUYzDZiT}mK$NgYND$JRrm&<4s<2&% z=TwD=4=A)I6+zo!7KXVz7%(RbK!q&?)H0NCtBYF-YDxNthF*l|^AvGQ#;qTZ8*Y!m ze6Xh?z~><^bOIs_HIWdAWj)$77X)z;ld1i2isF6)U`Yxz9jBhi@Io?(@YFVdMSvy| z$oiuX6lMY@<6{7hA^178vavbzi-tD*53C(+Qg;5pk6EPyS{Rt&uJjRYlcSexAxJb@ zlXP8V_C-ibi2Hs4T$8%#AhIMF089E9WhJNuOG2$3Zlh2OmIO3lNk9*lgxUz)f+ZC` zlv>in7s!$vz40y~ArecnevB*$m`qtx5kgZ81d1S61Pm0Uj$za&##)m>2xajr5GZ2k z;@2Ylg4==MlLi5!CPh3nh7lP!6emDw*pdp#z;(Rt+PZ)Ctt z&#~5mh6nOO+IIKoH7X&vl=Is-83ADNA1!j&MKEX{(k` z8M6o>8wSQtmI(BK=RhcG6a%5jAOy9*a{$N!&w*Y`5E6I}z*68jU&b!h!%Mc zc#?Py(GBc5fD(BQ5G3&&ldjlv08g3cc-P;t+f6<>^1u-uXV5Ar+e_4vI41&7(7v!A*AJJ2%$g`4xmg@r?3dH6!1MQ1zwWqvN;OZCWz<*2bcXa z834kNDCq}|Ew_Iw14Hdb?rr(_3^Rpt@R{lh0BU#YLeh=fg&o@(O$>3eIx|?%hV7TfZ2%2G^;qm5bu+d zUqFz@T(BDj41}M~VhqIeLfB}Widq52QUVX$$K%XOG$|n&(Kj5!Z47Sf!x@=wl4!yD z2&_ms595 z2zpIth3BV%cog8JD}We}W#&P85mJNFth*j4%{qXhI8Pq~7HA>OLaiNcqfiTJ7SKSN z1@w?+p*8}yD9zd-_0-`TP@2^OEdne+FG;gDVZZ`RrqV11%zw+P_dMLw;$;LPR-L1G z8Kt--<2D1gBGe8RV449Vv}_9l)Q;iyPAb4lT%~SE6@CmyCKIg(93>KhQi}xEzC8Hd z5n?j+12!|54_gB(y^WUv8v_gnur^>@U@_rX90~U)a2sfdXt-3g>|htv4!Pn!7$tl; zT>KJT#6}70feSYzL}HXVpspyyL}rvRupzM+Qxo8DNNSFtgqIIH^`T4%tXc@PWG8^z zs459a7wVA6S;B#lQh?G4(n?Fqm1H;@KmnDO@EH@J(vmZC0y+p1X=Eek=T~YDpu#BV z4^>1ZKkfT(pc&9C6dO>n1#J+N8eo8n&j}P7U<-tJv*?mO9_P9t8fG1$?=Hh_1#Y3( zA~}XoY=L74#TIl7q1XcSIEWfK149vJB`dZ-iD4gI+Z>bv#q2 z;f@@MU0s(HjEN!7r4&e=*lZh)0TFC3P)v-S7JfWqf(6m}l?hm`#E+Fc^Z_sN4)zPh zIK9UgF5FD+Tlma2h1K3kSV0qaJ+Y1yqnJGClCS(yo{}Ptc31Hh2Y0cKWO~s~#JNhIR zc^97t^s5}SXgO;}?veK8IuvRjqWY|)Fbl^c0;;K!LheXi%N10}m6mRC`NkZ)q?7{1 z3rUR^Lhkd+kgNa&loFh0OhC?c<^*JISQ?7HV_yU%lM(c7xdMe>$&DbK8M$jl?m5VF zNq&e@J4r?eiP#tI@osin0pRN5$0qT@wm5_`)L|G+Jlt#Lh%$ShJ z2K`G}N+uxV!O}c5AvthH!EBjBfcNZaVu9zZuRYly(xv2%KoIy_{pO+TD9jF zmOZ|8aK?X!M@Tk|Bd6`YySOKHt9S%{Za7!vt79rD@2OlUS7EGEIsc8yC0CWCrz#mM zRL-4KN!_J#@qtRZp-N(|N?M>w@_W` zOHE02NV(XYl5U=oSeTMVPe~T0WG=M->V$pDcKZvJ_KXeo=dajbqS_~Q*=MY_KX=hS zb&vhUpY79)?Gta>rv=+557}ofZ~7{}DP?EVh59DO=BD#MG+jE_l=Px0V{Oy9^rqAU zO&6aurCT;77B{6uG$l_qWzIMID%vb%o7sg5GlsU=`O9XP+|81>W*IBZ&Lx_q?l!yl z$Sj>`mY8prMl(wmm}M?4_$sy_Wk@o_== z)`G;l1!-Xg$>Rl?3ZY-cgr-=BUSNhYbVJWy3%%qOn#2#yP!2tp9Gbc>^kQRZx@lRHXKJ}z`>e=V1$5c~Kr>358NR4k!jW$n>EliD} zrzQv;jxKaKb;9BJc89o1hc7ocoVns~lIn1_%i-8+htn4wPV8}r|Jfnh*dg|&Lrkzk z!cg0j#7&s?LQ^rD~T(~l|BPbbq)?4!pw(xXl3vA5|lA@qb1;ZY^wsRZG1 zJ7FA0_$5JjCP#SEUwF1pcuZ4xI#YPUSs33YjHU=L;Ryk0x$@%Bkim6KGnlU@s~Jl?YMTxCVwV&%20%0sWpQ@x5xs`8qr!Uk&f4rs~+!M^?dUx)#9qxB2Nz;PtP@;R1;4xMV^N~&r^d(CG))IuJKs6#&g9Q zs^J>1#Yr9;lRQp5IUh_>nv`su$nox-* zUQ2IzY`W#Cc8j|8mY2d|58cC_%7>|@hrN`BJqW{|n!{AeFn^9>r`py4Po3RGV)N(ecdpmxy`kS_q|aZH-T6*7 z@6GHki){WPzfK)L-WtCy6FN4QPouYs+{>S<(Ya28w?d=KP=miXqjO^hPc@^< zJcGZ`sdIx9Z?#jGu@irJYv<-x-rCkK%U1q;a;G+#x02jNB=eV+bZ#o)sg-nXE#WIf zcIrm*lq0)LBl$|APJ)OBf22$ijg-$iA^AEX*c3BKLrkE+3wB`{Jal z{Ym#wE0+(f+&rvYKeKW_ROMn*<>p)E>R9C-uH*8)j+>i~>!&*IL04TqzUub3tFE74 zbwA?iV&&=P>FK)H)BP~d#g^yh$8&Y!xks*Xd2fxI>l)WxYup2qTs}&2^GY((n7C{+adS6u-EHDdyXCUumYdHl*ZsHL!w$Px zA9nLP?7Hu;d&sbh-LRYgu&eX1dz51PyNVqyicddLFFLP zZg6Y+&eo2DtxpfMawEv?+sGa6Riu_igU|hfP(dP5A|z%Aag1{$ay>$tsIu-Dk1xWwEMy zSmg{>WecmKh{c@HEBjLK{yn{W-|AI$=#`(>t9+eD7WLt>XEXC;3)3_*UHVWe)e0o$a~L?zxxUQ`Os3p3zg; z+EY=|!xX8P9aFzwu72+u^{S`p<>%BZAE;O4sxt=|Wv3bUYZ&*wXH-3Bl&3N(n;8{_ z45rYr?1ba}O2>Ow9ILt<%P%@s{_I$B(~&vUQWoEGzrN+(4=q(MTFTQ~Dxb7e6t^%Z zNoCQb`xT^nmq}GzQh6e&@)4;bpTratmBkj_uPeHDy{M|Ys64Hx@^Mkc-6H0Acv(#N zeP;N*YvEP=@bcvF%Es`D+u_WSiL!)=`<#h;ITKZV6Xlr`m2DFhr4vl?n!#B-CSN=j z8#>=WYrdgDw87iahGyFg-r8okrNTh3!jPnGuwL8H=(52(mkll44RqWMO}Ga7Tto6o zgLNwn4HFGECK{UWHrTM+(D;$T=0}E>L<4Q2Au-=zQ@-I=nt?9O&{SYR5ExRH-q2c_ zYY=2uxn#8)>(p`$ z({5}`%QfGBW5fPjdP93_F7fV-O?Pv*hTYH&%QYRpK^V`aCJ3efF3iDq#LkaH_+%>z&qChExZDByaG-50s8zvvU0#Wq2ewA_>qZTjiu(!TfvBhA%sI=d-(AMiT2||_jQ#qunopM1 z8!W2}SjKrRuKL3`*4uHlU&OIh;%e;U>doTnLgP3KY^p!7VZCKj>tVxQWmEH+P5l;| zx&s%RvuYe!^(0nZIEyn+uljvG)_T2KH$C#)7<;&9Xt@ZR}uko$f>sxQ)TX)!(qu5hz z+r!fDsrBn&YxLAO_0*Gl>LPnMbJeTgQ)jJHuXR;tuTZbqrCx8SUKgm&Sk}?O3zNvEJCRF4&Q?yrp_)3u|*r?ZFoI z+LoFFE%lZybrCI``K0P?B$hU*)}6#&NvhdRswa}_Xe7?kqUs$*tW8C=K1FP`qMH3h z^;?VT!iqQw;nmjREZy*0uW+_deiW_kZ_LDM77-ni!f2^Kf%_VsBxaCr%cpE zO*kJH>&}^P^!j|_Ytcq3(ZmJYj8<(UE~_w7uOQCTHhM#wxa6|Yo0o};+>O?_6BW5e z8eHPsl}0O85*H^LsU{K^?lxMzo4EXu(b`AE`9!0YMB>tXBei^@0?kO7MpP0QX$tb^ zEX{j;Y5r@mc`C8_3wGqK+L6DkE>FEKf8M6NH#X%jxt{mt_54LXd24*~6}$5^y7T9% z<*iW5U!0bwnwGzCf8Of-`O6>Yt$mz7e{0^#t@%su=BeGyR|v~f4$D^>&(j>I%~1$? zU4iymOprnFdGKs8V2^yh5%%TCgq5;Bprn6x{u3=#5vspAC-;0zp^w}^V-!Sn0Y#5L)8pt^{8wTWC6jBz6 zXD`&{8wa+JM*Z&&1Y)JCPq&C=cUc!!Nd6Q{_Ngv;oD<4;QEs>#9xA^6sq}ea)52Bg zh5!6-^;TSz&3#u@(vy8UPM?cj-TQjT?uGxNuchPve;-FpRqf`?mC2Ufy~#dzNtwgC zD-D6@wfSFf-@k2>&r;3M4{fjI&iv!DS1#v%|8Y_@FGq=S_m2cx{3@^<3B+Ia;Buv} zrQ?6{aa4JJ+}D5nRc5af|A+l)-0XFtoc| Date: Thu, 29 Jun 2023 15:39:13 +0200 Subject: [PATCH 11/49] format --- examples/fluid/rectangular_tank_2d.jl | 2 +- .../fluid/weakly_compressible_sph/rhs.jl | 98 +++++++----------- test.vtu | Bin 36342 -> 0 bytes 3 files changed, 41 insertions(+), 59 deletions(-) delete mode 100644 test.vtu diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl index b09a0f234..4a503221f 100644 --- a/examples/fluid/rectangular_tank_2d.jl +++ b/examples/fluid/rectangular_tank_2d.jl @@ -65,7 +65,7 @@ semi = Semidiscretization(fluid_system, boundary_system, tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) -info_callback = InfoCallback(interval=10) +info_callback = InfoCallback(interval=50) saving_callback = SolutionSavingCallback(dt=0.02) callbacks = CallbackSet(info_callback, saving_callback) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 56c89d110..371ea65c7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -14,66 +14,48 @@ function interact!(dv, v_particle_system, u_particle_system, # Example # debug_array = zeros(eltype(particle_system), ndims(particle_system), nparticles(particle_system)) - @threaded for particle in each_moving_particle(particle_system) - particle_coords = extract_svector(system_coords, particle_system, particle) - - #Note: don't move outside of @threaded block causes a bug (@2/2023) - NDIMS = ndims(particle_system) - - dv_viscosity = zero(SVector{NDIMS, eltype(particle_system)}) - dv_pressure = zero(SVector{NDIMS, eltype(particle_system)}) - #dv_surface_tension = zero(SVector{NDIMS, eltype(particle_system)}) - - for neighbor in eachneighbor(particle_coords, neighborhood_search) - neighbor_coords = extract_svector(neighbor_system_coords, neighbor_system, - neighbor) - - pos_diff = particle_coords - neighbor_coords - distance2 = dot(pos_diff, pos_diff) - - - if eps() < distance2 <= compact_support(particle_system, neighbor_system)^2 - distance = sqrt(distance2) - - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 - - - # Determine correction values - viscosity_correction, pressure_correction = free_surface_correction(correction, - particle_system, - rho_mean) - - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) - - dv_pressure += calc_pressure_eq(pressure_correction, m_b, particle, - particle_system, v_particle_system, - neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel) - - - dv_viscosity += viscosity_correction * - viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_mean) - - continuity_equation!(dv, density_calculator, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - particle_system, neighbor_system, grad_kernel) - end - end - - for i in 1:NDIMS + # Loop over all pairs of particles and neighbors within the kernel cutoff. + for_particle_neighbor(particle_system, neighbor_system, + system_coords, neighbor_system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + rho_a = particle_density(v_particle_system, particle_system, particle) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + + # Determine correction values + viscosity_correction, pressure_correction = free_surface_correction(correction, + particle_system, + rho_mean) + + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + + m_a = hydrodynamic_mass(particle_system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + dv_pressure = calc_pressure_eq(pressure_correction, m_b, particle, + particle_system, v_particle_system, + neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, + pos_diff, distance, grad_kernel) + + dv_viscosity = viscosity_correction * + viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) + + continuity_equation!(dv, density_calculator, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + particle_system, neighbor_system, grad_kernel) + + for i in 1:ndims(particle_system) dv[i, particle] += dv_pressure[i] + dv_viscosity[i] # example - # debug_array[i, particle] = dv_pressure[i] + # debug_array[i, particle] += dv_pressure[i] end end # example diff --git a/test.vtu b/test.vtu deleted file mode 100644 index 038b61c2d0b4751ea6b27ea01c87d023d1ce4653..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 36342 zcmeHv30RZI*1tP-gIX7$AfQyJ6$ubo4M|%ov=yP&s!)+VP%g+85h9xw6%~*w$|3>U z`l;=eg1}WlSw#p-NUF9h6#)q(dO^0x5|-pUGf61e zzcb;ID4_>^wGLB*Xx{#Qmg^1m4c2Q>{oMUMy#2f^*9RZ++-$aI^;{;?&sm{>bG9& z3#u!OI!oDnLYePx4^X=oxQl|8}$68c_HkSksK zHLZO8UGV~KT@S*yme5aQ3hVRtqhfKy^;(|(K?hwAS*~{t2%!3TP(dmFo}M%+j2NsZ z$fS@6Go_O7%4BbUXAM#nNmk(Q@8?H#KjeMb`w;#2*J1c)>oDA6Ks5e!E&ocIWcjqe zAjrtT%-rlx1?_?FWd@wFV$o^%OJYn-ji-qr5oRb;vcG6g|3tm7Z1G{LufO{rD%a@G zlxt*aW;kszzoXn&WO-2iXn!c`PxsNt)O?H4pQ`s2QP{e`SpUjani&}XKb#qFIQgM} zNsQr^X=0>i^6Tco05WsQPLiy=slFb+zMQUWsilz!GiXSJ-*XIDDulONwz zYy95vc-0#&iyfA%Sz}dIZNIh2C9P+3!Sof05QNLEeA1mx;8#k(Slkyrl z;o@hv#XXv$E=^HvRKInUVB~hGP+(oI+Yu$&8a1)6bc7$(;7eI8-X|VRrikyTdg(-t zhm}qc#8*|L^hNvAD3j5pBm2cA&3*4r^i7B|Cmy#=^g3@5IGnK)9eG)wDPnqNY#JyX z@6Hte#ERN1>WdWIjk3HfZc%Tvn$Q$mJP%8!S)ncWnfjZ5>kGS#w zmj5g$&tm>EqtivZ#G>c%oht*6Xy48}^GMB{ZNq4-9^iMpwCW96XY1b?o&;SsZ3A>C z8lkJvXaxk^CqU~ElWtagZm2odnet>yF!7pcUIy=_RfAB=*56|?0ceLSffO|JJidqP z0A0`rbnluMpW{bQ8i+~^50?1br1hyi;*D4RL-_A}+%?G?IpZyeaQ`K|EybkG^lJS` z{9k??w3-aNYI@!d7T*&24k*Hm2H#y(OjY&!uefKvyW<6J~tXBeeBssTYl$R<$4 z+QXoT1}{S7qI`o0Seg00ks3ENcT(%pagu|{pEFR}GQGAur5 z{CoIvTo$NHJxU+c6{(Dlkk>}&lJ=o))Mn@g9s~kX>kmL1+;kdrl36@D^l0a9Q_x4Z zH)sy*4Kycy8@kRcI5&+Sf%eE5mC)_m4&9P8&@{P8iQBCv;Eh!k6_yzKXZy3#{;iKc zD>*S!OdCyk^+WM=T7-M6W7}VT)c60X)bxn9@B>)*-!!4b;5 z0$qV4Xv1856S~ZspwGe3+ei_vx}b=N;f_W-)16O1A6eJHj)KWB%|Uz_Xd?{a+&Jr^ zBNVy@x~0mf>j#<^x=`YBV)ZX|e;doY^X%I~{X1V||2rS|M$ji|4;+;bY-!xCTp#>k zJMF-szxsFr98%;n`8cGK%sk)TYkFPrpl(0T_{KV3WNlvI)F5}}NZPND7fjD2^j*`N zHiG#4`uOQkM?@l(Mxi#h*g2St>`tnmydZn5RjXvq9C?M{mB(5zv}1BIwDB9)dE}mt zTfrndhl`Z%{;jtn`_lBmNV?|XpIc%M7>~H_tB+{fLz6t#pzkfN(jCkg& z1OIC7Z)F8XK2odEMSDnD4d%s21>@jyz-_ejy$_Zr)PpXA|1s+3K$nyqk1Xvj(B>&_ z1)J&4S1ygbs-X(sC$)MH10L?j#SD3R?(9&qt`ZhVT*XI#w#tcMPPjgnIYx9}> zyQnLpM?!v|_*IH?V7-~+zig3;#q2J#DROeunaWR0II(_GC&b@=yI%0^e){FnQTl@) zigv+4oE)S-I4CH3wihYH>Cj^m;m*MfO7{)tkUMA$$9nCt-1m zE}4H3X+R=ZI)6IFCjV8M6p0YDG6}nMyO&mjr@;gs+ksuEd+C5n7^R^je)n}a{p20c z|?}3nni_2qJ5p&?%u%3E139!zt9ZqRhNx`1);j z&2>nJ3<~s2`>yrCiQ@tXq)N9ea1XQoi;@}&YO}1o24lw2Z5TU znR))LQ_Nph4RPX&&V{c>!DBb51XZgk33VtA)C*EUVYe01>}P*$D*+hw&<7fs*buH?ug5D0$+_G6(5rj#W{L8BxE@|2uJ zaE}or0I1Dvv8Zbb;3E8nW>vSpdfJfbS^Rx(^T*Vfo|EyH!F?8qs&-**MHY&xavh?$ zZ-Z6mlma48++)#`fn}Il0u}FmM@;x0&EF3jGqmWduCYKrZ8sC2|?v!%g5avLCyH%fPP% z_=SdY>!<`H;YBnXAMbKQa_wPt6;)1wWN zPm*Sj!ja2oGv)OChe~Ygk2vc}LJF-oXLo*U<3|uuocg-ol7tk;zV4s0MuN|S9W{FU$<0MtWg(nM-5{|xE^QU-``U*8>Ovqz9D!np+DU%B4WNt zFL|~dY0CNKdZej5MMoL0K@!hm<^&7c+E*rM=18SBTW4z@@orEqu7

2=jzo1F9cC zs@;@VI^VeAmg#xg0k<&i+x$EChpo<51fPBpub$cc#$&ZY=j|bfF<|sX*(XiD*5Zu$brxf>jo=QS%mL%$TwyB>@RSYlMpwij47yN^DAtoxPeag|Aw2yg^bzwzm*ADX5r*#5G4aW z!3CeCC_j~jq61T+0~$NES|rGZKzqISlKOFg zG^1t;Po*Fxb&7%sGo@}-{rFuSOc6PGzbJb81dsI85E)C!gaIrUNh2HewP)A7SAdWZ zhypvs1?Tz21GghBzZP$nlUgFrK#D5dI1Kf-iB%A2UC{uSNU4Mj)d;mh1Ju$beLT)} zLp01fMBiP8+X~!zlG%2-=rlQa@%{sXZxrtugjQuiAka%xpbM--n3b%S{IV)1ZBisY zC!Y+}EaZ6%oQlF7luIpFH1Ym2U3(uU=%8`5QACnmX;ieQ`rtKFxc7qPzlb@kUwwJ%$G9U6wa8*3eVpH;*r~yt^m#-$P^&F#6YkI zumj*Bzzzr~jbb!4If!r+UMb>2@r%gXK>OuOnbZzLn_&g5gYiNCt?7Xio~0-kgRYXT9z|v%{6grhubLB z(j5`aU?`&ZSet^{2;71vrNiQID&`2w@8%RiJMKnCLa-y{;+V2`W`%nWwl+t4I5Nq?GLv=LXBbwb5FxwLYPa$ zucPpbzyZr6VcH4FXvjeU4S7V2BHCc0W+^-&XUYzDZiT}mK$NgYND$JRrm&<4s<2&% z=TwD=4=A)I6+zo!7KXVz7%(RbK!q&?)H0NCtBYF-YDxNthF*l|^AvGQ#;qTZ8*Y!m ze6Xh?z~><^bOIs_HIWdAWj)$77X)z;ld1i2isF6)U`Yxz9jBhi@Io?(@YFVdMSvy| z$oiuX6lMY@<6{7hA^178vavbzi-tD*53C(+Qg;5pk6EPyS{Rt&uJjRYlcSexAxJb@ zlXP8V_C-ibi2Hs4T$8%#AhIMF089E9WhJNuOG2$3Zlh2OmIO3lNk9*lgxUz)f+ZC` zlv>in7s!$vz40y~ArecnevB*$m`qtx5kgZ81d1S61Pm0Uj$za&##)m>2xajr5GZ2k z;@2Ylg4==MlLi5!CPh3nh7lP!6emDw*pdp#z;(Rt+PZ)Ctt z&#~5mh6nOO+IIKoH7X&vl=Is-83ADNA1!j&MKEX{(k` z8M6o>8wSQtmI(BK=RhcG6a%5jAOy9*a{$N!&w*Y`5E6I}z*68jU&b!h!%Mc zc#?Py(GBc5fD(BQ5G3&&ldjlv08g3cc-P;t+f6<>^1u-uXV5Ar+e_4vI41&7(7v!A*AJJ2%$g`4xmg@r?3dH6!1MQ1zwWqvN;OZCWz<*2bcXa z834kNDCq}|Ew_Iw14Hdb?rr(_3^Rpt@R{lh0BU#YLeh=fg&o@(O$>3eIx|?%hV7TfZ2%2G^;qm5bu+d zUqFz@T(BDj41}M~VhqIeLfB}Widq52QUVX$$K%XOG$|n&(Kj5!Z47Sf!x@=wl4!yD z2&_ms595 z2zpIth3BV%cog8JD}We}W#&P85mJNFth*j4%{qXhI8Pq~7HA>OLaiNcqfiTJ7SKSN z1@w?+p*8}yD9zd-_0-`TP@2^OEdne+FG;gDVZZ`RrqV11%zw+P_dMLw;$;LPR-L1G z8Kt--<2D1gBGe8RV449Vv}_9l)Q;iyPAb4lT%~SE6@CmyCKIg(93>KhQi}xEzC8Hd z5n?j+12!|54_gB(y^WUv8v_gnur^>@U@_rX90~U)a2sfdXt-3g>|htv4!Pn!7$tl; zT>KJT#6}70feSYzL}HXVpspyyL}rvRupzM+Qxo8DNNSFtgqIIH^`T4%tXc@PWG8^z zs459a7wVA6S;B#lQh?G4(n?Fqm1H;@KmnDO@EH@J(vmZC0y+p1X=Eek=T~YDpu#BV z4^>1ZKkfT(pc&9C6dO>n1#J+N8eo8n&j}P7U<-tJv*?mO9_P9t8fG1$?=Hh_1#Y3( zA~}XoY=L74#TIl7q1XcSIEWfK149vJB`dZ-iD4gI+Z>bv#q2 z;f@@MU0s(HjEN!7r4&e=*lZh)0TFC3P)v-S7JfWqf(6m}l?hm`#E+Fc^Z_sN4)zPh zIK9UgF5FD+Tlma2h1K3kSV0qaJ+Y1yqnJGClCS(yo{}Ptc31Hh2Y0cKWO~s~#JNhIR zc^97t^s5}SXgO;}?veK8IuvRjqWY|)Fbl^c0;;K!LheXi%N10}m6mRC`NkZ)q?7{1 z3rUR^Lhkd+kgNa&loFh0OhC?c<^*JISQ?7HV_yU%lM(c7xdMe>$&DbK8M$jl?m5VF zNq&e@J4r?eiP#tI@osin0pRN5$0qT@wm5_`)L|G+Jlt#Lh%$ShJ z2K`G}N+uxV!O}c5AvthH!EBjBfcNZaVu9zZuRYly(xv2%KoIy_{pO+TD9jF zmOZ|8aK?X!M@Tk|Bd6`YySOKHt9S%{Za7!vt79rD@2OlUS7EGEIsc8yC0CWCrz#mM zRL-4KN!_J#@qtRZp-N(|N?M>w@_W` zOHE02NV(XYl5U=oSeTMVPe~T0WG=M->V$pDcKZvJ_KXeo=dajbqS_~Q*=MY_KX=hS zb&vhUpY79)?Gta>rv=+557}ofZ~7{}DP?EVh59DO=BD#MG+jE_l=Px0V{Oy9^rqAU zO&6aurCT;77B{6uG$l_qWzIMID%vb%o7sg5GlsU=`O9XP+|81>W*IBZ&Lx_q?l!yl z$Sj>`mY8prMl(wmm}M?4_$sy_Wk@o_== z)`G;l1!-Xg$>Rl?3ZY-cgr-=BUSNhYbVJWy3%%qOn#2#yP!2tp9Gbc>^kQRZx@lRHXKJ}z`>e=V1$5c~Kr>358NR4k!jW$n>EliD} zrzQv;jxKaKb;9BJc89o1hc7ocoVns~lIn1_%i-8+htn4wPV8}r|Jfnh*dg|&Lrkzk z!cg0j#7&s?LQ^rD~T(~l|BPbbq)?4!pw(xXl3vA5|lA@qb1;ZY^wsRZG1 zJ7FA0_$5JjCP#SEUwF1pcuZ4xI#YPUSs33YjHU=L;Ryk0x$@%Bkim6KGnlU@s~Jl?YMTxCVwV&%20%0sWpQ@x5xs`8qr!Uk&f4rs~+!M^?dUx)#9qxB2Nz;PtP@;R1;4xMV^N~&r^d(CG))IuJKs6#&g9Q zs^J>1#Yr9;lRQp5IUh_>nv`su$nox-* zUQ2IzY`W#Cc8j|8mY2d|58cC_%7>|@hrN`BJqW{|n!{AeFn^9>r`py4Po3RGV)N(ecdpmxy`kS_q|aZH-T6*7 z@6GHki){WPzfK)L-WtCy6FN4QPouYs+{>S<(Ya28w?d=KP=miXqjO^hPc@^< zJcGZ`sdIx9Z?#jGu@irJYv<-x-rCkK%U1q;a;G+#x02jNB=eV+bZ#o)sg-nXE#WIf zcIrm*lq0)LBl$|APJ)OBf22$ijg-$iA^AEX*c3BKLrkE+3wB`{Jal z{Ym#wE0+(f+&rvYKeKW_ROMn*<>p)E>R9C-uH*8)j+>i~>!&*IL04TqzUub3tFE74 zbwA?iV&&=P>FK)H)BP~d#g^yh$8&Y!xks*Xd2fxI>l)WxYup2qTs}&2^GY((n7C{+adS6u-EHDdyXCUumYdHl*ZsHL!w$Px zA9nLP?7Hu;d&sbh-LRYgu&eX1dz51PyNVqyicddLFFLP zZg6Y+&eo2DtxpfMawEv?+sGa6Riu_igU|hfP(dP5A|z%Aag1{$ay>$tsIu-Dk1xWwEMy zSmg{>WecmKh{c@HEBjLK{yn{W-|AI$=#`(>t9+eD7WLt>XEXC;3)3_*UHVWe)e0o$a~L?zxxUQ`Os3p3zg; z+EY=|!xX8P9aFzwu72+u^{S`p<>%BZAE;O4sxt=|Wv3bUYZ&*wXH-3Bl&3N(n;8{_ z45rYr?1ba}O2>Ow9ILt<%P%@s{_I$B(~&vUQWoEGzrN+(4=q(MTFTQ~Dxb7e6t^%Z zNoCQb`xT^nmq}GzQh6e&@)4;bpTratmBkj_uPeHDy{M|Ys64Hx@^Mkc-6H0Acv(#N zeP;N*YvEP=@bcvF%Es`D+u_WSiL!)=`<#h;ITKZV6Xlr`m2DFhr4vl?n!#B-CSN=j z8#>=WYrdgDw87iahGyFg-r8okrNTh3!jPnGuwL8H=(52(mkll44RqWMO}Ga7Tto6o zgLNwn4HFGECK{UWHrTM+(D;$T=0}E>L<4Q2Au-=zQ@-I=nt?9O&{SYR5ExRH-q2c_ zYY=2uxn#8)>(p`$ z({5}`%QfGBW5fPjdP93_F7fV-O?Pv*hTYH&%QYRpK^V`aCJ3efF3iDq#LkaH_+%>z&qChExZDByaG-50s8zvvU0#Wq2ewA_>qZTjiu(!TfvBhA%sI=d-(AMiT2||_jQ#qunopM1 z8!W2}SjKrRuKL3`*4uHlU&OIh;%e;U>doTnLgP3KY^p!7VZCKj>tVxQWmEH+P5l;| zx&s%RvuYe!^(0nZIEyn+uljvG)_T2KH$C#)7<;&9Xt@ZR}uko$f>sxQ)TX)!(qu5hz z+r!fDsrBn&YxLAO_0*Gl>LPnMbJeTgQ)jJHuXR;tuTZbqrCx8SUKgm&Sk}?O3zNvEJCRF4&Q?yrp_)3u|*r?ZFoI z+LoFFE%lZybrCI``K0P?B$hU*)}6#&NvhdRswa}_Xe7?kqUs$*tW8C=K1FP`qMH3h z^;?VT!iqQw;nmjREZy*0uW+_deiW_kZ_LDM77-ni!f2^Kf%_VsBxaCr%cpE zO*kJH>&}^P^!j|_Ytcq3(ZmJYj8<(UE~_w7uOQCTHhM#wxa6|Yo0o};+>O?_6BW5e z8eHPsl}0O85*H^LsU{K^?lxMzo4EXu(b`AE`9!0YMB>tXBei^@0?kO7MpP0QX$tb^ zEX{j;Y5r@mc`C8_3wGqK+L6DkE>FEKf8M6NH#X%jxt{mt_54LXd24*~6}$5^y7T9% z<*iW5U!0bwnwGzCf8Of-`O6>Yt$mz7e{0^#t@%su=BeGyR|v~f4$D^>&(j>I%~1$? zU4iymOprnFdGKs8V2^yh5%%TCgq5;Bprn6x{u3=#5vspAC-;0zp^w}^V-!Sn0Y#5L)8pt^{8wTWC6jBz6 zXD`&{8wa+JM*Z&&1Y)JCPq&C=cUc!!Nd6Q{_Ngv;oD<4;QEs>#9xA^6sq}ea)52Bg zh5!6-^;TSz&3#u@(vy8UPM?cj-TQjT?uGxNuchPve;-FpRqf`?mC2Ufy~#dzNtwgC zD-D6@wfSFf-@k2>&r;3M4{fjI&iv!DS1#v%|8Y_@FGq=S_m2cx{3@^<3B+Ia;Buv} zrQ?6{aa4JJ+}D5nRc5af|A+l)-0XFtoc| Date: Wed, 6 Sep 2023 12:10:04 +0200 Subject: [PATCH 12/49] fix --- src/schemes/fluid/viscosity.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 60d27ef92..13e403264 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -64,14 +64,9 @@ end v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, -<<<<<<< HEAD:src/schemes/fluid/weakly_compressible_sph/viscosity.jl distance, sound_speed, m_a, m_b, rho_mean) - @unpack smoothing_length = particle_system -======= - distance, sound_speed, m_a, m_b) (; smoothing_length) = particle_system ->>>>>>> main:src/schemes/fluid/viscosity.jl v_a = current_velocity(v_particle_system, particle_system, particle) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) From 5d47233e206c38a6ca092eb29299700baad4bf02 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 6 Sep 2023 13:55:43 +0200 Subject: [PATCH 13/49] remove @unpack --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 03f32f9c9..f5d6928a8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -80,7 +80,7 @@ end TotalLagrangianSPHSystem}, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) - @unpack boundary_model = neighbor_system + (;boundary_model) = neighbor_system return boundary_particle_impact(particle, neighbor, boundary_model, v_particle_system, v_neighbor_system, From bab4ca8400b530ff3082079dd1e705901dd3dd19 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 25 Sep 2023 15:40:50 +0200 Subject: [PATCH 14/49] format --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 951ce0ee4..f5b0afc95 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -83,7 +83,7 @@ end TotalLagrangianSPHSystem}, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel) - (;boundary_model) = neighbor_system + (; boundary_model) = neighbor_system return boundary_particle_impact(particle, neighbor, boundary_model, v_particle_system, v_neighbor_system, From f96c6a02f0829ce0b83e07a0c936fa2f5f74f495 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:02:45 +0200 Subject: [PATCH 15/49] review comments --- .../fluid/weakly_compressible_sph/rhs.jl | 16 +-- src/visualization/write2vtk.jl | 101 +++++------------- 2 files changed, 35 insertions(+), 82 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index f5b0afc95..963b8d2cc 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -14,7 +14,7 @@ function interact!(dv, v_particle_system, u_particle_system, neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) # Example - # debug_array = zeros(eltype(particle_system), ndims(particle_system), nparticles(particle_system)) + debug_array = zeros(eltype(particle_system), ndims(particle_system), nparticles(particle_system)) # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, @@ -49,19 +49,21 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_mean) + for i in 1:ndims(particle_system) + dv[i, particle] += dv_pressure[i] + dv_viscosity[i] + # example + debug_array[i, particle] += dv_pressure[i] + end + continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, particle_system, neighbor_system, grad_kernel) - for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity[i] - # example - # debug_array[i, particle] += dv_pressure[i] - end end # example - # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, custom_value=Dict("debug"=>debug_array), prefix="debug") + periodic_box = neighborhood_search.periodic_box + trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug") return dv end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ca6a444e7..4e2089659 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -1,24 +1,24 @@ """ trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", custom_quantities...) -Converts Trixi simulation data to VTK format. +Convert Trixi simulation data to VTK format. # Arguments -- `vu_ode::ODESolution`: Solution of the Trixi ODE system. -- `semi::SemiDiscretization`: Semi-discretization of the TrixiParticles simulation. -- `t::Float64`: Current time of the simulation. -- `iter::Union{Nothing, Int}` (optional): Iteration number. Defaults to `nothing`. -- `output_directory::AbstractString` (optional): Output directory path. Defaults to `"out"`. -- `prefix::AbstractString` (optional): Prefix for output files. Defaults to an empty string. -- `custom_quantities...`: Additional custom quantities to include in the VTK output. - -# Details -This function converts TrixiParticles simulation data to VTK format. -It iterates over each system in the semi-discretization and calls the overloaded `trixi2vtk` function for each system. +- `vu_ode`: Solution of the TrixiParticles ODE system at one time step. This expects an `ArrayPartition` as returned in the examples as `sol`. +- `semi`: Semidiscretization of the TrixiParticles simulation. +- `t`: Current time of the simulation. + +# Keywords +- `iter`: Iteration number when multiple iterations are to be stored in separate files. +- `output_directory`: Output directory path. Defaults to `"out"`. +- `prefix`: Prefix for output files. Defaults to an empty string. +- `custom_quantities...`: Additional custom quantities to include in the VTK output. TODO. + + # Example ```julia -trixi2vtk(vu, semi, 0.0, iter=1, output_directory="output", prefix="solution", velocity=compute_velocity) +trixi2vtk(sol[end], semi, 0.0, iter=1, output_directory="output", prefix="solution", velocity=compute_velocity) """ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", custom_quantities...) @@ -44,34 +44,9 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix end end -""" - trixi2vtk(v, u, t, system; output_directory="out", prefix="", iter=nothing, - system_name=vtkname(system), custom_quantities...) - -Converts a single Trixi system data to VTK format. - -# Arguments -- `v::AbstractVector`: Vector of solution variables. -- `u::AbstractVector`: Vector of unknowns. -- `t::Float64`: Current time of the simulation. -- `system::TrixiSystem`: Trixi system to convert to VTK. -- `output_directory::AbstractString` (optional): Output directory path. Defaults to "out". -- `prefix::AbstractString` (optional): Prefix for output files. Defaults to an empty string. -- `iter::Union{Nothing, Int}` (optional): Iteration number. Defaults to nothing. -- `system_name::AbstractString` (optional): Name of the system for the output file. Defaults to the VTK name of the system. -- `custom_quantities...`: Additional custom quantities to include in the VTK output. - -# Details -This function converts a single Trixi system's data to VTK format. It creates the necessary VTK file, -writes the solution variables and unknowns to the file, and includes additional custom quantities if provided. - -# Example -trixi2vtk(v, u, 0.0, fluid_system, output_directory="output", prefix="solution", velocity=compute_velocity) - -""" +# Convert a single Trixi system data to VTK format. function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix="", - iter=nothing, system_name=vtkname(system), custom_value=nothing, - custom_quantities...) + iter=nothing, system_name=vtkname(system), custom_quantities...) mkpath(output_directory) # handle "_" on optional pre/postfix strings @@ -91,50 +66,26 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix # Store particle index vtk["index"] = eachparticle(system) - if custom_value !== nothing - for (key, value) in custom_value - if axes(value, 1) == ndims(system) - vtk[string(key)] = view(value, 1:ndims(system), :) - else - vtk[string(key)] = value - end - end - end - # Extract custom quantities for this system - for (key, func) in custom_quantities - value = func(v, u, t, system) + for (key, quantity) in custom_quantities + value = custom_quantity(quantity, v, u, t, system) if value !== nothing - vtk[string(key)] = func(v, u, t, system) + vtk[string(key)] = value end end end end -""" - trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates") - -Converts coordinate data to VTK format. - -# Arguments -- `coordinates::AbstractMatrix`: Matrix of coordinate data. -- `output_directory::AbstractString` (optional): Output directory path. Defaults to `"out"`. -- `prefix::AbstractString` (optional): Prefix for the output file. Defaults to an empty string. -- `filename::AbstractString` (optional): Name of the output file. Defaults to `"coordinates"`. - -# Details -This function converts coordinate data to VTK format. -It creates a VTK file with the specified filename and saves the coordinate data as points in the VTK file. -Each coordinate is treated as a vertex cell in the VTK file. +function custom_quantity(quantity::AbstractArray, v, u, t, system) + return quantity +end -# Returns -- `file::AbstractString`: Path to the generated VTK file. +function custom_quantity(quantity, v, u, t, system) + # Assume `quantity` is a function of `v`, `u`, `t`, and `system` + return quantity(v, u, t, system) +end -# Example -```julia -coordinates = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] -vtk_file = trixi2vtk(coordinates, output_directory="output", prefix="data", filename="coords") -""" +# Convert coordinate data to VTK format. function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates", custom_quantities...) mkpath(output_directory) From a09e39acd100e1a633697b96175426767db71d49 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:20:30 +0200 Subject: [PATCH 16/49] rename --- .../fluid/weakly_compressible_sph/rhs.jl | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 963b8d2cc..b1c15ba6c 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -13,8 +13,10 @@ function interact!(dv, v_particle_system, u_particle_system, system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # Example - debug_array = zeros(eltype(particle_system), ndims(particle_system), nparticles(particle_system)) + # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code + # and the two other lines below that are marked as "debug example". + # + # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, @@ -37,7 +39,7 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - dv_pressure = calc_pressure_eq(pressure_correction, m_b, particle, + dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, @@ -51,8 +53,8 @@ function interact!(dv, v_particle_system, u_particle_system, for i in 1:ndims(particle_system) dv[i, particle] += dv_pressure[i] + dv_viscosity[i] - # example - debug_array[i, particle] += dv_pressure[i] + # Debug example + # debug_array[i, particle] += dv_pressure[i] end continuity_equation!(dv, density_calculator, @@ -62,13 +64,13 @@ function interact!(dv, v_particle_system, u_particle_system, end # example - periodic_box = neighborhood_search.periodic_box - trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug") + # periodic_box = neighborhood_search.periodic_box + # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug") return dv end -@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, +@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, @@ -79,7 +81,7 @@ end pressure_correction end -@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, +@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}, From 7ad06435ab9b715483ffe70046f03e27f3311d05 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:33:12 +0200 Subject: [PATCH 17/49] format --- .../fluid/weakly_compressible_sph/rhs.jl | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index b1c15ba6c..665afbb73 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -40,10 +40,10 @@ function interact!(dv, v_particle_system, u_particle_system, m_b = hydrodynamic_mass(neighbor_system, neighbor) dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, - particle_system, v_particle_system, - neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel) + particle_system, v_particle_system, + neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, + pos_diff, distance, grad_kernel) dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, @@ -61,7 +61,6 @@ function interact!(dv, v_particle_system, u_particle_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, particle_system, neighbor_system, grad_kernel) - end # example # periodic_box = neighborhood_search.periodic_box @@ -71,10 +70,10 @@ function interact!(dv, v_particle_system, u_particle_system, end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, - neighbor_system::WeaklyCompressibleSPHSystem, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel) + v_particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel) return (-m_b * (particle_system.pressure[particle] / rho_a^2 + neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) * @@ -82,11 +81,11 @@ end end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, - neighbor_system::Union{BoundarySPHSystem, - TotalLagrangianSPHSystem}, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel) + v_particle_system, neighbor, + neighbor_system::Union{BoundarySPHSystem, + TotalLagrangianSPHSystem}, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel) (; boundary_model) = neighbor_system return boundary_particle_impact(particle, neighbor, boundary_model, From 3add63907896b1d0a535bcf04013816266d4270e Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:35:55 +0200 Subject: [PATCH 18/49] revert --- src/visualization/write2vtk.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 4e2089659..e673a1a0d 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -85,7 +85,20 @@ function custom_quantity(quantity, v, u, t, system) return quantity(v, u, t, system) end -# Convert coordinate data to VTK format. +""" + trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates") + +Convert coordinate data to VTK format. + +# Arguments +- `coordinates`: Coordinates to be saved. +- `output_directory` (optional): Output directory path. Defaults to `"out"`. +- `prefix` (optional): Prefix for the output file. Defaults to an empty string. +- `filename` (optional): Name of the output file. Defaults to `"coordinates"`. + +# Returns +- `file::AbstractString`: Path to the generated VTK file. +""" function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates", custom_quantities...) mkpath(output_directory) From 960ba248ee5f75dbf4a56cc7b52cc8b2fd15e7f8 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:38:26 +0200 Subject: [PATCH 19/49] rename --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 2 +- src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- src/schemes/solid/total_lagrangian_sph/rhs.jl | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index efe5b3646..6356c022c 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -103,7 +103,7 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) print(io, ")") end -@inline function boundary_particle_impact(particle, boundary_particle, +@inline function calc_bnd_pressure(particle, boundary_particle, boundary_model::BoundaryModelDummyParticles, v_particle_system, v_boundary_system, particle_system, boundary_system, diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index e049b5f41..4f864aa5e 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -72,7 +72,7 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end -@inline function boundary_particle_impact(particle, boundary_particle, +@inline function calc_bnd_pressure(particle, boundary_particle, boundary_model::BoundaryModelMonaghanKajtar, v_particle_system, v_boundary_system, particle_system, boundary_system, diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 665afbb73..be0bc711d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -88,7 +88,7 @@ end grad_kernel) (; boundary_model) = neighbor_system - return boundary_particle_impact(particle, neighbor, boundary_model, + return calc_bnd_pressure(particle, neighbor, boundary_model, v_particle_system, v_neighbor_system, particle_system, neighbor_system, pos_diff, distance, m_b) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index f7078cbae..6b24f3a1d 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -95,7 +95,7 @@ function interact!(dv, v_particle_system, u_particle_system, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces - dv_boundary = boundary_particle_impact(neighbor, particle, boundary_model, + dv_boundary = calc_bnd_pressure(neighbor, particle, boundary_model, v_neighbor_system, v_particle_system, neighbor_system, particle_system, pos_diff, distance, m_a) From 56d4d8546c64002f171280615adc954e6b252990 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:38:42 +0200 Subject: [PATCH 20/49] format --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 8 ++++---- src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl | 8 ++++---- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 6 +++--- src/schemes/solid/total_lagrangian_sph/rhs.jl | 6 +++--- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 6356c022c..bbb60377e 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -104,10 +104,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) end @inline function calc_bnd_pressure(particle, boundary_particle, - boundary_model::BoundaryModelDummyParticles, - v_particle_system, v_boundary_system, - particle_system, boundary_system, - pos_diff, distance, m_b) + boundary_model::BoundaryModelDummyParticles, + v_particle_system, v_boundary_system, + particle_system, boundary_system, + pos_diff, distance, m_b) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index 4f864aa5e..fa2d5c49e 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -73,10 +73,10 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) end @inline function calc_bnd_pressure(particle, boundary_particle, - boundary_model::BoundaryModelMonaghanKajtar, - v_particle_system, v_boundary_system, - particle_system, boundary_system, - pos_diff, distance, m_b) + boundary_model::BoundaryModelMonaghanKajtar, + v_particle_system, v_boundary_system, + particle_system, boundary_system, + pos_diff, distance, m_b) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index be0bc711d..20acab45e 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -89,9 +89,9 @@ end (; boundary_model) = neighbor_system return calc_bnd_pressure(particle, neighbor, boundary_model, - v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - pos_diff, distance, m_b) + v_particle_system, v_neighbor_system, + particle_system, neighbor_system, + pos_diff, distance, m_b) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 6b24f3a1d..eb1926ce8 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -96,9 +96,9 @@ function interact!(dv, v_particle_system, u_particle_system, # Boundary forces dv_boundary = calc_bnd_pressure(neighbor, particle, boundary_model, - v_neighbor_system, v_particle_system, - neighbor_system, particle_system, - pos_diff, distance, m_a) + v_neighbor_system, v_particle_system, + neighbor_system, particle_system, + pos_diff, distance, m_a) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) From cc81ba93ad55c32ec50412294566e5f1a70eaea4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 12:39:55 +0200 Subject: [PATCH 21/49] add todo --- src/visualization/write2vtk.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index e673a1a0d..2349cbe3a 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -18,7 +18,9 @@ Convert Trixi simulation data to VTK format. # Example ```julia -trixi2vtk(sol[end], semi, 0.0, iter=1, output_directory="output", prefix="solution", velocity=compute_velocity) +trixi2vtk(sol[end], semi, 0.0, iter=1, output_directory="output", prefix="solution") + +TODO: example for custom_quantities """ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", custom_quantities...) From 6d9166971b5e6448db169ca783bf58661bd44de5 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 13:03:56 +0200 Subject: [PATCH 22/49] allow to save multiple files by using a local iter --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 10 ++++++---- src/util.jl | 10 ++++++++++ 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 20acab45e..ca4725d1f 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -1,3 +1,5 @@ +using .IteratorModule + # Computes the forces that particles in `particle_system` experience from particles # in `neighbor_system` and updates `dv` accordingly. # It takes into account pressure forces, viscosity, and for `ContinuityDensity` updates the density @@ -15,8 +17,7 @@ function interact!(dv, v_particle_system, u_particle_system, # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code # and the two other lines below that are marked as "debug example". - # - # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) + debug_array = zeros(ndims(particle_system), nparticles(particle_system)) # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, @@ -62,9 +63,10 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, particle_system, neighbor_system, grad_kernel) end - # example + # Debug example # periodic_box = neighborhood_search.periodic_box - # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug") + # Note: this saves a file every step within the integration scheme! + # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter()) return dv end diff --git a/src/util.jl b/src/util.jl index 797997044..583964835 100644 --- a/src/util.jl +++ b/src/util.jl @@ -299,3 +299,13 @@ macro autoinfiltrate(condition=true) lnn, esc(condition)) end + + +module IteratorModule + i = 0 + export iter + + function iter() + global i += 1 + end +end From efb1707d4de8f81cdb9d6a3d8e3ca110bc5b8d6b Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 13:05:54 +0200 Subject: [PATCH 23/49] add todo --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index ca4725d1f..d97c734d6 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -66,6 +66,7 @@ function interact!(dv, v_particle_system, u_particle_system, # Debug example # periodic_box = neighborhood_search.periodic_box # Note: this saves a file every step within the integration scheme! + # TODO: this call should use public API this requires some additional changes to simplify the calls # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter()) return dv From 2a9b4f64ffaa02105ac56dfba571c5eeaf0904c1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 13:06:03 +0200 Subject: [PATCH 24/49] format --- src/util.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/util.jl b/src/util.jl index 583964835..856c2aee1 100644 --- a/src/util.jl +++ b/src/util.jl @@ -300,12 +300,11 @@ macro autoinfiltrate(condition=true) esc(condition)) end - module IteratorModule - i = 0 - export iter +i = 0 +export iter - function iter() - global i += 1 - end +function iter() + global i += 1 +end end From 82c5705f79cb2ae58e74851bd7058c5278aa0850 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 14:45:43 +0200 Subject: [PATCH 25/49] rename functions --- .../boundary/dummy_particles/dummy_particles.jl | 11 ++++++----- .../boundary/monaghan_kajtar/monaghan_kajtar.jl | 10 +++++----- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 8 ++++---- src/schemes/solid/total_lagrangian_sph/rhs.jl | 9 +++++---- 4 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index bbb60377e..ae00b45da 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -103,11 +103,12 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) print(io, ")") end -@inline function calc_bnd_pressure(particle, boundary_particle, - boundary_model::BoundaryModelDummyParticles, - v_particle_system, v_boundary_system, - particle_system, boundary_system, - pos_diff, distance, m_b) +@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, + v_particle_system, boundary_particle, + boundary_system, + v_boundary_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, + boundary_model::BoundaryModelDummyParticles) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index fa2d5c49e..53087421e 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -72,11 +72,11 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar) print(io, ")") end -@inline function calc_bnd_pressure(particle, boundary_particle, - boundary_model::BoundaryModelMonaghanKajtar, - v_particle_system, v_boundary_system, - particle_system, boundary_system, - pos_diff, distance, m_b) +@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, + v_particle_system, neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, + boundary_model::BoundaryModelMonaghanKajtar) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index d97c734d6..b25ea04c1 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -91,10 +91,10 @@ end grad_kernel) (; boundary_model) = neighbor_system - return calc_bnd_pressure(particle, neighbor, boundary_model, - v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - pos_diff, distance, m_b) + return pressure_acceleration(pressure_correction, m_b, particle, particle_system, + v_particle_system, neighbor, neighbor_system, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, boundary_model) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index eb1926ce8..4f921917e 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -95,10 +95,11 @@ function interact!(dv, v_particle_system, u_particle_system, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces - dv_boundary = calc_bnd_pressure(neighbor, particle, boundary_model, - v_neighbor_system, v_particle_system, - neighbor_system, particle_system, - pos_diff, distance, m_a) + dv_boundary = pressure_acceleration(pressure_correction, m_b, particle, + particle_system, v_particle_system, neighbor, + neighbor_system, v_neighbor_system, rho_a, + rho_b, pos_diff, distance, + grad_kernel, boundary_model) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) From 7d7596748db162dee76711c900bcd71b702c29c5 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 15:27:53 +0200 Subject: [PATCH 26/49] fix --- src/schemes/solid/total_lagrangian_sph/rhs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4f921917e..73a9e62cd 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -95,7 +95,7 @@ function interact!(dv, v_particle_system, u_particle_system, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces - dv_boundary = pressure_acceleration(pressure_correction, m_b, particle, + dv_boundary = pressure_acceleration(1.0, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, From 9a231ade2f1389c9faa3ce2f490b66bb364f0ff4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 15:56:43 +0200 Subject: [PATCH 27/49] fix --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 2 -- src/schemes/solid/total_lagrangian_sph/rhs.jl | 2 ++ 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index ae00b45da..4dc6e82f9 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -112,8 +112,6 @@ end rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) - return -m_b * (particle_system.pressure[particle] / rho_a^2 + boundary_model.pressure[boundary_particle] / rho_b^2) * diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 73a9e62cd..ee351ec83 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -88,6 +88,8 @@ function interact!(dv, v_particle_system, u_particle_system, rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = (rho_a + rho_b) / 2 + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) + # use `m_a` to get the same viscosity as for the fluid-solid direction. dv_viscosity = viscosity(neighbor_system, particle_system, v_neighbor_system, v_particle_system, From c2237b069491fa963f3fb9996831342824c1c93c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 26 Sep 2023 21:06:17 +0200 Subject: [PATCH 28/49] mistake --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index b25ea04c1..15d2f0b32 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -17,7 +17,7 @@ function interact!(dv, v_particle_system, u_particle_system, # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code # and the two other lines below that are marked as "debug example". - debug_array = zeros(ndims(particle_system), nparticles(particle_system)) + # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, From fbe68d3b1f3b2f204c308cfff51f3733c7e4a715 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 01:24:56 +0200 Subject: [PATCH 29/49] fix --- src/schemes/solid/total_lagrangian_sph/rhs.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index ee351ec83..67520b1fb 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -97,9 +97,10 @@ function interact!(dv, v_particle_system, u_particle_system, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces - dv_boundary = pressure_acceleration(1.0, m_b, particle, - particle_system, v_particle_system, neighbor, - neighbor_system, v_neighbor_system, rho_a, + # Note: neighbor and particle are switched in this call + dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system, + v_neighbor_system, particle, + particle_system, v_particle_system, rho_a, rho_b, pos_diff, distance, grad_kernel, boundary_model) dv_particle = dv_boundary + dv_viscosity From 983e1284c51e6ed3fb666b3e7668239069dfd0ae Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 14:35:18 +0200 Subject: [PATCH 30/49] rename to viscosity_model --- src/schemes/boundary/system.jl | 2 +- src/schemes/fluid/entropically_damped_sph/rhs.jl | 6 +++--- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/system.jl | 2 +- src/schemes/solid/total_lagrangian_sph/system.jl | 2 +- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 1185fe10a..7c512d854 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -293,6 +293,6 @@ function restart_with!(system, model, ::ContinuityDensity, v, u) return system end -function system_viscosity(system::BoundarySPHSystem) +function viscosity_model(system::BoundarySPHSystem) return system.boundary_model.viscosity end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index cac6f3ed6..808b33187 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -4,7 +4,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::EntropicallyDampedSPHSystem, neighbor_system) (; sound_speed) = particle_system - viscosity = viscosity_function(neighbor_system) + viscosity = viscosity_model(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -80,5 +80,5 @@ end return dv end -@inline viscosity_function(system) = system.viscosity -@inline viscosity_function(system::BoundarySPHSystem) = system.boundary_model.viscosity +@inline viscosity_model(system) = system.viscosity +@inline viscosity_model(system::BoundarySPHSystem) = system.boundary_model.viscosity diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 15d2f0b32..14ad4a91e 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -11,7 +11,7 @@ function interact!(dv, v_particle_system, u_particle_system, (; density_calculator, state_equation, correction, smoothing_length) = particle_system (; sound_speed) = state_equation - viscosity = system_viscosity(neighbor_system) + viscosity = viscosity_model(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 2410c9fb7..54fea868b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -297,7 +297,7 @@ function restart_with!(system, ::ContinuityDensity, v, u) return system end -function system_viscosity(system::WeaklyCompressibleSPHSystem) +function viscosity_model(system::WeaklyCompressibleSPHSystem) return system.viscosity end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 6f2484bf4..5f44a430a 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -464,6 +464,6 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u) restart_with!(system, system.boundary_model, v, u) end -function system_viscosity(system::TotalLagrangianSPHSystem) +function viscosity_model(system::TotalLagrangianSPHSystem) return system.boundary_model.viscosity end From 36c82e42f16884fa698400f749b58a8fd81497ca Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 15:57:23 +0200 Subject: [PATCH 31/49] switch param call order --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 9 ++++----- src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl | 4 ++-- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 4 ++-- src/schemes/solid/total_lagrangian_sph/rhs.jl | 4 ++-- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 4dc6e82f9..ccf9775e3 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -104,11 +104,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, - boundary_system, - v_boundary_system, rho_a, rho_b, pos_diff, distance, - grad_kernel, - boundary_model::BoundaryModelDummyParticles) + v_particle_system, boundary_particle, boundary_system, + v_boundary_system, + boundary_model::BoundaryModelDummyParticles, rho_a, + rho_b, pos_diff, distance, grad_kernel) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index 53087421e..22f7a1b15 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -74,9 +74,9 @@ end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, + boundary_model::BoundaryModelMonaghanKajtar, v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel, - boundary_model::BoundaryModelMonaghanKajtar) + grad_kernel) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 14ad4a91e..57686ad62 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -93,8 +93,8 @@ end return pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel, boundary_model) + v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, distance, + grad_kernel) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 67520b1fb..ffd7aa062 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -100,9 +100,9 @@ function interact!(dv, v_particle_system, u_particle_system, # Note: neighbor and particle are switched in this call dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system, v_neighbor_system, particle, - particle_system, v_particle_system, rho_a, + particle_system, v_particle_system, boundary_model, rho_a, rho_b, pos_diff, distance, - grad_kernel, boundary_model) + grad_kernel) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) From 126a1cf9084864a8c27f4893b5005879d2e96317 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 16:03:16 +0200 Subject: [PATCH 32/49] review comments --- .../fluid/weakly_compressible_sph/rhs.jl | 9 ++++----- src/schemes/solid/total_lagrangian_sph/rhs.jl | 2 +- src/util.jl | 9 --------- src/visualization/write2vtk.jl | 19 +++++++++---------- 4 files changed, 14 insertions(+), 25 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 57686ad62..5559dfc3d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -1,5 +1,3 @@ -using .IteratorModule - # Computes the forces that particles in `particle_system` experience from particles # in `neighbor_system` and updates `dv` accordingly. # It takes into account pressure forces, viscosity, and for `ContinuityDensity` updates the density @@ -65,9 +63,10 @@ function interact!(dv, v_particle_system, u_particle_system, end # Debug example # periodic_box = neighborhood_search.periodic_box - # Note: this saves a file every step within the integration scheme! - # TODO: this call should use public API this requires some additional changes to simplify the calls - # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter()) + # Note: this saves a file in every stage of the integrator + # if !@isdefined iter; iter = 0; end + # TODO: This call should use public API. This requires some additional changes to simplify the calls. + # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) return dv end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index ffd7aa062..fdd019b41 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -97,7 +97,7 @@ function interact!(dv, v_particle_system, u_particle_system, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces - # Note: neighbor and particle are switched in this call + # Note: neighbor and particle are switched in this call and `pressure_correction` is set to `1.0` (no correction) dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system, v_neighbor_system, particle, particle_system, v_particle_system, boundary_model, rho_a, diff --git a/src/util.jl b/src/util.jl index 856c2aee1..797997044 100644 --- a/src/util.jl +++ b/src/util.jl @@ -299,12 +299,3 @@ macro autoinfiltrate(condition=true) lnn, esc(condition)) end - -module IteratorModule -i = 0 -export iter - -function iter() - global i += 1 -end -end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 2349cbe3a..a63770e84 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -5,17 +5,16 @@ Convert Trixi simulation data to VTK format. # Arguments - `vu_ode`: Solution of the TrixiParticles ODE system at one time step. This expects an `ArrayPartition` as returned in the examples as `sol`. -- `semi`: Semidiscretization of the TrixiParticles simulation. -- `t`: Current time of the simulation. +- `semi`: Semidiscretization of the TrixiParticles simulation. +- `t`: Current time of the simulation. # Keywords -- `iter`: Iteration number when multiple iterations are to be stored in separate files. -- `output_directory`: Output directory path. Defaults to `"out"`. -- `prefix`: Prefix for output files. Defaults to an empty string. +- `iter`: Iteration number when multiple iterations are to be stored in separate files. +- `output_directory`: Output directory path. Defaults to `"out"`. +- `prefix`: Prefix for output files. Defaults to an empty string. - `custom_quantities...`: Additional custom quantities to include in the VTK output. TODO. - # Example ```julia trixi2vtk(sol[end], semi, 0.0, iter=1, output_directory="output", prefix="solution") @@ -46,7 +45,7 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix end end -# Convert a single Trixi system data to VTK format. +# Convert data for a single TrixiParticle system to VTK format function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix="", iter=nothing, system_name=vtkname(system), custom_quantities...) mkpath(output_directory) @@ -93,10 +92,10 @@ end Convert coordinate data to VTK format. # Arguments -- `coordinates`: Coordinates to be saved. +- `coordinates`: Coordinates to be saved. - `output_directory` (optional): Output directory path. Defaults to `"out"`. -- `prefix` (optional): Prefix for the output file. Defaults to an empty string. -- `filename` (optional): Name of the output file. Defaults to `"coordinates"`. +- `prefix` (optional): Prefix for the output file. Defaults to an empty string. +- `filename` (optional): Name of the output file. Defaults to `"coordinates"`. # Returns - `file::AbstractString`: Path to the generated VTK file. From 530d504ebb170fde35919545f18127274cfdd47f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 16:03:28 +0200 Subject: [PATCH 33/49] format --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 5 +++-- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 3 ++- src/schemes/solid/total_lagrangian_sph/rhs.jl | 3 ++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index ccf9775e3..c135cfb04 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -104,9 +104,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles) end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, boundary_particle, boundary_system, + v_particle_system, boundary_particle, + boundary_system, v_boundary_system, - boundary_model::BoundaryModelDummyParticles, rho_a, + boundary_model::BoundaryModelDummyParticles, rho_a, rho_b, pos_diff, distance, grad_kernel) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 5559dfc3d..de03fdc07 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -92,7 +92,8 @@ end return pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, - v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, distance, + v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, + distance, grad_kernel) end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index fdd019b41..c220ac370 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -100,7 +100,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Note: neighbor and particle are switched in this call and `pressure_correction` is set to `1.0` (no correction) dv_boundary = pressure_acceleration(1.0, m_b, neighbor, neighbor_system, v_neighbor_system, particle, - particle_system, v_particle_system, boundary_model, rho_a, + particle_system, v_particle_system, + boundary_model, rho_a, rho_b, pos_diff, distance, grad_kernel) dv_particle = dv_boundary + dv_viscosity From 38657b3c285e5c77f64c9300f5444ba18a2ef5c3 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 27 Sep 2023 16:50:59 +0200 Subject: [PATCH 34/49] fix --- src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl | 6 +++--- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 3 +-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index 22f7a1b15..f45b1376d 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -74,9 +74,9 @@ end @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, - boundary_model::BoundaryModelMonaghanKajtar, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel) + v_neighbor_system, + boundary_model::BoundaryModelMonaghanKajtar, rho_a, + rho_b, pos_diff, distance, grad_kernel) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index de03fdc07..049bf5e04 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -93,8 +93,7 @@ end return pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, - distance, - grad_kernel) + distance, grad_kernel) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, From 18b1b36ba0c75b119afc274a1db85df1de01daee Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 14:11:16 +0200 Subject: [PATCH 35/49] move viscosity_model --- src/schemes/fluid/entropically_damped_sph/rhs.jl | 3 --- src/schemes/fluid/fluid.jl | 2 ++ src/schemes/fluid/weakly_compressible_sph/system.jl | 3 --- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 808b33187..b5e5266ec 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -79,6 +79,3 @@ end return dv end - -@inline viscosity_model(system) = system.viscosity -@inline viscosity_model(system::BoundarySPHSystem) = system.boundary_model.viscosity diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 258cb19f3..844e0a91d 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -17,6 +17,8 @@ function write_u0!(u0, system::FluidSystem) return u0 end +@inline viscosity_model(system::FluidSystem) = system.viscosity + include("viscosity.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 54fea868b..cd6a0c11a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -297,9 +297,6 @@ function restart_with!(system, ::ContinuityDensity, v, u) return system end -function viscosity_model(system::WeaklyCompressibleSPHSystem) - return system.viscosity -end @inline function smoothing_kernel_grad(system::WeaklyCompressibleSPHSystem, pos_diff, distance, particle) From f840bfe5ff4c10b8801b39604e7cf64010375c1b Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 14:33:01 +0200 Subject: [PATCH 36/49] fix --- src/schemes/fluid/entropically_damped_sph/rhs.jl | 3 ++- src/schemes/fluid/viscosity.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index b5e5266ec..aac412dc1 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -18,6 +18,7 @@ function interact!(dv, v_particle_system, u_particle_system, rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = 0.5 * (rho_a + rho_b) p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) @@ -39,7 +40,7 @@ function interact!(dv, v_particle_system, u_particle_system, dv_viscosity = viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b) + sound_speed, m_a, m_b, rho_mean) for i in 1:ndims(particle_system) dv[i, particle] += dv_pressure[i] + dv_viscosity[i] diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 13e403264..490aec9a1 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -153,7 +153,7 @@ end @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b) + distance, sound_speed, m_a, m_b, rho_mean) (; epsilon, nu) = viscosity (; smoothing_length) = particle_system diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 049bf5e04..32d639d80 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -26,7 +26,7 @@ function interact!(dv, v_particle_system, u_particle_system, rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 + rho_mean = 0.5 * (rho_a + rho_b) # Determine correction values viscosity_correction, pressure_correction = free_surface_correction(correction, From 82903edc516c2916a7719a87124b8e1402d780b9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 14:33:17 +0200 Subject: [PATCH 37/49] format --- src/schemes/fluid/weakly_compressible_sph/system.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index cd6a0c11a..b9390082c 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -297,7 +297,6 @@ function restart_with!(system, ::ContinuityDensity, v, u) return system end - @inline function smoothing_kernel_grad(system::WeaklyCompressibleSPHSystem, pos_diff, distance, particle) return corrected_kernel_grad(system.smoothing_kernel, pos_diff, distance, From 7e072e04163aabe2f02a64afa0d516ee77425ac9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 15:24:49 +0200 Subject: [PATCH 38/49] add comment --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 32d639d80..e1dcbdf7d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -109,6 +109,7 @@ end return dv end +# 'SummationDensity' is used so density is recalculated in wcsph/system.jl:compute_density!() @inline function continuity_equation!(dv, density_calculator::SummationDensity, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, From 407951b71d80e0d1cbec20effd413bfdb92a3e24 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 17:00:49 +0200 Subject: [PATCH 39/49] fix --- .../boundary/dummy_particles/dummy_particles.jl | 2 +- .../boundary/monaghan_kajtar/monaghan_kajtar.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 10 +++++----- src/schemes/solid/total_lagrangian_sph/rhs.jl | 4 ++-- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index c135cfb04..35f7c6d01 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -108,7 +108,7 @@ end boundary_system, v_boundary_system, boundary_model::BoundaryModelDummyParticles, rho_a, - rho_b, pos_diff, distance, grad_kernel) + rho_b, pos_diff, distance, grad_kernel, density_calculator) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index f45b1376d..162f822eb 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -76,7 +76,7 @@ end v_particle_system, neighbor, neighbor_system, v_neighbor_system, boundary_model::BoundaryModelMonaghanKajtar, rho_a, - rho_b, pos_diff, distance, grad_kernel) + rho_b, pos_diff, distance, grad_kernel, density_calculator) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index d7f0b1c3e..5412633f5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -41,8 +41,8 @@ function interact!(dv, v_particle_system, u_particle_system, dv_pressure = pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, - v_neighbor_system, rho_a, rho_b, - pos_diff, distance, grad_kernel) + v_neighbor_system, rho_a, rho_b, pos_diff, + distance, grad_kernel, density_calculator) dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, @@ -74,7 +74,7 @@ end # As shown in "Variational and momentum preservation aspects of Smooth # Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation # this form has to be used with ContinuityDensity -@inline function calc_pressure_eq(pressure_correction, m_b, particle, particle_system, +@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, v_neighbor_system, rho_a, rho_b, pos_diff, distance, @@ -104,13 +104,13 @@ end neighbor_system::Union{BoundarySPHSystem, TotalLagrangianSPHSystem}, v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel) + grad_kernel, density_calculator) (; boundary_model) = neighbor_system return pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system, v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff, - distance, grad_kernel) + distance, grad_kernel, density_calculator) end @inline function continuity_equation!(dv, density_calculator::ContinuityDensity, diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index c220ac370..163388519 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -61,7 +61,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem) (; boundary_model) = particle_system - (; state_equation, viscosity) = neighbor_system + (; density_calculator, state_equation, viscosity) = neighbor_system (; sound_speed) = state_equation system_coords = current_coordinates(u_particle_system, particle_system) @@ -103,7 +103,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, v_particle_system, boundary_model, rho_a, rho_b, pos_diff, distance, - grad_kernel) + grad_kernel, density_calculator) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) From b31f3a1e72db9aa118a73b94bc98e163481e3bc9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 17:01:01 +0200 Subject: [PATCH 40/49] format --- src/schemes/boundary/dummy_particles/dummy_particles.jl | 3 ++- src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl | 3 ++- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 8 ++++---- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 35f7c6d01..f9c2b990a 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -108,7 +108,8 @@ end boundary_system, v_boundary_system, boundary_model::BoundaryModelDummyParticles, rho_a, - rho_b, pos_diff, distance, grad_kernel, density_calculator) + rho_b, pos_diff, distance, grad_kernel, + density_calculator) rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle) diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index 162f822eb..a3795b39e 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl @@ -76,7 +76,8 @@ end v_particle_system, neighbor, neighbor_system, v_neighbor_system, boundary_model::BoundaryModelMonaghanKajtar, rho_a, - rho_b, pos_diff, distance, grad_kernel, density_calculator) + rho_b, pos_diff, distance, grad_kernel, + density_calculator) (; smoothing_length) = particle_system (; K, beta, boundary_particle_spacing) = boundary_model diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 5412633f5..413de32e1 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -75,10 +75,10 @@ end # Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation # this form has to be used with ContinuityDensity @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, - v_particle_system, neighbor, - neighbor_system::WeaklyCompressibleSPHSystem, - v_neighbor_system, rho_a, rho_b, pos_diff, distance, - grad_kernel, ::ContinuityDensity) + v_particle_system, neighbor, + neighbor_system::WeaklyCompressibleSPHSystem, + v_neighbor_system, rho_a, rho_b, pos_diff, distance, + grad_kernel, ::ContinuityDensity) return (-m_b * (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) / (rho_a * rho_b) * grad_kernel) * From 0d3ca812e661119e0c3346743c34665be7851a9d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 17:25:17 +0200 Subject: [PATCH 41/49] add error check --- src/general/semidiscretization.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index bbedb3bc9..2af5513ad 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -31,6 +31,16 @@ struct Semidiscretization{S, RU, RV, NS, DC} ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) + for sys in systems + if sys isa BoundarySPHSystem + for neighbor in systems + if neighbor isa WeaklyCompressibleSPHSystem && sys.boundary_model isa BoundaryModelDummyParticles && isnothing(sys.boundary_model.state_equation) + throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) + end + end + end + end + # Create (and initialize) a tuple of n neighborhood searches for each of the n systems # We will need one neighborhood search for each pair of systems. searches = Tuple(Tuple(create_neighborhood_search(system, neighbor, From cd9a8ad9821edcf9f7f3ab34d0843222ff4624a8 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 28 Sep 2023 17:25:28 +0200 Subject: [PATCH 42/49] format --- src/general/semidiscretization.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 2af5513ad..b6a6b126b 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -34,7 +34,9 @@ struct Semidiscretization{S, RU, RV, NS, DC} for sys in systems if sys isa BoundarySPHSystem for neighbor in systems - if neighbor isa WeaklyCompressibleSPHSystem && sys.boundary_model isa BoundaryModelDummyParticles && isnothing(sys.boundary_model.state_equation) + if neighbor isa WeaklyCompressibleSPHSystem && + sys.boundary_model isa BoundaryModelDummyParticles && + isnothing(sys.boundary_model.state_equation) throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) end end From 74c9b8a47233ca394cbaf7393ae36da5ae0b46e3 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Sep 2023 16:41:48 +0200 Subject: [PATCH 43/49] fix merge --- .../fluid/weakly_compressible_sph/rhs.jl | 25 +------------------ 1 file changed, 1 insertion(+), 24 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0ae68df78..cd086d2d8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,21 +6,13 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) - (; density_calculator, state_equation, correction, smoothing_length) = particle_system - neighbor_system) (; density_calculator, state_equation, correction, smoothing_length) = particle_system (; sound_speed) = state_equation - viscosity = viscosity_model(neighbor_system) viscosity = viscosity_model(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code - # and the two other lines below that are marked as "debug example". - # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) - neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # In order to visualize quantities like pressure forces or viscosity forces, uncomment the following code # and the two other lines below that are marked as "debug example". # debug_array = zeros(ndims(particle_system), nparticles(particle_system)) @@ -36,7 +28,6 @@ function interact!(dv, v_particle_system, u_particle_system, rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = 0.5 * (rho_a + rho_b) - rho_mean = 0.5 * (rho_a + rho_b) # Determine correction values viscosity_correction, pressure_correction = free_surface_correction(correction, @@ -45,8 +36,6 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) @@ -56,21 +45,17 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) - dv_viscosity = viscosity_correction * - viscosity(particle_system, neighbor_system, + dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_mean) - sound_speed, m_a, m_b, rho_mean) for i in 1:ndims(particle_system) dv[i, particle] += dv_pressure[i] + dv_viscosity[i] # Debug example # debug_array[i, particle] += dv_pressure[i] - # Debug example - # debug_array[i, particle] += dv_pressure[i] end continuity_equation!(dv, density_calculator, @@ -84,13 +69,6 @@ function interact!(dv, v_particle_system, u_particle_system, # if !@isdefined iter; iter = 0; end # TODO: This call should use public API. This requires some additional changes to simplify the calls. # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) - # Debug example - # periodic_box = neighborhood_search.periodic_box - # Note: this saves a file in every stage of the integrator - # if !@isdefined iter; iter = 0; end - # TODO: This call should use public API. This requires some additional changes to simplify the calls. - # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) - return dv end @@ -149,7 +127,6 @@ end return dv end -# 'SummationDensity' is used so density is recalculated in wcsph/system.jl:compute_density!() # 'SummationDensity' is used so density is recalculated in wcsph/system.jl:compute_density!() @inline function continuity_equation!(dv, density_calculator::SummationDensity, v_particle_system, v_neighbor_system, From e2d3ab7b0e7727350a1e8c3e45231fdcf9a76d2b Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Sep 2023 16:43:04 +0200 Subject: [PATCH 44/49] fix merge --- src/schemes/solid/total_lagrangian_sph/rhs.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index d1c678dcf..163388519 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -90,18 +90,11 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = (rho_a + rho_b) / 2 - - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) - # use `m_a` to get the same viscosity as for the fluid-solid direction. dv_viscosity = viscosity(neighbor_system, particle_system, v_neighbor_system, v_particle_system, neighbor, particle, pos_diff, distance, sound_speed, m_b, m_a, rho_mean) - pos_diff, distance, sound_speed, m_b, m_a, rho_mean) # Boundary forces # Note: neighbor and particle are switched in this call and `pressure_correction` is set to `1.0` (no correction) From 197b08ce7c9aba3fc90b863b7c1e79f74f137516 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Sep 2023 16:44:02 +0200 Subject: [PATCH 45/49] fix merge --- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index cd086d2d8..1b5021610 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -19,7 +19,6 @@ function interact!(dv, v_particle_system, u_particle_system, # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(particle_system, neighbor_system, - system_coords, neighbor_system_coords, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance # Only consider particles with a distance > 0. @@ -69,6 +68,7 @@ function interact!(dv, v_particle_system, u_particle_system, # if !@isdefined iter; iter = 0; end # TODO: This call should use public API. This requires some additional changes to simplify the calls. # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) + return dv end From c5d67b1550938e57bc9637a0ffe6d98f17823ea7 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Sep 2023 18:08:55 +0200 Subject: [PATCH 46/49] cleanup constructor --- src/general/semidiscretization.jl | 33 ++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index b6a6b126b..5618586d2 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -31,17 +31,7 @@ struct Semidiscretization{S, RU, RV, NS, DC} ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) - for sys in systems - if sys isa BoundarySPHSystem - for neighbor in systems - if neighbor isa WeaklyCompressibleSPHSystem && - sys.boundary_model isa BoundaryModelDummyParticles && - isnothing(sys.boundary_model.state_equation) - throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) - end - end - end - end + check_configuration(systems) # Create (and initialize) a tuple of n neighborhood searches for each of the n systems # We will need one neighborhood search for each pair of systems. @@ -491,3 +481,24 @@ function nhs_coords(system::BoundarySPHSystem, # Don't update return nothing end + +function check_configuration(systems) + for sys in systems + check_configuration(sys, systems) + end + +end + +check_configuration(sys, systems) = systems + + +function check_configuration(bnd_sys::BoundarySPHSystem, systems) + boundary_model = bnd_sys.boundary_model + for neighbor in systems + if neighbor isa WeaklyCompressibleSPHSystem && + boundary_model isa BoundaryModelDummyParticles && + isnothing(boundary_model.state_equation) + throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) + end + end +end From da5cee74df3b359187989efb979a55a130a05955 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 29 Sep 2023 18:09:17 +0200 Subject: [PATCH 47/49] format --- src/general/semidiscretization.jl | 6 ++---- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 1 - 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 5618586d2..7e0a3842e 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -486,18 +486,16 @@ function check_configuration(systems) for sys in systems check_configuration(sys, systems) end - end check_configuration(sys, systems) = systems - function check_configuration(bnd_sys::BoundarySPHSystem, systems) boundary_model = bnd_sys.boundary_model for neighbor in systems if neighbor isa WeaklyCompressibleSPHSystem && - boundary_model isa BoundaryModelDummyParticles && - isnothing(boundary_model.state_equation) + boundary_model isa BoundaryModelDummyParticles && + isnothing(boundary_model.state_equation) throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) end end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 1b5021610..413de32e1 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -44,7 +44,6 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, rho_a, rho_b, pos_diff, distance, grad_kernel, density_calculator) - dv_viscosity = viscosity_correction * viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, From dc1b4dbe14bc4161b74e63d85a84f6579289607e Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 2 Oct 2023 13:14:48 +0200 Subject: [PATCH 48/49] review comments --- src/callbacks/info.jl | 2 +- src/general/semidiscretization.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks/info.jl b/src/callbacks/info.jl index c197d7321..5e520db23 100644 --- a/src/callbacks/info.jl +++ b/src/callbacks/info.jl @@ -115,7 +115,7 @@ function initialize_info_callback(discrete_callback, u, t, integrator; show(io_context, MIME"text/plain"(), semi) println(io, "\n") systems = semi.systems - for system in systems + foreach_enumerate(systems) do (system_index, system) show(io_context, MIME"text/plain"(), system) println(io, "\n") end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 7e0a3842e..ac966f128 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -483,7 +483,7 @@ function nhs_coords(system::BoundarySPHSystem, end function check_configuration(systems) - for sys in systems + foreach_enumerate(systems) do (system_index, sys) check_configuration(sys, systems) end end @@ -492,7 +492,7 @@ check_configuration(sys, systems) = systems function check_configuration(bnd_sys::BoundarySPHSystem, systems) boundary_model = bnd_sys.boundary_model - for neighbor in systems + foreach_enumerate(systems) do (nghb_index, neighbor) if neighbor isa WeaklyCompressibleSPHSystem && boundary_model isa BoundaryModelDummyParticles && isnothing(boundary_model.state_equation) From 16584f5f1e8d4af52305f5860d74677e93dcb6e5 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 2 Oct 2023 15:45:20 +0200 Subject: [PATCH 49/49] fix review comments --- src/general/semidiscretization.jl | 17 ++++++++++------- .../fluid/weakly_compressible_sph/rhs.jl | 4 ++-- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index ac966f128..d91d0f8b2 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -31,6 +31,8 @@ struct Semidiscretization{S, RU, RV, NS, DC} ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) + # Check that the boundary systems are using a state equation if EDAC is not used. + # Other checks might be added here later. check_configuration(systems) # Create (and initialize) a tuple of n neighborhood searches for each of the n systems @@ -483,20 +485,21 @@ function nhs_coords(system::BoundarySPHSystem, end function check_configuration(systems) - foreach_enumerate(systems) do (system_index, sys) - check_configuration(sys, systems) + foreach_enumerate(systems) do (system_index, system) + check_configuration(system, systems) end end -check_configuration(sys, systems) = systems +check_configuration(system, systems) = nothing + +function check_configuration(boundary_system::BoundarySPHSystem, systems) + (; boundary_model) = boundary_system -function check_configuration(bnd_sys::BoundarySPHSystem, systems) - boundary_model = bnd_sys.boundary_model - foreach_enumerate(systems) do (nghb_index, neighbor) + foreach_enumerate(systems) do (neighbor_index, neighbor) if neighbor isa WeaklyCompressibleSPHSystem && boundary_model isa BoundaryModelDummyParticles && isnothing(boundary_model.state_equation) - throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a state_equation for boundary.")) + throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a `state_equation` for all boundary systems")) end end end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 413de32e1..ef5c2bbf3 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -73,7 +73,7 @@ end # As shown in "Variational and momentum preservation aspects of Smooth # Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation -# this form has to be used with ContinuityDensity +# this form has to be used with ContinuityDensity. @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem, @@ -87,7 +87,7 @@ end # As shown in "Variational and momentum preservation aspects of Smooth # Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation -# this form has to be used with SummationDensity +# this form has to be used with SummationDensity. @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system, v_particle_system, neighbor, neighbor_system::WeaklyCompressibleSPHSystem,