From 87470da113988698f84de2693ad7db42d1fbe10c Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 28 Oct 2024 11:28:11 +0100 Subject: [PATCH] Small modifications --- .JuliaFormatter.toml | 8 + .github/dependabot.yml | 7 + .github/workflows/FormatCheck.yml | 44 +++ .github/workflows/SpellCheck.yml | 13 + .github/workflows/ci.yml | 63 +++ .gitignore | 33 ++ .typos.toml | 4 + LICENSE.md | 21 + Project.toml | 15 + README.md | 47 ++- ...blast_wave_neuralnetwork_perssonperaire.jl | 115 ++++++ ...r_blast_wave_neuralnetwork_rayhesthaven.jl | 115 ++++++ ...ixir_euler_blast_wave_neuralnetwork_cnn.jl | 114 ++++++ ...blast_wave_neuralnetwork_perssonperaire.jl | 113 ++++++ ...r_blast_wave_neuralnetwork_rayhesthaven.jl | 116 ++++++ ...bility_amr_neuralnetwork_perssonperaire.jl | 140 +++++++ ...blast_wave_neuralnetwork_perssonperaire.jl | 132 +++++++ src/TrixiSmartShockFinder.jl | 19 + src/indicators.jl | 202 ++++++++++ src/indicators_1d.jl | 225 +++++++++++ src/indicators_2d.jl | 358 ++++++++++++++++++ test/Project.toml | 14 + test/runtests.jl | 7 + test/test_tree_1d_euler.jl | 40 ++ test/test_tree_2d_euler.jl | 134 +++++++ test/test_unit.jl | 34 ++ 26 files changed, 2132 insertions(+), 1 deletion(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/FormatCheck.yml create mode 100644 .github/workflows/SpellCheck.yml create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 .typos.toml create mode 100644 LICENSE.md create mode 100644 Project.toml create mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl create mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl create mode 100644 src/TrixiSmartShockFinder.jl create mode 100644 src/indicators.jl create mode 100644 src/indicators_1d.jl create mode 100644 src/indicators_2d.jl create mode 100644 test/Project.toml create mode 100644 test/runtests.jl create mode 100644 test/test_tree_1d_euler.jl create mode 100644 test/test_tree_2d_euler.jl create mode 100644 test/test_unit.jl diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..8518d20 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,8 @@ +# Use SciML style: https://github.com/SciML/SciMLStyle +style = "sciml" + +# Python style alignment. See https://github.com/domluna/JuliaFormatter.jl/pull/732. +yas_style_nesting = true + +# Align struct fields for better readability of large struct definitions +align_struct_field = true diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..d60f070 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "monthly" diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 0000000..fd91472 --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,44 @@ +name: format-check + +on: + push: + branches: + - 'main' + tags: '*' + pull_request: + +jobs: + check-format: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v4 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version = "0.13.0"))' + # + # TODO: Change the call below to + # format(".") + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(["examples", "src", "test"])' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml new file mode 100644 index 0000000..fb71fe4 --- /dev/null +++ b/.github/workflows/SpellCheck.yml @@ -0,0 +1,13 @@ +name: Spell Check + +on: [pull_request, workflow_dispatch] + +jobs: + typos-check: + name: Spell Check with Typos + runs-on: ubuntu-latest + steps: + - name: Checkout Actions Repository + uses: actions/checkout@v4 + - name: Check spelling + uses: crate-ci/typos@v1.16.21 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..4b5f38e --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,63 @@ +name: CI + +on: + push: + branches: + - main + paths-ignore: + - 'AUTHORS.md' + - 'LICENSE.md' + - 'README.md' + pull_request: + paths-ignore: + - 'AUTHORS.md' + - 'LICENSE.md' + - 'README.md' + workflow_dispatch: + +# Cancel redundant CI tests automatically +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + test: + if: "!contains(github.event.head_commit.message, 'skip ci')" + # We could also include the Julia version as in + # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }} + # to be more specific. However, that requires us updating the required CI tests whenever we update Julia. + name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.9' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)' + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@v1 + env: + PYTHON: "" + - name: Run tests with coverage + uses: julia-actions/julia-runtest@v1 + with: + coverage: true + - uses: julia-actions/julia-processcoverage@v1 + with: + directories: src,examples + - uses: coverallsapp/github-action@master + with: + github-token: ${{ secrets.GITHUB_TOKEN }} + flag-name: run-${{ matrix.trixi_test }}-${{ matrix.os }}-${{ matrix.version }}-${{ matrix.arch }}-${{ github.run_id }} + parallel: false + path-to-lcov: ./lcov.info diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3132b9a --- /dev/null +++ b/.gitignore @@ -0,0 +1,33 @@ +*.pdf +*.png +*.h5 +*.mp4 +*.mem +*.vtu +*.pvd +*.avi +*.ogv +*.mesh +*.bson +*.inp +**/Manifest.toml +out*/ +docs/build +docs/src/notebooks +docs/src/authors.md +docs/src/tutorials +public/ +coverage/ +coverage_report/ +**/*.jl.*.cov +**/.ipynb_checkpoints/ +**/solution_*.txt +**/restart_*.txt +.vscode/ + +.DS_Store + +run +run/* + +LocalPreferences.toml diff --git a/.typos.toml b/.typos.toml new file mode 100644 index 0000000..1b16de9 --- /dev/null +++ b/.typos.toml @@ -0,0 +1,4 @@ +[default.extend-words] +rcall = "rcall" +claus = "claus" +dum = "dum" diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..3dd3458 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023-present Michael Schlottke-Lakemper, Julia Odenthal + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..f4ddb2b --- /dev/null +++ b/Project.toml @@ -0,0 +1,15 @@ +name = "TrixiSmartShockFinder" +uuid = "6b42f25e-8cc0-41c2-936f-439385811e39" +authors = ["Michael Schlottke-Lakemper and contributors"] +version = "0.1.0" + +[deps] +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" + +[compat] +Flux = "0.13.15, 0.14" +MuladdMacro = "0.2.4" +# Trixi = "0.6" +julia = "1.9" diff --git a/README.md b/README.md index d5232c8..2951284 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,47 @@ # TrixiSmartShockFinder.jl -Spin-off repository of Trixi.jl with neural network-based shock indicators +[![Build Status](https://github.com/trixi-framework/TrixiSmartShockFinder.jl/workflows/CI/badge.svg)](https://github.com/trixi-framework/TrixiSmartShockFinder.jl/actions?query=workflow%3ACI) +[![Coveralls](https://coveralls.io/repos/github/trixi-framework/TrixiSmartShockFinder.jl/badge.svg?branch=main)](https://coveralls.io/github/trixi-framework/TrixiSmartShockFinder.jl?branch=main) +[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) + +Spin-off repository of [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) +with neural network-based shock indicators. + +**Note: This repository is currently not under development and has been archived. If you are +interested in using any of this, please get in touch with the developers of the +[Trixi.jl](https://github.com/trixi-framework/Trixi.jl) package. + + +## Usage +To run any of the elixirs with the neurl network-based indicators, you first need to install +all required auxiliary packages by running the following code in the Julia REPL: +```julia +julia> using Pkg + +julia> Pkg.add(["BSON", "Flux", "OrdinaryDiffEq", "Trixi"]) +``` + +Then, clone this repository +```shell +git clone git@github.com:trixi-framework/TrixiSmartShockFinder.jl.git +``` +enter the directory, and start Julia (tested with Julia v1.9) with the project set to +the clone directory +``` +cd TrixiSmartShockFinder.jl +julia --project=. +``` + +Now run one of the elixirs in the `examples` folder, e.g., +```julia +using TrixiSmartShockFinder +trixi_include("examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl") +``` + +## Authors +TrixiSmartShockFinder.jl was initiated by +[Michael Schlottke-Lakemper](https://lakemper.eu) (University of Augsburg, Germany) and +Julia Odenthal (University of Cologne, Germany). + + +## License and contributing +TrixiSmartShockFinder.jl is licensed under the MIT license (see [LICENSE.md](LICENSE.md)). diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl new file mode 100644 index 0000000..bf30825 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl @@ -0,0 +1,115 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelnnpp-0.97-0.0001.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.97-0.0001.bson", + network) +model1d = load(network, @__MODULE__)[:model1d] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations1D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) + +A medium blast wave 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::CompressibleEulerEquations1D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0) + x_norm = x[1] - inicenter[1] + r = abs(x_norm) + # The following code is equivalent to + # phi = atan(0.0, x_norm) + # cos_phi = cos(phi) + # in 1D but faster + cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = false, + alpha_amr = false, + variable = density_pressure, + network = model1d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.0,) +coordinates_max = (2.0,) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl new file mode 100644 index 0000000..94cc8ef --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl @@ -0,0 +1,115 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelnnrh-0.95-0.009.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnrh-0.95-0.009.bson", + network) +model1d = load(network, @__MODULE__)[:model1d] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations1D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) + +A medium blast wave 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::CompressibleEulerEquations1D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0) + x_norm = x[1] - inicenter[1] + r = abs(x_norm) + # The following code is equivalent to + # phi = atan(0.0, x_norm) + # cos_phi = cos(phi) + # in 1D but faster + cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkRayHesthaven(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = false, + alpha_amr = false, + variable = density_pressure, + network = model1d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.0,) +coordinates_max = (2.0,) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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_euler_blast_wave_neuralnetwork_cnn.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl new file mode 100644 index 0000000..5432688 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl @@ -0,0 +1,114 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelcnn-0.964-0.001.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelcnn-0.964-0.001.bson", + network) +model2dcnn = load(network, @__MODULE__)[:model2dcnn] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +NeuralNetworkCNN = TrixiSmartShockFinder.NeuralNetworkCNN +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave 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) -> "medium blast wave" + # 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 + rho = r > 0.5 ? 1.0 : 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-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_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkCNN(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable = density_pressure, + network = model2dcnn) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + 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 = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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.9) + +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_euler_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl new file mode 100644 index 0000000..33270b9 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl @@ -0,0 +1,113 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", + network) +model2d = load(network, @__MODULE__)[:model2d] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave 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) -> "medium blast wave" + # 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 + rho = r > 0.5 ? 1.0 : 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-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_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable = density_pressure, + network = model2d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + 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 = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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.9) + +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_euler_blast_wave_neuralnetwork_rayhesthaven.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl new file mode 100644 index 0000000..a37aee4 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl @@ -0,0 +1,116 @@ +using Downloads: download +using Flux +using Random +using BSON: load +network = joinpath(@__DIR__, "modelnnrhs-0.973-0.001.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnrhs-0.973-0.001.bson", + network) +model2d = load(network, @__MODULE__)[:model2d] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave 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) -> "medium blast wave" + # 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 + rho = r > 0.5 ? 1.0 : 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-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_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkRayHesthaven(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable = density_pressure, + network = model2d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + 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) +refinement_patches = () # To allow for specifying them via `trixi_include` +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + refinement_patches = refinement_patches, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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.9) + +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_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl new file mode 100644 index 0000000..61abaa7 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl @@ -0,0 +1,140 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", + network) +model2d = load(network, @__MODULE__)[:model2d] + +using Random: seed! +seed!(0) + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test. +""" +function initial_condition_kelvin_helmholtz_instability(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-0.5,0.5]^2 + dens0 = 1.0 # outside density + dens1 = 2.0 # inside density + velx0 = -0.5 # outside velocity + velx1 = 0.5 # inside velocity + slope = 50 # used for tanh instead of discontinuous initial condition + # pressure equilibrium + p = 2.5 + # y velocity v2 is only white noise + v2 = 0.01 * (rand(Float64, 1)[1] - 0.5) + # density + rho = dens0 + + (dens1 - dens0) * 0.5 * + (1 + (tanh(slope * (x[2] + 0.25)) - (tanh(slope * (x[2] - 0.25)) + 1))) + # x velocity is also augmented with noise + v1 = velx0 + + (velx1 - velx0) * 0.5 * + (1 + (tanh(slope * (x[2] + 0.25)) - (tanh(slope * (x[2] - 0.25)) + 1))) + + 0.01 * (rand(Float64, 1)[1] - 0.5) + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_kelvin_helmholtz_instability + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 0.002, + alpha_min = 0.0001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable = density_pressure, + network = model2d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-0.5, -0.5) +coordinates_max = (0.5, 0.5) +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, 5.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) + +amr_indicator = IndicatorNeuralNetwork(semi, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 1.0, + alpha_min = 0.0001, + alpha_smooth = false, + alpha_continuous = true, + alpha_amr = true, + variable = density_pressure, + network = model2d) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 4, + med_level = 6, med_threshold = 0.3, # med_level = current level + max_level = 7, max_threshold = 0.5) +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 1.3) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + amr_callback, 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_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl new file mode 100644 index 0000000..73a57ff --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl @@ -0,0 +1,132 @@ +using Downloads: download +using Flux +using BSON: load +network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") +download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", + network) +model2d = load(network, @__MODULE__)[:model2d] + +using OrdinaryDiffEq +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# This elixir was one of the setups used in the following master thesis: +# - Julia Odenthal (2021) +# Shock capturing with artificial neural networks +# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper +# This motivates the particular choice of fluxes, mesh resolution etc. + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +The Sedov blast wave setup based on Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 +""" +function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # 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) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + # r0 = 0.5 # = more reasonable setup + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + p0_outer = 1.0e-5 # = true Sedov setup + # p0_outer = 1.0e-3 # = more reasonable setup + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_sedov_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable = density_pressure, + network = model2d) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + 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 = 6, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +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) + +amr_indicator = IndicatorNeuralNetwork(semi, + indicator_type = NeuralNetworkPerssonPeraire(), + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = true, + variable = density_pressure, + network = model2d) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 4, + max_level = 6, max_threshold = 0.22) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + amr_callback, 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/src/TrixiSmartShockFinder.jl b/src/TrixiSmartShockFinder.jl new file mode 100644 index 0000000..0b2dfa2 --- /dev/null +++ b/src/TrixiSmartShockFinder.jl @@ -0,0 +1,19 @@ +module TrixiSmartShockFinder + +using MuladdMacro: @muladd +using Flux: params +using Trixi +using Trixi: AbstractIndicator, AbstractEquations, AbstractSemidiscretization, @threaded, + summary_box, eachdirection, get_node_vars, multiply_scalar_dimensionwise!, + apply_smoothing!, gauss_lobatto_nodes_weights, polynomial_interpolation_matrix, + has_any_neighbor, has_neighbor, has_children, trixi_include + +include("indicators.jl") +include("indicators_1d.jl") +include("indicators_2d.jl") + +export IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven, + NeuralNetworkCNN +export trixi_include + +end # module TrixiSmartShockFinder diff --git a/src/indicators.jl b/src/indicators.jl new file mode 100644 index 0000000..7784249 --- /dev/null +++ b/src/indicators.jl @@ -0,0 +1,202 @@ +# 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 + +""" + IndicatorNeuralNetwork + +Artificial neural network based indicator used for shock-capturing or AMR. +Depending on the indicator_type, different input values and corresponding trained networks are used. + +`indicator_type = NeuralNetworkPerssonPeraire()` +- Input: The energies in lower modes as well as nnodes(dg). + +`indicator_type = NeuralNetworkRayHesthaven()` +- 1d Input: Cell average of the cell and its neighboring cells as well as the interface values. +- 2d Input: Linear modal values of the cell and its neighboring cells. + +- Ray, Hesthaven (2018) + "An artificial neural network as a troubled-cell indicator" + [doi:10.1016/j.jcp.2018.04.029](https://doi.org/10.1016/j.jcp.2018.04.029) +- Ray, Hesthaven (2019) + "Detecting troubled-cells on two-dimensional unstructured grids using a neural network" + [doi:10.1016/j.jcp.2019.07.043](https://doi.org/10.1016/j.jcp.2019.07.043) + +`indicator_type = CNN (Only in 2d)` +- Based on convolutional neural network. +- 2d Input: Interpolation of the nodal values of the `indicator.variable` to the 4x4 LGL nodes. + +If `alpha_continuous == true` the continuous network output for troubled cells (`alpha > 0.5`) is considered. +If the cells are good (`alpha < 0.5`), `alpha` is set to `0`. +If `alpha_continuous == false`, the blending factor is set to `alpha = 0` for good cells and +`alpha = 1` for troubled cells. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +""" +struct IndicatorNeuralNetwork{IndicatorType, RealT <: Real, Variable, Chain, Cache} <: + AbstractIndicator + indicator_type::IndicatorType + alpha_max::RealT + alpha_min::RealT + alpha_smooth::Bool + alpha_continuous::Bool + alpha_amr::Bool + variable::Variable + network::Chain + cache::Cache +end + +# this method is used when the indicator is constructed as for shock-capturing volume integrals +function IndicatorNeuralNetwork(equations::AbstractEquations, basis; + indicator_type, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = false, + variable, + network) + alpha_max, alpha_min = promote(alpha_max, alpha_min) + IndicatorType = typeof(indicator_type) + cache = create_cache(IndicatorNeuralNetwork{IndicatorType}, equations, basis) + IndicatorNeuralNetwork{IndicatorType, typeof(alpha_max), typeof(variable), + typeof(network), typeof(cache)}(indicator_type, alpha_max, + alpha_min, alpha_smooth, + alpha_continuous, alpha_amr, + variable, + network, cache) +end + +# this method is used when the indicator is constructed as for AMR +function IndicatorNeuralNetwork(semi::AbstractSemidiscretization; + indicator_type, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + alpha_continuous = true, + alpha_amr = true, + variable, + network) + alpha_max, alpha_min = promote(alpha_max, alpha_min) + IndicatorType = typeof(indicator_type) + cache = create_cache(IndicatorNeuralNetwork{IndicatorType}, semi) + IndicatorNeuralNetwork{IndicatorType, typeof(alpha_max), typeof(variable), + typeof(network), typeof(cache)}(indicator_type, alpha_max, + alpha_min, alpha_smooth, + alpha_continuous, alpha_amr, + variable, + network, cache) +end + +function Base.show(io::IO, indicator::IndicatorNeuralNetwork) + @nospecialize indicator # reduce precompilation time + + print(io, "IndicatorNeuralNetwork(") + print(io, indicator.indicator_type) + print(io, ", alpha_max=", indicator.alpha_max) + print(io, ", alpha_min=", indicator.alpha_min) + print(io, ", alpha_smooth=", indicator.alpha_smooth) + print(io, ", alpha_continuous=", indicator.alpha_continuous) + print(io, indicator.variable) + print(io, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorNeuralNetwork) + @nospecialize indicator # reduce precompilation time + + if get(io, :compact, false) + show(io, indicator) + else + setup = [ + "indicator type" => indicator.indicator_type, + "max. α" => indicator.alpha_max, + "min. α" => indicator.alpha_min, + "smooth α" => (indicator.alpha_smooth ? "yes" : "no"), + "continuous α" => (indicator.alpha_continuous ? "yes" : "no"), + "indicator variable" => indicator.variable, + ] + summary_box(io, "IndicatorNeuralNetwork", setup) + end +end + +# Convert probability for troubled cell to indicator value for shockcapturing/AMR +@inline function probability_to_indicator(probability_troubled_cell, alpha_continuous, + alpha_amr, + alpha_min, alpha_max) + # Initialize indicator to zero + alpha_element = zero(probability_troubled_cell) + + if alpha_continuous && !alpha_amr + # Set good cells to 0 and troubled cells to continuous value of the network prediction + if probability_troubled_cell > 0.5 + alpha_element = probability_troubled_cell + else + alpha_element = zero(probability_troubled_cell) + end + + # Take care of the case close to pure FV + if alpha_element > 1 - alpha_min + alpha_element = one(alpha_element) + end + + # Scale the probability for a troubled cell (in [0,1]) to the maximum allowed alpha + alpha_element *= alpha_max + elseif !alpha_continuous && !alpha_amr + # Set good cells to 0 and troubled cells to 1 + if probability_troubled_cell > 0.5 + alpha_element = alpha_max + else + alpha_element = zero(alpha_max) + end + elseif alpha_amr + # The entire continuous output of the neural network is used for AMR + alpha_element = probability_troubled_cell + + # Scale the probability for a troubled cell (in [0,1]) to the maximum allowed alpha + alpha_element *= alpha_max + end + + return alpha_element +end + +""" + NeuralNetworkPerssonPeraire + +Indicator type for creating an `IndicatorNeuralNetwork` indicator. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`IndicatorNeuralNetwork`](@ref) +""" +struct NeuralNetworkPerssonPeraire end + +""" + NeuralNetworkRayHesthaven + +Indicator type for creating an `IndicatorNeuralNetwork` indicator. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`IndicatorNeuralNetwork`](@ref) +""" +struct NeuralNetworkRayHesthaven end + +""" + NeuralNetworkCNN + +Indicator type for creating an `IndicatorNeuralNetwork` indicator. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +See also: [`IndicatorNeuralNetwork`](@ref) +""" +struct NeuralNetworkCNN end +end # @muladd diff --git a/src/indicators_1d.jl b/src/indicators_1d.jl new file mode 100644 index 0000000..ee9ff54 --- /dev/null +++ b/src/indicators_1d.jl @@ -0,0 +1,225 @@ +# 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 indicator is constructed as for shock-capturing volume integrals +# empty cache is default +function create_cache(::Type{<:IndicatorNeuralNetwork}, + equations::AbstractEquations{1}, basis::LobattoLegendreBasis) + return NamedTuple() +end + +# cache for NeuralNetworkPerssonPeraire-type indicator +function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire}}, + equations::AbstractEquations{1}, basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + alpha_tmp = similar(alpha) + A = Array{real(basis), ndims(equations)} + + prototype = A(undef, nnodes(basis)) + indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + + return (; alpha, alpha_tmp, indicator_threaded, modal_threaded) +end + +# cache for NeuralNetworkRayHesthaven-type indicator +function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkRayHesthaven}}, + equations::AbstractEquations{1}, basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + alpha_tmp = similar(alpha) + A = Array{real(basis), ndims(equations)} + + prototype = A(undef, nnodes(basis)) + indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + neighbor_ids = Vector{Int}(undef, 2) + + return (; alpha, alpha_tmp, indicator_threaded, neighbor_ids) +end + +# this method is used when the indicator is constructed as for AMR +function create_cache(typ::Type{<:IndicatorNeuralNetwork}, + mesh, equations::AbstractEquations{1}, dg::DGSEM, cache) + create_cache(typ, equations, dg.basis) +end + +function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire})(u::AbstractArray{ + <:Any, + 3 + }, + mesh, + equations, + dg::DGSEM, + cache; + kwargs...) + @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann + + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_ann.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + + # Calculate indicator variables at Gauss-Lobatto nodes + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, element) + indicator[i] = indicator_ann.variable(u_local, equations) + end + + # Convert to modal representation + multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, + indicator) + + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for i in 1:nnodes(dg) + total_energy += modal[i]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 1) + total_energy_clip1 += modal[i]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 2) + total_energy_clip2 += modal[i]^2 + end + + # Calculate energy in highest modes + X1 = (total_energy - total_energy_clip1) / total_energy + X2 = (total_energy_clip1 - total_energy_clip2) / total_energy_clip1 + + # There are two versions of the network: + # The first one only takes the highest energy modes as input, the second one also the number of + # nodes. Automatically use the correct input by checking the number of inputs of the network. + if size(params(network)[1], 2) == 2 + network_input = SVector(X1, X2) + elseif size(params(network)[1], 2) == 3 + network_input = SVector(X1, X2, nnodes(dg)) + end + + # Scale input data + network_input = network_input / + max(maximum(abs, network_input), one(eltype(network_input))) + probability_troubled_cell = network(network_input)[1] + + # Compute indicator value + alpha[element] = probability_to_indicator(probability_troubled_cell, + alpha_continuous, + alpha_amr, alpha_min, alpha_max) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end + +function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})(u::AbstractArray{ + <:Any, + 3 + }, + mesh, + equations, + dg::DGSEM, + cache; + kwargs...) + @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann + + @unpack alpha, alpha_tmp, indicator_threaded, neighbor_ids = indicator_ann.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + c2e = zeros(Int, length(mesh.tree)) + for element in eachelement(dg, cache) + c2e[cache.elements.cell_ids[element]] = element + end + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + cell_id = cache.elements.cell_ids[element] + + for direction in eachdirection(mesh.tree) + if !has_any_neighbor(mesh.tree, cell_id, direction) + neighbor_ids[direction] = element + continue + end + if has_neighbor(mesh.tree, cell_id, direction) + neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id] + if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor + if direction == 1 + neighbor_ids[direction] = c2e[mesh.tree.child_ids[2, + neighbor_cell_id]] + else + neighbor_ids[direction] = c2e[mesh.tree.child_ids[1, + neighbor_cell_id]] + end + else # Cell has same refinement level neighbor + neighbor_ids[direction] = c2e[neighbor_cell_id] + end + else # Cell is small and has large neighbor + parent_id = mesh.tree.parent_ids[cell_id] + neighbor_cell_id = mesh.tree.neighbor_ids[direction, parent_id] + neighbor_ids[direction] = c2e[neighbor_cell_id] + end + end + + # Calculate indicator variables at Gauss-Lobatto nodes + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, element) + indicator[i] = indicator_ann.variable(u_local, equations) + end + + # Cell average and interface values of the cell + X2 = sum(indicator) / nnodes(dg) + X4 = indicator[1] + X5 = indicator[end] + + # Calculate indicator variables from left neighboring cell at Gauss-Lobatto nodes + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, neighbor_ids[1]) + indicator[i] = indicator_ann.variable(u_local, equations) + end + X1 = sum(indicator) / nnodes(dg) + + # Calculate indicator variables from right neighboring cell at Gauss-Lobatto nodes + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, neighbor_ids[2]) + indicator[i] = indicator_ann.variable(u_local, equations) + end + X3 = sum(indicator) / nnodes(dg) + network_input = SVector(X1, X2, X3, X4, X5) + + # Scale input data + network_input = network_input / + max(maximum(abs, network_input), one(eltype(network_input))) + probability_troubled_cell = network(network_input)[1] + + # Compute indicator value + alpha[element] = probability_to_indicator(probability_troubled_cell, + alpha_continuous, + alpha_amr, alpha_min, alpha_max) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end +end # @muladd diff --git a/src/indicators_2d.jl b/src/indicators_2d.jl new file mode 100644 index 0000000..42e2ad7 --- /dev/null +++ b/src/indicators_2d.jl @@ -0,0 +1,358 @@ +# 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 indicator is constructed as for shock-capturing volume integrals +# empty cache is default +function create_cache(::Type{IndicatorNeuralNetwork}, + equations::AbstractEquations{2}, basis::LobattoLegendreBasis) + return NamedTuple() +end + +# cache for NeuralNetworkPerssonPeraire-type indicator +function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire}}, + equations::AbstractEquations{2}, basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + alpha_tmp = similar(alpha) + A = Array{real(basis), ndims(equations)} + + @assert nnodes(basis)>=4 "Indicator only works for nnodes >= 4 (polydeg > 2)" + + prototype = A(undef, nnodes(basis), nnodes(basis)) + indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + modal_tmp1_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + + return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded) +end + +# cache for NeuralNetworkRayHesthaven-type indicator +function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkRayHesthaven}}, + equations::AbstractEquations{2}, basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + alpha_tmp = similar(alpha) + A = Array{real(basis), ndims(equations)} + + prototype = A(undef, nnodes(basis), nnodes(basis)) + indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + modal_tmp1_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + + network_input = Vector{Float64}(undef, 15) + neighbor_ids = Array{Int64}(undef, 8) + neighbor_mean = Array{Float64}(undef, 4, 3) + + return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded, + network_input, neighbor_ids, neighbor_mean) +end + +# cache for NeuralNetworkCNN-type indicator +function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkCNN}}, + equations::AbstractEquations{2}, basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + alpha_tmp = similar(alpha) + A = Array{real(basis), ndims(equations)} + + prototype = A(undef, nnodes(basis), nnodes(basis)) + indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + n_cnn = 4 + nodes, _ = gauss_lobatto_nodes_weights(nnodes(basis)) + cnn_nodes, _ = gauss_lobatto_nodes_weights(n_cnn) + vandermonde = polynomial_interpolation_matrix(nodes, cnn_nodes) + network_input = Array{Float32}(undef, n_cnn, n_cnn, 1, 1) + + return (; alpha, alpha_tmp, indicator_threaded, nodes, cnn_nodes, vandermonde, + network_input) +end + +# this method is used when the indicator is constructed as for AMR +function create_cache(typ::Type{<:IndicatorNeuralNetwork}, + mesh, equations::AbstractEquations{2}, dg::DGSEM, cache) + create_cache(typ, equations, dg.basis) +end + +function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire})(u, + mesh::TreeMesh{ + 2 + }, + equations, + dg::DGSEM, + cache; + kwargs...) + @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann + + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded = indicator_ann.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] + + # Calculate indicator variables at Gauss-Lobatto nodes + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + indicator[i, j] = indicator_ann.variable(u_local, equations) + end + + # Convert to modal representation + multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, + indicator, modal_tmp1) + + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for j in 1:nnodes(dg), i in 1:nnodes(dg) + total_energy += modal[i, j]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for j in 1:(nnodes(dg) - 1), i in 1:(nnodes(dg) - 1) + total_energy_clip1 += modal[i, j]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for j in 1:(nnodes(dg) - 2), i in 1:(nnodes(dg) - 2) + total_energy_clip2 += modal[i, j]^2 + end + total_energy_clip3 = zero(eltype(modal)) + for j in 1:(nnodes(dg) - 3), i in 1:(nnodes(dg) - 3) + total_energy_clip3 += modal[i, j]^2 + end + + # Calculate energy in higher modes and polynomial degree for the network input + X1 = (total_energy - total_energy_clip1) / total_energy + X2 = (total_energy_clip1 - total_energy_clip2) / total_energy_clip1 + X3 = (total_energy_clip2 - total_energy_clip3) / total_energy_clip2 + X4 = nnodes(dg) + network_input = SVector(X1, X2, X3, X4) + + # Scale input data + network_input = network_input / + max(maximum(abs, network_input), one(eltype(network_input))) + probability_troubled_cell = network(network_input)[1] + + # Compute indicator value + alpha[element] = probability_to_indicator(probability_troubled_cell, + alpha_continuous, + alpha_amr, alpha_min, alpha_max) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end + +function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})(u, + mesh::TreeMesh{ + 2 + }, + equations, + dg::DGSEM, + cache; + kwargs...) + @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann + + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded, network_input, neighbor_ids, neighbor_mean = indicator_ann.cache #X, network_input + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + c2e = zeros(Int, length(mesh.tree)) + for element in eachelement(dg, cache) + c2e[cache.elements.cell_ids[element]] = element + end + + X = Array{Float64}(undef, 3, nelements(dg, cache)) + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] + + # Calculate indicator variables at Gauss-Lobatto nodes + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + indicator[i, j] = indicator_ann.variable(u_local, equations) + end + + # Convert to modal representation + multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, + indicator, modal_tmp1) + # Save linear modal coefficients for the network input + X[1, element] = modal[1, 1] + X[2, element] = modal[1, 2] + X[3, element] = modal[2, 1] + end + + @threaded for element in eachelement(dg, cache) + cell_id = cache.elements.cell_ids[element] + + network_input[1] = X[1, element] + network_input[2] = X[2, element] + network_input[3] = X[3, element] + + for direction in eachdirection(mesh.tree) + if direction == 1 # -x + dir = 4 + elseif direction == 2 # +x + dir = 1 + elseif direction == 3 # -y + dir = 3 + elseif direction == 4 # +y + dir = 2 + end + + # Of no neighbor exists and current cell is not small + if !has_any_neighbor(mesh.tree, cell_id, direction) + network_input[3 * dir + 1] = X[1, element] + network_input[3 * dir + 2] = X[2, element] + network_input[3 * dir + 3] = X[3, element] + continue + end + + # Get Input data from neighbors + if has_neighbor(mesh.tree, cell_id, direction) + neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id] + if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor + # Mean over 4 neighbor cells + neighbor_ids[1] = mesh.tree.child_ids[1, neighbor_cell_id] + neighbor_ids[2] = mesh.tree.child_ids[2, neighbor_cell_id] + neighbor_ids[3] = mesh.tree.child_ids[3, neighbor_cell_id] + neighbor_ids[4] = mesh.tree.child_ids[4, neighbor_cell_id] + + for i in 1:4 + if has_children(mesh.tree, neighbor_ids[i]) + neighbor_ids5 = c2e[mesh.tree.child_ids[1, neighbor_ids[i]]] + neighbor_ids6 = c2e[mesh.tree.child_ids[2, neighbor_ids[i]]] + neighbor_ids7 = c2e[mesh.tree.child_ids[3, neighbor_ids[i]]] + neighbor_ids8 = c2e[mesh.tree.child_ids[4, neighbor_ids[i]]] + + neighbor_mean[i, 1] = (X[1, neighbor_ids5] + + X[1, neighbor_ids6] + + X[1, neighbor_ids7] + + X[1, neighbor_ids8]) / 4 + neighbor_mean[i, 2] = (X[2, neighbor_ids5] + + X[2, neighbor_ids6] + + X[2, neighbor_ids7] + + X[2, neighbor_ids8]) / 4 + neighbor_mean[i, 3] = (X[3, neighbor_ids5] + + X[3, neighbor_ids6] + + X[3, neighbor_ids7] + + X[3, neighbor_ids8]) / 4 + else + neighbor_id = c2e[neighbor_ids[i]] + neighbor_mean[i, 1] = X[1, neighbor_id] + neighbor_mean[i, 2] = X[2, neighbor_id] + neighbor_mean[i, 3] = X[3, neighbor_id] + end + end + network_input[3 * dir + 1] = (neighbor_mean[1, 1] + + neighbor_mean[2, 1] + + neighbor_mean[3, 1] + + neighbor_mean[4, 1]) / 4 + network_input[3 * dir + 2] = (neighbor_mean[1, 2] + + neighbor_mean[2, 2] + + neighbor_mean[3, 2] + + neighbor_mean[4, 2]) / 4 + network_input[3 * dir + 3] = (neighbor_mean[1, 3] + + neighbor_mean[2, 3] + + neighbor_mean[3, 3] + + neighbor_mean[4, 3]) / 4 + + else # Cell has same refinement level neighbor + neighbor_id = c2e[neighbor_cell_id] + network_input[3 * dir + 1] = X[1, neighbor_id] + network_input[3 * dir + 2] = X[2, neighbor_id] + network_input[3 * dir + 3] = X[3, neighbor_id] + end + else # Cell is small and has large neighbor + parent_id = mesh.tree.parent_ids[cell_id] + neighbor_id = c2e[mesh.tree.neighbor_ids[direction, parent_id]] + + network_input[3 * dir + 1] = X[1, neighbor_id] + network_input[3 * dir + 2] = X[2, neighbor_id] + network_input[3 * dir + 3] = X[3, neighbor_id] + end + end + + # Scale input data + network_input = network_input / + max(maximum(abs, network_input), one(eltype(network_input))) + probability_troubled_cell = network(network_input)[1] + + # Compute indicator value + alpha[element] = probability_to_indicator(probability_troubled_cell, + alpha_continuous, + alpha_amr, alpha_min, alpha_max) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end + +function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkCNN})(u, mesh::TreeMesh{2}, + equations, dg::DGSEM, + cache; kwargs...) + @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann + + @unpack alpha, alpha_tmp, indicator_threaded, nodes, cnn_nodes, vandermonde, network_input = indicator_ann.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + + # Calculate indicator variables at Gauss-Lobatto nodes + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + indicator[i, j] = indicator_ann.variable(u_local, equations) + end + + # Interpolate nodal data to 4x4 LGL nodes + for j in 1:4, i in 1:4 + acc = zero(eltype(indicator)) + for jj in eachnode(dg), ii in eachnode(dg) + acc += vandermonde[i, ii] * indicator[ii, jj] * vandermonde[j, jj] + end + network_input[i, j, 1, 1] = acc + end + + # Scale input data + network_input = network_input / + max(maximum(abs, network_input), one(eltype(network_input))) + probability_troubled_cell = network(network_input)[1] + + # Compute indicator value + alpha[element] = probability_to_indicator(probability_troubled_cell, + alpha_continuous, + alpha_amr, alpha_min, alpha_max) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end +end # @muladd diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..bdb6373 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,14 @@ +[deps] +BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" + +[compat] +BSON = "0.3.3" +Flux = "0.13.15, 0.14" +OrdinaryDiffEq = "6.49.1" +# Trixi = "0.6" diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..d9f4b99 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,7 @@ +using Test + +@time @testset "TrixiSmartShockFinder.jl tests" begin + include("test_unit.jl") + include("test_tree_1d_euler.jl") + include("test_tree_2d_euler.jl") +end diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl new file mode 100644 index 0000000..a09c8ab --- /dev/null +++ b/test/test_tree_1d_euler.jl @@ -0,0 +1,40 @@ +module TestExamples1DEuler + +using Test +using TrixiSmartShockFinder +using Trixi + +# Load testing functions from Trixi.jl +include(joinpath(pkgdir(TrixiSmartShockFinder.Trixi), "test", "test_trixi.jl")) + +EXAMPLES_DIR = pkgdir(TrixiSmartShockFinder, "examples", "tree_1d_dgsem") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "Compressible Euler 1D" begin +#! format: noindent + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), + l2=[0.21814833203212694, 0.2818328665444332, 0.5528379124720818], + linf=[1.5548653877320868, 1.4474018998129738, 2.071919577393772], + maxiters=30) +end + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), + l2=[0.22054468879127423, 0.2828269190680846, 0.5542369885642424], + linf=[ + 1.5623359741479623, + 1.4290121654488288, + 2.1040405133123072, + ], + maxiters=30) +end +end + +end # module diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl new file mode 100644 index 0000000..e5ee19e --- /dev/null +++ b/test/test_tree_2d_euler.jl @@ -0,0 +1,134 @@ +module TestExamples2DEuler + +using Test +using TrixiSmartShockFinder +using Trixi + +# Load testing functions from Trixi.jl +include(joinpath(pkgdir(TrixiSmartShockFinder.Trixi), "test", "test_trixi.jl")) + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "Compressible Euler 2D" begin +#! format: noindent + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), + l2=[ + 0.4758794741390833, + 0.21045415565179362, + 0.21045325630191866, + 0.7022517958549878, + ], + linf=[ + 1.710832148442441, + 0.9711663578827681, + 0.9703787873632452, + 2.9619758810532653, + ], + initial_refinement_level=4, + maxiters=50) +end + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), + l2=[ + 0.472445774440313, + 0.2090782039442978, + 0.20885558673697927, + 0.700569533591275, + ], + linf=[ + 1.7066492792835155, + 0.9856122336679919, + 0.9784316656930644, + 2.9372978989672873, + ], + initial_refinement_level=4, + maxiters=50) +end + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl with mortars" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), + l2=[ + 0.016486406327766923, + 0.03097329879894433, + 0.03101012918167401, + 0.15157175775429868, + ], + linf=[ + 0.27688647744873407, + 0.5653724536715139, + 0.565695523611447, + 2.513047611639946, + ], + refinement_patches=((type = "box", + coordinates_min = (-0.25, -0.25), + coordinates_max = (0.25, 0.25)), + (type = "box", + coordinates_min = (-0.125, -0.125), + coordinates_max = (0.125, 0.125))), + initial_refinement_level=4, + maxiters=5) +end + +@trixi_testset "elixir_euler_blast_wave_neuralnetwork_cnn.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_neuralnetwork_cnn.jl"), + l2=[ + 0.4795795496408325, + 0.2125148972465021, + 0.21311260934645868, + 0.7033388737692883, + ], + linf=[ + 1.8295385992182336, + 0.9687795218482794, + 0.9616033072376108, + 2.9513245978047133, + ], + initial_refinement_level=4, + maxiters=50, + rtol=1.0e-7) +end + +@trixi_testset "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl"), + l2=[ + 0.0845430093623868, + 0.09271459184623232, + 0.09271459184623232, + 0.4377291875101709, + ], + linf=[ + 1.3608553480069898, + 1.6822884847136004, + 1.6822884847135997, + 4.2201475428867035, + ], + maxiters=30, + coverage_override=(maxiters = 6,)) +end + +@trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl"), + # This stuff is experimental and annoying to test. In the future, we plan + # to move it to another repository. Thus, we save developer time right now + # and do not run these tests anymore. + # l2 = [0.0009823702998067061, 0.004943231496200673, 0.0048604522073091815, 0.00496983530893294], + # linf = [0.00855717053383187, 0.02087422420794427, 0.017121993783086185, 0.02720703869972585], + maxiters=30, + coverage_override=(maxiters = 2,)) +end +end + +end # module diff --git a/test/test_unit.jl b/test/test_unit.jl new file mode 100644 index 0000000..88a149a --- /dev/null +++ b/test/test_unit.jl @@ -0,0 +1,34 @@ +module TestUnit + +using Test +using TrixiSmartShockFinder +using Trixi + +# TODO: Remove once Trixi.jl is released without these exports +NeuralNetworkPerssonPeraire = TrixiSmartShockFinder.NeuralNetworkPerssonPeraire +NeuralNetworkRayHesthaven = TrixiSmartShockFinder.NeuralNetworkRayHesthaven +IndicatorNeuralNetwork = TrixiSmartShockFinder.IndicatorNeuralNetwork + +# Load testing functions from Trixi.jl +include(joinpath(pkgdir(TrixiSmartShockFinder.Trixi), "test", "test_trixi.jl")) + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +# Run various unit (= non-elixir-triggered) tests +@testset "Unit tests" begin +#! format: noindent + +@timed_testset "Printing indicators/controllers" begin + equations = CompressibleEulerEquations2D(1.4) + basis = LobattoLegendreBasis(3) + indicator_neuralnetwork = IndicatorNeuralNetwork(equations, basis, + indicator_type = NeuralNetworkPerssonPeraire(), + variable = density, + network = nothing) + @test_nowarn show(stdout, indicator_neuralnetwork) +end +end + +end #module