Skip to content

Commit

Permalink
Add AnalysisSurfaceIntegral (#1812)
Browse files Browse the repository at this point in the history 
* Add AnalysisSurfaceIntegral

* Minor change

* Update elixir_euler_subsonic_cylinder.jl

* Apply suggestions from code review

Co-authored-by: Daniel Doehring <[email protected]>

* Fix labels, add tests

* Attempt to fix CI

* Obtain indices from user chosen function

* Formatting

* Minor change

* Fix labels, add tests

* Attempt to fix CI

* Obtain indices from user chosen function

* Add new elixirs

* Add tests for NACA0012

* Fix tolerance, fix CI, add comments

* Run formatter

* Minor changes

* Run formatter

* Fix type instability and need for a vector

* Fix typo!

* Lower CFL to pass CI?

* Reduce amr_interval to pass CI?

* Maybe positivity limiter will fix CI?

* Remove AMR from CI

* That CI fail was on me!

* Test for write permit

* forces via names

* Boundary Names for other examples

* Shift to SSPRK54 to pass CI?

* docstrings

* alloc tests

* fmt

* Apply suggestions from code review

Co-authored-by: Daniel Doehring <[email protected]>
Co-authored-by: Joshua Lampert <[email protected]>

* Format

* non-allocating BC

* comments

* comments

* fmt

* Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl

* Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl

* comments

* make ready for review

* typos

* Apply suggestions from code review

Some key points

1) Change pre to p
2) Don't capitalize U
3) Don't append Swanson quantities with 'sw' 
4) linf -> l_inf (may be pending in some other places)
5) Other formatting improvements

Co-authored-by: Andrew Winters <[email protected]>

* Some fixes

* Format

* Update NEWS.md

* Update src/callbacks_step/analysis_surface_integral_2d.jl

* Polish

* Update examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl

* increase cfl

* fmt

* fix

---------

Co-authored-by: Daniel Doehring <[email protected]>
Co-authored-by: Daniel_Doehring <[email protected]>
Co-authored-by: Joshua Lampert <[email protected]>
Co-authored-by: Andrew Winters <[email protected]>
  • Loading branch information
@Arpit-Babbar @DanielDoehring
@JoshuaLampert @andrewwinters5000
5 people authored Mar 30, 2024
1 parent 73da92c commit 57f4b20
Show file tree
Hide file tree
Showing 11 changed files with 725 additions and 6 deletions.
3 changes: 2 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,10 @@ for human readability.
## Changes in the v0.7 lifecycle

#### Added
- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`.
- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
to 1D and 3D on `TreeMesh`.

- New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients.

## Changes when updating to v0.7 from v0.6.x

Expand Down
132 changes: 132 additions & 0 deletions examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

p_inf() = 1.0
rho_inf() = p_inf() / (1.0 * 287.87) # p_inf = 1.0, T = 1, R = 287.87
mach_inf() = 0.85
aoa() = pi / 180.0 # 1 Degree angle of attack
c_inf(equations) = sqrt(equations.gamma * p_inf() / rho_inf())
u_inf(equations) = mach_inf() * c_inf(equations)

@inline function initial_condition_mach085_flow(x, t,
equations::CompressibleEulerEquations2D)
v1 = u_inf(equations) * cos(aoa())
v2 = u_inf(equations) * sin(aoa())

prim = SVector(rho_inf(), v1, v2, p_inf())
return prim2cons(prim, equations)
end

initial_condition = initial_condition_mach085_flow

volume_flux = flux_ranocha_turbo
surface_flux = flux_lax_friedrichs

polydeg = 3
basis = LobattoLegendreBasis(polydeg)
shock_indicator = IndicatorHennemannGassner(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
volume_integral = volume_integral)

mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp",
joinpath(@__DIR__, "NACA0012.inp"))

mesh = P4estMesh{2}(mesh_file)

# The outer boundary is constant but subsonic, 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_mach085_flow(x, t, equations)

return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations)
end

boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant,
:Right => boundary_condition_subsonic_constant,
:Top => boundary_condition_subsonic_constant,
:Bottom => boundary_condition_subsonic_constant,
:AirfoilBottom => boundary_condition_slip_wall,
:AirfoilTop => boundary_condition_slip_wall)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

###############################################################################
# ODE solvers

# Run for a long time to reach a steady state
tspan = (0.0, 20.0)
ode = semidiscretize(semi, tspan)

# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 2000

l_inf = 1.0 # Length of airfoil

force_boundary_names = [:AirfoilBottom, :AirfoilTop]
drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
DragCoefficientPressure(aoa(), rho_inf(),
u_inf(equations), l_inf))

lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
LiftCoefficientPressure(aoa(), rho_inf(),
u_inf(equations), l_inf))

analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
output_directory = "out",
save_analysis = true,
analysis_integrals = (drag_coefficient,
lift_coefficient))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 500,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 1.0)

amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
base_level = 1,
med_level = 3, med_threshold = 0.05,
max_level = 4, max_threshold = 0.1)

amr_interval = 100
amr_callback = AMRCallback(semi, amr_controller,
interval = amr_interval,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
stepsize_callback, amr_callback)

###############################################################################
# run the simulation
sol = solve(ode, SSPRK54(thread = OrdinaryDiffEq.True()),
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
125 changes: 125 additions & 0 deletions examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

@inline function initial_condition_mach038_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_mach038_flow

volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used
surface_flux = flux_lax_friedrichs

polydeg = 3
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 = polydeg,
periodicity = (false, true))

# The boundary condition on the outer cylinder is constant but subsonic, so we cannot compute the
# boundary flux from 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_mach038_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 a steady 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
l_inf = 1.0 # Diameter of circle

drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
DragCoefficientPressure(aoa, rho_inf, u_inf,
l_inf))

lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
LiftCoefficientPressure(aoa, rho_inf, u_inf,
l_inf))

analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
output_directory = "out",
save_analysis = true,
analysis_integrals = (drag_coefficient,
lift_coefficient))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 500,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 2.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
stepsize_callback)

###############################################################################
# run the simulation
sol = solve(ode,
CarpenterKennedy2N54(williamson_condition = false;
thread = OrdinaryDiffEq.True()),
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
Loading

0 comments on commit 57f4b20

Please sign in to comment.