From 60cce54c37a1adfae896acb89ba3c254f78f5580 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 7 Jun 2022 21:02:19 +0200 Subject: [PATCH 01/71] first draft of container for Navier-Stokes constants and fluxes --- .../compressible_navier_stokes_2d.jl | 228 ++++++++++++++++++ src/equations/equations_parabolic.jl | 4 + 2 files changed, 232 insertions(+) create mode 100644 src/equations/compressible_navier_stokes_2d.jl diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl new file mode 100644 index 00000000000..cf8a8b17fc7 --- /dev/null +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -0,0 +1,228 @@ +@doc raw""" + CompressibleNaiverStokes2D(gamma, + Re, + Pr, + Ma_inf, + kappa, + equations) + +`CompressibleNaiverStokes2D` represents the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective from +the `CompressibleEulerEquations2D`. + +gamma: adiabatic constant, +Re: Reynolds number, +Pr: Prandtl number, +Ma_inf: free-stream Mach number +kappa: thermal diffusivity for Fick's law + +For the particular scaling the vicosity is set internally to be μ = 1/Re. +Further, the nondimensionalization takes the density-temperature-sound speed as +the principle quantities such that +rho_inf = 1.0 +T_ref = 1.0 +c_inf = 1.0 +p_inf = 1.0 / gamma +u_inf = Ma_inf +R = 1.0 / gamma + +Other normalization strategies exist, see the reference below for details. +- Marc Montagnac (2013) + Variable Normalization (nondimensionalization and scaling) for Navier-Stokes + equations: a practical guide + [CERFACS Technical report](https://www.cerfacs.fr/~montagna/TR-CFD-13-77.pdf) +The scaling used herein is Section 4.5 of the reference. + +In two spatial dimensions we require gradients for three quantities, e.g., +primitive quantities + grad(u), grad(v), and grad(T) +or the entropy variables + grad(w_2), grad(w_3), grad(w_4) +where + w_2 = rho v_1 / p, w_3 = rho v_2 / p, w_4 = -rho / p +""" +# TODO: +# 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from +# the Euler equations +# 2) Add more here and probably some equations +struct CompressibleNaiverStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + Re::RealT # Reynolds number + Pr::RealT # Prandtl number + Ma_inf::RealT # free-stream Mach number + kappa::RealT # thermal diffusivity for Fick's law + + p_inf::RealT # free-stream pressure + u_inf::RealT # free-stream velocity + R::RealT # gas constant (depends on nondimensional scaling!) + + equations::E # CompressibleEulerEquations2D +end + +function CompressibleNaiverStokes2D(gamma, Reynolds, Prandtl, Mach_frestream, kappa, equations) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + + # From the nondimensionalization discussed above set the remaining free-stream + # quantities + p_inf = 1.0 / γ + u_inf = Mach_frestream + R = 1.0 / γ + CompressibleNaiverStokes2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, + Reynolds, Prandtl, Mach_freestream, + kappa, p_inf, u_inf, R, + equations) +end + +# I was not sure what to do here to allow flexibility of selecting primitive or entropy +# grandient variables +varnames(variable_mapping, equations_parabolic::CompressibleNaiverStokes2D) = + varnames(variable_mapping, equations_parabolic.equations) + + +# no orientation specified since the flux is vector-valued +# This form of the difussive Navier-Stokes fluxes is taken from Section 2 +# of the paper by Svärd, Carpenter and Nordström +# "A stable high-order finite difference scheme for the compressible Navier–Stokes +# equations, far-field boundary conditions" +# Note, could be generalized to use Sutherland's law to get the molecular and thermal +# diffusivity +function flux(u, grad_u, equations::CompressibleNaiverStokes2D) + # Here grad_u is assumed to contain the gradients of the primitive variables (v1,v2,T) + # either computed directly or reverse engineered from the gradient of the entropy vairables + # by way of the `convert_gradient_variables` function + rho, rho_v1, rho_v2, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = temperature(u, equations) + + # I was not sure what shape this array has or or if it was a tuple + # or how to properly "unpack" it. So I just guessed... + dv1dx, dv1dy, dv2dx, dv2dy, dTdx, dTdy = grad_u + + # Components of viscous stress tensor + + # (4/3*(v1)_x - 2/3*(v2)_y) + tau_xx = ( 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy ) + # ((v1)_y + (v2)_x) + tau_xy = ( dv1dy + dv2dx ) + tau_yx = tau_xy + # (4/3*(v2)_y - 2/3*(v1)_x) + tau_yy = ( 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx ) + + # Fick's law q = -kappa*grad(T); constant is kappa*gamma/(Pr*(gamma-1)) + # Important note! Due to nondimensional scaling R = 1/gamma, so the + # temperature T in the gradient computation already contains a factor of gamma + q1 = ( equations.kappa * equations.inv_gamma_minus_one * dTdx ) / equations.Pr + q2 = ( equations.kappa * equations.inv_gamma_minus_one * dTdy ) / equations.Pr + + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_xx / equations.Re + f3 = tau_xy / equations.Re + f4 = ( v1 * tau_xx + v2 * tau_xy + q1) / equations.Re + + # viscous flux components in y-direction + g1 = zero(rho) + g2 = tau_yx / equations.Re + g3 = tau_yy / equations.Re + g4 = ( v1 * tau_yx + v2 * tau_yy + q2) / equations.Re + + # TODO: I was not sure how to return this properly. Right now it is a vector of vectors + return SVector( SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4) ) +end + + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNaiverStokes2D) + rho, rho_v1, rho_v2, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = temperature(u, equations) + + return SVector(v1, v2, T) +end + + +# Convert conservative variables to entropy +@inline function cons2entropy(u, equations::CompressibleNaiverStokes2D) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) + + rho_p = rho / p + + w2 = rho_p * v1 + w3 = rho_p * v2 + w4 = -rho_p + + return SVector(w2, w3, w4) +end + + +@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNaiverStokes2D) +# Takes the solution values `u` and gradient of the variables (w_2, w_3, w_4) and +# reverse engineers the gradients to be terms of the primitive vairables (u, v, T). +# Helpful because then the diffusive fluxes have the same form as on paper. + rho, rho_v1, rho_v2, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = temperature(u, equations) + + return SVector(equations.R * T * (grad_entropy_vars(1) + v1 * grad_entropy_vars(3)), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) + equations.R * T * (grad_entropy_vars(2) + v2 * grad_entropy_vars(3)), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) + equations.R * T * T * grad_entropy_vars(3) # grad(T) = R*T^2*grad(w_4)) + ) +end + + +@inline function temperature(u, equations::CompressibleNaiverStokes2D) + rho, rho_v1, rho_v2, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) + T = p / (equations.R * rho) + return T +end + + +# +# All this boundary conditions stuff I did not touch +# + +# TODO: should this remain in the equations file, be moved to solvers, or live in the elixir? +# The penalization depends on the solver, but also depends explicitly on physical parameters, +# and would probably need to be specialized for every different equation. +# function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) +# return dg.penalty_parameter * (u_outer - u_inner) * equations.diffusivity * inv_h +# end + +# # Dirichlet-type boundary condition for use with a parabolic solver in weak form +# @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +# x, t, operator_type::Gradient, +# equations::LaplaceDiffusion2D) +# return boundary_condition.boundary_value_function(x, t, equations) +# end + +# @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +# x, t, operator_type::Divergence, +# equations::LaplaceDiffusion2D) +# return u_inner +# end + +# @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +# x, t, operator_type::Divergence, +# equations::LaplaceDiffusion2D) +# return boundary_condition.boundary_normal_flux_function(x, t, equations) +# end + +# @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +# x, t, operator_type::Gradient, +# equations::LaplaceDiffusion2D) +# return flux_inner +# end diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index fec0eed6f63..df5f7077ad1 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -1,3 +1,7 @@ # Linear scalar diffusion for use in linear scalar advection-diffusion problems abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_2d.jl") + +# Compressible Navier-Stokes equations +abstract type AbstractCompressibleNavierStokesEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +include("compressible_navier_stokes_2d.jl") From 46ae5fb558fa96acc6a0cec1e4356ec60b765034 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 7 Jun 2022 21:03:33 +0200 Subject: [PATCH 02/71] remove unneeded temperature computation --- src/equations/compressible_navier_stokes_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index cf8a8b17fc7..b695ea02361 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -95,7 +95,6 @@ function flux(u, grad_u, equations::CompressibleNaiverStokes2D) v1 = rho_v1 / rho v2 = rho_v2 / rho - T = temperature(u, equations) # I was not sure what shape this array has or or if it was a tuple # or how to properly "unpack" it. So I just guessed... From 4f6558fa778f47c1ea104957fcad3fe4612fc4e1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 7 Jun 2022 21:47:45 +0200 Subject: [PATCH 03/71] draft of elixir with possible boundary condition structure --- .../elixir_navier_stokes_source_terms.jl | 153 ++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl new file mode 100644 index 00000000000..7d8c2397272 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -0,0 +1,153 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNaiverStokes2D(1.4, # gamma + 1000, # Reynolds number + 0.72, # Prandtl number + 0.5, # free-stream Mach number + 1.0, # thermal diffusivity + equations) + +# 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) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Define initial condition +# Note: If you change the diffusion parameter here, also change it in the parabolic equation definition +function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) + # Store translated coordinate for easy use of exact solution + x_trans = x - equation.advection_velocity * t + + # @unpack nu = equation + nu = 5.0e-2 + c = 1.0 + A = 0.5 + L = 2 + f = 1/L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) + return SVector(scalar) +end + +initial_condition = initial_condition_diffusive_convergence_test + + + +# TODO: Wasn't sure of the call structure, but this should be what we need +function boundary_condition_no_slip_adiabatic_wall(u_inner, orientation, direction, + x, t, surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + + +# TODO: Wasn't sure of the call structure, but this should be what we need +function boundary_condition_no_slip_adiabatic_wall_neumann(grad_u_inner, orientation, direction, + x, t, surface_flux_function, + equations::CompressibleNaiverStokes2D) + # Copy the inner gradients to an external state array + grad_u_boundary .= grad_u_inner + + # Project out the appropriate temperature gradient pieces. Wasn't sure of array shape + if orientation == 1 + grad_u_norm = grad_u[1,3] # temperature gradient in x-direction + u_x_tangent = grad_u[1,3] - grad_u_norm + u_y_tangent = grad_u[2,3] + + # Update the external state gradients + grad_u_boundary[1,3] = u_x_tangent - grad_u_norm + grad_u_boundary[2,3] = u_y_tangent + else # orientation == 2 + grad_u_norm = grad_u[2,3] # temperature gradient in y-direction + u_x_tangent = grad_u[1,3] + u_y_tangent = grad_u[2,3] - grad_u_norm + + # Update the external state gradients + grad_u_boundary[1,3] = u_x_tangent + grad_u_boundary[2,3] = u_y_tangent - grad_u_norm + end + + # Calculate boundary flux (just averages so has no orientation I think) + flux = surface_flux_function(grad_u_inner, grad_u_boundary, orientation, equations) + + return flux +end + + +# Below is my best guess as to how to set periodic in x direction and walls +# in the y direcitons +boundary_condition= ( + x_neg=boundary_condition_periodic, + x_pos=boundary_condition_periodic, + y_neg=boundary_condition_no_slip_adiabatic_wall, + y_pos=boundary_condition_no_slip_adiabatic_wall, + ) +boundary_conditions_parabolic = ( + x_neg=boundary_condition_periodic, + x_pos=boundary_condition_periodic, + y_neg=boundary_condition_no_slip_adiabatic_wall_neumann, + y_pos=boundary_condition_no_slip_adiabatic_wall_neumann, + ) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.5 +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +time_int_tol = 1.0e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + +# Print the timer summary +summary_callback() From dde54d318fe10da4ef8157d1291218a770cbc0d0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 7 Jun 2022 22:40:05 +0200 Subject: [PATCH 04/71] added manufactured solution and source term --- .../elixir_navier_stokes_source_terms.jl | 149 ++++++++++++++++-- 1 file changed, 134 insertions(+), 15 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index 7d8c2397272..f8cf3782759 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -26,24 +26,143 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max=30_000) # set maximum capacity of tree data structure # Define initial condition -# Note: If you change the diffusion parameter here, also change it in the parabolic equation definition -function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) - # Store translated coordinate for easy use of exact solution - x_trans = x - equation.advection_velocity * t - - # @unpack nu = equation - nu = 5.0e-2 - c = 1.0 - A = 0.5 - L = 2 - f = 1/L - omega = 2 * pi * f - scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) - return SVector(scalar) +# Note: If you change the parameters here, also change it in the corresponding source terms +function initial_condition_navier_stokes_convergence_test(x, t, equation::CompressibleEulerEquations2D) + + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) end -initial_condition = initial_condition_diffusive_convergence_test +initial_condition = initial_condition_navier_stokes_convergence_test + + +@inline function initial_condition_navier_stokes_convergence_test(u, x, t, + equations::CompressibleNaiverStokes2D) + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * equations.inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * equations.inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * equations.inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * equations.inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * equations.inv_gamma_minus_one * equations.kappa / equations.Pr + inv_Re = 1.0 / equations.Re + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # 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 + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * inv_Re + + 2.0 / 3.0 * v2_xy * inv_Re + - v1_yy * inv_Re + - v2_xy * inv_Re ) + # 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 + # stress tensor from y-direction + - v1_xy * inv_Re + - v2_xx * inv_Re + - 4.0 / 3.0 * v2_yy * inv_Re + + 2.0 / 3.0 * v1_xy * inv_Re ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * inv_Re + + 2.0 / 3.0 * v2_xy * v1 * inv_Re + - 4.0 / 3.0 * v1_x * v1_x * inv_Re + + 2.0 / 3.0 * v2_y * v1_x * inv_Re + - v1_xy * v2 * inv_Re + - v2_xx * v2 * inv_Re + - v1_y * v2_x * inv_Re + - v2_x * v2_x * inv_Re + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * inv_Re + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * inv_Re + - v2_xy * v1 * inv_Re + - v1_y * v1_y * inv_Re + - v2_x * v1_y * inv_Re + - 4.0 / 3.0 * v2_yy * v2 * inv_Re + + 2.0 / 3.0 * v1_xy * v2 * inv_Re + - 4.0 / 3.0 * v2_y * v2_y * inv_Re + + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * inv_Re ) + + return SVector(du1, du2, du3, du4) +end # TODO: Wasn't sure of the call structure, but this should be what we need From daaeac3a35539f6de9d1bf35b7150f19aeee288f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 7 Jun 2022 22:46:22 +0200 Subject: [PATCH 05/71] fix typo in function name for MMS --- examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index f8cf3782759..adbff3101aa 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -50,8 +50,8 @@ end initial_condition = initial_condition_navier_stokes_convergence_test -@inline function initial_condition_navier_stokes_convergence_test(u, x, t, - equations::CompressibleNaiverStokes2D) +@inline function source_terms_navier_stokes_convergence_test(u, x, t, + equations::CompressibleNaiverStokes2D) # Same settings as in `initial_condition` # Amplitude and shift A = 0.5 From 4a39883f81e3adf47f16dbf7c8008bca6061b925 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 8 Jun 2022 06:59:52 +0200 Subject: [PATCH 06/71] update variable names for consistency. improve comments --- .../compressible_navier_stokes_2d.jl | 36 +++++++++++-------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index b695ea02361..6433eabf0c7 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -81,10 +81,13 @@ varnames(variable_mapping, equations_parabolic::CompressibleNaiverStokes2D) = # no orientation specified since the flux is vector-valued -# This form of the difussive Navier-Stokes fluxes is taken from Section 2 +# Explicit formulas for the diffussive Navier-Stokes fluxes are avilable, e.g. in Section 2 # of the paper by Svärd, Carpenter and Nordström # "A stable high-order finite difference scheme for the compressible Navier–Stokes # equations, far-field boundary conditions" +# Although these authors use a different nondimensionalization so some constants are different +# particularly for Fick's law. +# # Note, could be generalized to use Sutherland's law to get the molecular and thermal # diffusivity function flux(u, grad_u, equations::CompressibleNaiverStokes2D) @@ -103,12 +106,12 @@ function flux(u, grad_u, equations::CompressibleNaiverStokes2D) # Components of viscous stress tensor # (4/3*(v1)_x - 2/3*(v2)_y) - tau_xx = ( 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy ) + tau_11 = ( 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy ) # ((v1)_y + (v2)_x) - tau_xy = ( dv1dy + dv2dx ) - tau_yx = tau_xy + tau_12 = ( dv1dy + dv2dx ) + tau_21 = tau_12 # (4/3*(v2)_y - 2/3*(v1)_x) - tau_yy = ( 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx ) + tau_22 = ( 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx ) # Fick's law q = -kappa*grad(T); constant is kappa*gamma/(Pr*(gamma-1)) # Important note! Due to nondimensional scaling R = 1/gamma, so the @@ -116,17 +119,20 @@ function flux(u, grad_u, equations::CompressibleNaiverStokes2D) q1 = ( equations.kappa * equations.inv_gamma_minus_one * dTdx ) / equations.Pr q2 = ( equations.kappa * equations.inv_gamma_minus_one * dTdy ) / equations.Pr + # molecular diffusivity is simply 1/Re for this nondimensionalization + mu = 1.0 / equations.Re + # viscous flux components in the x-direction f1 = zero(rho) - f2 = tau_xx / equations.Re - f3 = tau_xy / equations.Re - f4 = ( v1 * tau_xx + v2 * tau_xy + q1) / equations.Re + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = ( v1 * tau_11 + v2 * tau_12 + q1 ) * mu # viscous flux components in y-direction g1 = zero(rho) - g2 = tau_yx / equations.Re - g3 = tau_yy / equations.Re - g4 = ( v1 * tau_yx + v2 * tau_yy + q2) / equations.Re + g2 = tau_21 * mu + g3 = tau_22 * mu + g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * mu # TODO: I was not sure how to return this properly. Right now it is a vector of vectors return SVector( SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4) ) @@ -166,7 +172,7 @@ end @inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNaiverStokes2D) # Takes the solution values `u` and gradient of the variables (w_2, w_3, w_4) and -# reverse engineers the gradients to be terms of the primitive vairables (u, v, T). +# reverse engineers the gradients to be terms of the primitive vairables (v1, v2, T). # Helpful because then the diffusive fluxes have the same form as on paper. rho, rho_v1, rho_v2, _ = u @@ -174,9 +180,9 @@ end v2 = rho_v2 / rho T = temperature(u, equations) - return SVector(equations.R * T * (grad_entropy_vars(1) + v1 * grad_entropy_vars(3)), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) - equations.R * T * (grad_entropy_vars(2) + v2 * grad_entropy_vars(3)), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) - equations.R * T * T * grad_entropy_vars(3) # grad(T) = R*T^2*grad(w_4)) + return SVector(equations.R * T * (grad_entropy_vars[1] + v1 * grad_entropy_vars[3]), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) + equations.R * T * (grad_entropy_vars[2] + v2 * grad_entropy_vars[3]), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) + equations.R * T * T * grad_entropy_vars[3] # grad(T) = R*T^2*grad(w_4)) ) end From 3879a3fc65e4c30097842b49b831f48472be7b64 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 8 Jun 2022 13:41:40 +0200 Subject: [PATCH 07/71] fix dumb typos in new equation system name --- .../elixir_navier_stokes_source_terms.jl | 6 ++--- .../compressible_navier_stokes_2d.jl | 24 +++++++++---------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index adbff3101aa..7203c321e2a 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -7,7 +7,7 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNaiverStokes2D(1.4, # gamma +equations_parabolic = CompressibleNavierStokes2D(1.4, # gamma 1000, # Reynolds number 0.72, # Prandtl number 0.5, # free-stream Mach number @@ -51,7 +51,7 @@ initial_condition = initial_condition_navier_stokes_convergence_test @inline function source_terms_navier_stokes_convergence_test(u, x, t, - equations::CompressibleNaiverStokes2D) + equations::CompressibleNavierStokes2D) # Same settings as in `initial_condition` # Amplitude and shift A = 0.5 @@ -185,7 +185,7 @@ end # TODO: Wasn't sure of the call structure, but this should be what we need function boundary_condition_no_slip_adiabatic_wall_neumann(grad_u_inner, orientation, direction, x, t, surface_flux_function, - equations::CompressibleNaiverStokes2D) + equations::CompressibleNavierStokes2D) # Copy the inner gradients to an external state array grad_u_boundary .= grad_u_inner diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 6433eabf0c7..6ee74893c01 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,12 +1,12 @@ @doc raw""" - CompressibleNaiverStokes2D(gamma, + CompressibleNavierStokes2D(gamma, Re, Pr, Ma_inf, kappa, equations) -`CompressibleNaiverStokes2D` represents the diffusion (i.e. parabolic) terms applied +`CompressibleNavierStokes2D` represents the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective from the `CompressibleEulerEquations2D`. @@ -45,7 +45,7 @@ where # 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from # the Euler equations # 2) Add more here and probably some equations -struct CompressibleNaiverStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} +struct CompressibleNavierStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications Re::RealT # Reynolds number @@ -60,15 +60,15 @@ struct CompressibleNaiverStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierS equations::E # CompressibleEulerEquations2D end -function CompressibleNaiverStokes2D(gamma, Reynolds, Prandtl, Mach_frestream, kappa, equations) +function CompressibleNavierStokes2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations) γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) # From the nondimensionalization discussed above set the remaining free-stream # quantities p_inf = 1.0 / γ - u_inf = Mach_frestream + u_inf = Mach_freestream R = 1.0 / γ - CompressibleNaiverStokes2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, + CompressibleNavierStokes2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, Reynolds, Prandtl, Mach_freestream, kappa, p_inf, u_inf, R, equations) @@ -76,7 +76,7 @@ end # I was not sure what to do here to allow flexibility of selecting primitive or entropy # grandient variables -varnames(variable_mapping, equations_parabolic::CompressibleNaiverStokes2D) = +varnames(variable_mapping, equations_parabolic::CompressibleNavierStokes2D) = varnames(variable_mapping, equations_parabolic.equations) @@ -90,7 +90,7 @@ varnames(variable_mapping, equations_parabolic::CompressibleNaiverStokes2D) = # # Note, could be generalized to use Sutherland's law to get the molecular and thermal # diffusivity -function flux(u, grad_u, equations::CompressibleNaiverStokes2D) +function flux(u, grad_u, equations::CompressibleNavierStokes2D) # Here grad_u is assumed to contain the gradients of the primitive variables (v1,v2,T) # either computed directly or reverse engineered from the gradient of the entropy vairables # by way of the `convert_gradient_variables` function @@ -140,7 +140,7 @@ end # Convert conservative variables to primitive -@inline function cons2prim(u, equations::CompressibleNaiverStokes2D) +@inline function cons2prim(u, equations::CompressibleNavierStokes2D) rho, rho_v1, rho_v2, _ = u v1 = rho_v1 / rho @@ -152,7 +152,7 @@ end # Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleNaiverStokes2D) +@inline function cons2entropy(u, equations::CompressibleNavierStokes2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho @@ -170,7 +170,7 @@ end end -@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNaiverStokes2D) +@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNavierStokes2D) # Takes the solution values `u` and gradient of the variables (w_2, w_3, w_4) and # reverse engineers the gradients to be terms of the primitive vairables (v1, v2, T). # Helpful because then the diffusive fluxes have the same form as on paper. @@ -187,7 +187,7 @@ end end -@inline function temperature(u, equations::CompressibleNaiverStokes2D) +@inline function temperature(u, equations::CompressibleNavierStokes2D) rho, rho_v1, rho_v2, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) From 435e081ee11a2c9c1f0e2620d823f82f2e96bbca Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 8 Jun 2022 13:41:52 +0200 Subject: [PATCH 08/71] actually export new equations --- src/Trixi.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index e78a1bd9dc5..c176593e64a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -133,6 +133,7 @@ export AcousticPerturbationEquations2D, InviscidBurgersEquation1D, LaplaceDiffusion2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, + CompressibleNavierStokes2D, ShallowWaterEquations1D, ShallowWaterEquations2D export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, From 7b59fef47ef5cc6df80a194918af6ecf5dae9d5d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 8 Jun 2022 13:54:47 +0200 Subject: [PATCH 09/71] add comment near variable_mapping declaration. --- src/equations/compressible_navier_stokes_2d.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 6ee74893c01..4a364a3e9ab 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -74,8 +74,14 @@ function CompressibleNavierStokes2D(gamma, Reynolds, Prandtl, Mach_freestream, k equations) end + # I was not sure what to do here to allow flexibility of selecting primitive or entropy -# grandient variables +# gradient variables. I see that `transform_variables!` just copies data at the moment. + +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokes2D) = ("v1", "v2", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokes2D) = ("w2", "w3", "w4") + varnames(variable_mapping, equations_parabolic::CompressibleNavierStokes2D) = varnames(variable_mapping, equations_parabolic.equations) From 5dc830154b2ae77b4c1c456c7b5be9f1c2e65c05 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 9 Jun 2022 08:33:58 +0200 Subject: [PATCH 10/71] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl | 2 +- src/equations/compressible_navier_stokes_2d.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index 7203c321e2a..b8796315f4d 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -2,7 +2,7 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the linear advection-diffusion equation +# semidiscretization of the ideal compressible Navier-Stokes equations equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 4a364a3e9ab..e19a4eadded 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -65,9 +65,9 @@ function CompressibleNavierStokes2D(gamma, Reynolds, Prandtl, Mach_freestream, k # From the nondimensionalization discussed above set the remaining free-stream # quantities - p_inf = 1.0 / γ + p_inf = 1 / γ u_inf = Mach_freestream - R = 1.0 / γ + R = 1 / γ CompressibleNavierStokes2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, Reynolds, Prandtl, Mach_freestream, kappa, p_inf, u_inf, R, @@ -141,7 +141,7 @@ function flux(u, grad_u, equations::CompressibleNavierStokes2D) g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * mu # TODO: I was not sure how to return this properly. Right now it is a vector of vectors - return SVector( SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4) ) + return (SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4)) end From f316684fefa6d1da10fcecde4e6eb739d432a557 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 9 Jun 2022 08:28:13 +0200 Subject: [PATCH 11/71] parabolic equations now exported separately --- src/Trixi.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index c176593e64a..4037d41961c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -131,11 +131,12 @@ export AcousticPerturbationEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, InviscidBurgersEquation1D, - LaplaceDiffusion2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, - CompressibleNavierStokes2D, ShallowWaterEquations1D, ShallowWaterEquations2D +export LaplaceDiffusion2D, + CompressibleNavierStokes2D + export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, From 02c14e1ed9ff78c94f451cec23a308424b830a11 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 9 Jun 2022 08:30:52 +0200 Subject: [PATCH 12/71] change name to CompressibleNavierStokesEquations2D --- .../elixir_navier_stokes_source_terms.jl | 6 +-- .../compressible_navier_stokes_2d.jl | 44 +++++++++---------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index b8796315f4d..6dcf22b41fa 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -7,7 +7,7 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokes2D(1.4, # gamma +equations_parabolic = CompressibleNavierStokesEquations2D(1.4, # gamma 1000, # Reynolds number 0.72, # Prandtl number 0.5, # free-stream Mach number @@ -51,7 +51,7 @@ initial_condition = initial_condition_navier_stokes_convergence_test @inline function source_terms_navier_stokes_convergence_test(u, x, t, - equations::CompressibleNavierStokes2D) + equations::CompressibleNavierStokesEquations2D) # Same settings as in `initial_condition` # Amplitude and shift A = 0.5 @@ -185,7 +185,7 @@ end # TODO: Wasn't sure of the call structure, but this should be what we need function boundary_condition_no_slip_adiabatic_wall_neumann(grad_u_inner, orientation, direction, x, t, surface_flux_function, - equations::CompressibleNavierStokes2D) + equations::CompressibleNavierStokesEquations2D) # Copy the inner gradients to an external state array grad_u_boundary .= grad_u_inner diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index e19a4eadded..ffe6cf1d95b 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,12 +1,12 @@ @doc raw""" - CompressibleNavierStokes2D(gamma, - Re, - Pr, - Ma_inf, - kappa, - equations) - -`CompressibleNavierStokes2D` represents the diffusion (i.e. parabolic) terms applied + CompressibleNavierStokesEquations2D(gamma, + Re, + Pr, + Ma_inf, + kappa, + equations) + +`CompressibleNavierStokesEquations2D` represents the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective from the `CompressibleEulerEquations2D`. @@ -45,7 +45,7 @@ where # 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from # the Euler equations # 2) Add more here and probably some equations -struct CompressibleNavierStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} +struct CompressibleNavierStokesEquations2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications Re::RealT # Reynolds number @@ -60,7 +60,7 @@ struct CompressibleNavierStokes2D{RealT<:Real, E} <: AbstractCompressibleNavierS equations::E # CompressibleEulerEquations2D end -function CompressibleNavierStokes2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations) +function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations) γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) # From the nondimensionalization discussed above set the remaining free-stream @@ -68,10 +68,10 @@ function CompressibleNavierStokes2D(gamma, Reynolds, Prandtl, Mach_freestream, k p_inf = 1 / γ u_inf = Mach_freestream R = 1 / γ - CompressibleNavierStokes2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, - Reynolds, Prandtl, Mach_freestream, - kappa, p_inf, u_inf, R, - equations) + CompressibleNavierStokesEquations2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, + Reynolds, Prandtl, Mach_freestream, + kappa, p_inf, u_inf, R, + equations) end @@ -79,10 +79,10 @@ end # gradient variables. I see that `transform_variables!` just copies data at the moment. # This is the flexibility a user should have to select the different gradient variable types -# varnames(::typeof(cons2prim) , ::CompressibleNavierStokes2D) = ("v1", "v2", "T") -# varnames(::typeof(cons2entropy), ::CompressibleNavierStokes2D) = ("w2", "w3", "w4") +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesEquations2D) = ("v1", "v2", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesEquations2D) = ("w2", "w3", "w4") -varnames(variable_mapping, equations_parabolic::CompressibleNavierStokes2D) = +varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesEquations2D) = varnames(variable_mapping, equations_parabolic.equations) @@ -96,7 +96,7 @@ varnames(variable_mapping, equations_parabolic::CompressibleNavierStokes2D) = # # Note, could be generalized to use Sutherland's law to get the molecular and thermal # diffusivity -function flux(u, grad_u, equations::CompressibleNavierStokes2D) +function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) # Here grad_u is assumed to contain the gradients of the primitive variables (v1,v2,T) # either computed directly or reverse engineered from the gradient of the entropy vairables # by way of the `convert_gradient_variables` function @@ -146,7 +146,7 @@ end # Convert conservative variables to primitive -@inline function cons2prim(u, equations::CompressibleNavierStokes2D) +@inline function cons2prim(u, equations::CompressibleNavierStokesEquations2D) rho, rho_v1, rho_v2, _ = u v1 = rho_v1 / rho @@ -158,7 +158,7 @@ end # Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleNavierStokes2D) +@inline function cons2entropy(u, equations::CompressibleNavierStokesEquations2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho @@ -176,7 +176,7 @@ end end -@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNavierStokes2D) +@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNavierStokesEquations2D) # Takes the solution values `u` and gradient of the variables (w_2, w_3, w_4) and # reverse engineers the gradients to be terms of the primitive vairables (v1, v2, T). # Helpful because then the diffusive fluxes have the same form as on paper. @@ -193,7 +193,7 @@ end end -@inline function temperature(u, equations::CompressibleNavierStokes2D) +@inline function temperature(u, equations::CompressibleNavierStokesEquations2D) rho, rho_v1, rho_v2, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) From dd733bb4d05b69813bb2ed342465e388dc480c67 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 9 Jun 2022 08:37:36 +0200 Subject: [PATCH 13/71] export NS with proper name --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 4037d41961c..2b4ac68d1f5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -135,7 +135,7 @@ export AcousticPerturbationEquations2D, ShallowWaterEquations1D, ShallowWaterEquations2D export LaplaceDiffusion2D, - CompressibleNavierStokes2D + CompressibleNavierStokesEquations2D export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, From a1a48fe2b0627c55c2d7487997e540c7283639d3 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 10:17:27 +0200 Subject: [PATCH 14/71] explicitly require compressible Euler in the type parameter --- src/equations/compressible_navier_stokes_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index ffe6cf1d95b..291fdcf8dfe 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -45,7 +45,7 @@ where # 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from # the Euler equations # 2) Add more here and probably some equations -struct CompressibleNavierStokesEquations2D{RealT<:Real, E} <: AbstractCompressibleNavierStokesEquations{2, 3} +struct CompressibleNavierStokesEquations2D{RealT<:Real, E<:AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesEquations{2, 3} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications Re::RealT # Reynolds number @@ -60,7 +60,7 @@ struct CompressibleNavierStokesEquations2D{RealT<:Real, E} <: AbstractCompressib equations::E # CompressibleEulerEquations2D end -function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations) +function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations::CompressibleEulerEquations2D) γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) # From the nondimensionalization discussed above set the remaining free-stream From d8db2030d94f8aecf6af97d8a9d3825d77c5e5e4 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 10:19:58 +0200 Subject: [PATCH 15/71] name kinematic viscosity nu for consistency with Lattice-Boltzmann --- src/equations/compressible_navier_stokes_2d.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 291fdcf8dfe..dc11d84a82e 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -125,22 +125,21 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) q1 = ( equations.kappa * equations.inv_gamma_minus_one * dTdx ) / equations.Pr q2 = ( equations.kappa * equations.inv_gamma_minus_one * dTdy ) / equations.Pr - # molecular diffusivity is simply 1/Re for this nondimensionalization - mu = 1.0 / equations.Re + # kinematic viscosity is simply 1/Re for this nondimensionalization + nu = 1.0 / equations.Re # viscous flux components in the x-direction f1 = zero(rho) - f2 = tau_11 * mu - f3 = tau_12 * mu - f4 = ( v1 * tau_11 + v2 * tau_12 + q1 ) * mu + f2 = tau_11 * nu + f3 = tau_12 * nu + f4 = ( v1 * tau_11 + v2 * tau_12 + q1 ) * nu # viscous flux components in y-direction g1 = zero(rho) - g2 = tau_21 * mu - g3 = tau_22 * mu - g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * mu + g2 = tau_21 * nu + g3 = tau_22 * nu + g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * nu - # TODO: I was not sure how to return this properly. Right now it is a vector of vectors return (SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4)) end From 06e69b869ec215ed6a314f5d07442e12de051982 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 10:23:43 +0200 Subject: [PATCH 16/71] add promotion in constructor --- src/equations/compressible_navier_stokes_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index dc11d84a82e..236fa9298cf 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -61,7 +61,7 @@ struct CompressibleNavierStokesEquations2D{RealT<:Real, E<:AbstractCompressibleE end function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations::CompressibleEulerEquations2D) - γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + γ, inv_gamma_minus_one, Re, Pr, Ma, κ = promote(gamma, inv(gamma - 1), Reynolds, Prandtl, Mach_freestream, kappa) # From the nondimensionalization discussed above set the remaining free-stream # quantities @@ -69,8 +69,8 @@ function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_free u_inf = Mach_freestream R = 1 / γ CompressibleNavierStokesEquations2D{typeof(γ),typeof(equations)}(γ, inv_gamma_minus_one, - Reynolds, Prandtl, Mach_freestream, - kappa, p_inf, u_inf, R, + Re, Pr, Ma, κ, + p_inf, u_inf, R, equations) end From eda9737f042b05fe7a4461d4027630c84095ff66 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 10:38:11 +0200 Subject: [PATCH 17/71] make Reynolds, Prandtl, Mach, and kappa keyword arguments --- src/equations/compressible_navier_stokes_2d.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 236fa9298cf..73e5ad55c2d 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -60,8 +60,10 @@ struct CompressibleNavierStokesEquations2D{RealT<:Real, E<:AbstractCompressibleE equations::E # CompressibleEulerEquations2D end -function CompressibleNavierStokesEquations2D(gamma, Reynolds, Prandtl, Mach_freestream, kappa, equations::CompressibleEulerEquations2D) - γ, inv_gamma_minus_one, Re, Pr, Ma, κ = promote(gamma, inv(gamma - 1), Reynolds, Prandtl, Mach_freestream, kappa) +function CompressibleNavierStokesEquations2D(equations::CompressibleEulerEquations2D; Reynolds, Prandtl, Mach_freestream, kappa) + γ = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + Re, Pr, Ma, κ = promote(Reynolds, Prandtl, Mach_freestream, kappa) # From the nondimensionalization discussed above set the remaining free-stream # quantities From b469bcea07e6f7403d58a24bbccdcf0090e512b6 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 10:38:38 +0200 Subject: [PATCH 18/71] update constructor call in elixir --- .../elixir_navier_stokes_source_terms.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl index 6dcf22b41fa..bd85946f1b3 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl @@ -7,12 +7,12 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesEquations2D(1.4, # gamma - 1000, # Reynolds number - 0.72, # Prandtl number - 0.5, # free-stream Mach number - 1.0, # thermal diffusivity - equations) +equations_parabolic = CompressibleNavierStokesEquations2D(equations, + Reynolds=1000, + Prandtl=0.72, + Mach_freestream=0.5, + kappa=1.0 # thermal diffusivity + ) # 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) From fb646b0e765a4667a3b4da246eff463388e15a33 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 10 Jun 2022 12:07:30 +0200 Subject: [PATCH 19/71] reduce computation by exploiting stress tensor symmetry --- src/equations/compressible_navier_stokes_2d.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 73e5ad55c2d..4683d7d01ea 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -116,8 +116,8 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) # (4/3*(v1)_x - 2/3*(v2)_y) tau_11 = ( 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy ) # ((v1)_y + (v2)_x) - tau_12 = ( dv1dy + dv2dx ) - tau_21 = tau_12 + # stress tensor is symmetric + tau_12 = ( dv1dy + dv2dx ) # = tau_21 # (4/3*(v2)_y - 2/3*(v1)_x) tau_22 = ( 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx ) @@ -136,11 +136,13 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) f3 = tau_12 * nu f4 = ( v1 * tau_11 + v2 * tau_12 + q1 ) * nu - # viscous flux components in y-direction + # viscous flux components in the y-direction + # Note, symmetry is exploited for tau_12 = tau_21 g1 = zero(rho) - g2 = tau_21 * nu + g2 = f3 # tau_21 * nu g3 = tau_22 * nu - g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * nu + # g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * nu + g4 = ( v1 * tau_12 + v2 * tau_22 + q2 ) * nu return (SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4)) end From f6bca767de53032723ccbece8558c4d1892ccedc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 11:58:42 -0500 Subject: [PATCH 20/71] fix unpacking of flux --- src/equations/compressible_navier_stokes_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 4683d7d01ea..9adf98dc1b2 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -109,7 +109,8 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) # I was not sure what shape this array has or or if it was a tuple # or how to properly "unpack" it. So I just guessed... - dv1dx, dv1dy, dv2dx, dv2dy, dTdx, dTdy = grad_u + dv1dx, dv2dx, dTdx = grad_u[1] + dv1dy, dv2dy, dTdy = grad_u[2] # Components of viscous stress tensor From 9b5dd3596167dbc93463dfea76b7c8abec83f354 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 12:04:29 -0500 Subject: [PATCH 21/71] modifying parabolic cache creation in cache, we assume we take the gradient of all hyperbolic variables. since the number of parabolic variables can differ from the number of hyperbolic variables, we pass in the hyperbolic equations to `create_cache_parabolic` now --- .../semidiscretization_hyperbolic_parabolic.jl | 3 ++- src/solvers/dgmulti/dg_parabolic.jl | 7 +++++-- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 7 ++++--- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index ea96aa5c363..52a48d7aa03 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -102,7 +102,8 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache) - cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, solver_parabolic, RealT, uEltype)..., + cache_parabolic = (; create_cache_parabolic(mesh, equations, equations_parabolic, + solver, solver_parabolic, RealT, uEltype)..., initial_cache_parabolic...) SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 9578e3f9eff..cfb92fd1b4f 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,6 +1,9 @@ -function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsParabolic, +function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, dg::DGMulti, dg_parabolic, RealT, uEltype) - nvars = nvariables(equations) + # default to taking derivatives of all hyperbolic terms + # TODO: utilize "differentiated" parabolic variables in `equations_parabolic` + nvars = nvariables(equations_hyperbolic) @unpack M, Drst = dg.basis weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 01f0ead469f..1c332fc6f2e 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -438,14 +438,15 @@ 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{2}, equations_parabolic::AbstractEquationsParabolic, +function create_cache_parabolic(mesh::TreeMesh{2}, equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, dg::DG, dg_parabolic, 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_parabolic, dg.basis, RealT, uEltype) + elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) - n_vars = nvariables(equations_parabolic) + 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_elements) From 21e1c998e3b3c4f5cfc914bbeb45fade3b82cefb Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:29:21 -0500 Subject: [PATCH 22/71] comments --- src/solvers/dgmulti/dg_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index cfb92fd1b4f..561dfd2f1a4 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -248,7 +248,7 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh end # interpolates from solution coefficients to face quadrature points - viscous_flux_face_values = cache_parabolic.grad_u_face_values + viscous_flux_face_values = cache_parabolic.grad_u_face_values # reuse storage for dim in eachdim(mesh) prolong2interfaces!(viscous_flux_face_values[dim], viscous_flux[dim], mesh, equations, dg.surface_integral, dg, cache) From 19f004e6de741a921c1adb4716bf593168de09bd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:29:32 -0500 Subject: [PATCH 23/71] comments --- src/solvers/dgmulti/dg_parabolic.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 561dfd2f1a4..3ad91627f63 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,8 +1,9 @@ -function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, +function create_cache_parabolic(mesh::DGMultiMesh, + equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, dg::DGMulti, dg_parabolic, RealT, uEltype) # default to taking derivatives of all hyperbolic terms - # TODO: utilize "differentiated" parabolic variables in `equations_parabolic` + # TODO: utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache nvars = nvariables(equations_hyperbolic) @unpack M, Drst = dg.basis From de442bf8e830b46007eca3283a1033fcedfabbf9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:30:57 -0500 Subject: [PATCH 24/71] formatting and renaming equations to equations_hyperbolic formatting comments --- src/equations/compressible_navier_stokes_2d.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 9adf98dc1b2..e5c34687f19 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -45,7 +45,7 @@ where # 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from # the Euler equations # 2) Add more here and probably some equations -struct CompressibleNavierStokesEquations2D{RealT<:Real, E<:AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesEquations{2, 3} +struct CompressibleNavierStokesEquations2D{RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesEquations{2, 3} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications Re::RealT # Reynolds number @@ -57,7 +57,7 @@ struct CompressibleNavierStokesEquations2D{RealT<:Real, E<:AbstractCompressibleE u_inf::RealT # free-stream velocity R::RealT # gas constant (depends on nondimensional scaling!) - equations::E # CompressibleEulerEquations2D + equations_hyperbolic::E # CompressibleEulerEquations2D end function CompressibleNavierStokesEquations2D(equations::CompressibleEulerEquations2D; Reynolds, Prandtl, Mach_freestream, kappa) @@ -85,11 +85,11 @@ end # varnames(::typeof(cons2entropy), ::CompressibleNavierStokesEquations2D) = ("w2", "w3", "w4") varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesEquations2D) = - varnames(variable_mapping, equations_parabolic.equations) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) # no orientation specified since the flux is vector-valued -# Explicit formulas for the diffussive Navier-Stokes fluxes are avilable, e.g. in Section 2 +# Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2 # of the paper by Svärd, Carpenter and Nordström # "A stable high-order finite difference scheme for the compressible Navier–Stokes # equations, far-field boundary conditions" @@ -99,7 +99,7 @@ varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesEquation # Note, could be generalized to use Sutherland's law to get the molecular and thermal # diffusivity function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) - # Here grad_u is assumed to contain the gradients of the primitive variables (v1,v2,T) + # Here grad_u is assumed to contain the gradients of the primitive variables (v1, v2, T) # either computed directly or reverse engineered from the gradient of the entropy vairables # by way of the `convert_gradient_variables` function rho, rho_v1, rho_v2, _ = u From 298b726ad601e44079e7b85ce6b8feb7f15b2840 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:31:11 -0500 Subject: [PATCH 25/71] fix unpacking of gradients in flux --- src/equations/compressible_navier_stokes_2d.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index e5c34687f19..170a6de5acb 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -107,10 +107,9 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) v1 = rho_v1 / rho v2 = rho_v2 / rho - # I was not sure what shape this array has or or if it was a tuple - # or how to properly "unpack" it. So I just guessed... - dv1dx, dv2dx, dTdx = grad_u[1] - dv1dy, dv2dy, dTdy = grad_u[2] + # grad_u contains derivatives of each hyperbolic variable + _, dv1dx, dv2dx, dTdx = grad_u[1] + _, dv1dy, dv2dy, dTdy = grad_u[2] # Components of viscous stress tensor From 01f832517e90cd0358dbfc441dedce8ac9d54483 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:31:28 -0500 Subject: [PATCH 26/71] adding CNS BCs --- src/Trixi.jl | 3 +- .../compressible_navier_stokes_2d.jl | 65 +++++++++---------- 2 files changed, 33 insertions(+), 35 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 2b4ac68d1f5..76482b841c1 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -163,7 +163,8 @@ export boundary_condition_do_nothing, BoundaryConditionNeumann, boundary_condition_noslip_wall, boundary_condition_slip_wall, - boundary_condition_wall + boundary_condition_wall, + BoundaryConditionViscousWall, NoSlip, Adiabatic export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 170a6de5acb..6e302221fa9 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -204,39 +204,36 @@ end return T end +# abstract container for wall-type boundary conditions +struct BoundaryConditionViscousWall{V, H} + boundary_condition_velocity::V + boundary_condition_heat_flux::H +end -# -# All this boundary conditions stuff I did not touch -# +# no slip velocity BC +struct NoSlip{F} + boundary_value_function::F # value of the velocity vector on the boundary +end -# TODO: should this remain in the equations file, be moved to solvers, or live in the elixir? -# The penalization depends on the solver, but also depends explicitly on physical parameters, -# and would probably need to be specialized for every different equation. -# function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) -# return dg.penalty_parameter * (u_outer - u_inner) * equations.diffusivity * inv_h -# end - -# # Dirichlet-type boundary condition for use with a parabolic solver in weak form -# @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, -# x, t, operator_type::Gradient, -# equations::LaplaceDiffusion2D) -# return boundary_condition.boundary_value_function(x, t, equations) -# end - -# @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, -# x, t, operator_type::Divergence, -# equations::LaplaceDiffusion2D) -# return u_inner -# end - -# @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, -# x, t, operator_type::Divergence, -# equations::LaplaceDiffusion2D) -# return boundary_condition.boundary_normal_flux_function(x, t, equations) -# end - -# @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, -# x, t, operator_type::Gradient, -# equations::LaplaceDiffusion2D) -# return flux_inner -# end +# adiabatic temperature BC +struct Adiabatic{F} + boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn +end + +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesEquations2D) + rho = u_inner[1] # should not matter since the viscous terms don't depend on rho + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + return SVector(rho, v1, v2, u_inner[4]) +end + +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(flux_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesEquations2D) + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + tau_1n, tau_2n = flux_inner[2:3] + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux + return SVector(flux_inner[1:3]..., normal_energy_flux) +end From 3f04e7148ec1cd48b5f2c5b8364b3f8fc28978ae Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 22:31:35 -0500 Subject: [PATCH 27/71] adding lid-driven cavity elixir --- .../elixir_navier_stokes_lid_driven_cavity.jl | 70 +++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..63bb1ce29f8 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl @@ -0,0 +1,70 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=100, Prandtl=0.72, + Mach_freestream=0.1, kappa=1.0) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top(x, tol=50*eps()) = abs(x[2] - 1) < tol +rest_of_boundary(x, tol=50*eps()) = !top(x, tol) +is_on_boundary = Dict(:top => top, :rest_of_boundary => rest_of_boundary) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = 0.1 + rho = 1.0 + u, v = 0, 0 + p = 1.0 / (Ma^2 * equations.gamma) + return prim2cons(SVector(rho, u, v, p), equations) +end +initial_condition = initial_condition_cavity + +# BC types +velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0)) +velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_lid = BoundaryConditionViscousWall(velocity_bc_lid, heat_bc) +boundary_condition_cavity = BoundaryConditionViscousWall(velocity_bc_cavity, heat_bc) + +# define inviscid boundary conditions +boundary_conditions = (; :top => boundary_condition_slip_wall, + :rest_of_boundary => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top => boundary_condition_lid, + :rest_of_boundary => boundary_condition_cavity) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.5 +tspan = (0.0, 25.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 36e3b68708eefdea12b8f6b59930111a90ba18c6 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 14 Jun 2022 23:26:39 -0500 Subject: [PATCH 28/71] adding variable transform, editing cons2prim for CNS --- src/equations/compressible_navier_stokes_2d.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 6e302221fa9..18005104387 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -87,6 +87,13 @@ end varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesEquations2D) = varnames(variable_mapping, equations_parabolic.equations_hyperbolic) +# transform to primitive variables +# TODO: should we have this call a "gradient transformation" field? +function transform_variables!(u_transformed, u, equations_parabolic::CompressibleNavierStokesEquations2D) + @threaded for i in eachindex(u) + u_transformed[i] = cons2prim(u[i], equations_parabolic) + end +end # no orientation specified since the flux is vector-valued # Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2 @@ -156,7 +163,7 @@ end v2 = rho_v2 / rho T = temperature(u, equations) - return SVector(v1, v2, T) + return SVector(rho, v1, v2, T) end From 6bf91eb80eaf66642b1c1d283de29cbc2ad96e1f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 15 Jun 2022 10:18:38 -0500 Subject: [PATCH 29/71] add prim2cons for NS (inconsistent right now) --- src/equations/compressible_navier_stokes_2d.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 18005104387..ea360b11f19 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -166,6 +166,11 @@ end return SVector(rho, v1, v2, T) end +# TODO: make this consistent with cons2prim above and cons2prim for Euler! +@inline prim2cons(u, equations::CompressibleNavierStokesEquations2D) = + prim2cons(u, equations.equations_hyperbolic) + + # Convert conservative variables to entropy @inline function cons2entropy(u, equations::CompressibleNavierStokesEquations2D) From a0ea3f0958abb2eb8d1ddc7160a08199dcf5044b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 15 Jun 2022 10:18:55 -0500 Subject: [PATCH 30/71] add draft of DGMulti Navier Stokes convergence elixir --- .../elixir_navier_stokes_convergence.jl | 203 ++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_navier_stokes_convergence.jl diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl new file mode 100644 index 00000000000..21ab09aac06 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -0,0 +1,203 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=1000, Prandtl=0.72, + Mach_freestream=0.5, kappa=1.0) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Tri(), 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); periodicity=(true, false), is_on_boundary) + +# Define initial condition +# Note: If you change the parameters here, also change it in the corresponding source terms +function initial_condition_navier_stokes_convergence_test(x, t, equations) + + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + + y = x[2] + + # 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 + kappa = 1 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = 0.72 + Re = 100 + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one * kappa / Pr + inv_Re = 1.0 / Re + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # 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 + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * inv_Re + + 2.0 / 3.0 * v2_xy * inv_Re + - v1_yy * inv_Re + - v2_xy * inv_Re ) + # 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 + # stress tensor from y-direction + - v1_xy * inv_Re + - v2_xx * inv_Re + - 4.0 / 3.0 * v2_yy * inv_Re + + 2.0 / 3.0 * v1_xy * inv_Re ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * inv_Re + + 2.0 / 3.0 * v2_xy * v1 * inv_Re + - 4.0 / 3.0 * v1_x * v1_x * inv_Re + + 2.0 / 3.0 * v2_y * v1_x * inv_Re + - v1_xy * v2 * inv_Re + - v2_xx * v2 * inv_Re + - v1_y * v2_x * inv_Re + - v2_x * v2_x * inv_Re + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * inv_Re + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * inv_Re + - v2_xy * v1 * inv_Re + - v1_y * v1_y * inv_Re + - v2_x * v1_y * inv_Re + - 4.0 / 3.0 * v2_yy * v2 * inv_Re + + 2.0 / 3.0 * v1_xy * v2 * inv_Re + - 4.0 / 3.0 * v2_y * v2_y * inv_Re + + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * inv_Re ) + + return SVector(du1, du2, du3, du4) +end + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionViscousWall(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 from 0.0 to 1.5 +tspan = (0.0, .50) +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 c103e2ee7c6488f39d67122d35acb0a4e92d8476 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 10:10:07 -0500 Subject: [PATCH 31/71] converging solution using elixir for TreeMesh with BCs --- .../elixir_advection_diffusion_nonperiodic.jl | 88 +++++++++++ src/solvers/dgsem_tree/dg_2d_parabolic.jl | 145 +++++++++++++++++- 2 files changed, 230 insertions(+), 3 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..2dce083330e --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -0,0 +1,88 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +get_diffusivity() = 5.0e-2 +advection_velocity = (1.0, 0.0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) +# Note: If you change the diffusion parameter here, also change it in the initial condition +equations_parabolic = LaplaceDiffusion2D(get_diffusivity(), equations) + +# 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) + +coordinates_min = (-1.0, -0.5) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 0.0, 0.5) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# from "Robust DPG methods for transient convection-diffusion." +# Building bridges: connections and challenges in modern approaches to numerical partial differential equations. +# Springer, Cham, 2016. 179-203. Ellis, Truman, Jesse Chan, and Leszek Demkowicz." +function initial_condition_erikkson_johnson(x, t, equations) + l = 4 + epsilon = get_diffusivity() # TODO: this requires epsilon < .6 due to sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_erikkson_johnson + +# define periodic boundary conditions everywhere +boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition), + y_neg = BoundaryConditionDirichlet(initial_condition), + y_pos = BoundaryConditionDirichlet(initial_condition), + x_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.5 +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +time_int_tol = 1.0e-10 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + +# Print the timer summary +summary_callback() diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 1c332fc6f2e..fabc79ff46e 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -39,7 +39,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::Abstra end - +# TODO: dg_parabolic is not a DG type; it contains solver-specific information such as an LDG penalty parameter. function calc_divergence!(du, u, t, viscous_flux, mesh::TreeMesh{2}, equations_parabolic, boundary_conditions_parabolic, dg::DG, dg_parabolic, @@ -176,7 +176,8 @@ function calc_divergence!(du, u, t, viscous_flux, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - @assert boundary_conditions_parabolic == boundary_condition_periodic + calc_divergence_boundary_flux!(cache_parabolic, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, dg) end # Prolong solution to mortars @@ -250,6 +251,143 @@ function calc_viscous_fluxes!(viscous_flux, u, mesh::TreeMesh{2}, end end +# TODO: decide if we should keep this, and if so, extend to 3D. +function get_unsigned_normal_vector_2d(direction) + if direction > 4 + @warn "Direction = $direction; in 2D, direction should be 1, 2, 3, or 4." + end + if direction==1 || direction==2 + return SVector(1.0, 0.0) + else + return SVector(0.0, 1.0) + end +end + +function calc_gradient_boundary_flux!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{2}, 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_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_gradient_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, 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 i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, 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 + flux_inner = u_inner # TODO: revisit if we want to generalize beyond BR1, LDG. + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, boundary) + flux = boundary_condition(flux_inner, u_inner, get_unsigned_normal_vector_2d(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, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + +function calc_divergence_boundary_flux!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{2}, 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_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_divergence_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, 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` of `calc_divergence!`` + @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 i in eachnode(dg) + # Get viscous boundary fluxes + flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, 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, boundary) + + # TODO: we are passing in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information. + # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion2D and NoSlipWall and Adiabatic + # boundary conditions for CompressibleNavierStokesEquations2D as of 2022-6-27. It may not work with future implementations. + flux = boundary_condition(flux_inner, nothing, get_unsigned_normal_vector_2d(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, direction, neighbor] = flux[v] + end + end + end + + return nothing +end function calc_gradient!(u_grad, u, t, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic, @@ -374,7 +512,8 @@ function calc_gradient!(u_grad, u, t, # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin - @assert boundary_conditions_parabolic == boundary_condition_periodic + calc_gradient_boundary_flux!(cache_parabolic, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, dg) end # Prolong solution to mortars From 3d717a679bc49852b7f62ab41102c4e53ca630ba Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:07:20 -0500 Subject: [PATCH 32/71] fixing DGMulti advection diffusion elixir convergence --- .../elixir_advection_diffusion_nonperiodic.jl | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..086c3a995b7 --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl @@ -0,0 +1,75 @@ +using Trixi, OrdinaryDiffEq + +polydeg = 3 +r1D, w1D = StartUpDG.gauss_lobatto_quad(0, 0, polydeg) +rq, sq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D)) +wq = (x->x[1] .* x[2])(vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D))) + +dg = DGMulti(polydeg = polydeg, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm(); + quad_rule_vol=(rq,sq,wq), quad_rule_face=(r1D,w1D)) + +get_diffusivity() = 5.0e-2 + +equations = LinearScalarAdvectionEquation2D(1.0, 0.0) +equations_parabolic = LaplaceDiffusion2D(get_diffusivity(), equations) + +# from "Robust DPG methods for transient convection-diffusion." +# Building bridges: connections and challenges in modern approaches to numerical partial differential equations. +# Springer, Cham, 2016. 179-203. Ellis, Truman, Jesse Chan, and Leszek Demkowicz." +function initial_condition_erikkson_johnson(x, t, equations) + l = 4 + epsilon = get_diffusivity() # TODO: this requires epsilon < .6 due to sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_erikkson_johnson + +# tag different boundary segments +left(x, tol=50*eps()) = abs(x[1] + 1) < tol +right(x, tol=50*eps()) = abs(x[1]) < tol +bottom(x, tol=50*eps()) = abs(x[2] + .5) < tol +top(x, tol=50*eps()) = abs(x[2] - .5) < tol +entire_boundary(x, tol=50*eps()) = true +is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom, + :entire_boundary => entire_boundary) +mesh = DGMultiMesh(dg; coordinates_min=(-1.0, -0.5), coordinates_max=(0.0, 0.5), + cells_per_dimension=(16, 16), is_on_boundary) + +# BC types +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +# define inviscid boundary conditions +boundary_conditions = (; :left => boundary_condition, + :top => boundary_condition, + :bottom => boundary_condition, + :right => boundary_condition_do_nothing) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :entire_boundary => boundary_condition) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) + +tspan = (0.0, 1.5) +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 8471273b6a96bad5317e5d412bc25ca637970804 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:08:30 -0500 Subject: [PATCH 33/71] naming equations as parabolic/hyperbolic --- src/equations/laplace_diffusion_2d.jl | 38 +++++++++++++-------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index c171741716b..1baa08c8199 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -6,49 +6,49 @@ with diffusivity ``\kappa`` applied to each solution component defined by `equat """ struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N} diffusivity::T - equations::E + equations_hyperbolic::E end -LaplaceDiffusion2D(diffusivity, equations) = - LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations) +LaplaceDiffusion2D(diffusivity, equations_hyperbolic) = + LaplaceDiffusion2D{typeof(equations_hyperbolic), nvariables(equations_hyperbolic), typeof(diffusivity)}(diffusivity, equations_hyperbolic) varnames(variable_mapping, equations_parabolic::LaplaceDiffusion2D) = - varnames(variable_mapping, equations_parabolic.equations) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) # no orientation specified since the flux is vector-valued -function flux(u, grad_u, equations::LaplaceDiffusion2D) +function flux(u, grad_u, equations_parabolic::LaplaceDiffusion2D) dudx, dudy = grad_u - return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy) + return SVector(equations_parabolic.diffusivity * dudx, equations_parabolic.diffusivity * dudy) end # TODO: should this remain in the equations file, be moved to solvers, or live in the elixir? # The penalization depends on the solver, but also depends explicitly on physical parameters, # and would probably need to be specialized for every different equation. -function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) - return dg.penalty_parameter * (u_outer - u_inner) * equations.diffusivity * inv_h +function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity * inv_h end # Dirichlet-type boundary condition for use with a parabolic solver in weak form -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, - equations::LaplaceDiffusion2D) - return boundary_condition.boundary_value_function(x, t, equations) + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_value_function(x, t, equations_parabolic) end -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Divergence, - equations::LaplaceDiffusion2D) - return u_inner + equations_parabolic::LaplaceDiffusion2D) + return flux_inner end -@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Divergence, - equations::LaplaceDiffusion2D) - return boundary_condition.boundary_normal_flux_function(x, t, equations) + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic) end -@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, - equations::LaplaceDiffusion2D) + equations_parabolic::LaplaceDiffusion2D) return flux_inner end From 96667e655695b3b6f3694845ff4ad8941e7f7c06 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:08:51 -0500 Subject: [PATCH 34/71] generalizing transform_variables --- src/solvers/dgmulti/dg_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 3ad91627f63..a9279dde599 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -41,7 +41,7 @@ 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, equations) +function transform_variables!(u_transformed, u, extra_args...) @threaded for i in eachindex(u) u_transformed[i] = u[i] end From 73e4740c50c9940e926ef7bd8fbc365bc1b10785 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:10:32 -0500 Subject: [PATCH 35/71] add TODO more todos --- src/equations/compressible_navier_stokes_2d.jl | 5 ++++- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 14 ++++++++++---- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index ea360b11f19..68337c860ca 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -41,11 +41,14 @@ or the entropy variables where w_2 = rho v_1 / p, w_3 = rho v_2 / p, w_4 = -rho / p """ + # TODO: # 1) For now I save gamma and inv(gamma-1) again, but we could potentially reuse them from # the Euler equations # 2) Add more here and probably some equations -struct CompressibleNavierStokesEquations2D{RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesEquations{2, 3} + +# TODO: Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function +struct CompressibleNavierStokesEquations2D{RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesEquations{2, 4} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications Re::RealT # Reynolds number diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index fabc79ff46e..3d92cfa5744 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -307,7 +307,11 @@ function calc_gradient_boundary_flux_by_direction!(surface_flux_values::Abstract else # Element is on the right, boundary on the left u_inner = u_rr end - flux_inner = u_inner # TODO: revisit if we want to generalize beyond BR1, LDG. + + # 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, boundary) flux = boundary_condition(flux_inner, u_inner, get_unsigned_normal_vector_2d(direction), @@ -373,9 +377,11 @@ function calc_divergence_boundary_flux_by_direction!(surface_flux_values::Abstra x = get_node_coords(node_coordinates, equations_parabolic, dg, i, boundary) - # TODO: we are passing in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information. - # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion2D and NoSlipWall and Adiabatic - # boundary conditions for CompressibleNavierStokesEquations2D as of 2022-6-27. It may not work with future implementations. + # TODO: add a field in `cache.boundaries` for gradient information. + # 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 LaplaceDiffusion2D and + # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesEquations2D 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_2d(direction), x, t, Divergence(), equations_parabolic) From 2c05f6645fbb39a8b3fb046fb3b6de0695ea3b7f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:10:15 -0500 Subject: [PATCH 36/71] additional checks on get_unsigned_normal --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 3d92cfa5744..89461342ca5 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -253,7 +253,7 @@ end # TODO: decide if we should keep this, and if so, extend to 3D. function get_unsigned_normal_vector_2d(direction) - if direction > 4 + if direction > 4 || direction < 1 @warn "Direction = $direction; in 2D, direction should be 1, 2, 3, or 4." end if direction==1 || direction==2 From d5e7f9e9f30ad122dd1d329c0f2d77cfa95ac921 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:11:21 -0500 Subject: [PATCH 37/71] adding isothermal BC --- src/Trixi.jl | 2 +- .../compressible_navier_stokes_2d.jl | 30 +++++++++++++++---- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 76482b841c1..75507ec1fc4 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -164,7 +164,7 @@ export boundary_condition_do_nothing, boundary_condition_noslip_wall, boundary_condition_slip_wall, boundary_condition_wall, - BoundaryConditionViscousWall, NoSlip, Adiabatic + BoundaryConditionViscousWall, NoSlip, Adiabatic, Isothermal export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 68337c860ca..1078c759d7e 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -230,25 +230,45 @@ struct NoSlip{F} boundary_value_function::F # value of the velocity vector on the boundary end +# isothermal temperature BC +struct Isothermal{F} + boundary_value_function::F # value of the temperature on the boundary +end + # adiabatic temperature BC struct Adiabatic{F} boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn end -@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(u_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, equations::CompressibleNavierStokesEquations2D) - rho = u_inner[1] # should not matter since the viscous terms don't depend on rho v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) - return SVector(rho, v1, v2, u_inner[4]) + return SVector(zero(eltype(u_inner)), v1, v2, u_inner[4]) end -@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(flux_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Divergence, equations::CompressibleNavierStokesEquations2D) + # rho, v1, v2, _ = u_inner normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) - tau_1n, tau_2n = flux_inner[2:3] + tau_1n, tau_2n = flux_inner[2:3] # extract fluxes for 2nd and 3rd equations normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux return SVector(flux_inner[1:3]..., normal_energy_flux) end + +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesEquations2D) + v1, v2 = 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(zero(eltype(flux_inner)), v1, v2, T) +end + +@inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesEquations2D) + return flux_inner +end + From 89acae836971f9aa9819b0be78c55d9a8e775c58 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:11:36 -0500 Subject: [PATCH 38/71] commenting out unused CNS functions for now --- .../compressible_navier_stokes_2d.jl | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 1078c759d7e..48477e754d0 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -175,40 +175,40 @@ end -# Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleNavierStokesEquations2D) - rho, rho_v1, rho_v2, rho_e = u +# # Convert conservative variables to entropy +# @inline function cons2entropy(u, equations::CompressibleNavierStokesEquations2D) +# rho, rho_v1, rho_v2, rho_e = u - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) +# v1 = rho_v1 / rho +# v2 = rho_v2 / rho +# v_square = v1^2 + v2^2 +# p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) - rho_p = rho / p +# rho_p = rho / p - w2 = rho_p * v1 - w3 = rho_p * v2 - w4 = -rho_p +# w2 = rho_p * v1 +# w3 = rho_p * v2 +# w4 = -rho_p - return SVector(w2, w3, w4) -end +# return SVector(w2, w3, w4) +# end -@inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNavierStokesEquations2D) -# Takes the solution values `u` and gradient of the variables (w_2, w_3, w_4) and -# reverse engineers the gradients to be terms of the primitive vairables (v1, v2, T). -# Helpful because then the diffusive fluxes have the same form as on paper. - rho, rho_v1, rho_v2, _ = u +# @inline function convert_gradient_variables(u, grad_entropy_vars, equations::CompressibleNavierStokesEquations2D) +# # Takes the solution values `u` and gradient of the 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. +# rho, rho_v1, rho_v2, _ = u - v1 = rho_v1 / rho - v2 = rho_v2 / rho - T = temperature(u, equations) +# v1 = rho_v1 / rho +# v2 = rho_v2 / rho +# T = temperature(u, equations) - return SVector(equations.R * T * (grad_entropy_vars[1] + v1 * grad_entropy_vars[3]), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) - equations.R * T * (grad_entropy_vars[2] + v2 * grad_entropy_vars[3]), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) - equations.R * T * T * grad_entropy_vars[3] # grad(T) = R*T^2*grad(w_4)) - ) -end +# return SVector(equations.R * T * (grad_entropy_vars[1] + v1 * grad_entropy_vars[3]), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) +# equations.R * T * (grad_entropy_vars[2] + v2 * grad_entropy_vars[3]), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) +# equations.R * T * T * grad_entropy_vars[3] # grad(T) = R*T^2*grad(w_4)) +# ) +# end @inline function temperature(u, equations::CompressibleNavierStokesEquations2D) From f503b85b4693f51e3c79a75de710a3bd6efe927e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:12:46 -0500 Subject: [PATCH 39/71] fix call to transform_variables --- src/solvers/dgmulti/dg_parabolic.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index a9279dde599..be6fc066f73 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -101,7 +101,7 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, scalar_flux_face_values[idM] = 0.5 * (uP + uM) # TODO: use strong/weak formulation for curved meshes? end - calc_boundary_flux!(scalar_flux_face_values, nothing, t, Gradient(), boundary_conditions, + calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, mesh, equations, dg, cache, cache_parabolic) # compute surface contributions @@ -158,10 +158,9 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, # for both the gradient and the divergence, the boundary flux is scalar valued. # for the gradient, it is the solution; for divergence, it is the normal flux. - u_boundary = boundary_condition(flux_face_values[fid,e], - face_normal, face_coordinates, t, - operator_type, equations) - flux_face_values[fid,e] = u_boundary + flux_face_values[fid,e] = boundary_condition(flux_face_values[fid,e], u_face_values[fid,e], + face_normal, face_coordinates, t, + operator_type, equations) end end return nothing @@ -299,7 +298,7 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::Abstra reset_du!(du, dg) @unpack u_transformed, u_grad, viscous_flux = cache_parabolic - transform_variables!(u_transformed, u, equations_parabolic) + transform_variables!(u_transformed, u, mesh, dg, dg_parabolic, equations_parabolic) calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, cache_parabolic) From 41b41a9f37bb7c1f521dab7b9fbbf9734a41d9dc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:17:08 -0500 Subject: [PATCH 40/71] comments and cleanup --- src/solvers/dgmulti/dg_parabolic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index be6fc066f73..54f8c6e9924 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -265,14 +265,14 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh for dim in eachdim(mesh) uM = viscous_flux_face_values[dim][idM] uP = viscous_flux_face_values[dim][idP] - # TODO: use strong/weak formulation? + # TODO: use strong/weak formulation to ensure stability on curved meshes? flux_face_value = flux_face_value + 0.5 * (uP + uM) * nxyzJ[dim][face_node_index] end scalar_flux_face_values[idM] = flux_face_value end # TODO: decide what to pass in - calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(), + calc_boundary_flux!(scalar_flux_face_values, cache_parabolic.u_face_values, t, Divergence(), boundary_conditions, mesh, equations, dg, cache, cache_parabolic) calc_viscous_penalty!(scalar_flux_face_values, cache_parabolic.u_face_values, t, From 1107def502cf1f39233a1105b5a1c96c8b39011a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 15:58:20 -0500 Subject: [PATCH 41/71] changing default solver and Re for cavity --- .../dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl index 63bb1ce29f8..d93cff55c5f 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl @@ -7,13 +7,13 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=100, Prandtl=0.72, +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=1000, Prandtl=0.72, Mach_freestream=0.1, kappa=1.0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = Polynomial(), +dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(), surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), - volume_integral = VolumeIntegralWeakForm()) + volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) top(x, tol=50*eps()) = abs(x[2] - 1) < tol rest_of_boundary(x, tol=50*eps()) = !top(x, tol) @@ -52,7 +52,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers, callbacks etc. # Create ODE problem with time span from 0.0 to 1.5 -tspan = (0.0, 25.0) +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan); summary_callback = SummaryCallback() From f90afb97e3d783404998f355bca9ff26d794b749 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 16:47:37 -0500 Subject: [PATCH 42/71] adding more advection diffusion tests --- test/test_parabolic_2d.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 4a792af4840..8ac13d581ce 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -78,11 +78,19 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "elixir_advection_diffusion.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion.jl"), + @trixi_testset "elixir_advection_diffusion_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_nonperiodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), - l2 = [0.2485803335154642], - linf = [1.079606969242132] + l2 = [0.012561479036088107], + linf = [0.10067620068384618] + ) + end + + @trixi_testset "elixir_advection_diffusion_nonperiodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.0076468006081234705], + linf = [0.10067621027072088] ) end From 2e7019feb07b8a2b064c2151c3bb757df3c13aec Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 19:38:39 -0500 Subject: [PATCH 43/71] label tests --- test/test_parabolic_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 8ac13d581ce..a9d208b8fbc 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -70,7 +70,7 @@ isdir(outdir) && rm(outdir, recursive=true) @test getindex.(du, 1) ≈ 2 * y end - @trixi_testset "elixir_advection_diffusion_periodic.jl" begin + @trixi_testset "DGMulti: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_periodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.03180371984888462], @@ -78,7 +78,7 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "elixir_advection_diffusion_nonperiodic.jl" begin + @trixi_testset "DGMulti: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_nonperiodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.012561479036088107], @@ -86,7 +86,7 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "elixir_advection_diffusion_nonperiodic.jl" begin + @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), l2 = [0.0076468006081234705], From 063b602ebdb06a09be0f1cb34d5f24cdc86c6a83 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 19:39:27 -0500 Subject: [PATCH 44/71] add gradient_variables field to SemidiscretizationHyperbolicParabolic --- ...semidiscretization_hyperbolic_parabolic.jl | 33 +++++++++++++------ 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 52a48d7aa03..4c960066cf6 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -11,7 +11,7 @@ A struct containing everything needed to describe a spatial semidiscretization of a mixed hyperbolic-parabolic conservation law. """ -struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, +struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} <: AbstractSemidiscretization @@ -20,6 +20,8 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic equations::Equations equations_parabolic::EquationsParabolic + gradient_variables::GradientVariables + # This guy is a bit messy since we abuse it as some kind of "exact solution" # although this doesn't really exist... initial_condition::InitialCondition @@ -37,17 +39,18 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic performance_counter::PerformanceCounter - function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}( - mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, + function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}( + mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, + gradient_variables::GradientVariables, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic, - source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} + source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} @assert ndims(mesh) == ndims(equations) # Todo: assert nvariables(equations)==nvariables(equations_parabolic) performance_counter = PerformanceCounter() - new(mesh, equations, equations_parabolic, initial_condition, + new(mesh, equations, equations_parabolic, gradient_variables, initial_condition, boundary_conditions, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic, performance_counter) end @@ -55,6 +58,7 @@ end """ SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver; + gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), @@ -66,6 +70,7 @@ Construct a semidiscretization of a hyperbolic-parabolic PDE. """ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_condition, solver; + gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), @@ -79,7 +84,8 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_hyperbolic_cache, initial_cache_parabolic = initial_caches return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic, - initial_condition, solver; solver_parabolic, source_terms, + initial_condition, solver; + gradient_variables, solver_parabolic, source_terms, boundary_conditions=boundary_conditions_hyperbolic, boundary_conditions_parabolic=boundary_conditions_parabolic, RealT, uEltype, initial_cache=initial_hyperbolic_cache, @@ -88,6 +94,7 @@ end function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, solver; + gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, boundary_conditions=boundary_condition_periodic, @@ -106,10 +113,10 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo solver, solver_parabolic, RealT, uEltype)..., initial_cache_parabolic...) - SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), + SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(gradient_variables), typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic), typeof(source_terms), typeof(solver), typeof(solver_parabolic), typeof(cache), typeof(cache_parabolic)}( - mesh, equations, equations_parabolic, initial_condition, + mesh, equations, equations_parabolic, gradient_variables, initial_condition, _boundary_conditions, _boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic) end @@ -124,6 +131,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s mesh=semi.mesh, equations=semi.equations, equations_parabolic=semi.equations_parabolic, + gradient_variables=semi.gradient_variables, initial_condition=semi.initial_condition, solver=semi.solver, solver_parabolic=semi.solver_parabolic, @@ -135,7 +143,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s # special care if shock-capturing volume integrals are used (because of # the indicators and their own caches...). SemidiscretizationHyperbolicParabolic( - mesh, equations, equations_parabolic, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) + mesh, equations, equations_parabolic, gradient_variables, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) end function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) @@ -145,6 +153,7 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) print(io, semi.mesh) print(io, ", ", semi.equations) print(io, ", ", semi.equations_parabolic) + print(io, ", ", semi.gradient_variables) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) print(io, ", ", semi.boundary_conditions_parabolic) @@ -170,6 +179,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_line(io, "mesh", semi.mesh) summary_line(io, "hyperbolic equations", semi.equations |> typeof |> nameof) summary_line(io, "parabolic equations", semi.equations_parabolic |> typeof |> nameof) + summary_line(io, "gradient variables", semi.gradient_variables |> typeof |> nameof) summary_line(io, "initial condition", semi.initial_condition) # print_boundary_conditions(io, semi) @@ -243,8 +253,11 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi + @unpack mesh, equations, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi + # TODO: clean up. We wrap_array using `equations` and not `equations_parabolic`. + # This is because `u` is a hyperbolic solution and `equations_parabolic` can have a + # different number of variables for something like CompressibleNavierStokesEquations2D u = wrap_array(u_ode, mesh, equations_parabolic, solver, cache_parabolic) du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic) From 50d76df032bbf9c3c7614e81ebb823a4968a40be Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 27 Jun 2022 19:39:27 -0500 Subject: [PATCH 45/71] Revert "add gradient_variables field to SemidiscretizationHyperbolicParabolic" This reverts commit 063b602ebdb06a09be0f1cb34d5f24cdc86c6a83. --- ...semidiscretization_hyperbolic_parabolic.jl | 33 ++++++------------- 1 file changed, 10 insertions(+), 23 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 4c960066cf6..52a48d7aa03 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -11,7 +11,7 @@ A struct containing everything needed to describe a spatial semidiscretization of a mixed hyperbolic-parabolic conservation law. """ -struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, +struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} <: AbstractSemidiscretization @@ -20,8 +20,6 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic equations::Equations equations_parabolic::EquationsParabolic - gradient_variables::GradientVariables - # This guy is a bit messy since we abuse it as some kind of "exact solution" # although this doesn't really exist... initial_condition::InitialCondition @@ -39,18 +37,17 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic performance_counter::PerformanceCounter - function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}( - mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, - gradient_variables::GradientVariables, initial_condition::InitialCondition, + function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}( + mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic, - source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, GradientVariables, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} + source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} @assert ndims(mesh) == ndims(equations) # Todo: assert nvariables(equations)==nvariables(equations_parabolic) performance_counter = PerformanceCounter() - new(mesh, equations, equations_parabolic, gradient_variables, initial_condition, + new(mesh, equations, equations_parabolic, initial_condition, boundary_conditions, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic, performance_counter) end @@ -58,7 +55,6 @@ end """ SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver; - gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), @@ -70,7 +66,6 @@ Construct a semidiscretization of a hyperbolic-parabolic PDE. """ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_condition, solver; - gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), @@ -84,8 +79,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_hyperbolic_cache, initial_cache_parabolic = initial_caches return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic, - initial_condition, solver; - gradient_variables, solver_parabolic, source_terms, + initial_condition, solver; solver_parabolic, source_terms, boundary_conditions=boundary_conditions_hyperbolic, boundary_conditions_parabolic=boundary_conditions_parabolic, RealT, uEltype, initial_cache=initial_hyperbolic_cache, @@ -94,7 +88,6 @@ end function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, solver; - gradient_variables=cons2cons, solver_parabolic=default_parabolic_solver(), source_terms=nothing, boundary_conditions=boundary_condition_periodic, @@ -113,10 +106,10 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo solver, solver_parabolic, RealT, uEltype)..., initial_cache_parabolic...) - SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(gradient_variables), + SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic), typeof(source_terms), typeof(solver), typeof(solver_parabolic), typeof(cache), typeof(cache_parabolic)}( - mesh, equations, equations_parabolic, gradient_variables, initial_condition, + mesh, equations, equations_parabolic, initial_condition, _boundary_conditions, _boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic) end @@ -131,7 +124,6 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s mesh=semi.mesh, equations=semi.equations, equations_parabolic=semi.equations_parabolic, - gradient_variables=semi.gradient_variables, initial_condition=semi.initial_condition, solver=semi.solver, solver_parabolic=semi.solver_parabolic, @@ -143,7 +135,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s # special care if shock-capturing volume integrals are used (because of # the indicators and their own caches...). SemidiscretizationHyperbolicParabolic( - mesh, equations, equations_parabolic, gradient_variables, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) + mesh, equations, equations_parabolic, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) end function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) @@ -153,7 +145,6 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) print(io, semi.mesh) print(io, ", ", semi.equations) print(io, ", ", semi.equations_parabolic) - print(io, ", ", semi.gradient_variables) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) print(io, ", ", semi.boundary_conditions_parabolic) @@ -179,7 +170,6 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_line(io, "mesh", semi.mesh) summary_line(io, "hyperbolic equations", semi.equations |> typeof |> nameof) summary_line(io, "parabolic equations", semi.equations_parabolic |> typeof |> nameof) - summary_line(io, "gradient variables", semi.gradient_variables |> typeof |> nameof) summary_line(io, "initial condition", semi.initial_condition) # print_boundary_conditions(io, semi) @@ -253,11 +243,8 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, equations, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi + @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi - # TODO: clean up. We wrap_array using `equations` and not `equations_parabolic`. - # This is because `u` is a hyperbolic solution and `equations_parabolic` can have a - # different number of variables for something like CompressibleNavierStokesEquations2D u = wrap_array(u_ode, mesh, equations_parabolic, solver, cache_parabolic) du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic) From 549793ca9dbe0aaaddfab42e2a1c5816ec9e07fe Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 12:12:17 -0500 Subject: [PATCH 46/71] allowing for specialization in transform_variables adding a function `gradient_variable_transformation` which should get specialized if the gradient variables are not conservative variables --- .../compressible_navier_stokes_2d.jl | 13 +++--------- src/equations/equations_parabolic.jl | 4 ++++ src/solvers/dgmulti/dg_parabolic.jl | 8 +++++--- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 20 ++++++++++++++++++- 4 files changed, 31 insertions(+), 14 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 48477e754d0..217c3810d60 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -80,9 +80,6 @@ function CompressibleNavierStokesEquations2D(equations::CompressibleEulerEquatio end -# I was not sure what to do here to allow flexibility of selecting primitive or entropy -# gradient variables. I see that `transform_variables!` just copies data at the moment. - # This is the flexibility a user should have to select the different gradient variable types # varnames(::typeof(cons2prim) , ::CompressibleNavierStokesEquations2D) = ("v1", "v2", "T") # varnames(::typeof(cons2entropy), ::CompressibleNavierStokesEquations2D) = ("w2", "w3", "w4") @@ -90,13 +87,9 @@ end varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesEquations2D) = varnames(variable_mapping, equations_parabolic.equations_hyperbolic) -# transform to primitive variables -# TODO: should we have this call a "gradient transformation" field? -function transform_variables!(u_transformed, u, equations_parabolic::CompressibleNavierStokesEquations2D) - @threaded for i in eachindex(u) - u_transformed[i] = cons2prim(u[i], equations_parabolic) - end -end +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +gradient_variable_transformation(::CompressibleNavierStokesEquations2D, dg_parabolic) = cons2prim # no orientation specified since the flux is vector-valued # Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2 diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index df5f7077ad1..24f2b51503c 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -1,3 +1,7 @@ +# specify transformation of conservative variables prior to taking gradients. +# specialize this function to compute gradients e.g., of primitive variables instead of conservative +gradient_variable_transformation(::AbstractEquationsParabolic, dg_parabolic) = cons2cons + # Linear scalar diffusion for use in linear scalar advection-diffusion problems abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_2d.jl") diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 54f8c6e9924..41a190bdc1a 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -41,9 +41,10 @@ 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, extra_args...) +function transform_variables!(u_transformed, u, mesh, equations_parabolic::AbstractEquationsParabolic, + dg::DGMulti, dg_parabolic, cache, cache_parabolic) @threaded for i in eachindex(u) - u_transformed[i] = u[i] + u_transformed[i] = gradient_variable_transformation(equations_parabolic, dg_parabolic)(u[i], equations_parabolic) end end @@ -298,7 +299,8 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::Abstra reset_du!(du, dg) @unpack u_transformed, u_grad, viscous_flux = cache_parabolic - transform_variables!(u_transformed, u, mesh, dg, dg_parabolic, equations_parabolic) + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, dg_parabolic, cache, cache_parabolic) calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, cache_parabolic) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 89461342ca5..1896a5c1a9d 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -16,7 +16,9 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::Abstra @unpack u_transformed, u_grad = cache_parabolic @trixi_timeit timer() "transform variables" transform_variables!(u_transformed, u, - equations_parabolic) + mesh, equations_parabolic, + dg, dg_parabolic, + cache, cache_parabolic) @trixi_timeit timer() "calculate gradient" calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, @@ -39,6 +41,22 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::Abstra 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{2}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, dg_parabolic, cache, cache_parabolic) + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations_parabolic, dg, i, j, element) + u_transformed_node = gradient_variable_transformation(equations_parabolic, dg_parabolic)(u_node, equations_parabolic) + set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, element) + end + end +end + # TODO: dg_parabolic is not a DG type; it contains solver-specific information such as an LDG penalty parameter. function calc_divergence!(du, u, t, viscous_flux, mesh::TreeMesh{2}, equations_parabolic, From 3e3e41f2a1eab1845432997feb1280729f0f4c88 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:47:57 -0500 Subject: [PATCH 47/71] formatting and comments --- .../dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl index 086c3a995b7..3a263aa0ad1 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl @@ -20,7 +20,7 @@ equations_parabolic = LaplaceDiffusion2D(get_diffusivity(), equations) # Springer, Cham, 2016. 179-203. Ellis, Truman, Jesse Chan, and Leszek Demkowicz." function initial_condition_erikkson_johnson(x, t, equations) l = 4 - epsilon = get_diffusivity() # TODO: this requires epsilon < .6 due to sqrt + epsilon = get_diffusivity() # TODO: this requires epsilon < .6 due to the sqrt lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) @@ -46,10 +46,10 @@ mesh = DGMultiMesh(dg; coordinates_min=(-1.0, -0.5), coordinates_max=(0.0, 0.5), boundary_condition = BoundaryConditionDirichlet(initial_condition) # define inviscid boundary conditions -boundary_conditions = (; :left => boundary_condition, - :top => boundary_condition, +boundary_conditions = (; :left => boundary_condition, + :top => boundary_condition, :bottom => boundary_condition, - :right => boundary_condition_do_nothing) + :right => boundary_condition_do_nothing) # define viscous boundary conditions boundary_conditions_parabolic = (; :entire_boundary => boundary_condition) From 63471512ef9568abe98a11d2031cbe8e30cd03ce Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:48:04 -0500 Subject: [PATCH 48/71] reverting elixir --- .../elixir_advection_diffusion_nonperiodic.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl index 3a263aa0ad1..a8e9e3329af 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl @@ -1,14 +1,8 @@ using Trixi, OrdinaryDiffEq -polydeg = 3 -r1D, w1D = StartUpDG.gauss_lobatto_quad(0, 0, polydeg) -rq, sq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D)) -wq = (x->x[1] .* x[2])(vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D))) - -dg = DGMulti(polydeg = polydeg, element_type = Quad(), approximation_type = Polynomial(), +dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), - volume_integral = VolumeIntegralWeakForm(); - quad_rule_vol=(rq,sq,wq), quad_rule_face=(r1D,w1D)) + volume_integral = VolumeIntegralWeakForm()) get_diffusivity() = 5.0e-2 From eb98568296d02ec14e864c0b65b1add7a2b1e453 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:48:17 -0500 Subject: [PATCH 49/71] comments --- examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index 2dce083330e..c49ca70aadf 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -38,7 +38,6 @@ function initial_condition_erikkson_johnson(x, t, equations) end initial_condition = initial_condition_erikkson_johnson -# define periodic boundary conditions everywhere boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition), y_neg = BoundaryConditionDirichlet(initial_condition), y_pos = BoundaryConditionDirichlet(initial_condition), From 766cdc5b8e703937be7d2eda9b1a682563278dd0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:48:35 -0500 Subject: [PATCH 50/71] standardizing time tol --- .../tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index c49ca70aadf..c55459ae56b 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -79,7 +79,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -time_int_tol = 1.0e-10 +time_int_tol = 1.0e-8 sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false, callback=callbacks) From c865fb78d9c05df01f507cb5677ccc79c90bd94a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:49:23 -0500 Subject: [PATCH 51/71] minor fix to CNS boundary flux for convenience make it so that the density state is computed correctly, even though it's not used --- 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 217c3810d60..777300779e5 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -237,7 +237,7 @@ end x, t, operator_type::Gradient, equations::CompressibleNavierStokesEquations2D) v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) - return SVector(zero(eltype(u_inner)), v1, v2, u_inner[4]) + return SVector(u_inner[1], v1, v2, u_inner[4]) end @inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, From 08768a7ad29fc0b2ed3f8b7bb4eccdc7fc0cf6d5 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:49:33 -0500 Subject: [PATCH 52/71] formatting + comments --- src/equations/compressible_navier_stokes_2d.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 777300779e5..8527f5bbf9e 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -116,16 +116,16 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) # Components of viscous stress tensor - # (4/3*(v1)_x - 2/3*(v2)_y) + # (4/3 * (v1)_x - 2/3 * (v2)_y) tau_11 = ( 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy ) # ((v1)_y + (v2)_x) # stress tensor is symmetric tau_12 = ( dv1dy + dv2dx ) # = tau_21 - # (4/3*(v2)_y - 2/3*(v1)_x) + # (4/3 * (v2)_y - 2/3 * (v1)_x) tau_22 = ( 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx ) - # Fick's law q = -kappa*grad(T); constant is kappa*gamma/(Pr*(gamma-1)) - # Important note! Due to nondimensional scaling R = 1/gamma, so the + # Fick's law q = -kappa * grad(T); constant is kappa * gamma / (Pr * (gamma-1)) + # Important note! Due to nondimensional scaling R = 1 / gamma, so the # temperature T in the gradient computation already contains a factor of gamma q1 = ( equations.kappa * equations.inv_gamma_minus_one * dTdx ) / equations.Pr q2 = ( equations.kappa * equations.inv_gamma_minus_one * dTdy ) / equations.Pr @@ -144,7 +144,6 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) g1 = zero(rho) g2 = f3 # tau_21 * nu g3 = tau_22 * nu - # g4 = ( v1 * tau_21 + v2 * tau_22 + q2 ) * nu g4 = ( v1 * tau_12 + v2 * tau_22 + q2 ) * nu return (SVector(f1, f2, f3, f4) , SVector(g1, g2, g3, g4)) From 6c38bdbba6e120405acafebc10a58ab39af326ea Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:50:02 -0500 Subject: [PATCH 53/71] using primitive variables in viscous flux instead of conservative --- src/equations/compressible_navier_stokes_2d.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 8527f5bbf9e..7161e1f7783 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -105,10 +105,7 @@ function flux(u, grad_u, equations::CompressibleNavierStokesEquations2D) # Here grad_u is assumed to contain the gradients of the primitive variables (v1, v2, T) # either computed directly or reverse engineered from the gradient of the entropy vairables # by way of the `convert_gradient_variables` function - rho, rho_v1, rho_v2, _ = u - - v1 = rho_v1 / rho - v2 = rho_v2 / rho + rho, v1, v2, _ = u # grad_u contains derivatives of each hyperbolic variable _, dv1dx, dv2dx, dTdx = grad_u[1] From 0f67dd07bd902cab803583812c91f834cdf2e423 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 22:50:49 -0500 Subject: [PATCH 54/71] minor formatting --- .../dgmulti_2d/elixir_navier_stokes_convergence.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl index 21ab09aac06..a72ab18ab1b 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -4,10 +4,13 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations +get_Re() = 100 +get_Pr() = .72 + equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=1000, Prandtl=0.72, +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=get_Re(), Prandtl=get_Pr(), Mach_freestream=0.5, kappa=1.0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -50,8 +53,8 @@ initial_condition = initial_condition_navier_stokes_convergence_test # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 kappa = 1 inv_gamma_minus_one = inv(equations.gamma - 1) - Pr = 0.72 - Re = 100 + Pr = get_Pr() + Re = get_Re() # Same settings as in `initial_condition` # Amplitude and shift @@ -186,7 +189,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # Create ODE problem with time span from 0.0 to 1.5 tspan = (0.0, .50) -ode = semidiscretize(semi, tspan); +ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() alive_callback = AliveCallback(alive_interval=10) From 6a1d908bdba547e3d7e7b754f585cfe1e0c6491a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 23:27:15 -0500 Subject: [PATCH 55/71] add CNS tests --- test/test_parabolic_2d.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index a9d208b8fbc..2aae2ed279e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -71,7 +71,7 @@ isdir(outdir) && rm(outdir, recursive=true) end @trixi_testset "DGMulti: elixir_advection_diffusion_periodic.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_periodic.jl"), + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_periodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.03180371984888462], linf = [0.2136821621370909] @@ -79,13 +79,21 @@ isdir(outdir) && rm(outdir, recursive=true) end @trixi_testset "DGMulti: elixir_advection_diffusion_nonperiodic.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_nonperiodic.jl"), + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_nonperiodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.012561479036088107], linf = [0.10067620068384618] ) end + @trixi_testset "DGMulti: elixir_navier_stokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_convergence.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.00153550768125133, 0.0033843168272696357, 0.0036531858107444067, 0.009948436427519428], + linf = [0.005522560467190019, 0.013425258500731063, 0.013962115643483375, 0.027483102120516634] + ) + end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), From 6d5be6c52d2b65df8e402d4b1e57909cd273db08 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 12 Jul 2022 23:49:50 -0500 Subject: [PATCH 56/71] fix test --- test/test_parabolic_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 2aae2ed279e..94ab0235616 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -81,8 +81,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "DGMulti: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_nonperiodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), - l2 = [0.012561479036088107], - linf = [0.10067620068384618] + l2 = [0.002123168335604323], + linf = [0.00963640423513712] ) end From bc58754f578f2497c6b905af9756ba8611113606 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 13 Jul 2022 08:27:50 -0500 Subject: [PATCH 57/71] testing if scoping issues fix TreeMesh tests --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 1896a5c1a9d..e4ce2a3a775 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -470,7 +470,7 @@ function calc_gradient!(u_grad, u, t, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @unpack surface_flux_values = cache_parabolic.elements - @unpack u, neighbor_ids, orientations = cache_parabolic.interfaces + @unpack neighbor_ids, orientations = cache_parabolic.interfaces @threaded for interface in eachinterface(dg, cache_parabolic) # Get neighboring elements @@ -485,7 +485,8 @@ function calc_gradient!(u_grad, u, t, for i in eachnode(dg) # Call pointwise Riemann solver - u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, interface) + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, + equations_parabolic, dg, i, interface) flux = 0.5 * (u_ll + u_rr) # Copy flux to left and right element storage From f13045a8b7c637be1b9c8a7cc5856cdde6e4fa6f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 13 Jul 2022 11:31:25 -0500 Subject: [PATCH 58/71] decreasing timestep tol for TreeMesh parabolic test --- .../tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl | 2 +- test/test_parabolic_2d.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index c55459ae56b..74d0e83bdf0 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -79,7 +79,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -time_int_tol = 1.0e-8 +time_int_tol = 1.0e-11 sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false, callback=callbacks) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 94ab0235616..929072c9376 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -97,8 +97,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), - l2 = [0.0076468006081234705], - linf = [0.10067621027072088] + l2 = [0.007646800618485118], + linf = [0.10067621050468958] ) end From d29d583b1f75acaf0a8a9443a968f6fd3f65612f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 13 Jul 2022 16:05:07 -0500 Subject: [PATCH 59/71] enabling periodic BCs for TreeMesh + more tests --- .../elixir_advection_diffusion.jl | 2 +- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 18 ++++++++++++++++-- test/test_parabolic_2d.jl | 8 ++++++++ 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl index 1706d4d4184..d9ea702ea7c 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl @@ -77,7 +77,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -time_int_tol = 1.0e-8 +time_int_tol = 1.0e-11 sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false, callback=callbacks) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index e4ce2a3a775..a126da4afe9 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -122,7 +122,7 @@ function calc_divergence!(du, u, t, viscous_flux, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @unpack surface_flux_values = cache_parabolic.elements - @unpack u, neighbor_ids, orientations = cache_parabolic.interfaces + @unpack neighbor_ids, orientations = cache_parabolic.interfaces @threaded for interface in eachinterface(dg, cache_parabolic) # Get neighboring elements @@ -137,7 +137,8 @@ function calc_divergence!(du, u, t, viscous_flux, for i in eachnode(dg) # Call pointwise Riemann solver - u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, interface) + # TODO: should this u be cache_parabolic.interfaces.u?? + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, i, interface) flux = 0.5 * (u_ll + u_rr) # Copy flux to left and right element storage @@ -281,6 +282,18 @@ function get_unsigned_normal_vector_2d(direction) end end +function calc_gradient_boundary_flux!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_divergence_boundary_flux!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + function calc_gradient_boundary_flux!(cache, t, boundary_conditions_parabolic::NamedTuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG) @@ -485,6 +498,7 @@ function calc_gradient!(u_grad, u, t, for i in eachnode(dg) # Call pointwise Riemann solver + # TODO: should this be cache_parabolic.interfaces.u?? u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, i, interface) flux = 0.5 * (u_ll + u_rr) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 929072c9376..79084b35878 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -94,6 +94,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.4), polydeg=5, + l2 = [4.0915532997994255e-6], + linf = [2.3040850347877395e-5] + ) + end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), From abd020b943c447b3092ee2ac21fb30762e0909dd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 13 Jul 2022 16:05:26 -0500 Subject: [PATCH 60/71] fix density BC flux (unused, but could be useful) --- src/equations/compressible_navier_stokes_2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 7161e1f7783..1ae4ba87ced 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -163,7 +163,6 @@ end prim2cons(u, equations.equations_hyperbolic) - # # Convert conservative variables to entropy # @inline function cons2entropy(u, equations::CompressibleNavierStokesEquations2D) # rho, rho_v1, rho_v2, rho_e = u @@ -252,7 +251,7 @@ end equations::CompressibleNavierStokesEquations2D) v1, v2 = 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(zero(eltype(flux_inner)), v1, v2, T) + return SVector(u_inner[1], v1, v2, T) end @inline function (boundary_condition::BoundaryConditionViscousWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, From 663f4657ceb14ca7ddbccbfba3af12ffefea66bf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 15 Jul 2022 11:38:35 -0500 Subject: [PATCH 61/71] adding non-working TreeMesh elixirs --- .../elixir_navier_stokes_convergence.jl | 216 ++++++++++++++++++ .../elixir_navier_stokes_lid_driven_cavity.jl | 77 +++++++ 2 files changed, 293 insertions(+) create mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl new file mode 100644 index 00000000000..8909578982a --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -0,0 +1,216 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +get_Re() = 1000 +get_Pr() = .72 + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=get_Re(), Prandtl=get_Pr(), + Mach_freestream=0.5, kappa=1.0) + +# 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) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=(true, false), + n_cells_max=30_000) # set maximum capacity of tree data structure + + +# Define initial condition +# Note: If you change the parameters here, also change it in the corresponding source terms +function initial_condition_navier_stokes_convergence_test(x, t, equations) + + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + + y = x[2] + + # 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 + kappa = 1 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = get_Pr() + Re = get_Re() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one * kappa / Pr + inv_Re = 1.0 / Re + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # 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 + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * inv_Re + + 2.0 / 3.0 * v2_xy * inv_Re + - v1_yy * inv_Re + - v2_xy * inv_Re ) + # 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 + # stress tensor from y-direction + - v1_xy * inv_Re + - v2_xx * inv_Re + - 4.0 / 3.0 * v2_yy * inv_Re + + 2.0 / 3.0 * v1_xy * inv_Re ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * inv_Re + + 2.0 / 3.0 * v2_xy * v1 * inv_Re + - 4.0 / 3.0 * v1_x * v1_x * inv_Re + + 2.0 / 3.0 * v2_y * v1_x * inv_Re + - v1_xy * v2 * inv_Re + - v2_xx * v2 * inv_Re + - v1_y * v2_x * inv_Re + - v2_x * v2_x * inv_Re + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * inv_Re + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * inv_Re + - v2_xy * v1 * inv_Re + - v1_y * v1_y * inv_Re + - v2_x * v1_y * inv_Re + - 4.0 / 3.0 * v2_yy * v2 * inv_Re + + 2.0 / 3.0 * v1_xy * v2 * inv_Re + - 4.0 / 3.0 * v2_y * v2_y * inv_Re + + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * inv_Re ) + + return SVector(du1, du2, du3, du4) +end + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionViscousWall(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) + +# 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) + +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 from 0.0 to 1.5 +tspan = (0.0, .50) +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) + +############################################################################### +# 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_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..4784cb7f890 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl @@ -0,0 +1,77 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=1000, Prandtl=0.72, + Mach_freestream=0.1, kappa=1.0) + +# 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) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = 0.1 + rho = 1.0 + u, v = 0, 0 + p = 1.0 / (Ma^2 * equations.gamma) + return prim2cons(SVector(rho, u, v, p), equations) +end +initial_condition = initial_condition_cavity + +# BC types +velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0)) +velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_lid = BoundaryConditionViscousWall(velocity_bc_lid, heat_bc) +boundary_condition_cavity = BoundaryConditionViscousWall(velocity_bc_cavity, heat_bc) + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_slip_wall + +boundary_conditions_parabolic = (; x_neg = boundary_condition_cavity, + y_neg = boundary_condition_cavity, + y_pos = boundary_condition_lid, + x_pos = boundary_condition_cavity) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.5 +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) +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 5ae580b78255384ec077ee7b516ac04ad7b8db70 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 15 Jul 2022 15:38:01 -0500 Subject: [PATCH 62/71] adding AD checks --- .../elixir_navier_stokes_convergence.jl | 48 ++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl index 8909578982a..7ca891b85be 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -211,6 +211,52 @@ 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, +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 + +############################################################################### +# check accuracy of parabolic terms + +u_exact = (x, y, t) -> initial_condition_navier_stokes_convergence_test(SVector(x,y), t, equations) +dqdx = (x, y, t) -> ForwardDiff.derivative(x -> cons2prim(u_exact(x, y, t), equations_parabolic), x) +dqdy = (x, y, t) -> ForwardDiff.derivative(y -> cons2prim(u_exact(x, y, t), equations_parabolic), y) + +g = (x, y, t) -> Trixi.flux(cons2prim(u_exact(x,y,t), equations_parabolic), + (dqdx(x,y,t), dqdy(x,y,t)), equations_parabolic) +gx = (x, y, t) -> g(x, y, t)[1] +gy = (x, y, t) -> g(x, y, t)[2] + +dgdx = (x,y,t) -> ForwardDiff.derivative(x -> gx(x,y,t), x) +dgdy = (x,y,t) -> ForwardDiff.derivative(y -> gy(x,y,t), y) +div_f = (x,y,t) -> dgdx(x,y,t) + dgdy(x,y,t) + +x = vec(semi.cache.elements.node_coordinates[1, :, :, :]) +y = vec(semi.cache.elements.node_coordinates[2, :, :, :]) +t = tspan[end] + +trixi_wrap(u) = Trixi.wrap_array(vec(reinterpret(reshape, Float64, u)), semi) +dg = semi.solver +dg_parabolic = semi.solver_parabolic +@unpack cache, cache_parabolic = semi +boundary_conditions = semi.boundary_conditions_parabolic +u = u_exact.(x, y, t) +u = trixi_wrap(u) + +# check gradient +@unpack u_transformed, u_grad = cache_parabolic +Trixi.transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, dg_parabolic, cache, cache_parabolic) +Trixi.calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, cache_parabolic) + +println("error in parabolic gradients: $(maximum(abs.(u_grad[1] - trixi_wrap(dqdx.(x,y,t)))) + maximum(abs.(u_grad[2] - trixi_wrap(dqdy.(x,y,t)))))") + +# check fluxes +Trixi.calc_viscous_fluxes!(u_grad, u_transformed, mesh, equations_parabolic, dg, cache, cache_parabolic) +println("error in parabolic fluxes: $(maximum(abs.(u_grad[1] - trixi_wrap(gx.(x,y,t)))) + maximum(abs.(u_grad[2] - trixi_wrap(gy.(x,y,t)))))") + +# check divergence +du = similar(u) +Trixi.reset_du!(du, dg, cache) +Trixi.calc_divergence!(du, u_transformed, t, u_grad, mesh, equations_parabolic, + boundary_conditions, dg, dg_parabolic, cache, cache_parabolic) +println("error in parabolic divergence: $(maximum(abs.(du - trixi_wrap(div_f.(x,y,t)))))") \ No newline at end of file From bfffc0035c31013918dd53d85b1d234e90a30eed Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 14:26:17 -0500 Subject: [PATCH 63/71] standardizing parameters in convergence elixirs --- .../elixir_navier_stokes_convergence.jl | 4 +- .../elixir_navier_stokes_convergence.jl | 49 +------------------ 2 files changed, 3 insertions(+), 50 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl index a72ab18ab1b..f7d96c17c5f 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -14,13 +14,13 @@ equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=ge Mach_freestream=0.5, kappa=1.0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), +dg = DGMulti(polydeg = polydeg, element_type = Quad(), 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); periodicity=(true, false), is_on_boundary) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); periodicity=(true, false), is_on_boundary) # Define initial condition # Note: If you change the parameters here, also change it in the corresponding source terms diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl index 7ca891b85be..d731af84e32 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -4,7 +4,7 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -get_Re() = 1000 +get_Re() = 100 get_Pr() = .72 equations = CompressibleEulerEquations2D(1.4) @@ -25,7 +25,6 @@ mesh = TreeMesh(coordinates_min, coordinates_max, periodicity=(true, false), n_cells_max=30_000) # set maximum capacity of tree data structure - # Define initial condition # Note: If you change the parameters here, also change it in the corresponding source terms function initial_condition_navier_stokes_convergence_test(x, t, equations) @@ -193,7 +192,6 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), source_terms=source_terms_navier_stokes_convergence_test) - ############################################################################### # ODE solvers, callbacks etc. @@ -215,48 +213,3 @@ sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, dt = save_everystep=false, callback=callbacks) summary_callback() # print the timer summary -############################################################################### -# check accuracy of parabolic terms - -u_exact = (x, y, t) -> initial_condition_navier_stokes_convergence_test(SVector(x,y), t, equations) -dqdx = (x, y, t) -> ForwardDiff.derivative(x -> cons2prim(u_exact(x, y, t), equations_parabolic), x) -dqdy = (x, y, t) -> ForwardDiff.derivative(y -> cons2prim(u_exact(x, y, t), equations_parabolic), y) - -g = (x, y, t) -> Trixi.flux(cons2prim(u_exact(x,y,t), equations_parabolic), - (dqdx(x,y,t), dqdy(x,y,t)), equations_parabolic) -gx = (x, y, t) -> g(x, y, t)[1] -gy = (x, y, t) -> g(x, y, t)[2] - -dgdx = (x,y,t) -> ForwardDiff.derivative(x -> gx(x,y,t), x) -dgdy = (x,y,t) -> ForwardDiff.derivative(y -> gy(x,y,t), y) -div_f = (x,y,t) -> dgdx(x,y,t) + dgdy(x,y,t) - -x = vec(semi.cache.elements.node_coordinates[1, :, :, :]) -y = vec(semi.cache.elements.node_coordinates[2, :, :, :]) -t = tspan[end] - -trixi_wrap(u) = Trixi.wrap_array(vec(reinterpret(reshape, Float64, u)), semi) -dg = semi.solver -dg_parabolic = semi.solver_parabolic -@unpack cache, cache_parabolic = semi -boundary_conditions = semi.boundary_conditions_parabolic -u = u_exact.(x, y, t) -u = trixi_wrap(u) - -# check gradient -@unpack u_transformed, u_grad = cache_parabolic -Trixi.transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, dg_parabolic, cache, cache_parabolic) -Trixi.calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, cache_parabolic) - -println("error in parabolic gradients: $(maximum(abs.(u_grad[1] - trixi_wrap(dqdx.(x,y,t)))) + maximum(abs.(u_grad[2] - trixi_wrap(dqdy.(x,y,t)))))") - -# check fluxes -Trixi.calc_viscous_fluxes!(u_grad, u_transformed, mesh, equations_parabolic, dg, cache, cache_parabolic) -println("error in parabolic fluxes: $(maximum(abs.(u_grad[1] - trixi_wrap(gx.(x,y,t)))) + maximum(abs.(u_grad[2] - trixi_wrap(gy.(x,y,t)))))") - -# check divergence -du = similar(u) -Trixi.reset_du!(du, dg, cache) -Trixi.calc_divergence!(du, u_transformed, t, u_grad, mesh, equations_parabolic, - boundary_conditions, dg, dg_parabolic, cache, cache_parabolic) -println("error in parabolic divergence: $(maximum(abs.(du - trixi_wrap(div_f.(x,y,t)))))") \ No newline at end of file From 17eafbb876388c7a9f699a06ac6a9aab2da13f08 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 14:26:43 -0500 Subject: [PATCH 64/71] minor cleanup --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index a126da4afe9..78445a04b34 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -98,22 +98,16 @@ function calc_divergence!(du, u, t, viscous_flux, right_element = interfaces.neighbor_ids[2, interface] if orientations[interface] == 1 - # TODO Make this cleaner (remove the let) - let u = viscous_flux[1] - # interface in x-direction - for j in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces.u[1, v, j, interface] = u[v, nnodes(dg), j, left_element] - interfaces.u[2, v, j, interface] = u[v, 1, j, right_element] - end + # interface in x-direction + for j in eachnode(dg), v in eachvariable(equations_parabolic) + interfaces.u[1, v, j, interface] = viscous_flux[1][v, nnodes(dg), j, left_element] + interfaces.u[2, v, j, interface] = viscous_flux[1][v, 1, j, right_element] end else # if orientations[interface] == 2 - # TODO Make this cleaner (remove the let) - let u = viscous_flux[2] - # interface in y-direction - for i in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces.u[1, v, i, interface] = u[v, i, nnodes(dg), left_element] - interfaces.u[2, v, i, interface] = u[v, i, 1, right_element] - end + # interface in y-direction + for i in eachnode(dg), v in eachvariable(equations_parabolic) + interfaces.u[1, v, i, interface] = viscous_flux[2][v, i, nnodes(dg), left_element] + interfaces.u[2, v, i, interface] = viscous_flux[2][v, i, 1, right_element] end end end From 0c4e098bd26fe42851c9ec6c3f791c3fb8f9615e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 14:29:36 -0500 Subject: [PATCH 65/71] revert DGMulti CNS elixir setup back to the one used in tests --- examples/dgmulti_2d/elixir_navier_stokes_convergence.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl index f7d96c17c5f..b3af1e0f14f 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -14,7 +14,7 @@ equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=ge Mach_freestream=0.5, kappa=1.0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -dg = DGMulti(polydeg = polydeg, element_type = Quad(), approximation_type = Polynomial(), +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), volume_integral = VolumeIntegralWeakForm()) From b70bbc046c116624ace2d0ab46cbdc9a44ae0690 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 14:29:50 -0500 Subject: [PATCH 66/71] adding TreeMesh CNS convergence test --- test/test_parabolic_2d.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 79084b35878..fee3e83049a 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -110,6 +110,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.0021116725306635146, 0.0034322351490824465, 0.003874252819611102, 0.012469246082522416], + linf = [0.012006418939297214, 0.03552087120958058, 0.02451274749176294, 0.11191122588577151] + ) + end + end # Clean up afterwards: delete Trixi output directory From 13e0c8fa0f968c7da1bab33ba837afc8643a0c55 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 17:18:33 -0500 Subject: [PATCH 67/71] removing redundant elixir --- .../elixir_navier_stokes_source_terms.jl | 272 ------------------ 1 file changed, 272 deletions(-) delete mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl deleted file mode 100644 index bd85946f1b3..00000000000 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_source_terms.jl +++ /dev/null @@ -1,272 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the ideal compressible Navier-Stokes equations - -equations = CompressibleEulerEquations2D(1.4) -# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition -# I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesEquations2D(equations, - Reynolds=1000, - Prandtl=0.72, - Mach_freestream=0.5, - kappa=1.0 # thermal diffusivity - ) - -# 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) - -coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) -coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) - -# Create a uniformly refined mesh with periodic boundaries -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=4, - n_cells_max=30_000) # set maximum capacity of tree data structure - -# Define initial condition -# Note: If you change the parameters here, also change it in the corresponding source terms -function initial_condition_navier_stokes_convergence_test(x, t, equation::CompressibleEulerEquations2D) - - # Amplitude and shift - A = 0.5 - c = 2.0 - - # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t - - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0)) ) * cos(pi_t) - v2 = v1 - p = rho^2 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end - - -initial_condition = initial_condition_navier_stokes_convergence_test - - -@inline function source_terms_navier_stokes_convergence_test(u, x, t, - equations::CompressibleNavierStokesEquations2D) - # Same settings as in `initial_condition` - # Amplitude and shift - A = 0.5 - c = 2.0 - - # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t - - # compute the manufactured solution and all necessary derivatives - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) - rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) - rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) - rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - - v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) - v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) - v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) - v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) - - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) - - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) - v2 = v1 - v2_t = v1_t - v2_x = v1_x - v2_y = v1_y - v2_xx = v1_xx - v2_xy = v1_xy - v2_yy = v1_yy - - p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_y = 2.0 * rho * rho_y - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x - p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y - - # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * equations.inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) - E_t = p_t * equations.inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t - E_x = p_x * equations.inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - E_y = p_y * equations.inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y - - # Some convenience constants - T_const = equations.gamma * equations.inv_gamma_minus_one * equations.kappa / equations.Pr - inv_Re = 1.0 / equations.Re - inv_rho_cubed = 1.0 / (rho^3) - - # compute the source terms - # density equation - du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y - - # 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 - # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * inv_Re - + 2.0 / 3.0 * v2_xy * inv_Re - - v1_yy * inv_Re - - v2_xy * inv_Re ) - # 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 - # stress tensor from y-direction - - v1_xy * inv_Re - - v2_xx * inv_Re - - 4.0 / 3.0 * v2_yy * inv_Re - + 2.0 / 3.0 * v1_xy * inv_Re ) - # total energy equation - du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) - + v2_y * (E + p) + v2 * (E_y + p_y) - # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * inv_Re - + 2.0 / 3.0 * v2_xy * v1 * inv_Re - - 4.0 / 3.0 * v1_x * v1_x * inv_Re - + 2.0 / 3.0 * v2_y * v1_x * inv_Re - - v1_xy * v2 * inv_Re - - v2_xx * v2 * inv_Re - - v1_y * v2_x * inv_Re - - v2_x * v2_x * inv_Re - - T_const * inv_rho_cubed * ( p_xx * rho * rho - - 2.0 * p_x * rho * rho_x - + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * inv_Re - # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * inv_Re - - v2_xy * v1 * inv_Re - - v1_y * v1_y * inv_Re - - v2_x * v1_y * inv_Re - - 4.0 / 3.0 * v2_yy * v2 * inv_Re - + 2.0 / 3.0 * v1_xy * v2 * inv_Re - - 4.0 / 3.0 * v2_y * v2_y * inv_Re - + 2.0 / 3.0 * v1_x * v2_y * inv_Re - - T_const * inv_rho_cubed * ( p_yy * rho * rho - - 2.0 * p_y * rho * rho_y - + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * inv_Re ) - - return SVector(du1, du2, du3, du4) -end - - -# TODO: Wasn't sure of the call structure, but this should be what we need -function boundary_condition_no_slip_adiabatic_wall(u_inner, orientation, direction, - x, t, surface_flux_function, - equations::CompressibleEulerEquations2D) - u_boundary = SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) - - # Calculate boundary flux - if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary - flux = surface_flux_function(u_inner, u_boundary, orientation, equations) - else # u_boundary is "left" of boundary, u_inner is "right" of boundary - flux = surface_flux_function(u_boundary, u_inner, orientation, equations) - end - - return flux -end - - -# TODO: Wasn't sure of the call structure, but this should be what we need -function boundary_condition_no_slip_adiabatic_wall_neumann(grad_u_inner, orientation, direction, - x, t, surface_flux_function, - equations::CompressibleNavierStokesEquations2D) - # Copy the inner gradients to an external state array - grad_u_boundary .= grad_u_inner - - # Project out the appropriate temperature gradient pieces. Wasn't sure of array shape - if orientation == 1 - grad_u_norm = grad_u[1,3] # temperature gradient in x-direction - u_x_tangent = grad_u[1,3] - grad_u_norm - u_y_tangent = grad_u[2,3] - - # Update the external state gradients - grad_u_boundary[1,3] = u_x_tangent - grad_u_norm - grad_u_boundary[2,3] = u_y_tangent - else # orientation == 2 - grad_u_norm = grad_u[2,3] # temperature gradient in y-direction - u_x_tangent = grad_u[1,3] - u_y_tangent = grad_u[2,3] - grad_u_norm - - # Update the external state gradients - grad_u_boundary[1,3] = u_x_tangent - grad_u_boundary[2,3] = u_y_tangent - grad_u_norm - end - - # Calculate boundary flux (just averages so has no orientation I think) - flux = surface_flux_function(grad_u_inner, grad_u_boundary, orientation, equations) - - return flux -end - - -# Below is my best guess as to how to set periodic in x direction and walls -# in the y direcitons -boundary_condition= ( - x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_no_slip_adiabatic_wall, - y_pos=boundary_condition_no_slip_adiabatic_wall, - ) -boundary_conditions_parabolic = ( - x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_no_slip_adiabatic_wall_neumann, - y_pos=boundary_condition_no_slip_adiabatic_wall_neumann, - ) - -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolicParabolic(mesh, - (equations, equations_parabolic), - initial_condition, solver; - boundary_conditions=(boundary_conditions, - boundary_conditions_parabolic)) - - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0.0 to 1.5 -tspan = (0.0, 1.5) -ode = semidiscretize(semi, tspan); - -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) - -# The AliveCallback prints short status information in regular intervals -alive_callback = AliveCallback(analysis_interval=analysis_interval) - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) - - -############################################################################### -# run the simulation - -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -time_int_tol = 1.0e-8 -sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, - save_everystep=false, callback=callbacks) - -# Print the timer summary -summary_callback() From 8cecb5f47984b313bad84b8c8899934f9b9bd7c2 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 17:26:05 -0500 Subject: [PATCH 68/71] add more tests --- .../elixir_navier_stokes_lid_driven_cavity.jl | 4 ++-- test/test_parabolic_2d.jl | 16 ++++++++++++++++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl index 4784cb7f890..c42e7b605ab 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl @@ -57,11 +57,11 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers, callbacks etc. # Create ODE problem with time span from 0.0 to 1.5 -tspan = (0.0, 10.0) +tspan = (0.0, 25.0) ode = semidiscretize(semi, tspan); summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval=10) +alive_callback = AliveCallback(alive_interval=100) analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) callbacks = CallbackSet(summary_callback, alive_callback) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index fee3e83049a..41b8f76d6ad 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -94,6 +94,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "DGMulti: elixir_navier_stokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_lid_driven_cavity.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.5), + l2 = [0.0002215612522711349, 0.028318325921400257, 0.009509168701069035, 0.028267900513539248], + linf = [0.0015622789413053395, 0.14886653390741342, 0.07163235655334241, 0.19472785105216417] + ) + end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion.jl"), initial_refinement_level = 2, tspan=(0.0, 0.4), polydeg=5, @@ -118,6 +126,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "TreeMesh2D: elixir_navier_stokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_lid_driven_cavity.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.0001514457152968994, 0.018766076072331786, 0.007065070765651992, 0.020839900573430787], + linf = [0.0014523369373645734, 0.12366779944955876, 0.055324509971157544, 0.1609992780534526] + ) + end + end # Clean up afterwards: delete Trixi output directory From adb8eca2c14eb7b1f2dbfdfc140384c1d644990c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 28 Jul 2022 17:26:56 -0500 Subject: [PATCH 69/71] add more test --- test/test_parabolic_2d.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 41b8f76d6ad..42d6b1d3752 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -70,6 +70,14 @@ isdir(outdir) && rm(outdir, recursive=true) @test getindex.(du, 1) ≈ 2 * y end + @trixi_testset "DGMulti: elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.2485803335154642], + linf = [1.079606969242132] + ) + end + @trixi_testset "DGMulti: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_periodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), From 3281b49a070d232ba9d304143ba9e263fa72ef3a Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 29 Jul 2022 09:32:53 +0200 Subject: [PATCH 70/71] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../elixir_navier_stokes_convergence.jl | 6 +++--- .../elixir_navier_stokes_lid_driven_cavity.jl | 2 +- .../elixir_advection_diffusion_nonperiodic.jl | 2 +- .../elixir_navier_stokes_convergence.jl | 6 +++--- .../elixir_navier_stokes_lid_driven_cavity.jl | 2 +- src/equations/compressible_navier_stokes_2d.jl | 16 +++++++++------- 6 files changed, 18 insertions(+), 16 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl index b3af1e0f14f..e922a9cad12 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -5,7 +5,7 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations get_Re() = 100 -get_Pr() = .72 +get_Pr() = 0.72 equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition @@ -187,8 +187,8 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.5 -tspan = (0.0, .50) +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl index d93cff55c5f..bdcc3a60607 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl @@ -51,7 +51,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.5 +# Create ODE problem with time span `tspan` tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan); diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index 74d0e83bdf0..21e4866189b 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -56,7 +56,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.5 +# Create ODE problem with time span `tspan` tspan = (0.0, 1.5) ode = semidiscretize(semi, tspan); diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl index d731af84e32..ee2307a90e7 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -5,7 +5,7 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations get_Re() = 100 -get_Pr() = .72 +get_Pr() = 0.72 equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition @@ -195,8 +195,8 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.5 -tspan = (0.0, .50) +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl index c42e7b605ab..2a8fb635895 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl @@ -56,7 +56,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.5 +# Create ODE problem with time span `tspan` tspan = (0.0, 25.0) ode = semidiscretize(semi, tspan); diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 1ae4ba87ced..e9be251077a 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -8,23 +8,25 @@ `CompressibleNavierStokesEquations2D` represents the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective from -the `CompressibleEulerEquations2D`. +the [`CompressibleEulerEquations2D`](@ref). -gamma: adiabatic constant, -Re: Reynolds number, -Pr: Prandtl number, -Ma_inf: free-stream Mach number -kappa: thermal diffusivity for Fick's law +- `gamma`: adiabatic constant, +- `Re`: Reynolds number, +- `Pr`: Prandtl number, +- `Ma_inf`: free-stream Mach number +- `kappa`: thermal diffusivity for Fick's law For the particular scaling the vicosity is set internally to be μ = 1/Re. Further, the nondimensionalization takes the density-temperature-sound speed as the principle quantities such that +``` rho_inf = 1.0 T_ref = 1.0 c_inf = 1.0 p_inf = 1.0 / gamma u_inf = Ma_inf R = 1.0 / gamma +``` Other normalization strategies exist, see the reference below for details. - Marc Montagnac (2013) @@ -35,7 +37,7 @@ The scaling used herein is Section 4.5 of the reference. In two spatial dimensions we require gradients for three quantities, e.g., primitive quantities - grad(u), grad(v), and grad(T) + grad(v_1), grad(v_2), and grad(T) or the entropy variables grad(w_2), grad(w_3), grad(w_4) where From 98ea1bdc83e269062f933ba16dcc40e2f8fc78aa Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 29 Jul 2022 09:12:02 -0500 Subject: [PATCH 71/71] add docstrings --- .../compressible_navier_stokes_2d.jl | 34 ++++++++++++++++--- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index e9be251077a..3003de4fda2 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -209,23 +209,49 @@ end return T end -# abstract container for wall-type boundary conditions +""" + struct BoundaryConditionViscousWall{V, H} + +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. +""" struct BoundaryConditionViscousWall{V, H} boundary_condition_velocity::V boundary_condition_heat_flux::H end -# no slip velocity BC +""" + struct NoSlip{F} + +Creates a no-slip boundary condition with field `boundary_value_function`, which +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 -# isothermal temperature BC +""" + struct Isothermal{F} + +Creates an isothermal temperature boundary condition with field `boundary_value_function`, +which 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 -# adiabatic temperature BC +""" + struct Adiabatic{F} + +Creates an adiabatic temperature boundary condition with field `boundary_value_function`, +which should be a function with signature `boundary_value_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