Merge branch 'main' into velocity

trixi-framework · Oct 10, 2024 · cb06499 · cb06499
2 parents 82976ac + 0c2c367
commit cb06499
Show file tree

Hide file tree

Showing 46 changed files with 338 additions and 380 deletions.
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.24.3
+        uses: crate-ci/typos@v1.25.0
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,30 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.9 from v0.8.x
+
+#### Added
+
+- Boundary conditions are now supported on nonconservative terms ([#2062]).
+
+#### Changed
+
+- We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
+  `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`.
+  ([#2069])
+- In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
+  Users should not be affected by this.
+- Nonconservative terms depend only on `normal_direction_average` instead of both 
+  `normal_direction_average` and `normal_direction_ll`, such that the function signature is now 
+  identical with conservative fluxes. This required a change of the `normal_direction` in
+  `flux_nonconservative_powell` ([#2062]).
+
+#### Deprecated
+
+#### Removed
+
+
 ## Changes in the v0.8 lifecycle
 
 #### Changed

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.8.11-DEV"
+version = "0.9.1-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
@@ -35,7 +35,7 @@ plot!(Shape([(-1.7,-11), (1.7,-11), (1.7,-17), (-1.7,-17)]), linecolor=:blue, fi
 annotate!(1.5, -14, ("Trixi.jl", 12, :blue, :center))
 
 plot!(Shape([(-1.4,-14.5), (1.4,-14.5), (1.4,-16.5), (-1.4,-16.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
-annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, initial_condition, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
+annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
 
 
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
@@ -85,11 +85,11 @@ analysis_interval = 2000
 l_inf = 1.0 # Length of airfoil
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
@@ -89,11 +89,11 @@ rho_inf = 1.4
 u_inf = 0.38
 l_inf = 1.0 # Diameter of circle
 
-drag_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+drag_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            DragCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+lift_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            LiftCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -120,23 +120,23 @@ summary_callback = SummaryCallback()
 analysis_interval = 2000
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        DragCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        LiftCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -31,14 +31,6 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-# Extend the definition of the non-conservative Powell flux functions.
-import Trixi.flux_nonconservative_powell
-function flux_nonconservative_powell(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     equations::IdealGlmMhdEquations2D)
-    flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
-                                equations)
-end
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
 solver = DGSEM(polydeg = 3,
                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),

diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl
@@ -43,10 +43,7 @@ end
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-function inflow(x, t, equations::InviscidBurgersEquation1D)
-    return initial_condition_rarefaction(coordinate_min, t, equations)
-end
-boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,

diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
@@ -43,10 +43,7 @@ end
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-function inflow(x, t, equations::InviscidBurgersEquation1D)
-    return initial_condition_shock(coordinate_min, t, equations)
-end
-boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,

diff --git a/src/basic_types.jl b/src/basic_types.jl
@@ -71,14 +71,14 @@ struct BoundaryConditionDoNothing end
                                                 orientation_or_normal_direction,
                                                 direction::Integer, x, t, surface_flux,
                                                 equations)
-    return flux(u_inner, orientation_or_normal_direction, equations)
+    return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
 end
 
 # This version can be called by hyperbolic solvers on unstructured, curved meshes
 @inline function (::BoundaryConditionDoNothing)(u_inner,
                                                 outward_direction::AbstractVector,
                                                 x, t, surface_flux, equations)
-    return flux(u_inner, outward_direction, equations)
+    return surface_flux(u_inner, u_inner, outward_direction, equations)
 end
 
 # This version can be called by parabolic solvers

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
@@ -314,8 +314,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
         mpi_println(" #elements:      " *
                     @sprintf("% 14d", nelementsglobal(mesh, solver, cache)))
 
-        # Level information (only show for AMR)
-        print_amr_information(integrator.opts.callback, mesh, solver, cache)
+        # Level information (only for AMR and/or non-uniform `TreeMesh`es)
+        print_level_information(integrator.opts.callback, mesh, solver, cache)
         mpi_println()
 
         # Open file for appending and store time step and time information
@@ -489,21 +489,7 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
     return nothing
 end
 
-# Print level information only if AMR is enabled
-function print_amr_information(callbacks, mesh, solver, cache)
-
-    # Return early if there is nothing to print
-    uses_amr(callbacks) || return nothing
-
-    # Get global minimum and maximum level from the AMRController
-    min_level = max_level = 0
-    for cb in callbacks.discrete_callbacks
-        if cb.affect! isa AMRCallback
-            min_level = cb.affect!.controller.base_level
-            max_level = cb.affect!.controller.max_level
-        end
-    end
-
+function print_level_information(mesh, solver, cache, min_level, max_level)
     # Get local element count per level
     elements_per_level = get_elements_per_level(min_level, max_level, mesh, solver,
                                                 cache)
@@ -522,6 +508,45 @@ function print_amr_information(callbacks, mesh, solver, cache)
     return nothing
 end
 
+function print_level_information(callbacks, mesh::TreeMesh, solver, cache)
+    if uses_amr(callbacks)
+        # Get global minimum and maximum level from the AMRController
+        min_level = max_level = 0
+        for cb in callbacks.discrete_callbacks
+            if cb.affect! isa AMRCallback
+                min_level = cb.affect!.controller.base_level
+                max_level = cb.affect!.controller.max_level
+            end
+        end
+        print_level_information(mesh, solver, cache, min_level, max_level)
+        # Special check for `TreeMesh`es without AMR.
+        # These meshes may still be non-uniform due to `refinement_patches`, see 
+        # `refine_box!`, `coarsen_box!`, and `refine_sphere!`.
+    elseif minimum_level(mesh.tree) != maximum_level(mesh.tree)
+        min_level = minimum_level(mesh.tree)
+        max_level = maximum_level(mesh.tree)
+        print_level_information(mesh, solver, cache, min_level, max_level)
+    else # Uniform mesh
+        return nothing
+    end
+end
+
+function print_level_information(callbacks, mesh, solver, cache)
+    if uses_amr(callbacks)
+        # Get global minimum and maximum level from the AMRController
+        min_level = max_level = 0
+        for cb in callbacks.discrete_callbacks
+            if cb.affect! isa AMRCallback
+                min_level = cb.affect!.controller.base_level
+                max_level = cb.affect!.controller.max_level
+            end
+        end
+        print_level_information(mesh, solver, cache, min_level, max_level)
+    else # Uniform mesh
+        return nothing
+    end
+end
+
 function get_elements_per_level(min_level, max_level, mesh::P4estMesh, solver, cache)
     elements_per_level = zeros(P4EST_MAXLEVEL + 1)
 

diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -20,7 +20,6 @@ For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@
 drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
 name `:Airfoil` in 2D.
 
-- `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries
   where the quantity of interest is computed
 - `variable::Variable`: Quantity of interest, like lift or drag
@@ -29,8 +28,7 @@ struct AnalysisSurfaceIntegral{Variable, NBoundaries}
     variable::Variable # Quantity of interest, like lift or drag
     boundary_symbols::NTuple{NBoundaries, Symbol} # Name(s) of the boundary/boundaries
 
-    function AnalysisSurfaceIntegral(semi,
-                                     boundary_symbols::NTuple{NBoundaries, Symbol},
+    function AnalysisSurfaceIntegral(boundary_symbols::NTuple{NBoundaries, Symbol},
                                      variable) where {NBoundaries}
         return new{typeof(variable), NBoundaries}(variable, boundary_symbols)
     end

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
@@ -196,18 +196,15 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations2D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations2D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations2D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -254,8 +251,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations2D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -265,14 +261,9 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
@@ -300,8 +291,9 @@ of the [`IdealGlmMhdEquations2D`](@ref).
 
 This implementation uses a non-conservative term that can be written as the product
 of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm
-et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different
-results on non-conforming meshes(!).
+et al. [`flux_nonconservative_powell`](@ref) for conforming meshes but it yields different
+results on non-conforming meshes(!). On curvilinear meshes this formulation applies the
+local normal direction compared to the averaged one used in [`flux_nonconservative_powell`](@ref).
 
 The two other flux functions with the same name return either the local
 or symmetric portion of the non-conservative flux based on the type of the

diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
@@ -224,18 +224,15 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations3D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations3D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations3D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -292,8 +289,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -303,16 +299,11 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] +
-                 v3_ll * normal_direction_ll[3]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2] +
-                 B3_rr * normal_direction_average[3]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
+                 v3_ll * normal_direction[3]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2] +
+                 B3_rr * normal_direction[3]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})