Fix the validation setup based on Marrone2011

examples/fluid/dam_break_2d.jl

@@ -113,15 +113,9 @@ info_callback = InfoCallback(interval=100)
 solution_prefix = ""
 saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix)
 
-# Save at certain timepoints which allows comparison to the results of Marrone et al.,
-# i.e. (1.5, 2.36, 3.0, 5.7, 6.45).
-# Please note that the images in Marrone et al. are obtained at a particle_spacing = H/320,
-# which takes between 2 and 4 hours.
-saving_paper = SolutionSavingCallback(save_times=[0.0, 0.371, 0.584, 0.743, 1.411, 1.597],
-                                      prefix="marrone_times")
-
 # This can be overwritten with `trixi_include`
 extra_callback = nothing
+extra_callback2 = nothing
 
 use_reinit = false
 density_reinit_cb = use_reinit ?
@@ -130,7 +124,8 @@ density_reinit_cb = use_reinit ?
 stepsize_callback = StepsizeCallback(cfl=0.9)
 
 callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback,
-                        density_reinit_cb, saving_paper)
+                        extra_callback2,
+                        density_reinit_cb)
 
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
             dt=1.0, # This is overwritten by the stepsize callback

src/schemes/boundary/system.jl

@@ -460,7 +460,8 @@ function system_data(system::BoundarySPHSystem, v_ode, u_ode, semi)
     u = wrap_u(u_ode, system, semi)
 
     coordinates = current_coordinates(u, system)
-    velocity = current_velocity(v, system)
+    velocity = [current_velocity(v, system, system.movement, i)
+                for i in 1:nparticles(system)]
     density = current_density(v, system)
     pressure = current_pressure(v, system)
 

test/validation/validation.jl

@@ -50,25 +50,25 @@
         if Sys.ARCH === :aarch64
             # MacOS ARM produces slightly different pressure values than x86.
             # Note that pressure values are in the order of 1e5.
-            @test isapprox(error_edac_P1, 0, atol=5e-6)
-            @test isapprox(error_edac_P2, 0, atol=4e-11)
-            @test isapprox(error_wcsph_P1, 0, atol=400.0)
-            @test isapprox(error_wcsph_P2, 0, atol=0.03)
+            @test isapprox(error_edac_P1, 0, atol=25.9)
+            @test isapprox(error_edac_P2, 0, atol=7.3e-5)
+            @test isapprox(error_wcsph_P1, 0, atol=eps())
+            @test isapprox(error_wcsph_P2, 0, atol=eps())
         elseif VERSION < v"1.11"
             # 1.10 produces slightly different pressure values than 1.11.
             # This is most likely due to muladd and FMA instructions in the
             # density diffusion update (inside the StaticArrays matrix-vector product).
             # Note that pressure values are in the order of 1e5.
-            @test isapprox(error_edac_P1, 0, atol=eps())
-            @test isapprox(error_edac_P2, 0, atol=eps())
-            @test isapprox(error_wcsph_P1, 0, atol=8.0)
-            @test isapprox(error_wcsph_P2, 0, atol=5e-4)
-        else
-            # Reference values are computed with 1.11
-            @test isapprox(error_edac_P1, 0, atol=eps())
-            @test isapprox(error_edac_P2, 0, atol=eps())
+            @test isapprox(error_edac_P1, 0, atol=30.9)
+            @test isapprox(error_edac_P2, 0, atol=0.00016)
             @test isapprox(error_wcsph_P1, 0, atol=eps())
             @test isapprox(error_wcsph_P2, 0, atol=eps())
+        else
+            # Reference values are computed with 1.11
+            @test isapprox(error_edac_P1, 0, atol=29.8)
+            @test isapprox(error_edac_P2, 0, atol=0.00015)
+            @test isapprox(error_wcsph_P1, 0, atol=0.037)
+            @test isapprox(error_wcsph_P2, 0, atol=1.5e-7)
         end
 
         # Ignore method redefinitions from duplicate `include("../validation_util.jl")`
