From 247a5f18ef11be9a82c1c289bca91308120acdf7 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 16:26:53 +0530 Subject: [PATCH] Add AnalysisSurfaceIntegral --- .../elixir_euler_subsonic_cylinder.jl | 122 ++++++++++++++++++ src/callbacks_step/analysis.jl | 1 + src/callbacks_step/analysis_surface_2d.jl | 116 +++++++++++++++++ 3 files changed, 239 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl create mode 100644 src/callbacks_step/analysis_surface_2d.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl new file mode 100644 index 00000000000..e6b438944f3 --- /dev/null +++ b/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 diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index ba232032951..2c3351272fa 100644 --- a/src/callbacks_step/analysis.jl +++ b/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") diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl new file mode 100644 index 00000000000..f984849448b --- /dev/null +++ b/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