diff --git a/NEWS.md b/NEWS.md index 5dd911391a6..265979c3508 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ for human readability. ## Changes when updating to v0.6 from v0.5.x #### Added +- AMR for hyperbolic-parabolic equations on 2D `P4estMesh` #### Changed diff --git a/Project.toml b/Project.toml index 17c02a2adac..4e2a89dc133 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.6.2-pre" +version = "0.6.3-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" @@ -66,7 +66,7 @@ Makie = "0.19" MuladdMacro = "0.2.2" Octavian = "0.3.5" OffsetArrays = "1.3" -P4est = "0.4" +P4est = "0.4.9" Polyester = "0.7.5" PrecompileTools = "1.1" Printf = "1" @@ -84,7 +84,7 @@ StaticArrays = "1" StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" -T8code = "0.4.1" +T8code = "0.4.3" TimerOutputs = "0.5" Triangulate = "2.0" TriplotBase = "0.1" diff --git a/docs/src/meshes/dgmulti_mesh.md b/docs/src/meshes/dgmulti_mesh.md index fc086bba146..efa71334bf8 100644 --- a/docs/src/meshes/dgmulti_mesh.md +++ b/docs/src/meshes/dgmulti_mesh.md @@ -21,7 +21,7 @@ around Jonathan Shewchuk's [Triangle](https://www.cs.cmu.edu/~quake/triangle.htm ## The `DGMulti` solver type -Trixi.jl solvers on simplicial meshes use the `[DGMulti](@ref)` solver type, which allows users to specify +Trixi.jl solvers on simplicial meshes use the [`DGMulti`](@ref) solver type, which allows users to specify `element_type` and `approximation_type` in addition to `polydeg`, `surface_flux`, `surface_integral`, and `volume_integral`. diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md index fa6fc1a5d32..f599eb5fafe 100644 --- a/docs/src/parallelization.md +++ b/docs/src/parallelization.md @@ -53,30 +53,43 @@ a system-provided MPI installation with Trixi.jl can be found in the following s ### [Using a system-provided MPI installation](@id parallel_system_MPI) -When using Trixi.jl with a system-provided MPI backend the underlying -[`p4est`](https://github.com/cburstedde/p4est) and [`t8code`](https://github.com/DLR-AMR/t8code) -libraries need to be compiled with the same MPI installation. Therefore, you also need to -use system-provided `p4est` and `t8code` installations (for notes on how to install `p4est` -and `t8code` see e.g. [here](https://github.com/cburstedde/p4est/blob/master/README) and -[here](https://github.com/DLR-AMR/t8code/wiki/Installation), use the configure option -`--enable-mpi`). Note that `t8code` already comes with a `p4est` installation, so it suffices -to install `t8code`. In addition, [P4est.jl](https://github.com/trixi-framework/P4est.jl) and -[T8code.jl](https://github.com/DLR-AMR/T8code.jl) need to be configured to use the custom -installations. Follow the steps described -[here](https://github.com/DLR-AMR/T8code.jl/blob/main/README.md#installation) and -[here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation) for the -configuration. The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be -the same for P4est.jl and T8code.jl. This could e.g. be `libp4est.so` that usually can be found -in `lib/` or `local/lib/` in the installation directory of `t8code`. -In total, in your active Julia project you should have a LocalPreferences.toml file with sections -`[MPIPreferences]`, `[T8code]` and `[P4est]` as well as an entry `MPIPreferences` in your -Project.toml to use a custom MPI installation. A `LocalPreferences.toml` file +When using Trixi.jl with a system-provided MPI backend, the underlying +[`p4est`](https://github.com/cburstedde/p4est), [`t8code`](https://github.com/DLR-AMR/t8code) +and [`HDF5`](https://github.com/HDFGroup/hdf5) libraries need to be compiled with the same MPI +installation. If you want to use `p4est` (via the `P4estMesh`) or `t8code` (via the `T8codeMesh`) +from Trixi.jl, you also need to use system-provided `p4est` or `t8code` installations +(for notes on how to install `p4est` and `t8code` see, e.g., [here](https://github.com/cburstedde/p4est/blob/master/README) +and [here](https://github.com/DLR-AMR/t8code/wiki/Installation), use the configure option +`--enable-mpi`). Otherwise, there will be warnings that no preference is set for P4est.jl and +T8code.jl that can be ignored if you do not use these libraries from Trixi.jl. Note that +`t8code` already comes with a `p4est` installation, so it suffices to install `t8code`. +In order to use system-provided `p4est` and `t8code` installations, [P4est.jl](https://github.com/trixi-framework/P4est.jl) +and [T8code.jl](https://github.com/DLR-AMR/T8code.jl) need to be configured to use the custom +installations. Follow the steps described [here](https://github.com/DLR-AMR/T8code.jl/blob/main/README.md#installation) and +[here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation) for the configuration. +The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be +the same for P4est.jl and T8code.jl. This could, e.g., be `libp4est.so` that usually can be found +in `lib/` or `local/lib/` in the installation directory of `t8code`. Note that the `T8codeMesh`, however, +does not support MPI yet. +The preferences for [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) always need to be set, even if you +do not want to use `HDF5` from Trixi.jl, see also [issue #1079 in HDF5.jl](https://github.com/JuliaIO/HDF5.jl/issues/1079). +To set the preferences for HDF5.jl, follow the instructions described +[here](https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#Using-parallel-input-and-output). + +In total, in your active Julia project you should have a `LocalPreferences.toml` file with sections +`[MPIPreferences]`, `[T8code]` (only needed if `T8codeMesh` is used), `[P4est]` (only needed if +`P4estMesh` is used), and `[HDF5]` as well as an entry `MPIPreferences` in your +`Project.toml` to use a custom MPI installation. A `LocalPreferences.toml` file created as described above might look something like the following: ```toml [HDF5] libhdf5 = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so" libhdf5_hl = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so" +[HDF5_jll] +libhdf5_hl_path = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so" +libhdf5_path = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so" + [MPIPreferences] __clear__ = ["preloads_env_switch"] _format = "1.0" @@ -97,6 +110,22 @@ libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so" libt8 = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libt8.so" ``` +This file is created with the following sequence of commands: +```julia +julia> using MPIPreferences +julia> MPIPreferences.use_system_binary() +``` +Restart the Julia REPL +```julia +julia> using P4est +julia> P4est.set_library_p4est!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so") +julia> P4est.set_library_sc!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so") +julia> using T8code +julia> T8code.set_libraries_path!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/") +julia> using HDF5 +julia> HDF5.API.set_libraries!("/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so", "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so") +``` +After the preferences are set, restart the Julia REPL again. ### [Usage](@id parallel_usage) @@ -218,7 +247,7 @@ julia> HDF5.API.set_libraries!("/path/to/your/libhdf5.so", "/path/to/your/libhdf ``` For more information see also the [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). In total, you should -have a file called LocalPreferences.toml in the project directory that contains a section +have a file called `LocalPreferences.toml` in the project directory that contains a section `[MPIPreferences]`, a section `[HDF5]` with entries `libhdf5` and `libhdf5_hl`, a section `[P4est]` with the entry `libp4est` as well as a section `[T8code]` with the entries `libt8`, `libp4est` and `libsc`. diff --git a/docs/src/performance.md b/docs/src/performance.md index bbe3a3390b7..df66f451b79 100644 --- a/docs/src/performance.md +++ b/docs/src/performance.md @@ -34,7 +34,7 @@ Hence, you should at least investigate the performance roughly by comparing the timings of several elixirs. Deeper investigations and micro-benchmarks should usually use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl). For example, the following steps were used to benchmark the changes introduced in -https://github.com/trixi-framework/Trixi.jl/pull/256. +[PR #256](https://github.com/trixi-framework/Trixi.jl/pull/256). 1. `git checkout e7ebf3846b3fd62ee1d0042e130afb50d7fe8e48` (new version) 2. Start `julia --threads=1 --check-bounds=no`. diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl new file mode 100644 index 00000000000..f87c0e056ca --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -0,0 +1,98 @@ +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)) + +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) + +# 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)) + +# 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, 0.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 = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 0.9, + max_level = 3, max_threshold = 1.0) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 50) + +# 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, amr_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, dt = 1e-7, 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_advection_diffusion_periodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl new file mode 100644 index 00000000000..fdecb05f7bd --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.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)) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max) + +# 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 + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.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) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 1.25, + max_level = 3, max_threshold = 1.45) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 20) + +# 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, amr_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_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index e6566edb18a..bc28ae6ffb3 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -41,7 +41,6 @@ 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, diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl new file mode 100644 index 00000000000..898366969a4 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl @@ -0,0 +1,94 @@ +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 = (6, 6) +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) + +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 = 2000) +analysis_interval = 2000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.005, + max_level = 2, max_threshold = 0.01) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 50, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, amr_callback) +# 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/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index c855f448fd4..1b268422e1f 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -7,11 +7,97 @@ using Trixi gamma = 1.4 equations = CompressibleEulerEquations2D(gamma) -initial_condition = Trixi.initial_condition_double_mach_reflection +""" + initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D) + +Compressible Euler setup for a double Mach reflection problem. +Involves strong shock interactions as well as steady / unsteady flow structures. +Also exercises special boundary conditions along the bottom of the domain that is a mixture of +Dirichlet and slip wall. +See Section IV c on the paper below for details. + +- Paul Woodward and Phillip Colella (1984) + The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks. + [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6) +""" +@inline function initial_condition_double_mach_reflection(x, t, + equations::CompressibleEulerEquations2D) + if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3) + phi = pi / 6 + sin_phi, cos_phi = sincos(phi) + + rho = 8.0 + v1 = 8.25 * cos_phi + v2 = -8.25 * sin_phi + p = 116.5 + else + rho = 1.4 + v1 = 0.0 + v2 = 0.0 + p = 1.0 + end + + prim = SVector(rho, v1, v2, p) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_double_mach_reflection boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition) -boundary_conditions = (y_neg = Trixi.boundary_condition_mixed_dirichlet_wall, +# Special mixed boundary condition type for the :y_neg of the domain. +# It is charachteristic-based when x < 1/6 and a slip wall when x >= 1/6 +# Note: Only for StructuredMesh +@inline function boundary_condition_mixed_characteristic_wall(u_inner, + normal_direction::AbstractVector, + direction, + x, t, surface_flux_function, + equations::CompressibleEulerEquations2D) + if x[1] < 1 / 6 + # From the BoundaryConditionCharacteristic + # get the external state of the solution + u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, + u_inner, + normal_direction, + direction, x, t, + equations) + # Calculate boundary flux + flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + else # x[1] >= 1 / 6 + # Use the free slip wall BC otherwise + flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations) + end + + return flux +end + +# Note: Only for StructuredMesh +@inline function Trixi.get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), + normal_direction::AbstractVector, direction, + equations, + dg, indices...) + x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) + if x[1] < 1 / 6 # BoundaryConditionCharacteristic + u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, + u_inner, + normal_direction, + direction, x, t, + equations) + + else # if x[1] >= 1 / 6 # boundary_condition_slip_wall + u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations) + + u_outer = SVector(u_inner[1], + u_inner[2] - 2.0 * u_rotate[2], + u_inner[3] - 2.0 * u_rotate[3], + u_inner[4]) + end + return u_outer +end + +boundary_conditions = (y_neg = boundary_condition_mixed_characteristic_wall, y_pos = boundary_condition_inflow_outflow, x_pos = boundary_condition_inflow_outflow, x_neg = boundary_condition_inflow_outflow) diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl index 6a04b39ec4f..ad537b49684 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -7,11 +7,97 @@ using Trixi gamma = 1.4 equations = CompressibleEulerEquations2D(gamma) -initial_condition = Trixi.initial_condition_double_mach_reflection +""" + initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D) + +Compressible Euler setup for a double Mach reflection problem. +Involves strong shock interactions as well as steady / unsteady flow structures. +Also exercises special boundary conditions along the bottom of the domain that is a mixture of +Dirichlet and slip wall. +See Section IV c on the paper below for details. + +- Paul Woodward and Phillip Colella (1984) + The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks. + [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6) +""" +@inline function initial_condition_double_mach_reflection(x, t, + equations::CompressibleEulerEquations2D) + if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3) + phi = pi / 6 + sin_phi, cos_phi = sincos(phi) + + rho = 8.0 + v1 = 8.25 * cos_phi + v2 = -8.25 * sin_phi + p = 116.5 + else + rho = 1.4 + v1 = 0.0 + v2 = 0.0 + p = 1.0 + end + + prim = SVector(rho, v1, v2, p) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_double_mach_reflection boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition) -boundary_conditions = (y_neg = Trixi.boundary_condition_mixed_dirichlet_wall, +# Special mixed boundary condition type for the :y_neg of the domain. +# It is charachteristic-based when x < 1/6 and a slip wall when x >= 1/6 +# Note: Only for StructuredMesh +@inline function boundary_condition_mixed_characteristic_wall(u_inner, + normal_direction::AbstractVector, + direction, + x, t, surface_flux_function, + equations::CompressibleEulerEquations2D) + if x[1] < 1 / 6 + # From the BoundaryConditionCharacteristic + # get the external state of the solution + u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, + u_inner, + normal_direction, + direction, x, t, + equations) + # Calculate boundary flux + flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + else # x[1] >= 1 / 6 + # Use the free slip wall BC otherwise + flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations) + end + + return flux +end + +# Note: Only for StructuredMesh +@inline function Trixi.get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), + normal_direction::AbstractVector, direction, + equations, + dg, indices...) + x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) + if x[1] < 1 / 6 # BoundaryConditionCharacteristic + u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, + u_inner, + normal_direction, + direction, x, t, + equations) + + else # if x[1] >= 1 / 6 # boundary_condition_slip_wall + u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations) + + u_outer = SVector(u_inner[1], + u_inner[2] - 2.0 * u_rotate[2], + u_inner[3] - 2.0 * u_rotate[3], + u_inner[4]) + end + return u_outer +end + +boundary_conditions = (y_neg = boundary_condition_mixed_characteristic_wall, y_pos = boundary_condition_inflow_outflow, x_pos = boundary_condition_inflow_outflow, x_neg = boundary_condition_inflow_outflow) diff --git a/examples/structured_2d_dgsem/elixir_euler_shock_upstream_MCL.jl b/examples/structured_2d_dgsem/elixir_euler_shock_upstream_MCL.jl index cc12990ab10..2f2bc26f219 100644 --- a/examples/structured_2d_dgsem/elixir_euler_shock_upstream_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_shock_upstream_MCL.jl @@ -101,7 +101,7 @@ mapping_as_string = "a = sqrt(5.9^2 - 3.85^2); alpha = acos(3.85 / 5.9); l = (pi "f4(s) = SVector(0.0, -(a - 1.0) * 0.5 * (s + 1.0) + a); " * "faces = (f1, f2, f3, f4); mapping = Trixi.transfinite_mapping(faces)" -cells_per_dimension = (24, 36) +cells_per_dimension = (8, 12) mesh = StructuredMesh(cells_per_dimension, mapping_bow, mapping_as_string = mapping_as_string, periodicity = false) @@ -112,7 +112,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 10.0) +tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl b/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl new file mode 100644 index 00000000000..95f71f38130 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl @@ -0,0 +1,55 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler multicomponent equations +equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, 0.4)) + +initial_condition = initial_condition_convergence_test + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +cells_per_dimension = (16, 16) +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index 7bd1ec0c647..70d76fc9075 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -40,7 +40,6 @@ 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, diff --git a/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl index 968af339cbd..0b826254129 100644 --- a/src/auxiliary/p4est.jl +++ b/src/auxiliary/p4est.jl @@ -13,14 +13,20 @@ This function will check if `p4est` is already initialized and if yes, do nothing, thus it is safe to call it multiple times. """ function init_p4est() - p4est_package_id = P4est.package_id() - if p4est_package_id >= 0 - return nothing + # Only initialize p4est if P4est.jl can be used + if P4est.preferences_set_correctly() + p4est_package_id = P4est.package_id() + if p4est_package_id >= 0 + return nothing + end + + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations + p4est_init(C_NULL, SC_LP_ERROR) + else + @warn "Preferences for P4est.jl are not set correctly. Until fixed, using `P4estMesh` will result in a crash. " * + "See also https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#parallel_system_MPI" end - # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations - p4est_init(C_NULL, SC_LP_ERROR) - return nothing end diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index 37cb782bb93..bd781b21c1e 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -7,34 +7,40 @@ is already initialized and if yes, do nothing, thus it is safe to call it multiple times. """ function init_t8code() - t8code_package_id = t8_get_package_id() - if t8code_package_id >= 0 - return nothing - end + # Only initialize t8code if T8code.jl can be used + if T8code.preferences_set_correctly() + t8code_package_id = t8_get_package_id() + if t8code_package_id >= 0 + return nothing + end - # Initialize the sc library, has to happen before we initialize t8code. - let catch_signals = 0, print_backtrace = 0, log_handler = C_NULL - T8code.Libt8.sc_init(mpi_comm(), catch_signals, print_backtrace, log_handler, - T8code.Libt8.SC_LP_ERROR) - end + # Initialize the sc library, has to happen before we initialize t8code. + let catch_signals = 0, print_backtrace = 0, log_handler = C_NULL + T8code.Libt8.sc_init(mpi_comm(), catch_signals, print_backtrace, log_handler, + T8code.Libt8.SC_LP_ERROR) + end - if T8code.Libt8.p4est_is_initialized() == 0 - # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations - T8code.Libt8.p4est_init(C_NULL, T8code.Libt8.SC_LP_ERROR) - end + if T8code.Libt8.p4est_is_initialized() == 0 + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations + T8code.Libt8.p4est_init(C_NULL, T8code.Libt8.SC_LP_ERROR) + end - # Initialize t8code with log level ERROR to prevent a lot of output in AMR simulations. - t8_init(T8code.Libt8.SC_LP_ERROR) - - if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") - # Normally, `sc_finalize` should always be called during shutdown of an - # application. It checks whether there is still un-freed memory by t8code - # and/or T8code.jl and throws an exception if this is the case. For - # production runs this is not mandatory, but is helpful during - # development. Hence, this option is only activated when environment - # variable TRIXI_T8CODE_SC_FINALIZE exists. - @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." - MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) + # Initialize t8code with log level ERROR to prevent a lot of output in AMR simulations. + t8_init(T8code.Libt8.SC_LP_ERROR) + + if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") + # Normally, `sc_finalize` should always be called during shutdown of an + # application. It checks whether there is still un-freed memory by t8code + # and/or T8code.jl and throws an exception if this is the case. For + # production runs this is not mandatory, but is helpful during + # development. Hence, this option is only activated when environment + # variable TRIXI_T8CODE_SC_FINALIZE exists. + @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." + MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) + end + else + @warn "Preferences for T8code.jl are not set correctly. Until fixed, using `T8codeMesh` will result in a crash. " * + "See also https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#parallel_system_MPI" end return nothing diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index ba840ff9675..5854c8617c3 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -528,6 +528,103 @@ function copy_to_quad_iter_volume(info, user_data) return nothing end +# specialized callback which includes the `cache_parabolic` argument +function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, + equations, dg::DG, cache, cache_parabolic, + semi, + t, iter; + only_refine = false, only_coarsen = false, + passive_args = ()) + @unpack controller, adaptor = amr_callback + + u = wrap_array(u_ode, mesh, equations, dg, cache) + lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache, + t = t, iter = iter) + + @boundscheck begin + @assert axes(lambda)==(Base.OneTo(ncells(mesh)),) ("Indicator array (axes = $(axes(lambda))) and mesh cells (axes = $(Base.OneTo(ncells(mesh)))) have different axes") + end + + # Copy controller value of each quad to the quad's user data storage + iter_volume_c = cfunction(copy_to_quad_iter_volume, Val(ndims(mesh))) + + # The pointer to lambda will be interpreted as Ptr{Int} below + @assert lambda isa Vector{Int} + iterate_p4est(mesh.p4est, lambda; iter_volume_c = iter_volume_c) + + @trixi_timeit timer() "refine" if !only_coarsen + # Refine mesh + refined_original_cells = @trixi_timeit timer() "mesh" refine!(mesh) + + # Refine solver + @trixi_timeit timer() "solver" refine!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + refined_original_cells) + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" refine!(p_u_ode, adaptor, p_mesh, + p_equations, + p_dg, p_cache, + refined_original_cells) + end + else + # If there is nothing to refine, create empty array for later use + refined_original_cells = Int[] + end + + @trixi_timeit timer() "coarsen" if !only_refine + # Coarsen mesh + coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh) + + # coarsen solver + @trixi_timeit timer() "solver" coarsen!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + coarsened_original_cells) + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" coarsen!(p_u_ode, adaptor, p_mesh, + p_equations, + p_dg, p_cache, + coarsened_original_cells) + end + else + # If there is nothing to coarsen, create empty array for later use + coarsened_original_cells = Int[] + end + + # Store whether there were any cells coarsened or refined and perform load balancing + has_changed = !isempty(refined_original_cells) || !isempty(coarsened_original_cells) + # Check if mesh changed on other processes + if mpi_isparallel() + has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[] + end + + if has_changed # TODO: Taal decide, where shall we set this? + # don't set it to has_changed since there can be changes from earlier calls + mesh.unsaved_changes = true + + if mpi_isparallel() && amr_callback.dynamic_load_balancing + @trixi_timeit timer() "dynamic load balancing" begin + global_first_quadrant = unsafe_wrap(Array, + unsafe_load(mesh.p4est).global_first_quadrant, + mpi_nranks() + 1) + old_global_first_quadrant = copy(global_first_quadrant) + partition!(mesh) + rebalance_solver!(u_ode, mesh, equations, dg, cache, + old_global_first_quadrant) + end + end + + reinitialize_boundaries!(semi.boundary_conditions, cache) + # if the semidiscretization also stores parabolic boundary conditions, + # reinitialize them after each refinement step as well. + if hasproperty(semi, :boundary_conditions_parabolic) + reinitialize_boundaries!(semi.boundary_conditions_parabolic, cache) + end + end + + # Return true if there were any cells coarsened or refined, otherwise false + return has_changed +end + # 2D function cfunction(::typeof(copy_to_quad_iter_volume), ::Val{2}) @cfunction(copy_to_quad_iter_volume, Cvoid, diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 6395a9f348f..969f9c564f3 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -136,8 +136,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM return nothing end -# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes -function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, +function refine!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_refine) # Call `refine!` for the hyperbolic part, which does the heavy lifting of @@ -298,8 +298,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, return nothing end -# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes -function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, +function coarsen!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_remove) # Call `coarsen!` for the hyperbolic part, which does the heavy lifting of diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index f1cec04d282..182ce73d253 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -496,71 +496,6 @@ end return cons end -""" - initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D) - -Compressible Euler setup for a double Mach reflection problem. -Involves strong shock interactions as well as steady / unsteady flow structures. -Also exercises special boundary conditions along the bottom of the domain that is a mixture of -Dirichlet and slip wall. -See Section IV c on the paper below for details. - -- Paul Woodward and Phillip Colella (1984) - The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks. - [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6) -""" -@inline function initial_condition_double_mach_reflection(x, t, - equations::CompressibleEulerEquations2D) - if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3) - phi = pi / 6 - sin_phi, cos_phi = sincos(phi) - - rho = 8.0 - v1 = 8.25 * cos_phi - v2 = -8.25 * sin_phi - p = 116.5 - else - rho = 1.4 - v1 = 0.0 - v2 = 0.0 - p = 1.0 - end - - prim = SVector(rho, v1, v2, p) - return prim2cons(prim, equations) -end - -# Special mixed boundary condition type for the :Bottom of the domain. -# It is charachteristic when x < 1/6 and a slip wall when x >= 1/6 -@inline function boundary_condition_mixed_dirichlet_wall(u_inner, - normal_direction::AbstractVector, - direction, - x, t, surface_flux_function, - equations::CompressibleEulerEquations2D) - # Note: Only for StructuredMesh - if x[1] < 1 / 6 - # # From the BoundaryConditionDirichlet - # # get the external value of the solution - # u_boundary = initial_condition_double_mach_reflection(x, t, equations) - - # From the BoundaryConditionCharacteristic - # get the external state of the solution - u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, - u_inner, - normal_direction, - direction, x, t, - equations) - # Calculate boundary flux - flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) - else # x[1] >= 1 / 6 - # Use the free slip wall BC otherwise - flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, - surface_flux_function, equations) - end - - return flux -end - # Calculate 2D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 7b437f4a1b4..ecd3bc80c0a 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -270,6 +270,29 @@ end return vcat(f_other, f_rho) end +# Calculate 1D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + gamma = totalgamma(u, equations) + p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + + f_rho = densities(u, v_normal, equations) + f1 = rho_v1 * v_normal + p * normal_direction[1] + f2 = rho_v2 * v_normal + p * normal_direction[2] + f3 = (rho_e + p) * v_normal + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + """ flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) @@ -446,6 +469,76 @@ See also return vcat(f_other, f_rho) end +@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, cv = equations + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # helpful variables + help1_ll = zero(v1_ll) + help1_rr = zero(v1_rr) + enth_ll = zero(v1_ll) + enth_rr = zero(v1_rr) + for i in eachcomponent(equations) + enth_ll += u_ll[i + 3] * gas_constants[i] + enth_rr += u_rr[i + 3] * gas_constants[i] + help1_ll += u_ll[i + 3] * cv[i] + help1_rr += u_rr[i + 3] * cv[i] + end + + # temperature and pressure + T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + (v_dot_n_ll + v_dot_n_rr) + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] + f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) @@ -491,6 +584,50 @@ end return (abs(v1) + c, abs(v2) + c) end +@inline function rotate_to_x(u, normal_vector, + equations::CompressibleEulerMulticomponentEquations2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ n_1 n_2 0 0; + # t_1 t_2 0 0; + # 0 0 1 0 + # 0 0 0 1] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] + s * u[2], + -s * u[1] + c * u[2], + u[3], + densities...) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function rotate_from_x(u, normal_vector, + equations::CompressibleEulerMulticomponentEquations2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ n_1 t_1 0 0; + # n_2 t_2 0 0; + # 0 0 1 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] - s * u[2], + s * u[1] + c * u[2], + u[3], + densities...) +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 71782644b17..43be04f745d 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -61,6 +61,23 @@ struct FluxRotated{NumericalFlux} numerical_flux::NumericalFlux end +# Rotated surface flux computation (2D version) +@inline function (flux_rotated::FluxRotated)(u, + normal_direction::AbstractVector, + equations::AbstractEquations{2}) + @unpack numerical_flux = flux_rotated + + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal_vector = normal_direction / norm_ + + u_rotated = rotate_to_x(u, normal_vector, equations) + + f = numerical_flux(u_rotated, 1, equations) + + return rotate_from_x(f, normal_vector, equations) * norm_ +end + # Rotated surface flux computation (2D version) @inline function (flux_rotated::FluxRotated)(u_ll, u_rr, normal_direction::AbstractVector, diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 0176f5c6346..5fe68e06710 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -429,13 +429,17 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) @unpack boundaries = cache resize!(boundaries, required.boundaries) - # resize mortars container - @unpack mortars = cache - resize!(mortars, required.mortars) - - # re-initialize containers together to reduce - # the number of iterations over the mesh in `p4est` - init_surfaces!(interfaces, mortars, boundaries, mesh) + # re-initialize mortars container + if hasproperty(cache, :mortars) # cache_parabolic does not carry mortars + @unpack mortars = cache + resize!(mortars, required.mortars) + + # re-initialize containers together to reduce + # the number of iterations over the mesh in `p4est` + init_surfaces!(interfaces, mortars, boundaries, mesh) + else + init_surfaces!(interfaces, nothing, boundaries, mesh) + end end # A helper struct used in initialization methods below diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index cf07645b949..9dd10df16ae 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -94,8 +94,19 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{2}, dg.surface_integral, dg) end - # TODO: parabolic; extend to mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars (specialized for AbstractEquationsParabolic) + # !!! NOTE: we reuse the hyperbolic cache here since it contains "mortars" and "u_threaded". See https://github.com/trixi-framework/Trixi.jl/issues/1674 for a discussion + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes (specialized for AbstractEquationsParabolic) + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -176,14 +187,15 @@ function calc_gradient!(gradients, u_transformed, t, end end - # Prolong solution to interfaces + # Prolong solution to interfaces. + # This reuses `prolong2interfaces` for the purely hyperbolic case. @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) end - # Calculate interface fluxes for the gradient. This reuses P4est `calc_interface_flux!` along with a - # specialization for AbstractEquationsParabolic. + # Calculate interface fluxes for the gradient. + # This reuses `calc_interface_flux!` for the purely hyperbolic case. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, False(), # False() = no nonconservative terms @@ -202,8 +214,21 @@ function calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, dg.surface_integral, dg) end - # TODO: parabolic; mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars. This resues the hyperbolic version of `prolong2mortars` + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes. This reuses the hyperbolic version of `calc_mortar_flux`, + # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for + # AbstractEquationsParabolic. + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, + mesh, False(), # False() = no nonconservative terms + equations_parabolic, + dg.mortar, dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -303,6 +328,64 @@ function calc_gradient!(gradients, u_transformed, t, return nothing end +# This version is called during `calc_gradients!` and must be specialized because the +# flux for the gradient is {u}, which doesn't depend on the outward normal. Thus, +# you don't need to scale by 2 (e.g., the scaling factor in the normals (and in the +# contravariant vectors) along large/small elements across a non-conforming +# interface in 2D) and flip the sign when storing the mortar fluxes back +# into `surface_flux_values`. +@inline function mortar_fluxes_to_elements!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + dg::DGSEM, cache, mortar, fstar, u_buffer) + @unpack neighbor_ids, node_indices = cache.mortars + # Copy solution small to small + small_indices = node_indices[1, mortar] + small_direction = indices2direction(small_indices) + + for position in 1:2 + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, small_direction, element] = fstar[position][v, + i] + end + end + end + + # Project small fluxes to large element. + multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_upper, fstar[2], + mortar_l2.reverse_lower, fstar[1]) + + # Copy interpolated flux values from buffer to large element face in the + # correct orientation. + # Note that the index of the small sides will always run forward but + # the index of the large side might need to run backwards for flipped sides. + large_element = neighbor_ids[3, mortar] + large_indices = node_indices[2, mortar] + large_direction = indices2direction(large_indices) + + if :i_backward in large_indices + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, end + 1 - i, large_direction, large_element] = u_buffer[v, + i] + end + end + else + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, large_direction, large_element] = u_buffer[v, + i] + end + end + end + + return nothing +end + # This version is used for parabolic gradient computations @inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::False, @@ -321,7 +404,7 @@ end flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux - # Note that we don't flip the sign on the secondondary flux. This is because for parabolic terms, + # Note that we don't flip the sign on the secondary flux. This is because for parabolic terms, # the normals are not embedded in `flux_` for the parabolic gradient computations. for v in eachvariable(equations) surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = flux_[v] @@ -520,6 +603,163 @@ function calc_interface_flux!(surface_flux_values, return nothing end +function prolong2mortars_divergence!(cache, flux_viscous::Vector{Array{uEltype, 4}}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) where {uEltype <: Real} + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + flux_viscous_x, flux_viscous_y = flux_viscous + + @threaded for mortar in eachmortar(dg, cache) + # Copy solution data from the small elements using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + small_indices = node_indices[1, mortar] + direction_index = indices2direction(small_indices) + + i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], + index_range) + j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], + index_range) + + for position in 1:2 + i_small = i_small_start + j_small = j_small_start + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, + contravariant_vectors, + i_small, j_small, element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_small, j_small, element], + flux_viscous_y[v, i_small, j_small, element]) + + cache.mortars.u[1, v, position, i, mortar] = dot(flux_viscous, + normal_direction) + end + i_small += i_small_step + j_small += j_small_step + end + end + + # Buffer to copy solution values of the large element in the correct orientation + # before interpolating + u_buffer = cache.u_threaded[Threads.threadid()] + + # Copy solution of large element face to buffer in the + # correct orientation + large_indices = node_indices[2, mortar] + direction_index = indices2direction(large_indices) + + i_large_start, i_large_step = index_to_start_step_2d(large_indices[1], + index_range) + j_large_start, j_large_step = index_to_start_step_2d(large_indices[2], + index_range) + + i_large = i_large_start + j_large = j_large_start + element = neighbor_ids[3, mortar] + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, contravariant_vectors, + i_large, j_large, element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, element], + flux_viscous_y[v, i_large, j_large, element]) + + # We prolong the viscous flux dotted with respect the outward normal + # on the small element. We scale by -1/2 here because the normal + # direction on the large element is negative 2x that of the small + # element (these normal directions are "scaled" by the surface Jacobian) + u_buffer[v, i] = -0.5 * dot(flux_viscous, normal_direction) + end + i_large += i_large_step + j_large += j_large_step + end + + # Interpolate large element face data from buffer to small face locations + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, mortar), + mortar_l2.forward_lower, + u_buffer) + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, mortar), + mortar_l2.forward_upper, + u_buffer) + end + + return nothing +end + +# We specialize `calc_mortar_flux!` for the divergence part of +# the parabolic terms. +function calc_mortar_flux_divergence!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + @unpack fstar_upper_threaded, fstar_lower_threaded = cache + index_range = eachnode(dg) + + @threaded for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar = (fstar_lower_threaded[Threads.threadid()], + fstar_upper_threaded[Threads.threadid()]) + + for position in 1:2 + for node in eachnode(dg) + for v in eachvariable(equations) + viscous_flux_normal_ll = cache.mortars.u[1, v, position, node, mortar] + viscous_flux_normal_rr = cache.mortars.u[2, v, position, node, mortar] + + # TODO: parabolic; only BR1 at the moment + fstar[position][v, node] = 0.5 * (viscous_flux_normal_ll + + viscous_flux_normal_rr) + end + end + end + + # Buffer to interpolate flux values of the large element to before + # copying in the correct orientation + u_buffer = cache.u_threaded[Threads.threadid()] + + # this reuses the hyperbolic version of `mortar_fluxes_to_elements!` + mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar, u_buffer) + end + + return nothing +end + +# We structure `calc_interface_flux!` similarly to "calc_mortar_flux!" for +# hyperbolic equations with no nonconservative terms. +# The reasoning is that parabolic fluxes are treated like conservative +# terms (e.g., we compute a viscous conservative "flux") and thus no +# non-conservative terms are present. +@inline function calc_mortar_flux!(fstar, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + nonconservative_terms::False, + equations::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + mortar_index, position_index, normal_direction, + node_index) + @unpack u = cache.mortars + @unpack surface_flux = surface_integral + + u_ll, u_rr = get_surface_node_vars(u, equations, dg, position_index, node_index, + mortar_index) + + # TODO: parabolic; only BR1 at the moment + flux_ = 0.5 * (u_ll + u_rr) + + # Copy flux to buffer + set_node_vars!(fstar[position_index], flux_, equations, dg, node_index) +end + # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::P4estMesh{2}, diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 86b1da7c13a..7d997d8389a 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -85,7 +85,7 @@ function calc_volume_integral!(du, u, cache) if limiter.smoothness_indicator - @unpack element_ids_dg, element_ids_dgfv = cache + (; element_ids_dg, element_ids_dgfv) = cache # Calculate element-wise blending factors α alpha_element = @trixi_timeit timer() "element-wise blending factors" limiter.IndicatorHG(u, mesh, @@ -181,11 +181,11 @@ end nonconservative_terms::False, equations, volume_integral, limiter::SubcellLimiterMCL, dg::DGSEM, cache) - @unpack inverse_weights = dg.basis - @unpack volume_flux_dg, volume_flux_fv = volume_integral + (; inverse_weights) = dg.basis + (; volume_flux_dg, volume_flux_fv) = volume_integral # high-order DG fluxes - @unpack fhat1_L_threaded, fhat1_R_threaded, fhat2_L_threaded, fhat2_R_threaded = cache + (; fhat1_L_threaded, fhat1_R_threaded, fhat2_L_threaded, fhat2_R_threaded) = cache fhat1_L = fhat1_L_threaded[Threads.threadid()] fhat1_R = fhat1_R_threaded[Threads.threadid()] fhat2_L = fhat2_L_threaded[Threads.threadid()] @@ -195,7 +195,7 @@ end cache) # low-order FV fluxes - @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache + (; fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded) = cache fstar1_L = fstar1_L_threaded[Threads.threadid()] fstar2_L = fstar2_L_threaded[Threads.threadid()] fstar1_R = fstar1_R_threaded[Threads.threadid()] @@ -215,7 +215,7 @@ end limiter, dg, element, cache, fstar1_L, fstar2_L) - @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) du[v, i, j, element] += inverse_weights[i] * @@ -585,7 +585,7 @@ end u, mesh, nonconservative_terms::False, equations, limiter::SubcellLimiterMCL, dg, element, cache) - @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes for j in eachnode(dg), i in 2:nnodes(dg) for v in eachvariable(equations) @@ -622,7 +622,7 @@ end u, mesh, nonconservative_terms::True, equations, limiter::SubcellLimiterMCL, dg, element, cache) - @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes for j in eachnode(dg), i in 2:nnodes(dg) for v in eachvariable(equations) @@ -661,7 +661,7 @@ end if limiter isa SubcellLimiterIDP && !limiter.bar_states return nothing end - @unpack lambda1, lambda2, bar_states1, bar_states2 = limiter.cache.container_bar_states + (; lambda1, lambda2, bar_states1, bar_states2) = limiter.cache.container_bar_states # Calc lambdas and bar states inside elements @threaded for element in eachelement(dg, cache) @@ -860,8 +860,8 @@ end if !limiter.bar_states return nothing end - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - @unpack bar_states1, bar_states2 = limiter.cache.container_bar_states + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + (; bar_states1, bar_states2) = limiter.cache.container_bar_states # state variables if limiter.local_minmax @@ -911,27 +911,34 @@ end for j in eachnode(dg), i in eachnode(dg) s_min[i, j, element] = typemax(eltype(s_min)) end + # FV solution at node (i, j) for j in eachnode(dg), i in eachnode(dg) s = entropy_spec(get_node_vars(u, equations, dg, i, j, element), equations) s_min[i, j, element] = min(s_min[i, j, element], s) - # TODO: Add source? - # - xi direction + # TODO: Add source term! + end + # xi direction: subcell face between (i-1, j) and (i, j) + for j in eachnode(dg), i in 1:(nnodes(dg) + 1) s = entropy_spec(get_node_vars(bar_states1, equations, dg, i, j, element), equations) - s_min[i, j, element] = min(s_min[i, j, element], s) - # + xi direction - s = entropy_spec(get_node_vars(bar_states1, equations, dg, i + 1, j, - element), equations) - s_min[i, j, element] = min(s_min[i, j, element], s) - # - eta direction + if i <= nnodes(dg) + s_min[i, j, element] = min(s_min[i, j, element], s) + end + if i > 1 + s_min[i - 1, j, element] = min(s_min[i - 1, j, element], s) + end + end + # eta direction: subcell face between (i, j-1) and (i, j) + for j in 1:(nnodes(dg) + 1), i in eachnode(dg) s = entropy_spec(get_node_vars(bar_states2, equations, dg, i, j, element), equations) - s_min[i, j, element] = min(s_min[i, j, element], s) - # + eta direction - s = entropy_spec(get_node_vars(bar_states2, equations, dg, i, j + 1, - element), equations) - s_min[i, j, element] = min(s_min[i, j, element], s) + if j <= nnodes(dg) + s_min[i, j, element] = min(s_min[i, j, element], s) + end + if j > 1 + s_min[i, j - 1, element] = min(s_min[i, j - 1, element], s) + end end end end @@ -942,26 +949,34 @@ end for j in eachnode(dg), i in eachnode(dg) s_max[i, j, element] = typemin(eltype(s_max)) end + # FV solution at node (i, j) for j in eachnode(dg), i in eachnode(dg) s = entropy_math(get_node_vars(u, equations, dg, i, j, element), equations) s_max[i, j, element] = max(s_max[i, j, element], s) - # - xi direction + # TODO: Add source term! + end + # xi direction: subcell face between (i-1, j) and (i, j) + for j in eachnode(dg), i in 1:(nnodes(dg) + 1) s = entropy_math(get_node_vars(bar_states1, equations, dg, i, j, element), equations) - s_max[i, j, element] = max(s_max[i, j, element], s) - # + xi direction - s = entropy_math(get_node_vars(bar_states1, equations, dg, i + 1, j, - element), equations) - s_max[i, j, element] = max(s_max[i, j, element], s) - # - eta direction + if i <= nnodes(dg) + s_max[i, j, element] = max(s_max[i, j, element], s) + end + if i > 1 + s_max[i - 1, j, element] = max(s_max[i - 1, j, element], s) + end + end + # eta direction: subcell face between (i, j-1) and (i, j) + for j in 1:(nnodes(dg) + 1), i in eachnode(dg) s = entropy_math(get_node_vars(bar_states2, equations, dg, i, j, element), equations) - s_max[i, j, element] = max(s_max[i, j, element], s) - # + eta direction - s = entropy_math(get_node_vars(bar_states2, equations, dg, i, j + 1, - element), equations) - s_max[i, j, element] = max(s_max[i, j, element], s) + if j <= nnodes(dg) + s_max[i, j, element] = max(s_max[i, j, element], s) + end + if j > 1 + s_max[i, j - 1, element] = max(s_max[i, j - 1, element], s) + end end end end @@ -971,8 +986,8 @@ end @inline function calc_variable_bounds!(u, mesh, nonconservative_terms, equations, limiter::SubcellLimiterMCL, dg, cache) - @unpack var_min, var_max = limiter.cache.subcell_limiter_coefficients - @unpack bar_states1, bar_states2, lambda1, lambda2 = limiter.cache.container_bar_states + (; var_min, var_max) = limiter.cache.subcell_limiter_coefficients + (; bar_states1, bar_states2) = limiter.cache.container_bar_states @threaded for element in eachelement(dg, cache) for v in eachvariable(equations) @@ -1111,13 +1126,13 @@ end equations, limiter, dg, element, cache, fstar1, fstar2) - @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes - @unpack var_min, var_max = limiter.cache.subcell_limiter_coefficients - @unpack bar_states1, bar_states2, lambda1, lambda2 = limiter.cache.container_bar_states + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; var_min, var_max) = limiter.cache.subcell_limiter_coefficients + (; bar_states1, bar_states2, lambda1, lambda2) = limiter.cache.container_bar_states if limiter.Plotting - @unpack alpha, alpha_pressure, alpha_entropy, - alpha_mean, alpha_mean_pressure, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_pressure, alpha_entropy, alpha_mean, + alpha_mean_pressure, alpha_mean_entropy) = limiter.cache.subcell_limiter_coefficients for j in eachnode(dg), i in eachnode(dg) alpha_mean[:, i, j, element] .= zero(eltype(alpha_mean)) alpha[:, i, j, element] .= one(eltype(alpha)) @@ -1170,7 +1185,7 @@ end end if limiter.Plotting - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[1, i - 1, j, element] = min(alpha[1, i - 1, j, element], coefficient) alpha[1, i, j, element] = min(alpha[1, i, j, element], coefficient) @@ -1222,7 +1237,7 @@ end end if limiter.Plotting - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[1, i, j - 1, element] = min(alpha[1, i, j - 1, element], coefficient) alpha[1, i, j, element] = min(alpha[1, i, j, element], coefficient) @@ -1285,7 +1300,7 @@ end (g_limited + sign(g_limited) * eps()) / (g + sign(g_limited) * eps())) end - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[v, i - 1, j, element] = min(alpha[v, i - 1, j, element], coefficient) alpha[v, i, j, element] = min(alpha[v, i, j, element], coefficient) @@ -1336,7 +1351,7 @@ end (g_limited + sign(g_limited) * eps()) / (g + sign(g_limited) * eps())) end - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[v, i, j - 1, element] = min(alpha[v, i, j - 1, element], coefficient) alpha[v, i, j, element] = min(alpha[v, i, j, element], coefficient) @@ -1380,7 +1395,7 @@ end (antidiffusive_flux1_L[v, i, j, element] + sign(flux_limited) * eps())) end - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[v, i - 1, j, element] = min(alpha[v, i - 1, j, element], coefficient) alpha[v, i, j, element] = min(alpha[v, i, j, element], coefficient) @@ -1420,7 +1435,7 @@ end (antidiffusive_flux2_L[v, i, j, element] + sign(flux_limited) * eps())) end - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[v, i, j - 1, element] = min(alpha[v, i, j - 1, element], coefficient) alpha[v, i, j, element] = min(alpha[v, i, j, element], coefficient) @@ -1459,7 +1474,7 @@ end end if limiter.Plotting - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[1, i - 1, j, element] = min(alpha[1, i - 1, j, element], coefficient) alpha[1, i, j, element] = min(alpha[1, i, j, element], coefficient) @@ -1507,7 +1522,7 @@ end end if limiter.Plotting - @unpack alpha, alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha, alpha_mean) = limiter.cache.subcell_limiter_coefficients alpha[1, i, j - 1, element] = min(alpha[1, i, j - 1, element], coefficient) alpha[1, i, j, element] = min(alpha[1, i, j, element], coefficient) @@ -1534,7 +1549,7 @@ end # Divide alpha_mean by number of additions if limiter.Plotting - @unpack alpha_mean = limiter.cache.subcell_limiter_coefficients + (; alpha_mean) = limiter.cache.subcell_limiter_coefficients # Interfaces contribute with 1.0 if limiter.density_limiter || limiter.positivity_limiter_density for i in eachnode(dg) @@ -1564,7 +1579,7 @@ end # Limit pressure à la Kuzmin if limiter.positivity_limiter_pressure - @unpack alpha_pressure, alpha_mean_pressure = limiter.cache.subcell_limiter_coefficients + (; alpha_pressure, alpha_mean_pressure) = limiter.cache.subcell_limiter_coefficients for j in eachnode(dg), i in 2:nnodes(dg) bar_state_velocity = bar_states1[2, i, j, element]^2 + bar_states1[3, i, j, element]^2 @@ -1675,7 +1690,7 @@ end end end if limiter.Plotting - @unpack alpha_mean_pressure = limiter.cache.subcell_limiter_coefficients + (; alpha_mean_pressure) = limiter.cache.subcell_limiter_coefficients # Interfaces contribute with 1.0 for i in eachnode(dg) alpha_mean_pressure[i, 1, element] += 1.0 @@ -1732,7 +1747,7 @@ end end end if limiter.Plotting - @unpack alpha_entropy, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients + (; alpha_entropy, alpha_mean_entropy) = limiter.cache.subcell_limiter_coefficients alpha_entropy[i - 1, j, element] = min(alpha_entropy[i - 1, j, element], alpha) alpha_entropy[i, j, element] = min(alpha_entropy[i, j, element], alpha) @@ -1780,7 +1795,7 @@ end end end if limiter.Plotting - @unpack alpha_entropy, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients + (; alpha_entropy, alpha_mean_entropy) = limiter.cache.subcell_limiter_coefficients alpha_entropy[i, j - 1, element] = min(alpha_entropy[i, j - 1, element], alpha) alpha_entropy[i, j, element] = min(alpha_entropy[i, j, element], alpha) @@ -1789,7 +1804,7 @@ end end end if limiter.Plotting - @unpack alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients + (; alpha_mean_entropy) = limiter.cache.subcell_limiter_coefficients # Interfaces contribute with 1.0 for i in eachnode(dg) alpha_mean_entropy[i, 1, element] += 1.0 @@ -1816,44 +1831,29 @@ end return nothing end +@inline function get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_slip_wall), + orientation::Integer, direction, + equations, dg, indices...) + return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) +end + +@inline function get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_slip_wall), + normal_direction::AbstractVector, + direction, equations, dg, indices...) + u_rotate = rotate_to_x(u_inner, normal_direction, equations) + + return SVector(u_inner[1], + u_inner[2] - 2.0 * u_rotate[2], + u_inner[3] - 2.0 * u_rotate[3], + u_inner[4]) +end + +# Default implementation of `get_boundary_outer_state` returns inner value. @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition, orientation_or_normal, direction, equations, dg, indices...) - if boundary_condition == boundary_condition_slip_wall #boundary_condition_reflecting_euler_wall - if orientation_or_normal isa AbstractArray - u_rotate = rotate_to_x(u_inner, orientation_or_normal, equations) - - return SVector(u_inner[1], - u_inner[2] - 2.0 * u_rotate[2], - u_inner[3] - 2.0 * u_rotate[3], - u_inner[4]) - else # orientation_or_normal isa Integer - return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) - end - elseif boundary_condition == boundary_condition_mixed_dirichlet_wall - x = get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) - if x[1] < 1 / 6 # BoundaryConditionCharacteristic - u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, - u_inner, - orientation_or_normal, - direction, x, t, - equations) - - return u_outer - else # x[1] >= 1 / 6 # boundary_condition_slip_wall - if orientation_or_normal isa AbstractArray - u_rotate = rotate_to_x(u_inner, orientation_or_normal, equations) - - return SVector(u_inner[1], - u_inner[2] - 2.0 * u_rotate[2], - u_inner[3] - 2.0 * u_rotate[3], - u_inner[4]) - else # orientation_or_normal isa Integer - return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) - end - end - end - return u_inner end diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index c1c31a7a48f..e73e475a73f 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -270,7 +270,7 @@ function calc_normal_directions!(container_bar_states, mesh::TreeMesh, equations end function calc_normal_directions!(container_bar_states, - mesh::Union{StructuredMesh, P4estMesh}, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, equations, dg, cache) (; weights, derivative_matrix) = dg.basis (; contravariant_vectors) = cache.elements diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 152ca52a6ca..6632cd0bb27 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -541,6 +541,38 @@ end end end +@trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_advection_diffusion_periodic_amr.jl"), + tspan=(0.0, 0.01), + l2=[0.014715887539773128], + linf=[0.2285802791900049]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_advection_diffusion_nonperiodic_amr.jl"), + tspan=(0.0, 0.01), + l2=[0.00793438523666649], + linf=[0.11030633127144573]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +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"), @@ -635,6 +667,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_lid_driven_cavity_amr.jl"), + tspan=(0.0, 1.0), + l2=[ + 0.0005323841980601085, 0.07892044543547208, + 0.02909671646389337, 0.11717468256112017, + ], + linf=[ + 0.006045292737899444, 0.9233292581786228, + 0.7982129977236198, 1.6864546235292153, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 56d8e0dcbfa..50f032e9e0a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -247,6 +247,33 @@ end end end +@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 1.5123651627525257e-5, + 1.51236516273878e-5, + 2.4544918394022538e-5, + 5.904791661362391e-6, + 1.1809583322724782e-5, + ], + linf=[ + 8.393471747591974e-5, + 8.393471748258108e-5, + 0.00015028562494778797, + 3.504466610437795e-5, + 7.00893322087559e-5, + ]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_convergence_wavingflag_IDP.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_wavingflag_IDP.jl"), @@ -615,34 +642,42 @@ end @trixi_testset "elixir_euler_double_mach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach.jl"), l2=[ - 0.8955457632754655, - 6.8117495933240235, - 3.2697118944675716, - 77.5174041919109, + 0.8895735118233061, + 6.7732502383794655, + 3.2537120458330913, + 77.15211199077615, ], linf=[ - 10.16165871096883, - 133.2522870057006, - 38.23157147773949, - 1470.3950960145828, + 10.595204890391603, + 133.14647204580373, + 39.59427444784427, + 1467.649170200717, ], initial_refinement_level=3, tspan=(0.0, 0.05)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 10000 + end end @trixi_testset "elixir_euler_double_mach_MCL.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach_MCL.jl"), l2=[ - 0.9266313242695542, - 7.071517579972717, - 3.2627078543492787, - 80.24631724351916, + 0.9230564188859353, + 7.048067525695328, + 3.2513897966486693, + 80.02373234685636, ], linf=[ - 14.244598580563007, - 138.4745277257612, - 38.69633620234036, - 1574.6686216469134, + 13.994085024959638, + 138.49691741864086, + 38.68709126443809, + 1574.7770582279284, ], initial_refinement_level=3, tspan=(0.0, 0.05)) @@ -660,16 +695,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shock_upstream_sc_subcell.jl"), l2=[ - 1.2351468819080416, - 1.1269856428294935, - 1.7239124305681928, - 11.715260007491556, + 1.2351650580585378, + 1.1268555757641623, + 1.723421391322672, + 11.715156856828612, ], linf=[ - 5.385493056976312, - 6.575446146030286, - 10.06523457762742, - 51.00903155017642, + 5.378363125399197, + 6.5624165241817245, + 10.008377369949162, + 51.24533429408666, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations @@ -685,18 +720,17 @@ end @trixi_testset "elixir_euler_shock_upstream_MCL.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shock_upstream_MCL.jl"), l2=[ - 1.2607430289877726, - 1.1565837325291355, - 1.7791790302458714, - 11.891223800389232, + 1.2607430289878367, + 1.1565837325291748, + 1.779179030245946, + 11.891223800389836, ], linf=[ - 5.68876088477983, - 8.165554425950146, - 10.859100194836538, - 50.25822408989214, + 5.68876088478026, + 8.16555442595031, + 10.859100194836678, + 50.258224089906975, ], - cells_per_dimension=(8, 12), tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index af91a365458..f95f52a031b 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -296,9 +296,11 @@ end end Trixi.move_connectivity!(c::MyContainer, first, last, destination) = c Trixi.delete_connectivity!(c::MyContainer, first, last) = c - Trixi.reset_data_structures!(c::MyContainer) = (c.data = Vector{Int}(undef, - c.capacity + 1); - c) + function Trixi.reset_data_structures!(c::MyContainer) + (c.data = Vector{Int}(undef, + c.capacity + 1); + c) + end function Base.:(==)(c1::MyContainer, c2::MyContainer) return (c1.capacity == c2.capacity && c1.length == c2.length && @@ -617,6 +619,18 @@ end @test_throws ArgumentError TimeSeriesCallback(semi, [1.0 1.0 1.0; 2.0 2.0 2.0]) end +@timed_testset "Consistency check for single point flux: CEMCE" begin + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, 0.4)) + u = SVector(0.1, -0.5, 1.0, 1.0, 2.0) + + orientations = [1, 2] + for orientation in orientations + @test flux(u, orientation, equations) ≈ + flux_ranocha(u, u, orientation, equations) + end +end + @timed_testset "Consistency check for HLL flux (naive): CEE" begin flux_hll = FluxHLL(min_max_speed_naive) @@ -1227,6 +1241,29 @@ end end @testset "FluxRotated vs. direct implementation" begin + @timed_testset "CompressibleEulerMulticomponentEquations2D" begin + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, + 0.4)) + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0), + SVector(-0.1, -0.3, 1.2, 1.3, 1.4)] + + f_std = flux + f_rot = FluxRotated(f_std) + println(typeof(f_std)) + println(typeof(f_rot)) + for u in u_values, + normal_direction in normal_directions + + @test f_rot(u, normal_direction, equations) ≈ + f_std(u, normal_direction, equations) + end + end + @timed_testset "CompressibleEulerEquations2D" begin equations = CompressibleEulerEquations2D(1.4) normal_directions = [SVector(1.0, 0.0),