@@ -78,7 +78,10 @@
             r"WARNING: Method definition linear_interpolation.*\n",
             r"WARNING: Method definition interpolated_mse.*\n",
             r"WARNING: Method definition extract_number_from_filename.*\n",
-            r"WARNING: Method definition extract_resolution_from_filename.*\n"
+            r"WARNING: Method definition extract_resolution_from_filename.*\n",
+            r"WARNING: Method definition pressure_probe.*\n",
+            r"WARNING: Method definition interpolated_pressure.*\n",
+            r"WARNING: Method definition max_x_coord.*\n"
         ]
         # Verify number of plots
         @test length(axs_edac[1].scene.plots) >= 2

validation/dam_break_2d/plot_dam_break_results.jl

@@ -10,6 +10,10 @@ using Glob
 using Printf
 using TrixiParticles
 
+# Set to save figures as SVG
+save_figures = false
+include_sim_results = false
+
 # Initial width of the fluid
 H = 0.6
 W = 2 * H
@@ -21,20 +25,23 @@ case_dir = joinpath(validation_dir(), "dam_break_2d")
 
 edac_reference_files = glob("validation_reference_edac*.json",
                             case_dir)
-edac_sim_files = glob("validation_result_dam_break_edac*.json",
-                      "out/")
+edac_sim_files = include_sim_results ?
+                 glob("validation_result_dam_break_edac*.json",
+                      "out/") : []
 
 merged_files = vcat(edac_reference_files, edac_sim_files)
 edac_files = sort(merged_files, by=extract_number_from_filename)
 
 wcsph_reference_files = glob("validation_reference_wcsph*.json",
                              case_dir)
-wcsph_sim_files = glob("validation_result_dam_break_wcsph*.json",
-                       "out/")
+wcsph_sim_files = include_sim_results ?
+                  glob("validation_result_dam_break_wcsph*.json",
+                       "out/") : []
 
 merged_files = vcat(wcsph_reference_files, wcsph_sim_files)
 wcsph_files = sort(merged_files, by=extract_number_from_filename)
 
+# Read in external reference data
 surge_front = CSV.read(joinpath(case_dir, "exp_surge_front.csv"), DataFrame)
 
 exp_P1 = CSV.read(joinpath(case_dir, "exp_pressure_sensor_P1.csv"), DataFrame)
@@ -43,48 +50,77 @@ exp_P2 = CSV.read(joinpath(case_dir, "exp_pressure_sensor_P2.csv"), DataFrame)
 sim_P1 = CSV.read(joinpath(case_dir, "sim_pressure_sensor_P1.csv"), DataFrame)
 sim_P2 = CSV.read(joinpath(case_dir, "sim_pressure_sensor_P2.csv"), DataFrame)
 
-n_sensors = 2
-fig = Figure(size=(1200, 1200))
-axs_edac = [Axis(fig[1, i], title="Sensor P$i with EDAC") for i in 1:n_sensors]
-axs_wcsph = [Axis(fig[3, i], title="Sensor P$i with WCSPH") for i in 1:n_sensors]
-ax_max_x_edac = Axis(fig[5, 1], title="Surge Front with EDAC")
-ax_max_x_wcsph = Axis(fig[5, 2], title="Surge Front with WCSPH")
-
-function plot_results(axs, ax_max, files)
+function plot_sensor_results(axs, files)
     for ax in axs
-        ax.xlabel = "Time"
-        ax.ylabel = "Pressure"
-        xlims!(ax, 0.0, 8.0)
-        ylims!(ax, -0.1, 1.0)
+        ax.xlabel = "Time [s]"
+        ax.ylabel = "P/(ρ g H)"
+        xlims!(ax, 2, 8)
+        ylims!(ax, -0.2, 1.0)
     end
 
-    # Define a regex to extract the sensor number from the key names
-    sensor_number_regex = r"pressure_P(\d+)_fluid_\d+"
+    for (idx, json_file) in enumerate(files)
+        println("Processing file: $json_file")
 
