diff --git a/validation/dam_break_2d/validation_dam_break_2d.jl b/validation/dam_break_2d/validation_dam_break_2d.jl index bfe216fde..4bf57d59e 100644 --- a/validation/dam_break_2d/validation_dam_break_2d.jl +++ b/validation/dam_break_2d/validation_dam_break_2d.jl @@ -115,12 +115,12 @@ run_file_edac_name = joinpath("out", reference_data = JSON.parsefile(reference_file_edac_name) run_data = JSON.parsefile(run_file_edac_name) -error_edac_P1 = interpolated_mse(reference_data["pressure_P1_fluid_1"]["time"], +error_edac_P1 = interpolated_mre(reference_data["pressure_P1_fluid_1"]["time"], reference_data["pressure_P1_fluid_1"]["values"], run_data["pressure_P1_fluid_1"]["time"], run_data["pressure_P1_fluid_1"]["values"]) -error_edac_P2 = interpolated_mse(reference_data["pressure_P2_fluid_1"]["time"], +error_edac_P2 = interpolated_mre(reference_data["pressure_P2_fluid_1"]["time"], reference_data["pressure_P2_fluid_1"]["values"], run_data["pressure_P2_fluid_1"]["time"], run_data["pressure_P2_fluid_1"]["values"]) @@ -150,12 +150,12 @@ run_file_wcsph_name = joinpath("out", reference_data = JSON.parsefile(reference_file_wcsph_name) run_data = JSON.parsefile(run_file_wcsph_name) -error_wcsph_P1 = interpolated_mse(reference_data["pressure_P1_fluid_1"]["time"], +error_wcsph_P1 = interpolated_mre(reference_data["pressure_P1_fluid_1"]["time"], reference_data["pressure_P1_fluid_1"]["values"], run_data["pressure_P1_fluid_1"]["time"], run_data["pressure_P1_fluid_1"]["values"]) -error_wcsph_P2 = interpolated_mse(reference_data["pressure_P2_fluid_1"]["time"], +error_wcsph_P2 = interpolated_mre(reference_data["pressure_P2_fluid_1"]["time"], reference_data["pressure_P2_fluid_1"]["values"], run_data["pressure_P2_fluid_1"]["time"], run_data["pressure_P2_fluid_1"]["values"]) diff --git a/validation/validation_util.jl b/validation/validation_util.jl index 01db69635..56433f0f0 100644 --- a/validation/validation_util.jl +++ b/validation/validation_util.jl @@ -1,3 +1,9 @@ +# Perform linear interpolation to find a value at `interpolation_point` using arrays `x` and `y`. +# +# Arguments: +# - `x` : The array of abscissas (e.g., time points). +# - `y` : The array of ordinates (e.g., data points corresponding to `x`). +# - `interpolation_point` : The point at which to interpolate the data. function linear_interpolation(x, y, interpolation_point) if !(first(x) <= interpolation_point <= last(x)) throw(ArgumentError("`interpolation_point` at $interpolation_point is outside the interpolation range")) @@ -12,6 +18,15 @@ function linear_interpolation(x, y, interpolation_point) return y[i] + slope * (interpolation_point - x[i]) end + +# Calculate the mean squared error (MSE) between interpolated simulation values and reference values +# over a common time range. +# +# Arguments: +# - `reference_time` : Time points for the reference data. +# - `reference_values` : Data points for the reference data. +# - `simulation_time` : Time points for the simulation data. +# - `simulation_values` : Data points for the simulation data. function interpolated_mse(reference_time, reference_values, simulation_time, simulation_values) if last(simulation_time) > last(reference_time) @@ -34,6 +49,37 @@ function interpolated_mse(reference_time, reference_values, simulation_time, return mse end +# Calculate the mean relative error (MRE) between interpolated simulation values and reference values +# over a common time range. +# +# Arguments: +# - `reference_time` : Time points for the reference data. +# - `reference_values` : Data points for the reference data. +# - `simulation_time` : Time points for the simulation data. +# - `simulation_values` : Data points for the simulation data. +function interpolated_mre(reference_time, reference_values, simulation_time, simulation_values) + if last(simulation_time) > last(reference_time) + @warn "simulation time range is larger than reference time range. " * + "Only checking values within reference time range." + end + + # Remove reference time points outside the simulation time + start = searchsortedfirst(reference_time, first(simulation_time)) + end_ = searchsortedlast(reference_time, last(simulation_time)) + common_time_range = reference_time[start:end_] + + # Interpolate simulation data at the common time points + interpolated_values = [linear_interpolation(simulation_time, simulation_values, t) + for t in common_time_range] + + filtered_values = reference_values[start:end_] + + # Calculate MRE only over the common time range + mre = sum(abs.(interpolated_values .- filtered_values) ./ abs.(filtered_values)) / length(common_time_range) + return mre +end + + function extract_number_from_filename(filename) # This regex matches the last sequence of digits in the filename m = match(r"(\d+)(?!.*\d)", filename)