Handle obstacle from .celltape files

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "1bf254ffd9f7a81cec2510312af1c565bbc487e0"
+project_hash = "40441ff4083380f812414e462d16bb6e3dc99fb1"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"

Project.toml

@@ -15,13 +15,15 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ShockwaveProperties = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Documenter = "1.4.1"
 LinearAlgebra = "1"
 Plots = "1.40.4"
 Unitful = "1.20.0"
+Statistics = "1.10.0"
 julia = "1.10"
 
 [extras]

examples/obstacle.jl

@@ -1,5 +1,11 @@
 using ShockwaveDetection
 using ShockwaveProperties
-using Euler2D
+using Euler2D:Euler2D
 
-load_cell_sim("examples/data/obstacle/circular_obstacle_radius_1.celltape")
+flow_data = FlowData("examples/data/obstacle/circular_obstacle_radius_1.celltape")
+point_detect_algo = ImageProcessingShockDetectionAlgo(0.5, :prewitt)
+dbscan_algo = DBSCANAlgo(0.25, 3, 10)
+
+detection = detect(flow_data, point_detect_algo, dbscan_algo)
+
+plot_shock_fits_over_time(flow_data, detection, false)

src/ShockwaveDetection.jl

@@ -1,6 +1,6 @@
 module ShockwaveDetection
 
-using Base.Threads
+using Base.Threads: @threads
 
 export write_output
 export read_output_file, FlowData

src/fitting.jl

@@ -37,6 +37,13 @@ function line_model(xy, p)
     return m*x .+ b .- y
 end
 
+
+function hline_model(xy, p)
+    b = p
+    y = xy[:, 2]
+    return y .- b
+end
+
 function vline_model(xy, p)
     c = p
     x = xy[:, 1]
@@ -58,9 +65,9 @@ function fit_shock_cluster(cluster)
 
 
     xy = cluster_to_data_points(cluster)
-    models = [vline_model, line_model, circle_model, parabola_model] # Use only these three firsts
+    models = [vline_model, hline_model, line_model, circle_model, parabola_model] # Use only these three firsts
     #TODO: better parameter initialization from boundary conditions or information about the cluster??
-    p0s = [[1.0], [1.0, 1.0], [0.0, 0.0, 1.0], [1.0, 1.0, 1.0]]  # Initial parameters for each model
+    p0s = [[1.0], [1.0], [1.0, 1.0], [0.0, 0.0, 1.0], [1.0, 1.0, 1.0]]  # Initial parameters for each model
 
     best_fit = nothing
     least_error = Inf
