From c9756dae91f28a6bc68e753dc6f863c07b7bfd18 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Thu, 30 Jun 2022 16:30:50 +0200 Subject: [PATCH 01/21] Add node variables --- src/convert.jl | 12 +++++++++++- src/interpolate.jl | 43 +++++++++++++++++++++++++++++++++++++++++++ src/io.jl | 17 ++++++++++++++++- 3 files changed, 70 insertions(+), 2 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 0a6b66a..5535eb1 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -107,7 +107,7 @@ function trixi2vtk(filename::AbstractString...; if is_datafile verbose && println("| Reading data file...") @timeit "read data" (labels, data, n_elements, n_nodes, - element_variables, time) = read_datafile(filename) + element_variables, node_variables, time) = read_datafile(filename) assert_cells_elements(n_elements, mesh, filename, meshfile) @@ -154,6 +154,16 @@ function trixi2vtk(filename::AbstractString...; verbose && println("| | Element variable: $label...") @timeit label vtk_celldata[label] = variable end + + # Add node variables + for (label, variable) in node_variables + verbose && println("| | Node variable: $label...") + @timeit "interpolate cell data" interpolated_cell_data = interpolate_cell_data(Val(format), + variable, mesh, + n_visnodes, verbose) + # Add the "interpolated" cell_data to celldata, not node_data + @timeit label vtk_nodedata[label] = interpolated_cell_data + end end end end diff --git a/src/interpolate.jl b/src/interpolate.jl index 50dd902..3c56f05 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -8,6 +8,49 @@ function interpolate_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, v return raw2interpolated(input_data, nodes_out) end +function interpolate_cell_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, verbose) + # Calculate equidistant output nodes (visualization nodes) + nodes_out = collect(range(-1, 1, length=n_visnodes)) + + # Extract number of spatial dimensions + ndims_ = ndims(mesh) + + # Extract data shape information + n_nodes_in = size(input_data, 1) + n_nodes_out = length(nodes_out) + n_elements = size(input_data, ndims_ + 1) + + # Get node coordinates for DG locations on reference element + nodes_in, weights_in = gauss_lobatto_nodes_weights(n_nodes_in) + + if ndims_ == 2 + # Create output data structure + data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_elements) + + # Interpolate node data for visualization nodes to piecewise constant values and store to global data structure + for element_id in 1:n_elements + index_j = 1 + for j in 1:n_nodes_out + index_i = 1 + for i in 1:n_nodes_out + while -1.0 + sum(weights_in[1:index_i]) < nodes_out[i] + index_i += 1 + end + while -1.0 + sum(weights_in[1:index_j]) < nodes_out[j] + index_j += 1 + end + data_vis[i, j, element_id] = input_data[index_i, index_j, element_id] + end + end + end + else + error("Unsupported number of spatial dimensions: ", ndims_) + end + + # Return as one 1D array + return reshape(data_vis, n_nodes_out^ndims_ * n_elements) +end + # Interpolate data from input format to desired output format (StructuredMesh or UnstructuredMesh2D version) function interpolate_data(::Val{:vtu}, input_data, diff --git a/src/io.jl b/src/io.jl index 54b79fd..997cdf5 100644 --- a/src/io.jl +++ b/src/io.jl @@ -76,6 +76,21 @@ function read_datafile(filename::String) index +=1 end - return labels, data, n_elements, n_nodes, element_variables, time + # Extract node variable arrays + node_variables = Dict{String, Union{Array{Float64}, Array{Int}}}() + index = 1 + while haskey(file, "node_variables_$index") + varname = read(attributes(file["node_variables_$index"])["name"]) + nodedata = read(file["node_variables_$index"]) + if ndims_ == 2 + node_variables[varname] = Array{Float64}(undef, n_nodes, n_nodes, n_elements) + @views node_variables[varname][:, :, :] .= nodedata + else + error("Unsupported number of spatial dimensions: ", ndims_) + end + index +=1 + end + + return labels, data, n_elements, n_nodes, element_variables, node_variables, time end end From 2cef2a3bdb373050c9f9799f91c0f27c1f5dd2a2 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Tue, 22 Nov 2022 14:20:32 +0100 Subject: [PATCH 02/21] Add StructuredMesh --- src/interpolate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolate.jl b/src/interpolate.jl index 3c56f05..ec97697 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -8,7 +8,7 @@ function interpolate_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, v return raw2interpolated(input_data, nodes_out) end -function interpolate_cell_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, verbose) +function interpolate_cell_data(::Val{:vtu}, input_data, mesh::Union{TreeMesh, StructuredMesh}, n_visnodes, verbose) # Calculate equidistant output nodes (visualization nodes) nodes_out = collect(range(-1, 1, length=n_visnodes)) From 0f9202cd491f8c8b4cf8c7ce11d2df24df1b3e26 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 18 Jan 2023 09:29:56 +0100 Subject: [PATCH 03/21] Add branch to allow reinterpolate=false --- .gitignore | 1 + src/convert.jl | 75 +++++-- src/interpolate.jl | 251 ++++++++++----------- src/vtktools.jl | 218 +++++++++++++++--- test/Project.toml | 3 +- test/test_2d.jl | 490 +++++++++++++++++++++++++++++------------ test/test_3d.jl | 319 +++++++++++++++++++++------ test/test_manual.jl | 6 + test/test_trixi2vtk.jl | 151 ++++++++++++- 9 files changed, 1116 insertions(+), 398 deletions(-) diff --git a/.gitignore b/.gitignore index f1677ce..821d95b 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ **/Manifest.toml out/ artifacts/ +reference_files/ docs/build public/ coverage/ diff --git a/src/convert.jl b/src/convert.jl index 3e228de..a21527f 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -1,7 +1,8 @@ """ trixi2vtk(filename::AbstractString...; format=:vtu, verbose=false, hide_progress=false, pvd=nothing, - output_directory=".", nvisnodes=nothing) + output_directory=".", nvisnodes=nothing, save_celldata=true, + reinterpolate=true, data_is_uniform=false) Convert Trixi-generated output files to VTK files (VTU or VTI). @@ -21,6 +22,13 @@ Convert Trixi-generated output files to VTK files (VTU or VTI). A value of `0` (zero) uses the number of nodes in the DG elements. - `save_celldata`: Boolean value to determine if cell-based data should be saved. (default: `true`) +- `reinterpolate`: Boolean value to determine if data should be reinterpolated + onto uniform points. When `false` the raw data at the compute nodes + is copied into the appropriate format. + (default: `true`) +- `data_is_uniform`: Boolean to indicate if the data to be converted is from a finite difference + method on a uniform grid of points. + (default: `false`) # Examples ```julia @@ -30,7 +38,8 @@ julia> trixi2vtk("out/solution_000*.h5") """ function trixi2vtk(filename::AbstractString...; format=:vtu, verbose=false, hide_progress=false, pvd=nothing, - output_directory=".", nvisnodes=nothing, save_celldata=true) + output_directory=".", nvisnodes=nothing, save_celldata=true, + reinterpolate=true, data_is_uniform=false) # Reset timer reset_timer!() @@ -103,6 +112,9 @@ function trixi2vtk(filename::AbstractString...; verbose && println("| Reading mesh file...") @timeit "read mesh" mesh = Trixi.load_mesh_serial(meshfile; n_cells_max=0, RealT=Float64) + # Create an empty `node_set` such that converting a mesh.h5 file also works + node_set = [] + # Read data only if it is a data file if is_datafile verbose && println("| Reading data file...") @@ -113,24 +125,53 @@ function trixi2vtk(filename::AbstractString...; # Determine resolution for data interpolation n_visnodes = get_default_nvisnodes_solution(nvisnodes, n_nodes, mesh) + + # If a user requests that no reinterpolation is done automatically set + # `n_visnodes` to be the same as the number of nodes in the raw data. + if !reinterpolate + n_visnodes = n_nodes + end else # If file is a mesh file, do not interpolate data as detailed n_visnodes = get_default_nvisnodes_mesh(nvisnodes, mesh) end + # Check if the raw data is uniform (finite difference) or not (dg) + # and create the corresponding node set for reinterpolation / copying. + if (reinterpolate & !data_is_uniform) | (!reinterpolate & data_is_uniform) + # (1) Default settings; presumably the most common + # (2) Finite difference data + node_set = range(-1, 1, length=n_visnodes) + elseif !reinterpolate & !data_is_uniform + # raw data is on a set of LGL nodes + node_set, _ = gauss_lobatto_nodes_weights(n_visnodes) + else # reinterpolate & data_is_uniform + error("uniform data should not be reinterpolated! Set reinterpolate=false and try again.") + end + # Create output directory if it does not exist mkpath(output_directory) # Build VTK grids - vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, n_visnodes, verbose, - output_directory, is_datafile, filename) + vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, node_set, n_visnodes, verbose, + output_directory, is_datafile, filename, Val(reinterpolate)) # Interpolate data if is_datafile verbose && println("| Interpolating data...") - @timeit "interpolate data" interpolated_data = interpolate_data(Val(format), - data, mesh, - n_visnodes, verbose) + if reinterpolate + @timeit "interpolate data" interpolated_data = interpolate_data(Val(format), + data, mesh, + n_visnodes, verbose) + else # Copy the raw solution data; only works for `vtu` format + # Extract data shape information + ndims_ = ndims(data) - 2 + n_variables = length(labels) + # Save raw data as one 1D array for each variable + @timeit "interpolate data" interpolated_data = reshape(data, + n_visnodes^ndims_ * n_elements, + n_variables) + end end # Add data to file @@ -158,9 +199,14 @@ function trixi2vtk(filename::AbstractString...; # Add node variables for (label, variable) in node_variables verbose && println("| | Node variable: $label...") - @timeit "interpolate cell data" interpolated_cell_data = interpolate_cell_data(Val(format), + if reinterpolate + @timeit "interpolate cell data" interpolated_cell_data = interpolate_cell_data(Val(format), variable, mesh, n_visnodes, verbose) + else + @timeit "interpolate cell data" interpolated_cell_data = reshape(variable, + n_visnodes^ndims_ * n_elements) + end # Add the "interpolated" cell_data to celldata, not node_data @timeit label vtk_nodedata[label] = interpolated_cell_data end @@ -282,19 +328,10 @@ end # default number of visualization nodes if only the mesh should be visualized -function get_default_nvisnodes_mesh(nvisnodes, mesh::TreeMesh) - if nvisnodes === nothing - # for a Cartesian mesh, we do not need to interpolate - return 1 - else - return nvisnodes - end -end - function get_default_nvisnodes_mesh(nvisnodes, - mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh}) + mesh::Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh}) if nvisnodes === nothing - # for curved meshes, we need to get at least the vertices + # we need to get at least the vertices return 2 else return nvisnodes diff --git a/src/interpolate.jl b/src/interpolate.jl index ec97697..fd9d83b 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -62,130 +62,132 @@ function interpolate_data(::Val{:vtu}, input_data, return raw2interpolated(input_data, nodes_out) end - -# Interpolate data from input format to desired output format (vti version) -function interpolate_data(::Val{:vti}, input_data, mesh, n_visnodes, verbose) - coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) - - # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² - normalized_coordinates = similar(coordinates) - for element_id in axes(coordinates, 2) - @views normalized_coordinates[:, element_id] .= ( - (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 )) - end - - # Determine level-wise resolution - max_level = maximum(levels) - resolution = n_visnodes * 2^max_level - - # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate - # level-0-cell) that contains the number of visualization nodes for any - # refinement level to visualize on an equidistant grid - nvisnodes_per_level = [2^(max_level - level)*n_visnodes for level in 0:max_level] - - # Interpolate unstructured DG data to structured data - structured_data = unstructured2structured(input_data, normalized_coordinates, levels, - resolution, nvisnodes_per_level) - - return structured_data -end - - -# Interpolate unstructured DG data to structured data (cell-centered) -function unstructured2structured(unstructured_data::AbstractArray{Float64}, - normalized_coordinates::AbstractArray{Float64}, - levels::AbstractArray{Int}, resolution::Int, - nvisnodes_per_level::AbstractArray{Int}) - # Extract number of spatial dimensions - ndims_ = size(normalized_coordinates, 1) - - # Extract data shape information - n_nodes_in, _, n_elements, n_variables = size(unstructured_data) - - # Get node coordinates for DG locations on reference element - nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) - - # Calculate interpolation vandermonde matrices for each level - max_level = length(nvisnodes_per_level) - 1 - vandermonde_per_level = [] - for l in 0:max_level - n_nodes_out = nvisnodes_per_level[l + 1] - dx = 2 / n_nodes_out - nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out)) - push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out)) - end - - # For each element, calculate index position at which to insert data in global data structure - lower_left_index = element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level) - - # Create output data structure - structured = Array{Float64}(undef, resolution, resolution, n_variables) - - # For each variable, interpolate element data and store to global data structure - for v in 1:n_variables - # Reshape data array for use in interpolate_nodes function - reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements) - - for element_id in 1:n_elements - # Extract level for convenience - level = levels[element_id] - - # Determine target indices - n_nodes_out = nvisnodes_per_level[level + 1] - first = lower_left_index[:, element_id] - last = first .+ (n_nodes_out - 1) - - # Interpolate data - vandermonde = vandermonde_per_level[level + 1] - structured[first[1]:last[1], first[2]:last[2], v] .= ( - reshape(interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1), - n_nodes_out, n_nodes_out)) - end - end - - # Return as one 1D array for each variable - return reshape(structured, resolution^ndims_, n_variables) -end - - -# For a given normalized element coordinate, return the index of its lower left -# contribution to the global data structure -function element2index(normalized_coordinates::AbstractArray{Float64}, levels::AbstractArray{Int}, - resolution::Int, nvisnodes_per_level::AbstractArray{Int}) - # Extract number of spatial dimensions - ndims_ = size(normalized_coordinates, 1) - - n_elements = length(levels) - - # First, determine lower left coordinate for all cells - dx = 2 / resolution - lower_left_coordinate = Array{Float64}(undef, ndims_, n_elements) - for element_id in 1:n_elements - nvisnodes = nvisnodes_per_level[levels[element_id] + 1] - lower_left_coordinate[1, element_id] = ( - normalized_coordinates[1, element_id] - (nvisnodes - 1)/2 * dx) - lower_left_coordinate[2, element_id] = ( - normalized_coordinates[2, element_id] - (nvisnodes - 1)/2 * dx) - end - - # Then, convert coordinate to global index - indices = coordinate2index(lower_left_coordinate, resolution) - - return indices -end - - -# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1]) -function coordinate2index(coordinate, resolution::Integer) - # Calculate 1D normalized coordinates - dx = 2/resolution - mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution)) - - # Find index - id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], lt=(x,y)->x .< y .- dx/2) - id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], lt=(x,y)->x .< y .- dx/2) - return transpose(hcat(id_x, id_y)) -end +# +# # TODO: remove the `vti` option and all its helper routines +# +# # Interpolate data from input format to desired output format (vti version) +# function interpolate_data(::Val{:vti}, input_data, mesh, n_visnodes, verbose) +# coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + +# # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² +# normalized_coordinates = similar(coordinates) +# for element_id in axes(coordinates, 2) +# @views normalized_coordinates[:, element_id] .= ( +# (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 )) +# end + +# # Determine level-wise resolution +# max_level = maximum(levels) +# resolution = n_visnodes * 2^max_level + +# # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate +# # level-0-cell) that contains the number of visualization nodes for any +# # refinement level to visualize on an equidistant grid +# nvisnodes_per_level = [2^(max_level - level)*n_visnodes for level in 0:max_level] + +# # Interpolate unstructured DG data to structured data +# structured_data = unstructured2structured(input_data, normalized_coordinates, levels, +# resolution, nvisnodes_per_level) + +# return structured_data +# end + + +# # Interpolate unstructured DG data to structured data (cell-centered) +# function unstructured2structured(unstructured_data::AbstractArray{Float64}, +# normalized_coordinates::AbstractArray{Float64}, +# levels::AbstractArray{Int}, resolution::Int, +# nvisnodes_per_level::AbstractArray{Int}) +# # Extract number of spatial dimensions +# ndims_ = size(normalized_coordinates, 1) + +# # Extract data shape information +# n_nodes_in, _, n_elements, n_variables = size(unstructured_data) + +# # Get node coordinates for DG locations on reference element +# nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + +# # Calculate interpolation vandermonde matrices for each level +# max_level = length(nvisnodes_per_level) - 1 +# vandermonde_per_level = [] +# for l in 0:max_level +# n_nodes_out = nvisnodes_per_level[l + 1] +# dx = 2 / n_nodes_out +# nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out)) +# push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out)) +# end + +# # For each element, calculate index position at which to insert data in global data structure +# lower_left_index = element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level) + +# # Create output data structure +# structured = Array{Float64}(undef, resolution, resolution, n_variables) + +# # For each variable, interpolate element data and store to global data structure +# for v in 1:n_variables +# # Reshape data array for use in interpolate_nodes function +# reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements) + +# for element_id in 1:n_elements +# # Extract level for convenience +# level = levels[element_id] + +# # Determine target indices +# n_nodes_out = nvisnodes_per_level[level + 1] +# first = lower_left_index[:, element_id] +# last = first .+ (n_nodes_out - 1) + +# # Interpolate data +# vandermonde = vandermonde_per_level[level + 1] +# structured[first[1]:last[1], first[2]:last[2], v] .= ( +# reshape(interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1), +# n_nodes_out, n_nodes_out)) +# end +# end + +# # Return as one 1D array for each variable +# return reshape(structured, resolution^ndims_, n_variables) +# end + + +# # For a given normalized element coordinate, return the index of its lower left +# # contribution to the global data structure +# function element2index(normalized_coordinates::AbstractArray{Float64}, levels::AbstractArray{Int}, +# resolution::Int, nvisnodes_per_level::AbstractArray{Int}) +# # Extract number of spatial dimensions +# ndims_ = size(normalized_coordinates, 1) + +# n_elements = length(levels) + +# # First, determine lower left coordinate for all cells +# dx = 2 / resolution +# lower_left_coordinate = Array{Float64}(undef, ndims_, n_elements) +# for element_id in 1:n_elements +# nvisnodes = nvisnodes_per_level[levels[element_id] + 1] +# lower_left_coordinate[1, element_id] = ( +# normalized_coordinates[1, element_id] - (nvisnodes - 1)/2 * dx) +# lower_left_coordinate[2, element_id] = ( +# normalized_coordinates[2, element_id] - (nvisnodes - 1)/2 * dx) +# end + +# # Then, convert coordinate to global index +# indices = coordinate2index(lower_left_coordinate, resolution) + +# return indices +# end + + +# # Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1]) +# function coordinate2index(coordinate, resolution::Integer) +# # Calculate 1D normalized coordinates +# dx = 2/resolution +# mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution)) + +# # Find index +# id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], lt=(x,y)->x .< y .- dx/2) +# id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], lt=(x,y)->x .< y .- dx/2) +# return transpose(hcat(id_x, id_y)) +# end # Interpolate to specified output nodes @@ -254,6 +256,7 @@ function interpolate_nodes(data_in::AbstractArray{T, 3}, interpolate_nodes!(data_out, data_in, vandermonde, n_vars) end + function interpolate_nodes!(data_out::AbstractArray{T, 3}, data_in::AbstractArray{T, 3}, vandermonde, n_vars) where T n_nodes_out = size(vandermonde, 1) diff --git a/src/vtktools.jl b/src/vtktools.jl index 5ee500b..7b863e8 100644 --- a/src/vtktools.jl +++ b/src/vtktools.jl @@ -1,6 +1,7 @@ -# Create and return VTK grids that are ready to be filled with data (vtu version) -function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, n_visnodes, verbose, - output_directory, is_datafile, filename) +# Create and return VTK grids on enriched uniform nodes +# that are ready to be filled with reinterpolated data (vtu version) +function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, nodes, n_visnodes, verbose, + output_directory, is_datafile, filename, reinterpolate::Val{true}) coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) # Extract number of spatial dimensions @@ -43,18 +44,27 @@ function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, n_visnodes, verbose, end -# Create and return VTK grids that are ready to be filled with data (vti version) -function build_vtk_grids(::Val{:vti}, mesh, n_visnodes, verbose, - output_directory, is_datafile, filename) +# Create and return VTK grids with nodes on which the raw was created +# ready to be filled with raw data (vtu version) +function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, + nodes, n_visnodes, verbose, + output_directory, is_datafile, filename, reinterpolate::Val{false}) - coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, nodes, n_visnodes) - # Extract number of spatial dimensions - ndims_ = size(coordinates, 1) + # Calculate VTK points and cells + verbose && println("| Preparing VTK cells...") + if is_datafile + @timeit "prepare VTK cells (node data)" begin + vtk_points, vtk_cells = calc_vtk_points_cells(node_coordinates) + end + end # Prepare VTK points and cells for celldata file - @timeit "prepare VTK cells" vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells( - Val(ndims_), coordinates, levels, center_level_0, length_level_0, 1) + @timeit "prepare VTK cells (cell data)" begin + vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells(node_coordinates) + end + # Determine output file names base, _ = splitext(splitdir(filename)[2]) @@ -64,35 +74,74 @@ function build_vtk_grids(::Val{:vti}, mesh, n_visnodes, verbose, # Open VTK files verbose && println("| Building VTK grid...") if is_datafile - # Determine level-wise resolution - max_level = maximum(levels) - resolution = n_visnodes * 2^max_level - - Nx = Ny = resolution + 1 - dx = dy = length_level_0/resolution - origin = center_level_0 .- 1/2 * length_level_0 - spacing = (dx, dy) - @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, Nx, Ny, - origin=tuple(origin...), - spacing=spacing) + @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, vtk_points, + vtk_cells) else vtk_nodedata = nothing end @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename, - vtk_celldata_points, - vtk_celldata_cells) + vtk_celldata_points, + vtk_celldata_cells) return vtk_nodedata, vtk_celldata end +# # Create and return VTK grids on enriched uniform nodes +# # that are ready to be filled with reinterpolated data (vti version) +# # TODO: Decide whether or not to keep this format. Appears to be broken... +# # OBS! Requires modifications to structure and docs here as well as the Trixi docs +# function build_vtk_grids(::Val{:vti}, mesh, nodes, n_visnodes, verbose, +# output_directory, is_datafile, filename, reinterpolate::Val{true}) + +# coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + +# # Extract number of spatial dimensions +# ndims_ = size(coordinates, 1) + +# # Prepare VTK points and cells for celldata file +# @timeit "prepare VTK cells" vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells( +# Val(ndims_), coordinates, levels, center_level_0, length_level_0, 1) + +# # Determine output file names +# base, _ = splitext(splitdir(filename)[2]) +# vtk_filename = joinpath(output_directory, base) +# vtk_celldata_filename = vtk_filename * "_celldata" + +# # Open VTK files +# verbose && println("| Building VTK grid...") +# if is_datafile +# # Determine level-wise resolution +# max_level = maximum(levels) +# resolution = n_visnodes * 2^max_level + +# Nx = Ny = resolution + 1 +# dx = dy = length_level_0/resolution +# origin = center_level_0 .- 1/2 * length_level_0 +# spacing = (dx, dy) +# @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, Nx, Ny, +# origin=tuple(origin...), +# spacing=spacing) +# else +# vtk_nodedata = nothing +# end +# @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename, +# vtk_celldata_points, +# vtk_celldata_cells) + +# return vtk_nodedata, vtk_celldata +# end + + # Create and return VTK grids that are ready to be filled with data -# (StructuredMesh/UnstructuredMesh2D/P4estMesh version) +# (StructuredMesh/UnstructuredMesh2D/P4estMesh version). +# Routine is agnostic with respect to reinterpolation. function build_vtk_grids(::Val{:vtu}, mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh}, - n_visnodes, verbose, output_directory, is_datafile, filename) + nodes, n_visnodes, verbose, output_directory, is_datafile, filename, + reinterpolate::Union{Val{true}, Val{false}}) - @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, n_visnodes) + @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, nodes, n_visnodes) # Calculate VTK points and cells verbose && println("| Preparing VTK cells...") @@ -128,12 +177,24 @@ function build_vtk_grids(::Val{:vtu}, end -function calc_node_coordinates(mesh::StructuredMesh, n_visnodes) +function calc_node_coordinates(mesh::TreeMesh, nodes, n_visnodes) + coordinates, levels, _, _ = extract_mesh_information(mesh) + + # Extract number of spatial dimensions + ndims_ = size(coordinates, 1) + n_elements = length(levels) + + node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., n_elements) + + return calc_node_coordinates!(node_coordinates, nodes, mesh) +end + + +function calc_node_coordinates(mesh::StructuredMesh, nodes, n_visnodes) # Extract number of spatial dimensions ndims_ = ndims(mesh) n_elements = prod(size(mesh)) - nodes = range(-1, 1, length=n_visnodes) f(args...; kwargs...) = Base.invokelatest(mesh.mapping, args...; kwargs...) node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., n_elements) @@ -142,14 +203,11 @@ function calc_node_coordinates(mesh::StructuredMesh, n_visnodes) end -function calc_node_coordinates(mesh::UnstructuredMesh2D, n_visnodes) +function calc_node_coordinates(mesh::UnstructuredMesh2D, nodes, n_visnodes) # Extract number of spatial dimensions ndims_ = ndims(mesh) n_elements = length(mesh) - # get the uniform nodes - nodes = range(-1, 1, length=n_visnodes) - # intialize the container for the node coordinates node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., n_elements) @@ -173,12 +231,10 @@ function calc_node_coordinates(mesh::UnstructuredMesh2D, n_visnodes) end -function calc_node_coordinates(mesh::P4estMesh, n_visnodes) +function calc_node_coordinates(mesh::P4estMesh, nodes, n_visnodes) # Extract number of spatial dimensions ndims_ = ndims(mesh) - nodes = range(-1, 1, length=n_visnodes) - node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., Trixi.ncells(mesh)) @@ -187,6 +243,96 @@ function calc_node_coordinates(mesh::P4estMesh, n_visnodes) end +# Calculation of the node coordinates for `TreeMesh` in 2D +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, nodes, mesh) + _, levels, _, length_level_0 = extract_mesh_information(mesh) + + # Extract number of spatial dimensions + n_elements = length(levels) + + # Set the reference length + # TODO: Is this ever different from 2? + reference_length = 2.0 + # Compute the offset of the midpoint of the 1D reference interval + # (its difference from zero) + reference_offset = first(nodes) + reference_length / 2 + + # Recompute the cell ids + cell_ids = Trixi.local_leaf_cells(mesh.tree) + + # Calculate inverse Jacobian and node coordinates + for element in 1:n_elements + # Get cell id + cell_id = cell_ids[element] + + # Get cell length + dx = length_level_0 / 2^levels[element] + + # Get Jacobian value + jacobian = dx / reference_length + + # Calculate node coordinates + # Note that the `tree_coordinates` are the midpoints of the cells. + # Hence, we need to add an offset for `nodes` with a midpoint + # different from zero. + for j in eachindex(nodes), i in eachindex(nodes) + node_coordinates[1, i, j, element] = ( + mesh.tree.coordinates[1, cell_id] + jacobian * (nodes[i] - reference_offset)) + node_coordinates[2, i, j, element] = ( + mesh.tree.coordinates[2, cell_id] + jacobian * (nodes[j] - reference_offset)) + end + end + + return node_coordinates +end + + +# Calculation of the node coordinates for `TreeMesh` in 3D +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, nodes, mesh) + _, levels, _, length_level_0 = extract_mesh_information(mesh) + + # Extract number of spatial dimensions + n_elements = length(levels) + + # Set the reference length + # TODO: Is this ever different from 2? + reference_length = 2.0 + # Compute the offset of the midpoint of the 1D reference interval + # (its difference from zero) + reference_offset = first(nodes) + reference_length / 2 + + # Recompute the cell ids + cell_ids = Trixi.local_leaf_cells(mesh.tree) + + # Calculate inverse Jacobian and node coordinates + for element in 1:n_elements + # Get cell id + cell_id = cell_ids[element] + + # Get cell length + dx = length_level_0 / 2^levels[element] + + # Get Jacobian value + jacobian = dx / reference_length + + # Calculate node coordinates + # Note that the `tree_coordinates` are the midpoints of the cells. + # Hence, we need to add an offset for `nodes` with a midpoint + # different from zero. + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + node_coordinates[1, i, j, k, element] = ( + mesh.tree.coordinates[1, cell_id] + jacobian * (nodes[i] - reference_offset)) + node_coordinates[2, i, j, k, element] = ( + mesh.tree.coordinates[2, cell_id] + jacobian * (nodes[j] - reference_offset)) + node_coordinates[3, i, j, k, element] = ( + mesh.tree.coordinates[3, cell_id] + jacobian * (nodes[k] - reference_offset)) + end + end + + return node_coordinates +end + + # Calculation of the node coordinates for `StructuredMesh` in 2D function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, f, nodes, mesh) linear_indices = LinearIndices(size(mesh)) diff --git a/test/Project.toml b/test/Project.toml index b2d49ec..05a0bad 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,11 +2,12 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" +ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] Documenter = "0.27" OrdinaryDiffEq = "6" +ReadVTK = "0.1" Trixi = "0.4" diff --git a/test/test_2d.jl b/test/test_2d.jl index a12ef06..9492dab 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -9,175 +9,383 @@ include("test_trixi2vtk.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) +# Windows github runners encounter memory issues saving all the output files. +# So artifacts are only save on Ubuntu and Mac # Create artifacts directory where all files that should be preserved will be stored artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") if !isdir(artifacts_dir) mkdir(artifacts_dir) end - @testset "2D" begin @testset "TreeMesh" begin isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_extended.jl"), maxiters=1) - - @timed_testset "uniform mesh" begin - if Sys.iswindows() - # This test fails on Windows due to globbing not working - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"), - ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")]) - outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu") - - else - test_trixi2vtk("solution_00000*.h5", outdir, - hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"), - ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830"), - ("solution_00000.pvd", "b9e6742dc2b397b14d8d3964e90dcfadea5d98cb"), - ("solution_00000_celldata.pvd", "c12004714a980581450cd4bad16a2541c5ec8f26")]) - outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu", - "solution_00000.pvd", "solution_00000_celldata.pvd") - end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - testname = "2d-tree-mesh-uniform-mesh" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) - end + run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave.jl"), maxiters=10) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000010.h5"), output_directory=outdir) + outfilename = "mesh_000010_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_amr_mesh_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_amr_mesh_10.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir) + outfilename = "solution_000010_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_amr_celldata_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_amr_celldata_10.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_amr_reinterp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_amr_reinterp_10.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_amr_no_reinterp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_amr_no_reinterp_10.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_amr_no_reinterp_uniform_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_amr_no_reinterp_uniform_10.vtu", remote_filename) + compare_point_info(out_file, ref_file) end - @timed_testset "uniform mesh with vti output" begin - if Sys.isapple() - # This test fails on MacOS due to differing binary VTI files (even though the contents match) - test_trixi2vtk("restart_000001.h5", outdir, - format=:vti) - else - test_trixi2vtk("restart_000001.h5", outdir, - hashes=[("restart_000001.vti", "9ade067c71f1f6492242a8aa215bd0d633caf9bc"), - ("restart_000001_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")], - format=:vti) - end - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("restart_000001.vti", "restart_000001_celldata.vtu") - testname = "2d-tree-mesh-uniform-mesh-with-vti-output" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + @timed_testset "attempt reinterpolate with uniform data" begin + # Purposely request a bad configuration and check that an error message gets thrown + # OBS! Only needs tested once across all mesh types and dimensions + let err = nothing + try + trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, data_is_uniform=true) + catch err + end + @test err isa Exception + @test sprint(showerror, err) == "uniform data should not be reinterpolated! Set reinterpolate=false and try again." end end end - @testset "StructuredMesh" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_waving_flag.jl"), maxiters=1) - - @timed_testset "waving flag" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "a93dbd393647627a861d890568e65598be0062f9")]) - else - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "564701ed0a9a90230f3a67f8bddd0616c818319b")]) - end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "2d-curved-mesh-waving-flag" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + + if !Sys.iswindows() + # OBS! Only `TreeMesh` results are tested on Windows runners due to memory limits. + # All remaining mesh types are tested on Ubuntu and Mac + @testset "StructuredMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_waving_flag.jl"), maxiters=1) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"), output_directory=outdir) + outfilename = "mesh_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-structuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/structuredmesh/dgsem_adv_mesh_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_mesh_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) end - end - @timed_testset "waving flag (supersampling)" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, nvisnodes=6, - hashes=[("solution_000000.vtu", "c6d74ab831bf4b6de2ba8cf537b6653ad611cfe7")]) - else - test_trixi2vtk("solution_000000.h5", outdir, nvisnodes=6, - hashes=[("solution_000000.vtu", "ae8c059c110aaabe2ed7dcfa8516d336c15ba618")]) - end - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "2d-curved-mesh-waving-flag-supersampling" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + @timed_testset "solution celldata" begin + # create the output file to be tested + # OBS! This exercises passing multiple files (in this case 2 files) to `trixi2vtk` + # that only needs to be tested once. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_00000*.h5"), output_directory=outdir) + + outfilename = "solution_000001_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-structuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/structuredmesh/dgsem_adv_celldata_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_celldata_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-structuredmesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/structuredmesh/dgsem_adv_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-structuredmesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/structuredmesh/dgsem_adv_no_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_no_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-structuredmesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/structuredmesh/dgsem_adv_no_reinterp_uniform_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_no_reinterp_uniform_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end - end - @testset "UnstructuredMesh2D" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "unstructured_2d_dgsem", "elixir_euler_basic.jl"), maxiters=1) - - @timed_testset "basic" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "0daedeea99d03d53b925ce5691bd7924abe88861")]) - else - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "acc03f295d2b7daee2cb0b4e29b90035c5e92bd7")]) - end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "2d-unstructured-quad-basic" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + @testset "UnstructuredMesh2D" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "unstructured_2d_dgsem", "elixir_shallowwater_ec.jl"), maxiters=1) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"), output_directory=outdir) + outfilename = "mesh_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-unstructuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/unstructuredmesh/dgsem_swe_mesh_01.vtu" + ref_file = get_test_reference_file("dgsem_swe_mesh_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir) + outfilename = "solution_000001_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-unstructuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/unstructuredmesh/dgsem_swe_celldata_01.vtu" + ref_file = get_test_reference_file("dgsem_swe_celldata_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-unstructuredmesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/unstructuredmesh/dgsem_swe_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_swe_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-unstructuredmesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/unstructuredmesh/dgsem_swe_no_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_swe_no_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-unstructuredmesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/unstructuredmesh/dgsem_swe_no_reinterp_uniform_01.vtu" + ref_file = get_test_reference_file("dgsem_swe_no_reinterp_uniform_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end - end - @testset "P4estMesh" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), initial_refinement_level=0, maxiters=1) - - @timed_testset "nonperiodic" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "476eaaa9c13c3a8dba826b614a68a9d3e7979e6b")]) - else - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "a80aadb353ce6ec40baa1b94d278c480f17d0419")]) - end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "2d-p4est-mesh-nonperiodic" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + @testset "P4estMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_mhd_rotor.jl"), maxiters=5) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000005.h5"), output_directory=outdir) + outfilename = "mesh_000005_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-p4estmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/p4estmesh/dgsem_rotor_amr_mesh_05.vtu" + ref_file = get_test_reference_file("dgsem_rotor_amr_mesh_05.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir) + outfilename = "solution_000005_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-p4estmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/p4estmesh/dgsem_rotor_amr_celldata_05.vtu" + ref_file = get_test_reference_file("dgsem_rotor_amr_celldata_05.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir) + outfilename = "solution_000005.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-p4estmesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/p4estmesh/dgsem_rotor_amr_reinterp_05.vtu" + ref_file = get_test_reference_file("dgsem_rotor_amr_reinterp_05.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000005.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-p4estmesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/p4estmesh/dgsem_rotor_amr_no_reinterp_05.vtu" + ref_file = get_test_reference_file("dgsem_rotor_amr_no_reinterp_05.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000005.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-p4estmesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/p4estmesh/dgsem_rotor_amr_no_reinterp_uniform_05.vtu" + ref_file = get_test_reference_file("dgsem_rotor_amr_no_reinterp_uniform_05.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end end end -# Clean up afterwards: delete Trixi output directory +# Clean up afterwards: delete Trixi output directory and reference file directory @test_nowarn rm(outdir, recursive=true) +@test_nowarn rm(TEST_REFERENCE_DIR, recursive=true) end diff --git a/test/test_3d.jl b/test/test_3d.jl index 76b9a68..5e288b6 100644 --- a/test/test_3d.jl +++ b/test/test_3d.jl @@ -9,93 +9,280 @@ include("test_trixi2vtk.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) +# Windows github runners encounter memory issues saving all the output files. +# So artifacts are only save on Ubuntu and Mac # Create artifacts directory where all files that should be preserved will be stored artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") if !isdir(artifacts_dir) mkdir(artifacts_dir) end - @testset "3D" begin @testset "TreeMesh" begin isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_extended.jl"), maxiters=1) - - @timed_testset "uniform mesh" begin - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "6ab3aa525851187ee0839e1d670a254a66be4ad7"), - ("solution_000000_celldata.vtu", "fdfee2d4200ecdad08067b37908412813016f4e7")]) - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu") - testname = "3d-tree-mesh-uniform-mesh" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) - end + run_trixi(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_blob_amr.jl"), maxiters=4) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000004.h5"), output_directory=outdir) + outfilename = "mesh_000004_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-treemesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/treemesh/dgsem_blob_amr_mesh_04.vtu" + ref_file = get_test_reference_file("dgsem_blob_amr_mesh_04.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir) + outfilename = "solution_000004_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-treemesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/treemesh/dgsem_blob_amr_celldata_04.vtu" + ref_file = get_test_reference_file("dgsem_blob_amr_celldata_04.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir) + outfilename = "solution_000004.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-treemesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/treemesh/dgsem_blob_amr_reinterp_04.vtu" + ref_file = get_test_reference_file("dgsem_blob_amr_reinterp_04.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000004.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-treemesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/treemesh/dgsem_blob_amr_no_reinterp_04.vtu" + ref_file = get_test_reference_file("dgsem_blob_amr_no_reinterp_04.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000004.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-treemesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/treemesh/dgsem_blob_amr_no_reinterp_uniform_04.vtu" + ref_file = get_test_reference_file("dgsem_blob_amr_no_reinterp_uniform_04.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end - @testset "StructuredMesh" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "structured_3d_dgsem", "elixir_advection_basic.jl"), maxiters=1) - - @timed_testset "basic" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "58e07f981fd6c005ea17e47054bd509c2c66d771")]) - else - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "58e07f981fd6c005ea17e47054bd509c2c66d771")]) + if !Sys.iswindows() + # OBS! Only `TreeMesh` results are tested on Windows runners due to memory limits. + # All remaining mesh types are tested on Ubuntu and Mac + @testset "StructuredMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "structured_3d_dgsem", "elixir_advection_nonperiodic_curved.jl"), maxiters=1) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"), output_directory=outdir) + outfilename = "mesh_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-structuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/structuredmesh/dgsem_adv_mesh_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_mesh_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "3d-curved-mesh-basic" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir) + outfilename = "solution_000001_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-structuredmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/structuredmesh/dgsem_adv_celldata_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_celldata_01.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-structuredmesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/structuredmesh/dgsem_adv_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-structuredmesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/structuredmesh/dgsem_adv_no_reinterp_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_no_reinterp_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000001.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-structuredmesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/structuredmesh/dgsem_adv_no_reinterp_uniform_01.vtu" + ref_file = get_test_reference_file("dgsem_adv_no_reinterp_uniform_01.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end - end - @testset "P4estMesh" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "p4est_3d_dgsem", "elixir_advection_amr_unstructured_curved.jl"), maxiters=1) - - @timed_testset "unstructured curved" begin - if Sys.isapple() - # This file has a different hash on macOS for some reason - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "fe0f45809ef6f3f0c1b7f5331198585f406923c9")]) - else - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "0fa5a099378d153aa3a1bb7dcf3559ea5d6bf9c5")]) + @testset "P4estMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "p4est_3d_dgsem", "elixir_advection_cubed_sphere.jl"), maxiters=2) + + @timed_testset "mesh data" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"), output_directory=outdir) + outfilename = "mesh_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-p4estmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/p4estmesh/dgsem_adv_sphere_mesh_02.vtu" + ref_file = get_test_reference_file("dgsem_adv_sphere_mesh_02.vtu", remote_filename) + compare_cell_info(out_file, ref_file) + end + + @timed_testset "solution celldata" begin + # create the output file to be tested + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir) + outfilename = "solution_000002_celldata.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-p4estmesh" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/p4estmesh/dgsem_adv_sphere_celldata_02.vtu" + ref_file = get_test_reference_file("dgsem_adv_sphere_celldata_02.vtu", remote_filename) + compare_cell_info(out_file, ref_file) end - @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5")) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu",) - testname = "3d-p4est-mesh-unstructured-curved" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + + @timed_testset "reinterpolate with nonuniform data" begin + # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`) + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir) + outfilename = "solution_000002.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-p4estmesh-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/p4estmesh/dgsem_adv_sphere_reinterp_02.vtu" + ref_file = get_test_reference_file("dgsem_adv_sphere_reinterp_02.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with nonuniform data" begin + # Create and test output without reinterpolation on LGL nodes + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000002.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-p4estmesh-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/p4estmesh/dgsem_adv_sphere_no_reinterp_02.vtu" + ref_file = get_test_reference_file("dgsem_adv_sphere_no_reinterp_02.vtu", remote_filename) + compare_point_info(out_file, ref_file) + end + + @timed_testset "do not reinterpolate with uniform data" begin + # Create and test output without reinterpolation on uniform nodes + # OBS! This is a dummy test just to exercise code. The resulting plot will look weird. + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true) + outfilename = "solution_000002.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "3d-p4estmesh-no-reinterp-uniform" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "3d/p4estmesh/dgsem_adv_sphere_no_reinterp_uniform_02.vtu" + ref_file = get_test_reference_file("dgsem_adv_sphere_no_reinterp_uniform_02.vtu", remote_filename) + compare_point_info(out_file, ref_file) end end end end -# Clean up afterwards: delete Trixi output directory -@test_skip rm(outdir, recursive=true) +# Clean up afterwards: delete Trixi output directory and reference file directory +@test_nowarn rm(outdir, recursive=true) +@test_nowarn rm(TEST_REFERENCE_DIR, recursive=true) end - diff --git a/test/test_manual.jl b/test/test_manual.jl index bb05478..a051ee4 100644 --- a/test/test_manual.jl +++ b/test/test_manual.jl @@ -39,6 +39,12 @@ isdir(outdir) && rm(outdir, recursive=true) @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"); output_directory=outdir) end + @testset "trixi2vtk set number of output nodes" begin + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000000.h5"); nvisnodes=0) + + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000000.h5"); nvisnodes=5) + end + @timed_testset "pvd_filenames" begin @test Trixi2Vtk.pvd_filenames("", "manual", "out") == (joinpath("out", "manual"), joinpath("out", "manual_celldata")) @test_throws ErrorException Trixi2Vtk.pvd_filenames(("a", "b"), nothing, "out") diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl index 6b1b3b9..34c2aab 100644 --- a/test/test_trixi2vtk.jl +++ b/test/test_trixi2vtk.jl @@ -1,7 +1,8 @@ using Test: @test_nowarn, @test, @testset, @test_skip -using SHA +using Downloads: download using Trixi using Trixi2Vtk +using ReadVTK function run_trixi(elixir; kwargs...) @@ -10,25 +11,153 @@ function run_trixi(elixir; kwargs...) end -function sha1file(filename) - open(filename) do f - bytes2hex(sha1(f)) +""" + get_reference_file(filename, remotename; head="main", output_directory=".", force=false) + +Retrieve an reference file from the +[`Trixi2Vtk_reference_files` repository](https://github.com/trixi-framework/Trixi2Vtk_reference_files) +at commit/branch `head` and store it in the `output_directory`. If the file already +exists locally, do not download the file again unless `force` is true. +Provide the `remotename` because the file structure in the `Trixi2Vtk_reference_files` repository +is different than what will be created locally. Return the local path to the downloaded file. +""" +function get_reference_file(filename, remotename; head="main", output_directory=".", force=false) + filepath = joinpath(output_directory, filename) + if !isfile(filepath) || force + url = ("https://github.com/trixi-framework/Trixi2Vtk_reference_files/raw/" + * head + * "/reference_files/" + * remotename) + download(url, filepath) end + + return filepath end +# Commit in the reference file repository for which the test files will be downloaded +# Note: The purpose of using a specific commit hash (instead of `main`) is to be able to tie a given +# version of Trixi2Vtk to a specific version of the test file repository. This way, also tests +# for older Trixi2Vtk releases should continue to work. +TEST_REFERENCE_COMMIT = "1fb95363322c32d2391267e06671273ad2889144" + +# Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`! +TEST_REFERENCE_DIR = "reference_files" + + +get_test_reference_file(filename, remotename) = get_reference_file(filename, remotename, + head=TEST_REFERENCE_COMMIT, + output_directory=TEST_REFERENCE_DIR) + + +# Start with a clean environment: remove reference file directory if it exists +isdir(TEST_REFERENCE_DIR) && rm(TEST_REFERENCE_DIR, recursive=true) +mkpath(TEST_REFERENCE_DIR) + + +""" + compare_cell_info(out_filename, ref_filename; atol=500*eps(), rtol=sqrt(eps())) + +Test values from the VTK file header and acutal (possibly interpolated) cell data. Uses +`out_filename` created during testing and compares against `ref_filename` that comes +from the +[`Trixi2Vtk_reference_files` repository](https://github.com/trixi-framework/Trixi2Vtk_reference_files). +""" +function compare_cell_info(out_filename, ref_filename; atol=500*eps(), rtol=sqrt(eps())) + ref_vtk = VTKFile(ref_filename) + vtk = VTKFile(out_filename) + + # check that the number of cells and points match + @test vtk.n_cells == ref_vtk.n_cells + @test vtk.n_points == ref_vtk.n_points + + # Compare data information + + # extract the data sets + out_cell_data = get_cell_data(vtk) + ref_cell_data = get_cell_data(ref_vtk) -function test_trixi2vtk(filenames, outdir; hashes=nothing, kwargs...) - @test_nowarn trixi2vtk(joinpath(outdir, filenames); output_directory=outdir, kwargs...) + # check that the number of variables and their names match + @test length(out_cell_data.names) == length(ref_cell_data.names) + @test out_cell_data.names == ref_cell_data.names - if !isnothing(hashes) - for (filename, hash_expected) in hashes - hash_measured = sha1file(joinpath(outdir, filename)) - @test_skip hash_expected == hash_measured - end + # Note!! + # Occassionally, for the last equation variable there is an issue + # that the size of the data extracted with + # [`ReadVTK.jl`](https://github.com/trixi-framework/ReadVTK.jl) + # is one byte larger than expected. Somehow this is related to + # the zlib format and/or how + # [`WriteVTK.jl`](https://github.com/jipolanco/WriteVTK.jl) created + # the file in the first place. Even though there is this access error, + # the `vtu` file is valid and can be used with ParaView or VisIt without issue. + # Therefore, we check all the variable data except for the last one. + + eqn_check_number = length(ref_cell_data.names) == 1 ? length(ref_cell_data.names) : length(ref_cell_data.names)-1 + + # check that the actual plot data is (approximately) the same + for i in 1:eqn_check_number + variable_name = ref_cell_data.names[i] + out_data = get_data(out_cell_data[variable_name]) + ref_data = get_data(ref_cell_data[variable_name]) + @test isapprox(out_data, ref_data, atol=atol, rtol=rtol) end end +""" + compare_point_info(out_filename, ref_filename; atol=500*eps(), rtol=sqrt(eps())) + +Test values from the VTK file header and acutal (possibly interpolated) point data. Uses +`out_filename` created during testing and compares against `ref_filename` that comes +from the +[`Trixi2Vtk_reference_files` repository](https://github.com/trixi-framework/Trixi2Vtk_reference_files). +""" +function compare_point_info(out_filename, ref_filename; atol=500*eps(), rtol=sqrt(eps())) + # Load the data from both files + ref_vtk = VTKFile(ref_filename) + vtk = VTKFile(out_filename) + + # Compare header information + @test vtk.byte_order == ref_vtk.byte_order + @test vtk.compressor == ref_vtk.compressor + @test vtk.file_type == ref_vtk.file_type + @test vtk.version == ref_vtk.version + + # check that the number of cells and points match + @test vtk.n_cells == ref_vtk.n_cells + @test vtk.n_points == ref_vtk.n_points + + # Compare data information + + # extract the data sets + out_point_data = get_point_data(vtk) + ref_point_data = get_point_data(ref_vtk) + + # check that the number of variables and their names match + @test length(out_point_data.names) == length(ref_point_data.names) + @test out_point_data.names == ref_point_data.names + + # Note!! + # Occassionally, for the last equation variable there is an issue + # that the size of the data extracted with + # [`ReadVTK.jl`](https://github.com/trixi-framework/ReadVTK.jl) + # is one byte larger than expected. Somehow this is related to + # the zlib format and/or how + # [`WriteVTK.jl`](https://github.com/jipolanco/WriteVTK.jl) created + # the file in the first place. Even though there is this access error, + # the `vtu` file is valid and can be used with ParaView or VisIt without issue. + # Therefore, we check all the variable data except for the last one. + + eqn_check_number = length(ref_point_data.names) == 1 ? length(ref_point_data.names) : length(ref_point_data.names)-1 + + # check that the actual plot data is (approximately) the same + for i in 1:eqn_check_number + variable_name = ref_point_data.names[i] + out_data = get_data(out_point_data[variable_name]) + ref_data = get_data(ref_point_data[variable_name]) + @test isapprox(out_data, ref_data, atol=atol, rtol=rtol) + end +end + """ @timed_testset "name of the testset" #= code to test #= From 2f0fb880789100ba09dcf0797200c86b893d1ac7 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 18 Jan 2023 10:21:02 +0100 Subject: [PATCH 04/21] Adapt to main --- src/convert.jl | 26 +++++++++++--------------- src/interpolate.jl | 2 -- test/test_2d.jl | 2 -- test/test_3d.jl | 2 -- test/test_trixi2vtk.jl | 2 +- 5 files changed, 12 insertions(+), 22 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 3524dfc..466c270 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -153,19 +153,6 @@ function trixi2vtk(filename::AbstractString...; node_set = Array{Float64}(undef, n_visnodes) end - # Check if the raw data is uniform (finite difference) or not (dg) - # and create the corresponding node set for reinterpolation / copying. - if (reinterpolate & !data_is_uniform) | (!reinterpolate & data_is_uniform) - # (1) Default settings; presumably the most common - # (2) Finite difference data - node_set = range(-1, 1, length=n_visnodes) - elseif !reinterpolate & !data_is_uniform - # raw data is on a set of LGL nodes - node_set, _ = gauss_lobatto_nodes_weights(n_visnodes) - else # reinterpolate & data_is_uniform - error("uniform data should not be reinterpolated! Set reinterpolate=false and try again.") - end - # Create output directory if it does not exist mkpath(output_directory) @@ -345,10 +332,19 @@ end # default number of visualization nodes if only the mesh should be visualized +function get_default_nvisnodes_mesh(nvisnodes, mesh::TreeMesh) + if nvisnodes === nothing + # for a Cartesian mesh, we do not need to interpolate + return 1 + else + return nvisnodes + end +end + function get_default_nvisnodes_mesh(nvisnodes, - mesh::Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh}) + mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh}) if nvisnodes === nothing - # we need to get at least the vertices + # for curved meshes, we need to get at least the vertices return 2 else return nvisnodes diff --git a/src/interpolate.jl b/src/interpolate.jl index 5225773..7187015 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -120,13 +120,11 @@ function unstructured2structured(unstructured_data::AbstractArray{Float64}, # Create output data structure structured = Array{Float64}(undef, resolution, resolution, n_variables) - # For each variable, interpolate element data and store to global data structure for v in 1:n_variables # Reshape data array for use in interpolate_nodes function reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements) - for element_id in 1:n_elements # Extract level for convenience level = levels[element_id] diff --git a/test/test_2d.jl b/test/test_2d.jl index b1e3e52..8c0a0ce 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -9,8 +9,6 @@ include("test_trixi2vtk.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -# Windows github runners encounter memory issues saving all the output files. -# So artifacts are only save on Ubuntu and Mac # Create artifacts directory where all files that should be preserved will be stored artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") if !isdir(artifacts_dir) diff --git a/test/test_3d.jl b/test/test_3d.jl index 06bf282..e7217a8 100644 --- a/test/test_3d.jl +++ b/test/test_3d.jl @@ -9,8 +9,6 @@ include("test_trixi2vtk.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -# Windows github runners encounter memory issues saving all the output files. -# So artifacts are only save on Ubuntu and Mac # Create artifacts directory where all files that should be preserved will be stored artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") if !isdir(artifacts_dir) diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl index efa7a90..b9aa8e6 100644 --- a/test/test_trixi2vtk.jl +++ b/test/test_trixi2vtk.jl @@ -1,4 +1,4 @@ -using Test: @test_nowarn, @test, @testset, @test_skip +using Test: @test_nowarn, @test, @testset using Downloads: download using Trixi using Trixi2Vtk From 0aee544cdecf7beeab8bfa04df660192691a190d Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 17 May 2023 14:46:13 +0200 Subject: [PATCH 05/21] Fix interpolation of node variables --- src/convert.jl | 10 +++++----- src/interpolate.jl | 43 ------------------------------------------- 2 files changed, 5 insertions(+), 48 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 57f6084..94a0bde 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -204,12 +204,12 @@ function trixi2vtk(filename::AbstractString...; for (label, variable) in node_variables verbose && println("| | Node variable: $label...") if reinterpolate - @timeit "interpolate cell data" interpolated_cell_data = interpolate_cell_data(Val(format), - variable, mesh, - n_visnodes, verbose) + @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format), + reshape(variable, size(variable)..., 1), + mesh, n_visnodes, verbose) else - @timeit "interpolate cell data" interpolated_cell_data = reshape(variable, - n_visnodes^ndims_ * n_elements) + @timeit "interpolate data" interpolated_cell_data = reshape(variable, + n_visnodes^ndims_ * n_elements) end # Add the "interpolated" cell_data to celldata, not node_data @timeit label vtk_nodedata[label] = interpolated_cell_data diff --git a/src/interpolate.jl b/src/interpolate.jl index 7187015..528ece6 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -8,49 +8,6 @@ function interpolate_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, v return raw2interpolated(input_data, nodes_out) end -function interpolate_cell_data(::Val{:vtu}, input_data, mesh::Union{TreeMesh, StructuredMesh}, n_visnodes, verbose) - # Calculate equidistant output nodes (visualization nodes) - nodes_out = collect(range(-1, 1, length=n_visnodes)) - - # Extract number of spatial dimensions - ndims_ = ndims(mesh) - - # Extract data shape information - n_nodes_in = size(input_data, 1) - n_nodes_out = length(nodes_out) - n_elements = size(input_data, ndims_ + 1) - - # Get node coordinates for DG locations on reference element - nodes_in, weights_in = gauss_lobatto_nodes_weights(n_nodes_in) - - if ndims_ == 2 - # Create output data structure - data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_elements) - - # Interpolate node data for visualization nodes to piecewise constant values and store to global data structure - for element_id in 1:n_elements - index_j = 1 - for j in 1:n_nodes_out - index_i = 1 - for i in 1:n_nodes_out - while -1.0 + sum(weights_in[1:index_i]) < nodes_out[i] - index_i += 1 - end - while -1.0 + sum(weights_in[1:index_j]) < nodes_out[j] - index_j += 1 - end - data_vis[i, j, element_id] = input_data[index_i, index_j, element_id] - end - end - end - else - error("Unsupported number of spatial dimensions: ", ndims_) - end - - # Return as one 1D array - return reshape(data_vis, n_nodes_out^ndims_ * n_elements) -end - # Interpolate data from input format to desired output format (StructuredMesh or UnstructuredMesh2D version) function interpolate_data(::Val{:vtu}, input_data, From 1e93ad3877670f444d5908ef0bb2614b72b431eb Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 21 Aug 2023 13:28:06 +0200 Subject: [PATCH 06/21] Add tests for subcell limiting coefficients --- test/test_2d.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_2d.jl b/test/test_2d.jl index d752644..31c5ffb 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -395,6 +395,18 @@ end end end end + + @testset "Subcell limiting coefficients" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_shockcapturing_subcell.jl"), maxiters=10) + + @timed_testset "do not reinterpolate" begin + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) + end + @timed_testset "do reinterpolate" begin + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) + end + end end # Clean up afterwards: delete Trixi output directory and reference file directory From cdc6d236d90576c33f6cdeb83060ffa96235652b Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 23 Aug 2023 14:24:37 +0200 Subject: [PATCH 07/21] Add data comparison for subcell shockcapturing to tests --- test/test_2d.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/test_2d.jl b/test/test_2d.jl index 31c5ffb..d85826c 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -401,10 +401,35 @@ end run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_shockcapturing_subcell.jl"), maxiters=10) @timed_testset "do not reinterpolate" begin + # Create and test output without reinterpolation @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-shockcapturing-subcell-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_blast_shockcapturing_subcell_no_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_blast_shockcapturing_subcell_no_interp_10.vtu", remote_filename) + compare_point_data(out_file, ref_file) end + @timed_testset "do reinterpolate" begin + # Create and test output without reinterpolation @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) + + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-shockcapturing-subcell-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_blast_shockcapturing_subcell_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_blast_shockcapturing_subcell_interp_10.vtu", remote_filename) + compare_point_data(out_file, ref_file) end end end From 2f4ec01b62616059cbfd02f45ebe3c45672c5335 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 25 Oct 2023 10:26:37 +0200 Subject: [PATCH 08/21] Add info when interpolating alphas --- src/convert.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/convert.jl b/src/convert.jl index 29c2d18..e4a3f80 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -207,6 +207,9 @@ function trixi2vtk(filename::AbstractString...; for (label, variable) in node_variables verbose && println("| | Node variable: $label...") if reinterpolate + @info "The limiting coefficients are no continuous field but happens to be" * + "represented by a piecewise-constant approximation. Thus, reinterpolation does" * + "not give a meaningful representation" @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format), reshape(variable, size(variable)..., 1), mesh, n_visnodes, verbose) From 2140c285bdffe083ec508c9a031e3f4732c17d28 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 25 Oct 2023 10:34:03 +0200 Subject: [PATCH 09/21] Use `println` instead of `info` to not interrupt the test --- src/convert.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index e4a3f80..267fa2b 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -207,9 +207,9 @@ function trixi2vtk(filename::AbstractString...; for (label, variable) in node_variables verbose && println("| | Node variable: $label...") if reinterpolate - @info "The limiting coefficients are no continuous field but happens to be" * + println("The limiting coefficients are no continuous field but happens to be" * "represented by a piecewise-constant approximation. Thus, reinterpolation does" * - "not give a meaningful representation" + "not give a meaningful representation") @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format), reshape(variable, size(variable)..., 1), mesh, n_visnodes, verbose) From 7cdb7dc52cf2807ab6e1e631102e5a23379dbcc4 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 6 Dec 2023 10:39:56 +0100 Subject: [PATCH 10/21] Print reinterpolate warning only once --- src/convert.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 267fa2b..2620af1 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -207,9 +207,11 @@ function trixi2vtk(filename::AbstractString...; for (label, variable) in node_variables verbose && println("| | Node variable: $label...") if reinterpolate - println("The limiting coefficients are no continuous field but happens to be" * - "represented by a piecewise-constant approximation. Thus, reinterpolation does" * - "not give a meaningful representation") + if index == 1 + println("WARNING: The limiting coefficients are no continuous field but happens " * + "to be represented by a piecewise-constant approximation. Thus, reinterpolation " * + "does not give a meaningful representation.") + end @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format), reshape(variable, size(variable)..., 1), mesh, n_visnodes, verbose) From 7a39d91efe252794d158b54c657f007d6447562d Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Tue, 23 Apr 2024 14:16:12 +0200 Subject: [PATCH 11/21] Adapt test structure --- test/test_2d.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_2d.jl b/test/test_2d.jl index d85826c..50655a2 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -398,7 +398,7 @@ end @testset "Subcell limiting coefficients" begin isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_shockcapturing_subcell.jl"), maxiters=10) + run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), maxiters=10) @timed_testset "do not reinterpolate" begin # Create and test output without reinterpolation @@ -411,8 +411,8 @@ end cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) # remote file path is actually a URL so it always has the same path structure - remote_filename = "2d/treemesh/dgsem_blast_shockcapturing_subcell_no_interp_10.vtu" - ref_file = get_test_reference_file("dgsem_blast_shockcapturing_subcell_no_interp_10.vtu", remote_filename) + remote_filename = "2d/treemesh/dgsem_sedov_subcell_no_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_subcell_no_interp_10.vtu", remote_filename) compare_point_data(out_file, ref_file) end @@ -427,9 +427,9 @@ end cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) # remote file path is actually a URL so it always has the same path structure - remote_filename = "2d/treemesh/dgsem_blast_shockcapturing_subcell_interp_10.vtu" - ref_file = get_test_reference_file("dgsem_blast_shockcapturing_subcell_interp_10.vtu", remote_filename) - compare_point_data(out_file, ref_file) + remote_filename = "2d/treemesh/dgsem_sedov_subcell_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_subcell_interp_10.vtu", remote_filename) + compare_cell_data(out_file, ref_file) end end end From 5a1ab1ca5e552dd241d7312ec8ab37f2fe24c8ee Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 24 Apr 2024 10:18:33 +0200 Subject: [PATCH 12/21] Update test set names --- test/test_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_2d.jl b/test/test_2d.jl index 50655a2..61aacab 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -400,7 +400,7 @@ end isdir(outdir) && rm(outdir, recursive=true) run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), maxiters=10) - @timed_testset "do not reinterpolate" begin + @timed_testset "without reinterpolation" begin # Create and test output without reinterpolation @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) outfilename = "solution_000010.vtu" @@ -416,7 +416,7 @@ end compare_point_data(out_file, ref_file) end - @timed_testset "do reinterpolate" begin + @timed_testset "with reinterpolation" begin # Create and test output without reinterpolation @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) outfilename = "solution_000010.vtu" From 526b7cec93696b3416d83e173da4b27a0e8a91b1 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Fri, 21 Jun 2024 11:34:21 +0200 Subject: [PATCH 13/21] Adapt comment (suggestion); Adapt commit hash of repo with reference files --- src/convert.jl | 2 +- test/test_trixi2vtk.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 2620af1..ea6b66d 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -219,7 +219,7 @@ function trixi2vtk(filename::AbstractString...; @timeit "interpolate data" interpolated_cell_data = reshape(variable, n_visnodes^ndims_ * n_elements) end - # Add the "interpolated" cell_data to celldata, not node_data + # Add to node_data @timeit label vtk_nodedata[label] = interpolated_cell_data end end diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl index 15e59ad..40e3dc5 100644 --- a/test/test_trixi2vtk.jl +++ b/test/test_trixi2vtk.jl @@ -41,7 +41,7 @@ end if VERSION < v"1.8" const TEST_REFERENCE_COMMIT = "c0a966b06489f9b2ee3aefeb0a5c0dae733df36f" else - const TEST_REFERENCE_COMMIT = "86a43fe8dc254130345754fb512268204cf2233c" + const TEST_REFERENCE_COMMIT = "114c4d4dbb428ca31b6127935c2b30870436c605" end # Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`! From aeb6802cd7a8c58db19b8fc3f85afa6ec686ca9e Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Fri, 21 Jun 2024 15:00:26 +0200 Subject: [PATCH 14/21] Adapt if clause for printing warning --- src/convert.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/convert.jl b/src/convert.jl index ea6b66d..1e0f0e2 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -86,6 +86,10 @@ function trixi2vtk(filename::AbstractString...; barlen = 40) end + # Show warning when reinterpolating node-level data of subcell limiting + # Auxiliary variable to show warning only once + has_warned_about_interpolation = false + # Iterate over input files for (index, filename) in enumerate(filenames) verbose && println("Processing file $filename ($(index)/$(length(filenames)))...") @@ -207,10 +211,12 @@ function trixi2vtk(filename::AbstractString...; for (label, variable) in node_variables verbose && println("| | Node variable: $label...") if reinterpolate - if index == 1 + # Show warning if node-level data of subcell limiting are reinterpolated. + if label == "limiting_coefficient" && !has_warned_about_interpolation println("WARNING: The limiting coefficients are no continuous field but happens " * "to be represented by a piecewise-constant approximation. Thus, reinterpolation " * "does not give a meaningful representation.") + has_warned_about_interpolation = true end @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format), reshape(variable, size(variable)..., 1), From 458b0fb1cb2b8471f103f1f8e92e37d7fe2ad2bf Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 24 Jun 2024 10:42:31 +0200 Subject: [PATCH 15/21] Subcell tests only for julia v1.8+ --- test/test_2d.jl | 58 +++++++++++++++++++++++++------------------------ 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/test/test_2d.jl b/test/test_2d.jl index 61aacab..7ba0211 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -396,40 +396,42 @@ end end end - @testset "Subcell limiting coefficients" begin - isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), maxiters=10) + if VERSION >= v"1.8" + @testset "Subcell limiting coefficients" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), maxiters=10) - @timed_testset "without reinterpolation" begin - # Create and test output without reinterpolation - @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) - outfilename = "solution_000010.vtu" - out_file = joinpath(outdir, outfilename) + @timed_testset "without reinterpolation" begin + # Create and test output without reinterpolation + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) - # save output file to `artifacts` to facilitate debugging of failing tests - testname = "2d-treemesh-shockcapturing-subcell-no-reinterp" - cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-shockcapturing-subcell-no-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) - # remote file path is actually a URL so it always has the same path structure - remote_filename = "2d/treemesh/dgsem_sedov_subcell_no_interp_10.vtu" - ref_file = get_test_reference_file("dgsem_sedov_subcell_no_interp_10.vtu", remote_filename) - compare_point_data(out_file, ref_file) - end + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_subcell_no_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_subcell_no_interp_10.vtu", remote_filename) + compare_point_data(out_file, ref_file) + end - @timed_testset "with reinterpolation" begin - # Create and test output without reinterpolation - @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) - outfilename = "solution_000010.vtu" - out_file = joinpath(outdir, outfilename) + @timed_testset "with reinterpolation" begin + # Create and test output without reinterpolation + @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) + outfilename = "solution_000010.vtu" + out_file = joinpath(outdir, outfilename) - # save output file to `artifacts` to facilitate debugging of failing tests - testname = "2d-treemesh-shockcapturing-subcell-reinterp" - cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) + # save output file to `artifacts` to facilitate debugging of failing tests + testname = "2d-treemesh-shockcapturing-subcell-reinterp" + cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true) - # remote file path is actually a URL so it always has the same path structure - remote_filename = "2d/treemesh/dgsem_sedov_subcell_interp_10.vtu" - ref_file = get_test_reference_file("dgsem_sedov_subcell_interp_10.vtu", remote_filename) - compare_cell_data(out_file, ref_file) + # remote file path is actually a URL so it always has the same path structure + remote_filename = "2d/treemesh/dgsem_sedov_subcell_interp_10.vtu" + ref_file = get_test_reference_file("dgsem_sedov_subcell_interp_10.vtu", remote_filename) + compare_cell_data(out_file, ref_file) + end end end end From 7dec585a60857b55b88023f54560c2ea2eb3b5af Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 24 Jun 2024 10:50:27 +0200 Subject: [PATCH 16/21] Update hash of newest commit --- test/test_trixi2vtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl index 40e3dc5..dc3e867 100644 --- a/test/test_trixi2vtk.jl +++ b/test/test_trixi2vtk.jl @@ -41,7 +41,7 @@ end if VERSION < v"1.8" const TEST_REFERENCE_COMMIT = "c0a966b06489f9b2ee3aefeb0a5c0dae733df36f" else - const TEST_REFERENCE_COMMIT = "114c4d4dbb428ca31b6127935c2b30870436c605" + const TEST_REFERENCE_COMMIT = "ece6f552cf759f5b0f1416f108a15baedc3ad887" end # Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`! From cc95a2f8c6267b3d93b5e828db99fbd802be2951 Mon Sep 17 00:00:00 2001 From: bbolm Date: Thu, 10 Oct 2024 14:08:05 +0200 Subject: [PATCH 17/21] Add leading zeros to test filenames --- test/test_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_2d.jl b/test/test_2d.jl index 97db5d7..879008d 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -403,8 +403,8 @@ end @timed_testset "without reinterpolation" begin # Create and test output without reinterpolation - @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false) - outfilename = "solution_000010.vtu" + @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=false) + outfilename = "solution_" * LEADING_ZEROS * "000010.vtu" out_file = joinpath(outdir, outfilename) # save output file to `artifacts` to facilitate debugging of failing tests @@ -419,8 +419,8 @@ end @timed_testset "with reinterpolation" begin # Create and test output without reinterpolation - @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=true) - outfilename = "solution_000010.vtu" + @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=true) + outfilename = "solution_" * LEADING_ZEROS * "000010.vtu" out_file = joinpath(outdir, outfilename) # save output file to `artifacts` to facilitate debugging of failing tests From 61ece047f017be14d0edd4b944f072ce6f81e053 Mon Sep 17 00:00:00 2001 From: bbolm Date: Thu, 10 Oct 2024 14:26:26 +0200 Subject: [PATCH 18/21] Reduce initial refinement level in test --- test/test_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_2d.jl b/test/test_2d.jl index 879008d..d6f0519 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -399,7 +399,8 @@ end if VERSION >= v"1.8" @testset "Subcell limiting coefficients" begin isdir(outdir) && rm(outdir, recursive=true) - run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), maxiters=10) + run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + maxiters=10, initial_refinement_level=4) @timed_testset "without reinterpolation" begin # Create and test output without reinterpolation From b4f61abf43c4c8151c8a92d6818226f78159400b Mon Sep 17 00:00:00 2001 From: bbolm Date: Wed, 16 Oct 2024 10:47:13 +0200 Subject: [PATCH 19/21] Update commit hash for RefRepo --- test/test_trixi2vtk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl index be234db..3729341 100644 --- a/test/test_trixi2vtk.jl +++ b/test/test_trixi2vtk.jl @@ -56,7 +56,7 @@ end if VERSION < v"1.8" const TEST_REFERENCE_COMMIT = "c0a966b06489f9b2ee3aefeb0a5c0dae733df36f" else - const TEST_REFERENCE_COMMIT = "ece6f552cf759f5b0f1416f108a15baedc3ad887" + const TEST_REFERENCE_COMMIT = "7bb59953a7f2a6511fede9549fe8b062011ed788" end # Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`! From 8d7f021333ac146b35af39d75a8d2607eb9e7631 Mon Sep 17 00:00:00 2001 From: bbolm Date: Thu, 17 Oct 2024 16:10:18 +0200 Subject: [PATCH 20/21] Run macOS tests on apple chip `arm64` --- .github/workflows/ci.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 45a7f91..0e48049 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,11 +37,13 @@ jobs: # - 'nightly' os: - ubuntu-latest - - macOS-latest - windows-latest arch: - x64 include: + - version: '1.9' + os: macOS-latest + arch: arm64 - version: '1.7' os: ubuntu-latest arch: x64 From 48ec21317ce6a0a9c42f5b40c13f9c9d98710ff7 Mon Sep 17 00:00:00 2001 From: bbolm Date: Mon, 28 Oct 2024 11:12:17 +0100 Subject: [PATCH 21/21] Add comment about julia version --- test/test_2d.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_2d.jl b/test/test_2d.jl index d6f0519..dd397b2 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -397,6 +397,8 @@ end end if VERSION >= v"1.8" + # Julia v1.7 heavily downgrades Trixi.jl. Subcell limiting is not yet supported. + # Therefore, only perform tests with Julia v1.8 or newer. @testset "Subcell limiting coefficients" begin isdir(outdir) && rm(outdir, recursive=true) run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"),