Multiphase

Project.toml

@@ -1,8 +1,38 @@
+authors = ["erik.faulhaber <[email protected]>"]
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
 version = "0.4.2"
 
+[compat]
+Adapt = "4"
+CSV = "0.10"
+DataFrames = "1.6"
+DelimitedFiles = "1"
+DiffEqCallbacks = "4"
+FastPow = "0.1"
+FileIO = "1"
+ForwardDiff = "1"
+GPUArraysCore = "0.2"
+JSON = "0.21"
+KernelAbstractions = "0.9"
+MuladdMacro = "0.2"
+OrdinaryDiffEq = "6.91"
+OrdinaryDiffEqCore = "1"
+PointNeighbors = "0.6.3"
+Polyester = "0.7.10"
+ReadVTK = "0.2"
+RecipesBase = "1"
+Reexport = "1"
+SciMLBase = "2"
+StaticArrays = "1"
+Statistics = "1"
+StrideArraysCore = "0.5.7"
+TimerOutputs = "0.5.25"
+TrixiBase = "0.1.5"
+WriteVTK = "1.21.2"
+julia = "1.10"
+
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -33,38 +63,12 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
-[weakdeps]
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-
 [extensions]
 TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"]
 
-[compat]
-Adapt = "4"
-CSV = "0.10"
-DataFrames = "1.6"
-DelimitedFiles = "1"
-DiffEqCallbacks = "4"
-FastPow = "0.1"
-FileIO = "1"
-ForwardDiff = "1"
-GPUArraysCore = "0.2"
-JSON = "1"
-KernelAbstractions = "0.9"
-MuladdMacro = "0.2"
-OrdinaryDiffEq = "6.91"
-OrdinaryDiffEqCore = "1"
-PointNeighbors = "0.6.3"
-Polyester = "0.7.10"
-ReadVTK = "0.2"
-RecipesBase = "1"
-Reexport = "1"
-SciMLBase = "2"
-StaticArrays = "1"
-Statistics = "1"
-StrideArraysCore = "0.5.7"
-TimerOutputs = "0.5.25"
-TrixiBase = "0.1.5"
-WriteVTK = "1.21.2"
-julia = "1.10"
+[extras]
+CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
+
+[weakdeps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"

examples/fluid/dam_break_2d_gpu.jl

@@ -36,9 +36,11 @@ trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               neighborhood_search=neighborhood_search,
               fluid_particle_spacing=fluid_particle_spacing,
-              tspan=tspan, smoothing_length=smoothing_length,
+              tspan=tspan,
+              smoothing_length=smoothing_length,
               density_diffusion=density_diffusion,
-              boundary_layers=boundary_layers, spacing_ratio=spacing_ratio,
+              boundary_layers=boundary_layers,
+              spacing_ratio=spacing_ratio,
               boundary_model=boundary_model,
               parallelization_backend=PolyesterBackend(),
               boundary_density_calculator=boundary_density_calculator)

examples/fluid/dam_break_ds.jl

@@ -0,0 +1,38 @@
+using TrixiParticles
+using CUDA
+
+fluid_particle_spacing = 0.001
+
+# Load setup from dam break example
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+              fluid_particle_spacing=fluid_particle_spacing,
+              tank_size=(4.0, 3.0), W=1.0, H=2.0,
+              spacing_ratio=1, boundary_layers=1,
+              sol=nothing, ode=nothing)
+
+tank.fluid.coordinates .+= 0.005
+tank.boundary.coordinates .+= 0.005
+
+# Define a GPU-compatible neighborhood search
+min_corner = minimum(tank.boundary.coordinates, dims=2) .- 1
+max_corner = maximum(tank.boundary.coordinates, dims=2) .+ 1
+cell_list = FullGridCellList(; min_corner, max_corner, max_points_per_cell=30)
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list,
+                                                update_strategy=ParallelUpdate())
+
+# Run the dam break simulation with this neighborhood search
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+              tank=tank,
+              smoothing_length=1.414216 * fluid_particle_spacing,
+              time_integration_scheme=SymplecticPositionVerlet(),
+              boundary_density_calculator=ContinuityDensity(),
+              fluid_particle_spacing=fluid_particle_spacing,
+              tank_size=(4.0, 3.0), W=1.0, H=2.0,
+              spacing_ratio=1, boundary_layers=1,
+              tspan=(0.0, 0.05), cfl=0.2,
+              neighborhood_search=neighborhood_search,
+              viscosity_wall=viscosity_fluid,
+              dt=1e-5, stepsize_callback=nothing, saving_callback=nothing,
+              parallelization_backend=CUDABackend())

