diff --git a/Project.toml b/Project.toml index d318389a6d2..4efbc40b098 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.45-pre" +version = "0.5.46-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl new file mode 100644 index 00000000000..935f132ba4b --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl @@ -0,0 +1,215 @@ +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=(false, 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) + +boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test) + +# define inviscid boundary conditions +boundary_conditions = Dict(:x_neg => boundary_condition_left_right, + :x_pos => boundary_condition_left_right, + :y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right, + :x_pos => boundary_condition_left_right, + :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/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index fad42b11098..a3e49007fb6 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -297,7 +297,7 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi) mpi_println(" " * " " * " " * " PID: " * @sprintf("%10.8e s", performance_index)) - mpi_println(" #DOF: " * @sprintf("% 14d", ndofs(semi)) * + mpi_println(" #DOFs per field:" * @sprintf("% 14d", ndofs(semi)) * " " * " alloc'd memory: " * @sprintf("%14.3f MiB", memory_use)) mpi_println(" #elements: " * diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index b10ffa3b9d3..80857999017 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -456,4 +456,33 @@ end }) return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4]) end + +# Dirichlet Boundary Condition for P4est mesh + +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion2D{ + GradientVariablesPrimitive + }) + # BCs are usually specified as conservative variables so we convert them to primitive variables + # because the gradients are assumed to be with respect to the primitive variables + u_boundary = boundary_condition.boundary_value_function(x, t, equations) + + return cons2prim(u_boundary, equations) +end + +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, + u_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion2D{ + GradientVariablesPrimitive + }) + # for Dirichlet boundary conditions, we do not impose any conditions on the viscous fluxes + return flux_inner +end end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b7adff78425..49763b12e5d 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -71,7 +71,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled) summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms) summary_line(increment_indent(io), "solver", solver |> typeof |> nameof) end - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 50b2c21c14e..9d465bfcc5f 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -243,7 +243,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_line(io, "source terms", semi.source_terms) summary_line(io, "solver", semi.solver |> typeof |> nameof) - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index b12ecadb58b..35340d65b1e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -228,7 +228,7 @@ function Base.show(io::IO, ::MIME"text/plain", summary_line(io, "source terms", semi.source_terms) summary_line(io, "solver", semi.solver |> typeof |> nameof) summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof) - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl index ba02d812041..d434d3146ce 100644 --- a/src/solvers/dgmulti/sbp.jl +++ b/src/solvers/dgmulti/sbp.jl @@ -36,312 +36,6 @@ function DGMulti(element_type::AbstractElemShape, surface_integral = surface_integral, volume_integral = volume_integral) end -function construct_1d_operators(D::AbstractDerivativeOperator, tol) - nodes_1d = collect(grid(D)) - M = SummationByPartsOperators.mass_matrix(D) - if M isa UniformScaling - weights_1d = M * ones(Bool, length(nodes_1d)) - else - weights_1d = diag(M) - end - - # StartUpDG assumes nodes from -1 to +1. Thus, we need to re-scale everything. - # We can adjust the grid spacing as follows. - xmin = SummationByPartsOperators.xmin(D) - xmax = SummationByPartsOperators.xmax(D) - factor = 2 / (xmax - xmin) - @. nodes_1d = factor * (nodes_1d - xmin) - 1 - @. weights_1d = factor * weights_1d - - D_1d = droptol!(inv(factor) * sparse(D), tol) - I_1d = Diagonal(ones(Bool, length(nodes_1d))) - - return nodes_1d, weights_1d, D_1d, I_1d -end - -function StartUpDG.RefElemData(element_type::Line, - D::AbstractDerivativeOperator; - tol = 100 * eps()) - approximation_type = D - N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree - - # 1D operators - nodes_1d, weights_1d, D_1d = construct_1d_operators(D, tol) - - # volume - rq = r = nodes_1d - wq = weights_1d - Dr = D_1d - M = Diagonal(wq) - Pq = LinearAlgebra.I - Vq = LinearAlgebra.I - - VDM = nothing # unused generalized Vandermonde matrix - - rst = (r,) - rstq = (rq,) - Drst = (Dr,) - - # face - face_vertices = StartUpDG.face_vertices(element_type) - face_mask = [1, length(nodes_1d)] - - rf = [-1.0; 1.0] - nrJ = [-1.0; 1.0] - wf = [1.0; 1.0] - if D isa AbstractPeriodicDerivativeOperator - # we do not need any face stuff for periodic operators - Vf = spzeros(length(wf), length(wq)) - else - Vf = sparse([1, 2], [1, length(nodes_1d)], [1.0, 1.0]) - end - LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) - - rstf = (rf,) - nrstJ = (nrJ,) - - # low order interpolation nodes - r1 = StartUpDG.nodes(element_type, 1) - V1 = StartUpDG.vandermonde(element_type, 1, r) / - StartUpDG.vandermonde(element_type, 1, r1) - - return RefElemData(element_type, approximation_type, N, - face_vertices, V1, - rst, VDM, face_mask, - rst, LinearAlgebra.I, # plotting - rstq, wq, Vq, # quadrature - rstf, wf, Vf, nrstJ, # faces - M, Pq, Drst, LIFT) -end - -function StartUpDG.RefElemData(element_type::Quad, - D::AbstractDerivativeOperator; - tol = 100 * eps()) - approximation_type = D - N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree - - # 1D operators - nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) - - # volume - s, r = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d)) # this is to match - # ordering of nrstJ - rq = r - sq = s - wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d)) - wq = wr .* ws - Dr = kron(I_1d, D_1d) - Ds = kron(D_1d, I_1d) - M = Diagonal(wq) - Pq = LinearAlgebra.I - Vq = LinearAlgebra.I - - VDM = nothing # unused generalized Vandermonde matrix - - rst = (r, s) - rstq = (rq, sq) - Drst = (Dr, Ds) - - # face - face_vertices = StartUpDG.face_vertices(element_type) - face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s)...) - - rf, sf, wf, nrJ, nsJ = StartUpDG.init_face_data(element_type, - quad_rule_face = (nodes_1d, weights_1d)) - if D isa AbstractPeriodicDerivativeOperator - # we do not need any face stuff for periodic operators - Vf = spzeros(length(wf), length(wq)) - else - Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask))) - end - LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) - - rstf = (rf, sf) - nrstJ = (nrJ, nsJ) - - # low order interpolation nodes - r1, s1 = StartUpDG.nodes(element_type, 1) - V1 = StartUpDG.vandermonde(element_type, 1, r, s) / - StartUpDG.vandermonde(element_type, 1, r1, s1) - - return RefElemData(element_type, approximation_type, N, - face_vertices, V1, - rst, VDM, face_mask, - rst, LinearAlgebra.I, # plotting - rstq, wq, Vq, # quadrature - rstf, wf, Vf, nrstJ, # faces - M, Pq, Drst, LIFT) -end - -function StartUpDG.RefElemData(element_type::Hex, - D::AbstractDerivativeOperator; - tol = 100 * eps()) - approximation_type = D - N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree - - # 1D operators - nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) - - # volume - # to match ordering of nrstJ - s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d)) - rq = r - sq = s - tq = t - wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d)) - wq = wr .* ws .* wt - Dr = kron(I_1d, I_1d, D_1d) - Ds = kron(I_1d, D_1d, I_1d) - Dt = kron(D_1d, I_1d, I_1d) - M = Diagonal(wq) - Pq = LinearAlgebra.I - Vq = LinearAlgebra.I - - VDM = nothing # unused generalized Vandermonde matrix - - rst = (r, s, t) - rstq = (rq, sq, tq) - Drst = (Dr, Ds, Dt) - - # face - face_vertices = StartUpDG.face_vertices(element_type) - face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s, t)...) - - rf, sf, tf, wf, nrJ, nsJ, ntJ = let - rf, sf = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d)) - wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d)) - wf = wr .* ws - StartUpDG.init_face_data(element_type, quad_rule_face = (rf, sf, wf)) - end - Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask))) - LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) - - rstf = (rf, sf, tf) - nrstJ = (nrJ, nsJ, ntJ) - - # low order interpolation nodes - r1, s1, t1 = StartUpDG.nodes(element_type, 1) - V1 = StartUpDG.vandermonde(element_type, 1, r, s, t) / - StartUpDG.vandermonde(element_type, 1, r1, s1, t1) - - return RefElemData(element_type, approximation_type, N, - face_vertices, V1, - rst, VDM, face_mask, - rst, LinearAlgebra.I, # plotting - rstq, wq, Vq, # quadrature - rstf, wf, Vf, nrstJ, # faces - M, Pq, Drst, LIFT) -end - -# specialized Hex constructor in 3D to reduce memory usage. -function StartUpDG.RefElemData(element_type::Hex, - D::AbstractPeriodicDerivativeOperator; - tol = 100 * eps()) - approximation_type = D - N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree - - # 1D operators - nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) - - # volume - # to match ordering of nrstJ - s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d)) - rq = r - sq = s - tq = t - wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d)) - wq = wr .* ws .* wt - Dr = kron(I_1d, I_1d, D_1d) - Ds = kron(I_1d, D_1d, I_1d) - Dt = kron(D_1d, I_1d, I_1d) - M = Diagonal(wq) - Pq = LinearAlgebra.I - Vq = LinearAlgebra.I - - VDM = nothing # unused generalized Vandermonde matrix - - rst = (r, s, t) - rstq = (rq, sq, tq) - Drst = (Dr, Ds, Dt) - - # face - # We do not need any face data for periodic operators. Thus, we just - # pass `nothing` to save memory. - face_vertices = ntuple(_ -> nothing, 3) - face_mask = nothing - wf = nothing - rstf = ntuple(_ -> nothing, 3) - nrstJ = ntuple(_ -> nothing, 3) - Vf = nothing - LIFT = nothing - - # low order interpolation nodes - V1 = nothing # do not need to store V1, since we specialize StartUpDG.MeshData to avoid using it. - - return RefElemData(element_type, approximation_type, N, - face_vertices, V1, - rst, VDM, face_mask, - rst, LinearAlgebra.I, # plotting - rstq, wq, Vq, # quadrature - rstf, wf, Vf, nrstJ, # faces - M, Pq, Drst, LIFT) -end - -function Base.show(io::IO, mime::MIME"text/plain", - rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {NDIMS, - ElementType <: - StartUpDG.AbstractElemShape, - ApproximationType <: - AbstractDerivativeOperator - } - @nospecialize rd - print(io, "RefElemData for an approximation using an ") - show(IOContext(io, :compact => true), rd.approximation_type) - print(io, " on $(rd.element_type) element") -end - -function Base.show(io::IO, - rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {NDIMS, - ElementType <: - StartUpDG.AbstractElemShape, - ApproximationType <: - AbstractDerivativeOperator - } - @nospecialize rd - print(io, "RefElemData{", summary(rd.approximation_type), ", ", rd.element_type, "}") -end - -function StartUpDG.inverse_trace_constant(rd::RefElemData{NDIMS, ElementType, - ApproximationType}) where {NDIMS, - ElementType <: - Union{ - Line, - Quad, - Hex - }, - ApproximationType <: - AbstractDerivativeOperator - } - D = rd.approximation_type - - # the inverse trace constant is the maximum eigenvalue corresponding to - # M_f * v = λ * M * v - # where M_f is the face mass matrix and M is the volume mass matrix. - # Since M is diagonal and since M_f is just the boundary "mask" matrix - # (which extracts the first and last entries of a vector), the maximum - # eigenvalue is the inverse of the first or last mass matrix diagonal. - left_weight = SummationByPartsOperators.left_boundary_weight(D) - right_weight = SummationByPartsOperators.right_boundary_weight(D) - max_eigenvalue = max(inv(left_weight), inv(right_weight)) - - # For tensor product elements, the trace constant for higher dimensional - # elements is the one-dimensional trace constant multiplied by `NDIMS`. See - # "GPU-accelerated discontinuous Galerkin methods on hybrid meshes." - # Chan, Jesse, et al (2016), https://doi.org/10.1016/j.jcp.2016.04.003 - # for more details (specifically, Appendix A.1, Theorem A.4). - return NDIMS * max_eigenvalue -end - # type alias for specializing on a periodic SBP operator const DGMultiPeriodicFDSBP{NDIMS, ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxType, @@ -450,30 +144,6 @@ end @muladd begin #! format: noindent -# This is used in `estimate_dt`. `estimate_h` uses that `Jf / J = O(h^{NDIMS-1}) / O(h^{NDIMS}) = O(1/h)`. -# However, since we do not initialize `Jf` for periodic FDSBP operators, we specialize `estimate_h` -# based on the reference grid provided by SummationByPartsOperators.jl and information about the domain size -# provided by `md::MeshData``. -function StartUpDG.estimate_h(e, rd::RefElemData{NDIMS, ElementType, ApproximationType}, - md::MeshData) where {NDIMS, - ElementType <: - StartUpDG.AbstractElemShape, - ApproximationType <: - SummationByPartsOperators.AbstractPeriodicDerivativeOperator - } - D = rd.approximation_type - x = grid(D) - - # we assume all SummationByPartsOperators.jl reference grids are rescaled to [-1, 1] - xmin = SummationByPartsOperators.xmin(D) - xmax = SummationByPartsOperators.xmax(D) - factor = 2 / (xmax - xmin) - - # If the domain has size L^NDIMS, then `minimum(md.J)^(1 / NDIMS) = L`. - # WARNING: this is not a good estimate on anisotropic grids. - return minimum(diff(x)) * factor * minimum(md.J)^(1 / NDIMS) -end - # specialized for DGMultiPeriodicFDSBP since there are no face nodes # and thus no inverse trace constant for periodic domains. function estimate_dt(mesh::DGMultiMesh, dg::DGMultiPeriodicFDSBP) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 97447c0f438..22bc99692f2 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -296,6 +296,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence_nonperiodic.jl"), + initial_refinement_level = 1, tspan=(0.0, 0.2), + l2 = [0.00040364962558511795, 0.0005869762481506936, 0.00091488537427274, 0.0011984191566376762], + linf = [0.0024993634941723464, 0.009487866203944725, 0.004505829506628117, 0.011634902776245681] + ) + 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),