From 59d7c0b22d08b11274cc910798eeb95f0d5fbefc Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 16 Oct 2024 10:48:48 +0200 Subject: [PATCH 01/43] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index f01673791..dad92016a 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" From 335d3f03a2ae2bb25e0470c3b3c5dddcddf7413f Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Sat, 9 Nov 2024 14:38:41 +0100 Subject: [PATCH 02/43] vtk2trixi --- examples/vtk-reader/read_two_points.jl | 39 +++++++++++++++++++ examples/vtk-reader/write_n_points.jl | 51 +++++++++++++++++++++++++ examples/vtk-reader/write_one_point.jl | 28 ++++++++++++++ examples/vtk-reader/write_two_points.jl | 28 ++++++++++++++ examples/vtk-reader/write_with_trixi.jl | 8 ++++ 5 files changed, 154 insertions(+) create mode 100644 examples/vtk-reader/read_two_points.jl create mode 100644 examples/vtk-reader/write_n_points.jl create mode 100644 examples/vtk-reader/write_one_point.jl create mode 100644 examples/vtk-reader/write_two_points.jl create mode 100644 examples/vtk-reader/write_with_trixi.jl diff --git a/examples/vtk-reader/read_two_points.jl b/examples/vtk-reader/read_two_points.jl new file mode 100644 index 000000000..0baf20b2c --- /dev/null +++ b/examples/vtk-reader/read_two_points.jl @@ -0,0 +1,39 @@ +using ReadVTK + +vtk = VTKFile("out_vtk/two_points.vtu") + +""" +# Fields of vtk + +- `filename`: original path to the VTK file that has been read in +- `xml_file`: object that represents the XML file +- `file_type`: currently only `"UnstructuredGrid"` or `"ImageData"` are supported +- `version`: VTK XML file format version +- `byte_order`: can be `LittleEndian` or `BigEndian` and must currently be the same as the system's +- `compressor`: can be empty (no compression) or `vtkZLibDataCompressor` +- `appended_data`: in case of appended data (see XML documentation), the data is stored here for + convenient retrieval (otherwise it is empty) +- `n_points`: number of points in the VTK file +- `n_cells`: number of cells in the VTK file` +""" + +cell_data = get_cell_data(vtk) + +point_data = get_point_data(vtk) +# create an vector with the names of the point data like "velocity" +names = point_data.names + +# create an vector with the information about the Data like "offset" +data_arrays = point_data.data_arrays + +# information about the file like "UnstructedGrid" +vtk_file = point_data.vtk_file + +# gets the data of the vtk-file +appended_data = vtk.appended_data + +# gets the coords of all the points +coords = get_points(vtk) + +# gets the information about the cells +cells = get_cells(vtk) diff --git a/examples/vtk-reader/write_n_points.jl b/examples/vtk-reader/write_n_points.jl new file mode 100644 index 000000000..a8f5e198b --- /dev/null +++ b/examples/vtk-reader/write_n_points.jl @@ -0,0 +1,51 @@ +using WriteVTK + +# Function to create n points and their associated data +function create_vtk_file(n::Int, coords::Matrix{Float64}, scalars::Vector{Float64}, + velocities::Matrix{Float64}, filename::String) + # Ensure the coordinates matrix has the correct size (3 x n) + @assert size(coords)==(3, n) "coords should be a 3 x n matrix, where n is the number of points" + # Ensure scalars is a vector of length n + @assert length(scalars)==n "scalars should be a vector of length n" + # Ensure velocities is a matrix of size (3 x n) + @assert size(velocities)==(3, n) "velocities should be a 3 x n matrix, where n is the number of points" + + output_directory = "out_vtk" + mkpath(output_directory) + + file = joinpath(output_directory, filename) + + # Create VTK_VERTEX cells for all points + cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in 1:n] + + # Initialize the VTK grid with the points and cells + vtk_file = vtk_grid(file, coords, cells) + + # Assign scalar data to the points + vtk_file["scalar"] = scalars + + # Flatten velocities into a 1D array for VTK + vtk_file["velocity"] = vec(velocities) + + # Save the VTK file + vtk_save(vtk_file) +end + +# Example usage with n points +n = 3 # Number of points + +# Define coordinates for n points (3 rows for x, y, z; n columns for the points) +coords = [1.0 2.0 3.0; # x-coordinates + 1.0 2.0 3.0; # y-coordinates + 0.0 0.0 0.0] # z-coordinates (all on the z=0 plane) + +# Define scalar values for each point +scalars = [100.0, 200.0, 300.0] + +# Define velocity vectors for each point (3 rows for vx, vy, vz; n columns for the points) +velocities = [10.0 0.0 5.0; # x-components + 0.0 5.0 2.0; # y-components + 0.0 0.0 0.0] # z-components (all on the z=0 plane) + +# Create the VTK file +create_vtk_file(n, coords, scalars, velocities, "n_points") diff --git a/examples/vtk-reader/write_one_point.jl b/examples/vtk-reader/write_one_point.jl new file mode 100644 index 000000000..a66f3abd6 --- /dev/null +++ b/examples/vtk-reader/write_one_point.jl @@ -0,0 +1,28 @@ +using WriteVTK + +output_directory = "out_vtk" +mkpath(output_directory) + +file = joinpath(output_directory, "one_point") + +# Define the points array in the expected format: 3 rows (x, y, z), 1 column (for 1 point) +points = [1.0 # x-coordinate + 2.0 # y-coordinate + 0.0] # z-coordinate + +# Create a VTK_VERTEX cell for the single point +cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,))] # Cell for the single point + +# Initialize the VTK grid with the point and cell +vtk_file = vtk_grid(file, points, cells) + +# Assign scalar data to the point +vtk_file["scalar"] = [100.0] # Scalar value for the point + +# Assign vector data (e.g., velocity) to the point +vtk_file["velocity"] = [10.0, 0.0, 0.0] # Velocity for the point + +# Save the VTK file +vtk_save(vtk_file) + +# Surprisingly, the code works not for one point, but for >=2 points. \ No newline at end of file diff --git a/examples/vtk-reader/write_two_points.jl b/examples/vtk-reader/write_two_points.jl new file mode 100644 index 000000000..02bf9e453 --- /dev/null +++ b/examples/vtk-reader/write_two_points.jl @@ -0,0 +1,28 @@ +using WriteVTK + +output_directory = "out_vtk" +mkpath(output_directory) + +file = joinpath(output_directory, "two_points") + +# Define the points array in the expected format: 3 rows (x, y, z), 2 columns (for 2 points) +points = [1.0 3.0 # x-coordinates + 2.0 4.0 # y-coordinates + 0.0 0.0] # z-coordinates (both points on the z=0 plane) + +# Create VTK_VERTEX cells for both points +cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,)), # Cell for first point + MeshCell(VTKCellTypes.VTK_VERTEX, (2,))] # Cell for second point + +# Initialize the VTK grid with the points and cells +vtk_file = vtk_grid(file, points, cells) + +# Assign scalar data to the two points +vtk_file["scalar"] = [100.0, 200.0] # Scalar values for the two points + +# Assign vector data (e.g., velocity) to the two points +vtk_file["velocity"] = [10.0, 0.0, 0.0, # Velocity for first point + 0.0, 5.0, 0.0] # Velocity for second point + +# Save the VTK file +vtk_save(vtk_file) diff --git a/examples/vtk-reader/write_with_trixi.jl b/examples/vtk-reader/write_with_trixi.jl new file mode 100644 index 000000000..bf12d26ed --- /dev/null +++ b/examples/vtk-reader/write_with_trixi.jl @@ -0,0 +1,8 @@ +using TrixiParticles + +rectangle = RectangularShape(0.1, (1, 2), (0.0, 0.0), velocity=[10.0, 0.0], density=1000.0) + +custom_quantity = ones(nparticles(rectangle)) + +trixi2vtk(rectangle; filename="rectangle_of_water", output_directory="out_vtk", + custom_quantity=custom_quantity) \ No newline at end of file From 41696389f4502a7e8642e99930e8ec9ee46d845d Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Sat, 9 Nov 2024 14:38:54 +0100 Subject: [PATCH 03/43] vtk2trixi --- Project.toml | 5 ++ examples/readvtk/trixi2vtk.jl | 8 +++ examples/vtk-reader/particle_saving.jl | 61 +++++++++++++++++++ .../vtk-reader/read_two_points_dynamic.jl | 23 +++++++ ...wo_points.jl => read_two_points_static.jl} | 18 +++++- examples/vtk-reader/write_two_points.jl | 4 +- examples/vtk-reader/write_with_trixi.jl | 2 +- src/TrixiParticles.jl | 2 +- src/preprocessing/preprocessing.jl | 1 + src/preprocessing/readvtk/vtk2trixi.jl | 36 +++++++++++ 10 files changed, 153 insertions(+), 7 deletions(-) create mode 100644 examples/readvtk/trixi2vtk.jl create mode 100644 examples/vtk-reader/particle_saving.jl create mode 100644 examples/vtk-reader/read_two_points_dynamic.jl rename examples/vtk-reader/{read_two_points.jl => read_two_points_static.jl} (80%) create mode 100644 src/preprocessing/readvtk/vtk2trixi.jl diff --git a/Project.toml b/Project.toml index dad92016a..a29efeea3 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.2.3-dev" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -19,9 +20,11 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" @@ -29,6 +32,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" +VTKDataIO = "c6703add-1d23-52c6-9943-3ad88652b9b2" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] @@ -44,6 +48,7 @@ GPUArraysCore = "0.1" JSON = "0.21" KernelAbstractions = "0.9" MuladdMacro = "0.2" +OrdinaryDiffEq = "6.87.0" PointNeighbors = "0.4.2" Polyester = "0.7.10" RecipesBase = "1" diff --git a/examples/readvtk/trixi2vtk.jl b/examples/readvtk/trixi2vtk.jl new file mode 100644 index 000000000..74d8f734d --- /dev/null +++ b/examples/readvtk/trixi2vtk.jl @@ -0,0 +1,8 @@ +# Example File to read a vtk file and convert it to InitialCondition +using TrixiParticles +using OrdinaryDiffEq + +ic = vtk2trixi("out/fluid_1_1.vtu") + +trixi2vtk(ic; filename="trixi2vtk_test", output_directory="out_vtk", + custom_quantity=nothing) diff --git a/examples/vtk-reader/particle_saving.jl b/examples/vtk-reader/particle_saving.jl new file mode 100644 index 000000000..5d1058a23 --- /dev/null +++ b/examples/vtk-reader/particle_saving.jl @@ -0,0 +1,61 @@ +using ReadVTK + +struct Particle + index::Int64 + coordinates::Tuple{Float64, Float64, Float64} # Assuming 3D coordinates + density::Float64 + velocity::Tuple{Float64, Float64, Float64} # Assuming 3D velocity + pressure::Float64 +end + +vtk_file = VTKFile("out/fluid_1_1.vtu") +#vtk_file = VTKFile("out_vtk/rectangle_of_water.vtu") + +# Retrieve particle coordinates +coords = get_points(vtk_file) + +# Retrieve point data fields (e.g., pressure, velocity, ...) +vtk_point_data = get_point_data(vtk_file) + +# Dynamically get all available field names from point_data +field_data = Dict() +for field_name in keys(vtk_point_data) + field_data[field_name] = get_data(vtk_point_data[field_name]) +end + +# Create an array of Particle instances +particles = Vector{Particle}(undef, size(coords, 2)) + +for i in 1:size(coords, 2) + # Retrieve required field "index" + index = field_data["index"][i] + + # Coordinates + coordinates = (coords[1, i], coords[2, i], coords[3, i]) + + # Retrieve each required field directly, assuming all are present + density = field_data["density"][i] + pressure = field_data["pressure"][i] + + velocity = if size(field_data["velocity"], 1) == 2 + # If velocity is 2D, add 0.0 for the z component + (field_data["velocity"][1, i], field_data["velocity"][2, i], 0.0) + else + # If velocity is 3D, use all three components + (field_data["velocity"][1, i], field_data["velocity"][2, i], + field_data["velocity"][3, i]) + end + + # Create a new Particle instance + particles[i] = Particle(index, coordinates, density, velocity, pressure) +end + +# Display some particles for verification +for particle in particles[1:2] + println(particle) +end + +println("Coords: ", size(coords)) +println("velocity ", size(field_data["velocity"])) + +# TODO FieldData \ No newline at end of file diff --git a/examples/vtk-reader/read_two_points_dynamic.jl b/examples/vtk-reader/read_two_points_dynamic.jl new file mode 100644 index 000000000..607dc584c --- /dev/null +++ b/examples/vtk-reader/read_two_points_dynamic.jl @@ -0,0 +1,23 @@ +using ReadVTK + +#vtk = VTKFile("out/fluid_1_0.vtu") +#vtk = VTKFile("out_vtk/two_points.vtu") +vtk = VTKFile("out_vtk/rectangle_of_water.vtu") + +point_data_file = get_point_data(vtk) + +point_names = keys(point_data_file) + +for point_name in point_names + point_array = point_data_file[point_name] + + point_data = get_data(point_array) + + println(point_name) + println(point_data) +end + +coords = get_points(vtk) +for i in 1:size(coords, 2) + println("Coords $i: (", coords[1, i], ", ", coords[2, i], ", ", coords[3, i], ")") +end diff --git a/examples/vtk-reader/read_two_points.jl b/examples/vtk-reader/read_two_points_static.jl similarity index 80% rename from examples/vtk-reader/read_two_points.jl rename to examples/vtk-reader/read_two_points_static.jl index 0baf20b2c..db3927f03 100644 --- a/examples/vtk-reader/read_two_points.jl +++ b/examples/vtk-reader/read_two_points_static.jl @@ -1,6 +1,7 @@ using ReadVTK vtk = VTKFile("out_vtk/two_points.vtu") +#vtk = VTKFile("out/fluid_1_0.vtu") """ # Fields of vtk @@ -17,7 +18,9 @@ vtk = VTKFile("out_vtk/two_points.vtu") - `n_cells`: number of cells in the VTK file` """ -cell_data = get_cell_data(vtk) +""" Examples + +#cell_data = get_cell_data(vtk) point_data = get_point_data(vtk) # create an vector with the names of the point data like "velocity" @@ -33,7 +36,16 @@ vtk_file = point_data.vtk_file appended_data = vtk.appended_data # gets the coords of all the points -coords = get_points(vtk) +#coords = get_points(vtk) # gets the information about the cells -cells = get_cells(vtk) +#cells = get_cells(vtk) +""" +point_data = get_point_data(vtk) + +velocity_array = point_data["velocity"] + +velocity_data = get_data(velocity_array) + +println("Velocity Data:") +println(velocity_data) \ No newline at end of file diff --git a/examples/vtk-reader/write_two_points.jl b/examples/vtk-reader/write_two_points.jl index 02bf9e453..a4d7162b6 100644 --- a/examples/vtk-reader/write_two_points.jl +++ b/examples/vtk-reader/write_two_points.jl @@ -21,8 +21,8 @@ vtk_file = vtk_grid(file, points, cells) vtk_file["scalar"] = [100.0, 200.0] # Scalar values for the two points # Assign vector data (e.g., velocity) to the two points -vtk_file["velocity"] = [10.0, 0.0, 0.0, # Velocity for first point - 0.0, 5.0, 0.0] # Velocity for second point +vtk_file["velocity"] = [1.0, 2.0, 3.0, # Velocity for first point + 4.0, 5.0, 6.0] # Velocity for second point # Save the VTK file vtk_save(vtk_file) diff --git a/examples/vtk-reader/write_with_trixi.jl b/examples/vtk-reader/write_with_trixi.jl index bf12d26ed..ac8260133 100644 --- a/examples/vtk-reader/write_with_trixi.jl +++ b/examples/vtk-reader/write_with_trixi.jl @@ -5,4 +5,4 @@ rectangle = RectangularShape(0.1, (1, 2), (0.0, 0.0), velocity=[10.0, 0.0], dens custom_quantity = ones(nparticles(rectangle)) trixi2vtk(rectangle; filename="rectangle_of_water", output_directory="out_vtk", - custom_quantity=custom_quantity) \ No newline at end of file + custom_quantity=nothing) \ No newline at end of file diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d6e0c2996..58d6b56e3 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -71,7 +71,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx PressureMirroring, PressureZeroing, BoundaryModelLastiwka export BoundaryMovement export examples_dir, validation_dir, trixi_include -export trixi2vtk +export trixi2vtk, vtk2trixi export RectangularTank, RectangularShape, SphereShape, ComplexShape export WindingNumberHormann, WindingNumberJacobson export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry diff --git a/src/preprocessing/preprocessing.jl b/src/preprocessing/preprocessing.jl index f575048b4..353349ff8 100644 --- a/src/preprocessing/preprocessing.jl +++ b/src/preprocessing/preprocessing.jl @@ -1,2 +1,3 @@ include("geometries/geometries.jl") include("point_in_poly/point_in_poly.jl") +include("readvtk/vtk2trixi.jl") diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl new file mode 100644 index 000000000..ef35536b7 --- /dev/null +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -0,0 +1,36 @@ +using ReadVTK + +""" + Convert data from VTK-file to InitialCondition + + # FluidSystem data only +""" + +function vtk2trixi(filename) + vtk_file = VTKFile(filename) + + # Retrieve particle coordinates + coordinates = get_points(vtk_file) + + # Retrieve point data fields (e.g., pressure, velocity, ...) + vtk_point_data = get_point_data(vtk_file) + + # create field data arrays + density = get_data(vtk_point_data["density"]) + + pressure = get_data(vtk_point_data["pressure"]) + + velocity = get_data(vtk_point_data["velocity"]) + if size(velocity, 1) == 2 + # If velocity is 2D, add 0.0 for the z component + velocity = vcat(velocity, zeros(1, size(velocity, 2))) + end + + mass = ones(size(coordinates, 2)) + + return InitialCondition(; coordinates, velocity, mass, density, pressure) +end + +# TODO: edit the mass array --> InitialCondition needs a mass or a particle_spacing +# TODO: example file in folder examples/readvtk +# TODO: make it work with 2D velocity --> In ParaView the velocity vecor is 3D after trixi2vtk. It should be 2D if the initial data is 2D. \ No newline at end of file From 550b066fe2e5af032fa501933596850bfe1a729f Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Sat, 9 Nov 2024 14:43:52 +0100 Subject: [PATCH 04/43] example_file --- examples/readvtk/trixi2vtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/readvtk/trixi2vtk.jl b/examples/readvtk/trixi2vtk.jl index 74d8f734d..272c3a417 100644 --- a/examples/readvtk/trixi2vtk.jl +++ b/examples/readvtk/trixi2vtk.jl @@ -2,7 +2,7 @@ using TrixiParticles using OrdinaryDiffEq -ic = vtk2trixi("out/fluid_1_1.vtu") +ic = vtk2trixi("examples/readvtk/fluid_1_1.vtu") trixi2vtk(ic; filename="trixi2vtk_test", output_directory="out_vtk", custom_quantity=nothing) From caf94397f84e0e1ba0f55e31d7468e5a95bf6ee5 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Fri, 15 Nov 2024 12:06:37 +0100 Subject: [PATCH 05/43] add readvtk --- Project.toml | 5 -- examples/readvtk/readvtk.jl | 10 +++ examples/readvtk/trixi2vtk.jl | 8 --- examples/vtk-reader/particle_saving.jl | 61 ------------------- .../vtk-reader/read_two_points_dynamic.jl | 23 ------- examples/vtk-reader/read_two_points_static.jl | 51 ---------------- examples/vtk-reader/write_n_points.jl | 51 ---------------- examples/vtk-reader/write_one_point.jl | 28 --------- examples/vtk-reader/write_two_points.jl | 28 --------- examples/vtk-reader/write_with_trixi.jl | 8 --- src/TrixiParticles.jl | 1 + src/preprocessing/readvtk/vtk2trixi.jl | 34 +++++------ 12 files changed, 26 insertions(+), 282 deletions(-) create mode 100644 examples/readvtk/readvtk.jl delete mode 100644 examples/readvtk/trixi2vtk.jl delete mode 100644 examples/vtk-reader/particle_saving.jl delete mode 100644 examples/vtk-reader/read_two_points_dynamic.jl delete mode 100644 examples/vtk-reader/read_two_points_static.jl delete mode 100644 examples/vtk-reader/write_n_points.jl delete mode 100644 examples/vtk-reader/write_one_point.jl delete mode 100644 examples/vtk-reader/write_two_points.jl delete mode 100644 examples/vtk-reader/write_with_trixi.jl diff --git a/Project.toml b/Project.toml index 0c9ea908c..7122dda9a 100644 --- a/Project.toml +++ b/Project.toml @@ -5,13 +5,11 @@ version = "0.2.4-dev" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -20,7 +18,6 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -33,7 +30,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" -VTKDataIO = "c6703add-1d23-52c6-9943-3ad88652b9b2" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] @@ -49,7 +45,6 @@ GPUArraysCore = "0.1, 0.2" JSON = "0.21" KernelAbstractions = "0.9" MuladdMacro = "0.2" -OrdinaryDiffEq = "6.87.0" PointNeighbors = "0.4.2" Polyester = "0.7.10" RecipesBase = "1" diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl new file mode 100644 index 000000000..06c9f4201 --- /dev/null +++ b/examples/readvtk/readvtk.jl @@ -0,0 +1,10 @@ +# Example File to read a vtk file and convert it to InitialCondition +using TrixiParticles + +filename = "read_vtk_example" +file = joinpath("examples", "preprocessing", "data", filename * ".vtu") + +ic = vtk2trixi(file) + +trixi2vtk(ic; filename="readvtk", output_directory="out", + custom_quantity=nothing) diff --git a/examples/readvtk/trixi2vtk.jl b/examples/readvtk/trixi2vtk.jl deleted file mode 100644 index 272c3a417..000000000 --- a/examples/readvtk/trixi2vtk.jl +++ /dev/null @@ -1,8 +0,0 @@ -# Example File to read a vtk file and convert it to InitialCondition -using TrixiParticles -using OrdinaryDiffEq - -ic = vtk2trixi("examples/readvtk/fluid_1_1.vtu") - -trixi2vtk(ic; filename="trixi2vtk_test", output_directory="out_vtk", - custom_quantity=nothing) diff --git a/examples/vtk-reader/particle_saving.jl b/examples/vtk-reader/particle_saving.jl deleted file mode 100644 index 5d1058a23..000000000 --- a/examples/vtk-reader/particle_saving.jl +++ /dev/null @@ -1,61 +0,0 @@ -using ReadVTK - -struct Particle - index::Int64 - coordinates::Tuple{Float64, Float64, Float64} # Assuming 3D coordinates - density::Float64 - velocity::Tuple{Float64, Float64, Float64} # Assuming 3D velocity - pressure::Float64 -end - -vtk_file = VTKFile("out/fluid_1_1.vtu") -#vtk_file = VTKFile("out_vtk/rectangle_of_water.vtu") - -# Retrieve particle coordinates -coords = get_points(vtk_file) - -# Retrieve point data fields (e.g., pressure, velocity, ...) -vtk_point_data = get_point_data(vtk_file) - -# Dynamically get all available field names from point_data -field_data = Dict() -for field_name in keys(vtk_point_data) - field_data[field_name] = get_data(vtk_point_data[field_name]) -end - -# Create an array of Particle instances -particles = Vector{Particle}(undef, size(coords, 2)) - -for i in 1:size(coords, 2) - # Retrieve required field "index" - index = field_data["index"][i] - - # Coordinates - coordinates = (coords[1, i], coords[2, i], coords[3, i]) - - # Retrieve each required field directly, assuming all are present - density = field_data["density"][i] - pressure = field_data["pressure"][i] - - velocity = if size(field_data["velocity"], 1) == 2 - # If velocity is 2D, add 0.0 for the z component - (field_data["velocity"][1, i], field_data["velocity"][2, i], 0.0) - else - # If velocity is 3D, use all three components - (field_data["velocity"][1, i], field_data["velocity"][2, i], - field_data["velocity"][3, i]) - end - - # Create a new Particle instance - particles[i] = Particle(index, coordinates, density, velocity, pressure) -end - -# Display some particles for verification -for particle in particles[1:2] - println(particle) -end - -println("Coords: ", size(coords)) -println("velocity ", size(field_data["velocity"])) - -# TODO FieldData \ No newline at end of file diff --git a/examples/vtk-reader/read_two_points_dynamic.jl b/examples/vtk-reader/read_two_points_dynamic.jl deleted file mode 100644 index 607dc584c..000000000 --- a/examples/vtk-reader/read_two_points_dynamic.jl +++ /dev/null @@ -1,23 +0,0 @@ -using ReadVTK - -#vtk = VTKFile("out/fluid_1_0.vtu") -#vtk = VTKFile("out_vtk/two_points.vtu") -vtk = VTKFile("out_vtk/rectangle_of_water.vtu") - -point_data_file = get_point_data(vtk) - -point_names = keys(point_data_file) - -for point_name in point_names - point_array = point_data_file[point_name] - - point_data = get_data(point_array) - - println(point_name) - println(point_data) -end - -coords = get_points(vtk) -for i in 1:size(coords, 2) - println("Coords $i: (", coords[1, i], ", ", coords[2, i], ", ", coords[3, i], ")") -end diff --git a/examples/vtk-reader/read_two_points_static.jl b/examples/vtk-reader/read_two_points_static.jl deleted file mode 100644 index db3927f03..000000000 --- a/examples/vtk-reader/read_two_points_static.jl +++ /dev/null @@ -1,51 +0,0 @@ -using ReadVTK - -vtk = VTKFile("out_vtk/two_points.vtu") -#vtk = VTKFile("out/fluid_1_0.vtu") - -""" -# Fields of vtk - -- `filename`: original path to the VTK file that has been read in -- `xml_file`: object that represents the XML file -- `file_type`: currently only `"UnstructuredGrid"` or `"ImageData"` are supported -- `version`: VTK XML file format version -- `byte_order`: can be `LittleEndian` or `BigEndian` and must currently be the same as the system's -- `compressor`: can be empty (no compression) or `vtkZLibDataCompressor` -- `appended_data`: in case of appended data (see XML documentation), the data is stored here for - convenient retrieval (otherwise it is empty) -- `n_points`: number of points in the VTK file -- `n_cells`: number of cells in the VTK file` -""" - -""" Examples - -#cell_data = get_cell_data(vtk) - -point_data = get_point_data(vtk) -# create an vector with the names of the point data like "velocity" -names = point_data.names - -# create an vector with the information about the Data like "offset" -data_arrays = point_data.data_arrays - -# information about the file like "UnstructedGrid" -vtk_file = point_data.vtk_file - -# gets the data of the vtk-file -appended_data = vtk.appended_data - -# gets the coords of all the points -#coords = get_points(vtk) - -# gets the information about the cells -#cells = get_cells(vtk) -""" -point_data = get_point_data(vtk) - -velocity_array = point_data["velocity"] - -velocity_data = get_data(velocity_array) - -println("Velocity Data:") -println(velocity_data) \ No newline at end of file diff --git a/examples/vtk-reader/write_n_points.jl b/examples/vtk-reader/write_n_points.jl deleted file mode 100644 index a8f5e198b..000000000 --- a/examples/vtk-reader/write_n_points.jl +++ /dev/null @@ -1,51 +0,0 @@ -using WriteVTK - -# Function to create n points and their associated data -function create_vtk_file(n::Int, coords::Matrix{Float64}, scalars::Vector{Float64}, - velocities::Matrix{Float64}, filename::String) - # Ensure the coordinates matrix has the correct size (3 x n) - @assert size(coords)==(3, n) "coords should be a 3 x n matrix, where n is the number of points" - # Ensure scalars is a vector of length n - @assert length(scalars)==n "scalars should be a vector of length n" - # Ensure velocities is a matrix of size (3 x n) - @assert size(velocities)==(3, n) "velocities should be a 3 x n matrix, where n is the number of points" - - output_directory = "out_vtk" - mkpath(output_directory) - - file = joinpath(output_directory, filename) - - # Create VTK_VERTEX cells for all points - cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in 1:n] - - # Initialize the VTK grid with the points and cells - vtk_file = vtk_grid(file, coords, cells) - - # Assign scalar data to the points - vtk_file["scalar"] = scalars - - # Flatten velocities into a 1D array for VTK - vtk_file["velocity"] = vec(velocities) - - # Save the VTK file - vtk_save(vtk_file) -end - -# Example usage with n points -n = 3 # Number of points - -# Define coordinates for n points (3 rows for x, y, z; n columns for the points) -coords = [1.0 2.0 3.0; # x-coordinates - 1.0 2.0 3.0; # y-coordinates - 0.0 0.0 0.0] # z-coordinates (all on the z=0 plane) - -# Define scalar values for each point -scalars = [100.0, 200.0, 300.0] - -# Define velocity vectors for each point (3 rows for vx, vy, vz; n columns for the points) -velocities = [10.0 0.0 5.0; # x-components - 0.0 5.0 2.0; # y-components - 0.0 0.0 0.0] # z-components (all on the z=0 plane) - -# Create the VTK file -create_vtk_file(n, coords, scalars, velocities, "n_points") diff --git a/examples/vtk-reader/write_one_point.jl b/examples/vtk-reader/write_one_point.jl deleted file mode 100644 index a66f3abd6..000000000 --- a/examples/vtk-reader/write_one_point.jl +++ /dev/null @@ -1,28 +0,0 @@ -using WriteVTK - -output_directory = "out_vtk" -mkpath(output_directory) - -file = joinpath(output_directory, "one_point") - -# Define the points array in the expected format: 3 rows (x, y, z), 1 column (for 1 point) -points = [1.0 # x-coordinate - 2.0 # y-coordinate - 0.0] # z-coordinate - -# Create a VTK_VERTEX cell for the single point -cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,))] # Cell for the single point - -# Initialize the VTK grid with the point and cell -vtk_file = vtk_grid(file, points, cells) - -# Assign scalar data to the point -vtk_file["scalar"] = [100.0] # Scalar value for the point - -# Assign vector data (e.g., velocity) to the point -vtk_file["velocity"] = [10.0, 0.0, 0.0] # Velocity for the point - -# Save the VTK file -vtk_save(vtk_file) - -# Surprisingly, the code works not for one point, but for >=2 points. \ No newline at end of file diff --git a/examples/vtk-reader/write_two_points.jl b/examples/vtk-reader/write_two_points.jl deleted file mode 100644 index a4d7162b6..000000000 --- a/examples/vtk-reader/write_two_points.jl +++ /dev/null @@ -1,28 +0,0 @@ -using WriteVTK - -output_directory = "out_vtk" -mkpath(output_directory) - -file = joinpath(output_directory, "two_points") - -# Define the points array in the expected format: 3 rows (x, y, z), 2 columns (for 2 points) -points = [1.0 3.0 # x-coordinates - 2.0 4.0 # y-coordinates - 0.0 0.0] # z-coordinates (both points on the z=0 plane) - -# Create VTK_VERTEX cells for both points -cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (1,)), # Cell for first point - MeshCell(VTKCellTypes.VTK_VERTEX, (2,))] # Cell for second point - -# Initialize the VTK grid with the points and cells -vtk_file = vtk_grid(file, points, cells) - -# Assign scalar data to the two points -vtk_file["scalar"] = [100.0, 200.0] # Scalar values for the two points - -# Assign vector data (e.g., velocity) to the two points -vtk_file["velocity"] = [1.0, 2.0, 3.0, # Velocity for first point - 4.0, 5.0, 6.0] # Velocity for second point - -# Save the VTK file -vtk_save(vtk_file) diff --git a/examples/vtk-reader/write_with_trixi.jl b/examples/vtk-reader/write_with_trixi.jl deleted file mode 100644 index ac8260133..000000000 --- a/examples/vtk-reader/write_with_trixi.jl +++ /dev/null @@ -1,8 +0,0 @@ -using TrixiParticles - -rectangle = RectangularShape(0.1, (1, 2), (0.0, 0.0), velocity=[10.0, 0.0], density=1000.0) - -custom_quantity = ones(nparticles(rectangle)) - -trixi2vtk(rectangle; filename="rectangle_of_water", output_directory="out_vtk", - custom_quantity=nothing) \ No newline at end of file diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 737256bcb..9e7217dca 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -18,6 +18,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch using Printf: @printf, @sprintf +using ReadVTK: VTKFile, get_data, get_point_data, get_points using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index ef35536b7..6352cb726 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -1,18 +1,13 @@ -using ReadVTK - """ - Convert data from VTK-file to InitialCondition - - # FluidSystem data only + Convert data from VTK-file to InitialCondition """ +# TODO: write documentation +# TODO: write tests -function vtk2trixi(filename) - vtk_file = VTKFile(filename) +function vtk2trixi(file) + vtk_file = VTKFile(file) - # Retrieve particle coordinates - coordinates = get_points(vtk_file) - - # Retrieve point data fields (e.g., pressure, velocity, ...) + # retrieve point data fields (e.g., pressure, velocity, ...) vtk_point_data = get_point_data(vtk_file) # create field data arrays @@ -21,16 +16,17 @@ function vtk2trixi(filename) pressure = get_data(vtk_point_data["pressure"]) velocity = get_data(vtk_point_data["velocity"]) - if size(velocity, 1) == 2 - # If velocity is 2D, add 0.0 for the z component - velocity = vcat(velocity, zeros(1, size(velocity, 2))) - end + + # retrieve particle coordinates + # point coordinates are stored in a 3xN matrix, but velocity can be stored either as 3xN or 2xN matrix + coordinates = get_points(vtk_file)[axes(velocity, 1), :] mass = ones(size(coordinates, 2)) + # TODO: read out mass as soon as mass is written out in vtu-file by Trixi + + # TODO: get custom_quantities from vtk file + + # TODO: include the cases, that flieds like velocity are not stored in the vtk file return InitialCondition(; coordinates, velocity, mass, density, pressure) end - -# TODO: edit the mass array --> InitialCondition needs a mass or a particle_spacing -# TODO: example file in folder examples/readvtk -# TODO: make it work with 2D velocity --> In ParaView the velocity vecor is 3D after trixi2vtk. It should be 2D if the initial data is 2D. \ No newline at end of file From 5f37b55b8638111467a940b4246dfa753032704d Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 20 Nov 2024 11:17:53 +0100 Subject: [PATCH 06/43] boundary-fluid-concept --- examples/readvtk/readvtk.jl | 8 +- examples/readvtk/readvtk_testfile.jl | 117 +++++++++++++++++++++++++ src/preprocessing/readvtk/vtk2trixi.jl | 68 ++++++++++---- 3 files changed, 174 insertions(+), 19 deletions(-) create mode 100644 examples/readvtk/readvtk_testfile.jl diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 06c9f4201..c23661229 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,10 +1,14 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -filename = "read_vtk_example" +# the example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with an added 'custom_quantity' +filename = "readvtk_boundary_example" file = joinpath("examples", "preprocessing", "data", filename * ".vtu") -ic = vtk2trixi(file) +# Read the vtk file and convert it to an 'InitialCondition' +# 'system_type' must be "fluid" or "boundary" +system_type = "boundary" +ic = vtk2trixi(file, system_type) trixi2vtk(ic; filename="readvtk", output_directory="out", custom_quantity=nothing) diff --git a/examples/readvtk/readvtk_testfile.jl b/examples/readvtk/readvtk_testfile.jl new file mode 100644 index 000000000..c602c7140 --- /dev/null +++ b/examples/readvtk/readvtk_testfile.jl @@ -0,0 +1,117 @@ +# 2D dam break simulation based on +# +# S. Marrone, M. Antuono, A. Colagrossi, G. Colicchio, D. le Touzé, G. Graziani. +# "δ-SPH model for simulating violent impact flows". +# In: Computer Methods in Applied Mechanics and Engineering, Volume 200, Issues 13–16 (2011), pages 1526–1542. +# https://doi.org/10.1016/J.CMA.2010.12.016 + +using TrixiParticles +using OrdinaryDiffEq + +# Size parameters +H = 0.6 +W = 2 * H + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = H / 40 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +boundary_particle_spacing = fluid_particle_spacing / spacing_ratio + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 + +tspan = (0.0, 1.4) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (W, H) +tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_spacing, 4.0) + +fluid_density = 1000.0 +sound_speed = 20 * sqrt(gravity * H) +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, clip_negative_pressure=false) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + acceleration=(0.0, -gravity), state_equation=state_equation) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +# nu = 0.02 * smoothing_length * sound_speed/8 +# viscosity = ViscosityMorris(nu=nu) +# viscosity = ViscosityAdami(nu=nu) +# Alternatively the density diffusion model by Molteni & Colagrossi can be used, +# which will run faster. +# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity), correction=nothing, + surface_tension=nothing) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + correction=nothing) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) + +# ========================================================================================== +# ==== Simulation +# `nothing` will automatically choose the best update strategy. This is only to be able +# to change this with `trixi_include`. +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +info_callback = InfoCallback(interval=100) + +# custom quantity +kinetic_energy = ones(nparticles(tank.fluid)) +potential_energy = zeros(nparticles(tank.fluid)) + +solution_prefix = "" +saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix, + custom_quantity_1=kinetic_energy, + custom_quantity_2=potential_energy) + +# Save at certain timepoints which allows comparison to the results of Marrone et al., +# i.e. (1.5, 2.36, 3.0, 5.7, 6.45). +# Please note that the images in Marrone et al. are obtained at a particle_spacing = H/320, +# which takes between 2 and 4 hours. +saving_paper = SolutionSavingCallback(save_times=[0.0, 0.371, 0.584, 0.743, 1.411, 1.597], + prefix="marrone_times") + +# This can be overwritten with `trixi_include` +extra_callback = nothing + +use_reinit = false +density_reinit_cb = use_reinit ? + DensityReinitializationCallback(semi.systems[1], interval=10) : + nothing +stepsize_callback = StepsizeCallback(cfl=0.9) + +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback, + density_reinit_cb, saving_paper) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 6352cb726..e8ab10984 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -1,32 +1,66 @@ """ Convert data from VTK-file to InitialCondition """ -# TODO: write documentation -# TODO: write tests +# TODO: Write documentation +# TODO: Write tests + +function vtk2trixi(file, system_type) + # Validate the system type + if !(system_type in ["fluid", "boundary"]) + throw(ArgumentError("Invalid system_type argument. Must be 'fluid' or 'boundary'.")) + end -function vtk2trixi(file) vtk_file = VTKFile(file) - # retrieve point data fields (e.g., pressure, velocity, ...) - vtk_point_data = get_point_data(vtk_file) + # Retrieve point data fields (e.g., pressure, velocity, ...) + point_data = get_point_data(vtk_file) + + coordinates = get_points(vtk_file) + # Check for 2D or 3D coordinates + if all(coordinates[3, :] .== 0) + # If the third row is all zeros, reduce to 2xN + coordinates = coordinates[1:2, :] + end + + if system_type == "fluid" - # create field data arrays - density = get_data(vtk_point_data["density"]) + # Required fields + required_fields = ["density", "pressure", "velocity"] #, "mass"] + missing_fields = [field for field in required_fields if field ∉ keys(point_data)] - pressure = get_data(vtk_point_data["pressure"]) + # Throw an error if any required field is missing + if !isempty(missing_fields) + throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) + end - velocity = get_data(vtk_point_data["velocity"]) + # Retrieve required field + density = get_data(point_data["density"]) + pressure = get_data(point_data["pressure"]) + velocity = get_data(point_data["velocity"]) + # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - # retrieve particle coordinates - # point coordinates are stored in a 3xN matrix, but velocity can be stored either as 3xN or 2xN matrix - coordinates = get_points(vtk_file)[axes(velocity, 1), :] + return InitialCondition(; coordinates, velocity, mass=ones(size(coordinates, 2)), + density, + pressure) - mass = ones(size(coordinates, 2)) - # TODO: read out mass as soon as mass is written out in vtu-file by Trixi + else + # Required fields + required_fields = ["hydrodynamic_density", "pressure"] #, "mass"] + missing_fields = [field for field in required_fields if field ∉ keys(point_data)] - # TODO: get custom_quantities from vtk file + # Throw an error if any required field is missing + if !isempty(missing_fields) + throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) + end - # TODO: include the cases, that flieds like velocity are not stored in the vtk file + # Retrieve required field + hydrodynamic_density = get_data(point_data["hydrodynamic_density"]) + pressure = get_data(point_data["pressure"]) + # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - return InitialCondition(; coordinates, velocity, mass, density, pressure) + return InitialCondition(; coordinates, velocity=zeros(size(coordinates)), + mass=ones(size(coordinates, 2)), + density=hydrodynamic_density, + pressure) + end end From 1d007c012eacf3c9ba3a56394e6a91609cb37c29 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 20 Nov 2024 15:29:53 +0100 Subject: [PATCH 07/43] Fix callback `show` and remove hacky workaround (#669) * Fix callback `show` and remove hacky workaround * Reforamt code * Fix tests * Make it work for all versions --- src/callbacks/solution_saving.jl | 86 +++++++++++++++++++++++-------- test/callbacks/solution_saving.jl | 2 +- 2 files changed, 65 insertions(+), 23 deletions(-) diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index c43cb16da..91dc03ff7 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -106,7 +106,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0, -1, Ref("UnknownVersion")) if length(save_times) > 0 - # See the large comment below for an explanation why we use `finalize` here + # See the large comment below for an explanation why we use `finalize` here. + # When support for Julia 1.9 is dropped, the `finalize` argument can be removed. return PresetTimeCallback(save_times, solution_callback, finalize=solution_callback) elseif dt > 0 # Add a `tstop` every `dt`, and save the final solution @@ -191,31 +192,61 @@ function (solution_callback::SolutionSavingCallback)(integrator) return nothing end -# `finalize` -# This is a hack to make dispatch on a `PresetTimeCallback` possible. -# -# The type of the returned `DiscreteCallback` is +# The type of the `DiscreteCallback` returned by the constructor is # `DiscreteCallback{typeof(condition), typeof(affect!), typeof(initialize), typeof(finalize)}`. -# For the `PeriodicCallback`, `typeof(affect!)` contains the type of the -# `SolutionSavingCallback`. The `PresetTimeCallback` uses anonymous functions as `condition` -# and `affect!`, so this doesn't work here. # -# This hacky workaround makes use of the type parameter `typeof(finalize)` above. -# It's set to `FINALIZE_DEFAULT` by default in the `PresetTimeCallback`, which is a function -# that just returns `nothing`. -# Instead, we pass the `SolutionSavingCallback` as `finalize`, and define it to also just -# return `nothing` when called as `initialize`. +# When `interval` is used, this is +# `DiscreteCallback{<:SolutionSavingCallback, +# <:SolutionSavingCallback, +# typeof(TrixiParticles.initialize_save_cb!), +# typeof(SciMLBase.FINALIZE_DEFAULT)}`. +# +# When `dt` is used, this is +# `DiscreteCallback{DiffEqCallbacks.var"#99#103"{...}, +# DiffEqCallbacks.PeriodicCallbackAffect{<:SolutionSavingCallback}, +# DiffEqCallbacks.var"#100#104"{...} +# typeof(SciMLBase.FINALIZE_DEFAULT)}`. +# +# When `save_times` is used, this is +# `DiscreteCallback{DiffEqCallbacks.var"#115#117"{...}, +# <:SolutionSavingCallback, +# DiffEqCallbacks.var"#116#118"{...}, +# typeof(SciMLBase.FINALIZE_DEFAULT)}}`. +# +# So we can unambiguously dispatch on +# - `DiscreteCallback{<:SolutionSavingCallback, <:SolutionSavingCallback}`, +# - `DiscreteCallback{<:Any, <:PeriodicCallbackAffect{<:SolutionSavingCallback}}`, +# - `DiscreteCallback{<:Any, <:SolutionSavingCallback}`. +# +# WORKAROUND FOR JULIA 1.9: +# When `save_times` is used, the `affect!` is also wrapped in an anonymous function: +# `DiscreteCallback{DiffEqCallbacks.var"#110#113"{...}, +# DiffEqCallbacks.var"#111#114"{<:SolutionSavingCallback}, +# DiffEqCallbacks.var"#116#118"{...}, +# typeof(SciMLBase.FINALIZE_DEFAULT)}}`. +# +# To dispatch here, we set `finalize` to the callback itself, so that the fourth parameter +# becomes `<:SolutionSavingCallback`. This is only used in Julia 1.9. 1.10 and later uses +# a newer version of DiffEqCallbacks.jl that does not have this issue. +# +# To use the callback as `finalize`, we have to define the following function. +# `finalize` is set to `FINALIZE_DEFAULT` by default in the `PresetTimeCallback`, +# which is a function that just returns `nothing`. +# We define the `SolutionSavingCallback` to do the same when called with these arguments. function (finalize::SolutionSavingCallback)(c, u, t, integrator) return nothing end -function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SolutionSavingCallback}) +# With `interval` +function Base.show(io::IO, + cb::DiscreteCallback{<:SolutionSavingCallback, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time solution_saving = cb.affect! print(io, "SolutionSavingCallback(interval=", solution_saving.interval, ")") end +# With `dt` function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:PeriodicCallbackAffect{<:SolutionSavingCallback}}) @@ -225,16 +256,23 @@ function Base.show(io::IO, print(io, "SolutionSavingCallback(dt=", solution_saving.interval, ")") end +# With `save_times`, also working in Julia 1.9. +# When support for Julia 1.9 is dropped, this can be changed to +# `DiscreteCallback{<:Any, <:SolutionSavingCallback}`, and the `finalize` argument +# in the constructor of `SolutionSavingCallback` can be removed. function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:Any, <:Any, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time + # This has to be changed to `cb.affect!` when support for Julia 1.9 is dropped + # and finalize is removed from the constructor of `SolutionSavingCallback`. solution_saving = cb.finalize print(io, "SolutionSavingCallback(save_times=", solution_saving.save_times, ")") end +# With `interval` function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:SolutionSavingCallback}) + cb::DiscreteCallback{<:SolutionSavingCallback, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time if get(io, :compact, false) @@ -257,18 +295,20 @@ function Base.show(io::IO, ::MIME"text/plain", end end +# With `dt` function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:Any, <:Any, <:SolutionSavingCallback}) + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:SolutionSavingCallback}}) @nospecialize cb # reduce precompilation time if get(io, :compact, false) show(io, cb) else - solution_saving = cb.finalize + solution_saving = cb.affect!.affect! cq = collect(solution_saving.custom_quantities) setup = [ - "save_times" => solution_saving.save_times, + "dt" => solution_saving.interval, "custom quantities" => isempty(cq) ? nothing : cq, "save initial solution" => solution_saving.save_initial_solution ? "yes" : "no", @@ -281,19 +321,21 @@ function Base.show(io::IO, ::MIME"text/plain", end end +# With `save_times`. See comments above. function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, - <:PeriodicCallbackAffect{<:SolutionSavingCallback}}) + cb::DiscreteCallback{<:Any, <:Any, <:Any, <:SolutionSavingCallback}) @nospecialize cb # reduce precompilation time if get(io, :compact, false) show(io, cb) else - solution_saving = cb.affect!.affect! + # This has to be changed to `cb.affect!` when support for Julia 1.9 is dropped + # and finalize is removed from the constructor of `SolutionSavingCallback`. + solution_saving = cb.finalize cq = collect(solution_saving.custom_quantities) setup = [ - "dt" => solution_saving.interval, + "save_times" => solution_saving.save_times, "custom quantities" => isempty(cq) ? nothing : cq, "save initial solution" => solution_saving.save_initial_solution ? "yes" : "no", diff --git a/test/callbacks/solution_saving.jl b/test/callbacks/solution_saving.jl index 3cc765b65..7b01c217b 100644 --- a/test/callbacks/solution_saving.jl +++ b/test/callbacks/solution_saving.jl @@ -44,7 +44,7 @@ @test repr("text/plain", callback) == show_box end - @testset verbose=true "interval" begin + @testset verbose=true "save_times" begin callback = SolutionSavingCallback(save_times=[1.0, 2.0, 3.0], prefix="test", output_directory=out) From 07fa6f03c8adbe333ed7dfb86e03b25312f16ab6 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 20 Nov 2024 16:42:26 +0100 Subject: [PATCH 08/43] Add ideal gas equation (#607) * set test up for 1.11 * add basic struct * add export * make equation available * format * typo * fix docs * add exmaple * format * fix * Increase errors for 1.11 * Fix invalidations * Fix tests * Update ci.yml * revert * Update ci.yml * Update test/validation/validation.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * format * review fixes * add test * fix example * fix test * review comments * fix docs * Update src/schemes/fluid/weakly_compressible_sph/state_equations.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * forgot to edit the other doc --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- examples/fluid/dam_break_2phase_2d.jl | 89 +++++++++++++++++++ examples/fluid/dam_break_oil_film_2d.jl | 21 +++-- src/TrixiParticles.jl | 2 +- .../state_equations.jl | 55 +++++++++++- src/visualization/write2vtk.jl | 11 ++- test/examples/examples.jl | 12 +++ .../weakly_compressible_sph/state_equation.jl | 33 +++++++ 7 files changed, 210 insertions(+), 13 deletions(-) create mode 100644 examples/fluid/dam_break_2phase_2d.jl diff --git a/examples/fluid/dam_break_2phase_2d.jl b/examples/fluid/dam_break_2phase_2d.jl new file mode 100644 index 000000000..82b6e343e --- /dev/null +++ b/examples/fluid/dam_break_2phase_2d.jl @@ -0,0 +1,89 @@ +# 2D dam break simulation with an air layer on top + +using TrixiParticles +using OrdinaryDiffEq + +# Size parameters +H = 0.6 +W = 2 * H + +# ========================================================================================== +# ==== Resolution + +# Note: The resolution is very coarse. A better result is obtained with H/60 or higher (which takes over 1 hour) +fluid_particle_spacing = H / 20 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 2.0) + +# Numerical settings +smoothing_length = 3.5 * fluid_particle_spacing +sound_speed = 100.0 +# when using the Ideal gas equation +# sound_speed = 343.0 + +# physical values +nu_water = 8.9E-7 +nu_air = 1.544E-5 +nu_ratio = nu_water / nu_air + +nu_sim_air = 0.02 * smoothing_length * sound_speed +nu_sim_water = nu_ratio * nu_sim_air + +air_viscosity = ViscosityMorris(nu=nu_sim_air) +water_viscosity = ViscosityMorris(nu=nu_sim_water) + +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + sol=nothing, fluid_particle_spacing=fluid_particle_spacing, + viscosity=water_viscosity, smoothing_length=smoothing_length, + gravity=gravity, tspan=tspan, density_diffusion=nothing, + sound_speed=sound_speed, exponent=7, + tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing, + 2.6 * H)) + +# ========================================================================================== +# ==== Setup air_system layer + +air_size = (tank_size[1], 1.5 * H) +air_size2 = (tank_size[1] - W, H) +air_density = 1.0 + +# Air above the initial water +air_system = RectangularShape(fluid_particle_spacing, + round.(Int, air_size ./ fluid_particle_spacing), + zeros(length(air_size)), density=air_density) + +# Air for the rest of the empty volume +air_system2 = RectangularShape(fluid_particle_spacing, + round.(Int, air_size2 ./ fluid_particle_spacing), + (W, 0.0), density=air_density) + +# move on top of the water +for i in axes(air_system.coordinates, 2) + air_system.coordinates[:, i] .+= [0.0, H] +end + +air_system = union(air_system, air_system2) + +air_eos = StateEquationCole(; sound_speed, reference_density=air_density, exponent=1, + clip_negative_pressure=false) +#air_eos = StateEquationIdealGas(; sound_speed, reference_density=air_density, gamma=1.4) + +air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator, + air_eos, smoothing_kernel, smoothing_length, + viscosity=air_viscosity, + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, air_system_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index ecfe66481..b221d8cd6 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -17,15 +17,20 @@ fluid_particle_spacing = H / 60 smoothing_length = 3.5 * fluid_particle_spacing sound_speed = 20 * sqrt(gravity * H) -nu = 0.02 * smoothing_length * sound_speed / 8 -oil_viscosity = ViscosityMorris(nu=10 * nu) +# Physical values +nu_water = 8.9E-7 +nu_oil = 6E-5 +nu_ratio = nu_water / nu_oil + +nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) +nu_sim_water = nu_ratio * nu_sim_oil + +oil_viscosity = ViscosityMorris(nu=nu_sim_oil) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity=ViscosityMorris(nu=nu), smoothing_length=smoothing_length, - gravity=gravity, - density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1), - sound_speed=sound_speed) + viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, + gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed) # ========================================================================================== # ==== Setup oil layer @@ -42,8 +47,10 @@ for i in axes(oil.coordinates, 2) oil.coordinates[:, i] .+= [0.0, H] end +oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density, + exponent=1, clip_negative_pressure=false) oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, - state_equation, smoothing_kernel, + oil_state_equation, smoothing_kernel, smoothing_length, viscosity=oil_viscosity, density_diffusion=density_diffusion, acceleration=(0.0, -gravity), diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ad4e3c478..d85096860 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,7 +65,7 @@ export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel -export StateEquationCole +export StateEquationCole, StateEquationIdealGas export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono diff --git a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl index ea0f7d2ef..a0f9ffdb9 100644 --- a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl +++ b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl @@ -6,10 +6,11 @@ Equation of state to describe the relationship between pressure and density of water up to high pressures. # Keywords -- `sound_speed`: Artificial speed of sound. -- `reference_density`: Reference density of the fluid. -- `exponent`: A value of `7` is usually used for most simulations. +- `sound_speed`: Artificial speed of sound. +- `reference_density`: Reference density of the fluid. +- `exponent`: A value of `7` is usually used for most simulations. - `background_pressure=0.0`: Background pressure. +- `clip_negative_pressure=false`: Negative pressure values are clipped to 0, which prevents spurious surface tension with `SummationDensity` but allows unphysical rarefaction of the fluid. """ struct StateEquationCole{ELTYPE, CLIP} # Boolean to clip negative pressure sound_speed :: ELTYPE @@ -49,3 +50,51 @@ function inverse_state_equation(state_equation::StateEquationCole, pressure) return reference_density * tmp^(1 / exponent) end + +@doc raw""" + StateEquationIdealGas( ;sound_speed, reference_density, gamma, background_pressure=0.0, + clip_negative_pressure=false) + +Equation of state to describe the relationship between pressure and density +of a gas using the Ideal Gas Law. + +# Keywords +- `sound_speed` : Artificial speed of sound. +- `reference_density` : Reference density of the fluid. +- `gamma` : Heat-capacity ratio +- `background_pressure=0.0` : Background pressure. +- `clip_negative_pressure=false`: Negative pressure values are clipped to 0, which prevents spurious surface tension with `SummationDensity` but allows unphysical rarefaction of the fluid. +""" +struct StateEquationIdealGas{ELTYPE, CLIP} + sound_speed :: ELTYPE + reference_density :: ELTYPE + gamma :: ELTYPE + background_pressure :: ELTYPE + + function StateEquationIdealGas(; sound_speed, reference_density, gamma, + background_pressure=0.0, clip_negative_pressure=false) + new{typeof(sound_speed), clip_negative_pressure}(sound_speed, reference_density, + gamma, background_pressure) + end +end + +clip_negative_pressure(::StateEquationIdealGas{<:Any, CLIP}) where {CLIP} = CLIP + +function (state_equation::StateEquationIdealGas)(density) + (; reference_density, sound_speed, gamma, background_pressure) = state_equation + pressure = (density - reference_density) * sound_speed^2 / gamma + background_pressure + + # This is determined statically and has therefore no overhead + if clip_negative_pressure(state_equation) + return max(0.0, pressure) + end + + return pressure +end + +function inverse_state_equation(state_equation::StateEquationIdealGas, pressure) + (; reference_density, sound_speed, gamma, background_pressure) = state_equation + density = (pressure - background_pressure) * gamma / sound_speed^2 + reference_density + + return density +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ef5f68cc6..9f657f733 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -260,12 +260,19 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) if system.correction isa AkinciFreeSurfaceCorrection vtk["correction_rho0"] = system.correction.rho0 end + + if system.state_equation isa StateEquationCole + vtk["state_equation_exponent"] = system.state_equation.exponent + end + + if system.state_equation isa StateEquationIdealGas + vtk["state_equation_gamma"] = system.state_equation.gamma + end + vtk["state_equation"] = type2string(system.state_equation) vtk["state_equation_rho0"] = system.state_equation.reference_density vtk["state_equation_pa"] = system.state_equation.background_pressure vtk["state_equation_c"] = system.state_equation.sound_speed - vtk["state_equation_exponent"] = system.state_equation.exponent - vtk["solver"] = "WCSPH" else vtk["solver"] = "EDAC" diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 51e1dd030..b80e18353 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -203,6 +203,18 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/dam_break_2phase_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2phase_2d.jl"), + tspan=(0.0, 0.05)) [ + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n" + ] + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/dam_break_3d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/schemes/fluid/weakly_compressible_sph/state_equation.jl b/test/schemes/fluid/weakly_compressible_sph/state_equation.jl index 5ece451b8..f86abd47f 100644 --- a/test/schemes/fluid/weakly_compressible_sph/state_equation.jl +++ b/test/schemes/fluid/weakly_compressible_sph/state_equation.jl @@ -139,4 +139,37 @@ end end end + @testset verbose=true "StateEquationIdealGas" begin + # This is the classical ideal gas equation giving a linear relationship + # between density and pressure. + @testset "Physical Properties of Water" begin + # Standard speed of sound for air at 20°C and 1 atm + sound_speed = 343.3 + + # Density of air at 20°C and 1 atm + rest_density = 1.205 + + # Heat capacity ratio of air + gamma = 1.4 + + # Work with pressures in ATM + ATM = 101_325.0 + + state_equation = StateEquationIdealGas(; sound_speed, gamma=gamma, + reference_density=rest_density, + background_pressure=1ATM) + + # Results by manual calculation + @test TrixiParticles.inverse_state_equation(state_equation, 1ATM) == + rest_density + @test TrixiParticles.inverse_state_equation(state_equation, 100ATM) == + 120.36547777058719 + @test TrixiParticles.inverse_state_equation(state_equation, 500ATM) == + 601.8219536113436 + + @test state_equation(rest_density) == 1ATM + @test state_equation(120.36547777058719) == 100ATM + @test state_equation(601.8219536113436) == 500ATM + end + end end From 06a93526f4f02a437d90131cfff2bc9e9c490fa0 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 20 Nov 2024 17:47:14 +0100 Subject: [PATCH 09/43] Add tutorial for setting up a simulation (#514) * Add tutorial for setting up a simulation * Add section about custom smoothing kernels * Implement suggestions * Add kernel plot * Try to remove warnings * Implement suggestions * Make example a dam break and plot initial condition * Fix rendering * Fix tank size * Fix range for kernel plot * Fix descriptions after changing the example * Add equation for kernel * Fix minor errors --- docs/make.jl | 5 + docs/src/examples.md | 21 +- docs/src/systems/weakly_compressible_sph.md | 2 +- docs/src/tutorials_template/tut_setup.md | 317 +++++++++++++++++- .../density_diffusion.jl | 2 +- 5 files changed, 330 insertions(+), 17 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index ac541b0a1..1fb16e935 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -100,6 +100,11 @@ copy_file("NEWS.md") # Define module-wide setups such that the respective modules are available in doctests DocMeta.setdocmeta!(TrixiParticles, :DocTestSetup, :(using TrixiParticles); recursive=true) +# Define environment variables to create plots without warnings +# https://discourse.julialang.org/t/test-plots-on-travis-gks-cant-open-display/9465/2 +ENV["PLOTS_TEST"] = "true" +ENV["GKSwstype"] = "100" + bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib")) makedocs(sitename="TrixiParticles.jl", diff --git a/docs/src/examples.md b/docs/src/examples.md index 57ac2eb7c..f19e623e9 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -1,4 +1,4 @@ -# Examples +# [Examples](@id examples) ## Fluid @@ -11,17 +11,17 @@ ```@raw html ``` - + ### Dam Break 3D (`fluid/dam_break_3d.jl`) ```@raw html ``` - + ### Falling Water Column (`fluid/falling_water_column_2d.jl`) ```@raw html ``` - + ### Hydrostatic Water Column (`fluid/hydrostatic_water_column_*.jl`) ```@raw html @@ -31,38 +31,37 @@ ```@raw html ``` - + ### Oscillating Drop (`fluid/oscillating_drop_2d.jl`) ```@raw html ``` - + ### Periodic Channel (`fluid/periodic_channel_2d.jl`) ```@raw html ``` - + ## Fluid Structure Interaction ### Dam Break with Elastic Plate (`fsi/dam_break_plate_2d.jl`) ```@raw html ``` - + ### Falling Sphere 2D (`fsi/falling_sphere_2d.jl`) ```@raw html ``` - + ### Falling Spheres 2D (`fsi/falling_spheres_2d.jl`) ```@raw html ``` - + ## Structure Mechanics ### Oscillating Beam (`solid/oscillating_beam_2d.jl`) ```@raw html ``` - \ No newline at end of file diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index b4968211a..e561d2cb9 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -53,7 +53,7 @@ Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "viscosity.jl")] ``` -## Density Diffusion +## [Density Diffusion](@id density_diffusion) Density diffusion can be used with [`ContinuityDensity`](@ref) to remove the noise in the pressure field. It is highly recommended to use density diffusion when using WCSPH. diff --git a/docs/src/tutorials_template/tut_setup.md b/docs/src/tutorials_template/tut_setup.md index 0edb01d82..99873dcdf 100644 --- a/docs/src/tutorials_template/tut_setup.md +++ b/docs/src/tutorials_template/tut_setup.md @@ -1,13 +1,322 @@ # Setting up your simulation from scratch -## Hydrostatic tank +In this tutorial, we will guide you through the general structure of simulation files. +We will set up a simulation similar to the example simulation +[`examples/fluid/dam_break_2d.jl`](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/examples/fluid/dam_break_2d.jl), +which is one of our simplest example simulations. +In the second part of this tutorial, we will show how to replace components +of TrixiParticles.jl by custom implementations from within a simulation file, +without ever cloning the repository. +For different setups and physics, have a look at [our other example files](@ref examples). +First, we import TrixiParticles.jl and +[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), which we will +use at the very end for the time integration. +```@example tut_setup +using TrixiParticles +using OrdinaryDiffEq +``` + +## Resolution + +Now, we define the particle spacing, which is our numerical resolution. +For a fluid, we usually call the variable `fluid_particle_spacing`, so that we can easily change +the resolution of an example file by overwriting this variable with [`trixi_include`](@ref). +In 2D, the number of particles will grow quadratically, in 3D cubically with the spacing. + +We also set the number of boundary layers, which need to be sufficiently +large, depending on the smoothing kernel and smoothing length, so that +the compact support of the smoothing kernel is fully sampled with particles +for a fluid particle close to a boundary. +In particular, we require the boundary thickness `boundary_layers * fluid_particle_spacing` +to be larger than the compact support of the kernel. The compact support of each kernel +can be found [in the smoothing kernel overview](@ref smoothing_kernel). +```@example tut_setup +fluid_particle_spacing = 0.02 +boundary_layers = 3 +nothing # hide +``` + +## Experiment setup + +We want to simulate a small dam break problem inside a rectangular tank. +![Experiment Setup](https://github.com/user-attachments/assets/862a1189-b758-4bad-b1e2-6abc42870bb2) +First, we define physical parameters like gravitational acceleration, simulation time, +initial fluid size, tank size and fluid density. +```@example tut_setup +gravity = 9.81 +tspan = (0.0, 1.0) +initial_fluid_size = (1.0, 0.5) +tank_size = (2.0, 1.0) +fluid_density = 1000.0 +nothing # hide +``` + +In order to have the initial particle mass and density correspond to the +hydrostatic pressure gradient, we need to define a [state equation](@ref equation_of_state), +which relates the fluid density to pressure. +Note that we could also skip this part here and define the state equation +later when we define the fluid system, but then the fluid would be initialized +with constant density, which would cause it to oscillate under gravity. +```@example tut_setup +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) +nothing # hide +``` +The speed of sound here is numerical and not physical. +We artificially lower the speed of sound, since the physical speed of sound +in water would lead to prohibitively small time steps. +The speed of sound in Weakly Compressible SPH should be chosen as small as +possible for numerical efficiency, but large enough to limit density fluctuations +to about 1%. + +TrixiParticles.jl requires the initial particle positions and quantities in +form of an [`InitialCondition`](@ref). +Instead of manually defining particle positions, you can work with our +pre-defined setups. +Among others, we provide setups for rectangular shapes, circles, and spheres. +Initial conditions can also be combined with common set operations. +See [this page](@ref initial_condition) for a list of pre-defined setups +and details on set operations on initial conditions. + +Here, we use the [`RectangularTank`](@ref) setup, which generates a rectangular +fluid inside a rectangular tank, and supports a hydrostatic pressure gradient +by passing a gravitational acceleration and a state equation (see above). +```@example tut_setup +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, + fluid_density, n_layers=boundary_layers, + acceleration=(0.0, -gravity), state_equation=state_equation) +nothing # hide +``` +A `RectangularTank` consists of two [`InitialCondition`](@ref)s, `tank.fluid` and `tank.boundary`. +We can plot these initial conditions to visualize the initial setup. +```@example tut_setup +using Plots +plot(tank.fluid, tank.boundary, labels=["fluid" "boundary"]) +plot!(dpi=200); savefig("tut_setup_plot_tank.png"); nothing # hide +``` +![plot tank](tut_setup_plot_tank.png) + +## Fluid system + +To model the water column, we use the [Weakly Compressible Smoothed Particle +Hydrodynamics (WCSPH) method](@ref wcsph). +This method requires a smoothing kernel and a corresponding smoothing length, +which should be chosen in relation to the particle spacing. +```@example tut_setup +smoothing_length = 1.2 * fluid_particle_spacing +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +nothing # hide +``` +You can find an overview over smoothing kernels and corresponding smoothing +lengths [here](@ref smoothing_kernel). + +For stability, we need numerical dissipation in form of an artificial viscosity term. +[Other viscosity models](@ref viscosity_wcsph) offer a physical approach +based on the kinematic viscosity of the fluid. +```@example tut_setup +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +nothing # hide +``` +We choose the parameters as small as possible to avoid non-physical behavior, +but as large as possible to stabilize the simulation. + +The WCSPH method can either compute the particle density directly with a kernel summation +over all neighboring particles (see [`SummationDensity`](@ref)) or by making +the particle density a variable in the ODE system and integrating its change over time. +We choose the latter approach here by using the density calculator +[`ContinuityDensity`](@ref), which is more efficient and handles free surfaces +without the need for additional correction terms. +The simulation quality greatly benefits from using [density diffusion](@ref density_diffusion). +```@example tut_setup +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) +nothing # hide +``` + +## Boundary system + +To model the boundary, we use particle-based boundary conditions, in which particles +are sampled in the boundary that interact with the fluid particles to avoid penetration. +In order to define a boundary system, we first have to choose a boundary model, +which defines how the fluid interacts with boundary particles. +We will use the [`BoundaryModelDummyParticles`](@ref) with [`AdamiPressureExtrapolation`](@ref). +See [here](@ref boundary_models) for a comprehensive overview over boundary models. +```@example tut_setup +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) +nothing # hide +``` + +## Semidiscretization + +The key component of every simulation is the [`Semidiscretization`](@ref), +which couples all systems of the simulation. +All simulation methods in TrixiParticles.jl are semidiscretizations, which discretize +the equations in space to provide an ordinary differential equation that still +has to be solved in time. +By providing a simulation time span, we can call [`semidiscretize`](@ref), +which returns an `ODEProblem` that can be solved with a time integration method. +```@example tut_setup +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) +nothing # hide +``` +## Time integration -# Example file -```julia -!!include:examples/fluid/hydrostatic_water_column_2d.jl!! +We use the methods provided by +[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), +but note that other packages or custom implementations can also be used. +OrdinaryDiffEq.jl supports callbacks, which are executed during the simulation. +For this simulation, we use the [`InfoCallback`](@ref), which prints +information about the simulation setup at the beginning of the simulation, +information about the current simulation time and runtime during the simulation, +and a performance summary at the end of the simulation. +We also want to save the current solution in regular intervals in terms of +simulation time as VTK, so that we can [look at the solution in ParaView](@ref Visualization). +The [`SolutionSavingCallback`](@ref) provides this functionality. +To pass the callbacks to OrdinaryDiffEq.jl, we have to bundle them into a +`CallbackSet`. +```@example tut_setup +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02) + +callbacks = CallbackSet(info_callback, saving_callback) +nothing # hide +``` + +Finally, we can start the simulation by solving the `ODEProblem`. +We use the method `RDPK3SpFSAL35` of OrdinaryDiffEq.jl, which is a Runge-Kutta +method with automatic (error based) time step size control. +This method is usually a good choice for prototyping, since we do not have to +worry about choosing a stable step size and can just run the simulation. +For better performance, it might be beneficial to tweak the tolerances +of this method or choose a different method that is more efficient for the +respective simulation. +You can find both approaches in our [example files](@ref examples). +Here, we just use the method with the default parameters, and only disable +`save_everystep` to avoid expensive saving of the solution in every time step. +```@example tut_setup +sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); +nothing # hide +``` + +See [Visualization](@ref) for how to visualize the final solution. +For the simplest visualization, we can use [Plots.jl](https://docs.juliaplots.org/stable/): +```@example tut_setup +using Plots +plot(sol) +plot!(dpi=200); savefig("tut_setup_plot.png"); nothing # hide +``` +![plot](tut_setup_plot.png) + +## Replacing components with custom implementations + +If we would like to use an implementation of a component that is not available +in TrixiParticles.jl, we can implement it ourselves within the simulation file, +without ever cloning the TrixiParticles.jl repository. +A good starting point is to check out the available implementations in +TrixiParticles.jl, then copy the relevant functions to the simulation file +and modify them as needed. + +### Custom smoothing kernel + +To implement a custom smoothing kernel, we define a struct extending +`TrixiParticles.SmoothingKernel`. +This abstract struct has a type parameter for the number of dimensions, +which we set to 2 in this case. + +```@example tut_setup +struct MyGaussianKernel <: TrixiParticles.SmoothingKernel{2} end ``` +This kernel is going to be an implementation of the Gaussian kernel with +a cutoff for compact support, which reads +```math +W(r, h) = +\begin{cases} +\frac{1}{\pi h^2} \exp(-(r/h)^2) & \text{for } r < 2h\\ +0 & \text{for } r \geq 2h. +\end{cases} +``` +Note that the same kernel in a more optimized version and with a cutoff at ``3`` +is already implemented in TrixiParticles.jl as [`GaussianKernel`](@ref). + +In order to use our new kernel, we have to define three functions. +`TrixiParticles.kernel`, which is the kernel function itself, +`TrixiParticles.kernel_deriv`, which is the derivative of the kernel function, +and `TrixiParticles.compact_support`, which defines the compact support of the +kernel in relation to the smoothing length. +The latter is relevant for determining the search radius of the neighborhood search. + +```@example tut_setup +function TrixiParticles.kernel(kernel::MyGaussianKernel, r, h) + q = r / h + if q < 2 + return 1 / (pi * h^2) * exp(-q^2) + end + + return 0.0 +end + +function TrixiParticles.kernel_deriv(kernel::MyGaussianKernel, r, h) + q = r / h + + if q < 2 + return 1 / (pi * h^2) * (-2 * q) * exp(-q^2) / h + end + + return 0.0 +end + +TrixiParticles.compact_support(::MyGaussianKernel, h) = 2 * h +``` + +For this kernel, we use a different smoothing length, which yields a similar kernel +to the `SchoenbergCubicSplineKernel` that we used earlier. +```@example tut_setup +smoothing_length_gauss = 1.0 * fluid_particle_spacing +nothing # hide +``` +We can compare these kernels in a plot. +```@example tut_setup +using Plots +x = range(-0.05, 0.05, length=500) +plot(x, r -> TrixiParticles.kernel(SchoenbergCubicSplineKernel{2}(), abs(r), smoothing_length), + label="SchoenbergCubicSplineKernel", xlabel="r") +plot!(x, r -> TrixiParticles.kernel(MyGaussianKernel(), abs(r), smoothing_length_gauss), + label="MyGaussianKernel") +plot!(dpi=200); savefig("tut_setup_plot2.png"); nothing # hide +``` +![plot](tut_setup_plot2.png) + +This is all we need to use our custom kernel implementation in a simulation. +We only need to replace the definition above by +```@example tut_setup +smoothing_kernel = MyGaussianKernel() +nothing # hide +``` +and run the simulation file again. + +In order to use our kernel in a pre-defined example file, we can use the function +[`trixi_include`](@ref) to replace the definition of the variable `smoothing_kernel`. +The following will run the example simulation +`examples/fluid/hydrostatic_water_column_2d.jl` with our custom kernel and the corresponding +smoothing length. +```@example tut_setup +trixi_include(joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), + smoothing_kernel=MyGaussianKernel(), smoothing_length=smoothing_length_gauss); +nothing # hide +``` diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c10740ab8..8a04c88bf 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -11,7 +11,7 @@ Currently, the following formulations are available: | [`DensityDiffusionFerrari`](@ref) | ❌ | ✅ | | [`DensityDiffusionAntuono`](@ref) | ✅ | ❌ | -See [Density Diffusion](@ref) for a comparison and more details. +See [Density Diffusion](@ref density_diffusion) for a comparison and more details. """ abstract type DensityDiffusion end From f13fef2ff49436d9fe96b0544ef3578c943b202d Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Fri, 22 Nov 2024 11:32:12 +0100 Subject: [PATCH 10/43] compact-concept --- examples/readvtk/readvtk.jl | 9 +- examples/readvtk/readvtk_testfile.jl | 117 ------------------------- src/preprocessing/readvtk/vtk2trixi.jl | 61 +++++-------- 3 files changed, 29 insertions(+), 158 deletions(-) delete mode 100644 examples/readvtk/readvtk_testfile.jl diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index c23661229..8c1693e17 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,13 +1,16 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -# the example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with an added 'custom_quantity' -filename = "readvtk_boundary_example" +# The example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with added 'custom_quantities' +filename = "readvtk_fluid_example" file = joinpath("examples", "preprocessing", "data", filename * ".vtu") +# filename = "boundary_1_1" +# file = joinpath("out", filename * ".vtu") + # Read the vtk file and convert it to an 'InitialCondition' # 'system_type' must be "fluid" or "boundary" -system_type = "boundary" +system_type = "fluid" ic = vtk2trixi(file, system_type) trixi2vtk(ic; filename="readvtk", output_directory="out", diff --git a/examples/readvtk/readvtk_testfile.jl b/examples/readvtk/readvtk_testfile.jl deleted file mode 100644 index c602c7140..000000000 --- a/examples/readvtk/readvtk_testfile.jl +++ /dev/null @@ -1,117 +0,0 @@ -# 2D dam break simulation based on -# -# S. Marrone, M. Antuono, A. Colagrossi, G. Colicchio, D. le Touzé, G. Graziani. -# "δ-SPH model for simulating violent impact flows". -# In: Computer Methods in Applied Mechanics and Engineering, Volume 200, Issues 13–16 (2011), pages 1526–1542. -# https://doi.org/10.1016/J.CMA.2010.12.016 - -using TrixiParticles -using OrdinaryDiffEq - -# Size parameters -H = 0.6 -W = 2 * H - -# ========================================================================================== -# ==== Resolution -fluid_particle_spacing = H / 40 - -# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model -boundary_layers = 4 -spacing_ratio = 1 - -boundary_particle_spacing = fluid_particle_spacing / spacing_ratio - -# ========================================================================================== -# ==== Experiment Setup -gravity = 9.81 - -tspan = (0.0, 1.4) - -# Boundary geometry and initial fluid particle positions -initial_fluid_size = (W, H) -tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_spacing, 4.0) - -fluid_density = 1000.0 -sound_speed = 20 * sqrt(gravity * H) -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=1, clip_negative_pressure=false) - -tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - n_layers=boundary_layers, spacing_ratio=spacing_ratio, - acceleration=(0.0, -gravity), state_equation=state_equation) - -# ========================================================================================== -# ==== Fluid -smoothing_length = 3.5 * fluid_particle_spacing -smoothing_kernel = WendlandC2Kernel{2}() - -fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) -# nu = 0.02 * smoothing_length * sound_speed/8 -# viscosity = ViscosityMorris(nu=nu) -# viscosity = ViscosityAdami(nu=nu) -# Alternatively the density diffusion model by Molteni & Colagrossi can be used, -# which will run faster. -# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) -density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) - -fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, - state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, - density_diffusion=density_diffusion, - acceleration=(0.0, -gravity), correction=nothing, - surface_tension=nothing) - -# ========================================================================================== -# ==== Boundary -boundary_density_calculator = AdamiPressureExtrapolation() -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - state_equation=state_equation, - boundary_density_calculator, - smoothing_kernel, smoothing_length, - correction=nothing) - -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) - -# ========================================================================================== -# ==== Simulation -# `nothing` will automatically choose the best update strategy. This is only to be able -# to change this with `trixi_include`. -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing)) -ode = semidiscretize(semi, tspan, data_type=nothing) - -info_callback = InfoCallback(interval=100) - -# custom quantity -kinetic_energy = ones(nparticles(tank.fluid)) -potential_energy = zeros(nparticles(tank.fluid)) - -solution_prefix = "" -saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix, - custom_quantity_1=kinetic_energy, - custom_quantity_2=potential_energy) - -# Save at certain timepoints which allows comparison to the results of Marrone et al., -# i.e. (1.5, 2.36, 3.0, 5.7, 6.45). -# Please note that the images in Marrone et al. are obtained at a particle_spacing = H/320, -# which takes between 2 and 4 hours. -saving_paper = SolutionSavingCallback(save_times=[0.0, 0.371, 0.584, 0.743, 1.411, 1.597], - prefix="marrone_times") - -# This can be overwritten with `trixi_include` -extra_callback = nothing - -use_reinit = false -density_reinit_cb = use_reinit ? - DensityReinitializationCallback(semi.systems[1], interval=10) : - nothing -stepsize_callback = StepsizeCallback(cfl=0.9) - -callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback, - density_reinit_cb, saving_paper) - -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # This is overwritten by the stepsize callback - save_everystep=false, callback=callbacks); diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index e8ab10984..ad3875a7a 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -22,45 +22,30 @@ function vtk2trixi(file, system_type) coordinates = coordinates[1:2, :] end - if system_type == "fluid" - - # Required fields - required_fields = ["density", "pressure", "velocity"] #, "mass"] - missing_fields = [field for field in required_fields if field ∉ keys(point_data)] - - # Throw an error if any required field is missing - if !isempty(missing_fields) - throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) - end - - # Retrieve required field - density = get_data(point_data["density"]) - pressure = get_data(point_data["pressure"]) - velocity = get_data(point_data["velocity"]) - # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - - return InitialCondition(; coordinates, velocity, mass=ones(size(coordinates, 2)), - density, - pressure) - - else - # Required fields - required_fields = ["hydrodynamic_density", "pressure"] #, "mass"] - missing_fields = [field for field in required_fields if field ∉ keys(point_data)] + # Define required fields based on system type + required_fields, density_key = system_type == "fluid" ? + (["density", "pressure", "velocity"], "density") : + (["hydrodynamic_density", "pressure"], + "hydrodynamic_density") + + # Check for missing fields + missing_fields = [field for field in required_fields if field ∉ keys(point_data)] + if !isempty(missing_fields) + throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) + end - # Throw an error if any required field is missing - if !isempty(missing_fields) - throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) - end + # Retrieve fields + density = get_data(point_data[density_key]) + pressure = get_data(point_data["pressure"]) + velocity = system_type == "fluid" ? get_data(point_data["velocity"]) : + zeros(size(coordinates)) - # Retrieve required field - hydrodynamic_density = get_data(point_data["hydrodynamic_density"]) - pressure = get_data(point_data["pressure"]) - # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi + # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - return InitialCondition(; coordinates, velocity=zeros(size(coordinates)), - mass=ones(size(coordinates, 2)), - density=hydrodynamic_density, - pressure) - end + return InitialCondition( + ; coordinates, + velocity, + mass=ones(size(coordinates, 2)), + density, + pressure) end From 231e2e020cf484ccb46fc28d0eee39c4490aef91 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Fri, 22 Nov 2024 15:13:20 +0100 Subject: [PATCH 11/43] Validation TVF (#640) * add `density_calculator` to example * add perturbation option * fix typo * implement suggestions --------- Co-authored-by: Sven Berger --- examples/fluid/taylor_green_vortex_2d.jl | 7 ++++++- .../validation_taylor_green_vortex_2d.jl | 15 +++++++++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 551863ccb..fc4ad9e63 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -58,13 +58,18 @@ background_pressure = sound_speed^2 * fluid_density smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() +# To be set via `trixi_include` +perturb_coordinates = true fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), - coordinates_perturbation=0.2, # To avoid stagnant streamlines when not using TVF. + # Perturb particle coordinates to avoid stagnant streamlines without TVF + coordinates_perturbation=perturb_coordinates ? 0.2 : nothing, # To avoid stagnant streamlines when not using TVF. density=fluid_density, pressure=initial_pressure_function, velocity=initial_velocity_function) +density_calculator = SummationDensity() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, + density_calculator=density_calculator, transport_velocity=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index cd9d18159..99b83ca03 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -9,6 +9,9 @@ particle_spacings = [0.02, 0.01, 0.005] tspan = (0.0, 5.0) reynolds_number = 100.0 +density_calculators = [SummationDensity(), ContinuityDensity()] +perturb_coordinates = [false, true] + function compute_l1v_error(v, u, t, system) v_analytical_avg = 0.0 v_avg = 0.0 @@ -60,10 +63,16 @@ function diff_p_loc_p_avg(v, u, t, system) return v[end, :] .- p_avg_tot end -for particle_spacing in particle_spacings +for density_calculator in density_calculators, perturbation in perturb_coordinates, + particle_spacing in particle_spacings + n_particles_xy = round(Int, 1.0 / particle_spacing) - output_directory = joinpath("out_tgv", + name_density_calculator = density_calculator isa SummationDensity ? + "summation_density" : "continuity_density" + name_perturbation = perturbation ? "perturbed" : "" + + output_directory = joinpath("out_tgv", "$(name_density_calculator)_$name_perturbation", "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory, @@ -79,6 +88,8 @@ for particle_spacing in particle_spacings # Import variables into scope trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "taylor_green_vortex_2d.jl"), + density_calculator=density_calculator, + perturb_coordinates=perturbation, particle_spacing=particle_spacing, reynolds_number=reynolds_number, tspan=tspan, saving_callback=saving_callback, pp_callback=pp_callback) end From e28b79b3e85ba4b36585e332cd3c123f2cc6ad87 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 08:33:19 +0100 Subject: [PATCH 12/43] integrate try and catch --- examples/readvtk/readvtk.jl | 26 ++++++---- src/TrixiParticles.jl | 2 +- src/preprocessing/readvtk/vtk2trixi.jl | 66 ++++++++++++++------------ src/visualization/write2vtk.jl | 1 + 4 files changed, 56 insertions(+), 39 deletions(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 8c1693e17..8ff760ea3 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,17 +1,27 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -# The example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with added 'custom_quantities' -filename = "readvtk_fluid_example" -file = joinpath("examples", "preprocessing", "data", filename * ".vtu") +coordinates = [0.0 0.0 0.0; 1.0 0.0 0.0] +initial_condition = InitialCondition(; coordinates, density=1.0, particle_spacing=0.1) +trixi2vtk(initial_condition; filename="initial_condition", output_directory="out", + custom_quantity=nothing) + +filename = "initial_condition" +file = joinpath("out", filename * ".vtu") -# filename = "boundary_1_1" -# file = joinpath("out", filename * ".vtu") +#rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.0) +#trixi2vtk(rectangle; filename="rectangle", output_directory="out", +# custom_quantity=nothing) + +#filename = "rectangle" +#file = joinpath("out", filename * ".vtu") + +# The example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with added 'custom_quantities' +#filename = "fluid_1_6" +#file = joinpath("out", filename * ".vtu") # Read the vtk file and convert it to an 'InitialCondition' -# 'system_type' must be "fluid" or "boundary" -system_type = "fluid" -ic = vtk2trixi(file, system_type) +ic = vtk2trixi(file) trixi2vtk(ic; filename="readvtk", output_directory="out", custom_quantity=nothing) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 1b4ee8ded..67b6bfa0b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -19,7 +19,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch using Printf: @printf, @sprintf -using ReadVTK: VTKFile, get_data, get_point_data, get_points +using ReadVTK: VTKFile, get_data, get_field_data, get_point_data, get_points using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index ad3875a7a..1f96d8aee 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -1,46 +1,52 @@ +# TODO: Write documentation +# TODO: Write tests """ Convert data from VTK-file to InitialCondition """ -# TODO: Write documentation -# TODO: Write tests - -function vtk2trixi(file, system_type) - # Validate the system type - if !(system_type in ["fluid", "boundary"]) - throw(ArgumentError("Invalid system_type argument. Must be 'fluid' or 'boundary'.")) - end - +function vtk2trixi(file; iter=0, input_directory="out", prefix="") vtk_file = VTKFile(file) # Retrieve point data fields (e.g., pressure, velocity, ...) point_data = get_point_data(vtk_file) - - coordinates = get_points(vtk_file) - # Check for 2D or 3D coordinates - if all(coordinates[3, :] .== 0) - # If the third row is all zeros, reduce to 2xN - coordinates = coordinates[1:2, :] + meta_data = get_field_data(vtk_file) + # TODO: Shapes created directly with write2vtk do not have 'meta_data'. Add this feature + + NDIMS = get_data(meta_data["ndims"]) + coordinates = get_points(vtk_file)[1:NDIMS[1], :] + + # Retrieve fields + density = try + get_data(point_data["density"]) + catch + try + get_data(point_data["hydrodynamic_density"]) + catch + zeros(size(coordinates, 2)) + end end - # Define required fields based on system type - required_fields, density_key = system_type == "fluid" ? - (["density", "pressure", "velocity"], "density") : - (["hydrodynamic_density", "pressure"], - "hydrodynamic_density") - - # Check for missing fields - missing_fields = [field for field in required_fields if field ∉ keys(point_data)] - if !isempty(missing_fields) - throw(ArgumentError("The following required fields are missing in the VTK file: $(missing_fields)")) + pressure = try + get_data(point_data["pressure"]) + catch + zeros(size(coordinates, 2)) end - # Retrieve fields - density = get_data(point_data[density_key]) - pressure = get_data(point_data["pressure"]) - velocity = system_type == "fluid" ? get_data(point_data["velocity"]) : - zeros(size(coordinates)) + velocity = try + get_data(point_data["velocity"]) + catch + try + get_data(point_data["initial_velocity"]) + catch + zeros(size(coordinates)) + end + end # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi + # mass = try + # get_data(point_data["mass"]) + # catch + # zeros(size(coordinates, 2)) + # end return InitialCondition( ; coordinates, diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 9f657f733..82ad2ee7e 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -123,6 +123,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre # Store particle index vtk["index"] = active_particles(system) vtk["time"] = t + vtk["ndims"] = ndims(system) if write_meta_data vtk["solver_version"] = git_hash From 499e2859343cb130b00043e7351c94a773e3c9da Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 10:13:51 +0100 Subject: [PATCH 13/43] added wall_velocity case --- examples/readvtk/readvtk.jl | 12 ++---------- src/preprocessing/readvtk/vtk2trixi.jl | 25 +++++++++++-------------- src/visualization/write2vtk.jl | 2 +- 3 files changed, 14 insertions(+), 25 deletions(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 8ff760ea3..cebcb3d55 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,14 +1,6 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -coordinates = [0.0 0.0 0.0; 1.0 0.0 0.0] -initial_condition = InitialCondition(; coordinates, density=1.0, particle_spacing=0.1) -trixi2vtk(initial_condition; filename="initial_condition", output_directory="out", - custom_quantity=nothing) - -filename = "initial_condition" -file = joinpath("out", filename * ".vtu") - #rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.0) #trixi2vtk(rectangle; filename="rectangle", output_directory="out", # custom_quantity=nothing) @@ -17,8 +9,8 @@ file = joinpath("out", filename * ".vtu") #file = joinpath("out", filename * ".vtu") # The example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with added 'custom_quantities' -#filename = "fluid_1_6" -#file = joinpath("out", filename * ".vtu") +filename = "boundary_1_6" +file = joinpath("out", filename * ".vtu") # Read the vtk file and convert it to an 'InitialCondition' ic = vtk2trixi(file) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 1f96d8aee..1433beaec 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -9,26 +9,18 @@ function vtk2trixi(file; iter=0, input_directory="out", prefix="") # Retrieve point data fields (e.g., pressure, velocity, ...) point_data = get_point_data(vtk_file) meta_data = get_field_data(vtk_file) - # TODO: Shapes created directly with write2vtk do not have 'meta_data'. Add this feature + # Meta data only written out in simulations NDIMS = get_data(meta_data["ndims"]) coordinates = get_points(vtk_file)[1:NDIMS[1], :] - # Retrieve fields + # Retrieve fields + pressure = get_data(point_data["pressure"]) + density = try get_data(point_data["density"]) catch - try - get_data(point_data["hydrodynamic_density"]) - catch - zeros(size(coordinates, 2)) - end - end - - pressure = try - get_data(point_data["pressure"]) - catch - zeros(size(coordinates, 2)) + get_data(point_data["hydrodynamic_density"]) end velocity = try @@ -37,7 +29,12 @@ function vtk2trixi(file; iter=0, input_directory="out", prefix="") try get_data(point_data["initial_velocity"]) catch - zeros(size(coordinates)) + try + get_data(point_data["wall_velocity"]) + catch + # Case is only used for 'boundary_systems' in simulations where velocity is not written out + zeros(size(coordinates)) + end end end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 82ad2ee7e..aa696579c 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -123,7 +123,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre # Store particle index vtk["index"] = active_particles(system) vtk["time"] = t - vtk["ndims"] = ndims(system) + vtk["ndmis"] = ndims(system) if write_meta_data vtk["solver_version"] = git_hash From e827bcace0acdef7f04c9fcc6923a57c4778fd4d Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 11:28:17 +0100 Subject: [PATCH 14/43] correct ndims --- src/visualization/write2vtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index aa696579c..82ad2ee7e 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -123,7 +123,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre # Store particle index vtk["index"] = active_particles(system) vtk["time"] = t - vtk["ndmis"] = ndims(system) + vtk["ndims"] = ndims(system) if write_meta_data vtk["solver_version"] = git_hash From a0e042e688c60e145cd2bb6c7efc9982c27bc9b0 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 14:35:44 +0100 Subject: [PATCH 15/43] stable v1 --- src/preprocessing/readvtk/vtk2trixi.jl | 26 +++++++------------------- src/visualization/write2vtk.jl | 4 +++- 2 files changed, 10 insertions(+), 20 deletions(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 1433beaec..9d8f67d52 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -8,14 +8,10 @@ function vtk2trixi(file; iter=0, input_directory="out", prefix="") # Retrieve point data fields (e.g., pressure, velocity, ...) point_data = get_point_data(vtk_file) - meta_data = get_field_data(vtk_file) - - # Meta data only written out in simulations - NDIMS = get_data(meta_data["ndims"]) - coordinates = get_points(vtk_file)[1:NDIMS[1], :] # Retrieve fields pressure = get_data(point_data["pressure"]) + # mass = get_data(point_data["mass"]) # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi density = try get_data(point_data["density"]) @@ -27,28 +23,20 @@ function vtk2trixi(file; iter=0, input_directory="out", prefix="") get_data(point_data["velocity"]) catch try - get_data(point_data["initial_velocity"]) + get_data(point_data["wall_velocity"]) catch - try - get_data(point_data["wall_velocity"]) - catch - # Case is only used for 'boundary_systems' in simulations where velocity is not written out - zeros(size(coordinates)) - end + get_data(point_data["initial_velocity"]) end end - # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - # mass = try - # get_data(point_data["mass"]) - # catch - # zeros(size(coordinates, 2)) - # end + # Coordinates are stored as a 3xN array + # Adjust the dimension of coordinates to match the velocity field + coordinates = get_points(vtk_file)[axes(velocity, 1), :] return InitialCondition( ; coordinates, velocity, - mass=ones(size(coordinates, 2)), + mass=zeros(size(coordinates, 2)), density, pressure) end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 82ad2ee7e..ec296b7d6 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -123,7 +123,6 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre # Store particle index vtk["index"] = active_particles(system) vtk["time"] = t - vtk["ndims"] = ndims(system) if write_meta_data vtk["solver_version"] = git_hash @@ -385,6 +384,9 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; if model.viscosity isa ViscosityAdami vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :) + else + vtk["velocity"] = [current_velocity(v, system, particle) + for particle in active_particles(system)] end end From 1dbec171c52007e858fa771b97b52fa6bd8bdf94 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 15:05:01 +0100 Subject: [PATCH 16/43] delete using get_field_data() --- src/TrixiParticles.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 67b6bfa0b..1b4ee8ded 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -19,7 +19,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch using Printf: @printf, @sprintf -using ReadVTK: VTKFile, get_data, get_field_data, get_point_data, get_points +using ReadVTK: VTKFile, get_data, get_point_data, get_points using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, From 642b417c38b362ff8fcad4e1006f36e2bca40df5 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 25 Nov 2024 15:44:43 +0100 Subject: [PATCH 17/43] update example file --- examples/readvtk/readvtk.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index cebcb3d55..7adc53b90 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,15 +1,11 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -#rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.0) -#trixi2vtk(rectangle; filename="rectangle", output_directory="out", -# custom_quantity=nothing) - -#filename = "rectangle" -#file = joinpath("out", filename * ".vtu") +rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.0) +trixi2vtk(rectangle; filename="rectangle", output_directory="out", + custom_quantity=nothing) -# The example file is from the simulation “dam_break_2d.jl” at time step 70 (1.4s) with added 'custom_quantities' -filename = "boundary_1_6" +filename = "rectangle" file = joinpath("out", filename * ".vtu") # Read the vtk file and convert it to an 'InitialCondition' From 83e5c3bdeb34c07fa9cdd5b64dfaab7d6491aeb2 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 1 Dec 2024 07:36:23 +0100 Subject: [PATCH 18/43] Bump crate-ci/typos from 1.26.8 to 1.28.1 (#678) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.26.8 to 1.28.1. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.26.8...v1.28.1) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 27a81105f..3d5579157 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.26.8 + uses: crate-ci/typos@v1.28.1 From a1cc9d611b26e353ee100323a3ac95ebb7966627 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 2 Dec 2024 10:19:31 +0100 Subject: [PATCH 19/43] Bump codecov/codecov-action from 4 to 5 (#679) Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v4...v5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 06c868258..cae43ee26 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -95,7 +95,7 @@ jobs: - name: Upload coverage report to Codecov # Only run coverage in one Job (Ubuntu and latest Julia version) if: matrix.os == 'ubuntu-latest' && matrix.version == '1' - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: files: lcov.info fail_ci_if_error: true From c93ab39970eacd67ca8e98a484db6dc814380473 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Mon, 2 Dec 2024 11:53:35 +0100 Subject: [PATCH 20/43] Plot recipes: support automatic marker size with `xlims` and `ylims` (#677) * add `xlims` and `ylmis` * implement suggestions --- src/visualization/recipes_plots.jl | 34 +++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index 53b588e34..a14de2948 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -10,7 +10,8 @@ RecipesBase.@recipe function f(sol::TrixiParticlesODESolution) end RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; - size=(600, 400)) # Default size + size=(600, 400), # Default size + xlims=(-Inf, Inf), ylims=(-Inf, Inf)) systems_data = map(semi.systems) do system u = wrap_u(u_ode, system, semi) coordinates = active_coordinates(u, system) @@ -25,6 +26,14 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; x_min, y_min = minimum(coordinates, dims=2) .- 0.5particle_spacing x_max, y_max = maximum(coordinates, dims=2) .+ 0.5particle_spacing + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, label=timer_name(system)) end @@ -32,7 +41,8 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; return (semi, systems_data...) end -RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) +RecipesBase.@recipe function f((initial_conditions::InitialCondition)...; + xlims=(Inf, Inf), ylims=(Inf, Inf)) idx = 0 ics = map(initial_conditions) do ic x = collect(ic.coordinates[1, :]) @@ -46,6 +56,14 @@ RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) x_min, y_min = minimum(ic.coordinates, dims=2) .- 0.5particle_spacing x_max, y_max = maximum(ic.coordinates, dims=2) .+ 0.5particle_spacing + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + idx += 1 return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, @@ -56,12 +74,22 @@ RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) end RecipesBase.@recipe function f(::Union{InitialCondition, Semidiscretization}, - data...; zcolor=nothing, size=(600, 400), colorbar_title="") + data...; zcolor=nothing, size=(600, 400), colorbar_title="", + xlims=(Inf, Inf), ylims=(Inf, Inf)) x_min = minimum(obj.x_min for obj in data) x_max = maximum(obj.x_max for obj in data) + y_min = minimum(obj.y_min for obj in data) y_max = maximum(obj.y_max for obj in data) + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + # Note that this assumes the plot area to be ~10% smaller than `size`, # which is the case when showing a single plot with the legend inside. # With the legend outside, this is no longer the case, so the `markersize` has to be From cda17e7779eec8fc0b52adf856c2b5ec8994d54e Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Tue, 3 Dec 2024 15:53:35 +0100 Subject: [PATCH 21/43] use Regular Expressions --- examples/readvtk/readvtk.jl | 5 ++- src/TrixiParticles.jl | 2 +- src/preprocessing/readvtk/vtk2trixi.jl | 56 ++++++++++++++++---------- src/visualization/write2vtk.jl | 4 +- 4 files changed, 40 insertions(+), 27 deletions(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 7adc53b90..763a0a0dd 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -1,7 +1,10 @@ # Example File to read a vtk file and convert it to InitialCondition using TrixiParticles -rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.0) +# Create a example rectangle +rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0), + pressure=1000.0) + trixi2vtk(rectangle; filename="rectangle", output_directory="out", custom_quantity=nothing) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 1b4ee8ded..67b6bfa0b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -19,7 +19,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch using Printf: @printf, @sprintf -using ReadVTK: VTKFile, get_data, get_point_data, get_points +using ReadVTK: VTKFile, get_data, get_field_data, get_point_data, get_points using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 9d8f67d52..829e0e336 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -1,37 +1,49 @@ -# TODO: Write documentation -# TODO: Write tests """ - Convert data from VTK-file to InitialCondition + vtk2trixi(file) + +Convert data from VTK-file to InitialCondition """ -function vtk2trixi(file; iter=0, input_directory="out", prefix="") +function vtk2trixi(file) vtk_file = VTKFile(file) # Retrieve point data fields (e.g., pressure, velocity, ...) point_data = get_point_data(vtk_file) # Retrieve fields - pressure = get_data(point_data["pressure"]) - # mass = get_data(point_data["mass"]) # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi + velocity = nothing + coordinates = nothing + try + velocity = get_data(point_data[first(k + for k in keys(point_data) + if occursin(r"(?i)velocity", k))]) - density = try - get_data(point_data["density"]) + # Coordinates are stored as a 3xN array + # Adjust the dimension of coordinates to match the velocity field + coordinates = get_points(vtk_file)[axes(velocity, 1), :] catch - get_data(point_data["hydrodynamic_density"]) - end + field_data = get_field_data(vtk_file) + ndims = get_data(field_data["ndims"]) - velocity = try - get_data(point_data["velocity"]) - catch - try - get_data(point_data["wall_velocity"]) - catch - get_data(point_data["initial_velocity"]) - end + # If no velocity field was found, adjust the coordinates to the dimensions + coordinates = get_points(vtk_file)[1:ndims[1], :] + + # Make sure that the 'InitialCondition' has a velocity field + @warn "No 'velocity' field found in VTK file. Velocity is set to zero." + velocity = zeros(size(coordinates)) end - # Coordinates are stored as a 3xN array - # Adjust the dimension of coordinates to match the velocity field - coordinates = get_points(vtk_file)[axes(velocity, 1), :] + density = get_data(point_data[first(k + for k in keys(point_data) + if occursin(r"(?i)density", k))]) + + pressure = get_data(point_data[first(k + for k in keys(point_data) + if occursin(r"(?i)pressure", k))]) + + # mass = get_data(point_data[first(k + # for k in keys(point_data) + # if occursin(r"(?i)mass", k))]) + # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi return InitialCondition( ; coordinates, @@ -39,4 +51,4 @@ function vtk2trixi(file; iter=0, input_directory="out", prefix="") mass=zeros(size(coordinates, 2)), density, pressure) -end +end \ No newline at end of file diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ec296b7d6..82ad2ee7e 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -123,6 +123,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre # Store particle index vtk["index"] = active_particles(system) vtk["time"] = t + vtk["ndims"] = ndims(system) if write_meta_data vtk["solver_version"] = git_hash @@ -384,9 +385,6 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; if model.viscosity isa ViscosityAdami vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :) - else - vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] end end From d3531245df2e3f715359bf5a2a0c86b92c7acae8 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 9 Dec 2024 13:51:46 +0100 Subject: [PATCH 22/43] tesset begin --- examples/fluid/hydrostatic_water_column_3d.jl | 6 +---- test/preprocessing/preprocessing.jl | 1 + test/preprocessing/readvtk/vtk2trixi.jl | 25 +++++++++++++++++++ 3 files changed, 27 insertions(+), 5 deletions(-) create mode 100644 test/preprocessing/readvtk/vtk2trixi.jl diff --git a/examples/fluid/hydrostatic_water_column_3d.jl b/examples/fluid/hydrostatic_water_column_3d.jl index 7168b0557..7fddfebc5 100644 --- a/examples/fluid/hydrostatic_water_column_3d.jl +++ b/examples/fluid/hydrostatic_water_column_3d.jl @@ -2,8 +2,4 @@ using TrixiParticles trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), - fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9), - tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81), - smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0), - maxiters=10^5, fluid_density_calculator=ContinuityDensity(), - clip_negative_pressure=false) + tspan=(0.0, 0.1)) diff --git a/test/preprocessing/preprocessing.jl b/test/preprocessing/preprocessing.jl index d2a994f93..223044b74 100644 --- a/test/preprocessing/preprocessing.jl +++ b/test/preprocessing/preprocessing.jl @@ -1 +1,2 @@ include("geometries/geometries.jl") +include("readvtk/readvtk.jl") diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl new file mode 100644 index 000000000..9b1919e08 --- /dev/null +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -0,0 +1,25 @@ +@testset verbose=true "vtk2trixi" begin + @testset verbose=true "Basic Functionality - Inital File" begin + min_coordinates = [(0, 0), (0, 0, 0)] + n_particles_per_dimension = [(10, 20), (10, 20, 30)] + velocity = [(1.0, 2.0), (1.0, 2.0, 3.0)] + + for i in 1:2 # 2d and 3d case + saved_ic = RectangularShape(0.1, n_particles_per_dimension[i], + min_coordinates[i], density=1.5, + velocity=velocity[i], pressure=1000.0) + filename = "is_write_out" + file = trixi2vtk(saved_ic; filename=filename) + + loaded_ic = vtk2trixi(joinpath("out", filename * ".vtu")) + + @test saved_ic.coordinates == loaded_ic.coordinates + @test saved_ic.velocity == loaded_ic.velocity + @test saved_ic.density == loaded_ic.density + @test saved_ic.pressure == loaded_ic.pressure + #@test saved_ic.mass == loaded_ic.mass + end + end + + @testset verbose=true "Basic Functionality - Simulation File" begin end +end \ No newline at end of file From f051b55922aceab34ebe436ecf90a8a0b9d5e6d6 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 11 Dec 2024 06:39:53 +0100 Subject: [PATCH 23/43] Improve Test Coverage (#674) * refactoring and improve coverage * add test * cleanup * format * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * fix --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- test/callbacks/postprocess.jl | 12 + .../dummy_particles/dummy_particles.jl | 273 +++++++++++------- 2 files changed, 173 insertions(+), 112 deletions(-) diff --git a/test/callbacks/postprocess.jl b/test/callbacks/postprocess.jl index 372a7b453..551392000 100644 --- a/test/callbacks/postprocess.jl +++ b/test/callbacks/postprocess.jl @@ -1,4 +1,16 @@ @testset verbose=true "PostprocessCallback" begin + @testset verbose=true "errors" begin + error_str1 = "`funcs` cannot be empty" + @test_throws ArgumentError(error_str1) PostprocessCallback(interval=10, + write_file_interval=0) + + error_str2 = "setting both `interval` and `dt` is not supported" + @test_throws ArgumentError(error_str2) PostprocessCallback(interval=10, + write_file_interval=0, + dt=0.1, + another_function=(v, u, t, system) -> 1) + end + @testset verbose=true "show" begin function example_function(v, u, t, system) return 0 diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index a72ac409c..dc8491cbb 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -1,142 +1,191 @@ @testset verbose=true "Dummy Particles" begin @testset "show" begin - boundary_model = BoundaryModelDummyParticles([1000.0], [1.0], + boundary_model = BoundaryModelDummyParticles([1000.0], + [1.0], SummationDensity(), - SchoenbergCubicSplineKernel{2}(), 0.1) + SchoenbergCubicSplineKernel{2}(), + 0.1) - show_compact = "BoundaryModelDummyParticles(SummationDensity, Nothing)" - @test repr(boundary_model) == show_compact + expected_repr = "BoundaryModelDummyParticles(SummationDensity, Nothing)" + @test repr(boundary_model) == expected_repr end - @testset "Viscosity Adami: Wall Velocity" begin + @testset "Viscosity Adami/Bernoulli: Wall Velocity" begin particle_spacing = 0.1 # Boundary particles in fluid compact support - boundary_1 = RectangularShape(particle_spacing, (10, 1), (0.0, 0.2), density=257.0) - boundary_2 = RectangularShape(particle_spacing, (10, 1), (0.0, 0.1), density=257.0) + boundary_1 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.2), + density=257.0) + boundary_2 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.1), + density=257.0) # Boundary particles out of fluid compact support - boundary_3 = RectangularShape(particle_spacing, (10, 1), (0, 0), density=257.0) + boundary_3 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.0), + density=257.0) boundary = union(boundary_1, boundary_2, boundary_3) particles_in_compact_support = length(boundary_1.mass) + length(boundary_2.mass) - fluid = RectangularShape(particle_spacing, (16, 5), (-0.3, 0.3), density=257.0, + # Define fluid particles + fluid = RectangularShape(particle_spacing, + (16, 5), + (-0.3, 0.3), + density=257.0, loop_order=:x_first) + # Simulation parameters smoothing_kernel = SchoenbergCubicSplineKernel{2}() smoothing_length = 1.2 * particle_spacing viscosity = ViscosityAdami(nu=1e-6) - state_equation = StateEquationCole(sound_speed=10, reference_density=257, + state_equation = StateEquationCole(sound_speed=10.0, + reference_density=257.0, exponent=7) - boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, - state_equation=state_equation, - AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length, - viscosity=viscosity) - - boundary_system = BoundarySPHSystem(boundary, boundary_model) - - fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), - state_equation, - smoothing_kernel, smoothing_length) - - neighborhood_search = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius=1.0, - eachpoint=TrixiParticles.eachparticle(fluid_system)) - - velocities = [[0; -1], [1; 1], [-1; 0], [0.7; 0.2], [0.3; 0.8]] - - @testset "Wall Velocity $v_fluid" for v_fluid in velocities - viscosity = boundary_system.boundary_model.viscosity - volume = boundary_system.boundary_model.cache.volume - - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - boundary_system.boundary_model.viscosity) - TrixiParticles.boundary_pressure_extrapolation!(boundary_model, - boundary_system, - fluid_system, - boundary.coordinates, - fluid.coordinates, v_fluid, - v_fluid .* - ones(size(fluid.coordinates)), - neighborhood_search) - - for particle in TrixiParticles.eachparticle(boundary_system) - if volume[particle] > eps() - TrixiParticles.compute_wall_velocity!(viscosity, boundary_system, - boundary.coordinates, particle) + # Define pressure extrapolation methods to test + pressure_extrapolations = [ + AdamiPressureExtrapolation(), + BernoulliPressureExtrapolation() + ] + + for pressure_extrapolation in pressure_extrapolations + @testset "Pressure Extrapolation: $(typeof(pressure_extrapolation))" begin + # Create boundary and fluid systems + boundary_model = BoundaryModelDummyParticles(boundary.density, + boundary.mass, + state_equation=state_equation, + pressure_extrapolation, + smoothing_kernel, + smoothing_length, + viscosity=viscosity) + boundary_system = BoundarySPHSystem(boundary, boundary_model) + fluid_system = WeaklyCompressibleSPHSystem(fluid, + SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length) + + neighborhood_search = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius=1.0, + eachpoint=TrixiParticles.eachparticle(fluid_system)) + + velocities = [ + [0.0; -1.0], + [1.0; 1.0], + [-1.0; 0.0], + [0.7; 0.2], + [0.3; 0.8] + ] + + @testset "Wall Velocity $v_fluid" for v_fluid in velocities + viscosity = boundary_system.boundary_model.viscosity + volume = boundary_system.boundary_model.cache.volume + + # Reset cache and perform pressure extrapolation + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + boundary_system.boundary_model.viscosity) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, + v_fluid, + v_fluid .* + ones(size(fluid.coordinates)), + neighborhood_search) + + # Compute wall velocities + for particle in TrixiParticles.eachparticle(boundary_system) + if volume[particle] > eps() + TrixiParticles.compute_wall_velocity!(viscosity, + boundary_system, + boundary.coordinates, + particle) + end + end + + # Expected wall velocities + v_wall = zeros(size(boundary.coordinates)) + v_wall[:, 1:particles_in_compact_support] .= -v_fluid + + @test isapprox(boundary_system.boundary_model.cache.wall_velocity, + v_wall) end - end - - v_wall = zeros(size(boundary.coordinates)) - v_wall[:, 1:particles_in_compact_support] .= -v_fluid - - @test isapprox(boundary_system.boundary_model.cache.wall_velocity, v_wall) - end - - scale_v = [1, 0.5, 0.7, 1.8, 67.5] - @testset "Wall Velocity Staggerd: Factor $scale" for scale in scale_v - viscosity = boundary_system.boundary_model.viscosity - volume = boundary_system.boundary_model.cache.volume - - # For a constant velocity profile (each fluid particle has the same velocity) - # the wall velocity is `v_wall = -v_fluid` (see eq. 22 in Adami_2012). - # Thus, generate a staggered velocity profile to test the smoothed velocity field - # for a variable velocity profile. - v_fluid = zeros(size(fluid.coordinates)) - for i in TrixiParticles.eachparticle(fluid_system) - if mod(i, 2) == 1 - v_fluid[:, i] .= scale - end - end - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - boundary_system.boundary_model.viscosity) - TrixiParticles.boundary_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, - boundary.coordinates, - fluid.coordinates, v_fluid, - v_fluid, - neighborhood_search) - - for particle in TrixiParticles.eachparticle(boundary_system) - if volume[particle] > eps() - TrixiParticles.compute_wall_velocity!(viscosity, boundary_system, - boundary.coordinates, particle) + scale_factors = [1.0, 0.5, 0.7, 1.8, 67.5] + + # For a constant velocity profile (each fluid particle has the same velocity), + # the wall velocity is `v_wall = -v_fluid` (see eq. 22 in Adami_2012). + # With a staggered velocity profile, we can test the smoothed velocity field + # for a variable velocity profile. + @testset "Wall Velocity Staggered: Factor $scale" for scale in scale_factors + viscosity = boundary_system.boundary_model.viscosity + volume = boundary_system.boundary_model.cache.volume + + # Create a staggered velocity profile + v_fluid = zeros(size(fluid.coordinates)) + for i in TrixiParticles.eachparticle(fluid_system) + if mod(i, 2) == 1 + v_fluid[:, i] .= scale + end + end + + # Reset cache and perform pressure extrapolation + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + boundary_system.boundary_model.viscosity) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, + v_fluid, + v_fluid, + neighborhood_search) + + # Compute wall velocities + for particle in TrixiParticles.eachparticle(boundary_system) + if volume[particle] > eps() + TrixiParticles.compute_wall_velocity!(viscosity, + boundary_system, + boundary.coordinates, + particle) + end + end + + # Expected wall velocities + v_wall = zeros(size(boundary.coordinates)) + + # First boundary row + for i in 1:length(boundary_1.mass) + if mod(i, 2) == 1 + # Particles with a diagonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.42040669416720744 * scale + else + # Particles with an orthogonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.5795933058327924 * scale + end + end + + # Second boundary row + for i in (length(boundary_1.mass) + 1):particles_in_compact_support + if mod(i, 2) == 1 + # Particles with a diagonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.12101100073462243 * scale + else + # Particles with an orthogonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.8789889992653775 * scale + end + end + + @test isapprox(boundary_system.boundary_model.cache.wall_velocity, + v_wall) end end - - v_wall = zeros(size(boundary.coordinates)) - - # First boundary row - for i in 1:length(boundary_1.mass) - if mod(i, 2) == 1 - - # Particles with a diagonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.42040669416720744 * scale - else - - # Particles with a orthogonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.5795933058327924 * scale - end - end - - # Second boundary row - for i in (length(boundary_1.mass) + 1):particles_in_compact_support - if true == mod(i, 2) - - # Particles with a diagonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.12101100073462243 * scale - else - - # Particles with a orthogonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.8789889992653775 * scale - end - end - - @test isapprox(boundary_system.boundary_model.cache.wall_velocity, v_wall) end end end From 91025f137806a9fe85d7c8dde02d46016a350256 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 11 Dec 2024 16:25:33 +0100 Subject: [PATCH 24/43] use match --- examples/readvtk/readvtk.jl | 2 +- src/preprocessing/readvtk/vtk2trixi.jl | 68 ++++++++++++-------------- src/visualization/write2vtk.jl | 2 +- 3 files changed, 32 insertions(+), 40 deletions(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 763a0a0dd..0ca097eea 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -8,7 +8,7 @@ rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, trixi2vtk(rectangle; filename="rectangle", output_directory="out", custom_quantity=nothing) -filename = "rectangle" +filename = "open_boundary_1_0" file = joinpath("out", filename * ".vtu") # Read the vtk file and convert it to an 'InitialCondition' diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 829e0e336..45fae688a 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -6,49 +6,41 @@ Convert data from VTK-file to InitialCondition function vtk2trixi(file) vtk_file = VTKFile(file) - # Retrieve point data fields (e.g., pressure, velocity, ...) + # Retrieve data fields (e.g., pressure, velocity, ...) point_data = get_point_data(vtk_file) + field_data = get_field_data(vtk_file) # Retrieve fields - velocity = nothing - coordinates = nothing - try - velocity = get_data(point_data[first(k - for k in keys(point_data) - if occursin(r"(?i)velocity", k))]) - - # Coordinates are stored as a 3xN array - # Adjust the dimension of coordinates to match the velocity field - coordinates = get_points(vtk_file)[axes(velocity, 1), :] - catch - field_data = get_field_data(vtk_file) - ndims = get_data(field_data["ndims"]) - - # If no velocity field was found, adjust the coordinates to the dimensions - coordinates = get_points(vtk_file)[1:ndims[1], :] - - # Make sure that the 'InitialCondition' has a velocity field - @warn "No 'velocity' field found in VTK file. Velocity is set to zero." - velocity = zeros(size(coordinates)) + ndims = get_data(field_data["ndims"]) + coordinates = get_points(vtk_file)[1:ndims[1], :] + + fields = ["velocity", "density", "pressure", "mass"] + results = Dict{String, Array{Float64}}() + + for field in fields + found = false + for k in keys(point_data) + if match(Regex("$field"), k) !== nothing + results[field] = get_data(point_data[k]) + found = true + break + end + end + if !found + # Set fields to zero if not found + if field in ["density", "pressure", "mass"] + results[field] = zeros(size(coordinates, 2)) + else + results[field] = zeros(size(coordinates)) + end + @info "No '$field' field found in VTK file. $field is set to zero." + end end - density = get_data(point_data[first(k - for k in keys(point_data) - if occursin(r"(?i)density", k))]) - - pressure = get_data(point_data[first(k - for k in keys(point_data) - if occursin(r"(?i)pressure", k))]) - - # mass = get_data(point_data[first(k - # for k in keys(point_data) - # if occursin(r"(?i)mass", k))]) - # TODO: Read out mass correctly as soon as mass is written out in vtu-file by Trixi - return InitialCondition( ; coordinates, - velocity, - mass=zeros(size(coordinates, 2)), - density, - pressure) + velocity=results["velocity"], + mass=results["mass"], + density=results["density"], + pressure=results["pressure"]) end \ No newline at end of file diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 82ad2ee7e..4029cf586 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -195,6 +195,7 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo vtk_grid(file, points, cells) do vtk # Store particle index. vtk["index"] = [i for i in axes(coordinates, 2)] + vtk["ndims"] = size(coordinates, 1) # Extract custom quantities for this system. for (key, quantity) in custom_quantities @@ -203,7 +204,6 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo end end end - return file end From 09c0a52d0eba6bc76c29f33628f0b9ba299a6932 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 11 Dec 2024 16:26:37 +0100 Subject: [PATCH 25/43] correct example file --- examples/readvtk/readvtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/readvtk/readvtk.jl b/examples/readvtk/readvtk.jl index 0ca097eea..763a0a0dd 100644 --- a/examples/readvtk/readvtk.jl +++ b/examples/readvtk/readvtk.jl @@ -8,7 +8,7 @@ rectangle = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, trixi2vtk(rectangle; filename="rectangle", output_directory="out", custom_quantity=nothing) -filename = "open_boundary_1_0" +filename = "rectangle" file = joinpath("out", filename * ".vtu") # Read the vtk file and convert it to an 'InitialCondition' From 3b008ee229f39638ae56cc583a1a23d5996cd66f Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Fri, 13 Dec 2024 14:58:26 +0100 Subject: [PATCH 26/43] Update src/preprocessing/readvtk/vtk2trixi.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 45fae688a..a8e7462a4 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -1,5 +1,5 @@ """ - vtk2trixi(file) + vtk2trixi(file::String) Convert data from VTK-file to InitialCondition """ From 46c8beae5bf58a08a23075f7e80351756490a863 Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Fri, 13 Dec 2024 14:59:05 +0100 Subject: [PATCH 27/43] Update src/preprocessing/readvtk/vtk2trixi.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/preprocessing/readvtk/vtk2trixi.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index a8e7462a4..46d4c1e91 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -2,6 +2,9 @@ vtk2trixi(file::String) Convert data from VTK-file to InitialCondition + +# Arguments +- `file`: Name of the file to be loaded. """ function vtk2trixi(file) vtk_file = VTKFile(file) From 716f488c342ea15ac15354f2abdb926fb1bf5843 Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Fri, 13 Dec 2024 14:59:35 +0100 Subject: [PATCH 28/43] Update src/preprocessing/readvtk/vtk2trixi.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 46d4c1e91..dcfcb6840 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -23,7 +23,7 @@ function vtk2trixi(file) for field in fields found = false for k in keys(point_data) - if match(Regex("$field"), k) !== nothing + if !isnothing(match(Regex("$field"), k)) results[field] = get_data(point_data[k]) found = true break From 599c97c5d1de690da8edf9e8010a64ad75703717 Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Fri, 13 Dec 2024 15:00:10 +0100 Subject: [PATCH 29/43] Update src/preprocessing/readvtk/vtk2trixi.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index dcfcb6840..fca9ce21c 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -34,7 +34,7 @@ function vtk2trixi(file) if field in ["density", "pressure", "mass"] results[field] = zeros(size(coordinates, 2)) else - results[field] = zeros(size(coordinates)) + results[field] = zero(coordinates) end @info "No '$field' field found in VTK file. $field is set to zero." end From eae739573d587c05042d0cf38b92452a9712397a Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Fri, 13 Dec 2024 15:02:42 +0100 Subject: [PATCH 30/43] info update --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index fca9ce21c..8c0db60cd 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -36,7 +36,7 @@ function vtk2trixi(file) else results[field] = zero(coordinates) end - @info "No '$field' field found in VTK file. $field is set to zero." + @info "No '$field' field found in VTK file. Is set to zero." end end From 9c1a50a4c8b103be824c040d9504ee5d4b1b33f5 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Fri, 13 Dec 2024 15:15:52 +0100 Subject: [PATCH 31/43] update info --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 8c0db60cd..580f72421 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -36,7 +36,7 @@ function vtk2trixi(file) else results[field] = zero(coordinates) end - @info "No '$field' field found in VTK file. Is set to zero." + @info "No '$field' field found in VTK file. Will be set to zero." end end From 61e88fe67fdaed3d46f79e7bbb84f88db6903ca9 Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Mon, 16 Dec 2024 11:03:51 +0100 Subject: [PATCH 32/43] Update hydrostatic_water_column_3d.jl --- examples/fluid/hydrostatic_water_column_3d.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/fluid/hydrostatic_water_column_3d.jl b/examples/fluid/hydrostatic_water_column_3d.jl index 7fddfebc5..de918b563 100644 --- a/examples/fluid/hydrostatic_water_column_3d.jl +++ b/examples/fluid/hydrostatic_water_column_3d.jl @@ -1,5 +1,9 @@ using TrixiParticles trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), + joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9), + tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81), + smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0), + maxiters=10^5, fluid_density_calculator=ContinuityDensity(), + clip_negative_pressure=false) tspan=(0.0, 0.1)) From e1a59725d7b262a98d2c87107a5a4d9396b7312c Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 16 Dec 2024 11:05:27 +0100 Subject: [PATCH 33/43] changes --- test/preprocessing/readvtk/vtk2trixi.jl | 44 ++++++++++++++----------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl index 9b1919e08..8552e2762 100644 --- a/test/preprocessing/readvtk/vtk2trixi.jl +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -1,25 +1,31 @@ @testset verbose=true "vtk2trixi" begin - @testset verbose=true "Basic Functionality - Inital File" begin - min_coordinates = [(0, 0), (0, 0, 0)] - n_particles_per_dimension = [(10, 20), (10, 20, 30)] - velocity = [(1.0, 2.0), (1.0, 2.0, 3.0)] + @testset verbose=true "Functionality Check - 2D" begin - for i in 1:2 # 2d and 3d case - saved_ic = RectangularShape(0.1, n_particles_per_dimension[i], - min_coordinates[i], density=1.5, - velocity=velocity[i], pressure=1000.0) - filename = "is_write_out" - file = trixi2vtk(saved_ic; filename=filename) + # 'InitialCondition'-Files + saved_ic = RectangularShape(0.1, (10, 20), + (0, 0), density=1.5, + velocity=(1.0, 2.0), pressure=1000.0) + filename = "is_write_out" + file = trixi2vtk(saved_ic; filename=filename) - loaded_ic = vtk2trixi(joinpath("out", filename * ".vtu")) + loaded_ic = vtk2trixi(joinpath("out", filename * ".vtu")) - @test saved_ic.coordinates == loaded_ic.coordinates - @test saved_ic.velocity == loaded_ic.velocity - @test saved_ic.density == loaded_ic.density - @test saved_ic.pressure == loaded_ic.pressure - #@test saved_ic.mass == loaded_ic.mass - end - end + @test saved_ic.coordinates == loaded_ic.coordinates + @test saved_ic.velocity == loaded_ic.velocity + @test saved_ic.density == loaded_ic.density + @test saved_ic.pressure == loaded_ic.pressure + @test saved_ic.mass == loaded_ic.mass + + # 'Fluidsystem'-Files + state_equation = StateEquationCole(; sound_speed=1.0, reference_density=1.0, + exponent=1, clip_negative_pressure=false) + fluid_system = WeaklyCompressibleSPHSystem(saved_ic, + ContinuityDensity(), state_equation, + WendlandC2Kernel{2}(), 0.2) - @testset verbose=true "Basic Functionality - Simulation File" begin end + trixi2vtk(zeros(), zeros(), 0.0, fluid_system, nothing;) + # trixi2vtk(sol.u[end], semi, 0.0, iter=1, output_directory="output", + # prefix="solution") + + end end \ No newline at end of file From a51cf2ffad59faef70448c24fdafb40484f6a228 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 16 Dec 2024 11:07:54 +0100 Subject: [PATCH 34/43] repair --- examples/fluid/hydrostatic_water_column_3d.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/fluid/hydrostatic_water_column_3d.jl b/examples/fluid/hydrostatic_water_column_3d.jl index de918b563..d8b53494f 100644 --- a/examples/fluid/hydrostatic_water_column_3d.jl +++ b/examples/fluid/hydrostatic_water_column_3d.jl @@ -1,9 +1,10 @@ using TrixiParticles trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9), + joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), + fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9), tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81), smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0), maxiters=10^5, fluid_density_calculator=ContinuityDensity(), - clip_negative_pressure=false) - tspan=(0.0, 0.1)) + clip_negative_pressure=false, + tspan=(0.0, 0.1)) \ No newline at end of file From 3ba15d8a0f9d9bb18fd8571eb65cfd7f4be9ceb8 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 16 Dec 2024 11:12:07 +0100 Subject: [PATCH 35/43] tests --- test/preprocessing/readvtk/vtk2trixi.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl index 8552e2762..7dc1d922b 100644 --- a/test/preprocessing/readvtk/vtk2trixi.jl +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -2,7 +2,7 @@ @testset verbose=true "Functionality Check - 2D" begin # 'InitialCondition'-Files - saved_ic = RectangularShape(0.1, (10, 20), + saved_ic = RectangularShape(0.1, (2, 2), (0, 0), density=1.5, velocity=(1.0, 2.0), pressure=1000.0) filename = "is_write_out" @@ -23,7 +23,7 @@ ContinuityDensity(), state_equation, WendlandC2Kernel{2}(), 0.2) - trixi2vtk(zeros(), zeros(), 0.0, fluid_system, nothing;) + trixi2vtk(saved_ic.velocity, saved_ic.coordinates, 0.0, fluid_system, nothing;) # trixi2vtk(sol.u[end], semi, 0.0, iter=1, output_directory="output", # prefix="solution") From 85792cda77d0e83ef408f79f3391909e44cdcf59 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Mon, 16 Dec 2024 11:16:14 +0100 Subject: [PATCH 36/43] tests --- examples/fluid/hydrostatic_water_column_3d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/fluid/hydrostatic_water_column_3d.jl b/examples/fluid/hydrostatic_water_column_3d.jl index d8b53494f..7168b0557 100644 --- a/examples/fluid/hydrostatic_water_column_3d.jl +++ b/examples/fluid/hydrostatic_water_column_3d.jl @@ -6,5 +6,4 @@ trixi_include(@__MODULE__, tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81), smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0), maxiters=10^5, fluid_density_calculator=ContinuityDensity(), - clip_negative_pressure=false, - tspan=(0.0, 0.1)) \ No newline at end of file + clip_negative_pressure=false) From b31753f7bcaf725cb8a6a8bcdec867a94ab995cb Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Wed, 18 Dec 2024 15:30:53 +0100 Subject: [PATCH 37/43] Update src/TrixiParticles.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/TrixiParticles.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 67b6bfa0b..c4017fccf 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -19,7 +19,7 @@ using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch using Printf: @printf, @sprintf -using ReadVTK: VTKFile, get_data, get_field_data, get_point_data, get_points +using ReadVTK: ReadVTK using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, From 5b272ad593342537f8181e5b6b4f18471b6cd581 Mon Sep 17 00:00:00 2001 From: Marcel Schurer <153067321+marcelschurer@users.noreply.github.com> Date: Wed, 18 Dec 2024 15:31:03 +0100 Subject: [PATCH 38/43] Update src/preprocessing/readvtk/vtk2trixi.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/preprocessing/readvtk/vtk2trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 580f72421..1e4f90165 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -10,7 +10,7 @@ function vtk2trixi(file) vtk_file = VTKFile(file) # Retrieve data fields (e.g., pressure, velocity, ...) - point_data = get_point_data(vtk_file) + point_data = ReadVTK.get_point_data(vtk_file) field_data = get_field_data(vtk_file) # Retrieve fields From c59f3705e59eed8a0c87cfebb6d1b83cb33f3e56 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Thu, 19 Dec 2024 10:33:02 +0100 Subject: [PATCH 39/43] implement suggestions --- test/preprocessing/readvtk/vtk2trixi.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl index 7dc1d922b..50c251fdf 100644 --- a/test/preprocessing/readvtk/vtk2trixi.jl +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -2,19 +2,21 @@ @testset verbose=true "Functionality Check - 2D" begin # 'InitialCondition'-Files - saved_ic = RectangularShape(0.1, (2, 2), + saved_ic = RectangularShape(0.1, (10, 20), (0, 0), density=1.5, velocity=(1.0, 2.0), pressure=1000.0) - filename = "is_write_out" - file = trixi2vtk(saved_ic; filename=filename) - loaded_ic = vtk2trixi(joinpath("out", filename * ".vtu")) + file = trixi2vtk(saved_ic; output_directory="test/preprocessing/readvtk", + filename="initial_condition") - @test saved_ic.coordinates == loaded_ic.coordinates - @test saved_ic.velocity == loaded_ic.velocity - @test saved_ic.density == loaded_ic.density - @test saved_ic.pressure == loaded_ic.pressure - @test saved_ic.mass == loaded_ic.mass + loaded_ic = vtk2trixi(joinpath("test/preprocessing/readvtk", + "initial_condition.vtu")) + + @test isapprox(saved_ic.coordinates, loaded_ic.coordinates) + @test isapprox(saved_ic.velocity, loaded_ic.velocity) + @test isapprox(saved_ic.density, loaded_ic.density) + @test isapprox(saved_ic.pressure, loaded_ic.pressure) + @test isapprox(saved_ic.mass, loaded_ic.mass) # 'Fluidsystem'-Files state_equation = StateEquationCole(; sound_speed=1.0, reference_density=1.0, From 3c0b0057ea2249e7f7880ef8b93d494db2e2ac30 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Thu, 19 Dec 2024 10:34:11 +0100 Subject: [PATCH 40/43] read out simulation file --- test/preprocessing/readvtk/vtk2trixi.jl | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl index 50c251fdf..532b936ed 100644 --- a/test/preprocessing/readvtk/vtk2trixi.jl +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -22,12 +22,19 @@ state_equation = StateEquationCole(; sound_speed=1.0, reference_density=1.0, exponent=1, clip_negative_pressure=false) fluid_system = WeaklyCompressibleSPHSystem(saved_ic, - ContinuityDensity(), state_equation, + SummationDensity(), state_equation, WendlandC2Kernel{2}(), 0.2) - trixi2vtk(saved_ic.velocity, saved_ic.coordinates, 0.0, fluid_system, nothing;) - # trixi2vtk(sol.u[end], semi, 0.0, iter=1, output_directory="output", - # prefix="solution") + saved_fs = trixi2vtk(saved_ic.velocity, saved_ic.coordinates, 0.0, fluid_system, + nothing; + output_directory="test/preprocessing/readvtk", iter=1) + loaded_fs = vtk2trixi(joinpath("test/preprocessing/readvtk", "fluid_1.vtu")) + + @test isapprox(saved_fs.coordinates, loaded_fs.coordinates) + @test isapprox(saved_fs.velocity, loaded_fs.velocity) + @test isapprox(saved_fs.density, loaded_fs.density) + @test isapprox(saved_fs.pressure, loaded_fs.pressure) + @test isapprox(saved_fs.mass, loaded_fs.mass) end end \ No newline at end of file From 5bf34d437b3bad9b30e3baa4598b74f1434b38b2 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Thu, 19 Dec 2024 10:34:48 +0100 Subject: [PATCH 41/43] implement suggestions --- src/preprocessing/readvtk/vtk2trixi.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 1e4f90165..6ea697a12 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -7,24 +7,24 @@ Convert data from VTK-file to InitialCondition - `file`: Name of the file to be loaded. """ function vtk2trixi(file) - vtk_file = VTKFile(file) + vtk_file = ReadVTK.VTKFile(file) # Retrieve data fields (e.g., pressure, velocity, ...) point_data = ReadVTK.get_point_data(vtk_file) - field_data = get_field_data(vtk_file) + field_data = ReadVTK.get_field_data(vtk_file) # Retrieve fields - ndims = get_data(field_data["ndims"]) - coordinates = get_points(vtk_file)[1:ndims[1], :] + ndims = ReadVTK.get_data(field_data["ndims"]) + coordinates = ReadVTK.get_points(vtk_file)[1:ndims[1], :] fields = ["velocity", "density", "pressure", "mass"] results = Dict{String, Array{Float64}}() for field in fields found = false - for k in keys(point_data) + for k in ReadVTK.keys(point_data) if !isnothing(match(Regex("$field"), k)) - results[field] = get_data(point_data[k]) + results[field] = ReadVTK.get_data(point_data[k]) found = true break end From 1363075e85ea40ffd4edd4edffbe5f20bb0a1829 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 1 Jan 2025 11:53:39 +0100 Subject: [PATCH 42/43] Reading out 'particle_spacing' --- src/preprocessing/readvtk/vtk2trixi.jl | 4 +++- src/visualization/write2vtk.jl | 4 ++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/preprocessing/readvtk/vtk2trixi.jl b/src/preprocessing/readvtk/vtk2trixi.jl index 6ea697a12..4a2c2afa1 100644 --- a/src/preprocessing/readvtk/vtk2trixi.jl +++ b/src/preprocessing/readvtk/vtk2trixi.jl @@ -15,6 +15,8 @@ function vtk2trixi(file) # Retrieve fields ndims = ReadVTK.get_data(field_data["ndims"]) + particle_spacing = first(ReadVTK.get_data(field_data["particle_spacing"])) + coordinates = ReadVTK.get_points(vtk_file)[1:ndims[1], :] fields = ["velocity", "density", "pressure", "mass"] @@ -41,7 +43,7 @@ function vtk2trixi(file) end return InitialCondition( - ; coordinates, + ; coordinates, particle_spacing, velocity=results["velocity"], mass=results["mass"], density=results["density"], diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 4029cf586..230668a47 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -124,6 +124,7 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre vtk["index"] = active_particles(system) vtk["time"] = t vtk["ndims"] = ndims(system) + vtk["particle_spacing"] = system.initial_condition.particle_spacing if write_meta_data vtk["solver_version"] = git_hash @@ -184,6 +185,7 @@ Convert coordinate data to VTK format. - `file::AbstractString`: Path to the generated VTK file. """ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coordinates", + particle_spacing=-1, custom_quantities...) mkpath(output_directory) file = prefix === "" ? joinpath(output_directory, filename) : @@ -196,6 +198,7 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo # Store particle index. vtk["index"] = [i for i in axes(coordinates, 2)] vtk["ndims"] = size(coordinates, 1) + vtk["particle_spacing"] = particle_spacing # Extract custom quantities for this system. for (key, quantity) in custom_quantities @@ -231,6 +234,7 @@ function trixi2vtk(initial_condition::InitialCondition; output_directory="out", return trixi2vtk(coordinates; output_directory, prefix, filename, density=density, initial_velocity=velocity, mass=mass, + particle_spacing=initial_condition.particle_spacing, pressure=pressure, custom_quantities...) end From d56757324295c82440472ccee0bb1e111ef98d67 Mon Sep 17 00:00:00 2001 From: marcelschurer Date: Wed, 1 Jan 2025 15:10:02 +0100 Subject: [PATCH 43/43] add test for 'fluid_system' and 'boundary_system' --- test/preprocessing/readvtk/vtk2trixi.jl | 131 +++++++++++++++++++----- 1 file changed, 104 insertions(+), 27 deletions(-) diff --git a/test/preprocessing/readvtk/vtk2trixi.jl b/test/preprocessing/readvtk/vtk2trixi.jl index 532b936ed..1c5d7b5a7 100644 --- a/test/preprocessing/readvtk/vtk2trixi.jl +++ b/test/preprocessing/readvtk/vtk2trixi.jl @@ -2,9 +2,7 @@ @testset verbose=true "Functionality Check - 2D" begin # 'InitialCondition'-Files - saved_ic = RectangularShape(0.1, (10, 20), - (0, 0), density=1.5, - velocity=(1.0, 2.0), pressure=1000.0) + saved_ic = rectangular_patch(0.1, (10, 20)) file = trixi2vtk(saved_ic; output_directory="test/preprocessing/readvtk", filename="initial_condition") @@ -12,29 +10,108 @@ loaded_ic = vtk2trixi(joinpath("test/preprocessing/readvtk", "initial_condition.vtu")) - @test isapprox(saved_ic.coordinates, loaded_ic.coordinates) - @test isapprox(saved_ic.velocity, loaded_ic.velocity) - @test isapprox(saved_ic.density, loaded_ic.density) - @test isapprox(saved_ic.pressure, loaded_ic.pressure) - @test isapprox(saved_ic.mass, loaded_ic.mass) - - # 'Fluidsystem'-Files - state_equation = StateEquationCole(; sound_speed=1.0, reference_density=1.0, - exponent=1, clip_negative_pressure=false) - fluid_system = WeaklyCompressibleSPHSystem(saved_ic, - SummationDensity(), state_equation, - WendlandC2Kernel{2}(), 0.2) - - saved_fs = trixi2vtk(saved_ic.velocity, saved_ic.coordinates, 0.0, fluid_system, - nothing; - output_directory="test/preprocessing/readvtk", iter=1) - - loaded_fs = vtk2trixi(joinpath("test/preprocessing/readvtk", "fluid_1.vtu")) - - @test isapprox(saved_fs.coordinates, loaded_fs.coordinates) - @test isapprox(saved_fs.velocity, loaded_fs.velocity) - @test isapprox(saved_fs.density, loaded_fs.density) - @test isapprox(saved_fs.pressure, loaded_fs.pressure) - @test isapprox(saved_fs.mass, loaded_fs.mass) + @test isapprox(saved_ic.coordinates, loaded_ic.coordinates, rtol=1e-5) + @test isapprox(saved_ic.velocity, loaded_ic.velocity, rtol=1e-5) + @test isapprox(saved_ic.density, loaded_ic.density, rtol=1e-5) + @test isapprox(saved_ic.pressure, loaded_ic.pressure, rtol=1e-5) + #@test isapprox(saved_ic.mass, loaded_ic.mass, rtol=1e-5) #TODO: wait until mass is written out with 'write2vtk' + + # Simulations-Files + particle_spacing = 0.05 + smoothing_length = 0.15 + boundary_layers = 3 + open_boundary_layers = 6 + + density = 1000.0 + pressure = 1.0 + domain_size = (1.0, 0.4) + + flow_direction = [1.0, 0.0] + reynolds_number = 100 + prescribed_velocity = 2.0 + sound_speed = 20.0 + + boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2]) + + state_equation = StateEquationCole(; sound_speed=sound_speed, + reference_density=density, + exponent=7, background_pressure=pressure) + + sim_ic = RectangularTank(particle_spacing, domain_size, boundary_size, density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true)) + + sim_ic.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers + + # ==== Fluid System + smoothing_kernel = WendlandC2Kernel{2}() + fluid_density_calculator = SummationDensity() + + kinematic_viscosity = 0.008 + n_buffer_particles = 0 + + viscosity = ViscosityAdami(nu=kinematic_viscosity) + + fluid_system = EntropicallyDampedSPHSystem(sim_ic.fluid, smoothing_kernel, + smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + + fluid_system.cache.density .= sim_ic.fluid.density + # Write out 'fluid_system' Simulation-File + trixi2vtk(sim_ic.fluid.velocity, sim_ic.fluid.coordinates, 0.0, fluid_system, + nothing; + output_directory="test/preprocessing/readvtk", iter=1) + + # Load 'fluid_system' Simulation-File + loaded_fs = vtk2trixi(joinpath("test/preprocessing/readvtk", + "fluid_1.vtu")) + + @test isapprox(sim_ic.fluid.coordinates, loaded_fs.coordinates, rtol=1e-5) + @test isapprox(sim_ic.fluid.velocity, loaded_fs.velocity, rtol=1e-5) + @test isapprox(sim_ic.fluid.density, loaded_fs.density, rtol=1e-5) + @test isapprox(sim_ic.fluid.pressure, loaded_fs.pressure, rtol=1e-5) + #@test isapprox(sim_ic.fluid.mass, loaded_fs.mass, rtol=1e-5) #TODO: wait until mass is written out with 'write2vtk' + + # ==== Boundary System + boundary_model = BoundaryModelDummyParticles(sim_ic.boundary.density, + sim_ic.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + smoothing_kernel, smoothing_length) + + boundary_system = BoundarySPHSystem(sim_ic.boundary, boundary_model) + + trixi2vtk(sim_ic.boundary.velocity, sim_ic.boundary.coordinates, 0.0, + boundary_system, + nothing; + output_directory="test/preprocessing/readvtk", iter=1) + + # Load 'boundary_system' Simulation-File + loaded_bs = vtk2trixi(joinpath("test/preprocessing/readvtk", + "boundary_1.vtu")) + + @test isapprox(sim_ic.boundary.coordinates, loaded_bs.coordinates, rtol=1e-5) + @test isapprox(sim_ic.boundary.velocity, loaded_bs.velocity, rtol=1e-5) + @test isapprox(sim_ic.boundary.density, loaded_bs.density, rtol=1e-5) + @test isapprox(sim_ic.boundary.pressure, loaded_bs.pressure, rtol=1e-5) + #@test isapprox(sim_ic.boundary.mass, loaded_bs.mass, rtol=1e-5) #TODO: wait until mass is written out with 'write2vtk' + + # ==== Open Boundary System + plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) + inflow = InFlow(; plane=plane_in, flow_direction, + open_boundary_layers, density=density, particle_spacing) + + function velocity_function2d(pos, t) + return SVector(prescribed_velocity, 0.0) + end + open_boundary = OpenBoundarySPHSystem(inflow; fluid_system, + boundary_model=BoundaryModelLastiwka(), + buffer_size=n_buffer_particles, + reference_density=density, + reference_pressure=pressure, + reference_velocity=velocity_function2d) end end \ No newline at end of file