From 3f48f64d07ab90cc42c324a430e3b372162d91e4 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 2 Feb 2023 16:13:40 +0100 Subject: [PATCH] Compressible Navier-Stokes on TreeMesh3D (#1239) * initial commit of 3d approximation for CNS * typo fix in CNS 2D file * add TGV elixir for compressible Navier-Stokes * add DGMulti CNS 3D elixir * need TreeMesh versions of slip wall BC in 3D compressible Euler for convergence testing * update comments in CNS 3D file * add convergence test elixir for 3D compressible Navier-Stokes * remove debug statements from new elixir * adjust maximum tree size in new elixir * add convergence test for DGMulti 3D * add battery of 3D parabolic tests. All pass locally * remove old comments and TODOs * bug fix in divergence B analysis routine * add vorticity and enstrophy functions * add first version of an enstrophy callback * add enstrophy callback to TGV elixir * update TGV test values since switched surface flux to flux_hllc * adding draft DGMulti analysis callback * adding draft DGMulti analysis callback forgot the return statement * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * make sure Experimental code text renders properly in docs * adding extra analysis integrals * remove unused file path from parabolic test scripts * remove threading for reduction * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper * add designation to the 2D test set * remove specific types in new integrate routine * update docstrings for CNS to match the constructor * Update examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl --------- Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> Co-authored-by: Michael Schlottke-Lakemper --- .../elixir_navierstokes_convergence.jl | 253 +++++++ ...elixir_navierstokes_taylor_green_vortex.jl | 71 ++ .../elixir_navierstokes_convergence.jl | 267 +++++++ ...elixir_navierstokes_taylor_green_vortex.jl | 76 ++ src/Trixi.jl | 5 +- src/callbacks_step/analysis.jl | 18 + src/callbacks_step/analysis_dg3d.jl | 20 +- src/callbacks_step/analysis_dgmulti.jl | 37 + src/equations/compressible_euler_3d.jl | 46 ++ .../compressible_navier_stokes_2d.jl | 12 +- .../compressible_navier_stokes_3d.jl | 380 ++++++++++ src/equations/equations_parabolic.jl | 1 + src/solvers/dgsem_tree/dg.jl | 1 + src/solvers/dgsem_tree/dg_3d_parabolic.jl | 667 ++++++++++++++++++ test/runtests.jl | 1 + test/test_parabolic_2d.jl | 3 +- test/test_parabolic_3d.jl | 86 +++ 17 files changed, 1932 insertions(+), 12 deletions(-) create mode 100644 examples/dgmulti_3d/elixir_navierstokes_convergence.jl create mode 100644 examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl create mode 100644 examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl create mode 100644 examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl create mode 100644 src/equations/compressible_navier_stokes_3d.jl create mode 100644 src/solvers/dgsem_tree/dg_3d_parabolic.jl create mode 100644 test/test_parabolic_3d.jl diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..b08eda5cf8c --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl @@ -0,0 +1,253 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) +mesh = DGMultiMesh(dg, cells_per_dimension=(8, 8, 8); periodicity=(true, false, true), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl new file mode 100644 index 00000000000..9d0d926338a --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -0,0 +1,71 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), + Prandtl=prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + +The classical inviscid Taylor-Green vortex. +""" +function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3]) + v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3]) + v3 = 0.0 + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3])) + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = GaussSBP(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = ( 1.0, 1.0, 1.0) .* pi +mesh = DGMultiMesh(dg; coordinates_min, coordinates_max, + cells_per_dimension=(8, 8, 8), + periodicity=(true, true, true)) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, dg) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg), + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)) +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_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..9199174652a --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,267 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = ( 1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + periodicity=(true, false, true), + n_cells_max=50_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_top_bottom, + y_pos = boundary_condition_top_bottom, + z_neg = boundary_condition_periodic, + z_pos = boundary_condition_periodic) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary + diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl new file mode 100644 index 00000000000..d3f1bcbd3bd --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -0,0 +1,76 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), + Prandtl=prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + +The classical inviscid Taylor-Green vortex. +""" +function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3]) + v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3]) + v3 = 0.0 + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3])) + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = ( 1.0, 1.0, 1.0) .* pi +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=100_000) + + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval,) + +callbacks = CallbackSet(summary_callback, + analysis_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/src/Trixi.jl b/src/Trixi.jl index 81bdee58d58..e3c4638ea5c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -136,7 +136,7 @@ export AcousticPerturbationEquations2D, LinearizedEulerEquations2D export LaplaceDiffusion2D, - CompressibleNavierStokesDiffusion2D + CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D export GradientVariablesPrimitive, GradientVariablesEntropy @@ -181,7 +181,8 @@ export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_c export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure -export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity +export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, + enstrophy export lake_at_rest_error export ncomponents, eachcomponent diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index b6880509d56..46a6b6c3fbe 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -549,6 +549,21 @@ pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) +# Special analyze for `SemidiscretizationHyperbolicParabolic` such that +# precomputed gradients are available. For now only implemented for the `enstrophy` +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. +function analyze(quantity::typeof(enstrophy), du, u, t, semi::SemidiscretizationHyperbolicParabolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + cache_parabolic = semi.cache_parabolic + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, cache_parabolic) +end +function analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, cache_parabolic) + integrate(quantity, u, mesh, equations, equations_parabolic, solver, cache, cache_parabolic, normalize=true) +end + + function entropy_timederivative end pretty_form_utf(::typeof(entropy_timederivative)) = "∑∂S/∂U ⋅ Uₜ" pretty_form_ascii(::typeof(entropy_timederivative)) = "dsdu_ut" @@ -573,6 +588,9 @@ pretty_form_ascii(::typeof(energy_magnetic)) = "e_magnetic" pretty_form_utf(::typeof(cross_helicity)) = "∑v⋅B" pretty_form_ascii(::typeof(cross_helicity)) = "v_dot_B" +pretty_form_utf(::typeof(enstrophy)) = "∑enstrophy" +pretty_form_ascii(::typeof(enstrophy)) = "enstrophy" + pretty_form_utf(::Val{:l2_divb}) = "L2 ∇⋅B" pretty_form_ascii(::Val{:l2_divb}) = "l2_divb" diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 69bfb291472..77cf1f819ea 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -199,6 +199,22 @@ function integrate(func::Func, u, end +function integrate(func::Func, u, + mesh::TreeMesh{3}, + equations, equations_parabolic, + dg::DGSEM, + cache, cache_parabolic; normalize=true) where {Func} + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, k, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, k, element) + gradients_1_local = get_node_vars(gradients_x, equations_parabolic, dg, i, j, k, element) + gradients_2_local = get_node_vars(gradients_y, equations_parabolic, dg, i, j, k, element) + gradients_3_local = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k, element) + return func(u_local, (gradients_1_local, gradients_2_local, gradients_3_local), equations_parabolic) + end +end + + function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, equations, dg::DG, cache) @@ -220,7 +236,7 @@ function analyze(::Val{:l2_divb}, du, u, t, for l in eachnode(dg) divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] + derivative_matrix[j, l] * u[7, i, l, k, element] + - derivative_matrix[k, l] * u[7, i, j, l, element] ) + derivative_matrix[k, l] * u[8, i, j, l, element] ) end divb *= cache.elements.inverse_jacobian[element] divb^2 @@ -262,7 +278,7 @@ function analyze(::Val{:linf_divb}, du, u, t, for l in eachnode(dg) divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] + derivative_matrix[j, l] * u[7, i, l, k, element] + - derivative_matrix[k, l] * u[7, i, j, l, element] ) + derivative_matrix[k, l] * u[8, i, j, l, element] ) end divb *= cache.elements.inverse_jacobian[element] linf_divb = max(linf_divb, abs(divb)) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index c67dc892338..9489c8bc753 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -135,6 +135,43 @@ function analyze(::Val{:linf_divb}, du, u, t, return linf_divB end +function integrate(func::typeof(enstrophy), u, + mesh::DGMultiMesh, + equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, + dg::DGMulti, + cache, cache_parabolic; normalize=true) + + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + + # allocate local storage for gradients. + # TODO: can we avoid allocating here? + local_gradient_quadrature_values = ntuple(_ -> similar(cache_parabolic.local_u_values_threaded), 3) + + integral = zero(eltype(u)) + for e in eachelement(mesh, dg) + u_quadrature_values = cache_parabolic.local_u_values_threaded[Threads.threadid()] + gradient_x_quadrature_values = local_gradient_quadrature_values[1][Threads.threadid()] + gradient_y_quadrature_values = local_gradient_quadrature_values[2][Threads.threadid()] + gradient_z_quadrature_values = local_gradient_quadrature_values[3][Threads.threadid()] + + # interpolate to quadrature on each element + apply_to_each_field(mul_by(dg.basis.Vq), u_quadrature_values, view(u, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e)) + apply_to_each_field(mul_by(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e)) + + # integrate over the element + for i in eachindex(u_quadrature_values) + gradients_i = SVector(gradient_x_quadrature_values[i], + gradient_y_quadrature_values[i], + gradient_z_quadrature_values[i]) + integral += mesh.md.wJq[i, e] * func(u_quadrature_values[i], gradients_i, equations) + end + end + return integral +end + + function create_cache_analysis(analyzer, mesh::DGMultiMesh, equations, dg::DGMulti, cache, RealT, uEltype) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index a795a235164..e66fac5ab4c 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -323,6 +323,52 @@ Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of zero(eltype(u_inner))) * norm_ end +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations3D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(1.0, 0.0, 0.0) + elseif orientation == 2 + normal_direction = SVector(0.0, 1.0, 0.0) + else # orientation == 3 + normal_direction = SVector(0.0, 0.0, 1.0) + end + + # compute and return the flux using `boundary_condition_slip_wall` routine above + return boundary_condition_slip_wall(u_inner, normal_direction, direction, + x, t, surface_flux_function, equations) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations3D) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations3D) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x, -y, and -z sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, equations) + end + + return boundary_flux +end # Calculate 1D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations3D) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 9f7758aa476..c509015ba7b 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,15 +1,14 @@ @doc raw""" - CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, equations, + CompressibleNavierStokesDiffusion2D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) -These equations contain the diffusion (i.e. parabolic) terms applied +Contains the diffusion (i.e. parabolic) terms applied to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations2D`](@ref). -- `gamma`: adiabatic constant, +- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `mu`: dynamic viscosity, - `Pr`: Prandtl number, -- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. @@ -37,7 +36,8 @@ where the system is closed with the ideal gas assumption giving ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` -as the pressure. The terms on the right hand side of the system above +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations2D`](@ref). +The terms on the right hand side of the system above are built from the viscous stress tensor ```math \underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} @@ -142,7 +142,7 @@ gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientV # Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 -# of the paper by Rueda-Ramíreza, Hennemann, Hindenlang, Winters, and Gassner +# of the paper by Rueda-Ramírez, Hennemann, Hindenlang, Winters, and Gassner # "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive # MHD Equations. Part II: Subcell Finite Volume Shock Capturing" # where one sets the magnetic field components equal to 0. diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl new file mode 100644 index 00000000000..d4821c55e58 --- /dev/null +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -0,0 +1,380 @@ +@doc raw""" + CompressibleNavierStokesDiffusion3D(equations; mu, Pr, + gradient_variables=GradientVariablesPrimitive()) + +Contains the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective terms from +the [`CompressibleEulerEquations3D`](@ref). + +- `equations`: instance of the [`CompressibleEulerEquations3D`](@ref) +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `gradient_variables`: which variables the gradients are taken with respect to. + Defaults to `GradientVariablesPrimitive()`. + +Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., +[``\mu``] = kg m⁻¹ s⁻¹. + +The particular form of the compressible Navier-Stokes implemented is +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho \mathbf{v} \\ \rho e +\end{pmatrix} ++ +\nabla \cdot +\begin{pmatrix} + \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e + p) \mathbf{v} +\end{pmatrix} += +\nabla \cdot +\begin{pmatrix} +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +\end{pmatrix} +``` +where the system is closed with the ideal gas assumption giving +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right) +``` +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations2D`](@ref). +The terms on the right hand side of the system above +are built from the viscous stress tensor +```math +\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} +``` +where ``\underline{I}`` is the ``3\times 3`` identity matrix and the heat flux is +```math +\nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} +``` +where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fick's law. +Under the assumption that the gas has a constant Prandtl number, +the thermal conductivity is +```math +\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}. +``` +From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see +that the gas constant `R` cancels and the heat flux becomes +```math +\nabla q = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) +``` +which is the form implemented below in the [`flux`](@ref) function. + +In two spatial dimensions we require gradients for three quantities, e.g., +primitive quantities +```math +\nabla v_1,\, \nabla v_2,\, \nabla v_3,\, \nabla T +``` +or the entropy variables +```math +\nabla w_2,\, \nabla w_3,\, \nabla w_4\, \nabla w_5 +``` +where +```math +w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} +``` + +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. +""" +struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3}} <: AbstractCompressibleNavierStokesDiffusion{3, 5} + # TODO: parabolic + # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations + # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + mu::RealT # viscosity + Pr::RealT # Prandtl number + kappa::RealT # thermal diffusivity for Fick's law + + equations_hyperbolic::E # CompressibleEulerEquations3D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +# default to primitive gradient variables +function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquations3D; + mu, Prandtl, + gradient_variables = GradientVariablesPrimitive()) + gamma = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + μ, Pr = promote(mu, Prandtl) + + # Under the assumption of constant Prandtl number the thermal conductivity + # constant is kappa = gamma μ / ((gamma-1) Pr). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * inv_gamma_minus_one / Pr + + CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, + μ, Pr, kappa, + equations, gradient_variables) +end + +# TODO: parabolic +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion3D) = ("v1", "v2", "v3", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion3D) = ("w2", "w3", "w4", "w5") + +varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusion3D) = + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) + +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) = cons2prim +gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) = cons2entropy + + +# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramírez, Hennemann, Hindenlang, Winters, and Gassner +# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive +# MHD Equations. Part II: Subcell Finite Volume Shock Capturing" +# where one sets the magnetic field components equal to 0. +function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion3D) + # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. + rho, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T) + # either computed directly or reverse engineered from the gradient of the entropy vairables + # by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, dv3dx, dTdx = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, dv3dy, dTdy = convert_derivative_to_primitive(u, gradients[2], equations) + _, dv1dz, dv2dz, dv3dz, dTdz = convert_derivative_to_primitive(u, gradients[3], equations) + + # Components of viscous stress tensor + + # Diagonal parts + # (4/3 * (v1)_x - 2/3 * ((v2)_y + (v3)_z) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * (dv2dy + dv3dz) + # (4/3 * (v2)_y - 2/3 * ((v1)_x + (v3)_z) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * (dv1dx + dv3dz) + # (4/3 * (v3)_z - 2/3 * ((v1)_x + (v2)_y) + tau_33 = 4.0 / 3.0 * dv3dz - 2.0 / 3.0 * (dv1dx + dv2dy) + + # Off diagonal parts, exploit that stress tensor is symmetric + # ((v1)_y + (v2)_x) + tau_12 = dv1dy + dv2dx # = tau_21 + # ((v1)_z + (v3)_x) + tau_13 = dv1dz + dv3dx # = tau_31 + # ((v2)_z + (v3)_y) + tau_23 = dv2dz + dv3dy # = tau_32 + + # Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho)) + # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) + # Note, the gas constant cancels under this formulation, so it is not present + # in the implementation + q1 = equations.kappa * dTdx + q2 = equations.kappa * dTdy + q3 = equations.kappa * dTdz + + # Constant dynamic viscosity is copied to a variable for readibility. + # Offers flexibility for dynamic viscosity via Sutherland's law where it depends + # on temperature and reference values, Ts and Tref such that mu(T) + mu = equations.mu + + if orientation == 1 + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = tau_13 * mu + f5 = ( v1 * tau_11 + v2 * tau_12 + v3 * tau_13 + q1 ) * mu + + return SVector(f1, f2, f3, f4, f5) + elseif orientation == 2 + # viscous flux components in the y-direction + # Note, symmetry is exploited for tau_12 = tau_21 + g1 = zero(rho) + g2 = tau_12 * mu # tau_21 * mu + g3 = tau_22 * mu + g4 = tau_23 * mu + g5 = ( v1 * tau_12 + v2 * tau_22 + v3 * tau_23 + q2 ) * mu + + return SVector(g1, g2, g3, g4, g5) + else # if orientation == 3 + # viscous flux components in the z-direction + # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32 + h1 = zero(rho) + h2 = tau_13 * mu # tau_31 * mu + h3 = tau_23 * mu # tau_32 * mu + h4 = tau_33 * mu + h5 = ( v1 * tau_13 + v2 * tau_23 + v3 * tau_33 + q3 ) * mu + + return SVector(h1, h2, h3, h4, h5) + end +end + + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + T = temperature(u, equations) + + return SVector(rho, v1, v2, v3, T) +end + +# Convert conservative variables to entropy +# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms +# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion2D`, +# but this may be confusing to new users. +cons2entropy(u, equations::CompressibleNavierStokesDiffusion3D) = cons2entropy(u, equations.equations_hyperbolic) +entropy2cons(w, equations::CompressibleNavierStokesDiffusion3D) = entropy2cons(w, equations.equations_hyperbolic) + +# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. +# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# variables into primitive variables. +@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + return u_transformed +end + +# TODO: parabolic. Make this more efficient! +@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + # note: this uses CompressibleNavierStokesDiffusion3D versions of cons2prim and entropy2cons + return cons2prim(entropy2cons(u_transformed, equations), equations) +end + + +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4, w_5) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, v3, T). +# Helpful because then the diffusive fluxes have the same form as on paper. +# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. +# TODO: parabolic; entropy stable viscous terms +@inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + return gradient +end + +# the first argument is always the "transformed" variables. +@inline function convert_derivative_to_primitive(w, gradient_entropy_vars, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + + # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back. + # We can fix this if we directly compute v1, v2, v3, T from the entropy variables + u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion3D + rho, rho_v1, rho_v2, rho_v3, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + T = temperature(u, equations) + + return SVector(gradient_entropy_vars[1], + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[5]), # grad(u) = T*(grad(w_2)+v1*grad(w_5)) + T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_3)+v2*grad(w_5)) + T * (gradient_entropy_vars[4] + v3 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_4)+v3*grad(w_5)) + T * T * gradient_entropy_vars[5] # grad(T) = T^2*grad(w_5)) + ) +end + + +# This routine is required because `prim2cons` is called in `initial_condition`, which +# is called with `equations::CompressibleEulerEquations3D`. This means it is inconsistent +# with `cons2prim(..., ::CompressibleNavierStokesDiffusion3D)` as defined above. +# TODO: parabolic. Is there a way to clean this up? +@inline prim2cons(u, equations::CompressibleNavierStokesDiffusion3D) = + prim2cons(u, equations.equations_hyperbolic) + + +@inline function temperature(u, equations::CompressibleNavierStokesDiffusion3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) + T = p / rho + return T +end + + +@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D) + # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v + + omega = vorticity(u, gradients, equations) + return 0.5 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2) +end + + +@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D) + # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, dv3dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, dv3dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) + _, dv1dz, dv2dz, dv3dz, _ = convert_derivative_to_primitive(u, gradients[3], equations) + + return SVector(dv3dy - dv2dz , dv1dz - dv3dx , dv2dx - dv1dy) +end + + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + return SVector(u_inner[1], v1, v2, v3, u_inner[5]) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + # rho, v1, v2, v3, _ = u_inner + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], normal_energy_flux) +end + + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) + return SVector(u_inner[1], v1, v2, v3, T) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) + return flux_inner +end + +# specialized BC impositions for GradientVariablesEntropy. + +# This should return a SVector containing the boundary values of entropy variables. +# Here, `w_inner` are the transformed variables (e.g., entropy variables). +# +# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions +# for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022. +# DOI: 10.1016/j.jcp.2021.110723 +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + negative_rho_inv_p = w_inner[5] # w_5 = -rho / p + return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p, -v3 * negative_rho_inv_p, negative_rho_inv_p) +end + +# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness. +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) + + # the entropy variables w2 = rho * v1 / p = v1 / T = -v1 * w5. Similarly for w3 and w4 + w5 = -1 / T + return SVector(w_inner[1], -v1 * w5, -v2 * w5, -v3 * w5, w5) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy}) + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], flux_inner[5]) +end diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 340cf0e4294..4fcdfff5f94 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -9,3 +9,4 @@ include("laplace_diffusion_2d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("compressible_navier_stokes_2d.jl") +include("compressible_navier_stokes_3d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 0ffaec2eced..e0066c2d2d8 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -64,6 +64,7 @@ include("dg_2d_parabolic.jl") # 3D DG implementation include("dg_3d.jl") +include("dg_3d_parabolic.jl") # Auxiliary functions that are specialized on this solver # as well as specialized implementations used to improve performance diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl new file mode 100644 index 00000000000..d3a47cb06be --- /dev/null +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -0,0 +1,667 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + +# This file collects all methods that have been updated to work with parabolic systems of equations +# +# assumptions: parabolic terms are of the form div(f(u, grad(u))) and +# will be discretized first order form as follows: +# 1. compute grad(u) +# 2. compute f(u, grad(u)) +# 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call) +# boundary conditions will be applied to both grad(u) and div(f(u, grad(u))). +function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + initial_condition, boundary_conditions_parabolic, source_terms, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @unpack u_transformed, gradients, flux_viscous = cache_parabolic + + # Convert conservative variables to a form more suitable for viscous flux calculations + @trixi_timeit timer() "transform variables" transform_variables!( + u_transformed, u, mesh, equations_parabolic, dg, parabolic_scheme, cache, cache_parabolic) + + # Compute the gradients of the transformed variables + @trixi_timeit timer() "calculate gradient" calc_gradient!( + gradients, u_transformed, t, mesh, equations_parabolic, boundary_conditions_parabolic, dg, + cache, cache_parabolic) + + # Compute and store the viscous fluxes + @trixi_timeit timer() "calculate viscous fluxes" calc_viscous_fluxes!( + flux_viscous, gradients, u_transformed, mesh, equations_parabolic, dg, cache, cache_parabolic) + + # The remainder of this function is essentially a regular rhs! for parabolic equations (i.e., it + # computes the divergence of the viscous fluxes) + # + # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have + # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # `interfaces` and `boundaries` containers in `cache_parabolic` to interpolate and store the + # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we + # do not need to recreate the existing data structure only with a different name, and c) we do not + # need to interpolate solutions *and* gradients to the surfaces. + + # TODO: parabolic; reconsider current data structure reuse strategy + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" calc_volume_integral!( + du, flux_viscous, mesh, equations_parabolic, dg, cache) + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" calc_interface_flux!( + cache_parabolic.elements.surface_flux_values, mesh, equations_parabolic, dg, cache_parabolic) + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" calc_boundary_flux_divergence!( + cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic, + dg.surface_integral, dg) + + # TODO: parabolic; extend to mortars + @assert nmortars(dg, cache) == 0 + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" calc_surface_integral!( + du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic) + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!( + du, mesh, equations_parabolic, dg, cache_parabolic) + + return nothing +end + +# Transform solution variables prior to taking the gradient +# (e.g., conservative to primitive variables). Defaults to doing nothing. +# TODO: can we avoid copying data? +function transform_variables!(u_transformed, u, mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations_parabolic, dg, i, j, k, element) + u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, equations_parabolic) + set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, k, element) + end + end +end + +# This is the version used when calculating the divergence of the viscous fluxes +function calc_volume_integral!(du, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + flux_1_node = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, k, element) + flux_2_node = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, k, element) + flux_3_node = get_node_vars(flux_viscous_z, equations_parabolic, dg, i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux_1_node, equations_parabolic, dg, ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[jj, j], flux_2_node, equations_parabolic, dg, i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[kk, k], flux_3_node, equations_parabolic, dg, i, j, kk, element) + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +# We pass the `surface_integral` argument solely for dispatch +function prolong2interfaces!(cache_parabolic, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache) + @unpack interfaces = cache_parabolic + @unpack orientations = interfaces + + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for interface in eachinterface(dg, cache) + left_element = interfaces.neighbor_ids[1, interface] + right_element = interfaces.neighbor_ids[2, interface] + + if orientations[interface] == 1 + # interface in x-direction + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j, k, left_element] + interfaces.u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k, right_element] + end + elseif orientations[interface] == 2 + # interface in y-direction + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg), k, left_element] + interfaces.u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k, right_element] + end + else # if orientations[interface] == 3 + # interface in z-direction + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + interfaces.u[1, v, i, j, interface] = flux_viscous_z[v, i, j, nnodes(dg), left_element] + interfaces.u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1, right_element] + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +function calc_interface_flux!(surface_flux_values, + mesh::TreeMesh{3}, equations_parabolic, + dg::DG, cache_parabolic) + @unpack neighbor_ids, orientations = cache_parabolic.interfaces + + @threaded for interface in eachinterface(dg, cache_parabolic) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + # orientation = 3: left -> 6, right -> 5 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for j in eachnode(dg), i in eachnode(dg) + # Get precomputed fluxes at interfaces + flux_ll, flux_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, + dg, i, j, interface) + + # Compute interface flux as mean of left and right viscous fluxes + # TODO: parabolic; only BR1 at the moment + flux = 0.5 * (flux_ll + flux_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, left_direction, left_id] = flux[v] + surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + end + end + end + + return nothing +end + + +# This is the version used when calculating the divergence of the viscous fluxes +function prolong2boundaries!(cache_parabolic, flux_viscous, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache) + @unpack boundaries = cache_parabolic + @unpack orientations, neighbor_sides = boundaries + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for boundary in eachboundary(dg, cache_parabolic) + element = boundaries.neighbor_ids[boundary] + + if orientations[boundary] == 1 + # boundary in x-direction + if neighbor_sides[boundary] == 1 + # element in -x direction of boundary + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg), j, k, element] + end + else # Element in +x direction of boundary + for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k, element] + end + end + elseif orientations[boundary] == 2 + # boundary in y-direction + if neighbor_sides[boundary] == 1 + # element in -y direction of boundary + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, i, k, boundary] = flux_viscous_y[v, i, nnodes(dg), k, element] + end + else + # element in +y direction of boundary + for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k, element] + end + end + else # if orientations[boundary] == 3 + # boundary in z-direction + if neighbor_sides[boundary] == 1 + # element in -z direction of boundary + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[1, v, i, j, boundary] = flux_viscous_z[v, i, j, nnodes(dg), element] + end + else + # element in +z direction of boundary + for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + boundaries.u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1, element] + end + end + end + end + + return nothing +end + + +function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache, cache_parabolic) + gradients_x, gradients_y, gradients_z = gradients + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays + + @threaded for element in eachelement(dg, cache) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + # Get solution and gradients + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, element) + gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, k, element) + gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, k, element) + gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k, element) + + # Calculate viscous flux and store each component for later use + flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 1, equations_parabolic) + flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 2, equations_parabolic) + flux_viscous_node_z = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 3, equations_parabolic) + set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, i, j, k, element) + set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, i, j, k, element) + set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg, i, j, k, element) + end + end +end + + +# TODO: parabolic; decide if we should keep this. +function get_unsigned_normal_vector_3d(direction) + if direction > 6 || direction < 1 + error("Direction = $direction; in 3D, direction should be 1, 2, 3, 4, 5, or 6.") + end + if direction == 1 || direction == 2 + return SVector(1.0, 0.0, 0.0) + elseif direction == 3 || direction == 4 + return SVector(0.0, 1.0, 0.0) + else + return SVector(0.0, 0.0, 1.0) + end +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[5], + equations_parabolic, surface_integral, dg, cache, + 5, firsts[5], lasts[5]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[6], + equations_parabolic, surface_integral, dg, cache, + 6, firsts[6], lasts[6]) +end + + +function calc_boundary_flux_by_direction_gradient!(surface_flux_values::AbstractArray{<:Any,5}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, j, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + end + + # TODO: revisit if we want more general boundary treatments. + # This assumes the gradient numerical flux at the boundary is the gradient variable, + # which is consistent with BR1, LDG. + flux_inner = u_inner + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, j, boundary) + flux = boundary_condition(flux_inner, u_inner, get_unsigned_normal_vector_3d(direction), + x, t, Gradient(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[5], + equations_parabolic, surface_integral, dg, cache, + 5, firsts[5], lasts[5]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[6], + equations_parabolic, surface_integral, dg, cache, + 6, firsts[6], lasts[6]) +end +function calc_boundary_flux_by_direction_divergence!(surface_flux_values::AbstractArray{<:Any,5}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + + # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") + # of the viscous flux, as computed in `prolong2boundaries!` + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Get viscous boundary fluxes + flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, j, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + flux_inner = flux_ll + else # Element is on the right, boundary on the left + flux_inner = flux_rr + end + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, j, boundary) + + # TODO: add a field in `cache.boundaries` for gradient information. UPDATE THIS COMMENT + # Here, we pass in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information. + # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion3D and + # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesDiffusion3D as of 2022-6-27. + # It will not work with implementations which utilize `u_inner` to impose boundary conditions. + flux = boundary_condition(flux_inner, nothing, get_unsigned_normal_vector_3d(direction), + x, t, Divergence(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + + +# Calculate the gradient of the transformed variables +function calc_gradient!(gradients, u_transformed, t, + mesh::TreeMesh{3}, equations_parabolic, + boundary_conditions_parabolic, dg::DG, cache, cache_parabolic) + + gradients_x, gradients_y, gradients_z = gradients + + # Reset du + @trixi_timeit timer() "reset gradients" begin + reset_du!(gradients_x, dg, cache) + reset_du!(gradients_y, dg, cache) + reset_du!(gradients_z, dg, cache) + end + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + @unpack derivative_dhat = dg.basis + @threaded for element in eachelement(dg, cache) + + # Calculate volume terms in one element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], u_node, equations_parabolic, dg, ii, j, k, element) + end + + for jj in eachnode(dg) + multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], u_node, equations_parabolic, dg, i, jj, k, element) + end + + for kk in eachnode(dg) + multiply_add_to_node_vars!(gradients_z, derivative_dhat[kk, k], u_node, equations_parabolic, dg, i, j, kk, element) + end + end + end + end + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!( + cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" begin + @unpack surface_flux_values = cache_parabolic.elements + @unpack neighbor_ids, orientations = cache_parabolic.interfaces + + @threaded for interface in eachinterface(dg, cache_parabolic) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + # orientation = 3: left -> 6, right -> 5 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for j in eachnode(dg), i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, + equations_parabolic, dg, i, j, interface) + flux = 0.5 * (u_ll + u_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, left_direction, left_id] = flux[v] + surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + end + end + end + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!( + cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!( + cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic, + dg.surface_integral, dg) + + # TODO: parabolic; mortars + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache_parabolic.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + factor_1 = boundary_interpolation[1, 1] + factor_2 = boundary_interpolation[nnodes(dg), 2] + @threaded for element in eachelement(dg, cache) + for m in eachnode(dg), l in eachnode(dg) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients_x[v, 1, l, m, element] = ( + gradients_x[v, 1, l, m, element] - surface_flux_values[v, l, m, 1, element] * factor_1) + + # surface at +x + gradients_x[v, nnodes(dg), l, m, element] = ( + gradients_x[v, nnodes(dg), l, m, element] + surface_flux_values[v, l, m, 2, element] * factor_2) + + # surface at -y + gradients_y[v, l, 1, m, element] = ( + gradients_y[v, l, 1, m, element] - surface_flux_values[v, l, m, 3, element] * factor_1) + + # surface at +y + gradients_y[v, l, nnodes(dg), m, element] = ( + gradients_y[v, l, nnodes(dg), m, element] + surface_flux_values[v, l, m, 4, element] * factor_2) + + # surface at -z + gradients_z[v, l, m, 1, element] = ( + gradients_z[v, l, m, 1, element] - surface_flux_values[v, l, m, 5, element] * factor_1) + + # surface at +z + gradients_z[v, l, m, nnodes(dg), element] = ( + gradients_z[v, l, m, nnodes(dg), element] + surface_flux_values[v, l, m, 6, element] * factor_2) + end + end + end + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" begin + apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic) + end + + return nothing +end + + +# This method is called when a SemidiscretizationHyperbolic is constructed. +# It constructs the basic `cache` used throughout the simulation to compute +# the RHS etc. +function create_cache_parabolic(mesh::TreeMesh{3}, equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, RealT, uEltype) + # Get cells for which an element needs to be created (i.e. all leaf cells) + leaf_cell_ids = local_leaf_cells(mesh.tree) + + elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) + + n_vars = nvariables(equations_hyperbolic) + n_nodes = nnodes(elements) + n_elements = nelements(elements) + u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) + gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) + flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) + + interfaces = init_interfaces(leaf_cell_ids, mesh, elements) + + boundaries = init_boundaries(leaf_cell_ids, mesh, elements) + + # mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar) + + # cache = (; elements, interfaces, boundaries, mortars) + cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + + # Add specialized parts of the cache required to compute the mortars etc. + # cache = (;cache..., create_cache(mesh, equations_parabolic, dg.mortar, uEltype)...) + + return cache +end + + +# Needed to *not* flip the sign of the inverse Jacobian. +# This is because the parabolic fluxes are assumed to be of the form +# `du/dt + df/dx = dg/dx + source(x,t)`, +# where f(u) is the inviscid flux and g(u) is the viscous flux. +function apply_jacobian!(du, mesh::TreeMesh{3}, + equations::AbstractEquationsParabolic, dg::DG, cache) + + @threaded for element in eachelement(dg, cache) + factor = cache.elements.inverse_jacobian[element] + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + du[v, i, j, k, element] *= factor + end + end + end + + return nothing +end + +end # @muladd diff --git a/test/runtests.jl b/test/runtests.jl index 79a8e54e791..d16ec13256b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -88,6 +88,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" include("test_parabolic_2d.jl") + include("test_parabolic_3d.jl") end @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index c08dc57be51..fbb23695c38 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -5,13 +5,12 @@ using Trixi include("test_trixi.jl") -EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") # Start with a clean environment: remove Trixi output directory if it exists outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -@testset "SemidiscretizationHyperbolicParabolic" begin +@testset "SemidiscretizationHyperbolicParabolic (2D)" begin @trixi_testset "DGMulti 2D rhs_parabolic!" begin diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl new file mode 100644 index 00000000000..c0728bfbb8c --- /dev/null +++ b/test/test_parabolic_3d.jl @@ -0,0 +1,86 @@ +module TestExamplesParabolic3D + +using Test +using Trixi + +include("test_trixi.jl") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +@testset "SemidiscretizationHyperbolicParabolic (3D)" begin + + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), + l2 = [0.0005532847115849239, 0.000659263490965341, 0.0007776436127362806, 0.0006592634909662951, 0.0038073628897809185], + linf = [0.0017039861523615585, 0.002628561703560073, 0.003531057425112172, 0.0026285617036090336, 0.015587829540351095] + ) + end + + @trixi_testset "DGMulti: elixir_navierstokes_taylor_green_vortex.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_taylor_green_vortex.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.25), + l2 = [0.0001825713444029892, 0.015589736382772248, 0.015589736382771884, 0.021943924667273653, 0.01927370280244222], + linf = [0.0006268463584697681, 0.03218881662749007, 0.03218881662697948, 0.053872495395614256, 0.05183822000984151] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.0019582188528512257, 0.002653449504302844, 0.002898264205184629, 0.002653449504302853, 0.009511572365085706], + linf = [0.013680656759085918, 0.0356910450154318, 0.023526343547736236, 0.035691045015431855, 0.11482570604041165] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.00195468651965362, 0.0026554367897028506, 0.002892730402724066, 0.002655436789702817, 0.009596351796609566], + linf = [0.013680508110645473, 0.035673446359424356, 0.024024936779729028, 0.03567344635942474, 0.11839497110809383] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (Entropy gradient variables)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + l2 = [0.0019770444875099307, 0.0026524750946399327, 0.00290860030832445, 0.0026524750946399396, 0.009509568981439294], + linf = [0.01387936112914212, 0.03526260609304053, 0.023554197097368997, 0.035262606093040896, 0.11719963716509518] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.001974631423398113, 0.002654768259143932, 0.002907031063651286, 0.002654768259143901, 0.009587792882971452], + linf = [0.01387919380137137, 0.035244084526358944, 0.02398614622061363, 0.03524408452635828, 0.12005056512506407] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_convergence.jl (flux differencing)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + volume_integral=VolumeIntegralFluxDifferencing(flux_central), + l2 = [0.0019582188528180213, 0.002653449504301736, 0.0028982642051960006, 0.0026534495043017384, 0.009511572364811033], + linf = [0.013680656758949583, 0.035691045015224444, 0.02352634354676752, 0.035691045015223424, 0.11482570603751441] + ) + end + + @trixi_testset "TreeMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.25), + l2 = [0.00024173250389635442, 0.015684268393762454, 0.01568426839376248, 0.021991909545192333, 0.02825413672911425], + linf = [0.0008410587892853094, 0.04740176181772552, 0.04740176181772507, 0.07483494924031157, 0.150181591534448] + ) + end + +end + +# Clean up afterwards: delete Trixi output directory +@test_nowarn isdir(outdir) && rm(outdir, recursive=true) + +end # module