Validation Oscillating Beam

examples/fsi/falling_water_column_2d.jl

@@ -58,7 +58,8 @@ boundary_model = BoundaryModelMonaghanKajtar(k, spacing_ratio, particle_spacing,
 
 solid_system = TotalLagrangianSPHSystem(solid,
                                         smoothing_kernel, smoothing_length,
-                                        E, nu, boundary_model=boundary_model,
+                                        material.E, material.nu,
+                                        boundary_model=boundary_model,
                                         n_fixed_particles=nparticles(fixed_particles),
                                         acceleration=(0.0, -gravity))
 

examples/postprocessing/postprocessing.jl

@@ -13,15 +13,15 @@ function hello(v, u, t, system)
     # Value stored for output in the postprocessing output file
     return 2 * t
 end
-example_cb = PostprocessCallback(; interval=10, hello)
+example_cb = PostprocessCallback(; interval=10, hello, write_file_interval=0)
 
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"),
               extra_callback=example_cb, tspan=(0.0, 0.1));
 
 # Lets write the average pressure and kinetic energy every 0.01s
 pp = PostprocessCallback(; dt=0.005, filename="example_pressure_ekin", avg_pressure,
-                         kinetic_energy)
+                         kinetic_energy, write_file_interval=0)
 
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"),

examples/solid/oscillating_beam_2d.jl

@@ -10,36 +10,32 @@ n_particles_y = 5
 gravity = 2.0
 tspan = (0.0, 5.0)
 
-length_beam = 0.35
-thickness = 0.02
+elastic_beam = (length=0.35, thickness=0.02)
+material = (density=1000.0, E=1.4e6, nu=0.4)
 clamp_radius = 0.05
-material_density = 1000.0
-
-# Young's modulus and Poisson ratio
-E = 1.4e6
-nu = 0.4
 
 # The structure starts at the position of the first particle and ends
 # at the position of the last particle.
-particle_spacing = thickness / (n_particles_y - 1)
+particle_spacing = elastic_beam.thickness / (n_particles_y - 1)
 
 # Add particle_spacing/2 to the clamp_radius to ensure that particles are also placed on the radius
 fixed_particles = SphereShape(particle_spacing, clamp_radius + particle_spacing / 2,
-                              (0.0, thickness / 2), material_density,
-                              cutout_min=(0.0, 0.0), cutout_max=(clamp_radius, thickness),
+                              (0.0, elastic_beam.thickness / 2), material.density,
+                              cutout_min=(0.0, 0.0),
+                              cutout_max=(clamp_radius, elastic_beam.thickness),
                               tlsph=true)
 
 n_particles_clamp_x = round(Int, clamp_radius / particle_spacing)
 
 # Beam and clamped particles
-n_particles_per_dimension = (round(Int, length_beam / particle_spacing) +
+n_particles_per_dimension = (round(Int, elastic_beam.length / particle_spacing) +
                              n_particles_clamp_x + 1, n_particles_y)
 
 # Note that the `RectangularShape` puts the first particle half a particle spacing away
 # from the boundary, which is correct for fluids, but not for solids.
 # We therefore need to pass `tlsph=true`.
 beam = RectangularShape(particle_spacing, n_particles_per_dimension,
-                        (0.0, 0.0), density=material_density, tlsph=true)
+                        (0.0, 0.0), density=material.density, tlsph=true)
 
 solid = union(beam, fixed_particles)
 
@@ -51,9 +47,10 @@ smoothing_length = 2 * sqrt(2) * particle_spacing
 smoothing_kernel = WendlandC2Kernel{2}()
 
 solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length,
-                                        E, nu,
+                                        material.E, material.nu,
                                         n_fixed_particles=nparticles(fixed_particles),
-                                        acceleration=(0.0, -gravity))
+                                        acceleration=(0.0, -gravity),
+                                        penalty_force=nothing)
 
 # ==========================================================================================
 # ==== Simulation
@@ -63,26 +60,22 @@ ode = semidiscretize(semi, tspan)
 info_callback = InfoCallback(interval=100)
 
 # Track the position of the particle in the middle of the tip of the beam.
-particle_id = Int(n_particles_per_dimension[1] * (n_particles_per_dimension[2] + 1) / 2)
-
-shift_x = beam.coordinates[1, particle_id]
-shift_y = beam.coordinates[2, particle_id]
+middle_particle_id = Int(n_particles_per_dimension[1] * (n_particles_per_dimension[2] + 1) /
+                         2)
+startposition_x = beam.coordinates[1, middle_particle_id]
+startposition_y = beam.coordinates[2, middle_particle_id]
 
