Add capability to handle with node variables (#67)

* Add node variables

* Add StructuredMesh

* Add branch to allow reinterpolate=false

* Adapt to main

* Fix interpolation of node variables

* Add tests for subcell limiting coefficients

* Add data comparison for subcell shockcapturing to tests

* Add info when interpolating alphas

* Use `println` instead of `info` to not interrupt the test

* Print reinterpolate warning only once

* Adapt test structure

* Update test set names

* Adapt comment (suggestion); Adapt commit hash of repo with reference files

* Adapt if clause for printing warning

* Subcell tests only for julia v1.8+

* Update hash of newest commit

* Add leading zeros to test filenames

* Reduce initial refinement level in test

* Update commit hash for RefRepo

* Run macOS tests on apple chip `arm64`

* Add comment about julia version

.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

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)))...")
@@ -124,7 +128,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)
 
@@ -202,6 +206,28 @@ 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...")
+            if reinterpolate
+              # 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),
+                                                                    mesh, n_visnodes, verbose)
+            else
+              @timeit "interpolate data" interpolated_cell_data = reshape(variable,
+                                                                          n_visnodes^ndims_ * n_elements)
+            end
+            # Add to node_data
+            @timeit label vtk_nodedata[label] = interpolated_cell_data
+          end
         end
       end
     end

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

test/test_2d.jl

@@ -395,6 +395,48 @@ end
       end
     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"),
+                maxiters=10, initial_refinement_level=4)
+
+      @timed_testset "without reinterpolation" begin
+        # Create and test output without reinterpolation
+        @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
+        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
+
+      @timed_testset "with reinterpolation" begin
+        # Create and test output without reinterpolation
+        @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
+        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)
+      end
+    end
+  end
 end
 
 # Clean up afterwards: delete Trixi output directory and reference file directory

test/test_trixi2vtk.jl

@@ -56,7 +56,7 @@ end
 if VERSION < v"1.8"
   const TEST_REFERENCE_COMMIT = "c0a966b06489f9b2ee3aefeb0a5c0dae733df36f"
 else
-  const TEST_REFERENCE_COMMIT = "86a43fe8dc254130345754fb512268204cf2233c"
+  const TEST_REFERENCE_COMMIT = "7bb59953a7f2a6511fede9549fe8b062011ed788"
 end
 
 # Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`!