refactor everything

Project.toml

@@ -4,7 +4,10 @@ authors = ["Warisa <[email protected]> and contributors"]
 version = "1.0.0-DEV"
 
 [deps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ShockwaveProperties = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 julia = "1.9"

examples/1D_shock_scenario1_detection.jl

@@ -1 +1,11 @@
-using ShockwaveDetection
+using ShockwaveDetection
+
+x0_xmax, t_values, u_values, dims_u = read_output_file("examples/data/euler_scenario_1.out")
+density_field, velocity_field, pressure_field = convert_to_primitive(u_values)
+
+# Plot the animation of the properties over time
+anim = create_wave_animation(x0_xmax, t_values, density_field, velocity_field, pressure_field)
+
+#Plot the shock create_wave_animation
+shock_positions_over_time = detect_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values)
+anim_with_shock = create_wave_animation_with_shock(x0_xmax, t_values, density_field, velocity_field, pressure_field, shock_positions_over_time)

examples/1D_shock_scenario2_detection.jl

@@ -1 +1,11 @@
-using ShockwaveDetection
+using ShockwaveDetection
+
+x0_xmax, t_values, u_values, dims_u = read_output_file("examples/data/euler_scenario_2.out")
+density_field, velocity_field, pressure_field = convert_to_primitive(u_values)
+
+# Plot the animation of the properties over time
+anim = create_wave_animation(x0_xmax, t_values, density_field, velocity_field, pressure_field)
+
+#Plot the shock create_wave_animation
+shock_positions_over_time = detect_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values)
+anim_with_shock = create_wave_animation_with_shock(x0_xmax, t_values, density_field, velocity_field, pressure_field, shock_positions_over_time)

src/ShockwaveDetection.jl

@@ -1,5 +1,16 @@
 module ShockwaveDetection
 
-# Write your package code here.
+export write_output
+export read_output_file
+export convert_to_primitive
+export detect_shock
+export create_wave_animation, create_wave_animation_with_shock
+
+
+include("dummy.jl")
+include("input.jl")
+include("variable_utils.jl")
+include("detect_shock.jl")
+include("visualize.jl")
 
 end

src/detect_shock.jl

@@ -1,35 +1,25 @@
-include("visualize.jl")
-
 # Use Interpolations.jl to perform np.gradient-like operations
-using Interpolations
-
-using StatsPlots
+using Interpolations: interpolate, Gridded, Linear, gradient
 
 # 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))
+    grad = only.(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)
+function detect_shocks_at_timestep(density_at_t, velocity_at_t, pressure_at_t, x, threshold)
+
     # 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)
@@ -40,38 +30,40 @@ function detect_shocks_at_timestep(density_at_t, velocity_at_t, pressure_at_t, x
     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_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values; threshold=0.5)
 
-# Detect shocks for each time step and store the positions
-shock_positions_over_time = []
+Detects shock positions over time based on given density, velocity, and pressure fields.
 
-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
+# Arguments
+- `density_field::Matrix`: Matrix representing the density field over space and time.
+- `velocity_field::Matrix`: Matrix representing the velocity field over space and time.
+- `pressure_field::Matrix`: Matrix representing the pressure field over space and time.
+- `x0_xmax::Tuple`: Tuple containing the start and end positions of the spatial domain.
+- `t_values::Vector`: Vector containing time values.
+- `threshold::Float64`: Threshold value for detecting shocks. Default is `0.5`.
 
-# 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))
+# Returns
+- `shock_positions_over_time::Vector`: Vector containing shock positions at each time step.
+
+# Details
+This function iterates over each time step and detects shock positions using the `detect_shocks_at_timestep` function.
+
+"""
+function detect_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values; threshold = 0.5)
+    len_x = size(density_field, 1)
+    x = range(x0_xmax[1], stop=x0_xmax[2], length=len_x)
+
+    shock_positions_over_time = []
 
-    # 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)
+    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, threshold)
+        push!(shock_positions_over_time, shock_positions)
     end
 
-    plot(p1, p2, p3, layout = (3, 1))
-end
-
-# Save the animation as a gif
-gif(anim, "shock_over_time.gif", fps = 10)
+    return shock_positions_over_time
+end

src/dummy.jl

@@ -0,0 +1,3 @@
+function write_output()
+    return "Hello from Julia!"
+end

src/input.jl

@@ -0,0 +1,48 @@
+"""
+    read_output_file(filename)
+
+Reads data from an output file.
+
+# Arguments
+- `filename::String`: The path to the output file.
+
+# Returns
+- `x0_xmax::Vector{Float64}`: A vector containing the first and last x values.
+- `t_values::Vector{Float64}`: A vector containing the time values.
+- `u_values::Array{Float64, 3}`: A 3D array containing the u values.
+- `dims_u::Vector{Int}`: A vector containing the dimensions of the u values (N_u x N_x).
+
+This function reads data from the specified output file generated from 1d_plots.jl that use Euler2D.jl to solve, which is assumed to have a specific format. 
+It reads the dimensions of the u values, the first and last x values, the number of time steps, the time values, and the u values themselves. 
+It reshapes the u values into a 3D array and returns all the read data.
+"""
+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
+

src/variable_utils.jl

@@ -0,0 +1,42 @@
+using ShockwaveProperties: primitive_state_vector, pressure, speed_of_sound, DRY_AIR
+using Unitful: ustrip
+
+"""
+    convert_to_primitive(u_values)
+
+Converts conserved state variables to primitive variables.
+
+# Arguments
+- `u_values::Array{Float64, 3}`: A 3D array containing the conserved state variables (density, momentum, total energy) at each point x and time t.
+
+# Returns
+- `density_field::Array{Float64, 2}`: A 2D array containing the density field.
+- `velocity_field::Array{Float64, 2}`: A 2D array containing the velocity field.
+- `pressure_field::Array{Float64, 2}`: A 2D array containing the pressure field.
+
+Euler equations are typically represented in a "conserved" form, where the vector u contains (density, momentum, total energy) at each point x and time t. This function converts these conserved state variables to primitive variables (density, velocity, pressure) in order to calculate delta1 and delta2 according to [6].
+
+The `u_values` parameter is a 3D array representing the conserved state variables at each point x and time t. This function iterates over each time step and calculates the primitive state vector for density, velocity, and pressure using appropriate transformations. The resulting primitive variables are stored in separate arrays `density_field`, `velocity_field`, and `pressure_field`, and returned.
+"""
+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