L1 errors

examples/dgmulti_2d/elixir_navierstokes_convergence.jl

@@ -194,7 +194,7 @@ summary_callback = SummaryCallback()
 alive_callback = AliveCallback(alive_interval=10)
 analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
-callbacks = CallbackSet(summary_callback, alive_callback)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 
 ###############################################################################
 # run the simulation

examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl

@@ -62,7 +62,7 @@ summary_callback = SummaryCallback()
 alive_callback = AliveCallback(alive_interval=10)
 analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
-callbacks = CallbackSet(summary_callback, alive_callback)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 
 ###############################################################################
 # run the simulation

examples/structured_2d_dgsem/elixir_advection_extended.jl

@@ -48,6 +48,7 @@ summary_callback = SummaryCallback()
 # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
 analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:l1_error,),
                                      extra_analysis_integrals=(entropy, energy_total))
 
 # The AliveCallback prints short status information in regular intervals

src/callbacks_step/analysis.jl

@@ -149,6 +149,11 @@ function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrat
           @printf(io, "   %-14s", "linf_" * v)
         end
       end
+      if :l1_error in analysis_errors
+        for v in varnames(cons2cons, equations)
+          @printf(io, "   %-14s", "l1_" * v)
+        end
+      end
       if :conservation_error in analysis_errors
         for v in varnames(cons2cons, equations)
           @printf(io, "   %-14s", "cons_" * v)
@@ -169,6 +174,11 @@ function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrat
           @printf(io, "   %-14s", "linf_" * v)
         end
       end
+      if :l1_error_primitive in analysis_errors
+        for v in varnames(cons2prim, equations)
+          @printf(io, "   %-14s", "l1_" * v)
+        end
+      end
 
       for quantity in analysis_integrals
         @printf(io, "   %-14s", pretty_form_ascii(quantity))
@@ -292,7 +302,7 @@ function (analysis_callback::AnalysisCallback)(integrator)
     @notimeit timer() integrator.f(du_ode, integrator.u, semi, t)
     u  = wrap_array(integrator.u, mesh, equations, solver, cache)
     du = wrap_array(du_ode,       mesh, equations, solver, cache)
-    l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)
+    l2_error, linf_error, l1_error = analysis_callback(io, du, u, integrator.u, t, semi)
 
     mpi_println("─"^100)
     mpi_println()
@@ -315,7 +325,7 @@ function (analysis_callback::AnalysisCallback)(integrator)
   analysis_callback.ncalls_rhs_last_analysis = ncalls(semi.performance_counter)
 
   # Return errors for EOC analysis
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end
 
 
@@ -329,16 +339,16 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
   # Calculate and print derived quantities (error norms, entropy etc.)
   # Variable names required for L2 error, Linf error, and conservation error
   if any(q in analysis_errors for q in
-         (:l2_error, :linf_error, :conservation_error, :residual)) && mpi_isroot()
+         (:l2_error, :linf_error, :l1_error, :conservation_error, :residual)) && mpi_isroot()
     print(" Variable:    ")
     for v in eachvariable(equations)
       @printf("   %-14s", varnames(cons2cons, equations)[v])
     end
     println()
   end
 
-  # Calculate L2/Linf errors, which are also returned
-  l2_error, linf_error = calc_error_norms(u_ode, t, analyzer, semi, cache_analysis)
+  # Calculate L2, L1, Linf errors, which are also returned
+  l2_error, linf_error, l1_error = calc_error_norms(u_ode, t, analyzer, semi, cache_analysis)
 
   if mpi_isroot()
     # L2 error
@@ -360,6 +370,16 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
       end
       println()
     end
+
+    # L1 error
+    if :l1_error in analysis_errors
+      print(" L1 error:    ")
+      for v in eachvariable(equations)
+        @printf("  % 10.8e", l1_error[v])
+        @printf(io, "  % 10.8e", l1_error[v])
+      end
+      println()
+    end
   end
 
 
@@ -400,9 +420,10 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
     mpi_println()
   end
 
-  # L2/L∞ errors of the primitive variables
-  if :l2_error_primitive in analysis_errors || :linf_error_primitive in analysis_errors
-    l2_error_prim, linf_error_prim = calc_error_norms(cons2prim, u_ode, t, analyzer, semi, cache_analysis)
+  # L2, L1, L∞ errors of the primitive variables
+  if :l2_error_primitive in analysis_errors || :linf_error_primitive in analysis_errors || 
+    :l1_error_primitive in analysis_errors
+    l2_error_prim, linf_error_prim, l1_error_prim = calc_error_norms(cons2prim, u_ode, t, analyzer, semi, cache_analysis)
 
     if mpi_isroot()
       print(" Variable:    ")
@@ -430,13 +451,23 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
         end
         println()
       end
+
+      # L1 error
+      if :l1_error_primitive in analysis_errors
+        print(" L1 error prim.: ")
+        for v in eachvariable(equations)
+          @printf("%10.8e   ", l1_error_prim[v])
+          @printf(io, "  % 10.8e", l1_error_prim[v])
+        end
+        println()
+      end
     end
   end
 
   # additional integrals
   analyze_integrals(analysis_integrals, io, du, u, t, semi)
 
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end
 
 
@@ -530,8 +561,8 @@ function (cb::DiscreteCallback{Condition,Affect!})(sol) where {Condition, Affect
   @unpack analyzer = analysis_callback
   cache_analysis = analysis_callback.cache
 
-  l2_error, linf_error = calc_error_norms(sol.u[end], sol.t[end], analyzer, semi, cache_analysis)
-  (; l2=l2_error, linf=linf_error)
+  l2_error, linf_error, l1_error = calc_error_norms(sol.u[end], sol.t[end], analyzer, semi, cache_analysis)
+  (; l2=l2_error, linf=linf_error, l1=l1_error)
 end
 
 

src/callbacks_step/analysis_dg1d.jl

@@ -51,6 +51,7 @@ function calc_error_norms(func, u, t, analyzer,
   # Set up data structures
   l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
   linf_error = copy(l2_error)
+  l1_error   = copy(l2_error)
   total_volume = zero(real(mesh))
 
   # Iterate over all elements for error calculations
@@ -66,16 +67,20 @@ function calc_error_norms(func, u, t, analyzer,
     for i in eachnode(analyzer)
       u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t, equations)
       diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i), equations)
