Merge branch 'main' into fsi-check-configuration

.gitignore

@@ -18,3 +18,4 @@ run
 *.vtu
 *.pvd
 *.mp4
+*.csv

docs/src/install.md

@@ -30,11 +30,28 @@ The advantage of using a separate `run` directory is that you can also add other
 related packages (e.g., OrdinaryDiffEq.jl, see above) to the project in the `run` folder
 and always have a reproducible environment at hand to share with others.
 
-## Optional Software/Packages
-- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- a Julia package of ordinary differential equation solvers that is used in the examples
+## Optional software/packages
+- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- A Julia package of ordinary differential equation solvers that is used in the examples
 - [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia Plotting library that is used in some examples
 - [PythonPlot.jl](https://github.com/JuliaPy/PythonPlot.jl) -- Plotting library that can be used instead of Plots.jl
 - [ParaView](https://www.paraview.org/) -- Software that can be used for visualization of results
 
 ## Common issues
 
+If you followed the [installation instructions for developers](@ref for-developers) and you
+run into any problems with packages when pulling the latest version of TrixiParticles.jl,
+start Julia with the project in the `run` folder,
+```bash
+   julia --project=run
+```
+update all packages in that project, resolve all conflicts in the project, and install all
+new dependencies:
+```julia
+julia> using Pkg
+
+julia> Pkg.update()
+
+julia> Pkg.resolve()
+
+julia> Pkg.instantiate()
+```

examples/fluid/oscillating_drop_2d.jl

@@ -0,0 +1,93 @@
+# Oscillating drop in a central force field, based on
+#
+# J. J. Monaghan, Ashkan Rafiee.
+# "A Simple SPH Algorithm for Multi-Fluid Flow with High Density Ratios."
+# In: International Journal for Numerical Methods in Fluids 71, no. 5 (2013), pages 537-61.
+# https://doi.org/10.1002/fld.3671.
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+fluid_particle_spacing = 0.05
+
+# ==========================================================================================
+# ==== Experiment Setup
+radius = 1.0
+sigma = 0.5
+# Make this a constant because global variables in the source terms are slow
+const OMEGA = 1.0
+
+source_terms = (coords, velocity, density, pressure) -> -OMEGA^2 * coords
+
+# 1 period in the exact solution as computed below (but integrated with a small timestep)
+period = 4.567375
+tspan = (0.0, 1period)
+
+fluid_density = 1000.0
+sound_speed = 10.0
+state_equation = StateEquationCole(; sound_speed, exponent=7,
+                                   reference_density=fluid_density)
+
+# Equation A.19 in the paper rearranged.
+# sigma^2 = Q / rho * (a^2 + b^2) / (a^2 b^2) - OMEGA^2.
+Q = (sigma^2 + OMEGA^2) * fluid_density / 2
+pressure = coords -> Q * (1 - coords[1]^2 - coords[2]^2)
+density = coords -> TrixiParticles.inverse_state_equation(state_equation, pressure(coords))
+
+fluid = SphereShape(fluid_particle_spacing, radius, (0.0, 0.0),
+                    density, pressure=pressure,
+                    sphere_type=RoundSphere(),
+                    velocity=coords -> sigma .* (coords[1], -coords[2]))
+
+# ==========================================================================================
+# ==== Fluid
+smoothing_length = 3.0 * fluid_particle_spacing
+smoothing_kernel = WendlandC2Kernel{2}()
+
+fluid_density_calculator = ContinuityDensity()
+viscosity = ArtificialViscosityMonaghan(alpha=0.01, beta=0.0)
+
+density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1)
+fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           smoothing_length, viscosity=viscosity,
+                                           density_diffusion=density_diffusion,
+                                           source_terms=source_terms)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(fluid_system)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+saving_callback = SolutionSavingCallback(dt=0.04, prefix="")
+
+callbacks = CallbackSet(info_callback, saving_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control.
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-7, # Default abstol is 1e-6 (may need to be tuned to prevent intabilities)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent intabilities)
+            save_everystep=false, callback=callbacks);
+
+@inline function exact_solution_rhs(u, p, t)
+    sigma, A, B = u
+
+    dsigma = (sigma^2 + OMEGA^2) * ((B^2 - A^2) / (B^2 + A^2))
+    dA = sigma * A
+    dB = -sigma * B
+
+    return SVector(dsigma, dA, dB)
+end
+
+exact_u0 = SVector(sigma, radius, radius)
+exact_solution_ode = ODEProblem(exact_solution_rhs, exact_u0, tspan)
+
+# Use the same time integrator to avoid compilation of another integrator in CI
+sol_exact = solve(exact_solution_ode, RDPK3SpFSAL49(), save_everystep=false)
+
+# Error in the semi-major axis of the elliptical drop
+error_A = maximum(sol.u[end].x[2]) + 0.5fluid_particle_spacing -
+          maximum(sol_exact.u[end][2:3])

examples/fsi/dam_break_2d.jl

@@ -88,8 +88,8 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Solid
-solid_smoothing_length = sqrt(2) * solid_particle_spacing
-solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing
+solid_smoothing_kernel = WendlandC2Kernel{2}()
 
 # For the FSI we need the hydrodynamic masses and densities in the solid boundary model
 hydrodynamic_densites = fluid_density * ones(size(solid.density))

examples/fsi/dam_break_gate_2d.jl

@@ -117,8 +117,8 @@ boundary_system_gate = BoundarySPHSystem(gate, boundary_model_gate, movement=gat
 
 # ==========================================================================================
 # ==== Solid
-solid_smoothing_length = sqrt(2) * solid_particle_spacing
-solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing
+solid_smoothing_kernel = WendlandC2Kernel{2}()
 
 # For the FSI we need the hydrodynamic masses and densities in the solid boundary model
 hydrodynamic_densites = fluid_density * ones(size(solid.density))

examples/solid/oscillating_beam_2d.jl

@@ -45,12 +45,12 @@ solid = union(beam, fixed_particles)
 
 # ==========================================================================================
 # ==== Solid
+# The kernel in the reference uses a differently scaled smoothing length,
+# so this is equivalent to the smoothing length of `sqrt(2) * particle_spacing` used in the paper.
+smoothing_length = 2 * sqrt(2) * particle_spacing
+smoothing_kernel = WendlandC2Kernel{2}()
 
-smoothing_length = sqrt(2) * particle_spacing
-smoothing_kernel = SchoenbergCubicSplineKernel{2}()
-
-solid_system = TotalLagrangianSPHSystem(solid,
-                                        smoothing_kernel, smoothing_length,
+solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length,
                                         E, nu,
                                         n_fixed_particles=nparticles(fixed_particles),
                                         acceleration=(0.0, -gravity))
@@ -61,7 +61,28 @@ semi = Semidiscretization(solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
-saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+# 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]
+
+function x_deflection(v, u, t, system)
+    particle_position = TrixiParticles.current_coords(u, system, particle_id)
+
+    return particle_position[1] - shift_x
+end
+
+function y_deflection(v, u, t, system)
+    particle_position = TrixiParticles.current_coords(u, system, particle_id)
+
+    return particle_position[2] - shift_y
+end
+
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="",
+                                         x_deflection=x_deflection,
+                                         y_deflection=y_deflection)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

src/general/initial_condition.jl

@@ -100,9 +100,6 @@ struct InitialCondition{ELTYPE}
 
         if velocity isa AbstractMatrix
             velocities = velocity
-            if size(coordinates) != size(velocities)
-                throw(ArgumentError("`coordinates` and `velocities` must be of the same size"))
-            end
         else
             # Assuming `velocity` is a scalar or a function
             velocity_fun = wrap_function(velocity, Val(NDIMS))
@@ -111,7 +108,10 @@ struct InitialCondition{ELTYPE}
                                     "for $NDIMS-dimensional `coordinates`"))
             end
             velocities_svector = velocity_fun.(coordinates_svector)
