diff --git a/docs/src/general/smoothing_kernels.md b/docs/src/general/smoothing_kernels.md index ea36cbf8aa..bdd368edae 100644 --- a/docs/src/general/smoothing_kernels.md +++ b/docs/src/general/smoothing_kernels.md @@ -13,6 +13,7 @@ The following smoothing kernels are currently available: | [`Poly6Kernel`](@ref) | $[0, 1h]$ | $1.5$ to $2.5$ | Academic | + | | [`SpikyKernel`](@ref) | $[0, 1h]$ | $1.5$ to $3.0$ | Academic | + | | [`LaguerreGaussKernel`](@ref) | $[0, 2h]$ | $1.3$ to $1.5$ | General | ++++ | +| [`ParabolicKernel`](@ref) | $[0, 1h]$ | $2.1$ to $3.0$ | TLSPH only | - | Any Kernel with a stability rating of more than '+++' doesn't suffer from pairing-instability. @@ -40,6 +41,7 @@ other_kernels = [ ("Gaussian", GaussianKernel{2}()), ("Poly6", Poly6Kernel{2}()), ("Laguerre-Gauss", LaguerreGaussKernel{2}()), + ("Parabolic", ParabolicKernel{2}()), ] spiky_kernel_group = [ diff --git a/docs/src/refs.bib b/docs/src/refs.bib index c479e59210..70fbd17fc5 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -342,6 +342,19 @@ @Article{Ganzenmueller2015 publisher = {Elsevier {BV}}, } +@Article{Ganzenmueller2015b, + author = {Ganzenm{\"{u}}ller, Georg C. and Hiermaier, Stefan and May, Michael}, + title = {On the similarity of meshless discretizations of {Peridynamics} and {Smooth-Particle} {Hydrodynamics}}, + journal = {Computers \& Structures}, + year = {2015}, + volume = {150}, + pages = {71--78}, + month = apr, + issn = {0045-7949}, + doi = {10.1016/j.compstruc.2014.12.011}, + publisher = {Elsevier {BV}}, +} + @Book{Gasser2021, author = {Gasser, T. Christian}, title = {Vascular Biomechanics: Concepts, Models, and Applications}, diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md index b6dc4e6f04..ae17586312 100644 --- a/docs/src/systems/total_lagrangian_sph.md +++ b/docs/src/systems/total_lagrangian_sph.md @@ -57,6 +57,57 @@ are the Lamé coefficients, where $E$ is the Young's modulus and $\nu$ is the Po The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGanzenmueller`](@ref). +### Bond-Associated Quadrature + +The traditional formulation above is selected with `model=StandardTLSPHModel()`. +With `model=BondAssociatedTLSPHModel()`, the corrected gradient is interpreted as the +`:C1` reproducing-kernel gradient weight used by the bond-associated correspondence +formulation in Peridynamics.jl. For +$\bm{\xi}_{ab} = \bm{X}_b - \bm{X}_a$ and particle volume +$V_b = m_{0b}/\rho_{0b}$, define +```math +\omega_{ab} = +-\frac{\nabla_{0a} W(\bm{X}_{ab}) \cdot \bm{X}_{ab}} + {\lVert\bm{X}_{ab}\rVert^2}, \qquad +w_a = \sum_b \omega_{ab} V_b, +``` +and +```math +\bm{\Phi}_{ab} = V_b \bm{L}_{0a} \nabla_{0a} W(\bm{X}_{ab}). +``` +These are algebraically the `:C1` moment matrix and gradient weights. + +The deformation gradient associated with the center of bond $ab$ is +```math +\bm{F}_{ab} += \frac{\bm{F}_a + \bm{F}_b}{2} ++ \left(\bm{x}_b - \bm{x}_a + - \frac{\bm{F}_a + \bm{F}_b}{2}\bm{\xi}_{ab}\right) + \frac{\bm{\xi}_{ab}^T}{\lVert\bm{\xi}_{ab}\rVert^2}. +``` +Let $\bm{P}_{ab}$ be the PK1 stress evaluated from $\bm{F}_{ab}$ with the material +parameters of particle $a$. The transverse stress integral is +```math +\bm{Q}_a = \sum_b \omega_{ab} +\left(\frac{1}{2w_a} + \frac{1}{2w_b}\right) V_b +\bm{P}_{ab} +\left(\bm{I} - \frac{\bm{\xi}_{ab}\bm{\xi}_{ab}^T} +{\lVert\bm{\xi}_{ab}\rVert^2}\right). +``` +The directed force state and acceleration are then +```math +\bm{t}_{ab} += \frac{\omega_{ab}}{w_a\lVert\bm{\xi}_{ab}\rVert^2} + \bm{P}_{ab}\bm{\xi}_{ab} ++ \frac{\bm{Q}_a\bm{\Phi}_{ab}}{V_b}, +\qquad +\frac{\mathrm{d}\bm{v}_a}{\mathrm{d}t} += \frac{1}{\rho_{0a}}\sum_b V_b(\bm{t}_{ab} - \bm{t}_{ba}). +``` +This is the bond-associated quadrature update implemented in +[`Peridynamics.jl`](https://github.com/kaipartmann/Peridynamics.jl/blob/main/src/physics/rk_correspondence.jl), +adapted to the existing TLSPH kernels and two- or three-dimensional particle systems. + ```@autodocs Modules = [TrixiParticles] Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "system.jl")] diff --git a/examples/structure/oscillating_beam_2d.jl b/examples/structure/oscillating_beam_2d.jl index 784692001c..ec19fae0be 100644 --- a/examples/structure/oscillating_beam_2d.jl +++ b/examples/structure/oscillating_beam_2d.jl @@ -67,6 +67,7 @@ structure_system = TotalLagrangianSPHSystem(structure; smoothing_kernel, smoothi acceleration=(0.0, -gravity), penalty_force=nothing, viscosity=nothing, clamped_particles_motion=nothing, + model=StandardTLSPHModel(), self_interaction_nhs=:default) # ========================================================================================== diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 0bc56370f2..0c345bfd37 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -67,6 +67,7 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition, apply_angular_velocity export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, + StandardTLSPHModel, BondAssociatedTLSPHModel, RigidBodySystem, WallBoundarySystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySystem, ImplicitIncompressibleSPHSystem @@ -82,7 +83,7 @@ export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechni ContinuityEquationTermSun2019, MomentumEquationTermSun2019, VelocityAveraging export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, - WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel + WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel, ParabolicKernel export StateEquationCole, StateEquationIdealGas, StateEquationAdaptiveCole export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, ViscosityMorrisSGS, ViscosityCarreauYasuda diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index 456d55a13f..47341ef946 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -820,3 +820,68 @@ end # C' = 1 / (4 * pi * h^3 * C) = C_ / (h^3) return C_ * h_inv^3 end + +@doc raw""" + ParabolicKernel{NDIMS}() + +Parabolic smoothing kernel given by +```math + W(r, h) = \frac{1}{h^d} w(r/h) +``` +with +```math +w(q) = \sigma \begin{cases} + 1 - q^2 & \text{if } 0 \leq q < 1, \\ + 0 & \text{if } q \geq 1, +\end{cases} +``` +where ``d`` is the number of dimensions and ``\sigma`` is a normalization factor. +The normalization factor ``\sigma`` is ``3 / 4`` in one dimension, ``2 / \pi`` +in two dimensions, and ``15 / (8 \pi)`` in three dimensions. + +This kernel function has a compact support of ``[0, h]``. + +The parabolic kernel is intended for Total Lagrangian SPH (TLSPH), where the +gradient is evaluated with respect to initial coordinates. With +``X_{ij} = X_i - X_j``, this kernel yields +```math + \nabla W_{ij} = -c X_{ij}, +``` +where ``c`` is constant for fixed smoothing length. This is the kernel that makes +the basic Peridynamics discretization equivalent to TLSPH, as shown by +[Ganzenmüller et al. (2015)](@cite Ganzenmueller2015b). + +!!! warning + This kernel is not useful for fluid SPH. Its simple gradient causes severe + particle clustering in fluid simulations. Use it only for TLSPH applications. + +For general information and usage see [Smoothing Kernels](@ref smoothing_kernel). +""" +struct ParabolicKernel{NDIMS} <: AbstractSmoothingKernel{NDIMS} end + +@inline function kernel_unsafe(kernel::ParabolicKernel, r::Real, h) + h_inv = 1 / h + q = r * h_inv + + return normalization_factor(kernel, h_inv) * (1 - q^2) +end + +@inline function kernel_deriv_div_r_unsafe(kernel::ParabolicKernel, r::Real, h) + h_inv = 1 / h + + return -2 * normalization_factor(kernel, h_inv) * h_inv^2 +end + +@inline compact_support(::ParabolicKernel, h) = h + +@inline function normalization_factor(::ParabolicKernel{1}, h_inv) + return oftype(h_inv, 3 / 4) * h_inv +end + +@inline function normalization_factor(::ParabolicKernel{2}, h_inv) + return oftype(h_inv, 2 / pi) * h_inv^2 +end + +@inline function normalization_factor(::ParabolicKernel{3}, h_inv) + return oftype(h_inv, 15 / (8 * pi)) * h_inv^3 +end diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index d3ee370237..dfdb804a2d 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -18,6 +18,13 @@ end # Function barrier without dispatch for unit testing @inline function interact_structure_structure!(dv, v_system, system, semi; eachparticle=each_integrated_particle(system)) + interact_structure_structure!(dv, v_system, system, tlsph_model(system), semi; + eachparticle) +end + +@inline function interact_structure_structure!(dv, v_system, system, + ::StandardTLSPHModel, semi; + eachparticle=each_integrated_particle(system)) (; penalty_force) = system # Everything here is done in the initial coordinates @@ -101,6 +108,95 @@ end return dv end +@inline function interact_structure_structure!(dv, v_system, system, + ::BondAssociatedTLSPHModel, semi; + eachparticle=each_integrated_particle(system)) + (; penalty_force, mass, material_density) = system + (; weighted_volume) = system.cache + + system_coords = initial_coordinates(system) + neighborhood_search = get_neighborhood_search(system, semi) + backend = semi.parallelization_backend + h = initial_smoothing_length(system) + almostzero = sqrt(eps(h^2)) + + @threaded semi for particle in eachparticle + m_a = @inbounds mass[particle] + rho_a = @inbounds material_density[particle] + volume_a = div_fast(m_a, rho_a) + current_coords_a = @inbounds current_coords(system, particle) + F_a = @inbounds deformation_gradient(system, particle) + stress_integral_a = @inbounds bond_associated_stress_integral(system, particle) + weighted_volume_a = @inbounds weighted_volume[particle] + + dv_particle = Ref(zero(current_coords_a)) + + @inbounds foreach_neighbor(system_coords, system_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + initial_distance < almostzero && return + + rho_b = material_density[neighbor] + m_b = mass[neighbor] + volume_b = div_fast(m_b, rho_b) + current_coords_b = current_coords(system, neighbor) + F_b = deformation_gradient(system, neighbor) + + current_pos_diff_ = current_coords_a - current_coords_b + current_pos_diff = convert.(eltype(system), current_pos_diff_) + current_distance = norm(current_pos_diff) + initial_bond = -initial_pos_diff + current_bond = -current_pos_diff + + F_ab = bond_deformation_gradient(F_a, F_b, initial_bond, current_bond, + initial_distance) + P_ab = pk1_stress_tensor(F_ab, system, particle) + P_ba = pk1_stress_tensor(F_ab, system, neighbor) + + grad_kernel_a = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + weight_a = bond_weight(grad_kernel_a, initial_pos_diff, initial_distance) + gradient_weight_a = volume_b * correction_matrix(system, particle) * + grad_kernel_a + + neighbor_pos_diff = -initial_pos_diff + grad_kernel_b = smoothing_kernel_grad_unsafe(system, neighbor_pos_diff, + initial_distance, neighbor) + weight_b = bond_weight(grad_kernel_b, neighbor_pos_diff, initial_distance) + gradient_weight_b = volume_a * correction_matrix(system, neighbor) * + grad_kernel_b + + stress_integral_b = bond_associated_stress_integral(system, neighbor) + weighted_volume_b = weighted_volume[neighbor] + + force_state_ab = weight_a / (weighted_volume_a * initial_distance^2) * + (P_ab * initial_bond) + + stress_integral_a * gradient_weight_a / volume_b + force_state_ba = weight_b / (weighted_volume_b * initial_distance^2) * + (P_ba * -initial_bond) + + stress_integral_b * gradient_weight_b / volume_a + + dv_particle[] += volume_b / rho_a * (force_state_ab - force_state_ba) + + @inbounds dv_penalty_force!(dv_particle, penalty_force, particle, neighbor, + initial_pos_diff, initial_distance, + current_pos_diff, current_distance, + system, m_a, m_b, rho_a, rho_b, F_a, F_b) + + @inbounds dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor, + current_pos_diff, current_distance, + m_a, m_b, rho_a, rho_b, F_a, grad_kernel_a) + end + + for i in 1:ndims(system) + @inbounds dv[i, particle] += dv_particle[][i] + end + end + + return dv +end + # Structure-fluid interaction function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index fbd4dc0333..b0821120f5 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -1,6 +1,24 @@ +abstract type AbstractTLSPHModel end + +""" + StandardTLSPHModel() + +Standard total Lagrangian SPH correspondence formulation. +""" +struct StandardTLSPHModel <: AbstractTLSPHModel end + +""" + BondAssociatedTLSPHModel() + +Total Lagrangian SPH correspondence formulation with bond-associated quadrature +integration at the center of each bond. +""" +struct BondAssociatedTLSPHModel <: AbstractTLSPHModel end + @doc raw""" TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing_length, young_modulus, poisson_ratio, + model=StandardTLSPHModel(), clamped_particles=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), @@ -26,6 +44,10 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. See [Smoothing Kernels](@ref smoothing_kernel). - `young_modulus`: Young's modulus. - `poisson_ratio`: Poisson ratio. +- `model`: Discretization of the correspondence formulation. + [`StandardTLSPHModel`](@ref) (default) uses the traditional TLSPH + discretization. [`BondAssociatedTLSPHModel`](@ref) uses + bond-associated quadrature integration. - `clamped_particles`: Indices specifying the clamped particles that are fixed and not integrated to clamp the structure. - `clamped_particles_motion`: Prescribed motion of the clamped particles. @@ -74,7 +96,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. `structure`, so their indices are `1:nparticles(clamped_particles)`. """ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D, - YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS, VA, + YM, PR, LL, LM, K, MOD, PF, V, ST, M, IM, NHS, VA, C} <: AbstractStructureSystem{NDIMS} initial_condition :: IC initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle] @@ -94,6 +116,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, smoothing_length :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} boundary_model :: BM + model :: MOD penalty_force :: PF viscosity :: V source_terms :: ST @@ -105,7 +128,9 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, end function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing_length, - young_modulus, poisson_ratio, clamped_particles=Int[], + young_modulus, poisson_ratio, + model::AbstractTLSPHModel=StandardTLSPHModel(), + clamped_particles=Int[], clamped_particles_motion=nothing, acceleration=ntuple(_ -> zero(eltype(initial_condition)), ndims(smoothing_kernel)), @@ -164,7 +189,8 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing initialize_prescribed_motion!(clamped_particles_motion, initial_condition_sorted, n_clamped_particles) - cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)..., + cache = (; create_cache_tlsph_model(model, initial_condition_sorted)..., + create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)..., create_cache_tlsph(velocity_averaging, initial_condition_sorted)...) return TotalLagrangianSPHSystem(initial_condition_sorted, initial_coordinates, @@ -174,7 +200,7 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing poisson_ratio_sorted, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, - penalty_force, viscosity, source_terms, + model, penalty_force, viscosity, source_terms, clamped_particles_motion, ismoving, self_interaction_nhs, velocity_averaging, cache) end @@ -234,7 +260,7 @@ function initialize_self_interaction_nhs(system::TotalLagrangianSPHSystem, system.poisson_ratio, system.lame_lambda, system.lame_mu, system.smoothing_kernel, system.smoothing_length, system.acceleration, - system.boundary_model, system.penalty_force, + system.boundary_model, system.model, system.penalty_force, system.viscosity, system.source_terms, system.clamped_particles_motion, system.clamped_particles_moving, @@ -247,6 +273,19 @@ extract_periodic_box(nhs) = nhs.periodic_box create_cache_tlsph(::Nothing, initial_condition) = (;) +create_cache_tlsph_model(::StandardTLSPHModel, initial_condition) = (;) + +function create_cache_tlsph_model(::BondAssociatedTLSPHModel, initial_condition) + n_particles = nparticles(initial_condition) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + weighted_volume = Vector{ELTYPE}(undef, n_particles) + stress_integral = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) + + return (; weighted_volume, stress_integral) +end + function create_cache_tlsph(::PrescribedMotion, initial_condition) velocity = zero(initial_condition.velocity) acceleration = zero(initial_condition.velocity) @@ -358,6 +397,10 @@ end extract_smatrix(system.pk1_rho2, system, particle) end +@propagate_inbounds function bond_associated_stress_integral(system, particle) + extract_smatrix(system.cache.stress_integral, system, particle) +end + @propagate_inbounds function young_modulus(system::TotalLagrangianSPHSystem, particle) return young_modulus(system, system.young_modulus, particle) end @@ -388,6 +431,9 @@ end return system.velocity_averaging end +@inline tlsph_model(system::TotalLagrangianSPHSystem) = system.model +@inline tlsph_model(system) = StandardTLSPHModel() + function initialize!(system::TotalLagrangianSPHSystem, semi) (; correction_matrix) = system @@ -398,6 +444,41 @@ function initialize!(system::TotalLagrangianSPHSystem, semi) # Calculate correction matrix compute_gradient_correction_matrix!(correction_matrix, system, initial_coords, density_fun, semi) + initialize_tlsph_model!(system, system.model, semi) +end + +@inline initialize_tlsph_model!(system, ::StandardTLSPHModel, semi) = system + +function initialize_tlsph_model!(system, ::BondAssociatedTLSPHModel, semi) + (; weighted_volume) = system.cache + (; mass, material_density) = system + + initial_coords = initial_coordinates(system) + neighborhood_search = get_neighborhood_search(system, semi) + backend = semi.parallelization_backend + h = initial_smoothing_length(system) + almostzero = sqrt(eps(h^2)) + + @threaded semi for particle in eachparticle(system) + result = Ref(zero(eltype(system))) + + @inbounds foreach_neighbor(initial_coords, initial_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + initial_distance < almostzero && return + + grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + volume = div_fast(mass[neighbor], material_density[neighbor]) + result[] += bond_weight(grad_kernel, initial_pos_diff, initial_distance) * + volume + end + + @inbounds weighted_volume[particle] = result[] + end + + return system end function update_positions!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t) @@ -452,7 +533,29 @@ end function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t) # Precompute PK1 stress tensor - @trixi_timeit timer() "stress tensor" compute_pk1_corrected!(system, semi) + @trixi_timeit timer() "stress tensor" compute_stress!(system, system.model, semi) + + return system +end + +@inline compute_stress!(system, ::StandardTLSPHModel, semi) = + compute_pk1_corrected!(system, semi) + +@inline function compute_stress!(system, ::BondAssociatedTLSPHModel, semi) + (; deformation_grad, pk1_rho2, material_density) = system + + calc_deformation_grad!(deformation_grad, system, semi) + + @threaded semi for particle in eachparticle(system) + pk1_particle = @inbounds pk1_stress_tensor(system, particle) + rho2_inv = 1 / @inbounds material_density[particle]^2 + + for j in 1:ndims(system), i in 1:ndims(system) + @inbounds pk1_rho2[i, j, particle] = pk1_particle[i, j] * rho2_inv + end + end + + calc_bond_associated_stress_integral!(system, semi) return system end @@ -556,11 +659,80 @@ end return deformation_grad end +@inline function bond_weight(grad_kernel, initial_pos_diff, initial_distance) + return -dot(grad_kernel, initial_pos_diff) / initial_distance^2 +end + +@inline function bond_deformation_gradient(F_a, F_b, initial_bond, current_bond, + initial_distance) + F_average = (F_a + F_b) / 2 + F_correction = (current_bond - F_average * initial_bond) * + (initial_bond' / initial_distance^2) + return F_average + F_correction +end + +function calc_bond_associated_stress_integral!(system, semi) + (; stress_integral, weighted_volume) = system.cache + (; mass, material_density) = system + + initial_coords = initial_coordinates(system) + neighborhood_search = get_neighborhood_search(system, semi) + backend = semi.parallelization_backend + h = initial_smoothing_length(system) + almostzero = sqrt(eps(h^2)) + + @threaded semi for particle in eachparticle(system) + F_a = @inbounds deformation_gradient(system, particle) + current_coords_a = @inbounds current_coords(system, particle) + weighted_volume_a = @inbounds weighted_volume[particle] + result = Ref(zero(F_a)) + + @inbounds foreach_neighbor(initial_coords, initial_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + initial_distance < almostzero && return + + grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + weight = bond_weight(grad_kernel, initial_pos_diff, initial_distance) + volume_b = div_fast(mass[neighbor], material_density[neighbor]) + + initial_bond = -initial_pos_diff + current_coords_b = current_coords(system, neighbor) + current_bond_ = current_coords_b - current_coords_a + current_bond = convert.(eltype(system), current_bond_) + + F_b = deformation_gradient(system, neighbor) + F_ab = bond_deformation_gradient(F_a, F_b, initial_bond, current_bond, + initial_distance) + P_ab = pk1_stress_tensor(F_ab, system, particle) + transverse_projection = I - + initial_bond * initial_bond' / initial_distance^2 + weighted_volume_b = weighted_volume[neighbor] + quadrature_weight = weight * + (inv(2 * weighted_volume_a) + + inv(2 * weighted_volume_b)) * volume_b + + result[] += quadrature_weight * (P_ab * transverse_projection) + end + + for j in 1:ndims(system), i in 1:ndims(system) + @inbounds stress_integral[i, j, particle] = result[][i, j] + end + end + + return stress_integral +end + # First Piola-Kirchhoff stress tensor @propagate_inbounds function pk1_stress_tensor(system, particle) - (; lame_lambda, lame_mu) = system - F = deformation_gradient(system, particle) + return pk1_stress_tensor(F, system, particle) +end + +@propagate_inbounds function pk1_stress_tensor(F, system, particle) + (; lame_lambda, lame_mu) = system S = pk2_stress_tensor(F, lame_lambda, lame_mu, particle) return F * S @@ -640,6 +812,10 @@ end # The von-Mises stress is one form of equivalent stress, where sigma is the deviatoric stress. # See pages 32 and 123. function von_mises_stress(system) + return von_mises_stress(system, system.model) +end + +function von_mises_stress(system, ::StandardTLSPHModel) von_mises_stress_vector = zeros(eltype(system.pk1_rho2), nparticles(system)) @threaded default_backend(von_mises_stress_vector) for particle in @@ -650,6 +826,20 @@ function von_mises_stress(system) return von_mises_stress_vector end +function von_mises_stress(system, ::BondAssociatedTLSPHModel) + cauchy_stress_tensors = cauchy_stress(system) + von_mises_stress_vector = zeros(eltype(system.pk1_rho2), nparticles(system)) + + @threaded default_backend(von_mises_stress_vector) for particle in + each_integrated_particle(system) + sigma = extract_smatrix(cauchy_stress_tensors, system, particle) + s = sigma - (1.0 / 3.0) * tr(sigma) * I + von_mises_stress_vector[particle] = sqrt(3.0 / 2.0 * sum(s .^ 2)) + end + + return von_mises_stress_vector +end + # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl # with `@batch` (`@threaded`). # Otherwise, `@threaded` does not work here with Julia ARM on macOS. @@ -671,6 +861,10 @@ end # See here page 473 for the relation between the `pk1`, the first Piola-Kirchhoff tensor, # and the Cauchy stress. function cauchy_stress(system::TotalLagrangianSPHSystem) + return cauchy_stress(system, system.model) +end + +function cauchy_stress(system, ::StandardTLSPHModel) NDIMS = ndims(system) cauchy_stress_tensors = zeros(eltype(system.pk1_rho2), NDIMS, NDIMS, @@ -688,6 +882,55 @@ function cauchy_stress(system::TotalLagrangianSPHSystem) return cauchy_stress_tensors end +function cauchy_stress(system, ::BondAssociatedTLSPHModel) + NDIMS = ndims(system) + cauchy_stress_tensors = zeros(eltype(system.pk1_rho2), NDIMS, NDIMS, + nparticles(system)) + (; mass, material_density) = system + (; weighted_volume) = system.cache + + initial_coords = initial_coordinates(system) + neighborhood_search = system.self_interaction_nhs + backend = default_backend(cauchy_stress_tensors) + h = initial_smoothing_length(system) + almostzero = sqrt(eps(h^2)) + + @threaded backend for particle in each_integrated_particle(system) + F_a = @inbounds deformation_gradient(system, particle) + current_coords_a = @inbounds current_coords(system, particle) + weighted_volume_a = @inbounds weighted_volume[particle] + sigma_a = Ref(zero(F_a)) + + @inbounds foreach_neighbor(initial_coords, initial_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + initial_distance < almostzero && return + + grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + weight = bond_weight(grad_kernel, initial_pos_diff, initial_distance) + volume_b = div_fast(mass[neighbor], material_density[neighbor]) + initial_bond = -initial_pos_diff + current_bond_ = current_coords(system, neighbor) - current_coords_a + current_bond = convert.(eltype(system), current_bond_) + F_b = deformation_gradient(system, neighbor) + F_ab = bond_deformation_gradient(F_a, F_b, initial_bond, current_bond, + initial_distance) + P_ab = pk1_stress_tensor(F_ab, system, particle) + sigma_ab = P_ab * F_ab' / det(F_ab) + + sigma_a[] += weight / weighted_volume_a * volume_b * sigma_ab + end + + for j in 1:NDIMS, i in 1:NDIMS + @inbounds cauchy_stress_tensors[i, j, particle] = sigma_a[][i, j] + end + end + + return cauchy_stress_tensors +end + function calculate_dt(v_ode, u_ode, cfl_number, system::TotalLagrangianSPHSystem, semi) # TODO variable smoothing length smoothing_length_ = initial_smoothing_length(system) @@ -759,6 +1002,7 @@ function Base.show(io::IO, system::TotalLagrangianSPHSystem) print(io, "", system.smoothing_kernel) print(io, ", ", system.acceleration) print(io, ", ", system.boundary_model) + show_tlsph_model(io, system.model) print(io, ", ", system.penalty_force) print(io, ", ", system.viscosity) print(io, ") with ", nparticles(system), " particles") @@ -790,12 +1034,24 @@ function Base.show(io::IO, ::MIME"text/plain", system::TotalLagrangianSPHSystem) summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "acceleration", system.acceleration) summary_line(io, "boundary model", system.boundary_model) + show_tlsph_model_summary(io, system.model) summary_line(io, "penalty force", system.penalty_force) summary_line(io, "viscosity", system.viscosity) summary_footer(io) end end +@inline show_tlsph_model(io, ::StandardTLSPHModel) = nothing +@inline show_tlsph_model_summary(io, ::StandardTLSPHModel) = nothing + +@inline function show_tlsph_model(io, model::BondAssociatedTLSPHModel) + print(io, ", ", model) +end + +@inline function show_tlsph_model_summary(io, model::BondAssociatedTLSPHModel) + summary_line(io, "model", model |> typeof |> nameof) +end + function check_configuration(system::TotalLagrangianSPHSystem, systems, nhs) (; boundary_model) = system diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index e112d9562b..c5d0137295 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -4,38 +4,40 @@ # All smoothing kernels should integrate to something close to 1. # We integrate slightly beyond the compact support to verify that the kernel is # correctly evaluating to zero there. - function integrate_kernel_2d(smk) + function integrate_kernel_2d(kernel, h) integral_2d_radial, - _ = quadgk(r -> r * TrixiParticles.kernel(smk, r, 1.0), 0, - TrixiParticles.compact_support(smk, 1.0) * 1.1, + _ = quadgk(r -> r * TrixiParticles.kernel(kernel, r, h), 0, + TrixiParticles.compact_support(kernel, h) * 1.1, rtol=1e-15) return 2 * pi * integral_2d_radial end - function integrate_kernel_1d(smk) + function integrate_kernel_1d(kernel, h) integral_1d_half, - _ = quadgk(r -> TrixiParticles.kernel(smk, r, 1.0), 0, - TrixiParticles.compact_support(smk, 1.0) * 1.1, + _ = quadgk(r -> TrixiParticles.kernel(kernel, r, h), 0, + TrixiParticles.compact_support(kernel, h) * 1.1, rtol=1e-15) return 2 * integral_1d_half end - function integrate_kernel_3d(smk) + function integrate_kernel_3d(kernel, h) integral_3d_radial, - _ = quadgk(r -> r^2 * TrixiParticles.kernel(smk, r, 1.0), 0, - TrixiParticles.compact_support(smk, 1.0) * 1.1, + _ = quadgk(r -> r^2 * TrixiParticles.kernel(kernel, r, h), 0, + TrixiParticles.compact_support(kernel, h) * 1.1, rtol=1e-15) return 4 * pi * integral_3d_radial end # Treat the truncated Gaussian kernel separately @testset "GaussianKernel" begin - error_2d = abs(integrate_kernel_2d(GaussianKernel{2}()) - 1.0) - error_3d = abs(integrate_kernel_3d(GaussianKernel{3}()) - 1.0) + for h in [0.5, 1.0] + error_2d = abs(integrate_kernel_2d(GaussianKernel{2}(), h) - 1.0) + error_3d = abs(integrate_kernel_3d(GaussianKernel{3}(), h) - 1.0) - # Large error due to truncation - @test 1e-4 < error_2d < 1e-3 - @test 1e-4 < error_3d < 1e-3 + # Large error due to truncation + @test 1e-4 < error_2d < 1e-3 + @test 1e-4 < error_3d < 1e-3 + end end # All other kernels @@ -48,27 +50,31 @@ WendlandC6Kernel, SpikyKernel, Poly6Kernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] kernels_1d = [ SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] @testset "$kernel" for kernel in kernels - # The integral should be 1 for all kernels - error_2d = abs(integrate_kernel_2d(kernel{2}()) - 1.0) - error_3d = abs(integrate_kernel_3d(kernel{3}()) - 1.0) + for h in [0.5, 1.0] + # The integral should be 1 for all kernels + error_2d = abs(integrate_kernel_2d(kernel{2}(), h) - 1.0) + error_3d = abs(integrate_kernel_3d(kernel{3}(), h) - 1.0) - @test error_2d <= 1e-15 - @test error_3d <= 1e-15 + @test error_2d <= 2e-15 + @test error_3d <= 3e-15 - if kernel in kernels_1d - error_1d = abs(integrate_kernel_1d(kernel{1}()) - 1.0) - @test error_1d <= 1e-15 + if kernel in kernels_1d + error_1d = abs(integrate_kernel_1d(kernel{1}(), h) - 1.0) + @test error_1d <= 1e-15 + end end end end @@ -86,14 +92,16 @@ WendlandC6Kernel, SpikyKernel, Poly6Kernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] kernels_1d = [ SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] # Test 4 different smoothing lengths @@ -132,6 +140,33 @@ end end + @testset "ParabolicKernel Gradient Linearity" begin + # The `ParabolicKernel` is designed to have a linear gradient: ∇Wᵢⱼ = -c * rᵢⱼ. + smoothing_lengths = (0.5, 1.0) + + for ndims in 1:3 + kernel_ = ParabolicKernel{ndims}() + + @testset "ParabolicKernel{$ndims}" begin + for h in smoothing_lengths + h_inv = inv(h) + c = 2 * TrixiParticles.normalization_factor(kernel_, h_inv) * h_inv^2 + + for factor in (0.25, 0.5, 0.75) + pos_diff = SVector(ntuple(i -> i * factor * h / (2 * ndims), + ndims)) + distance = norm(pos_diff) + + gradient = TrixiParticles.kernel_grad(kernel_, pos_diff, + distance, h) + + @test gradient ≈ -c .* pos_diff + end + end + end + end + end + @testset verbose=false "Return Type" begin # Test that the return type of the kernel and kernel derivative preserve # the input type. We don't want to return `Float64` when working with `Float32`. @@ -145,14 +180,16 @@ WendlandC6Kernel, SpikyKernel, Poly6Kernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] kernels_1d = [ SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, - LaguerreGaussKernel + LaguerreGaussKernel, + ParabolicKernel ] # Test different smoothing length types diff --git a/test/schemes/structure/total_lagrangian_sph/rhs.jl b/test/schemes/structure/total_lagrangian_sph/rhs.jl index e0abb48e4a..598756a8be 100644 --- a/test/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/test/schemes/structure/total_lagrangian_sph/rhs.jl @@ -202,4 +202,179 @@ end end end + + @testset verbose=true "Bond-Associated Quadrature" begin + coordinate_range = 0.1:0.1:0.5 + coordinates = hcat(collect.(Iterators.product(coordinate_range, + coordinate_range))...) + n_particles = size(coordinates, 2) + mass = fill(10.0, n_particles) + density = fill(1000.0, n_particles) + young_modulus = collect(range(2.0, 3.0, length=n_particles)) + poisson_ratio = fill(0.25, n_particles) + initial_condition = InitialCondition(; coordinates, mass, density) + system = TotalLagrangianSPHSystem(initial_condition; + smoothing_kernel=SchoenbergCubicSplineKernel{2}(), + smoothing_length=0.12, young_modulus, + poisson_ratio, + model=BondAssociatedTLSPHModel()) + semi = Semidiscretization(system, parallelization_backend=SerialBackend()) + system = semi.systems[1] + ode = semidiscretize(semi, (0.0, 1.0)) + v_ode = copy(ode.u0.x[1]) + u_ode = copy(ode.u0.x[2]) + u = TrixiParticles.wrap_u(u_ode, system, semi) + + for particle in axes(u, 2) + x, y = coordinates[:, particle] + u[:, particle] .= (1.05 * x + 0.1 * y + 0.03 * x * y, + -0.04 * x + 0.95 * y + 0.02 * x^2) + end + + dv_ode = zero(v_ode) + TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) + dv = TrixiParticles.wrap_v(dv_ode, system, semi) + + # Reconstruct the C1 reproducing-kernel quantities from Peridynamics.jl directly. + volumes = system.mass ./ system.material_density + weights = zeros(n_particles, n_particles) + gradient_weights = zeros(2, n_particles, n_particles) + weighted_volume_reference = zeros(n_particles) + moment_matrices = zeros(2, 2, n_particles) + system_coords = TrixiParticles.initial_coordinates(system) + neighborhood_search = TrixiParticles.get_neighborhood_search(system, semi) + + for particle in TrixiParticles.eachparticle(system) + TrixiParticles.foreach_neighbor(system_coords, system_coords, + neighborhood_search, SerialBackend(), + particle) do _, neighbor, initial_pos_diff, + initial_distance + initial_distance < sqrt(eps(system.smoothing_length^2)) && return + + grad_kernel = TrixiParticles.smoothing_kernel_grad_unsafe(system, + initial_pos_diff, + initial_distance, + particle) + initial_bond = -initial_pos_diff + weight = -dot(grad_kernel, initial_pos_diff) / initial_distance^2 + weights[neighbor, particle] = weight + weighted_volume_reference[particle] += weight * volumes[neighbor] + moment_matrices[:, :, particle] .+= weight * volumes[neighbor] * + initial_bond * initial_bond' + end + + moment_matrix_inv = inv(moment_matrices[:, :, particle]) + for neighbor in TrixiParticles.eachparticle(system) + initial_bond = coordinates[:, neighbor] - coordinates[:, particle] + gradient_weights[:, neighbor, particle] .= weights[neighbor, particle] * + volumes[neighbor] * + moment_matrix_inv * initial_bond + end + end + + deformation_gradients = zeros(2, 2, n_particles) + for particle in TrixiParticles.eachparticle(system) + deformation_gradients[:, :, particle] .= Matrix{Float64}(I, 2, 2) + for neighbor in TrixiParticles.eachparticle(system) + initial_bond = coordinates[:, neighbor] - coordinates[:, particle] + current_bond = u[:, neighbor] - u[:, particle] + displacement_bond = current_bond - initial_bond + deformation_gradients[:, :, particle] .+= + displacement_bond * gradient_weights[:, neighbor, particle]' + end + end + + stress_integrals = zeros(2, 2, n_particles) + bond_stresses = zeros(2, 2, n_particles, n_particles) + for particle in TrixiParticles.eachparticle(system) + for neighbor in TrixiParticles.eachparticle(system) + weight = weights[neighbor, particle] + iszero(weight) && continue + + initial_bond = coordinates[:, neighbor] - coordinates[:, particle] + current_bond = u[:, neighbor] - u[:, particle] + initial_distance = norm(initial_bond) + F_average = (deformation_gradients[:, :, particle] + + deformation_gradients[:, :, neighbor]) / 2 + F_bond = F_average + + (current_bond - F_average * initial_bond) * + (initial_bond' / initial_distance^2) + P_bond = TrixiParticles.pk1_stress_tensor(F_bond, system, particle) + bond_stresses[:, :, neighbor, particle] .= P_bond + transverse_projection = I - + initial_bond * initial_bond' / initial_distance^2 + quadrature_weight = weight * + (1 / (2 * weighted_volume_reference[particle]) + + 1 / (2 * weighted_volume_reference[neighbor])) * + volumes[neighbor] + stress_integrals[:, :, particle] .+= + quadrature_weight * P_bond * transverse_projection + end + end + + @test system.cache.weighted_volume ≈ weighted_volume_reference + for particle in TrixiParticles.eachparticle(system) + @test TrixiParticles.deformation_gradient(system, particle) ≈ + deformation_gradients[:, :, particle] + @test TrixiParticles.bond_associated_stress_integral(system, particle) ≈ + stress_integrals[:, :, particle] + end + + # Evaluate the directed-bond force updates from Peridynamics.jl directly. + force_density = zeros(size(dv)) + for particle in TrixiParticles.eachparticle(system) + for neighbor in TrixiParticles.eachparticle(system) + weight = weights[neighbor, particle] + iszero(weight) && continue + + initial_bond = coordinates[:, neighbor] - coordinates[:, particle] + initial_distance = norm(initial_bond) + force_state = weight / + (weighted_volume_reference[particle] * initial_distance^2) * + (bond_stresses[:, :, neighbor, particle] * initial_bond) + + stress_integrals[:, :, particle] * + gradient_weights[:, neighbor, particle] / volumes[neighbor] + + force_density[:, particle] += force_state * volumes[neighbor] + force_density[:, neighbor] -= force_state * volumes[particle] + end + end + + dv_reference = force_density ./ reshape(system.material_density, 1, :) + @test dv ≈ dv_reference rtol=2e-13 atol=2e-13 + @test vec(sum(dv .* reshape(system.mass, 1, :), dims=2)) ≈ zeros(2) atol=1e-13 + + gpu_system = TotalLagrangianSPHSystem(initial_condition; + smoothing_kernel=SchoenbergCubicSplineKernel{2}(), + smoothing_length=0.12, young_modulus, + poisson_ratio, + model=BondAssociatedTLSPHModel()) + gpu_backend = TrixiParticles.KernelAbstractions.CPU() + gpu_semi = Semidiscretization(gpu_system; parallelization_backend=gpu_backend) + gpu_ode = semidiscretize(gpu_semi, (0.0, 1.0)) + gpu_v_ode = gpu_ode.u0.x[1] + gpu_u = zeros(size(coordinates)) + for particle in axes(gpu_u, 2) + x, y = coordinates[:, particle] + gpu_u[:, particle] .= (1.05 * x + 0.1 * y + 0.03 * x * y, + -0.04 * x + 0.95 * y + 0.02 * x^2) + end + gpu_u_ode = vec(gpu_u) + gpu_dv_ode = zero(gpu_v_ode) + TrixiParticles.kick!(gpu_dv_ode, gpu_v_ode, gpu_u_ode, gpu_ode.p, 0.0) + gpu_dv = TrixiParticles.wrap_v(gpu_dv_ode, gpu_semi.systems[1], gpu_semi) + @test gpu_dv ≈ dv_reference rtol=2e-13 atol=2e-13 + + # A uniform deformation must be reproduced exactly in the interior. + affine_map = [1.1 0.2; -0.1 0.9] + for particle in axes(u, 2) + u[:, particle] .= affine_map * coordinates[:, particle] + end + TrixiParticles.update_tlsph_positions!(system, u, semi) + TrixiParticles.compute_stress!(system, system.model, semi) + @test TrixiParticles.deformation_gradient(system, 13) ≈ affine_map + P = TrixiParticles.pk1_stress_tensor(affine_map, system, 13) + sigma = TrixiParticles.cauchy_stress(system)[:, :, 13] + @test sigma ≈ P * affine_map' / det(affine_map) + end end; diff --git a/test/systems/tlsph_system.jl b/test/systems/tlsph_system.jl index b39f0a6b2a..321f93415e 100644 --- a/test/systems/tlsph_system.jl +++ b/test/systems/tlsph_system.jl @@ -43,6 +43,15 @@ @test TrixiParticles.initial_smoothing_length(system) == smoothing_length @test system.acceleration == [0.0 for _ in 1:NDIMS] @test system.boundary_model == boundary_model + @test system.model isa StandardTLSPHModel + + baq_system = TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, young_modulus=E, + poisson_ratio=nu, boundary_model, + model=BondAssociatedTLSPHModel()) + @test baq_system.model isa BondAssociatedTLSPHModel + @test size(baq_system.cache.weighted_volume) == (2,) + @test size(baq_system.cache.stress_integral) == (NDIMS, NDIMS, 2) end end diff --git a/validation/oscillating_beam_2d/plot_oscillating_beam_results.jl b/validation/oscillating_beam_2d/plot_oscillating_beam_results.jl index 1e24c633f0..1d0e687a91 100644 --- a/validation/oscillating_beam_2d/plot_oscillating_beam_results.jl +++ b/validation/oscillating_beam_2d/plot_oscillating_beam_results.jl @@ -1,81 +1,122 @@ include("../validation_util.jl") -# Activate for interactive plot +# Activate for interactive plots # using GLMakie using CairoMakie using CSV using DataFrames using JSON using Glob -using Printf using TrixiParticles -elastic_plate = (length=0.35, thickness=0.02) +const RESOLUTIONS = (5, 9, 17, 33, 65) +const COLORS = ("#ff0000", "#f5a000", "#800080", "#3caf75") +const TIME_LIMITS = (0.35, 0.55) +const DEFLECTION_LIMITS = (-0.235, -0.105) -# 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", "out") -merged_files = vcat(reference_files, simulation_files) -input_files = sort(merged_files, by=extract_number_from_filename) - -# Regular expressions for matching keys -key_pattern_x = r"deflection_x_structure_\d+" -key_pattern_y = r"deflection_y_structure_\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") +function load_tip_deflection(file_name) json_data = JSON.parsefile(file_name) + data = json_data["deflection_y_structure_1"] - resolution = parse(Int, split(split(file_name, "_")[end], ".")[1]) - particle_spacing_ = elastic_plate.thickness / (resolution - 1) + return Float64.(data["time"]), Float64.(data["values"]) +end - 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)))) +function file_for_resolution(files, resolution) + index = findfirst(file -> extract_number_from_filename(file) == resolution, files) + return isnothing(index) ? nothing : files[index] +end - if isempty(matching_keys_x) - error("No matching keys found in: $file_name") +reference_directory = joinpath(validation_dir(), "oscillating_beam_2d") +reference_files = glob("validation_reference_*.json", reference_directory) +simulation_files = glob("validation_run_oscillating_beam_2d_*.json", "out") +turek = CSV.read(joinpath(reference_directory, "reference_turek.csv"), DataFrame) + +set_theme!(Theme(fontsize=28, fonts=(; regular="TeX Gyre Termes"))) + +fig = Figure(size=(1080, 760), figure_padding=(20, 25, 15, 15)) +ax = Axis(fig[1, 1]; + xlabel="Time (s)", + ylabel="Tip Y - Deflection (m)", + limits=(TIME_LIMITS, DEFLECTION_LIMITS), + xticks=0.35:0.05:0.55, + yticks=-0.13:0.01:-0.11, + xlabelsize=34, + ylabelsize=34, + xticklabelsize=30, + yticklabelsize=30, + spinewidth=2.5, + xtickwidth=2.5, + ytickwidth=2.5, + xticksize=10, + yticksize=10) + +line_plots = [] +marker_plots = [] + +for (resolution, color) in zip(RESOLUTIONS, COLORS) + reference_file = file_for_resolution(reference_files, resolution) + isnothing(reference_file) && + error("Reference data for resolution $resolution not found") + + reference_time, reference_deflection = load_tip_deflection(reference_file) + marker_plot = scatter!(ax, reference_time, reference_deflection; + color=:transparent, + strokecolor=color, + strokewidth=3, + marker=:circle, + markersize=15) + push!(marker_plots, marker_plot) + + simulation_file = file_for_resolution(simulation_files, resolution) + if isnothing(simulation_file) + @warn "Simulation data for resolution $resolution not found; plotting the reference data as the solid curve" + simulation_time = reference_time + simulation_deflection = reference_deflection + else + simulation_time, simulation_deflection = load_tip_deflection(simulation_file) 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) - end - end + line_plot = lines!(ax, simulation_time, simulation_deflection; + color, linewidth=4) + push!(line_plots, line_plot) 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) - -# save("validation_osc_beam.png", fig) +turek_plot = lines!(ax, turek.time, turek.Uy; + color=:black, linestyle=:dash, linewidth=3) + +# Four empty entries make the two reference labels occupy the first two rows +# of the right legend column. +empty_entry = LineElement(color=:transparent) +legend_elements = [line_plots..., + MarkerElement(marker=:circle, color=:transparent, + strokecolor=:black, strokewidth=3, + markersize=15), + turek_plot, + empty_entry, + empty_entry] +legend_labels = [L"t_s/dp = 8", + L"t_s/dp = 16", + L"t_s/dp = 32", + L"t_s/dp = 64", + "O'Connor and Rogers, 2021", + "Turek & Hron, 2007", + "", + ""] + +Legend(fig[1, 1], legend_elements, legend_labels; + nbanks=2, + tellwidth=false, + tellheight=false, + halign=:center, + valign=:top, + margin=(20, 20, 8, 8), + padding=(8, 8, 5, 5), + patchsize=(35, 20), + rowgap=1, + colgap=16, + framevisible=true, + framewidth=1, + backgroundcolor=:white) + +# save("validation_oscillating_beam_2d.png", fig, px_per_unit=2) fig