Add FVO2 for curved 3D meshes

examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded_fvO2.jl

@@ -0,0 +1,95 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+initial_condition = initial_condition_constant
+
+boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))
+
+# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
+# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
+# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
+# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
+# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
+# `StepsizeCallback` (CFL-Condition) and less diffusion.
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
+
+volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                      volume_flux_fv = surface_flux,
+                                                      reconstruction_mode = reconstruction_O2_full,
+                                                      slope_limiter = monotonized_central)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D.
+# This particular mesh is unstructured in the yz-plane, but extruded in x-direction.
+# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded
+# in x-direction to ensure free stream preservation on a non-conforming mesh.
+# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
+function mapping(xi, eta_, zeta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    eta = 1.5 * eta_ + 1.5
+    zeta = 1.5 * zeta_ + 1.5
+
+    z = zeta +
+        1 / 6 * (cos(1.5 * pi * (2 * eta - 3) / 3) *
+                 cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    y = eta + 1 / 6 * (cos(0.5 * pi * (2 * eta - 3) / 3) *
+                       cos(2 * pi * (2 * z - 3) / 3))
+
+    return SVector(xi, y, z)
+end
+
+# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3
+mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
+                           joinpath(@__DIR__, "cube_unstructured_2.inp"))
+
+mesh = P4estMesh{3}(mesh_file, polydeg = polydeg,
+                    mapping = mapping)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.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)
+
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        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);

examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl

@@ -0,0 +1,88 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:x_neg => boundary_condition,
+                           :x_pos => boundary_condition,
+                           :y_neg => boundary_condition,
+                           :y_pos => boundary_condition,
+                           :z_neg => boundary_condition,
+                           :z_pos => boundary_condition)
+
+polydeg = 2
+basis = LobattoLegendreBasis(polydeg)
+surface_flux = flux_hll
+
+volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                      volume_flux_fv = surface_flux,
+                                                      reconstruction_mode = reconstruction_O2_full,
+                                                      slope_limiter = monotonized_central)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+function mapping(xi, eta, zeta)
+    # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
+    # xi = 1.5 * xi_ + 1.5
+    # eta = 1.5 * eta_ + 1.5
+    # zeta = 1.5 * zeta_ + 1.5
+
+    y = eta +
+        1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+         cos(0.5 * pi * (2 * eta - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    x = xi +
+        1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+         cos(2 * pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    z = zeta +
+        1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) *
+         cos(pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (2, 2, 2)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg, mapping = mapping,
+                 initial_refinement_level = 1,
+                 periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions,
+                                    source_terms = source_terms)
+###############################################################################
+# 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 = 1.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+###############################################################################
+# run the simulation
+
+sol = solve(ode, ParsaniKetchesonDeconinck3S82();
+            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/calc_volume_integral.jl

@@ -146,7 +146,8 @@ function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{1}, StructuredMesh{1},
                                            TreeMesh{2}, StructuredMesh{2}, P4estMesh{2},
                                            UnstructuredMesh2D, T8codeMesh{2},
-                                           TreeMesh{3}},
+                                           TreeMesh{3}, StructuredMesh{3},
+                                           P4estMesh{3}, T8codeMesh{3}},
                                have_nonconservative_terms, equations,
                                volume_integral::VolumeIntegralPureLGLFiniteVolumeO2,
                                dg::DGSEM, cache)

src/solvers/dgsem_structured/dg_3d.jl

@@ -528,6 +528,92 @@ end
     return nothing
 end
 
+@inline function calcflux_fvO2!(fstar1_L, fstar1_R, fstar2_L, fstar2_R,
+                                fstar3_L, fstar3_R, u,
+                                mesh::Union{StructuredMesh{3}, P4estMesh{3},
+                                            T8codeMesh{3}},
+                                have_nonconservative_terms::False,
+                                equations,
+                                volume_flux_fv, dg::DGSEM, element, cache,
+                                sc_interface_coords, reconstruction_mode, slope_limiter)
+    @unpack normal_vectors_1, normal_vectors_2, normal_vectors_3 = cache.normal_vectors
+
+    for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
+        u_ll = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), j, k, element),
+                         equations)
+        u_lr = cons2prim(get_node_vars(u, equations, dg, i - 1, j, k, element),
+                         equations)
+        u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element),
+                         equations)
+
+        u_rr = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1), j, k,
+                                       element), equations)
+
+        u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr,
+                                       sc_interface_coords, i,
+                                       slope_limiter, dg)
+
+        normal_direction = get_normal_vector(normal_vectors_1, i - 1, j, k, element)
+
+        contravariant_flux = volume_flux_fv(prim2cons(u_l, equations),
+                                            prim2cons(u_r, equations),
+                                            normal_direction, equations)
+
+        set_node_vars!(fstar1_L, contravariant_flux, equations, dg, i, j, k)
+        set_node_vars!(fstar1_R, contravariant_flux, equations, dg, i, j, k)
+    end
+
+    for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
+        u_ll = cons2prim(get_node_vars(u, equations, dg, i, max(1, j - 2), k, element),
+                         equations)
+        u_lr = cons2prim(get_node_vars(u, equations, dg, i, j - 1, k, element),
+                         equations)
+        u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element),
+                         equations)
+        u_rr = cons2prim(get_node_vars(u, equations, dg, i, min(nnodes(dg), j + 1), k,
+                                       element), equations)
+
+        u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr,
+                                       sc_interface_coords, j,
+                                       slope_limiter, dg)
+
+        normal_direction = get_normal_vector(normal_vectors_2, i, j - 1, k, element)
+
+        contravariant_flux = volume_flux_fv(prim2cons(u_l, equations),
+                                            prim2cons(u_r, equations),
+                                            normal_direction, equations)
+
+        set_node_vars!(fstar2_L, contravariant_flux, equations, dg, i, j, k)
+        set_node_vars!(fstar2_R, contravariant_flux, equations, dg, i, j, k)
+    end
+
+    for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
+        u_ll = cons2prim(get_node_vars(u, equations, dg, i, j, max(1, k - 2), element),
+                         equations)
+        u_lr = cons2prim(get_node_vars(u, equations, dg, i, j, k - 1, element),
+                         equations)
+        u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element),
+                         equations)
+        u_rr = cons2prim(get_node_vars(u, equations, dg, i, j, min(nnodes(dg), k + 1),
+                                       element), equations)
+
+        u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr,
+                                       sc_interface_coords, k,
+                                       slope_limiter, dg)
+
+        normal_direction = get_normal_vector(normal_vectors_3, i, j, k - 1, element)
+
+        contravariant_flux = volume_flux_fv(prim2cons(u_l, equations),
+                                            prim2cons(u_r, equations),
+                                            normal_direction, equations)
+
+        set_node_vars!(fstar3_L, contravariant_flux, equations, dg, i, j, k)
+        set_node_vars!(fstar3_R, contravariant_flux, equations, dg, i, j, k)
+    end
+
+    return nothing
+end
+
 function calc_interface_flux!(cache, u, mesh::StructuredMesh{3},
                               have_nonconservative_terms, # can be True/False
                               equations, surface_integral, dg::DG)