diff --git a/NEWS.md b/NEWS.md index 953b15c627d..c8f30e0d3c6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,10 @@ for human readability. the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) and nonlinear variables (using `positivity_variables_nonlinear`) is supported. `BoundsCheckCallback` is not supported in 3D yet. +- Support for user-defined RHS splitting for IMEX methods via `SemidiscretizationHyperbolicSplit` ([#2518]). + The splitting follows the form `y_t = f_1(y) + f_2(y)`, allowing users to define separate solvers + for the stiff (`f_1`) and non-stiff (`f_2`) parts of the right-hand side. Boundary conditions + and source terms can be specified independently for the stiff and non-stiff parts. - Optimized 2D and 3D kernels for nonconservative fluxes with `P4estMesh` were added ([#2653], [#2663]). The optimized kernel can be enabled via the trait `Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)`. When the trait is set to `Trixi.True()`, a single method has to be defined, that computes and returns the tuple diff --git a/examples/p4est_2d_dgsem/elixir_euler_imex_warm_bubble.jl b/examples/p4est_2d_dgsem/elixir_euler_imex_warm_bubble.jl new file mode 100644 index 00000000000..a7e158c8e07 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_imex_warm_bubble.jl @@ -0,0 +1,251 @@ +# This elixir demonstrates how an implicit-explicit (IMEX) time integration scheme can be applied to the stiff and non-stiff parts of a right hand side, respectively. +# We define separate solvers, boundary conditions, and source terms, and create a `SemidiscretizationHyperbolicSplit`, which will return a `SplitODEProblem` compatible with `OrdinaryDiffEqBDF`, cf. https://docs.sciml.ai/OrdinaryDiffEq/stable/implicit/SDIRK/#IMEX-DIRK . +# Note: This is currently more of a proof of concept and not particularly useful in practice, as fully explicit methods are still faster at the moment. + +using Trixi +using OrdinaryDiffEqBDF +using ADTypes # This is needed to set 'autodiff = AutoFiniteDiff()' in the ODE solver. +using LinearSolve + +function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) + g = 9.81 + c_p = 1004.0 + c_v = 717.0 + + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +# The standard Trixi implementation of the slip wall boundary condition is not directly +# compatible with this general splitting approach, since it is based on Toro's Riemann solver +# which always returns boundary condition values for the entire right-hand side. +# This function computes the boundary condition based on the surface flux function of the +# explicit and implicit parts, where the splitting has been defined and thus accounts for it. +@inline function boundary_condition_slip_wall_2(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + # normalize the outward pointing direction + normal = normal_direction / Trixi.norm(normal_direction) + + # compute the normal momentum from + # u = (rho, rho v1, rho v2, rho e) + u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3] + + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + u_inner[2] - 2 * u_normal * normal[1], + u_inner[3] - 2 * u_normal * normal[2], + u_inner[4]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + + return flux +end + +# The total flux is split into: +# - Fast (implicit/stiff) part: Contains all pressure-related terms responsible for acoustic waves. +# Uses LMARS for surface fluxes and Kennedy-Gruber for volume fluxes. +# - Slow (explicit/non-stiff) part: Contains convective terms (advection). +@inline function flux_lmars_fast(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + a = 340.0 + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = Trixi.norm(normal_direction) + + rho = 0.5f0 * (rho_ll + rho_rr) + + p_interface = 0.5f0 * (p_ll + p_rr) - 0.5f0 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f4 = p_ll * v_interface + else + f4 = p_rr * v_interface + end + + return SVector(zero(eltype(u_ll)), + p_interface * normal_direction[1], + p_interface * normal_direction[2], + f4) +end + +@inline function flux_lmars_slow(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + a = 340.0 + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = Trixi.norm(normal_direction) + + rho = 0.5f0 * (rho_ll + rho_rr) + + v_interface = 0.5f0 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f1, f2, f3, f4 = u_ll * v_interface + else + f1, f2, f3, f4 = u_rr * v_interface + end + + return SVector(f1, f2, f3, f4) +end + +@inline function flux_kennedy_gruber_slow(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Unpack left and right state + rho_e_ll = last(u_ll) + rho_e_rr = last(u_rr) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Average each factor of products in flux + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + p_avg = 0.5f0 * (p_ll + p_rr) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_avg * v_dot_n_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * e_avg + + return SVector(f1, f2, f3, f4) +end + +@inline function flux_kennedy_gruber_fast(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Average each factor of products in flux + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + p_avg = 0.5f0 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f2 = p_avg * normal_direction[1] + f3 = p_avg * normal_direction[2] + f4 = p_avg * v_dot_n_avg + + return SVector(zero(eltype(u_ll)), f2, f3, f4) +end + +@inline function source_terms_gravity(u, x, t, equations::CompressibleEulerEquations2D) + g = 9.81 + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +gamma = 1004 / 717 +equations = CompressibleEulerEquations2D(gamma) + +polydeg = 2 +basis = LobattoLegendreBasis(polydeg) + +volume_integral_explicit = VolumeIntegralFluxDifferencing(flux_kennedy_gruber_slow) +solver_explicit = DGSEM(basis, flux_lmars_slow, volume_integral_explicit) + +volume_integral_implicit = VolumeIntegralFluxDifferencing(flux_kennedy_gruber_fast) +solver_implicit = DGSEM(basis, flux_lmars_fast, volume_integral_implicit) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) +trees_per_dimension = (16, 8) +mesh = P4estMesh(trees_per_dimension; polydeg = polydeg, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = (true, false), initial_refinement_level = 0) + +boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall_2, + :y_pos => boundary_condition_slip_wall_2) + +initial_condition = initial_condition_warm_bubble + +semi = SemidiscretizationHyperbolicSplit(mesh, + (equations, equations), + initial_condition, + (solver_implicit, solver_explicit); + boundary_conditions = (boundary_conditions, + boundary_conditions), + source_terms = (nothing, source_terms_gravity),) +tspan = (0.0, 1000.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_callback = AnalysisCallback(semi, interval = 1000) + +alive_callback = AliveCallback(analysis_interval = 1000) + +save_solution = SaveSolutionCallback(interval = 1000, solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, alive_callback) + +############################################################################### +# run the simulation + +# Tolerances for GMRES residual, see https://jso.dev/Krylov.jl/stable/solvers/unsymmetric/#Krylov.gmres +atol_lin_solve = 1e-5 +rtol_lin_solve = 1e-5 + +# Jacobian-free Newton-Krylov (GMRES) solver, see +# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov +linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve) + +# Use second-order implicit Runge-Kutta BDF, see +# https://docs.sciml.ai/OrdinaryDiffEq/stable/imex/IMEXBDF/ +ode_alg = SBDF2(autodiff = AutoFiniteDiff(), linsolve = linsolve) + +atol_ode_solve = 1e-4 +rtol_ode_solve = 1e-4 +sol = solve(ode, ode_alg; + dt = 0.5, + abstol = atol_ode_solve, reltol = rtol_ode_solve, + ode_default_options()..., callback = callbacks); diff --git a/src/Trixi.jl b/src/Trixi.jl index 21a73867718..db081d8266c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -148,6 +148,7 @@ include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("semidiscretization/semidiscretization_coupled.jl") +include("semidiscretization/semidiscretization_split.jl") include("time_integration/time_integration.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") @@ -285,6 +286,8 @@ export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integ export SemidiscretizationHyperbolicParabolic +export SemidiscretizationHyperbolicSplit + export SemidiscretizationEulerAcoustics export SemidiscretizationEulerGravity, ParametersEulerGravity, diff --git a/src/semidiscretization/semidiscretization_split.jl b/src/semidiscretization/semidiscretization_split.jl new file mode 100644 index 00000000000..b8f2cbfb0b8 --- /dev/null +++ b/src/semidiscretization/semidiscretization_split.jl @@ -0,0 +1,236 @@ +# 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 +#! format: noindent + +""" + SemidiscretizationHyperbolicSplit + +A struct containing everything needed to describe a spatial semidiscretization +of a split-rhs corresponding to a hyperbolic conservation/balance law. +""" +struct SemidiscretizationHyperbolicSplit{Mesh, EquationsStiff, EquationsNonStiff, + InitialCondition, + BoundaryConditionsStiff, + BoundaryConditionsNonStiff, + SourceTermsStiff, SourceTermsNonStiff, + SolverStiff, SolverNonStiff, + CacheStiff, CacheNonStiff} <: + AbstractSemidiscretization + mesh::Mesh + + equations_stiff::EquationsStiff + equations_nonstiff::EquationsNonStiff + + initial_condition::InitialCondition + + boundary_conditions_stiff::BoundaryConditionsStiff + boundary_conditions_nonstiff::BoundaryConditionsNonStiff + + source_terms_stiff::SourceTermsStiff + source_terms_nonstiff::SourceTermsNonStiff + + solver_stiff::SolverStiff + solver_nonstiff::SolverNonStiff + + cache_stiff::CacheStiff + cache_nonstiff::CacheNonStiff + + performance_counter::PerformanceCounterList{2} +end + +""" + SemidiscretizationHyperbolicSplit(mesh, equations::Tuple, + initial_condition, + solver_stiff, solver_nonstiff; + source_terms=(nothing, nothing), + boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), + RealT=real(solver), uEltype=RealT, + initial_caches=(NamedTuple(), NamedTuple())) + +Construct a semidiscretization of a hyperbolic-split PDE. +""" +function SemidiscretizationHyperbolicSplit(mesh, equations::Tuple, + initial_condition, solvers::Tuple; + source_terms = (nothing, nothing), + boundary_conditions = (boundary_condition_periodic, + boundary_condition_periodic), + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT = real(first(solvers)), + uEltype = RealT, + initial_caches = (NamedTuple(), + NamedTuple())) + solver_stiff, solver_nonstiff = solvers + equations_stiff, equations_nonstiff = equations + @assert ndims(mesh) == ndims(equations_stiff) + @assert ndims(mesh) == ndims(equations_nonstiff) + @assert nvariables(equations_stiff) == nvariables(equations_nonstiff) + + boundary_conditions_stiff, boundary_conditions_nonstiff = boundary_conditions + initial_cache_stiff, initial_cache_nonstiff = initial_caches + source_terms_stiff, source_terms_nonstiff = source_terms + + cache_stiff = (; + create_cache(mesh, equations_stiff, solver_stiff, RealT, uEltype)..., + initial_cache_stiff...) + cache_nonstiff = (; + create_cache(mesh, equations_nonstiff, solver_nonstiff, RealT, + uEltype)..., + initial_cache_nonstiff...) + + _boundary_conditions_stiff = digest_boundary_conditions(boundary_conditions_stiff, + mesh, + solver_stiff, + cache_stiff) + _boundary_conditions_nonstiff = digest_boundary_conditions(boundary_conditions_nonstiff, + mesh, solver_nonstiff, + cache_nonstiff) + + check_periodicity_mesh_boundary_conditions(mesh, _boundary_conditions_stiff) + + performance_counter = PerformanceCounterList{2}(false) + + SemidiscretizationHyperbolicSplit{typeof(mesh), + typeof(equations_stiff), + typeof(equations_nonstiff), + typeof(initial_condition), + typeof(_boundary_conditions_stiff), + typeof(_boundary_conditions_nonstiff), + typeof(source_terms_stiff), + typeof(source_terms_nonstiff), + typeof(solver_stiff), typeof(solver_nonstiff), + typeof(cache_stiff), + typeof(cache_nonstiff)}(mesh, equations_stiff, + equations_nonstiff, + initial_condition, + _boundary_conditions_stiff, + _boundary_conditions_nonstiff, + source_terms_stiff, + source_terms_nonstiff, + solver_stiff, + solver_nonstiff, + cache_stiff, + cache_nonstiff, + performance_counter) +end + +function Base.show(io::IO, ::MIME"text/plain", + semi::SemidiscretizationHyperbolicSplit) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationHyperbolicSplit") + summary_line(io, "#spatial dimensions", ndims(semi.equations_stiff)) + summary_line(io, "mesh", semi.mesh) + summary_line(io, "hyperbolic equations stiff", + semi.equations_stiff |> typeof |> nameof) + summary_line(io, "hyperbolic equations nonstiff", + semi.equations_nonstiff |> typeof |> nameof) + summary_line(io, "initial condition", semi.initial_condition) + + summary_line(io, "source terms stiff", semi.source_terms_stiff) + summary_line(io, "source terms nonstiff", semi.source_terms_nonstiff) + summary_line(io, "solver stiff", semi.solver_stiff |> typeof |> nameof) + summary_line(io, "solver nonstiff", semi.solver_nonstiff |> typeof |> nameof) + summary_line(io, "total #DOFs per field", ndofs(semi)) + summary_footer(io) + end +end + +@inline Base.ndims(semi::SemidiscretizationHyperbolicSplit) = ndims(semi.mesh) + +# retain dispatch on hyperbolic non-stiff equations only +@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicSplit) + @unpack mesh, equations_nonstiff, solver_nonstiff, cache_nonstiff = semi + return mesh, equations_nonstiff, solver_nonstiff, cache_nonstiff +end + +function compute_coefficients(t, semi::SemidiscretizationHyperbolicSplit) + # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` + compute_coefficients(semi.initial_condition, t, semi) +end + +""" + semidiscretize(semi::SemidiscretizationHyperbolicSplit, tspan) + +Wrap the semidiscretization `semi` as a split ODE problem in the time interval `tspan` +that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). +The stiff hyperbolic right-hand side is the first function of the split ODE problem and +will be used by default by the implicit part of IMEX methods from the +SciML ecosystem. +""" +function Trixi.semidiscretize(semi::SemidiscretizationHyperbolicSplit, tspan; + reset_threads = true) + # Optionally reset Polyester.jl threads. See + # https://github.com/trixi-framework/Trixi.jl/issues/1583 + # https://github.com/JuliaSIMD/Polyester.jl/issues/30 + if reset_threads + Trixi.Polyester.reset_threads!() + end + u0_ode = compute_coefficients(first(tspan), semi) + # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using + # mpi_isparallel() && MPI.Barrier(mpi_comm()) + # See https://github.com/trixi-framework/Trixi.jl/issues/328 + iip = true # is-inplace, i.e., we modify a vector when calling rhs_parabolic!, rhs! + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer function first. + return SplitODEProblem{iip}(rhs_stiff!, rhs!, u0_ode, tspan, semi) +end + +function rhs_stiff!(du_ode, u_ode, semi::SemidiscretizationHyperbolicSplit, t) + @unpack mesh, equations_stiff, initial_condition, boundary_conditions_stiff, source_terms_stiff, solver_stiff, cache_stiff = semi + + u = wrap_array(u_ode, mesh, equations_stiff, solver_stiff, cache_stiff) + du = wrap_array(du_ode, mesh, equations_stiff, solver_stiff, cache_stiff) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "rhs! stiff" rhs!(du, u, t, mesh, equations_stiff, + boundary_conditions_stiff, + source_terms_stiff, + solver_stiff, + cache_stiff) + runtime = time_ns() - time_start + put!(semi.performance_counter.counters[1], runtime) + + return nothing +end + +# nonstiff `rhs!` +function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicSplit, t) + @unpack mesh, equations_nonstiff, initial_condition, boundary_conditions_nonstiff, source_terms_nonstiff, solver_nonstiff, cache_nonstiff = semi + + u = wrap_array(u_ode, mesh, equations_nonstiff, solver_nonstiff, cache_nonstiff) + du = wrap_array(du_ode, mesh, equations_nonstiff, solver_nonstiff, cache_nonstiff) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "rhs! nonstiff" rhs!(du, u, t, mesh, equations_nonstiff, + boundary_conditions_nonstiff, + source_terms_nonstiff, + solver_nonstiff, + cache_nonstiff) + runtime = time_ns() - time_start + put!(semi.performance_counter.counters[2], runtime) + + return nothing +end + +# Here we only pass the nonstiff solver and cache to the calc_error_norms function, +# since they are needed only for auxiliary functions, such as get_node_vars, etc. +function calc_error_norms(func, u_ode, t, analyzer, + semi::SemidiscretizationHyperbolicSplit, + cache_analysis) + @unpack mesh, equations_nonstiff, initial_condition, solver_nonstiff, cache_nonstiff = semi + u = wrap_array(u_ode, mesh, equations_nonstiff, solver_nonstiff, cache_nonstiff) + + calc_error_norms(func, u, t, analyzer, mesh, equations_nonstiff, initial_condition, + solver_nonstiff, cache_nonstiff, cache_analysis) +end +end # @muladd diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ca6af4c24ff..7d75f3dca47 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -891,6 +891,26 @@ end # (e.g., from type instabilities) @test_allocations(Trixi.rhs!, semi, sol, 1000) end + +@trixi_testset "elixir_euler_imex_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_imex_warm_bubble.jl"), + l2=[ + 2.5906685915145214e-5, + 0.001327868084047019, + 0.0013043939534472446, + 2.1316157151842505 + ], + linf=[ + 0.0003013316385283016, + 0.013751491038817676, + 0.011904250984857088, + 19.010794903791975 + ], + tspan=(0.0, 0.1)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_unit.jl b/test/test_unit.jl index 7937d0b7fce..bca2bb97c26 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -2846,6 +2846,42 @@ end # the sparsity pattern of the parabolic part of a hyperbolic-parabolic problem always includes the hyperbolic one @test jac_prototype_parabolic == jac_prototype_hyperbolic_parabolic end -end +@testset "ndims function for SemidiscretizaionHyperbolicSplit" begin + gamma = 1004 / 717 + equations = CompressibleEulerEquations2D(gamma) + + polydeg = 2 + basis = LobattoLegendreBasis(polydeg) + + volume_integral_explicit = VolumeIntegralFluxDifferencing(flux_ranocha) + solver_explicit = DGSEM(basis, flux_ranocha, volume_integral_explicit) + + volume_integral_implicit = VolumeIntegralFluxDifferencing(flux_ranocha) + solver_implicit = DGSEM(basis, flux_ranocha, volume_integral_implicit) + + coordinates_min = (0.0, 0.0) + coordinates_max = (20_000.0, 10_000.0) + trees_per_dimension = (16, 8) + mesh = P4estMesh(trees_per_dimension; polydeg = polydeg, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + periodicity = (true, false), initial_refinement_level = 0) + + boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall) + + initial_condition = initial_condition_convergence_test + + semi = SemidiscretizationHyperbolicSplit(mesh, + (equations, equations), + initial_condition, + (solver_implicit, solver_explicit); + boundary_conditions = (boundary_conditions, + boundary_conditions), + source_terms = (nothing, nothing),) + + @test Trixi.ndims(semi) == 2 +end +end end #module