diff --git a/README.md b/README.md index e1bf6bb4..8b1411cd 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,7 @@ Example scripts live in `example/` and read geometry from `input/`. The solver c src/ ├── AuxiliaryFunctions.jl # Small helper utilities ├── OpenExternalPrograms.jl # Convenience wrappers for logs and ParaView +├── PointMeasurements.jl # Point probe configuration helpers ├── PreProcess.jl # Load inputs and allocate arrays ├── ProduceHDFVTK.jl # Write simulation data in HDF5/VTK format ├── SPHCellList.jl # Custom neighbour search and time stepping @@ -89,6 +90,14 @@ To color exported cell grids by particle counts, set `ParticleNeighborsPerCell` array that includes each cell's particle count minus one plus the particles in its neighbor stencil. +To record a property at a specific point, populate `PointMeasures` in the +`SimulationMetaData` configuration. Each `PointMeasure` names a location and a +list of variables to interpolate within the neighbor stencil. These values are +written into a dedicated `*_PointMeasures.vtkhdf` file so you can load probe +points alongside the particle data. Set `ExportMultiBlockVTKHDF=true` to emit a +single `*_MultiBlock.vtkhdf` file that groups particle data and point measures +as separate blocks for ParaView. + ## Help Questions or issues can be posted on the GitHub issue tracker. Response times may vary but all feedback is welcome. diff --git a/example/StillWedgeMDBC.jl b/example/StillWedgeMDBC.jl index 76795e98..1833ee1f 100644 --- a/example/StillWedgeMDBC.jl +++ b/example/StillWedgeMDBC.jl @@ -1,4 +1,5 @@ using SPHExample +using StaticArrays let Dimensions = 2 @@ -15,6 +16,13 @@ let VisualizeInParaview=true, ExportSingleVTKHDF=true, ExportGridCells=true, + PointMeasures=[ + PointMeasure( + "Middle", + SVector{Dimensions, FloatType}(0.5, 0.25), + ["Pressure", "Velocity"], + ), + ], OpenLogFile=true, ) @@ -98,4 +106,3 @@ let # display(plt) end - diff --git a/src/OpenExternalPrograms.jl b/src/OpenExternalPrograms.jl index fdb51e4e..756ec2dc 100644 --- a/src/OpenExternalPrograms.jl +++ b/src/OpenExternalPrograms.jl @@ -70,7 +70,11 @@ function AutoOpenParaview(SimMetaData::SimulationMetaData, OutputVariableNames; if SimMetaData.ExportSingleVTKHDF ParaViewStateFileName = joinpath(SimMetaData.SaveLocation, SimMetaData.SimulationName) * "_SingleVTKHDFStateFile.py" - py_regex = "$(SimMetaData.SimulationName).vtkhdf" + if SimMetaData.ExportMultiBlockVTKHDF + py_regex = "$(SimMetaData.SimulationName)_MultiBlock.vtkhdf" + else + py_regex = "$(SimMetaData.SimulationName).vtkhdf" + end else ParaViewStateFileName = joinpath(SimMetaData.SaveLocation, SimMetaData.SimulationName) * "_StateFile.py" py_regex = "^$(SimMetaData.SimulationName)_(\\d+).vtk" #^ means to anchor the regex to the start of the string diff --git a/src/PointMeasurements.jl b/src/PointMeasurements.jl new file mode 100644 index 00000000..8cdfb227 --- /dev/null +++ b/src/PointMeasurements.jl @@ -0,0 +1,237 @@ +module PointMeasurements + +using LinearAlgebra +using StaticArrays + +using ..SPHKernels: Wᵢⱼ +using ..SPHNeighborList: MapFloor + +export PointMeasure, PointMeasureFieldNames, PointMeasureOutputVariableNames, + FillPointMeasureData!, FillPointMeasureSnapshot! + +""" + PointMeasure(Name, Position, Variables) + +Define a probe point at `Position` for measuring `Variables` during output. +""" +struct PointMeasure{Dimensions, FloatType} + Name::String + Position::SVector{Dimensions, FloatType} + Variables::Vector{String} +end + +@inline function To3DPoint(::Val{2}, Point::SVector{2, T}) where {T} + return SVector{3, T}(Point[1], Point[2], zero(T)) +end + +@inline function To3DPoint(::Val{3}, Point::SVector{3, T}) where {T} + return Point +end + +function PointMeasureFieldNames(PointMeasures::Vector{<:PointMeasure}) + FieldNames = String[] + for Measure in PointMeasures + push!(FieldNames, "PointMeasure_$(Measure.Name)_Position") + for Variable in Measure.Variables + push!(FieldNames, "PointMeasure_$(Measure.Name)_$(Variable)") + end + end + return FieldNames +end + +function PointMeasureOutputVariableNames(PointMeasures::Vector{<:PointMeasure}) + OutputNames = String[] + Seen = Set{String}() + for Measure in PointMeasures + for Variable in Measure.Variables + if !(Variable in Seen) + push!(OutputNames, Variable) + push!(Seen, Variable) + end + end + end + return OutputNames +end + +function FillPointMeasureData!(FieldData, PointMeasures::Vector{<:PointMeasure}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields, + neighbor_data = nothing) + resize!(FieldData, 0) + if isempty(PointMeasures) + return FieldData + end + if neighbor_data === nothing + error("PointMeasure evaluation requires neighbor data; pass `neighbor_data` from the cell list.") + end + + for Measure in PointMeasures + AppendPointMeasureResults!(FieldData, Measure, SimParticles, SimKernel, Dimensions; + field_map = field_map, + vector_fields = vector_fields, + neighbor_data = neighbor_data) + end + + return FieldData +end + +function AppendPointMeasureResults!(FieldData, Measure::PointMeasure{D, T}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields, + neighbor_data) where {D, T} + Positions = SimParticles.Position + NumberOfPoints = length(Positions) + if NumberOfPoints == 0 + error("PointMeasure cannot evaluate without particles.") + end + + Position = Measure.Position + push!(FieldData, To3DPoint(Val(Dimensions), Position)) + + Variables = Measure.Variables + Accumulators = Vector{Any}(undef, length(Variables)) + InterpolateMask = BitVector(undef, length(Variables)) + Sources = Vector{Any}(undef, length(Variables)) + + for (Index, Variable) in pairs(Variables) + Property = get(field_map, Variable, Symbol(Variable)) + if !hasproperty(SimParticles, Property) + error("PointMeasure includes $(Variable) but SimParticles has no field $(Property).") + end + Source = getproperty(SimParticles, Property) + Sources[Index] = Source + if Variable in vector_fields + Accumulators[Index] = zero(eltype(Source)) + InterpolateMask[Index] = eltype(eltype(Source)) <: AbstractFloat + else + Accumulators[Index] = zero(eltype(Source)) + InterpolateMask[Index] = eltype(Source) <: AbstractFloat + end + end + + NearestIndex = 1 + NearestDistance = typemax(T) + WeightSum = zero(T) + + cell_dict = neighbor_data.cell_dict + particle_ranges = neighbor_data.particle_ranges + full_stencil = neighbor_data.full_stencil + cell_index = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(Position))) + + HasCandidates = false + @inbounds for offset in full_stencil + neighbor_cell = cell_index + offset + neighbor_index = get(cell_dict, neighbor_cell, 0) + if neighbor_index != 0 + start_index = particle_ranges[neighbor_index] + end_index = particle_ranges[neighbor_index + 1] - 1 + if start_index <= end_index + HasCandidates = true + for particle_index in start_index:end_index + Offset = Position - Positions[particle_index] + Distance = norm(Offset) + if Distance < NearestDistance + NearestDistance = Distance + NearestIndex = particle_index + end + q = clamp(Distance * SimKernel.h⁻¹, zero(T), T(2)) + if q <= T(2) + Weight = Wᵢⱼ(SimKernel, q) + if Weight != zero(T) + WeightSum += Weight + for VariableIndex in eachindex(Variables) + if InterpolateMask[VariableIndex] + Accumulators[VariableIndex] += Sources[VariableIndex][particle_index] * Weight + end + end + end + end + end + end + end + end + + if !HasCandidates + for (Index, Variable) in pairs(Variables) + if Variable in vector_fields + push!(FieldData, MissingVectorValue(Val(Dimensions), T)) + else + Source = Sources[Index] + push!(FieldData, MissingScalarValue(eltype(Source))) + end + end + return FieldData + end + + for (Index, Variable) in pairs(Variables) + Source = Sources[Index] + Value = if InterpolateMask[Index] && WeightSum > zero(T) + Accumulators[Index] / WeightSum + else + Source[NearestIndex] + end + + if Variable in vector_fields + push!(FieldData, To3DPoint(Val(Dimensions), Value)) + else + push!(FieldData, Value) + end + end + + return FieldData +end + +@inline function MissingScalarValue(::Type{T}) where {T} + return T <: AbstractFloat ? T(NaN) : zero(T) +end + +@inline function MissingVectorValue(::Val{2}, ::Type{T}) where {T} + return SVector{3, T}(MissingScalarValue(T), MissingScalarValue(T), MissingScalarValue(T)) +end + +@inline function MissingVectorValue(::Val{3}, ::Type{T}) where {T} + return SVector{3, T}(MissingScalarValue(T), MissingScalarValue(T), MissingScalarValue(T)) +end + +function FillPointMeasureSnapshot!(snapshot_positions, snapshot_output_data, + PointMeasures::Vector{<:PointMeasure}, + OutputVariables::Vector{String}, + SimParticles, SimKernel, Dimensions; + field_map, vector_fields, + neighbor_data = nothing) + if isempty(PointMeasures) + return snapshot_positions, snapshot_output_data + end + if neighbor_data === nothing + error("PointMeasure evaluation requires neighbor data; pass `neighbor_data` from the cell list.") + end + + resize!(snapshot_positions, length(PointMeasures)) + for Output in snapshot_output_data + resize!(Output, length(PointMeasures)) + fill!(Output, zero(eltype(Output))) + end + + output_index = Dict(name => idx for (idx, name) in pairs(OutputVariables)) + + for (MeasureIndex, Measure) in pairs(PointMeasures) + snapshot_positions[MeasureIndex] = To3DPoint(Val(Dimensions), Measure.Position) + + buffer = Any[] + AppendPointMeasureResults!(buffer, Measure, SimParticles, SimKernel, Dimensions; + field_map = field_map, + vector_fields = vector_fields, + neighbor_data = neighbor_data) + + data_index = 2 + for Variable in Measure.Variables + index = output_index[Variable] + snapshot_output_data[index][MeasureIndex] = buffer[data_index] + data_index += 1 + end + end + + return snapshot_positions, snapshot_output_data +end + +end diff --git a/src/ProduceHDFVTK.jl b/src/ProduceHDFVTK.jl index b1be0eb5..de9e566b 100644 --- a/src/ProduceHDFVTK.jl +++ b/src/ProduceHDFVTK.jl @@ -23,20 +23,27 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, using StaticArrays using ..AuxiliaryFunctions: to_3d! + using ..PointMeasurements: PointMeasureOutputVariableNames, FillPointMeasureSnapshot! const idType = Int64 const fType = Float64 - struct ParticleSnapshot{P, V} + struct ParticleSnapshot{P, V, M} positions::P output_data::V + point_measure_snapshot::M end - struct ParticleWriteJob{P, V} + struct ParticleWriteJob{P, V, M} iteration::Int time::AbstractFloat - snapshot::ParticleSnapshot{P, V} + snapshot::ParticleSnapshot{P, V, M} + end + + struct PointMeasureSnapshot{P, V} + positions::P + output_data::V end struct GridWriteJob{N} @@ -56,21 +63,21 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.write_attribute(attr, dtype, value) end - """ - Compute points and connectivity information for an unstructured grid given - `UniqueCells` from the SPH cell list. Returns `(points, connectivity, - offsets, cell_types, cell_ids, dims)` where `dims` is 2 or 3. - """ - function compute_grid_geometry(SimKernel, UniqueCells) - ExtractDimensionality(::AbstractVector{CartesianIndex{N}}) where N = N - - dims = ExtractDimensionality(UniqueCells) + function create_ordered_group(parent, name::String) + gcpl = HDF5.API.h5p_create(HDF5.API.H5P_GROUP_CREATE) + flags = HDF5.API.H5P_CRT_ORDER_TRACKED | HDF5.API.H5P_CRT_ORDER_INDEXED + HDF5.API.h5pset_link_creation_order(gcpl, flags) + HDF5.API.h5pset_attr_creation_order(gcpl, flags) + group_id = HDF5.API.h5gcreate(parent.id, name, HDF5.API.H5P_DEFAULT, gcpl, HDF5.API.H5P_DEFAULT) + HDF5.API.h5p_close(gcpl) + return HDF5.Group(group_id) + end - dx = dy = dz = SimKernel.H + function create_soft_link(target_path::String, parent, name::String) + HDF5.API.h5l_create_soft(target_path, parent.id, name, HDF5.API.H5P_DEFAULT, HDF5.API.H5P_DEFAULT) + return nothing + end - if dims == 2 - minx, maxx = minimum(ci -> ci[1], UniqueCells), maximum(ci -> ci[1], UniqueCells) - miny, maxy = minimum(ci -> ci[2], UniqueCells), maximum(ci -> ci[2], UniqueCells) nx = maxx - minx + 1 elseif dims == 3 minx, maxx = minimum(ci -> ci[1], UniqueCells), maximum(ci -> ci[1], UniqueCells) @@ -231,6 +238,8 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(pData, var_name, arg_val_type, ((0,),(-1,)), chunk=(chunk_size,)) end end + + FieldData = HDF5.create_group(root, "FieldData") elseif vtk_file_type == "UnstructuredGrid" HDF5.create_dataset(root, "Connectivity" , idType , ((0,),(-1,)), chunk=(chunk_size,)) HDF5.create_dataset(root, "NumberOfCells" , idType , ((0,),(-1,)), chunk=(chunk_size,)) @@ -238,7 +247,7 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, HDF5.create_dataset(root, "Offsets" , idType , ((0,),(-1,)), chunk=(chunk_size,)) HDF5.create_dataset(root, "Types" , UInt8 , ((0,),(-1,)), chunk=(chunk_size,)) #Must be UInt8 - FieldData = HDF5.create_group(root, "FieldData") #Currently just empty group + FieldData = HDF5.create_group(root, "FieldData") CellData = HDF5.create_group(root, "CellData") for name in cell_data_names @@ -359,7 +368,6 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, root["PointData"][var_name][PointsStartIndex:(PointsStartIndex + PositionLength - 1)] = arg end end - connectivities = ["Vertices", "Lines", "Polygons", "Strips"] for connect in connectivities @@ -546,10 +554,42 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, # File naming functions particle_filename = (iter) -> "$(particle_savepath)_$(lpad(iter,6,"0")).vtkhdf" grid_filename = (iter) -> "$(grid_savepath)_$(lpad(iter,6,"0")).vtkhdf" - + point_filename = (iter) -> "$(particle_savepath)_PointMeasures_$(lpad(iter,6,"0")).vtkhdf" + multiblock_filename = "$(particle_savepath)_MultiBlock.vtkhdf" + output_vars = SimMetaData.OutputVariables - + point_measure_vars = PointMeasureOutputVariableNames(SimMetaData.PointMeasures) + + vector_fields = Set([ + "KernelGradient", + "Velocity", + "Acceleration", + "GhostPoints", + "GhostNormals", + ]) + field_map = Dict( + "Kernel" => :Kernel, + "KernelGradient" => :KernelGradient, + "Density" => :Density, + "Pressure" => :Pressure, + "Velocity" => :Velocity, + "Acceleration" => :Acceleration, + "BoundaryBool" => :BoundaryBool, + "ID" => :ID, + "Type" => :Type, + "GroupMarker" => :GroupMarker, + "GhostPoints" => :GhostPoints, + "GhostNormals" => :GhostNormals, + ) + + T = eltype(eltype(SimParticles.Position)) # Initialize storage for file handles + if SimMetaData.ExportMultiBlockVTKHDF && !SimMetaData.ExportSingleVTKHDF + error("ExportMultiBlockVTKHDF requires ExportSingleVTKHDF=true.") + end + + export_multiblock = SimMetaData.ExportMultiBlockVTKHDF + file_handles = if !SimMetaData.ExportSingleVTKHDF # Multi-file mode: vector for particle files n_outputs = if SimMetaData.OutputTimes isa AbstractVector @@ -560,26 +600,17 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, ( particle_files = Vector{HDF5.File}(undef, n_outputs), grid_files = nothing, + point_files = isempty(point_measure_vars) ? nothing : Vector{HDF5.File}(undef, n_outputs), ) else # Single-file mode: handles for both files - OutputVTKHDF = h5open("$(particle_savepath).vtkhdf", "w") - root = HDF5.create_group(OutputVTKHDF, "VTKHDF") - - field_map = Dict( - "Kernel" => :Kernel, - "KernelGradient" => :KernelGradient, - "Density" => :Density, - "Pressure" => :Pressure, - "Velocity" => :Velocity, - "Acceleration" => :Acceleration, - "BoundaryBool" => :BoundaryBool, - "ID" => :ID, - "Type" => :Type, - "GroupMarker" => :GroupMarker, - "GhostPoints" => :GhostPoints, - "GhostNormals" => :GhostNormals, - ) + OutputVTKHDF = h5open(export_multiblock ? multiblock_filename : "$(particle_savepath).vtkhdf", "w") + root = export_multiblock ? create_ordered_group(OutputVTKHDF, "VTKHDF") : HDF5.create_group(OutputVTKHDF, "VTKHDF") + root_particles = create_ordered_group(root, "Particles") + root_points = create_ordered_group(root, "PointMeasures") + assembly = create_ordered_group(root, "Assembly") + create_soft_link("/VTKHDF/PointMeasures", assembly, "PointMeasures") + end output_data_init = Vector{Any}(undef, length(output_vars)) for (i, name) in pairs(output_vars) prop = get(field_map, name, nothing) @@ -593,8 +624,28 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - GenerateGeometryStructure(root, output_vars, output_data_init...; chunk_size=1024) - GenerateStepStructure(root, output_vars, output_data_init...) + root_particles = root + root_points = nothing + + if export_multiblock + HDF5.attrs(root)["Version"] = Int32.([2, 3]) + write_ascii_attribute(root, "Type", "MultiBlock") + blocks = HDF5.create_group(root, "Blocks") + root_particles = HDF5.create_group(blocks, "Particles") + root_points = HDF5.create_group(blocks, "PointMeasures") + end + + GenerateGeometryStructure( + root_particles, + output_vars, + output_data_init...; + chunk_size = 1024, + ) + GenerateStepStructure( + root_particles, + output_vars, + output_data_init...; + ) # Initialize grid file if needed if SimMetaData.ExportGridCells @@ -616,35 +667,56 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, cell_data_names=cell_data_names, ) - (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid) + point_handle = nothing + if !isempty(point_measure_vars) + point_root = root_points + if !export_multiblock + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + end + point_output_init = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_init[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_init[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) + GenerateStepStructure(point_root, point_measure_vars, point_output_init...) + point_handle = export_multiblock ? (file = OutputVTKHDF, root = point_root, shared = true) : (file = OutputVTKHDFPoints, root = point_root, shared = false) + end + + (particle_files = OutputVTKHDF, grid_files = OutputVTKHDFGrid, point_files = point_handle, particle_root = root_particles) else - (particle_files = OutputVTKHDF, grid_files = nothing) + point_handle = nothing + if !isempty(point_measure_vars) + point_root = root_points + if !export_multiblock + OutputVTKHDFPoints = h5open("$(particle_savepath)_PointMeasures.vtkhdf", "w") + point_root = HDF5.create_group(OutputVTKHDFPoints, "VTKHDF") + end + point_output_init = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_init[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_init[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + GenerateGeometryStructure(point_root, point_measure_vars, point_output_init...; chunk_size = 1024) + GenerateStepStructure(point_root, point_measure_vars, point_output_init...) + point_handle = export_multiblock ? (file = OutputVTKHDF, root = point_root, shared = true) : (file = OutputVTKHDFPoints, root = point_root, shared = false) + end + + (particle_files = OutputVTKHDF, grid_files = nothing, point_files = point_handle, particle_root = root_particles) end end - vector_fields = Set([ - "KernelGradient", - "Velocity", - "Acceleration", - "GhostPoints", - "GhostNormals", - ]) - field_map = Dict( - "Kernel" => :Kernel, - "KernelGradient" => :KernelGradient, - "Density" => :Density, - "Pressure" => :Pressure, - "Velocity" => :Velocity, - "Acceleration" => :Acceleration, - "BoundaryBool" => :BoundaryBool, - "ID" => :ID, - "Type" => :Type, - "GroupMarker" => :GroupMarker, - "GhostPoints" => :GhostPoints, - "GhostNormals" => :GhostNormals, - ) - - T = eltype(eltype(SimParticles.Position)) function allocate_particle_snapshot() n = length(SimParticles.Position) positions = Vector{SVector{3, T}}(undef, n) @@ -663,10 +735,22 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, output_data[i] = similar(src, n) end end - return ParticleSnapshot(positions, output_data) + point_positions = Vector{SVector{3, T}}(undef, length(SimMetaData.PointMeasures)) + point_output_data = Vector{Any}(undef, length(point_measure_vars)) + for (index, name) in pairs(point_measure_vars) + if name in vector_fields + point_output_data[index] = fill(zero(SVector{3, T}), length(SimMetaData.PointMeasures)) + else + prop = field_map[name] + src = getproperty(SimParticles, prop) + point_output_data[index] = fill(zero(eltype(src)), length(SimMetaData.PointMeasures)) + end + end + point_snapshot = PointMeasureSnapshot(point_positions, point_output_data) + return ParticleSnapshot(positions, output_data, point_snapshot) end - function fill_particle_snapshot!(snapshot) + function fill_particle_snapshot!(snapshot; neighbor_data = nothing) if Dimensions == 2 to_3d!(snapshot.positions, SimParticles.Position) else @@ -700,6 +784,18 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, copy!(buf, src) end end + FillPointMeasureSnapshot!( + snapshot.point_measure_snapshot.positions, + snapshot.point_measure_snapshot.output_data, + SimMetaData.PointMeasures, + point_measure_vars, + SimParticles, + SimKernel, + Dimensions; + field_map = field_map, + vector_fields = vector_fields, + neighbor_data = neighbor_data, + ) return snapshot end @@ -719,17 +815,37 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, particle_filename(job.iteration), snapshot.positions, output_vars, - snapshot.output_data..., + snapshot.output_data...; ) else AppendVTKHDFData( - root, + file_handles.particle_root, job.time, snapshot.positions, output_vars, - snapshot.output_data..., + snapshot.output_data...; ) end + if !isempty(point_measure_vars) + if !SimMetaData.ExportSingleVTKHDF + SaveVTKHDF( + file_handles.point_files, + job.iteration, + point_filename(job.iteration), + snapshot.point_measure_snapshot.positions, + point_measure_vars, + snapshot.point_measure_snapshot.output_data..., + ) + else + AppendVTKHDFData( + file_handles.point_files.root, + job.time, + snapshot.point_measure_snapshot.positions, + point_measure_vars, + snapshot.point_measure_snapshot.output_data..., + ) + end + end put!(buffer_pool, snapshot) elseif job isa GridWriteJob if !SimMetaData.ExportSingleVTKHDF @@ -754,9 +870,9 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, end end - function enqueue_particle_data(iteration) + function enqueue_particle_data(iteration; neighbor_data = nothing) snapshot = take!(buffer_pool) - fill_particle_snapshot!(snapshot) + fill_particle_snapshot!(snapshot; neighbor_data = neighbor_data) job = ParticleWriteJob(iteration, SimMetaData.TotalTime, snapshot) put!(job_channel, job) end @@ -796,12 +912,22 @@ export SaveVTKHDF, GenerateGeometryStructure, GenerateStepStructure, for f in file_handles.particle_files isopen(f) && close(f) end + if file_handles.point_files !== nothing + for f in file_handles.point_files + isopen(f) && close(f) + end + end else # Close single-file handles isopen(file_handles.particle_files) && close(file_handles.particle_files) if file_handles.grid_files !== nothing isopen(file_handles.grid_files) && close(file_handles.grid_files) end + if file_handles.point_files !== nothing + if !file_handles.point_files.shared + isopen(file_handles.point_files.file) && close(file_handles.point_files.file) + end + end end end diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index eb55a820..094a3076 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -826,7 +826,14 @@ using LinearAlgebra # Save initial state, use 1 else this cannot be used to index fid vector SimMetaData.OutputIterationCounter = 1 - output.enqueue_particles(SimMetaData.OutputIterationCounter) + output.enqueue_particles( + SimMetaData.OutputIterationCounter; + neighbor_data = ( + particle_ranges = ParticleRanges, + cell_dict = CellDict, + full_stencil = FullStencil, + ), + ) if SimMetaData.IndexCounter > 0 unique_cells_view = view(UniqueCells, 1:SimMetaData.IndexCounter) cell_particle_counts = nothing @@ -884,7 +891,14 @@ using LinearAlgebra end @timeit SimMetaData.HourGlass "13 Save Particle Data" begin - output.enqueue_particles(SimMetaData.OutputIterationCounter) + output.enqueue_particles( + SimMetaData.OutputIterationCounter; + neighbor_data = ( + particle_ranges = ParticleRanges, + cell_dict = CellDict, + full_stencil = FullStencil, + ), + ) output.enqueue_grid(SimMetaData.OutputIterationCounter, UniqueCellsView, cell_particle_counts=cell_particle_counts, cell_neighbor_counts=cell_neighbor_counts) end diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 28ef3ae7..2bf4a2c2 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -3,6 +3,8 @@ module SPHExample include("AuxiliaryFunctions.jl"); include("SPHKernels.jl") include("SPHViscosityModels.jl") + include("SPHNeighborList.jl") + include("PointMeasurements.jl") include("ProduceHDFVTK.jl") include("SimulationGeometry.jl") include("SimulationMetaDataConfiguration.jl"); @@ -13,7 +15,6 @@ module SPHExample include("PreProcess.jl"); include("OpenExternalPrograms.jl") include("SPHDensityDiffusionModels.jl") - include("SPHNeighborList.jl") include("SPHCellList.jl") #Must be last # Re-export desired functions from each submodule @@ -23,6 +24,9 @@ module SPHExample using .SPHKernels export SPHKernel, SPHKernelInstance, WendlandC2, CubicSpline, Wᵢⱼ, ∇Wᵢⱼ, tensile_correction + using .PointMeasurements + export PointMeasure, PointMeasureFieldNames, FillPointMeasureData! + using .SPHViscosityModels export SPHViscosity, ZeroViscosity, ArtificialViscosity, Laminar, LaminarSPS, compute_viscosity diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index c0ce65df..ba107d55 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -3,6 +3,8 @@ module SimulationMetaDataConfiguration using Parameters using TimerOutputs +using ..PointMeasurements: PointMeasure + export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, @@ -46,6 +48,7 @@ struct StoreLog <: LogMode end IndexCounter::Int = 0 VisualizeInParaview::Bool = true ExportSingleVTKHDF::Bool = true + ExportMultiBlockVTKHDF::Bool = false ExportGridCells::Bool = false ExportGridCellParticleCounts::Bool = false OutputVariables::Vector{String} = [ @@ -56,6 +59,7 @@ struct StoreLog <: LogMode end "Type", "GroupMarker", ] + PointMeasures::Vector{PointMeasure} = PointMeasure[] OpenLogFile::Bool = true Δx::FloatType = zero(FloatType) end