From 4aad0732f9d7c2f9bf435691f523356f81315a17 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 31 Jan 2024 15:53:50 +0100 Subject: [PATCH 01/12] add example --- examples/fsi/hydrostatic_water_column_2d.jl | 127 ++++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 examples/fsi/hydrostatic_water_column_2d.jl diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl new file mode 100644 index 0000000000..bdb7494c0e --- /dev/null +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -0,0 +1,127 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +solid_particle_spacing = 0.01 +fluid_particle_spacing = solid_particle_spacing + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 3 +spacing_ratio = 1 + +boundary_particle_spacing = fluid_particle_spacing / spacing_ratio + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (1.0, 2.0) +plate_size = (1.0, 0.05) + +fluid_density = 1000.0 +solid_density = 2700.0 + +# Young's modulus and Poisson ratio +E = 67.5e9 +nu = 0.34 + +sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +boundary_density_calculator = AdamiPressureExtrapolation() + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, initial_fluid_size, + min_coordinates=(0.0, fluid_particle_spacing / 2), + fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, false, false), acceleration=(0.0, -gravity), + state_equation=state_equation) + +# Beam and clamped particles +n_particles_plate_x = round(Int, plate_size[1] / solid_particle_spacing + 1) +n_particles_plate_y = round(Int, plate_size[2] / solid_particle_spacing + 1) +n_particles_per_dimension = (n_particles_plate_x, n_particles_plate_y) + +plate = RectangularShape(solid_particle_spacing, n_particles_per_dimension, + (0.0, -plate_size[2]), density=solid_density, tlsph=true) + +fixed_particles_1 = RectangularShape(solid_particle_spacing, (3, n_particles_plate_y), + (-3solid_particle_spacing, -plate_size[2]), + density=solid_density, tlsph=true) +fixed_particles_2 = RectangularShape(solid_particle_spacing, (3, n_particles_plate_y), + (plate_size[1] + solid_particle_spacing, + -plate_size[2]), + density=solid_density, tlsph=true) +fixed_particles = union(fixed_particles_1, fixed_particles_2) + +solid = union(plate, fixed_particles) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 1.2 * fluid_particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) + +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Solid +smoothing_length_solid = sqrt(2) * solid_particle_spacing + +# For the FSI we need the hydrodynamic masses and densities in the solid boundary model +hydrodynamic_densites = fluid_density * ones(size(solid.density)) +hydrodynamic_masses = hydrodynamic_densites * solid_particle_spacing^ndims(solid) + +boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites, + hydrodynamic_masses, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_solid) + +solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length_solid, + E, nu, boundary_model_solid, + n_fixed_particles=nparticles(fixed_particles), + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(solid_system, fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +const PARTICLE_ID = Int(prod(n_particles_per_dimension) - + (n_particles_per_dimension[1] - 1) / 2) + +y_deflection(v, u, t, system) = nothing + +function y_deflection(v, u, t, system::TotalLagrangianSPHSystem) + return TrixiParticles.current_coords(u, system, PARTICLE_ID)[2] +end + +saving_callback = SolutionSavingCallback(dt=0.0005, prefix="", y_deflection=y_deflection) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(), save_everystep=false, callback=callbacks); From 3224b304ac55969c3760836700b7a5a1366c7f74 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 31 Jan 2024 15:54:00 +0100 Subject: [PATCH 02/12] fix viscosity --- src/schemes/solid/total_lagrangian_sph/rhs.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index e3994db6b5..f14229ab6c 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -60,9 +60,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::TotalLagrangianSPHSystem, neighbor_system::WeaklyCompressibleSPHSystem) - (; boundary_model) = particle_system - (; density_calculator, state_equation, viscosity, correction) = neighbor_system - (; sound_speed) = state_equation + (; sound_speed) = neighbor_system.state_equation + viscosity = viscosity_model(particle_system) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -109,7 +108,8 @@ function interact!(dv, v_particle_system, u_particle_system, # `pressure_correction` is set to `1.0` (no correction). dv_boundary = pressure_acceleration(neighbor_system, particle_system, particle, m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, - distance, grad_kernel, 1.0, correction) + distance, grad_kernel, 1.0, + neighbor_system.correction) dv_particle = dv_boundary + dv_viscosity for i in 1:ndims(particle_system) From 1288c183114f504d594c33febb6e60a6be5610ef Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 1 Feb 2024 17:07:00 +0100 Subject: [PATCH 03/12] add analytical solution --- examples/fsi/hydrostatic_water_column_2d.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index bdb7494c0e..3cd3279260 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -124,4 +124,9 @@ saving_callback = SolutionSavingCallback(dt=0.0005, prefix="", y_deflection=y_de callbacks = CallbackSet(info_callback, saving_callback) # Use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL49(), save_everystep=false, callback=callbacks); +sol = nothing#solve(ode, RDPK3SpFSAL49(), save_everystep=false, callback=callbacks); + +# The pseudostatic midpoint deflection of a 2-D plate: +D = E*plate_size[2]^3/(12*(1-nu^2)) +analytical_sol = 0.0026 * gravity * + (fluid_density * initial_fluid_size[2] + solid_density * plate_size[2]) / D From 3b9b8a9599d10b34994159c637b33371851ee644 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 12:08:45 +0100 Subject: [PATCH 04/12] add backup for pp callback --- src/callbacks/post_process.jl | 36 +++++++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 33d544dc29..725e9f24a2 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -28,6 +28,7 @@ a fixed interval of simulation time (`dt`). - `write_csv=true`: If set to `true`, write a csv file. - `write_json=true`: If set to `true`, write a json file. - `append_timestep=false`: If set to `true`, the current timestamp will be added to the filename. +- `backup_period=0`: Write a backup for each multiple of `interval` or `dt`. # Examples ```julia @@ -44,6 +45,7 @@ postprocess_callback = PostprocessCallback(example_function, dt=0.1) """ struct PostprocessCallback{I, F} interval :: I + backup_period :: Int data :: Dict{String, Vector{Any}} times :: Array{Float64, 1} exclude_boundary :: Bool @@ -53,12 +55,13 @@ struct PostprocessCallback{I, F} append_timestamp :: Bool write_csv :: Bool write_json :: Bool + write_backup :: Bool end function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=true, output_directory="out", filename="values", append_timestamp=false, write_csv=true, write_json=true, - funcs...) + backup_period::Integer=0, funcs...) if isempty(funcs) throw(ArgumentError("`funcs` cannot be empty")) end @@ -71,9 +74,13 @@ function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=tru interval = Float64(dt) end - post_callback = PostprocessCallback(interval, Dict{String, Vector{Any}}(), Float64[], + write_backup = backup_period > 0 ? true : false + + post_callback = PostprocessCallback(interval, backup_period, + Dict{String, Vector{Any}}(), Float64[], exclude_boundary, funcs, filename, output_directory, - append_timestamp, write_csv, write_json) + append_timestamp, write_csv, write_json, + write_backup) if dt > 0 # Add a `tstop` every `dt`, and save the final solution. return PeriodicCallback(post_callback, dt, @@ -208,16 +215,26 @@ function (pp::PostprocessCallback)(integrator) push!(pp.times, t) end - if isfinished(integrator) - write_postprocess_callback(pp) + if isfinished(integrator) || (pp.write_backup && backup_condition(pp, integrator)) + write_postprocess_callback(pp, integrator) end # Tell OrdinaryDiffEq that u has not been modified u_modified!(integrator, false) end +@inline function backup_condition(cb::PostprocessCallback{Int}, integrator) + @autoinfiltrate + + return integrator.stats.naccept % cb.backup_period == 0 +end + +@inline function backup_condition(cb::PostprocessCallback, integrator) + return round(Int, integrator.t / cb.interval) % cb.backup_period == 0 +end + # After the simulation has finished, this function is called to write the data to a JSON file -function write_postprocess_callback(pp::PostprocessCallback) +function write_postprocess_callback(pp::PostprocessCallback, integrator) if isempty(pp.data) return end @@ -230,8 +247,11 @@ function write_postprocess_callback(pp::PostprocessCallback) if pp.append_timestamp time_stamp = string("_", Dates.format(now(), "YY-mm-ddTHHMMSS")) end - filename_json = pp.filename * time_stamp * ".json" - filename_csv = pp.filename * time_stamp * ".csv" + + backup_suffix = isfinished(integrator) ? "" : "_backup" + + filename_json = pp.filename * time_stamp * backup_suffix * ".json" + filename_csv = pp.filename * time_stamp * backup_suffix * ".csv" if pp.write_json abs_file_path = joinpath(abspath(pp.output_directory), filename_json) From 09560d64db02806cab130ccf37cc7a1407bbd3da Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 12:11:41 +0100 Subject: [PATCH 05/12] remove autoinfiltrator --- src/callbacks/post_process.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 725e9f24a2..f3bdd8a1c4 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -74,7 +74,7 @@ function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=tru interval = Float64(dt) end - write_backup = backup_period > 0 ? true : false + write_backup = backup_period > 0 ? true : false post_callback = PostprocessCallback(interval, backup_period, Dict{String, Vector{Any}}(), Float64[], @@ -224,8 +224,6 @@ function (pp::PostprocessCallback)(integrator) end @inline function backup_condition(cb::PostprocessCallback{Int}, integrator) - @autoinfiltrate - return integrator.stats.naccept % cb.backup_period == 0 end From db4dcae9edd4ad84f83116e368123a73e14fca8a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 12:37:42 +0100 Subject: [PATCH 06/12] fix bug --- src/callbacks/post_process.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index f3bdd8a1c4..ff8a1dabd2 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -224,7 +224,7 @@ function (pp::PostprocessCallback)(integrator) end @inline function backup_condition(cb::PostprocessCallback{Int}, integrator) - return integrator.stats.naccept % cb.backup_period == 0 + return round(integrator.stats.naccept / cb.interval) % cb.backup_period == 0 end @inline function backup_condition(cb::PostprocessCallback, integrator) From 558484bcbd8911e750eb8a8ab629f2c38aebbf97 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 13:00:11 +0100 Subject: [PATCH 07/12] allow `nothing` --- src/callbacks/post_process.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index ff8a1dabd2..96ba9c443e 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -206,8 +206,10 @@ function (pp::PostprocessCallback)(integrator) u = wrap_u(u_ode, system, semi) for (key, f) in pp.func result = f(v, u, t, system) - add_entry!(pp, string(key), t, result, filenames[system_index]) - new_data = true + if result !== nothing + add_entry!(pp, string(key), t, result, filenames[system_index]) + new_data = true + end end end @@ -224,11 +226,13 @@ function (pp::PostprocessCallback)(integrator) end @inline function backup_condition(cb::PostprocessCallback{Int}, integrator) - return round(integrator.stats.naccept / cb.interval) % cb.backup_period == 0 + return integrator.stats.naccept > 0 && + round(integrator.stats.naccept / cb.interval) % cb.backup_period == 0 end @inline function backup_condition(cb::PostprocessCallback, integrator) - return round(Int, integrator.t / cb.interval) % cb.backup_period == 0 + return integrator.stats.naccept > 0 && + round(Int, integrator.t / cb.interval) % cb.backup_period == 0 end # After the simulation has finished, this function is called to write the data to a JSON file From 144b34da442a45176d6073d8b07c6b96f1c9da74 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 14:59:53 +0100 Subject: [PATCH 08/12] implement suggestions --- src/callbacks/post_process.jl | 18 +++++++++-------- test/callbacks/postprocess.jl | 37 +++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 8 deletions(-) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 96ba9c443e..6b9d42ec03 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -1,7 +1,7 @@ """ PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=true, filename="values", output_directory="out", append_timestamp=false, write_csv=true, - write_json=true, funcs...) + write_json=true, backup_period=0, funcs...) Create a callback to post-process simulation data at regular intervals. This callback allows for the execution of a user-defined function `func` at specified @@ -28,7 +28,9 @@ a fixed interval of simulation time (`dt`). - `write_csv=true`: If set to `true`, write a csv file. - `write_json=true`: If set to `true`, write a json file. - `append_timestep=false`: If set to `true`, the current timestamp will be added to the filename. -- `backup_period=0`: Write a backup for each multiple of `interval` or `dt`. +- `backup_period=0`: Specifies that a backup should be created after every `backup_period` + number of postprocessing execution steps. A value of 0 indicates that + backups should not be automatically generated during postprocessing. # Examples ```julia @@ -55,7 +57,6 @@ struct PostprocessCallback{I, F} append_timestamp :: Bool write_csv :: Bool write_json :: Bool - write_backup :: Bool end function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=true, @@ -74,13 +75,10 @@ function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=tru interval = Float64(dt) end - write_backup = backup_period > 0 ? true : false - post_callback = PostprocessCallback(interval, backup_period, Dict{String, Vector{Any}}(), Float64[], exclude_boundary, funcs, filename, output_directory, - append_timestamp, write_csv, write_json, - write_backup) + append_timestamp, write_csv, write_json) if dt > 0 # Add a `tstop` every `dt`, and save the final solution. return PeriodicCallback(post_callback, dt, @@ -123,6 +121,8 @@ function Base.show(io::IO, ::MIME"text/plain", callback = cb.affect! setup = [ "interval" => string(callback.interval), + "write backup" => callback.backup_period > 0 ? + "every $(callback.backup_period) * interval" : "no", "exclude boundary" => callback.exclude_boundary ? "yes" : "no", "filename" => callback.filename, "output directory" => callback.output_directory, @@ -148,6 +148,8 @@ function Base.show(io::IO, ::MIME"text/plain", callback = cb.affect!.affect! setup = [ "dt" => string(callback.interval), + "write backup" => callback.backup_period > 0 ? + "every $(callback.backup_period) * dt" : "no", "exclude boundary" => callback.exclude_boundary ? "yes" : "no", "filename" => callback.filename, "output directory" => callback.output_directory, @@ -217,7 +219,7 @@ function (pp::PostprocessCallback)(integrator) push!(pp.times, t) end - if isfinished(integrator) || (pp.write_backup && backup_condition(pp, integrator)) + if isfinished(integrator) || (pp.backup_period > 0 && backup_condition(pp, integrator)) write_postprocess_callback(pp, integrator) end diff --git a/test/callbacks/postprocess.jl b/test/callbacks/postprocess.jl index 1d0dd95153..a3ffca9b93 100644 --- a/test/callbacks/postprocess.jl +++ b/test/callbacks/postprocess.jl @@ -15,6 +15,7 @@ │ PostprocessCallback │ │ ═══════════════════ │ │ interval: ……………………………………………………… 10 │ + │ write backup: …………………………………………… no │ │ exclude boundary: ………………………………… yes │ │ filename: ……………………………………………………… values │ │ output directory: ………………………………… out │ @@ -36,6 +37,42 @@ │ PostprocessCallback │ │ ═══════════════════ │ │ dt: ……………………………………………………………………… 0.1 │ + │ write backup: …………………………………………… no │ + │ exclude boundary: ………………………………… yes │ + │ filename: ……………………………………………………… values │ + │ output directory: ………………………………… out │ + │ append timestamp: ………………………………… no │ + │ write json file: …………………………………… yes │ + │ write csv file: ……………………………………… yes │ + │ function1: …………………………………………………… example_function │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + + callback = PostprocessCallback(; dt=0.1, example_function, backup_period=3) + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ PostprocessCallback │ + │ ═══════════════════ │ + │ dt: ……………………………………………………………………… 0.1 │ + │ write backup: …………………………………………… every 3 * dt │ + │ exclude boundary: ………………………………… yes │ + │ filename: ……………………………………………………… values │ + │ output directory: ………………………………… out │ + │ append timestamp: ………………………………… no │ + │ write json file: …………………………………… yes │ + │ write csv file: ……………………………………… yes │ + │ function1: …………………………………………………… example_function │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + + callback = PostprocessCallback(; interval=23, example_function, backup_period=4) + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ PostprocessCallback │ + │ ═══════════════════ │ + │ interval: ……………………………………………………… 23 │ + │ write backup: …………………………………………… every 4 * interval │ │ exclude boundary: ………………………………… yes │ │ filename: ……………………………………………………… values │ │ output directory: ………………………………… out │ From 0ab92a4257ed4dd5c5834fe28b796c4fe903f092 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 16:22:24 +0100 Subject: [PATCH 09/12] add `PostProcessCallback` in example --- examples/fsi/hydrostatic_water_column_2d.jl | 79 ++++++++++++--------- 1 file changed, 44 insertions(+), 35 deletions(-) diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index 3cd3279260..c2324313a0 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -1,17 +1,16 @@ using TrixiParticles using OrdinaryDiffEq +ddt = false + # ========================================================================================== # ==== Resolution -solid_particle_spacing = 0.01 -fluid_particle_spacing = solid_particle_spacing +n_particles_plate_y = 5 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model boundary_layers = 3 spacing_ratio = 1 -boundary_particle_spacing = fluid_particle_spacing / spacing_ratio - # ========================================================================================== # ==== Experiment Setup gravity = 9.81 @@ -32,7 +31,9 @@ sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7, clip_negative_pressure=false) -boundary_density_calculator = AdamiPressureExtrapolation() +solid_particle_spacing = plate_size[2] / (n_particles_plate_y - 1) +fluid_particle_spacing = solid_particle_spacing +boundary_particle_spacing = fluid_particle_spacing / spacing_ratio tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, initial_fluid_size, min_coordinates=(0.0, fluid_particle_spacing / 2), @@ -42,7 +43,6 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, initial_fluid # Beam and clamped particles n_particles_plate_x = round(Int, plate_size[1] / solid_particle_spacing + 1) -n_particles_plate_y = round(Int, plate_size[2] / solid_particle_spacing + 1) n_particles_per_dimension = (n_particles_plate_x, n_particles_plate_y) plate = RectangularShape(solid_particle_spacing, n_particles_per_dimension, @@ -59,29 +59,16 @@ fixed_particles = union(fixed_particles_1, fixed_particles_2) solid = union(plate, fixed_particles) -# ========================================================================================== -# ==== Fluid -smoothing_length = 1.2 * fluid_particle_spacing -smoothing_kernel = SchoenbergQuinticSplineKernel{2}() - -fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) - -density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) -fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, - state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, - density_diffusion=density_diffusion, - acceleration=(0.0, -gravity)) - # ========================================================================================== # ==== Solid -smoothing_length_solid = sqrt(2) * solid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() +smoothing_length_solid = 2 * sqrt(2) * solid_particle_spacing # For the FSI we need the hydrodynamic masses and densities in the solid boundary model hydrodynamic_densites = fluid_density * ones(size(solid.density)) hydrodynamic_masses = hydrodynamic_densites * solid_particle_spacing^ndims(solid) +boundary_density_calculator = AdamiPressureExtrapolation() boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites, hydrodynamic_masses, state_equation=state_equation, @@ -89,10 +76,24 @@ boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites, smoothing_kernel, smoothing_length_solid) solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length_solid, - E, nu, boundary_model_solid, + E, nu, boundary_model=boundary_model_solid, n_fixed_particles=nparticles(fixed_particles), acceleration=(0.0, -gravity)) +# ========================================================================================== +# ==== Fluid +smoothing_length_fluid = 2 * sqrt(2) * fluid_particle_spacing + +fluid_density_calculator = ContinuityDensity() +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) + +density_diffusion = ddt ? DensityDiffusionMolteniColagrossi(delta=0.1) : nothing +fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length_fluid, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + # ========================================================================================== # ==== Boundary @@ -110,23 +111,31 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) -const PARTICLE_ID = Int(prod(n_particles_per_dimension) - - (n_particles_per_dimension[1] - 1) / 2) +# Track the position of the particle in the center of the of the beam. +particle_id = Int(n_particles_per_dimension[1] * (n_particles_plate_y + 1) / 2 - + (n_particles_per_dimension[1] + 1) / 2 + 1) +function y_deflection(v, u, t, system::TotalLagrangianSPHSystem) + return TrixiParticles.current_coords(u, system, particle_id)[2] + plate_size[2] / 2 +end y_deflection(v, u, t, system) = nothing -function y_deflection(v, u, t, system::TotalLagrangianSPHSystem) - return TrixiParticles.current_coords(u, system, PARTICLE_ID)[2] +# The pseudostatic midpoint deflection of a 2-D plate +function analytical_sol(v, u, t, system::TotalLagrangianSPHSystem) + # Flexural rigidity of the plate + D = E * plate_size[2]^3 / (12 * (1 - nu^2)) + + return -0.0026 * gravity * + (fluid_density * initial_fluid_size[2] + solid_density * plate_size[2]) / D end +analytical_sol(v, u, t, system) = nothing -saving_callback = SolutionSavingCallback(dt=0.0005, prefix="", y_deflection=y_deflection) +saving_callback = SolutionSavingCallback(dt=0.005, prefix="") -callbacks = CallbackSet(info_callback, saving_callback) +pp = PostprocessCallback(; interval=100, filename="hydrostatic_water_column_2d", + y_deflection, analytical_sol, kinetic_energy, backup_period=10) -# Use a Runge-Kutta method with automatic (error based) time step size control -sol = nothing#solve(ode, RDPK3SpFSAL49(), save_everystep=false, callback=callbacks); +callbacks = CallbackSet(info_callback, saving_callback, pp) -# The pseudostatic midpoint deflection of a 2-D plate: -D = E*plate_size[2]^3/(12*(1-nu^2)) -analytical_sol = 0.0026 * gravity * - (fluid_density * initial_fluid_size[2] + solid_density * plate_size[2]) / D +# Use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(), save_everystep=false, callback=callbacks); From 6f7402cb9d933e50ced6e51c31b87fff1025ae92 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 16:48:16 +0100 Subject: [PATCH 10/12] fix edac fsi --- src/schemes/fluid/entropically_damped_sph/rhs.jl | 14 ++++++++++++++ .../fluid/entropically_damped_sph/system.jl | 3 +++ .../fluid/weakly_compressible_sph/system.jl | 2 ++ src/schemes/solid/total_lagrangian_sph/rhs.jl | 6 +++--- 4 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index aac412dc19..76cb224033 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -80,3 +80,17 @@ end return dv end + +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + neighbor_system, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, pressure_correction, correction) + volume_a = m_a / rho_a + volume_b = m_b / rho_b + volume_term = (volume_a^2 + volume_b^2) / m_a + + # Inter-particle averaged pressure + pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b) + + return -volume_term * pressure_tilde * W_a +end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 9a9aea3719..2f81403595 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -68,6 +68,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, ST} <: viscosity :: V nu_edac :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} + correction :: Nothing source_terms :: ST function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, @@ -139,6 +140,8 @@ end @inline v_nvariables(system::EntropicallyDampedSPHSystem) = ndims(system) + 1 +@inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed + function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) summation_density!(system, semi, u, u_ode, system.density) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 1d23b0db8a..21aed1aaaa 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -201,6 +201,8 @@ end return system.pressure[particle] end +@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed + function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, density_diffusion, correction) = system diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index f14229ab6c..35ec50cb19 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -59,8 +59,8 @@ end function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::TotalLagrangianSPHSystem, - neighbor_system::WeaklyCompressibleSPHSystem) - (; sound_speed) = neighbor_system.state_equation + neighbor_system::FluidSystem) + sound_speed = system_sound_speed(neighbor_system) viscosity = viscosity_model(particle_system) system_coords = current_coordinates(u_particle_system, particle_system) @@ -133,7 +133,7 @@ end particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, particle_system::TotalLagrangianSPHSystem, - neighbor_system::WeaklyCompressibleSPHSystem, + neighbor_system::FluidSystem, grad_kernel) return dv end From fcfa4b358f932a0b0d318478242d54dc1a4ddd75 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 19 Feb 2024 17:22:56 +0100 Subject: [PATCH 11/12] add edac system --- examples/fsi/hydrostatic_water_column_2d.jl | 30 +++++++++++++-------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index c2324313a0..b18bdb9d6b 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -2,6 +2,7 @@ using TrixiParticles using OrdinaryDiffEq ddt = false +edac = true # ========================================================================================== # ==== Resolution @@ -28,7 +29,8 @@ E = 67.5e9 nu = 0.34 sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, +state_equation = edac ? nothing : + StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7, clip_negative_pressure=false) solid_particle_spacing = plate_size[2] / (n_particles_plate_y - 1) @@ -84,15 +86,21 @@ solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_lengt # ==== Fluid smoothing_length_fluid = 2 * sqrt(2) * fluid_particle_spacing -fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) - -density_diffusion = ddt ? DensityDiffusionMolteniColagrossi(delta=0.1) : nothing -fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, - state_equation, smoothing_kernel, - smoothing_length_fluid, viscosity=viscosity, - density_diffusion=density_diffusion, - acceleration=(0.0, -gravity)) +if edac + fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, + smoothing_length_fluid, + sound_speed, + acceleration=(0.0, -gravity)) +else + fluid_density_calculator = ContinuityDensity() + + density_diffusion = ddt ? DensityDiffusionMolteniColagrossi(delta=0.1) : nothing + fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length_fluid, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) +end # ========================================================================================== # ==== Boundary @@ -100,7 +108,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, - smoothing_kernel, smoothing_length) + smoothing_kernel, smoothing_length_fluid) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) From 37b13e67498bdc3d33fc90489a221388a37a336c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 14 Nov 2024 15:51:56 +0100 Subject: [PATCH 12/12] fix example --- examples/fsi/hydrostatic_water_column_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index b18bdb9d6b..9aba8a8eb0 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -6,7 +6,7 @@ edac = true # ========================================================================================== # ==== Resolution -n_particles_plate_y = 5 +n_particles_plate_y = 3 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model boundary_layers = 3 @@ -18,7 +18,7 @@ gravity = 9.81 tspan = (0.0, 1.0) # Boundary geometry and initial fluid particle positions -initial_fluid_size = (1.0, 2.0) +initial_fluid_size = (1.0, 1.0) plate_size = (1.0, 0.05) fluid_density = 1000.0 @@ -141,7 +141,7 @@ analytical_sol(v, u, t, system) = nothing saving_callback = SolutionSavingCallback(dt=0.005, prefix="") pp = PostprocessCallback(; interval=100, filename="hydrostatic_water_column_2d", - y_deflection, analytical_sol, kinetic_energy, backup_period=10) + y_deflection, analytical_sol, kinetic_energy, write_file_interval=10) callbacks = CallbackSet(info_callback, saving_callback, pp)