Merge branch 'dev' of github.com:trixi-framework/TrixiParticles.jl into dev

Author: efaulhaber

NEWS.md

@@ -28,6 +28,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
   It is now possible to pass multiple `BoundaryZone`s to a single `OpenBoundarySystem`.
   Reference values are now assigned individually to each `BoundaryZone`. (#866)
 
+- Rename keyword arguments `plane` and `plane_normal` for `BoundaryZone` to `boundary_face` and `face_normal` (#597).
+
 - The argument of `TransportVelocityAdami` is now a keyword argument.
   `TransportVelocityAdami(1000.0)` now becomes
   `TransportVelocityAdami(background_pressure=1000.0)` (#884).

Project.toml

@@ -27,6 +27,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StrideArraysCore = "7792a7ef-975c-4747-a70f-980b88e8d1da"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"

docs/src/systems/boundary.md

@@ -233,6 +233,11 @@ Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "boundary", "open_boundary", "boundary_zones.jl")]
 ```
 
+```@autodocs
+Modules = [TrixiParticles]
+Filter = t -> typeof(t) === typeof(TrixiParticles.planar_geometry_to_face)
+```
+
 # [Open Boundary Models](@id open_boundary_models)
 We offer two models for open boundaries, with the choice depending on the specific problem and flow characteristics near the boundary:
 1. [**Method of characteristics**](@ref method_of_characteristics): The method of characteristics is typically used in problems where tracking of wave propagation

examples/fluid/pipe_flow_2d.jl

@@ -113,9 +113,9 @@ reference_velocity_in = velocity_function2d
 reference_pressure_in = nothing
 reference_density_in = nothing
 boundary_type_in = InFlow()
-plane_in = ([0.0, 0.0], [0.0, domain_size[2]])
-inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers,
-                      density=fluid_density, particle_spacing,
+face_in = ([0.0, 0.0], [0.0, domain_size[2]])
+inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction,
+                      open_boundary_layers, density=fluid_density, particle_spacing,
                       reference_density=reference_density_in,
                       reference_pressure=reference_pressure_in,
                       reference_velocity=reference_velocity_in,
@@ -125,8 +125,8 @@ reference_velocity_out = nothing
 reference_pressure_out = nothing
 reference_density_out = nothing
 boundary_type_out = OutFlow()
-plane_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]])
-outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction),
+face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]])
+outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction),
                        open_boundary_layers, density=fluid_density, particle_spacing,
                        reference_density=reference_density_out,
                        reference_pressure=reference_pressure_out,

examples/fluid/pipe_flow_3d.jl

@@ -40,7 +40,7 @@ trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"),
               tspan=tspan, prescribed_velocity=prescribed_velocity,
               open_boundary_layers=open_boundary_layers, min_coords_inlet=min_coords_inlet,
               min_coords_outlet=min_coords_outlet,
-              plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0],
-                        [0.0, 0.0, domain_size[3]]),
-              plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0],
-                         [domain_size[1], 0.0, domain_size[3]]))
+              face_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0],
+                       [0.0, 0.0, domain_size[3]]),
+              face_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0],
+                        [domain_size[1], 0.0, domain_size[3]]))

src/TrixiParticles.jl

@@ -27,6 +27,7 @@ using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!
                  get_du
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
+using Statistics: Statistics
 using StrideArraysCore: PtrArray, StaticInt
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!, @notimeit
 using TrixiBase: @trixi_timeit, timer, timeit_debug_enabled,
@@ -92,7 +93,7 @@ export RectangularTank, RectangularShape, SphereShape, ComplexShape
 export ParticlePackingSystem, SignedDistanceField
 export WindingNumberHormann, WindingNumberJacobson
 export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry,
-       sample_boundary
+       sample_boundary, planar_geometry_to_face
 export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection

src/callbacks/post_process.jl

@@ -252,8 +252,10 @@ function (pp::PostprocessCallback)(integrator)
             system_index = system_indices(system, semi)
 
             for (key, f) in pp.func
-                result = custom_quantity(f, system, dv_ode, du_ode, v_ode, u_ode, semi, t)
-                if result !== nothing
+                result_ = custom_quantity(f, system, dv_ode, du_ode, v_ode, u_ode, semi, t)
+                if result_ !== nothing
+                    # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
+                    result = transfer2cpu(result_)
                     add_entry!(pp, string(key), t, result, filenames[system_index])
                     new_data = true
                 end

src/general/interpolation.jl

@@ -628,6 +628,11 @@ end
     return (; velocity, jacobian, von_mises_stress, cauchy_stress)
 end
 
+function interpolate_system!(cache, v, neighbor_system,
+                             point, neighbor, volume_b, W_ab, clip_negative_pressure)
+    return cache
+end
+
 @inline function interpolate_system!(cache, v, system::AbstractFluidSystem,
                                      point, neighbor, volume_b, W_ab,
                                      clip_negative_pressure)

src/io/io.jl

@@ -135,12 +135,8 @@ function add_system_data!(system_data, system::OpenBoundarySystem)
     system_data["system_type"] = type2string(system)
     system_data["fluid_system_index"] = system.fluid_system_index[]
     system_data["smoothing_length"] = system.smoothing_length
+    system_data["number_of_boundary_zones"] = length(system.boundary_zones)
     add_system_data!(system_data, system.boundary_model)
-
-    system_data["boundary_zones"] = Dict{String, Any}()
-    for (indice, boundary_zone) in enumerate(system.boundary_zones)
-        add_system_data!(system_data["boundary_zones"], boundary_zone, indice)
-    end
 end
 
 function add_system_data!(system_data, system::ParticlePackingSystem)
@@ -287,7 +283,7 @@ function add_system_data!(system_data, boundary_zone::BoundaryZone, indice)
     system_data[zone_name]["zone_origin"] = boundary_zone.zone_origin
     system_data[zone_name]["zone_width"] = boundary_zone.zone_width
     system_data[zone_name]["flow_direction"] = boundary_zone.flow_direction
-    system_data[zone_name]["plane_normal"] = boundary_zone.plane_normal
+    system_data[zone_name]["face_normal"] = boundary_zone.face_normal
     system_data[zone_name]["reference_values"] = boundary_zone.reference_values
     system_data[zone_name]["average_inflow_velocity"] = boundary_zone.average_inflow_velocity
     system_data[zone_name]["prescribed_density"] = boundary_zone.prescribed_density
@@ -299,8 +295,6 @@ function add_system_data!(system_data, motion::PrescribedMotion)
     system_data["prescribed_motion"] = Dict{String, Any}()
     system_data["prescribed_motion"]["model"] = type2string(motion)
     system_data["prescribed_motion"]["movement_function"] = type2string(motion.movement_function)
-    system_data["prescribed_motion"]["is_moving"] = type2string(motion.is_moving)
-    system_data["prescribed_motion"]["moving_particles"] = motion.moving_particles
 end
 
 function add_system_data!(system_data, penalty_force::PenaltyForceGanzenmueller)

src/io/write_vtk.jl

@@ -140,8 +140,8 @@ function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
 
         # Extract custom quantities for this system
         if !isempty(custom_quantities)
-            dv_ode, du_ode = dvdu_ode_.x
-            dv_ode, du_ode = transfer2cpu(dv_ode, du_ode)
+            dv_ode_, du_ode_ = dvdu_ode_.x
+            dv_ode, du_ode = transfer2cpu(dv_ode_, du_ode_)
 
             for (key, quantity) in custom_quantities
                 value = custom_quantity(quantity, system, dv_ode, du_ode, v_ode, u_ode,
@@ -173,8 +173,8 @@ function transfer2cpu(v_, u_, system_, semi_)
 end
 
 function transfer2cpu(v_::AbstractGPUArray, u_)
-    v = Adapt.adapt(Array, v_)
-    u = Adapt.adapt(Array, u_)
+    v = transfer2cpu(v_)
+    u = transfer2cpu(u_)
 
     return v, u
 end
@@ -183,6 +183,14 @@ function transfer2cpu(v_, u_)
     return v_, u_
 end
 
+function transfer2cpu(a_::AbstractGPUArray)
+    return Adapt.adapt(Array, a_)
+end
+
+function transfer2cpu(a_)
+    return a_
+end
+
 function custom_quantity(quantity::AbstractArray, system, dv_ode, du_ode, v_ode, u_ode,
                          semi, t)
     return quantity