-    for (file_number, json_file) in enumerate(files)
-        json_data = JSON.parsefile(json_file)
-        time = json_data["pressure_P1_fluid_1"]["time"] .* normalization_factor_time
-        pressure_P1 = json_data["pressure_P1_fluid_1"]["values"] ./
+        jd = JSON.parsefile(json_file)
+        t = jd["interpolated_pressure_P1_fluid_1"]["time"] .* normalization_factor_time
+        pressure_P1 = jd["interpolated_pressure_P1_fluid_1"]["values"] /
                       normalization_factor_pressure
-        pressure_P2 = json_data["pressure_P2_fluid_1"]["values"] ./
+        pressure_P2 = jd["interpolated_pressure_P2_fluid_1"]["values"] /
                       normalization_factor_pressure
+        probe_P1 = jd["particle_pressure_P1_boundary_1"]["values"] /
+                   normalization_factor_pressure
+        probe_P2 = jd["particle_pressure_P2_boundary_1"]["values"] /
+                   normalization_factor_pressure
+        lab = occursin("reference", json_file) ? "Ref. " : ""
+        res = extract_resolution_from_filename(json_file)
+
+        lines!(axs[1], t, pressure_P1; label="$lab dp=$res",
+               color=idx, colormap=:tab10, colorrange=(1, 10))
+        lines!(axs[2], t, pressure_P2; label="$lab dp=$res",
+               color=idx, colormap=:tab10, colorrange=(1, 10))
+        lines!(axs[3], t, probe_P1; label="$lab dp=$res",
+               color=idx, colormap=:tab10, colorrange=(1, 10))
+        lines!(axs[4], t, probe_P2; label="$lab dp=$res",
+               color=idx, colormap=:tab10, colorrange=(1, 10))
+    end
+end
 
-        label_prefix = occursin("reference", json_file) ? "Reference " : ""
-
-        lines!(axs[1], time, pressure_P1,
-               label="$(label_prefix)dp=$(extract_resolution_from_filename(json_file))",
-               color=file_number, colormap=:tab10, colorrange=(1, 10))
-        lines!(axs[2], time, pressure_P2,
-               label="$(label_prefix)dp=$(extract_resolution_from_filename(json_file))",
-               color=file_number, colormap=:tab10, colorrange=(1, 10))
-        value = json_data["max_x_coord_fluid_1"]
-        lines!(ax_max, value["time"] .* sqrt(9.81), Float64.(value["values"]) ./ W,
-               label="$(label_prefix)dp=$(extract_resolution_from_filename(json_file))")
+function plot_surge_results(ax, files)
+    ax.xlabel = "Time [s]"
+    ax.ylabel = "x / W"
+    xlims!(ax, 0, 3)
+    ylims!(ax, 1, 3)
+
+    for (idx, json_file) in enumerate(files)
+        jd = JSON.parsefile(json_file)
+        lab = occursin("reference", json_file) ? "Ref. " : ""
+        val = jd["max_x_coord_fluid_1"]
+        lines!(ax, val["time"] .* sqrt(9.81),
+               Float64.(val["values"]) ./ W;
+               label="$lab dp=$(extract_resolution_from_filename(json_file))",
+               color=idx, colormap=:tab10, colorrange=(1, 10))
     end
+
+    # experimental reference
+    scatter!(ax, surge_front.time, surge_front.surge_front;
+             color=:black, marker=:utriangle, markersize=6,
+             label="Martin & Moyce 1952 (exp)")
 end
 
