diff --git a/.github/review-checklist.md b/.github/review-checklist.md new file mode 100644 index 00000000000..2d8a24f1971 --- /dev/null +++ b/.github/review-checklist.md @@ -0,0 +1,38 @@ +### Review checklist + +This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging. + +#### Purpose and scope +- [ ] The PR has a single goal that is clear from the PR title and/or description. +- [ ] All code changes represent a single set of modifications that logically belong together. +- [ ] No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs. + +#### Code quality +- [ ] The code can be understood easily. +- [ ] Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions. +- [ ] There are no redundancies that can be removed by simple modularization/refactoring. +- [ ] There are no leftover debug statements or commented code sections. +- [ ] The code adheres to our [conventions](https://trixi-framework.github.io/Trixi.jl/stable/conventions/) and [style guide](https://trixi-framework.github.io/Trixi.jl/stable/styleguide/), and to the [Julia guidelines](https://docs.julialang.org/en/v1/manual/style-guide/). + +#### Documentation +- [ ] New functions and types are documented with a docstring or top-level comment. +- [ ] Relevant publications are referenced in docstrings (see [example](https://github.com/trixi-framework/Trixi.jl/blob/7f83a1a938eecd9b841efe215a6e482e67cfdcc1/src/equations/compressible_euler_2d.jl#L601-L615) for formatting). +- [ ] Inline comments are used to document longer or unusual code sections. +- [ ] Comments describe intent ("why?") and not just functionality ("what?"). +- [ ] If the PR introduces a significant change or new feature, it is documented in `NEWS.md`. + +#### Testing +- [ ] The PR passes all tests. +- [ ] New or modified lines of code are covered by tests. +- [ ] New or modified tests run in less then 10 seconds. + +#### Performance +- [ ] There are no type instabilities or memory allocations in performance-critical parts. +- [ ] If the PR intent is to improve performance, before/after [time measurements](https://trixi-framework.github.io/Trixi.jl/stable/performance/#Manual-benchmarking) are posted in the PR. + +#### Verification +- [ ] The correctness of the code was verified using appropriate tests. +- [ ] If new equations/methods are added, a convergence test has been run and the results + are posted in the PR. + +*Created with :heart: by the Trixi.jl community.* \ No newline at end of file diff --git a/.github/workflows/ReviewChecklist.yml b/.github/workflows/ReviewChecklist.yml new file mode 100644 index 00000000000..959a04752d7 --- /dev/null +++ b/.github/workflows/ReviewChecklist.yml @@ -0,0 +1,20 @@ +name: Add review checklist + +on: + pull_request_target: + types: [opened] + +permissions: + pull-requests: write + +jobs: + add-review-checklist: + runs-on: ubuntu-latest + steps: + - name: Check out repository + uses: actions/checkout@v3 + - name: Add review checklist + uses: trixi-framework/add-pr-review-checklist@v1 + with: + file: '.github/review-checklist.md' + no-checklist-keyword: '[no checklist]' \ No newline at end of file diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 6ebb288ea30..a06121e7ca1 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v3 - name: Check spelling - uses: crate-ci/typos@v1.16.5 + uses: crate-ci/typos@v1.16.9 diff --git a/.zenodo.json b/.zenodo.json index 95879af1e90..905c0170ab9 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -15,7 +15,7 @@ "orcid": "0000-0002-1752-1158" }, { - "affiliation": "Applied Mathematics, University of Hamburg, Germany", + "affiliation": "Numerical Mathematics, Johannes Gutenberg University Mainz, Germany", "name": "Ranocha, Hendrik", "orcid": "0000-0002-3456-2277" }, diff --git a/AUTHORS.md b/AUTHORS.md index 74bfaa9c852..f1debf8ba76 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -12,7 +12,7 @@ provided substantial additions or modifications. Together, these two groups form * [Gregor Gassner](https://www.mi.uni-koeln.de/NumSim/gregor-gassner), University of Cologne, Germany * [Hendrik Ranocha](https://ranocha.de), - University of Hamburg, Germany + Johannes Gutenberg University Mainz, Germany * [Andrew Winters](https://liu.se/en/employee/andwi94), Linköping University, Sweden * [Jesse Chan](https://jlchan.github.io), diff --git a/NEWS.md b/NEWS.md index 10125c40d17..4b96e1e2834 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,7 @@ for human readability. - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations. - Capability to set truly discontinuous initial conditions in 1D. - Wetting and drying feature and examples for 1D and 2D shallow water equations +- Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` #### Changed diff --git a/Project.toml b/Project.toml index dd937ed213b..d37c0548a6a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.39-pre" +version = "0.5.42-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" @@ -56,7 +56,7 @@ DiffEqCallbacks = "2.25" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.18" -HDF5 = "0.14, 0.15, 0.16" +HDF5 = "0.14, 0.15, 0.16, 0.17" IfElse = "0.1" LinearMaps = "2.7, 3.0" LoopVectorization = "0.12.118" diff --git a/README.md b/README.md index 7eaee8750dd..63540b1f640 100644 --- a/README.md +++ b/README.md @@ -247,7 +247,7 @@ Schlottke-Lakemper](https://lakemper.eu) (RWTH Aachen University/High-Performance Computing Center Stuttgart (HLRS), Germany) and [Gregor Gassner](https://www.mi.uni-koeln.de/NumSim/gregor-gassner) (University of Cologne, Germany). Together with [Hendrik Ranocha](https://ranocha.de) -(University of Hamburg, Germany), [Andrew Winters](https://liu.se/en/employee/andwi94) +(Johannes Gutenberg University Mainz, Germany), [Andrew Winters](https://liu.se/en/employee/andwi94) (Linköping University, Sweden), and [Jesse Chan](https://jlchan.github.io) (Rice University, US), they are the principal developers of Trixi.jl. The full list of contributors can be found in [AUTHORS.md](AUTHORS.md). diff --git a/docs/literate/src/files/differentiable_programming.jl b/docs/literate/src/files/differentiable_programming.jl index ecc09d05dcf..5c5a7cd7440 100644 --- a/docs/literate/src/files/differentiable_programming.jl +++ b/docs/literate/src/files/differentiable_programming.jl @@ -128,7 +128,7 @@ condition_number = cond(V) # you can compute the gradient of an entropy-dissipative semidiscretization with respect to the # ideal gas constant of the compressible Euler equations as described in the following. This example # is also available as the elixir -# [examples/special\_elixirs/elixir\_euler\_ad.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/special_elixirs/elixir_euler_ad.jl) +# [`examples/special_elixirs/elixir_euler_ad.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/special_elixirs/elixir_euler_ad.jl) # First, we create a semidiscretization of the compressible Euler equations. diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index 5b669881502..0c8de66bf42 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -116,7 +116,7 @@ # ## Examples in Trixi.jl # Trixi.jl already contains several more coding examples, the so-called `elixirs`. You can find them -# in the folder [`examples`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/). +# in the folder [`examples/`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/). # They are structured by the underlying mesh type and the respective number of spatial dimensions. # The name of an elixir is composed of the underlying system of conservation equations (for instance # `advection` or `euler`) and other special characteristics like the initial condition diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md index a85f8e8191b..1d3e5e34b51 100644 --- a/docs/src/callbacks.md +++ b/docs/src/callbacks.md @@ -15,7 +15,7 @@ control, adaptive mesh refinement, I/O, and more. ### CFL-based time step control Time step control can be performed with a [`StepsizeCallback`](@ref). An example making use -of this can be found at [examples/tree_2d_dgsem/elixir\_advection\_basic.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_basic.jl) +of this can be found at [`examples/tree_2d_dgsem/elixir_advection_basic.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_basic.jl) ### Adaptive mesh refinement Trixi.jl uses a hierarchical Cartesian mesh which can be locally refined in a solution-adaptive way. @@ -24,12 +24,12 @@ passing an [`AMRCallback`](@ref) to the ODE solver. The `AMRCallback` requires a [`ControllerThreeLevel`](@ref) or [`ControllerThreeLevelCombined`](@ref) to tell the AMR algorithm which cells to refine/coarsen. -An example elixir using AMR can be found at [examples/tree_2d_dgsem/elixir\_advection\_amr.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_amr.jl). +An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_advection_amr.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_amr.jl). ### Analyzing the numerical solution The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate errors or user-specified integrals, and print the results to the screen. The results can also be -saved in a file. An example can be found at [examples/tree_2d_dgsem/elixir\_euler\_vortex.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl). +saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl). In [Performance metrics of the `AnalysisCallback`](@ref) you can find a detailed description of the different performance metrics the `AnalysisCallback` computes. @@ -38,15 +38,15 @@ description of the different performance metrics the `AnalysisCallback` computes #### Solution and restart files To save the solution in regular intervals you can use a [`SaveSolutionCallback`](@ref). It is also possible to create restart files using the [`SaveRestartCallback`](@ref). An example making use -of these can be found at [examples/tree_2d_dgsem/elixir\_advection\_extended.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_extended.jl). +of these can be found at [`examples/tree_2d_dgsem/elixir_advection_extended.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_extended.jl). An example showing how to restart a simulation from a restart file can be found at -[examples/tree_2d_dgsem/elixir\_advection\_restart.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_restart.jl). +[`examples/tree_2d_dgsem/elixir_advection_restart.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_restart.jl). #### Time series Sometimes it is useful to record the evaluations of state variables over time at a given set of points. This can be achieved by the [`TimeSeriesCallback`](@ref), which is used, e.g., in -[examples/tree_2d_dgsem/elixir\_acoustics\_gaussian\_source.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl). +[`examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl). The `TimeSeriesCallback` constructor expects a semidiscretization and a list of points at which the solution should be recorded in regular time step intervals. After the last time step, the entire record is stored in an HDF5 file. @@ -113,12 +113,12 @@ will yield the following plot: Some callbacks provided by Trixi.jl implement specific features for certain equations: * The [`LBMCollisionCallback`](@ref) implements the Lattice-Boltzmann method (LBM) collision operator and should only be used when solving the Lattice-Boltzmann equations. See e.g. - [examples/tree_2d_dgsem/elixir\_lbm\_constant.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_lbm_constant.jl) + [`examples/tree_2d_dgsem/elixir_lbm_constant.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_lbm_constant.jl) * The [`SteadyStateCallback`](@ref) terminates the time integration when the residual steady state falls below a certain threshold. This checks the convergence of the potential ``\phi`` for - hyperbolic diffusion. See e.g. [examples/tree_2d_dgsem/elixir\_hypdiff\_nonperiodic.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_hypdiff_nonperiodic.jl). + hyperbolic diffusion. See e.g. [`examples/tree_2d_dgsem/elixir_hypdiff_nonperiodic.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_hypdiff_nonperiodic.jl). * The [`GlmSpeedCallback`](@ref) updates the divergence cleaning wave speed `c_h` for the ideal - GLM-MHD equations. See e.g. [examples/tree_2d_dgsem/elixir\_mhd\_alfven\_wave.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_mhd_alfven_wave.jl). + GLM-MHD equations. See e.g. [`examples/tree_2d_dgsem/elixir_mhd_alfven_wave.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_mhd_alfven_wave.jl). ## Usage of step callbacks Step callbacks are passed to the `solve` method from the ODE solver via the keyword argument @@ -152,7 +152,7 @@ more callbacks, you need to turn them into a `CallbackSet` first by calling ## Stage callbacks [`PositivityPreservingLimiterZhangShu`](@ref) is a positivity-preserving limiter, used to enforce physical constraints. An example elixir using this feature can be found at -[examples/tree_2d_dgsem/elixir\_euler\_positivity.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_positivity.jl). +[`examples/tree_2d_dgsem/elixir_euler_positivity.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_positivity.jl). ## Implementing new callbacks Since Trixi.jl is compatible with [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), @@ -162,4 +162,4 @@ Step callbacks are just called [callbacks](https://diffeq.sciml.ai/latest/featur Stage callbacks are called [`stage_limiter!`](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)). An example elixir showing how to implement a new simple stage callback and a new simple step -callback can be found at [examples/tree_2d_dgsem/elixir\_advection\_callbacks.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_callbacks.jl). +callback can be found at [`examples/tree_2d_dgsem/elixir_advection_callbacks.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_callbacks.jl). diff --git a/docs/src/index.md b/docs/src/index.md index 3af785bc681..bb2afd1019f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -324,7 +324,7 @@ Schlottke-Lakemper](https://lakemper.eu) (RWTH Aachen University/High-Performance Computing Center Stuttgart (HLRS), Germany) and [Gregor Gassner](https://www.mi.uni-koeln.de/NumSim/gregor-gassner) (University of Cologne, Germany). Together with [Hendrik Ranocha](https://ranocha.de) -(University of Hamburg, Germany) and [Andrew Winters](https://liu.se/en/employee/andwi94) +(Johannes Gutenberg University Mainz, Germany) and [Andrew Winters](https://liu.se/en/employee/andwi94) (Linköping University, Sweden), and [Jesse Chan](https://jlchan.github.io) (Rice University, US), they are the principal developers of Trixi.jl. The full list of contributors can be found under [Authors](@ref). diff --git a/docs/src/meshes/dgmulti_mesh.md b/docs/src/meshes/dgmulti_mesh.md index e07ba70a80a..fc086bba146 100644 --- a/docs/src/meshes/dgmulti_mesh.md +++ b/docs/src/meshes/dgmulti_mesh.md @@ -81,16 +81,20 @@ type, but will be more efficient at high orders of approximation. ## Trixi.jl elixirs on simplicial and tensor product element meshes Example elixirs with triangular, quadrilateral, and tetrahedral meshes can be found in -the `examples/dgmulti_2d` and `examples/dgmulti_3d` folders. Some key elixirs to look at: +the [`examples/dgmulti_2d/`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_2d/) +and [`examples/dgmulti_3d/`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_3d/) +folders. Some key elixirs to look at: -* `examples/dgmulti_2d/elixir_euler_weakform.jl`: basic weak form DG discretization on a uniform triangular mesh. +* [`examples/dgmulti_2d/elixir_euler_weakform.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_2d/elixir_euler_weakform.jl): + basic weak form DG discretization on a uniform triangular mesh. Changing `element_type = Quad()` or `approximation_type = SBP()` will switch to a quadrilateral mesh or an SBP-type discretization. Changing `surface_integral = SurfaceIntegralWeakForm(flux_ec)` and `volume_integral = VolumeIntegralFluxDifferencing(flux_ec)` for some entropy conservative flux (e.g., [`flux_chandrashekar`](@ref) or [`flux_ranocha`](@ref)) will switch to an entropy conservative formulation. -* `examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl`: uses an unstructured mesh generated by - [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl). -* `examples/dgmulti_3d/elixir_euler_weakform.jl`: basic weak form DG discretization on a uniform tetrahedral mesh. +* [`examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl): + uses an unstructured mesh generated by [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl). +* [`examples/dgmulti_3d/elixir_euler_weakform.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_3d/elixir_euler_weakform.jl): + ´basic weak form DG discretization on a uniform tetrahedral mesh. Changing `element_type = Hex()` will switch to a hexahedral mesh. Changing `surface_integral = SurfaceIntegralWeakForm(flux_ec)` and `volume_integral = VolumeIntegralFluxDifferencing(flux_ec)` for some entropy conservative flux diff --git a/docs/src/overview.md b/docs/src/overview.md index 519ec2ca424..46bc28b6025 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -5,7 +5,7 @@ conservation laws. Thus, it is not a monolithic PDE solver that is configured at via parameter files, as it is often found in classical numerical simulation codes. Instead, each simulation is configured by pure Julia code. Many examples of such simulation setups, called *elixirs* in Trixi.jl, are provided in the -[examples](https://github.com/trixi-framework/Trixi.jl/blob/main/examples) +[`examples/`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples) folder. Trixi.jl uses the method of lines, i.e., the full space-time discretization is separated into two steps; @@ -77,7 +77,7 @@ Further information can be found in the ## Next steps We explicitly encourage people interested in Trixi.jl to have a look at the -[examples](https://github.com/trixi-framework/Trixi.jl/blob/main/examples) +[`examples/`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples) bundled with Trixi.jl to get an impression of what is possible and the general look and feel of Trixi.jl. Before doing that, it is usually good to get an idea of diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md index 08470fd064a..d56777c9af4 100644 --- a/docs/src/parallelization.md +++ b/docs/src/parallelization.md @@ -53,16 +53,24 @@ 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) -library needs to be compiled with the same MPI installation. Therefore, you also need to use -a system-provided `p4est` installation (for notes on how to install `p4est` see e.g. -[here](https://github.com/cburstedde/p4est/blob/master/README), use the configure option -`--enable-mpi`). In addition, [P4est.jl](https://github.com/trixi-framework/P4est.jl) needs to -be configured to use the custom `p4est` installation. Follow the steps described -[here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md) for the configuration. +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]` and `[P4est]` as well as an entry `MPIPreferences` in your Project.toml to -use a custom MPI installation. +`[MPIPreferences]`, `[T8code]` and `[P4est]` as well as an entry `MPIPreferences` in your +Project.toml to use a custom MPI installation. ### [Usage](@id parallel_usage) @@ -158,17 +166,36 @@ section, specifically at the descriptions of the performance index (PID). ### Using error-based step size control with MPI -If you use error-based step size control (see also the section on [error-based adaptive step sizes](@ref adaptive_step_sizes)) -together with MPI you need to pass `internalnorm=ode_norm` and you should pass -`unstable_check=ode_unstable_check` to OrdinaryDiffEq's [`solve`](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/), +If you use error-based step size control (see also the section on +[error-based adaptive step sizes](@ref adaptive_step_sizes)) together with MPI you need to pass +`internalnorm=ode_norm` and you should pass `unstable_check=ode_unstable_check` to +OrdinaryDiffEq's [`solve`](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/), which are both included in [`ode_default_options`](@ref). ### Using parallel input and output -Trixi.jl allows parallel I/O using MPI by leveraging parallel HDF5.jl. To enable this, you first need -to use a system-provided MPI library, see also [here](@ref parallel_system_MPI) and you need to tell -[HDF5.jl](https://github.com/JuliaIO/HDF5.jl) to use this library. -To do so, set the environment variable `JULIA_HDF5_PATH` to the local path -that contains the `libhdf5.so` shared object file and build HDF5.jl by executing `using Pkg; Pkg.build("HDF5")`. -For more information see also the [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). - -If you do not perform these steps to use parallel HDF5 or if the HDF5 is not MPI-enabled, Trixi.jl will fall back on a less efficient I/O mechanism. In that case, all disk I/O is performed only on rank zero and data is distributed to/gathered from the other ranks using regular MPI communication. +Trixi.jl allows parallel I/O using MPI by leveraging parallel HDF5.jl. On most systems, this is +enabled by default. Additionally, you can also use a local installation of the HDF5 library +(with MPI support). For this, you first need to use a system-provided MPI library, see also +[here](@ref parallel_system_MPI) and you need to tell [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) +to use this library. To do so with HDF5.jl v0.17 and newer, set the preferences `libhdf5` and +`libhdf5_hl` to the local paths of the libraries `libhdf5` and `libhdf5_hl`, which can be done by +```julia +julia> using Preferences, UUIDs +julia> set_preferences!( + UUID("f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"), # UUID of HDF5.jl + "libhdf5" => "/path/to/your/libhdf5.so", + "libhdf5_hl" => "/path/to/your/libhdf5_hl.so", force = true) +``` +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 +`[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`. +If you use HDF5.jl v0.16 or older, instead of setting the preferences for HDF5.jl, you need to set +the environment variable `JULIA_HDF5_PATH` to the path, where the HDF5 binaries are located and +then call `]build HDF5` from Julia. + +If HDF5 is not MPI-enabled, Trixi.jl will fall back on a less efficient I/O mechanism. In that +case, all disk I/O is performed only on rank zero and data is distributed to/gathered from the +other ranks using regular MPI communication. diff --git a/docs/src/restart.md b/docs/src/restart.md index d24d93cb297..767269ff27d 100644 --- a/docs/src/restart.md +++ b/docs/src/restart.md @@ -18,7 +18,7 @@ save_restart = SaveRestartCallback(interval=100, Make this part of your `CallbackSet`. An example is -[```examples/examples/structured_2d_dgsem/elixir_advection_extended.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl). +[`examples/examples/structured_2d_dgsem/elixir_advection_extended.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl). ## [Perform the simulation restart](@id restart_perform) @@ -26,7 +26,7 @@ Since all of the information about the simulation can be obtained from the last snapshot, the restart can be done with relatively few lines in an extra elixir file. However, some might prefer to keep everything in one elixir and -conditionals like ```if restart``` with a boolean variable ```restart``` that is user defined. +conditionals like `if restart` with a boolean variable `restart` that is user defined. First we need to define from which file we want to restart, e.g. ```julia @@ -50,7 +50,7 @@ time the one form the snapshot: tspan = (load_time(restart_filename), 2.0) ``` -We now also take the last ```dt```, so that our solver does not need to first find +We now also take the last `dt`, so that our solver does not need to first find one to fulfill the CFL condition: ```julia dt = load_dt(restart_filename) @@ -63,7 +63,7 @@ ode = semidiscretize(semi, tspan, restart_filename) You should now define a [`SaveSolutionCallback`](@ref) similar to the [original simulation](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl), -but with ```save_initial_solution=false```, otherwise our initial snapshot will be overwritten. +but with `save_initial_solution=false`, otherwise our initial snapshot will be overwritten. If you are using one file for the original simulation and the restart you can reuse your [`SaveSolutionCallback`](@ref), but need to set ```julia @@ -86,4 +86,4 @@ Now we can compute the solution: sol = solve!(integrator) ``` -An example is in `[``examples/structured_2d_dgsem/elixir_advection_restart.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_restart.jl). +An example is in [`examples/structured_2d_dgsem/elixir_advection_restart.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_restart.jl). diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 8f72bb4b1c6..36a7e8f5ac8 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -375,7 +375,7 @@ During the simulation, the visualization callback creates and displays visualizations of the current solution in regular intervals. This can be useful to, e.g., monitor the validity of a long-running simulation or for illustrative purposes. An example for how to create a `VisualizationCallback` can be found in -[examples/tree_2d_dgsem/elixir\_advection\_amr\_visualization.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl): +[`examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl): ```julia [...] diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl new file mode 100644 index 00000000000..6b69e4db563 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -0,0 +1,92 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave (modified to lower density and higher pressure) taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> modified to lower density, higher pressure + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables "normal" medium blast wave + rho = r > 0.5 ? 0.1 : 0.2691 # rho = r > 0.5 ? 1 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0E-1 : 1.245 # p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons=[1], + positivity_correction_factor=0.5) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.0, -2.0) +coordinates_max = ( 2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=5, + n_cells_max=100_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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.6) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); + 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_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl new file mode 100644 index 00000000000..a67eaeb5b2b --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl @@ -0,0 +1,140 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler multicomponent equations + +# 1) Dry Air 2) Helium + 28% Air +equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), + gas_constants = (0.287, 1.578)) + +""" + initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2}) + +A shock-bubble testcase for multicomponent Euler equations +- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman + Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations + [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972) +""" +function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2}) + # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972 + # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf + # typical domain is rectangular, we change it to a square, as Trixi can only do squares + @unpack gas_constants = equations + + # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving + delta = 0.03 + + # Region I + rho1_1 = delta + rho2_1 = 1.225 * gas_constants[1]/gas_constants[2] - delta + v1_1 = zero(delta) + v2_1 = zero(delta) + p_1 = 101325 + + # Region II + rho1_2 = 1.225-delta + rho2_2 = delta + v1_2 = zero(delta) + v2_2 = zero(delta) + p_2 = 101325 + + # Region III + rho1_3 = 1.6861 - delta + rho2_3 = delta + v1_3 = -113.5243 + v2_3 = zero(delta) + p_3 = 159060 + + # Set up Region I & II: + inicenter = SVector(zero(delta), zero(delta)) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + if (x[1] > 0.50) + # Set up Region III + rho1 = rho1_3 + rho2 = rho2_3 + v1 = v1_3 + v2 = v2_3 + p = p_3 + elseif (r < 0.25) + # Set up Region I + rho1 = rho1_1 + rho2 = rho2_1 + v1 = v1_1 + v2 = v2_1 + p = p_1 + else + # Set up Region II + rho1 = rho1_2 + rho2 = rho2_2 + v1 = v1_2 + v2 = v2_2 + p = p_2 + end + + return prim2cons(SVector(v1, v2, p, rho1, rho2), equations) +end +initial_condition = initial_condition_shock_bubble + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) + +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons=[(i+3 for i in eachcomponent(equations))...]) + +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.25, -2.225) +coordinates_max = ( 2.20, 2.225) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=1_000_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.01) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 300 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_integrals=(Trixi.density,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=300, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); + 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 \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index 78ddaa3ca7f..ec4d20558e5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -121,10 +121,10 @@ include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("semidiscretization/semidiscretization_coupled.jl") +include("time_integration/time_integration.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") -include("time_integration/time_integration.jl") # `trixi_include` and special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") @@ -229,6 +229,9 @@ export DG, SurfaceIntegralUpwind, MortarL2 +export VolumeIntegralSubcellLimiting, + SubcellLimiterIDP, SubcellLimiterIDPCorrection + export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index ab0f34efb78..976af327e6f 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -6,6 +6,7 @@ #! format: noindent include("positivity_zhang_shu.jl") +include("subcell_limiter_idp_correction.jl") # TODO: TrixiShallowWater: move specific limiter file include("positivity_shallow_water.jl") end # @muladd diff --git a/src/callbacks_stage/subcell_limiter_idp_correction.jl b/src/callbacks_stage/subcell_limiter_idp_correction.jl new file mode 100644 index 00000000000..69125ebecd9 --- /dev/null +++ b/src/callbacks_stage/subcell_limiter_idp_correction.jl @@ -0,0 +1,69 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + SubcellLimiterIDPCorrection() + +Perform antidiffusive correction stage for the a posteriori IDP limiter [`SubcellLimiterIDP`](@ref) +called with [`VolumeIntegralSubcellLimiting`](@ref). + +!!! note + This callback and the actual limiter [`SubcellLimiterIDP`](@ref) only work together. + This is not a replacement but a necessary addition. + +## References + +- Rueda-Ramírez, Pazner, Gassner (2022) + Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods + [DOI: 10.1016/j.compfluid.2022.105627](https://doi.org/10.1016/j.compfluid.2022.105627) +- Pazner (2020) + Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting + [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876) + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +struct SubcellLimiterIDPCorrection end + +function (limiter!::SubcellLimiterIDPCorrection)(u_ode, + integrator::Trixi.SimpleIntegratorSSP, + stage) + semi = integrator.p + limiter!(u_ode, semi, integrator.t, integrator.dt, + semi.solver.volume_integral) +end + +function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt, + volume_integral::VolumeIntegralSubcellLimiting) + @trixi_timeit timer() "a posteriori limiter" limiter!(u_ode, semi, t, dt, + volume_integral.limiter) +end + +function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt, + limiter::SubcellLimiterIDP) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + + u = wrap_array(u_ode, mesh, equations, solver, cache) + + # Calculate blending factor alpha in [0,1] + # f_ij = alpha_ij * f^(FV)_ij + (1 - alpha_ij) * f^(DG)_ij + # = f^(FV)_ij + (1 - alpha_ij) * f^(antidiffusive)_ij + @trixi_timeit timer() "blending factors" solver.volume_integral.limiter(u, semi, + solver, t, + dt) + + perform_idp_correction!(u, dt, mesh, equations, solver, cache) + + return nothing +end + +init_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing + +finalize_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing + +include("subcell_limiter_idp_correction_2d.jl") +end # @muladd diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl new file mode 100644 index 00000000000..f6b91444578 --- /dev/null +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -0,0 +1,44 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) + @unpack inverse_weights = dg.basis + @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients + + @threaded for element in eachelement(dg, cache) + # Sign switch as in apply_jacobian! + inverse_jacobian = -cache.elements.inverse_jacobian[element] + + for j in eachnode(dg), i in eachnode(dg) + # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} + alpha_flux1 = (1 - alpha1[i, j, element]) * + get_node_vars(antidiffusive_flux1, equations, dg, i, j, + element) + alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * + get_node_vars(antidiffusive_flux1, equations, dg, i + 1, + j, element) + alpha_flux2 = (1 - alpha2[i, j, element]) * + get_node_vars(antidiffusive_flux2, equations, dg, i, j, + element) + alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * + get_node_vars(antidiffusive_flux2, equations, dg, i, + j + 1, element) + + for v in eachvariable(equations) + u[v, i, j, element] += dt * inverse_jacobian * + (inverse_weights[i] * + (alpha_flux1_ip1[v] - alpha_flux1[v]) + + inverse_weights[j] * + (alpha_flux2_jp1[v] - alpha_flux2[v])) + end + end + end + + return nothing +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 90b2cd62191..570a25cece9 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -75,8 +75,14 @@ end @inline Base.ndims(::AbstractEquations{NDIMS}) where {NDIMS} = NDIMS -# equations act like scalars in broadcasting -Base.broadcastable(equations::AbstractEquations) = Ref(equations) +# Equations act like scalars in broadcasting. +# Using `Ref(equations)` would be more convenient in some circumstances. +# However, this does not work with Julia v1.9.3 correctly due to a (performance) +# bug in Julia, see +# - https://github.com/trixi-framework/Trixi.jl/pull/1618 +# - https://github.com/JuliaLang/julia/issues/51118 +# Thus, we use the workaround below. +Base.broadcastable(equations::AbstractEquations) = (equations,) """ flux(u, orientation_or_normal, equations) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 495e0ffc4a4..36bbc6de361 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -174,6 +174,46 @@ function Base.show(io::IO, ::MIME"text/plain", end end +""" + VolumeIntegralSubcellLimiting(limiter; + volume_flux_dg, volume_flux_fv) + +A subcell limiting volume integral type for DG methods based on subcell blending approaches +with a low-order FV method. Used with limiter [`SubcellLimiterIDP`](@ref). + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +struct VolumeIntegralSubcellLimiting{VolumeFluxDG, VolumeFluxFV, Limiter} <: + AbstractVolumeIntegral + volume_flux_dg::VolumeFluxDG + volume_flux_fv::VolumeFluxFV + limiter::Limiter +end + +function VolumeIntegralSubcellLimiting(limiter; volume_flux_dg, + volume_flux_fv) + VolumeIntegralSubcellLimiting{typeof(volume_flux_dg), typeof(volume_flux_fv), + typeof(limiter)}(volume_flux_dg, volume_flux_fv, + limiter) +end + +function Base.show(io::IO, mime::MIME"text/plain", + integral::VolumeIntegralSubcellLimiting) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + summary_header(io, "VolumeIntegralSubcellLimiting") + summary_line(io, "volume flux DG", integral.volume_flux_dg) + summary_line(io, "volume flux FV", integral.volume_flux_fv) + summary_line(io, "limiter", integral.limiter |> typeof |> nameof) + show(increment_indent(io), mime, integral.limiter) + summary_footer(io) + end +end + # TODO: FD. Should this definition live in a different file because it is # not strictly a DG method? """ diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index d80522d42fd..9148b936312 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -77,7 +77,7 @@ end eachelement(elements::ElementContainer2D) Return an iterator over the indices that specify the location in relevant data structures -for the elements in `elements`. +for the elements in `elements`. In particular, not the elements themselves are returned. """ @inline eachelement(elements::ElementContainer2D) = Base.OneTo(nelements(elements)) @@ -1254,4 +1254,138 @@ function init_mpi_mortars!(mpi_mortars, elements, mesh::TreeMesh2D) return nothing end + +# Container data structure (structure-of-arrays style) for FCT-type antidiffusive fluxes +# (i, j+1) +# | +# flux2(i, j+1) +# | +# (i-1, j) ---flux1(i, j)--- (i, j) ---flux1(i+1, j)--- (i+1, j) +# | +# flux2(i, j) +# | +# (i, j-1) +mutable struct ContainerAntidiffusiveFlux2D{uEltype <: Real} + antidiffusive_flux1::Array{uEltype, 4} # [variables, i, j, elements] + antidiffusive_flux2::Array{uEltype, 4} # [variables, i, j, elements] + # internal `resize!`able storage + _antidiffusive_flux1::Vector{uEltype} + _antidiffusive_flux2::Vector{uEltype} +end + +function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables, + n_nodes) where {uEltype <: Real} + nan_uEltype = convert(uEltype, NaN) + + # Initialize fields with defaults + _antidiffusive_flux1 = fill(nan_uEltype, + n_variables * (n_nodes + 1) * n_nodes * capacity) + antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), + (n_variables, n_nodes + 1, n_nodes, capacity)) + + _antidiffusive_flux2 = fill(nan_uEltype, + n_variables * n_nodes * (n_nodes + 1) * capacity) + antidiffusive_flux2 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), + (n_variables, n_nodes, n_nodes + 1, capacity)) + + return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1, + antidiffusive_flux2, + _antidiffusive_flux1, + _antidiffusive_flux2) +end + +nvariables(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1, 1) +nnodes(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1, 3) + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity) + n_nodes = nnodes(fluxes) + n_variables = nvariables(fluxes) + + @unpack _antidiffusive_flux1, _antidiffusive_flux2 = fluxes + + resize!(_antidiffusive_flux1, n_variables * (n_nodes + 1) * n_nodes * capacity) + fluxes.antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), + (n_variables, n_nodes + 1, n_nodes, + capacity)) + resize!(_antidiffusive_flux2, n_variables * n_nodes * (n_nodes + 1) * capacity) + fluxes.antidiffusive_flux2 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), + (n_variables, n_nodes, n_nodes + 1, + capacity)) + + return nothing +end + +# Container data structure (structure-of-arrays style) for variables used for IDP limiting +mutable struct ContainerSubcellLimiterIDP2D{uEltype <: Real} + alpha::Array{uEltype, 3} # [i, j, element] + alpha1::Array{uEltype, 3} + alpha2::Array{uEltype, 3} + variable_bounds::Vector{Array{uEltype, 3}} + # internal `resize!`able storage + _alpha::Vector{uEltype} + _alpha1::Vector{uEltype} + _alpha2::Vector{uEltype} + _variable_bounds::Vector{Vector{uEltype}} +end + +function ContainerSubcellLimiterIDP2D{uEltype}(capacity::Integer, n_nodes, + length) where {uEltype <: Real} + nan_uEltype = convert(uEltype, NaN) + + # Initialize fields with defaults + _alpha = fill(nan_uEltype, n_nodes * n_nodes * capacity) + alpha = unsafe_wrap(Array, pointer(_alpha), (n_nodes, n_nodes, capacity)) + _alpha1 = fill(nan_uEltype, (n_nodes + 1) * n_nodes * capacity) + alpha1 = unsafe_wrap(Array, pointer(_alpha1), (n_nodes + 1, n_nodes, capacity)) + _alpha2 = fill(nan_uEltype, n_nodes * (n_nodes + 1) * capacity) + alpha2 = unsafe_wrap(Array, pointer(_alpha2), (n_nodes, n_nodes + 1, capacity)) + + _variable_bounds = Vector{Vector{uEltype}}(undef, length) + variable_bounds = Vector{Array{uEltype, 3}}(undef, length) + for i in 1:length + _variable_bounds[i] = fill(nan_uEltype, n_nodes * n_nodes * capacity) + variable_bounds[i] = unsafe_wrap(Array, pointer(_variable_bounds[i]), + (n_nodes, n_nodes, capacity)) + end + + return ContainerSubcellLimiterIDP2D{uEltype}(alpha, alpha1, alpha2, + variable_bounds, + _alpha, _alpha1, _alpha2, + _variable_bounds) +end + +nnodes(container::ContainerSubcellLimiterIDP2D) = size(container.alpha, 1) + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(container::ContainerSubcellLimiterIDP2D, capacity) + n_nodes = nnodes(container) + + @unpack _alpha, _alpha1, _alpha2 = container + resize!(_alpha, n_nodes * n_nodes * capacity) + container.alpha = unsafe_wrap(Array, pointer(_alpha), (n_nodes, n_nodes, capacity)) + resize!(_alpha1, (n_nodes + 1) * n_nodes * capacity) + container.alpha1 = unsafe_wrap(Array, pointer(_alpha1), + (n_nodes + 1, n_nodes, capacity)) + resize!(_alpha2, n_nodes * (n_nodes + 1) * capacity) + container.alpha2 = unsafe_wrap(Array, pointer(_alpha2), + (n_nodes, n_nodes + 1, capacity)) + + @unpack _variable_bounds = container + for i in 1:length(_variable_bounds) + resize!(_variable_bounds[i], n_nodes * n_nodes * capacity) + container.variable_bounds[i] = unsafe_wrap(Array, pointer(_variable_bounds[i]), + (n_nodes, n_nodes, capacity)) + end + + return nothing +end end # @muladd diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index cb28dad968c..6e02bc1d94a 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -71,4 +71,9 @@ include("dg_3d_parabolic.jl") # as well as specialized implementations used to improve performance include("dg_2d_compressible_euler.jl") include("dg_3d_compressible_euler.jl") + +# Subcell limiters +include("subcell_limiters.jl") +include("subcell_limiters_2d.jl") +include("dg_2d_subcell_limiters.jl") end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl new file mode 100644 index 00000000000..70ff346740d --- /dev/null +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -0,0 +1,193 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function create_cache(mesh::TreeMesh{2}, equations, + volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) + cache = create_cache(mesh, equations, + VolumeIntegralPureLGLFiniteVolume(volume_integral.volume_flux_fv), + dg, uEltype) + + A3dp1_x = Array{uEltype, 3} + A3dp1_y = Array{uEltype, 3} + A3d = Array{uEltype, 3} + + fhat1_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, + nnodes(dg)) for _ in 1:Threads.nthreads()] + fhat2_threaded = A3dp1_y[A3dp1_y(undef, nvariables(equations), nnodes(dg), + nnodes(dg) + 1) for _ in 1:Threads.nthreads()] + flux_temp_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), nnodes(dg)) + for _ in 1:Threads.nthreads()] + + antidiffusive_fluxes = Trixi.ContainerAntidiffusiveFlux2D{uEltype}(0, + nvariables(equations), + nnodes(dg)) + + return (; cache..., antidiffusive_fluxes, fhat1_threaded, fhat2_threaded, + flux_temp_threaded) +end + +function calc_volume_integral!(du, u, + mesh::TreeMesh{2}, + nonconservative_terms, equations, + volume_integral::VolumeIntegralSubcellLimiting, + dg::DGSEM, cache) + @unpack limiter = volume_integral + + @threaded for element in eachelement(dg, cache) + subcell_limiting_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + volume_integral, limiter, + dg, cache) + end +end + +@inline function subcell_limiting_kernel!(du, u, + element, mesh::TreeMesh{2}, + nonconservative_terms::False, equations, + volume_integral, limiter::SubcellLimiterIDP, + dg::DGSEM, cache) + @unpack inverse_weights = dg.basis + @unpack volume_flux_dg, volume_flux_fv = volume_integral + + # high-order DG fluxes + @unpack fhat1_threaded, fhat2_threaded = cache + + fhat1 = fhat1_threaded[Threads.threadid()] + fhat2 = fhat2_threaded[Threads.threadid()] + calcflux_fhat!(fhat1, fhat2, u, mesh, + nonconservative_terms, equations, volume_flux_dg, dg, element, cache) + + # low-order FV fluxes + @unpack 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()] + fstar2_R = fstar2_R_threaded[Threads.threadid()] + calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, + nonconservative_terms, equations, volume_flux_fv, dg, element, cache) + + # antidiffusive flux + calcflux_antidiffusive!(fhat1, fhat2, fstar1_L, fstar2_L, u, mesh, + nonconservative_terms, equations, limiter, dg, element, + cache) + + # Calculate volume integral contribution of low-order FV flux + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + du[v, i, j, element] += inverse_weights[i] * + (fstar1_L[v, i + 1, j] - fstar1_R[v, i, j]) + + inverse_weights[j] * + (fstar2_L[v, i, j + 1] - fstar2_R[v, i, j]) + end + end + + return nothing +end + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**without non-conservative terms**). +# +# See also `flux_differencing_kernel!`. +@inline function calcflux_fhat!(fhat1, fhat2, u, + mesh::TreeMesh{2}, nonconservative_terms::False, + equations, + volume_flux, dg::DGSEM, element, cache) + @unpack weights, derivative_split = dg.basis + @unpack flux_temp_threaded = cache + + flux_temp = flux_temp_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of the `volume_flux` to save half of the possible two-point flux + # computations. + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, element) + flux1 = volume_flux(u_node, u_node_ii, 1, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], flux1, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, + equations, dg, ii, j) + end + end + + # FV-form flux `fhat` in x direction + fhat1[:, 1, :] .= zero(eltype(fhat1)) + fhat1[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1)) + + for j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations) + fhat1[v, i + 1, j] = fhat1[v, i, j] + weights[i] * flux_temp[v, i, j] + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, element) + flux2 = volume_flux(u_node, u_node_jj, 2, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], flux2, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, + equations, dg, i, jj) + end + end + + # FV-form flux `fhat` in y direction + fhat2[:, :, 1] .= zero(eltype(fhat2)) + fhat2[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2)) + + for j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations) + fhat2[v, i, j + 1] = fhat2[v, i, j] + weights[j] * flux_temp[v, i, j] + end + + return nothing +end + +# Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar`. +@inline function calcflux_antidiffusive!(fhat1, fhat2, fstar1, fstar2, u, mesh, + nonconservative_terms, equations, + limiter::SubcellLimiterIDP, dg, element, cache) + @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + + for j in eachnode(dg), i in 2:nnodes(dg) + for v in eachvariable(equations) + antidiffusive_flux1[v, i, j, element] = fhat1[v, i, j] - fstar1[v, i, j] + end + end + for j in 2:nnodes(dg), i in eachnode(dg) + for v in eachvariable(equations) + antidiffusive_flux2[v, i, j, element] = fhat2[v, i, j] - fstar2[v, i, j] + end + end + + antidiffusive_flux1[:, 1, :, element] .= zero(eltype(antidiffusive_flux1)) + antidiffusive_flux1[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1)) + + antidiffusive_flux2[:, :, 1, element] .= zero(eltype(antidiffusive_flux2)) + antidiffusive_flux2[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2)) + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl new file mode 100644 index 00000000000..3a707de3bc7 --- /dev/null +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -0,0 +1,103 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +abstract type AbstractSubcellLimiter end + +function create_cache(typ::Type{LimiterType}, + semi) where {LimiterType <: AbstractSubcellLimiter} + create_cache(typ, mesh_equations_solver_cache(semi)...) +end + +""" + SubcellLimiterIDP(equations::AbstractEquations, basis; + positivity_variables_cons = [], + positivity_correction_factor = 0.1) + +Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) +including: +- positivity limiting for conservative variables (`positivity_variables_cons`) + +The bounds are calculated using the low-order FV solution. The positivity limiter uses +`positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. + +!!! note + This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together. + Without the callback, no limiting takes place, leading to a standard flux-differencing DGSEM scheme. + +## References + +- Rueda-Ramírez, Pazner, Gassner (2022) + Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods + [DOI: 10.1016/j.compfluid.2022.105627](https://doi.org/10.1016/j.compfluid.2022.105627) +- Pazner (2020) + Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting + [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876) + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter + positivity::Bool + positivity_variables_cons::Vector{Int} # Positivity for conservative variables + positivity_correction_factor::RealT + cache::Cache +end + +# this method is used when the indicator is constructed as for shock-capturing volume integrals +function SubcellLimiterIDP(equations::AbstractEquations, basis; + positivity_variables_cons = [], + positivity_correction_factor = 0.1) + positivity = (length(positivity_variables_cons) > 0) + number_bounds = length(positivity_variables_cons) + + cache = create_cache(SubcellLimiterIDP, equations, basis, number_bounds) + + SubcellLimiterIDP{typeof(positivity_correction_factor), typeof(cache)}(positivity, + positivity_variables_cons, + positivity_correction_factor, + cache) +end + +function Base.show(io::IO, limiter::SubcellLimiterIDP) + @nospecialize limiter # reduce precompilation time + @unpack positivity = limiter + + print(io, "SubcellLimiterIDP(") + if !(positivity) + print(io, "No limiter selected => pure DG method") + else + print(io, "limiter=(") + positivity && print(io, "positivity") + print(io, "), ") + end + print(io, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) + @nospecialize limiter # reduce precompilation time + @unpack positivity = limiter + + if get(io, :compact, false) + show(io, limiter) + else + if !(positivity) + setup = ["limiter" => "No limiter selected => pure DG method"] + else + setup = ["limiter" => ""] + if positivity + string = "positivity with conservative variables $(limiter.positivity_variables_cons)" + setup = [setup..., "" => string] + setup = [ + setup..., + "" => " positivity correction factor = $(limiter.positivity_correction_factor)", + ] + end + end + summary_box(io, "SubcellLimiterIDP", setup) + end +end +end # @muladd diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl new file mode 100644 index 00000000000..09ab84ed11a --- /dev/null +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -0,0 +1,114 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# this method is used when the limiter is constructed as for shock-capturing volume integrals +function create_cache(indicator::Type{SubcellLimiterIDP}, + equations::AbstractEquations{2}, + basis::LobattoLegendreBasis, number_bounds) + subcell_limiter_coefficients = Trixi.ContainerSubcellLimiterIDP2D{real(basis) + }(0, + nnodes(basis), + number_bounds) + + cache = (; subcell_limiter_coefficients) + + return cache +end + +function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t, + dt; + kwargs...) + @unpack alpha = limiter.cache.subcell_limiter_coefficients + alpha .= zero(eltype(alpha)) + + if limiter.positivity + @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, + semi) + end + + # Calculate alpha1 and alpha2 + @unpack alpha1, alpha2 = limiter.cache.subcell_limiter_coefficients + @threaded for element in eachelement(dg, semi.cache) + for j in eachnode(dg), i in 2:nnodes(dg) + alpha1[i, j, element] = max(alpha[i - 1, j, element], alpha[i, j, element]) + end + for j in 2:nnodes(dg), i in eachnode(dg) + alpha2[i, j, element] = max(alpha[i, j - 1, element], alpha[i, j, element]) + end + alpha1[1, :, element] .= zero(eltype(alpha1)) + alpha1[nnodes(dg) + 1, :, element] .= zero(eltype(alpha1)) + alpha2[:, 1, element] .= zero(eltype(alpha2)) + alpha2[:, nnodes(dg) + 1, element] .= zero(eltype(alpha2)) + end + + return nothing +end + +@inline function idp_positivity!(alpha, limiter, u, dt, semi) + # Conservative variables + for (index, variable) in enumerate(limiter.positivity_variables_cons) + idp_positivity!(alpha, limiter, u, dt, semi, variable, index) + end + + return nothing +end + +@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable, index) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + @unpack inverse_weights = dg.basis + @unpack positivity_correction_factor = limiter + + @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients + + var_min = variable_bounds[index] + + @threaded for element in eachelement(dg, semi.cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + var = u[variable, i, j, element] + if var < 0 + error("Safe $variable is not safe. element=$element, node: $i $j, value=$var") + end + + # Compute bound + var_min[i, j, element] = positivity_correction_factor * var + + # Real one-sided Zalesak-type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2[variable, i, j + 1, element] + + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100) + + # Calculate alpha + alpha[i, j, element] = max(alpha[i, j, element], 1 - Qm) + end + end + + return nothing +end +end # @muladd diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl new file mode 100644 index 00000000000..8ecad69748b --- /dev/null +++ b/src/time_integration/methods_SSP.jl @@ -0,0 +1,241 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# Abstract base type for time integration schemes of explicit strong stability-preserving (SSP) +# Runge-Kutta (RK) methods. They are high-order time discretizations that guarantee the TVD property. +abstract type SimpleAlgorithmSSP end + +""" + SimpleSSPRK33(; stage_callbacks=()) + +The third-order SSP Runge-Kutta method of Shu and Osher. + +## References + +- Shu, Osher (1988) + "Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes" (Eq. 2.18) + [DOI: 10.1016/0021-9991(88)90177-5](https://doi.org/10.1016/0021-9991(88)90177-5) + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP + a::SVector{3, Float64} + b::SVector{3, Float64} + c::SVector{3, Float64} + stage_callbacks::StageCallbacks + + function SimpleSSPRK33(; stage_callbacks = ()) + a = SVector(0.0, 3 / 4, 1 / 3) + b = SVector(1.0, 1 / 4, 2 / 3) + c = SVector(0.0, 1.0, 1 / 2) + + # Butcher tableau + # c | a + # 0 | + # 1 | 1 + # 1/2 | 1/4 1/4 + # -------------------- + # b | 1/6 1/6 2/3 + + new{typeof(stage_callbacks)}(a, b, c, stage_callbacks) + end +end + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 +mutable struct SimpleIntegratorSSPOptions{Callback} + callback::Callback # callbacks; used in Trixi + adaptive::Bool # whether the algorithm is adaptive; ignored + dtmax::Float64 # ignored + maxiters::Int # maximal number of time steps + tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored +end + +function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) + SimpleIntegratorSSPOptions{typeof(callback)}(callback, false, Inf, maxiters, + [last(tspan)]) +end + +# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77 +# This implements the interface components described at +# https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1 +# which are used in Trixi. +mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, + SimpleIntegratorSSPOptions} + u::uType + du::uType + r0::uType + t::RealT + dt::RealT # current time step + dtcache::RealT # ignored + iter::Int # current number of time steps (iteration) + p::Params # will be the semidiscretization from Trixi + sol::Sol # faked + f::F + alg::Alg + opts::SimpleIntegratorSSPOptions + finalstep::Bool # added for convenience +end + +# Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) +function Base.getproperty(integrator::SimpleIntegratorSSP, field::Symbol) + if field === :stats + return (naccept = getfield(integrator, :iter),) + end + # general fallback + return getfield(integrator, field) +end + +""" + solve(ode, alg; dt, callbacks, kwargs...) + +The following structures and methods provide the infrastructure for SSP Runge-Kutta methods +of type `SimpleAlgorithmSSP`. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; + dt, callback = nothing, kwargs...) + u = copy(ode.u0) + du = similar(u) + r0 = similar(u) + t = first(ode.tspan) + iter = 0 + integrator = SimpleIntegratorSSP(u, du, r0, t, dt, zero(dt), iter, ode.p, + (prob = ode,), ode.f, alg, + SimpleIntegratorSSPOptions(callback, ode.tspan; + kwargs...), false) + + # resize container + resize!(integrator.p, nelements(integrator.p.solver, integrator.p.cache)) + + # initialize callbacks + if callback isa CallbackSet + for cb in callback.continuous_callbacks + error("unsupported") + end + for cb in callback.discrete_callbacks + cb.initialize(cb, integrator.u, integrator.t, integrator) + end + elseif !isnothing(callback) + error("unsupported") + end + + for stage_callback in alg.stage_callbacks + init_callback(stage_callback, integrator.p) + end + + solve!(integrator) +end + +function solve!(integrator::SimpleIntegratorSSP) + @unpack prob = integrator.sol + @unpack alg = integrator + t_end = last(prob.tspan) + callbacks = integrator.opts.callback + + integrator.finalstep = false + while !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end + + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end + + @. integrator.r0 = integrator.u + for stage in eachindex(alg.c) + t_stage = integrator.t + integrator.dt * alg.c[stage] + # compute du + integrator.f(integrator.du, integrator.u, integrator.p, t_stage) + + # perform forward Euler step + @. integrator.u = integrator.u + integrator.dt * integrator.du + + for stage_callback in alg.stage_callbacks + stage_callback(integrator.u, integrator, stage) + end + + # perform convex combination + @. integrator.u = alg.a[stage] * integrator.r0 + alg.b[stage] * integrator.u + end + + integrator.iter += 1 + integrator.t += integrator.dt + + # handle callbacks + if callbacks isa CallbackSet + for cb in callbacks.discrete_callbacks + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) + end + end + end + + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end + end + + for stage_callback in alg.stage_callbacks + finalize_callback(stage_callback, integrator.p) + end + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), prob) +end + +# get a cache where the RHS can be stored +get_du(integrator::SimpleIntegratorSSP) = integrator.du +get_tmp_cache(integrator::SimpleIntegratorSSP) = (integrator.r0,) + +# some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u +u_modified!(integrator::SimpleIntegratorSSP, ::Bool) = false + +# used by adaptive timestepping algorithms in DiffEq +function set_proposed_dt!(integrator::SimpleIntegratorSSP, dt) + integrator.dt = dt +end + +# stop the time integration +function terminate!(integrator::SimpleIntegratorSSP) + integrator.finalstep = true + empty!(integrator.opts.tstops) +end + +# used for AMR +function Base.resize!(integrator::SimpleIntegratorSSP, new_size) + resize!(integrator.u, new_size) + resize!(integrator.du, new_size) + resize!(integrator.r0, new_size) + + # Resize container + resize!(integrator.p, new_size) +end + +function Base.resize!(semi::AbstractSemidiscretization, new_size) + resize!(semi, semi.solver.volume_integral, new_size) +end + +Base.resize!(semi, volume_integral::AbstractVolumeIntegral, new_size) = nothing + +function Base.resize!(semi, volume_integral::VolumeIntegralSubcellLimiting, new_size) + # Resize container antidiffusive_fluxes + resize!(semi.cache.antidiffusive_fluxes, new_size) + + # Resize container subcell_limiter_coefficients + @unpack limiter = volume_integral + resize!(limiter.cache.subcell_limiter_coefficients, new_size) +end +end # @muladd diff --git a/src/time_integration/time_integration.jl b/src/time_integration/time_integration.jl index 539e00ff700..c1e53527121 100644 --- a/src/time_integration/time_integration.jl +++ b/src/time_integration/time_integration.jl @@ -15,4 +15,5 @@ end include("methods_2N.jl") include("methods_3Sstar.jl") +include("methods_SSP.jl") end # @muladd diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index cafa17edd4c..1e5aeac1786 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -102,7 +102,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), l2 = [0.17979210479598923, 1.2377495706611434, 6.289818963361573e-8], linf = [0.845938394800688, 3.3740800777086575, 4.4541473087633676e-7], - tspan = (0.0, 0.05)) + tspan = (0.0, 0.05), + atol = 3e-10) # see https://github.com/trixi-framework/Trixi.jl/issues/1617 end @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 6de380288db..e1e3ad32e7d 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -63,6 +63,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [0.18527440131928286, 0.2404798030563736, 0.23269573860381076, 0.6874012187446894]) end + @trixi_testset "elixir_euler_shockcapturing_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_subcell.jl"), + l2 = [0.08508147906199143, 0.04510299017724501, 0.045103019801950375, 0.6930704343869766], + linf = [0.31123546471463326, 0.5616274869594462, 0.5619692712224448, 2.88670199345138]) + end + @trixi_testset "elixir_euler_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"), l2 = [0.14170569763947993, 0.11647068900798814, 0.11647072556898294, 0.3391989213659599], diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 800dc31f84f..606afca1034 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -19,6 +19,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") tspan = (0.0, 0.001)) end + @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), + l2 = [81.52845664909304, 2.5455678559421346, 63229.190712645846, 0.19929478404550321, 0.011068604228443425], + linf = [249.21708417382013, 40.33299887640794, 174205.0118831558, 0.6881458768113586, 0.11274401158173972], + initial_refinement_level = 3, + tspan = (0.0, 0.001)) + end + @trixi_testset "elixir_eulermulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), l2 = [0.050182236154087095, 0.050189894464434635, 0.2258715597305131, 0.06175171559771687], diff --git a/test/test_trixi.jl b/test/test_trixi.jl index ddace6b4fbe..f2cd0cab94d 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -5,7 +5,7 @@ import Trixi # inside an elixir. """ @test_trixi_include(elixir; l2=nothing, linf=nothing, - atol=10*eps(), rtol=0.001, + atol=500*eps(), rtol=sqrt(eps()), parameters...) Test Trixi by calling `trixi_include(elixir; parameters...)`. diff --git a/test/test_unit.jl b/test/test_unit.jl index e70a9be6a4a..5c5291c2430 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -402,6 +402,9 @@ isdir(outdir) && rm(outdir, recursive=true) indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) + limiter_idp = SubcellLimiterIDP(true, [1], 0.1, "cache") + @test_nowarn show(stdout, limiter_idp) + # TODO: TrixiShallowWater: move unit test indicator_hg_swe = IndicatorHennemannGassnerShallowWater(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg_swe) @@ -637,7 +640,7 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for HLL flux (naive): SWE" begin @@ -674,7 +677,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] - orientations = [1, 2] + orientations = [1, 2] u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2),] @@ -704,7 +707,7 @@ isdir(outdir) && rm(outdir, recursive=true) for u in u_values, normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: CEE" begin @@ -718,7 +721,7 @@ isdir(outdir) && rm(outdir, recursive=true) for orientation in orientations @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations) end - + equations = CompressibleEulerEquations2D(1.4) u = SVector(1.1, -0.5, 2.34, 5.5) @@ -734,7 +737,7 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end equations = CompressibleEulerEquations3D(1.4) u = SVector(1.1, -0.5, 2.34, 2.4, 5.5) @@ -752,7 +755,7 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin @@ -815,7 +818,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] - orientations = [1, 2] + orientations = [1, 2] u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2),] @@ -845,7 +848,7 @@ isdir(outdir) && rm(outdir, recursive=true) for u in u_values, normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for HLLE flux: CEE" begin @@ -873,7 +876,7 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end equations = CompressibleEulerEquations3D(1.4) u = SVector(1.1, -0.5, 2.34, 2.4, 5.5) @@ -891,7 +894,7 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for HLLE flux: SWE" begin @@ -907,7 +910,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] - orientations = [1, 2] + orientations = [1, 2] u = SVector(1, 0.5, 0.5, 0.0) @@ -937,7 +940,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] - orientations = [1, 2] + orientations = [1, 2] u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2),] @@ -956,7 +959,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(0.0, 0.0, 1.0), SVector(0.5, -0.5, 0.2), SVector(-1.2, 0.3, 1.4)] - orientations = [1, 2, 3] + orientations = [1, 2, 3] u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2),] @@ -967,7 +970,7 @@ isdir(outdir) && rm(outdir, recursive=true) for u in u_values, normal_direction in normal_directions @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) - end + end end @timed_testset "Consistency check for Godunov flux" begin @@ -1137,7 +1140,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(-1.2, 0.3)] u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0),] - fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, + fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, flux_hll, FluxHLL(min_max_speed_davis)] for f_std in fluxes @@ -1157,7 +1160,7 @@ isdir(outdir) && rm(outdir, recursive=true) SVector(-1.2, 0.3, 1.4)] u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0),] - fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), + fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis)] for f_std in fluxes @@ -1173,11 +1176,11 @@ isdir(outdir) && rm(outdir, recursive=true) normal_directions = [SVector(1.0, 0.0), SVector(0.0, 1.0), SVector(0.5, -0.5), - SVector(-1.2, 0.3)] + SVector(-1.2, 0.3)] u = SVector(1, 0.5, 0.5, 0.0) - fluxes = [flux_central, flux_fjordholm_etal, flux_wintermeyer_etal, + fluxes = [flux_central, flux_fjordholm_etal, flux_wintermeyer_etal, flux_hll, FluxHLL(min_max_speed_davis), FluxHLL(min_max_speed_einfeldt)] end