diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl new file mode 100644 index 00000000000..ab5f5c5ab81 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -0,0 +1,188 @@ +using OrdinaryDiffEqLowStorageRK +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 + +# In Margolin et al., the Navier-Stokes equations are given for an +# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu +mu_isotropic() = 0.15 +mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity + +rho_0() = 1 +v() = 1 # Shock speed + +domain_length = 4.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_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale + +""" + 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 + + # Coordinate transformation. See eq. (33) in Margolin et al. (2017) + chi = 2 * exp(y / (2 * l())) + + w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2)) + + rho = rho_0() / w + u = v() * (1 - w) + p = p_0() * 1 / 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_bar() * 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) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + periodicity = (false, true), + n_cells_max = 100_000) + +### Inviscid boundary conditions ### + +# Prescribe pure influx based on initial conditions +function boundary_condition_inflow(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_cons = initial_condition_viscous_shock(x, t, equations) + return flux(u_cons, orientation, equations) +end + +# Completely free outflow +function boundary_condition_outflow(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + # Calculate the boundary flux entirely from the internal solution state + return flux(u_inner, orientation, equations) +end + +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +### Viscous boundary conditions ### +# For the viscous BCs, we use the known analytical solution +velocity_bc = NoSlip() do x, t, equations_parabolic + Trixi.velocity(initial_condition_viscous_shock(x, + t, + equations_parabolic), + equations_parabolic) +end + +heat_bc = Isothermal() do x, t, equations_parabolic + Trixi.temperature(initial_condition_viscous_shock(x, + t, + equations_parabolic), + equations_parabolic) +end + +boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) + +boundary_conditions_parabolic = (x_neg = boundary_condition_parabolic, + x_pos = boundary_condition_parabolic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +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) + +# Admissible stepsize is governed by the diffusive CFL condition. +# Unless the advective cfl number `cfl` is not reduced to e.g. `0.1` +# (which is overly restrictive for this problem), +# the diffusive CFL restricts the timestep for this problem. +stepsize_callback = StepsizeCallback(cfl = 0.2, + cfl_diffusive = 0.2) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# Use time integrator tailored to compressible Navier-Stokes +sol = solve(ode, CKLLSRK95_4S(), adaptive = false, dt = 1.0, + save_everystep = false, callback = callbacks); diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl new file mode 100644 index 00000000000..d52bb84295d --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -0,0 +1,190 @@ +using OrdinaryDiffEqLowStorageRK +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 + +# In Margolin et al., the Navier-Stokes equations are given for an +# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu +mu_isotropic() = 0.15 +mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity + +rho_0() = 1 +v() = 1 # Shock speed + +domain_length = 4.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_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale + +""" + 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 + + # Coordinate transformation. See eq. (33) in Margolin et al. (2017) + chi = 2 * exp(y / (2 * l())) + + w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2)) + + rho = rho_0() / w + u = v() * (1 - w) + p = p_0() * 1 / 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_bar() * 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) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + periodicity = (false, true, true), + n_cells_max = 100_000) + +### Inviscid boundary conditions ### + +# Prescribe pure influx based on initial conditions +function boundary_condition_inflow(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + u_cons = initial_condition_viscous_shock(x, t, equations) + return flux(u_cons, orientation, equations) +end + +# Completely free outflow +function boundary_condition_outflow(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + # Calculate the boundary flux entirely from the internal solution state + return flux(u_inner, orientation, equations) +end + +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) + +### Viscous boundary conditions ### +# For the viscous BCs, we use the known analytical solution +velocity_bc = NoSlip() do x, t, equations_parabolic + Trixi.velocity(initial_condition_viscous_shock(x, + t, + equations_parabolic), + equations_parabolic) +end + +heat_bc = Isothermal() do x, t, equations_parabolic + Trixi.temperature(initial_condition_viscous_shock(x, + t, + equations_parabolic), + equations_parabolic) +end + +boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) + +boundary_conditions_parabolic = (x_neg = boundary_condition_parabolic, + x_pos = boundary_condition_parabolic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) + +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) + +# For this setup, both advective and diffusive time step restrictions are relevant, i.e., +# may not be increased beyond the given values. +stepsize_callback = StepsizeCallback(cfl = 0.4, + cfl_diffusive = 0.2) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# Use time integrator tailored to compressible Navier-Stokes +sol = solve(ode, CKLLSRK95_4S(), adaptive = false, dt = 1.0, + save_everystep = false, callback = callbacks); diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index 61a9eb2e7f5..ef3a4a50e4c 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -33,19 +33,18 @@ function max_dt(u, t, mesh::TreeMesh{1}, constant_diffusivity::False, equations, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) - # to avoid a division by zero if the speed vanishes everywhere, - # e.g. for steady-state linear advection + # to avoid a division by zero if the diffusivity vanishes everywhere max_scaled_speed = nextfloat(zero(t)) @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) - max_lambda1 = zero(max_scaled_speed) + max_diffusivity_ = zero(max_scaled_speed) for i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, element) - lambda1, = max_diffusivity(u_node, equations_parabolic) - max_lambda1 = max(max_lambda1, lambda1) + diffusivity = max_diffusivity(u_node, equations_parabolic) + max_diffusivity_ = max(max_diffusivity_, diffusivity) end inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx - max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_lambda1) + max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_) end # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2 @@ -72,18 +71,20 @@ function max_dt(u, t, mesh::TreeMesh{1}, return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::TreeMesh{1}, +function max_dt(u, t, mesh::TreeMesh, constant_diffusivity::True, equations, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) - # to avoid a division by zero if the speed vanishes everywhere, - # e.g. for steady-state linear advection + # to avoid a division by zero if the diffusivity vanishes everywhere max_scaled_speed = nextfloat(zero(t)) - max_lambda1, = max_diffusivity(equations_parabolic) + max_lambda1 = max_diffusivity(equations_parabolic) @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx + # Note: For the currently supported parabolic equations + # Diffusion & Navier-Stokes, we only have one diffusivity, + # so this is valid for 1D, 2D and 3D. max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_lambda1) end diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index d0bd3d49f32..28268be65bd 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -29,6 +29,30 @@ function max_dt(u, t, mesh::TreeMesh{2}, return 2 / (nnodes(dg) * max_scaled_speed) end +function max_dt(u, t, mesh::TreeMesh{2}, + constant_diffusivity::False, equations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) + # to avoid a division by zero if the diffusivity vanishes everywhere + max_scaled_speed = nextfloat(zero(t)) + + @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) + max_diffusivity_ = zero(max_scaled_speed) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + # Note: For the currently supported parabolic equations + # Diffusion & Navier-Stokes, we only have one diffusivity. + diffusivity = max_diffusivity(u_node, equations_parabolic) + max_diffusivity_ = max(max_diffusivity_, diffusivity) + end + inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx + max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_) + end + + # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2 + return 4 / (nnodes(dg) * max_scaled_speed) +end + function max_dt(u, t, mesh::TreeMesh{2}, constant_speed::True, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, @@ -82,6 +106,10 @@ function max_dt(u, t, mesh::ParallelTreeMesh{2}, return dt end +# Hyperbolic-parabolic simulations are not yet supported on MPI-parallel meshes. +# Thus, there is no `max_dt` function for `ParallelTreeMesh{2}` and +# `equations_parabolic::AbstractEquationsParabolic` implemented. + function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}, StructuredMeshView{2}}, diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index dc1e7af6be8..dba95999d3a 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -31,6 +31,30 @@ function max_dt(u, t, mesh::TreeMesh{3}, return 2 / (nnodes(dg) * max_scaled_speed) end +function max_dt(u, t, mesh::TreeMesh{3}, + constant_diffusivity::False, equations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) + # to avoid a division by zero if the diffusivity vanishes everywhere + max_scaled_speed = nextfloat(zero(t)) + + @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) + max_diffusivity_ = zero(max_scaled_speed) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + # Note: For the currently supported parabolic equations + # Diffusion & Navier-Stokes, we only have one diffusivity. + diffusivity = max_diffusivity(u_node, equations_parabolic) + max_diffusivity_ = max(max_diffusivity_, diffusivity) + end + inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx + max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_) + end + + # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2 + return 4 / (nnodes(dg) * max_scaled_speed) +end + function max_dt(u, t, mesh::TreeMesh{3}, constant_speed::True, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 96f00c866e7..315dd9c5d82 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -96,6 +96,7 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, Mu, mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law + max_4over3_kappa::RealT # max(4/3, kappa) used for diffusive CFL => `max_diffusivity` equations_hyperbolic::E # CompressibleEulerEquations2D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy @@ -118,6 +119,7 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, mu, Prandtl, kappa, + max(4 / 3, kappa), equations, gradient_variables) end @@ -200,6 +202,34 @@ function flux(u, gradients, orientation::Integer, end end +@inline function max_diffusivity(u, + equations_parabolic::CompressibleNavierStokesDiffusion2D) + # For the diffusive estimate we use the eigenvalues of the diffusivity matrix, + # as suggested in Section 3.5 of + # + # FLEXI: A high order discontinuous Galerkin framework for hyperbolic–parabolic conservation laws + # https://doi.org/10.1016/j.camwa.2020.05.004 + # + # For details on the derivation of eigenvalues of the diffusivity matrix + # for the compressible Navier-Stokes equations see for instance + # + # - Richard P. Dwight (2006) + # Efficiency improvements of RANS-based analysis and optimization using implicit and adjoint methods on unstructured grids + # PhD Thesis, University of Manchester + # https://elib.dlr.de/50794/1/rdwight-PhDThesis-ImplicitAndAdjoint.pdf + # See especially equations (2.79), (3.24), and (3.25) from Chapter 3.2.3 + # + # The eigenvalues of the diffusivity matrix in 2D are + # -mu/rho .* {0, 4/3, 1, kappa} + # + # See for instance also the computation in FLUXO: + # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L122-L128 + # + # Accordingly, the spectral radius/largest absolute eigenvalue can be computed as: + return dynamic_viscosity(u, equations_parabolic) / u[1] * + equations_parabolic.max_4over3_kappa +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion2D) rho, rho_v1, rho_v2, _ = u diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 6c615a11ced..0f96022ecea 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -96,6 +96,7 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu, mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law + max_4over3_kappa::RealT # max(4/3, kappa) used for diffusive CFL => `max_diffusivity` equations_hyperbolic::E # CompressibleEulerEquations3D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy @@ -118,6 +119,7 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, mu, Prandtl, kappa, + max(4 / 3, kappa), equations, gradient_variables) end @@ -225,6 +227,34 @@ function flux(u, gradients, orientation::Integer, end end +@inline function max_diffusivity(u, + equations_parabolic::CompressibleNavierStokesDiffusion3D) + # For the diffusive estimate we use the eigenvalues of the diffusivity matrix, + # as suggested in Section 3.5 of + # + # FLEXI: A high order discontinuous Galerkin framework for hyperbolic–parabolic conservation laws + # https://doi.org/10.1016/j.camwa.2020.05.004 + # + # For details on the derivation of eigenvalues of the diffusivity matrix + # for the compressible Navier-Stokes equations see for instance + # + # - Richard P. Dwight (2006) + # Efficiency improvements of RANS-based analysis and optimization using implicit and adjoint methods on unstructured grids + # PhD Thesis, University of Manchester + # https://elib.dlr.de/50794/1/rdwight-PhDThesis-ImplicitAndAdjoint.pdf + # See especially equations (2.79), (3.24), and (3.25) from Chapter 3.2.3 + # + # The eigenvalues of the diffusivity matrix in 3D are + # -mu/rho .* {0, 4/3, 1, 1, kappa} + # + # See for instance also the computation in FLUXO: + # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L122-L128 + # + # Accordingly, the spectral radius/largest absolute eigenvalue can be computed as: + return dynamic_viscosity(u, equations_parabolic) / u[1] * + equations_parabolic.max_4over3_kappa +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion3D) rho, rho_v1, rho_v2, rho_v3, _ = u diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 2dc8f8828f9..5888bb837f9 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -488,6 +488,23 @@ end tspan=(0.0, 1.0)) end +@trixi_testset "TreeMesh2D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + l2=[ + 2.817640352994614e-5, + 1.3827801939742e-5, + 3.1001993851549174e-17, + 1.7535689010948764e-5 + ], + linf=[ + 0.0002185837290411552, + 0.00013405261969601234, + 1.8273738729889617e-16, + 0.00015782934605046428 + ]) +end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"), diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 34a7cfd337e..a346671127b 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -434,6 +434,25 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "TreeMesh3D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_3d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + l2=[ + 0.0002576235757916208, + 0.00014336914033937938, + 3.361746364570895e-17, + 3.1399702631471645e-17, + 0.00017369856897561295 + ], + linf=[ + 0.0016731997562187129, + 0.0010638567566626511, + 1.733671234158084e-16, + 1.9060786274399122e-16, + 0.001149518946967798 + ]) +end + @trixi_testset "P4estMesh3D: elixir_navierstokes_blast_wave_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_navierstokes_blast_wave_amr.jl"),