-plot_results(axs_edac, ax_max_x_edac, edac_files)
-plot_results(axs_wcsph, ax_max_x_wcsph, wcsph_files)
+# ------------------------------------------------------------
+# 1) Pressure-sensor figure
+# ------------------------------------------------------------
+n_sensors = 4
+fig_sensors = Figure(size=(2400, 1000))
+axs_edac = [Axis(fig_sensors[1, i],
+                 title=(i>2) ? "Boundary values at P$(i-2) (EDAC)" : "Sensor P$i (EDAC)")
+            for i in 1:n_sensors]
+axs_wcsph = [Axis(fig_sensors[3, i],
+                  title=(i>2) ? "Boundary values at P$(i-2) (WCSPH)" : "Sensor P$i (WCSPH)")
+             for i in 1:n_sensors]
+
+plot_sensor_results(axs_edac, edac_files)
+plot_sensor_results(axs_wcsph, wcsph_files)
 
 # Plot reference values
 function plot_experiment(ax, time, data, label, color=:black, marker=:utriangle,
@@ -98,37 +134,62 @@ end
 
 # Plot for Pressure Sensor P1
 plot_experiment(axs_edac[1], exp_P1.time, exp_P1.P1, "Buchner 2002 (exp)")
-plot_simulation(axs_edac[1], sim_P1.time, sim_P1.h320, "Marrone et al. 2011 (sim)")
+plot_simulation(axs_edac[1], sim_P1.time, sim_P1.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
 plot_experiment(axs_wcsph[1], exp_P1.time, exp_P1.P1, "Buchner 2002 (exp)")
-plot_simulation(axs_wcsph[1], sim_P1.time, sim_P1.h320, "Marrone et al. 2011 (sim)")
+plot_simulation(axs_wcsph[1], sim_P1.time, sim_P1.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
+plot_experiment(axs_edac[3], exp_P1.time, exp_P1.P1, "Buchner 2002 (exp)")
+plot_simulation(axs_edac[3], sim_P1.time, sim_P1.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
+plot_experiment(axs_wcsph[3], exp_P1.time, exp_P1.P1, "Buchner 2002 (exp)")
+plot_simulation(axs_wcsph[3], sim_P1.time, sim_P1.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
 
 # Plot for Pressure Sensor P2
 plot_experiment(axs_edac[2], exp_P2.time, exp_P2.P2, "Buchner 2002 (exp)")
-plot_simulation(axs_edac[2], sim_P2.time, sim_P2.h320, "Marrone et al. 2011 (sim)")
+plot_simulation(axs_edac[2], sim_P2.time, sim_P2.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
 plot_experiment(axs_wcsph[2], exp_P2.time, exp_P2.P2, "Buchner 2002 (exp)")
-plot_simulation(axs_wcsph[2], sim_P2.time, sim_P2.h320, "Marrone et al. 2011 (sim)")
-
-# Plot for Surge Front
-for ax_max in [ax_max_x_edac, ax_max_x_wcsph]
-    ax_max.xlabel = "Time"
-    ax_max.ylabel = "Surge Front"
-    xlims!(ax_max, 0.0, 3.0)
-    ylims!(ax_max, 1, 3.0)
-    plot_experiment(ax_max, surge_front.time, surge_front.surge_front,
-                    "Martin and Moyce 1952 (exp)")
-end
+plot_simulation(axs_wcsph[2], sim_P2.time, sim_P2.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
+plot_experiment(axs_edac[4], exp_P2.time, exp_P2.P2, "Buchner 2002 (exp)")
+plot_simulation(axs_edac[4], sim_P2.time, sim_P2.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
+plot_experiment(axs_wcsph[4], exp_P2.time, exp_P2.P2, "Buchner 2002 (exp)")
+plot_simulation(axs_wcsph[4], sim_P2.time, sim_P2.h320,
+                "Marrone et al. 2011 dp=0.001875 (sim)")
 
 for (i, ax) in enumerate(axs_edac)
-    Legend(fig[2, i], ax; tellwidth=false, orientation=:horizontal, valign=:top, nbanks=3)
+    Legend(fig_sensors[2, i], ax; tellwidth=false, orientation=:horizontal, valign=:top,
+           nbanks=3)
 end
 for (i, ax) in enumerate(axs_wcsph)
-    Legend(fig[4, i], ax; tellwidth=false, orientation=:horizontal, valign=:top, nbanks=3)
+    Legend(fig_sensors[4, i], ax; tellwidth=false, orientation=:horizontal, valign=:top,
+           nbanks=3)
+end
+
+if save_figures
+    save("dam_break_pressure.svg", fig_sensors)
+else
+    display(fig_sensors)
+end
+
+# ------------------------------------------------------------
+# 2) Surge-front figure
+# ------------------------------------------------------------
+fig_surge = Figure(size=(1200, 400))
+ax_surge_edac = Axis(fig_surge[1, 1], title="Surge Front – EDAC")
+ax_surge_wcsph = Axis(fig_surge[1, 2], title="Surge Front – WCSPH")
+
+plot_surge_results(ax_surge_edac, edac_files)
+plot_surge_results(ax_surge_wcsph, wcsph_files)
+
+Legend(fig_surge[2, 1], ax_surge_edac; orientation=:horizontal, valign=:top, nbanks=3)
+Legend(fig_surge[2, 2], ax_surge_wcsph; orientation=:horizontal, valign=:top, nbanks=3)
+
+if save_figures
+    save("dam_break_surge_front.svg", fig_surge)
+else
+    display(fig_surge)
 end
-Legend(fig[6, 1], ax_max_x_edac, tellwidth=false, orientation=:horizontal, valign=:top,
-       nbanks=2)
-Legend(fig[6, 2], ax_max_x_wcsph, tellwidth=false, orientation=:horizontal, valign=:top,
-       nbanks=2)
-
-fig
-# Uncomment to save the figure
-# save("dam_break_validation.svg", fig)

validation/dam_break_2d/sensors.jl

@@ -0,0 +1,69 @@
+function max_x_coord(system, data, t)
+    return maximum(j -> data.coordinates[1, j], axes(data.coordinates, 2))
+end
+
+function interpolated_pressure(coord_top, coord_bottom, v_ode, u_ode, t, system, semi) end
+
+function interpolated_pressure(coord_top, coord_bottom, v_ode, u_ode, t,
+                               system::TrixiParticles.FluidSystem, semi)
+    # use at least 5 interpolation points for low resolution simulations
+    # otherwise use at least the number of particles present
+    n_interpolation_points = min(5, Int(ceil(sensor_size / particle_spacing)))
+    interpolated_values = interpolate_line(coord_top, coord_bottom,
+                                           n_interpolation_points, semi, system, v_ode,
+                                           u_ode,
+                                           smoothing_length=2.0 *
+                                                            TrixiParticles.initial_smoothing_length(system),
+                                           clip_negative_pressure=true)
+    return sum(map(x -> isnan(x) ? 0.0 : x, interpolated_values.pressure)) /
+           n_interpolation_points
+end
+
+function pressure_probe(coord_top, coord_bottom, v_ode, u_ode, t, system, semi) end
+
+function pressure_probe(coord_top, coord_bottom, v_ode, u_ode, t,
+                        system::TrixiParticles.BoundarySystem, semi)
+    # The sensor is at the right wall, so its x-coordinate is the same for top and bottom.
+    x_sensor = coord_top[1]
+
+    # Use the initial particle spacing as a reference for the thickness of the averaging region.
+    # A thickness of one or two particle spacings is usually a good choice.
+    region_thickness = 2.0 * particle_spacing
+
+    # Define the rectangular region for averaging
+    x_min = x_sensor - region_thickness
+    x_max = x_sensor + region_thickness
+    y_min = coord_bottom[2]
+    y_max = coord_top[2]
+
+    sum_of_pressures = 0.0
+    num_particles_in_region = 0
+
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
+    # Iterate over each particle in the specified fluid system
+    for particle_idx in TrixiParticles.eachparticle(system)
+        pc = TrixiParticles.current_coords(u, system, particle_idx)
+
+        # Get coordinates for the current particle from the 1D vector
+        px = pc[1] # x-coordinate
+        py = pc[2] # y-coordinate
+
+        # Check if the particle is inside the sensor's rectangular region
+        if (x_min <= px <= x_max) && (y_min <= py <= y_max)
+            # Add its pressure to the sum and increment the count
+            sum_of_pressures += TrixiParticles.current_pressure(v, system, particle_idx)
+            num_particles_in_region += 1
+        end
+    end
+
+    # If no particles are in the region (e.g., before the water hits the wall),
+    # the pressure is zero.
+    if num_particles_in_region == 0
+        return 0.0
+    end
+
+    # Return the calculated average pressure
+    return sum_of_pressures / num_particles_in_region
+end