-      l2_error += diff.^2 * (weights[i] * jacobian_local[i])
+
+      l2_error  += diff.^2 * (weights[i] * jacobian_local[i])
       linf_error = @. max(linf_error, abs(diff))
+      l1_error  += abs.(diff) * (weights[i] * jacobian_local[i])
+
       total_volume += weights[i] * jacobian_local[i]
     end
   end
 
-  # For L2 error, divide by total volume
-  l2_error = @. sqrt(l2_error / total_volume)
+  # For L2/L1 error, divide by total volume
+  l2_error  = @. sqrt(l2_error / total_volume)
+  l1_error /= total_volume
 
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end
 
 
@@ -89,6 +94,7 @@ function calc_error_norms(func, u, t, analyzer,
   # Set up data structures
   l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
   linf_error = copy(l2_error)
+  l1_error   = copy(l2_error)
 
   # Iterate over all elements for error calculations
   for element in eachelement(dg, cache)
@@ -102,16 +108,19 @@ function calc_error_norms(func, u, t, analyzer,
     for i in eachnode(analyzer)
       u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t, equations)
       diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i), equations)
-      l2_error += diff.^2 * (weights[i] * volume_jacobian_)
+
+      l2_error  += diff.^2 * (weights[i] * volume_jacobian_)
       linf_error = @. max(linf_error, abs(diff))
+      l1_error  += abs.(diff) * (weights[i] * volume_jacobian_)
     end
   end
 
-  # For L2 error, divide by total volume
+  # For L2/L1 error, divide by total volume
   total_volume_ = total_volume(mesh)
-  l2_error = @. sqrt(l2_error / total_volume_)
+  l2_error  = @. sqrt(l2_error / total_volume_)
+  l1_error /= total_volume_
 
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end
 
 

src/callbacks_step/analysis_dg2d.jl

@@ -61,15 +61,17 @@ function calc_error_norms(func, u, t, analyzer,
   # Set up data structures
   l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations))
   linf_error = copy(l2_error)
+  l1_error   = copy(l2_error)
 
   # Iterate over all elements for error calculations
   # Accumulate L2 error on the element first so that the order of summation is the
   # same as in the parallel case to ensure exact equality. This facilitates easier parallel
   # development and debugging (see
   # https://github.com/trixi-framework/Trixi.jl/pull/850#pullrequestreview-757463943 for details).
   for element in eachelement(dg, cache)
-    # Set up data structures for local element L2 error
+    # Set up data structures for local element L2/L1 error
     l2_error_local = zero(l2_error)
+    l1_error_local = zero(l1_error)
 
     # Interpolate solution and node locations to analysis nodes
     multiply_dimensionwise!(u_local, vandermonde, view(u,                :, :, :, element), u_tmp1)
@@ -81,17 +83,21 @@ function calc_error_norms(func, u, t, analyzer,
     for j in eachnode(analyzer), i in eachnode(analyzer)
       u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations)
       diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations)
+
       l2_error_local += diff.^2 * (weights[i] * weights[j] * volume_jacobian_)
-      linf_error = @. max(linf_error, abs(diff))
+      linf_error      = @. max(linf_error, abs(diff))
+      l1_error_local += abs.(diff) * (weights[i] * weights[j] * volume_jacobian_)
     end
     l2_error += l2_error_local
+    l1_error += l1_error_local
   end
 
-  # For L2 error, divide by total volume
+  # For L2/L1 error, divide by total volume
   total_volume_ = total_volume(mesh)
-  l2_error = @. sqrt(l2_error / total_volume_)
+  l2_error  = @. sqrt(l2_error / total_volume_)
+  l1_error /= total_volume_
 
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end
 
 
@@ -105,6 +111,7 @@ function calc_error_norms(func, u, t, analyzer,
   # Set up data structures
   l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations))
   linf_error = copy(l2_error)
+  l1_error   = copy(l2_error)
   total_volume = zero(real(mesh))
 
   # Iterate over all elements for error calculations
@@ -120,16 +127,20 @@ function calc_error_norms(func, u, t, analyzer,
     for j in eachnode(analyzer), i in eachnode(analyzer)
       u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations)
       diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations)
+
       l2_error += diff.^2 * (weights[i] * weights[j] * jacobian_local[i, j])
       linf_error = @. max(linf_error, abs(diff))
+      l1_error += abs.(diff) * (weights[i] * weights[j] * jacobian_local[i, j])
+
       total_volume += weights[i] * weights[j] * jacobian_local[i, j]
     end
   end
 
-  # For L2 error, divide by total volume
-  l2_error = @. sqrt(l2_error / total_volume)
+  # For L2/L1 error, divide by total volume
+  l2_error  = @. sqrt(l2_error / total_volume)
+  l1_error /= total_volume
 
-  return l2_error, linf_error
+  return l2_error, linf_error, l1_error
 end