examples/paper/dam_break_2d_val_600.jl

@@ -0,0 +1,42 @@
+# This file computes the pressure sensor data of the dam break setup described in
+#
+# J. J. De Courcy, T. C. S. Rendall, L. Constantin, B. Titurus, J. E. Cooper.
+# "Incompressible δ-SPH via artificial compressibility".
+# In: Computer Methods in Applied Mechanics and Engineering, Volume 420 (2024),
+# https://doi.org/10.1016/j.cma.2023.116700
+
+using TrixiParticles
+using TrixiParticles.JSON
+using CUDA
+
+# When using data center CPUs with large numbers of cores, especially on multi-socket
+# systems with multiple NUMA nodes, pinning threads to cores can significantly
+# improve performance, even for low resolutions.
+# using ThreadPinning
+# pinthreads(:numa)
+
+# `resolution` in this case is set relative to `H`, the initial height of the fluid.
+# Use 40, 80 or 400 for validation.
+# Note: 400 takes about 30 minutes on a large data center CPU (much longer with serial update)
+resolution = 600
+
+min_corner = (-1.5, -1.5)
+max_corner = (6.5, 5.0)
+cell_list = FullGridCellList(; min_corner, max_corner, max_points_per_cell=30)
+
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, update_strategy=ParallelUpdate())
+# neighborhood_search = GridNeighborhoodSearch{2}(; cell_list)
+
+
+# ==========================================================================================
+# ==== WCSPH simulation
+trixi_include_changeprecision(Float32, @__MODULE__,
+              joinpath(validation_dir(), "dam_break_2d",
+                       "setup_marrone_2011.jl"),
+              use_edac=false,
+              particles_per_height=resolution,
+              sound_speed=50 * sqrt(9.81 * 0.6), # This is used by De Courcy et al. (2024)
+              alpha=0.01, # This is used by De Courcy et al. (2024)
+              tspan=(0.0, 7 / sqrt(9.81 / 0.6)), # This is used by De Courcy et al. (2024)
+              parallelization_backend=CUDABackend(),
+              neighborhood_search=neighborhood_search)

examples/paper/dam_break_2d_val_600_phys_viscosity.jl

@@ -0,0 +1,49 @@
+# This file computes the pressure sensor data of the dam break setup described in
+#
+# J. J. De Courcy, T. C. S. Rendall, L. Constantin, B. Titurus, J. E. Cooper.
+# "Incompressible δ-SPH via artificial compressibility".
+# In: Computer Methods in Applied Mechanics and Engineering, Volume 420 (2024),
+# https://doi.org/10.1016/j.cma.2023.116700
+
+using TrixiParticles
+using TrixiParticles.JSON
+using CUDA
+
+# When using data center CPUs with large numbers of cores, especially on multi-socket
+# systems with multiple NUMA nodes, pinning threads to cores can significantly
+# improve performance, even for low resolutions.
+# using ThreadPinning
+# pinthreads(:numa)
+
+# `resolution` in this case is set relative to `H`, the initial height of the fluid.
+# Use 40, 80 or 400 for validation.
+# Note: 400 takes about 30 minutes on a large data center CPU (much longer with serial update)
+resolution = 600
+
+min_corner = (-1.5, -1.5)
+max_corner = (6.5, 5.0)
+cell_list = FullGridCellList(; min_corner, max_corner, max_points_per_cell=30)
+
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, update_strategy=ParallelUpdate())
+# neighborhood_search = GridNeighborhoodSearch{2}(; cell_list)
+
+# switch to physical viscosity model
+viscosity_fluid = ViscosityMorris(nu = 8.9E-7)
+# viscosity_fluid = ViscosityAdami(nu = 8.9E-7)
+#  viscosity_fluid = nothing
+
+
+# ==========================================================================================
+# ==== WCSPH simulation
+trixi_include(@__MODULE__,
+              joinpath(validation_dir(), "dam_break_2d",
+                       "setup_marrone_2011.jl"),
+              use_edac=false,
+              extra_string = "_phys_viscosity",
+              viscosity_fluid=viscosity_fluid,
+              particles_per_height=resolution,
+              sound_speed=50 * sqrt(9.81 * 0.6), # This is used by De Courcy et al. (2024)
+              tspan=(0.0, 7 / sqrt(9.81 / 0.6)), # This is used by De Courcy et al. (2024)
+            #   tspan=(0.0, 0.0001), # This is used by De Courcy et al. (2024)
+              parallelization_backend=CUDABackend(),
+              neighborhood_search=neighborhood_search)

