diff --git a/Manifest.toml b/Manifest.toml index 79f86bf..f349bb7 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,6 +2,110 @@ julia_version = "1.9.4" manifest_format = "2.0" -project_hash = "d61a61014ab7dc6e0a8e034fba489fc14f7fd619" +project_hash = "a52433223847df22e7a48359808fd13703c14ae3" -[deps] +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra"] +git-tree-sha1 = "71acdbf594aab5bbb2cec89b208c41b4c411e49f" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.24.0" + + [deps.ChainRulesCore.extensions] + ChainRulesCoreSparseArraysExt = "SparseArrays" + + [deps.ChainRulesCore.weakdeps] + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.Compat]] +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.15.0" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.5+0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.21+4" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.ShockwaveProperties]] +deps = ["ChainRulesCore", "LinearAlgebra", "Unitful", "UnitfulChainRules"] +git-tree-sha1 = "b5c1b1cc410447f01e3a16ccd174af56a39d2310" +repo-rev = "master" +repo-url = "https://github.com/STCE-at-RWTH/ShockwaveProperties.jl" +uuid = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50" +version = "0.1.11" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.Unitful]] +deps = ["Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "dd260903fdabea27d9b6021689b3cd5401a57748" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.20.0" + + [deps.Unitful.extensions] + ConstructionBaseUnitfulExt = "ConstructionBase" + InverseFunctionsUnitfulExt = "InverseFunctions" + + [deps.Unitful.weakdeps] + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.UnitfulChainRules]] +deps = ["ChainRulesCore", "LinearAlgebra", "Unitful"] +git-tree-sha1 = "24a349fa6f4fbcdc12830c5640baa6eae532f83f" +uuid = "f31437dd-25a7-4345-875f-756556e6935d" +version = "0.1.2" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.8.0+0" diff --git a/Project.toml b/Project.toml index 44def77..5be1f28 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,9 @@ uuid = "a75b99cc-be56-4638-99bc-d89fa43b9ca1" authors = ["Warisa and contributors"] version = "1.0.0-DEV" +[deps] +ShockwaveProperties = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50" + [compat] julia = "1.9" diff --git a/examples/1D_shock_scenario1_detection.jl b/examples/1D_shock_scenario1_detection.jl new file mode 100644 index 0000000..ba2f2e0 --- /dev/null +++ b/examples/1D_shock_scenario1_detection.jl @@ -0,0 +1 @@ +using ShockwaveDetection \ No newline at end of file diff --git a/examples/1D_shock_scenario2_detection.jl b/examples/1D_shock_scenario2_detection.jl new file mode 100644 index 0000000..ba2f2e0 --- /dev/null +++ b/examples/1D_shock_scenario2_detection.jl @@ -0,0 +1 @@ +using ShockwaveDetection \ No newline at end of file diff --git a/examples/data/euler_scenario_1.out b/examples/data/euler_scenario_1.out new file mode 100644 index 0000000..6ba030d Binary files /dev/null and b/examples/data/euler_scenario_1.out differ diff --git a/examples/data/euler_scenario_1.tape b/examples/data/euler_scenario_1.tape new file mode 100644 index 0000000..fd05074 Binary files /dev/null and b/examples/data/euler_scenario_1.tape differ diff --git a/examples/data/euler_scenario_2.out b/examples/data/euler_scenario_2.out new file mode 100644 index 0000000..ad92493 Binary files /dev/null and b/examples/data/euler_scenario_2.out differ diff --git a/examples/data/euler_scenario_2.tape b/examples/data/euler_scenario_2.tape new file mode 100644 index 0000000..8881368 Binary files /dev/null and b/examples/data/euler_scenario_2.tape differ diff --git a/src/detect_shock.jl b/src/detect_shock.jl new file mode 100644 index 0000000..81ad868 --- /dev/null +++ b/src/detect_shock.jl @@ -0,0 +1,77 @@ +include("visualize.jl") + +# Use Interpolations.jl to perform np.gradient-like operations +using Interpolations + +using StatsPlots + +# Function to compute gradients using Interpolations.jl +# without explicit function definition +function compute_gradients(arr::AbstractVector, x) + # Convert StepRangeLen to LinRange because that's what Interpolations.jl expects + x_linrange = LinRange(first(x), last(x), length(x)) + itp = interpolate((x_linrange,), arr, Gridded(Linear())) + grad = only.(Interpolations.gradient.(Ref(itp), x_linrange)) + return grad +end + +# Function to compute zero crossings, detecting points where the second derivative changes sign\ +# suggesting a potential shock +function zero_crossings(arr) + return [i for i in 2:length(arr) if arr[i-1] * arr[i] < 0] +end + +# Function to detect shock locations at a given time step +function detect_shocks_at_timestep(density_at_t, velocity_at_t, pressure_at_t, x) + # Compute first gradients + density_grad = compute_gradients(density_at_t, x) + velocity_grad = compute_gradients(velocity_at_t, x) + pressure_grad = compute_gradients(pressure_at_t, x) + + # Find points where the gradient exceeds a certain threshold + threshold = 0.5 # Adjust this threshold as needed + shock_location_density = findall(gradient -> abs(gradient) > threshold, density_grad) + shock_location_velocity = findall(gradient -> abs(gradient) > threshold, velocity_grad) + shock_location_pressure = findall(gradient -> abs(gradient) > threshold, pressure_grad) + + # Combine detections (common shock location across variables) + shock_locations = intersect(intersect(shock_location_density, shock_location_velocity), shock_location_pressure) + + return shock_locations +end + +# Read the data +x0_xmax, t_values, u_values, dims_u = read_output_file("C:/Users/user/Documents/School/Sem4/softwareentwicklungspraktikum/shock_wave_detection/ShockwaveProperties.jl/example/data/euler_scenario_2.out") +x = range(x0_xmax[1], stop=x0_xmax[2], length=dims_u[2]) +density_field, velocity_field, pressure_field = convert_to_primitive(u_values) + +# Detect shocks for each time step and store the positions +shock_positions_over_time = [] + +for t_step in 1: length(t_values) + density_field_t = density_field[:, t_step] + velocity_field_t = velocity_field[:, t_step] + pressure_field_t = pressure_field[:, t_step] + shock_positions = detect_shocks_at_timestep(density_field_t,velocity_field_t,pressure_field_t, x) + push!(shock_positions_over_time, shock_positions) +end + +# Create an animation +anim = @animate for (t_step, t) in enumerate(t_values) + p1 = plot(x, density_field[:, t_step], title="Density at Time $t", xlabel="x", ylabel="Density", label = "Density across x", size=(800, 600)) + p2 = plot(x, velocity_field[:, t_step], title="Velocity at Time $t", xlabel="x", ylabel="Velocity", label = "Velocity across x", size=(800, 600)) + p3 = plot(x, pressure_field[:, t_step], title="Pressure at Time $t", xlabel="x", ylabel="Pressure", label = "Pressure across x", size=(800, 600)) + + # Add markers for the shock positions + shock_positions_t = shock_positions_over_time[t_step] + for pos in shock_positions_t + scatter!(p1, [x[pos]], [density_field[pos, t_step]], color=:red, label=false) + scatter!(p2, [x[pos]], [velocity_field[pos, t_step]], color=:red, label=false) + scatter!(p3, [x[pos]], [pressure_field[pos, t_step]], color=:red, label=false) + end + + plot(p1, p2, p3, layout = (3, 1)) +end + +# Save the animation as a gif +gif(anim, "shock_over_time.gif", fps = 10) \ No newline at end of file diff --git a/src/visualize.jl b/src/visualize.jl new file mode 100644 index 0000000..36d8b74 --- /dev/null +++ b/src/visualize.jl @@ -0,0 +1,86 @@ +using Plots +using FFMPEG +using ShockwaveProperties: primitive_state_vector, pressure, speed_of_sound, DRY_AIR +using Unitful: Pa, ustrip + +function read_output_file(filename) + open(filename, "r") do f + dims_u = Vector{Int}(undef, 2) + read!(f, dims_u) + + # Read the first and last x values + x0_xmax = Vector{Float64}(undef, 2) + read!(f, x0_xmax) + + # Read the number of time steps + num_timesteps = Vector{Int}(undef, 1) + read!(f, num_timesteps) + + # Read the time values + t_values = Vector{Float64}(undef, num_timesteps[1]) + read!(f, t_values) + + # Read the u values + u_values = Vector{Float64}(undef, prod(dims_u)*num_timesteps[1]) + read!(f, u_values) + + + # Reshape u_values to a 3D array + u_values = reshape(u_values, dims_u[1], dims_u[2], num_timesteps[1]) # N_u x N_x x N_t as u a vector of 3 is written in range of x according to each time step + + + return x0_xmax, t_values, u_values, dims_u + end +end + +# Euler is in a "conserved" form and the vector u contains (density, momentum, total energy) at each point x and time t +# So we want to convert these to (density, velocity, pressure) to calculate delta1 and delta2 according to [6] + +# Read the data and create the animation +# x0_xmax, t_values, u_values, dims_u = read_output_file("C:/Users/user/Documents/School/Sem4/softwareentwicklungspraktikum/shock_wave_detection/ShockwaveProperties.jl/example/data/euler_scenario_2.out") + +# Convert u_values to primitive variables. It's convert to density, velocity, pressure +function convert_to_primitive(u_values) + u_prim = zeros(size(u_values)) + for i in 1:size(u_values, 2) + for j in 1:size(u_values, 3) + # primitive_state_vector returns value without units + u_p_M_T = primitive_state_vector(u_values[:, i, j]; gas=DRY_AIR) + p_u = pressure(u_p_M_T[1], u_p_M_T[3]; gas=DRY_AIR) + # Store density + u_prim[1, i, j] = u_p_M_T[1] + # Convert Mach to m/s using speed_of_sound + u_prim[2, i, j] = u_p_M_T[2] * ustrip(speed_of_sound(u_p_M_T[3]; gas=DRY_AIR)) + # Strip the unit of pressure so that it can be stored in an empty array + u_prim[3, i, j] = ustrip(p_u) + end + end + + # Extract the density, velocity, and pressure fields + density_field = u_prim[1, :, :] + velocity_field = u_prim[2, :, :] + pressure_field = u_prim[3, :, :] + return density_field, velocity_field, pressure_field +end + +function create_animation(x0_xmax, t_values, density_field, velocity_field, pressure_field) + # Create a range of x values + x = range(x0_xmax[1], stop=x0_xmax[2], length=dims_u[2]) + + # Create a range of t values + t = t_values + + # Create an animation + anim = @animate for (i, t) in enumerate(t_values) + p1 = Plots.plot(x, density_field[:, i], title="Density at Time $t", xlabel="x", ylabel="Density(kg/m^3)", label = "Density across x", size=(800, 600)) + p2 = Plots.plot(x, velocity_field[:, i], title="Velocity at Time $t", xlabel="x", ylabel="Velocity(m/s)", label = "Velocity across x", size=(800, 600)) + p3 = Plots.plot(x, pressure_field[:, i], title="Pressure at Time $t", xlabel="x", ylabel="Pressure(Pa)", label = "Pressure across x", size=(800, 600)) + plot(p1, p2, p3, layout = (3, 1)) + end + # Save the animation as a gif + gif(anim, "density_velocity_pressure_over_time.gif", fps = 10) +end + +#density_field, velocity_field, pressure_field = convert_to_primitive(u_values) + +#create_animation(x0_xmax, t_values, density_field, velocity_field, pressure_field)