-            velocities = reinterpret(reshape, ELTYPE, velocities_svector)
+            velocities = stack(velocities_svector)
+        end
+        if size(coordinates) != size(velocities)
+            throw(ArgumentError("`coordinates` and `velocities` must be of the same size"))
         end
 
         if density isa AbstractVector

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -33,23 +33,23 @@ The zero subscript on quantities denotes that the quantity is to be measured in
 The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$,
 the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$.
 
-For the computation of the PK1 stress tensor, the deformation gradient $\bm{J}$ is computed per particle as
+For the computation of the PK1 stress tensor, the deformation gradient $\bm{F}$ is computed per particle as
 ```math
-\bm{J}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\
+\bm{F}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\
     \qquad  = -\left(\sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ab} (\nabla_{0a} W(\bm{X}_{ab}))^T \right) \bm{L}_{0a}^T
 ```
 with $1 \leq i,j \leq d$.
 From the deformation gradient, the Green-Lagrange strain
 ```math
-\bm{E} = \frac{1}{2}(\bm{J}^T\bm{J} - \bm{I})
+\bm{E} = \frac{1}{2}(\bm{F}^T\bm{F} - \bm{I})
 ```
 and the second Piola-Kirchhoff stress tensor
 ```math
 \bm{S} = \lambda \operatorname{tr}(\bm{E}) \bm{I} + 2\mu \bm{E}
 ```
 are computed to obtain the PK1 stress tensor as
 ```math
-\bm{P} = \bm{J}\bm{S}.
+\bm{P} = \bm{F}\bm{S}.
 ```
 
 Here,
@@ -285,8 +285,8 @@ end
     calc_deformation_grad!(deformation_grad, neighborhood_search, system)
 
     @threaded for particle in eachparticle(system)