@@ -203,6 +210,24 @@ function calculate_normal_vector(fit::Fitting, evenly_spaced_range, flow_data, t
         normals_x = [end_x - start_x]
         normals_y = [0]
 
+    elseif fit.model == hline_model
+        b = fit.parameters[1]
+        start_y = b
+        # Find where b is in the y direction
+        y = range(bounds[2][1], bounds[2][2], length=ncells[2])
+        differences = abs.(y .- b)
+        index_of_b = argmin(differences)
+
+        # Check if the density is increasing or decreasing in the y direction
+        # to identify which side is ahead of the shock (low density) and which side is behind the shock (high density)
+        if density_field_t[round(Int, ncells[1] / 2), index_of_b-5] > density_field_t[round(Int, ncells[1] / 2), index_of_b+5]
+            end_y = b + 1
+        else
+            end_y = b - 1
+        end
+        normals_x = [0]
+        normals_y = [end_y - start_y]
+
     elseif fit.model == circle_model
         angles = evenly_spaced_range
 

src/input.jl

@@ -4,16 +4,18 @@ using LinearAlgebra
 using ShockwaveProperties
 using Unitful
 
-struct FlowData{N, NAXES, T}
+
+struct FlowData{N, T}
     ncells::NTuple{N, Int}
     nsteps::Int
     bounds::NTuple{N, Tuple{T, T}}
     tsteps::Vector{T}
-    u::Array{T, NAXES}
-    density_field::Array{T} # Dimension NAXES-1 but Julia doesn't allow this
+    u::Union{Array{T}, Nothing} # Nothing for .celltape data since all other important data are stored in other attributes already
+    density_field::Array{T}
     velocity_field::Array{T}
     pressure_field::Array{T}
     mach_to_m_s::Bool
+    cell_ids::Union{Matrix{Int64}, Nothing}
 end
 
 
@@ -25,11 +27,21 @@ function FlowData(file_path::String, mach_to_m_s=true)
         bounds = sim_data.bounds
         tsteps = sim_data.tsteps
         u = sim_data.u
+        cell_ids = nothing
         density_field, velocity_field, pressure_field = convert_to_primitive(u, ncells, nsteps, mach_to_m_s)
+    elseif endswith(file_path, ".celltape")
+        sim_data = Euler2D.load_cell_sim(file_path)
+        ncells = sim_data.ncells
+        nsteps = sim_data.nsteps
+        bounds = sim_data.bounds
+        tsteps = sim_data.tsteps
+        cell_ids = sim_data.cell_ids
+        u = nothing
+        density_field, velocity_field, pressure_field = convert_to_primitive(sim_data, nsteps, mach_to_m_s)
     else
-        error("Unsupported file type. Please provide a .tape file.")
+        error("Unsupported file type. Please provide a .tape  or .celltape file.")
     end
-    return FlowData(ncells, nsteps, bounds, tsteps, u, density_field, velocity_field, pressure_field, mach_to_m_s)
+    return FlowData(ncells, nsteps, bounds, tsteps, u, density_field, velocity_field, pressure_field, mach_to_m_s, cell_ids)
 end
 
 # TODO: if the initial condition is given like the scripts, convert EulerSim to FlowData

src/pipeline.jl

@@ -8,7 +8,11 @@ end
 
 # Pipeline is only for 2D shock detection now
 function detect(flow_data::FlowData, shock_point_algo::Abstract2DShockDetectionAlgo, cluster_algo::DBSCANAlgo)
-    shock_positions_over_time = detect_shock_points(flow_data, shock_point_algo)
+    has_obstacle = false
+    if isnothing(flow_data.u)
+        has_obstacle = true
+    end
+    shock_positions_over_time = detect_shock_points(flow_data, shock_point_algo, has_obstacle)
     shock_clusters_over_time = cluster_shock_points(cluster_algo, shock_positions_over_time, flow_data)
     shock_fits_over_time = fit_shock_clusters_over_time(shock_clusters_over_time)
     return ShockDetectionResult2D(shock_positions_over_time, shock_clusters_over_time, shock_fits_over_time)

src/shock_point_detectors/2D/image_processing.jl

@@ -1,4 +1,32 @@
 using Images
+using Statistics: mean
+
+# Function to replace NaNs with the mean of neighboring values for cells near obstacles and obstacles cells
+function replace_nan_with_mean!(matrix)
+    for i in 1:size(matrix, 1)
+        for j in 1:size(matrix, 2)
+            if isnan(matrix[i, j])
+                # Get the neighborhood (3x3 window) around the NaN value since the kernel size is 3x3
+                imin = max(1, i-1)
+                imax = min(size(matrix, 1), i+1)
+                jmin = max(1, j-1)
+                jmax = min(size(matrix, 2), j+1)
+
+                # Extract the valid (non-NaN) neighbors
+                neighbors = matrix[imin:imax, jmin:jmax]
+                valid_neighbors = neighbors[.!isnan.(neighbors)]
+
+                # Replace NaN with the mean of valid neighbors
+                if !isempty(valid_neighbors)
+                    matrix[i, j] = mean(valid_neighbors)
+                else # If all neighbors are NaN (in the obstacle zone)
+                    matrix[i, j] = 0.0
+                end
+            end
+        end
+    end
+end
+
 
 struct ImageProcessingShockDetectionAlgo <: Abstract2DShockDetectionAlgo
     threshold::Float64
@@ -58,7 +86,7 @@ function detect_discon_at_timestep(density_at_t, velocity_at_t, pressure_at_t, a
     return high_gradient_intersection
 end
 
-function detect_shock_points(flow_data::FlowData, alg::ImageProcessingShockDetectionAlgo)
+function detect_shock_points(flow_data::FlowData, alg::ImageProcessingShockDetectionAlgo, has_obstacle)
     # Unpack all the values from the detector
     nsteps = flow_data.nsteps
     density_field = flow_data.density_field
@@ -71,8 +99,16 @@ function detect_shock_points(flow_data::FlowData, alg::ImageProcessingShockDetec
         density_field_t = density_field[:, :, t_step]
         velocity_field_x_t = velocity_field[1, :, :, t_step]
         velocity_field_y_t = velocity_field[2, :, :, t_step]
-        velocity_field_t = sqrt.(velocity_field_x_t.^2 + velocity_field_y_t.^2)
         pressure_field_t = pressure_field[:, :, t_step]
+
+        if has_obstacle
+            replace_nan_with_mean!(density_field_t)
+            replace_nan_with_mean!(velocity_field_x_t)
+            replace_nan_with_mean!(velocity_field_y_t)
+            replace_nan_with_mean!(pressure_field_t)
+        end
+
+        velocity_field_t = sqrt.(velocity_field_x_t.^2 + velocity_field_y_t.^2) # Calculate magnitude after handling replacement of NaNs
 
         # Use the simple shock detection algorithm to detect the normal shock
         shock_positions_over_time[t_step] = detect_discon_at_timestep(density_field_t, velocity_field_t, pressure_field_t, alg)

src/variable_utils.jl

@@ -1,4 +1,5 @@
 using ShockwaveProperties: primitive_state_vector, pressure, speed_of_sound, DRY_AIR
+using Euler2D: Euler2D
 using Unitful: ustrip
 
 """
@@ -77,6 +78,38 @@ function convert_to_primitive(u_values::Array{T, 4}, ncells, nsteps, mach_to_m_s
     return density_field, velocity_field, pressure_field
 end
 
+function convert_to_primitive(sim_data, nsteps, mach_to_m_s=false)
+    density_field = []
+    velocity_field = []
+    pressure_field = []
+
+
+    for t in 1:nsteps
+        density_t = [x !== nothing ? ustrip(x) : NaN for x in Euler2D.density_field(sim_data, t)]
+        pressure_t = Euler2D.pressure_field(sim_data, t, DRY_AIR)
+        velocity_t = [x !== nothing ? ustrip(x) : NaN for x in Euler2D.velocity_field(sim_data, t)]
+
+       # TODO: find a way to make this work properly. Since speed of sound needs density or pressure and we have no access to temperature 
+        if mach_to_m_s
+            speed_of_sound_t = [x !== nothing ? ustrip(speed_of_sound(ustrip(x), DRY_AIR)) : NaN for x in pressure_t]
+            #velocity_t = velocity_t .* speed_of_sound_t
+        end
+
+        pressure_t = [x !== nothing ? ustrip(x) : NaN for x in pressure_t]
+
+        push!(density_field, density_t)
+        push!(pressure_field, pressure_t)
+        push!(velocity_field, velocity_t)
+    end
+
+    # Convert lists of lists into multidimensional arrays
+    density_field_array = cat(density_field..., dims=3)
+    velocity_field_array = cat(velocity_field..., dims=4)
+    pressure_field_array = cat(pressure_field..., dims=3)
+
+    return density_field_array, velocity_field_array, pressure_field_array
+end
+
 """
     cartesian_index_to_xy(shock_positions_t, x, y) -> Matrix