diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl new file mode 100644 index 00000000000..3890846fd63 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -0,0 +1,176 @@ +using OrdinaryDiffEq +using Trixi + +# This is the classic 1D viscous shock wave problem with analytical solution +# for a special value of the Prandtl number. +# The original references are: +# +# - R. Becker (1922) +# Stoßwelle und Detonation. +# [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605) +# +# English translations: +# Impact waves and detonation. Part I. +# https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf +# +# Impact waves and detonation. Part II. +# https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf +# +# - M. Morduchow, P. A. Libby (1949) +# On a Complete Solution of the One-Dimensional Flow Equations +# of a Viscous, Head-Conducting, Compressible Gas +# [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882) +# +# +# The particular problem considered here is described in +# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) +# Entropy in self-similar shock profiles +# [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003) + +### Fixed parameters ### + +# Special value for which nonlinear solver can be omitted +# Corresponds essentially to fixing the Mach number +alpha = 0.5 +# We want kappa = cp * mu = mu_bar to ensure constant enthalpy +prandtl_number() = 3 / 4 + +### Free choices: ### +gamma = 5 / 3 +rho_0 = 1 +# In Margolin et al., the Navier-Stokes equations are given for an +# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu +mu_isotropic = 0.1 +v = 1 # Shock speed + +domain_length = 5.0 + +### Derived quantities ### + +Ma = 2 / sqrt(3 - gamma) # Mach number for alpha = 0.5 +c_0 = v / Ma # Speed of sound ahead of the shock + +# From constant enthalpy condition +p_0 = c_0^2 * rho_0 / gamma + +l = mu_isotropic / (rho_0 * v) * 2 * gamma / (gamma + 1) # Appropriate length scale + +# Helper function for coordinate transformation. See eq. (33) in Margolin et al. (2017) +chi_of_y(y) = 2 * exp(y / (2 * l)) + +""" + initial_condition_viscous_shock(x, t, equations) + +Classic 1D viscous shock wave problem with analytical solution +for a special value of the Prandtl number. +The version implemented here is described in +- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) + Entropy in self-similar shock profiles + [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003) +""" +function initial_condition_viscous_shock(x, t, equations) + y = x[1] - v * t # Translated coordinate + + chi = chi_of_y(y) + w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2)) + + rho = rho_0 / w + u = v * (1 - w) + p = p_0 / w * (1 + (gamma - 1) / 2 * Ma^2 * (1 - w^2)) + + return prim2cons(SVector(rho, u, 0, p), equations) +end +initial_condition = initial_condition_viscous_shock + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +equations = CompressibleEulerEquations2D(gamma) + +# Trixi implements the stress tensor in deviatoric form, thus we need to +# convert the "isotropic viscosity" to the "deviatoric viscosity" +mu_deviatoric() = mu_isotropic * 3 / 4 +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu_deviatoric(), + Prandtl = prandtl_number(), + gradient_variables = GradientVariablesPrimitive()) + +solver = DGSEM(polydeg = 3, surface_flux = flux_hlle) + +coordinates_min = (-domain_length / 2, -domain_length / 2) +coordinates_max = (domain_length / 2, domain_length / 2) + +trees_per_dimension = (8, 2) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 0, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = (false, true)) + +### Inviscid boundary conditions ### + +# Prescribe pure influx based on initial conditions +function boundary_condition_inflow(u_inner, normal_direction, x, t, + surface_flux_function, equations) + u_cons = initial_condition(x, t, equations) + flux = Trixi.flux(u_cons, normal_direction, equations) + + return flux +end + +# Completely free outflow +function boundary_condition_outflow(u_inner, normal_direction, x, t, + surface_flux_function, equations) + # Calculate the boundary flux entirely from the internal solution state + flux = Trixi.flux(u_inner, normal_direction, equations) + + return flux +end + +boundary_conditions = Dict(:x_neg => boundary_condition_inflow, + :x_pos => boundary_condition_outflow) + +### Viscous boundary conditions ### +# For the viscous BCs, we use the known analytical solution +velocity_bc = NoSlip((x, t, equations) -> Trixi.velocity(initial_condition(x, + t, + equations), + equations_parabolic)) + +heat_bc = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition(x, + t, + equations), + equations_parabolic)) + +boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) + +boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic, + :x_pos => boundary_condition_parabolic) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +alive_callback = AliveCallback(alive_interval = 10) + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + dt = 1e-3, ode_default_options()..., callback = callbacks) + +summary_callback() # print the timer summary diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl new file mode 100644 index 00000000000..34efc064476 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -0,0 +1,176 @@ +using OrdinaryDiffEq +using Trixi + +# This is the classic 1D viscous shock wave problem with analytical solution +# for a special value of the Prandtl number. +# The original references are: +# +# - R. Becker (1922) +# Stoßwelle und Detonation. +# [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605) +# +# English translations: +# Impact waves and detonation. Part I. +# https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf +# +# Impact waves and detonation. Part II. +# https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf +# +# - M. Morduchow, P. A. Libby (1949) +# On a Complete Solution of the One-Dimensional Flow Equations +# of a Viscous, Head-Conducting, Compressible Gas +# [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882) +# +# +# The particular problem considered here is described in +# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) +# Entropy in self-similar shock profiles +# [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003) + +### Fixed parameters ### + +# Special value for which nonlinear solver can be omitted +# Corresponds essentially to fixing the Mach number +alpha = 0.5 +# We want kappa = cp * mu = mu_bar to ensure constant enthalpy +prandtl_number() = 3 / 4 + +### Free choices: ### +gamma = 5 / 3 +rho_0 = 1 +# In Margolin et al., the Navier-Stokes equations are given for an +# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu +mu_isotropic = 0.1 +v = 1 # Shock speed + +domain_length = 5.0 + +### Derived quantities ### + +Ma = 2 / sqrt(3 - gamma) # Mach number for alpha = 0.5 +c_0 = v / Ma # Speed of sound ahead of the shock + +# From constant enthalpy condition +p_0 = c_0^2 * rho_0 / gamma + +l = mu_isotropic / (rho_0 * v) * 2 * gamma / (gamma + 1) # Appropriate length scale + +# Helper function for coordinate transformation. See eq. (33) in Margolin et al. (2017) +chi_of_y(y) = 2 * exp(y / (2 * l)) + +""" + initial_condition_viscous_shock(x, t, equations) + +Classic 1D viscous shock wave problem with analytical solution +for a special value of the Prandtl number. +The version implemented here is described in +- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) + Entropy in self-similar shock profiles + [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003) +""" +function initial_condition_viscous_shock(x, t, equations) + y = x[1] - v * t # Translated coordinate + + chi = chi_of_y(y) + w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2)) + + rho = rho_0 / w + u = v * (1 - w) + p = p_0 / w * (1 + (gamma - 1) / 2 * Ma^2 * (1 - w^2)) + + return prim2cons(SVector(rho, u, 0, 0, p), equations) +end +initial_condition = initial_condition_viscous_shock + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +equations = CompressibleEulerEquations3D(gamma) + +# Trixi implements the stress tensor in deviatoric form, thus we need to +# convert the "isotropic viscosity" to the "deviatoric viscosity" +mu_deviatoric() = mu_isotropic * 3 / 4 +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu_deviatoric(), + Prandtl = prandtl_number(), + gradient_variables = GradientVariablesPrimitive()) + +solver = DGSEM(polydeg = 3, surface_flux = flux_hlle) + +coordinates_min = (-domain_length / 2, -domain_length / 2, -domain_length / 2) +coordinates_max = (domain_length / 2, domain_length / 2, domain_length / 2) + +trees_per_dimension = (8, 2, 2) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 0, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = (false, true, true)) + +### Inviscid boundary conditions ### + +# Prescribe pure influx based on initial conditions +function boundary_condition_inflow(u_inner, normal_direction, x, t, + surface_flux_function, equations) + u_cons = initial_condition(x, t, equations) + flux = Trixi.flux(u_cons, normal_direction, equations) + + return flux +end + +# Completely free outflow +function boundary_condition_outflow(u_inner, normal_direction, x, t, + surface_flux_function, equations) + # Calculate the boundary flux entirely from the internal solution state + flux = Trixi.flux(u_inner, normal_direction, equations) + + return flux +end + +boundary_conditions = Dict(:x_neg => boundary_condition_inflow, + :x_pos => boundary_condition_outflow) + +### Viscous boundary conditions ### +# For the viscous BCs, we use the known analytical solution +velocity_bc = NoSlip((x, t, equations) -> Trixi.velocity(initial_condition(x, + t, + equations), + equations_parabolic)) + +heat_bc = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition(x, + t, + equations), + equations_parabolic)) + +boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) + +boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic, + :x_pos => boundary_condition_parabolic) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +alive_callback = AliveCallback(alive_interval = 10) + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + dt = 1e-3, ode_default_options()..., callback = callbacks) + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl index 13e6a7d1d19..d33daceea1a 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -51,7 +51,7 @@ c_0 = v / Ma # Speed of sound ahead of the shock # From constant enthalpy condition p_0 = c_0^2 * rho_0 / gamma -l = mu() / (rho_0 * v) * 2 * gamma / (gamma + 1) # Appropriate lenght scale +l = mu() / (rho_0 * v) * 2 * gamma / (gamma + 1) # Appropriate length scale # Helper function for coordinate transformation. See eq. (33) in Margolin et al. (2017) chi_of_y(y) = 2 * exp(y / (2 * l)) @@ -78,6 +78,7 @@ function initial_condition_viscous_shock(x, t, equations) return prim2cons(SVector(rho, u, p), equations) end +initial_condition = initial_condition_viscous_shock ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations @@ -85,7 +86,7 @@ end equations = CompressibleEulerEquations1D(gamma) equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu(), Prandtl = prandtl_number(), - gradient_variables = GradientVariablesEntropy()) + gradient_variables = GradientVariablesPrimitive()) solver = DGSEM(polydeg = 3, surface_flux = flux_hlle) @@ -93,19 +94,17 @@ coordinates_min = -domain_length / 2 coordinates_max = domain_length / 2 mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, + initial_refinement_level = 3, periodicity = false, n_cells_max = 30_000) -initial_condition = initial_condition_viscous_shock - ### Inviscid boundary conditions ### # Prescribe pure influx based on initial conditions function boundary_condition_inflow(u_inner, orientation, normal_direction, x, t, surface_flux_function, equations) u_cons = initial_condition(x, t, equations) - flux = Trixi.flux(u_cons, normal_direction, equations) + flux = Trixi.flux(u_cons, orientation, equations) return flux end @@ -114,7 +113,7 @@ end function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t, surface_flux_function, equations) # Calculate the boundary flux entirely from the internal solution state - flux = Trixi.flux(u_inner, normal_direction, equations) + flux = Trixi.flux(u_inner, orientation, equations) return flux end @@ -127,7 +126,7 @@ boundary_conditions = (; x_neg = boundary_condition_inflow, velocity_bc = NoSlip((x, t, equations) -> Trixi.velocity(initial_condition(x, t, equations), - equations)) + equations_parabolic)) heat_bc = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition(x, t, @@ -165,6 +164,6 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) time_int_tol = 1e-8 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, - dt = 1e-5, ode_default_options()..., callback = callbacks) + dt = 1e-3, ode_default_options()..., callback = callbacks) summary_callback() # print the timer summary diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 024ce2d7e87..fe79a7e0590 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -148,7 +148,7 @@ end function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion1D) # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. - rho, v1, _ = convert_transformed_to_primitive(u, equations) + _, v1, _ = convert_transformed_to_primitive(u, equations) # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T) # either computed directly or reverse engineered from the gradient of the entropy variables # by way of the `convert_gradient_variables` function. @@ -257,6 +257,12 @@ end return T end +@inline function velocity(u, equations::CompressibleNavierStokesDiffusion1D) + rho = u[1] + v1 = u[2] / rho + return v1 +end + @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 5708abb4e00..977dae18a43 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -148,7 +148,7 @@ end function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion2D) # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. - rho, v1, v2, _ = convert_transformed_to_primitive(u, equations) + _, v1, v2, _ = convert_transformed_to_primitive(u, equations) # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T) # either computed directly or reverse engineered from the gradient of the entropy variables # by way of the `convert_gradient_variables` function. @@ -281,6 +281,13 @@ end return T end +@inline function velocity(u, equations::CompressibleNavierStokesDiffusion2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + return SVector(v1, v2) +end + @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 9622882fc1a..551462e7595 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -148,7 +148,7 @@ end function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion3D) # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. - rho, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations) + _, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations) # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T) # either computed directly or reverse engineered from the gradient of the entropy variables # by way of the `convert_gradient_variables` function. @@ -309,6 +309,14 @@ end return T end +@inline function velocity(u, equations::CompressibleNavierStokesDiffusion3D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D) # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 48b2b9ce4d6..a04009a2908 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -228,6 +228,29 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "TreeMesh1D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + l2=[ + 0.0037611419514324567, + 0.0031130623550597895, + 0.0032646075737842403 + ], + linf=[ + 0.019311368481275126, + 0.013898566510286758, + 0.01496339059101559 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi output directory diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index ceefb65e99b..b607d81b025 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -756,6 +756,31 @@ end @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13) @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13) end + +@trixi_testset "P4estMesh2D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + l2=[ + 0.003761142833904016, + 0.003113063164051044, + 8.268420784251351e-17, + 0.0032646082100107336 + ], + linf=[ + 0.019311388566492393, + 0.013898567062704204, + 4.776829432039756e-16, + 0.014963390663296883 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 2690a08cbb9..0158faa13f7 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -505,6 +505,33 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "P4estMesh3D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + l2=[ + 0.0037611432649265466, + 0.003113063576315902, + 1.78506109857939e-16, + 1.6328152981807113e-16, + 0.0032646085623784424 + ], + linf=[ + 0.01931139690659478, + 0.01389856746278273, + 3.7278358922962055e-15, + 3.903521486285084e-15, + 0.014963390916563846 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory