From 4f319f8919ab52baed5f02c78862659e179823d4 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 29 Jan 2024 03:45:16 +0000 Subject: [PATCH 01/11] Bump ZedThree/clang-tidy-review from 0.17.0 to 0.17.1 (#511) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bumps [ZedThree/clang-tidy-review](https://github.com/zedthree/clang-tidy-review) from 0.17.0 to 0.17.1.
Release notes

Sourced from ZedThree/clang-tidy-review's releases.

v0.17.1

What's Changed

Full Changelog: https://github.com/ZedThree/clang-tidy-review/compare/v0.17.0...v0.17.1

Commits

[![Dependabot compatibility score](https://dependabot-badges.githubapp.com/badges/compatibility_score?dependency-name=ZedThree/clang-tidy-review&package-manager=github_actions&previous-version=0.17.0&new-version=0.17.1)](https://docs.github.com/en/github/managing-security-vulnerabilities/about-dependabot-security-updates#about-compatibility-scores) Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/clang-tidy-comments.yml | 2 +- .github/workflows/clang-tidy.yml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/clang-tidy-comments.yml b/.github/workflows/clang-tidy-comments.yml index 76b671405..512ce7ee2 100644 --- a/.github/workflows/clang-tidy-comments.yml +++ b/.github/workflows/clang-tidy-comments.yml @@ -14,7 +14,7 @@ jobs: runs-on: ubuntu-latest steps: - - uses: ZedThree/clang-tidy-review/post@v0.17.0 + - uses: ZedThree/clang-tidy-review/post@v0.17.1 # lgtm_comment_body, max_comments, and annotations need to be set on the posting workflow in a split setup with: # adjust options as necessary diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml index c64592588..3e2fd32e8 100644 --- a/.github/workflows/clang-tidy.yml +++ b/.github/workflows/clang-tidy.yml @@ -17,7 +17,7 @@ jobs: submodules: true fetch-depth: 0 - - uses: ZedThree/clang-tidy-review@v0.17.0 + - uses: ZedThree/clang-tidy-review@v0.17.1 id: review with: config_file: src/.clang-tidy @@ -27,7 +27,7 @@ jobs: split_workflow: true # Uploads an artefact containing clang_fixes.json - - uses: ZedThree/clang-tidy-review/upload@v0.17.0 + - uses: ZedThree/clang-tidy-review/upload@v0.17.1 # If there are any comments, fail the check - if: steps.review.outputs.total_comments > 0 From 089ba268926590db08c30b5242941b7e3de99b21 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 1 Feb 2024 03:10:31 +1100 Subject: [PATCH 02/11] Make ComputeRadPressure return F and S (#512) ### Description Make ComputeRadPressure in radiation_system.hpp return F and S instead of modifying them at place. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [ ] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- src/radiation_system.hpp | 59 +++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 25 deletions(-) diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index 14b1af6aa..d0710a6f9 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -63,6 +63,12 @@ template struct RadSystem_Traits { static constexpr amrex::GpuArray::nGroups + 1> radBoundaries = {0., inf}; }; +// A struct to hold the results of the ComputeRadPressure function. +struct RadPressureResult { + quokka::valarray F; // components of radiation pressure tensor + double S; // maximum wavespeed for the radiation system +}; + /// Class for the radiation moment equations /// template class RadSystem : public HyperbolicSystem @@ -189,9 +195,9 @@ template class RadSystem : public HyperbolicSystem &cons); - template - AMREX_GPU_DEVICE static void ComputeRadPressure(TARRAY &F_L, double &S_L, double erad_L, double Fx_L, double Fy_L, double Fz_L, double fx_L, - double fy_L, double fz_L); + template + AMREX_GPU_DEVICE static auto ComputeRadPressure(double erad_L, double Fx_L, double Fy_L, double Fz_L, double fx_L, double fy_L, double fz_L) + -> RadPressureResult; }; // Compute radiation energy fractions for each photon group from a Planck function, given nGroups, radBoundaries, and temperature @@ -694,10 +700,12 @@ AMREX_GPU_DEVICE auto RadSystem::ComputeCellOpticalDepth(const quokka } template -template -AMREX_GPU_DEVICE void RadSystem::ComputeRadPressure(TARRAY &F, double &S, const double erad, const double Fx, const double Fy, const double Fz, - const double fx, const double fy, const double fz) +template +AMREX_GPU_DEVICE auto RadSystem::ComputeRadPressure(const double erad, const double Fx, const double Fy, const double Fz, const double fx, + const double fy, const double fz) -> RadPressureResult { + // Compute the radiation pressure tensor and the maximum signal speed and return them as a struct. + // check that states are physically admissible AMREX_ASSERT(erad > 0.0); // AMREX_ASSERT(f < 1.0); // there is sometimes a small (<1%) flux @@ -786,9 +794,11 @@ AMREX_GPU_DEVICE void RadSystem::ComputeRadPressure(TARRAY &F, double AMREX_ASSERT(Pny != NAN); AMREX_ASSERT(Pnz != NAN); - F = {Fn, Pnx, Pny, Pnz}; + RadPressureResult result{}; + result.F = {Fn, Pnx, Pny, Pnz}; + result.S = std::max(0.1, std::sqrt(Tnormal)); - S = std::max(0.1, std::sqrt(Tnormal)); + return result; } template @@ -879,17 +889,10 @@ void RadSystem::ComputeFluxes(array_t &x1Flux_in, array_t &x1FluxDiff f_R = std::sqrt(fx_R * fx_R + fy_R * fy_R + fz_R * fz_R); } - // compute radiation pressure F_L and F_R using ComputeRadPressure - quokka::valarray F_L{}; - quokka::valarray F_R{}; - double S_L = NAN; - double S_R = NAN; - // ComputeRadPressure will modify F_L and S_L - ComputeRadPressure(F_L, S_L, erad_L, Fx_L, Fy_L, Fz_L, fx_L, fy_L, fz_L); + // ComputeRadPressure returns F_L_and_S_L or F_R_and_S_R + auto [F_L, S_L] = ComputeRadPressure(erad_L, Fx_L, Fy_L, Fz_L, fx_L, fy_L, fz_L); S_L *= -1.; // speed sign is -1 - // ComputeRadPressure will modify F_R and S_R - ComputeRadPressure(F_R, S_R, erad_R, Fx_R, Fy_R, Fz_R, fx_R, fy_R, fz_R); - // speed sign is +1, S_R unchanged. + auto [F_R, S_R] = ComputeRadPressure(erad_R, Fx_R, Fy_R, Fz_R, fx_R, fy_R, fz_R); // correct for reduced speed of light F_L[0] *= c_hat_ / c_light_; @@ -1240,13 +1243,19 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne auto fy = Fy / (c_light_ * erad); auto fz = Fz / (c_light_ * erad); - std::array, 3> P{}; - double S = NAN; - // loop DIR over {FluxDir::X1, FluxDir::X2, FluxDir::X3} - // ComputeRadPressure will modify P[DIR] and S - ComputeRadPressure(P[0], S, erad, Fx, Fy, Fz, fx, fy, fz); - ComputeRadPressure(P[1], S, erad, Fx, Fy, Fz, fx, fy, fz); - ComputeRadPressure(P[2], S, erad, Fx, Fy, Fz, fx, fy, fz); + std::array, 3> P{}; + { + auto [F, S] = ComputeRadPressure(erad, Fx, Fy, Fz, fx, fy, fz); + P[0] = F; + } + { + auto [F, S] = ComputeRadPressure(erad, Fx, Fy, Fz, fx, fy, fz); + P[1] = F; + } + { + auto [F, S] = ComputeRadPressure(erad, Fx, Fy, Fz, fx, fy, fz); + P[2] = F; + } // loop over spatial dimensions for (int n = 0; n < 3; ++n) { From 11d743650c45629029aa86b30d30fc7dff79d023 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Wed, 31 Jan 2024 14:33:59 -0500 Subject: [PATCH 03/11] update macOS CI to use M1 runner (#513) ### Description macOS GitHub actions can now run on M1 (https://github.blog/changelog/2024-01-30-github-actions-introducing-the-new-m1-macos-runner-available-to-open-source/), so we update it to do so. ### Related issues N/A ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [ ] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --- .github/workflows/cmake-macos.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/cmake-macos.yml b/.github/workflows/cmake-macos.yml index f1b3c037b..a9be31917 100644 --- a/.github/workflows/cmake-macos.yml +++ b/.github/workflows/cmake-macos.yml @@ -17,7 +17,7 @@ env: jobs: build: - runs-on: macos-13 + runs-on: macos-14 steps: - uses: actions/checkout@v4 @@ -60,8 +60,7 @@ jobs: working-directory: ${{runner.workspace}}/build shell: bash # Execute the build. You can specify a specific target with "--target " - # Restrict to 1 build process at a time to avoid OOM kills. - run: cmake --build . --config $BUILD_TYPE --parallel 1 + run: cmake --build . --config $BUILD_TYPE --parallel 2 - name: Create test output directory run: cmake -E make_directory $GITHUB_WORKSPACE/tests From c0344d5e07a19dcce207ddf29be2a06898c99598 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 1 Feb 2024 03:53:08 -0500 Subject: [PATCH 04/11] add HWASAN CI run (#510) ### Description This adds a CI run on Oracle Cloud Infrastructure (with Ampere aarch64 processors) with HWASAN enabled. Processor info: ``` Architecture: aarch64 CPU op-mode(s): 32-bit, 64-bit Byte Order: Little Endian CPU(s): 4 On-line CPU(s) list: 0-3 Vendor ID: ARM Model name: Neoverse-N1 Model: 1 Thread(s) per core: 1 Core(s) per socket: 4 Socket(s): 1 Stepping: r3p1 BogoMIPS: 50.00 Flags: fp asimd evtstrm aes pmull sha1 sha2 crc32 atomics fphp asimdhp cpuid asimdrdm lrcpc dcpop asimddp ssbs ``` ### Related issues Closes https://github.com/quokka-astro/quokka/issues/509. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [x] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [ ] I have requested a reviewer for this PR. --------- Co-authored-by: Piyush Sharda <34922596+psharda@users.noreply.github.com> --- .ci/azure-pipelines-aarch64-hwasan.yml | 41 ++++++++++++++++++++++++++ .ci/azure-pipelines-asan.yml | 40 ------------------------- .ci/azure-pipelines-build-llvm.yml | 35 ---------------------- .ci/azure-pipelines-debug.yml | 40 ------------------------- CMakeLists.txt | 1 + src/CMakeLists.txt | 7 +++++ 6 files changed, 49 insertions(+), 115 deletions(-) create mode 100644 .ci/azure-pipelines-aarch64-hwasan.yml delete mode 100644 .ci/azure-pipelines-asan.yml delete mode 100644 .ci/azure-pipelines-build-llvm.yml delete mode 100644 .ci/azure-pipelines-debug.yml diff --git a/.ci/azure-pipelines-aarch64-hwasan.yml b/.ci/azure-pipelines-aarch64-hwasan.yml new file mode 100644 index 000000000..9dfbe7186 --- /dev/null +++ b/.ci/azure-pipelines-aarch64-hwasan.yml @@ -0,0 +1,41 @@ +# Starter pipeline +# Start with a minimal pipeline that you can customize to build and deploy your code. +# Add steps that build, run tests, deploy, and more: +# https://aka.ms/yaml + +trigger: +- development + +pr: + autoCancel: true + branches: + include: + - development + +jobs: +- job: BuildAndTest + timeoutInMinutes: 180 # how long to run the job before automatically cancelling + pool: oracle-cloud + steps: + - task: CMake@1 + displayName: 'Configure CMake' + inputs: + cmakeArgs: '.. -DCMAKE_BUILD_TYPE=Release -DENABLE_HWASAN=ON -DCMAKE_CXX_FLAGS="-fPIC" -DCMAKE_EXE_LINKER_FLAGS="-fuse-ld=lld" -DCMAKE_SHARED_LINKER_FLAGS="-fuse-ld=lld" -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DAMReX_SPACEDIM=1' + + - task: CMake@1 + displayName: 'Build Quokka' + inputs: + cmakeArgs: '--build .' + + - task: CMake@1 + displayName: 'Run CTest' + inputs: + cmakeArgs: '-E chdir . ctest -j 2 -T Test --output-on-failure' + + - task: PublishTestResults@2 + inputs: + testResultsFormat: cTest + testResultsFiles: build/Testing/*/Test.xml + testRunTitle: $(Agent.JobName) + condition: succeededOrFailed() + displayName: Publish test results diff --git a/.ci/azure-pipelines-asan.yml b/.ci/azure-pipelines-asan.yml deleted file mode 100644 index e5b97c68c..000000000 --- a/.ci/azure-pipelines-asan.yml +++ /dev/null @@ -1,40 +0,0 @@ -# Starter pipeline -# Start with a minimal pipeline that you can customize to build and deploy your code. -# Add steps that build, run tests, deploy, and more: -# https://aka.ms/yaml - -trigger: -- development - -pr: - autoCancel: true - branches: - include: - - development - -pool: avatar - -steps: - -- task: CMake@1 - displayName: 'Configure CMake' - inputs: - cmakeArgs: '.. -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DAMReX_SPACEDIM=1 -DAMReX_GPU_BACKEND=CUDA -DQUOKKA_PYTHON=OFF -DENABLE_ASAN=ON -G Ninja' - -- task: CMake@1 - displayName: 'Build Quokka' - inputs: - cmakeArgs: '--build .' - -- task: CMake@1 - displayName: 'Run CTest' - inputs: - cmakeArgs: '-E chdir . ctest -E "MatterEnergyExchange*" -T Test --output-on-failure' - -- task: PublishTestResults@2 - inputs: - testResultsFormat: cTest - testResultsFiles: build/Testing/*/Test.xml - testRunTitle: $(Agent.JobName) - condition: succeededOrFailed() - displayName: Publish test results diff --git a/.ci/azure-pipelines-build-llvm.yml b/.ci/azure-pipelines-build-llvm.yml deleted file mode 100644 index bf45ccf99..000000000 --- a/.ci/azure-pipelines-build-llvm.yml +++ /dev/null @@ -1,35 +0,0 @@ -trigger: -- development - -pr: - autoCancel: true - branches: - include: - - development - -pool: avatar - -steps: - -- task: CMake@1 - displayName: 'Configure CMake' - inputs: - cmakeArgs: '.. -DCMAKE_BUILD_TYPE=Release -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CUDA_COMPILER=clang -DAMReX_SPACEDIM=3 -DAMReX_GPU_BACKEND=CUDA -DCMAKE_CUDA_ARCHITECTURES=70 -DQUOKKA_PYTHON=OFF -G Ninja' - -- task: CMake@1 - displayName: 'Build Quokka' - inputs: - cmakeArgs: '--build .' - -- task: CMake@1 - displayName: 'Run CTest' - inputs: - cmakeArgs: '-E chdir . ctest -E "MatterEnergyExchange*" -T Test --output-on-failure' - -- task: PublishTestResults@2 - inputs: - testResultsFormat: cTest - testResultsFiles: build/Testing/*/Test.xml - testRunTitle: $(Agent.JobName) - condition: succeededOrFailed() - displayName: Publish test results \ No newline at end of file diff --git a/.ci/azure-pipelines-debug.yml b/.ci/azure-pipelines-debug.yml deleted file mode 100644 index 1ef6b68b6..000000000 --- a/.ci/azure-pipelines-debug.yml +++ /dev/null @@ -1,40 +0,0 @@ -# Starter pipeline -# Start with a minimal pipeline that you can customize to build and deploy your code. -# Add steps that build, run tests, deploy, and more: -# https://aka.ms/yaml - -trigger: -- development - -pr: - autoCancel: true - branches: - include: - - development - -pool: avatar - -steps: - -- task: CMake@1 - displayName: 'Configure CMake' - inputs: - cmakeArgs: '.. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DAMReX_SPACEDIM=1 -DAMReX_GPU_BACKEND=CUDA -DQUOKKA_PYTHON=OFF -G Ninja' - -- task: CMake@1 - displayName: 'Build Quokka' - inputs: - cmakeArgs: '--build .' - -- task: CMake@1 - displayName: 'Run CTest' - inputs: - cmakeArgs: '-E chdir . ctest -E "MatterEnergyExchange*" -T Test --output-on-failure' - -- task: PublishTestResults@2 - inputs: - testResultsFormat: cTest - testResultsFiles: build/Testing/*/Test.xml - testRunTitle: $(Agent.JobName) - condition: succeededOrFailed() - displayName: Publish test results \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 9bbb5cebc..d0b58a3cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,6 +30,7 @@ set(AMReX_TINY_PROFILE ON CACHE BOOL "" FORCE) option(QUOKKA_PYTHON "Compile with Python support (on/off)" ON) option(DISABLE_FMAD "Disable fused multiply-add instructions on GPU (on/off)" ON) option(ENABLE_ASAN "Enable AddressSanitizer and UndefinedBehaviorSanitizer" OFF) +option(ENABLE_HWASAN "Enable HWAddressSanitizer" OFF) option(WARNINGS_AS_ERRORS "Treat compiler warnings as errors" OFF) option(QUOKKA_OPENPMD "Enable OpenPMD output (on/off)" OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4d007bdd8..e38854230 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -60,6 +60,13 @@ if(ENABLE_ASAN) add_link_options(-fsanitize=address -fsanitize=undefined) endif(ENABLE_ASAN) +if(ENABLE_HWASAN) + # enable HWAddressSanitizer for debugging + message(STATUS "Compiling with HWAddressSanitizer *enabled*") + add_compile_options(-fsanitize=hwaddress) + add_link_options(-fsanitize=hwaddress) +endif(ENABLE_HWASAN) + if(DISABLE_FMAD) message(STATUS "Fused multiply-add (FMAD) is disabled for device code.") add_compile_options($<$:-ffp-contract=off>) From b4222b7e92d41231fc115e0298c75faf7d6300ca Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 1 Feb 2024 03:53:48 -0500 Subject: [PATCH 05/11] re-enable CodeQL (#514) ### Description This re-enables CodeQL static analysis based on the settings used by AMReX. ### Related issues Closes https://github.com/quokka-astro/quokka/issues/508. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [ ] I have tested this PR on my local computer and all tests pass. - [ ] I have manually triggered the GPU tests with the magic comment `/azp run`. - [ ] I have requested a reviewer for this PR. --------- Co-authored-by: Piyush Sharda <34922596+psharda@users.noreply.github.com> --- .github/workflows/codeql.yml | 110 +++++++++++++++++++++ .github/workflows/codeql/codeql-config.yml | 10 ++ 2 files changed, 120 insertions(+) create mode 100644 .github/workflows/codeql.yml create mode 100644 .github/workflows/codeql/codeql-config.yml diff --git a/.github/workflows/codeql.yml b/.github/workflows/codeql.yml new file mode 100644 index 000000000..ae613abeb --- /dev/null +++ b/.github/workflows/codeql.yml @@ -0,0 +1,110 @@ +name: CodeQL + +on: + push: + branches: [ "development" ] + pull_request: + branches: [ "development" ] + schedule: + - cron: "27 3 * * 0" + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-codeql + cancel-in-progress: true + +jobs: + analyze: + if: ${{ github.repository == 'quokka-astro/quokka' || github.event_name != 'schedule' }} + name: Analyze + runs-on: ubuntu-latest + permissions: + actions: read + contents: read + security-events: write + + strategy: + fail-fast: false + matrix: + language: [ python, cpp ] + + steps: + - name: Checkout + uses: actions/checkout@v4 + with: + submodules: true + fetch-depth: 0 + + - name: Install Packages (C++) + if: ${{ matrix.language == 'cpp' }} + run: | + sudo apt-get update + sudo apt-get install --yes cmake openmpi-bin libopenmpi-dev libhdf5-openmpi-dev + .github/workflows/dependencies/dependencies_ccache.sh + sudo ln -s /usr/local/bin/ccache /usr/local/bin/g++ + + - name: Set Up Cache + if: ${{ matrix.language == 'cpp' }} + uses: actions/cache@v4 + with: + path: ~/.cache/ccache + key: ccache-${{ github.workflow }}-${{ github.job }}-git-${{ github.sha }} + restore-keys: | + ccache-${{ github.workflow }}-${{ github.job }}-git- + + - name: Configure (C++) + if: ${{ matrix.language == 'cpp' }} + run: | + cmake -S . -B build \ + -DQUOKKA_PYTHON=OFF \ + -DCMAKE_VERBOSE_MAKEFILE=ON \ + -DCMAKE_CXX_COMPILER="/usr/local/bin/g++" + + - name: Initialize CodeQL + uses: github/codeql-action/init@v3 + with: + languages: ${{ matrix.language }} + queries: +security-and-quality + config-file: ./.github/workflows/codeql/codeql-config.yml + + - name: Build (py) + uses: github/codeql-action/autobuild@v3 + if: ${{ matrix.language == 'python' }} + + - name: Build (C++) + if: ${{ matrix.language == 'cpp' }} + run: | + export CCACHE_COMPRESS=1 + export CCACHE_COMPRESSLEVEL=10 + export CCACHE_MAXSIZE=30M + ccache -z + + cmake --build build -j 4 + + ccache -s + du -hs ~/.cache/ccache + + # Make sure CodeQL has something to do + touch src/main.cpp + export CCACHE_DISABLE=1 + cd build + make -j 4 + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v3 + with: + category: "/language:${{ matrix.language }}" + + save_pr_number: + if: github.event_name == 'pull_request' + runs-on: ubuntu-latest + steps: + - name: Save PR number + env: + PR_NUMBER: ${{ github.event.number }} + run: | + echo $PR_NUMBER > pr_number.txt + - uses: actions/upload-artifact@v4 + with: + name: pr_number + path: pr_number.txt + retention-days: 1 diff --git a/.github/workflows/codeql/codeql-config.yml b/.github/workflows/codeql/codeql-config.yml new file mode 100644 index 000000000..90eaf7ad5 --- /dev/null +++ b/.github/workflows/codeql/codeql-config.yml @@ -0,0 +1,10 @@ +query-filters: + - exclude: + id: + - cpp/commented-out-code + - cpp/complex-condition + - cpp/equality-on-floats + - cpp/fixme-comment + - cpp/path-injection + - cpp/poorly-documented-function + - cpp/use-of-goto From d03f52ad638516d77b4a4de02ea767dd842a9104 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:00:01 -0500 Subject: [PATCH 06/11] Bump extern/openPMD-api from `a3fe9b7` to `d9fce66` (#520) Bumps [extern/openPMD-api](https://github.com/openPMD/openPMD-api) from `a3fe9b7` to `d9fce66`.
Commits
  • d9fce66 Don't allow Container insertion in READ_LINEAR (#1590)
  • 8faadde Initial support for HDF5 subfiling (#1580)
  • 17e757f Refactor: Extract ADIOS2 BufferedActions struct to own file, rename to ADIOS2...
  • 3a8a220 Fix: Strings with a single char in Python (#1585)
  • d5524ec Fix availableChunks for READ_LINEAR in ADIOS2 (#1586)
  • 7bb2948 Fix docstring for MyPath (#1587)
  • See full diff in compare view

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- extern/openPMD-api | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/openPMD-api b/extern/openPMD-api index a3fe9b75d..d9fce66dd 160000 --- a/extern/openPMD-api +++ b/extern/openPMD-api @@ -1 +1 @@ -Subproject commit a3fe9b75d81aa1df569b9fb73ad67b772638961f +Subproject commit d9fce66ddf256beb527c3d3e98e212725dae7f47 From 2479420183b2c756d89c2b14eccc293fb3c032e0 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:00:14 -0500 Subject: [PATCH 07/11] Bump extern/fmt from `a5bacf3` to `6321a97` (#518) Bumps [extern/fmt](https://github.com/fmtlib/fmt) from `a5bacf3` to `6321a97`.
Commits

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- extern/fmt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/fmt b/extern/fmt index a5bacf3fe..6321a97d6 160000 --- a/extern/fmt +++ b/extern/fmt @@ -1 +1 @@ -Subproject commit a5bacf3fefcbe927b8aafac2ab17d5c6895cabe2 +Subproject commit 6321a97d6bd1096274b4e1b22f6210cc9031cab4 From d96715f1e9318e412c0d2aac9e23958ea2a7112d Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:00:26 -0500 Subject: [PATCH 08/11] Bump extern/yaml-cpp from `94710bb` to `4afd53b` (#516) Bumps [extern/yaml-cpp](https://github.com/jbeder/yaml-cpp) from `94710bb` to `4afd53b`.
Commits
  • 4afd53b Bump the github-actions group with 1 update
  • 96f5c88 assign fallback value
  • c67d701 Fix indentation of empty sequences and add test
  • 9eb1142 Fix GIT_TAG field in cmake integration example
  • c28295b Add CMake integration example
  • See full diff in compare view

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- extern/yaml-cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/yaml-cpp b/extern/yaml-cpp index 94710bb22..4afd53b0d 160000 --- a/extern/yaml-cpp +++ b/extern/yaml-cpp @@ -1 +1 @@ -Subproject commit 94710bb2213cb6793121eaaa4edfb15abcca7af9 +Subproject commit 4afd53b0d3140ba2eb85816bfced209cf37ff096 From 16395247a9eb6330f6ed9e619ac8ef696475e7d7 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 17:01:08 +0000 Subject: [PATCH 09/11] Bump extern/amrex from `b96b731` to `689144d` (#519) Bumps [extern/amrex](https://github.com/AMReX-Codes/amrex) from `b96b731` to `689144d`.
Commits

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- extern/amrex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/amrex b/extern/amrex index b96b731d4..689144d15 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit b96b731d4a6fa496f80e54c2c0476a50b9dd4e86 +Subproject commit 689144d157a0106faf3d0ae89f8d90b0250cf975 From 4e15bbc51b442a9462f3d98ce70f3a384b789014 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 19:27:52 -0500 Subject: [PATCH 10/11] Bump extern/Microphysics from `7de2124` to `fbeddb3` (#517) Bumps [extern/Microphysics](https://github.com/psharda/Microphysics) from `7de2124` to `fbeddb3`.
Commits
  • fbeddb3 update pynucastro nets to have bion/mion constexpr (#1437)
  • 3700590 constexpr some helm EOS table parameters (#1440)
  • 6fceeec update tols for Urca test_react so it runs (#1438)
  • 6dc51ad fix NSE+SDC when we are doing drive_initial_convection (#1431)
  • d99cee0 update github actions to latest versions (#1429)
  • a0deaae add struct support for runtime_parameters.py (#1422)
  • 9d0655b update for 24.01
  • 3185974 optional quantum correction terms for chabrier1998 (#1428)
  • a83fe00 [WIP] fix powi script and bug in chugunov2009 screening (#1432)
  • 32fa759 Replace std::pow with integer powers to amrex::Math::powi (#1426)
  • Additional commits viewable in compare view

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Piyush Sharda <34922596+psharda@users.noreply.github.com> --- extern/Microphysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/Microphysics b/extern/Microphysics index 7de212422..fbeddb367 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit 7de212422501a88c1c1f536b93217721e3dc6701 +Subproject commit fbeddb367c6ef7d71052e23c20deda164243e643 From 3455fc8ef2ce93033df9d7be11651e427a1f49a1 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 3 Feb 2024 07:54:24 +1100 Subject: [PATCH 11/11] Multigroup advecting pulse test (#486) ### Description Define the multigroup, variable opacity advecting pulse test. See Sec. 4.1.6 of the CASTRO III paper (Zhang et al. 2013). The code perfectly reproduces Figure 6 of the paper. We use the original Newton-Raphson method on the (Egas, Erad) basis. ### Related issues N/A ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [ ] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [x] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .ci/azure-pipelines-aarch64-debug.yml | 2 +- src/CMakeLists.txt | 1 + src/RadhydroPulse/CMakeLists.txt | 2 +- src/RadhydroPulseMG/CMakeLists.txt | 9 + .../test_radhydro_pulse_MG.cpp | 523 ++++++++++++++++++ .../test_radhydro_pulse_MG.hpp | 21 + .../test_radhydro_shock_multigroup.cpp | 3 +- src/planck_integral.hpp | 3 + src/radiation_system.hpp | 172 +++--- 9 files changed, 662 insertions(+), 74 deletions(-) create mode 100644 src/RadhydroPulseMG/CMakeLists.txt create mode 100644 src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp create mode 100644 src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp diff --git a/.ci/azure-pipelines-aarch64-debug.yml b/.ci/azure-pipelines-aarch64-debug.yml index 3d53e9d87..6da3152a7 100644 --- a/.ci/azure-pipelines-aarch64-debug.yml +++ b/.ci/azure-pipelines-aarch64-debug.yml @@ -30,7 +30,7 @@ jobs: - task: CMake@1 displayName: 'Run CTest' inputs: - cmakeArgs: '-E chdir . ctest -j 2 -T Test --output-on-failure' + cmakeArgs: '-E chdir . ctest -j 2 -T Test -E "RadhydroPulseMG*" --output-on-failure' - task: PublishTestResults@2 inputs: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e38854230..689728792 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -172,6 +172,7 @@ add_subdirectory(RadhydroShockCGS) add_subdirectory(RadhydroShockMultigroup) add_subdirectory(RadhydroUniformAdvecting) add_subdirectory(RadhydroPulse) +add_subdirectory(RadhydroPulseMG) add_subdirectory(ODEIntegration) add_subdirectory(Cooling) diff --git a/src/RadhydroPulse/CMakeLists.txt b/src/RadhydroPulse/CMakeLists.txt index c77a58d81..0cc17b65f 100644 --- a/src/RadhydroPulse/CMakeLists.txt +++ b/src/RadhydroPulse/CMakeLists.txt @@ -6,4 +6,4 @@ if (AMReX_SPACEDIM EQUAL 1) endif(AMReX_GPU_BACKEND MATCHES "CUDA") add_test(NAME RadhydroPulse COMMAND test_radhydro_pulse RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) -endif() \ No newline at end of file +endif() diff --git a/src/RadhydroPulseMG/CMakeLists.txt b/src/RadhydroPulseMG/CMakeLists.txt new file mode 100644 index 000000000..e62068b4f --- /dev/null +++ b/src/RadhydroPulseMG/CMakeLists.txt @@ -0,0 +1,9 @@ +if (AMReX_SPACEDIM EQUAL 1) + add_executable(test_radhydro_pulse_MG test_radhydro_pulse_MG.cpp ../fextract.cpp ${QuokkaObjSources}) + + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_pulse_MG) + endif(AMReX_GPU_BACKEND MATCHES "CUDA") + + add_test(NAME RadhydroPulseMG COMMAND test_radhydro_pulse_MG RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp new file mode 100644 index 000000000..23ca8a5c0 --- /dev/null +++ b/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp @@ -0,0 +1,523 @@ +/// \file test_radhydro_pulse_MG.cpp +/// \brief Defines a test problem for multigroup radiation in the diffusion regime with advection by gas. +/// + +#include "test_radhydro_pulse_MG.hpp" +#include "AMReX_BC_TYPES.H" +#include "AMReX_Print.H" +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "physics_info.hpp" + +struct PulseProblem { +}; // dummy type to allow compile-type polymorphism via template specialization +struct AdvPulseProblem { +}; + +constexpr int isconst = 0; +// constexpr int n_groups_ = 1; +// constexpr amrex::GpuArray rad_boundaries_{0., inf}; +constexpr int n_groups_ = 2; +constexpr amrex::GpuArray rad_boundaries_{1e15, 1e17, 1e19}; +// constexpr int n_groups_ = 8; +// constexpr amrex::GpuArray rad_boundaries_{1e15, 3.16e15, 1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19}; +// constexpr int n_groups_ = 64; +// constexpr amrex::GpuArray rad_boundaries_{1.00000000e+15, 1.15478198e+15, 1.33352143e+15, 1.53992653e+15, +// 1.77827941e+15, 2.05352503e+15, 2.37137371e+15, 2.73841963e+15, +// 3.16227766e+15, 3.65174127e+15, 4.21696503e+15, 4.86967525e+15, +// 5.62341325e+15, 6.49381632e+15, 7.49894209e+15, 8.65964323e+15, +// 1.00000000e+16, 1.15478198e+16, 1.33352143e+16, 1.53992653e+16, +// 1.77827941e+16, 2.05352503e+16, 2.37137371e+16, 2.73841963e+16, +// 3.16227766e+16, 3.65174127e+16, 4.21696503e+16, 4.86967525e+16, +// 5.62341325e+16, 6.49381632e+16, 7.49894209e+16, 8.65964323e+16, +// 1.00000000e+17, 1.15478198e+17, 1.33352143e+17, 1.53992653e+17, +// 1.77827941e+17, 2.05352503e+17, 2.37137371e+17, 2.73841963e+17, +// 3.16227766e+17, 3.65174127e+17, 4.21696503e+17, 4.86967525e+17, +// 5.62341325e+17, 6.49381632e+17, 7.49894209e+17, 8.65964323e+17, +// 1.00000000e+18, 1.15478198e+18, 1.33352143e+18, 1.53992653e+18, +// 1.77827941e+18, 2.05352503e+18, 2.37137371e+18, 2.73841963e+18, +// 3.16227766e+18, 3.65174127e+18, 4.21696503e+18, 4.86967525e+18, +// 5.62341325e+18, 6.49381632e+18, 7.49894209e+18, 8.65964323e+18, +// 1.00000000e+19}; + +constexpr double kappa0 = 180.; // cm^2 g^-1 +constexpr double kappa_const = 100.; // cm^2 g^-1 +constexpr double T0 = 1.0e7; // K (temperature) +constexpr double T1 = 2.0e7; // K (temperature) +constexpr double rho0 = 1.2; // g cm^-3 (matter density) +constexpr double a_rad = C::a_rad; +constexpr double c = C::c_light; // speed of light (cgs) +constexpr double chat = c; +constexpr double width = 24.0; // cm, width of the pulse +constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-10; +constexpr double mu = 2.33 * C::m_u; +constexpr double h_planck = C::hplanck; +constexpr double k_B = C::k_B; +constexpr double v0_nonadv = 0.; // non-advecting pulse + +// static diffusion: (for single group) tau = 2e3, beta = 3e-5, beta tau = 6e-2 +constexpr double v0_adv = 1.0e6; // advecting pulse +constexpr double max_time = 2.4e-5; // max_time = 1.0 * width / v1; +// constexpr double max_time = 4.8e-5; // max_time = 2.0 * width / v1; + +// dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 +// constexpr double kappa0 = 1000.; // cm^2 g^-1 +// constexpr double v0_adv = 1.0e8; // advecting pulse +// constexpr double max_time = 1.2e-4; // max_time = 2.0 * width / v1; + +constexpr double T_ref = T0; +constexpr double nu_ref = 1.0e18; // Hz +constexpr double coeff_ = h_planck * nu_ref / (k_B * T_ref); + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_groups_; +}; +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_groups_; +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr bool compute_v_over_c_terms = true; + static constexpr double energy_unit = h_planck; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; +}; +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr bool compute_v_over_c_terms = true; + static constexpr double energy_unit = h_planck; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; +}; + +AMREX_GPU_HOST_DEVICE +auto compute_initial_Tgas(const double x) -> double +{ + // compute temperature profile for Gaussian radiation pulse + const double sigma = width; + return T0 + (T1 - T0) * std::exp(-x * x / (2.0 * sigma * sigma)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_exact_rho(const double x) -> double +{ + // compute density profile for Gaussian radiation pulse + auto T = compute_initial_Tgas(x); + return rho0 * T0 / T + (a_rad * mu / 3. / k_B) * (std::pow(T0, 4) / T - std::pow(T, 3)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_kappa(const double nu, const double Tgas) -> double +{ + // cm^-1 + auto T_ = Tgas / T_ref; + auto nu_ = nu / nu_ref; + return kappa0 * std::pow(T_, -0.5) * std::pow(nu_, -3) * (1.0 - std::exp(-coeff_ * nu_ / T_)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_repres_nu() -> quokka::valarray +{ + // return the geometrical mean as the representative frequency for each group + quokka::valarray nu_rep{}; + if constexpr (n_groups_ == 1) { + nu_rep[0] = nu_ref; + } else { + for (int g = 0; g < n_groups_; ++g) { + nu_rep[g] = std::sqrt(rad_boundaries_[g] * rad_boundaries_[g + 1]); + } + } + return nu_rep; +} + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + quokka::valarray kappaPVec{}; + auto nu_rep = compute_repres_nu(); + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = compute_kappa(nu_rep[g], Tgas) / rho; + } + if constexpr (isconst == 1) { + kappaPVec.fillin(kappa_const); + } + return kappaPVec; +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return RadSystem::ComputePlanckOpacity(rho, Tgas); +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) + -> quokka::valarray +{ + quokka::valarray opacity_deriv{}; + const auto nu_rep = compute_repres_nu(); + const auto T = Tgas / T0; + for (int g = 0; g < nGroups_; ++g) { + const auto nu = nu_rep[g] / nu_ref; + opacity_deriv[g] = + 1. / T0 * kappa0 * (-0.5 * std::pow(T, -1.5) * (1. - std::exp(-coeff_ * nu / T)) - coeff_ * std::pow(T, -2.5) * std::exp(-coeff_ * nu / T)); + opacity_deriv[g] /= rho; + } + if constexpr (isconst == 1) { + opacity_deriv.fillin(0.0); + } + return opacity_deriv; +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) + -> quokka::valarray +{ + return RadSystem::ComputePlanckOpacityTempDerivative(rho, Tgas); +} + +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double +{ + return (1. / 3.); // Eddington approximation +} +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double +{ + return (1. / 3.); // Eddington approximation +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + amrex::GpuArray prob_hi = grid_elem.prob_hi_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); + + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + const double Trad = compute_initial_Tgas(x - x0); + const double rho = compute_exact_rho(x - x0); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double v0 = v0_nonadv; + + auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + amrex::GpuArray prob_hi = grid_elem.prob_hi_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); + + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + const double Trad = compute_initial_Tgas(x - x0); + const double rho = compute_exact_rho(x - x0); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double v0 = v0_adv; + + auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +auto problem_main() -> int +{ + // This problem is a test of radiation diffusion plus advection by gas. + // This makes this problem a stringent test of the radiation advection + // in the diffusion limit. + + // Problem parameters + const int64_t max_timesteps = 1e8; + const double CFL_number = 0.8; + // const int nx = 32; + + const double max_dt = 1e-3; // t_cr = 2 cm / cs = 7e-8 s + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + // periodic boundary condition in the x-direction will not work + BCs_cc[n].setLo(0, amrex::BCType::foextrap); // extrapolate + BCs_cc[n].setHi(0, amrex::BCType::foextrap); + for (int i = 1; i < AMREX_SPACEDIM; ++i) { + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); + } + } + + // Problem 1: non-advecting pulse + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + sim.radiationCflNumber_ = CFL_number; + sim.cflNumber_ = CFL_number; + sim.maxDt_ = max_dt; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); + const int nx = static_cast(position.size()); + amrex::GpuArray prob_lo = sim.geom[0].ProbLoArray(); + amrex::GpuArray prob_hi = sim.geom[0].ProbHiArray(); + + std::vector xs(nx); + std::vector Trad(nx); + std::vector Tgas(nx); + std::vector Vgas(nx); + std::vector rhogas(nx); + + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + // const auto Erad_t = values.at(RadSystem::radEnergy_index)[i]; + double Erad_t = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + const auto rho_t = values.at(RadSystem::gasDensity_index)[i]; + const auto v_t = values.at(RadSystem::x1GasMomentum_index)[i] / rho_t; + const auto Egas = values.at(RadSystem::gasInternalEnergy_index)[i]; + rhogas.at(i) = rho_t; + Trad.at(i) = Trad_t; + Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); + Vgas.at(i) = 1e-5 * v_t; + } + // END OF PROBLEM 1 + + // Problem 2: advecting pulse + + // Problem initialization + RadhydroSimulation sim2(BCs_cc); + + sim2.radiationReconstructionOrder_ = 3; // PPM + sim2.stopTime_ = max_time; + sim2.radiationCflNumber_ = CFL_number; + sim2.maxDt_ = max_dt; + sim2.maxTimesteps_ = max_timesteps; + sim2.plotfileInterval_ = -1; + + // initialize + sim2.setInitialConditions(); + + // evolve + sim2.evolve(); + + // read output variables + auto [position2, values2] = fextract(sim2.state_new_cc_[0], sim2.Geom(0), 0, 0.0); + prob_lo = sim2.geom[0].ProbLoArray(); + prob_hi = sim2.geom[0].ProbHiArray(); + // compute the pixel size + const double dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx); + const double move = v0_adv * sim2.tNew_[0]; + const int n_p = static_cast(move / dx); + const int half = static_cast(nx / 2.0); + const double drift = move - static_cast(n_p) * dx; + const int shift = n_p - static_cast((n_p + half) / nx) * nx; + + std::vector xs2(nx); + std::vector Trad2(nx); + std::vector Tgas2(nx); + std::vector Vgas2(nx); + std::vector rhogas2(nx); + + for (int i = 0; i < nx; ++i) { + int index_ = 0; + if (shift >= 0) { + if (i < shift) { + index_ = nx - shift + i; + } else { + index_ = i - shift; + } + } else { + if (i <= nx - 1 + shift) { + index_ = i - shift; + } else { + index_ = i - (nx + shift); + } + } + const amrex::Real x = position2[i]; + // const auto Erad_t = values2.at(RadSystem::radEnergy_index)[i]; + double Erad_t = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values2.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + const auto rho_t = values2.at(RadSystem::gasDensity_index)[i]; + const auto v_t = values2.at(RadSystem::x1GasMomentum_index)[i] / rho_t; + const auto Egas = values2.at(RadSystem::gasInternalEnergy_index)[i]; + xs2.at(i) = x - drift; + rhogas2.at(index_) = rho_t; + Trad2.at(index_) = Trad_t; + Tgas2.at(index_) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); + Vgas2.at(index_) = 1e-5 * (v_t - v0_adv); + } + // END OF PROBLEM 2 + + // compute error norm + double err_norm = 0.; + double sol_norm = 0.; + for (size_t i = 0; i < xs2.size(); ++i) { + err_norm += std::abs(Tgas[i] - Trad[i]); + err_norm += std::abs(Trad2[i] - Trad[i]); + err_norm += std::abs(Tgas2[i] - Trad[i]); + sol_norm += std::abs(Trad[i]) * 3.0; + } + const double error_tol = 0.008; + const double rel_error = err_norm / sol_norm; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + +#ifdef HAVE_PYTHON + // plot temperature + matplotlibcpp::clf(); + std::map Trad_args; + std::map Tgas_args; + Trad_args["label"] = "Trad (non-advecting)"; + Trad_args["linestyle"] = "-."; + Tgas_args["label"] = "Tgas (non-advecting)"; + Tgas_args["linestyle"] = "--"; + matplotlibcpp::plot(xs, Trad, Trad_args); + matplotlibcpp::plot(xs, Tgas, Tgas_args); + Trad_args["label"] = "Trad (advecting)"; + Tgas_args["label"] = "Tgas (advecting)"; + matplotlibcpp::plot(xs2, Trad2, Trad_args); + matplotlibcpp::plot(xs2, Tgas2, Tgas_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("temperature (K)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_temperature.pdf"); + + // plot gas density profile + matplotlibcpp::clf(); + std::map rho_args; + rho_args["label"] = "gas density (non-advecting)"; + rho_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, rhogas, rho_args); + rho_args["label"] = "gas density (advecting))"; + matplotlibcpp::plot(xs2, rhogas2, rho_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("density (g cm^-3)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_density.pdf"); + + // plot gas velocity profile + matplotlibcpp::clf(); + std::map vgas_args; + vgas_args["label"] = "gas velocity (non-advecting)"; + vgas_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, Vgas, vgas_args); + vgas_args["label"] = "gas velocity (advecting)"; + matplotlibcpp::plot(xs2, Vgas2, vgas_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("velocity (km s^-1)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_velocity.pdf"); + +#endif + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + return status; +} diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp new file mode 100644 index 000000000..eed311866 --- /dev/null +++ b/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp @@ -0,0 +1,21 @@ +#ifndef TEST_RADHYDRO_PULSE_MG_HPP_ // NOLINT +#define TEST_RADHYDRO_PULSE_MG_HPP_ +/// \file test_radiation_marshak.hpp +/// \brief Defines a test problem for radiation in the static diffusion regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// function definitions + +#endif // TEST_RADHYDRO_PULSE_MG_HPP_ diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index 1640fb193..1a21b2dd1 100644 --- a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -40,6 +40,7 @@ constexpr double v1 = 1.73e7; // cm s^-1 [1.72875e7] constexpr double chat = 10.0 * (v0 + c_s0); // reduced speed of light constexpr double Erad0 = a_rad * (T0 * T0 * T0 * T0); // erg cm^-3 +constexpr double Erad_floor_ = Erad0 * 1e-12; constexpr double Egas0 = rho0 * c_v * T0; // erg cm^-3 constexpr double Erad1 = a_rad * (T1 * T1 * T1 * T1); // erg cm^-3 constexpr double Egas1 = rho1 * c_v * T1; // erg cm^-3 @@ -63,7 +64,7 @@ template <> struct RadSystem_Traits { static constexpr double c_light = c; static constexpr double c_hat = chat; static constexpr double radiation_constant = a_rad; - static constexpr double Erad_floor = 0.; + static constexpr double Erad_floor = Erad_floor_; static constexpr bool compute_v_over_c_terms = true; static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{ diff --git a/src/planck_integral.hpp b/src/planck_integral.hpp index deb22e8ac..56d53a337 100644 --- a/src/planck_integral.hpp +++ b/src/planck_integral.hpp @@ -245,8 +245,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto integrate_planck_from_0_to_x(const // y = x * x * x / 3.0; // 1st order y = (-4 + x) * x + 8 * std::log((2 + x) / 2); // 2nd order // Y_INTERP_MIN is the minimum value returned from interpolate_planck_integral. To ensure y is monotonic with respect to x: + // AMREX_ASSERT(y <= Y_INTERP_MIN); if (y > Y_INTERP_MIN) { y = Y_INTERP_MIN; + } else if (y < 0.) { + y = 0.; } } else if (logx >= LOG_X_MAX) { return 1.0; diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index d0710a6f9..4c6ce87b9 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -110,7 +110,6 @@ template class RadSystem : public HyperbolicSystem::c_hat; static constexpr double radiation_constant_ = RadSystem_Traits::radiation_constant; - static constexpr double Erad_floor_ = RadSystem_Traits::Erad_floor; static constexpr bool compute_v_over_c_terms_ = RadSystem_Traits::compute_v_over_c_terms; static constexpr int nGroups_ = Physics_Traits::nGroups; @@ -122,6 +121,7 @@ template class RadSystem : public HyperbolicSystem::Erad_floor / nGroups_; static constexpr double mean_molecular_mass_ = quokka::EOS_Traits::mean_molecular_mass; static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; @@ -174,6 +174,7 @@ template class RadSystem : public HyperbolicSystem double; AMREX_GPU_HOST_DEVICE static auto ComputePlanckOpacity(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeFluxMeanOpacity(double rho, double Tgas) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto ComputeEnergyMeanOpacity(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputePlanckOpacityTempDerivative(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeEintFromEgas(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Etot) -> double; AMREX_GPU_HOST_DEVICE static auto ComputeEgasFromEint(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Eint) -> double; @@ -186,6 +187,13 @@ template class RadSystem : public HyperbolicSystem const &boundaries, amrex::Real temperature) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray; + + AMREX_GPU_HOST_DEVICE static auto ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) + -> quokka::valarray; + template AMREX_GPU_DEVICE static auto ComputeCellOpticalDepth(const quokka::Array4View &consVar, amrex::GpuArray dx, int i, int j, int k) @@ -218,18 +226,41 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(am radEnergyFractions[g] = y - previous; previous = y; } - // A hack to avoid zeros in radEnergyFractions - // TODO(CCH): use energy_floor - for (int g = 0; g < nGroups_; ++g) { - if (radEnergyFractions[g] < Erad_fraction_floor) { - radEnergyFractions[g] = Erad_fraction_floor; - } - } - radEnergyFractions /= sum(radEnergyFractions); + auto tote = sum(radEnergyFractions); + // AMREX_ALWAYS_ASSERT(tote <= 1.0); + // AMREX_ALWAYS_ASSERT(tote > 0.9999); + radEnergyFractions /= tote; return radEnergyFractions; } } +// define ComputeThermalRadiation, returns the thermal radiation power for each photon group. = a_r * T^4 * radEnergyFractions +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + double power = radiation_constant_ * std::pow(temperature, 4); + auto Erad_g = power * radEnergyFractions; + // set floor + for (int g = 0; g < nGroups_; ++g) { + if (Erad_g[g] < Erad_floor_) { + Erad_g[g] = Erad_floor_; + } + } + return Erad_g; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + // by default, d emission/dT = 4 emission / T + auto erad = ComputeThermalRadiation(temperature, boundaries); + return 4. * erad / temperature; +} + // Linear equation solver for matrix with non-zeros at the first row, first column, and diagonal only. // solve the linear system // [a00 a0i] [x0] = [y0] @@ -919,6 +950,8 @@ void RadSystem::ComputeFluxes(array_t &x1Flux_in, array_t &x1FluxDiff // const quokka::valarray epsilon = {S_corr, S_corr, S_corr, S_corr}; // Jiang et al. (2013) quokka::valarray epsilon = {S_corr * S_corr, S_corr, S_corr, S_corr}; // this code + // fix odd-even instability that appears in the asymptotic diffusion limit + // [for details, see section 3.1: https://ui.adsabs.harvard.edu/abs/2022MNRAS.512.1499R/abstract] if ((i + j + k) % 2 == 1) { // revert to more diffusive flux (has no effect in optically-thin limit) epsilon = {S_corr, 1.0, 1.0, 1.0}; @@ -974,6 +1007,12 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const do return kappaFVec; } +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEnergyMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} + template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double /* rho */, const double /* Tgas */) -> quokka::valarray @@ -1043,7 +1082,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (int g = 0; g < nGroups_; ++g) { Erad0Vec[g] = consPrev(i, j, k, radEnergy_index + numRadVars_ * g); } - AMREX_ASSERT(min(Erad0Vec) >= 0.0); + AMREX_ASSERT(min(Erad0Vec) > 0.0); const double Erad0 = sum(Erad0Vec); // load radiation energy source term @@ -1053,12 +1092,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne Src[g] = dt * (chat * radEnergySource(i, j, k, g)); } - const double a_rad = radiation_constant_; double Egas0 = NAN; double Ekin0 = NAN; double Egas_guess = NAN; double T_gas = NAN; - double fourPiB = NAN; + quokka::valarray fourPiB{}; quokka::valarray EradVec_guess{}; quokka::valarray kappaPVec{}; quokka::valarray kappaEVec{}; @@ -1093,14 +1131,14 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne double deltaEgas = NAN; double Rtot = NAN; double dRtot_dEgas = NAN; - quokka::valarray Rvec; - quokka::valarray dRtot_dErad; - quokka::valarray dRvec_dEgas; - quokka::valarray dFG_dErad; - quokka::valarray dFR_dEgas; - quokka::valarray dFR_i_dErad_i; - quokka::valarray deltaErad; - quokka::valarray F_R; + quokka::valarray Rvec{}; + quokka::valarray dRtot_dErad{}; + quokka::valarray dRvec_dEgas{}; + quokka::valarray dFG_dErad{}; + quokka::valarray dFR_dEgas{}; + quokka::valarray dFR_i_dErad_i{}; + quokka::valarray deltaErad{}; + quokka::valarray F_R{}; const double resid_tol = 1.0e-13; // 1.0e-15; const int maxIter = 400; @@ -1111,30 +1149,29 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(T_gas >= 0.); // compute opacity, emissivity - fourPiB = chat * a_rad * std::pow(T_gas, 4); - auto B = fourPiB / (4.0 * M_PI); + fourPiB = chat * ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); - kappaEVec = kappaFVec; // TODO(CCH): define ComputeEnergyMeanOpacity + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); // compute derivatives w/r/t T_gas - const double dB_dTgas = (4.0 * B) / T_gas; + const auto dfourPiB_dTgas = chat * ComputeThermalRadiationTempDerivative(T_gas, radBoundaries_g_copy); // compute residuals - if constexpr (nGroups_ > 1) { - radEnergyFractions = ComputePlanckEnergyFractions(radBoundaries_g_copy, T_gas); - AMREX_ASSERT(min(radEnergyFractions) > 0.); - } else { - radEnergyFractions[0] = 1.0; - } - - Rvec = dt * rho * (kappaPVec * fourPiB * radEnergyFractions - chat * kappaEVec * EradVec_guess); + const auto absorption = chat * kappaEVec * EradVec_guess; + const auto emission = kappaPVec * fourPiB; + const auto Ediff = emission - absorption; + // auto tot_ediff = sum(abs(Ediff)); + // if ((tot_ediff <= Erad_floor_) || (tot_ediff / (sum(absorption) + sum(emission)) < resid_tol)) { + // break; + // } + Rvec = dt * rho * Ediff; Rtot = sum(Rvec); F_G = Egas_guess - Egas0 + (c / chat) * Rtot; F_R = EradVec_guess - Erad0Vec - (Rvec + Src); - // check if converged - if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * max(abs(F_R)) / Etot0 < resid_tol)) { + // check relative convergence of the residuals + if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * sum(abs(F_R)) / Etot0 < resid_tol)) { break; } @@ -1143,20 +1180,17 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(!std::isnan(sum(dkappaP_dTgas))); // prepare to compute Jacobian elements - const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); + const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T dRtot_dErad = -dt * rho * kappaEVec * chat; - // TODO(CCH): dB_dTgas * radEnergyFractions is a small-dT approximation. Define ComputePlanckEnergyFractionsTempDerivative. Note - // this is accurate in 1-group case. - dRvec_dEgas = dt * rho / c_v * - (4 * M_PI * kappaPVec * dB_dTgas * radEnergyFractions + dkappaP_dTgas * fourPiB * radEnergyFractions - - chat * dkappaE_dTgas * EradVec_guess); - dRtot_dEgas = sum(dRvec_dEgas); - // TODO(CCH): Consider Ben's suggestion of using a global Newton method. + dRvec_dEgas = dt * rho / c_v * (kappaPVec * dfourPiB_dTgas + dkappaP_dTgas * fourPiB - chat * dkappaE_dTgas * EradVec_guess); // Equivalent of eta = (eta > 0.0) ? eta : 0.0, to ensure possitivity of Egas_guess - if (dRtot_dEgas < 0.0) { - dRtot_dEgas = 0.0; + for (int g = 0; g < nGroups_; ++g) { + // AMREX_ASSERT(dRvec_dEgas[g] >= 0.0); + if (dRvec_dEgas[g] < 0.0) { + dRvec_dEgas[g] = 0.0; + } } - // AMREX_ASSERT(dRtot_dEgas >= 0.0); + dRtot_dEgas = sum(dRvec_dEgas); // compute Jacobian elements dFG_dEgas = 1.0 + (c / chat) * dRtot_dEgas; @@ -1173,25 +1207,22 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(!std::isnan(deltaEgas)); AMREX_ASSERT(!deltaErad.hasnan()); - // check relative and absolute convergence of E_r - if ((max(abs(deltaErad / Erad0Vec)) < resid_tol) || (max(abs(deltaErad)) < 0.1 * Erad_floor_)) { - break; - } - EradVec_guess += deltaErad; Egas_guess += deltaEgas; AMREX_ASSERT(min(EradVec_guess) >= 0.); AMREX_ASSERT(Egas_guess > 0.); - // set negative values in EradVec_guess to Erad_floor_ and return the excess to Egas_guess - for (int g = 0; g < nGroups_; ++g) { - if (EradVec_guess[g] < 0.0) { - Egas_guess += (EradVec_guess[g] - Erad_floor_) * (c / chat); - EradVec_guess[g] = Erad_floor_; - } + + // check relative and absolute convergence of E_r + // if ((sum(abs(deltaEgas / Egas_guess)) < 1e-7) || (sum(abs(deltaErad)) <= Erad_floor_)) { + // break; + // } + if (std::abs(deltaEgas / Egas_guess) < 1e-7) { + break; } } // END NEWTON-RAPHSON LOOP AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); + // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); } // endif gamma != 1.0 @@ -1209,24 +1240,17 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::GpuArray Frad_t0{}; amrex::GpuArray Frad_t1{}; amrex::GpuArray dMomentum{0., 0., 0.}; - std::array gasMtm{}; - double realFourPiB = NAN; T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note the T_gas is allowed to be NAN - + quokka::valarray realFourPiB{}; if constexpr (gamma_ != 1.0) { - if constexpr (nGroups_ > 1) { - radEnergyFractions = ComputePlanckEnergyFractions(radBoundaries_g_copy, T_gas); - AMREX_ASSERT(min(radEnergyFractions) > 0.); - } else { - radEnergyFractions[0] = 1.0; - } - kappaPVec = ComputePlanckOpacity(rho, T_gas); - realFourPiB = c * a_rad * std::pow(T_gas, 4); - gasMtm = {x1GasMom0, x2GasMom0, x3GasMom0}; + realFourPiB = c * ComputeThermalRadiation(T_gas, radBoundaries_g_copy); } + kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); + kappaPVec = ComputePlanckOpacity(rho, T_gas); + std::array gasMtm = {x1GasMom0, x2GasMom0, x3GasMom0}; + for (int g = 0; g < nGroups_; ++g) { Frad_t0[0] = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); @@ -1259,7 +1283,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // loop over spatial dimensions for (int n = 0; n < 3; ++n) { - double lastTwoTerms = gasMtm[n] * kappaPVec[g] * realFourPiB * radEnergyFractions[g] * chat / c; + double lastTwoTerms = gasMtm[n] * kappaPVec[g] * realFourPiB[g] * chat / c; // loop over the second rank of the radiation pressure tensor for (int z = 0; z < 3; ++z) { lastTwoTerms += chat * kappaFVec[g] * gasMtm[z] * P[n][z + 1]; @@ -1327,6 +1351,12 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } for (int g = 0; g < nGroups_; ++g) { auto radEnergyNew = EradVec_guess[g] + dErad_work * energyLossFractions[g]; + // AMREX_ASSERT(radEnergyNew > 0.0); + if (radEnergyNew < Erad_floor_) { + // return energy to Egas_guess + Egas_guess -= (Erad_floor_ - radEnergyNew) * (c / chat); + radEnergyNew = Erad_floor_; + } consNew(i, j, k, radEnergy_index + numRadVars_ * g) = radEnergyNew; } } else {