Merge branch 'add-APR-struct' into additional-terms

src/TrixiParticles.jl

@@ -40,6 +40,7 @@ include("general/system.jl")
 # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded`
 include("util.jl")
 include("preprocessing/preprocessing.jl")
+include("multi_resolution/multi_resolution.jl")
 include("callbacks/callbacks.jl")
 include("general/general.jl")
 include("setups/setups.jl")

src/multi_resolution/multi_resolution.jl

@@ -0,0 +1 @@
+include("particle_refinement.jl")

src/multi_resolution/particle_refinement.jl

@@ -0,0 +1 @@
+struct ParticleRefinement end

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -63,55 +63,55 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
 end
 
 function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
-                                         smoothing_length, sound_speed;
-                                         pressure_acceleration=inter_particle_averaged_pressure,
-                                         density_calculator=SummationDensity(),
-                                         transport_velocity=nothing,
-                                         alpha=0.5, viscosity=nothing,
-                                         acceleration=ntuple(_ -> 0.0,
-                                                             ndims(smoothing_kernel)),
-                                         particle_refinement=nothing,
-                                         source_terms=nothing, buffer_size=nothing)
-        buffer = isnothing(buffer_size) ? nothing :
-                 SystemBuffer(nparticles(initial_condition), buffer_size)
-
-        initial_condition = allocate_buffer(initial_condition, buffer)
-
-        NDIMS = ndims(initial_condition)
-        ELTYPE = eltype(initial_condition)
-
-        mass = copy(initial_condition.mass)
-
-        if ndims(smoothing_kernel) != NDIMS
-            throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
-        end
-
-        acceleration_ = SVector(acceleration...)
-        if length(acceleration_) != NDIMS
-            throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
-        end
+                                     smoothing_length, sound_speed;
+                                     pressure_acceleration=inter_particle_averaged_pressure,
+                                     density_calculator=SummationDensity(),
+                                     transport_velocity=nothing,
+                                     alpha=0.5, viscosity=nothing,
+                                     acceleration=ntuple(_ -> 0.0,
+                                                         ndims(smoothing_kernel)),
+                                     particle_refinement=nothing,
+                                     source_terms=nothing, buffer_size=nothing)
+    buffer = isnothing(buffer_size) ? nothing :
+             SystemBuffer(nparticles(initial_condition), buffer_size)
+
+    initial_condition = allocate_buffer(initial_condition, buffer)
+
+    NDIMS = ndims(initial_condition)
+    ELTYPE = eltype(initial_condition)
+
+    mass = copy(initial_condition.mass)
+
+    if ndims(smoothing_kernel) != NDIMS
+        throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
+    end
 
-        pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
-                                                                         density_calculator,
-                                                                         NDIMS, ELTYPE,
-                                                                         nothing)
-
-        nu_edac = (alpha * smoothing_length * sound_speed) /
-                  (2 * ndims(initial_condition) + 4)
-
-        cache = create_cache_density(initial_condition, density_calculator)
-        cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
-        cache = (;
-             create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS,
-                                     nparticles(initial_condition))...,
-             cache...)
-
-        return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, smoothing_kernel,
-            sound_speed, viscosity, nu_edac, acceleration_,
-            nothing, pressure_acceleration, transport_velocity, source_terms,
-            buffer, particle_refinement, cache)
+    acceleration_ = SVector(acceleration...)
+    if length(acceleration_) != NDIMS
+        throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
+    pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
+                                                                     density_calculator,
+                                                                     NDIMS, ELTYPE,
+                                                                     nothing)
+
+    nu_edac = (alpha * smoothing_length * sound_speed) /
+              (2 * ndims(initial_condition) + 4)
+
+    cache = create_cache_density(initial_condition, density_calculator)
+    cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
+    cache = (;
+             create_cache_refinement(initial_condition, particle_refinement,
+                                     smoothing_length)..., cache...)
+
+    return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator,
+                                       smoothing_kernel, sound_speed, viscosity, nu_edac,
+                                       acceleration_, nothing, pressure_acceleration,
+                                       transport_velocity, source_terms, buffer,
+                                       particle_refinement, cache)
+end
+
 function Base.show(io::IO, system::EntropicallyDampedSPHSystem)
     @nospecialize system # reduce precompilation time
 
@@ -148,18 +148,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst
     end
 end
 
-function smoothing_length(system::EntropicallyDampedSPHSystem, particle)
-    return smoothing_length(system, system.particle_refinement, particle)
-end
-
-function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle)
-    return system.cache.smoothing_length
-end
-
-function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle)
-    return system.cache.smoothing_length[particle]
-end
-
 create_cache_edac(initial_condition, ::Nothing) = (;)
 
 function create_cache_edac(initial_condition, ::TransportVelocityAdami)

src/schemes/fluid/fluid.jl

@@ -13,8 +13,28 @@ function create_cache_density(ic, ::ContinuityDensity)
     return (;)
 end
 
+function create_cache_refinement(initial_condition, ::Nothing, smoothing_length)
+    return (; smoothing_length)
+end
+
+function create_cache_refinement(initial_condition, refinement, smoothing_length)
+    return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)))
+end
+
 @inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle]
 
+@inline function smoothing_length(system::FluidSystem, particle)
+    return smoothing_length(system, system.particle_refinement, particle)
+end
+
+@inline function smoothing_length(system::FluidSystem, ::Nothing, particle)
+    return system.cache.smoothing_length
+end
+
+@inline function smoothing_length(system::FluidSystem, refinement, particle)
+    return system.cache.smoothing_length[particle]
+end
+
 function write_u0!(u0, system::FluidSystem)
     (; initial_condition) = system
 

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -44,8 +44,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method.
 
 
 """
-struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K,
-                                   V, DD, COR, PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC}
+struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR,
+                                   PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC}
     initial_condition                 :: IC
     mass                              :: MA     # Array{ELTYPE, 1}
     pressure                          :: P      # Array{ELTYPE, 1}
@@ -118,8 +118,8 @@ function WeaklyCompressibleSPHSystem(initial_condition,
              create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)...,
              cache...)
     cache = (;
-             create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, n_particles)...,
-             cache...)
+             create_cache_refinement(initial_condition, particle_refinement,
+                                     smoothing_length)..., cache...)
 
     return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure,
                                        density_calculator, state_equation,
@@ -160,26 +160,6 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles)
     return (; surface_normal)
 end
 
-function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles)
-    return (; smoothing_length)
-end
-
-function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles)
-    return (; smoothing_length=smoothing_length * ones(n_particles))
-end
-
-function smoothing_length(system::WeaklyCompressibleSPHSystem, particle)
-    return smoothing_length(system, system.particle_refinement, particle)
-end
-
-function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle)
-    return system.cache.smoothing_length
-end
-
-function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle)
-    return system.cache.smoothing_length[particle]
-end
-
 function Base.show(io::IO, system::WeaklyCompressibleSPHSystem)
     @nospecialize system # reduce precompilation time