AMR-Compat SurfaceIntegrals

NEWS.md

@@ -4,6 +4,19 @@ 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.8 from v0.7.x
+
+#### Added
+
+#### Changed
+
+- The specification of boundary names on which `AnalysisSurfaceIntegral`s are computed (such as drag and lift coefficients) has changed from `Symbol` and `Vector{Symbol}` to `NTuple{Symbol}`.
+Thus, for one boundary the syntax changes from `:boundary` to `(:boundary,)` and for `Vector`s `[:boundary1, :boundary2]` to `(:boundary1, boundary2)`.
+
+#### Deprecated
+
+#### Removed
+
 ## Changes in the v0.7 lifecycle
 
 #### Added

examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl

@@ -84,7 +84,7 @@ analysis_interval = 2000
 
 l_inf = 1.0 # Length of airfoil
 
-force_boundary_names = [:AirfoilBottom, :AirfoilTop]
+force_boundary_names = (:AirfoilBottom, :AirfoilTop)
 drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))

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(semi, (:x_neg,),
                                            DragCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
+lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
                                            LiftCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 

examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl

@@ -119,7 +119,7 @@ summary_callback = SummaryCallback()
 
 analysis_interval = 2000
 
-force_boundary_names = [:AirfoilBottom, :AirfoilTop]
+force_boundary_names = (:AirfoilBottom, :AirfoilTop)
 drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),

src/callbacks_step/analysis.jl

@@ -657,6 +657,15 @@ include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
 
+# This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires
+# `semi` to be passed along to retrieve the current boundary indices, which are non-static 
+# in the case of AMR.
+function analyze(quantity::AnalysisSurfaceIntegral, du, u, t,
+                 semi::AbstractSemidiscretization)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, semi)
+end
+
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
 # Note that this needs to be included after `analysis_surface_integral_2d.jl` to
@@ -669,6 +678,6 @@ function analyze(quantity::AnalysisSurfaceIntegral{Variable},
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     equations_parabolic = semi.equations_parabolic
     cache_parabolic = semi.cache_parabolic
-    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
             cache_parabolic)
 end

src/callbacks_step/analysis_surface_integral_2d.jl

@@ -21,30 +21,18 @@ drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the b
 name `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
-- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries
+- `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
 """
-struct AnalysisSurfaceIntegral{Variable}
-    indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
+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_symbol, variable)
-        @unpack boundary_symbol_indices = semi.boundary_conditions
-        indices = boundary_symbol_indices[boundary_symbol]
-
-        return new{typeof(variable)}(indices, variable)
-    end
-
-    function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable)
-        @unpack boundary_symbol_indices = semi.boundary_conditions
-        indices = Vector{Int}()
-        for name in boundary_symbols
-            append!(indices, boundary_symbol_indices[name])
-        end
-        sort!(indices)
-
-        return new{typeof(variable)}(indices, variable)
+    function AnalysisSurfaceIntegral(semi,
+                                     boundary_symbols::NTuple{NBoundaries, Symbol},
+                                     variable) where {NBoundaries}
+        return new{typeof(variable), NBoundaries}(variable, boundary_symbols)
     end
 end
 
@@ -255,14 +243,26 @@ function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x,
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
+function get_boundary_indices(boundary_symbols, boundary_symbol_indices)
+    indices = Vector{Int}()
+    for name in boundary_symbols
+        append!(indices, boundary_symbol_indices[name])
+    end
+    sort!(indices) # Try to achieve some data locality by sorting
+
+    return indices
+end
+
 function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
                  mesh::P4estMesh{2},
-                 equations, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache, semi)
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
 
-    @unpack indices, variable = surface_variable
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
 
     surface_integral = zero(eltype(u))
     index_range = eachnode(dg)
@@ -308,13 +308,15 @@ end
 function analyze(surface_variable::AnalysisSurfaceIntegral{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
-                 dg::DGSEM, cache,
+                 dg::DGSEM, cache, semi,
                  cache_parabolic) where {Variable <: VariableViscous}
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
 
-    @unpack indices, variable = surface_variable
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
 
     # Additions for parabolic
     @unpack viscous_container = cache_parabolic

test/test_p4est_2d.jl

@@ -604,9 +604,9 @@ end
     u = Trixi.wrap_array(u_ode, semi)
     du = Trixi.wrap_array(du_ode, semi)
     drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver,
-                         semi.cache)
+                         semi.cache, semi)
     lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver,
-                         semi.cache)
+                         semi.cache, semi)
 
     @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13)
     @test isapprox(drag, 2.588589856781827, atol = 1e-13)
@@ -617,13 +617,14 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_NACA0012airfoil_mach085.jl"),
                         l2=[
-                            5.371568111383228e-7, 6.4158131303956445e-6,
-                            1.0324346542348325e-5, 0.0006348064933187732,
+                            5.634402680811982e-7, 6.748066107517321e-6,
+                            1.091879472416885e-5, 0.0006686372064029146,
                         ],
                         linf=[
-                            0.0016263400091978443, 0.028471072159724428,
-                            0.02986133204785877, 1.9481060511014872,
+                            0.0021456247890772823, 0.03957142889488085,
+                            0.03832024233032798, 2.6628739573358495,
                         ],
+                        amr_interval=1,
                         base_level=0, med_level=1, max_level=1,
                         tspan=(0.0, 0.0001),
                         adapt_initial_condition=false,
@@ -647,12 +648,12 @@ end
     u = Trixi.wrap_array(u_ode, semi)
     du = Trixi.wrap_array(du_ode, semi)
     drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver,
-                         semi.cache)
+                         semi.cache, semi)
     lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver,
-                         semi.cache)
+                         semi.cache, semi)
 
-    @test isapprox(lift, 0.0262382560809345, atol = 1e-13)
-    @test isapprox(drag, 0.10898248971932244, atol = 1e-13)
+    @test isapprox(lift, 0.029076443678087403, atol = 1e-13)
+    @test isapprox(drag, 0.13564720009197903, atol = 1e-13)
 end
 end
 

test/test_parabolic_2d.jl

@@ -739,16 +739,16 @@ end
     du = Trixi.wrap_array(du_ode, semi)
 
     drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver,
-                           semi.cache)
+                           semi.cache, semi)
     lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver,
-                           semi.cache)
+                           semi.cache, semi)
 
     drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh,
                            equations, equations_parabolic, solver,
-                           semi.cache, semi.cache_parabolic)
+                           semi.cache, semi, semi.cache_parabolic)
     lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh,
                            equations, equations_parabolic, solver,
-                           semi.cache, semi.cache_parabolic)
+                           semi.cache, semi, semi.cache_parabolic)
 
     @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13)
     @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13)