From c1db7ac97963fd362e9807ba1f042bd951ddeb64 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 28 Jul 2023 13:17:50 +0200 Subject: [PATCH 01/18] Remove doubled implementations --- src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 54 ---------------------- 1 file changed, 54 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 5370c927e05..6439cad69bb 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -563,60 +563,6 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -# # Function barrier for type stability -# !!! TODO: Figure out why this cannot removed eventhough it exists in the dg_2d_parabolic.jl file -function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh, - equations, surface_integral, dg::DG) - (; boundary_condition_types, boundary_indices) = boundary_conditions - - calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, - Gradient(), mesh, equations, surface_integral, dg) - return nothing -end - -function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh, - equations, surface_integral, dg::DG) - (; boundary_condition_types, boundary_indices) = boundary_conditions - - calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, - Divergence(), mesh, equations, surface_integral, dg) - return nothing -end - -# Iterate over tuples of boundary condition types and associated indices -# in a type-stable way using "lispy tuple programming". -function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, - BC_indices::NTuple{N, Vector{Int}}, - operator_type, - mesh::P4estMesh, - equations, surface_integral, dg::DG) where {N} - # Extract the boundary condition type and index vector - boundary_condition = first(BCs) - boundary_condition_indices = first(BC_indices) - # Extract the remaining types and indices to be processed later - remaining_boundary_conditions = Base.tail(BCs) - remaining_boundary_condition_indices = Base.tail(BC_indices) - - # process the first boundary condition type - calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, - operator_type, mesh, equations, surface_integral, dg) - - # recursively call this method with the unprocessed boundary types - calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions, - remaining_boundary_condition_indices, - operator_type, - mesh, equations, surface_integral, dg) - - return nothing -end - -# terminate the type-stable iteration over tuples -function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, - operator_type, mesh::P4estMesh, equations, - surface_integral, dg::DG) - nothing -end - function calc_boundary_flux!(cache, t, boundary_condition_parabolic, # works with Dict types boundary_condition_indices, From be3b0675ff98e908ef58413a0fd5b41dfb695ac4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 28 Jul 2023 13:20:16 +0200 Subject: [PATCH 02/18] kepp main updated with true main --- src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 54 ++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 6439cad69bb..5370c927e05 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -563,6 +563,60 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end +# # Function barrier for type stability +# !!! TODO: Figure out why this cannot removed eventhough it exists in the dg_2d_parabolic.jl file +function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Gradient(), mesh, equations, surface_integral, dg) + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Divergence(), mesh, equations, surface_integral, dg) + return nothing +end + +# Iterate over tuples of boundary condition types and associated indices +# in a type-stable way using "lispy tuple programming". +function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, + BC_indices::NTuple{N, Vector{Int}}, + operator_type, + mesh::P4estMesh, + equations, surface_integral, dg::DG) where {N} + # Extract the boundary condition type and index vector + boundary_condition = first(BCs) + boundary_condition_indices = first(BC_indices) + # Extract the remaining types and indices to be processed later + remaining_boundary_conditions = Base.tail(BCs) + remaining_boundary_condition_indices = Base.tail(BC_indices) + + # process the first boundary condition type + calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, + operator_type, mesh, equations, surface_integral, dg) + + # recursively call this method with the unprocessed boundary types + calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions, + remaining_boundary_condition_indices, + operator_type, + mesh, equations, surface_integral, dg) + + return nothing +end + +# terminate the type-stable iteration over tuples +function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, + operator_type, mesh::P4estMesh, equations, + surface_integral, dg::DG) + nothing +end + function calc_boundary_flux!(cache, t, boundary_condition_parabolic, # works with Dict types boundary_condition_indices, From 06a419612e785b6bf62f8ac63f845a0300f6bf39 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 31 Jul 2023 16:50:56 +0200 Subject: [PATCH 03/18] Avoid allocations in parabolic boundary fluxes --- examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl | 3 ++- examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl | 5 ++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 36a9f52e39d..4208929f018 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -170,7 +170,8 @@ end initial_condition = initial_condition_navier_stokes_convergence_test # BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +velocity_bc_top_bottom = NoSlip((x, t, equations) -> SVector(initial_condition_navier_stokes_convergence_test(x, t, equations)[2], + initial_condition_navier_stokes_convergence_test(x, t, equations)[3])) heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index b32355c48df..7ca7ab35121 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -220,7 +220,10 @@ end initial_condition = initial_condition_navier_stokes_convergence_test # BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +velocity_bc_top_bottom = NoSlip((x, t, equations) -> SVector(initial_condition_navier_stokes_convergence_test(x, t, equations)[2], + initial_condition_navier_stokes_convergence_test(x, t, equations)[3], + initial_condition_navier_stokes_convergence_test(x, t, equations)[4]) + ) heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) From 40b2796b78dbd9ef3de3e4216c81a603a3a779fd Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 3 Aug 2023 09:11:18 +0200 Subject: [PATCH 04/18] Correct shear layer IC --- examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl index a7cb2fc89f1..dd26fd8097b 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl @@ -14,14 +14,16 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number()) function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) + # Shear layer parameters k = 80 delta = 0.05 u0 = 1.0 + Ms = 0.1 # maximum Mach number rho = 1.0 - v1 = x[2] <= 0.5 ? u0*tanh(k*(x[2]*0.5 - 0.25)) : tanh(k*(0.75 -x[2]*0.5)) - v2 = u0*delta * sin(2*pi*(x[1]*0.5 + 0.25)) + v1 = x[2] <= 0.5 ? u0 * tanh(k*(x[2]*0.5 - 0.25)) : u0 * tanh(k*(0.75 -x[2]*0.5)) + v2 = u0 * delta * sin(2*pi*(x[1]*0.5 + 0.25)) p = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms return prim2cons(SVector(rho, v1, v2, p), equations) From 80519b1a024243eec3ca4a635d5650d2d8db57de Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 09:03:52 +0200 Subject: [PATCH 05/18] Whitespaces --- src/solvers/dgsem_tree/containers_2d.jl | 6 +++--- src/solvers/dgsem_tree/containers_3d.jl | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 5cf256d3499..d80522d42fd 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -764,10 +764,10 @@ end # Container data structure (structure-of-arrays style) for DG MPI interfaces mutable struct MPIInterfaceContainer2D{uEltype <: Real} <: AbstractContainer - u::Array{uEltype, 4} # [leftright, variables, i, interfaces] + u::Array{uEltype, 4} # [leftright, variables, i, interfaces] local_neighbor_ids::Vector{Int} # [interfaces] - orientations::Vector{Int} # [interfaces] - remote_sides::Vector{Int} # [interfaces] + orientations::Vector{Int} # [interfaces] + remote_sides::Vector{Int} # [interfaces] # internal `resize!`able storage _u::Vector{uEltype} end diff --git a/src/solvers/dgsem_tree/containers_3d.jl b/src/solvers/dgsem_tree/containers_3d.jl index 0318946e34d..5fc027ad001 100644 --- a/src/solvers/dgsem_tree/containers_3d.jl +++ b/src/solvers/dgsem_tree/containers_3d.jl @@ -520,14 +520,14 @@ end # Left and right are used *both* for the numbering of the mortar faces *and* for the position of the # elements with respect to the axis orthogonal to the mortar. mutable struct L2MortarContainer3D{uEltype <: Real} <: AbstractContainer - u_upper_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] + u_upper_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_upper_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] - u_lower_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] + u_lower_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_lower_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] - neighbor_ids::Array{Int, 2} # [position, mortars] + neighbor_ids::Array{Int, 2} # [position, mortars] # Large sides: left -> 1, right -> 2 large_sides::Vector{Int} # [mortars] - orientations::Vector{Int} # [mortars] + orientations::Vector{Int} # [mortars] # internal `resize!`able storage _u_upper_left::Vector{uEltype} _u_upper_right::Vector{uEltype} From a2632899338cc07e14616321e71e4cceafd48920 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 16:11:40 +0200 Subject: [PATCH 06/18] Restore main --- examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl | 3 +-- examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl | 6 ++---- examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl | 5 +---- src/solvers/dgsem_tree/containers_2d.jl | 6 +++--- src/solvers/dgsem_tree/containers_3d.jl | 8 ++++---- 5 files changed, 11 insertions(+), 17 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 4208929f018..36a9f52e39d 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -170,8 +170,7 @@ end initial_condition = initial_condition_navier_stokes_convergence_test # BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> SVector(initial_condition_navier_stokes_convergence_test(x, t, equations)[2], - initial_condition_navier_stokes_convergence_test(x, t, equations)[3])) +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl index dd26fd8097b..a7cb2fc89f1 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl @@ -14,16 +14,14 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number()) function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) - # Shear layer parameters k = 80 delta = 0.05 u0 = 1.0 - Ms = 0.1 # maximum Mach number rho = 1.0 - v1 = x[2] <= 0.5 ? u0 * tanh(k*(x[2]*0.5 - 0.25)) : u0 * tanh(k*(0.75 -x[2]*0.5)) - v2 = u0 * delta * sin(2*pi*(x[1]*0.5 + 0.25)) + v1 = x[2] <= 0.5 ? u0*tanh(k*(x[2]*0.5 - 0.25)) : tanh(k*(0.75 -x[2]*0.5)) + v2 = u0*delta * sin(2*pi*(x[1]*0.5 + 0.25)) p = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms return prim2cons(SVector(rho, v1, v2, p), equations) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index 7ca7ab35121..b32355c48df 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -220,10 +220,7 @@ end initial_condition = initial_condition_navier_stokes_convergence_test # BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> SVector(initial_condition_navier_stokes_convergence_test(x, t, equations)[2], - initial_condition_navier_stokes_convergence_test(x, t, equations)[3], - initial_condition_navier_stokes_convergence_test(x, t, equations)[4]) - ) +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index d80522d42fd..5cf256d3499 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -764,10 +764,10 @@ end # Container data structure (structure-of-arrays style) for DG MPI interfaces mutable struct MPIInterfaceContainer2D{uEltype <: Real} <: AbstractContainer - u::Array{uEltype, 4} # [leftright, variables, i, interfaces] + u::Array{uEltype, 4} # [leftright, variables, i, interfaces] local_neighbor_ids::Vector{Int} # [interfaces] - orientations::Vector{Int} # [interfaces] - remote_sides::Vector{Int} # [interfaces] + orientations::Vector{Int} # [interfaces] + remote_sides::Vector{Int} # [interfaces] # internal `resize!`able storage _u::Vector{uEltype} end diff --git a/src/solvers/dgsem_tree/containers_3d.jl b/src/solvers/dgsem_tree/containers_3d.jl index 5fc027ad001..2d8d1e63a2c 100644 --- a/src/solvers/dgsem_tree/containers_3d.jl +++ b/src/solvers/dgsem_tree/containers_3d.jl @@ -520,13 +520,13 @@ end # Left and right are used *both* for the numbering of the mortar faces *and* for the position of the # elements with respect to the axis orthogonal to the mortar. mutable struct L2MortarContainer3D{uEltype <: Real} <: AbstractContainer - u_upper_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] + u_upper_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_upper_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] - u_lower_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] + u_lower_left::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_lower_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] - neighbor_ids::Array{Int, 2} # [position, mortars] + neighbor_ids::Array{Int, 2} # [position, mortars] # Large sides: left -> 1, right -> 2 - large_sides::Vector{Int} # [mortars] + large_sides::Vector{Int} # [mortars] orientations::Vector{Int} # [mortars] # internal `resize!`able storage _u_upper_left::Vector{uEltype} From 520332ddfb7ba06cc72e9fe8004af9114aa4b90b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 16:12:32 +0200 Subject: [PATCH 07/18] restore main --- src/solvers/dgsem_tree/containers_3d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/containers_3d.jl b/src/solvers/dgsem_tree/containers_3d.jl index 2d8d1e63a2c..0318946e34d 100644 --- a/src/solvers/dgsem_tree/containers_3d.jl +++ b/src/solvers/dgsem_tree/containers_3d.jl @@ -526,8 +526,8 @@ mutable struct L2MortarContainer3D{uEltype <: Real} <: AbstractContainer u_lower_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] neighbor_ids::Array{Int, 2} # [position, mortars] # Large sides: left -> 1, right -> 2 - large_sides::Vector{Int} # [mortars] - orientations::Vector{Int} # [mortars] + large_sides::Vector{Int} # [mortars] + orientations::Vector{Int} # [mortars] # internal `resize!`able storage _u_upper_left::Vector{uEltype} _u_upper_right::Vector{uEltype} From dbfbb977c2239480d33c164d6a8120a7d1030f1c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 16:26:50 +0200 Subject: [PATCH 08/18] 1D Navier Stokes --- .../elixir_navierstokes_convergence.jl | 136 +++++ src/Trixi.jl | 3 +- .../compressible_navier_stokes_1d.jl | 482 ++++++++++++++++++ .../compressible_navier_stokes_2d.jl | 71 --- src/equations/equations_parabolic.jl | 1 + test/test_parabolic_1d.jl | 8 +- 6 files changed, 628 insertions(+), 73 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl create mode 100644 src/equations/compressible_navier_stokes_1d.jl diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..3f72d319b0b --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,136 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number()) + +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +# (Simplified version of the 2D) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_t) + v1 = sin(pi_x) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_navier_stokes_convergence_test + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_t) + + v1 = sin(pi_x) * cos(pi_t) + v1_t = -pi * sin(pi_x) * sin(pi_t) + v1_x = pi * cos(pi_x) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * cos(pi_t) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + # stress tensor from x-direction + - v1_xx * mu_) + + # total energy equation + du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + # stress tensor and temperature gradient terms from x-direction + - v1_xx * v1 * mu_ + - v1_x * v1_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_) + + return SVector(du1, du2, du3) +end + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=100_000) + + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver, + source_terms = source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval,) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-9 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index 990c33f3c94..78ddaa3ca7f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -152,7 +152,8 @@ export AcousticPerturbationEquations2D, LinearizedEulerEquations2D export LaplaceDiffusion1D, LaplaceDiffusion2D, - CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D + CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, + CompressibleNavierStokesDiffusion3D export GradientVariablesPrimitive, GradientVariablesEntropy diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl new file mode 100644 index 00000000000..ad167fc2783 --- /dev/null +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -0,0 +1,482 @@ +@doc raw""" + CompressibleNavierStokesDiffusion1D(equations; mu, Pr, + gradient_variables=GradientVariablesPrimitive()) + +Contains the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective terms from +the [`CompressibleEulerEquations1D`](@ref). + +- `equations`: instance of the [`CompressibleEulerEquations1D`](@ref) +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `gradient_variables`: which variables the gradients are taken with respect to. + Defaults to `GradientVariablesPrimitive()`. + +Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., +[``\mu``] = kg m⁻¹ s⁻¹. + +The particular form of the compressible Navier-Stokes implemented is +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho v \\ \rho e +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v \\ \rho v^2 + p \\ (\rho e + p) v +\end{pmatrix} += +\frac{\partial}{\partial x} +\begin{pmatrix} +0 \\ \tau \\ \tau v - q +\end{pmatrix} +``` +where the system is closed with the ideal gas assumption giving +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho v^2 \right) +``` +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations1D`](@ref). +The terms on the right hand side of the system above +are built from the viscous stress +```math +\tau = \mu \frac{\partial}{\partial x} v +``` +where the heat flux is +```math +q = -\kappa \frac{\partial}{\partial x} \left(T\right),\quad T = \frac{p}{R\rho} +``` +where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fick's law. +Under the assumption that the gas has a constant Prandtl number, +the thermal conductivity is +```math +\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}. +``` +From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see +that the gas constant `R` cancels and the heat flux becomes +```math +q = -\kappa \frac{\partial}{\partial x} \left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}} \frac{\partial}{\partial x} \left(\frac{p}{\rho}\right) +``` +which is the form implemented below in the [`flux`](@ref) function. + +In one spatial dimensions we require gradients for two quantities, e.g., +primitive quantities +```math +\frac{\partial}{\partial x} v,\, \frac{\partial}{\partial x} T +``` +or the entropy variables +```math +\frac{\partial}{\partial x} w_2,\, \frac{\partial}{\partial x} w_3 +``` +where +```math +w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} +``` + +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. +""" +struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, + E <: AbstractCompressibleEulerEquations{1}} <: + AbstractCompressibleNavierStokesDiffusion{1, 3} + # TODO: parabolic + # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations + # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + mu::RealT # viscosity + Pr::RealT # Prandtl number + kappa::RealT # thermal diffusivity for Fick's law + + equations_hyperbolic::E # CompressibleEulerEquations1D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +""" +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. + +`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters +for `CompressibleNavierStokesDiffusion1D`. By default, the gradient variables are set to be +`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable +formulation from +- Hughes, Mallet, Franca (1986) + A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the + compressible Euler and Navier-Stokes equations and the second law of thermodynamics. + [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + +Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. +""" +struct GradientVariablesPrimitive end +struct GradientVariablesEntropy end + +# default to primitive gradient variables +function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquations1D; + mu, Prandtl, + gradient_variables = GradientVariablesPrimitive()) + gamma = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + μ, Pr = promote(mu, Prandtl) + + # Under the assumption of constant Prandtl number the thermal conductivity + # constant is kappa = gamma μ / ((gamma-1) Pr). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * inv_gamma_minus_one / Pr + + CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma), + typeof(equations)}(gamma, inv_gamma_minus_one, + μ, Pr, kappa, + equations, gradient_variables) +end + +# TODO: parabolic +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion1D) = ("v1", "v2", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion1D) = ("w2", "w3", "w4") + +function varnames(variable_mapping, + equations_parabolic::CompressibleNavierStokesDiffusion1D) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) +end + +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +function gradient_variable_transformation(::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + cons2prim +end +function gradient_variable_transformation(::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + cons2entropy +end + +# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramírez, Hennemann, Hindenlang, Winters, and Gassner +# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive +# MHD Equations. Part II: Subcell Finite Volume Shock Capturing" +# where one sets the magnetic field components equal to 0. +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) + # 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. + _, dv1dx, dTdx = convert_derivative_to_primitive(u, gradients, equations) + + # Viscous stress (tensor) + tau_11 = dv1dx + + # Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho)) + # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) + # Note, the gas constant cancels under this formulation, so it is not present + # in the implementation + q1 = equations.kappa * dTdx + + # Constant dynamic viscosity is copied to a variable for readability. + # Offers flexibility for dynamic viscosity via Sutherland's law where it depends + # on temperature and reference values, Ts and Tref such that mu(T) + mu = equations.mu + + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = (v1 * tau_11 + q1) * mu + + return SVector(f1, f2, f3) +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion1D) + rho, rho_v1, _ = u + + v1 = rho_v1 / rho + T = temperature(u, equations) + + return SVector(rho, v1, T) +end + +# Convert conservative variables to entropy +# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms +# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion1D`, +# but this may be confusing to new users. +function cons2entropy(u, equations::CompressibleNavierStokesDiffusion1D) + cons2entropy(u, equations.equations_hyperbolic) +end +function entropy2cons(w, equations::CompressibleNavierStokesDiffusion1D) + entropy2cons(w, equations.equations_hyperbolic) +end + +# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. +# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# variables into primitive variables. +@inline function convert_transformed_to_primitive(u_transformed, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return u_transformed +end + +# TODO: parabolic. Make this more efficient! +@inline function convert_transformed_to_primitive(u_transformed, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + # note: this uses CompressibleNavierStokesDiffusion1D versions of cons2prim and entropy2cons + return cons2prim(entropy2cons(u_transformed, equations), equations) +end + +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, T). +# Helpful because then the diffusive fluxes have the same form as on paper. +# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. +# TODO: parabolic; entropy stable viscous terms +@inline function convert_derivative_to_primitive(u, gradient, + ::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return gradient +end + +# the first argument is always the "transformed" variables. +@inline function convert_derivative_to_primitive(w, gradient_entropy_vars, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + + # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back. + # We can fix this if we directly compute v1, v2, T from the entropy variables + u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion1D + rho, rho_v1, _ = u + + v1 = rho_v1 / rho + T = temperature(u, equations) + + return SVector(gradient_entropy_vars[1], + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[3]), # grad(u) = T*(grad(w_2)+v1*grad(w_3)) + T * T * gradient_entropy_vars[3]) +end + +# This routine is required because `prim2cons` is called in `initial_condition`, which +# is called with `equations::CompressibleEulerEquations1D`. This means it is inconsistent +# with `cons2prim(..., ::CompressibleNavierStokesDiffusion1D)` as defined above. +# TODO: parabolic. Is there a way to clean this up? +@inline function prim2cons(u, equations::CompressibleNavierStokesDiffusion1D) + prim2cons(u, equations.equations_hyperbolic) +end + +@inline function temperature(u, equations::CompressibleNavierStokesDiffusion1D) + rho, rho_v1, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1^2 / rho) + T = p / rho + return T +end + +#= +@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion1D) + # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v + + omega = vorticity(u, gradients, equations) + return 0.5 * u[1] * omega^2 +end + +@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion1D) + # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. + _, dv1dx, _ = convert_derivative_to_primitive(u, gradients, equations) + + return dv1dx +end +=# + +# TODO: can we generalize this to MHD? +""" + struct BoundaryConditionNavierStokesWall + +Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. +The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended +to be boundary condition types such as the `NoSlip` velocity boundary condition and the +`Adiabatic` or `Isothermal` heat boundary condition. + +!!! warning "Experimental feature" + This is an experimental feature and may change in future releases. +""" +struct BoundaryConditionNavierStokesWall{V, H} + boundary_condition_velocity::V + boundary_condition_heat_flux::H +end + +""" + struct NoSlip + +Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` +should be a function with signature `boundary_value_function(x, t, equations)` +and should return a `SVector{NDIMS}` whose entries are the velocity vector at a +point `x` and time `t`. +""" +struct NoSlip{F} + boundary_value_function::F # value of the velocity vector on the boundary +end + +""" + struct Isothermal + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_function` should be a function with signature +`boundary_value_function(x, t, equations)` and return a scalar value for the +temperature at point `x` and time `t`. +""" +struct Isothermal{F} + boundary_value_function::F # value of the temperature on the boundary +end + +""" + struct Adiabatic + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_normal_flux_function` should be a function with signature +`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the +normal heat flux at point `x` and time `t`. +""" +struct Adiabatic{F} + boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + return SVector(u_inner[1], v1, u_inner[3]) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + # rho, v1, v2, _ = u_inner + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, + t, + equations) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + _, tau_1n, _ = flux_inner # extract fluxes for 2nd and 3rd equations + normal_energy_flux = v1 * tau_1n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, + equations) + return SVector(u_inner[1], v1, T) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return flux_inner +end + +# specialized BC impositions for GradientVariablesEntropy. + +# This should return a SVector containing the boundary values of entropy variables. +# Here, `w_inner` are the transformed variables (e.g., entropy variables). +# +# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions +# for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022. +# DOI: 10.1016/j.jcp.2021.110723 +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + w_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + negative_rho_inv_p = w_inner[3] # w_3 = -rho / p + return SVector(w_inner[1], -v1 * negative_rho_inv_p, negative_rho_inv_p) +end + +# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness. +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + w_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, + t, + equations) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + _, tau_1n, _ = flux_inner # extract fluxes for 2nd and 3rd equations + normal_energy_flux = v1 * tau_1n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + w_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, + equations) + + # the entropy variables w2 = rho * v1 / p = v1 / T = -v1 * w3. + w3 = -1 / T + return SVector(w_inner[1], -v1 * w3, w3) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + w_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + return SVector(flux_inner[1], flux_inner[2], flux_inner[3]) +end diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index a1f11717e69..2f7de151688 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -93,24 +93,6 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy end -""" -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. - -`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters -for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be -`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable -formulation from -- Hughes, Mallet, Franca (1986) - A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the - compressible Euler and Navier-Stokes equations and the second law of thermodynamics. - [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) - -Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. -""" -struct GradientVariablesPrimitive end -struct GradientVariablesEntropy end - # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; mu, Prandtl, @@ -315,59 +297,6 @@ end return dv2dx - dv1dy end -# TODO: can we generalize this to MHD? -""" - struct BoundaryConditionNavierStokesWall - -Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. -The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended -to be boundary condition types such as the `NoSlip` velocity boundary condition and the -`Adiabatic` or `Isothermal` heat boundary condition. - -!!! warning "Experimental feature" - This is an experimental feature and may change in future releases. -""" -struct BoundaryConditionNavierStokesWall{V, H} - boundary_condition_velocity::V - boundary_condition_heat_flux::H -end - -""" - struct NoSlip - -Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` -should be a function with signature `boundary_value_function(x, t, equations)` -and should return a `SVector{NDIMS}` whose entries are the velocity vector at a -point `x` and time `t`. -""" -struct NoSlip{F} - boundary_value_function::F # value of the velocity vector on the boundary -end - -""" - struct Isothermal - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_function` should be a function with signature -`boundary_value_function(x, t, equations)` and return a scalar value for the -temperature at point `x` and time `t`. -""" -struct Isothermal{F} - boundary_value_function::F # value of the temperature on the boundary -end - -""" - struct Adiabatic - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_normal_flux_function` should be a function with signature -`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the -normal heat flux at point `x` and time `t`. -""" -struct Adiabatic{F} - boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn -end - @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 6c0be43798a..194fc824afd 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -11,5 +11,6 @@ include("laplace_diffusion_2d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +include("compressible_navier_stokes_1d.jl") include("compressible_navier_stokes_2d.jl") include("compressible_navier_stokes_3d.jl") diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 1aaf23d576a..563099c541b 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -19,7 +19,13 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [2.847421658558336e-05] ) end - + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence.jl"), + l2 = [0.0001133835907077494, 6.226282245610444e-5, 0.0002820171699999139], + linf = [0.0006255102377159538, 0.00036195501456059986, 0.0016147729485886941] + ) + end end # Clean up afterwards: delete Trixi output directory From a0339cf5fac0cf5462017fa0f3cd1710648c70cc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 16:32:22 +0200 Subject: [PATCH 09/18] Conventional notation for heat flux --- src/equations/compressible_navier_stokes_2d.jl | 6 +++--- src/equations/compressible_navier_stokes_3d.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 2f7de151688..f762fe5d5ee 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -29,7 +29,7 @@ The particular form of the compressible Navier-Stokes implemented is = \nabla \cdot \begin{pmatrix} -0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix} ``` where the system is closed with the ideal gas assumption giving @@ -44,7 +44,7 @@ are built from the viscous stress tensor ``` where ``\underline{I}`` is the ``2\times 2`` identity matrix and the heat flux is ```math -\nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} +\mathbf{q} = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} ``` where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fick's law. Under the assumption that the gas has a constant Prandtl number, @@ -55,7 +55,7 @@ the thermal conductivity is From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see that the gas constant `R` cancels and the heat flux becomes ```math -\nabla q = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) +\mathbf{q} = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) ``` which is the form implemented below in the [`flux`](@ref) function. diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 0b770dff1ca..166b53bf615 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -29,7 +29,7 @@ The particular form of the compressible Navier-Stokes implemented is = \nabla \cdot \begin{pmatrix} -0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix} ``` where the system is closed with the ideal gas assumption giving @@ -44,7 +44,7 @@ are built from the viscous stress tensor ``` where ``\underline{I}`` is the ``3\times 3`` identity matrix and the heat flux is ```math -\nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} +\mathbf{q} = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} ``` where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fick's law. Under the assumption that the gas has a constant Prandtl number, @@ -55,7 +55,7 @@ the thermal conductivity is From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see that the gas constant `R` cancels and the heat flux becomes ```math -\nabla q = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) +\mathbf{q} = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) ``` which is the form implemented below in the [`flux`](@ref) function. From 1b7e55a264414d519f45444bb0f925ce8f0ee680 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 7 Aug 2023 16:35:05 +0200 Subject: [PATCH 10/18] remove multi-dim artefacts --- .../compressible_navier_stokes_1d.jl | 20 ++----------------- 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index ad167fc2783..396c0e14e2e 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -276,22 +276,6 @@ end return T end -#= -@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion1D) - # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v - - omega = vorticity(u, gradients, equations) - return 0.5 * u[1] * omega^2 -end - -@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion1D) - # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. - _, dv1dx, _ = convert_derivative_to_primitive(u, gradients, equations) - - return dv1dx -end -=# - # TODO: can we generalize this to MHD? """ struct BoundaryConditionNavierStokesWall @@ -374,7 +358,7 @@ end equations) v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) - _, tau_1n, _ = flux_inner # extract fluxes for 2nd and 3rd equations + _, tau_1n, _ = flux_inner # extract fluxes for 2nd equation normal_energy_flux = v1 * tau_1n + normal_heat_flux return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) end @@ -445,7 +429,7 @@ end equations) v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) - _, tau_1n, _ = flux_inner # extract fluxes for 2nd and 3rd equations + _, tau_1n, _ = flux_inner # extract fluxes for 2nd equation normal_energy_flux = v1 * tau_1n + normal_heat_flux return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) end From 7e5a7a5698b185e8581c089d7682647280a2df85 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 09:18:38 +0200 Subject: [PATCH 11/18] Move general part into own file --- src/equations/compressible_navier_stokes.jl | 71 +++++++++++++++++++ .../compressible_navier_stokes_1d.jl | 71 ------------------- src/equations/equations_parabolic.jl | 1 + 3 files changed, 72 insertions(+), 71 deletions(-) create mode 100644 src/equations/compressible_navier_stokes.jl diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl new file mode 100644 index 00000000000..becc8073fc0 --- /dev/null +++ b/src/equations/compressible_navier_stokes.jl @@ -0,0 +1,71 @@ +# TODO: can we generalize this to MHD? +""" + struct BoundaryConditionNavierStokesWall + +Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. +The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended +to be boundary condition types such as the `NoSlip` velocity boundary condition and the +`Adiabatic` or `Isothermal` heat boundary condition. + +!!! warning "Experimental feature" + This is an experimental feature and may change in future releases. +""" +struct BoundaryConditionNavierStokesWall{V, H} + boundary_condition_velocity::V + boundary_condition_heat_flux::H +end + +""" + struct NoSlip + +Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` +should be a function with signature `boundary_value_function(x, t, equations)` +and should return a `SVector{NDIMS}` whose entries are the velocity vector at a +point `x` and time `t`. +""" +struct NoSlip{F} + boundary_value_function::F # value of the velocity vector on the boundary +end + +""" + struct Isothermal + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_function` should be a function with signature +`boundary_value_function(x, t, equations)` and return a scalar value for the +temperature at point `x` and time `t`. +""" +struct Isothermal{F} + boundary_value_function::F # value of the temperature on the boundary +end + +""" + struct Adiabatic + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_normal_flux_function` should be a function with signature +`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the +normal heat flux at point `x` and time `t`. +""" +struct Adiabatic{F} + boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn +end + + +""" +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. + +`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters +for `CompressibleNavierStokesDiffusion1D`. By default, the gradient variables are set to be +`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable +formulation from +- Hughes, Mallet, Franca (1986) + A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the + compressible Euler and Navier-Stokes equations and the second law of thermodynamics. + [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + +Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. +""" +struct GradientVariablesPrimitive end +struct GradientVariablesEntropy end \ No newline at end of file diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 396c0e14e2e..5f938dd754d 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -93,24 +93,6 @@ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy end -""" -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. - -`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters -for `CompressibleNavierStokesDiffusion1D`. By default, the gradient variables are set to be -`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable -formulation from -- Hughes, Mallet, Franca (1986) - A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the - compressible Euler and Navier-Stokes equations and the second law of thermodynamics. - [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) - -Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. -""" -struct GradientVariablesPrimitive end -struct GradientVariablesEntropy end - # default to primitive gradient variables function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquations1D; mu, Prandtl, @@ -276,59 +258,6 @@ end return T end -# TODO: can we generalize this to MHD? -""" - struct BoundaryConditionNavierStokesWall - -Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. -The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended -to be boundary condition types such as the `NoSlip` velocity boundary condition and the -`Adiabatic` or `Isothermal` heat boundary condition. - -!!! warning "Experimental feature" - This is an experimental feature and may change in future releases. -""" -struct BoundaryConditionNavierStokesWall{V, H} - boundary_condition_velocity::V - boundary_condition_heat_flux::H -end - -""" - struct NoSlip - -Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` -should be a function with signature `boundary_value_function(x, t, equations)` -and should return a `SVector{NDIMS}` whose entries are the velocity vector at a -point `x` and time `t`. -""" -struct NoSlip{F} - boundary_value_function::F # value of the velocity vector on the boundary -end - -""" - struct Isothermal - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_function` should be a function with signature -`boundary_value_function(x, t, equations)` and return a scalar value for the -temperature at point `x` and time `t`. -""" -struct Isothermal{F} - boundary_value_function::F # value of the temperature on the boundary -end - -""" - struct Adiabatic - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_normal_flux_function` should be a function with signature -`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the -normal heat flux at point `x` and time `t`. -""" -struct Adiabatic{F} - boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn -end - @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 194fc824afd..66214025044 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -11,6 +11,7 @@ include("laplace_diffusion_2d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +include("compressible_navier_stokes.jl") include("compressible_navier_stokes_1d.jl") include("compressible_navier_stokes_2d.jl") include("compressible_navier_stokes_3d.jl") From feebac5f2d7d8c56a0879878c48128b93b2e0751 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 10:15:43 +0200 Subject: [PATCH 12/18] Slip Wall BC for 1D Compressible Euler --- src/equations/compressible_euler_1d.jl | 51 ++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index e4fd0997eae..146bdc6ea8c 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -198,6 +198,57 @@ function initial_condition_eoc_test_coupled_euler_gravity(x, t, return prim2cons(SVector(rho, v1, p), equations) end +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations1D) +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + + Should be used together with [`TreeMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations1D) + # compute the primitive variables + rho_local, v_normal, p_local = cons2prim(u_inner, equations) + + if isodd(direction) # flip sign of normal to make it outward pointing + v_normal *= -1 + end + + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0.0 + A = 2 / ((equations.gamma + 1) * rho_local) + B = p_local * (equations.gamma - 1) / (equations.gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner)), + p_star, + zero(eltype(u_inner))) +end + # Calculate 1D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u From 43b41f12d9f75a0f149b084ffb7ef86c4fd3cb95 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 10:15:59 +0200 Subject: [PATCH 13/18] Correct arguments for 1D BCs --- .../compressible_navier_stokes_1d.jl | 24 ++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 5f938dd754d..7ccad4ea53a 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -261,7 +261,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ @@ -275,7 +276,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ @@ -295,7 +297,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ @@ -311,7 +314,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ @@ -331,7 +335,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ @@ -347,7 +352,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ @@ -366,7 +372,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ @@ -385,7 +392,8 @@ end @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, - normal::AbstractVector, + orientation::Integer, + direction, x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ From f088c13a6c80e28b3240ec73f8a25c0979ab2b84 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 10:28:36 +0200 Subject: [PATCH 14/18] format --- src/equations/compressible_euler_1d.jl | 14 +++++++------- src/equations/compressible_navier_stokes.jl | 3 +-- src/equations/compressible_navier_stokes_1d.jl | 2 +- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 146bdc6ea8c..9204989e8be 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -232,21 +232,21 @@ are available in the paper: if v_normal <= 0.0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * - equations.gamma * - equations.inv_gamma_minus_one) + (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) else # v_normal > 0.0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + - 0.5 * v_normal / A * - (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end # For the slip wall we directly set the flux as the normal velocity is zero return SVector(zero(eltype(u_inner)), - p_star, - zero(eltype(u_inner))) + p_star, + zero(eltype(u_inner))) end # Calculate 1D flux for a single point diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl index becc8073fc0..af7897d4586 100644 --- a/src/equations/compressible_navier_stokes.jl +++ b/src/equations/compressible_navier_stokes.jl @@ -51,7 +51,6 @@ struct Adiabatic{F} boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn end - """ !!! warning "Experimental code" This code is experimental and may be changed or removed in any future release. @@ -68,4 +67,4 @@ formulation from Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. """ struct GradientVariablesPrimitive end -struct GradientVariablesEntropy end \ No newline at end of file +struct GradientVariablesEntropy end diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 7ccad4ea53a..dca846cac1e 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -262,7 +262,7 @@ end <:Adiabatic})(flux_inner, u_inner, orientation::Integer, - direction, + direction, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ From 8346da7f3a90e20d2526eff07a815161b7301e46 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 10:28:49 +0200 Subject: [PATCH 15/18] Add convergence test with walls --- ...ixir_navierstokes_convergence_periodic.jl} | 0 .../elixir_navierstokes_convergence_walls.jl | 154 ++++++++++++++++++ test/test_parabolic_1d.jl | 11 +- 3 files changed, 163 insertions(+), 2 deletions(-) rename examples/tree_1d_dgsem/{elixir_navierstokes_convergence.jl => elixir_navierstokes_convergence_periodic.jl} (100%) create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl similarity index 100% rename from examples/tree_1d_dgsem/elixir_navierstokes_convergence.jl rename to examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl new file mode 100644 index 00000000000..cf19bd7c897 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -0,0 +1,154 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = -1.0 +coordinates_max = 1.0 + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion1D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion1D`) +# and by the initial condition (which passes in `CompressibleEulerEquations1D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * cos(pi_x) * cos(pi_t) + v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0)) ) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + x = x[1] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * cos(pi_x) * cos(pi_t) + rho_t = -pi * A * cos(pi_x) * sin(pi_t) + rho_x = -pi * A * sin(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) + v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) + v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) + v1_xx = (( 2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) + - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # y-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + # stress tensor from y-direction + - v1_xx * mu_) + + # total energy equation + du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + # stress tensor and temperature gradient terms from x-direction + - v1_xx * v1 * mu_ + - v1_x * v1_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ ) + + return SVector(du1, du2, du3) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2]) +heat_bc_left_right = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_left_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left_right) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_left_right, + x_pos = boundary_condition_left_right) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +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-6 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 563099c541b..d5bf4eba8e2 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -20,12 +20,19 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence.jl"), + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), l2 = [0.0001133835907077494, 6.226282245610444e-5, 0.0002820171699999139], linf = [0.0006255102377159538, 0.00036195501456059986, 0.0016147729485886941] ) end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), + l2 = [0.000321596075226452, 0.0003182128802451823, 0.0014390613437572657], + linf = [0.0012224119045414206, 0.0026168886703672534, 0.012041112066421888] + ) + end end # Clean up afterwards: delete Trixi output directory From 4f0180ddef3aeaf7c7b41dc3faa3398e7377ff1c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Aug 2023 17:00:38 +0200 Subject: [PATCH 16/18] Test gradient with entropy variables --- test/test_parabolic_1d.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index d5bf4eba8e2..5295fc4bf78 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -27,6 +27,16 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl: GradientVariablesEntropy" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number(), + gradient_variables = GradientVariablesEntropy()), + l2 = [0.00011310615871043463, 6.216495207074201e-5, 0.00028195843110817814], + linf = [0.0006240837363233886, 0.0003616694320713876, 0.0016147339542413874] + ) + end + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), l2 = [0.000321596075226452, 0.0003182128802451823, 0.0014390613437572657], From 6aeeeea6b1582d93639e80923cb20c8e40ae5cb9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 9 Aug 2023 08:49:00 +0200 Subject: [PATCH 17/18] Test isothermal BC, test gradient in entropy vars --- .../elixir_navierstokes_convergence_walls.jl | 18 ++++++++++++------ test/test_parabolic_1d.jl | 14 ++++++++++++-- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl index cf19bd7c897..181a2cb209f 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -117,16 +117,22 @@ initial_condition = initial_condition_navier_stokes_convergence_test # BC types velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2]) -heat_bc_left_right = Adiabatic((x, t, equations) -> 0.0) -boundary_condition_left_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left_right) + +heat_bc_left = Isothermal((x, t, equations) -> + Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), + equations_parabolic)) +heat_bc_right = Adiabatic((x, t, equations) -> 0.0) + +boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left) +boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_right) # define inviscid boundary conditions boundary_conditions = (; x_neg = boundary_condition_slip_wall, x_pos = boundary_condition_slip_wall) # define viscous boundary conditions -boundary_conditions_parabolic = (; x_neg = boundary_condition_left_right, - x_pos = boundary_condition_left_right) +boundary_conditions_parabolic = (; x_neg = boundary_condition_left, + x_pos = boundary_condition_right) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), @@ -148,7 +154,7 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) ############################################################################### # run the simulation -time_int_tol = 1e-6 +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) -summary_callback() # print the timer summary +summary_callback() # print the timer summary \ No newline at end of file diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 5295fc4bf78..7c45a0e4829 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -39,8 +39,18 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), - l2 = [0.000321596075226452, 0.0003182128802451823, 0.0014390613437572657], - linf = [0.0012224119045414206, 0.0026168886703672534, 0.012041112066421888] + l2 = [0.00047023310868269237, 0.00032181736027057234, 0.0014966266486095025], + linf = [0.002996375101363302, 0.002863904256059634, 0.012691132946258676] + ) + end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl: GradientVariablesEntropy" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number(), + gradient_variables = GradientVariablesEntropy()), + l2 = [0.0004702396543424139, 0.00032222144750623604, 0.0014999908210611393], + linf = [0.002994783988792271, 0.002864514179705944, 0.012717662985073375] ) end end From 955a95eb3b43b3b0c9152a7c75cf9d28302f148c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 9 Aug 2023 09:40:30 +0200 Subject: [PATCH 18/18] Correct test data --- test/test_parabolic_1d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 7c45a0e4829..06a55100d62 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -49,8 +49,8 @@ isdir(outdir) && rm(outdir, recursive=true) equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number(), gradient_variables = GradientVariablesEntropy()), - l2 = [0.0004702396543424139, 0.00032222144750623604, 0.0014999908210611393], - linf = [0.002994783988792271, 0.002864514179705944, 0.012717662985073375] + l2 = [0.0004608500483647771, 0.00032431091222851285, 0.0015159733360626845], + linf = [0.002754803146635787, 0.0028567714697580906, 0.012941794048176192] ) end end