Merge pull request #172 from ranocha/calc_error

extend calc_error_norms(func, dg, t)

Author: ranocha

src/equations/equations.jl

@@ -80,6 +80,9 @@ DG method (except floating point errors).
 end
 
 
+@inline cons2cons(u, equation::AbstractEquation) = u
+
+
 ####################################################################################################
 # Include files with actual implementations for different systems of equations.
 

src/solvers/dg/2d/dg.jl

@@ -706,7 +706,7 @@ integrate(u, dg::Dg2D; normalize=true) = integrate(identity, u, dg; normalize=no
 
 
 # Calculate L2/Linf error norms based on "exact solution"
-function calc_error_norms(dg::Dg2D, t)
+function calc_error_norms(func, dg::Dg2D, t)
   # Gather necessary information
   equation = equations(dg)
   n_nodes_analysis = size(dg.analysis_vandermonde, 1)
@@ -722,8 +722,8 @@ function calc_error_norms(dg::Dg2D, t)
             2, size(dg.analysis_vandermonde, 1), size(dg.analysis_vandermonde, 2))
 
   # Set up data structures
-  l2_error   = @SVector zeros(nvariables(equation))
-  linf_error = @SVector zeros(nvariables(equation))
+  l2_error   = zero(func(get_node_vars(dg.elements.u, dg, 1, 1, 1), equation))
+  linf_error = zero(func(get_node_vars(dg.elements.u, dg, 1, 1, 1), equation))
 
   # Iterate over all elements for error calculations
   for element_id in 1:dg.n_elements
@@ -736,7 +736,7 @@ function calc_error_norms(dg::Dg2D, t)
     jacobian_volume = inv(dg.elements.inverse_jacobian[element_id])^ndims(dg)
     for j in 1:n_nodes_analysis, i in 1:n_nodes_analysis
       u_exact = dg.initial_conditions(get_node_coords(x, dg, i, j), t, equation)
-      diff = u_exact - get_node_vars(u, dg, i, j)
+      diff = func(u_exact, equation) - func(get_node_vars(u, dg, i, j), equation)
       l2_error += diff.^2 * (weights[i] * weights[j] * jacobian_volume)
       linf_error = @. max(linf_error, abs(diff))
     end