examples/paper/dam_break_2d_val_600_phys_viscosity_gpu.jl

@@ -0,0 +1,16 @@
+# This file computes the pressure sensor data of the dam break setup described in
+#
+# J. J. De Courcy, T. C. S. Rendall, L. Constantin, B. Titurus, J. E. Cooper.
+# "Incompressible δ-SPH via artificial compressibility".
+# In: Computer Methods in Applied Mechanics and Engineering, Volume 420 (2024),
+# https://doi.org/10.1016/j.cma.2023.116700
+
+using TrixiParticles
+using TrixiParticles.JSON
+using CUDA
+
+# ==========================================================================================
+# ==== WCSPH simulation
+trixi_include_changeprecision(Float32, @__MODULE__,
+              joinpath(examples_dir(), "paper",
+                       "dam_break_2d_val_600_phys_viscosity.jl"))

src/callbacks/post_process.jl

@@ -230,20 +230,20 @@ end
 function (pp::PostprocessCallback)(integrator; from_initialize=false)
     @trixi_timeit timer() "apply postprocess cb" begin
         vu_ode = integrator.u
-        if from_initialize
+        # if from_initialize
             # Avoid calling `get_du` here, since it will call the RHS function
             # if it is called before the first time step.
             # This would cause problems with `semi.update_callback_used`,
             # which might not yet be set to `true` at this point if the `UpdateCallback`
             # comes AFTER the `PostprocessCallback` in the `CallbackSet`.
             dv_ode, du_ode = zero(vu_ode).x
-        else
-            # Depending on the time integration scheme, this might call the RHS function
-            @trixi_timeit timer() "update du" begin
-                # Don't create sub-timers here to avoid cluttering the timer output
-                @notimeit timer() dv_ode, du_ode = get_du(integrator).x
-            end
-        end
+        # else
+        #     # Depending on the time integration scheme, this might call the RHS function
+        #     @trixi_timeit timer() "update du" begin
+        #         # Don't create sub-timers here to avoid cluttering the timer output
+        #         @notimeit timer() dv_ode, du_ode = get_du(integrator).x
+        #     end
+        # end
         v_ode, u_ode = vu_ode.x
         semi = integrator.p
         t = integrator.t

src/callbacks/solution_saving.jl

@@ -155,16 +155,17 @@ function (solution_callback::SolutionSavingCallback)(integrator; from_initialize
      prefix, latest_saved_iter, max_coordinates) = solution_callback
 
     vu_ode = integrator.u
-    if from_initialize
+    # TODO doesn't work when `v_ode` and `u_ode` have different eltypes
+    # if from_initialize
         # Avoid calling `get_du` here, since it will call the RHS function
         # if it is called before the first time step.
         # This would cause problems with `semi.update_callback_used`,
         # which might not yet be set to `true` at this point if the `UpdateCallback`
         # comes AFTER the `SolutionSavingCallback` in the `CallbackSet`.
         dvdu_ode = zero(vu_ode)
-    else
-        dvdu_ode = get_du(integrator)
-    end
+    # else
+    #     dvdu_ode = get_du(integrator)
+    # end
     semi = integrator.p
     iter = get_iter(interval, integrator)
 

src/general/abstract_system.jl

@@ -97,6 +97,8 @@ end
 # This can be dispatched by system type
 @inline initial_coordinates(system) = system.initial_condition.coordinates
 
+@inline coordinates_eltype(system::AbstractSystem) = eltype(initial_coordinates(system))
+
 @propagate_inbounds function current_velocity(v, system, particle)
     return extract_svector(current_velocity(v, system), system, particle)
 end

