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_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_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/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index a465571989b..eca06c8c203 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -285,7 +285,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, f_ll = flux(u_ll, orientation, equations) f_rr = flux(u_rr, orientation, equations) - SsL, SsR = min_max_speed_naive(u_ll, u_rr, orientation, equations) + SsL, SsR = min_max_speed_einfeldt(u_ll, u_rr, orientation, equations) sMu_L = SsL - v1_ll sMu_R = SsR - v1_rr if SsL >= 0 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/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 96202e00f58..1addc29e3e6 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -247,6 +247,32 @@ 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_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), # Expected errors are exactly the same as with TreeMesh! diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index 447572eee88..8895fe30e8b 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -210,15 +210,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"), surface_flux=flux_hllc, l2=[ - 0.4574266553239646, 0.4794143154876439, 0.3407079689595056, - 0.44797768430829343, 0.9206916204424165, - 1.3216517820475193e-16, 0.2889748702415378, - 0.25529778018020927, + 0.45738965718253993, 0.479402222862685, 0.34069729746967664, + 0.44795514335568865, 0.9206813325913135, + 1.3216517820475193e-16, 0.2889672868491632, + 0.2552794220777942, ], linf=[ - 1.217943947570543, 0.8868438459815245, 0.878215340656725, - 0.9710882819266371, 1.6742759645320984, - 2.220446049250313e-16, 0.704710220504591, 0.6562122176458641, + 1.2181099854251536, 0.8869319941747589, 0.8763562906332134, + 0.9712221036087284, 1.6734231113527818, + 2.220446049250313e-16, 0.7035011427822779, + 0.6562884129650286, ]) # 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 707315bc353..080f387e8fc 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 && @@ -612,6 +614,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) @@ -1242,6 +1256,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),