diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 51b13c681..eb04b4e61 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -51,14 +51,32 @@ results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, so (density = ...) ``` """ -function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, sol; +function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, + sol::ODESolution; smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false) # Filter out particles without neighbors filter_no_neighbors = true + v_ode, u_ode = sol.u[end].x + + results, _, _ = interpolate_plane_2d(min_corner, max_corner, resolution, + semi, ref_system, v_ode, u_ode, + filter_no_neighbors, smoothing_length, cut_off_bnd, + clip_negative_pressure) + + return results +end + +function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, + v_ode, u_ode; + smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, clip_negative_pressure=false) + # Filter out particles without neighbors + filter_no_neighbors = true + results, _, _ = interpolate_plane_2d(min_corner, max_corner, resolution, - semi, ref_system, sol, filter_no_neighbors, - smoothing_length, cut_off_bnd, + semi, ref_system, v_ode, u_ode, + filter_no_neighbors, smoothing_length, cut_off_bnd, clip_negative_pressure) return results @@ -112,14 +130,28 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref), results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, sol) ``` """ -function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system, sol; +function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system, + sol::ODESolution; clip_negative_pressure=false, + smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, + output_directory="out", filename="plane") + v_ode = sol.u[end].x[1] + u_ode = sol.u[end].x[2] + + interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system, + v_ode, u_ode; clip_negative_pressure, + smoothing_length, cut_off_bnd, output_directory, filename) +end + +function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system, + v_ode, u_ode; smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false, output_directory="out", filename="plane") # Don't filter out particles without neighbors to keep 2D grid structure filter_no_neighbors = false results, x_range, y_range = interpolate_plane_2d(min_corner, max_corner, resolution, - semi, ref_system, sol, + semi, ref_system, v_ode, u_ode, filter_no_neighbors, smoothing_length, cut_off_bnd, clip_negative_pressure) @@ -135,9 +167,9 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ end end -function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, sol, - filter_no_neighbors, smoothing_length, cut_off_bnd, - clip_negative_pressure) +function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, + v_ode, u_ode, filter_no_neighbors, smoothing_length, + cut_off_bnd, clip_negative_pressure) dims = length(min_corner) if dims != 2 || length(max_corner) != 2 throw(ArgumentError("function is intended for 2D coordinates only")) @@ -157,7 +189,7 @@ function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_syst # Generate points within the plane points_coords = [SVector(x, y) for x in x_range, y in y_range] - results = interpolate_point(points_coords, semi, ref_system, sol, + results = interpolate_point(points_coords, semi, ref_system, v_ode, u_ode, smoothing_length=smoothing_length, cut_off_bnd=cut_off_bnd, clip_negative_pressure=clip_negative_pressure) @@ -226,9 +258,21 @@ results = interpolate_plane_3d([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0] (density = ...) ``` """ -function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol; +function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, + sol::ODESolution; smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false) + v_ode = sol.u[end].x[1] + u_ode = sol.u[end].x[2] + + interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, + v_ode, u_ode; smoothing_length, cut_off_bnd, + clip_negative_pressure) +end + +function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, + v_ode, u_ode; smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, clip_negative_pressure=false) # Verify that points are in 3D space if length(point1) != 3 || length(point2) != 3 || length(point3) != 3 throw(ArgumentError("all points must be 3D coordinates")) @@ -270,7 +314,7 @@ function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_syst end # Interpolate using the generated points - results = interpolate_point(points_coords, semi, ref_system, sol, + results = interpolate_point(points_coords, semi, ref_system, v_ode, u_ode, smoothing_length=smoothing_length, cut_off_bnd=cut_off_bnd, clip_negative_pressure=clip_negative_pressure) @@ -335,8 +379,17 @@ results = interpolate_line([1.0, 0.0], [1.0, 1.0], 5, semi, ref_system, sol) (density = ...) ``` """ -function interpolate_line(start, end_, n_points, semi, ref_system, sol; endpoint=true, - smoothing_length=ref_system.smoothing_length, +function interpolate_line(start, end_, n_points, semi, ref_system, sol::ODESolution; + endpoint=true, smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, clip_negative_pressure=false) + v_ode = sol.u[end].x[1] + u_ode = sol.u[end].x[2] + + interpolate_line(start, end_, n_points, semi, ref_system, v_ode, u_ode; + endpoint, smoothing_length, cut_off_bnd, clip_negative_pressure) +end +function interpolate_line(start, end_, n_points, semi, ref_system, v_ode, u_ode; + endpoint=true, smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false) start_svector = SVector{ndims(ref_system)}(start) end_svector = SVector{ndims(ref_system)}(end_) @@ -346,10 +399,9 @@ function interpolate_line(start, end_, n_points, semi, ref_system, sol; endpoint points_coords = points_coords[2:(end - 1)] end - return interpolate_point(points_coords, semi, ref_system, sol, + return interpolate_point(points_coords, semi, ref_system, v_ode, u_ode; smoothing_length=smoothing_length, - cut_off_bnd=cut_off_bnd, - clip_negative_pressure=clip_negative_pressure) + cut_off_bnd=cut_off_bnd, clip_negative_pressure) end @doc raw""" @@ -411,9 +463,19 @@ results = interpolate_point(points, semi, ref_system, sol) - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as accurate as a real surface reconstruction. """ -function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi, ref_system, - sol; smoothing_length=ref_system.smoothing_length, - cut_off_bnd=true, clip_negative_pressure=false) +@inline function interpolate_point(point_coords, semi, ref_system, sol::ODESolution; + smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, clip_negative_pressure=false) + v_ode = sol.u[end].x[1] + u_ode = sol.u[end].x[2] + interpolate_point(point_coords, semi, ref_system, v_ode, u_ode; + smoothing_length, cut_off_bnd, clip_negative_pressure) +end + +@inline function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi, + ref_system, v_ode, u_ode; + smoothing_length=ref_system.smoothing_length, + cut_off_bnd=true, clip_negative_pressure=false) num_points = length(points_coords) coords = similar(points_coords) velocities = similar(points_coords) @@ -421,14 +483,13 @@ function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi, pressures = Vector{Float64}(undef, num_points) neighbor_counts = Vector{Int}(undef, num_points) - neighborhood_searches = process_neighborhood_searches(semi, sol, ref_system, + neighborhood_searches = process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length) for (i, point) in enumerate(points_coords) - result = interpolate_point(SVector{ndims(ref_system)}(point), semi, ref_system, sol, - neighborhood_searches, smoothing_length=smoothing_length, - cut_off_bnd=cut_off_bnd, - clip_negative_pressure=clip_negative_pressure) + result = interpolate_point(SVector{ndims(ref_system)}(point), semi, ref_system, + v_ode, u_ode, neighborhood_searches; + smoothing_length, cut_off_bnd, clip_negative_pressure) densities[i] = result.density neighbor_counts[i] = result.neighbor_count coords[i] = result.coord @@ -440,28 +501,27 @@ function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi, velocity=velocities, pressure=pressures) end -function interpolate_point(point_coords, semi, ref_system, sol; +function interpolate_point(point_coords, semi, ref_system, v_ode, u_ode; smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false) - neighborhood_searches = process_neighborhood_searches(semi, sol, ref_system, + neighborhood_searches = process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length) return interpolate_point(SVector{ndims(ref_system)}(point_coords), semi, ref_system, - sol, neighborhood_searches, smoothing_length=smoothing_length, - cut_off_bnd=cut_off_bnd, - clip_negative_pressure=clip_negative_pressure) + v_ode, u_ode, neighborhood_searches; + smoothing_length, cut_off_bnd, clip_negative_pressure) end -function process_neighborhood_searches(semi, sol, ref_system, smoothing_length) +function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length) if isapprox(smoothing_length, ref_system.smoothing_length) # Update existing NHS - update_nhs(sol.u[end].x[2], semi) + update_nhs(u_ode, semi) neighborhood_searches = semi.neighborhood_searches[system_indices(ref_system, semi)] else ref_smoothing_kernel = ref_system.smoothing_kernel search_radius = compact_support(ref_smoothing_kernel, smoothing_length) neighborhood_searches = map(semi.systems) do system - u = wrap_u(sol.u[end].x[2], system, semi) + u = wrap_u(u_ode, system, semi) system_coords = current_coordinates(u, system) old_nhs = get_neighborhood_search(ref_system, system, semi) nhs = copy_neighborhood_search(old_nhs, search_radius, system_coords) @@ -472,7 +532,7 @@ function process_neighborhood_searches(semi, sol, ref_system, smoothing_length) return neighborhood_searches end -@inline function interpolate_point(point_coords, semi, ref_system, sol, +@inline function interpolate_point(point_coords, semi, ref_system, v_ode, u_ode, neighborhood_searches; smoothing_length=ref_system.smoothing_length, cut_off_bnd=true, clip_negative_pressure=false) @@ -498,8 +558,8 @@ end nhs = neighborhood_searches[system_id] (; search_radius, periodic_box) = nhs - v = wrap_v(sol.u[end].x[1], system, semi) - u = wrap_u(sol.u[end].x[2], system, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) system_coords = current_coordinates(u, system) diff --git a/test/general/interpolation.jl b/test/general/interpolation.jl index f45534f16..899e60528 100644 --- a/test/general/interpolation.jl +++ b/test/general/interpolation.jl @@ -115,16 +115,12 @@ # Density is integrated with `ContinuityDensity` v_no_bnd = vcat(fluid.velocity, fluid.density') - sol_no_boundary = (; u=[(; x=(v_no_bnd, u_no_bnd))]) - u_bnd = hcat(fluid.coordinates, bnd.coordinates) v_bnd_velocity = hcat(fluid.velocity, bnd.velocity) v_bnd_density = vcat(fluid.density, bnd.density) v_bnd = vcat(v_bnd_velocity, v_bnd_density') - sol_boundary = (; u=[(; x=(v_bnd, u_bnd))]) - semi_no_boundary = Semidiscretization(fluid_system, neighborhood_search=GridNeighborhoodSearch) semi_boundary = Semidiscretization(fluid_system, boundary_system, @@ -139,7 +135,8 @@ interpolation_walldistance(y) = TrixiParticles.interpolate_point([0.0, y], semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, + u_no_bnd, cut_off_bnd=cut_off_bnd) # top outside @@ -205,7 +202,7 @@ result_multipoint = TrixiParticles.interpolate_point(multi_point_coords, semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, u_no_bnd, cut_off_bnd=cut_off_bnd) expected_multi = (density=[666.0, 666.0000000000001, 666.0], @@ -224,13 +221,13 @@ result_endpoint = TrixiParticles.interpolate_line([1.0, -0.05], [1.0, 1.0], 5, semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, u_no_bnd, endpoint=true, cut_off_bnd=cut_off_bnd) result = TrixiParticles.interpolate_line([1.0, -0.05], [1.0, 1.0], 5, semi_no_boundary, - fluid_system, sol_no_boundary, + fluid_system, v_no_bnd, u_no_bnd, endpoint=false, cut_off_bnd=cut_off_bnd) @@ -282,7 +279,7 @@ result = interpolate_plane_2d(interpolation_start, interpolation_end, resolution, semi_no_boundary, - fluid_system, sol_no_boundary) + fluid_system, v_no_bnd, u_no_bnd) expected_res = (density=[ 666.0, 666.0, 666.0, 666.0, 666.0, 666.0, 666.0, @@ -403,7 +400,7 @@ result = interpolate_plane_2d(interpolation_start, interpolation_end, resolution, semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, u_no_bnd, smoothing_length=0.5 * smoothing_length) expected_res = (density=[ 666.0, 666.0, 666.0, 666.0, 666.0, 666.0, 666.0, 666.0, @@ -490,7 +487,8 @@ interpolation_walldistance(y) = TrixiParticles.interpolate_point([0.0, y], semi_boundary, fluid_system, - sol_boundary, + v_bnd, + u_bnd, cut_off_bnd=cut_off_bnd) # top outside @@ -578,7 +576,7 @@ result_multipoint = TrixiParticles.interpolate_point(multi_point_coords, semi_boundary, fluid_system, - sol_boundary, + v_bnd, u_bnd, cut_off_bnd=cut_off_bnd) if cut_off_bnd expected_multi = (density=[666.0, 666.0000000000001, 666.0], @@ -618,13 +616,13 @@ result_endpoint = TrixiParticles.interpolate_line([1.0, -0.05], [1.0, 1.0], 5, semi_boundary, fluid_system, - sol_boundary, + v_bnd, u_bnd, endpoint=true, cut_off_bnd=cut_off_bnd) result = TrixiParticles.interpolate_line([1.0, -0.05], [1.0, 1.0], 5, semi_no_boundary, - fluid_system, sol_no_boundary, + fluid_system, v_no_bnd, u_no_bnd, endpoint=false, cut_off_bnd=cut_off_bnd) if cut_off_bnd @@ -823,16 +821,12 @@ # Density is integrated with `ContinuityDensity` v_no_bnd = vcat(fluid.velocity, fluid.density') - sol_no_boundary = (; u=[(; x=(v_no_bnd, u_no_bnd))]) - u_bnd = hcat(fluid.coordinates, bnd.coordinates) v_bnd_velocity = hcat(fluid.velocity, bnd.velocity) v_bnd_density = vcat(fluid.density, bnd.density) v_bnd = vcat(v_bnd_velocity, v_bnd_density') - sol_boundary = (; u=[(; x=(v_bnd, u_bnd))]) - semi_no_boundary = Semidiscretization(fluid_system, neighborhood_search=GridNeighborhoodSearch) semi_boundary = Semidiscretization(fluid_system, boundary_system, @@ -851,7 +845,8 @@ ], semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, + u_no_bnd, cut_off_bnd=cut_off_bnd) # top outside @@ -920,7 +915,7 @@ result_multipoint = TrixiParticles.interpolate_point(multi_point_coords, semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, u_no_bnd, cut_off_bnd=cut_off_bnd) expected_multi = (density=[666.0, 666.0, 666.0], neighbor_count=[4, 4, 9], @@ -944,7 +939,8 @@ ], semi_boundary, fluid_system, - sol_boundary, + v_no_bnd, + u_no_bnd, cut_off_bnd=cut_off_bnd) # top outside @@ -1031,7 +1027,7 @@ result_multipoint = TrixiParticles.interpolate_point(multi_point_coords, semi_no_boundary, fluid_system, - sol_no_boundary, + v_no_bnd, u_no_bnd, cut_off_bnd=cut_off_bnd) expected_multi = (density=[666.0, 666.0, 666.0], neighbor_count=[4, 4, 9], @@ -1053,7 +1049,7 @@ result = interpolate_plane_3d(p1, p2, p3, resolution, semi_no_boundary, - fluid_system, sol_no_boundary) + fluid_system, v_no_bnd, u_no_bnd) expected_res = (density=[ 666.0, diff --git a/validation/dam_break_2d/validation_dam_break_2d.jl b/validation/dam_break_2d/validation_dam_break_2d.jl index baf2c41b4..bfe216fde 100644 --- a/validation/dam_break_2d/validation_dam_break_2d.jl +++ b/validation/dam_break_2d/validation_dam_break_2d.jl @@ -65,10 +65,9 @@ function max_x_coord(v, u, t, system) end function interpolated_pressure(coord_top, coord_bottom, v, u, t, system) - sol = (; u=[(; x=(v, u))]) n_interpolation_points = 10 interpolated_values = interpolate_line(coord_top, coord_bottom, - n_interpolation_points, semi, system, sol, + n_interpolation_points, semi, system, v, u, smoothing_length=2.0 * system.smoothing_length, clip_negative_pressure=true) return sum(map(x -> isnan(x) ? 0.0 : x, interpolated_values.pressure)) /