-function x_deflection(v, u, t, system)
-    particle_position = TrixiParticles.current_coords(u, system, particle_id)
-
-    return particle_position[1] - shift_x
+function deflection_x(v, u, t, system)
+    return system.current_coordinates[1, middle_particle_id] - startposition_x
 end
 
-function y_deflection(v, u, t, system)
-    particle_position = TrixiParticles.current_coords(u, system, particle_id)
-
-    return particle_position[2] - shift_y
+function deflection_y(v, u, t, system)
+    return system.current_coordinates[2, middle_particle_id] - startposition_y
 end
 
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="",
-                                         x_deflection=x_deflection,
-                                         y_deflection=y_deflection)
+                                         deflection_x=deflection_x,
+                                         deflection_y=deflection_y)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

src/TrixiParticles.jl

@@ -58,7 +58,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx
        PressureMirroring, PressureZeroing
 export BoundaryMovement
 export GridNeighborhoodSearch, TrivialNeighborhoodSearch
-export examples_dir, trixi_include
+export examples_dir, validation_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape
 export VoxelSphere, RoundSphere, reset_wall!

src/schemes/solid/total_lagrangian_sph/penalty_force.jl

@@ -52,16 +52,16 @@ struct PenaltyForceGanzenmueller{ELTYPE}
 end
 
 @inline function calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                                     initial_distance, system,
+                                     initial_distance, system, m_a, m_b, rho_a, rho_b,
                                      penalty_force::PenaltyForceGanzenmueller)
-    (; mass, material_density, young_modulus) = system
+    (; young_modulus) = system
 
     current_pos_diff = current_coords(system, particle) -
                        current_coords(system, neighbor)
     current_distance = norm(current_pos_diff)
 
-    volume_particle = mass[particle] / material_density[particle]
-    volume_neighbor = mass[neighbor] / material_density[neighbor]
+    volume_particle = m_a / rho_a
+    volume_neighbor = m_b / rho_b
 
     kernel_weight = smoothing_kernel(system, initial_distance)
 
@@ -76,9 +76,9 @@ end
         kernel_weight / initial_distance^2 * young_modulus * delta_sum *
         current_pos_diff / current_distance
 
-    for i in 1:ndims(system)
+    @inbounds for i in 1:ndims(system)
         # Divide force by mass to obtain acceleration
-        dv[i, particle] += f[i] / mass[particle]
+        dv[i, particle] += f[i] / m_a
     end
 
     return dv

src/schemes/solid/total_lagrangian_sph/rhs.jl

@@ -35,19 +35,21 @@ end
         grad_kernel = smoothing_kernel_grad(particle_system, initial_pos_diff,
                                             initial_distance)
 
+        m_a = particle_system.mass[particle]
         m_b = neighbor_system.mass[neighbor]
 
         dv_particle = m_b *
                       (pk1_corrected(particle_system, particle) / rho_a^2 +
                        pk1_corrected(neighbor_system, neighbor) / rho_b^2) *
                       grad_kernel
 
-        for i in 1:ndims(particle_system)
+        @inbounds for i in 1:ndims(particle_system)
             dv[i, particle] += dv_particle[i]
         end
 
         calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                            initial_distance, particle_system, penalty_force)
+                            initial_distance, particle_system, m_a, m_b, rho_a, rho_b,
+                            penalty_force)
 
         # TODO continuity equation?
     end

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -383,7 +383,8 @@ end
 end
 
 @inline function calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                                     initial_distance, system, ::Nothing)
+                                     initial_distance, system, m_a, m_b, rho_a, rho_b,
+                                     ::Nothing)
     return dv
 end
 

src/util.jl

@@ -140,6 +140,22 @@ readdir(examples_dir())
 """
 examples_dir() = pkgdir(TrixiParticles, "examples")
 
+"""
+    validation_dir()
+
+Return the directory where the validation files provided with TrixiParticles.jl are located. If TrixiParticles is
+installed as a regular package (with `]add TrixiParticles`), these files are read-only and should *not* be
+modified. To find out which files are available, use, e.g., `readdir`.
+
+Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).
+
+# Examples
+```@example
+readdir(validation_dir())
+```
+"""
+validation_dir() = pkgdir(TrixiParticles, "validation")
+
 """
     @autoinfiltrate
     @autoinfiltrate condition::Bool

validation/investigate_relaxation.jl → validation/general/investigate_relaxation.jl

