|
| 1 | +using Trixi |
| 2 | +using OrdinaryDiffEqSSPRK |
| 3 | + |
| 4 | +############################################################################### |
| 5 | +# Geometry & boundary conditions |
| 6 | + |
| 7 | +# Mapping to create a "close-up" mesh around the second quadrant of a cylinder, |
| 8 | +# implemented by Georgii Oblapenko. If you use this in your own work, please cite: |
| 9 | +# |
| 10 | +# - G. Oblapenko and A. Tarnovskiy (2024) |
| 11 | +# Reproducibility Repository for the paper: |
| 12 | +# Entropy-stable fluxes for high-order Discontinuous Galerkin simulations of high-enthalpy flows. |
| 13 | +# [DOI: 10.5281/zenodo.13981615](https://doi.org/10.5281/zenodo.13981615) |
| 14 | +# [GitHub](https://github.com/knstmrd/paper_ec_trixi_chem) |
| 15 | +# |
| 16 | +# as well as the corresponding paper: |
| 17 | +# - G. Oblapenko and M. Torrilhon (2025) |
| 18 | +# Entropy-conservative high-order methods for high-enthalpy gas flows. |
| 19 | +# Computers & Fluids, 2025. |
| 20 | +# [DOI: 10.1016/j.compfluid.2025.106640](https://doi.org/10.1016/j.compfluid.2025.106640) |
| 21 | +# |
| 22 | +# The mapping produces the following geometry & shock (indicated by the asterisks `* `): |
| 23 | +# ____x_neg____ |
| 24 | +# | | |
| 25 | +# | | |
| 26 | +# | | |
| 27 | +# | * | |
| 28 | +# | * y |
| 29 | +# | Inflow * _ |
| 30 | +# | state * p |
| 31 | +# x * o |
| 32 | +# _ * s |
| 33 | +# n * | |
| 34 | +# e * | |
| 35 | +# g Shock . |
| 36 | +# | * . |
| 37 | +# | * . <- x_pos |
| 38 | +# | * . |
| 39 | +# | * . (Cylinder) |
| 40 | +# |_______y_neg_______. |
| 41 | +function mapping_cylinder_shock_fitted(xi_, eta_, |
| 42 | + cylinder_radius, spline_points) |
| 43 | + shock_shape = [ |
| 44 | + (spline_points[1], 0.0), # Shock position on the stagnation line (`y_neg`, y = 0) |
| 45 | + (spline_points[2], spline_points[2]), # Shock position at -45° angle |
| 46 | + (0.0, spline_points[3]) # Shock position at outflow (`y_pos`, x = x_max) |
| 47 | + ] # 3 points that define the geometry of the mesh which follows the shape of the shock (known a-priori) |
| 48 | + R = [sqrt(shock_shape[i][1]^2 + shock_shape[i][2]^2) for i in 1:3] # 3 radii |
| 49 | + |
| 50 | + # Construct spline with form R[1] + c2 * eta_01^2 + c3 * eta_01^3, |
| 51 | + # chosen such that derivative w.r.t eta_01 is 0 at eta_01 = 0 such that |
| 52 | + # we have symmetry along the stagnation line (`y_neg`, y = 0). |
| 53 | + # |
| 54 | + # A single cubic spline doesn't fit the shock perfectly, |
| 55 | + # but is the simplest curve that does a reasonable job and it also can be easily computed analytically. |
| 56 | + # The choice of points on the stagnation line and outflow region is somewhat self-evident |
| 57 | + # (capture the minimum and maximum extent of the shock stand-off), |
| 58 | + # and the point at the 45 degree angle seemed the most logical to add |
| 59 | + # since it only requires one additional value (and not two), |
| 60 | + # simplifies the math a bit, and the angle lies exactly in between the other angles. |
| 61 | + spline_matrix = [1.0 1.0; 0.25 0.125] |
| 62 | + spline_RHS = [R[3] - R[1], R[2] - R[1]] |
| 63 | + spline_coeffs = spline_matrix \ spline_RHS # c2, c3 |
| 64 | + |
| 65 | + eta_01 = (eta_ + 1) / 2 # Transform `eta_` in [-1, 1] to `eta_01` in [0, 1] |
| 66 | + # "Flip" `xi_` in [-1, 1] to `xi_01` in [0, 1] since |
| 67 | + # shock positions where originally for first quadrant, here we use second quadrant |
| 68 | + xi_01 = (-xi_ + 1) / 2 |
| 69 | + |
| 70 | + R_outer = R[1] + spline_coeffs[1] * eta_01^2 + spline_coeffs[2] * eta_01^3 |
| 71 | + |
| 72 | + angle = -π / 4 + eta_ * π / 4 # Angle runs from -90° to 0° |
| 73 | + r = (cylinder_radius + xi_01 * (R_outer - cylinder_radius)) |
| 74 | + |
| 75 | + return SVector(round(r * sin(angle); digits = 8), round(r * cos(angle); digits = 8)) |
| 76 | +end |
| 77 | + |
| 78 | +@inline function initial_condition_mach6_flow(x, t, |
| 79 | + equations::NonIdealCompressibleEulerEquations2D) |
| 80 | + RealT = eltype(x) |
| 81 | + eos = equations.equation_of_state |
| 82 | + |
| 83 | + # Freestream conditions at 40 km altitude |
| 84 | + rho = 0.00385101 # [kg/m^3] |
| 85 | + V = inv(rho) # [m^3/kg] |
| 86 | + |
| 87 | + p = 277.522 # [Pa] |
| 88 | + |
| 89 | + # invert for temperature given p, V |
| 90 | + T = temperature_given_Vp(V, p, eos; initial_T = 251.05, # [K] |
| 91 | + tol = 100 * eps(RealT), maxiter = 100) |
| 92 | + |
| 93 | + a = Trixi.speed_of_sound(V, T, eos) |
| 94 | + M = 6.0 # [1] |
| 95 | + v1 = M * a # [m/s] |
| 96 | + v2 = 0.0 # [m/s] |
| 97 | + |
| 98 | + return thermo2cons(SVector(V, v1, v2, T), equations) |
| 99 | +end |
| 100 | + |
| 101 | +@inline function boundary_condition_supersonic_inflow(u_inner, |
| 102 | + normal_direction::AbstractVector, |
| 103 | + x, t, |
| 104 | + surface_flux_function, |
| 105 | + equations) |
| 106 | + u_boundary = initial_condition_mach6_flow(x, t, equations) |
| 107 | + return flux(u_boundary, normal_direction, equations) |
| 108 | +end |
| 109 | + |
| 110 | +# For physical significance of boundary conditions, see sketch at `mapping_cylinder_shock_fitted` |
| 111 | +boundary_conditions = (; x_neg = boundary_condition_supersonic_inflow, # Supersonic inflow |
| 112 | + y_neg = boundary_condition_slip_wall, # Induce symmetry by slip wall |
| 113 | + y_pos = boundary_condition_do_nothing, # Free outflow |
| 114 | + x_pos = boundary_condition_slip_wall) # Cylinder |
| 115 | + |
| 116 | +############################################################################### |
| 117 | +# Equations, mesh and solver |
| 118 | + |
| 119 | +# The default values correspond to air with temperature bounds/intervals [200.0, 1000.0, 6000.0], see |
| 120 | +# https://ntrs.nasa.gov/api/citations/20020085330/downloads/20020085330.pdf |
| 121 | +eos = ThermallyPerfectGas9PolyFit() |
| 122 | + |
| 123 | +equations = NonIdealCompressibleEulerEquations2D(eos) |
| 124 | + |
| 125 | +# Reduce tolerance to speed things up (otherwise, `eos_newton_maxiter` would need to be increased) |
| 126 | +Trixi.eos_newton_tol(eos::ThermallyPerfectGas9PolyFit) = 1e-5 |
| 127 | + |
| 128 | +polydeg = 3 |
| 129 | +basis = LobattoLegendreBasis(polydeg) |
| 130 | + |
| 131 | +surface_flux = flux_lax_friedrichs |
| 132 | +volume_flux = flux_terashima_etal |
| 133 | + |
| 134 | +shock_indicator = IndicatorHennemannGassner(equations, basis, |
| 135 | + alpha_max = 1.0, |
| 136 | + alpha_min = 0.001, |
| 137 | + alpha_smooth = true, |
| 138 | + variable = density_pressure) |
| 139 | +volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; |
| 140 | + volume_flux_dg = volume_flux, |
| 141 | + volume_flux_fv = surface_flux) |
| 142 | + |
| 143 | +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, |
| 144 | + volume_integral = volume_integral) |
| 145 | + |
| 146 | +trees_per_dimension = (20, 16) |
| 147 | + |
| 148 | +cylinder_radius = 0.5 |
| 149 | +# Follow from a-priori known shock shape, originally for first qaudrant, |
| 150 | +# here transformed to second quadrant, see `mapping_cylinder_shock_fitted`. |
| 151 | +spline_points = 0.6 .* [1.32, 1.05, 2.25] |
| 152 | +cylinder_mapping = (xi, eta) -> mapping_cylinder_shock_fitted(xi, eta, |
| 153 | + cylinder_radius, |
| 154 | + spline_points) |
| 155 | + |
| 156 | +mesh = P4estMesh(trees_per_dimension, |
| 157 | + polydeg = polydeg, |
| 158 | + mapping = cylinder_mapping, |
| 159 | + periodicity = false) |
| 160 | + |
| 161 | +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_mach6_flow, solver; |
| 162 | + boundary_conditions = boundary_conditions) |
| 163 | + |
| 164 | +############################################################################### |
| 165 | +# Semidiscretization & callbacks |
| 166 | + |
| 167 | +tspan = (0.0, 1e-3) |
| 168 | +ode = semidiscretize(semi, tspan) |
| 169 | + |
| 170 | +summary_callback = SummaryCallback() |
| 171 | + |
| 172 | +analysis_callback = AnalysisCallback(semi, interval = 5000) |
| 173 | +alive_callback = AliveCallback(alive_interval = 200) |
| 174 | + |
| 175 | +# Add `:gamma` to `extra_node_variables` tuple ... |
| 176 | +extra_node_variables = (:gamma,) |
| 177 | + |
| 178 | +# ... and specify the function `get_node_variable` for this symbol, |
| 179 | +# with first argument matching the symbol (turned into a type via `Val`) for dispatching. |
| 180 | +function Trixi.get_node_variable(::Val{:gamma}, u, mesh, equations, dg, cache) |
| 181 | + n_nodes = nnodes(dg) |
| 182 | + n_elements = nelements(dg, cache) |
| 183 | + # By definition, the variable must be provided at every node of every element! |
| 184 | + # Otherwise, the `SaveSolutionCallback` will crash. |
| 185 | + gamma_array = zeros(eltype(cache.elements), |
| 186 | + n_nodes, n_nodes, # equivalent: `ntuple(_ -> n_nodes, ndims(mesh))...,` |
| 187 | + n_elements) |
| 188 | + |
| 189 | + eos = equations.equation_of_state |
| 190 | + |
| 191 | + # We can accelerate the computation by thread-parallelizing the loop over elements |
| 192 | + # by using the `@threaded` macro. |
| 193 | + Trixi.@threaded for element in eachelement(dg, cache) |
| 194 | + for j in eachnode(dg), i in eachnode(dg) |
| 195 | + u_node = get_node_vars(u, equations, dg, i, j, element) |
| 196 | + |
| 197 | + # Get temperature |
| 198 | + u_node_thermo = cons2thermo(u_node, equations) |
| 199 | + T = u_node_thermo[4] |
| 200 | + |
| 201 | + gamma_array[i, j, element] = Trixi.gamma(T, eos) |
| 202 | + end |
| 203 | + end |
| 204 | + |
| 205 | + return gamma_array |
| 206 | +end |
| 207 | + |
| 208 | +save_solution = SaveSolutionCallback(interval = 5000, |
| 209 | + solution_variables = cons2thermo, |
| 210 | + extra_node_variables = extra_node_variables) |
| 211 | + |
| 212 | +amr_controller = ControllerThreeLevel(semi, shock_indicator; |
| 213 | + base_level = 0, |
| 214 | + med_level = 1, med_threshold = 0.175, |
| 215 | + max_level = 2, max_threshold = 0.35) |
| 216 | + |
| 217 | +amr_callback = AMRCallback(semi, amr_controller, |
| 218 | + interval = 25, |
| 219 | + adapt_initial_condition = true, |
| 220 | + adapt_initial_condition_only_refine = true) |
| 221 | + |
| 222 | +callbacks = CallbackSet(summary_callback, |
| 223 | + analysis_callback, alive_callback, |
| 224 | + save_solution, amr_callback) |
| 225 | + |
| 226 | +############################################################################### |
| 227 | +# Run the simulation |
| 228 | + |
| 229 | +sol = solve(ode, SSPRK43(; thread = Trixi.Threaded()); |
| 230 | + dt = 1e-7, abstol = 1e-4, reltol = 1e-4, |
| 231 | + ode_default_options()..., callback = callbacks); |
0 commit comments