From 370c560361ca9ba08a44a4ae0ed7df9e756029d9 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 11 Oct 2022 12:02:56 -0500 Subject: [PATCH] Incorporate parabolic terms into `main` (#1149) * Parabolic terms for DGMulti (#1136) * minor formatting * adding multi-term semidiscretization * fixing bug and adding test * changing name * changing name, importing SplitODEFunction * fixing test by adding norm * trying another test fix * adding norm back * adjusting tol * try to fix test again * syntax error * removing norm frm test * adding cache to parabolic terms (enable reuse of inviscid cache) * generalize dgmulti iterators * fix typo * generalizing reset_du! for dgmulti * WIP gradient and divergence computations * fix name * more generalization * divergence + remove repeated code * adding cache * remove test * update diffusion test * update diffusion rhs * moving routines around * fix multiple includes * removing time dependence for now from viscous rhs * export ScalarDiffusion2D * fix test * working periodic advection RHS * cleanup and comments * making periodic advection-diffusion example * trim whitespace * adding t back in to parabolic rhs add t back in * refactoring boundary condition code * update elixir with BCs * change default DGMulti interface flux behavior same as in https://github.com/trixi-framework/Trixi.jl/pull/1139 * cleanup * renaming variables * BC prepping * cleanup of elixir * more name fixing * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Hendrik Ranocha * adding Neumann BCs * switching to weak form weak form more weak form weak form * fixing cache variable names * simplified parabolic boundary flux treatement * updated elixir with different types of BCs * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper * rhs! -> rhs_parabolic! * move Gradient/Divergence types * add "no boundary condition" BC * rename ScalarDiffusion -> LaplaceDiffusion r r * rhs -> rhs_parabolic! again * comments * moving tests around * scalarDiffusion -> LaplaceDiffusion * tuple-based constructor * updating examples * add ScalarPlotData2D docstring (#1145) Co-authored-by: Jesse Chan * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha * renaming for consistency parabolic_equations -> equations_parabolic no_boundary_condition -> boundary_condition_do_nothing * move Neumann BCs into equations.jl * renaming, e.g., ParabolicEquations -> EquationsParabolic, etc * rename test module * generalize LaplaceDiffusion2D to any number of variables * removing saveat times and standardizing `tol` names * parabolic_cache -> cache_parabolic cache * renaming `create_cache` -> `create_cache_parabolic` avoids conflicts when using dg::DGMultiFluxDiff * Update examples/dgmulti_2d/elixir_advection_diffusion.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper * parabolic_boundary_conditions -> boundary_conditions_parabolic similarly for hyperbolic d * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper * add equations as field to LaplaceDiffusion2D * moving Neumann BCs into equations.jl and adding default BC behavior * running parabolic tests, adding integration test * adding parabolic tests to ci * missing comma * Update src/equations/laplace_diffusion_2d.jl Co-authored-by: Hendrik Ranocha * fix test * fix test typo * fix tests (for real this time) * improve test coverage Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha Co-authored-by: Michael Schlottke-Lakemper * Additional work on parabolic terms (#1148) * make "do nothing BC" a struct, move it to basic_types * add varnames for LaplaceDiffusion * add forgotten @threaded * update comment * Revert "add forgotten @threaded" This reverts commit 1564501a64ecba856ef9cc28bda1a34278ff1eb6. * add @threaded (without enumerate) * remove periodic advection diffusion elixir (redundant) * adding test * one more test * Update src/basic_types.jl Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha * Reorganization of boundary conditions (#1152) * moving "default" boundary conditions back into equations file * add periodic advection_diffusion back in * updating test * update test * Apply suggestions from code review Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha * Adding parabolic solver option (#1153) * adding parabolic solver types * incorporating parabolic solver types into dg_parabolic penalties * specializing parabolic solver options in LaplaceDiffusion2D * moving things around to fix CI * dropped arg in test * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * renaming files and ViscousFlux solvers renaming files rename rename * Added TODO note for penalty term * ViscousFlux___ -> ViscousFormulation___ * move function def around Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha * rearrange include order * 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 * Add commentary on reuse of interfaces/boundaries data structures (#1199) * Add comments on storing flux values in `u` * Rearrange function names for consistency * 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: Andrew Winters Co-authored-by: Hendrik Ranocha Co-authored-by: Andrew Winters * Adding parabolic term description to NEWS.md (#1207) * Add Literate docs for advection-diffusion (#1208) * adding literate docs for advection-diffusion * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * rename "adding parabolic terms" -> "parabolic terms" Co-authored-by: Hendrik Ranocha * Adding a tutorial on adding new parabolic terms (#1209) * Viscous terms from entropy gradients (#1202) * first draft of entropy stable gradients * update TODOs * added elixir with discontinuous solution for ES testing * forgot one TODO marker * adding GradientVariables**** type parameter * fix bug (missing type param in constructor) * add elixir with periodic CNS manufactured solution * fix gradient conversion * fix entropy variables w/BCs * updating elixirs with `gradient_variables=GradientVariablesPrimitive()` * Update src/equations/compressible_navier_stokes_2d.jl Co-authored-by: Michael Schlottke-Lakemper * removing unused routines * removing specialized entropy2cons and cons2entropy routines * adding more tests * removing unused elixirs * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * adding docs/warning for gradient variable type * extracting w_inner[4] into a more clearly named variable Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan Co-authored-by: Jesse Chan Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> Co-authored-by: Hendrik Ranocha * adding some minor documentation (#1214) * fix parabolic PID (#1224) * fix parabolic PID * activate analysis_callback in Navier-Stokes convergence test * fix typo * Rescale compressible Navier-Stokes (#1230) * remove the nondimensionalization and update the docstring for CNS * update TreeMesh elixirs and test values * update the elixirs and test values for DGMulti * fix unbalanced parentheses * file renaming for consistency * update file names in the parabolic test script * remove the gas constant as a parameter and update comments * rename viscosity to dynamic_viscosity * add comment about units in the docstring * Update src/equations/compressible_navier_stokes_2d.jl Co-authored-by: Michael Schlottke-Lakemper * unify notation to call the viscosity mu * fix undefined variable Co-authored-by: Michael Schlottke-Lakemper * improve type stability of semidiscretize for parabolic semi (#1231) * Apply suggestions from code review Co-authored-by: Hendrik Ranocha Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha Co-authored-by: Andrew Winters --- .github/workflows/ci.yml | 1 + NEWS.md | 16 + .../src/files/adding_new_parabolic_terms.jl | 160 +++++ docs/literate/src/files/parabolic_terms.jl | 88 +++ docs/make.jl | 2 + .../dgmulti_2d/elixir_advection_diffusion.jl | 56 ++ .../elixir_advection_diffusion_nonperiodic.jl | 72 +++ .../elixir_advection_diffusion_periodic.jl | 35 + .../elixir_navierstokes_convergence.jl | 205 ++++++ .../elixir_navierstokes_lid_driven_cavity.jl | 73 +++ .../elixir_advection_diffusion.jl | 83 +++ .../elixir_advection_diffusion_nonperiodic.jl | 89 +++ .../elixir_navierstokes_convergence.jl | 213 ++++++ .../elixir_navierstokes_lid_driven_cavity.jl | 79 +++ src/Trixi.jl | 21 +- src/auxiliary/auxiliary.jl | 32 +- src/basic_types.jl | 14 + src/equations/compressible_euler_1d.jl | 4 +- src/equations/compressible_euler_2d.jl | 7 +- src/equations/compressible_euler_3d.jl | 8 +- .../compressible_euler_multicomponent_1d.jl | 10 +- .../compressible_euler_multicomponent_2d.jl | 12 +- .../compressible_navier_stokes_2d.jl | 407 ++++++++++++ src/equations/equations.jl | 20 + src/equations/equations_parabolic.jl | 11 + src/equations/laplace_diffusion_2d.jl | 58 ++ ...semidiscretization_hyperbolic_parabolic.jl | 264 ++++++++ src/solvers/dgmulti.jl | 3 + src/solvers/dgmulti/dg.jl | 42 +- src/solvers/dgmulti/dg_parabolic.jl | 322 ++++++++++ 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 | 605 ++++++++++++++++++ src/solvers/dgsem_tree/dg_3d.jl | 1 + src/solvers/dgsem_unstructured/dg_2d.jl | 1 + src/solvers/solvers.jl | 3 +- src/solvers/solvers_parabolic.jl | 31 + test/runtests.jl | 4 + test/test_dgmulti_2d.jl | 1 + test/test_parabolic_2d.jl | 185 ++++++ 44 files changed, 3196 insertions(+), 49 deletions(-) create mode 100644 docs/literate/src/files/adding_new_parabolic_terms.jl create mode 100644 docs/literate/src/files/parabolic_terms.jl create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion.jl create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl create mode 100644 examples/dgmulti_2d/elixir_navierstokes_convergence.jl create mode 100644 examples/dgmulti_2d/elixir_navierstokes_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_navierstokes_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl create mode 100644 src/equations/compressible_navier_stokes_2d.jl create mode 100644 src/equations/equations_parabolic.jl create mode 100644 src/equations/laplace_diffusion_2d.jl create mode 100644 src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl create mode 100644 src/solvers/dgmulti/dg_parabolic.jl create mode 100644 src/solvers/dgsem_tree/dg_2d_parabolic.jl create mode 100644 src/solvers/solvers_parabolic.jl create mode 100644 test/test_parabolic_2d.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f257da07a44..6aba503691e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -69,6 +69,7 @@ jobs: - p4est_part1 - p4est_part2 - unstructured_dgmulti + - parabolic - paper_self_gravitating_gas_dynamics - misc_part1 - misc_part2 diff --git a/NEWS.md b/NEWS.md index 134173068b9..c816bdfe561 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,22 @@ for human readability. #### Added +- Experimental support for 2D parabolic diffusion terms has been added. + * `LaplaceDiffusion2D` and `CompressibleNavierStokesDiffusion2D` can be used to add + diffusion to systems. `LaplaceDiffusion2D` can be used to add scalar diffusion to each + equation of a system, while `CompressibleNavierStokesDiffusion2D` can be used to add + Navier-Stokes diffusion to `CompressibleEulerEquations2D`. + * Parabolic boundary conditions can be imposed as well. For `LaplaceDiffusion2D`, both + `Dirichlet` and `Neumann` conditions are supported. For `CompressibleNavierStokesDiffusion2D`, + viscous no-slip velocity boundary conditions are supported, along with adiabatic and isothermal + temperature boundary conditions. See the boundary condition container + `BoundaryConditionNavierStokesWall` and boundary condition types `NoSlip`, `Adiabatic`, and + `Isothermal` for more information. + * `CompressibleNavierStokesDiffusion2D` can utilize both primitive variables (which are not + guaranteed to provably dissipate entropy) and entropy variables (which provably dissipate + entropy at the semi-discrete level). + * Please check the `examples` directory for further information about the supported setups. + Further documentation will be added later. - Numerical fluxes `flux_shima_etal_turbo` and `flux_ranocha_turbo` that are equivalent to their non-`_turbo` counterparts but may enable specialized methods making use of SIMD instructions to increase runtime efficiency diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl new file mode 100644 index 00000000000..a6cecb351be --- /dev/null +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -0,0 +1,160 @@ +#src # Adding new parabolic terms. + +# This demo illustrates the steps involved in adding new parabolic terms for the scalar +# advection equation. In particular, we will add an anisotropic diffusion. We begin by +# defining the hyperbolic (advection) part of the advection-diffusion equation. + +using OrdinaryDiffEq +using Trixi + + +advection_velocity = (1.0, 1.0) +equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity); + +# ## Define a new parabolic equation type +# +# Next, we define a 2D parabolic diffusion term type. This is similar to [`LaplaceDiffusion2D`](@ref) +# except that the `diffusivity` field refers to a spatially constant diffusivity matrix now. Note that +# `ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have +# information about the hyperbolic system available to the parabolic part so that we can reuse +# functions defined for hyperbolic equations (such as `varnames`). + +struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1} + diffusivity::T + equations_hyperbolic::E +end + +varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) = + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) + +# Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system +# is of the form +# ```math +# \partial_t u(t,x) + \partial_x (f_1(u) - g_1(u, \nabla u)) +# + \partial_y (f_2(u) - g_2(u, \nabla u)) = 0 +# ``` +# where ``f_1(u)``, ``f_2(u)`` are the hyperbolic fluxes and ``g_1(u, \nabla u)``, ``g_2(u, \nabla u)`` denote +# the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components +# of the matrix-vector product involving `diffusivity` and the gradient vector. +# +# Here, we specialize the flux to our new parabolic equation type `ConstantAnisotropicDiffusion2D`. + +function Trixi.flux(u, gradients, orientation::Integer, equations_parabolic::ConstantAnisotropicDiffusion2D) + @unpack diffusivity = equations_parabolic + dudx, dudy = gradients + if orientation == 1 + return SVector(diffusivity[1, 1] * dudx + diffusivity[1, 2] * dudy) + else # if orientation == 2 + return SVector(diffusivity[2, 1] * dudx + diffusivity[2, 2] * dudy) + end +end + +# ## Defining boundary conditions + +# Trixi.jl's implementation of parabolic terms discretizes both the gradient and divergence +# using weak formulation. In other words, we discretize the system +# ```math +# \begin{aligned} +# \bm{q} &= \nabla u \\ +# \bm{\sigma} &= \begin{pmatrix} g_1(u, \bm{q}) \\ g_2(u, \bm{q}) \end{pmatrix} \\ +# \text{viscous contribution } &= \nabla \cdot \bm{\sigma} +# \end{aligned} +# ``` +# +# Boundary data must be specified for all spatial derivatives, e.g., for both the gradient +# equation ``\bm{q} = \nabla u`` and the divergence of the viscous flux +# ``\nabla \cdot \bm{\sigma}``. We account for this by introducing internal `Gradient` +# and `Divergence` types which are used to dispatch on each type of boundary condition. +# +# As an example, let us introduce a Dirichlet boundary condition with constant boundary data. + +struct BoundaryConditionConstantDirichlet{T <: Real} + boundary_value::T +end + +# This boundary condition contains only the field `boundary_value`, which we assume to be some +# real-valued constant which we will impose as the Dirichlet data on the boundary. +# +# Boundary conditions have generally been defined as "callable structs" (also known as "functors"). +# For each boundary condition, we need to specify the appropriate boundary data to return for both +# the `Gradient` and `Divergence`. Since the gradient is operating on the solution `u`, the boundary +# data should be the value of `u`, and we can directly impose Dirichlet data. + +@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Gradient, + equations_parabolic::ConstantAnisotropicDiffusion2D) + return boundary_condition.boundary_value +end + +# While the gradient acts on the solution `u`, the divergence acts on the viscous flux ``\bm{\sigma}``. +# Thus, we have to supply boundary data for the `Divergence` operator that corresponds to ``\bm{\sigma}``. +# However, we've already imposed boundary data on `u` for a Dirichlet boundary condition, and imposing +# boundary data for ``\bm{\sigma}`` might overconstrain our problem. +# +# Thus, for the `Divergence` boundary data under a Dirichlet boundary condition, we simply return +# `flux_inner`, which is boundary data for ``\bm{\sigma}`` computed using the "inner" or interior solution. +# This way, we supply boundary data for the divergence operation without imposing any additional conditions. + +@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Divergence, + equations_parabolic::ConstantAnisotropicDiffusion2D) + return flux_inner +end + +# ### A note on the choice of gradient variables +# +# It is often simpler to transform the solution variables (and solution gradients) to another set of +# variables prior to computing the viscous fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref) +# for an example of this). If this is done, then the boundary condition for the `Gradient` operator +# should be modified accordingly as well. +# +# ## Putting things together +# +# Finally, we can instantiate our new parabolic equation type, define boundary conditions, +# and run a simulation. The specific anisotropic diffusion matrix we use produces more +# dissipation in the direction ``(1, -1)`` as an isotropic diffusion. +# +# For boundary conditions, we impose that ``u=1`` on the left wall, ``u=2`` on the bottom +# wall, and ``u = 0`` on the outflow walls. The initial condition is taken to be ``u = 0``. +# Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary +# conditions, since we have not defined its behavior for the hyperbolic part. + +using Trixi: SMatrix +diffusivity = 5.0e-2 * SMatrix{2, 2}([2 -1; -1 2]) +equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); + +boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)), + y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(2.0)), + y_pos = boundary_condition_do_nothing, + x_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(1.0), + y_neg = BoundaryConditionConstantDirichlet(2.0), + y_pos = BoundaryConditionConstantDirichlet(0.0), + x_pos = BoundaryConditionConstantDirichlet(0.0)); + +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)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure + +initial_condition = (x, t, equations) -> SVector(0.0) + +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations_hyperbolic, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions_hyperbolic, + boundary_conditions_parabolic)) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) +callbacks = CallbackSet(SummaryCallback()) +time_int_tol = 1.0e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks); + +using Plots +plot(sol) + diff --git a/docs/literate/src/files/parabolic_terms.jl b/docs/literate/src/files/parabolic_terms.jl new file mode 100644 index 00000000000..aeceb7b7e6f --- /dev/null +++ b/docs/literate/src/files/parabolic_terms.jl @@ -0,0 +1,88 @@ +#src # Parabolic terms (advection-diffusion). + +# Experimental support for parabolic diffusion terms is available in Trixi.jl. +# This demo illustrates parabolic terms for the advection-diffusion equation. + +using OrdinaryDiffEq +using Trixi + +# ## Splitting a system into hyperbolic and parabolic parts. + +# For a mixed hyperbolic-parabolic system, we represent the hyperbolic and parabolic +# parts of the system separately. We first define the hyperbolic (advection) part of +# the advection-diffusion equation. + +advection_velocity = (1.5, 1.0) +equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity); + +# Next, we define the parabolic diffusion term. The constructor requires knowledge of +# `equations_hyperbolic` to be passed in because the [`LaplaceDiffusion2D`](@ref) applies +# diffusion to every variable of the hyperbolic system. + +diffusivity = 5.0e-2 +equations_parabolic = LaplaceDiffusion2D(diffusivity, equations_hyperbolic); + +# ## Boundary conditions + +# As with the equations, we define boundary conditions separately for the hyperbolic and +# parabolic part of the system. For this example, we impose inflow BCs for the hyperbolic +# system (no condition is imposed on the outflow), and we impose Dirichlet boundary conditions +# for the parabolic equations. Both `BoundaryConditionDirichlet` and `BoundaryConditionNeumann` +# are defined for `LaplaceDiffusion2D`. +# +# The hyperbolic and parabolic boundary conditions are assumed to be consistent with each other. + +boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) + +boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), + y_neg = boundary_condition_zero_dirichlet, + y_pos = boundary_condition_do_nothing, + x_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), + y_neg = boundary_condition_zero_dirichlet, + y_pos = boundary_condition_zero_dirichlet, + x_pos = boundary_condition_zero_dirichlet); + +# ## Defining the solver and mesh + +# The process of creating the DG solver and mesh is the same as for a purely +# hyperbolic system of equations. + +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)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure + +initial_condition = (x, t, equations) -> SVector(0.0); + +# ## Semidiscretizing and solving + +# To semidiscretize a hyperbolic-parabolic system, we create a [`SemidiscretizationHyperbolicParabolic`](@ref). +# This differs from a [`SemidiscretizationHyperbolic`](@ref) in that we pass in a `Tuple` containing both the +# hyperbolic and parabolic equation, as well as a `Tuple` containing the hyperbolic and parabolic +# boundary conditions. + +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations_hyperbolic, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions_hyperbolic, + boundary_conditions_parabolic)) + +# The rest of the code is identical to the hyperbolic case. We create a system of ODEs through +# `semidiscretize`, defining callbacks, and then passing the system to OrdinaryDiffEq.jl. + +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan) +callbacks = CallbackSet(SummaryCallback()) +time_int_tol = 1.0e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks); + +# We can now visualize the solution, which develops a boundary layer at the outflow boundaries. + +using Plots +plot(sol) + diff --git a/docs/make.jl b/docs/make.jl index 57321677480..d8d1298fba6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -39,6 +39,8 @@ files = [ # Topic: equations "Adding a new scalar conservation law" => "adding_new_scalar_equations.jl", "Adding a non-conservative equation" => "adding_nonconservative_equation.jl", + "Parabolic terms" => "parabolic_terms.jl", + "Adding new parabolic terms" => "adding_new_parabolic_terms.jl", # Topic: meshes "Adaptive mesh refinement" => "adaptive_mesh_refinement.jl", "Structured mesh with curvilinear mapping" => "structured_mesh_mapping.jl", diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl new file mode 100644 index 00000000000..1273e25859e --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -0,0 +1,56 @@ +using Trixi, OrdinaryDiffEq + +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +equations = LinearScalarAdvectionEquation2D(1.5, 1.0) +equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations) + +initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) +initial_condition = initial_condition_zero + +# tag different boundary segments +left(x, tol=50*eps()) = abs(x[1] + 1) < tol +right(x, tol=50*eps()) = abs(x[1] - 1) < tol +bottom(x, tol=50*eps()) = abs(x[2] + 1) < tol +top(x, tol=50*eps()) = abs(x[2] - 1) < tol +is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) + +# BC types +boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.1 * x[2])) +boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) +boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0)) + +# define inviscid boundary conditions +boundary_conditions = (; :left => boundary_condition_left, + :bottom => boundary_condition_zero, + :top => boundary_condition_do_nothing, + :right => boundary_condition_do_nothing) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :left => boundary_condition_left, + :bottom => boundary_condition_zero, + :top => boundary_condition_zero, + :right => boundary_condition_neumann_zero) + +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-6 +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_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..bb22f494374 --- /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] + 0.5) < tol +top(x, tol=50*eps()) = abs(x[2] - 0.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_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl new file mode 100644 index 00000000000..8058e156969 --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -0,0 +1,35 @@ +using Trixi, OrdinaryDiffEq + +dg = DGMulti(polydeg = 1, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +equations = LinearScalarAdvectionEquation2D(0.0, 0.0) +equations_parabolic = LaplaceDiffusion2D(5.0e-1, equations) + +function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) + return SVector(exp(-100 * (x[1]^2 + x[2]^2))) +end +initial_condition = initial_condition_sharp_gaussian + +mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, dg) + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + dt = time_int_tol, save_everystep=false, callback=callbacks) + +summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..f7039627238 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -0,0 +1,205 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +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, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = 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() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_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_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 * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) + # 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 * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) + # 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 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ + - 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 ) * mu_ ) + + 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_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..0837cc2fd55 --- /dev/null +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -0,0 +1,73 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 0.001 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), + Prandtl=prandtl_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 = 0.1 + rho = 1.0 + u, v = 0.0, 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..2ef0ffe766f --- /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_eriksson_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_eriksson_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_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..1145c22fcaf --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,213 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = (-1.0, -1.0) # 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() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_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_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 * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) + # 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 * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) + # 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 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ + - 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 ) * mu_ ) + + 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, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary + diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..d94187eade1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -0,0 +1,79 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 0.001 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), + Prandtl=prandtl_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 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = 0.1 + rho = 1.0 + u, v = 0.0, 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 7142436edd3..5811a950505 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -26,7 +26,8 @@ using SparseArrays: AbstractSparseMatrix, AbstractSparseMatrixCSC, sparse, dropt using Reexport: @reexport using SciMLBase: CallbackSet, DiscreteCallback, - ODEProblem, ODESolution, ODEFunction + ODEProblem, ODESolution, ODEFunction, + SplitODEProblem import SciMLBase: get_du, get_tmp_cache, u_modified!, AbstractODEIntegrator, init, step!, check_error, get_proposed_dt, set_proposed_dt!, @@ -104,8 +105,10 @@ include("auxiliary/p4est.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") +include("equations/equations_parabolic.jl") # these depend on parabolic solver types include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") +include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") @@ -131,6 +134,11 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D +export LaplaceDiffusion2D, + CompressibleNavierStokesDiffusion2D + +export GradientVariablesPrimitive, GradientVariablesEntropy + 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, @@ -151,11 +159,14 @@ export initial_condition_constant, initial_condition_density_wave, initial_condition_weak_blast_wave -export boundary_condition_periodic, +export boundary_condition_do_nothing, + boundary_condition_periodic, BoundaryConditionDirichlet, + 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 @@ -185,6 +196,8 @@ export nelements, nnodes, nvariables, export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate +export SemidiscretizationHyperbolicParabolic + export SemidiscretizationEulerAcoustics export SemidiscretizationEulerGravity, ParametersEulerGravity, @@ -215,6 +228,8 @@ export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure export DGMulti, estimate_dt, DGMultiMesh, GaussSBP export VertexMappedMesh # TODO: DGMulti, v0.5. Remove deprecated VertexMappedMesh in next release +export ViscousFormulationBassiRebay1, ViscousFormulationLocalDG + # Visualization-related exports export PlotData1D, PlotData2D, ScalarPlotData2D, getmesh, adapt_to_mesh_level!, adapt_to_mesh_level diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 2152727940c..bc23e831a45 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -24,9 +24,9 @@ timer() = main_timer """ - PerformanceCounter + PerformanceCounter() -A `PerformanceCounter` be used to track the runtime performance of some calls. +A `PerformanceCounter` can be used to track the runtime performance of some calls. Add a new runtime measurement via `put!(counter, runtime)` and get the averaged runtime of all measurements added so far via `take!(counter)`, resetting the `counter`. @@ -51,6 +51,34 @@ function Base.put!(counter::PerformanceCounter, runtime::Real) end +""" + PerformanceCounterList{N}() + +A `PerformanceCounterList{N}` can be used to track the runtime performance of +calls to multiple functions, adding them up. +Add a new runtime measurement via `put!(counter.counters[i], runtime)` and get +the averaged runtime of all measurements added so far via `take!(counter)`, +resetting the `counter`. +""" +struct PerformanceCounterList{N} + counters::NTuple{N, PerformanceCounter} +end + +function PerformanceCounterList{N}() where {N} + counters = ntuple(_ -> PerformanceCounter(), Val{N}()) + return PerformanceCounterList{N}(counters) +end + +function Base.take!(counter::PerformanceCounterList) + time_per_call = 0.0 + for c in counter.counters + time_per_call += take!(c) + end + return time_per_call +end + + + """ examples_dir() diff --git a/src/basic_types.jl b/src/basic_types.jl index d6ae7f9fceb..1116f3e6b66 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -70,5 +70,19 @@ const boundary_condition_periodic = BoundaryConditionPeriodic() Base.show(io::IO, ::BoundaryConditionPeriodic) = print(io, "boundary_condition_periodic") +""" + boundary_condition_do_nothing = BoundaryConditionDoNothing() + +Imposing no boundary condition just evaluates the flux at the inner state. +""" +struct BoundaryConditionDoNothing end + +@inline function (boundary_condition::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) + return inner_flux_or_state +end + +const boundary_condition_do_nothing = BoundaryConditionDoNothing() + +Base.show(io::IO, ::BoundaryConditionDoNothing) = print(io, "boundary_condition_do_nothing") end # @muladd 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..7bda9999605 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} @@ -36,7 +36,6 @@ Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` the pressure. - """ struct CompressibleEulerEquations2D{RealT<:Real} <: AbstractCompressibleEulerEquations{2, 4} gamma::RealT # ratio of specific heats 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..9f7758aa476 --- /dev/null +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -0,0 +1,407 @@ +@doc raw""" + CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, equations, + gradient_variables=GradientVariablesPrimitive()) + +These equations contain the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective terms from +the [`CompressibleEulerEquations2D`](@ref). + +- `gamma`: adiabatic constant, +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) +- `gradient_variables`: which variables the gradients are taken with respect to. + Defaults to `GradientVariablesPrimitive()`. + +Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., +[``\mu``] = kg m⁻¹ s⁻¹. + +The particular form of the compressible Navier-Stokes implemented is +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho \mathbf{v} \\ \rho e +\end{pmatrix} ++ +\nabla \cdot +\begin{pmatrix} + \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e + p) \mathbf{v} +\end{pmatrix} += +\nabla \cdot +\begin{pmatrix} +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +\end{pmatrix} +``` +where the system is closed with the ideal gas assumption giving +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \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}}. +``` +From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see +that the gas constant `R` cancels and the heat flux becomes +```math +\nabla q = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) +``` +which is the form implemented below in the [`flux`](@ref) function. + +In two spatial dimensions we require gradients for three quantities, e.g., +primitive quantities +```math +\nabla v_1,\, \nabla v_2,\, \nabla 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} +``` + +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. +""" +struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4} + # TODO: parabolic + # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations + # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + mu::RealT # viscosity + Pr::RealT # Prandtl number + kappa::RealT # thermal diffusivity for Fick's law + + equations_hyperbolic::E # CompressibleEulerEquations2D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +""" +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. + +`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters +for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be +`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable +formulation from +- Hughes, Mallet, Franca (1986) + A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the + compressible Euler and Navier-Stokes equations and the second law of thermodynamics. + [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + +Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. +""" +struct GradientVariablesPrimitive end +struct GradientVariablesEntropy end + +# default to primitive gradient variables +function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; + mu, Prandtl, + gradient_variables = GradientVariablesPrimitive()) + gamma = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + μ, Pr = promote(mu, Prandtl) + + # Under the assumption of constant Prandtl number the thermal conductivity + # constant is kappa = gamma μ / ((gamma-1) Pr). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * inv_gamma_minus_one / Pr + + CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, + μ, Pr, kappa, + equations, gradient_variables) +end + +# TODO: parabolic +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::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{GradientVariablesPrimitive}) = cons2prim +gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) = cons2entropy + + +# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramíreza, Hennemann, Hindenlang, Winters, and Gassner +# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive +# MHD Equations. Part II: Subcell Finite Volume Shock Capturing" +# where one sets the magnetic field components equal to 0. +function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion2D) + # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. + rho, v1, v2, _ = convert_transformed_to_primitive(u, equations) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T) + # either computed directly or reverse engineered from the gradient of the entropy vairables + # by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, dTdx = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, dTdy = convert_derivative_to_primitive(u, gradients[2], equations) + + # 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) = -kappa * grad(p / (R rho)) + # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) + # Note, the gas constant cancels under this formulation, so it is not present + # in the implementation + q1 = equations.kappa * dTdx + q2 = equations.kappa * dTdy + + # Constant dynamic viscosity is copied to a variable for readibility. + # Offers flexibility for dynamic viscosity via Sutherland's law where it depends + # on temperature and reference values, Ts and Tref such that mu(T) + mu = equations.mu + + if orientation == 1 + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = ( 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 + +# Convert conservative variables to entropy +# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms +# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion2D`, +# but this may be confusing to new users. +cons2entropy(u, equations::CompressibleNavierStokesDiffusion2D) = cons2entropy(u, equations.equations_hyperbolic) +entropy2cons(w, equations::CompressibleNavierStokesDiffusion2D) = entropy2cons(w, equations.equations_hyperbolic) + +# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. +# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# variables into primitive variables. +@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) + return u_transformed +end + +# TODO: parabolic. Make this more efficient! +@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + # note: this uses CompressibleNavierStokesDiffusion2D versions of cons2prim and entropy2cons + return cons2prim(entropy2cons(u_transformed, equations), equations) +end + + +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, T). +# Helpful because then the diffusive fluxes have the same form as on paper. +# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. +# TODO: parabolic; entropy stable viscous terms +@inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) + return gradient +end + +# the first argument is always the "transformed" variables. +@inline function convert_derivative_to_primitive(w, gradient_entropy_vars, + equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + + # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back. + # We can fix this if we directly compute v1, v2, T from the entropy variables + u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion2D + rho, rho_v1, rho_v2, _ = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = temperature(u, equations) + + return SVector(gradient_entropy_vars[1], + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[4]), # grad(u) = T*(grad(w_2)+v1*grad(w_4)) + T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[4]), # grad(v) = T*(grad(w_3)+v2*grad(w_4)) + T * T * gradient_entropy_vars[4] # grad(T) = T^2*grad(w_4)) + ) +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 / 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{GradientVariablesPrimitive}) + 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{GradientVariablesPrimitive}) + # 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{GradientVariablesPrimitive}) + 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{GradientVariablesPrimitive}) + return flux_inner +end + +# specialized BC impositions for GradientVariablesEntropy. + +# This should return a SVector containing the boundary values of entropy variables. +# Here, `w_inner` are the transformed variables (e.g., entropy variables). +# +# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions +# for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022. +# DOI: 10.1016/j.jcp.2021.110723 +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) + negative_rho_inv_p = w_inner[4] # w_4 = -rho / p + return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p, negative_rho_inv_p) +end + +# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness. +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + 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, w_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + 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) + + # the entropy variables w2 = rho * v1 / p = v1 / T = -v1 * w4. Similarly for w3 + w4 = -1 / T + return SVector(w_inner[1], -v1 * w4, -v2 * w4, w4) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) + return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4]) +end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index dc3bfffba30..f5fbd8cb411 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -175,6 +175,25 @@ end return flux end +# operator types used for dispatch on parabolic boundary fluxes +struct Gradient end +struct Divergence end + +""" + BoundaryConditionNeumann(boundary_normal_flux_function) + +Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic +equations that uses the function `boundary_normal_flux_function` to specify the values of the normal +flux at the boundary. +The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as +```julia +boundary_normal_flux_function(x, t, equations) +``` +where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations. +""" +struct BoundaryConditionNeumann{B} + boundary_normal_flux_function::B +end # set sensible default values that may be overwritten by specific equations """ @@ -329,5 +348,6 @@ include("lattice_boltzmann_3d.jl") abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("acoustic_perturbation_2d.jl") +abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end end # @muladd diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl new file mode 100644 index 00000000000..340cf0e4294 --- /dev/null +++ b/src/equations/equations_parabolic.jl @@ -0,0 +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 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 new file mode 100644 index 00000000000..2f1afe25a6d --- /dev/null +++ b/src/equations/laplace_diffusion_2d.jl @@ -0,0 +1,58 @@ +@doc raw""" + LaplaceDiffusion2D(diffusivity, equations) + +`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} <: AbstractLaplaceDiffusion{2, N} + diffusivity::T + equations_hyperbolic::E +end + +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_hyperbolic) + +# no orientation specified since the flux is vector-valued +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: 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_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)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_value_function(x, t, equations_parabolic) +end + +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations_parabolic::LaplaceDiffusion2D) + return flux_inner +end + +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations_parabolic::LaplaceDiffusion2D) + return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic) +end + +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations_parabolic::LaplaceDiffusion2D) + return flux_inner +end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl new file mode 100644 index 00000000000..c1f9534162a --- /dev/null +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -0,0 +1,264 @@ +# 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 + + +""" + SemidiscretizationHyperbolicParabolic + +A struct containing everything needed to describe a spatial semidiscretization +of a mixed hyperbolic-parabolic conservation law. +""" +struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, + BoundaryConditions, BoundaryConditionsParabolic, + SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} <: AbstractSemidiscretization + + mesh::Mesh + + equations::Equations + equations_parabolic::EquationsParabolic + + # This guy is a bit messy since we abuse it as some kind of "exact solution" + # although this doesn't really exist... + initial_condition::InitialCondition + + boundary_conditions::BoundaryConditions + boundary_conditions_parabolic::BoundaryConditionsParabolic + + source_terms::SourceTerms + + solver::Solver + solver_parabolic::SolverParabolic + + cache::Cache + cache_parabolic::CacheParabolic + + performance_counter::PerformanceCounterList{2} + + function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}( + mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, + boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic, + source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} + @assert ndims(mesh) == ndims(equations) + + # Todo: assert nvariables(equations)==nvariables(equations_parabolic) + + performance_counter = PerformanceCounterList{2}() + + new(mesh, equations, equations_parabolic, initial_condition, + boundary_conditions, boundary_conditions_parabolic, + source_terms, solver, solver_parabolic, cache, cache_parabolic, performance_counter) + end +end + +""" + SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver; + solver_parabolic=default_parabolic_solver(), + source_terms=nothing, + both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), + RealT=real(solver), + uEltype=RealT, + both_initial_caches=(NamedTuple(), NamedTuple())) + +Construct a semidiscretization of a hyperbolic-parabolic PDE. +""" +function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, + initial_condition, solver; + solver_parabolic=default_parabolic_solver(), + source_terms=nothing, + boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT=real(solver), uEltype=RealT, + initial_caches=(NamedTuple(), NamedTuple())) + + equations_hyperbolic, equations_parabolic = equations + boundary_conditions_hyperbolic, boundary_conditions_parabolic = boundary_conditions + initial_hyperbolic_cache, initial_cache_parabolic = initial_caches + + return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic, + initial_condition, solver; solver_parabolic, source_terms, + boundary_conditions=boundary_conditions_hyperbolic, + boundary_conditions_parabolic=boundary_conditions_parabolic, + RealT, uEltype, initial_cache=initial_hyperbolic_cache, + initial_cache_parabolic=initial_cache_parabolic) +end + +function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, + initial_condition, solver; + solver_parabolic=default_parabolic_solver(), + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + boundary_conditions_parabolic=boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT=real(solver), uEltype=RealT, + initial_cache=NamedTuple(), + initial_cache_parabolic=NamedTuple()) + + cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) + _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, equations_parabolic, + solver, solver_parabolic, RealT, uEltype)..., + initial_cache_parabolic...) + + SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), + typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic), + typeof(source_terms), typeof(solver), typeof(solver_parabolic), typeof(cache), typeof(cache_parabolic)}( + mesh, equations, equations_parabolic, initial_condition, + _boundary_conditions, _boundary_conditions_parabolic, source_terms, + solver, solver_parabolic, cache, cache_parabolic) +end + + +# Create a new semidiscretization but change some parameters compared to the input. +# `Base.similar` follows a related concept but would require us to `copy` the `mesh`, +# which would impact the performance. Instead, `SciMLBase.remake` has exactly the +# semantics we want to use here. In particular, it allows us to re-use mutable parts, +# e.g. `remake(semi).mesh === semi.mesh`. +function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.solver), + mesh=semi.mesh, + equations=semi.equations, + equations_parabolic=semi.equations_parabolic, + initial_condition=semi.initial_condition, + solver=semi.solver, + solver_parabolic=semi.solver_parabolic, + source_terms=semi.source_terms, + boundary_conditions=semi.boundary_conditions, + boundary_conditions_parabolic=semi.boundary_conditions_parabolic + ) + # TODO: Which parts do we want to `remake`? At least the solver needs some + # special care if shock-capturing volume integrals are used (because of + # the indicators and their own caches...). + SemidiscretizationHyperbolicParabolic( + mesh, equations, equations_parabolic, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) +end + +function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) + @nospecialize semi # reduce precompilation time + + print(io, "SemidiscretizationHyperbolicParabolic(") + print(io, semi.mesh) + print(io, ", ", semi.equations) + print(io, ", ", semi.equations_parabolic) + print(io, ", ", semi.initial_condition) + print(io, ", ", semi.boundary_conditions) + print(io, ", ", semi.boundary_conditions_parabolic) + print(io, ", ", semi.source_terms) + print(io, ", ", semi.solver) + print(io, ", ", semi.solver_parabolic) + print(io, ", cache(") + for (idx,key) in enumerate(keys(semi.cache)) + idx > 1 && print(io, " ") + print(io, key) + end + print(io, "))") +end + +function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperbolicParabolic) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationHyperbolicParabolic") + summary_line(io, "#spatial dimensions", ndims(semi.equations)) + summary_line(io, "mesh", semi.mesh) + summary_line(io, "hyperbolic equations", semi.equations |> typeof |> nameof) + summary_line(io, "parabolic equations", semi.equations_parabolic |> typeof |> nameof) + summary_line(io, "initial condition", semi.initial_condition) + + # print_boundary_conditions(io, semi) + + summary_line(io, "source terms", semi.source_terms) + summary_line(io, "solver", semi.solver |> typeof |> nameof) + summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof) + summary_line(io, "total #DOFs", ndofs(semi)) + summary_footer(io) + end +end + +@inline Base.ndims(semi::SemidiscretizationHyperbolicParabolic) = ndims(semi.mesh) + +@inline nvariables(semi::SemidiscretizationHyperbolicParabolic) = nvariables(semi.equations) + +@inline Base.real(semi::SemidiscretizationHyperbolicParabolic) = real(semi.solver) + +# retain dispatch on hyperbolic equations only +@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic) + @unpack mesh, equations, solver, cache = semi + return mesh, equations, solver, cache +end + + +function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) + @unpack mesh, equations, initial_condition, solver, cache = semi + u = wrap_array(u_ode, mesh, equations, solver, cache) + + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) +end + + +function compute_coefficients(t, semi::SemidiscretizationHyperbolicParabolic) + # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` + compute_coefficients(semi.initial_condition, t, semi) +end + +function compute_coefficients!(u_ode, t, semi::SemidiscretizationHyperbolicParabolic) + compute_coefficients!(u_ode, semi.initial_condition, t, semi) +end + +""" + semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan) + +Wrap the semidiscretization `semi` as a Split ODE problem in the time interval `tspan` +that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). +""" +function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan) + u0_ode = compute_coefficients(first(tspan), semi) + # TODO: MPI, do we want to synchonize loading and print debug statements, e.g. using + # mpi_isparallel() && MPI.Barrier(mpi_comm()) + # See https://github.com/trixi-framework/Trixi.jl/issues/328 + iip = true # is-inplace, i.e., we modify a vector when calling rhs!, rhs_parabolic! + return SplitODEProblem{iip}(rhs!, rhs_parabolic!, u0_ode, tspan, semi) +end + +function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) + @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi + + u = wrap_array(u_ode, mesh, equations, solver, cache) + du = wrap_array(du_ode, mesh, equations, solver, cache) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "rhs!" rhs!(du, u, t, mesh, equations, initial_condition, + boundary_conditions, source_terms, solver, cache) + runtime = time_ns() - time_start + put!(semi.performance_counter.counters[1], runtime) + + return nothing +end + +function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) + @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi + + u = wrap_array(u_ode, mesh, equations_parabolic, solver, cache_parabolic) + du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition, + boundary_conditions_parabolic, source_terms, + solver, solver_parabolic, cache, cache_parabolic) + runtime = time_ns() - time_start + put!(semi.performance_counter.counters[2], runtime) + + return nothing +end + + +end # @muladd diff --git a/src/solvers/dgmulti.jl b/src/solvers/dgmulti.jl index d5c8f12d395..318a11b678e 100644 --- a/src/solvers/dgmulti.jl +++ b/src/solvers/dgmulti.jl @@ -9,3 +9,6 @@ include("dgmulti/sbp.jl") # specialization of DGMulti to specific equations include("dgmulti/flux_differencing_compressible_euler.jl") + +# parabolic terms for DGMulti solvers +include("dgmulti/dg_parabolic.jl") \ No newline at end of file diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 870533420a6..25135b784f2 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -31,33 +31,32 @@ mul_by_accum!(A::UniformScaling) = MulByAccumUniformScaling() # solution storage formats. @inline apply_to_each_field(f::MulByUniformScaling, out, x, args...) = copy!(out, x) @inline function apply_to_each_field(f::MulByAccumUniformScaling, out, x, args...) - # TODO: DGMulti speed up using threads - for (i, x_i) in enumerate(x) - out[i] = out[i] + x_i + @threaded for i in eachindex(x) + out[i] = out[i] + x[i] end end @inline eachdim(mesh) = Base.OneTo(ndims(mesh)) # iteration over all elements in a mesh -@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, cache) = dg.basis.Np * mesh.md.num_elements -@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(mesh.md.num_elements) +@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, other_args...) = dg.basis.Np * mesh.md.num_elements +@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(mesh.md.num_elements) # iteration over quantities in a single element @inline nnodes(basis::RefElemData) = basis.Np -@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq) -@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq) +@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq) +@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq) # iteration over quantities over the entire mesh (dofs, quad nodes, face nodes). -@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(ndofs(mesh, dg, cache)) -@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements) -@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements) +@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(ndofs(mesh, dg, other_args...)) +@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements) +@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements) # interface with semidiscretization_hyperbolic wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode -function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, mesh::DGMultiMesh, - dg::DGMulti, cache) where {Keys,ValueTypes<:NTuple{N,Any}} where {N} +function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueTypes}, mesh::DGMultiMesh, + dg::DGMulti, cache) where {Keys, ValueTypes<:NTuple{N, Any}} where {N} return boundary_conditions end @@ -67,7 +66,7 @@ function allocate_nested_array(uEltype, nvars, array_dimensions, dg) return StructArray{SVector{nvars, uEltype}}(ntuple(_->zeros(uEltype, array_dimensions...), nvars)) end -function reset_du!(du, dg::DGMulti, cache) +function reset_du!(du, dg::DGMulti, other_args...) @threaded for i in eachindex(du) du[i] = zero(eltype(du)) end @@ -182,6 +181,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 @@ -231,13 +231,9 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm, # inner (idM -> minus) and outer (idP -> plus) indices idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] - - # compute flux if node is not a boundary node - if idM != idP - uP = u_face_values[idP] - normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM] - flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM] - end + uP = u_face_values[idP] + normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM] + flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM] end end @@ -287,6 +283,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 @@ -376,9 +373,10 @@ end # Todo: DGMulti. Specialize for modal DG on curved meshes using WADG -function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache) +# inverts Jacobian and scales by -1.0 +function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1) @threaded for i in each_dof_global(mesh, dg, cache) - du[i] *= -cache.invJ[i] + du[i] *= scaling * cache.invJ[i] end end diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl new file mode 100644 index 00000000000..50cfd8ab17d --- /dev/null +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -0,0 +1,322 @@ +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) + + # u_transformed stores "transformed" variables for computing the gradient + @unpack md = mesh + u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) + 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) + 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_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 + inv_h = similar(mesh.md.Jf) + J = dg.basis.Vf * mesh.md.J # interp to face nodes + for e in eachelement(mesh, dg) + for i in each_face_node(mesh, dg) + inv_h[i, e] = mesh.md.Jf[i, e] / J[i, e] + end + end + + return (; u_transformed, gradients, flux_viscous, + weak_differentiation_matrices, inv_h, + 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, mesh, equations_parabolic::AbstractEquationsParabolic, + dg::DGMulti, parabolic_scheme, cache, cache_parabolic) + @threaded for i in eachindex(u) + 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(gradients, u, scalar_flux_face_values, + mesh, equations::AbstractEquationsParabolic, + dg::DGMulti, cache, cache_parabolic) + @unpack local_flux_face_values_threaded = cache_parabolic + @threaded for e in eachelement(mesh, dg) + local_flux_values = local_flux_face_values_threaded[Threads.threadid()] + for dim in eachdim(mesh) + for i in eachindex(local_flux_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(gradients[dim], :, e), local_flux_values) + end + end +end + +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(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: DGMulti. Assumes mesh is affine here. + apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), + view(gradients[i], :, e), view(u, :, e)) + end + end + + @unpack u_face_values = cache_parabolic + prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + + # compute fluxes at interfaces + @unpack scalar_flux_face_values = cache_parabolic + @unpack mapM, mapP, Jf = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg) + idM, idP = mapM[face_node_index], mapP[face_node_index] + uM = u_face_values[idM] + uP = u_face_values[idP] + 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, u_face_values, t, Gradient(), boundary_conditions, + mesh, equations, dg, cache, cache_parabolic) + + # compute surface contributions + calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, + mesh, equations, dg, cache, cache_parabolic) + + for dim in eachdim(mesh) + invert_jacobian!(gradients[dim], mesh, equations, dg, cache; scaling=1.0) + end +end + +# do nothing for periodic domains +function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, + mesh, equations::AbstractEquationsParabolic, dg::DGMulti, + cache, cache_parabolic) + return nothing +end + +# "lispy tuple programming" instead of for loop for type stability +function calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions, + mesh, equations, dg::DGMulti, cache, cache_parabolic) + + # peel off first boundary condition + calc_single_boundary_flux!(flux, u, t, operator_type, first(boundary_conditions), first(keys(boundary_conditions)), + mesh, equations, dg, cache, cache_parabolic) + + # recurse on the remainder of the boundary conditions + calc_boundary_flux!(flux, u, t, operator_type, Base.tail(boundary_conditions), + mesh, equations, dg, cache, cache_parabolic) +end + +# terminate recursion +calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, + mesh, equations, dg::DGMulti, cache, cache_parabolic) = nothing + +# TODO: DGMulti. Decide if we want to use the input `u_face_values` (currently unused) +function calc_single_boundary_flux!(flux_face_values, u_face_values, t, + operator_type, boundary_condition, boundary_key, + mesh, equations, dg::DGMulti{NDIMS}, cache, cache_parabolic) where {NDIMS} + rd = dg.basis + md = mesh.md + + num_pts_per_face = rd.Nfq ÷ rd.Nfaces + @unpack xyzf, nxyzJ, Jf = md + for f in mesh.boundary_faces[boundary_key] + for i in Base.OneTo(num_pts_per_face) + + # reverse engineer element + face node indices (avoids reshaping arrays) + e = ((f-1) ÷ rd.Nfaces) + 1 + fid = i + ((f-1) % rd.Nfaces) * num_pts_per_face + + face_normal = SVector{NDIMS}(getindex.(nxyzJ, fid, e)) / Jf[fid,e] + face_coordinates = SVector{NDIMS}(getindex.(xyzf, fid, e)) + + # 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. + 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!(flux_viscous, u, gradients, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, + dg::DGMulti, cache, cache_parabolic) + + for dim in eachdim(mesh) + reset_du!(flux_viscous[dim], dg) + end + + @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_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_flux_viscous[dim], zero(eltype(local_flux_viscous[dim]))) + end + + # 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_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] + 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(flux_viscous[dim], :, e), local_flux_viscous[dim]) + end + end +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, + 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, + parabolic_scheme, cache, cache_parabolic) + # compute fluxes at interfaces + @unpack scalar_flux_face_values, inv_h = cache_parabolic + @unpack mapM, mapP = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg) + 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, parabolic_scheme) + end + return nothing +end + + +function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, + boundary_conditions, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) + + @unpack weak_differentiation_matrices = cache_parabolic + + reset_du!(du, dg) + + # compute volume contributions to divergence + @threaded for e in eachelement(mesh, dg) + 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(flux_viscous[i], :, e)) + end + end + + # interpolates from solution coefficients to face quadrature points + flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage + for dim in eachdim(mesh) + prolong2interfaces!(flux_viscous_face_values[dim], flux_viscous[dim], mesh, equations, + dg.surface_integral, dg, cache) + end + + # compute fluxes at interfaces + @unpack scalar_flux_face_values = cache_parabolic + @unpack mapM, mapP, nxyzJ = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg, cache, cache_parabolic) + idM, idP = mapM[face_node_index], mapP[face_node_index] + + # compute f(u, ∇u) ⋅ n + flux_face_value = zero(eltype(scalar_flux_face_values)) + for dim in eachdim(mesh) + 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 + + 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, 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 + +# 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(u) +# 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, parabolic_scheme, cache, cache_parabolic) + + reset_du!(du, dg) + + @unpack u_transformed, gradients, flux_viscous = cache_parabolic + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, parabolic_scheme, cache, cache_parabolic) + + calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, + boundary_conditions, dg, cache, cache_parabolic) + + calc_viscous_fluxes!(flux_viscous, u_transformed, gradients, + mesh, equations_parabolic, dg, cache, cache_parabolic) + + calc_divergence!(du, u_transformed, t, flux_viscous, mesh, equations_parabolic, + boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) + + return nothing + +end diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl index 6f73d7f44a2..d1b40702024 100644 --- a/src/solvers/dgmulti/sbp.jl +++ b/src/solvers/dgmulti/sbp.jl @@ -435,6 +435,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 c7f8ed1ddd7..469bc072a31 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -363,6 +363,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..ca6394172ad --- /dev/null +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -0,0 +1,605 @@ +# 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) + # + # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have + # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # `interfaces` and `boundaries` containers in `cache_parabolic` to interpolate and store the + # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we + # do not need to recreate the existing data structure only with a different name, and c) we do not + # need to interpolate solutions *and* gradients to the surfaces. + + # TODO: parabolic; reconsider current data structure reuse strategy + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" calc_volume_integral!( + du, flux_viscous, mesh, equations_parabolic, dg, cache) + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" calc_interface_flux!( + cache_parabolic.elements.surface_flux_values, mesh, equations_parabolic, dg, cache_parabolic) + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!( + cache_parabolic, flux_viscous, mesh, equations_parabolic, dg.surface_integral, dg, cache) + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" calc_boundary_flux_divergence!( + cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic, + dg.surface_integral, dg) + + # TODO: parabolic; extend to mortars + @assert nmortars(dg, cache) == 0 + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" calc_surface_integral!( + du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic) + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!( + du, mesh, equations_parabolic, dg, cache_parabolic) + + return nothing +end + +# Transform solution variables prior to taking the gradient +# (e.g., conservative to primitive variables). Defaults to doing nothing. +# TODO: can we avoid copying data? +function transform_variables!(u_transformed, u, mesh::TreeMesh{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) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + 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) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + 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) + # Get precomputed fluxes at interfaces + flux_ll, flux_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, + dg, i, interface) + + # Compute interface flux as mean of left and right viscous fluxes + # TODO: parabolic; only BR1 at the moment + flux = 0.5 * (flux_ll + flux_rr) + + # Copy flux to left and right element storage + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, 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) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + 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) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + 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) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + 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) + # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! + 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_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_gradient!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_boundary_flux_by_direction_gradient!(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_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[1], + equations_parabolic, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[2], + equations_parabolic, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[3], + equations_parabolic, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction_divergence!(surface_flux_values, t, boundary_conditions_parabolic[4], + equations_parabolic, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) +end +function calc_boundary_flux_by_direction_divergence!(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!` + @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/src/solvers/solvers.jl b/src/solvers/solvers.jl index a52f717a612..465e051482c 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -4,9 +4,10 @@ # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin +# define types for parabolic solvers +include("solvers_parabolic.jl") include("dg.jl") include("dgmulti.jl") - end # @muladd diff --git a/src/solvers/solvers_parabolic.jl b/src/solvers/solvers_parabolic.jl new file mode 100644 index 00000000000..f253cdbd03d --- /dev/null +++ b/src/solvers/solvers_parabolic.jl @@ -0,0 +1,31 @@ +""" + ViscousFormulationBassiRebay1() + +The classical BR1 flux from + +- F. Bassi, S. Rebay (1997) + A High-Order Accurate Discontinuous Finite Element Method for + the Numerical Solution of the Compressible Navier-Stokes Equations + [DOI: 10.1006/jcph.1996.5572](https://doi.org/10.1006/jcph.1996.5572) +""" +struct ViscousFormulationBassiRebay1 end + +""" + ViscousFormulationLocalDG(penalty_parameter) + +The local DG (LDG) flux from "The Local Discontinuous Galerkin Method for Time-Dependent +Convection-Diffusion Systems" by Cockburn and Shu (1998). + +Note that, since this implementation does not involve the parabolic "upwinding" vector, +the LDG solver is equivalent to [`ViscousFormulationBassiRebay1`](@ref) with an LDG-type penalization. + +- Cockburn and Shu (1998). + The Local Discontinuous Galerkin Method for Time-Dependent + Convection-Diffusion Systems + [DOI: 10.1137/S0036142997316712](https://doi.org/10.1137/S0036142997316712) +""" +struct ViscousFormulationLocalDG{P} + penalty_parameter::P +end + +default_parabolic_solver() = ViscousFormulationBassiRebay1() \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index e294d10746b..79a8e54e791 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -86,6 +86,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_dgmulti_3d.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" + include("test_parabolic_2d.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" include("test_unit.jl") include("test_visualization.jl") diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 2651ae06ee6..d3aa1ab0f7f 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -12,6 +12,7 @@ outdir = "out" isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin + @trixi_testset "elixir_euler_weakform.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), cells_per_dimension = (4, 4), diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl new file mode 100644 index 00000000000..d9c3967315a --- /dev/null +++ b/test/test_parabolic_2d.jl @@ -0,0 +1,185 @@ +module TestExamplesParabolic2D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +@testset "SemidiscretizationHyperbolicParabolic" begin + + @trixi_testset "DGMulti 2D rhs_parabolic!" begin + + dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralWeakForm()) + mesh = DGMultiMesh(dg, cells_per_dimension=(2, 2)) + + # test with polynomial initial condition x^2 * y + # test if we recover the exact second derivative + initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) + + equations = LinearScalarAdvectionEquation2D(1.0, 1.0) + equations_parabolic = LaplaceDiffusion2D(1.0, equations) + + semi = SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, dg) + @test_nowarn_debug show(stdout, semi) + @test_nowarn_debug show(stdout, MIME"text/plain"(), semi) + @test_nowarn_debug show(stdout, boundary_condition_do_nothing) + + @test nvariables(semi)==nvariables(equations) + @test Base.ndims(semi)==Base.ndims(mesh) + @test Base.real(semi)==Base.real(dg) + + ode = semidiscretize(semi, (0.0, 0.01)) + u0 = similar(ode.u0) + Trixi.compute_coefficients!(u0, 0.0, semi) + @test u0 ≈ ode.u0 + + # test "do nothing" BC just returns first argument + @test boundary_condition_do_nothing(u0, nothing) == u0 + + @unpack cache, cache_parabolic, equations_parabolic = semi + @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!(gradients, ode.u0, t, mesh, equations_parabolic, + boundary_condition_periodic, dg, cache, cache_parabolic) + @unpack x, y = mesh.md + @test getindex.(gradients[1], 1) ≈ 2 * x .* y + @test getindex.(gradients[2], 1) ≈ x.^2 + + u_flux = similar.(gradients) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, gradients, mesh, equations_parabolic, + dg, cache, cache_parabolic) + @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, + dg, semi.solver_parabolic, cache, cache_parabolic) + @test getindex.(du, 1) ≈ 2 * y + end + + @trixi_testset "DGMulti: elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.2485803335154642], + linf = [1.079606969242132] + ) + end + + @trixi_testset "DGMulti: elixir_advection_diffusion_periodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_advection_diffusion_periodic.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.03180371984888462], + linf = [0.2136821621370909] + ) + end + + @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.002123168335604323], + linf = [0.00963640423513712] + ) + end + + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.0015355076812510957, 0.0033843168272696756, 0.0036531858107443434, 0.009948436427519214], + linf = [0.005522560467190019, 0.013425258500730508, 0.013962115643482154, 0.027483102120502423] + ) + end + + @trixi_testset "DGMulti: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_lid_driven_cavity.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.5), + l2 = [0.00022156125227115747, 0.028318325921401, 0.009509168701070296, 0.028267900513550506], + linf = [0.001562278941298234, 0.14886653390744856, 0.0716323565533752, 0.19472785105241996] + ) + 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_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + l2 = [0.002111672530658797, 0.0034322351490857846, 0.0038742528195910416, 0.012469246082568561], + linf = [0.012006418939223495, 0.035520871209746126, 0.024512747492231427, 0.11191122588756564] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.002103629650383915, 0.003435843933396454, 0.00386735987813341, 0.012670355349235728], + linf = [0.012006261793147788, 0.03550212518982032, 0.025107947319661185, 0.11647078036571124] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (Entropy gradient variables)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + l2 = [0.0021403742517389513, 0.0034258287094908572, 0.0038915122886898517, 0.012506862343013842], + linf = [0.012244412004628336, 0.03507559186162224, 0.024580892345558894, 0.11425600758350107] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), + l2 = [0.0021349737347844907, 0.0034301388278203033, 0.0038928324474291572, 0.012693611436230873], + linf = [0.01224423627586213, 0.035054066314102905, 0.025099598504931965, 0.11795616324751634] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (flux differencing)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.1), + volume_integral=VolumeIntegralFluxDifferencing(flux_central), + l2 = [0.0021116725306633594, 0.0034322351490827557, 0.0038742528196093542, 0.012469246082526909], + linf = [0.012006418939291663, 0.035520871209594115, 0.024512747491801577, 0.11191122588591007] + ) + end + + @trixi_testset "TreeMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.00015144571529699053, 0.018766076072331623, 0.007065070765652574, 0.0208399005734258], + linf = [0.0014523369373669048, 0.12366779944955864, 0.05532450997115432, 0.16099927805328207] + ) + end + +end + +# Clean up afterwards: delete Trixi output directory +@test_nowarn isdir(outdir) && rm(outdir, recursive=true) + +end # module