@@ -6,11 +6,12 @@ using DataFrames
 using Printf
 
 pp_cb = PostprocessCallback(; ekin=kinetic_energy, max_pressure, avg_density, dt=0.025,
-                            filename="relaxation", write_csv=false)
+                            filename="relaxation", write_csv=false, write_file_interval=0)
 
 pp_damped_cb = PostprocessCallback(; ekin=kinetic_energy, max_pressure, avg_density,
                                    dt=0.025,
-                                   filename="relaxation_damped", write_csv=false)
+                                   filename="relaxation_damped", write_csv=false,
+                                   write_file_interval=0)
 
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"),

validation/oscillating_beam_2d/README.md

@@ -0,0 +1,9 @@
+The following files are provided here:
+
+1) Validation simulation: validation_oscillating_beam_2d.jl
+2) Comparison with reference: compare_reference_oscillating_beam_results.jl
+3) Comparison with TrixiParticles.jl and literature reference: plot_oscillating_beam_results.jl
+4) TrixiParticles.jl reference files: validation_reference_[5, 9, 21, 35].j (5 is the default CI resolution, note that resolution 21 takes about 30-50 minutes while 35 takes about 4-6 hours)
+5) Reference file reference_turek.csv extracted from the file csm3_l4_t0p005.point here:
+   https://wwwold.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/fsi_benchmark/fsi_tests/fsi_csm_tests.html
+

validation/oscillating_beam_2d/plot_oscillating_beam_results.jl

@@ -0,0 +1,80 @@
+include("../validation_util.jl")
+
+# Activate for interactive plot
+# using GLMakie
+using CairoMakie
+using CSV
+using DataFrames
+using JSON
+using Glob
+using Printf
+using TrixiParticles
+
+elastic_plate = (length=0.35, thickness=0.02)
+
+# Load the reference simulation data
+ref = CSV.read(joinpath(validation_dir(), "oscillating_beam_2d/reference_turek.csv"),
+               DataFrame)
+
+# Get the list of JSON files
+reference_files = glob("validation_reference_*.json",
+                       joinpath(validation_dir(), "oscillating_beam_2d"))
+simulation_files = glob("validation_run_oscillating_beam_2d_*.json",
+                        joinpath(pkgdir(TrixiParticles), "out"))
+merged_files = vcat(reference_files, simulation_files)
+input_files = sort(merged_files, by=extract_number)
+
+# Regular expressions for matching keys
+key_pattern_x = r"deflection_x_solid_\d+"
+key_pattern_y = r"deflection_y_solid_\d+"
+
+# Setup for Makie plotting
+fig = Figure(size=(1200, 800))
+ax1 = Axis(fig, title="X-Axis Displacement", xlabel="Time [s]", ylabel="X Displacement")
+ax2 = Axis(fig, title="Y-Axis Displacement", xlabel="Time [s]", ylabel="Y Displacement")
+fig[1, 1] = ax1
+fig[2, 1] = ax2
+
+for file_name in input_files
+    println("Loading the input file: $file_name")
+    json_data = JSON.parsefile(file_name)
+
+    resolution = parse(Int, split(split(file_name, "_")[end], ".")[1])
+    particle_spacing = elastic_plate.thickness / (resolution - 1)
+
+    matching_keys_x = sort(collect(filter(key -> occursin(key_pattern_x, key),
+                                          keys(json_data))))
+    matching_keys_y = sort(collect(filter(key -> occursin(key_pattern_y, key),
+                                          keys(json_data))))
+
+    if isempty(matching_keys_x)
+        error("No matching keys found in: $file_name")
+    end
+
+    label_prefix = occursin("reference", file_name) ? "Reference" : ""
+
+    for (matching_keys, ax) in ((matching_keys_x, ax1), (matching_keys_y, ax2))
+        for key in matching_keys
+            data = json_data[key]
+            times = Float64.(data["time"])
+            displacements = Float64.(data["values"])
+
+            mse_results = occursin(key_pattern_x, key) ?
+                          interpolated_mse(ref.time, ref.Ux, data["time"], displacements) :
+                          interpolated_mse(ref.time, ref.Uy, data["time"], displacements)
+
+            label = "$label_prefix dp = $(@sprintf("%.8f", particle_spacing)) mse=$(@sprintf("%.8f", mse_results))"
+            lines!(ax, times, displacements, label=label)
+        end
+    end
+end
+
+# Plot reference data
+lines!(ax1, ref.time, ref.Ux, color=:black, linestyle=:dash,
+       label="Turek and Hron 2006")
+lines!(ax2, ref.time, ref.Uy, color=:black, linestyle=:dash,
+       label="Turek and Hron 2006")
+
+legend_ax1 = Legend(fig[1, 2], ax1)
+legend_ax2 = Legend(fig[2, 2], ax2)
+fig