From afac6d8421a269901b9c2220ec767d166c05cd35 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 1 Mar 2024 16:25:00 +0100 Subject: [PATCH 01/30] add functionality for TimeSeries callback on UnstructuredMesh2D --- .../elixir_euler_basic.jl | 6 +- src/callbacks_step/time_series_dg.jl | 4 +- src/callbacks_step/time_series_dg2d.jl | 182 +++++++++++++++++- 3 files changed, 185 insertions(+), 7 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index f8976120d53..9892178c874 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -58,12 +58,16 @@ save_solution = SaveSolutionCallback(interval = 10, stepsize_callback = StepsizeCallback(cfl = 0.9) +time_series = TimeSeriesCallback(semi, [(1.4, -0.9), (2.0, -0.5), (2.8, -0.1)]; + interval = 10) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart, save_solution, - stepsize_callback) + stepsize_callback, + time_series) ############################################################################### # run the simulation diff --git a/src/callbacks_step/time_series_dg.jl b/src/callbacks_step/time_series_dg.jl index 1b63979d579..628b1e880b3 100644 --- a/src/callbacks_step/time_series_dg.jl +++ b/src/callbacks_step/time_series_dg.jl @@ -6,7 +6,9 @@ #! format: noindent # Store time series file for a TreeMesh with a DG solver -function save_time_series_file(time_series_callback, mesh::TreeMesh, equations, dg::DG) +function save_time_series_file(time_series_callback, + mesh::Union{TreeMesh, UnstructuredMesh2D}, + equations, dg::DG) @unpack (interval, solution_variables, variable_names, output_directory, filename, point_coordinates, point_data, time, step, time_series_cache) = time_series_callback diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index c15945d6e16..62a722ccd56 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -6,7 +6,9 @@ #! format: noindent # Creates cache for time series callback -function create_cache_time_series(point_coordinates, mesh::TreeMesh{2}, dg, cache) +function create_cache_time_series(point_coordinates, + mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, + dg, cache) # Determine element ids for point coordinates element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache) @@ -68,6 +70,83 @@ function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, return element_ids end +function get_elements_by_coordinates!(element_ids, coordinates, + mesh::UnstructuredMesh2D, + dg, cache) + if length(element_ids) != size(coordinates, 2) + throw(DimensionMismatch("storage length for element ids does not match the number of coordinates")) + end + + # Reset element ids - 0 indicates "not (yet) found" + element_ids .= 0 + found_elements = 0 + + # Iterate over all elements + for element in eachelement(dg, cache) + + # Iterate over coordinates + for index in 1:length(element_ids) + # Skip coordinates for which an element has already been found + if element_ids[index] > 0 + continue + end + + # Construct point + x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) + + # Skip if point is not in the current element + if !is_point_in_quad(mesh, x, element) + continue + end + + # Otherwise point is in the current element + element_ids[index] = element + found_elements += 1 + + if mesh.element_is_curved[element] && mesh.polydeg > 1 + @warn "A time series point inside a curved element may contain errors" + end + end + + # Exit loop if all elements have already been found + if found_elements == length(element_ids) + break + end + end + + return element_ids +end + +# For the `UnstructuredMesh2D` this uses a simple method assuming a convex +# quadrilateral. It simply checks that the point lies on the "correct" side +# of each of the quadrilateral's edges (basically a ray casting strategy). +# OBS! One possibility for a more robust (and maybe faster) algorithm would +# be to replace this function with `inpolygon` from PolygonOps.jl +function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) + # Helper array for the current quadrilateral element + corners = zeros(eltype(mesh.corners), 2, 4) + + # Grab the four corners + for j in 1:2, i in 1:4 + # pull the (x,y) values of these corners out of the global corners array + corners[j, i] = mesh.corners[j, mesh.element_node_ids[i, element]] + end + + if cross_product(corners[:, 2] .- corners[:, 1], point .- corners[:, 1]) > 0 && + cross_product(corners[:, 3] .- corners[:, 2], point .- corners[:, 2]) > 0 && + cross_product(corners[:, 4] .- corners[:, 3], point .- corners[:, 3]) > 0 && + cross_product(corners[:, 1] .- corners[:, 4], point .- corners[:, 4]) > 0 + return true + else + return false + end +end + +# 2D cross product +function cross_product(u, v) + return u[1] * v[2] - u[2] * v[1] +end + function get_elements_by_coordinates(coordinates, mesh, dg, cache) element_ids = Vector{Int}(undef, size(coordinates, 2)) get_elements_by_coordinates!(element_ids, coordinates, mesh, dg, cache) @@ -106,8 +185,101 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, return interpolating_polynomials end -function calc_interpolating_polynomials(coordinates, element_ids, mesh::TreeMesh, dg, - cache) +function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, + element_ids, + mesh::UnstructuredMesh2D, dg::DGSEM, cache) + @unpack nodes = dg.basis + + wbary = barycentric_weights(nodes) + + for index in 1:length(element_ids) + # Construct point + x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) + + # Convert to unit coordinates + unit_coordinates = inverse_bilinear_interpolation(mesh, x, element_ids[index]) + + println("point ", x, " has unit coordinates ", unit_coordinates, " in element ", + element_ids[index]) + + # Calculate interpolating polynomial for each dimension, making use of tensor product structure + for d in 1:ndims(mesh) + interpolating_polynomials[:, d, index] .= lagrange_interpolating_polynomials(unit_coordinates[d], + nodes, + wbary) + end + end + + return interpolating_polynomials +end + +# Given a `point` within the current convex quadrilateral `element` with +# counter-clockwise ordering of the corners +# p4 ------------ p3 +# | | +# | | +# | x `point` | +# | | +# | | +# p1 ------------ p2 +# we use a Newton iteration to determine the computational coordinates +# (xi, eta) of the `point` that is given in physical coordinates by inverting +# a bi-linear interpolation. +# The residual function for the Newton iteration is +# r = p1*(1-s)*(1-t) + p2*s*(1-t) + p3*s*t + p4*(1-s)*t - point +# and the Jacobian entries are computed accorindaly. This implementation exploits +# the 2x2 nature of the problem and directly computes the matrix inverse to make things faster. +# The implementation below is a slight modifition of that given on Stack Overflow +# https://stackoverflow.com/a/18332009 where the author explicitly states that their +# code is released to the public domain. +# Originally the bi-linear interpolant was on the interval [0,1]^2 so a final mapping +# is applied at the end to obtain the point in [-1,1]^2. +function inverse_bilinear_interpolation(mesh::UnstructuredMesh2D, point, element) + # Grab the corners of the current element + p1x = mesh.corners[1, mesh.element_node_ids[1, element]] + p1y = mesh.corners[2, mesh.element_node_ids[1, element]] + p2x = mesh.corners[1, mesh.element_node_ids[2, element]] + p2y = mesh.corners[2, mesh.element_node_ids[2, element]] + p3x = mesh.corners[1, mesh.element_node_ids[3, element]] + p3y = mesh.corners[2, mesh.element_node_ids[3, element]] + p4x = mesh.corners[1, mesh.element_node_ids[4, element]] + p4y = mesh.corners[2, mesh.element_node_ids[4, element]] + + # Initial guess for the point + s = 0.5 + t = 0.5 + for k in 1:7 # Newton's method should converge quickly + # Compute residuals for each x and y coordinate + r1 = p1x * (1 - s) * (1 - t) + p2x * s * (1 - t) + p3x * s * t + + p4x * (1 - s) * t - point[1] + r2 = p1y * (1 - s) * (1 - t) + p2y * s * (1 - t) + p3y * s * t + + p4y * (1 - s) * t - point[2] + + # Elements of the Jacobian matrix J = (dr1/ds, dr1/dt; dr2/ds, dr2/dt) + J11 = -p1x * (1 - t) + p2x * (1 - t) + p3x * t - p4x * t + J21 = -p1y * (1 - t) + p2y * (1 - t) + p3y * t - p4y * t + J12 = -p1x * (1 - s) - p2x * s + p3x * s + p4x * (1 - s) + J22 = -p1y * (1 - s) - p2y * s + p3y * s + p4y * (1 - s) + + inv_detJ = 1 / (J11 * J22 - J12 * J21) + + s = s - inv_detJ * (J22 * r1 - J12 * r2) + t = t - inv_detJ * (-J21 * r1 + J11 * r2) + + s = min(max(s, 0), 1) + t = min(max(t, 0), 1) + end + + # Map results from the interval [0,1] to the interval [-1,1] + s = 2 * s - 1 + t = 2 * t - 1 + + return SVector(s, t) +end + +function calc_interpolating_polynomials(coordinates, element_ids, + mesh::Union{TreeMesh, UnstructuredMesh2D}, + dg, cache) interpolating_polynomials = Array{real(dg), 3}(undef, nnodes(dg), ndims(mesh), length(element_ids)) @@ -121,8 +293,8 @@ end # Record the solution variables at each given point function record_state_at_points!(point_data, u, solution_variables, n_solution_variables, - mesh::TreeMesh{2}, equations, dg::DG, - time_series_cache) + mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, + equations, dg::DG, time_series_cache) @unpack element_ids, interpolating_polynomials = time_series_cache old_length = length(first(point_data)) new_length = old_length + n_solution_variables From 89ba0e865c75c6f1c4f48135aec3c02da1c76743 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 1 Mar 2024 21:04:34 +0100 Subject: [PATCH 02/30] Update src/callbacks_step/time_series_dg.jl Co-authored-by: Hendrik Ranocha --- src/callbacks_step/time_series_dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/time_series_dg.jl b/src/callbacks_step/time_series_dg.jl index 628b1e880b3..ae394afbbfd 100644 --- a/src/callbacks_step/time_series_dg.jl +++ b/src/callbacks_step/time_series_dg.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent -# Store time series file for a TreeMesh with a DG solver +# Store time series file for a DG solver function save_time_series_file(time_series_callback, mesh::Union{TreeMesh, UnstructuredMesh2D}, equations, dg::DG) From 970ddaabfa38bc6df5162841ce3743d601fa2c66 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 07:26:45 +0100 Subject: [PATCH 03/30] Apply suggestions from code review Co-authored-by: Daniel Doehring --- src/callbacks_step/time_series_dg2d.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 62a722ccd56..6fd88622381 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -132,10 +132,10 @@ function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) corners[j, i] = mesh.corners[j, mesh.element_node_ids[i, element]] end - if cross_product(corners[:, 2] .- corners[:, 1], point .- corners[:, 1]) > 0 && - cross_product(corners[:, 3] .- corners[:, 2], point .- corners[:, 2]) > 0 && - cross_product(corners[:, 4] .- corners[:, 3], point .- corners[:, 3]) > 0 && - cross_product(corners[:, 1] .- corners[:, 4], point .- corners[:, 4]) > 0 + if cross_product_2d(corners[:, 2] .- corners[:, 1], point .- corners[:, 1]) > 0 && + cross_product_2d(corners[:, 3] .- corners[:, 2], point .- corners[:, 2]) > 0 && + cross_product_2d(corners[:, 4] .- corners[:, 3], point .- corners[:, 3]) > 0 && + cross_product_2d(corners[:, 1] .- corners[:, 4], point .- corners[:, 4]) > 0 return true else return false @@ -143,7 +143,7 @@ function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) end # 2D cross product -function cross_product(u, v) +@inline function cross_product_2d(u, v) return u[1] * v[2] - u[2] * v[1] end @@ -245,7 +245,7 @@ function inverse_bilinear_interpolation(mesh::UnstructuredMesh2D, point, element p4x = mesh.corners[1, mesh.element_node_ids[4, element]] p4y = mesh.corners[2, mesh.element_node_ids[4, element]] - # Initial guess for the point + # Initial guess for the point (center of the element) s = 0.5 t = 0.5 for k in 1:7 # Newton's method should converge quickly From ef91125d60e840f2fdbcc2a82f1b20b761e15d12 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 11:47:16 +0100 Subject: [PATCH 04/30] add strategy to correctly locate a gauge point within a curvilinear element --- src/callbacks_step/time_series_dg2d.jl | 142 ++++++++++++++----------- 1 file changed, 81 insertions(+), 61 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 6fd88622381..f95f3b65a0c 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -120,9 +120,10 @@ end # For the `UnstructuredMesh2D` this uses a simple method assuming a convex # quadrilateral. It simply checks that the point lies on the "correct" side # of each of the quadrilateral's edges (basically a ray casting strategy). +# Does not account for element curvature. # OBS! One possibility for a more robust (and maybe faster) algorithm would # be to replace this function with `inpolygon` from PolygonOps.jl -function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) +@inline function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) # Helper array for the current quadrilateral element corners = zeros(eltype(mesh.corners), 2, 4) @@ -192,13 +193,27 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, wbary = barycentric_weights(nodes) + # Helper array for a straight-sided quadrilateral element + corners = zeros(eltype(mesh.corners), 4, 2) + for index in 1:length(element_ids) # Construct point x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) - # Convert to unit coordinates - unit_coordinates = inverse_bilinear_interpolation(mesh, x, element_ids[index]) + # Convert to unit coordinates; procedure differs for straight-sided + # versus curvilinear elements + element = element_ids[index] + if !mesh.element_is_curved[element] + for j in 1:2, i in 1:4 + # Pull the (x,y) values of the element corners from the global corners array + corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]] + end + unit_coordinates = invert_bilinear_interpolation(mesh, x, corners) + else # mesh.element_is_curved[element] + unit_coordinates = invert_transfinite_interpolation(mesh, x, view(mesh.surface_curves, :, element)) + end + # TODO: debug statment for removal println("point ", x, " has unit coordinates ", unit_coordinates, " in element ", element_ids[index]) @@ -213,68 +228,73 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, return interpolating_polynomials end -# Given a `point` within the current convex quadrilateral `element` with -# counter-clockwise ordering of the corners -# p4 ------------ p3 -# | | -# | | -# | x `point` | -# | | -# | | -# p1 ------------ p2 -# we use a Newton iteration to determine the computational coordinates -# (xi, eta) of the `point` that is given in physical coordinates by inverting -# a bi-linear interpolation. +# Use a Newton iteration to determine the computational coordinates +# (xi, eta) of given (x,y) `point` that is given in physical coordinates +# by inverting the transformation. For straight-sided elements this +# amounts to inverting a bi-linear interpolation. For curved +# elements we invert the transfinite interpolation with linear blending. # The residual function for the Newton iteration is -# r = p1*(1-s)*(1-t) + p2*s*(1-t) + p3*s*t + p4*(1-s)*t - point -# and the Jacobian entries are computed accorindaly. This implementation exploits -# the 2x2 nature of the problem and directly computes the matrix inverse to make things faster. -# The implementation below is a slight modifition of that given on Stack Overflow -# https://stackoverflow.com/a/18332009 where the author explicitly states that their -# code is released to the public domain. -# Originally the bi-linear interpolant was on the interval [0,1]^2 so a final mapping -# is applied at the end to obtain the point in [-1,1]^2. -function inverse_bilinear_interpolation(mesh::UnstructuredMesh2D, point, element) - # Grab the corners of the current element - p1x = mesh.corners[1, mesh.element_node_ids[1, element]] - p1y = mesh.corners[2, mesh.element_node_ids[1, element]] - p2x = mesh.corners[1, mesh.element_node_ids[2, element]] - p2y = mesh.corners[2, mesh.element_node_ids[2, element]] - p3x = mesh.corners[1, mesh.element_node_ids[3, element]] - p3y = mesh.corners[2, mesh.element_node_ids[3, element]] - p4x = mesh.corners[1, mesh.element_node_ids[4, element]] - p4y = mesh.corners[2, mesh.element_node_ids[4, element]] - - # Initial guess for the point (center of the element) - s = 0.5 - t = 0.5 - for k in 1:7 # Newton's method should converge quickly - # Compute residuals for each x and y coordinate - r1 = p1x * (1 - s) * (1 - t) + p2x * s * (1 - t) + p3x * s * t + - p4x * (1 - s) * t - point[1] - r2 = p1y * (1 - s) * (1 - t) + p2y * s * (1 - t) + p3y * s * t + - p4y * (1 - s) * t - point[2] - - # Elements of the Jacobian matrix J = (dr1/ds, dr1/dt; dr2/ds, dr2/dt) - J11 = -p1x * (1 - t) + p2x * (1 - t) + p3x * t - p4x * t - J21 = -p1y * (1 - t) + p2y * (1 - t) + p3y * t - p4y * t - J12 = -p1x * (1 - s) - p2x * s + p3x * s + p4x * (1 - s) - J22 = -p1y * (1 - s) - p2y * s + p3y * s + p4y * (1 - s) - - inv_detJ = 1 / (J11 * J22 - J12 * J21) - - s = s - inv_detJ * (J22 * r1 - J12 * r2) - t = t - inv_detJ * (-J21 * r1 + J11 * r2) - - s = min(max(s, 0), 1) - t = min(max(t, 0), 1) +# r(xi,eta) = X(xi,eta) - point +# and the Jacobian entries are computed accordingly from either +# `straight_side_quad_map_metrics` or `transfinite_quad_map_metrics`. +# We exploit the 2x2 nature of the problem and directly compute the matrix +# inverse to make things faster. The implementations below are inspired by +# an answer on Stack Overflow (https://stackoverflow.com/a/18332009) where +# the author explicitly states that their code is released to the public domain. +@inline function invert_bilinear_interpolation(mesh::UnstructuredMesh2D, point, element_corners) + # Initial guess for the point (center of the reference element) + xi = zero(eltype(point)) + eta = zero(eltype(point)) + for k in 1:5 # Newton's method should converge quickly + # Compute current x and y coordinate and the Jacobian matrix + # J = (X_xi, X_eta; Y_xi, Y_eta) + x, y = straight_side_quad_map(xi, eta, element_corners) + J11, J12, J21, J22 = straight_side_quad_map_metrics(xi, eta, element_corners) + + # Compute residuals for the Newton teration for the current (x, y) coordinate + r1 = x - point[1] + r2 = y - point[2] + + # Newton update that directly applies the inverse of the 2x2 Jacobian matrix + inv_detJ = inv(J11 * J22 - J12 * J21) + + xi = xi - inv_detJ * (J22 * r1 - J12 * r2) + eta = eta - inv_detJ * (-J21 * r1 + J11 * r2) + + # Ensure updated point is in the reference element + xi = min(max(xi, -1), 1) + eta = min(max(eta, -1), 1) end - # Map results from the interval [0,1] to the interval [-1,1] - s = 2 * s - 1 - t = 2 * t - 1 + return SVector(xi, eta) +end + +@inline function invert_transfinite_interpolation(mesh::UnstructuredMesh2D, point, surface_curves::AbstractVector{<:CurvedSurface}) + # Initial guess for the point (center of the reference element) + xi = zero(eltype(point)) + eta = zero(eltype(point)) + for k in 1:5 # Newton's method should converge quickly + # Compute current x and y coordinate and the Jacobian matrix + # J = (X_xi, X_eta; Y_xi, Y_eta) + x, y = transfinite_quad_map(xi, eta, surface_curves) + J11, J12, J21, J22 = transfinite_quad_map_metrics(xi, eta, surface_curves) + + # Compute residuals for the Newton teration for the current (x,y) coordinate + r1 = x - point[1] + r2 = y - point[2] + + # Newton update that directly applies the inverse of the 2x2 Jacobian matrix + inv_detJ = inv(J11 * J22 - J12 * J21) + + xi = xi - inv_detJ * (J22 * r1 - J12 * r2) + eta = eta - inv_detJ * (-J21 * r1 + J11 * r2) + + # Ensure updated point is in the reference element + xi = min(max(xi, -1), 1) + eta = min(max(eta, -1), 1) + end - return SVector(s, t) + return SVector(xi, eta) end function calc_interpolating_polynomials(coordinates, element_ids, From b6840481c15a26bf499ae8e7a0b18b065ebee495 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 12:02:38 +0100 Subject: [PATCH 05/30] add sanity check that the Newton solution is correct --- src/callbacks_step/time_series_dg2d.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index f95f3b65a0c..aa24d112a6d 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -102,10 +102,6 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Otherwise point is in the current element element_ids[index] = element found_elements += 1 - - if mesh.element_is_curved[element] && mesh.polydeg > 1 - @warn "A time series point inside a curved element may contain errors" - end end # Exit loop if all elements have already been found @@ -209,8 +205,20 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]] end unit_coordinates = invert_bilinear_interpolation(mesh, x, corners) + + # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` + x_check = straight_side_quad_map(unit_coordinates[1], unit_coordinates[2], corners) + if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) + error("failed to compute computational coordinates for the time series point $(x)") + end else # mesh.element_is_curved[element] unit_coordinates = invert_transfinite_interpolation(mesh, x, view(mesh.surface_curves, :, element)) + + # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` + x_check = transfinite_quad_map(unit_coordinates[1], unit_coordinates[2], view(mesh.surface_curves, :, element)) + if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) + error("failed to compute computational coordinates for the time series point $(x)") + end end # TODO: debug statment for removal From bcbcffcebc8b22477df9045a6aebe58fe85edda0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 12:04:54 +0100 Subject: [PATCH 06/30] run formatter --- src/callbacks_step/time_series_dg2d.jl | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index aa24d112a6d..e4ec49ab29e 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -207,16 +207,22 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, unit_coordinates = invert_bilinear_interpolation(mesh, x, corners) # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` - x_check = straight_side_quad_map(unit_coordinates[1], unit_coordinates[2], corners) - if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) + x_check = straight_side_quad_map(unit_coordinates[1], unit_coordinates[2], + corners) + if !isapprox(x[1], x_check[1], atol = 1e-13) || + !isapprox(x[2], x_check[2], atol = 1e-13) error("failed to compute computational coordinates for the time series point $(x)") end else # mesh.element_is_curved[element] - unit_coordinates = invert_transfinite_interpolation(mesh, x, view(mesh.surface_curves, :, element)) + unit_coordinates = invert_transfinite_interpolation(mesh, x, + view(mesh.surface_curves, + :, element)) # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` - x_check = transfinite_quad_map(unit_coordinates[1], unit_coordinates[2], view(mesh.surface_curves, :, element)) - if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) + x_check = transfinite_quad_map(unit_coordinates[1], unit_coordinates[2], + view(mesh.surface_curves, :, element)) + if !isapprox(x[1], x_check[1], atol = 1e-13) || + !isapprox(x[2], x_check[2], atol = 1e-13) error("failed to compute computational coordinates for the time series point $(x)") end end @@ -249,7 +255,8 @@ end # inverse to make things faster. The implementations below are inspired by # an answer on Stack Overflow (https://stackoverflow.com/a/18332009) where # the author explicitly states that their code is released to the public domain. -@inline function invert_bilinear_interpolation(mesh::UnstructuredMesh2D, point, element_corners) +@inline function invert_bilinear_interpolation(mesh::UnstructuredMesh2D, point, + element_corners) # Initial guess for the point (center of the reference element) xi = zero(eltype(point)) eta = zero(eltype(point)) @@ -277,7 +284,8 @@ end return SVector(xi, eta) end -@inline function invert_transfinite_interpolation(mesh::UnstructuredMesh2D, point, surface_curves::AbstractVector{<:CurvedSurface}) +@inline function invert_transfinite_interpolation(mesh::UnstructuredMesh2D, point, + surface_curves::AbstractVector{<:CurvedSurface}) # Initial guess for the point (center of the reference element) xi = zero(eltype(point)) eta = zero(eltype(point)) From 25071a35590de73aa1b0b7096653d61dd99f5837 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 19:15:01 +0100 Subject: [PATCH 07/30] implement a more general approach that also works on curved element without issue --- .../elixir_euler_basic.jl | 2 +- src/callbacks_step/time_series_dg2d.jl | 84 ++++++------------- 2 files changed, 26 insertions(+), 60 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index 9892178c874..191f5540ab5 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -58,7 +58,7 @@ save_solution = SaveSolutionCallback(interval = 10, stepsize_callback = StepsizeCallback(cfl = 0.9) -time_series = TimeSeriesCallback(semi, [(1.4, -0.9), (2.0, -0.5), (2.8, -0.1)]; +time_series = TimeSeriesCallback(semi, [(2.0, -0.5), (0.28, 1/3), (1.87, 2/3)]; interval = 10) callbacks = CallbackSet(summary_callback, diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index e4ec49ab29e..52e3afc0fba 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -70,6 +70,11 @@ function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, return element_ids end +# Elements on an `UnstructuredMesh2D` are possibly curved. Assume that each +# element is convex, i.e., all interior angles are less than 180 degrees. +# We use the barycenter of each element to determine if a given coordinate (x,y) +# lies within said element. The shortest distance between the point and all the +# barycenters "wins". function get_elements_by_coordinates!(element_ids, coordinates, mesh::UnstructuredMesh2D, dg, cache) @@ -79,69 +84,34 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Reset element ids - 0 indicates "not (yet) found" element_ids .= 0 - found_elements = 0 - # Iterate over all elements - for element in eachelement(dg, cache) + # Compute and save the barycentric coordinate on each element + bary_centers = zeros(eltype(mesh.corners), 2, mesh.n_elements) + distances = zeros(eltype(mesh.corners), mesh.n_elements) + calc_bary_centers!(bary_centers, dg, cache) - # Iterate over coordinates - for index in 1:length(element_ids) - # Skip coordinates for which an element has already been found - if element_ids[index] > 0 - continue - end - - # Construct point - x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) - - # Skip if point is not in the current element - if !is_point_in_quad(mesh, x, element) - continue - end - - # Otherwise point is in the current element - element_ids[index] = element - found_elements += 1 + # Iterate over coordinates + for index in 1:length(element_ids) + # Compute the distance between the current point and all the bary centers. + for element in eachelement(dg, cache) + distances[element] = norm(coordinates[:, index] .- bary_centers[:, element]) end - # Exit loop if all elements have already been found - if found_elements == length(element_ids) - break - end + # Locate the minimal distance as this gives the element containing the point + element_ids[index] = argmin(distances) end return element_ids end -# For the `UnstructuredMesh2D` this uses a simple method assuming a convex -# quadrilateral. It simply checks that the point lies on the "correct" side -# of each of the quadrilateral's edges (basically a ray casting strategy). -# Does not account for element curvature. -# OBS! One possibility for a more robust (and maybe faster) algorithm would -# be to replace this function with `inpolygon` from PolygonOps.jl -@inline function is_point_in_quad(mesh::UnstructuredMesh2D, point, element) - # Helper array for the current quadrilateral element - corners = zeros(eltype(mesh.corners), 2, 4) - - # Grab the four corners - for j in 1:2, i in 1:4 - # pull the (x,y) values of these corners out of the global corners array - corners[j, i] = mesh.corners[j, mesh.element_node_ids[i, element]] - end - - if cross_product_2d(corners[:, 2] .- corners[:, 1], point .- corners[:, 1]) > 0 && - cross_product_2d(corners[:, 3] .- corners[:, 2], point .- corners[:, 2]) > 0 && - cross_product_2d(corners[:, 4] .- corners[:, 3], point .- corners[:, 3]) > 0 && - cross_product_2d(corners[:, 1] .- corners[:, 4], point .- corners[:, 4]) > 0 - return true - else - return false +# Use the avaiable `node_coordinates` on each element to compute and save the barycenter. +@inline function calc_bary_centers!(bary_centers, dg, cache) + n = nnodes(dg) + for element in eachelement(dg, cache) + bary_centers[1, element] = sum(cache.elements.node_coordinates[1,:,:,element]) / n^2 + bary_centers[2, element] = sum(cache.elements.node_coordinates[2,:,:,element]) / n^2 end -end - -# 2D cross product -@inline function cross_product_2d(u, v) - return u[1] * v[2] - u[2] * v[1] + return nothing end function get_elements_by_coordinates(coordinates, mesh, dg, cache) @@ -211,7 +181,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, corners) if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x)") + error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") end else # mesh.element_is_curved[element] unit_coordinates = invert_transfinite_interpolation(mesh, x, @@ -223,14 +193,10 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, view(mesh.surface_curves, :, element)) if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x)") + error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") end end - # TODO: debug statment for removal - println("point ", x, " has unit coordinates ", unit_coordinates, " in element ", - element_ids[index]) - # Calculate interpolating polynomial for each dimension, making use of tensor product structure for d in 1:ndims(mesh) interpolating_polynomials[:, d, index] .= lagrange_interpolating_polynomials(unit_coordinates[d], From 735ad18a4d81b408c044cdf9d57962952b0d2030 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 19:22:29 +0100 Subject: [PATCH 08/30] run formatter --- src/callbacks_step/time_series_dg2d.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 52e3afc0fba..ee6c0ad89c7 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -94,7 +94,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, for index in 1:length(element_ids) # Compute the distance between the current point and all the bary centers. for element in eachelement(dg, cache) - distances[element] = norm(coordinates[:, index] .- bary_centers[:, element]) + distances[element] = norm(coordinates[:, index] .- bary_centers[:, element]) end # Locate the minimal distance as this gives the element containing the point @@ -104,12 +104,14 @@ function get_elements_by_coordinates!(element_ids, coordinates, return element_ids end -# Use the avaiable `node_coordinates` on each element to compute and save the barycenter. +# Use the available `node_coordinates` on each element to compute and save the barycenter. @inline function calc_bary_centers!(bary_centers, dg, cache) n = nnodes(dg) for element in eachelement(dg, cache) - bary_centers[1, element] = sum(cache.elements.node_coordinates[1,:,:,element]) / n^2 - bary_centers[2, element] = sum(cache.elements.node_coordinates[2,:,:,element]) / n^2 + bary_centers[1, element] = sum(cache.elements.node_coordinates[1, :, :, + element]) / n^2 + bary_centers[2, element] = sum(cache.elements.node_coordinates[2, :, :, + element]) / n^2 end return nothing end From 935d01edbf077d240ffac6c52690c874eba74814 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 7 Mar 2024 19:24:59 +0100 Subject: [PATCH 09/30] forgot to format the examples --- examples/unstructured_2d_dgsem/elixir_euler_basic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index 191f5540ab5..92a33f64371 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -58,7 +58,7 @@ save_solution = SaveSolutionCallback(interval = 10, stepsize_callback = StepsizeCallback(cfl = 0.9) -time_series = TimeSeriesCallback(semi, [(2.0, -0.5), (0.28, 1/3), (1.87, 2/3)]; +time_series = TimeSeriesCallback(semi, [(2.0, -0.5), (0.28, 1 / 3), (1.87, 2 / 3)]; interval = 10) callbacks = CallbackSet(summary_callback, From b8b99dbc1077290ef87ae7d43e5ad09e4c8c79e4 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 8 Mar 2024 09:50:42 +0100 Subject: [PATCH 10/30] Apply suggestions from code review Co-authored-by: Hendrik Ranocha Co-authored-by: Daniel Doehring --- .../elixir_euler_basic.jl | 4 +-- src/callbacks_step/time_series_dg2d.jl | 29 ++++++++++++------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index 92a33f64371..4ebf5e29780 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -66,8 +66,8 @@ callbacks = CallbackSet(summary_callback, alive_callback, save_restart, save_solution, - stepsize_callback, - time_series) + time_series, + stepsize_callback) ############################################################################### # run the simulation diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index ee6c0ad89c7..b72575b2aed 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -87,18 +87,24 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Compute and save the barycentric coordinate on each element bary_centers = zeros(eltype(mesh.corners), 2, mesh.n_elements) - distances = zeros(eltype(mesh.corners), mesh.n_elements) calc_bary_centers!(bary_centers, dg, cache) # Iterate over coordinates for index in 1:length(element_ids) - # Compute the distance between the current point and all the bary centers. - for element in eachelement(dg, cache) - distances[element] = norm(coordinates[:, index] .- bary_centers[:, element]) + point = SVector(coordinates[1, index], + coordinates[2, index]) + + # Search for the element with minimal distance between the point + # we are interested in and the barycenter of the element. + element_ids[index] = argmin(eachelement(dg, cache)) do element + bary_center = SVector(bary_centers[1, element], + bary_centers[2, element]) + # Compute the squared Euclidean distance between the point we are + # interested in and the barycenter of the element. + # We compute the squared norm to avoid computing computationally + # more expensive square roots. + sum(abs2, point - bary_center) end - - # Locate the minimal distance as this gives the element containing the point - element_ids[index] = argmin(distances) end return element_ids @@ -107,7 +113,7 @@ end # Use the available `node_coordinates` on each element to compute and save the barycenter. @inline function calc_bary_centers!(bary_centers, dg, cache) n = nnodes(dg) - for element in eachelement(dg, cache) + @views for element in eachelement(dg, cache) bary_centers[1, element] = sum(cache.elements.node_coordinates[1, :, :, element]) / n^2 bary_centers[2, element] = sum(cache.elements.node_coordinates[2, :, :, @@ -176,6 +182,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, # Pull the (x,y) values of the element corners from the global corners array corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]] end + # Compute coordinates in reference system unit_coordinates = invert_bilinear_interpolation(mesh, x, corners) # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` @@ -183,7 +190,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, corners) if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") + error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)") end else # mesh.element_is_curved[element] unit_coordinates = invert_transfinite_interpolation(mesh, x, @@ -195,7 +202,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, view(mesh.surface_curves, :, element)) if !isapprox(x[1], x_check[1], atol = 1e-13) || !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") + error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)") end end @@ -241,6 +248,7 @@ end # Newton update that directly applies the inverse of the 2x2 Jacobian matrix inv_detJ = inv(J11 * J22 - J12 * J21) + # Update with explicitly inverted Jacobian xi = xi - inv_detJ * (J22 * r1 - J12 * r2) eta = eta - inv_detJ * (-J21 * r1 + J11 * r2) @@ -270,6 +278,7 @@ end # Newton update that directly applies the inverse of the 2x2 Jacobian matrix inv_detJ = inv(J11 * J22 - J12 * J21) + # Update with explicitly inverted Jacobian xi = xi - inv_detJ * (J22 * r1 - J12 * r2) eta = eta - inv_detJ * (-J21 * r1 + J11 * r2) From 818edbfcaeca177524d1b3e3268a6c82bfbf444e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 8 Mar 2024 23:19:42 +0100 Subject: [PATCH 11/30] working version of the element finding routine --- src/callbacks_step/time_series_dg2d.jl | 111 ++++++++++++++++++++----- 1 file changed, 88 insertions(+), 23 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index b72575b2aed..09029af1b88 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -72,9 +72,13 @@ end # Elements on an `UnstructuredMesh2D` are possibly curved. Assume that each # element is convex, i.e., all interior angles are less than 180 degrees. -# We use the barycenter of each element to determine if a given coordinate (x,y) -# lies within said element. The shortest distance between the point and all the -# barycenters "wins". +# This routine computes the shortest distance from a given point to each element +# surface in the mesh. These distances then indicate possible candidate elements. +# From these candidates we (essentially) apply a ray casting strategy and identify +# the element in which the point lies by comparing the ray formed by the point to +# the nearest boundary to the rays cast by the candidate element barycenters to the +# boundary. If these rays point in the same direction, then we have identified the +# desired element location. function get_elements_by_coordinates!(element_ids, coordinates, mesh::UnstructuredMesh2D, dg, cache) @@ -90,20 +94,50 @@ function get_elements_by_coordinates!(element_ids, coordinates, calc_bary_centers!(bary_centers, dg, cache) # Iterate over coordinates + distances = zeros(eltype(mesh.corners), mesh.n_elements) + indices = zeros(Int, mesh.n_elements, 2) for index in 1:length(element_ids) - point = SVector(coordinates[1, index], + # Grab the current point for which the element needs found + point = SVector(coordinates[1, index], coordinates[2, index]) - - # Search for the element with minimal distance between the point - # we are interested in and the barycenter of the element. - element_ids[index] = argmin(eachelement(dg, cache)) do element - bary_center = SVector(bary_centers[1, element], - bary_centers[2, element]) - # Compute the squared Euclidean distance between the point we are - # interested in and the barycenter of the element. - # We compute the squared norm to avoid computing computationally - # more expensive square roots. - sum(abs2, point - bary_center) + + # Compute the minimum distance between the `point` and all the element surfaces + # saved into `distances`. The point in `node_coordinates` that gives said minimum + # distance on each element is saved in `indices` + distances, indices = calc_minimum_surface_distance(point, + cache.elements.node_coordinates, + dg, mesh) + + # Get the candidate elements where the `point` might live + candidates = findall( abs.(minimum(distances) .- distances) .< 500 * eps(eltype(point)) ) + + # The minimal surface point is on a boundary so it plays no role which candidate + # we use to grab it. So just use the first one + surface_point = SVector(cache.elements.node_coordinates[1, indices[candidates[1], 1], indices[candidates[1], 2], candidates[1]], + cache.elements.node_coordinates[2, indices[candidates[1], 1], indices[candidates[1], 2], candidates[1]]) + + # Compute the vector pointing from the current `point` toward the surface + P = surface_point .- point + + # If the vector `P` is the zero vector then this `point` is at an element corner or + # on a surface. So any candidate element works, just use the first one + if sum(P .* P) < 500 * eps(eltype(point)) + element_ids[index] = candidates[1] + continue + end + + # Loop through all the element candidates until we find a vector from the barycenter + # to the surface that points in the same direction as the current `point` vector. + # This then gives us the correct element. + for element = 1:length(candidates) + bary_center = SVector(bary_centers[1, candidates[element]], + bary_centers[2, candidates[element]]) + # Vector pointing from the barycenter toward the minimal `surface_point` + B = surface_point .- bary_center + if sum(P .* B) > zero(eltype(bary_center)) + element_ids[index] = candidates[element] + break + end end end @@ -122,6 +156,39 @@ end return nothing end +# Compute the shortest distance from a `point` to the surface of each element +# using the available `node_coordinates`. Also return the index pair of this +# minimum surface point location. We compute the squared norm to avoid +# computing computationally more expensive square roots. +# OBS! Could be made more accuracte if the `node_coordinates` were super-sampled +# and reinterpolated onto a higher polynomial degree before this computation. +function calc_minimum_surface_distance(point, node_coordinates, + dg, mesh::UnstructuredMesh2D) + n = nnodes(dg) + hmin2 = Inf * ones(eltype(mesh.corners), length(mesh)) + indices = zeros(Int, length(mesh), 2) + for k = 1:length(mesh) + # used to ensure that only boundary points are used + on_surface = [false, false] + for j = 1:n + on_surface[2] = (j == 1) || (j == n) + for i = 1:n + on_surface[1] = (i == 1) || (i == n) + if !any(on_surface) + continue + end + if sum((node_coordinates[:, i, j, k] - point) .* (node_coordinates[:, i, j, k] - point)) < hmin2[k] + hmin2[k] = sum((node_coordinates[:, i, j, k] - point) .* (node_coordinates[:, i, j, k] - point)) + indices[k, 1] = i + indices[k, 2] = j + end + end + end + end + + return hmin2, indices + end + function get_elements_by_coordinates(coordinates, mesh, dg, cache) element_ids = Vector{Int}(undef, size(coordinates, 2)) get_elements_by_coordinates!(element_ids, coordinates, mesh, dg, cache) @@ -188,9 +255,8 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` x_check = straight_side_quad_map(unit_coordinates[1], unit_coordinates[2], corners) - if !isapprox(x[1], x_check[1], atol = 1e-13) || - !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)") + if !isapprox(x[1], x_check[1]) || !isapprox(x[2], x_check[2]) + error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") end else # mesh.element_is_curved[element] unit_coordinates = invert_transfinite_interpolation(mesh, x, @@ -200,9 +266,8 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x` x_check = transfinite_quad_map(unit_coordinates[1], unit_coordinates[2], view(mesh.surface_curves, :, element)) - if !isapprox(x[1], x_check[1], atol = 1e-13) || - !isapprox(x[2], x_check[2], atol = 1e-13) - error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)") + if !isapprox(x[1], x_check[1]) || !isapprox(x[2], x_check[2]) + error("failed to compute computational coordinates for the time series point $(x), closet candidate was $(x_check)") end end @@ -218,12 +283,12 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, end # Use a Newton iteration to determine the computational coordinates -# (xi, eta) of given (x,y) `point` that is given in physical coordinates +# (xi, eta) of an (x,y) `point` that is given in physical coordinates # by inverting the transformation. For straight-sided elements this # amounts to inverting a bi-linear interpolation. For curved # elements we invert the transfinite interpolation with linear blending. # The residual function for the Newton iteration is -# r(xi,eta) = X(xi,eta) - point +# r(xi, eta) = X(xi, eta) - point # and the Jacobian entries are computed accordingly from either # `straight_side_quad_map_metrics` or `transfinite_quad_map_metrics`. # We exploit the 2x2 nature of the problem and directly compute the matrix From 6f5e24e0fc8adfa4a88ddd85d6810f6068939ee1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 8 Mar 2024 23:22:11 +0100 Subject: [PATCH 12/30] run formatter --- src/callbacks_step/time_series_dg2d.jl | 35 ++++++++++++++++++-------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 09029af1b88..0e554daf836 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -109,12 +109,23 @@ function get_elements_by_coordinates!(element_ids, coordinates, dg, mesh) # Get the candidate elements where the `point` might live - candidates = findall( abs.(minimum(distances) .- distances) .< 500 * eps(eltype(point)) ) + candidates = findall(abs.(minimum(distances) .- distances) .< + 500 * eps(eltype(point))) # The minimal surface point is on a boundary so it plays no role which candidate # we use to grab it. So just use the first one - surface_point = SVector(cache.elements.node_coordinates[1, indices[candidates[1], 1], indices[candidates[1], 2], candidates[1]], - cache.elements.node_coordinates[2, indices[candidates[1], 1], indices[candidates[1], 2], candidates[1]]) + surface_point = SVector(cache.elements.node_coordinates[1, + indices[candidates[1], + 1], + indices[candidates[1], + 2], + candidates[1]], + cache.elements.node_coordinates[2, + indices[candidates[1], + 1], + indices[candidates[1], + 2], + candidates[1]]) # Compute the vector pointing from the current `point` toward the surface P = surface_point .- point @@ -129,7 +140,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Loop through all the element candidates until we find a vector from the barycenter # to the surface that points in the same direction as the current `point` vector. # This then gives us the correct element. - for element = 1:length(candidates) + for element in 1:length(candidates) bary_center = SVector(bary_centers[1, candidates[element]], bary_centers[2, candidates[element]]) # Vector pointing from the barycenter toward the minimal `surface_point` @@ -163,22 +174,24 @@ end # OBS! Could be made more accuracte if the `node_coordinates` were super-sampled # and reinterpolated onto a higher polynomial degree before this computation. function calc_minimum_surface_distance(point, node_coordinates, - dg, mesh::UnstructuredMesh2D) + dg, mesh::UnstructuredMesh2D) n = nnodes(dg) hmin2 = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) - for k = 1:length(mesh) + for k in 1:length(mesh) # used to ensure that only boundary points are used on_surface = [false, false] - for j = 1:n + for j in 1:n on_surface[2] = (j == 1) || (j == n) - for i = 1:n + for i in 1:n on_surface[1] = (i == 1) || (i == n) if !any(on_surface) continue end - if sum((node_coordinates[:, i, j, k] - point) .* (node_coordinates[:, i, j, k] - point)) < hmin2[k] - hmin2[k] = sum((node_coordinates[:, i, j, k] - point) .* (node_coordinates[:, i, j, k] - point)) + if sum((node_coordinates[:, i, j, k] - point) .* + (node_coordinates[:, i, j, k] - point)) < hmin2[k] + hmin2[k] = sum((node_coordinates[:, i, j, k] - point) .* + (node_coordinates[:, i, j, k] - point)) indices[k, 1] = i indices[k, 2] = j end @@ -187,7 +200,7 @@ function calc_minimum_surface_distance(point, node_coordinates, end return hmin2, indices - end +end function get_elements_by_coordinates(coordinates, mesh, dg, cache) element_ids = Vector{Int}(undef, size(coordinates, 2)) From 00d9b836f0cfb24672e243253b8d7c2a7549731e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 9 Mar 2024 00:04:27 +0100 Subject: [PATCH 13/30] add new elixir for the time series callback --- .../elixir_euler_basic.jl | 4 - .../elixir_euler_time_series.jl | 116 ++++++++++++++++++ 2 files changed, 116 insertions(+), 4 deletions(-) create mode 100644 examples/unstructured_2d_dgsem/elixir_euler_time_series.jl diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index 4ebf5e29780..f8976120d53 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -58,15 +58,11 @@ save_solution = SaveSolutionCallback(interval = 10, stepsize_callback = StepsizeCallback(cfl = 0.9) -time_series = TimeSeriesCallback(semi, [(2.0, -0.5), (0.28, 1 / 3), (1.87, 2 / 3)]; - interval = 10) - callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart, save_solution, - time_series, stepsize_callback) ############################################################################### diff --git a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl new file mode 100644 index 00000000000..c8b930fc64a --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl @@ -0,0 +1,116 @@ +# An elixir that has an alternative convergence test that uses +# the `TimeSeriesCallback` on several gauge points. Many of the +# gauge points are selected as "stress tests" for the element +# identification, e.g., a gauge point that lies on an +# element corner of a curvlinear mesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +# Modify the manufactured solution test to use `L = sqrt(2)` +# in the initial condition and source terms +function initial_condition_convergence_shifted(x, t, + equations::CompressibleEulerEquations2D) + c = 2 + A = 0.1 + L = sqrt(2) + f = 1 / L + ω = 2 * pi * f + ini = c + A * sin(ω * (x[1] + x[2] - t)) + + rho = ini + rho_v1 = ini + rho_v2 = ini + rho_e = ini^2 + + return SVector(rho, rho_v1, rho_v2, rho_e) +end + +@inline function source_terms_convergence_shifted(u, x, t, + equations::CompressibleEulerEquations2D) + # Same settings as in `initial_condition` + c = 2 + A = 0.1 + L = sqrt(2) + f = 1 / L + ω = 2 * pi * f + γ = equations.gamma + + x1, x2 = x + si, co = sincos(ω * (x1 + x2 - t)) + rho = c + A * si + rho_x = ω * A * co + # Note that d/dt rho = -d/dx rho = -d/dy rho. + + tmp = (2 * rho - 1) * (γ - 1) + + du1 = rho_x + du2 = rho_x * (1 + tmp) + du3 = du2 + du4 = 2 * rho_x * (rho + tmp) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_convergence_shifted + +source_term = source_terms_convergence_shifted + +############################################################################### +# Get the DG approximation space + +solver = DGSEM(polydeg = 6, surface_flux = flux_lax_friedrichs) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh", + joinpath(@__DIR__, "mesh_multiple_flips.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_term) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +time_series = TimeSeriesCallback(semi, + [(0.75, 0.7), (1.23, 0.302), (0.8, 1.0), + (0.353553390593274, 0.353553390593274), + (0.505, 1.125), (1.37, 0.89), (0.349, 0.7153), + (0.883883476483184, 0.406586401289607), + (sqrt(2), sqrt(2))]; + interval = 10) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + time_series, + alive_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-6, reltol = 1.0e-6, + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary From 7f52beda44e1db482a8a9aa5d9d8452af8b4b549 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 9 Mar 2024 00:13:38 +0100 Subject: [PATCH 14/30] add additional test for the time series callback on an unstructured mesh --- test/test_unstructured_2d.jl | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 8a62dcaec3c..0f72ce03818 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -198,6 +198,36 @@ end end end +@trixi_testset "elixir_euler_time_series.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_time_series.jl"), + l2=[ + 9.44589994864801e-5, + 7.375567345161402e-5, + 8.469037537040469e-5, + 0.0001638067441092813, + ], + linf=[ + 0.0006797328621439558, + 0.00048723709312858965, + 0.0006701229224337357, + 0.0011917195110209278, + ], + tspan=(0.0, 0.5)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + # Extra tests to make sure the `TimeSeriesCallback` made correct data + point_data_1 = time_series.affect!.point_data[1] + @test all(isapprox.(point_data_1[1:4], + [1.9548629504179071, 1.9548895017660788, + 1.954889292850916, 3.8217607623030623])) +end + @trixi_testset "elixir_acoustics_gauss_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss_wall.jl"), l2=[0.029330394861252995, 0.029345079728907965, From 76d0dcb28fd0f1e6832e44cdcade15c977b5daa0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 9 Mar 2024 00:20:26 +0100 Subject: [PATCH 15/30] add appropriate test --- examples/unstructured_2d_dgsem/elixir_euler_time_series.jl | 1 - test/test_unstructured_2d.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl index c8b930fc64a..670b17041d8 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl @@ -103,7 +103,6 @@ time_series = TimeSeriesCallback(semi, callbacks = CallbackSet(summary_callback, analysis_callback, - save_solution, time_series, alive_callback) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 0f72ce03818..a2f5c80883f 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -221,7 +221,7 @@ end du_ode = similar(u_ode) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end - # Extra tests to make sure the `TimeSeriesCallback` made correct data + # Extra test to make sure the `TimeSeriesCallback` made correct data point_data_1 = time_series.affect!.point_data[1] @test all(isapprox.(point_data_1[1:4], [1.9548629504179071, 1.9548895017660788, From ce44fe6935aac3bbb4176f61efda640495b7cbd8 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 9 Mar 2024 00:33:06 +0100 Subject: [PATCH 16/30] update docstring --- src/callbacks_step/time_series.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl index 7baa6b9c5a1..e8a3cecc211 100644 --- a/src/callbacks_step/time_series.jl +++ b/src/callbacks_step/time_series.jl @@ -24,7 +24,8 @@ The real data type `RealT` and data type for solution variables `uEltype` defaul types used in the solver and the cache. !!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. + This is an experimental feature and may change in future releases. Currently this callback + is only implemented for [`TreeMesh`](@ref) in 2D and [`UnstructuredMesh2D`](@ref). """ mutable struct TimeSeriesCallback{RealT <: Real, uEltype <: Real, SolutionVariables, VariableNames, Cache} From 04ffc1ac42960c807f8df98d63c825367a320cad Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 9 Mar 2024 00:42:21 +0100 Subject: [PATCH 17/30] add comment about the barycenter computation --- src/callbacks_step/time_series_dg2d.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 0e554daf836..d4c69ee1a90 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -156,6 +156,9 @@ function get_elements_by_coordinates!(element_ids, coordinates, end # Use the available `node_coordinates` on each element to compute and save the barycenter. +# In essence, the barycenter is like an average where all the x and y node coordinates are +# summed and then we divide by the total number of degrees of freedom on the element, i.e., +# the value of `n^2` in two spatial dimensions. @inline function calc_bary_centers!(bary_centers, dg, cache) n = nnodes(dg) @views for element in eachelement(dg, cache) From 0682471bbda28dbfce9137ce1b0be19a8d3af7b2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 08:06:33 +0100 Subject: [PATCH 18/30] add simplifications and comments from code review --- .../elixir_euler_time_series.jl | 2 +- src/callbacks_step/time_series_dg2d.jl | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl index 670b17041d8..e5af6ee1dad 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl @@ -2,7 +2,7 @@ # the `TimeSeriesCallback` on several gauge points. Many of the # gauge points are selected as "stress tests" for the element # identification, e.g., a gauge point that lies on an -# element corner of a curvlinear mesh +# element corner of a curvilinear mesh using OrdinaryDiffEq using Trixi diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index d4c69ee1a90..2b1fe9b302f 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -131,7 +131,9 @@ function get_elements_by_coordinates!(element_ids, coordinates, P = surface_point .- point # If the vector `P` is the zero vector then this `point` is at an element corner or - # on a surface. So any candidate element works, just use the first one + # on a surface. In this case the choice of a candidate element is ambiguous and + # we just use the first candidate. However, solutions might differ at discontinuous + # interfaces such that this choice may influence the result. if sum(P .* P) < 500 * eps(eltype(point)) element_ids[index] = candidates[1] continue @@ -174,12 +176,12 @@ end # using the available `node_coordinates`. Also return the index pair of this # minimum surface point location. We compute the squared norm to avoid # computing computationally more expensive square roots. -# OBS! Could be made more accuracte if the `node_coordinates` were super-sampled +# Note! Could be made more accurate if the `node_coordinates` were super-sampled # and reinterpolated onto a higher polynomial degree before this computation. function calc_minimum_surface_distance(point, node_coordinates, dg, mesh::UnstructuredMesh2D) n = nnodes(dg) - hmin2 = Inf * ones(eltype(mesh.corners), length(mesh)) + min_distance_squared = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) for k in 1:length(mesh) # used to ensure that only boundary points are used @@ -191,10 +193,8 @@ function calc_minimum_surface_distance(point, node_coordinates, if !any(on_surface) continue end - if sum((node_coordinates[:, i, j, k] - point) .* - (node_coordinates[:, i, j, k] - point)) < hmin2[k] - hmin2[k] = sum((node_coordinates[:, i, j, k] - point) .* - (node_coordinates[:, i, j, k] - point)) + if sum((node_coordinates[:, i, j, k] - point).^2) < min_distance_squared[k] + min_distance_squared[k] = sum((node_coordinates[:, i, j, k] - point).^2) indices[k, 1] = i indices[k, 2] = j end @@ -202,7 +202,7 @@ function calc_minimum_surface_distance(point, node_coordinates, end end - return hmin2, indices + return min_distance_squared, indices end function get_elements_by_coordinates(coordinates, mesh, dg, cache) From 5435e1b7bf576383c3bd5f1e0d6a9f8c4a826ffa Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 13:19:19 +0100 Subject: [PATCH 19/30] adjust variable name to avoid ugly formatting --- src/callbacks_step/time_series_dg2d.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 2b1fe9b302f..a74e4ef9bc7 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -174,14 +174,14 @@ end # Compute the shortest distance from a `point` to the surface of each element # using the available `node_coordinates`. Also return the index pair of this -# minimum surface point location. We compute the squared norm to avoid -# computing computationally more expensive square roots. +# minimum surface point location. We compute and store in `min_distance` +# the squared norm to avoid computing computationally more expensive square roots. # Note! Could be made more accurate if the `node_coordinates` were super-sampled # and reinterpolated onto a higher polynomial degree before this computation. function calc_minimum_surface_distance(point, node_coordinates, dg, mesh::UnstructuredMesh2D) n = nnodes(dg) - min_distance_squared = Inf * ones(eltype(mesh.corners), length(mesh)) + min_distance = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) for k in 1:length(mesh) # used to ensure that only boundary points are used @@ -193,8 +193,8 @@ function calc_minimum_surface_distance(point, node_coordinates, if !any(on_surface) continue end - if sum((node_coordinates[:, i, j, k] - point).^2) < min_distance_squared[k] - min_distance_squared[k] = sum((node_coordinates[:, i, j, k] - point).^2) + if sum((node_coordinates[:, i, j, k] - point) .^ 2) < min_distance[k] + min_distance[k] = sum((node_coordinates[:, i, j, k] - point) .^ 2) indices[k, 1] = i indices[k, 2] = j end From b002920b02f885ec7c5bd8eb192b564a3ac7c069 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 15:06:10 +0100 Subject: [PATCH 20/30] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/callbacks_step/time_series_dg2d.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index a74e4ef9bc7..26efcb286de 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -128,7 +128,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, candidates[1]]) # Compute the vector pointing from the current `point` toward the surface - P = surface_point .- point + P = surface_point - point # If the vector `P` is the zero vector then this `point` is at an element corner or # on a surface. In this case the choice of a candidate element is ambiguous and @@ -146,7 +146,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, bary_center = SVector(bary_centers[1, candidates[element]], bary_centers[2, candidates[element]]) # Vector pointing from the barycenter toward the minimal `surface_point` - B = surface_point .- bary_center + B = surface_point - bary_center if sum(P .* B) > zero(eltype(bary_center)) element_ids[index] = candidates[element] break @@ -185,7 +185,7 @@ function calc_minimum_surface_distance(point, node_coordinates, indices = zeros(Int, length(mesh), 2) for k in 1:length(mesh) # used to ensure that only boundary points are used - on_surface = [false, false] + on_surface = MVector(false, false) for j in 1:n on_surface[2] = (j == 1) || (j == n) for i in 1:n @@ -193,8 +193,11 @@ function calc_minimum_surface_distance(point, node_coordinates, if !any(on_surface) continue end - if sum((node_coordinates[:, i, j, k] - point) .^ 2) < min_distance[k] - min_distance[k] = sum((node_coordinates[:, i, j, k] - point) .^ 2) + node = SVector(node_coordinates[1, i, j, k], + node_coordinates[2, i, j, k]) + distance2 = sum(abs2, node - point) + if distance2 < min_distance[k] + min_distance[k] = distance2 indices[k, 1] = i indices[k, 2] = j end From 7147a106cf151b74b836969c853dfeeae1c16b1c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 15:17:25 +0100 Subject: [PATCH 21/30] fix variable name --- src/callbacks_step/time_series_dg2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl index 26efcb286de..ad7c6851c80 100644 --- a/src/callbacks_step/time_series_dg2d.jl +++ b/src/callbacks_step/time_series_dg2d.jl @@ -181,7 +181,7 @@ end function calc_minimum_surface_distance(point, node_coordinates, dg, mesh::UnstructuredMesh2D) n = nnodes(dg) - min_distance = Inf * ones(eltype(mesh.corners), length(mesh)) + min_distance2 = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) for k in 1:length(mesh) # used to ensure that only boundary points are used @@ -196,8 +196,8 @@ function calc_minimum_surface_distance(point, node_coordinates, node = SVector(node_coordinates[1, i, j, k], node_coordinates[2, i, j, k]) distance2 = sum(abs2, node - point) - if distance2 < min_distance[k] - min_distance[k] = distance2 + if distance2 < min_distance2[k] + min_distance2[k] = distance2 indices[k, 1] = i indices[k, 2] = j end @@ -205,7 +205,7 @@ function calc_minimum_surface_distance(point, node_coordinates, end end - return min_distance_squared, indices + return min_distance2, indices end function get_elements_by_coordinates(coordinates, mesh, dg, cache) From fe4f4313888e6b1b7f74fa0811fdd467c21acddf Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 19:39:57 +0100 Subject: [PATCH 22/30] remove Experimental status from the TimeSeriesCallback --- src/callbacks_step/time_series.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl index e8a3cecc211..8ca8fe4db8a 100644 --- a/src/callbacks_step/time_series.jl +++ b/src/callbacks_step/time_series.jl @@ -23,9 +23,8 @@ After the last time step, the results are stored in an HDF5 file `filename` in d The real data type `RealT` and data type for solution variables `uEltype` default to the respective types used in the solver and the cache. -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. Currently this callback - is only implemented for [`TreeMesh`](@ref) in 2D and [`UnstructuredMesh2D`](@ref). +Currently this callback is only implemented for [`TreeMesh`](@ref) in 2D +and [`UnstructuredMesh2D`](@ref). """ mutable struct TimeSeriesCallback{RealT <: Real, uEltype <: Real, SolutionVariables, VariableNames, Cache} From fa5daa800d30f4a1c8dd4b62479d065982f429f2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 20:00:51 +0100 Subject: [PATCH 23/30] move new TimeSeries test into the unit testing --- test/test_unit.jl | 11 +++++++++++ test/test_unstructured_2d.jl | 5 ----- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 1907a281718..68b99c488c2 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -600,6 +600,7 @@ end end @timed_testset "TimeSeriesCallback" begin + # Test the 2D TreeMesh version of the callback and some warnings @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_acoustics_gaussian_source.jl"), @@ -615,6 +616,16 @@ end @test_nowarn show(stdout, time_series) @test_throws ArgumentError TimeSeriesCallback(semi, [(1.0, 1.0)]; interval = -1) @test_throws ArgumentError TimeSeriesCallback(semi, [1.0 1.0 1.0; 2.0 2.0 2.0]) + # Test the 2D unstructured mesh version of the callback + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "unstructured_2d_dgsem", + "elixir_euler_time_series.jl"), + tspan = (0.0, 0.2)) + + point_data_1 = time_series.affect!.point_data[1] + @test all(isapprox.(point_data_1[1:4], + [ 1.9548629504390056, 1.9548895017871935, + 1.9548892928720158, 3.821760762385351])) end @timed_testset "Consistency check for single point flux: CEMCE" begin diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index a2f5c80883f..deb3377f6bd 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -221,11 +221,6 @@ end du_ode = similar(u_ode) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end - # Extra test to make sure the `TimeSeriesCallback` made correct data - point_data_1 = time_series.affect!.point_data[1] - @test all(isapprox.(point_data_1[1:4], - [1.9548629504179071, 1.9548895017660788, - 1.954889292850916, 3.8217607623030623])) end @trixi_testset "elixir_acoustics_gauss_wall.jl" begin From 771a3c5670f6fd7e918bd230e12205289b669909 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 20:08:57 +0100 Subject: [PATCH 24/30] add output_directory creation if not already done. Necessary if this callback is used without the SaveSolution callback --- src/callbacks_step/time_series.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl index 8ca8fe4db8a..f6d76f0fb15 100644 --- a/src/callbacks_step/time_series.jl +++ b/src/callbacks_step/time_series.jl @@ -96,6 +96,11 @@ function TimeSeriesCallback(mesh, equations, solver, cache, point_coordinates; throw(ArgumentError("`point_coordinates` must be a matrix of size n_points × ndims")) end + # create the output folder if it does not exist already + if mpi_isroot() && !isdir(output_directory) + mkpath(output_directory) + end + # Transpose point_coordinates to our usual format [ndims, n_points] # Note: They are accepted in a different format to allow direct input from `readdlm` point_coordinates_ = permutedims(point_coordinates) From e65207938838fae43b5484614265ac3be227a4dd Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 12 Mar 2024 20:14:01 +0100 Subject: [PATCH 25/30] formatting --- test/test_unit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 68b99c488c2..54de750adc5 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -624,7 +624,7 @@ end point_data_1 = time_series.affect!.point_data[1] @test all(isapprox.(point_data_1[1:4], - [ 1.9548629504390056, 1.9548895017871935, + [1.9548629504390056, 1.9548895017871935, 1.9548892928720158, 3.821760762385351])) end From 823891bc2918ab0bee62cc26349a3a1d76ea52ef Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Mar 2024 09:45:29 +0100 Subject: [PATCH 26/30] update test mesh to have one straight-sided element to trigger inverse bi-linear interpolation --- examples/unstructured_2d_dgsem/elixir_euler_time_series.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl index e5af6ee1dad..13233cdadbc 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl @@ -69,7 +69,7 @@ solver = DGSEM(polydeg = 6, surface_flux = flux_lax_friedrichs) ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) -mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh", +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/55c916cd8c0294a2d4a836e960dac7247b7c8ccf/mesh_multiple_flips.mesh", joinpath(@__DIR__, "mesh_multiple_flips.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) From afc6859d6189fed02a1fb9667c8d66558146e03f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Mar 2024 10:11:38 +0100 Subject: [PATCH 27/30] update test values --- test/test_unit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 54de750adc5..1d36c29011c 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -624,8 +624,8 @@ end point_data_1 = time_series.affect!.point_data[1] @test all(isapprox.(point_data_1[1:4], - [1.9548629504390056, 1.9548895017871935, - 1.9548892928720158, 3.821760762385351])) + [1.9546882708551676, 1.9547149531788077, + 1.9547142161310154, 3.821066781119142])) end @timed_testset "Consistency check for single point flux: CEMCE" begin From e73e3e7a30756d7178a6b034a5abfc168abef9fe Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Mar 2024 10:15:13 +0100 Subject: [PATCH 28/30] add news item --- NEWS.md | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index d70504d8c85..5b08d51ab89 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,13 +4,19 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes in the v0.7 lifecycle + +#### Added +- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. + + ## Changes when updating to v0.7 from v0.6.x #### Added #### Changed -- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` +- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` instead of `min_max_speed_naive`. #### Deprecated @@ -26,7 +32,7 @@ for human readability. #### Added - AMR for hyperbolic-parabolic equations on 3D `P4estMesh` - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` -- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, +- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` - Added Lighthill-Whitham-Richards (LWR) traffic model @@ -40,7 +46,7 @@ for human readability. #### Changed - The wave speed estimates for `flux_hll`, `FluxHLL()` are now consistent across equations. - In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now + In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now conceptually identical across equations. Users, who have been using `flux_hll` for MHD have now to use `flux_hlle` in order to use the Einfeldt wave speed estimate. From 099050dcf9e6c0190341b717b5694f179fde3452 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Mar 2024 10:51:43 +0100 Subject: [PATCH 29/30] forgot to update all new test values on the new mesh --- test/test_unstructured_2d.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index deb3377f6bd..66d11afd014 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -201,16 +201,16 @@ end @trixi_testset "elixir_euler_time_series.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_time_series.jl"), l2=[ - 9.44589994864801e-5, - 7.375567345161402e-5, - 8.469037537040469e-5, - 0.0001638067441092813, + 9.860036213913957e-5, + 7.767842271271052e-5, + 8.956109781622046e-5, + 0.0001771544151538529, ], linf=[ - 0.0006797328621439558, - 0.00048723709312858965, - 0.0006701229224337357, - 0.0011917195110209278, + 0.000656954971669732, + 0.000494381111073805, + 0.0006563610876975101, + 0.0011651838899400246, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations From d81b79ee7a8ffacc24de217033f022437a66a63f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Mar 2024 20:34:59 +0100 Subject: [PATCH 30/30] update tests and use coverage override to avoid redundancy --- test/test_unit.jl | 10 ---------- test/test_unstructured_2d.jl | 26 +++++++++++++++++--------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 1d36c29011c..03a78f6918a 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -616,16 +616,6 @@ end @test_nowarn show(stdout, time_series) @test_throws ArgumentError TimeSeriesCallback(semi, [(1.0, 1.0)]; interval = -1) @test_throws ArgumentError TimeSeriesCallback(semi, [1.0 1.0 1.0; 2.0 2.0 2.0]) - # Test the 2D unstructured mesh version of the callback - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "unstructured_2d_dgsem", - "elixir_euler_time_series.jl"), - tspan = (0.0, 0.2)) - - point_data_1 = time_series.affect!.point_data[1] - @test all(isapprox.(point_data_1[1:4], - [1.9546882708551676, 1.9547149531788077, - 1.9547142161310154, 3.821066781119142])) end @timed_testset "Consistency check for single point flux: CEMCE" begin diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 66d11afd014..6814250dd47 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -201,18 +201,26 @@ end @trixi_testset "elixir_euler_time_series.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_time_series.jl"), l2=[ - 9.860036213913957e-5, - 7.767842271271052e-5, - 8.956109781622046e-5, - 0.0001771544151538529, + 6.984024099236519e-5, + 6.289022520363763e-5, + 6.550951878107466e-5, + 0.00016222767700879948, ], linf=[ - 0.000656954971669732, - 0.000494381111073805, - 0.0006563610876975101, - 0.0011651838899400246, + 0.0005367823248620951, + 0.000671293180158461, + 0.0005656680962440319, + 0.0013910024779804075, ], - tspan=(0.0, 0.5)) + tspan=(0.0, 0.2), + # With the default `maxiters = 1` in coverage tests, + # there would be no time series to check against. + coverage_override=(maxiters = 20,)) + # Extra test that the `TimeSeries` callback creates reasonable data + point_data_1 = time_series.affect!.point_data[1] + @test all(isapprox.(point_data_1[1:4], + [1.9546882708551676, 1.9547149531788077, + 1.9547142161310154, 3.821066781119142])) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let