Fix StructuredMesh1D nonuniform meshes (#2620)

* Start fix `StructuredMesh1D` non-uniform

* Fix `StructuredMesh1D` nonuniform

* bb

* comms

* test & fmt

* fix test vals

Author: DanielDoehring

examples/structured_1d_dgsem/elixir_advection_nonuniform.jl

@@ -0,0 +1,41 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = 1.0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_godunov)
+
+cells_per_dimension = (24,)
+
+# This mapping converts [-1, 1] to [1, 9] with a non-uniform distribution of cells
+mapping(xi) = (xi + 2)^2
+
+mesh = StructuredMesh(cells_per_dimension, mapping)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+ode = semidiscretize(semi, (0.0, 0.5))
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0,
+            ode_default_options()..., callback = callbacks);

examples/structured_1d_dgsem/elixir_euler_convergence_nonuniform.jl

@@ -0,0 +1,49 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)
+solver = DGSEM(polydeg = 3, surface_flux = flux_hllc,
+               volume_integral = volume_integral)
+
+cells_per_dimension = (8,)
+
+# This mapping converts [-1, 1] to [0, 2] with a non-uniform distribution of cells
+mapping(xi) = ((xi + 2)^2 - 1) / 4
+
+mesh = StructuredMesh(cells_per_dimension, mapping)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/solvers/dgsem_structured/containers_2d.jl

@@ -102,7 +102,7 @@ function calc_jacobian_matrix!(jacobian_matrix, element,
     return jacobian_matrix
 end
 
-# Calculate contravarant vectors, multiplied by the Jacobian determinant J of the transformation mapping.
+# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping.
 # Those are called Ja^i in Kopriva's blue book.
 function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 5},
                                      element, jacobian_matrix)

src/solvers/dgsem_structured/dg_1d.jl

@@ -12,6 +12,7 @@ function calc_interface_flux!(cache, u, mesh::StructuredMesh{1},
 
     @threaded for element in eachelement(dg, cache)
         left_element = cache.elements.left_neighbors[1, element]
+        # => `element` is the right element of the interface
 
         if left_element > 0 # left_element = 0 at boundaries
             u_ll = get_node_vars(u, equations, dg, nnodes(dg), left_element)
@@ -67,4 +68,21 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple,
 
     return nothing
 end
+
+function apply_jacobian!(du, mesh::StructuredMesh{1},
+                         equations, dg::DG, cache)
+    @unpack inverse_jacobian = cache.elements
+
+    @threaded for element in eachelement(dg, cache)
+        for i in eachnode(dg)
+            factor = -inverse_jacobian[i, element]
+
+            for v in eachvariable(equations)
+                du[v, i, element] *= factor
+            end
+        end
+    end
+
+    return nothing
+end
 end # @muladd

src/solvers/dgsem_tree/dg_1d.jl

@@ -641,7 +641,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1
     return nothing
 end
 
-function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+function apply_jacobian!(du, mesh::TreeMesh{1},
                          equations, dg::DG, cache)
     @unpack inverse_jacobian = cache.elements
 

test/test_structured_1d.jl

@@ -33,6 +33,14 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_advection_nonuniform.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonuniform.jl"),
+                        l2=[0.0006665846145698006], linf=[0.00643334347367408])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_advection_shockcapturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_shockcapturing.jl"),
                         l2=[0.08015029105233593],
@@ -81,13 +89,31 @@ end
                         save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8), # Adding a small epsilon to avoid floating-point precision issues
                         callbacks=CallbackSet(summary_callback, save_solution,
                                               analysis_callback, alive_callback),
-                        l2=[5.726144824784944e-7],
-                        linf=[3.43073006914274e-6])
+                        l2=[5.725028892495733e-7],
+                        linf=[3.4292579200734252e-6])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     @test_allocations(Trixi.rhs!, semi, sol, 8000)
 end
 
+@trixi_testset "elixir_euler_convergence_nonuniform.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_convergence_nonuniform.jl"),
+                        l2=[
+                            6.0145954087568086e-5,
+                            4.865396929677369e-5,
+                            0.00012525608158029655
+                        ],
+                        linf=[
+                            0.00014050847320445925,
+                            0.00038307746603916115,
+                            0.0006872324329503243
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_euler_sedov.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
                         l2=[3.67478226e-01, 3.49491179e-01, 8.08910759e-01],