Merge dev to main

NEWS.md

@@ -4,6 +4,14 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.5
+
+### API Changes
+
+- `DensityDiffusionAntuono` now only takes only one kwarg `delta` (#1142).
+- Return type of `vtk2trixi` changed to `NamedTuple` including an optional
+  `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
+
 ## Version 0.4.4
 
 ### API Changes
@@ -54,7 +62,6 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Added GPU and FP32 support for DEM (#979).
 
-
 ### Performance
 - Improved GPU performance with shifting up to a factor of 10x (#974, #993).
 

examples/fluid/dam_break_2d.jl

@@ -63,7 +63,7 @@ viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
 # free surface in long-running simulations, but is significantly faster than the model
 # by Antuono. This simulation is short enough to use the faster model.
 density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
-# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
+# density_diffusion = DensityDiffusionAntuono(delta=0.1)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,

examples/fluid/falling_water_spheres_2d.jl

@@ -57,7 +57,7 @@ fluid_density_calculator = ContinuityDensity()
 nu = 0.005
 alpha = 8 * nu / (fluid_smoothing_length * sound_speed)
 viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
-density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1)
+density_diffusion = DensityDiffusionAntuono(delta=0.1)
 
 sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
                                                      fluid_smoothing_length,

examples/fluid/oscillating_drop_2d.jl

@@ -55,7 +55,7 @@ smoothing_kernel = WendlandC2Kernel{2}()
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 
-density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1)
+density_diffusion = DensityDiffusionAntuono(delta=0.1)
 fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
                                            smoothing_length, viscosity=viscosity,

examples/fluid/periodic_array_of_cylinders_2d.jl

@@ -63,7 +63,7 @@ smoothing_kernel = WendlandC2Kernel{2}()
 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                    exponent=1, clip_negative_pressure=false)
 
-density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1)
+density_diffusion = DensityDiffusionAntuono(delta=0.1)
 fluid_system = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
                                            density_diffusion=density_diffusion,

examples/fsi/hydrostatic_water_column_2d.jl

@@ -65,7 +65,7 @@ right_wall = RectangularShape(structure_particle_spacing, (3, n_particles_plate_
                               place_on_shell=true)
 fixed_particles = union(left_wall, right_wall)
 
-structure_geometry = union(plate, fixed_particles)
+structure_geometry = union(fixed_particles, plate)
 
 # ==========================================================================================
 # ==== Smoothing Kernel, Boundary, and Related Quantities
@@ -103,7 +103,7 @@ if use_edac
 else
     fluid_density_calculator = ContinuityDensity()
     density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
-    # density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
+    # density_diffusion = DensityDiffusionAntuono(delta=0.1)
     fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                                state_equation, smoothing_kernel,
                                                smoothing_length_fluid,
@@ -127,7 +127,7 @@ boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densities,
 structure_system = TotalLagrangianSPHSystem(structure_geometry, smoothing_kernel,
                                             smoothing_length_structure,
                                             E, nu, boundary_model=boundary_model_structure,
-                                            n_clamped_particles=nparticles(fixed_particles),
+                                            clamped_particles=1:nparticles(fixed_particles),
                                             acceleration=(0.0, -gravity))
 
 min_corner = min.(minimum(structure_geometry.coordinates, dims=2),

src/io/read_vtk.jl

@@ -1,7 +1,11 @@
 """
-	vtk2trixi(file::String; element_type=nothing, coordinates_eltype=nothing)
+	vtk2trixi(file::String; element_type=nothing, coordinates_eltype=nothing,
+              create_initial_condition=true)
 
-Load VTK file and convert data to an [`InitialCondition`](@ref).
+Read a VTK file and return a `NamedTuple` with keys
+`:coordinates`, `:velocity`, `:density`, `:pressure`, `:particle_spacing`, `:time`, `:initial_condition`,
+plus any custom quantities.
+Missing fields are zero-filled; `:particle_spacing` is scalar if constant, otherwise per-particle.
 
 # Arguments
 - `file`: Name of the VTK file to be loaded.
@@ -13,35 +17,32 @@ Load VTK file and convert data to an [`InitialCondition`](@ref).
 - `coordinates_eltype`: Element type for particle coordinates. By default, the type
                         stored in the VTK file is used.
                         Otherwise, data is converted to the specified type.
+- `create_initial_condition`: If `true`, an `InitialCondition` object is created
+                              and included in the returned `NamedTuple` under
+                              the key `:initial_condition`. Default is `true`.
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 
 # Example
-```jldoctest; output = false
+```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]"
 # Create a rectangular shape
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
 
-# Write the `InitialCondition` to a vtk file
-trixi2vtk(rectangular; filename="rectangular", output_directory="out")
+# Write the `InitialCondition` with custom quantity to a VTK file
+trixi2vtk(rectangular; filename="rectangular", output_directory="out",
+          my_custom_quantity=3.0)
 
-# Read the vtk file and convert it to `InitialCondition`
-ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
+# Read the VTK file and convert the data to a `NamedTuple`
+data = vtk2trixi(joinpath("out", "rectangular.vtu"))
 
 # output
-┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
-│ InitialCondition                                                                                 │
-│ ════════════════                                                                                 │
-│ #dimensions: ……………………………………………… 2                                                                │
-│ #particles: ………………………………………………… 100                                                              │
-│ particle spacing: ………………………………… 0.1                                                              │
-│ eltype: …………………………………………………………… Float64                                                          │
-│ coordinate eltype: ……………………………… Float64                                                          │
-└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], mass = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...], initial_condition = InitialCondition{Float64, Float64}())
 ```
 """
-function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing)
+function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
+                   create_initial_condition=true)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -53,38 +54,67 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing)
     ELTYPE = isnothing(element_type) ?
              eltype(first(ReadVTK.get_data(point_data["pressure"]))) : element_type
 
+    results = Dict{Symbol, Any}()
+
+    # Tracking used keys like `initial_velocity`
+    used_keys = String[]
+
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
     coordinates = convert.(cELTYPE, point_coords[1:ndims, :])
 
-    fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
-    results = Dict{String, Array{Float64}}()
-
+    fields = [:velocity, :density, :pressure, :particle_spacing]
     for field in fields
         # First look for an exact key match, then fall back to substring matching.
         all_keys = keys(point_data)
         idx = findfirst(k -> k == field, all_keys)
         if idx === nothing
-            idx = findfirst(k -> occursin(field, k), all_keys)
+            idx = findfirst(k -> occursin(string(field), k), all_keys)
         end
         if idx !== nothing
             results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
+            push!(used_keys, all_keys[idx])
         else
             # Use zeros as default values when a field is missing
-            results[field] = field in ["mass"] ?
-                             zeros(ELTYPE, size(coordinates, 2)) : zero(coordinates)
+            results[field] = string(field) in ["velocity"] ?
+                             zero(coordinates) : zeros(ELTYPE, size(coordinates, 2))
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
 
-    particle_spacing = allequal(results["particle_spacing"]) ?
-                       first(results["particle_spacing"]) :
-                       results["particle_spacing"]
+    results[:particle_spacing] = allequal(results[:particle_spacing]) ?
+                                 first(results[:particle_spacing]) :
+                                 results[:particle_spacing]
+    results[:coordinates] = coordinates
+    results[:time] = "time" in keys(field_data) ?
+                     first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
+
+    append!(used_keys, ["index", "ndims"])
+    # Load any custom quantities
+    for key in keys(point_data)
+        if !(key in used_keys)
+            results[Symbol(key)] = convert.(ELTYPE, ReadVTK.get_data(point_data[key]))
+        end
+    end
+
+    for key in keys(field_data)
+        if !(key in used_keys)
+            data = ReadVTK.get_data(field_data[key])
+            conv_data = convert.(ELTYPE, data)
+            results[Symbol(key)] = length(conv_data) == 1 ? first(conv_data) : conv_data
+        end
+    end
+
+    results = NamedTuple(results)
+
+    if create_initial_condition
+        ic = InitialCondition(; coordinates=results.coordinates,
+                              particle_spacing=results.particle_spacing,
+                              velocity=results.velocity, density=results.density,
+                              pressure=results.pressure)
 
-    return InitialCondition(; coordinates,
-                            particle_spacing=convert(ELTYPE, particle_spacing),
-                            velocity=results["velocity"],
-                            mass=results["mass"],
-                            density=results["density"],
-                            pressure=results["pressure"])
+        return (; results..., initial_condition=ic)
+    else
+        return results
+    end
 end

src/schemes/boundary/open_boundary/system.jl

@@ -178,14 +178,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
         # as it was already verified in `allocate_buffer` that the density array is constant.
         density_rest = first(initial_condition.density)
 
-        if density_diffusion isa DensityDiffusionAntuono
-            density_diffusion_ = DensityDiffusionAntuono(initial_condition;
-                                                         delta=density_diffusion.delta)
-        else
-            density_diffusion_ = density_diffusion
-        end
-
-        cache = (; density_diffusion=density_diffusion_,
+        cache = (; density_diffusion=density_diffusion,
+                 create_cache_density_diffusion(initial_condition, density_diffusion)...,
                  pressure_boundary=pressure_boundary,
                  density_rest=density_rest, cache...)
 

src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl

@@ -88,7 +88,7 @@ end
 end
 
 @doc raw"""
-    DensityDiffusionAntuono(initial_condition; delta)
+    DensityDiffusionAntuono(; delta)
 
 The commonly used density diffusion terms by [Antuono (2010)](@cite Antuono2010), also referred to as
 δ-SPH. The density diffusion term by [Molteni (2009)](@cite Molteni2009) is extended by a second
@@ -115,33 +115,31 @@ where ``d`` is the number of dimensions.
 See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion)
 for an overview and comparison of implemented density diffusion terms.
 """
-struct DensityDiffusionAntuono{NDIMS, ELTYPE, ARRAY2D, ARRAY3D} <: AbstractDensityDiffusion
-    delta                       :: ELTYPE
-    correction_matrix           :: ARRAY3D # Array{ELTYPE, 3}: [i, j, particle]
-    normalized_density_gradient :: ARRAY2D # Array{ELTYPE, 2}: [i, particle]
-
-    function DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient)
-        new{size(correction_matrix, 1), typeof(delta),
-            typeof(normalized_density_gradient),
-            typeof(correction_matrix)}(delta, correction_matrix,
-                                       normalized_density_gradient)
+struct DensityDiffusionAntuono{ELTYPE} <: AbstractDensityDiffusion
+    delta::ELTYPE
+
+    function DensityDiffusionAntuono(; delta)
+        new{typeof(delta)}(delta)
     end
 end
 
-function DensityDiffusionAntuono(initial_condition; delta)
+create_cache_density_diffusion(initial_condition, density_diffusion) = (;)
+
+function create_cache_density_diffusion(initial_condition,
+                                        ::DensityDiffusionAntuono)
     NDIMS = ndims(initial_condition)
     ELTYPE = eltype(initial_condition)
-    correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS,
-                                         nparticles(initial_condition))
+    n_particles = nparticles(initial_condition)
 
-    normalized_density_gradient = Array{ELTYPE, 2}(undef, NDIMS,
-                                                   nparticles(initial_condition))
+    density_diffusion_correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS,
+                                                           n_particles)
+    density_diffusion_normalized_density_gradient = Array{ELTYPE, 2}(undef, NDIMS,
+                                                                     n_particles)
 
-    return DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient)
+    return (; density_diffusion_correction_matrix,
+            density_diffusion_normalized_density_gradient)
 end
 
-@inline Base.ndims(::DensityDiffusionAntuono{NDIMS}) where {NDIMS} = NDIMS
-
 function Base.show(io::IO, density_diffusion::DensityDiffusionAntuono)
     @nospecialize density_diffusion # reduce precompilation time
 
@@ -154,22 +152,19 @@ function allocate_buffer(initial_condition, density_diffusion, buffer)
     return allocate_buffer(initial_condition, buffer), density_diffusion
 end
 
-function allocate_buffer(ic, dd::DensityDiffusionAntuono, buffer::SystemBuffer)
-    initial_condition = allocate_buffer(ic, buffer)
-    return initial_condition, DensityDiffusionAntuono(initial_condition; delta=dd.delta)
-end
-
 @propagate_inbounds function density_diffusion_psi(density_diffusion::DensityDiffusionAntuono,
                                                    rho_a, rho_b, pos_diff, distance, system,
                                                    particle, neighbor)
-    (; normalized_density_gradient) = density_diffusion
+    (; density_diffusion_normalized_density_gradient) = system.cache
 
     # First term by Molteni & Colagrossi
     result = 2 * (rho_a - rho_b)
 
     # Second correction term
-    normalized_gradient_a = extract_svector(normalized_density_gradient, system, particle)
-    normalized_gradient_b = extract_svector(normalized_density_gradient, system, neighbor)
+    normalized_gradient_a = extract_svector(density_diffusion_normalized_density_gradient,
+                                            system, particle)
+    normalized_gradient_b = extract_svector(density_diffusion_normalized_density_gradient,
+                                            system, neighbor)
     result -= dot(normalized_gradient_a + normalized_gradient_b, pos_diff)
 
     # Since this is one of the most performance critical functions, using fast divisions
@@ -179,13 +174,14 @@ end
 end
 
 function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi)
-    (; normalized_density_gradient) = density_diffusion
+    density_diffusion_correction_matrix = system.cache.density_diffusion_correction_matrix
+    normalized_density_gradient = system.cache.density_diffusion_normalized_density_gradient
 
     # Compute correction matrix
     density_fun = @propagate_inbounds(particle->current_density(v, system, particle))
     system_coords = current_coordinates(u, system)
 
-    compute_gradient_correction_matrix!(density_diffusion.correction_matrix,
+    compute_gradient_correction_matrix!(density_diffusion_correction_matrix,
                                         system, system_coords, density_fun, semi)
 
     # Compute normalized density gradient
@@ -198,7 +194,8 @@ function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi)
         rho_b = @inbounds current_density(v, system, neighbor)
 
         grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
-        L = @inbounds correction_matrix(density_diffusion, particle)
+        L = @inbounds extract_smatrix(density_diffusion_correction_matrix, system,
+                                      particle)
 
         m_b = @inbounds hydrodynamic_mass(system, neighbor)
         volume_b = m_b / rho_b

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -152,6 +152,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, stat
                                           n_particles)...,
              create_cache_refinement(initial_condition, particle_refinement,
                                      smoothing_length)...,
+             create_cache_density_diffusion(initial_condition, density_diffusion)...,
              create_cache_shifting(initial_condition, shifting_technique)...,
              # Per-system color tag for colorfield surface-normal logic and VTK output.
              color=Int(color_value))