From a090746fff16d6b83e98462ad195463b3abf9e3e Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 9 Aug 2022 22:10:39 +0200 Subject: [PATCH] Add support for hyperbolic-parabolic systems to TreeMesh2D (#1156) * Add prototype for advection-diffusion elixir for TreeMesh2D * Initial implementation of `calc_gradient!` for TreeMesh{2} * First complete implementation of the parabolic rhs operator for TreeMesh * Make it work * Add first elixir that works * Clean up elixir * first draft of container for Navier-Stokes constants and fluxes * remove unneeded temperature computation * draft of elixir with possible boundary condition structure * added manufactured solution and source term * fix typo in function name for MMS * update variable names for consistency. improve comments * fix dumb typos in new equation system name * actually export new equations * add comment near variable_mapping declaration. * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * parabolic equations now exported separately * change name to CompressibleNavierStokesEquations2D * export NS with proper name * explicitly require compressible Euler in the type parameter * name kinematic viscosity nu for consistency with Lattice-Boltzmann * add promotion in constructor * make Reynolds, Prandtl, Mach, and kappa keyword arguments * update constructor call in elixir * reduce computation by exploiting stress tensor symmetry * fix unpacking of flux * modifying parabolic cache creation in cache, we assume we take the gradient of all hyperbolic variables. since the number of parabolic variables can differ from the number of hyperbolic variables, we pass in the hyperbolic equations to `create_cache_parabolic` now * comments * comments * formatting and renaming equations to equations_hyperbolic formatting comments * fix unpacking of gradients in flux * adding CNS BCs * adding lid-driven cavity elixir * adding variable transform, editing cons2prim for CNS * add prim2cons for NS (inconsistent right now) * add draft of DGMulti Navier Stokes convergence elixir * converging solution using elixir for TreeMesh with BCs * fixing DGMulti advection diffusion elixir convergence * naming equations as parabolic/hyperbolic * generalizing transform_variables * add TODO more todos * additional checks on get_unsigned_normal * adding isothermal BC * commenting out unused CNS functions for now * fix call to transform_variables * comments and cleanup * changing default solver and Re for cavity * adding more advection diffusion tests * label tests * add gradient_variables field to SemidiscretizationHyperbolicParabolic * Revert "add gradient_variables field to SemidiscretizationHyperbolicParabolic" This reverts commit 063b602ebdb06a09be0f1cb34d5f24cdc86c6a83. * allowing for specialization in transform_variables adding a function `gradient_variable_transformation` which should get specialized if the gradient variables are not conservative variables * formatting and comments * reverting elixir * comments * standardizing time tol * minor fix to CNS boundary flux for convenience make it so that the density state is computed correctly, even though it's not used * formatting + comments * using primitive variables in viscous flux instead of conservative * minor formatting * add CNS tests * fix test * testing if scoping issues fix TreeMesh tests * decreasing timestep tol for TreeMesh parabolic test * enabling periodic BCs for TreeMesh + more tests * fix density BC flux (unused, but could be useful) * adding non-working TreeMesh elixirs * adding AD checks * standardizing parameters in convergence elixirs * minor cleanup * revert DGMulti CNS elixir setup back to the one used in tests * adding TreeMesh CNS convergence test * removing redundant elixir * add more tests * add more test * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * set version to v0.4.43 * set development version to v0.4.44-pre * add YouTube links to JuliaCon presentations (#1192) * add YouTube links to JuliaCon presentations * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Michael Schlottke-Lakemper * improved readability of equation docstrings in compressible Euler files * Navier-Stokes in 2D on DGMulti and TreeMesh (#1165) * first draft of container for Navier-Stokes constants and fluxes * remove unneeded temperature computation * draft of elixir with possible boundary condition structure * added manufactured solution and source term * fix typo in function name for MMS * update variable names for consistency. improve comments * fix dumb typos in new equation system name * actually export new equations * add comment near variable_mapping declaration. * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * parabolic equations now exported separately * change name to CompressibleNavierStokesEquations2D * export NS with proper name * explicitly require compressible Euler in the type parameter * name kinematic viscosity nu for consistency with Lattice-Boltzmann * add promotion in constructor * make Reynolds, Prandtl, Mach, and kappa keyword arguments * update constructor call in elixir * reduce computation by exploiting stress tensor symmetry * fix unpacking of flux * modifying parabolic cache creation in cache, we assume we take the gradient of all hyperbolic variables. since the number of parabolic variables can differ from the number of hyperbolic variables, we pass in the hyperbolic equations to `create_cache_parabolic` now * comments * comments * formatting and renaming equations to equations_hyperbolic formatting comments * fix unpacking of gradients in flux * adding CNS BCs * adding lid-driven cavity elixir * adding variable transform, editing cons2prim for CNS * add prim2cons for NS (inconsistent right now) * add draft of DGMulti Navier Stokes convergence elixir * converging solution using elixir for TreeMesh with BCs * fixing DGMulti advection diffusion elixir convergence * naming equations as parabolic/hyperbolic * generalizing transform_variables * add TODO more todos * additional checks on get_unsigned_normal * adding isothermal BC * commenting out unused CNS functions for now * fix call to transform_variables * comments and cleanup * changing default solver and Re for cavity * adding more advection diffusion tests * label tests * add gradient_variables field to SemidiscretizationHyperbolicParabolic * Revert "add gradient_variables field to SemidiscretizationHyperbolicParabolic" This reverts commit 063b602ebdb06a09be0f1cb34d5f24cdc86c6a83. * allowing for specialization in transform_variables adding a function `gradient_variable_transformation` which should get specialized if the gradient variables are not conservative variables * formatting and comments * reverting elixir * comments * standardizing time tol * minor fix to CNS boundary flux for convenience make it so that the density state is computed correctly, even though it's not used * formatting + comments * using primitive variables in viscous flux instead of conservative * minor formatting * add CNS tests * fix test * testing if scoping issues fix TreeMesh tests * decreasing timestep tol for TreeMesh parabolic test * enabling periodic BCs for TreeMesh + more tests * fix density BC flux (unused, but could be useful) * adding non-working TreeMesh elixirs * adding AD checks * standardizing parameters in convergence elixirs * minor cleanup * revert DGMulti CNS elixir setup back to the one used in tests * adding TreeMesh CNS convergence test * removing redundant elixir * add more tests * add more test * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * add docstrings Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan * add docstring for CNS equations * removing some addressed TODOs * adjust definition of the kappa constant * avoid overwriting u_grad with viscous_flux results * remove old TODOs * clarify TODOs for later * removing let statements * rearrange main docstring. Fix typos in math environment causing it to not parse * update TODO notes * Apply suggestions from code review Co-authored-by: Andrew Winters * update TODO notes * Apply suggestions from code review Co-authored-by: Andrew Winters * update TODO notes * changing diffusivity names Auto stash before merge of "parabolic-treemesh" and "origin/parabolic-treemesh" * Reynolds/Prandtl number name changes for consistency * Apply suggestions from code review Co-authored-by: Hendrik Ranocha Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Andrew Winters * fixing bug introduced by GH code review * moving manufactured solution into the CNS equations file * fix allocation issue * get_diffusivity -> diffusivity * removing cruft code, adding comment on messy prim2cons routine * moving manufactured solution back into elixirs * Update examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl * TODO -> Todo, and .6 -> 0.6 d * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl Co-authored-by: Hendrik Ranocha * adding note to Adiabatic/Isothermal BCs * replacing CompressibleNavierStokesEquations with ...Diffusion * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl Co-authored-by: Hendrik Ranocha * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl Co-authored-by: Hendrik Ranocha * documenting Jacobian sign flip * remove extra arg to `gradient_variable_transformation` * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper * change warning to error * fix bug introduced by search/replace * dg_parabolic -> parabolic_scheme * BoundaryConditionViscousWall -> BoundaryConditionNavierStokesWall * u_grad -> gradients * fix test * renaming * update comment * converting to signature `flux(u, gradients, orientation, equations)` * using prolong2interfaces/boundaries * unpacking `gradients` and `flux_viscous` * using calc_surface_integral! * adding @muladd * Apply suggestions from code review * update comments * Rename grad_u -> gradients * Update src/equations/compressible_navier_stokes_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Rename remaining grad_u/u_grad -> gradients * Refactor into prolong2interfaces! * Refactor calc_volume_integral! * Refactor calc_interface_flux * Refactor prolong2boundaries! * Refactor calc_divergence! * Formatting consistency * Remove dispatch on volume integral type * Add comment on why surface_integral is passed to prolong2interfaces! * Explicitly use weak form volume integral to allow easy overriding in test * Add flux differencing test for CNS test * Remove erroneous P4estMesh{2} dispatch * Remove non-cons terms from parabolic solver Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan Co-authored-by: Jesse Chan --- .../dgmulti_2d/elixir_advection_diffusion.jl | 2 +- .../elixir_advection_diffusion_nonperiodic.jl | 72 +++ .../elixir_navier_stokes_convergence.jl | 206 ++++++ .../elixir_navier_stokes_lid_driven_cavity.jl | 74 +++ .../elixir_advection_diffusion.jl | 83 +++ .../elixir_advection_diffusion_nonperiodic.jl | 89 +++ .../elixir_navier_stokes_convergence.jl | 214 +++++++ .../elixir_navier_stokes_lid_driven_cavity.jl | 81 +++ src/Trixi.jl | 7 +- src/equations/compressible_euler_1d.jl | 4 +- src/equations/compressible_euler_2d.jl | 6 +- src/equations/compressible_euler_3d.jl | 8 +- .../compressible_euler_multicomponent_1d.jl | 10 +- .../compressible_euler_multicomponent_2d.jl | 12 +- .../compressible_navier_stokes_2d.jl | 308 +++++++++ src/equations/equations_parabolic.jl | 10 +- src/equations/laplace_diffusion_2d.jl | 48 +- ...semidiscretization_hyperbolic_parabolic.jl | 3 +- src/solvers/dgmulti/dg.jl | 2 + src/solvers/dgmulti/dg_parabolic.jl | 129 ++-- src/solvers/dgmulti/sbp.jl | 1 + src/solvers/dgsem_p4est/dg_2d.jl | 1 + src/solvers/dgsem_p4est/dg_3d.jl | 1 + src/solvers/dgsem_tree/dg.jl | 1 + src/solvers/dgsem_tree/dg_1d.jl | 1 + src/solvers/dgsem_tree/dg_2d.jl | 2 + src/solvers/dgsem_tree/dg_2d_parabolic.jl | 585 ++++++++++++++++++ src/solvers/dgsem_tree/dg_3d.jl | 1 + src/solvers/dgsem_unstructured/dg_2d.jl | 1 + test/test_parabolic_2d.jl | 97 ++- 30 files changed, 1937 insertions(+), 122 deletions(-) create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl create mode 100644 examples/dgmulti_2d/elixir_navier_stokes_convergence.jl create mode 100644 examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl create mode 100644 examples/tree_2d_dgsem/elixir_advection_diffusion.jl create mode 100644 examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl create mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl create mode 100644 src/equations/compressible_navier_stokes_2d.jl create mode 100644 src/solvers/dgsem_tree/dg_2d_parabolic.jl diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 21effda986e..1273e25859e 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -5,7 +5,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.5, 1.0) -equations_parabolic = LaplaceDiffusion2D(5e-2, equations) +equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..b6adf242299 --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl @@ -0,0 +1,72 @@ +using Trixi, OrdinaryDiffEq + +dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +diffusivity() = 5.0e-2 + +equations = LinearScalarAdvectionEquation2D(1.0, 0.0) +equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) + +# Example setup taken from +# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016). +# Robust DPG methods for transient convection-diffusion. +# In: Building bridges: connections and challenges in modern approaches +# to numerical partial differential equations. +# [DOI](https://doi.org/10.1007/978-3-319-41640-3_6). +function initial_condition_erikkson_johnson(x, t, equations) + l = 4 + epsilon = diffusivity() # Note: this requires epsilon < 0.6 due to the sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_erikkson_johnson + +# tag different boundary segments +left(x, tol=50*eps()) = abs(x[1] + 1) < tol +right(x, tol=50*eps()) = abs(x[1]) < tol +bottom(x, tol=50*eps()) = abs(x[2] + .5) < tol +top(x, tol=50*eps()) = abs(x[2] - .5) < tol +entire_boundary(x, tol=50*eps()) = true +is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom, + :entire_boundary => entire_boundary) +mesh = DGMultiMesh(dg; coordinates_min=(-1.0, -0.5), coordinates_max=(0.0, 0.5), + cells_per_dimension=(16, 16), is_on_boundary) + +# BC types +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +# define inviscid boundary conditions, enforce "do nothing" boundary condition at the outflow +boundary_conditions = (; :left => boundary_condition, + :top => boundary_condition, + :bottom => boundary_condition, + :right => boundary_condition_do_nothing) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :entire_boundary => boundary_condition) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) + +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl new file mode 100644 index 00000000000..4995acccb78 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -0,0 +1,206 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +reynolds_number() = 100 +prandtl_number() = 0.72 + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(), + Mach_freestream=0.5) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); periodicity=(true, false), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + y = x[2] + + # 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() + Re = reynolds_number() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_Re = 1.0 / Re + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * inv_Re + + 2.0 / 3.0 * v2_xy * inv_Re + - v1_yy * inv_Re + - v2_xy * inv_Re ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + # stress tensor from y-direction + - v1_xy * inv_Re + - v2_xx * inv_Re + - 4.0 / 3.0 * v2_yy * inv_Re + + 2.0 / 3.0 * v1_xy * inv_Re ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * inv_Re + + 2.0 / 3.0 * v2_xy * v1 * inv_Re + - 4.0 / 3.0 * v1_x * v1_x * inv_Re + + 2.0 / 3.0 * v2_y * v1_x * inv_Re + - v1_xy * v2 * inv_Re + - v2_xx * v2 * inv_Re + - v1_y * v2_x * inv_Re + - v2_x * v2_x * inv_Re + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * inv_Re + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * inv_Re + - v2_xy * v1 * inv_Re + - v1_y * v1_y * inv_Re + - v2_x * v1_y * inv_Re + - 4.0 / 3.0 * v2_yy * v2 * inv_Re + + 2.0 / 3.0 * v1_xy * v2 * inv_Re + - 4.0 / 3.0 * v2_y * v2_y * inv_Re + + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * inv_Re ) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +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, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..be8ad1a0730 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl @@ -0,0 +1,74 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +reynolds_number() = 1000.0 +prandtl_number() = 0.72 +mach_number() = 0.1 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), + Prandtl=prandtl_number(), + Mach_freestream=mach_number()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) + +top(x, tol=50*eps()) = abs(x[2] - 1) < tol +rest_of_boundary(x, tol=50*eps()) = !top(x, tol) +is_on_boundary = Dict(:top => top, :rest_of_boundary => rest_of_boundary) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = mach_number() + rho = 1.0 + u, v = 0, 0 + p = 1.0 / (Ma^2 * equations.gamma) + return prim2cons(SVector(rho, u, v, p), equations) +end +initial_condition = initial_condition_cavity + +# BC types +velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0)) +velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc) +boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc) + +# define inviscid boundary conditions +boundary_conditions = (; :top => boundary_condition_slip_wall, + :rest_of_boundary => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top => boundary_condition_lid, + :rest_of_boundary => boundary_condition_cavity) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl new file mode 100644 index 00000000000..e538f701919 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl @@ -0,0 +1,83 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +advection_velocity = (1.5, 1.0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) +diffusivity() = 5.0e-2 +equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=true, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Define initial condition +function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) + # Store translated coordinate for easy use of exact solution + x_trans = x - equation.advection_velocity * t + + nu = diffusivity() + c = 1.0 + A = 0.5 + L = 2 + f = 1/L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) + return SVector(scalar) +end +initial_condition = initial_condition_diffusive_convergence_test + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_periodic +boundary_conditions_parabolic = boundary_condition_periodic + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.5 +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +time_int_tol = 1.0e-11 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + +# Print the timer summary +summary_callback() diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..1f0bfb788ad --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -0,0 +1,89 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +diffusivity() = 5.0e-2 +advection_velocity = (1.0, 0.0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) +equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -0.5) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 0.0, 0.5) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Example setup taken from +# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016). +# Robust DPG methods for transient convection-diffusion. +# In: Building bridges: connections and challenges in modern approaches +# to numerical partial differential equations. +# [DOI](https://doi.org/10.1007/978-3-319-41640-3_6). +function initial_condition_erikkson_johnson(x, t, equations) + l = 4 + epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_erikkson_johnson + +boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition), + y_neg = BoundaryConditionDirichlet(initial_condition), + y_pos = BoundaryConditionDirichlet(initial_condition), + x_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +time_int_tol = 1.0e-11 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + +# Print the timer summary +summary_callback() diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl new file mode 100644 index 00000000000..173c01365c1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -0,0 +1,214 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +reynolds_number() = 100 +prandtl_number() = 0.72 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(), + Mach_freestream=0.5) + +# 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) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=(true, false), + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + y = x[2] + + # 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() + Re = reynolds_number() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_Re = 1.0 / Re + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * inv_Re + + 2.0 / 3.0 * v2_xy * inv_Re + - v1_yy * inv_Re + - v2_xy * inv_Re ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + # stress tensor from y-direction + - v1_xy * inv_Re + - v2_xx * inv_Re + - 4.0 / 3.0 * v2_yy * inv_Re + + 2.0 / 3.0 * v1_xy * inv_Re ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * inv_Re + + 2.0 / 3.0 * v2_xy * v1 * inv_Re + - 4.0 / 3.0 * v1_x * v1_x * inv_Re + + 2.0 / 3.0 * v2_y * v1_x * inv_Re + - v1_xy * v2 * inv_Re + - v2_xx * v2 * inv_Re + - v1_y * v2_x * inv_Re + - v2_x * v2_x * inv_Re + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * inv_Re + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * inv_Re + - v2_xy * v1 * inv_Re + - v1_y * v1_y * inv_Re + - v2_x * v1_y * inv_Re + - 4.0 / 3.0 * v2_yy * v2 * inv_Re + + 2.0 / 3.0 * v1_xy * v2 * inv_Re + - 4.0 / 3.0 * v2_y * v2_y * inv_Re + + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * inv_Re ) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +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) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_top_bottom, + y_pos = boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +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, dt = 1e-5, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary + diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..bf40aaa41df --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl @@ -0,0 +1,81 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +reynolds_number() = 1000.0 +prandtl_number() = 0.72 +mach_number() = 0.1 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), + Prandtl=prandtl_number(), + Mach_freestream=mach_number()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = mach_number() + rho = 1.0 + u, v = 0, 0 + p = 1.0 / (Ma^2 * equations.gamma) + return prim2cons(SVector(rho, u, v, p), equations) +end +initial_condition = initial_condition_cavity + +# BC types +velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0)) +velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc) +boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc) + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_slip_wall + +boundary_conditions_parabolic = (; x_neg = boundary_condition_cavity, + y_neg = boundary_condition_cavity, + y_pos = boundary_condition_lid, + x_pos = boundary_condition_cavity) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 25.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=100) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary + + diff --git a/src/Trixi.jl b/src/Trixi.jl index 2533f71ea9e..e6e4e6da692 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -131,10 +131,12 @@ export AcousticPerturbationEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, InviscidBurgersEquation1D, - LaplaceDiffusion2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D +export LaplaceDiffusion2D, + CompressibleNavierStokesDiffusion2D + export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, @@ -161,7 +163,8 @@ export boundary_condition_do_nothing, BoundaryConditionNeumann, boundary_condition_noslip_wall, boundary_condition_slip_wall, - boundary_condition_wall + boundary_condition_wall, + BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 89ec07558bc..400650ed885 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -10,12 +10,12 @@ The compressible Euler equations ```math -\partial t +\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho e \end{pmatrix} + -\partial x +\frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ (\rho e +p) v_1 \end{pmatrix} diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 9c3d94cffd8..35852c0d1f9 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -10,17 +10,17 @@ The compressible Euler equations ```math -\partial t +\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \end{pmatrix} + -\partial x +\frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \end{pmatrix} + -\partial y +\frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \end{pmatrix} diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index d919c600112..1b93674c8a2 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -10,22 +10,22 @@ The compressible Euler equations ```math -\partial t +\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho e \end{pmatrix} + -\partial x +\frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ \rho v_1 v_3 \\ ( \rho e +p) v_1 \end{pmatrix} + -\partial y +\frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ \rho v_1 v_3 \\ ( \rho e +p) v_2 \end{pmatrix} + -\partial z +\frac{\partial}{\partial z} \begin{pmatrix} \rho v_3 \\ \rho v_1 v_3 \\ \rho v_2 v_3 \\ \rho v_3^2 + p \\ ( \rho e +p) v_3 \end{pmatrix} diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 39d8eb34f24..c5a3579ab3e 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -10,12 +10,12 @@ Multicomponent version of the compressible Euler equations ```math -\partial t +\frac{\partial}{\partial t} \begin{pmatrix} \rho v_1 \\ \rho e \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n} \end{pmatrix} + -\partial x +\frac{\partial}{\partial x} \begin{pmatrix} \rho v_1^2 + p \\ (\rho e +p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1 \end{pmatrix} @@ -25,15 +25,15 @@ Multicomponent version of the compressible Euler equations 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} ``` -for calorically perfect gas in one space dimension. -Here, ``\rho_i`` is the density of component ``i``, ``\rho=\sum_{i=1}^n\rho_i`` the sum of the individual ``\rho_i``, +for calorically perfect gas in one space dimension. +Here, ``\rho_i`` is the density of component ``i``, ``\rho=\sum_{i=1}^n\rho_i`` the sum of the individual ``\rho_i``, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho v_1^2 \right) ``` the pressure, ```math -\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}} +\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}} ``` total heat capacity ratio, ``\gamma_i`` heat capacity ratio of component ``i``, ```math diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index ad90969cf80..bb91cfbcb4e 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -10,17 +10,17 @@ Multicomponent version of the compressible Euler equations ```math -\partial t +\frac{\partial}{\partial t} \begin{pmatrix} \rho v_1 \\ \rho v_2 \\ \rho e \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n} \end{pmatrix} + -\partial x +\frac{\partial}{\partial x} \begin{pmatrix} \rho v_1^2 + p \\ \rho v_1 v_2 \\ ( \rho e +p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1 \end{pmatrix} + -\partial y +\frac{\partial}{\partial y} \begin{pmatrix} \rho v_1 v_2 \\ \rho v_2^2 + p \\ ( \rho e +p) v_2 \\ \rho_1 v_2 \\ \rho_2 v_2 \\ \vdots \\ \rho_{n} v_2 \end{pmatrix} @@ -29,15 +29,15 @@ Multicomponent version of the compressible Euler equations 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} ``` -for calorically perfect gas in two space dimensions. -Here, ``\rho_i`` is the density of component ``i``, ``\rho=\sum_{i=1}^n\rho_i`` the sum of the individual ``\rho_i``, +for calorically perfect gas in two space dimensions. +Here, ``\rho_i`` is the density of component ``i``, ``\rho=\sum_{i=1}^n\rho_i`` the sum of the individual ``\rho_i``, ``v_1``, ``v_2`` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2 + v_2^2) \right) ``` the pressure, ```math -\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}} +\gamma=\frac{\sum_{i=1}^n\rho_i C_{v,i}\gamma_i}{\sum_{i=1}^n\rho_i C_{v,i}} ``` total heat capacity ratio, ``\gamma_i`` heat capacity ratio of component ``i``, ```math diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl new file mode 100644 index 00000000000..683ff1f8514 --- /dev/null +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -0,0 +1,308 @@ +@doc raw""" + CompressibleNavierStokesDiffusion2D(gamma, Re, Pr, Ma_inf, equations) + +These equations contain 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, +- `Re`: Reynolds number, +- `Pr`: Prandtl number, +- `Ma_inf`: free-stream Mach number +- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) + +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) \right) +``` +as the pressure. The terms on the right hand side of the system above +are built from the viscous stress tensor +```math +\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I} +``` +where ``\underline{I}`` is the ``2\times 2`` 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}} +``` + +In two spatial dimensions we require gradients for three quantities, e.g., +primitive quantities +```math +\nabla v_1,\, \nabla v_2,\, \nabla T +``` +or the entropy variables +```math +\nabla w_2,\, \nabla w_3,\, \nabla w_4 +``` +where +```math +w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} +``` + +For this particular scaling the vicosity is set internally to be μ = 1/Re. +Further, the nondimensionalization takes the density-temperature-sound speed as +the principle quantities such that +``` +rho_inf = 1.0 +T_ref = 1.0 +c_inf = 1.0 +p_inf = 1.0 / gamma +u_inf = Ma_inf +R = 1.0 / gamma +``` + +Other normalization strategies exist, see the reference below for details. +- Marc Montagnac (2013) + Variable Normalization (nondimensionalization and scaling) for Navier-Stokes + equations: a practical guide + [CERFACS Technical report](https://www.cerfacs.fr/~montagna/TR-CFD-13-77.pdf) +The scaling used herein is Section 4.5 of the reference. +""" +struct CompressibleNavierStokesDiffusion2D{RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4} + # 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 + Re::RealT # Reynolds number + Pr::RealT # Prandtl number + Ma_inf::RealT # free-stream Mach number + kappa::RealT # thermal diffusivity for Fick's law + + p_inf::RealT # free-stream pressure + u_inf::RealT # free-stream velocity + R::RealT # gas constant (depends on nondimensional scaling!) + + equations_hyperbolic::E # CompressibleEulerEquations2D +end + +function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; Reynolds, Prandtl, Mach_freestream) + gamma = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + Re, Pr, Ma = promote(Reynolds, Prandtl, Mach_freestream) + + # Under the assumption of constant Prandtl number the thermal conductivity + # constant is kappa = gamma μ R / ((gamma-1) Pr). + # Important note! Due to nondimensional scaling R = 1 / gamma, this constant + # simplifies slightly. Also, the factor of μ is accounted for later. + kappa = inv_gamma_minus_one / Pr + + # From the nondimensionalization discussed above set the remaining free-stream + # quantities + p_inf = 1 / gamma + u_inf = Mach_freestream + R = 1 / gamma + CompressibleNavierStokesDiffusion2D{typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, + Re, Pr, Ma, kappa, + p_inf, u_inf, R, + equations) +end + + +# TODO: parabolic +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion2D) = ("v1", "v2", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion2D) = ("w2", "w3", "w4") + +varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusion2D) = + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) + +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D) = cons2prim + +# Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2 +# of the paper by Svärd, Carpenter and Nordström +# "A stable high-order finite difference scheme for the compressible Navier–Stokes +# equations, far-field boundary conditions" +# Although these authors use a different nondimensionalization so some constants are different +# particularly for Fick's law. +# +# Note, could be generalized to use Sutherland's law to get the molecular and thermal +# diffusivity +function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion2D) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T) + # either computed directly or reverse engineered from the gradient of the entropy vairables + # by way of the `convert_gradient_variables` function + rho, v1, v2, _ = u + + # gradients contains derivatives of each hyperbolic variable + _, dv1dx, dv2dx, dTdx = gradients[1] + _, dv1dy, dv2dy, dTdy = gradients[2] + + # Components of viscous stress tensor + + # (4/3 * (v1)_x - 2/3 * (v2)_y) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy + # ((v1)_y + (v2)_x) + # stress tensor is symmetric + tau_12 = dv1dy + dv2dx # = tau_21 + # (4/3 * (v2)_y - 2/3 * (v1)_x) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx + + # Fick's law q = -kappa * grad(T); constant is kappa = gamma μ R / ((gamma-1) Pr) + # Important note! Due to nondimensional scaling R = 1 / gamma, so the + # temperature T in the gradient computation already contains a factor of gamma + q1 = equations.kappa * dTdx + q2 = equations.kappa * dTdy + + # kinematic viscosity is simply 1/Re for this nondimensionalization + mu = 1.0 / equations.Re + + if orientation == 1 + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = ( v1 * tau_11 + v2 * tau_12 + q1 ) * mu + + return SVector(f1, f2, f3, f4) + else # if 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 = ( v1 * tau_12 + v2 * tau_22 + q2 ) * mu + + return SVector(g1, g2, g3, g4) + end +end + + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion2D) + rho, rho_v1, rho_v2, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = temperature(u, equations) + + return SVector(rho, v1, v2, T) +end + +# This routine is required because `prim2cons` is called in `initial_condition`, which +# is called with `equations::CompressibleEulerEquations2D`. This means it is inconsistent +# with `cons2prim(..., ::CompressibleNavierStokesDiffusion2D)` as defined above. +# TODO: parabolic. Is there a way to clean this up? +@inline prim2cons(u, equations::CompressibleNavierStokesDiffusion2D) = + prim2cons(u, equations.equations_hyperbolic) + + +@inline function temperature(u, equations::CompressibleNavierStokesDiffusion2D) + rho, rho_v1, rho_v2, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) + T = p / (equations.R * rho) + return T +end + +# TODO: can we generalize this to MHD? +""" + struct BoundaryConditionNavierStokesWall + +Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. +The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended +to be boundary condition types such as the `NoSlip` velocity boundary condition and the +`Adiabatic` or `Isothermal` heat boundary condition. + +!!! warning "Experimental feature" + This is an experimental feature and may change in future releases. +""" +struct BoundaryConditionNavierStokesWall{V, H} + boundary_condition_velocity::V + boundary_condition_heat_flux::H +end + +""" + struct NoSlip + +Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` +should be a function with signature `boundary_value_function(x, t, equations)` +and should return a `SVector{NDIMS}` whose entries are the velocity vector at a +point `x` and time `t`. +""" +struct NoSlip{F} + boundary_value_function::F # value of the velocity vector on the boundary +end + +""" + struct Isothermal + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_function` should be a function with signature +`boundary_value_function(x, t, equations)` and return a scalar value for the +temperature at point `x` and time `t`. +""" +struct Isothermal{F} + boundary_value_function::F # value of the temperature on the boundary +end + +""" + struct Adiabatic + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_normal_flux_function` should be a function with signature +`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the +normal heat flux at point `x` and time `t`. +""" +struct Adiabatic{F} + boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion2D) + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + return SVector(u_inner[1], v1, v2, u_inner[4]) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion2D) + # rho, v1, v2, _ = u_inner + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + _, tau_1n, tau_2n, _ = flux_inner # extract fluxes for 2nd and 3rd equations + normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], normal_energy_flux) +end + + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion2D) + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) + return SVector(u_inner[1], v1, v2, T) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion2D) + return flux_inner +end + diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index fec0eed6f63..340cf0e4294 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -1,3 +1,11 @@ +# specify transformation of conservative variables prior to taking gradients. +# specialize this function to compute gradients e.g., of primitive variables instead of conservative +gradient_variable_transformation(::AbstractEquationsParabolic) = cons2cons + # Linear scalar diffusion for use in linear scalar advection-diffusion problems -abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_2d.jl") + +# Compressible Navier-Stokes equations +abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +include("compressible_navier_stokes_2d.jl") diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index c171741716b..2f1afe25a6d 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -4,51 +4,55 @@ `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` with diffusivity ``\kappa`` applied to each solution component defined by `equations`. """ -struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N} +struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusion{2, N} diffusivity::T - equations::E + equations_hyperbolic::E end -LaplaceDiffusion2D(diffusivity, equations) = - LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations) +LaplaceDiffusion2D(diffusivity, equations_hyperbolic) = + LaplaceDiffusion2D{typeof(equations_hyperbolic), nvariables(equations_hyperbolic), typeof(diffusivity)}(diffusivity, equations_hyperbolic) varnames(variable_mapping, equations_parabolic::LaplaceDiffusion2D) = - varnames(variable_mapping, equations_parabolic.equations) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) # no orientation specified since the flux is vector-valued -function flux(u, grad_u, equations::LaplaceDiffusion2D) - dudx, dudy = grad_u - return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy) +function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion2D) + dudx, dudy = gradients + if orientation == 1 + return SVector(equations_parabolic.diffusivity * dudx) + else # if orientation == 2 + return SVector(equations_parabolic.diffusivity * dudy) + end end -# TODO: should this remain in the equations file, be moved to solvers, or live in the elixir? +# TODO: parabolic; should this remain in the equations file, be moved to solvers, or live in the elixir? # The penalization depends on the solver, but also depends explicitly on physical parameters, # and would probably need to be specialized for every different equation. -function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) - return dg.penalty_parameter * (u_outer - u_inner) * equations.diffusivity * inv_h +function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity * inv_h end # Dirichlet-type boundary condition for use with a parabolic solver in weak form -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, - equations::LaplaceDiffusion2D) - return boundary_condition.boundary_value_function(x, t, equations) + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_value_function(x, t, equations_parabolic) end -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Divergence, - equations::LaplaceDiffusion2D) - return u_inner + equations_parabolic::LaplaceDiffusion2D) + return flux_inner end -@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Divergence, - equations::LaplaceDiffusion2D) - return boundary_condition.boundary_normal_flux_function(x, t, equations) + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic) end -@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, x, t, operator_type::Gradient, - equations::LaplaceDiffusion2D) + equations_parabolic::LaplaceDiffusion2D) return flux_inner end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index ea96aa5c363..52a48d7aa03 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -102,7 +102,8 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache) - cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, solver_parabolic, RealT, uEltype)..., + cache_parabolic = (; create_cache_parabolic(mesh, equations, equations_parabolic, + solver, solver_parabolic, RealT, uEltype)..., initial_cache_parabolic...) SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index c4abe2e91b7..e6caf61c168 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -158,6 +158,7 @@ function max_dt(u, t, mesh::DGMultiMesh, end # interpolates from solution coefficients to face quadrature points +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::DGMultiMesh, equations, surface_integral, dg::DGMulti) rd = dg.basis @@ -259,6 +260,7 @@ function calc_surface_integral!(du, u, surface_integral::SurfaceIntegralWeakForm end # Specialize for nodal SBP discretizations. Uses that Vf*u = u[Fmask,:] +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::DGMultiMesh, equations, surface_integral, dg::DGMultiSBP) rd = dg.basis diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 9578e3f9eff..50cfd8ab17d 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,6 +1,10 @@ -function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - dg::DGMulti, dg_parabolic, RealT, uEltype) - nvars = nvariables(equations) +function create_cache_parabolic(mesh::DGMultiMesh, + equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, + dg::DGMulti, parabolic_scheme, RealT, uEltype) + # default to taking derivatives of all hyperbolic terms + # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache + nvars = nvariables(equations_hyperbolic) @unpack M, Drst = dg.basis weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst) @@ -8,15 +12,15 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsP # u_transformed stores "transformed" variables for computing the gradient @unpack md = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) - u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) - viscous_flux = similar.(u_grad) + gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) + flux_viscous = similar.(gradients) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) scalar_flux_face_values = similar(u_face_values) - grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) + gradients_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()] - local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] + local_flux_viscous_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] # precompute 1 / h for penalty terms @@ -28,28 +32,30 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsP end end - return (; u_transformed, u_grad, viscous_flux, + return (; u_transformed, gradients, flux_viscous, weak_differentiation_matrices, inv_h, - u_face_values, grad_u_face_values, scalar_flux_face_values, - local_u_values_threaded, local_viscous_flux_threaded, local_flux_face_values_threaded) + u_face_values, gradients_face_values, scalar_flux_face_values, + local_u_values_threaded, local_flux_viscous_threaded, local_flux_face_values_threaded) end # Transform solution variables prior to taking the gradient # (e.g., conservative to primitive variables). Defaults to doing nothing. # TODO: can we avoid copying data? -function transform_variables!(u_transformed, u, equations) +function transform_variables!(u_transformed, u, mesh, equations_parabolic::AbstractEquationsParabolic, + dg::DGMulti, parabolic_scheme, cache, cache_parabolic) @threaded for i in eachindex(u) - u_transformed[i] = u[i] + u_transformed[i] = gradient_variable_transformation(equations_parabolic)(u[i], equations_parabolic) end end # interpolates from solution coefficients to face quadrature points +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, surface_integral, dg::DGMulti, cache) apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end -function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, +function calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) @unpack local_flux_face_values_threaded = cache_parabolic @@ -60,27 +66,27 @@ function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, # compute flux * (nx, ny, nz) local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end - apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values) + apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(gradients[dim], :, e), local_flux_values) end end end -function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, +function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, boundary_conditions, dg::DGMulti, cache, cache_parabolic) @unpack weak_differentiation_matrices = cache_parabolic - for dim in eachindex(u_grad) - reset_du!(u_grad[dim], dg) + for dim in eachindex(gradients) + reset_du!(gradients[dim], dg) end # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) - dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: assumes mesh is affine + dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: DGMulti. Assumes mesh is affine here. apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), - view(u_grad[i], :, e), view(u, :, e)) + view(gradients[i], :, e), view(u, :, e)) end end @@ -97,15 +103,15 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, scalar_flux_face_values[idM] = 0.5 * (uP + uM) # TODO: use strong/weak formulation for curved meshes? end - calc_boundary_flux!(scalar_flux_face_values, nothing, t, Gradient(), boundary_conditions, + calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, mesh, equations, dg, cache, cache_parabolic) # compute surface contributions - calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, + calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, mesh, equations, dg, cache, cache_parabolic) for dim in eachdim(mesh) - invert_jacobian!(u_grad[dim], mesh, equations, dg, cache; scaling=1.0) + invert_jacobian!(gradients[dim], mesh, equations, dg, cache; scaling=1.0) end end @@ -154,52 +160,53 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, # for both the gradient and the divergence, the boundary flux is scalar valued. # for the gradient, it is the solution; for divergence, it is the normal flux. - u_boundary = boundary_condition(flux_face_values[fid,e], - face_normal, face_coordinates, t, - operator_type, equations) - flux_face_values[fid,e] = u_boundary + flux_face_values[fid,e] = boundary_condition(flux_face_values[fid,e], u_face_values[fid,e], + face_normal, face_coordinates, t, + operator_type, equations) end end return nothing end -function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, +function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) for dim in eachdim(mesh) - reset_du!(viscous_flux[dim], dg) + reset_du!(flux_viscous[dim], dg) end - @unpack local_viscous_flux_threaded, local_u_values_threaded = cache_parabolic + @unpack local_flux_viscous_threaded, local_u_values_threaded = cache_parabolic @threaded for e in eachelement(mesh, dg) # reset local storage for each element - local_viscous_flux = local_viscous_flux_threaded[Threads.threadid()] + local_flux_viscous = local_flux_viscous_threaded[Threads.threadid()] local_u_values = local_u_values_threaded[Threads.threadid()] fill!(local_u_values, zero(eltype(local_u_values))) for dim in eachdim(mesh) - fill!(local_viscous_flux[dim], zero(eltype(local_viscous_flux[dim]))) + fill!(local_flux_viscous[dim], zero(eltype(local_flux_viscous[dim]))) end - # interpolate u and gradient to quadrature points, store in `local_viscous_flux` - apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) # TODO: can we avoid this when we don't need it? + # interpolate u and gradient to quadrature points, store in `local_flux_viscous` + apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) # TODO: DGMulti. Specialize for nodal collocation methods (SBP, GaussSBP) for dim in eachdim(mesh) - apply_to_each_field(mul_by!(dg.basis.Vq), local_viscous_flux[dim], view(u_grad[dim], :, e)) + apply_to_each_field(mul_by!(dg.basis.Vq), local_flux_viscous[dim], view(gradients[dim], :, e)) end # compute viscous flux at quad points for i in eachindex(local_u_values) u_i = local_u_values[i] - u_grad_i = getindex.(local_viscous_flux, i) # TODO: check if this allocates. Shouldn't for tuples or SVector... - viscous_flux_i = flux(u_i, u_grad_i, equations) - setindex!.(local_viscous_flux, viscous_flux_i, i) + gradients_i = getindex.(local_flux_viscous, i) + for dim in eachdim(mesh) + flux_viscous_i = flux(u_i, gradients_i, dim, equations) + setindex!(local_flux_viscous[dim], flux_viscous_i, i) + end end # project back to the DG approximation space for dim in eachdim(mesh) - apply_to_each_field(mul_by!(dg.basis.Pq), view(viscous_flux[dim], :, e), local_viscous_flux[dim]) + apply_to_each_field(mul_by!(dg.basis.Pq), view(flux_viscous[dim], :, e), local_flux_viscous[dim]) end end end @@ -207,13 +214,13 @@ end # no penalization for a BR1 parabolic solver function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, boundary_conditions, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, - dg_parabolic::ViscousFormulationBassiRebay1, cache, cache_parabolic) + parabolic_scheme::ViscousFormulationBassiRebay1, cache, cache_parabolic) return nothing end function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, boundary_conditions, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, - dg_parabolic, cache, cache_parabolic) + parabolic_scheme, cache, cache_parabolic) # compute fluxes at interfaces @unpack scalar_flux_face_values, inv_h = cache_parabolic @unpack mapM, mapP = mesh.md @@ -221,15 +228,15 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, bounda idM, idP = mapM[face_node_index], mapP[face_node_index] uM, uP = u_face_values[idM], u_face_values[idP] inv_h_face = inv_h[face_node_index] - scalar_flux_face_values[idM] = scalar_flux_face_values[idM] + penalty(uP, uM, inv_h_face, equations, dg_parabolic) + scalar_flux_face_values[idM] = scalar_flux_face_values[idM] + penalty(uP, uM, inv_h_face, equations, parabolic_scheme) end return nothing end -function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh, +function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - boundary_conditions, dg::DGMulti, dg_parabolic, cache, cache_parabolic) + boundary_conditions, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) @unpack weak_differentiation_matrices = cache_parabolic @@ -240,14 +247,14 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), - view(du, :, e), view(viscous_flux[i], :, e)) + view(du, :, e), view(flux_viscous[i], :, e)) end end # interpolates from solution coefficients to face quadrature points - viscous_flux_face_values = cache_parabolic.grad_u_face_values + flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage for dim in eachdim(mesh) - prolong2interfaces!(viscous_flux_face_values[dim], viscous_flux[dim], mesh, equations, + prolong2interfaces!(flux_viscous_face_values[dim], flux_viscous[dim], mesh, equations, dg.surface_integral, dg, cache) end @@ -260,25 +267,28 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh # compute f(u, ∇u) ⋅ n flux_face_value = zero(eltype(scalar_flux_face_values)) for dim in eachdim(mesh) - uM = viscous_flux_face_values[dim][idM] - uP = viscous_flux_face_values[dim][idP] - # TODO: use strong/weak formulation? + uM = flux_viscous_face_values[dim][idM] + uP = flux_viscous_face_values[dim][idP] + # TODO: use strong/weak formulation to ensure stability on curved meshes? flux_face_value = flux_face_value + 0.5 * (uP + uM) * nxyzJ[dim][face_node_index] end scalar_flux_face_values[idM] = flux_face_value end - # TODO: decide what to pass in - calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(), + calc_boundary_flux!(scalar_flux_face_values, cache_parabolic.u_face_values, t, Divergence(), boundary_conditions, mesh, equations, dg, cache, cache_parabolic) calc_viscous_penalty!(scalar_flux_face_values, cache_parabolic.u_face_values, t, - boundary_conditions, mesh, equations, dg, dg_parabolic, + boundary_conditions, mesh, equations, dg, parabolic_scheme, cache, cache_parabolic) # surface contributions apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values) + # Note: we do not flip the sign of the geometric Jacobian here. + # 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. invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end @@ -290,21 +300,22 @@ end # boundary conditions will be applied to both grad(u) and div(u). function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions, source_terms, - dg::DGMulti, dg_parabolic, cache, cache_parabolic) + dg::DGMulti, parabolic_scheme, cache, cache_parabolic) reset_du!(du, dg) - @unpack u_transformed, u_grad, viscous_flux = cache_parabolic - transform_variables!(u_transformed, u, equations_parabolic) + @unpack u_transformed, gradients, flux_viscous = cache_parabolic + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, parabolic_scheme, cache, cache_parabolic) - calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, + calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, cache_parabolic) - calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad, + calc_viscous_fluxes!(flux_viscous, u_transformed, gradients, mesh, equations_parabolic, dg, cache, cache_parabolic) - calc_divergence!(du, u_transformed, t, viscous_flux, mesh, equations_parabolic, - boundary_conditions, dg, dg_parabolic, cache, cache_parabolic) + calc_divergence!(du, u_transformed, t, flux_viscous, mesh, equations_parabolic, + boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) return nothing diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl index 49bad9d8547..357e6550660 100644 --- a/src/solvers/dgmulti/sbp.jl +++ b/src/solvers/dgmulti/sbp.jl @@ -433,6 +433,7 @@ function estimate_dt(mesh::DGMultiMesh, dg::DGMultiPeriodicFDSBP) end # do nothing for interface terms if using a periodic operator +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::DGMultiMesh, equations, surface_integral, dg::DGMultiPeriodicFDSBP) @assert nelements(mesh, dg, cache) == 1 diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index e4103909646..68ff171b44f 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -56,6 +56,7 @@ end end end +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::P4estMesh{2}, equations, surface_integral, dg::DG) diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index 8723c5c70e0..66df88f0e7c 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -83,6 +83,7 @@ end return (i1, i2) end +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::P4estMesh{3}, equations, surface_integral, dg::DG) diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 2b947e15944..0ffaec2eced 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -60,6 +60,7 @@ include("dg_1d.jl") # 2D DG implementation include("dg_2d.jl") include("dg_2d_parallel.jl") +include("dg_2d_parabolic.jl") # 3D DG implementation include("dg_3d.jl") diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 16495000561..17d1c779364 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -371,6 +371,7 @@ end end +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack interfaces = cache diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 5215db6edc0..d725a8d828e 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -490,6 +490,7 @@ end end +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack interfaces = cache @@ -961,6 +962,7 @@ end end +# we pass in the hyperbolic `dg.surface_integral` as a dummy argument for dispatch function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DG, cache) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl new file mode 100644 index 00000000000..8099d5d6a6e --- /dev/null +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -0,0 +1,585 @@ +# 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{2}, 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) + + # 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{2}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations_parabolic, dg, i, j, element) + u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, equations_parabolic) + set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, 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{2}, equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + flux_viscous_x, flux_viscous_y = flux_viscous + + @threaded for element in eachelement(dg, cache) + # Calculate volume terms in one element + for j in eachnode(dg), i in eachnode(dg) + flux_1_node = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, element) + flux_2_node = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux_1_node, equations_parabolic, dg, ii, j, 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, 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{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache) + @unpack interfaces = cache_parabolic + @unpack orientations = interfaces + + flux_viscous_x, flux_viscous_y = 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 j in eachnode(dg), v in eachvariable(equations_parabolic) + interfaces.u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j, left_element] + interfaces.u[2, v, j, interface] = flux_viscous_x[v, 1, j, right_element] + end + else # if orientations[interface] == 2 + # interface in y-direction + for i in eachnode(dg), v in eachvariable(equations_parabolic) + interfaces.u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg), left_element] + interfaces.u[2, v, i, interface] = flux_viscous_y[v, i, 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{2}, 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 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, i, interface) + flux = 0.5 * (u_ll + u_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, left_direction, left_id] = flux[v] + surface_flux_values[v, i, 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{2}, 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 + + @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 l in eachnode(dg), v in eachvariable(equations_parabolic) + boundaries.u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l, element] + end + else # Element in +x direction of boundary + for l in eachnode(dg), v in eachvariable(equations_parabolic) + boundaries.u[2, v, l, boundary] = flux_viscous_x[v, 1, l, element] + end + end + else # if orientations[boundary] == 2 + # boundary in y-direction + if neighbor_sides[boundary] == 1 + # element in -y direction of boundary + for l in eachnode(dg), v in eachvariable(equations_parabolic) + boundaries.u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg), element] + end + else + # element in +y direction of boundary + for l in eachnode(dg), v in eachvariable(equations_parabolic) + boundaries.u[2, v, l, boundary] = flux_viscous_y[v, l, 1, element] + end + end + end + end + + return nothing +end + + +function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{2}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache, cache_parabolic) + gradients_x, gradients_y = gradients + flux_viscous_x, flux_viscous_y = flux_viscous # output arrays + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + # Get solution and gradients + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, element) + gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, element) + gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, element) + + # Calculate viscous flux and store each component for later use + flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node), 1, equations_parabolic) + flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node), 2, equations_parabolic) + set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, i, j, element) + set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, i, j, element) + end + end +end + + +# TODO: parabolic; decide if we should keep this, and if so, extend to 3D. +function get_unsigned_normal_vector_2d(direction) + if direction > 4 || direction < 1 + error("Direction = $direction; in 2D, direction should be 1, 2, 3, or 4.") + end + if direction == 1 || direction == 2 + return SVector(1.0, 0.0) + else + return SVector(0.0, 1.0) + end +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + return nothing +end + +function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_gradient_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_gradient_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, i, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + end + + # TODO: revisit if we want more general boundary treatments. + # This assumes the gradient numerical flux at the boundary is the gradient variable, + # which is consistent with BR1, LDG. + flux_inner = u_inner + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, boundary) + flux = boundary_condition(flux_inner, u_inner, get_unsigned_normal_vector_2d(direction), + x, t, Gradient(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::NamedTuple, + mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_divergence_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_divergence_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, + boundary_condition, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + + # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") + # of the viscous flux, as computed in `prolong2boundaries` of `calc_divergence!`` + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for i in eachnode(dg) + # Get viscous boundary fluxes + flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + flux_inner = flux_ll + else # Element is on the right, boundary on the left + flux_inner = flux_rr + end + + x = get_node_coords(node_coordinates, equations_parabolic, dg, i, boundary) + + # TODO: add a field in `cache.boundaries` for gradient information. + # Here, we pass in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information. + # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion2D and + # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesDiffusion2D as of 2022-6-27. + # It will not work with implementations which utilize `u_inner` to impose boundary conditions. + flux = boundary_condition(flux_inner, nothing, get_unsigned_normal_vector_2d(direction), + x, t, Divergence(), equations_parabolic) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + + +# Calculate the gradient of the transformed variables +function calc_gradient!(gradients, u_transformed, t, + mesh::TreeMesh{2}, equations_parabolic, + boundary_conditions_parabolic, dg::DG, cache, cache_parabolic) + + gradients_x, gradients_y = gradients + + # Reset du + @trixi_timeit timer() "reset gradients" begin + reset_du!(gradients_x, dg, cache) + reset_du!(gradients_y, 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 j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, element) + + for ii in eachnode(dg) + multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], u_node, equations_parabolic, dg, ii, j, 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, 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 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, + equations_parabolic, dg, i, interface) + flux = 0.5 * (u_ll + u_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, left_direction, left_id] = flux[v] + surface_flux_values[v, i, 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 l in eachnode(dg) + for v in eachvariable(equations_parabolic) + # surface at -x + gradients_x[v, 1, l, element] = ( + gradients_x[v, 1, l, element] - surface_flux_values[v, l, 1, element] * factor_1) + + # surface at +x + gradients_x[v, nnodes(dg), l, element] = ( + gradients_x[v, nnodes(dg), l, element] + surface_flux_values[v, l, 2, element] * factor_2) + + # surface at -y + gradients_y[v, l, 1, element] = ( + gradients_y[v, l, 1, element] - surface_flux_values[v, l, 3, element] * factor_1) + + # surface at +y + gradients_y[v, l, nnodes(dg), element] = ( + gradients_y[v, l, nnodes(dg), element] + surface_flux_values[v, l, 4, 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) + 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{2}, 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_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{2}, + equations::AbstractEquationsParabolic, dg::DG, cache) + + @threaded for element in eachelement(dg, cache) + factor = cache.elements.inverse_jacobian[element] + + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + du[v, i, j, element] *= factor + end + end + end + + return nothing +end + +end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 611691f9f85..acd1b31d646 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -544,6 +544,7 @@ end end +# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack interfaces = cache diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 4128e1c6295..f2872652c8b 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -81,6 +81,7 @@ end # prolong the solution into the convenience array in the interior interface container +# We pass the `surface_integral` argument solely for dispatch # Note! this routine is for quadrilateral elements with "right-handed" orientation function prolong2interfaces!(cache, u, mesh::UnstructuredMesh2D, diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 4a792af4840..c61e5cf918c 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -45,24 +45,24 @@ isdir(outdir) && rm(outdir, recursive=true) @test boundary_condition_do_nothing(u0, nothing) == u0 @unpack cache, cache_parabolic, equations_parabolic = semi - @unpack u_grad = cache_parabolic - for dim in eachindex(u_grad) - fill!(u_grad[dim], zero(eltype(u_grad[dim]))) + @unpack gradients = cache_parabolic + for dim in eachindex(gradients) + fill!(gradients[dim], zero(eltype(gradients[dim]))) end t = 0.0 # pass in `boundary_condition_periodic` to skip boundary flux/integral evaluation - Trixi.calc_gradient!(u_grad, ode.u0, t, mesh, equations_parabolic, + Trixi.calc_gradient!(gradients, ode.u0, t, mesh, equations_parabolic, boundary_condition_periodic, dg, cache, cache_parabolic) @unpack x, y = mesh.md - @test getindex.(u_grad[1], 1) ≈ 2 * x .* y - @test getindex.(u_grad[2], 1) ≈ x.^2 + @test getindex.(gradients[1], 1) ≈ 2 * x .* y + @test getindex.(gradients[2], 1) ≈ x.^2 - u_flux = similar.(u_grad) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, + u_flux = similar.(gradients) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, gradients, mesh, equations_parabolic, dg, cache, cache_parabolic) - @test u_flux[1] ≈ u_grad[1] - @test u_flux[2] ≈ u_grad[2] + @test u_flux[1] ≈ gradients[1] + @test u_flux[2] ≈ gradients[2] du = similar(ode.u0) Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, boundary_condition_periodic, @@ -70,19 +70,84 @@ isdir(outdir) && rm(outdir, recursive=true) @test getindex.(du, 1) ≈ 2 * y end - @trixi_testset "elixir_advection_diffusion_periodic.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion_periodic.jl"), + @trixi_testset "DGMulti: elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.2485803335154642], + linf = [1.079606969242132] + ) + end + + @trixi_testset "DGMulti: elixir_advection_diffusion_periodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_periodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.03180371984888462], linf = [0.2136821621370909] ) end - @trixi_testset "elixir_advection_diffusion.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion.jl"), + @trixi_testset "DGMulti: elixir_advection_diffusion_nonperiodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_nonperiodic.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), - l2 = [0.2485803335154642], - linf = [1.079606969242132] + l2 = [0.002123168335604323], + linf = [0.00963640423513712] + ) + end + + @trixi_testset "DGMulti: elixir_navier_stokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_convergence.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.00153550768125133, 0.0033843168272696357, 0.0036531858107444067, 0.009948436427519428], + linf = [0.005522560467190019, 0.013425258500731063, 0.013962115643483375, 0.027483102120516634] + ) + end + + @trixi_testset "DGMulti: elixir_navier_stokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_lid_driven_cavity.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.5), + l2 = [0.0002215612522711349, 0.028318325921400257, 0.009509168701069035, 0.028267900513539248], + linf = [0.0015622789413053395, 0.14886653390741342, 0.07163235655334241, 0.19472785105216417] + ) + end + + @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.4), polydeg=5, + l2 = [4.0915532997994255e-6], + linf = [2.3040850347877395e-5] + ) + end + + @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.007646800618485118], + linf = [0.10067621050468958] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.0021116725306635146, 0.0034322351490824465, 0.003874252819611102, 0.012469246082522416], + linf = [0.012006418939297214, 0.03552087120958058, 0.02451274749176294, 0.11191122588577151] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (flux differencing)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + volume_integral=VolumeIntegralFluxDifferencing(flux_central), + l2 = [0.0021116725306635146, 0.0034322351490824465, 0.003874252819611102, 0.012469246082522416], + linf = [0.012006418939297214, 0.03552087120958058, 0.02451274749176294, 0.11191122588577151] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navier_stokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_lid_driven_cavity.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.0001514457152968994, 0.018766076072331786, 0.007065070765651992, 0.020839900573430787], + linf = [0.0014523369373645734, 0.12366779944955876, 0.055324509971157544, 0.1609992780534526] ) end