diff --git a/Manifest.toml b/Manifest.toml index 21a0fb8..cf0b65b 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index 84d875e..b743756 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ 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] @@ -22,6 +23,7 @@ Documenter = "1.4.1" LinearAlgebra = "1" Plots = "1.40.4" Unitful = "1.20.0" +Statistics = "1.10.0" julia = "1.10" [extras] diff --git a/examples/obstacle.jl b/examples/obstacle.jl index c242093..8e55805 100644 --- a/examples/obstacle.jl +++ b/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") \ No newline at end of file +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) diff --git a/src/ShockwaveDetection.jl b/src/ShockwaveDetection.jl index 3d03677..6d871fd 100644 --- a/src/ShockwaveDetection.jl +++ b/src/ShockwaveDetection.jl @@ -1,6 +1,6 @@ module ShockwaveDetection -using Base.Threads +using Base.Threads: @threads export write_output export read_output_file, FlowData diff --git a/src/fitting.jl b/src/fitting.jl index 2df1042..ce1682a 100644 --- a/src/fitting.jl +++ b/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 diff --git a/src/input.jl b/src/input.jl index fa6065c..f4e2575 100644 --- a/src/input.jl +++ b/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 \ No newline at end of file diff --git a/src/pipeline.jl b/src/pipeline.jl index 463da2d..ebb87bb 100644 --- a/src/pipeline.jl +++ b/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) diff --git a/src/shock_point_detectors/2D/image_processing.jl b/src/shock_point_detectors/2D/image_processing.jl index aaf58c8..2a74456 100644 --- a/src/shock_point_detectors/2D/image_processing.jl +++ b/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) diff --git a/src/variable_utils.jl b/src/variable_utils.jl index 35f7f36..b769ae9 100644 --- a/src/variable_utils.jl +++ b/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