From d5cf820b7c0e1904537846426b861d775f8f35ea Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 18 Oct 2022 19:25:46 +0200 Subject: [PATCH 01/29] initial commit of 3d approximation for CNS --- src/Trixi.jl | 2 +- .../compressible_navier_stokes_3d.jl | 420 +++++++++++ src/equations/equations_parabolic.jl | 1 + src/solvers/dgsem_tree/dg.jl | 1 + src/solvers/dgsem_tree/dg_3d_parabolic.jl | 667 ++++++++++++++++++ 5 files changed, 1090 insertions(+), 1 deletion(-) create mode 100644 src/equations/compressible_navier_stokes_3d.jl create mode 100644 src/solvers/dgsem_tree/dg_3d_parabolic.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 5811a950505..33138545e25 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -135,7 +135,7 @@ export AcousticPerturbationEquations2D, ShallowWaterEquations1D, ShallowWaterEquations2D export LaplaceDiffusion2D, - CompressibleNavierStokesDiffusion2D + CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D export GradientVariablesPrimitive, GradientVariablesEntropy diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl new file mode 100644 index 00000000000..e682fd1ca4c --- /dev/null +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -0,0 +1,420 @@ +@doc raw""" + CompressibleNavierStokesDiffusion3D(gamma, mu, Pr, equations, + gradient_variables=GradientVariablesPrimitive()) + +These equations contain the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective terms from +the [`CompressibleEulerEquations3D`](@ref). + +- `gamma`: adiabatic constant, +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `equations`: instance of the [`CompressibleEulerEquations3D`](@ref) +- `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 \mathbf{v} \\ \rho e +\end{pmatrix} ++ +\nabla \cdot +\begin{pmatrix} + \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e + p) \mathbf{v} +\end{pmatrix} += +\nabla \cdot +\begin{pmatrix} +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla 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_1^2+v_2^2+v_3^2) \right) +``` +as the pressure. The terms on the right hand side of the system above +are built from the viscous stress tensor +```math +\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} +``` +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} +``` +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 +\nabla 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. + +In two spatial dimensions we require gradients for three quantities, e.g., +primitive quantities +```math +\nabla v_1,\, \nabla v_2,\, \nabla v_3,\, \nabla T +``` +or the entropy variables +```math +\nabla w_2,\, \nabla w_3,\, \nabla w_4\, \nabla w_5 +``` +where +```math +w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} +``` + +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. +""" +struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3}} <: AbstractCompressibleNavierStokesDiffusion{3, 5} + # 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 # CompressibleEulerEquations3D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +# TODO: Delete. I do not think these need declared again because they are in the 2D file and are +# dimension agnostic. +# struct GradientVariablesPrimitive end +# struct GradientVariablesEntropy end + +# default to primitive gradient variables +function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquations3D; + 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 + + CompressibleNavierStokesDiffusion3D{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) , ::CompressibleNavierStokesDiffusion3D) = ("v1", "v2", "v3", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion3D) = ("w2", "w3", "w4", "w5") + +varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusion3D) = + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) + +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) = cons2prim +gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) = cons2entropy + + +# 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::CompressibleNavierStokesDiffusion3D) + # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. + rho, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T) + # either computed directly or reverse engineered from the gradient of the entropy vairables + # by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, dv3dx, dTdx = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, dv3dy, dTdy = convert_derivative_to_primitive(u, gradients[2], equations) + _, dv1dz, dv2dz, dv3dz, dTdz = convert_derivative_to_primitive(u, gradients[3], equations) + + # Components of viscous stress tensor + + # Diagonal parts + # (4/3 * (v1)_x - 2/3 * ((v2)_y + (v3)_z) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * (dv2dy + dv3dz) + # (4/3 * (v2)_y - 2/3 * ((v1)_x + (v3)_z) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * (dv1dx + dv3dz) + # (4/3 * (v3)_z - 2/3 * ((v1)_x + (v2)_y) + tau_33 = 4.0 / 3.0 * dv3dz - 2.0 / 3.0 * (dv1dx + dv2dy) + + # Off diagonal parts, exploit that stress tensor is symmetric + # ((v1)_y + (v2)_x) + tau_12 = dv1dy + dv2dx # = tau_21 + # ((v1)_z + (v3)_x) + tau_13 = dv1dz + dv3dx # = tau_31 + # ((v2)_z + (v3)_y) + tau_23 = dv2dz + dv3dy # = tau_32 + + # 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 + q2 = equations.kappa * dTdy + q3 = equations.kappa * dTdz + + # Constant dynamic viscosity is copied to a variable for readibility. + # 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 + + if orientation == 1 + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = tau_13 * mu + f5 = ( v1 * tau_11 + v2 * tau_12 + v3 * tau_13 + q1 ) * mu + + return SVector(f1, f2, f3, f4, f5) + elseif orientation == 2 + # viscous flux components in the y-direction + # Note, symmetry is exploited for tau_12 = tau_21 + g1 = zero(rho) + g2 = tau_12 * mu # tau_21 * mu + g3 = tau_22 * mu + g4 = tau_23 * mu + g5 = ( v1 * tau_12 + v2 * tau_22 + v3 * tau_23 + q2 ) * mu + + return SVector(g1, g2, g3, g4, g5) + else # if orientation == 3 + # viscous flux components in the z-direction + # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32 + h1 = zero(rho) + h2 = tau_13 * mu # tau_31 * mu + h3 = tau_23 * mu # tau_32 * mu + h4 = tau_33 * mu + h5 = ( v1 * tau_13 + v2 * tau_23 + v3 * tau_33 + q3 ) * mu + + return SVector(h1, h2, h3, h4, h5) + end +end + + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + T = temperature(u, equations) + + return SVector(rho, v1, v2, v3, 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 `CompressibleNavierStokesDiffusion2D`, +# but this may be confusing to new users. +cons2entropy(u, equations::CompressibleNavierStokesDiffusion3D) = cons2entropy(u, equations.equations_hyperbolic) +entropy2cons(w, equations::CompressibleNavierStokesDiffusion3D) = entropy2cons(w, equations.equations_hyperbolic) + +# 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::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + return u_transformed +end + +# TODO: parabolic. Make this more efficient! +@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + # note: this uses CompressibleNavierStokesDiffusion3D 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, ::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + return gradient +end + +# the first argument is always the "transformed" variables. +@inline function convert_derivative_to_primitive(w, gradient_entropy_vars, + equations::CompressibleNavierStokesDiffusion3D{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 CompressibleNavierStokesDiffusion3D + rho, rho_v1, rho_v2, rho_v3, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + T = temperature(u, equations) + + return SVector(gradient_entropy_vars[1], + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[5]), # grad(u) = T*(grad(w_2)+v1*grad(w_5)) + T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_3)+v2*grad(w_5)) + T * (gradient_entropy_vars[4] + v3 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_4)+v3*grad(w_5)) + T * T * gradient_entropy_vars[5] # grad(T) = T^2*grad(w_5)) + ) +end + + +# This routine is required because `prim2cons` is called in `initial_condition`, which +# is called with `equations::CompressibleEulerEquations3D`. This means it is inconsistent +# with `cons2prim(..., ::CompressibleNavierStokesDiffusion3D)` as defined above. +# TODO: parabolic. Is there a way to clean this up? +@inline prim2cons(u, equations::CompressibleNavierStokesDiffusion3D) = + prim2cons(u, equations.equations_hyperbolic) + + +@inline function temperature(u, equations::CompressibleNavierStokesDiffusion3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) + T = p / rho + return T +end + +# TODO: can we generalize this to MHD? +# TODO: These structs only need declared once. +# """ +# 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::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + return SVector(u_inner[1], v1, v2, v3, u_inner[5]) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + # rho, v1, v2, v3, _ = u_inner + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], normal_energy_flux) +end + + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + v1, v2, v3 = 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, v2, v3, T) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{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::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + negative_rho_inv_p = w_inner[5] # w_5 = -rho / p + return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p, -v3 * 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::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + v1, v2, v3 = 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 * w5. Similarly for w3 and w4 + w5 = -1 / T + return SVector(w_inner[1], -v1 * w5, -v2 * w5, -v3 * w5, w5) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], flux_inner[5]) +end diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 340cf0e4294..4fcdfff5f94 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -9,3 +9,4 @@ include("laplace_diffusion_2d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("compressible_navier_stokes_2d.jl") +include("compressible_navier_stokes_3d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 0ffaec2eced..e0066c2d2d8 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -64,6 +64,7 @@ include("dg_2d_parabolic.jl") # 3D DG implementation include("dg_3d.jl") +include("dg_3d_parabolic.jl") # Auxiliary functions that are specialized on this solver # as well as specialized implementations used to improve performance diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl new file mode 100644 index 00000000000..d3a47cb06be --- /dev/null +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -0,0 +1,667 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + +# This file collects all methods that have been updated to work with parabolic systems of equations +# +# assumptions: parabolic terms are of the form div(f(u, grad(u))) and +# will be discretized first order form as follows: +# 1. compute grad(u) +# 2. compute f(u, grad(u)) +# 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call) +# boundary conditions will be applied to both grad(u) and div(f(u, grad(u))). +function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + initial_condition, boundary_conditions_parabolic, source_terms, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @unpack u_transformed, gradients, flux_viscous = cache_parabolic + + # Convert conservative variables to a form more suitable for viscous flux calculations + @trixi_timeit timer() "transform variables" transform_variables!( + u_transformed, u, mesh, equations_parabolic, dg, parabolic_scheme, cache, cache_parabolic) + + # Compute the gradients of the transformed variables + @trixi_timeit timer() "calculate gradient" calc_gradient!( + gradients, u_transformed, t, mesh, equations_parabolic, boundary_conditions_parabolic, dg, + cache, cache_parabolic) + + # Compute and store the viscous fluxes + @trixi_timeit timer() "calculate viscous fluxes" calc_viscous_fluxes!( + flux_viscous, gradients, u_transformed, mesh, equations_parabolic, dg, cache, cache_parabolic) + + # The remainder of this function is essentially a regular rhs! for parabolic equations (i.e., it + # computes the divergence of the viscous fluxes) + # + # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have + # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # `interfaces` and `boundaries` containers in `cache_parabolic` to interpolate and store the + # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we + # do not need to recreate the existing data structure only with a different name, and c) we do not + # need to interpolate solutions *and* gradients to the surfaces. + + # TODO: parabolic; reconsider current data structure reuse strategy + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" calc_volume_integral!( + du, flux_viscous, mesh, equations_parabolic, dg, cache) + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" calc_interface_flux!( + cache_parabolic.elements.surface_flux_values, mesh, equations_parabolic, dg, cache_parabolic) + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" calc_boundary_flux_divergence!( + cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic, + dg.surface_integral, dg) + + # TODO: parabolic; extend to mortars + @assert nmortars(dg, cache) == 0 + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" calc_surface_integral!( + du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic) + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!( + du, mesh, equations_parabolic, dg, cache_parabolic) + + return nothing +end + +# Transform solution variables prior to taking the gradient +# (e.g., conservative to primitive variables). Defaults to doing nothing. +# TODO: can we avoid copying data? +function transform_variables!(u_transformed, u, mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations_parabolic, dg, i, j, k, element) + u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, equations_parabolic) + set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, k, element) + end + end +end + +# This is the version used when calculating the divergence of the viscous fluxes +function calc_volume_integral!(du, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + flux_1_node = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, k, element) + flux_2_node = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, k, element) + flux_3_node = get_node_vars(flux_viscous_z, equations_parabolic, dg, i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux_1_node, equations_parabolic, dg, ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[jj, j], flux_2_node, equations_parabolic, dg, i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[kk, k], flux_3_node, equations_parabolic, dg, i, j, kk, element) + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +# We pass the `surface_integral` argument solely for dispatch +function prolong2interfaces!(cache_parabolic, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache) + @unpack interfaces = cache_parabolic + @unpack orientations = interfaces + + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for interface in eachinterface(dg, cache) + left_element = interfaces.neighbor_ids[1, interface] + right_element = interfaces.neighbor_ids[2, interface] + + if orientations[interface] == 1 + # interface in x-direction + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j, k, left_element] + interfaces.u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k, right_element] + end + elseif orientations[interface] == 2 + # interface in y-direction + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg), k, left_element] + interfaces.u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k, right_element] + end + else # if orientations[interface] == 3 + # interface in z-direction + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, i, j, interface] = flux_viscous_z[v, i, j, nnodes(dg), left_element] + interfaces.u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1, right_element] + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +function calc_interface_flux!(surface_flux_values, + mesh::TreeMesh{3}, equations_parabolic, + dg::DG, cache_parabolic) + @unpack neighbor_ids, orientations = cache_parabolic.interfaces + + @threaded for interface in eachinterface(dg, cache_parabolic) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + # orientation = 3: left -> 6, right -> 5 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for j in eachnode(dg), i in eachnode(dg) + # Get precomputed fluxes at interfaces + flux_ll, flux_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, + dg, i, j, interface) + + # Compute interface flux as mean of left and right viscous fluxes + # TODO: parabolic; only BR1 at the moment + flux = 0.5 * (flux_ll + flux_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, left_direction, left_id] = flux[v] + surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +function prolong2boundaries!(cache_parabolic, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache) + @unpack boundaries = cache_parabolic + @unpack orientations, neighbor_sides = boundaries + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for boundary in eachboundary(dg, cache_parabolic) + element = boundaries.neighbor_ids[boundary] + + if orientations[boundary] == 1 + # boundary in x-direction + if neighbor_sides[boundary] == 1 + # element in -x direction of boundary + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg), j, k, element] + end + else # Element in +x direction of boundary + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k, element] + end + end + elseif orientations[boundary] == 2 + # boundary in y-direction + if neighbor_sides[boundary] == 1 + # element in -y direction of boundary + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, i, k, boundary] = flux_viscous_y[v, i, nnodes(dg), k, element] + end + else + # element in +y direction of boundary + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k, element] + end + end + else # if orientations[boundary] == 3 + # boundary in z-direction + if neighbor_sides[boundary] == 1 + # element in -z direction of boundary + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, i, j, boundary] = flux_viscous_z[v, i, j, nnodes(dg), element] + end + else + # element in +z direction of boundary + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1, element] + end + end + end + end + + return nothing +end + + +function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache, cache_parabolic) + gradients_x, gradients_y, gradients_z = gradients + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays + + @threaded for element in eachelement(dg, cache) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + # Get solution and gradients + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, element) + gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, k, element) + gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, k, element) + gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k, element) + + # Calculate viscous flux and store each component for later use + flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 1, equations_parabolic) + flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 2, equations_parabolic) + flux_viscous_node_z = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 3, equations_parabolic) + set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, i, j, k, element) + set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, i, j, k, element) + set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg, i, j, k, element) + end + end +end + + +# TODO: parabolic; decide if we should keep this. +function get_unsigned_normal_vector_3d(direction) + if direction > 6 || direction < 1 + error("Direction = $direction; in 3D, direction should be 1, 2, 3, 4, 5, or 6.") + end + if direction == 1 || direction == 2 + return SVector(1.0, 0.0, 0.0) + elseif direction == 3 || direction == 4 + return SVector(0.0, 1.0, 0.0) + else + return SVector(0.0, 0.0, 1.0) + end +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[5], + equations_parabolic, surface_integral, dg, cache, + 5, firsts[5], lasts[5]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[6], + equations_parabolic, surface_integral, dg, cache, + 6, firsts[6], lasts[6]) +end + + +function calc_boundary_flux_by_direction_gradient!(surface_flux_values::AbstractArray{<:Any,5}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, j, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + end + + # TODO: revisit if we want more general boundary treatments. + # This assumes the gradient numerical flux at the boundary is the gradient variable, + # which is consistent with BR1, LDG. + flux_inner = u_inner + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, j, boundary) + flux = boundary_condition(flux_inner, u_inner, get_unsigned_normal_vector_3d(direction), + x, t, Gradient(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[5], + equations_parabolic, surface_integral, dg, cache, + 5, firsts[5], lasts[5]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[6], + equations_parabolic, surface_integral, dg, cache, + 6, firsts[6], lasts[6]) +end +function calc_boundary_flux_by_direction_divergence!(surface_flux_values::AbstractArray{<:Any,5}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + + # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") + # of the viscous flux, as computed in `prolong2boundaries!` + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Get viscous boundary fluxes + flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, j, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + flux_inner = flux_ll + else # Element is on the right, boundary on the left + flux_inner = flux_rr + end + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, j, boundary) + + # TODO: add a field in `cache.boundaries` for gradient information. UPDATE THIS COMMENT + # Here, we pass in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information. + # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion3D and + # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesDiffusion3D as of 2022-6-27. + # It will not work with implementations which utilize `u_inner` to impose boundary conditions. + flux = boundary_condition(flux_inner, nothing, get_unsigned_normal_vector_3d(direction), + x, t, Divergence(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + + +# Calculate the gradient of the transformed variables +function calc_gradient!(gradients, u_transformed, t, + mesh::TreeMesh{3}, equations_parabolic, + boundary_conditions_parabolic, dg::DG, cache, cache_parabolic) + + gradients_x, gradients_y, gradients_z = gradients + + # Reset du + @trixi_timeit timer() "reset gradients" begin + reset_du!(gradients_x, dg, cache) + reset_du!(gradients_y, dg, cache) + reset_du!(gradients_z, dg, cache) + end + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + @unpack derivative_dhat = dg.basis + @threaded for element in eachelement(dg, cache) + + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], u_node, equations_parabolic, dg, ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], u_node, equations_parabolic, dg, i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], u_node, equations_parabolic, dg, i, j, kk, element) + end + end + end + end + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!( + cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" begin + @unpack surface_flux_values = cache_parabolic.elements + @unpack neighbor_ids, orientations = cache_parabolic.interfaces + + @threaded for interface in eachinterface(dg, cache_parabolic) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + # orientation = 3: left -> 6, right -> 5 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for j in eachnode(dg), i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, + equations_parabolic, dg, i, j, interface) + flux = 0.5 * (u_ll + u_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, left_direction, left_id] = flux[v] + surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + end + end + end + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!( + cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!( + cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic, + dg.surface_integral, dg) + + # TODO: parabolic; mortars + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for m in eachnode(dg), l in eachnode(dg) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients_x[v, 1, l, m, element] = ( + gradients_x[v, 1, l, m, element] - surface_flux_values[v, l, m, 1, element] * factor_1) + + # surface at +x + gradients_x[v, nnodes(dg), l, m, element] = ( + gradients_x[v, nnodes(dg), l, m, element] + surface_flux_values[v, l, m, 2, element] * factor_2) + + # surface at -y + gradients_y[v, l, 1, m, element] = ( + gradients_y[v, l, 1, m, element] - surface_flux_values[v, l, m, 3, element] * factor_1) + + # surface at +y + gradients_y[v, l, nnodes(dg), m, element] = ( + gradients_y[v, l, nnodes(dg), m, element] + surface_flux_values[v, l, m, 4, element] * factor_2) + + # surface at -z + gradients_z[v, l, m, 1, element] = ( + gradients_z[v, l, m, 1, element] - surface_flux_values[v, l, m, 5, element] * factor_1) + + # surface at +z + gradients_z[v, l, m, nnodes(dg), element] = ( + gradients_z[v, l, m, nnodes(dg), element] + surface_flux_values[v, l, m, 6, element] * factor_2) + end + end + end + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" begin + apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic) + end + + return nothing +end + + +# This method is called when a SemidiscretizationHyperbolic is constructed. +# It constructs the basic `cache` used throughout the simulation to compute +# the RHS etc. +function create_cache_parabolic(mesh::TreeMesh{3}, equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, RealT, uEltype) + # Get cells for which an element needs to be created (i.e. all leaf cells) + leaf_cell_ids = local_leaf_cells(mesh.tree) + + elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) + + n_vars = nvariables(equations_hyperbolic) + n_nodes = nnodes(elements) + n_elements = nelements(elements) + u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) + gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) + flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) + + interfaces = init_interfaces(leaf_cell_ids, mesh, elements) + + boundaries = init_boundaries(leaf_cell_ids, mesh, elements) + + # mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar) + + # cache = (; elements, interfaces, boundaries, mortars) + cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + + # Add specialized parts of the cache required to compute the mortars etc. + # cache = (;cache..., create_cache(mesh, equations_parabolic, dg.mortar, uEltype)...) + + return cache +end + + +# Needed to *not* flip the sign of the inverse Jacobian. +# This is because the parabolic fluxes are assumed to be of the form +# `du/dt + df/dx = dg/dx + source(x,t)`, +# where f(u) is the inviscid flux and g(u) is the viscous flux. +function apply_jacobian!(du, mesh::TreeMesh{3}, + equations::AbstractEquationsParabolic, dg::DG, cache) + + @threaded for element in eachelement(dg, cache) + factor = cache.elements.inverse_jacobian[element] + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + du[v, i, j, k, element] *= factor + end + end + end + + return nothing +end + +end # @muladd From b5e13545f7b1d76e0630403e283f721b32277316 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 18 Oct 2022 19:26:05 +0200 Subject: [PATCH 02/29] typo fix in CNS 2D file --- src/equations/compressible_navier_stokes_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 9f7758aa476..1234a02b4a8 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -142,7 +142,7 @@ gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientV # Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 -# of the paper by Rueda-Ramíreza, Hennemann, Hindenlang, Winters, and Gassner +# 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. From ab9329c6bd98ec7d9682a245cb1c451174bd47e3 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 18 Oct 2022 19:26:41 +0200 Subject: [PATCH 03/29] add TGV elixir for compressible Navier-Stokes --- ...elixir_navierstokes_taylor_green_vortex.jl | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl new file mode 100644 index 00000000000..e660fcaa951 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -0,0 +1,78 @@ + +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 = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), + Prandtl=prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + +The classical inviscid Taylor-Green vortex. +""" +function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3]) + v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3]) + v3 = 0.0 + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3])) + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = ( 1.0, 1.0, 1.0) .* pi +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=100_000) + + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary From 82369ec0e65e9939968959cb35ca7a413a4175ab Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 18 Oct 2022 19:58:57 -0500 Subject: [PATCH 04/29] add DGMulti CNS 3D elixir --- ...elixir_navierstokes_taylor_green_vortex.jl | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl new file mode 100644 index 00000000000..c39e0438849 --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -0,0 +1,69 @@ + +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 = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), + Prandtl=prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + +The classical inviscid Taylor-Green vortex. +""" +function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3]) + v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3]) + v3 = 0.0 + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3])) + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = GaussSBP(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = ( 1.0, 1.0, 1.0) .* pi +mesh = DGMultiMesh(dg; coordinates_min, coordinates_max, + cells_per_dimension=(8, 8, 8), + periodicity=(true, true, true)) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, dg) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary From 5929572717c3db108659c057f5e83492b95fca0f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 16:14:48 +0200 Subject: [PATCH 05/29] need TreeMesh versions of slip wall BC in 3D compressible Euler for convergence testing --- src/equations/compressible_euler_3d.jl | 46 ++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 1b93674c8a2..eafaaff7376 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -323,6 +323,52 @@ Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of zero(eltype(u_inner))) * norm_ end +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations3D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(1.0, 0.0, 0.0) + elseif orientation == 2 + normal_direction = SVector(0.0, 1.0, 0.0) + else # orientation == 3 + normal_direction = SVector(0.0, 0.0, 1.0) + end + + # compute and return the flux using `boundary_condition_slip_wall` routine above + return boundary_condition_slip_wall(u_inner, normal_direction, direction, + x, t, surface_flux_function, equations) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations3D) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x, -y, and -z sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, equations) + end + + return boundary_flux +end # Calculate 1D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations3D) From 9267369b9ac99ef03ac52ace8cba7eaa16f84419 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 16:15:18 +0200 Subject: [PATCH 06/29] update comments in CNS 3D file --- .../compressible_navier_stokes_3d.jl | 61 ++----------------- 1 file changed, 4 insertions(+), 57 deletions(-) diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index e682fd1ca4c..0cbf24d3f6d 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -116,7 +116,7 @@ end # TODO: parabolic # This is the flexibility a user should have to select the different gradient variable types -# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion3D) = ("v1", "v2", "v3", "T") +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion3D) = ("v1", "v2", "v3", "T") # varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion3D) = ("w2", "w3", "w4", "w5") varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusion3D) = @@ -240,8 +240,8 @@ end 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). +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4, w_5) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, v3, 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 @@ -254,7 +254,7 @@ end equations::CompressibleNavierStokesDiffusion3D{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 + # We can fix this if we directly compute v1, v2, v3, T from the entropy variables u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion3D rho, rho_v1, rho_v2, rho_v3, _ = u @@ -288,59 +288,6 @@ end return T end -# TODO: can we generalize this to MHD? -# TODO: These structs only need declared once. -# """ -# 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, From 01b3752d8b551c133794651f6efa4ffc20afcbbb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 16:16:18 +0200 Subject: [PATCH 07/29] add convergence test elixir for 3D compressible Navier-Stokes --- .../elixir_navierstokes_convergence.jl | 269 ++++++++++++++++++ 1 file changed, 269 insertions(+) create mode 100644 examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..d5c9951d0e4 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,269 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + # gradient_variables=GradientVariablesPrimitive()) + gradient_variables=GradientVariablesEntropy()) +# 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, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = ( 1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + periodicity=(true, false, true), + n_cells_max=100_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # TODO: cleanup constants + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 # sqrt(2.0 / 3.0) + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # 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() + + # TODO: cleanup constants + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 # sqrt(2.0 / 3.0) + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +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]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_top_bottom, + y_pos = boundary_condition_top_bottom, + 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), + 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-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary + From 2921f1d0e7f0158afff8c7c376aa98eace401571 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 16:17:04 +0200 Subject: [PATCH 08/29] remove debug statements from new elixir --- examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index d5c9951d0e4..8abb3a5e04d 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -9,8 +9,8 @@ mu() = 0.01 equations = CompressibleEulerEquations3D(1.4) equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), - # gradient_variables=GradientVariablesPrimitive()) - gradient_variables=GradientVariablesEntropy()) + 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()) From bfb7114b049998f89d2d3486591e7f18c8f6991b Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 16:17:21 +0200 Subject: [PATCH 09/29] adjust maximum tree size in new elixir --- examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index 8abb3a5e04d..f6cf23b88fc 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -22,7 +22,7 @@ coordinates_max = ( 1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max( mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=3, periodicity=(true, false, true), - n_cells_max=100_000) # set maximum capacity of tree data structure + n_cells_max=50_000) # set maximum capacity of tree data structure # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) From bfccb0746f39b53281ebfa1b1e3254108ff0eccb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 19:15:38 +0200 Subject: [PATCH 10/29] add convergence test for DGMulti 3D --- .../elixir_navierstokes_convergence.jl | 255 ++++++++++++++++++ .../elixir_navierstokes_convergence.jl | 4 +- 2 files changed, 257 insertions(+), 2 deletions(-) create mode 100644 examples/dgmulti_3d/elixir_navierstokes_convergence.jl diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..e128cb1f08e --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl @@ -0,0 +1,255 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(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 +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) +mesh = DGMultiMesh(dg, cells_per_dimension=(8, 8, 8); periodicity=(true, false, true), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # TODO: cleanup constants + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # 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() + + # TODO: cleanup constants + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +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]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + 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, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index f6cf23b88fc..a7dbd73886e 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -33,7 +33,7 @@ function initial_condition_navier_stokes_convergence_test(x, t, equations) # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 - A2 = 1.0 # sqrt(2.0 / 3.0) + A2 = 1.0 A3 = 0.5 # Convenience values for trig. functions @@ -63,7 +63,7 @@ end # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 - A2 = 1.0 # sqrt(2.0 / 3.0) + A2 = 1.0 A3 = 0.5 # Convenience values for trig. functions From 323e4cfc763803814f5b6c93226caa8cea51faef Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 19:45:13 +0200 Subject: [PATCH 11/29] add battery of 3D parabolic tests. All pass locally --- test/runtests.jl | 1 + test/test_parabolic_3d.jl | 88 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+) create mode 100644 test/test_parabolic_3d.jl diff --git a/test/runtests.jl b/test/runtests.jl index 79a8e54e791..d16ec13256b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -88,6 +88,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" include("test_parabolic_2d.jl") + include("test_parabolic_3d.jl") end @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl new file mode 100644 index 00000000000..5213685b6d7 --- /dev/null +++ b/test/test_parabolic_3d.jl @@ -0,0 +1,88 @@ +module TestExamplesParabolic3D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_3d") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +@testset "SemidiscretizationHyperbolicParabolic" begin + + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), + l2 = [0.0005532847115849239, 0.000659263490965341, 0.0007776436127362806, 0.0006592634909662951, 0.0038073628897809185], + linf = [0.0017039861523615585, 0.002628561703560073, 0.003531057425112172, 0.0026285617036090336, 0.015587829540351095] + ) + end + + @trixi_testset "DGMulti: elixir_navierstokes_taylor_green_vortex.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_taylor_green_vortex.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.25), + l2 = [0.0001825713444029892, 0.015589736382772248, 0.015589736382771884, 0.021943924667273653, 0.01927370280244222], + linf = [0.0006268463584697681, 0.03218881662749007, 0.03218881662697948, 0.053872495395614256, 0.05183822000984151] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.0019582188528512257, 0.002653449504302844, 0.002898264205184629, 0.002653449504302853, 0.009511572365085706], + linf = [0.013680656759085918, 0.0356910450154318, 0.023526343547736236, 0.035691045015431855, 0.11482570604041165] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.00195468651965362, 0.0026554367897028506, 0.002892730402724066, 0.002655436789702817, 0.009596351796609566], + linf = [0.013680508110645473, 0.035673446359424356, 0.024024936779729028, 0.03567344635942474, 0.11839497110809383] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (Entropy gradient variables)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + l2 = [0.0019770444875099307, 0.0026524750946399327, 0.00290860030832445, 0.0026524750946399396, 0.009509568981439294], + linf = [0.01387936112914212, 0.03526260609304053, 0.023554197097368997, 0.035262606093040896, 0.11719963716509518] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.001974631423398113, 0.002654768259143932, 0.002907031063651286, 0.002654768259143901, 0.009587792882971452], + linf = [0.01387919380137137, 0.035244084526358944, 0.02398614622061363, 0.03524408452635828, 0.12005056512506407] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (flux differencing)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + volume_integral=VolumeIntegralFluxDifferencing(flux_central), + l2 = [0.0019582188528180213, 0.002653449504301736, 0.0028982642051960006, 0.0026534495043017384, 0.009511572364811033], + linf = [0.013680656758949583, 0.035691045015224444, 0.02352634354676752, 0.035691045015223424, 0.11482570603751441] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.25), + l2 = [0.00024694280871739453, 0.01572456212714984, 0.01572456212714988, 0.02202557418601852, 0.030617743072153637], + linf = [0.0008556454642751898, 0.04269724800380153, 0.04269724800380213, 0.07517356785606294, 0.16689008447096398] + ) + end + +end + +# Clean up afterwards: delete Trixi output directory +@test_nowarn isdir(outdir) && rm(outdir, recursive=true) + +end # module From 1929552942cabb3979890f103821c4c5c5f0fc4a Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 19 Oct 2022 19:54:46 +0200 Subject: [PATCH 12/29] remove old comments and TODOs --- examples/dgmulti_3d/elixir_navierstokes_convergence.jl | 2 -- examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl | 2 -- src/equations/compressible_navier_stokes_3d.jl | 5 ----- 3 files changed, 9 deletions(-) diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl index e128cb1f08e..b08eda5cf8c 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl @@ -25,7 +25,6 @@ mesh = DGMultiMesh(dg, cells_per_dimension=(8, 8, 8); periodicity=(true, false, # and by the initial condition (which passes in `CompressibleEulerEquations3D`). # This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) function initial_condition_navier_stokes_convergence_test(x, t, equations) - # TODO: cleanup constants # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 @@ -55,7 +54,6 @@ end Pr = prandtl_number() mu_ = mu() - # TODO: cleanup constants # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index a7dbd73886e..9199174652a 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -29,7 +29,6 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # and by the initial condition (which passes in `CompressibleEulerEquations3D`). # This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) function initial_condition_navier_stokes_convergence_test(x, t, equations) - # TODO: cleanup constants # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 @@ -59,7 +58,6 @@ end Pr = prandtl_number() mu_ = mu() - # TODO: cleanup constants # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` c = 2.0 A1 = 0.5 diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 0cbf24d3f6d..25ef2780b75 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -91,11 +91,6 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E < gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy end -# TODO: Delete. I do not think these need declared again because they are in the 2D file and are -# dimension agnostic. -# struct GradientVariablesPrimitive end -# struct GradientVariablesEntropy end - # default to primitive gradient variables function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquations3D; mu, Prandtl, From 559f47c86f34202b07b3f2ba732ea24d4879d822 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 20 Oct 2022 06:50:54 +0200 Subject: [PATCH 13/29] bug fix in divergence B analysis routine --- src/callbacks_step/analysis_dg3d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 5267960e2f0..e7737169b07 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -220,7 +220,7 @@ function analyze(::Val{:l2_divb}, du, u, t, for l in eachnode(dg) divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] + derivative_matrix[j, l] * u[7, i, l, k, element] + - derivative_matrix[k, l] * u[7, i, j, l, element] ) + derivative_matrix[k, l] * u[8, i, j, l, element] ) end divb *= cache.elements.inverse_jacobian[element] divb^2 @@ -262,7 +262,7 @@ function analyze(::Val{:linf_divb}, du, u, t, for l in eachnode(dg) divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] + derivative_matrix[j, l] * u[7, i, l, k, element] + - derivative_matrix[k, l] * u[7, i, j, l, element] ) + derivative_matrix[k, l] * u[8, i, j, l, element] ) end divb *= cache.elements.inverse_jacobian[element] linf_divb = max(linf_divb, abs(divb)) From 2f051402485547374326531765711f3f8d354654 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 1 Nov 2022 14:54:16 +0100 Subject: [PATCH 14/29] add vorticity and enstrophy functions --- src/Trixi.jl | 3 ++- src/equations/compressible_navier_stokes_3d.jl | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 33138545e25..15687e6ab62 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -176,7 +176,8 @@ export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_c export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure -export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity +export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, + enstrophy export lake_at_rest_error export ncomponents, eachcomponent diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 25ef2780b75..762bf2b9987 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -284,6 +284,24 @@ end end +@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D) + # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v + + omega = vorticity(u, gradients, equations) + return 0.5 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2) +end + + +@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D) + # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, dv3dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, dv3dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) + _, dv1dz, dv2dz, dv3dz, _ = convert_derivative_to_primitive(u, gradients[3], equations) + + return SVector(dv3dy - dv2dz , dv1dz - dv3dx , dv2dx - dv1dy) +end + + @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) From 54ff12a1673df9856b26edb5fed4a0917caa410d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 1 Nov 2022 14:54:52 +0100 Subject: [PATCH 15/29] add first version of an enstrophy callback --- src/callbacks_step/analysis.jl | 18 ++++++++++++++++++ src/callbacks_step/analysis_dg3d.jl | 16 ++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 59b2d2e9113..62aeffe2fcb 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -470,6 +470,21 @@ pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) +# Special analyze for `SemidiscretizationHyperbolicParabolic` such that +# precomputed gradients are available. For now only implemented for the `enstrophy` +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. +function analyze(quantity::typeof(enstrophy), du, u, t, semi::SemidiscretizationHyperbolicParabolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + cache_parabolic = semi.cache_parabolic + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, cache_parabolic) +end +function analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, cache_parabolic) + integrate(quantity, u, mesh, equations, equations_parabolic, solver, cache, cache_parabolic, normalize=true) +end + + function entropy_timederivative end pretty_form_utf(::typeof(entropy_timederivative)) = "∑∂S/∂U ⋅ Uₜ" pretty_form_ascii(::typeof(entropy_timederivative)) = "dsdu_ut" @@ -494,6 +509,9 @@ pretty_form_ascii(::typeof(energy_magnetic)) = "e_magnetic" pretty_form_utf(::typeof(cross_helicity)) = "∑v⋅B" pretty_form_ascii(::typeof(cross_helicity)) = "v_dot_B" +pretty_form_utf(::typeof(enstrophy)) = "∑enstrophy" +pretty_form_ascii(::typeof(enstrophy)) = "enstrophy" + pretty_form_utf(::Val{:l2_divb}) = "L2 ∇⋅B" pretty_form_ascii(::Val{:l2_divb}) = "l2_divb" diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index e7737169b07..f3a5b564c38 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -199,6 +199,22 @@ function integrate(func::Func, u, end +function integrate(func::typeof(enstrophy), u, + mesh::TreeMesh{3}, + equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, + dg::DGSEM, + cache, cache_parabolic; normalize=true) + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, k, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, k, element) + gradients_1_local = get_node_vars(gradients_x, equations_parabolic, dg, i, j, k, element) + gradients_2_local = get_node_vars(gradients_y, equations_parabolic, dg, i, j, k, element) + gradients_3_local = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k, element) + return func(u_local, (gradients_1_local, gradients_2_local, gradients_3_local), equations_parabolic) + end +end + + function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, equations, dg::DG, cache) From 91bb5a100978a14ae622f9284db29d078761d1e7 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 1 Nov 2022 14:55:23 +0100 Subject: [PATCH 16/29] add enstrophy callback to TGV elixir --- ...elixir_navierstokes_taylor_green_vortex.jl | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index e660fcaa951..d3f1bcbd3bd 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -34,7 +34,7 @@ end initial_condition = initial_condition_taylor_green_vortex volume_flux = flux_ranocha -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (-1.0, -1.0, -1.0) .* pi @@ -50,24 +50,22 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 10.0) +tspan = (0.0, 20.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)) -alive_callback = AliveCallback(analysis_interval=analysis_interval) - -save_solution = SaveSolutionCallback(interval=100, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2prim) +alive_callback = AliveCallback(analysis_interval=analysis_interval,) callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, - save_solution) + analysis_callback, + alive_callback) ############################################################################### # run the simulation From 104b8c8a700b5f02a7f5b7a1ba4438f707080b1d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 2 Nov 2022 08:30:18 +0100 Subject: [PATCH 17/29] update TGV test values since switched surface flux to flux_hllc --- test/test_parabolic_3d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 5213685b6d7..67c55dde6e9 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -75,8 +75,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "TreeMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"), initial_refinement_level = 2, tspan=(0.0, 0.25), - l2 = [0.00024694280871739453, 0.01572456212714984, 0.01572456212714988, 0.02202557418601852, 0.030617743072153637], - linf = [0.0008556454642751898, 0.04269724800380153, 0.04269724800380213, 0.07517356785606294, 0.16689008447096398] + l2 = [0.00024173250389635442, 0.015684268393762454, 0.01568426839376248, 0.021991909545192333, 0.02825413672911425], + linf = [0.0008410587892853094, 0.04740176181772552, 0.04740176181772507, 0.07483494924031157, 0.150181591534448] ) end From 841670528aa5dfcd549ee0c7d1798d63a0d17506 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 3 Nov 2022 11:46:54 -0500 Subject: [PATCH 18/29] adding draft DGMulti analysis callback --- src/callbacks_step/analysis_dgmulti.jl | 36 ++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index 3f1ea51e01c..65adc14137b 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -135,6 +135,42 @@ function analyze(::Val{:linf_divb}, du, u, t, return linf_divB end +function integrate(func::typeof(enstrophy), u, + mesh::DGMultiMesh, + equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, + dg::DGMulti, + cache, cache_parabolic; normalize=true) + + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + + # allocate local storage for gradients. + # TODO: can we avoid allocating here? + local_gradient_quadrature_values = ntuple(_ -> similar(cache_parabolic.local_u_values_threaded), 3) + + integral = zero(eltype(u)) + @threaded for e in eachelement(mesh, dg) + u_quadrature_values = cache_parabolic.local_u_values_threaded[Threads.threadid()] + gradient_x_quadrature_values = local_gradient_quadrature_values[1][Threads.threadid()] + gradient_y_quadrature_values = local_gradient_quadrature_values[2][Threads.threadid()] + gradient_z_quadrature_values = local_gradient_quadrature_values[3][Threads.threadid()] + + # interpolate to quadrature on each element + apply_to_each_field(mul_by(dg.basis.Vq), u_quadrature_values, view(u, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e)) + + # integrate over the element + for i in eachindex(u_quadrature_values) + gradients_i = SVector(gradient_x_quadrature_values[i], + gradient_y_quadrature_values[i], + gradient_z_quadrature_values[i]) + integral += mesh.md.wJq[i, e] * func(u_quadrature_values[i], gradients_i, equations) + end + end +end + + function create_cache_analysis(analyzer, mesh::DGMultiMesh, equations, dg::DGMulti, cache, RealT, uEltype) From 28c3eca81bfee3f0ad5c0e6fa73fe692f16c17fc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 3 Nov 2022 11:47:34 -0500 Subject: [PATCH 19/29] adding draft DGMulti analysis callback forgot the return statement --- src/callbacks_step/analysis_dgmulti.jl | 37 ++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index 3f1ea51e01c..f5ac3e4459d 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -135,6 +135,43 @@ function analyze(::Val{:linf_divb}, du, u, t, return linf_divB end +function integrate(func::typeof(enstrophy), u, + mesh::DGMultiMesh, + equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, + dg::DGMulti, + cache, cache_parabolic; normalize=true) + + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + + # allocate local storage for gradients. + # TODO: can we avoid allocating here? + local_gradient_quadrature_values = ntuple(_ -> similar(cache_parabolic.local_u_values_threaded), 3) + + integral = zero(eltype(u)) + @threaded for e in eachelement(mesh, dg) + u_quadrature_values = cache_parabolic.local_u_values_threaded[Threads.threadid()] + gradient_x_quadrature_values = local_gradient_quadrature_values[1][Threads.threadid()] + gradient_y_quadrature_values = local_gradient_quadrature_values[2][Threads.threadid()] + gradient_z_quadrature_values = local_gradient_quadrature_values[3][Threads.threadid()] + + # interpolate to quadrature on each element + apply_to_each_field(mul_by(dg.basis.Vq), u_quadrature_values, view(u, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e)) + + # integrate over the element + for i in eachindex(u_quadrature_values) + gradients_i = SVector(gradient_x_quadrature_values[i], + gradient_y_quadrature_values[i], + gradient_z_quadrature_values[i]) + integral += mesh.md.wJq[i, e] * func(u_quadrature_values[i], gradients_i, equations) + end + end + return integral +end + + function create_cache_analysis(analyzer, mesh::DGMultiMesh, equations, dg::DGMulti, cache, RealT, uEltype) From 5f9e40b08586a57402707a0c84ae1f99f57943bd Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 14 Nov 2022 10:55:57 +0100 Subject: [PATCH 20/29] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/equations/compressible_euler_3d.jl | 3 ++- src/equations/compressible_navier_stokes_3d.jl | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index eafaaff7376..64fbcad8a05 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -351,7 +351,8 @@ end boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, surface_flux_function, equations::CompressibleEulerEquations3D) -Should be used together with [`StructuredMesh`](@ref). +Should be used together with [`StructuredMesh`](@ref) and other possible curved +mesh types. """ @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, direction, x, t, diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 762bf2b9987..e3f9ca48ca5 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3}} <: AbstractCompressibleNavierStokesDiffusion{3, 5} # TODO: parabolic From c512038cb7ddcebd72c9439a828068fad1f5942e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 14 Nov 2022 10:58:36 +0100 Subject: [PATCH 21/29] make sure Experimental code text renders properly in docs --- src/equations/compressible_navier_stokes_2d.jl | 10 +++++----- src/equations/compressible_navier_stokes_3d.jl | 2 +- 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 1234a02b4a8..a5f0b132d81 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4} # TODO: parabolic @@ -92,9 +92,6 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E < 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 @@ -105,6 +102,9 @@ formulation from [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. + +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct GradientVariablesPrimitive end struct GradientVariablesEntropy end diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index e3f9ca48ca5..c3ee910c37c 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -66,7 +66,7 @@ primitive quantities ``` or the entropy variables ```math -\nabla w_2,\, \nabla w_3,\, \nabla w_4\, \nabla w_5 +\nabla w_2,\, \nabla w_3,\, \nabla w_4,\, \nabla w_5 ``` where ```math From 8c2e24255c84ebc21c5de80dce096f67793205f9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 14 Nov 2022 10:05:39 -0600 Subject: [PATCH 22/29] adding extra analysis integrals --- .../dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl index c39e0438849..b289762ca48 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -57,7 +57,10 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() alive_callback = AliveCallback(alive_interval=10) analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg), + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)) callbacks = CallbackSet(summary_callback, alive_callback) ############################################################################### From 8854647a9fc31ee9cf2ad4653b742c7e3a960164 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 16 Nov 2022 21:20:51 +0100 Subject: [PATCH 23/29] remove unused file path from parabolic test scripts --- test/test_parabolic_2d.jl | 2 +- test/test_parabolic_3d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index d9c3967315a..886a3740411 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -5,7 +5,7 @@ using Trixi include("test_trixi.jl") -EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") +# EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") # Start with a clean environment: remove Trixi output directory if it exists outdir = "out" diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 67c55dde6e9..ec3cb3f1dc7 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -5,7 +5,7 @@ using Trixi include("test_trixi.jl") -EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_3d") +# EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_3d") # Start with a clean environment: remove Trixi output directory if it exists outdir = "out" From 9575600f9d32bf0891d6ea631397250793416588 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 22 Nov 2022 09:30:38 -0600 Subject: [PATCH 24/29] remove threading for reduction --- src/callbacks_step/analysis_dgmulti.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index f5ac3e4459d..e2c6466398f 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -148,7 +148,7 @@ function integrate(func::typeof(enstrophy), u, local_gradient_quadrature_values = ntuple(_ -> similar(cache_parabolic.local_u_values_threaded), 3) integral = zero(eltype(u)) - @threaded for e in eachelement(mesh, dg) + for e in eachelement(mesh, dg) u_quadrature_values = cache_parabolic.local_u_values_threaded[Threads.threadid()] gradient_x_quadrature_values = local_gradient_quadrature_values[1][Threads.threadid()] gradient_y_quadrature_values = local_gradient_quadrature_values[2][Threads.threadid()] From 9f2f5e25347c1459b8faf51942b88451cf4c4578 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 29 Nov 2022 18:57:00 +0100 Subject: [PATCH 25/29] Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper --- test/test_parabolic_2d.jl | 1 - test/test_parabolic_3d.jl | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 886a3740411..83a6a0dbea7 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -5,7 +5,6 @@ using Trixi include("test_trixi.jl") -# EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") # Start with a clean environment: remove Trixi output directory if it exists outdir = "out" diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index ec3cb3f1dc7..c0728bfbb8c 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -5,13 +5,11 @@ using Trixi include("test_trixi.jl") -# EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_3d") - # Start with a clean environment: remove Trixi output directory if it exists outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -@testset "SemidiscretizationHyperbolicParabolic" begin +@testset "SemidiscretizationHyperbolicParabolic (3D)" begin @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence.jl"), From 29c5aa247e19b7b2597ed0b0a0b8a8b2f6ee0f84 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 29 Nov 2022 18:58:18 +0100 Subject: [PATCH 26/29] add designation to the 2D test set --- test/test_parabolic_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 886a3740411..338e4a0aa26 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -11,7 +11,7 @@ include("test_trixi.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -@testset "SemidiscretizationHyperbolicParabolic" begin +@testset "SemidiscretizationHyperbolicParabolic (2D)" begin @trixi_testset "DGMulti 2D rhs_parabolic!" begin From d8ecc1a383c36cd621c148b0e3e515b579922db0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 30 Jan 2023 07:20:44 +0100 Subject: [PATCH 27/29] remove specific types in new integrate routine --- src/callbacks_step/analysis_dg3d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index b23158b4d31..77cf1f819ea 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -199,11 +199,11 @@ function integrate(func::Func, u, end -function integrate(func::typeof(enstrophy), u, +function integrate(func::Func, u, mesh::TreeMesh{3}, - equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, + equations, equations_parabolic, dg::DGSEM, - cache, cache_parabolic; normalize=true) + cache, cache_parabolic; normalize=true) where {Func} gradients_x, gradients_y, gradients_z = cache_parabolic.gradients integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, k, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, k, element) From f159f69fb7806ba341c5b74c4d2189b8af800ac1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 30 Jan 2023 09:31:28 +0100 Subject: [PATCH 28/29] update docstrings for CNS to match the constructor --- src/equations/compressible_navier_stokes_2d.jl | 10 +++++----- src/equations/compressible_navier_stokes_3d.jl | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 1234a02b4a8..c509015ba7b 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,15 +1,14 @@ @doc raw""" - CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, equations, + CompressibleNavierStokesDiffusion2D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) -These equations contain the diffusion (i.e. parabolic) terms applied +Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations2D`](@ref). -- `gamma`: adiabatic constant, +- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `mu`: dynamic viscosity, - `Pr`: Prandtl number, -- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. @@ -37,7 +36,8 @@ where the system is closed with the ideal gas assumption giving ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` -as the pressure. The terms on the right hand side of the system above +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations2D`](@ref). +The terms on the right hand side of the system above are built from the viscous stress tensor ```math \underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 762bf2b9987..d4821c55e58 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -1,15 +1,14 @@ @doc raw""" - CompressibleNavierStokesDiffusion3D(gamma, mu, Pr, equations, + CompressibleNavierStokesDiffusion3D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) -These equations contain the diffusion (i.e. parabolic) terms applied +Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations3D`](@ref). -- `gamma`: adiabatic constant, +- `equations`: instance of the [`CompressibleEulerEquations3D`](@ref) - `mu`: dynamic viscosity, - `Pr`: Prandtl number, -- `equations`: instance of the [`CompressibleEulerEquations3D`](@ref) - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. @@ -37,7 +36,8 @@ where the system is closed with the ideal gas assumption giving ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right) ``` -as the pressure. The terms on the right hand side of the system above +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations2D`](@ref). +The terms on the right hand side of the system above are built from the viscous stress tensor ```math \underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} From d00329f98c666544a3e0310bb671e79845e2f706 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 2 Feb 2023 07:36:17 +0100 Subject: [PATCH 29/29] Update examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl --- examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl index b289762ca48..9d0d926338a 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -33,7 +33,6 @@ function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEule end initial_condition = initial_condition_taylor_green_vortex -volume_flux = flux_ranocha # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = GaussSBP(), surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),