add some (not yet functional) implementation of fankine hugoniot

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.4"
 manifest_format = "2.0"
-project_hash = "9923b5a75b91e93f71a63c5f54fd461141cfe947"
+project_hash = "fb97f759a59eb4572eab12e1bcd17a2202ee2103"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"

Project.toml

@@ -6,6 +6,7 @@ version = "1.0.0-DEV"
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ShockwaveProperties = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

README.md

@@ -2,3 +2,16 @@
 
 [![Build Status](https://github.com/warisa-r/ShockwaveDetection.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/warisa-r/ShockwaveDetection.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/warisa-r/ShockwaveDetection.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/warisa-r/ShockwaveDetection.jl)
+
+
+## Brief description
+
+This package aims at processing the data from Euler's equation's solver [Euler2D.jl](https://github.com/STCE-at-RWTH/ShockwaveProperties.jl) and utilizing [ShockwaveProperties.jl](https://github.com/STCE-at-RWTH/ShockwaveProperties.jl) and detect where the shock is.
+Currently, it can
+- Detect 1D shock by detecting large gradients across properties with a customizable paramamer `threshold`
+- Visualize the change of properties along with shock positions
+
+## Goals
+- Develop 'better' methods to detect the shock
+- Develop visualization of shock position in x axis against time
+- Process 2D data (Waiting for the materials)

docs/src/index.md

@@ -17,7 +17,7 @@ You can install ShockDetection.jl using the Julia package manager. From the Juli
 ```@docs
 read_output_file
 convert_to_primitive
-detect_shock
+detect_normal_shock
 create_wave_animation
 create_wave_animation_with_shock
 ```

examples/1D_shock_scenario1_detection.jl

@@ -4,8 +4,8 @@ x0_xmax, t_values, u_values, dims_u = read_output_file("examples/data/euler_scen
 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)
+#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)
+shock_positions_over_time = detect_normal_shock(u_values, 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

@@ -7,5 +7,5 @@ density_field, velocity_field, pressure_field = convert_to_primitive(u_values)
 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)
+shock_positions_over_time = detect_normal_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

@@ -3,7 +3,7 @@ module ShockwaveDetection
 export write_output
 export read_output_file
 export convert_to_primitive
-export detect_shock
+export detect_normal_shock
 export create_wave_animation, create_wave_animation_with_shock
 
 

src/detect_shock.jl

@@ -1,5 +1,8 @@
 # Use Interpolations.jl to perform np.gradient-like operations
+using LinearAlgebra: I
+using Unitful: ustrip
 using Interpolations: interpolate, Gridded, Linear, gradient
+using ShockwaveProperties: state_behind, ConservedProps, DRY_AIR, mach_number, shock_temperature_ratio, shock_density_ratio, shock_pressure_ratio
 
 # Function to compute gradients using Interpolations.jl
 # without explicit function definition
@@ -11,8 +14,43 @@ function compute_gradients(arr::AbstractVector, x)
     return grad
 end
 
+#TODO: make this work
+function rankine_hugonoit(u_values_t, shock_locations)
+    # Define normal and tangential vectors for 1D case because state_behind needs it
+    n = [1]
+    t = [0]
+
+    # flux for the euler equations from ShockwaveProperties.jl's test since it's not exported
+    function F(u::ConservedProps)
+        v = u.ρv / u.ρ # velocity
+        P = pressure(u; gas = DRY_AIR)
+        return vcat(u.ρv', (u.ρv .* v' + I * P), (v .* (u.ρE + P))') # stack row vectors
+    end
+
+    for shock_location in shock_locations
+        u_L = ConservedProps(u_values_t[:, shock_location])
+        if mach_number(u_L; gas = DRY_AIR)[1] <= 0.0
+            println("Mach number leq 0 at shock location $shock_location and mach number is $(mach_number(u_L; gas = DRY_AIR))")
+        else
+            println("Mach number gt 0 at shock location $shock_location and mach number is $(mach_number(u_L; gas = DRY_AIR))")
+            println(shock_temperature_ratio(mach_number(u_L; gas = DRY_AIR), n; gas = DRY_AIR))
+            println(shock_pressure_ratio(mach_number(u_L; gas = DRY_AIR), n; gas = DRY_AIR))
+            println(shock_density_ratio(mach_number(u_L; gas = DRY_AIR), n; gas = DRY_AIR))
+            u_R = state_behind(u_L, n, t; gas = DRY_AIR)
+            # Check Rankine-Hugonoit conditions
+            if all(isapprox.(ustrip.(F(u_L) .* n), ustrip.(F(u_R) .* n); atol=1e-6))
+                println("Rankine-Hugonoit conditions are satisfied at shock location $shock_location")
+            else
+                println("Rankine-Hugonoit conditions are not satisfied at shock location $shock_location")
+            end
+        end
+
+    end
+
+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, threshold)
+function detect_normal_shocks_at_timestep(u_values_t, density_at_t, velocity_at_t, pressure_at_t, x, threshold)
 
     # Compute first gradients
     density_grad = compute_gradients(density_at_t, x)
@@ -26,16 +64,20 @@ function detect_shocks_at_timestep(density_at_t, velocity_at_t, pressure_at_t, x
 
     # Combine detections (common shock location across variables)
     shock_locations = intersect(intersect(shock_location_density, shock_location_velocity), shock_location_pressure)
+
+    #Comment this out everything works for now
+    #rankine_hugonoit(u_values_t, shock_locations)
 
     return shock_locations
 end
 
 """
-    detect_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values; threshold=0.5)
+    detect_normal_shock(density_field, velocity_field, pressure_field, x0_xmax, t_values; threshold=0.5)
 
 Detects shock positions over time based on given density, velocity, and pressure fields.
 
 # Arguments
+- `u_values::Array`: Array containing the conserved state variables at each point x and time t.
 - `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.
@@ -50,7 +92,7 @@ Detects shock positions over time based on given density, velocity, and pressure
 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)
+function detect_normal_shock(u_values, 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)
 
@@ -60,10 +102,10 @@ function detect_shock(density_field, velocity_field, pressure_field, x0_xmax, t_
         density_field_t = density_field[:, t_step]
         velocity_field_t = velocity_field[:, t_step]
         pressure_field_t = pressure_field[:, t_step]
+        u_values_t = u_values[:, :, t_step]
 
-        shock_positions = detect_shocks_at_timestep(density_field_t, velocity_field_t, pressure_field_t, x, threshold)
+        shock_positions = detect_normal_shocks_at_timestep(u_values_t, density_field_t, velocity_field_t, pressure_field_t, x, threshold)
         push!(shock_positions_over_time, shock_positions)
     end
-
     return shock_positions_over_time
 end