add some files from prior developement

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"

Project.toml

@@ -3,6 +3,9 @@ uuid = "a75b99cc-be56-4638-99bc-d89fa43b9ca1"
 authors = ["Warisa <[email protected]> and contributors"]
 version = "1.0.0-DEV"
 
+[deps]
+ShockwaveProperties = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
+
 [compat]
 julia = "1.9"
 

examples/1D_shock_scenario1_detection.jl

@@ -0,0 +1 @@
+using ShockwaveDetection

examples/1D_shock_scenario2_detection.jl

@@ -0,0 +1 @@
+using ShockwaveDetection

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)

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)