diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl new file mode 100644 index 00000000000..55682f73fce --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl @@ -0,0 +1,96 @@ +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) + +# Example setup taken from +# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016). +# Robust DPG methods for transient convection-diffusion. +# In: Building bridges: connections and challenges in modern approaches +# to numerical partial differential equations. +# [DOI](https://doi.org/10.1007/978-3-319-41640-3_6). +function initial_condition_eriksson_johnson(x, t, equations) + l = 4 + epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_eriksson_johnson + +boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), + :y_neg => BoundaryConditionDirichlet(initial_condition), + :y_pos => BoundaryConditionDirichlet(initial_condition), + :x_pos => boundary_condition_do_nothing) + +boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), + :x_pos => BoundaryConditionDirichlet(initial_condition), + :y_neg => BoundaryConditionDirichlet(initial_condition), + :y_pos => BoundaryConditionDirichlet(initial_condition)) + +# 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) +coordinates_max = ( 0.0, 0.5) + +# This maps the domain [-1, 1]^2 to [-1, 0] x [-0.5, 0.5] while also +# introducing a curved warping to interior nodes. +function mapping(xi, eta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(0.5 * (1 + x) - 1, 0.5 * y) +end + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg=3, initial_refinement_level=2, + mapping=mapping, periodicity=(false, false)) + +# 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.0) +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, + ode_default_options()..., callback=callbacks) + +# Print the timer summary +summary_callback() diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 00000000000..8111df8251a --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,209 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg=3, initial_refinement_level=2, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity=(true, false)) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + y = x[2] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + # stress tensor from y-direction + - v1_xy * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * mu_ ) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = Dict(:y_neg => boundary_condition_top_bottom, + :y_pos => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +# ############################################################################### +# # ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary + diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl new file mode 100644 index 00000000000..051f4defe54 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -0,0 +1,82 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 0.001 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), + Prandtl=prandtl_number()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg=3, initial_refinement_level=2, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity=(false, false)) + +function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) + Ma = 0.1 + rho = 1.0 + u, v = 0.0, 0.0 + p = 1.0 / (Ma^2 * equations.gamma) + return prim2cons(SVector(rho, u, v, p), equations) +end +initial_condition = initial_condition_cavity + +# BC types +velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0)) +velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc) +boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc) + +# define periodic boundary conditions everywhere +boundary_conditions = Dict( :x_neg => boundary_condition_slip_wall, + :y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall, + :x_pos => boundary_condition_slip_wall) + +boundary_conditions_parabolic = Dict( :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, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary + + diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index e73a8cda9b8..73ac47ed1e3 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -365,6 +365,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, return nothing end +# This version is used for divergence flux computations function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, equations_parabolic, dg::DG, cache_parabolic) @@ -405,7 +406,7 @@ function calc_interface_flux!(surface_flux_values, end for node in eachnode(dg) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the viscous flux dotted with respect the outward normal on the # primary element. We assume a BR-1 type of flux. viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, @@ -446,6 +447,7 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # a start value and a step size to get the correct face and orientation. element = boundaries.neighbor_ids[boundary] node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) @@ -454,15 +456,12 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, j_node = j_node_start for i in eachnode(dg) # this is the outward normal direction on the primary element - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_node, j_node, primary_element) + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, - primary_element], - flux_viscous_y[v, i_primary, j_primary, - primary_element]) + flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], + flux_viscous_y[v, i_node, j_node, element]) boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction) end @@ -470,6 +469,124 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, j_node += j_node_step end end + return nothing +end + +function calc_boundary_flux_gradients!(cache, t, + boundary_condition::Union{BoundaryConditionPeriodic, + BoundaryConditionDoNothing + }, + mesh::P4estMesh, equations, surface_integral, dg::DG) + @assert isempty(eachboundary(dg, cache)) +end + +# Function barrier for type stability +function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Gradient(), mesh, equations, surface_integral, dg) return nothing end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Divergence(), mesh, equations, surface_integral, dg) + return nothing +end + +# Iterate over tuples of boundary condition types and associated indices +# in a type-stable way using "lispy tuple programming". +function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, + BC_indices::NTuple{N, Vector{Int}}, + operator_type, + mesh::P4estMesh, + equations, surface_integral, dg::DG) where {N} + # Extract the boundary condition type and index vector + boundary_condition = first(BCs) + boundary_condition_indices = first(BC_indices) + # Extract the remaining types and indices to be processed later + remaining_boundary_conditions = Base.tail(BCs) + remaining_boundary_condition_indices = Base.tail(BC_indices) + + # process the first boundary condition type + calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, + operator_type, mesh, equations, surface_integral, dg) + + # recursively call this method with the unprocessed boundary types + calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions, + remaining_boundary_condition_indices, + operator_type, + mesh, equations, surface_integral, dg) + + return nothing +end + +# terminate the type-stable iteration over tuples +function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, + operator_type, mesh::P4estMesh, equations, + surface_integral, dg::DG) + nothing +end + +function calc_boundary_flux!(cache, t, + boundary_condition_parabolic, # works with Dict types + boundary_condition_indices, + operator_type, mesh::P4estMesh{2}, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + (; boundaries) = cache + (; node_coordinates, surface_flux_values) = cache.elements + (; contravariant_vectors) = cache.elements + index_range = eachnode(dg) + + @threaded for local_index in eachindex(boundary_condition_indices) + # Use the local index to get the global boundary index from the pre-sorted list + boundary_index = boundary_condition_indices[local_index] + + # Get information on the adjacent element, compute the surface fluxes, + # and store them + element = boundaries.neighbor_ids[boundary_index] + node_indices = boundaries.node_indices[boundary_index] + direction_index = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in eachnode(dg) + # Extract solution data from boundary container + u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index, + boundary_index) + + # Outward-pointing normal direction (not normalized) + normal_direction = get_normal_direction(direction_index, contravariant_vectors, + i_node, j_node, element) + + # 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 + + # Coordinates at boundary node + x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, + element) + + flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, + x, t, operator_type, equations_parabolic) + + # Copy flux to element storage in the correct orientation + for v in eachvariable(equations_parabolic) + surface_flux_values[v, node_index, direction_index, element] = flux_[v] + end + + i_node += i_node_step + j_node += j_node_step + end + end +end diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index b0ac63d4ce9..471b976e990 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -200,6 +200,38 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"), + trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.012380458938507371], + linf = [0.10860506906472567] + ) + end + + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_curved.jl"), + trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.04933902988507035], + linf = [0.2550261714590271] + ) + end + + @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 1, tspan=(0.0, 0.2), + l2 = [0.0003811978985836709, 0.0005874314969169538, 0.0009142898787923481, 0.0011613918899727263], + linf = [0.0021633623982135752, 0.009484348274135372, 0.004231572066492217, 0.011661660275365193] + ) + end + + @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.00028716166408816073, 0.08101204560401647, 0.02099595625377768, 0.05008149754143295], + linf = [0.014804500261322406, 0.9513271652357098, 0.7223919625994717, 1.4846907331004786] + ) + end + end # Clean up afterwards: delete Trixi.jl output directory