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