@@ -920,6 +920,37 @@ function analyze_solution(dg::Dg2D, mesh::TreeMesh, time::Real, dt::Real, step::
     println()
   end
 
+  # L2/L∞ errors of the primitive variables
+  if :l2_error_primitive in dg.analysis_quantities || :linf_error_primitive in dg.analysis_quantities
+    l2_error_prim, linf_error_prim = calc_error_norms(cons2prim, dg, time)
+
+    print(" Variable:    ")
+    for v in 1:nvariables(equation)
+      @printf("   %-14s", varnames_prim(equation)[v])
+    end
+    println()
+
+    # L2 error
+    if :l2_error_primitive in dg.analysis_quantities
+      print(" L2 error prim.: ")
+      for v in 1:nvariables(equation)
+        @printf("%10.8e   ", l2_error_prim[v])
+        dg.save_analysis && @printf(f, "  % 10.8e", l2_error_prim[v])
+      end
+      println()
+    end
+
+    # L∞ error
+    if :linf_error_primitive in dg.analysis_quantities
+      print(" Linf error pri.:")
+      for v in 1:nvariables(equation)
+        @printf("%10.8e   ", linf_error_prim[v])
+        dg.save_analysis && @printf(f, "  % 10.8e", linf_error_prim[v])
+      end
+      println()
+    end
+  end
+
   # Entropy time derivative
   if :dsdu_ut in dg.analysis_quantities
     duds_ut = calc_entropy_timederivative(dg, time)
@@ -1088,6 +1119,16 @@ function save_analysis_header(filename, quantities, equation::AbstractEquation{2
         @printf(f, "   %-14s", "res_" * v)
       end
     end
+    if :l2_error_primitive in quantities
+      for v in varnames_prim(equation)
+        @printf(f, "   %-14s", "l2_" * v)
+      end
+    end
+    if :linf_error_primitive in quantities
+      for v in varnames_prim(equation)
+        @printf(f, "   %-14s", "linf_" * v)
+      end
+    end
     if :dsdu_ut in quantities
       @printf(f, "   %-14s", "dsdu_ut")
     end

src/solvers/dg/3d/dg.jl

@@ -750,7 +750,7 @@ integrate(u, dg::Dg3D; normalize=true) = integrate(identity, u, dg; normalize=no
 
 
 # Calculate L2/Linf error norms based on "exact solution"
-function calc_error_norms(dg::Dg3D, t)
+function calc_error_norms(func, dg::Dg3D, t)
   # Gather necessary information
   equation = equations(dg)
   n_nodes_analysis = size(dg.analysis_vandermonde, 1)
@@ -770,8 +770,8 @@ function calc_error_norms(dg::Dg3D, t)
             3, size(dg.analysis_vandermonde, 1), size(dg.analysis_vandermonde, 1), size(dg.analysis_vandermonde, 2))
 
   # Set up data structures
-  l2_error   = @SVector zeros(nvariables(equation))
-  linf_error = @SVector zeros(nvariables(equation))
+  l2_error   = zero(func(get_node_vars(dg.elements.u, dg, 1, 1, 1, 1), equation))
+  linf_error = zero(func(get_node_vars(dg.elements.u, dg, 1, 1, 1, 1), equation))
 
   # Iterate over all elements for error calculations
   for element_id in 1:dg.n_elements
@@ -784,7 +784,7 @@ function calc_error_norms(dg::Dg3D, t)
     jacobian_volume = inv(dg.elements.inverse_jacobian[element_id])^ndims(dg)
     for k in 1:n_nodes_analysis, j in 1:n_nodes_analysis, i in 1:n_nodes_analysis
       u_exact = dg.initial_conditions(get_node_coords(x, dg, i, j, k), t, equation)
-      diff = u_exact - get_node_vars(u, dg, i, j, k)
+      diff = func(u_exact, equation) - func(get_node_vars(u, dg, i, j, k), equation)
       l2_error += diff.^2 * (weights[i] * weights[j] * weights[k] * jacobian_volume)
       linf_error = @. max(linf_error, abs(diff))
     end
@@ -968,7 +968,7 @@ function analyze_solution(dg::Dg3D, mesh::TreeMesh, time::Real, dt::Real, step::
   end
 
   # Calculate L2/Linf errors, which are also returned by analyze_solution
-   l2_error, linf_error = calc_error_norms(dg, time)
+  l2_error, linf_error = calc_error_norms(dg, time)
 
   # L2 error
   if :l2_error in dg.analysis_quantities
@@ -1022,6 +1022,37 @@ function analyze_solution(dg::Dg3D, mesh::TreeMesh, time::Real, dt::Real, step::
     println()
   end
 
+  # L2/L∞ errors of the primitive variables
+  if :l2_error_primitive in dg.analysis_quantities || :linf_error_primitive in dg.analysis_quantities
+    l2_error_prim, linf_error_prim = calc_error_norms(cons2prim, dg, time)
+
+    print(" Variable:    ")
+    for v in 1:nvariables(equation)
+      @printf("   %-14s", varnames_prim(equation)[v])
+    end
+    println()
+
+    # L2 error
+    if :l2_error_primitive in dg.analysis_quantities
+      print(" L2 error prim.: ")
+      for v in 1:nvariables(equation)
+        @printf("%10.8e   ", l2_error_prim[v])
+        dg.save_analysis && @printf(f, "  % 10.8e", l2_error_prim[v])
+      end
+      println()
+    end
+
+    # L∞ error
+    if :linf_error_primitive in dg.analysis_quantities
+      print(" Linf error pri.:")
+      for v in 1:nvariables(equation)
+        @printf("%10.8e   ", linf_error_prim[v])
+        dg.save_analysis && @printf(f, "  % 10.8e", linf_error_prim[v])
+      end
+      println()
+    end
+  end
+
   # Entropy time derivative
   if :dsdu_ut in dg.analysis_quantities
     duds_ut = calc_entropy_timederivative(dg, time)
@@ -1181,6 +1212,16 @@ function save_analysis_header(filename, quantities, equation::AbstractEquation{3
         @printf(f, "   %-14s", "res_" * v)
       end
     end
+    if :l2_error_primitive in quantities
+      for v in varnames_prim(equation)
+        @printf(f, "   %-14s", "l2_" * v)
+      end
+    end
+    if :linf_error_primitive in quantities
+      for v in varnames_prim(equation)
+        @printf(f, "   %-14s", "linf_" * v)
+      end
+    end
     if :dsdu_ut in quantities
       @printf(f, "   %-14s", "dsdu_ut")
     end

src/solvers/dg/dg.jl

@@ -53,6 +53,7 @@ end
   return nothing
 end
 
+
 # Include utilities
 include("interpolation.jl")
 include("l2projection.jl")

src/solvers/solvers.jl

@@ -38,6 +38,17 @@ function make_solver(name::String, equations::AbstractEquation, mesh::TreeMesh;
 end
 
 
+"""
+    calc_error_norms([func=(u,equation)->u,] solver, t)
+
+Calculate the discrete L2 and L∞ errors of `func` applied to the conservative variables of
+the problem encapsulated by `solver` at time `t`, where `func` is called as `func(u, equation)`.
+"""
+function calc_error_norms end
+
+@inline calc_error_norms(solver::AbstractSolver, t) = calc_error_norms(cons2cons, solver, t)
+
+
 ####################################################################################################
 # Include files with actual implementations for different systems of equations.
 

test/test_examples_2d.jl

@@ -75,6 +75,12 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "
             l2   = [0.06159341742582756, 0.05012484425381723, 0.05013298724507752, 0.22537740506116724],
             linf = [0.29912627861573327, 0.30886767304359375, 0.3088108573487326, 1.0657556075017878])
   end
+  @testset "parameters_density_wave.toml" begin
+    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_density_wave.toml"),
+            l2   = [0.001060077845747576, 0.00010600778457107525, 0.00021201556914875742, 2.6501946139091318e-5],
+            linf = [0.0065356386867677085, 0.0006535638688170142, 0.0013071277374487877, 0.0001633909674296774],
+            extra_analysis_quantities=["l2_error_primitive", "linf_error_primitive"], t_end=0.5)
+  end
   @testset "parameters_ec_mhd.toml" begin
     test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec_mhd.toml"),
             l2   = [0.03607862694368351, 0.04281395008247395, 0.04280207686965749, 0.025746770192645763, 0.1611518499414067, 0.017455917249117023, 0.017456981264942977, 0.02688321120361229, 0.00015024027267648003],

test/test_examples_3d.jl

@@ -120,7 +120,8 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "
     test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_density_pulse.toml"),
             l2   = [0.05719652660597408, 0.0571965266059741, 0.05719652660597407, 0.05719652660597409, 0.08579478990896279],
             linf = [0.27375961853433606, 0.27375961853433517, 0.27375961853433384, 0.2737596185343343, 0.4106394278015033],
-            source_terms=Trixi.source_terms_harmonic)
+            source_terms=Trixi.source_terms_harmonic,
+            extra_analysis_quantities=["l2_error_primitive", "linf_error_primitive"])
   end
   @testset "parameters_taylor_green_vortex.toml" begin
     test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_taylor_green_vortex.toml"),