Add AnalysisSurfaceIntegral

examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl

@@ -0,0 +1,122 @@
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+# include("analysis_surface_2d.jl")
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+@inline function initial_condition_mach01_flow(x, t,
+                                               equations::CompressibleEulerEquations2D)
+    # set the freestream flow parameters
+    rho_freestream = 1.4
+    v1 = 0.38
+    v2 = 0.0
+    p_freestream = 1.0
+
+    prim = SVector(rho_freestream, v1, v2, p_freestream)
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_mach01_flow
+
+volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used
+surface_flux = flux_lax_friedrichs
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+function mapping2cylinder(xi, eta)
+    xi_, eta_ = 0.5 * (xi + 1), 0.5 * (eta + 1.0) # Map from [-1,1] to [0,1] for simplicity
+
+    R2 = 50.0 # Bigger circle
+    R1 = 0.5  # Smaller circle
+
+    # Ensure an isotropic mesh by using elements with smaller radial length near the inner circle
+
+    r = R1 * exp(xi_ * log(R2 / R1))
+    theta = 2.0 * pi * eta_
+
+    x = r * cos(theta)
+    y = r * sin(theta)
+    return (x, y)
+end
+
+cells_per_dimension = (64, 64)
+# xi = -1 maps to the inner circle and xi = +1 maps to the outer circle and we can specify boundary
+# conditions there. However, the image of eta = -1, +1 coincides at the line y = 0. There is no
+# physical boundary there so we specify periodicity = true there and the solver treats the
+# element across eta = -1, +1 as neighbours which is what we want
+mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = 3,
+                 periodicity = (false, true))
+
+# The boundary of the outer cylinder is constant but supersonic, so we cannot compute the
+# boundary flux for the external information alone. Thus, we use the numerical flux to distinguish
+# between inflow and outflow characteristics
+@inline function boundary_condition_subsonic_constant(u_inner,
+                                                      normal_direction::AbstractVector, x,
+                                                      t,
+                                                      surface_flux_function,
+                                                      equations::CompressibleEulerEquations2D)
+    u_boundary = initial_condition_mach01_flow(x, t, equations)
+
+    return surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall,
+                           :x_pos => boundary_condition_subsonic_constant)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers
+
+# Run for a long time to reach state state
+tspan = (0.0, 100.0)
+ode = semidiscretize(semi, tspan)
+
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 2000
+
+aoa = 0.0
+rho_inf = 1.4
+U_inf = 0.38
+linf = 1.0 # Diameter of circle
+
+indices = semi_ -> semi.boundary_conditions.boundary_indices[2]
+my_drag_force = Trixi.AnalysisSurfaceIntegral(indices,
+                                              Trixi.DragForcePressure(aoa, rho_inf, U_inf,
+                                                                      linf))
+
+my_lift_force = Trixi.AnalysisSurfaceIntegral(indices,
+                                              Trixi.LiftForcePressure(aoa, rho_inf, U_inf,
+                                                                      linf))
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_errors = Symbol[],
+                                     output_directory = "analysis_results",
+                                     save_analysis = true,
+                                     analysis_integrals = (my_drag_force, my_lift_force))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, SSPRK43();
+            ode_default_options()..., callback = callbacks);
+summary_callback() # print the timer summary

src/callbacks_step/analysis.jl

@@ -691,6 +691,7 @@ end # @muladd
 # specialized implementations specific to some solvers
 include("analysis_dg1d.jl")
 include("analysis_dg2d.jl")
+include("analysis_surface_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")

src/callbacks_step/analysis_surface_2d.jl

@@ -0,0 +1,116 @@
+using Trixi
+using Trixi: integrate_via_indices, norm, apply_jacobian_parabolic!, @threaded,
+             indices2direction,
+             index_to_start_step_2d, get_normal_direction, dot, get_node_coords
+import Trixi: analyze, pretty_form_ascii, pretty_form_utf
+
+struct AnalysisSurfaceIntegral{Indices, Variable}
+    indices::Indices
+    variable::Variable
+end
+
+struct ForceState{RealT <: Real}
+    Ψl::Tuple{RealT, RealT}
+    rhoinf::RealT
+    uinf::RealT
+    linf::RealT
+end
+
+# TODO - This should be a struct in ForceState
+struct FreeStreamVariables{RealT <: Real}
+    rhoinf::RealT
+    uinf::RealT
+    linf::RealT
+end
+
+struct LiftForcePressure{RealT <: Real}
+    force_state::ForceState{RealT}
+end
+
+struct DragForcePressure{RealT <: Real}
+    force_state::ForceState{RealT}
+end
+
+function LiftForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real)
+    Ψl = (-sin(aoa), cos(aoa))
+    force_state = ForceState(Ψl, rhoinf, uinf, linf)
+    return LiftForcePressure(force_state)
+end
+
+function DragForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real)
+    Ψd = (cos(aoa), sin(aoa))
+    return DragForcePressure(ForceState(Ψd, rhoinf, uinf, linf))
+end
+
+function (lift_force::LiftForcePressure)(u, normal_direction, equations)
+    p = pressure(u, equations)
+    @unpack Ψl, rhoinf, uinf, linf = lift_force.force_state
+    n = dot(normal_direction, Ψl) / norm(normal_direction)
+    return p * n / (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function (drag_force::DragForcePressure)(u, normal_direction, equations)
+    p = pressure(u, equations)
+    @unpack Ψl, rhoinf, uinf, linf = drag_force.force_state
+    n = dot(normal_direction, Ψl) / norm(normal_direction)
+    return p * n / (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
+                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}},
+                 equations, dg::DGSEM, cache)
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+    @unpack indices, variable = surface_variable
+    # TODO - Use initialize callbacks to move boundary_conditions to cache
+    indices_ = indices(cache)
+
+    surface_integral = zero(eltype(u))
+    index_range = eachnode(dg)
+    for local_index in eachindex(indices_)
+        # Use the local index to get the global boundary index from the pre-sorted list
+        boundary = indices_[local_index]
+
+        # Get information on the adjacent element, compute the surface fluxes,
+        # and store them
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in eachnode(dg)
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node,
+                                                    element)
+
+            # L2 norm of normal direction is the surface element
+            # 0.5 factor is NOT needed, the norm(normal_direction) is all the factor needed
+            dS = weights[node_index] * norm(normal_direction)
+            surface_integral += variable(u_node, normal_direction, equations) * dS
+
+            i_node += i_node_step
+            j_node += j_node_step
+        end
+    end
+    return surface_integral
+end
+
+function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:LiftForcePressure{<:Any}})
+    "Pressure_lift"
+end
+function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftForcePressure{<:Any}})
+    "Pressure_lift"
+end
+function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}})
+    "Pressure_drag"
+end
+function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}})
+    "Pressure_drag"
+end