-        J_particle = deformation_gradient(system, particle)
-        pk1_particle = pk1_stress_tensor(J_particle, system)
+        F_particle = deformation_gradient(system, particle)
+        pk1_particle = pk1_stress_tensor(F_particle, system)
         pk1_particle_corrected = pk1_particle * correction_matrix(system, particle)
 
         @inbounds for j in 1:ndims(system), i in 1:ndims(system)
@@ -332,18 +332,18 @@ end
 end
 
 # First Piola-Kirchhoff stress tensor
-@inline function pk1_stress_tensor(J, system)
-    S = pk2_stress_tensor(J, system)
+@inline function pk1_stress_tensor(F, system)
+    S = pk2_stress_tensor(F, system)
 
-    return J * S
+    return F * S
 end
 
 # Second Piola-Kirchhoff stress tensor
-@inline function pk2_stress_tensor(J, system)
+@inline function pk2_stress_tensor(F, system)
     (; lame_lambda, lame_mu) = system
 
     # Compute the Green-Lagrange strain
-    E = 0.5 * (transpose(J) * J - I)
+    E = 0.5 * (transpose(F) * F - I)
 
     return lame_lambda * tr(E) * I + 2 * lame_mu * E
 end
@@ -411,3 +411,47 @@ end
 function viscosity_model(system::TotalLagrangianSPHSystem)
     return system.boundary_model.viscosity
 end
+
+# An explanation of these equation can be found in
+# J. Lubliner, 2008. Plasticity theory.
+# See here below Equation 5.3.21 for the equation for the equivalent stress.
+# 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::TotalLagrangianSPHSystem)
+    von_mises_stress = zeros(eltype(system.pk1_corrected), nparticles(system))
+
+    @threaded for particle in each_moving_particle(system)
+        F = deformation_gradient(system, particle)
+        J = det(F)
+        P = pk1_corrected(system, particle)
+        sigma = (1.0 / J) * P * F'
+
+        # Calculate deviatoric stress tensor
+        s = sigma - (1.0 / 3.0) * tr(sigma) * I
+
+        von_mises_stress[particle] = sqrt(3.0 / 2.0 * sum(s .^ 2))
+    end
+
+    return von_mises_stress
+end
+
+# An explanation of these equation can be found in
+# J. Lubliner, 2008. Plasticity theory.
+# See here page 473 for the relation between the `pk1`, the first Piola-Kirchhoff tensor,
+# and the Cauchy stress.
+function cauchy_stress(system::TotalLagrangianSPHSystem)
+    NDIMS = ndims(system)
+
+    cauchy_stress_tensors = zeros(eltype(system.pk1_corrected), NDIMS, NDIMS,
+                                  nparticles(system))
+
+    @threaded for particle in each_moving_particle(system)
+        F = deformation_gradient(system, particle)
+        J = det(F)
+        P = pk1_corrected(system, particle)
+        sigma = (1.0 / J) * P * F'
+        cauchy_stress_tensors[:, :, particle] = sigma
+    end
+
+    return cauchy_stress_tensors
+end

src/visualization/write2vtk.jl

@@ -216,11 +216,25 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d
 
     vtk["velocity"] = hcat(view(v, 1:ndims(system), :),
                            zeros(ndims(system), n_fixed_particles))
+    vtk["jacobian"] = [det(deformation_gradient(system, particle))
+                       for particle in eachparticle(system)]
+
+    vtk["von_mises_stress"] = von_mises_stress(system)
+
+    sigma = cauchy_stress(system)
+    vtk["sigma_11"] = sigma[1, 1, :]
+    vtk["sigma_22"] = sigma[2, 2, :]
+    if ndims(system) == 3
+        vtk["sigma_33"] = sigma[3, 3, :]
+    end
+
     vtk["material_density"] = system.material_density
     vtk["young_modulus"] = system.young_modulus
     vtk["poisson_ratio"] = system.poisson_ratio
     vtk["lame_lambda"] = system.lame_lambda
     vtk["lame_mu"] = system.lame_mu
+    vtk["smoothing_kernel"] = type2string(system.smoothing_kernel)
+    vtk["smoothing_length"] = system.smoothing_length
 
     write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)
 end

test/examples/examples.jl

@@ -2,6 +2,15 @@
 # but without checking the correctness of the solution.
 @testset verbose=true "Examples" begin
     @testset verbose=true "Fluid" begin
+        @trixi_testset "fluid/oscillating_drop_2d.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid",
+                                                "oscillating_drop_2d.jl"))
+            @test sol.retcode == ReturnCode.Success
+            # This is the error on an Apple M2 Pro. We need this tolerance to make CI pass.
+            @test isapprox(error_A, 0.0001717690010767381, atol=1e-8)
+        end
+
         @trixi_testset "fluid/rectangular_tank_2d.jl" begin
             @test_nowarn trixi_include(@__MODULE__,
                                        joinpath(examples_dir(), "fluid",