src/general/initial_condition.jl

@@ -81,17 +81,17 @@ initial_condition = InitialCondition(; coordinates, velocity=x -> 2x, mass=1.0,
 
 # output
 ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
-│ InitialCondition{Float64}                                                                        │
-│ ═════════════════════════                                                                        │
+│ InitialCondition{Float64, Float64}                                                               │
+│ ═════════════════════════=========                                                               │
 │ #dimensions: ……………………………………………… 2                                                                │
 │ #particles: ………………………………………………… 3                                                                │
 │ particle spacing: ………………………………… -1.0                                                             │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-struct InitialCondition{ELTYPE, MATRIX, VECTOR}
+struct InitialCondition{ELTYPE, C, MATRIX, VECTOR}
     particle_spacing :: ELTYPE
-    coordinates      :: MATRIX # Array{ELTYPE, 2}
+    coordinates      :: C      # Array{coordinates_eltype, 2}
     velocity         :: MATRIX # Array{ELTYPE, 2}
     mass             :: VECTOR # Array{ELTYPE, 1}
     density          :: VECTOR # Array{ELTYPE, 1}
@@ -111,17 +111,28 @@ end
 # Function barrier to make `NDIMS` static and therefore SVectors type-stable
 function InitialCondition{NDIMS}(coordinates, velocity, mass, density,
                                  pressure, particle_spacing) where {NDIMS}
-    ELTYPE = eltype(coordinates)
+    if !(particle_spacing isa AbstractFloat)
+        throw(ArgumentError("particle spacing must be a floating point number. " *
+                            "The type of the particle spacing determines the eltype " *
+                            "of the `InitialCondition`."))
+    end
+
+    # The arguments can all be functions, so use `particle_spacing` for the eltype
+    ELTYPE = typeof(particle_spacing)
     n_particles = size(coordinates, 2)
 
+    # Coordinates can have a different eltype, e.g., Float64 if others are Float32
+    coordinates_eltype = eltype(coordinates)
+
     if n_particles == 0
         return InitialCondition(particle_spacing, coordinates, zeros(ELTYPE, NDIMS, 0),
                                 zeros(ELTYPE, 0), zeros(ELTYPE, 0), zeros(ELTYPE, 0))
     end
 
     # SVector of coordinates to pass to functions.
     # This will return a vector of SVectors in 2D and 3D, but an 1×n matrix in 1D.
-    coordinates_svector_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, coordinates)
+    coordinates_svector_ = reinterpret(reshape, SVector{NDIMS, coordinates_eltype},
+                                       coordinates)
     # In 1D, this will reshape the 1×n matrix to a vector, in 2D/3D it will do nothing
     coordinates_svector = reshape(coordinates_svector_, length(coordinates_svector_))
 
@@ -194,19 +205,21 @@ end
 function Base.show(io::IO, ic::InitialCondition)
     @nospecialize ic # reduce precompilation time
 
-    print(io, "InitialCondition{$(eltype(ic))}()")
+    print(io, "InitialCondition{$(eltype(ic)), $(coordinates_eltype(ic))}()")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", ic::InitialCondition)
     @nospecialize ic # reduce precompilation time
 
     if get(io, :compact, false)
-        show(io, system)
+        show(io, ic)
     else
-        summary_header(io, "InitialCondition{$(eltype(ic))}")
+        summary_header(io, "InitialCondition")
         summary_line(io, "#dimensions", "$(ndims(ic))")
         summary_line(io, "#particles", "$(nparticles(ic))")
         summary_line(io, "particle spacing", "$(first(ic.particle_spacing))")
+        summary_line(io, "eltype", "$(eltype(ic))")
+        summary_line(io, "coordinate eltype", "$(coordinates_eltype(ic))")
         summary_footer(io)
     end
 end
@@ -229,7 +242,9 @@ end
     return size(initial_condition.coordinates, 1)
 end
 
-@inline function Base.eltype(initial_condition::InitialCondition)
+@inline Base.eltype(::InitialCondition{ELTYPE}) where {ELTYPE} = ELTYPE
+
+@inline function coordinates_eltype(initial_condition::InitialCondition)
     return eltype(initial_condition.coordinates)
 end