From 596a44e648aef97b2567b151b2a8fdb3e1cd5354 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Mon, 20 Nov 2023 12:00:12 +0000 Subject: [PATCH 1/9] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index a7e6187bc..e6f8fa7fb 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,4 @@ artisoptions.h build compile_commands.json *.dSYM/ +.cache/ \ No newline at end of file From bae30ffac1407cdfdc778c518faa12d99ea9bf26 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Mon, 27 Nov 2023 17:00:42 +0000 Subject: [PATCH 2/9] Node-share corrphotoion array and shrink estimators to non-empty cells and included ions (#43) Linting: use trailing return types and avoid redundant void --- .clang-tidy | 12 +- .clangd | 5 + .github/workflows/ci-checks.yml | 110 ++++++-- .github/workflows/ci.yml | 21 +- .gitignore | 5 +- .pre-commit-config.yaml | 7 +- .vscode/extensions.json | 5 + .vscode/launch.json | 18 ++ .vscode/settings.json | 13 + .vscode/tasks.json | 21 ++ Makefile | 151 +++++++--- artisoptions_doc.md | 1 + atomic.cc | 56 ++-- atomic.h | 70 ++--- decay.cc | 2 - decay.h | 30 +- gammapkt.cc | 26 +- globals.cc | 1 - globals.h | 13 +- grid.cc | 260 +++++++++++------- grid.h | 83 +++--- input.cc | 27 +- input.h | 4 +- kpkt.h | 2 +- ltepop.cc | 4 +- ltepop.h | 19 +- macroatom.h | 38 +-- nltepop.h | 2 +- nonthermal.cc | 222 +++++++-------- nonthermal.h | 18 +- packet.h | 4 +- radfield.cc | 33 +-- radfield.h | 18 +- ratecoeff.cc | 8 +- ratecoeff.h | 32 +-- rpkt.cc | 20 +- rpkt.h | 18 +- scripts/clean.sh | 2 +- sn3d.cc | 138 ++++++---- sn3d.h | 90 +++--- stats.cc | 16 +- stats.h | 4 +- tests/.gitignore | 2 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- .../results_md5_final.txt | 4 +- .../results_md5_job0.txt | 4 +- thermalbalance.cc | 11 +- thermalbalance.h | 2 +- update_grid.cc | 42 ++- vectors.h | 4 +- vpkt.h | 2 +- 60 files changed, 1015 insertions(+), 725 deletions(-) create mode 100644 .clangd create mode 100644 .vscode/extensions.json create mode 100644 .vscode/launch.json create mode 100644 .vscode/settings.json create mode 100644 .vscode/tasks.json diff --git a/.clang-tidy b/.clang-tidy index 4771bae80..0036c2bb9 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,4 +1,5 @@ --- +InheritParentConfig: true Checks: > *, -*static-assert*, @@ -15,23 +16,28 @@ Checks: > -cert-err33-c, -cert-err34-c, -cert-err58-cpp, + -clang-analyzer-deadcode.DeadStores, + -clang-diagnostic-unused-macros, -cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-avoid-non-const-global-variables, + -cppcoreguidelines-macro-usage, -cppcoreguidelines-no-malloc, -cppcoreguidelines-owning-memory, -cppcoreguidelines-pro-bounds-array-to-pointer-decay, -cppcoreguidelines-pro-bounds-constant-array-index, -cppcoreguidelines-pro-bounds-pointer-arithmetic, + -cppcoreguidelines-pro-type-cstyle-cast, -cppcoreguidelines-pro-type-reinterpret-cast, -cppcoreguidelines-pro-type-vararg, - -clang-diagnostic-error, -fuchsia-*, + -google-objc*, -google-readability-todo, -google-runtime-int, -hicpp-vararg, -hicpp-no-array-decay, -hicpp-no-malloc, -hicpp-signed-bitwise, + -llvm-header-guard, -misc-use-anonymous-namespace, -misc-no-recursion, -misc-non-private-member-variables-in-classes, @@ -49,8 +55,10 @@ CheckOptions: - key: cppcoreguidelines-init-variables.MathHeader value: - key: cppcoreguidelines-narrowing-conversions.IgnoreConversionFromTypes - value: size_t;ptrdiff_t;size_type;difference_type;time_t;MPI_Aint;unsigned long + value: size_t;ptrdiff_t;size_t;difference_type;time_t;MPI_Aint - key: cppcoreguidelines-narrowing-conversions.WarnOnFloatingPointNarrowingConversion value: 'false' - key: cppcoreguidelines-narrowing-conversions.WarnOnIntegerToFloatingPointNarrowingConversion value: 'false' + - key: cppcoreguidelines-pro-type-member-init.IgnoreArrays + value: 'true' diff --git a/.clangd b/.clangd new file mode 100644 index 000000000..50b0c280d --- /dev/null +++ b/.clangd @@ -0,0 +1,5 @@ +CompileFlags: + # to make the compilation database, pip install compiledb, then run: + # compiledb -n --full-path make + CompilationDatabase: compile_commands.json + Add: [-DMPI_ON=true] diff --git a/.github/workflows/ci-checks.yml b/.github/workflows/ci-checks.yml index 9041d052f..37bb5f256 100644 --- a/.github/workflows/ci-checks.yml +++ b/.github/workflows/ci-checks.yml @@ -5,21 +5,21 @@ on: push: branches-ignore: - classic* - # pull_request: - # branches-ignore: - # - classic* + workflow_dispatch: + # pull_request: + # branches-ignore: + # - classic* jobs: cppcheck: - runs-on: macos-latest - + runs-on: macos-13 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: install dependencies run: | brew update - brew install gsl openmpi cppcheck + brew install gsl openmpi cppcheck || true cp artisoptions_nltenebular.h artisoptions.h - name: run cppcheck and check for errors @@ -34,9 +34,8 @@ jobs: clang-format: runs-on: ubuntu-22.04 - steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Run clang-format style check uses: jidicula/clang-format-action@v4.11.0 @@ -44,29 +43,87 @@ jobs: clang-format-version: '16' check-path: . + clang-tidy: + runs-on: macos-13 + steps: + - uses: actions/checkout@v4 + + - uses: maxim-lobanov/setup-xcode@v1 + with: + xcode-version: 'latest-stable' + + - name: install dependencies + run: | + brew update + brew install openmpi gsl || true + brew link openmpi + + - name: install llvm + run: | + brew install llvm || true + echo "/usr/local/opt/llvm/bin" >> $GITHUB_PATH + + echo "CXX=clang++" >> $GITHUB_ENV + echo "OMPI_CXX=clang++" >> $GITHUB_ENV + + echo "LDFLAGS=-L/usr/local/opt/llvm/lib" >> $GITHUB_ENV + echo "CPPFLAGS=-I/usr/local/opt/llvm/include" >> $GITHUB_ENV + + echo "LDFLAGS=-L/usr/local/opt/libomp/lib" >> $GITHUB_ENV + echo "CXXFLAGS=-I/usr/local/opt/libomp/include" >> $GITHUB_ENV + + echo "LIBRARY_PATH=/usr/local/opt/lib" >> $GITHUB_ENV + echo "LD_LIBRARY_PATH=/usr/local/opt/lib" >> $GITHUB_ENV + + echo "CPATH=/usr/local/opt/include" >> $GITHUB_ENV + + echo "compiledb" > requirements.txt + + - name: Set up Python + uses: actions/setup-python@v4 + with: + cache: pip + python-version: '3.12' + + - name: Generate compile_commands.json + run: | + python3 -m pip install compiledb + cp artisoptions_nltenebular.h artisoptions.h + make version.h + python3 -m compiledb -n --full-path make + + - name: cat compile_commands.json + run: cat compile_commands.json + + - name: run-clang-tidy + run: | + clang-tidy --version + run-clang-tidy + compile: runs-on: ubuntu-22.04 + env: + CXX: g++ strategy: matrix: - compiler: [{name: gcc, ver: 11}, {name: gcc, ver: 12}, {name: gcc, ver: 13}, {name: clang, ver: 14}, {name: clang, ver: 15}] + compiler: + [ + {name: gcc, ver: 11}, + {name: gcc, ver: 12}, + {name: gcc, ver: 13}, + {name: clang, ver: 14}, + {name: clang, ver: 15}, + ] mpi: [ON, OFF] openmp: [ON, OFF] - artisoptionsfile: [artisoptions_classic.h, artisoptions_nltenebular.h] exclude: - - compiler: {name: gcc, ver: 11} - openmp: ON - - compiler: {name: gcc, ver: 12} - openmp: ON - - compiler: {name: gcc, ver: 13} - openmp: ON - compiler: {name: clang, ver: 15} openmp: ON fail-fast: false name: ${{ matrix.compiler.name }}-${{ matrix.compiler.ver }}${{ matrix.mpi == 'ON' && ' MPI' || ''}}${{ matrix.openmp == 'ON' && ' OpenMP' || ''}} - ${{ matrix.artisoptionsfile }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: CPU type and core count id: cpu-count @@ -108,9 +165,9 @@ jobs: - name: install gsl run: sudo apt install libgsl-dev - # - name: Set compiler environment variables (MPI off) - # if: matrix.mpi == 'OFF' - # run: echo "CXX=${{ matrix.compiler.cmd }}" >> $GITHUB_ENV + # - name: Set compiler environment variables (MPI off) + # if: matrix.mpi == 'OFF' + # run: echo "CXX=${{ matrix.compiler.cmd }}" >> $GITHUB_ENV - name: Set compiler environment variables (MPI on) if: matrix.mpi == 'ON' @@ -123,7 +180,12 @@ jobs: run: | ${{ env.CXX }} --version - - name: Compile + - name: Compile classic mode + run: | + cp -v -p artisoptions_classic.h artisoptions.h + make MPI=${{matrix.mpi}} OPENMP=${{matrix.openmp}} -j${{ steps.cpu-count.outputs.count}} sn3d exspec + + - name: Compile nebular mode run: | - cp -v -p ${{ matrix.artisoptionsfile }} artisoptions.h + cp -v -p artisoptions_classic.h artisoptions.h make MPI=${{matrix.mpi}} OPENMP=${{matrix.openmp}} -j${{ steps.cpu-count.outputs.count}} sn3d exspec diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 66db0d6ba..47cec82bd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -7,6 +7,7 @@ on: - classic* schedule: - cron: 0 13 * * 1 + workflow_dispatch: # pull_request: # branches-ignore: # - classic* @@ -21,7 +22,16 @@ jobs: # os: ['ubuntu-latest', 'self-hosted'] os: [ubuntu-22.04] testmode: [OFF, ON] - testname: [classicmode_1d_3dgrid, classicmode_3d, kilonova_1d_1dgrid, kilonova_1d_3dgrid, kilonova_2d_2dgrid, kilonova_2d_3dgrid, nebularonezone_1d_3dgrid] + testname: + [ + classicmode_1d_3dgrid, + classicmode_3d, + kilonova_1d_1dgrid, + kilonova_1d_3dgrid, + kilonova_2d_2dgrid, + kilonova_2d_3dgrid, + nebularonezone_1d_3dgrid, + ] exclude: - os: self-hosted testmode: ON @@ -32,7 +42,7 @@ jobs: name: ${{ matrix.testname }}${{ matrix.testmode == 'ON' && ' testmode ON' || ''}} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: fetch-depth: 0 @@ -139,6 +149,7 @@ jobs: if: always() working-directory: tests/${{ matrix.testname }}_testrun/ run: | + rm *.tmp mkdir job1 ../../scripts/movefiles.sh job1 @@ -170,9 +181,6 @@ jobs: working-directory: tests/${{ matrix.testname }}_testrun/ run: | mkdir output - cp -r output_0-0.txt exspec.txt ratecoeff*.dat *.out output/ - rsync -av job0 output/ - rsync -av job1 output/ cat exspec.txt - name: Checksum job1 output files @@ -187,7 +195,6 @@ jobs: working-directory: tests/${{ matrix.testname }}_testrun run: | touch requirements.txt - rm output/packets*.* find . -name "*_res_*.out" -exec zstd -v -T0 --rm -f {} \; - name: Upload output files @@ -195,7 +202,7 @@ jobs: if: always() && matrix.os != 'selfhosted' && matrix.testmode == 'OFF' with: name: test-${{ matrix.testname }}-output - path: tests/${{ matrix.testname }}_testrun/output + path: tests/${{ matrix.testname }}_testrun - name: Upload checksum files uses: actions/upload-artifact@v3 diff --git a/.gitignore b/.gitignore index e6f8fa7fb..56d06209e 100644 --- a/.gitignore +++ b/.gitignore @@ -13,7 +13,7 @@ browse.VC.db sn3d.dSYM .tags* .tags_swap -.vscode +.vscode/c_cpp_properties.json build .DS_Store artisoptions.h @@ -21,4 +21,5 @@ artisoptions.h build compile_commands.json *.dSYM/ -.cache/ \ No newline at end of file +.cache +*.log diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1bfa2da8f..99fe3f88b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -53,11 +53,6 @@ repos: hooks: - id: make name: make - entry: make + entry: make OPTIMIZE=OFF language: system - types: [c++, makefile] pass_filenames: false - - repo: https://github.com/jumanjihouse/pre-commit-hook-yamlfmt - rev: 0.2.3 # or other specific tag - hooks: - - id: yamlfmt diff --git a/.vscode/extensions.json b/.vscode/extensions.json new file mode 100644 index 000000000..7c35db31e --- /dev/null +++ b/.vscode/extensions.json @@ -0,0 +1,5 @@ +{ + "recommendations": [ + "llvm-vs-code-extensions.vscode-clangd" + ] +} \ No newline at end of file diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 000000000..95fb01828 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,18 @@ +{ + "configurations": [ + { + "name": "C/C++: clang++ build and debug active file", + "type": "cppdbg", + "request": "launch", + "program": "${fileDirname}/${fileBasenameNoExtension}", + "args": [], + "stopAtEntry": false, + "cwd": "${fileDirname}", + "environment": [], + "externalConsole": false, + "MIMode": "lldb", + "preLaunchTask": "C/C++: clang++ build active file" + } + ], + "version": "2.0.0" +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..29a66a967 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,13 @@ +{ + "[cpp]": { + "editor.tabSize": 2, + "editor.codeActionsOnSave": { + "source.fixAll": true + }, + "editor.formatOnSave": true + }, + "[yaml]": { + "editor.tabSize": 4 + }, + "sonarlint.pathToCompileCommands": "${workspaceFolder}/compile_commands.json" +} \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 000000000..373c3f918 --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,21 @@ +{ + "version": "2.0.0", + "tasks": [ + { + "label": "make", + "type": "process", + "command": "make", + "problemMatcher": [], + "group": { + "kind": "build", + "isDefault": true + } + }, + { + "label": "clean", + "type": "process", + "command": "make clean", + "problemMatcher": [] + } + ] +} \ No newline at end of file diff --git a/Makefile b/Makefile index 0955d5017..ce25967df 100644 --- a/Makefile +++ b/Makefile @@ -3,13 +3,83 @@ # place in architecture folder, e.g. build/arm64 BUILD_DIR = build/$(shell uname -m) -CXXFLAGS += -std=c++20 -fstrict-aliasing -ftree-vectorize -flto=auto +CXXFLAGS += -std=c++20 -fstrict-aliasing -ftree-vectorize -flto=auto -Wno-error=unknown-pragmas + +ifeq ($(MPI),) + # MPI option not specified. set to true by default + MPI := ON +endif +ifeq ($(MPI),ON) + CXX = mpicxx + CXXFLAGS += -DMPI_ON=true + BUILD_DIR := $(BUILD_DIR)_mpi +else ifeq ($(MPI),OFF) +else +$(error bad value for MPI option. Should be ON or OFF) +endif + +ifeq ($(TESTMODE),ON) +else ifeq ($(TESTMODE),OFF) +else ifeq ($(TESTMODE),) +else +$(error bad value for testmode option. Should be ON or OFF) +endif + +COMPILER_VERSION := $(shell $(CXX) --version) + +ifneq '' '$(findstring clang,$(COMPILER_VERSION))' + COMPILER_IS_CLANG := TRUE +else ifneq '' '$(findstring g++,$(COMPILER_VERSION))' + COMPILER_IS_CLANG := FALSE +else + $(warning Unknown compiler) + COMPILER_IS_CLANG := FALSE +endif + +ifeq ($(OPENMP),ON) + BUILD_DIR := $(BUILD_DIR)_openmp + + ifeq ($(COMPILER_IS_CLANG),TRUE) + CXXFLAGS += -Xpreprocessor -fopenmp + LDFLAGS += -lomp + else + CXXFLAGS += -fopenmp + endif + +else ifeq ($(OPENMP),OFF) +else ifeq ($(OPENMP),) +else + $(error bad value for openmp option. Should be ON or OFF) +endif + +ifeq ($(STDPAR),ON) + ifeq ($(OPENMP),ON) + $(error cannot combine OPENMP and STDPAR) + endif + + CXXFLAGS += -DSTDPAR_ON=true + BUILD_DIR := $(BUILD_DIR)_stdpar + + ifeq ($(COMPILER_IS_CLANG),TRUE) + else + # CXXFLAGS += -Xlinker -debug_snapshot + LDFLAGS += -ltbb + endif +else ifeq ($(STDPAR),OFF) +else ifeq ($(STDPAR),) +else + $(error bad value for STDPAR option. Should be ON or OFF) +endif -# CXXFLAGS += -Wunreachable-code ifeq ($(shell uname -s),Darwin) # macOS + ifeq ($(COMPILER_IS_CLANG),FALSE) + # fixes linking on macOS with gcc + LDFLAGS += -Wl,-ld_classic + endif + ifeq ($(shell uname -m),arm64) # On Arm, -mcpu combines -march and -mtune CXXFLAGS += -mcpu=native @@ -31,6 +101,7 @@ ifeq ($(shell uname -s),Darwin) # add -lprofiler for gperftools # LDFLAGS += $(LIB) # LDFLAGS += -lprofiler + CXXFLAGS += $(shell pkg-config --cflags ompi) else ifeq ($(USER),localadmin_ccollins) # CXX = c++ @@ -75,47 +146,37 @@ CXXFLAGS += $(shell pkg-config --cflags gsl) CXXFLAGS += -DHAVE_INLINE -DGSL_C99_INLINE ifeq ($(TESTMODE),ON) - CXXFLAGS += -DTESTMODE=true -O3 -DLIBCXX_ENABLE_DEBUG_MODE + CXXFLAGS += -DTESTMODE=true -DLIBCXX_ENABLE_DEBUG_MODE # makes GitHub actions classic test run forever? - # CXXFLAGS += -D_GLIBCXX_DEBUG - CXXFLAGS += -fsanitize=address,undefined -fno-omit-frame-pointer -fno-common - BUILD_DIR := $(BUILD_DIR)_testmode -else - # skip array range checking for better performance and use optimizations - CXXFLAGS += -DTESTMODE=false -DGSL_RANGE_CHECK_OFF -O3 -endif + # CXXFLAGS += -D_GLIBCXX_DEBUG=1 + CXXFLAGS += -fno-omit-frame-pointer -CXXFLAGS += -Werror -Werror=undef -Winline -Wall -Wpedantic -Wredundant-decls -Wundef -Wno-unused-parameter -Wno-unused-function -Wstrict-aliasing -Wno-inline + ifeq ($(COMPILER_IS_CLANG),TRUE) + CXXFLAGS += -fsanitize=address,undefined,integer + else + CXXFLAGS += -fsanitize=address,undefined + endif -ifeq ($(MPI),) - # MPI option not specified. set to true by default - MPI := ON -endif -ifeq ($(MPI),ON) - CXX = mpicxx - CXXFLAGS += -DMPI_ON=true - BUILD_DIR := $(BUILD_DIR)_mpi -else ifeq ($(MPI),OFF) + BUILD_DIR := $(BUILD_DIR)_testmode else -$(error bad value for MPI option. Should be ON or OFF) + # skip GSL range checking for better performance + CXXFLAGS += -DTESTMODE=false -DGSL_RANGE_CHECK_OFF endif -ifeq ($(TESTMODE),ON) -else ifeq ($(TESTMODE),OFF) -else ifeq ($(TESTMODE),) +ifeq ($(OPTIMIZE),OFF) + BUILD_DIR := $(BUILD_DIR)_optimizeoff + CXXFLAGS += -O0 else -$(error bad value for testmode option. Should be ON or OFF) + ifeq ($(FASTMATH),ON) + BUILD_DIR := $(BUILD_DIR)_fastmath + CXXFLAGS += -Ofast -ffast-math -funsafe-math-optimizations -fno-finite-math-only + else + CXXFLAGS += -O3 + endif endif -ifeq ($(OPENMP),ON) - CXXFLAGS += -Xpreprocessor -fopenmp - LDFLAGS += -lomp - BUILD_DIR := $(BUILD_DIR)_openmp -else ifeq ($(OPENMP),OFF) -else ifeq ($(OPENMP),) -else -$(error bad value for testmode option. Should be ON or OFF) -endif +CXXFLAGS += -Werror -Werror=undef -Winline -Wall -Wpedantic -Wredundant-decls -Wundef -Wno-unused-parameter -Wno-unused-function -Wunused-macros -Wno-inline -Wsign-compare + ### use pg when you want to use gprof profiler #CXXFLAGS = -g -pg -Wall -I$(INCLUDE) @@ -133,21 +194,23 @@ exspec_dep = $(exspec_objects:%.o=%.d) all: sn3d exspec -sn3d: $(sn3d_objects) - $(CXX) $(CXXFLAGS) $(sn3d_objects) $(LDFLAGS) -o sn3d -# $(LINK.cpp) $(filter %.o,$^) -o $@ --include $(sn3d_dep) - -sn3dwhole: version.h - $(CXX) $(CXXFLAGS) -g $(sn3d_files) $(LDFLAGS) -o sn3d - $(BUILD_DIR)/%.o: %.cc artisoptions.h Makefile @mkdir -p $(@D) $(CXX) $(CXXFLAGS) -MD -MP -c $< -o $@ $(BUILD_DIR)/sn3d.o $(BUILD_DIR)/exspec.o: version.h artisoptions.h Makefile -exspec: $(exspec_objects) +check: $(sn3d_files) + run-clang-tidy $(sn3d_files) + +sn3d: $(sn3d_objects) artisoptions.h Makefile + $(CXX) $(CXXFLAGS) $(sn3d_objects) $(LDFLAGS) -o sn3d +-include $(sn3d_dep) + +sn3dwhole: version.h artisoptions.h Makefile + $(CXX) $(CXXFLAGS) -g $(sn3d_files) $(LDFLAGS) -o sn3d + +exspec: $(exspec_objects) artisoptions.h Makefile $(CXX) $(CXXFLAGS) $(exspec_objects) $(LDFLAGS) -o exspec -include $(exspec_dep) @@ -162,4 +225,4 @@ version.h: @echo "constexpr const char* GIT_STATUS = \"$(shell git status --short)\";" >> version.h clean: - rm -rf sn3d exspec build version.h *.o *.d + rm -rf sn3d exspec build *.o *.d diff --git a/artisoptions_doc.md b/artisoptions_doc.md index a3f56f342..f53ffff75 100644 --- a/artisoptions_doc.md +++ b/artisoptions_doc.md @@ -220,4 +220,5 @@ constexpr bool KEEP_ALL_RESTART_FILES; // multiply bound-free cooling coefficient by upper level population instead of the upper ion target level population constexpr bool BFCOOLING_USELEVELPOPNOTIONPOP; + ``` \ No newline at end of file diff --git a/atomic.cc b/atomic.cc index 693ad4746..219fd0608 100644 --- a/atomic.cc +++ b/atomic.cc @@ -1,17 +1,24 @@ #include "atomic.h" +#include #include -#include #include "artisoptions.h" #include "grid.h" #include "ltepop.h" #include "sn3d.h" -double last_phixs_nuovernuedge = - -1; // last photoion cross section point as a factor of nu_edge = last_phixs_nuovernuedge -int maxnions = 0; // highest number of ions for any element -int includedions = 0; // number of ions of any element +// last photoion cross section point as a factor of nu_edge = last_phixs_nuovernuedge +double last_phixs_nuovernuedge = -1; + +// highest number of ions for any element +int maxnions = 0; + +// number of ions of any element +int includedions = 0; + +// number of ions of any element excluding the highest ionisation stage of each element +int includedions_excludehighest = 0; std::array phixs_file_version_exists; auto get_continuumindex_phixstargetindex(const int element, const int ion, const int level, const int phixstargetindex) @@ -155,7 +162,7 @@ auto photoionization_crosssection_fromtable(const float *const photoion_xs, cons void set_nelements(const int nelements_in) { globals::elements.resize(nelements_in); } -auto get_nelements() -> int { return globals::elements.size(); } +auto get_nelements() -> int { return static_cast(globals::elements.size()); } auto get_atomicnumber(const int element) -> int /// Returns the atomic number associated with a given elementindex. @@ -184,9 +191,11 @@ auto get_elementindex(const int Z) -> int void update_includedions_maxnions() { includedions = 0; + includedions_excludehighest = 0; maxnions = 0; for (int element = 0; element < get_nelements(); element++) { includedions += get_nions(element); + includedions_excludehighest += get_nions(element) > 0 ? get_nions(element) - 1 : 0; maxnions = std::max(maxnions, get_nions(element)); } } @@ -197,6 +206,12 @@ auto get_includedions() -> int return includedions; } +auto get_includedions_excludehighest() -> int +// returns the number of ions of all elements combined +{ + return includedions_excludehighest; +} + void update_max_nions(const int nions) // Will ensure that maxnions is always greater than or equal to the number of nions // this is called at startup once per element with the number of ions @@ -275,28 +290,21 @@ auto get_uniqueionindex(const int element, const int ion) -> int { assert_testmodeonly(element < get_nelements()); assert_testmodeonly(ion < get_nions(element)); - int index = 0; - for (int e = 0; e < element; e++) { - index += get_nions(e); - } - index += ion; - assert_testmodeonly(index == globals::elements[element].ions[ion].uniqueionindex); - assert_testmodeonly(index < get_includedions()); - return index; + return globals::elements[element].uniqueionindexstart + ion; } -void get_ionfromuniqueionindex(const int allionsindex, int *element, int *ion) { +auto get_ionfromuniqueionindex(const int allionsindex) -> std::tuple { assert_testmodeonly(allionsindex < get_includedions()); - int allionsindex_thiselementfirstion = 0; - for (int e = 0; e < get_nelements(); e++) { - if ((allionsindex - allionsindex_thiselementfirstion) >= get_nions(e)) { - allionsindex_thiselementfirstion += get_nions(e); // skip this element - } else { - *element = e; - *ion = allionsindex - allionsindex_thiselementfirstion; - assert_testmodeonly(get_uniqueionindex(*element, *ion) == allionsindex); - return; + + for (int element = 0; element < get_nelements(); element++) { + if (get_nions(element) == 0) { + continue; + } + const int ion = allionsindex - globals::elements[element].uniqueionindexstart; + if (ion < get_nions(element)) { + assert_testmodeonly(get_uniqueionindex(element, ion) == allionsindex); + return {element, ion}; } } assert_always(false); // allionsindex too high to be valid diff --git a/atomic.h b/atomic.h index 170f4e620..a1c36a9a6 100644 --- a/atomic.h +++ b/atomic.h @@ -9,46 +9,46 @@ constexpr std::array phixsdata_filenames = {"version0ignore", " extern std::array phixs_file_version_exists; // first value in this array is not used but exists so the // indexes match those of the phixsdata_filenames array -int get_continuumindex_phixstargetindex(int element, int ion, int level, int phixstargetindex); -int get_continuumindex(int element, int ion, int level, int upperionlevel); -int get_phixtargetindex(int element, int ion, int level, int upperionlevel); -double get_tau_sobolev(int modelgridindex, int lineindex, double t_current); +auto get_continuumindex_phixstargetindex(int element, int ion, int level, int phixstargetindex) -> int; +auto get_continuumindex(int element, int ion, int level, int upperionlevel) -> int; +auto get_phixtargetindex(int element, int ion, int level, int upperionlevel) -> int; +auto get_tau_sobolev(int modelgridindex, int lineindex, double t_current) -> double; auto get_nnion_tot(int modelgridindex) -> double; -bool is_nlte(int element, int ion, int level); -bool level_isinsuperlevel(int element, int ion, int level); -double photoionization_crosssection_fromtable(const float *photoion_xs, double nu_edge, double nu); +auto is_nlte(int element, int ion, int level) -> bool; +auto level_isinsuperlevel(int element, int ion, int level) -> bool; +auto photoionization_crosssection_fromtable(const float *photoion_xs, double nu_edge, double nu) -> double; void set_nelements(int nelements_in); -int get_nelements(); -int get_atomicnumber(int element); -int get_elementindex(int Z); -int get_includedions(); +auto get_nelements() -> int; +auto get_atomicnumber(int element) -> int; +auto get_elementindex(int Z) -> int; +auto get_includedions() -> int; void update_includedions_maxnions(); -int get_max_nions(); -int get_nions(int element); -int get_ionstage(int element, int ion); -int get_nlevels(int element, int ion); -int get_nlevels_nlte(int element, int ion); -int get_nlevels_groundterm(int element, int ion); -int get_ionisinglevels(int element, int ion); -int get_uniqueionindex(int element, int ion); -void get_ionfromuniqueionindex(int allionsindex, int *element, int *ion); -int get_uniquelevelindex(int element, int ion, int level); +auto get_max_nions() -> int; +auto get_nions(int element) -> int; +auto get_ionstage(int element, int ion) -> int; +auto get_nlevels(int element, int ion) -> int; +auto get_nlevels_nlte(int element, int ion) -> int; +auto get_nlevels_groundterm(int element, int ion) -> int; +auto get_ionisinglevels(int element, int ion) -> int; +auto get_uniqueionindex(int element, int ion) -> int; +[[nodiscard]] auto get_ionfromuniqueionindex(int allionsindex) -> std::tuple; +auto get_uniquelevelindex(int element, int ion, int level) -> int; void get_levelfromuniquelevelindex(int alllevelsindex, int *element, int *ion, int *level); -double epsilon(int element, int ion, int level); -double stat_weight(int element, int ion, int level); -int get_maxrecombininglevel(int element, int ion); -bool ion_has_superlevel(int element, int ion); -int get_ndowntrans(int element, int ion, int level); -int get_nuptrans(int element, int ion, int level); -int get_nphixstargets(int element, int ion, int level); -int get_phixsupperlevel(int element, int ion, int level, int phixstargetindex); -double get_phixsprobability(int element, int ion, int level, int phixstargetindex); +auto epsilon(int element, int ion, int level) -> double; +auto stat_weight(int element, int ion, int level) -> double; +auto get_maxrecombininglevel(int element, int ion) -> int; +auto ion_has_superlevel(int element, int ion) -> bool; +auto get_ndowntrans(int element, int ion, int level) -> int; +auto get_nuptrans(int element, int ion, int level) -> int; +auto get_nphixstargets(int element, int ion, int level) -> int; +auto get_phixsupperlevel(int element, int ion, int level, int phixstargetindex) -> int; +auto get_phixsprobability(int element, int ion, int level, int phixstargetindex) -> double; void set_ndowntrans(int element, int ion, int level, int ndowntrans); void set_nuptrans(int element, int ion, int level, int nuptrans); -double einstein_spontaneous_emission(int lineindex); -double photoionization_crosssection(int element, int ion, int level, double nu_edge, double nu); -double get_phixs_threshold(int element, int ion, int level, int phixstargetindex); -bool elem_has_nlte_levels(int element); -bool elem_has_nlte_levels_search(int element); +auto einstein_spontaneous_emission(int lineindex) -> double; +auto photoionization_crosssection(int element, int ion, int level, double nu_edge, double nu) -> double; +auto get_phixs_threshold(int element, int ion, int level, int phixstargetindex) -> double; +auto elem_has_nlte_levels(int element) -> bool; +auto elem_has_nlte_levels_search(int element) -> bool; #endif // ATOMIC_H diff --git a/decay.cc b/decay.cc index 2015ebe5c..25adb348d 100644 --- a/decay.cc +++ b/decay.cc @@ -8,8 +8,6 @@ #include #include #include -#include -#include #include #include #include diff --git a/decay.h b/decay.h index 9818656ba..1ba263f6d 100644 --- a/decay.h +++ b/decay.h @@ -25,25 +25,25 @@ constexpr std::array all_decaytypes = { decaytypes::DECAYTYPE_BETAMINUS, decaytypes::DECAYTYPE_NONE}; void init_nuclides(const std::vector &zlist, const std::vector &alist); -int get_nucstring_z(const std::string &strnuc); -int get_nucstring_a(const std::string &strnuc); -int get_num_nuclides(); -const char *get_elname(int z); -int get_nuc_z(int nucindex); -int get_nuc_a(int nucindex); -int get_nucindex(int z, int a); -bool nuc_exists(int z, int a); -double nucdecayenergygamma(int nucindex); -double nucdecayenergygamma(int z, int a); +auto get_nucstring_z(const std::string &strnuc) -> int; +auto get_nucstring_a(const std::string &strnuc) -> int; +auto get_num_nuclides() -> int; +auto get_elname(int z) -> const char *; +auto get_nuc_z(int nucindex) -> int; +auto get_nuc_a(int nucindex) -> int; +auto get_nucindex(int z, int a) -> int; +auto nuc_exists(int z, int a) -> bool; +auto nucdecayenergygamma(int nucindex) -> double; +auto nucdecayenergygamma(int z, int a) -> double; void set_nucdecayenergygamma(int nucindex, double value); void update_abundances(int modelgridindex, int timestep, double t_current); -double get_endecay_per_ejectamass_t0_to_time_withexpansion(int modelgridindex, double tstart); -double get_modelcell_simtime_endecay_per_mass(int mgi); +auto get_endecay_per_ejectamass_t0_to_time_withexpansion(int modelgridindex, double tstart) -> double; +auto get_modelcell_simtime_endecay_per_mass(int mgi) -> double; void setup_decaypath_energy_per_mass(); void free_decaypath_energy_per_mass(); -double get_qdot_modelcell(int modelgridindex, double t, int decaytype); -double get_particle_injection_rate(int modelgridindex, double t, int decaytype); -double get_global_etot_t0_tinf(); +auto get_qdot_modelcell(int modelgridindex, double t, int decaytype) -> double; +auto get_particle_injection_rate(int modelgridindex, double t, int decaytype) -> double; +auto get_global_etot_t0_tinf() -> double; void fprint_nuc_abundances(FILE *estimators_file, int modelgridindex, double t_current, int element); void setup_radioactive_pellet(double e0, int mgi, struct packet *pkt_ptr); void cleanup(); diff --git a/gammapkt.cc b/gammapkt.cc index 545fb711a..3798d6adf 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -4,13 +4,11 @@ #include #include #include -#include #include #include #include "decay.h" #include "grid.h" -#include "nonthermal.h" #include "packet.h" #include "sn3d.h" #include "stats.h" @@ -202,17 +200,17 @@ void init_gamma_linelist() { void normalise_grey(int nts) { const double dt = globals::timesteps[nts].width; globals::timesteps[nts].gamma_dep_pathint = 0.; - for (int mgi = 0; mgi < grid::get_npts_model(); mgi++) { - if (grid::get_numassociatedcells(mgi) > 0) { - const double dV = grid::get_modelcell_assocvolume_tmin(mgi) * pow(globals::timesteps[nts].mid / globals::tmin, 3); + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); - globals::timesteps[nts].gamma_dep_pathint += globals::rpkt_emiss[mgi] / globals::nprocs; + const double dV = grid::get_modelcell_assocvolume_tmin(mgi) * pow(globals::timesteps[nts].mid / globals::tmin, 3); - globals::rpkt_emiss[mgi] = globals::rpkt_emiss[mgi] * ONEOVER4PI / dV / dt / globals::nprocs; + globals::timesteps[nts].gamma_dep_pathint += globals::rpkt_emiss[mgi] / globals::nprocs; - assert_testmodeonly(globals::rpkt_emiss[mgi] >= 0.); - assert_testmodeonly(isfinite(globals::rpkt_emiss[mgi])); - } + globals::rpkt_emiss[mgi] = globals::rpkt_emiss[mgi] * ONEOVER4PI / dV / dt / globals::nprocs; + + assert_testmodeonly(globals::rpkt_emiss[mgi] >= 0.); + assert_testmodeonly(isfinite(globals::rpkt_emiss[mgi])); } } @@ -327,13 +325,7 @@ static auto get_chi_compton_rf(const struct packet *pkt_ptr) -> double { // Use this to decide whether the Thompson limit is acceptable. - double sigma_cmf; - if (xx < THOMSON_LIMIT) { - sigma_cmf = SIGMA_T; - } else { - const double fmax = (1 + (2 * xx)); - sigma_cmf = sigma_compton_partial(xx, fmax); - } + const double sigma_cmf = (xx < THOMSON_LIMIT) ? SIGMA_T : sigma_compton_partial(xx, 1 + (2 * xx)); // Now need to multiply by the electron number density. const double chi_cmf = sigma_cmf * grid::get_nnetot(grid::get_cell_modelgridindex(pkt_ptr->where)); diff --git a/globals.cc b/globals.cc index 0e0fa01fb..369193f5e 100644 --- a/globals.cc +++ b/globals.cc @@ -1,6 +1,5 @@ #include "globals.h" -#include "sn3d.h" #ifdef MPI_ON #include #endif diff --git a/globals.h b/globals.h index d9cf9f6a5..81d2914fd 100644 --- a/globals.h +++ b/globals.h @@ -125,6 +125,7 @@ struct elementlist_entry { // cell /// Be aware that this must not be used outside of the update_grid routine /// and their daughters. Neither it will work with OpenMP threads. + int uniqueionindexstart{-1}; /// Index of the lowest ionisation stage of this element float abundance; /// float initstablemeannucmass; /// Atomic mass number in multiple of MH bool has_nlte_levels; @@ -141,9 +142,9 @@ struct linelist_entry { }; struct gslintegration_paras { - const double nu_edge; - const float T; - const float *const photoion_xs; + double nu_edge; + float T; + float *photoion_xs; }; struct rpkt_continuum_absorptioncoeffs { @@ -159,17 +160,17 @@ struct rpkt_continuum_absorptioncoeffs { }; template -struct _chphixstargets { +struct chphixstargets { double corrphotoioncoeff; }; template <> -struct _chphixstargets { +struct chphixstargets { double corrphotoioncoeff; double separatestimrecomb; }; -using chphixstargets_t = struct _chphixstargets; +using chphixstargets_t = struct chphixstargets; #include "macroatom.h" diff --git a/grid.cc b/grid.cc index fbae21906..d0db59a01 100644 --- a/grid.cc +++ b/grid.cc @@ -7,8 +7,8 @@ #include #include #include +#include #include -#include #include #include #include @@ -20,6 +20,10 @@ #include "decay.h" #include "globals.h" #include "input.h" +#include "rpkt.h" +#ifdef MPI_ON +#include "mpi.h" +#endif #include "nltepop.h" #include "nonthermal.h" #include "packet.h" @@ -62,6 +66,7 @@ double *totmassradionuclide = nullptr; /// total mass of each radionuclide in t #ifdef MPI_ON MPI_Win win_nltepops_allcells = MPI_WIN_NULL; MPI_Win win_initradioabund_allcells = MPI_WIN_NULL; +MPI_Win win_corrphotoionrenorm = MPI_WIN_NULL; #endif float *initradioabund_allcells = nullptr; @@ -307,6 +312,12 @@ void set_nnetot(int modelgridindex, float x) { } static void set_ffegrp(int modelgridindex, float x) { + if (!(x >= 0.)) { + printout("WARNING: Fe-group mass fraction %g is negative in cell %d\n", x, modelgridindex); + assert_always(x > -1e-6); + x = 0.; + } + assert_always(x >= 0); assert_always(x <= 1.001); modelgrid[modelgridindex].ffegrp = x; @@ -460,10 +471,17 @@ auto get_modelinitradioabund(const int modelgridindex, const int nucindex) -> fl return modelgrid[modelgridindex].initradioabund[nucindex]; } -static void set_modelinitradioabund(const int modelgridindex, const int nucindex, const float abund) { +static void set_modelinitradioabund(const int modelgridindex, const int nucindex, float abund) { // set the mass fraction of a nuclide in a model grid cell at t=t_model by nuclide index // initradioabund array is in node shared memory assert_always(nucindex >= 0); + if (!(abund >= 0.)) { + printout("WARNING: nuclear mass fraction for nucindex %d = %g is negative in cell %d\n", abund, nucindex, + modelgridindex); + assert_always(abund > -1e-6); + abund = 0.; + } + assert_always(abund >= 0.); assert_always(abund <= 1.); @@ -561,8 +579,8 @@ static void set_elem_stable_abund_from_total(const int mgi, const int element, c double massfracstable = elemabundance - isofracsum; if (massfracstable < 0.) { - if ((isofracsum / elemabundance - 1.) > 1e-4) //  allow some roundoff error - { + //  allow some roundoff error before we complain + if ((isofracsum - elemabundance - 1.) > 1e-4 && std::abs(isofracsum - elemabundance) > 1e-6) { printout("WARNING: cell %d Z=%d element abundance is less than the sum of its radioisotope abundances \n", mgi, atomic_number); printout(" massfrac(Z) %g massfrac_radioisotopes(Z) %g\n", elemabundance, isofracsum); @@ -662,7 +680,8 @@ static void calculate_kappagrey() { double kappa = 0.; if (globals::opacity_case == 0) { kappa = globals::GREY_OP; - } else if (globals::opacity_case == 1) { + } else if (globals::opacity_case == 1 || globals::opacity_case == 4) { + /// kappagrey used for initial grey approximation in case 4 kappa = ((0.9 * get_ffegrp(mgi)) + 0.1) * globals::GREY_OP / ((0.9 * mfeg / mtot_input) + 0.1); } else if (globals::opacity_case == 2) { const double opcase2_normal = globals::GREY_OP * rho_sum / ((0.9 * fe_sum) + (0.1 * (ngrid - empty_cells))); @@ -670,10 +689,6 @@ static void calculate_kappagrey() { } else if (globals::opacity_case == 3) { globals::opcase3_normal = globals::GREY_OP * rho_sum / opcase3_sum; kappa = get_kappagrey(mgi) * globals::opcase3_normal; - } else if (globals::opacity_case == 4) { - /// kappagrey used for initial grey approximation in this case - kappa = ((0.9 * get_ffegrp(mgi)) + 0.1) * globals::GREY_OP / ((0.9 * mfeg / mtot_input) + 0.1); - // kappa = SIGMA_T; } else if (globals::opacity_case == 5) { // electron-fraction-dependent opacities // values from table 1 of Tanaka et al. (2020). @@ -723,22 +738,52 @@ static void allocate_composition_cooling() { const int npts_nonempty = get_nonempty_npts_model(); // add one for the combined empty cell at the end - auto *initmassfracstable_allcells = static_cast(malloc(npts_nonempty * get_nelements() * sizeof(float))); - auto *elem_meanweight_allcells = static_cast(malloc(npts_nonempty * get_nelements() * sizeof(float))); +#ifdef MPI_ON + size_t my_rank_cells_nonempty = nonempty_npts_model / globals::node_nprocs; + // rank_in_node 0 gets any remainder + if (globals::rank_in_node == 0) { + my_rank_cells_nonempty += nonempty_npts_model - (my_rank_cells_nonempty * globals::node_nprocs); + } +#endif + + float *initmassfracstable_allcells = nullptr; + float *elem_meanweight_allcells = nullptr; + +#ifdef MPI_ON + { + MPI_Aint size = my_rank_cells_nonempty * get_nelements() * sizeof(float); + int disp_unit = sizeof(float); + MPI_Win mpiwin = MPI_WIN_NULL; + assert_always(MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, + &initmassfracstable_allcells, &mpiwin) == MPI_SUCCESS); + assert_always(MPI_Win_shared_query(mpiwin, 0, &size, &disp_unit, &initmassfracstable_allcells) == MPI_SUCCESS); + } + + { + MPI_Aint size = my_rank_cells_nonempty * get_nelements() * sizeof(float); + int disp_unit = sizeof(float); + MPI_Win mpiwin = MPI_WIN_NULL; + + assert_always(MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, + &elem_meanweight_allcells, &mpiwin) == MPI_SUCCESS); + assert_always(MPI_Win_shared_query(mpiwin, 0, &size, &disp_unit, &elem_meanweight_allcells) == MPI_SUCCESS); + MPI_Barrier(globals::mpi_comm_node); + } +#else + initmassfracstable_allcells = static_cast(malloc(npts_nonempty * get_nelements() * sizeof(float))); + elem_meanweight_allcells = static_cast(malloc(npts_nonempty * get_nelements() * sizeof(float))); +#endif double *nltepops_allcells = nullptr; if (globals::total_nlte_levels > 0) { #ifdef MPI_ON - int my_rank_cells = nonempty_npts_model / globals::node_nprocs; - // rank_in_node 0 gets any remainder - if (globals::rank_in_node == 0) { - my_rank_cells += nonempty_npts_model - (my_rank_cells * globals::node_nprocs); - } - MPI_Aint size = my_rank_cells * globals::total_nlte_levels * sizeof(double); + MPI_Aint size = my_rank_cells_nonempty * globals::total_nlte_levels * sizeof(double); int disp_unit = sizeof(double); assert_always(MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, &nltepops_allcells, &win_nltepops_allcells) == MPI_SUCCESS); assert_always(MPI_Win_shared_query(win_nltepops_allcells, 0, &size, &disp_unit, &nltepops_allcells) == MPI_SUCCESS); + + MPI_Barrier(globals::mpi_comm_node); #else nltepops_allcells = static_cast(malloc(npts_nonempty * globals::total_nlte_levels * sizeof(double))); #endif @@ -907,6 +952,53 @@ static void allocate_nonemptymodelcells() { allocate_composition_cooling(); + globals::rpkt_emiss = static_cast(calloc((get_npts_model() + 1), sizeof(double))); + + size_t ionestimsize = get_nonempty_npts_model() * get_includedions() * sizeof(double); + +#ifdef MPI_ON + if constexpr (USE_LUT_PHOTOION) { + auto my_rank_cells = get_nonempty_npts_model() / globals::node_nprocs; + // rank_in_node 0 gets any remainder + if (globals::rank_in_node == 0) { + my_rank_cells += get_nonempty_npts_model() - (my_rank_cells * globals::node_nprocs); + } + + auto size = static_cast(my_rank_cells * get_includedions() * sizeof(double)); + int disp_unit = sizeof(double); + assert_always(MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, + &globals::corrphotoionrenorm, &win_corrphotoionrenorm) == MPI_SUCCESS); + assert_always(MPI_Win_shared_query(win_corrphotoionrenorm, 0, &size, &disp_unit, &globals::corrphotoionrenorm) == + MPI_SUCCESS); + } +#else + if constexpr (USE_LUT_PHOTOION) { + globals::corrphotoionrenorm = static_cast(malloc(ionestimsize)); + } +#endif + + if constexpr (USE_LUT_BFHEATING) { + globals::bfheatingestimator = static_cast(malloc(ionestimsize)); +#ifdef DO_TITER + globals::bfheatingestimator_save = static_cast(malloc(ionestimsize)); +#endif + } + + if constexpr (USE_LUT_PHOTOION) { + globals::gammaestimator = static_cast(malloc(ionestimsize)); + +#ifdef DO_TITER + globals::gammaestimator_save = static_cast(malloc(ionestimsize)); +#endif + } + + globals::ffheatingestimator = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); + globals::colheatingestimator = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); +#ifdef DO_TITER + globals::ffheatingestimator_save = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); + globals::colheatingestimator_save = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); +#endif + #ifdef MPI_ON // barrier to make sure node master has set abundance values to node shared memory MPI_Barrier(MPI_COMM_WORLD); @@ -1016,14 +1108,20 @@ static void abundances_read() { // The abundances begin with hydrogen, helium, etc, going as far up the atomic numbers as required double normfactor = 0.; float abundances_in[150] = {0.}; + double abund_in = 0.; for (int anumber = 1; anumber <= 150; anumber++) { abundances_in[anumber - 1] = 0.; - if (!(ssline >> abundances_in[anumber - 1])) { - assert_always(anumber > 1); // at least one element (hydrogen) should have been specified + if (!(ssline >> abund_in)) { + // at least one element (hydrogen) should have been specified for nonempty cells + assert_always(anumber > 1 || get_numassociatedcells(mgi) == 0); break; } - assert_always(abundances_in[anumber - 1] >= 0.); + if (abund_in < 0. || abund_in < std::numeric_limits::min()) { + assert_always(abund_in > -1e-6); + abund_in = 0.; + } + abundances_in[anumber - 1] = static_cast(abund_in); normfactor += abundances_in[anumber - 1]; } @@ -1085,11 +1183,7 @@ static void parse_model_headerline(const std::string &line, std::vector &zl assert_always((columnindex == 4 && get_model_type() == GRID_CARTESIAN3D) || (columnindex == 3 && get_model_type() == GRID_CYLINDRICAL2D)); continue; - } else if (token == "X_Fegroup") { - colnames.push_back(token); - zlist.push_back(-1); - alist.push_back(-1); - } else if (token.starts_with("X_")) { + } else if (token.starts_with("X_") && token != "X_Fegroup") { colnames.push_back(token); const int z = decay::get_nucstring_z(token.substr(2)); // + 2 skips the 'X_' const int a = decay::get_nucstring_a(token.substr(2)); @@ -1138,7 +1232,6 @@ static void read_model_radioabundances(std::fstream &fmodel, std::istringstream assert_always(ssline >> valuein); // usually a mass fraction, but now can be anything if (nucindexlist[i] >= 0) { - assert_testmodeonly(valuein >= 0.); assert_testmodeonly(valuein <= 1.); set_modelinitradioabund(mgi, nucindexlist[i], valuein); } else if (colnames[i] == "X_Fegroup") { @@ -1460,8 +1553,10 @@ static void read_3d_model() for (int axis = 0; axis < 3; axis++) { const double cellwidth = 2 * xmax_tmodel / ncoordgrid[axis]; const double cellpos_expected = -xmax_tmodel + cellwidth * get_cellcoordpointnum(mgi, axis); - // printout("n %d coord %d expected %g found %g rmax %g get_cellcoordpointnum(n, axis) %d ncoordgrid %d\n", - // n, axis, cellpos_expected, cellpos_in[axis], xmax_tmodel, get_cellcoordpointnum(n, axis), ncoordgrid[axis]); + // printout("mgi %d coord %d expected %g found %g or %g rmax %g get_cellcoordpointnum(mgi, axis) %d ncoordgrid + // %d\n", + // mgi, axis, cellpos_expected, cellpos_in[axis], cellpos_in[2 - axis], xmax_tmodel, + // get_cellcoordpointnum(mgi, axis), ncoordgrid[axis]); if (fabs(cellpos_expected - cellpos_in[axis]) > 0.5 * cellwidth) { posmatch_xyz = false; } @@ -1498,7 +1593,7 @@ static void read_3d_model() abort(); } - assert_always(posmatch_zyx ^ posmatch_xyz); // xor because if both match then probably an infinity occurred + // assert_always(posmatch_zyx ^ posmatch_xyz); // xor because if both match then probably an infinity occurred if (posmatch_xyz) { printout("Cell positions in model.txt are consistent with calculated values when x-y-z column order is used.\n"); } @@ -1506,6 +1601,12 @@ static void read_3d_model() printout("Cell positions in model.txt are consistent with calculated values when z-y-x column order is used.\n"); } + if (!posmatch_xyz && !posmatch_zyx) { + printout( + "WARNING: Cell positions in model.txt are not consistent with calculated values in either x-y-z or z-y-x " + "order.\n"); + } + printout("min_den %g [g/cm3]\n", min_den); printout("Effectively used model grid cells: %d\n", nonemptymgi); } @@ -1601,36 +1702,6 @@ void read_ejecta_model() { printout("tmin %g [s] = %.2f [d]\n", globals::tmin, globals::tmin / 86400.); printout("rmax %g [cm] (at t=tmin)\n", globals::rmax); - globals::rpkt_emiss = static_cast(calloc((get_npts_model() + 1), sizeof(double))); - - if constexpr (USE_LUT_PHOTOION) { - globals::corrphotoionrenorm = - static_cast(malloc((get_npts_model() + 1) * get_nelements() * get_max_nions() * sizeof(double))); - globals::gammaestimator = - static_cast(malloc((get_npts_model() + 1) * get_nelements() * get_max_nions() * sizeof(double))); - -#ifdef DO_TITER - globals::gammaestimator_save = - static_cast(malloc((get_npts_model() + 1) * get_nelements() * get_max_nions() * sizeof(double))); -#endif - } - - if constexpr (USE_LUT_BFHEATING) { - globals::bfheatingestimator = - static_cast(malloc((get_npts_model() + 1) * get_nelements() * get_max_nions() * sizeof(double))); -#ifdef DO_TITER - globals::bfheatingestimator_save = - static_cast(malloc((get_npts_model() + 1) * get_nelements() * get_max_nions() * sizeof(double))); -#endif - } - - globals::ffheatingestimator = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); - globals::colheatingestimator = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); -#ifdef DO_TITER - globals::ffheatingestimator_save = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); - globals::colheatingestimator_save = static_cast(malloc((get_npts_model() + 1) * sizeof(double))); -#endif - calc_modelinit_totmassradionuclides(); read_possible_yefile(); @@ -1673,7 +1744,8 @@ static void read_grid_restart_data(const int timestep) { assert_always(fscanf(gridsave_file, "%d ", ×tep_in) == 1); assert_always(timestep_in == timestep); - for (int mgi = 0; mgi < get_npts_model(); mgi++) { + for (int nonemptymgi = 0; nonemptymgi < get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); int mgi_in = -1; float T_R = 0.; float T_e = 0.; @@ -1682,15 +1754,13 @@ static void read_grid_restart_data(const int timestep) { int thick = 0; double rpkt_emiss = 0.; - if (get_numassociatedcells(mgi) > 0) { - assert_always(fscanf(gridsave_file, "%d %a %a %a %a %d %la %a %a", &mgi_in, &T_R, &T_e, &W, &T_J, &thick, - &rpkt_emiss, &modelgrid[mgi].nne, &modelgrid[mgi].nnetot) == 9); + assert_always(fscanf(gridsave_file, "%d %a %a %a %a %d %la %a %a", &mgi_in, &T_R, &T_e, &W, &T_J, &thick, + &rpkt_emiss, &modelgrid[mgi].nne, &modelgrid[mgi].nnetot) == 9); - if (mgi_in != mgi) { - printout("[fatal] read_grid_restart_data: cell mismatch in reading input gridsave.dat ... abort\n"); - printout("[fatal] read_grid_restart_data: read cellnumber %d, expected cellnumber %d\n", mgi_in, mgi); - assert_always(mgi_in == mgi); - } + if (mgi_in != mgi) { + printout("[fatal] read_grid_restart_data: cell mismatch in reading input gridsave.dat ... abort\n"); + printout("[fatal] read_grid_restart_data: read cellnumber %d, expected cellnumber %d\n", mgi_in, mgi); + assert_always(mgi_in == mgi); } assert_always(T_R >= 0.); @@ -1708,8 +1778,7 @@ static void read_grid_restart_data(const int timestep) { if constexpr (USE_LUT_PHOTOION) { for (int element = 0; element < get_nelements(); element++) { - const int nions = get_nions(element); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < (get_nions(element) - 1); ion++) { const int estimindex = get_ionestimindex(mgi, element, ion); assert_always(fscanf(gridsave_file, " %la %la", &globals::corrphotoionrenorm[estimindex], &globals::gammaestimator[estimindex]) == 2); @@ -1752,22 +1821,19 @@ void write_grid_restart_data(const int timestep) { fprintf(gridsave_file, "%d ", timestep); - for (int mgi = 0; mgi < get_npts_model(); mgi++) { - const bool nonemptycell = (get_numassociatedcells(mgi) > 0); + for (int nonemptymgi = 0; nonemptymgi < get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); - if (nonemptycell) { - assert_always(globals::rpkt_emiss[mgi] >= 0.); - fprintf(gridsave_file, "%d %a %a %a %a %d %la %a %a", mgi, get_TR(mgi), get_Te(mgi), get_W(mgi), get_TJ(mgi), - modelgrid[mgi].thick, globals::rpkt_emiss[mgi], modelgrid[mgi].nne, modelgrid[mgi].nnetot); - } + assert_always(globals::rpkt_emiss[mgi] >= 0.); + fprintf(gridsave_file, "%d %a %a %a %a %d %la %a %a", mgi, get_TR(mgi), get_Te(mgi), get_W(mgi), get_TJ(mgi), + modelgrid[mgi].thick, globals::rpkt_emiss[mgi], modelgrid[mgi].nne, modelgrid[mgi].nnetot); if constexpr (USE_LUT_PHOTOION) { for (int element = 0; element < get_nelements(); element++) { - const int nions = get_nions(element); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < (get_nions(element) - 1); ion++) { const int estimindex = get_ionestimindex(mgi, element, ion); - fprintf(gridsave_file, " %la %la", (nonemptycell ? globals::corrphotoionrenorm[estimindex] : 0.), - (nonemptycell ? globals::gammaestimator[estimindex] : 0.)); + fprintf(gridsave_file, " %la %la", globals::corrphotoionrenorm[estimindex], + globals::gammaestimator[estimindex]); } } } @@ -1796,7 +1862,11 @@ static void assign_initial_temperatures() /// according to the local energy density resulting from the 56Ni decay. /// The dilution factor is W=1 in LTE. + printout("Assigning initial temperatures...\n"); + const double tstart = globals::timesteps[0].mid; + int cells_below_mintemp = 0; + int cells_above_maxtemp = 0; for (int nonempymgi = 0; nonempymgi < get_nonempty_npts_model(); nonempymgi++) { const int mgi = get_mgi_of_nonemptymgi(nonempymgi); @@ -1810,11 +1880,13 @@ static void assign_initial_temperatures() pow(CLIGHT / 4 / STEBO * pow(globals::tmin / tstart, 3) * get_rho_tmin(mgi) * decayedenergy_per_mass, 1. / 4.); if (T_initial < MINTEMP) { - printout("mgi %d: T_initial of %g is below MINTEMP %g K, setting to MINTEMP.\n", mgi, T_initial, MINTEMP); + // printout("mgi %d: T_initial of %g is below MINTEMP %g K, setting to MINTEMP.\n", mgi, T_initial, MINTEMP); T_initial = MINTEMP; + cells_below_mintemp++; } else if (T_initial > MAXTEMP) { - printout("mgi %d: T_initial of %g is above MAXTEMP %g K, setting to MAXTEMP.\n", mgi, T_initial, MAXTEMP); + // printout("mgi %d: T_initial of %g is above MAXTEMP %g K, setting to MAXTEMP.\n", mgi, T_initial, MAXTEMP); T_initial = MAXTEMP; + cells_above_maxtemp++; } else if (!std::isfinite(T_initial)) { printout("mgi %d: T_initial of %g is infinite!\n", mgi, T_initial); } @@ -1827,6 +1899,8 @@ static void assign_initial_temperatures() set_W(mgi, 1.); modelgrid[mgi].thick = 0; } + printout(" cells below MINTEMP %g: %d\n", MINTEMP, cells_below_mintemp); + printout(" cells above MAXTEMP %g: %d\n", MAXTEMP, cells_above_maxtemp); } static void setup_nstart_ndo() { @@ -2124,22 +2198,20 @@ void grid_init(int my_rank) } double totmassradionuclide_actual = 0.; - for (int mgi = 0; mgi < get_npts_model(); mgi++) { - if (get_numassociatedcells(mgi) > 0) { - totmassradionuclide_actual += - get_modelinitradioabund(mgi, nucindex) * get_rho_tmin(mgi) * get_modelcell_assocvolume_tmin(mgi); - } + for (int nonemptymgi = 0; nonemptymgi < get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); + totmassradionuclide_actual += + get_modelinitradioabund(mgi, nucindex) * get_rho_tmin(mgi) * get_modelcell_assocvolume_tmin(mgi); } if (totmassradionuclide_actual > 0.) { const double ratio = totmassradionuclide[nucindex] / totmassradionuclide_actual; // printout("nuclide %d ratio %g\n", nucindex, ratio); - for (int mgi = 0; mgi < get_npts_model(); mgi++) { - if (get_numassociatedcells(mgi) > 0) { - const double prev_abund = get_modelinitradioabund(mgi, nucindex); - const double new_abund = prev_abund * ratio; - set_modelinitradioabund(mgi, nucindex, new_abund); - } + for (int nonemptymgi = 0; nonemptymgi < get_nonempty_npts_model(); nonemptymgi++) { + const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi); + const double prev_abund = get_modelinitradioabund(mgi, nucindex); + const double new_abund = prev_abund * ratio; + set_modelinitradioabund(mgi, nucindex, new_abund); } } } diff --git a/grid.h b/grid.h index 80d0834b2..5cc060a2e 100644 --- a/grid.h +++ b/grid.h @@ -58,7 +58,7 @@ struct modelgrid_t { uint_fast8_t thick = 0; }; -constexpr int get_ngriddimensions(void) { +constexpr auto get_ngriddimensions() -> int { switch (GRID_TYPE) { case GRID_SPHERICAL1D: return 1; @@ -66,8 +66,9 @@ constexpr int get_ngriddimensions(void) { return 2; case GRID_CARTESIAN3D: return 3; + default: + assert_always(false); } - assert_always(false); } extern struct modelgrid_t *modelgrid; @@ -76,29 +77,29 @@ extern int ncoordgrid[3]; extern int ngrid; extern char coordlabel[3]; -int get_elements_uppermost_ion(int modelgridindex, int element); +auto get_elements_uppermost_ion(int modelgridindex, int element) -> int; void set_elements_uppermost_ion(int modelgridindex, int element, int newvalue); -double wid_init(int cellindex, int axis); -double get_modelcell_assocvolume_tmin(int modelgridindex); -double get_gridcell_volume_tmin(int cellindex); -double get_cellcoordmax(int cellindex, int axis); -double get_cellcoordmin(int cellindex, int axis); -int get_cellcoordpointnum(int cellindex, int axis); -int get_coordcellindexincrement(int axis); -float get_rho_tmin(int modelgridindex); -float get_rho(int modelgridindex); -float get_nne(int modelgridindex); -float get_nnetot(int modelgridindex); -float get_ffegrp(int modelgridindex); +auto wid_init(int cellindex, int axis) -> double; +auto get_modelcell_assocvolume_tmin(int modelgridindex) -> double; +auto get_gridcell_volume_tmin(int cellindex) -> double; +auto get_cellcoordmax(int cellindex, int axis) -> double; +auto get_cellcoordmin(int cellindex, int axis) -> double; +auto get_cellcoordpointnum(int cellindex, int axis) -> int; +auto get_coordcellindexincrement(int axis) -> int; +auto get_rho_tmin(int modelgridindex) -> float; +auto get_rho(int modelgridindex) -> float; +auto get_nne(int modelgridindex) -> float; +auto get_nnetot(int modelgridindex) -> float; +auto get_ffegrp(int modelgridindex) -> float; void set_elem_abundance(int modelgridindex, int element, float newabundance); -double get_elem_numberdens(int modelgridindex, int element); -double get_initelectronfrac(int modelgridindex); -double get_initenergyq(int modelgridindex); -float get_kappagrey(int modelgridindex); -float get_Te(int modelgridindex); -float get_TR(int modelgridindex); -float get_TJ(int modelgridindex); -float get_W(int modelgridindex); +auto get_elem_numberdens(int modelgridindex, int element) -> double; +auto get_initelectronfrac(int modelgridindex) -> double; +auto get_initenergyq(int modelgridindex) -> double; +auto get_kappagrey(int modelgridindex) -> float; +auto get_Te(int modelgridindex) -> float; +auto get_TR(int modelgridindex) -> float; +auto get_TJ(int modelgridindex) -> float; +auto get_W(int modelgridindex) -> float; void set_nne(int modelgridindex, float nne); void set_nnetot(int modelgridindex, float x); void set_kappagrey(int modelgridindex, float kappagrey); @@ -108,33 +109,33 @@ void set_TR(int modelgridindex, float TR); void set_TJ(int modelgridindex, float TJ); void set_W(int modelgridindex, float W); void grid_init(int my_rank); -float get_modelinitradioabund(int modelgridindex, int nucindex); -float get_stable_initabund(int mgi, int element); -float get_element_meanweight(int mgi, int element); +auto get_modelinitradioabund(int modelgridindex, int nucindex) -> float; +auto get_stable_initabund(int mgi, int element) -> float; +auto get_element_meanweight(int mgi, int element) -> float; void set_element_meanweight(int mgi, int element, float meanweight); -double get_electronfrac(int modelgridindex); -int get_numassociatedcells(int modelgridindex); -int get_modelcell_nonemptymgi(int mgi); -int get_mgi_of_nonemptymgi(int nonemptymgi); -enum gridtypes get_model_type(); +auto get_electronfrac(int modelgridindex) -> double; +auto get_numassociatedcells(int modelgridindex) -> int; +auto get_modelcell_nonemptymgi(int mgi) -> int; +auto get_mgi_of_nonemptymgi(int nonemptymgi) -> int; +auto get_model_type() -> enum gridtypes; void set_model_type(enum gridtypes model_type_value); -int get_npts_model(); -int get_nonempty_npts_model(); -double get_t_model(); -int get_cell_modelgridindex(int cellindex); +auto get_npts_model() -> int; +auto get_nonempty_npts_model() -> int; +auto get_t_model() -> double; +auto get_cell_modelgridindex(int cellindex) -> int; int get_cellindex_from_pos(std::span pos, double time); void read_ejecta_model(); void write_grid_restart_data(int timestep); -int get_maxndo(); -int get_nstart(int rank); -int get_ndo(int rank); -int get_ndo_nonempty(int rank); -double get_totmassradionuclide(int z, int a); +auto get_maxndo() -> int; +auto get_nstart(int rank) -> int; +auto get_ndo(int rank) -> int; +auto get_ndo_nonempty(int rank) -> int; +auto get_totmassradionuclide(int z, int a) -> double; double boundary_distance(std::span dir, std::span pos, double tstart, int cellindex, int *snext, enum cell_boundary *pkt_last_cross); void change_cell(struct packet *pkt_ptr, int snext); -static inline float get_elem_abundance(int modelgridindex, int element) +static inline auto get_elem_abundance(int modelgridindex, int element) -> float // mass fraction of an element (all isotopes combined) { return modelgrid[modelgridindex].composition[element].abundance; diff --git a/input.cc b/input.cc index fcaacbb9b..c39aee20b 100644 --- a/input.cc +++ b/input.cc @@ -10,19 +10,14 @@ #include #include #include -#include #include #include #include "atomic.h" #include "decay.h" -#include "exspec.h" -#include "gammapkt.h" #include "grid.h" #include "kpkt.h" -#include "nltepop.h" #include "ratecoeff.h" -#include "rpkt.h" #include "sn3d.h" #include "vpkt.h" @@ -497,8 +492,8 @@ static void add_transitions_to_unsorted_linelist(const int element, const int io struct transitions *transitions, int *lineindex, std::vector &temp_linelist) { const int lineindex_initial = *lineindex; - const int tottransitions = transitiontable.size(); - int totupdowntrans = 0; + const auto tottransitions = transitiontable.size(); + size_t totupdowntrans = 0; // pass 0 to get transition counts of each level // pass 1 to allocate and fill transition arrays for (int pass = 0; pass < 2; pass++) { @@ -511,7 +506,7 @@ static void add_transitions_to_unsorted_linelist(const int element, const int io MPI_Barrier(MPI_COMM_WORLD); MPI_Win win = MPI_WIN_NULL; - int my_rank_trans = totupdowntrans / globals::node_nprocs; + size_t my_rank_trans = totupdowntrans / globals::node_nprocs; // rank_in_node 0 gets any remainder if (globals::rank_in_node == 0) { my_rank_trans += totupdowntrans - (my_rank_trans * globals::node_nprocs); @@ -543,7 +538,7 @@ static void add_transitions_to_unsorted_linelist(const int element, const int io } totupdowntrans = 0; - for (int ii = 0; ii < tottransitions; ii++) { + for (size_t ii = 0; ii < tottransitions; ii++) { const int level = transitiontable[ii].upper; const int targetlevel = transitiontable[ii].lower; if (pass == 0) { @@ -765,6 +760,7 @@ static void read_atomicdata_files() { globals::elements[element].nions = nions; globals::elements[element].abundance = abundance; /// abundances are expected to be given by mass globals::elements[element].initstablemeannucmass = mass_amu * MH; + globals::elements[element].uniqueionindexstart = uniqueionindex; /// Initialize the elements ionlist globals::elements[element].ions = static_cast(calloc(nions, sizeof(ionlist_entry))); @@ -976,7 +972,7 @@ static void read_atomicdata_files() { #ifdef MPI_ON MPI_Win win = MPI_WIN_NULL; - int my_rank_lines = globals::nlines / globals::node_nprocs; + size_t my_rank_lines = globals::nlines / globals::node_nprocs; // rank_in_node 0 gets any remainder if (globals::rank_in_node == 0) { my_rank_lines += globals::nlines - (my_rank_lines * globals::node_nprocs); @@ -1136,7 +1132,7 @@ static auto search_groundphixslist(double nu_edge, int *index_in_groundlevelcont index = -1; *index_in_groundlevelcontestimator = -1; } else { - int i; + int i = 1; int element = -1; int ion = -1; for (i = 1; i < globals::nbfcontinua_ground; i++) { @@ -1178,7 +1174,7 @@ static auto search_groundphixslist(double nu_edge, int *index_in_groundlevelcont element = globals::groundcont[index].element; ion = globals::groundcont[index].ion; } - *index_in_groundlevelcontestimator = element * get_max_nions() + ion; + *index_in_groundlevelcontestimator = get_uniqueionindex(element, ion); } return index; @@ -1458,7 +1454,7 @@ static void setup_phixs_list() { static_cast(malloc(globals::nbfcontinua * sizeof(struct fullphixslist))); printout("[info] mem_usage: photoionisation list occupies %.3f MB\n", globals::nbfcontinua * (sizeof(fullphixslist)) / 1024. / 1024.); - int nbftables = 0; + size_t nbftables = 0; int allcontindex = 0; for (int element = 0; element < get_nelements(); element++) { const int nions = get_nions(element); @@ -1514,7 +1510,8 @@ static void setup_phixs_list() { #ifdef MPI_ON float *allphixsblock = nullptr; MPI_Win win_allphixsblock = MPI_WIN_NULL; - MPI_Aint size = (globals::rank_in_node == 0) ? nbftables * globals::NPHIXSPOINTS * sizeof(float) : 0; + auto size = + static_cast((globals::rank_in_node == 0) ? nbftables * globals::NPHIXSPOINTS * sizeof(float) : 0); int disp_unit = sizeof(linelist_entry); MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, &allphixsblock, &win_allphixsblock); @@ -1526,7 +1523,7 @@ static void setup_phixs_list() { #endif assert_always(allphixsblock != nullptr); - int nbftableschanged = 0; + size_t nbftableschanged = 0; for (int i = 0; i < globals::nbfcontinua; i++) { globals::allcont_nu_edge[i] = nonconstallcont[i].nu_edge; diff --git a/input.h b/input.h index 49f63025e..12d31ed34 100644 --- a/input.h +++ b/input.h @@ -10,9 +10,9 @@ void read_parameterfile(int rank); void update_parameterfile(int nts); void time_init(); void write_timestep_file(); -bool get_noncommentline(std::fstream &input, std::string &line); +auto get_noncommentline(std::fstream &input, std::string &line) -> bool; -static inline bool lineiscommentonly(const std::string &line) +static inline auto lineiscommentonly(const std::string &line) -> bool // return true for whitepace-only lines, and lines that are exclusively whitepace up to a '#' character { for (size_t i = 0; i < line.length(); i++) { diff --git a/kpkt.h b/kpkt.h index f7f46deda..28bb64b11 100644 --- a/kpkt.h +++ b/kpkt.h @@ -14,7 +14,7 @@ void calculate_cooling_rates(int modelgridindex, struct heatingcoolingrates *hea void do_kpkt_blackbody(struct packet *pkt_ptr); void do_kpkt(struct packet *pkt_ptr, double t2, int nts); -static inline int get_coolinglistoffset(int element, int ion) { +static inline auto get_coolinglistoffset(int element, int ion) -> int { return globals::elements[element].ions[ion].coolingoffset; } diff --git a/ltepop.cc b/ltepop.cc index c9dae7f07..376207534 100644 --- a/ltepop.cc +++ b/ltepop.cc @@ -4,17 +4,15 @@ #include #include -#include #include "atomic.h" #include "decay.h" #include "grid.h" -#include "macroatom.h" #include "nltepop.h" #include "nonthermal.h" #include "ratecoeff.h" +#include "rpkt.h" #include "sn3d.h" -#include "update_grid.h" struct nne_solution_paras { int modelgridindex; diff --git a/ltepop.h b/ltepop.h index 680e573ff..cddec48fa 100644 --- a/ltepop.h +++ b/ltepop.h @@ -7,16 +7,17 @@ #include "atomic.h" #include "sn3d.h" -double get_groundlevelpop(int modelgridindex, int element, int ion); -double calculate_levelpop(int modelgridindex, int element, int ion, int level); -double calculate_levelpop_lte(int modelgridindex, int element, int ion, int level); -double get_levelpop(int modelgridindex, int element, int ion, int level); -double calculate_sahafact(int element, int ion, int level, int upperionlevel, double T, double E_threshold); -double get_nnion(int modelgridindex, int element, int ion); +auto get_groundlevelpop(int modelgridindex, int element, int ion) -> double; +auto calculate_levelpop(int modelgridindex, int element, int ion, int level) -> double; +auto calculate_levelpop_lte(int modelgridindex, int element, int ion, int level) -> double; +auto get_levelpop(int modelgridindex, int element, int ion, int level) -> double; +[[nodiscard]] auto calculate_sahafact(int element, int ion, int level, int upperionlevel, double T, double E_threshold) + -> double; +[[nodiscard]] auto get_nnion(int modelgridindex, int element, int ion) -> double; void calculate_ion_balance_nne(int modelgridindex); void calculate_cellpartfuncts(int modelgridindex, int element); -[[nodiscard]] auto calculate_ionfractions(const int element, const int modelgridindex, const double nne, - const bool force_lte) -> std::vector; -void set_groundlevelpops(const int modelgridindex, const int element, const float nne, const bool force_lte); +[[nodiscard]] auto calculate_ionfractions(int element, int modelgridindex, double nne, bool use_phi_lte) + -> std::vector; +void set_groundlevelpops(int modelgridindex, int element, float nne, bool force_lte); #endif // LTEPOP_H diff --git a/macroatom.h b/macroatom.h index 1bc6ccf6c..b3df8b56b 100644 --- a/macroatom.h +++ b/macroatom.h @@ -27,7 +27,6 @@ enum ma_action { #include "atomic.h" #include "constants.h" -#include "globals.h" #include "packet.h" void macroatom_open_file(int my_rank); @@ -35,23 +34,24 @@ void macroatom_close_file(); void do_macroatom(struct packet *pkt_ptr, int timestep); -double rad_deexcitation_ratecoeff(int modelgridindex, int element, int ion, int upper, int lower, double epsilon_trans, - float A_ul, double upperstatweight, double t_current); -double rad_excitation_ratecoeff(int modelgridindex, int element, int ion, int lower, int uptransindex, - double epsilon_trans, int lineindex, double t_current); -double rad_recombination_ratecoeff(float T_e, float nne, int element, int upperion, int upperionlevel, - int lowerionlevel, int modelgridindex); -double stim_recombination_ratecoeff(float nne, int element, int upperion, int upper, int lower, int modelgridindex); - -double col_recombination_ratecoeff(int modelgridindex, int element, int upperion, int upper, int lower, - double epsilon_trans); -double col_ionization_ratecoeff(float T_e, float nne, int element, int ion, int lower, int phixstargetindex, - double epsilon_trans); - -double col_deexcitation_ratecoeff(float T_e, float nne, double epsilon_trans, int element, int ion, int upper, - const struct level_transition &downtransition); - -double col_excitation_ratecoeff(float T_e, float nne, int element, int ion, int lower, int uptransindex, - double epsilon_trans, double lowerstatweight); +auto rad_deexcitation_ratecoeff(int modelgridindex, int element, int ion, int upper, int lower, double epsilon_trans, + float A_ul, double upperstatweight, double t_current) -> double; +auto rad_excitation_ratecoeff(int modelgridindex, int element, int ion, int lower, int uptransindex, + double epsilon_trans, int lineindex, double t_current) -> double; +auto rad_recombination_ratecoeff(float T_e, float nne, int element, int upperion, int upperionlevel, int lowerionlevel, + int modelgridindex) -> double; +auto stim_recombination_ratecoeff(float nne, int element, int upperion, int upper, int lower, int modelgridindex) + -> double; + +auto col_recombination_ratecoeff(int modelgridindex, int element, int upperion, int upper, int lower, + double epsilon_trans) -> double; +auto col_ionization_ratecoeff(float T_e, float nne, int element, int ion, int lower, int phixstargetindex, + double epsilon_trans) -> double; + +auto col_deexcitation_ratecoeff(float T_e, float nne, double epsilon_trans, int element, int ion, int upper, + const struct level_transition &downtransition) -> double; + +auto col_excitation_ratecoeff(float T_e, float nne, int element, int ion, int lower, int uptransindex, + double epsilon_trans, double lowerstatweight) -> double; #endif // MACROATOM_H diff --git a/nltepop.h b/nltepop.h index 6ab5cda90..9b1d9aec4 100644 --- a/nltepop.h +++ b/nltepop.h @@ -4,7 +4,7 @@ #include void solve_nlte_pops_element(int element, int modelgridindex, int timestep, int nlte_iter); -double superlevel_boltzmann(int modelgridindex, int element, int ion, int level); +auto superlevel_boltzmann(int modelgridindex, int element, int ion, int level) -> double; void nltepop_write_to_file(int modelgridindex, int timestep); void nltepop_open_file(int my_rank); void nltepop_close_file(); diff --git a/nonthermal.cc b/nonthermal.cc index 48e0021b5..907f5b424 100644 --- a/nonthermal.cc +++ b/nonthermal.cc @@ -13,12 +13,10 @@ #include "atomic.h" #include "decay.h" #include "grid.h" -#include "kpkt.h" #include "ltepop.h" #include "macroatom.h" #include "sn3d.h" #include "stats.h" -#include "update_grid.h" namespace nonthermal { @@ -404,9 +402,7 @@ static void zero_all_effionpot(const int modelgridindex) { nt_solution[modelgridindex].ionenfrac_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a] = 0.; } - int element = 0; - int ion = 0; - get_ionfromuniqueionindex(uniqueionindex, &element, &ion); + const auto [element, ion] = get_ionfromuniqueionindex(uniqueionindex); assert_always(fabs(get_auger_probability(modelgridindex, element, ion, 0) - 1.0) < 1e-3); assert_always(fabs(get_ion_auger_enfrac(modelgridindex, element, ion, 0) - 1.0) < 1e-3); } @@ -1684,28 +1680,8 @@ auto nt_excitation_ratecoeff(const int modelgridindex, const int element, const return ratecoeffperdeposition * deposition_rate_density; } -static void select_nt_ionization(int modelgridindex, int *element, int *lowerion) -// select based on stored frac_deposition for each ion -{ - const double zrand = rng_uniform(); - double frac_deposition_ion_sum = 0.; - // zrand is between zero and frac_ionization - // keep subtracting off deposition fractions of ionizations transitions until we hit the right one - // e.g. if zrand was less than frac_dep_trans1, then use the first transition - // e.g. if zrand was between frac_dep_trans1 and frac_dep_trans2 then use the second transition, etc - for (int allionindex = 0; allionindex < get_includedions(); allionindex++) { - frac_deposition_ion_sum += nt_solution[modelgridindex].fracdep_ionization_ion[allionindex]; - if (frac_deposition_ion_sum >= zrand) { - get_ionfromuniqueionindex(allionindex, element, lowerion); - - return; - } - } - assert_always(false); // should not reach here -} - static auto ion_ntion_energyrate(int modelgridindex, int element, int lowerion) -> double { - // returns the energy rate [erg/s] going toward non-thermal ionisation of lowerion + // returns the energy rate [erg/cm3/s] going toward non-thermal ionisation of lowerion const double nnlowerion = get_nnion(modelgridindex, element, lowerion); double enrate = 0.; for (int upperion = lowerion + 1; upperion <= nt_ionisation_maxupperion(element, lowerion); upperion++) { @@ -1737,19 +1713,35 @@ static auto get_ntion_energyrate(int modelgridindex) -> double { return ratetotal; } -static void select_nt_ionization2(int modelgridindex, int *element, int *lowerion) { +static auto select_nt_ionization(int modelgridindex) -> std::tuple { + const double zrand = rng_uniform(); + + // // select based on stored frac_deposition for each ion + // double frac_deposition_ion_sum = 0.; + // // zrand is between zero and frac_ionization + // // keep subtracting off deposition fractions of ionizations transitions until we hit the right one + // // e.g. if zrand was less than frac_dep_trans1, then use the first transition + // // e.g. if zrand was between frac_dep_trans1 and frac_dep_trans2 then use the second transition, etc + // for (int allionindex = 0; allionindex < get_includedions(); allionindex++) { + // frac_deposition_ion_sum += nt_solution[modelgridindex].fracdep_ionization_ion[allionindex]; + // if (frac_deposition_ion_sum >= zrand) { + // get_ionfromuniqueionindex(allionindex, element, lowerion); + + // return; + // } + // } + // assert_always(false); // should not reach here + const double ratetotal = get_ntion_energyrate(modelgridindex); - const double zrand = rng_uniform(); + // select based on the calcuated energy going to ionisation for each ion double ratesum = 0.; for (int ielement = 0; ielement < get_nelements(); ielement++) { const int nions = get_nions(ielement); for (int ilowerion = 0; ilowerion < nions - 1; ilowerion++) { ratesum += ion_ntion_energyrate(modelgridindex, ielement, ilowerion); if (ratesum >= zrand * ratetotal) { - *element = ielement; - *lowerion = ilowerion; - return; + return {ielement, ilowerion}; } } } @@ -1777,10 +1769,7 @@ void do_ntlepton(struct packet *pkt_ptr) { // const double frac_ionization = 0.; if (zrand < frac_ionization) { - int element = -1; - int lowerion = -1; - // select_nt_ionization(modelgridindex, &element, &lowerion); - select_nt_ionization2(modelgridindex, &element, &lowerion); + const auto [element, lowerion] = select_nt_ionization(modelgridindex); const int upperion = nt_random_upperion(modelgridindex, element, lowerion, true); // const int upperion = lowerion + 1; @@ -1819,8 +1808,8 @@ void do_ntlepton(struct packet *pkt_ptr) { zrand -= frac_ionization; // now zrand is between zero and frac_excitation // the selection algorithm is the same as for the ionization transitions - const int frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); - for (int excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { + const auto frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); + for (size_t excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { const double frac_deposition_exc = nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition; if (zrand < frac_deposition_exc) { @@ -2050,8 +2039,8 @@ static void analyse_sf_solution(const int modelgridindex, const int timestep, co const auto T_e = grid::get_Te(modelgridindex); printout(" Top non-thermal excitation fractions (total excitations = %d):\n", nt_solution[modelgridindex].frac_excitations_list.size()); - int ntransdisplayed = nt_solution[modelgridindex].frac_excitations_list.size(); - ntransdisplayed = (ntransdisplayed > 50) ? 50 : ntransdisplayed; + int ntransdisplayed = std::min(50, static_cast(nt_solution[modelgridindex].frac_excitations_list.size())); + for (excitationindex = 0; excitationindex < ntransdisplayed; excitationindex++) { const double frac_deposition = nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition; if (frac_deposition > 0.) { @@ -2563,45 +2552,44 @@ void write_restart_data(FILE *gridsave_file) { fprintf(gridsave_file, "%d\n", 24724518); // special number marking the beginning of NT data fprintf(gridsave_file, "%d %la %la\n", SFPTS, SF_EMIN, SF_EMAX); - for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { - if (grid::get_numassociatedcells(modelgridindex) > 0) { - fprintf(gridsave_file, "%d %d %la ", modelgridindex, deposition_rate_density_timestep[modelgridindex], - deposition_rate_density[modelgridindex]); + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + const int modelgridindex = grid::get_mgi_of_nonemptymgi(nonemptymgi); + fprintf(gridsave_file, "%d %d %la ", modelgridindex, deposition_rate_density_timestep[modelgridindex], + deposition_rate_density[modelgridindex]); - if (NT_ON && NT_SOLVE_SPENCERFANO) { - check_auger_probabilities(modelgridindex); + if (NT_ON && NT_SOLVE_SPENCERFANO) { + check_auger_probabilities(modelgridindex); - fprintf(gridsave_file, "%a %a %a %a\n", nt_solution[modelgridindex].nneperion_when_solved, - nt_solution[modelgridindex].frac_heating, nt_solution[modelgridindex].frac_ionization, - nt_solution[modelgridindex].frac_excitation); + fprintf(gridsave_file, "%a %a %a %a\n", nt_solution[modelgridindex].nneperion_when_solved, + nt_solution[modelgridindex].frac_heating, nt_solution[modelgridindex].frac_ionization, + nt_solution[modelgridindex].frac_excitation); - for (int uniqueionindex = 0; uniqueionindex < get_includedions(); uniqueionindex++) { - fprintf(gridsave_file, "%la ", nt_solution[modelgridindex].fracdep_ionization_ion[uniqueionindex]); - fprintf(gridsave_file, "%a ", nt_solution[modelgridindex].eff_ionpot[uniqueionindex]); + for (int uniqueionindex = 0; uniqueionindex < get_includedions(); uniqueionindex++) { + fprintf(gridsave_file, "%la ", nt_solution[modelgridindex].fracdep_ionization_ion[uniqueionindex]); + fprintf(gridsave_file, "%a ", nt_solution[modelgridindex].eff_ionpot[uniqueionindex]); - for (int a = 0; a <= NT_MAX_AUGER_ELECTRONS; a++) { - fprintf(gridsave_file, "%a %a ", - nt_solution[modelgridindex].prob_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a], - nt_solution[modelgridindex].ionenfrac_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a]); - } + for (int a = 0; a <= NT_MAX_AUGER_ELECTRONS; a++) { + fprintf(gridsave_file, "%a %a ", + nt_solution[modelgridindex].prob_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a], + nt_solution[modelgridindex].ionenfrac_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a]); } + } - // write NT excitations - fprintf(gridsave_file, "%d\n", static_cast(nt_solution[modelgridindex].frac_excitations_list.size())); + // write NT excitations + fprintf(gridsave_file, "%d\n", static_cast(nt_solution[modelgridindex].frac_excitations_list.size())); - const int frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); - for (int excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { - fprintf(gridsave_file, "%la %la %d\n", - nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition, - nt_solution[modelgridindex].frac_excitations_list[excitationindex].ratecoeffperdeposition, - nt_solution[modelgridindex].frac_excitations_list[excitationindex].lineindex); - } + const auto frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); + for (size_t excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { + fprintf(gridsave_file, "%la %la %d\n", + nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition, + nt_solution[modelgridindex].frac_excitations_list[excitationindex].ratecoeffperdeposition, + nt_solution[modelgridindex].frac_excitations_list[excitationindex].lineindex); + } - // write non-thermal spectrum - if (STORE_NT_SPECTRUM) { - for (int s = 0; s < SFPTS; s++) { - fprintf(gridsave_file, "%la\n", nt_solution[modelgridindex].yfunc[s]); - } + // write non-thermal spectrum + if (STORE_NT_SPECTRUM) { + for (int s = 0; s < SFPTS; s++) { + fprintf(gridsave_file, "%la\n", nt_solution[modelgridindex].yfunc[s]); } } } @@ -2630,60 +2618,58 @@ void read_restart_data(FILE *gridsave_file) { abort(); } - for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { - if (grid::get_numassociatedcells(modelgridindex) > 0) { - int mgi_in = 0; - assert_always(fscanf(gridsave_file, "%d %d %la ", &mgi_in, &deposition_rate_density_timestep[modelgridindex], - &deposition_rate_density[modelgridindex]) == 3); + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + const int modelgridindex = grid::get_mgi_of_nonemptymgi(nonemptymgi); - if (NT_ON && NT_SOLVE_SPENCERFANO) { - assert_always(fscanf(gridsave_file, "%a %a %a %a\n", &nt_solution[modelgridindex].nneperion_when_solved, - &nt_solution[modelgridindex].frac_heating, &nt_solution[modelgridindex].frac_ionization, - &nt_solution[modelgridindex].frac_excitation) == 4); + int mgi_in = 0; + assert_always(fscanf(gridsave_file, "%d %d %la ", &mgi_in, &deposition_rate_density_timestep[modelgridindex], + &deposition_rate_density[modelgridindex]) == 3); - if (mgi_in != modelgridindex) { - printout("ERROR: expected data for cell %d but found cell %d\n", modelgridindex, mgi_in); - abort(); - } + if (NT_ON && NT_SOLVE_SPENCERFANO) { + assert_always(fscanf(gridsave_file, "%a %a %a %a\n", &nt_solution[modelgridindex].nneperion_when_solved, + &nt_solution[modelgridindex].frac_heating, &nt_solution[modelgridindex].frac_ionization, + &nt_solution[modelgridindex].frac_excitation) == 4); - for (int uniqueionindex = 0; uniqueionindex < get_includedions(); uniqueionindex++) { - assert_always( - fscanf(gridsave_file, "%la ", &nt_solution[modelgridindex].fracdep_ionization_ion[uniqueionindex]) == 1); - assert_always(fscanf(gridsave_file, "%a ", &nt_solution[modelgridindex].eff_ionpot[uniqueionindex]) == 1); + if (mgi_in != modelgridindex) { + printout("ERROR: expected data for cell %d but found cell %d\n", modelgridindex, mgi_in); + abort(); + } - for (int a = 0; a <= NT_MAX_AUGER_ELECTRONS; a++) { - assert_always( - fscanf(gridsave_file, "%a %a ", - &nt_solution[modelgridindex].prob_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a], - &nt_solution[modelgridindex] - .ionenfrac_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a]) == 2); - } + for (int uniqueionindex = 0; uniqueionindex < get_includedions(); uniqueionindex++) { + assert_always( + fscanf(gridsave_file, "%la ", &nt_solution[modelgridindex].fracdep_ionization_ion[uniqueionindex]) == 1); + assert_always(fscanf(gridsave_file, "%a ", &nt_solution[modelgridindex].eff_ionpot[uniqueionindex]) == 1); + + for (int a = 0; a <= NT_MAX_AUGER_ELECTRONS; a++) { + assert_always( + fscanf(gridsave_file, "%a %a ", + &nt_solution[modelgridindex].prob_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a], + &nt_solution[modelgridindex] + .ionenfrac_num_auger[uniqueionindex * (NT_MAX_AUGER_ELECTRONS + 1) + a]) == 2); } + } - check_auger_probabilities(modelgridindex); + check_auger_probabilities(modelgridindex); - // read NT excitations - int frac_excitations_list_size_in = 0; - assert_always(fscanf(gridsave_file, "%d\n", &frac_excitations_list_size_in) == 1); + // read NT excitations + int frac_excitations_list_size_in = 0; + assert_always(fscanf(gridsave_file, "%d\n", &frac_excitations_list_size_in) == 1); - if (static_cast(nt_solution[modelgridindex].frac_excitations_list.size()) != - frac_excitations_list_size_in) { - nt_solution[modelgridindex].frac_excitations_list.resize(frac_excitations_list_size_in); - } + if (static_cast(nt_solution[modelgridindex].frac_excitations_list.size()) != frac_excitations_list_size_in) { + nt_solution[modelgridindex].frac_excitations_list.resize(frac_excitations_list_size_in); + } - for (int excitationindex = 0; excitationindex < frac_excitations_list_size_in; excitationindex++) { - assert_always( - fscanf(gridsave_file, "%la %la %d\n", - &nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition, - &nt_solution[modelgridindex].frac_excitations_list[excitationindex].ratecoeffperdeposition, - &nt_solution[modelgridindex].frac_excitations_list[excitationindex].lineindex) == 3); - } + for (int excitationindex = 0; excitationindex < frac_excitations_list_size_in; excitationindex++) { + assert_always(fscanf(gridsave_file, "%la %la %d\n", + &nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition, + &nt_solution[modelgridindex].frac_excitations_list[excitationindex].ratecoeffperdeposition, + &nt_solution[modelgridindex].frac_excitations_list[excitationindex].lineindex) == 3); + } - // read non-thermal spectrum - if (STORE_NT_SPECTRUM) { - for (int s = 0; s < SFPTS; s++) { - assert_always(fscanf(gridsave_file, "%la\n", &nt_solution[modelgridindex].yfunc[s]) == 1); - } + // read non-thermal spectrum + if (STORE_NT_SPECTRUM) { + for (int s = 0; s < SFPTS; s++) { + assert_always(fscanf(gridsave_file, "%la\n", &nt_solution[modelgridindex].yfunc[s]) == 1); } } } @@ -2721,16 +2707,16 @@ void nt_MPI_Bcast(const int modelgridindex, const int root) { MPI_FLOAT, root, MPI_COMM_WORLD); // communicate NT excitations - const int frac_excitations_list_size_old = nt_solution[modelgridindex].frac_excitations_list.size(); - int frac_excitations_list_size_new = nt_solution[modelgridindex].frac_excitations_list.size(); + const auto frac_excitations_list_size_old = nt_solution[modelgridindex].frac_excitations_list.size(); + auto frac_excitations_list_size_new = nt_solution[modelgridindex].frac_excitations_list.size(); MPI_Bcast(&frac_excitations_list_size_new, 1, MPI_INT, root, MPI_COMM_WORLD); if (frac_excitations_list_size_new != frac_excitations_list_size_old) { nt_solution[modelgridindex].frac_excitations_list.resize(frac_excitations_list_size_new); } - const int frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); - for (int excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { + const auto frac_excitations_list_size = nt_solution[modelgridindex].frac_excitations_list.size(); + for (size_t excitationindex = 0; excitationindex < frac_excitations_list_size; excitationindex++) { MPI_Bcast(&nt_solution[modelgridindex].frac_excitations_list[excitationindex].frac_deposition, 1, MPI_DOUBLE, root, MPI_COMM_WORLD); MPI_Bcast(&nt_solution[modelgridindex].frac_excitations_list[excitationindex].ratecoeffperdeposition, 1, diff --git a/nonthermal.h b/nonthermal.h index fbf2bb382..18093c485 100644 --- a/nonthermal.h +++ b/nonthermal.h @@ -9,16 +9,16 @@ namespace nonthermal { void init(int my_rank, int ndo_nonempty); void close_file(); void solve_spencerfano(int modelgridindex, int timestep, int iteration); -double nt_ionization_ratecoeff(int modelgridindex, int element, int ion); -double nt_ionization_upperion_probability(int modelgridindex, int element, int lowerion, int upperion, - bool energyweighted); -int nt_ionisation_maxupperion(int element, int lowerion); -int nt_random_upperion(int modelgridindex, int element, int lowerion, bool energyweighted); +auto nt_ionization_ratecoeff(int modelgridindex, int element, int ion) -> double; +auto nt_ionization_upperion_probability(int modelgridindex, int element, int lowerion, int upperion, + bool energyweighted) -> double; +auto nt_ionisation_maxupperion(int element, int lowerion) -> int; +auto nt_random_upperion(int modelgridindex, int element, int lowerion, bool energyweighted) -> int; void calculate_deposition_rate_density(int modelgridindex, int timestep); -double get_deposition_rate_density(int modelgridindex); -float get_nt_frac_heating(int modelgridindex); -double nt_excitation_ratecoeff(int modelgridindex, int element, int ion, int lowerlevel, int uptransindex, - double epsilon_trans, int lineindex); +auto get_deposition_rate_density(int modelgridindex) -> double; +auto get_nt_frac_heating(int modelgridindex) -> float; +auto nt_excitation_ratecoeff(int modelgridindex, int element, int ion, int lowerlevel, int uptransindex, + double epsilon_trans, int lineindex) -> double; void do_ntlepton(struct packet *pkt_ptr); void write_restart_data(FILE *gridsave_file); void read_restart_data(FILE *gridsave_file); diff --git a/packet.h b/packet.h index 4c282b833..477da1cfe 100644 --- a/packet.h +++ b/packet.h @@ -79,7 +79,7 @@ struct packet { float trueemissionvelocity = -1; struct mastate mastate; - inline bool operator==(const packet &rhs) { + inline auto operator==(const packet &rhs) -> bool { return (number == rhs.number && type == rhs.type && (em_pos[0] == rhs.em_pos[0] && em_pos[1] == rhs.em_pos[1] && em_pos[2] == rhs.em_pos[2]) && nu_cmf == rhs.nu_cmf && where == rhs.where && prop_time == rhs.prop_time && @@ -92,6 +92,6 @@ void packet_init(struct packet *pkt); void write_packets(char filename[], const struct packet *pkt); void read_packets(const char filename[], struct packet *pkt); void read_temp_packetsfile(int timestep, int my_rank, struct packet *pkt); -bool verify_temp_packetsfile(int timestep, int my_rank, const struct packet *pkt); +auto verify_temp_packetsfile(int timestep, int my_rank, const struct packet *pkt) -> bool; #endif // PACKET_H diff --git a/radfield.cc b/radfield.cc index bb117a7fc..520246688 100644 --- a/radfield.cc +++ b/radfield.cc @@ -8,12 +8,9 @@ #include #include #include -#include #include "atomic.h" #include "grid.h" -#include "ltepop.h" -#include "rpkt.h" #include "sn3d.h" #include "vectors.h" @@ -398,28 +395,28 @@ void init(int my_rank, int ndo_nonempty) /// Initialise estimator arrays which hold the last time steps values (used to damp out /// fluctuations over timestep iterations if DO_TITER is defined) to -1. void initialise_prev_titer_photoionestimators() { - for (int n = 0; n < grid::get_npts_model(); n++) { #ifdef DO_TITER - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); - J_reduced_save[nonemptymgi] = -1.; -#endif -#ifdef DO_TITER - nuJ_reduced_save[n] = -1.; + for (int nonemptymgi = 0; n < grid::get_npts_model(); n++) { globals::ffheatingestimator_save[n] = -1.; globals::colheatingestimator_save[n] = -1.; -#endif - for (int element = 0; element < get_nelements(); element++) { - const int nions = get_nions(element); - for (int ion = 0; ion < nions - 1; ion++) { -#ifdef DO_TITER - globals::gammaestimator_save[get_ionestimindex(n, element, ion)] = -1.; - if constexpr (USE_LUT_BFHEATING) { - globals::bfheatingestimator_save[get_ionestimindex(n, element, ion)] = -1.; + if (grid::get_numassociatedcells(modelgridindex) > 0) { + const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + J_reduced_save[nonemptymgi] = -1.; + nuJ_reduced_save[nonemptymgi] = -1.; + for (int element = 0; element < get_nelements(); element++) { + const int nions = get_nions(element); + for (int ion = 0; ion < nions - 1; ion++) { + if constexpr (USE_LUT_PHOTOION) { + globals::gammaestimator_save[get_ionestimindex(n, element, ion)] = -1.; + } + if constexpr (USE_LUT_BFHEATING) { + globals::bfheatingestimator_save[get_ionestimindex(n, element, ion)] = -1.; + } } -#endif } } } +#endif } auto get_Jblueindex(const int lineindex) -> int diff --git a/radfield.h b/radfield.h index 690f842c6..21bd1e879 100644 --- a/radfield.h +++ b/radfield.h @@ -16,15 +16,15 @@ void write_to_file(int modelgridindex, int timestep); void close_file(); void update_estimators(int modelgridindex, double distance_e_cmf, double nu_cmf, const struct packet *pkt_ptr); void update_lineestimator(int modelgridindex, int lineindex, double increment); -double radfield(double nu, int modelgridindex); +[[nodiscard]] auto radfield(double nu, int modelgridindex) -> double; void fit_parameters(int modelgridindex, int timestep); void set_J_normfactor(int modelgridindex, double normfactor); void normalise_J(int modelgridindex, double estimator_normfactor_over4pi); void normalise_nuJ(int modelgridindex, double estimator_normfactor_over4pi); -double get_T_J_from_J(int modelgridindex); -int get_Jblueindex(int lineindex); -double get_Jb_lu(int modelgridindex, int jblueindex); -int get_Jb_lu_contribcount(int modelgridindex, int jblueindex); +[[nodiscard]] auto get_T_J_from_J(int modelgridindex) -> double; +[[nodiscard]] auto get_Jblueindex(int lineindex) -> int; +[[nodiscard]] auto get_Jb_lu(int modelgridindex, int jblueindex) -> double; +[[nodiscard]] auto get_Jb_lu_contribcount(int modelgridindex, int jblueindex) -> int; void titer_J(int modelgridindex); void titer_nuJ(int modelgridindex); void reduce_estimators(); @@ -32,14 +32,14 @@ void do_MPI_Bcast(int modelgridindex, int root, int root_node_id); void write_restart_data(FILE *gridsave_file); void read_restart_data(FILE *gridsave_file); void normalise_bf_estimators(int modelgridindex, double estimator_normfactor_over_H); -double get_bfrate_estimator(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex); +auto get_bfrate_estimator(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex) -> double; void print_bfrate_contributions(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex, double nnlowerlevel, double nnlowerion); void reset_bfrate_contributions(int modelgridindex); -int integrate(const gsl_function *f, double nu_a, double nu_b, double epsabs, double epsrel, size_t limit, int key, - gsl_integration_workspace *workspace, double *result, double *abserr); +auto integrate(const gsl_function *f, double nu_a, double nu_b, double epsabs, double epsrel, size_t limit, int key, + gsl_integration_workspace *workspace, double *result, double *abserr) -> int; -constexpr double dbb(double nu, auto T, auto W) +[[nodiscard]] constexpr auto dbb(double nu, auto T, auto W) -> double // returns J_nu [ergs/s/sr/cm2/Hz] for a dilute black body with temperature T and dilution factor W { return W * TWOHOVERCLIGHTSQUARED * std::pow(nu, 3) / std::expm1(HOVERKB * nu / T); diff --git a/ratecoeff.cc b/ratecoeff.cc index a81ec589d..bbdf776d8 100644 --- a/ratecoeff.cc +++ b/ratecoeff.cc @@ -2,10 +2,9 @@ #include -#include #include #include -#define D_POSIX_SOURCE +// #define D_POSIX_SOURCE #include #include @@ -13,11 +12,11 @@ #include "artisoptions.h" #include "atomic.h" #include "grid.h" -#include "input.h" #include "ltepop.h" #include "macroatom.h" #include "md5.h" #include "radfield.h" +#include "rpkt.h" #include "sn3d.h" // typedef struct gslintegration_ffheatingparas @@ -1322,13 +1321,14 @@ auto get_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex } else { const double W = grid::get_W(modelgridindex); const double T_R = grid::get_TR(modelgridindex); + const auto nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); gammacorr = W * interpolate_corrphotoioncoeff(element, ion, level, phixstargetindex, T_R); const int index_in_groundlevelcontestimator = globals::elements[element].ions[ion].levels[level].closestgroundlevelcont; if (index_in_groundlevelcontestimator >= 0) { gammacorr *= - globals::corrphotoionrenorm[get_ionestimindex(modelgridindex, 0, 0) + index_in_groundlevelcontestimator]; + globals::corrphotoionrenorm[nonemptymgi * get_includedions() + index_in_groundlevelcontestimator]; } } } diff --git a/ratecoeff.h b/ratecoeff.h index aa2df0d3b..46eaeb44e 100644 --- a/ratecoeff.h +++ b/ratecoeff.h @@ -3,30 +3,30 @@ void ratecoefficients_init(); -void setup_photoion_luts(void); +void setup_photoion_luts(); -double select_continuum_nu(int element, int lowerion, int lower, int upperionlevel, float T_e); +auto select_continuum_nu(int element, int lowerion, int lower, int upperionlevel, float T_e) -> double; -double interpolate_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, double T); +auto interpolate_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, double T) -> double; -double get_spontrecombcoeff(int element, int ion, int level, int phixstargetindex, float T_e); -double get_stimrecombcoeff(int element, int lowerion, int level, int phixstargetindex, int modelgridindex); +auto get_spontrecombcoeff(int element, int ion, int level, int phixstargetindex, float T_e) -> double; +auto get_stimrecombcoeff(int element, int lowerion, int level, int phixstargetindex, int modelgridindex) -> double; -double get_bfcoolingcoeff(int element, int ion, int level, int phixstargetindex, float T_e); +auto get_bfcoolingcoeff(int element, int ion, int level, int phixstargetindex, float T_e) -> double; -double get_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, int modelgridindex); -double get_corrphotoioncoeff_ana(int element, int ion, int level, int phixstargetindex, int modelgridindex); +auto get_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, int modelgridindex) -> double; +auto get_corrphotoioncoeff_ana(int element, int ion, int level, int phixstargetindex, int modelgridindex) -> double; -bool iongamma_is_zero(int modelgridindex, int element, int ion); +auto iongamma_is_zero(int modelgridindex, int element, int ion) -> bool; -double calculate_iongamma_per_gspop(int modelgridindex, int element, int ion); -double calculate_iongamma_per_ionpop(int modelgridindex, float T_e, int element, int lowerion, bool assume_lte, - bool collisional_not_radiative, bool printdebug, bool force_bfest, - bool force_bfintegral); +auto calculate_iongamma_per_gspop(int modelgridindex, int element, int ion) -> double; +auto calculate_iongamma_per_ionpop(int modelgridindex, float T_e, int element, int lowerion, bool assume_lte, + bool collisional_not_radiative, bool printdebug, bool force_bfest, + bool force_bfintegral) -> double; -double calculate_ionrecombcoeff(int modelgridindex, float T_e, int element, int upperion, bool assume_lte, - bool collisional_not_radiative, bool printdebug, bool lower_superlevel_only, - bool per_groundmultipletpop, bool stimonly); +auto calculate_ionrecombcoeff(int modelgridindex, float T_e, int element, int upperion, bool assume_lte, + bool collisional_not_radiative, bool printdebug, bool lower_superlevel_only, + bool per_groundmultipletpop, bool stimonly) -> double; // CUDA integration. might be worth revisting for SYCL support // template diff --git a/rpkt.cc b/rpkt.cc index 066980e99..3d06fe5f6 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -5,15 +5,14 @@ #include #include +#include "artisoptions.h" #include "atomic.h" +#include "globals.h" #include "grid.h" -#include "kpkt.h" #include "ltepop.h" -#include "macroatom.h" #include "radfield.h" #include "sn3d.h" #include "stats.h" -#include "update_grid.h" #include "vectors.h" #include "vpkt.h" @@ -535,13 +534,11 @@ static void rpkt_event_thickcell(struct packet *pkt_ptr) pkt_ptr->em_time = pkt_ptr->prop_time; } -static void update_estimators(const struct packet *pkt_ptr, const double distance) +static void update_estimators(const struct packet *pkt_ptr, const double distance, const int modelgridindex) /// Update the volume estimators J and nuJ /// This is done in another routine than move, as we sometimes move dummy /// packets which do not contribute to the radiation field. { - const int modelgridindex = grid::get_cell_modelgridindex(pkt_ptr->where); - /// Update only non-empty cells if (modelgridindex == grid::get_npts_model()) { return; @@ -556,8 +553,7 @@ static void update_estimators(const struct packet *pkt_ptr, const double distanc if constexpr (USE_LUT_PHOTOION || USE_LUT_BFHEATING) { const double distance_e_cmf_over_nu = distance_e_cmf / nu; - const int nelements = get_nelements(); - const int max_nions = get_max_nions(); + const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); for (int i = 0; i < globals::nbfcontinua_ground; i++) { const double nu_edge = globals::groundcont[i].nu_edge; @@ -568,7 +564,7 @@ static void update_estimators(const struct packet *pkt_ptr, const double distanc /// the estimators if (grid::get_elem_abundance(modelgridindex, element) > 0) { const int ion = globals::groundcont[i].ion; - const int ionestimindex = modelgridindex * nelements * max_nions + element * max_nions + ion; + const int ionestimindex = get_ionestimindex_nonemptymgi(nonemptymgi, element, ion); if constexpr (USE_LUT_PHOTOION) { safeadd(globals::gammaestimator[ionestimindex], @@ -695,7 +691,7 @@ static auto do_rpkt_step(struct packet *pkt_ptr, const double t2) -> bool if ((sdist < tdist) && (sdist < edist)) { // Move it into the new cell. move_pkt_withtime(pkt_ptr, sdist / 2.); - update_estimators(pkt_ptr, sdist); + update_estimators(pkt_ptr, sdist, mgi); move_pkt_withtime(pkt_ptr, sdist / 2.); if (snext != pkt_ptr->where) { @@ -712,7 +708,7 @@ static auto do_rpkt_step(struct packet *pkt_ptr, const double t2) -> bool if ((edist < sdist) && (edist < tdist)) { // bound-bound or continuum event move_pkt_withtime(pkt_ptr, edist / 2.); - update_estimators(pkt_ptr, edist); + update_estimators(pkt_ptr, edist, mgi); move_pkt_withtime(pkt_ptr, edist / 2.); // The previously selected and in pkt_ptr stored event occurs. Handling is done by rpkt_event @@ -732,7 +728,7 @@ static auto do_rpkt_step(struct packet *pkt_ptr, const double t2) -> bool if ((tdist < sdist) && (tdist < edist)) { // reaches end of timestep before cell boundary or interaction move_pkt_withtime(pkt_ptr, tdist / 2.); - update_estimators(pkt_ptr, tdist); + update_estimators(pkt_ptr, tdist, mgi); pkt_ptr->prop_time = t2; move_pkt(pkt_ptr, tdist / 2.); pkt_ptr->last_event = pkt_ptr->last_event + 1000; diff --git a/rpkt.h b/rpkt.h index 6b4db7593..55f5ea3e3 100644 --- a/rpkt.h +++ b/rpkt.h @@ -3,11 +3,13 @@ #include "artisoptions.h" #include "constants.h" +#include "grid.h" +#include "sn3d.h" void do_rpkt(struct packet *pkt_ptr, double t2); void emit_rpkt(struct packet *pkt_ptr); -int closest_transition(double nu_cmf, int next_trans); -double calculate_chi_bf_gammacontr(int modelgridindex, double nu); +auto closest_transition(double nu_cmf, int next_trans) -> int; +auto calculate_chi_bf_gammacontr(int modelgridindex, double nu) -> double; void calculate_chi_rpkt_cont(double nu_cmf, struct rpkt_continuum_absorptioncoeffs *chi_rpkt_cont_thisthread, int modelgridindex, bool usecellhistupdatephixslist); @@ -30,4 +32,16 @@ constexpr auto get_linedistance(const double prop_time, const double nu_cmf, con return CLIGHT * prop_time * (nu_cmf / nu_trans - 1); } + +[[nodiscard]] inline auto get_ionestimindex_nonemptymgi(const int nonemptymgi, const int element, const int ion) + -> int { + assert_testmodeonly(ion >= 0); + assert_testmodeonly(ion < get_nions(element) - 1); + return nonemptymgi * get_includedions() + get_uniqueionindex(element, ion); +} + +[[nodiscard]] inline auto get_ionestimindex(const int mgi, const int element, const int ion) -> int { + return get_ionestimindex_nonemptymgi(grid::get_modelcell_nonemptymgi(mgi), element, ion); +} + #endif // RPKT_H diff --git a/scripts/clean.sh b/scripts/clean.sh index 33537f161..adc9b1e61 100755 --- a/scripts/clean.sh +++ b/scripts/clean.sh @@ -1,6 +1,6 @@ #!/usr/bin/env bash -paths="*.tmp *.out *.out.* out.txt output_*-?.txt exspec.txt machine.file.* core.* *.slurm packets bflist.dat logfiles.tar*" +paths="*.tmp *.out *.out.* out.txt output_*-*.txt exspec.txt machine.file.* core.* *.slurm packets bflist.dat logfiles.tar*" # if [[ "$1" == "-d" ]]; then # echo 1 diff --git a/sn3d.cc b/sn3d.cc index 0125d0279..40c249a8f 100644 --- a/sn3d.cc +++ b/sn3d.cc @@ -14,6 +14,7 @@ #include +#include "artisoptions.h" #include "atomic.h" #include "decay.h" #include "gammapkt.h" @@ -24,6 +25,7 @@ #include "nonthermal.h" #include "radfield.h" #include "ratecoeff.h" +#include "rpkt.h" #include "spectrum.h" #include "stats.h" #include "update_grid.h" @@ -42,7 +44,7 @@ static time_t real_time_start = -1; static time_t time_timestep_start = -1; // this will be set after the first update of the grid and before packet prop static FILE *estimators_file = nullptr; -int mpi_grid_buffer_size = 0; +size_t mpi_grid_buffer_size = 0; char *mpi_grid_buffer = nullptr; static void initialise_linestat_file() { @@ -189,7 +191,7 @@ static void write_deposition_file(const int nts, const int my_rank, const int ns #ifdef MPI_ON static void mpi_communicate_grid_properties(const int my_rank, const int nprocs, const int nstart, const int ndo, - char *mpi_grid_buffer, const int mpi_grid_buffer_size) { + char *mpi_grid_buffer, const size_t mpi_grid_buffer_size) { int position = 0; for (int root = 0; root < nprocs; root++) { MPI_Barrier(MPI_COMM_WORLD); @@ -201,23 +203,35 @@ static void mpi_communicate_grid_properties(const int my_rank, const int nprocs, MPI_Bcast(&root_node_id, 1, MPI_INT, root, MPI_COMM_WORLD); for (int modelgridindex = root_nstart; modelgridindex < (root_nstart + root_ndo); modelgridindex++) { + if (grid::get_numassociatedcells(modelgridindex) < 1) { + continue; + } + radfield::do_MPI_Bcast(modelgridindex, root, root_node_id); - if (grid::get_numassociatedcells(modelgridindex) > 0) { - nonthermal::nt_MPI_Bcast(modelgridindex, root); - if (globals::total_nlte_levels > 0 && globals::rank_in_node == 0) { - MPI_Bcast(grid::modelgrid[modelgridindex].nlte_pops, globals::total_nlte_levels, MPI_DOUBLE, root_node_id, - globals::mpi_comm_internode); - } + nonthermal::nt_MPI_Bcast(modelgridindex, root); + if (globals::total_nlte_levels > 0 && globals::rank_in_node == 0) { + MPI_Bcast(grid::modelgrid[modelgridindex].nlte_pops, globals::total_nlte_levels, MPI_DOUBLE, root_node_id, + globals::mpi_comm_internode); + } - if constexpr (USE_LUT_PHOTOION) { - assert_always(globals::corrphotoionrenorm != nullptr); - MPI_Bcast(&globals::corrphotoionrenorm[get_ionestimindex(modelgridindex, 0, 0)], - get_nelements() * get_max_nions(), MPI_DOUBLE, root, MPI_COMM_WORLD); - assert_always(globals::gammaestimator != nullptr); - MPI_Bcast(&globals::gammaestimator[get_ionestimindex(modelgridindex, 0, 0)], - get_nelements() * get_max_nions(), MPI_DOUBLE, root, MPI_COMM_WORLD); + if constexpr (USE_LUT_PHOTOION) { + const auto nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + assert_always(globals::corrphotoionrenorm != nullptr); + if (globals::rank_in_node == 0) { + MPI_Bcast(&globals::corrphotoionrenorm[nonemptymgi * get_includedions()], get_includedions(), MPI_DOUBLE, + root_node_id, globals::mpi_comm_internode); } + + assert_always(globals::gammaestimator != nullptr); + MPI_Bcast(&globals::gammaestimator[nonemptymgi * get_includedions()], get_includedions(), MPI_DOUBLE, root, + MPI_COMM_WORLD); + } + + assert_always(grid::modelgrid[modelgridindex].elem_meanweight != nullptr); + if (globals::rank_in_node == 0) { + MPI_Bcast(grid::modelgrid[modelgridindex].elem_meanweight, get_nelements(), MPI_FLOAT, root_node_id, + globals::mpi_comm_internode); } } @@ -250,18 +264,20 @@ static void mpi_communicate_grid_properties(const int my_rank, const int nprocs, MPI_COMM_WORLD); for (int element = 0; element < get_nelements(); element++) { - MPI_Pack(grid::modelgrid[mgi].composition[element].groundlevelpop, get_nions(element), MPI_FLOAT, - mpi_grid_buffer, mpi_grid_buffer_size, &position, MPI_COMM_WORLD); - MPI_Pack(grid::modelgrid[mgi].composition[element].partfunct, get_nions(element), MPI_FLOAT, - mpi_grid_buffer, mpi_grid_buffer_size, &position, MPI_COMM_WORLD); - MPI_Pack(grid::modelgrid[mgi].cooling_contrib_ion[element], get_nions(element), MPI_DOUBLE, mpi_grid_buffer, - mpi_grid_buffer_size, &position, MPI_COMM_WORLD); + if (get_nions(element) > 0) { + MPI_Pack(grid::modelgrid[mgi].composition[element].groundlevelpop, get_nions(element), MPI_FLOAT, + mpi_grid_buffer, mpi_grid_buffer_size, &position, MPI_COMM_WORLD); + MPI_Pack(grid::modelgrid[mgi].composition[element].partfunct, get_nions(element), MPI_FLOAT, + mpi_grid_buffer, mpi_grid_buffer_size, &position, MPI_COMM_WORLD); + MPI_Pack(grid::modelgrid[mgi].cooling_contrib_ion[element], get_nions(element), MPI_DOUBLE, + mpi_grid_buffer, mpi_grid_buffer_size, &position, MPI_COMM_WORLD); + } } } } printout("[info] mem_usage: MPI_BUFFER: used %d of %d bytes allocated to mpi_grid_buffer\n", position, mpi_grid_buffer_size); - assert_always(position <= mpi_grid_buffer_size); + assert_always(static_cast(position) <= mpi_grid_buffer_size); } MPI_Barrier(MPI_COMM_WORLD); MPI_Bcast(mpi_grid_buffer, mpi_grid_buffer_size, MPI_PACKED, root, MPI_COMM_WORLD); @@ -297,14 +313,17 @@ static void mpi_communicate_grid_properties(const int my_rank, const int nprocs, MPI_COMM_WORLD); for (int element = 0; element < get_nelements(); element++) { - MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, - grid::modelgrid[mgi].composition[element].groundlevelpop, get_nions(element), MPI_FLOAT, - MPI_COMM_WORLD); - MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, - grid::modelgrid[mgi].composition[element].partfunct, get_nions(element), MPI_FLOAT, - MPI_COMM_WORLD); - MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, - grid::modelgrid[mgi].cooling_contrib_ion[element], get_nions(element), MPI_DOUBLE, MPI_COMM_WORLD); + if (get_nions(element) > 0) { + MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, + grid::modelgrid[mgi].composition[element].groundlevelpop, get_nions(element), MPI_FLOAT, + MPI_COMM_WORLD); + MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, + grid::modelgrid[mgi].composition[element].partfunct, get_nions(element), MPI_FLOAT, + MPI_COMM_WORLD); + MPI_Unpack(mpi_grid_buffer, mpi_grid_buffer_size, &position, + grid::modelgrid[mgi].cooling_contrib_ion[element], get_nions(element), MPI_DOUBLE, + MPI_COMM_WORLD); + } } } } @@ -320,14 +339,17 @@ static void mpi_reduce_estimators(int nts) { MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD); - const int arraylen = grid::get_npts_model() * get_nelements() * get_max_nions(); + const int arraylen = grid::get_nonempty_npts_model() * get_includedions(); + if constexpr (USE_LUT_PHOTOION) { MPI_Barrier(MPI_COMM_WORLD); assert_always(globals::gammaestimator != nullptr); MPI_Allreduce(MPI_IN_PLACE, globals::gammaestimator, arraylen, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } + if constexpr (USE_LUT_BFHEATING) { MPI_Barrier(MPI_COMM_WORLD); + assert_always(globals::bfheatingestimator != nullptr); MPI_Allreduce(MPI_IN_PLACE, globals::bfheatingestimator, arraylen, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } @@ -518,32 +540,36 @@ static void save_grid_and_packets(const int nts, const int my_rank, struct packe } static void zero_estimators() { - // printout("zero_estimators()"); - for (int n = 0; n < grid::get_npts_model(); n++) { - if (grid::get_numassociatedcells(n) > 0) { - radfield::zero_estimators(n); +#ifdef MPI_ON + MPI_Barrier(MPI_COMM_WORLD); +#endif + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + const auto modelgridindex = grid::get_mgi_of_nonemptymgi(nonemptymgi); + radfield::zero_estimators(modelgridindex); - globals::ffheatingestimator[n] = 0.; - globals::colheatingestimator[n] = 0.; + globals::ffheatingestimator[modelgridindex] = 0.; + globals::colheatingestimator[modelgridindex] = 0.; - if constexpr (TRACK_ION_STATS) { - stats::reset_ion_stats(n); - } + if constexpr (TRACK_ION_STATS) { + stats::reset_ion_stats(modelgridindex); + } - for (int element = 0; element < get_nelements(); element++) { - for (int ion = 0; ion < get_max_nions(); ion++) { - if constexpr (USE_LUT_PHOTOION) { - globals::gammaestimator[get_ionestimindex(n, element, ion)] = 0.; - } - if constexpr (USE_LUT_BFHEATING) { - globals::bfheatingestimator[get_ionestimindex(n, element, ion)] = 0.; - } + for (int element = 0; element < get_nelements(); element++) { + for (int ion = 0; ion < (get_nions(element) - 1); ion++) { + if constexpr (USE_LUT_PHOTOION) { + globals::gammaestimator[get_ionestimindex(modelgridindex, element, ion)] = 0.; + } + if constexpr (USE_LUT_BFHEATING) { + globals::bfheatingestimator[get_ionestimindex(modelgridindex, element, ion)] = 0.; } } - - globals::rpkt_emiss[n] = 0.0; } + + globals::rpkt_emiss[modelgridindex] = 0.; } +#ifdef MPI_ON + MPI_Barrier(MPI_COMM_WORLD); +#endif } static auto do_timestep(const int nts, const int titer, const int my_rank, const int nstart, const int ndo, @@ -599,7 +625,9 @@ static auto do_timestep(const int nts, const int titer, const int my_rank, const // and also the photoion and stimrecomb estimators zero_estimators(); - // MPI_Barrier(MPI_COMM_WORLD); +#ifdef MPI_ON + MPI_Barrier(MPI_COMM_WORLD); +#endif if ((nts < globals::timestep_finish) && do_this_full_loop) { /// Now process the packets. @@ -732,9 +760,9 @@ auto main(int argc, char *argv[]) -> int { globals::startofline[tid] = true; #ifdef _OPENMP - printout("OpenMP parallelisation active with %d threads (max %d)\n", get_num_threads(), get_max_threads()); + printout("OpenMP parallelisation is active with %d threads (max %d)\n", get_num_threads(), get_max_threads()); #else - printout("OpenMP is not enabled in this build (this is normal)\n"); + printout("OpenMP parallelisation is not enabled in this build (this is normal)\n"); #endif gslworkspace = gsl_integration_workspace_alloc(GSLWSIZE); @@ -749,7 +777,7 @@ auto main(int argc, char *argv[]) -> int { #endif int opt = 0; - while ((opt = getopt(argc, argv, "w:")) != -1) { + while ((opt = getopt(argc, argv, "w:")) != -1) { // NOLINT(concurrency-mt-unsafe) if (opt == 'w') { printout("Command line argument specifies wall time hours '%s', setting ", optarg); const float walltimehours = strtof(optarg, nullptr); @@ -873,7 +901,7 @@ auto main(int argc, char *argv[]) -> int { #ifdef MPI_ON MPI_Barrier(MPI_COMM_WORLD); - const int maxndo = grid::get_maxndo(); + const size_t maxndo = grid::get_maxndo(); /// Initialise the exchange buffer /// The factor 4 comes from the fact that our buffer should contain elements of 4 byte /// instead of 1 byte chars. But the MPI routines don't care about the buffers datatype diff --git a/sn3d.h b/sn3d.h index 09067bba3..610bf1547 100644 --- a/sn3d.h +++ b/sn3d.h @@ -15,9 +15,6 @@ #include #include -#include "artisoptions.h" -#include "globals.h" - // #define _OPENMP #ifdef _OPENMP #include @@ -69,35 +66,38 @@ extern gsl_integration_workspace *gslworkspace; } #endif -// #define printout(...) fprintf(output_file, __VA_ARGS__) +#include "artisoptions.h" +#include "globals.h" -extern int tid; +// #define printout(...) fprintf(output_file, __VA_ARGS__) template -static int printout(const char *format, Args... args) { +static auto printout(const char *format, Args... args) -> int { if (globals::startofline[tid]) { const time_t now_time = time(nullptr); char s[32] = ""; - strftime(s, 32, "%FT%TZ", gmtime(&now_time)); + struct tm buf {}; + strftime(s, 32, "%FT%TZ", gmtime_r(&now_time, &buf)); fprintf(output_file, "%s ", s); } globals::startofline[tid] = (format[strlen(format) - 1] == '\n'); return fprintf(output_file, format, args...); } -static int printout(const char *format) { +static auto printout(const char *format) -> int { if (globals::startofline[tid]) { const time_t now_time = time(nullptr); char s[32] = ""; - strftime(s, 32, "%FT%TZ", gmtime(&now_time)); + struct tm buf {}; + strftime(s, 32, "%FT%TZ", gmtime_r(&now_time, &buf)); fprintf(output_file, "%s ", s); } globals::startofline[tid] = (format[strlen(format) - 1] == '\n'); return fprintf(output_file, "%s", format); } -static inline int get_bflutindex(const int tempindex, const int element, const int ion, const int level, - const int phixstargetindex) { +static inline auto get_bflutindex(const int tempindex, const int element, const int ion, const int level, + const int phixstargetindex) -> int { const int contindex = -1 - globals::elements[element].ions[ion].levels[level].cont_index + phixstargetindex; const int bflutindex = tempindex * globals::nbfcontinua + contindex; @@ -105,13 +105,27 @@ static inline int get_bflutindex(const int tempindex, const int element, const i return bflutindex; } +template +inline void safeadd(T &var, T val) { #ifdef _OPENMP -#define safeadd(var, val) _Pragma("omp atomic update") var += val +#pragma omp atomic update + var += val; +#else +#ifdef STDPAR_ON +#ifdef __cpp_lib_atomic_ref + static_assert(std::atomic::is_always_lock_free); + std::atomic_ref(var).fetch_add(val, std::memory_order_relaxed); #else -#define safeadd(var, val) var = var + val + // this works on clang but not gcc for doubles. + __atomic_fetch_add(&var, val, __ATOMIC_RELAXED); #endif +#else + var += val; +#endif +#endif +} -#define safeincrement(var) safeadd(var, 1) +#define safeincrement(var) safeadd((var), 1) // #define DO_TITER @@ -123,7 +137,7 @@ static inline void gsl_error_handler_printout(const char *reason, const char *fi } } -static FILE *fopen_required(const std::string &filename, const char *mode) { +static auto fopen_required(const std::string &filename, const char *mode) -> FILE * { // look in the data folder first const std::string datafolderfilename = "data/" + filename; if (mode[0] == 'r' && std::filesystem::exists(datafolderfilename)) { @@ -139,7 +153,7 @@ static FILE *fopen_required(const std::string &filename, const char *mode) { return file; } -static std::fstream fstream_required(const std::string &filename, std::ios_base::openmode mode) { +static auto fstream_required(const std::string &filename, std::ios_base::openmode mode) -> std::fstream { const std::string datafolderfilename = "data/" + filename; if (mode == std::ios::in && std::filesystem::exists(datafolderfilename)) { return fstream_required(datafolderfilename, mode); @@ -152,7 +166,7 @@ static std::fstream fstream_required(const std::string &filename, std::ios_base: return file; } -static int get_timestep(const double time) { +static auto get_timestep(const double time) -> int { assert_always(time >= globals::tmin); assert_always(time < globals::tmax); for (int nts = 0; nts < globals::ntimesteps; nts++) { @@ -166,7 +180,7 @@ static int get_timestep(const double time) { return -1; } -inline int get_max_threads(void) { +inline auto get_max_threads() -> int { #if defined _OPENMP return omp_get_max_threads(); #else @@ -174,7 +188,7 @@ inline int get_max_threads(void) { #endif } -inline int get_num_threads(void) { +inline auto get_num_threads() -> int { #if defined _OPENMP return omp_get_num_threads(); #else @@ -182,7 +196,7 @@ inline int get_num_threads(void) { #endif } -inline int get_thread_num(void) { +inline auto get_thread_num() -> int { #if defined _OPENMP return omp_get_thread_num(); #else @@ -190,19 +204,19 @@ inline int get_thread_num(void) { #endif } -inline float rng_uniform(void) { - float zrand; - do { - zrand = std::generate_canonical::digits>(stdrng); - } while (zrand == 1.); +inline auto rng_uniform() -> float { + const auto zrand = std::generate_canonical::digits>(stdrng); + if (zrand == 1.) { + return rng_uniform(); + } return zrand; } -inline float rng_uniform_pos(void) { - float zrand = 0.; - do { - zrand = rng_uniform(); - } while (zrand <= 0.); +inline auto rng_uniform_pos() -> float { + const auto zrand = std::generate_canonical::digits>(stdrng); + if (zrand <= 0.) { + return rng_uniform_pos(); + } return zrand; } @@ -211,22 +225,20 @@ inline void rng_init(const uint_fast64_t zseed) { stdrng.seed(zseed); } -inline bool is_pid_running(pid_t pid) { - while (waitpid(-1, 0, WNOHANG) > 0) { +inline auto is_pid_running(pid_t pid) -> bool { + while (waitpid(-1, nullptr, WNOHANG) > 0) { // Wait for defunct.... } - if (0 == kill(pid, 0)) return true; // Process exists - - return false; + return (0 == kill(pid, 0)); } -inline void check_already_running(void) { +inline void check_already_running() { pid_t artispid = getpid(); if (std::filesystem::exists("artis.pid")) { auto pidfile = std::fstream("artis.pid", std::ios::in); - pid_t artispid_in; + pid_t artispid_in = 0; pidfile >> artispid_in; pidfile.close(); if (is_pid_running(artispid_in)) { @@ -244,8 +256,4 @@ inline void check_already_running(void) { pidfile.close(); } -inline int get_ionestimindex(const int mgi, const int element, const int ion) { - return mgi * get_nelements() * get_max_nions() + element * get_max_nions() + ion; -} - #endif // SN3D_H diff --git a/stats.cc b/stats.cc index f991c087d..f584e2ec8 100644 --- a/stats.cc +++ b/stats.cc @@ -1,20 +1,22 @@ #include "stats.h" #include -#include #include "atomic.h" #include "globals.h" #include "grid.h" #include "ltepop.h" #include "nonthermal.h" +#include "sn3d.h" namespace stats { static double *ionstats = nullptr; -static std::array, COUNTER_COUNT> eventstats; +static std::vector> eventstats; void init() { + eventstats.resize(get_max_threads(), {0}); + if constexpr (TRACK_ION_STATS) { ionstats = static_cast(malloc(grid::get_npts_model() * get_includedions() * ION_STAT_COUNT * sizeof(double))); @@ -158,12 +160,12 @@ void normalise_ion_estimators(const int mgi, const double deltat, const double d void increment(enum eventcounters i) { assert_testmodeonly(i >= 0); assert_testmodeonly(i < COUNTER_COUNT); - eventstats[i]++; + eventstats[tid][i]++; } void pkt_action_counters_reset() { for (int i = 0; i < COUNTER_COUNT; i++) { - eventstats[i] = 0; + eventstats[tid][i] = 0; } nonthermal::nt_reset_stats(); @@ -172,7 +174,11 @@ void pkt_action_counters_reset() { auto get_counter(enum eventcounters i) -> int { assert_always(i < COUNTER_COUNT); - return eventstats[i]; + int count = 0; + for (int t = 0; t < get_num_threads(); t++) { + count += eventstats[t][i]; + } + return count; } void pkt_action_counters_printout(const struct packet *const pkt, const int nts) { diff --git a/stats.h b/stats.h index 2649be66c..fce850f86 100644 --- a/stats.h +++ b/stats.h @@ -89,7 +89,7 @@ void increment_ion_stats(int modelgridindex, int element, int ion, enum ionstatt void increment_ion_stats_contabsorption(const struct packet *pkt_ptr, int modelgridindex, int element, int ion); -double get_ion_stats(int modelgridindex, int element, int ion, enum ionstattypes ionstattype); +auto get_ion_stats(int modelgridindex, int element, int ion, enum ionstattypes ionstattype) -> double; void set_ion_stats(int modelgridindex, int element, int ion, enum ionstattypes ionstattype, double newvalue); @@ -101,7 +101,7 @@ void increment(enum eventcounters); void pkt_action_counters_reset(); -int get_counter(enum eventcounters i); +auto get_counter(enum eventcounters i) -> int; void pkt_action_counters_printout(const struct packet *pkt, int nts); diff --git a/tests/.gitignore b/tests/.gitignore index 05570f648..01c4d6f08 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -5,7 +5,6 @@ adata.txt bflist.dat *.tmp ratecoeff.dat -phixsdata_v2.txt transitiondata.txt *.tar.* atomic_data_logs @@ -19,3 +18,4 @@ binding_energies.txt *_testrun atomicdata_* *.pdf +checksums.zip diff --git a/tests/classicmode_1d_3dgrid_inputfiles/results_md5_final.txt b/tests/classicmode_1d_3dgrid_inputfiles/results_md5_final.txt index 432d20e55..8dd304787 100644 --- a/tests/classicmode_1d_3dgrid_inputfiles/results_md5_final.txt +++ b/tests/classicmode_1d_3dgrid_inputfiles/results_md5_final.txt @@ -21,5 +21,5 @@ bbf9c836a1bbc64f7919d8c2dccb6b52 specpol.out b2a6c81cfd8605508b89ed47f1f83be6 vpkt_grid_1-0.out 8d19f1a64dd5cb69869048df175a8a02 vspecpol_0-0.out 804848d2d65244f3c5d481db0384213b vspecpol_1-0.out -263568b30ce0c2633fbec2157aaa1114 job1/estimators_0000.out -b2ebf78e08410d7c7562d0d06120c95e job1/estimators_0001.out +c17b719d51868938ea2569ee03526ed6 job1/estimators_0000.out +83f41b85cdb98c088c001619135cde42 job1/estimators_0001.out diff --git a/tests/classicmode_1d_3dgrid_inputfiles/results_md5_job0.txt b/tests/classicmode_1d_3dgrid_inputfiles/results_md5_job0.txt index 8e9004e13..adff299e0 100644 --- a/tests/classicmode_1d_3dgrid_inputfiles/results_md5_job0.txt +++ b/tests/classicmode_1d_3dgrid_inputfiles/results_md5_job0.txt @@ -17,5 +17,5 @@ c4e4e8d00846618f3931dc01cf52615f modelgridrankassignments.out 35ce3462e5b9493bde385b8bad9ea954 vpkt_grid_1-0.out 2ed5cc4ea567b624eb8a4cb9137a2e01 vspecpol_0-0.out 01fd53e3ca382c4e2b0f5b588edbd9e5 vspecpol_1-0.out -2c67a882e73f240ecfeffa05f04d8c5f job0/estimators_0000.out -33ea538a70082b2a2ae43ae71569fc3e job0/estimators_0001.out +3779b9203782f26c2568c7e9caaa17c1 job0/estimators_0000.out +bee1cab6351770f62701c9fe3f9b7c33 job0/estimators_0001.out diff --git a/tests/classicmode_3d_inputfiles/results_md5_final.txt b/tests/classicmode_3d_inputfiles/results_md5_final.txt index 4fc5b78cb..28605dfd4 100644 --- a/tests/classicmode_3d_inputfiles/results_md5_final.txt +++ b/tests/classicmode_3d_inputfiles/results_md5_final.txt @@ -520,5 +520,5 @@ d88b1007563f469de4c3fa6485354950 spec.out 0463bd17f8c10622f2f13c1d4da2fb9f specpol.out edffea9cc8b28c23286d9438a28f7924 specpol_res.out bc01f046eab9bf3802149e0ff63d5069 timesteps.out -6ea207f7cffa562e62f9d2311376ed30 job1/estimators_0000.out -676fc28b1c1145b219cf3d7a93b2d589 job1/estimators_0001.out +2856bee28f5dc1de72dbcd79767f6da1 job1/estimators_0000.out +bfe2e3027aa471cb6775e5d00bc1bd03 job1/estimators_0001.out diff --git a/tests/classicmode_3d_inputfiles/results_md5_job0.txt b/tests/classicmode_3d_inputfiles/results_md5_job0.txt index 20077b81c..cb25b6f06 100644 --- a/tests/classicmode_3d_inputfiles/results_md5_job0.txt +++ b/tests/classicmode_3d_inputfiles/results_md5_job0.txt @@ -13,5 +13,5 @@ fb8c9b73ab7c5bb1fac556612d74812d linestat.out 663b0393f95a92f6a0f0fd931957f5c9 packets00_0001.out 8aad0eabdebf722460b05185de9a1a8c spec.out bc01f046eab9bf3802149e0ff63d5069 timesteps.out -1a934a7e3ebea1474f6bf2c90fa8bb1e job0/estimators_0000.out -c30f3fc1c4e98443f219a528dfecc6c8 job0/estimators_0001.out +5015246ae58957425f1276a6796adb08 job0/estimators_0000.out +f755dd3d0ac134553559cfd12fe2f3aa job0/estimators_0001.out diff --git a/tests/kilonova_1d_1dgrid_inputfiles/results_md5_final.txt b/tests/kilonova_1d_1dgrid_inputfiles/results_md5_final.txt index f8bcbf0f0..cd5ceaa23 100644 --- a/tests/kilonova_1d_1dgrid_inputfiles/results_md5_final.txt +++ b/tests/kilonova_1d_1dgrid_inputfiles/results_md5_final.txt @@ -14,5 +14,5 @@ f5b097e63480ea2b54bb88459f72b19a modelgridrankassignments.out 0b33390bb2afa6eb72e3eb66f8c6173c packets00_0001.out 5fbbcd0c30de4611d0f6e8e3b7644174 spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -6a6ebdc0648b2b21aef452bf8175f133 job1/estimators_0000.out -bc4574bd849763bb85d1c3811f193e5e job1/estimators_0001.out +2c3b065fff996381b2e0a45c70738f1e job1/estimators_0000.out +0639426a93fe8a9349b24676e19e0da5 job1/estimators_0001.out diff --git a/tests/kilonova_1d_1dgrid_inputfiles/results_md5_job0.txt b/tests/kilonova_1d_1dgrid_inputfiles/results_md5_job0.txt index c0e23510d..af67dca2f 100644 --- a/tests/kilonova_1d_1dgrid_inputfiles/results_md5_job0.txt +++ b/tests/kilonova_1d_1dgrid_inputfiles/results_md5_job0.txt @@ -13,5 +13,5 @@ d9dcb4378a7ed332f79a5071bae91a9e packets00_0000.out b0d755d49c026bf4e67bba200191fd35 packets00_0001.out 8269008eb999358088439c0330a996fd spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -d791e9b02393ab142c5a3e616dd1b668 job0/estimators_0000.out -944041c0ef95bf74c720c433505efe86 job0/estimators_0001.out +13aed5caf0552d2edc24823c650fbb51 job0/estimators_0000.out +7db264f657940fe58fbaeddcd88d9b62 job0/estimators_0001.out diff --git a/tests/kilonova_1d_3dgrid_inputfiles/results_md5_final.txt b/tests/kilonova_1d_3dgrid_inputfiles/results_md5_final.txt index f8d0fb41b..42b4cff39 100644 --- a/tests/kilonova_1d_3dgrid_inputfiles/results_md5_final.txt +++ b/tests/kilonova_1d_3dgrid_inputfiles/results_md5_final.txt @@ -14,5 +14,5 @@ def7d4e67866038f5a075bc74d1d225b packets00_0000.out a3bfac14bde7b6e25811144bfeacec4f packets00_0001.out 821414204b3fb0fab174b9895cbc6e68 spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -5d2348e48bb2d1ab37357fe5a90e0fcc job1/estimators_0000.out -dcf178e682d15f7bf16633e49bbf9bea job1/estimators_0001.out +3c8106ae605fffd169219d901d7ab212 job1/estimators_0000.out +c7b390ad5ca2566dd0b4d6477f50a3e6 job1/estimators_0001.out diff --git a/tests/kilonova_1d_3dgrid_inputfiles/results_md5_job0.txt b/tests/kilonova_1d_3dgrid_inputfiles/results_md5_job0.txt index 4acf72fc3..3af2dec78 100644 --- a/tests/kilonova_1d_3dgrid_inputfiles/results_md5_job0.txt +++ b/tests/kilonova_1d_3dgrid_inputfiles/results_md5_job0.txt @@ -13,5 +13,5 @@ a84d3d4a9c7d1f61394f1b903d5bb244 packets00_0000.out 627178c8fb51ac34030cd547cc9d268d packets00_0001.out 6fe9dc86f38f3b31814186c336019313 spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -d77d94d656fc59fd894cefffc97aca2b job0/estimators_0000.out -6b0031f1c76b0152f43013f4defc7e10 job0/estimators_0001.out +103e8461578e4790aefd83f90befebba job0/estimators_0000.out +9bd6654f641a8bd1911caf09f57f9d4c job0/estimators_0001.out diff --git a/tests/kilonova_2d_2dgrid_inputfiles/results_md5_final.txt b/tests/kilonova_2d_2dgrid_inputfiles/results_md5_final.txt index 0e9cef292..ac62fe61b 100644 --- a/tests/kilonova_2d_2dgrid_inputfiles/results_md5_final.txt +++ b/tests/kilonova_2d_2dgrid_inputfiles/results_md5_final.txt @@ -316,5 +316,5 @@ a1c91f6d89797e1feaf21d4d6186545d packets00_0000.out 6de6558232711b33cd042ab495c496ca spec.out 0dd442806e68ef09fb6b1d1f6b4a98ff spec_res.out a351f1711fecd60c023d0ba7332092db timesteps.out -ade51ef62b573258ac4974bbea28f39f job1/estimators_0000.out -cead8ecae6e4fbeeb9f49b1824176379 job1/estimators_0001.out +df59315c5b70e4891f2144e5eb183e55 job1/estimators_0000.out +395cc0f7be2bd78a2ce13b02c3d3ff4a job1/estimators_0001.out diff --git a/tests/kilonova_2d_2dgrid_inputfiles/results_md5_job0.txt b/tests/kilonova_2d_2dgrid_inputfiles/results_md5_job0.txt index 456097991..80d0e04b6 100644 --- a/tests/kilonova_2d_2dgrid_inputfiles/results_md5_job0.txt +++ b/tests/kilonova_2d_2dgrid_inputfiles/results_md5_job0.txt @@ -10,5 +10,5 @@ fa80b4c76909e5984e9f17d0548a9a73 modelgridrankassignments.out bd4a751fc9d0dc1c7fe7f60865b41a71 packets00_0001.out 600fc00772d3c1ad16948b6332791a9a spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -5ce4fd7564baff28fdc53625e53a83b1 job0/estimators_0000.out -85a8cb83f6ae8bb6594abdfa7ccd060c job0/estimators_0001.out +3fcf2b82ba502978e91e69837d81ab4f job0/estimators_0000.out +39a7dbf36a436eaaf3ba4ddc36987b92 job0/estimators_0001.out diff --git a/tests/kilonova_2d_3dgrid_inputfiles/results_md5_final.txt b/tests/kilonova_2d_3dgrid_inputfiles/results_md5_final.txt index dcc5c4250..3a986fb73 100644 --- a/tests/kilonova_2d_3dgrid_inputfiles/results_md5_final.txt +++ b/tests/kilonova_2d_3dgrid_inputfiles/results_md5_final.txt @@ -316,5 +316,5 @@ fa80b4c76909e5984e9f17d0548a9a73 modelgridrankassignments.out f2bea0bcfa94e8edcc00ce440857aa26 spec.out bd96aa33323c64e02bb6b901f304bfbb spec_res.out a351f1711fecd60c023d0ba7332092db timesteps.out -ce0c00c2d37a2829bbc04f1860622864 job1/estimators_0000.out -a25884811cb9bcd68c7a342852de882c job1/estimators_0001.out +7efbbab22bfa60568079d4b935fed3dd job1/estimators_0000.out +f8a807c360c450e26069965bd41ed4ee job1/estimators_0001.out diff --git a/tests/kilonova_2d_3dgrid_inputfiles/results_md5_job0.txt b/tests/kilonova_2d_3dgrid_inputfiles/results_md5_job0.txt index be6e9d755..45724882f 100644 --- a/tests/kilonova_2d_3dgrid_inputfiles/results_md5_job0.txt +++ b/tests/kilonova_2d_3dgrid_inputfiles/results_md5_job0.txt @@ -10,5 +10,5 @@ b435eb715747a70d3372b1cc98250065 packets00_0000.out cfb5460a9316c79c660bf7d913ea8e08 packets00_0001.out f8cfd141dfa301309bd0fba60c6f0c71 spec.out a351f1711fecd60c023d0ba7332092db timesteps.out -b66700748a8a5585562388ebca38f851 job0/estimators_0000.out -566cb793642c4e4ac43b50de2d114193 job0/estimators_0001.out +3f72e6eb4b1e5657f2681eceb7469337 job0/estimators_0000.out +a5a3c6eb38412ad0a274a14f80593191 job0/estimators_0001.out diff --git a/thermalbalance.cc b/thermalbalance.cc index 8b3f775ce..41518c48d 100644 --- a/thermalbalance.cc +++ b/thermalbalance.cc @@ -16,7 +16,6 @@ #include "radfield.h" #include "ratecoeff.h" #include "sn3d.h" -#include "update_grid.h" struct Te_solution_paras { double t_current; @@ -140,6 +139,7 @@ static auto get_bfheatingcoeff(int element, int ion, int level) -> double } void calculate_bfheatingcoeffs(int modelgridindex) { + const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); const double minelfrac = 0.01; for (int element = 0; element < get_nelements(); element++) { if (grid::get_elem_abundance(modelgridindex, element) <= minelfrac && !USE_LUT_BFHEATING) { @@ -171,8 +171,8 @@ void calculate_bfheatingcoeffs(int modelgridindex) { const int index_in_groundlevelcontestimator = globals::elements[element].ions[ion].levels[level].closestgroundlevelcont; if (index_in_groundlevelcontestimator >= 0) { - bfheatingcoeff *= globals::bfheatingestimator[modelgridindex * get_nelements() * get_max_nions() + - index_in_groundlevelcontestimator]; + bfheatingcoeff *= + globals::bfheatingestimator[nonemptymgi * get_includedions() + index_in_groundlevelcontestimator]; } } } @@ -282,11 +282,6 @@ static void calculate_heating_rates(const int modelgridindex, const double T_e, // */ // } - /// Bound-free heating (from estimators) - /// ------------------------------------ - // if (ion < nions-1) bfheating += - // globals::bfheatingestimator[cellnumber*get_nelements()*get_max_nions()+element*get_max_nions()+ion]; - /// Bound-free heating (renormalised analytical calculation) /// -------------------------------------------------------- /// We allow bound free-transitions only if there is a higher ionisation stage diff --git a/thermalbalance.h b/thermalbalance.h index 598a6c341..c8dbd95b2 100644 --- a/thermalbalance.h +++ b/thermalbalance.h @@ -15,7 +15,7 @@ struct heatingcoolingrates { void call_T_e_finder(int modelgridindex, int timestep, double t_current, double T_min, double T_max, struct heatingcoolingrates *heatingcoolingrates); -double get_bfheatingcoeff_ana(int element, int ion, int level, int phixstargetindex, double T, double W); +auto get_bfheatingcoeff_ana(int element, int ion, int level, int phixstargetindex, double T, double W) -> double; void calculate_bfheatingcoeffs(int modelgridindex); #endif // THERMALBALANCE_H diff --git a/update_grid.cc b/update_grid.cc index 0712cacbb..945107d84 100644 --- a/update_grid.cc +++ b/update_grid.cc @@ -2,9 +2,11 @@ #include +#include "artisoptions.h" #include "atomic.h" #include "constants.h" #include "decay.h" +#include "globals.h" #include "grid.h" #include "kpkt.h" #include "ltepop.h" @@ -625,13 +627,13 @@ static void write_to_estimators_file(FILE *estimators_file, const int mgi, const if (USE_LUT_PHOTOION && globals::nbfcontinua > 0) { fprintf(estimators_file, "corrphotoionrenorm Z=%2d", get_atomicnumber(element)); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < nions - 1; ion++) { fprintf(estimators_file, " %d: %9.3e", get_ionstage(element, ion), globals::corrphotoionrenorm[get_ionestimindex(mgi, element, ion)]); } fprintf(estimators_file, "\n"); fprintf(estimators_file, "gammaestimator Z=%2d", get_atomicnumber(element)); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < nions - 1; ion++) { fprintf(estimators_file, " %d: %9.3e", get_ionstage(element, ion), globals::gammaestimator[get_ionestimindex(mgi, element, ion)]); } @@ -867,13 +869,13 @@ static void solve_Te_nltepops(const int n, const int nts, const int titer, } } -static void update_gamma_corrphotoionrenorm_bfheating_estimators(const int n, const double estimator_normfactor) { +static void update_gamma_corrphotoionrenorm_bfheating_estimators(const int mgi, const double estimator_normfactor) { assert_always(USE_LUT_PHOTOION || USE_LUT_BFHEATING); if constexpr (USE_LUT_PHOTOION) { for (int element = 0; element < get_nelements(); element++) { const int nions = get_nions(element); for (int ion = 0; ion < nions - 1; ion++) { - const int ionestimindex = get_ionestimindex(n, element, ion); + const int ionestimindex = get_ionestimindex(mgi, element, ion); globals::gammaestimator[ionestimindex] *= estimator_normfactor / H; #ifdef DO_TITER if (globals::gammaestimator_save[ionestimindex] >= 0) { @@ -884,14 +886,14 @@ static void update_gamma_corrphotoionrenorm_bfheating_estimators(const int n, co #endif globals::corrphotoionrenorm[ionestimindex] = - globals::gammaestimator[ionestimindex] / get_corrphotoioncoeff_ana(element, ion, 0, 0, n); + globals::gammaestimator[ionestimindex] / get_corrphotoioncoeff_ana(element, ion, 0, 0, mgi); if (!std::isfinite(globals::corrphotoionrenorm[ionestimindex])) { printout( "[fatal] about to set corrphotoionrenorm = NaN = gammaestimator / " "get_corrphotoioncoeff_ana(%d,%d,%d,%d,%d)=%g/%g", - element, ion, 0, 0, n, globals::gammaestimator[ionestimindex], - get_corrphotoioncoeff_ana(element, ion, 0, 0, n)); + element, ion, 0, 0, mgi, globals::gammaestimator[ionestimindex], + get_corrphotoioncoeff_ana(element, ion, 0, 0, mgi)); abort(); } } @@ -910,10 +912,10 @@ static void update_gamma_corrphotoionrenorm_bfheating_estimators(const int n, co /// Reuse the gammaestimator array as temporary storage of the Gamma values during /// the remaining part of the update_grid phase. Afterwards it is reset to record /// the next timesteps gamma estimators. - const int ionestimindex = get_ionestimindex(n, element, ion); + const int ionestimindex = get_ionestimindex(mgi, element, ion); if constexpr (USE_LUT_PHOTOION) { - globals::gammaestimator[ionestimindex] = calculate_iongamma_per_gspop(n, element, ion); + globals::gammaestimator[ionestimindex] = calculate_iongamma_per_gspop(mgi, element, ion); } if constexpr (USE_LUT_BFHEATING) { @@ -929,14 +931,15 @@ static void update_gamma_corrphotoionrenorm_bfheating_estimators(const int n, co /// get_bfheating in the remaining part of update_grid. Later on it's reset and new /// contributions are added up. - const double bfheatingcoeff_ana = get_bfheatingcoeff_ana(element, ion, 0, 0, grid::get_TR(n), grid::get_W(n)); + const double bfheatingcoeff_ana = + get_bfheatingcoeff_ana(element, ion, 0, 0, grid::get_TR(mgi), grid::get_W(mgi)); globals::bfheatingestimator[ionestimindex] = globals::bfheatingestimator[ionestimindex] / bfheatingcoeff_ana; if (!std::isfinite(globals::bfheatingestimator[ionestimindex])) { printout( "[fatal] about to set bfheatingestimator = NaN = bfheatingestimator / " "get_bfheatingcoeff_ana(%d,%d,%d,%d,%d)=%g/%g", - element, ion, 0, 0, n, globals::bfheatingestimator[ionestimindex], bfheatingcoeff_ana); + element, ion, 0, 0, mgi, globals::bfheatingestimator[ionestimindex], bfheatingcoeff_ana); abort(); } } @@ -961,8 +964,7 @@ static void titer_average_estimators(const int n) { static void zero_gammaestimator(const int modelgridindex) { assert_always(USE_LUT_PHOTOION); for (int element = 0; element < get_nelements(); element++) { - const int nions = get_nions(element); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < (get_nions(element) - 1); ion++) { globals::gammaestimator[get_ionestimindex(modelgridindex, element, ion)] = 0.; } } @@ -971,8 +973,7 @@ static void zero_gammaestimator(const int modelgridindex) { static void set_all_corrphotoionrenorm(const int modelgridindex, const double value) { assert_always(USE_LUT_PHOTOION); for (int element = 0; element < get_nelements(); element++) { - const int nions = get_nions(element); - for (int ion = 0; ion < nions; ion++) { + for (int ion = 0; ion < (get_nions(element) - 1); ion++) { globals::corrphotoionrenorm[get_ionestimindex(modelgridindex, element, ion)] = value; } } @@ -1184,11 +1185,6 @@ static void update_grid_cell(const int mgi, const int nts, const int nts_prev, c grid::set_TJ(mgi, 0.); grid::set_Te(mgi, 0.); grid::set_W(mgi, 0.); - - if constexpr (USE_LUT_PHOTOION) { - zero_gammaestimator(mgi); - set_all_corrphotoionrenorm(mgi, 0.); - } } } @@ -1208,9 +1204,7 @@ void update_grid(FILE *estimators_file, const int nts, const int nts_prev, const /// unless they have been read from file if ((!globals::simulation_continued_from_saved) || (nts - globals::timestep_initial != 0) || (titer != 0)) { printout("nts %d, titer %d: reset corr photoionrenorm\n", nts, titer); - for (int i = 0; i < grid::get_npts_model() * get_nelements() * get_max_nions(); i++) { - globals::corrphotoionrenorm[i] = 0.; - } + std::fill_n(globals::corrphotoionrenorm, grid::get_nonempty_npts_model() * get_includedions(), 0.); printout("after nts %d, titer %d: reset corr photoionrenorm\n", nts, titer); } } @@ -1274,7 +1268,7 @@ void update_grid(FILE *estimators_file, const int nts, const int nts_prev, const #endif { write_to_estimators_file(estimators_file, mgi, nts, titer, &heatingcoolingrates); } - } else { + } else if (grid::get_numassociatedcells(mgi) > 0) { /// else, only reset gammaestimator to zero. This allows us to do a global MPI /// communication after update_grid to synchronize gammaestimator /// and write a contiguous restart file with grid properties diff --git a/vectors.h b/vectors.h index 8b17aae12..3689104da 100644 --- a/vectors.h +++ b/vectors.h @@ -190,7 +190,7 @@ constexpr void move_pkt_withtime(struct packet *pkt_ptr, const double distance) pkt_ptr->nu_cmf = std::min(pkt_ptr->nu_cmf, nu_cmf_old); } -[[nodiscard]] [[gnu::pure]] constexpr double get_arrive_time(const struct packet *const pkt_ptr) +[[nodiscard]] [[gnu::pure]] constexpr auto get_arrive_time(const struct packet *const pkt_ptr) -> double /// We know that a packet escaped at "escape_time". However, we have /// to allow for travel time. Use the formula in Leon's paper. The extra /// distance to be travelled beyond the reference surface is ds = r_ref (1 - mu). @@ -198,7 +198,7 @@ constexpr void move_pkt_withtime(struct packet *pkt_ptr, const double distance) return pkt_ptr->escape_time - (dot(pkt_ptr->pos, pkt_ptr->dir) / CLIGHT_PROP); } -inline double get_arrive_time_cmf(const struct packet *pkt_ptr) { +inline auto get_arrive_time_cmf(const struct packet *pkt_ptr) -> double { return pkt_ptr->escape_time * std::sqrt(1. - (globals::vmax * globals::vmax / CLIGHTSQUARED)); } diff --git a/vpkt.h b/vpkt.h index dd385d828..95db78daa 100644 --- a/vpkt.h +++ b/vpkt.h @@ -14,7 +14,7 @@ void frame_transform(std::span n_rf, double *Q, double *U, std: void read_parameterfile_vpkt(); void vpkt_init(int nts, int my_rank, int tid, bool continued_from_saved); -void vpkt_call_estimators(struct packet *pkt_ptr, const enum packet_type); +void vpkt_call_estimators(struct packet *pkt_ptr, enum packet_type); void vpkt_write_timestep(int nts, int my_rank, int tid, bool is_final); void vpkt_remove_temp_file(int nts, int my_rank); From dc5b22f97fc183fed48f00c2fe481563651d670d Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Tue, 28 Nov 2023 05:01:26 +0000 Subject: [PATCH 3/9] OpenMP: use nonmonotonic schedule for packet updates --- radfield.cc | 8 +++----- radfield.h | 2 +- rpkt.cc | 4 ++-- update_packets.cc | 2 +- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/radfield.cc b/radfield.cc index 520246688..dfb5749bb 100644 --- a/radfield.cc +++ b/radfield.cc @@ -323,7 +323,7 @@ void init(int my_rank, int ndo_nonempty) if (globals::rank_in_node == 0) { my_rank_cells += nonempty_npts_model - (my_rank_cells * globals::node_nprocs); } - MPI_Aint size = my_rank_cells * RADFIELDBINCOUNT * sizeof(struct radfieldbin_solution); + auto size = static_cast(my_rank_cells * RADFIELDBINCOUNT * sizeof(struct radfieldbin_solution)); int disp_unit = sizeof(struct radfieldbin_solution); MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, &radfieldbin_solutions, &win_radfieldbin_solutions); @@ -354,7 +354,7 @@ void init(int my_rank, int ndo_nonempty) if (globals::rank_in_node == 0) { my_rank_cells += nonempty_npts_model - (my_rank_cells * globals::node_nprocs); } - MPI_Aint size = my_rank_cells * globals::nbfcontinua * sizeof(float); + auto size = static_cast(my_rank_cells * globals::nbfcontinua * sizeof(float)); int disp_unit = sizeof(float); MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, &prev_bfrate_normed, &win_prev_bfrate_normed); @@ -718,10 +718,8 @@ static void update_bfestimators(const int nonemptymgi, const double distance_e_c } } -void update_estimators(const int modelgridindex, const double distance_e_cmf, const double nu_cmf, +void update_estimators(const int nonemptymgi, const double distance_e_cmf, const double nu_cmf, const struct packet *const pkt_ptr) { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); - safeadd(J[nonemptymgi], distance_e_cmf); safeadd(nuJ[nonemptymgi], distance_e_cmf * nu_cmf); diff --git a/radfield.h b/radfield.h index 21bd1e879..052dae006 100644 --- a/radfield.h +++ b/radfield.h @@ -14,7 +14,7 @@ void init(int my_rank, int ndo_nonempty); void initialise_prev_titer_photoionestimators(); void write_to_file(int modelgridindex, int timestep); void close_file(); -void update_estimators(int modelgridindex, double distance_e_cmf, double nu_cmf, const struct packet *pkt_ptr); +void update_estimators(int nonemptymgi, double distance_e_cmf, double nu_cmf, const struct packet *pkt_ptr); void update_lineestimator(int modelgridindex, int lineindex, double increment); [[nodiscard]] auto radfield(double nu, int modelgridindex) -> double; void fit_parameters(int modelgridindex, int timestep); diff --git a/rpkt.cc b/rpkt.cc index 3d06fe5f6..1135deab5 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -545,7 +545,8 @@ static void update_estimators(const struct packet *pkt_ptr, const double distanc } const double distance_e_cmf = distance * pkt_ptr->e_cmf; const double nu = pkt_ptr->nu_cmf; - radfield::update_estimators(modelgridindex, distance_e_cmf, nu, pkt_ptr); + const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + radfield::update_estimators(nonemptymgi, distance_e_cmf, nu, pkt_ptr); /// ffheatingestimator does not depend on ion and element, so an array with gridsize is enough. /// quick and dirty solution: store info in element=ion=0, and leave the others untouched (i.e. zero) @@ -553,7 +554,6 @@ static void update_estimators(const struct packet *pkt_ptr, const double distanc if constexpr (USE_LUT_PHOTOION || USE_LUT_BFHEATING) { const double distance_e_cmf_over_nu = distance_e_cmf / nu; - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); for (int i = 0; i < globals::nbfcontinua_ground; i++) { const double nu_edge = globals::groundcont[i].nu_edge; diff --git a/update_packets.cc b/update_packets.cc index 27ea21191..e13ee7ca4 100644 --- a/update_packets.cc +++ b/update_packets.cc @@ -274,7 +274,7 @@ void update_packets(const int my_rank, const int nts, struct packet *packets) const int updatecellcounter_beforepass = stats::get_counter(stats::COUNTER_UPDATECELL); #ifdef _OPENMP -#pragma omp parallel for schedule(dynamic) +#pragma omp parallel for schedule(nonmonotonic : dynamic) #endif for (int n = 0; n < globals::npkts; n++) { struct packet *pkt_ptr = &packets[n]; From f00180438eab5559b8f52771cda6780fe6f52602 Mon Sep 17 00:00:00 2001 From: Fionntan Callan <56071719+fionntancallan@users.noreply.github.com> Date: Fri, 1 Dec 2023 16:59:36 +0000 Subject: [PATCH 4/9] Set vpkts negative line optical depths to zero (#48) -Set vpkts line optical depths to zero if they are negative. Line optical depths were never negative in ARTIS-CLASSIC as the level populations are coming from Boltzmann statistics but population inversions can occur when running in ARTIS-NLTE resulting in negative line optical depths. There was previously no check for this for vpkts. -Set line optical depths to zero in unused function get_tau_sobolev for consistency -Add slurm submit script for cosma8 HPC --- atomic.cc | 2 +- scripts/artis-cosma8.sh | 50 +++++++++++++++++++++++++++++++++++++++++ vpkt.cc | 4 ++-- 3 files changed, 53 insertions(+), 3 deletions(-) create mode 100755 scripts/artis-cosma8.sh diff --git a/atomic.cc b/atomic.cc index 219fd0608..d08dd9378 100644 --- a/atomic.cc +++ b/atomic.cc @@ -60,7 +60,7 @@ auto get_tau_sobolev(const int modelgridindex, const int lineindex, const double const double B_ul = CLIGHTSQUAREDOVERTWOH / pow(nu_trans, 3) * A_ul; const double B_lu = stat_weight(element, ion, upper) / stat_weight(element, ion, lower) * B_ul; - const double tau_sobolev = (B_lu * n_l - B_ul * n_u) * HCLIGHTOVERFOURPI * t_current; + const double tau_sobolev = std::max(0., (B_lu * n_l - B_ul * n_u) * HCLIGHTOVERFOURPI * t_current); return tau_sobolev; } diff --git a/scripts/artis-cosma8.sh b/scripts/artis-cosma8.sh new file mode 100755 index 000000000..2eb2391d8 --- /dev/null +++ b/scripts/artis-cosma8.sh @@ -0,0 +1,50 @@ +#!/bin/bash -l + +#SBATCH --ntasks 3072 +#SBATCH -J artis-cosma8.sh +## SBATCH -o standard_output_file.%J.out +## SBATCH -e standard_error_file.%J.err +#SBATCH -p cosma8 +#SBATCH -A dp033 +## SBATCH --exclusive +#SBATCH -t 70:00:00 +#SBATCH --mail-type=ALL # notifications for job done & fail +## SBATCH --mail-user=f.callan@qub.ac.uk + +module purge +#load the modules used to build your program. +module load cosma +module load gsl/2.4 +module load gnu_comp/13.1.0 +module load openmpi/4.1.4 +module load python/3.10.12 + +module list + +cd $SLURM_SUBMIT_DIR + +echo "CPU type: $(c++ -march=native -Q --help=target | grep -- '-march= ' | cut -f3)" + +hoursleft=$(python3 ./artis/scripts/slurmjobhoursleft.py ${SLURM_JOB_ID}) +echo "$(date): before srun sn3d. hours left: $hoursleft" +# time srun -- ./sn3d -w $hoursleft > out.txt +time mpirun ./sn3d -w $hoursleft > out.txt +hoursleftafter=$(python3 ./artis/scripts/slurmjobhoursleft.py ${SLURM_JOB_ID}) +echo "$(date): after srun sn3d finished. hours left: $hoursleftafter" +hourselapsed=$(python3 -c "print($hoursleft - $hoursleftafter)") +echo "hours of runtime: $hourselapsed" +cpuhrs=$(python3 -c "print($SLURM_NTASKS * $hourselapsed)") +echo "ntasks: $SLURM_NTASKS -> CPU core hrs: $cpuhrs" + +mkdir ${SLURM_JOB_ID}.slurm +./artis/scripts/movefiles.sh ${SLURM_JOB_ID}.slurm + +if grep -q "RESTART_NEEDED" "output_0-0.txt" +then + sbatch ./artis-cosma8.sh + # sbatch $SLURM_JOB_NAME +fi + +if [ -f packets00_0000.out ]; then + sbatch ./artis/scripts/exspec-gzip-cosma8.sh +fi diff --git a/vpkt.cc b/vpkt.cc index 27d98c803..0f5bb28fb 100644 --- a/vpkt.cc +++ b/vpkt.cc @@ -346,7 +346,7 @@ static void rlc_emiss_vpkt(const struct packet *const pkt_ptr, const double t_cu const auto n_u = calculate_levelpop(mgi, element, ion, upper); const auto n_l = calculate_levelpop(mgi, element, ion, lower); - const double tau_line = (B_lu * n_l - B_ul * n_u) * HCLIGHTOVERFOURPI * t_line; + const double tau_line = std::max(0., (B_lu * n_l - B_ul * n_u) * HCLIGHTOVERFOURPI * t_line); // Check on the element to exclude // NB: ldist before need to be computed anyway (I want to move the packets to the @@ -1120,4 +1120,4 @@ void frame_transform(std::span n_rf, double *Q, double *U, std: // Compute Stokes Parameters in the CMF *Q = cos(2 * theta_rot) * p; *U = sin(2 * theta_rot) * p; -} \ No newline at end of file +} From 8049b378c99018b23fbea2716a04bc9b4224498b Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Tue, 12 Dec 2023 12:54:24 +0000 Subject: [PATCH 5/9] Fix zseed input type (#49) --- input.cc | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/input.cc b/input.cc index c39aee20b..f4f19ac15 100644 --- a/input.cc +++ b/input.cc @@ -1731,33 +1731,31 @@ auto get_noncommentline(std::fstream &input, std::string &line) -> bool void read_parameterfile(int rank) /// Subroutine to read in input parameters from input.txt. { - uint_least64_t pre_zseed = 0; - auto file = fstream_required("input.txt", std::ios::in); std::string line; assert_always(get_noncommentline(file, line)); - uint_fast64_t zseed_input = 0; + int64_t zseed_input = -1; std::istringstream(line) >> zseed_input; - if (zseed_input > 0) { - pre_zseed = zseed_input; // random number seed - printout("using input.txt specified random number seed of %lu\n", pre_zseed); - } else { - pre_zseed = std::random_device{}(); - printout("randomly-generated random number seed is %lu\n", pre_zseed); - } - #ifdef _OPENMP #pragma omp parallel { #endif + uint64_t pre_zseed{}; + if (zseed_input > 0) { + pre_zseed = static_cast(zseed_input); // random number seed + printout("using input.txt specified random number seed of %lu\n", pre_zseed); + } else { + pre_zseed = std::random_device{}(); + printout("randomly-generated random number seed is %lu\n", pre_zseed); + } /// For MPI parallelisation, the random seed is changed based on the rank of the process /// For OpenMP parallelisation rng is a threadprivate variable and the seed changed according /// to the thread-ID tid. - const auto zseed = pre_zseed + 13 * (rank * get_num_threads() + tid); + const uint64_t zseed = pre_zseed + static_cast(13 * (rank * get_num_threads() + tid)); printout("rank %d: thread %d has zseed %lu\n", rank, tid, zseed); rng_init(zseed); /// call it a few times From dc427c990bdc04747daf24c963c2b840c90ddcbc Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Tue, 6 Feb 2024 16:19:37 +0000 Subject: [PATCH 6/9] Fix potential segfault --- .vscode/settings.json | 2 +- rpkt.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 29a66a967..b8bda0698 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -2,7 +2,7 @@ "[cpp]": { "editor.tabSize": 2, "editor.codeActionsOnSave": { - "source.fixAll": true + "source.fixAll": "explicit" }, "editor.formatOnSave": true }, diff --git a/rpkt.cc b/rpkt.cc index 1135deab5..41411caf3 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -44,7 +44,7 @@ auto closest_transition(const double nu_cmf, const int next_trans) -> int // will find the highest frequency (lowest index) line with nu_line <= nu_cmf // lower_bound matches the first element where the comparison function is false const auto *matchline = - std::lower_bound(&globals::linelist[next_trans], &globals::linelist[globals::nlines], nu_cmf, + std::lower_bound(globals::linelist, globals::linelist + globals::nlines, nu_cmf, [](const auto &line, const double nu_cmf) -> bool { return line.nu > nu_cmf; }); const int matchindex = std::distance(globals::linelist, matchline); if (matchindex >= globals::nlines) { From e941118ad01d9ba0c21e0afd29d7267fa37a321b Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Tue, 6 Feb 2024 16:34:59 +0000 Subject: [PATCH 7/9] Fix cppcheck warning by passing const reference --- rpkt.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rpkt.cc b/rpkt.cc index 41411caf3..6292aab64 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -363,7 +363,8 @@ static void electron_scatter_rpkt(struct packet *pkt_ptr) { } static void rpkt_event_continuum(struct packet *pkt_ptr, - struct rpkt_continuum_absorptioncoeffs chi_rpkt_cont_thisthread, int modelgridindex) { + const struct rpkt_continuum_absorptioncoeffs &chi_rpkt_cont_thisthread, + int modelgridindex) { const double nu = pkt_ptr->nu_cmf; const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt_ptr->pos, pkt_ptr->dir, pkt_ptr->prop_time); From 330867f327e5f5dcf8c84861917188e4f266996b Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Wed, 13 Mar 2024 13:41:38 +0000 Subject: [PATCH 8/9] JUWELS: Update job script with --hint=nomultithread for slurm 22.05 change --- scripts/artis-juwels.sh | 3 ++- scripts/exspec-gzip-juwels.sh | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/scripts/artis-juwels.sh b/scripts/artis-juwels.sh index 08d8b05c8..6b945551c 100755 --- a/scripts/artis-juwels.sh +++ b/scripts/artis-juwels.sh @@ -11,6 +11,7 @@ ##SBATCH --mail-user=luke.shingles@gmail.com module load Stages/2024 GCC ParaStationMPI GSL +module load UCX-settings/plain module list @@ -20,7 +21,7 @@ echo "CPU type: $(c++ -march=native -Q --help=target | grep -- '-march= ' | cut hoursleft=$(python3 ./artis/scripts/slurmjobhoursleft.py ${SLURM_JOB_ID}) echo "$(date): before srun sn3d. hours left: $hoursleft" -time srun -- ./sn3d -w $hoursleft > out.txt +time srun --hint=nomultithread -- ./sn3d -w $hoursleft > out.txt hoursleftafter=$(python3 ./artis/scripts/slurmjobhoursleft.py ${SLURM_JOB_ID}) echo "$(date): after srun sn3d finished. hours left: $hoursleftafter" hourselapsed=$(python3 -c "print($hoursleft - $hoursleftafter)") diff --git a/scripts/exspec-gzip-juwels.sh b/scripts/exspec-gzip-juwels.sh index 070bd74c4..e294fe84a 100755 --- a/scripts/exspec-gzip-juwels.sh +++ b/scripts/exspec-gzip-juwels.sh @@ -9,7 +9,7 @@ #SBATCH --mail-type=ALL ##SBATCH --mail-user=luke.shingles@gmail.com -module load Stages/2023 GCC ParaStationMPI GSL zstd/.1.5.2 +module load Stages/2024 ParaStationMPI GSL zstd/.1.5.5 cd $SLURM_SUBMIT_DIR From d1eec12fdb046b0dbcd37bde93af9b1b0a88e50b Mon Sep 17 00:00:00 2001 From: jpollin98 <152508849+jpollin98@users.noreply.github.com> Date: Mon, 25 Mar 2024 16:18:09 +0000 Subject: [PATCH 9/9] Add option to limit number of BF estimators (#47) Co-authored-by: Luke Shingles --- .github/workflows/ci-checks.yml | 2 +- .github/workflows/ci.yml | 1 + .idea/.gitignore | 8 ++ .idea/codeStyles/Project.xml | 7 ++ .idea/codeStyles/codeStyleConfig.xml | 5 + .idea/misc.xml | 18 +++ .idea/vcs.xml | 6 + .pre-commit-config.yaml | 2 +- artisoptions_christinenonthermal.h | 2 + artisoptions_classic.h | 2 + artisoptions_kilonova_lte.h | 2 + artisoptions_nltenebular.h | 8 +- artisoptions_nltewithoutnonthermal.h | 2 + atomic.cc | 8 +- decay.cc | 16 +-- globals.cc | 2 + globals.h | 3 + grid.cc | 4 +- input.cc | 15 ++- ltepop.cc | 8 +- ltepop.h | 8 +- macroatom.cc | 4 +- macroatom.h | 4 +- nltepop.cc | 10 +- nonthermal.cc | 12 +- radfield.cc | 115 ++++++++++-------- ratecoeff.cc | 4 +- rpkt.h | 4 +- stats.cc | 4 +- .../results_md5_final.txt | 20 +++ .../results_md5_job0.txt | 19 +++ ...tup_nebularonezone_1d_3dgrid_limitbfest.sh | 45 +++++++ thermalbalance.cc | 4 +- vpkt.cc | 4 +- 34 files changed, 274 insertions(+), 104 deletions(-) create mode 100644 .idea/.gitignore create mode 100644 .idea/codeStyles/Project.xml create mode 100644 .idea/codeStyles/codeStyleConfig.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/vcs.xml create mode 100644 tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_final.txt create mode 100644 tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_job0.txt create mode 100755 tests/setup_nebularonezone_1d_3dgrid_limitbfest.sh diff --git a/.github/workflows/ci-checks.yml b/.github/workflows/ci-checks.yml index 37bb5f256..c7058e23c 100644 --- a/.github/workflows/ci-checks.yml +++ b/.github/workflows/ci-checks.yml @@ -40,7 +40,7 @@ jobs: - name: Run clang-format style check uses: jidicula/clang-format-action@v4.11.0 with: - clang-format-version: '16' + clang-format-version: '18' check-path: . clang-tidy: diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 47cec82bd..d2fb353ff 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -31,6 +31,7 @@ jobs: kilonova_2d_2dgrid, kilonova_2d_3dgrid, nebularonezone_1d_3dgrid, + nebularonezone_1d_3dgrid_limitbfest, ] exclude: - os: self-hosted diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 000000000..13566b81b --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,8 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/.idea/codeStyles/Project.xml b/.idea/codeStyles/Project.xml new file mode 100644 index 000000000..f60388162 --- /dev/null +++ b/.idea/codeStyles/Project.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/codeStyles/codeStyleConfig.xml b/.idea/codeStyles/codeStyleConfig.xml new file mode 100644 index 000000000..79ee123c2 --- /dev/null +++ b/.idea/codeStyles/codeStyleConfig.xml @@ -0,0 +1,5 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 000000000..53624c9e1 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,18 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 000000000..35eb1ddfb --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 99fe3f88b..afd2ac7f4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -27,7 +27,7 @@ repos: rev: v1.3.5 hooks: - id: clang-format - args: [-i] + args: [--version=18, -i] # - id: clang-tidy # args: [-extra-arg=-std=c++20] # - id: oclint diff --git a/artisoptions_christinenonthermal.h b/artisoptions_christinenonthermal.h index 9c5cf6ead..53d0b6f92 100644 --- a/artisoptions_christinenonthermal.h +++ b/artisoptions_christinenonthermal.h @@ -81,6 +81,8 @@ constexpr bool DETAILED_BF_ESTIMATORS_ON = true; constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13; +constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) { return true; } + constexpr bool USE_LUT_PHOTOION = false; constexpr bool USE_LUT_BFHEATING = false; diff --git a/artisoptions_classic.h b/artisoptions_classic.h index 08eb7f280..4e395310f 100644 --- a/artisoptions_classic.h +++ b/artisoptions_classic.h @@ -75,6 +75,8 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false; constexpr bool DETAILED_BF_ESTIMATORS_ON = false; +constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) { return false; } + constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13; constexpr bool USE_LUT_PHOTOION = true; diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h index 02b8000fc..250824bad 100644 --- a/artisoptions_kilonova_lte.h +++ b/artisoptions_kilonova_lte.h @@ -74,6 +74,8 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false; constexpr bool DETAILED_BF_ESTIMATORS_ON = false; +constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) { return false; } + constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13; constexpr bool USE_LUT_PHOTOION = true; diff --git a/artisoptions_nltenebular.h b/artisoptions_nltenebular.h index 4de2507d9..ca2704692 100644 --- a/artisoptions_nltenebular.h +++ b/artisoptions_nltenebular.h @@ -82,6 +82,12 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false; constexpr bool DETAILED_BF_ESTIMATORS_ON = true; +constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) { + // To only BF estimators for NLTE levels: + // return LEVEL_IS_NLTE(element_z, ionstage, level); + return true; +} + constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13; constexpr bool USE_LUT_PHOTOION = false; @@ -146,4 +152,4 @@ constexpr bool KEEP_ALL_RESTART_FILES = false; constexpr bool BFCOOLING_USELEVELPOPNOTIONPOP = false; // NOLINTEND(modernize*,misc-unused-parameters) -#endif // ARTISOPTIONS_H \ No newline at end of file +#endif // ARTISOPTIONS_H diff --git a/artisoptions_nltewithoutnonthermal.h b/artisoptions_nltewithoutnonthermal.h index d42344dcb..f295413f4 100644 --- a/artisoptions_nltewithoutnonthermal.h +++ b/artisoptions_nltewithoutnonthermal.h @@ -78,6 +78,8 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false; constexpr bool DETAILED_BF_ESTIMATORS_ON = true; +constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) { return true; } + constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13; constexpr bool USE_LUT_PHOTOION = false; diff --git a/atomic.cc b/atomic.cc index d08dd9378..a6f01bebb 100644 --- a/atomic.cc +++ b/atomic.cc @@ -21,8 +21,8 @@ int includedions = 0; int includedions_excludehighest = 0; std::array phixs_file_version_exists; -auto get_continuumindex_phixstargetindex(const int element, const int ion, const int level, const int phixstargetindex) - -> int +auto get_continuumindex_phixstargetindex(const int element, const int ion, const int level, + const int phixstargetindex) -> int /// Returns the index of the continuum associated to the given level. { return globals::elements[element].ions[ion].levels[level].cont_index - phixstargetindex; @@ -90,8 +90,8 @@ auto level_isinsuperlevel(const int element, const int ion, const int level) -> return (!is_nlte(element, ion, level) && level != 0 && (get_nlevels_nlte(element, ion) > 0)); } -auto photoionization_crosssection_fromtable(const float *const photoion_xs, const double nu_edge, const double nu) - -> double +auto photoionization_crosssection_fromtable(const float *const photoion_xs, const double nu_edge, + const double nu) -> double /// Calculates the photoionisation cross-section at frequency nu out of the atomic data. /// Input: - edge frequency nu_edge of the desired bf-continuum /// - nu diff --git a/decay.cc b/decay.cc index 25adb348d..b2f4e3a46 100644 --- a/decay.cc +++ b/decay.cc @@ -751,8 +751,8 @@ static auto sample_decaytime(const int decaypathindex, const double tdecaymin, c } static constexpr auto calculate_decaychain(const double firstinitabund, const std::vector &lambdas, - const int num_nuclides, const double timediff, const bool useexpansionfactor) - -> double { + const int num_nuclides, const double timediff, + const bool useexpansionfactor) -> double { // calculate final number abundance from multiple decays, e.g., Ni56 -> Co56 -> Fe56 (nuc[0] -> nuc[1] -> nuc[2]) // the top nuclide initial abundance is set and the chain-end abundance is returned (all intermediates nuclides // are assumed to start with zero abundance) @@ -912,8 +912,8 @@ static auto get_endecay_to_tinf_per_ejectamass_at_time(const int modelgridindex, } auto get_endecay_per_ejectamass_t0_to_time_withexpansion_chain_numerical(const int modelgridindex, - const int decaypathindex, const double tstart) - -> double + const int decaypathindex, + const double tstart) -> double // just here as as check on the analytic result from get_endecay_per_ejectamass_t0_to_time_withexpansion() // this version does an Euler integration { @@ -979,8 +979,8 @@ auto get_endecay_per_ejectamass_t0_to_time_withexpansion(const int modelgridinde return tot_endecay; } -static auto get_endecay_per_ejectamass_between_times(const int mgi, const int decaypathindex, double tlow, double thigh) - -> double +static auto get_endecay_per_ejectamass_between_times(const int mgi, const int decaypathindex, double tlow, + double thigh) -> double // get decay energy per mass [erg/g] released by a decaypath between times tlow [s] and thigh [s] { assert_always(tlow <= thigh); @@ -1016,8 +1016,8 @@ static auto get_simtime_endecay_per_ejectamass(const int mgi, const int decaypat return chainendecay; } -static auto get_decaypath_power_per_ejectamass(const int decaypathindex, const int modelgridindex, const double time) - -> double +static auto get_decaypath_power_per_ejectamass(const int decaypathindex, const int modelgridindex, + const double time) -> double // total decay power per mass [erg/s/g] for a given decaypath { // only decays at the end of the chain contributed from the initial abundance of the top of the chain are counted diff --git a/globals.cc b/globals.cc index 369193f5e..82664375a 100644 --- a/globals.cc +++ b/globals.cc @@ -12,6 +12,8 @@ std::unique_ptr timesteps = nullptr; double *rpkt_emiss = nullptr; /// Volume estimator for the rpkt emissivity +int bfestimcount = 0; + // for USE_LUT_PHOTOION = true double *corrphotoionrenorm = nullptr; double *gammaestimator = nullptr; diff --git a/globals.h b/globals.h index 81d2914fd..0352f5b76 100644 --- a/globals.h +++ b/globals.h @@ -50,6 +50,7 @@ struct fullphixslist { float *photoion_xs; double probability; int index_in_groundphixslist; + int bfestimindex; }; struct groundphixslist { @@ -209,6 +210,8 @@ extern std::unique_ptr timesteps; extern double *rpkt_emiss; +extern int bfestimcount; + // for USE_LUT_PHOTOION = true extern double *corrphotoionrenorm; extern double *gammaestimator; diff --git a/grid.cc b/grid.cc index d0db59a01..7f5d20627 100644 --- a/grid.cc +++ b/grid.cc @@ -2425,8 +2425,8 @@ static auto get_coordboundary_distances_cylindrical2d(std::span } [[nodiscard]] auto boundary_distance(std::span dir, std::span pos, - const double tstart, int cellindex, int *snext, enum cell_boundary *pkt_last_cross) - -> double + const double tstart, int cellindex, int *snext, + enum cell_boundary *pkt_last_cross) -> double /// Basic routine to compute distance to a cell boundary. { if constexpr (FORCE_SPHERICAL_ESCAPE_SURFACE) { diff --git a/input.cc b/input.cc index f4f19ac15..e886db7a6 100644 --- a/input.cc +++ b/input.cc @@ -1118,8 +1118,8 @@ static void read_atomicdata_files() { printout("cont_index %d\n", cont_index); } -static auto search_groundphixslist(double nu_edge, int *index_in_groundlevelcontestimator, int el, int in, int ll) - -> int +static auto search_groundphixslist(double nu_edge, int *index_in_groundlevelcontestimator, int el, int in, + int ll) -> int /// Return the closest ground level continuum index to the given edge /// frequency. If the given edge frequency is redder than the reddest /// continuum return -1. @@ -1456,6 +1456,7 @@ static void setup_phixs_list() { globals::nbfcontinua * (sizeof(fullphixslist)) / 1024. / 1024.); size_t nbftables = 0; int allcontindex = 0; + globals::bfestimcount = 0; for (int element = 0; element < get_nelements(); element++) { const int nions = get_nions(element); for (int ion = 0; ion < nions - 1; ion++) { @@ -1482,6 +1483,13 @@ static void setup_phixs_list() { nonconstallcont[allcontindex].probability = get_phixsprobability(element, ion, level, phixstargetindex); nonconstallcont[allcontindex].upperlevel = get_phixsupperlevel(element, ion, level, phixstargetindex); + if (LEVEL_HAS_BFEST(get_atomicnumber(element), get_ionstage(element, ion), level)) { + nonconstallcont[allcontindex].bfestimindex = globals::bfestimcount; + globals::bfestimcount++; + } else { + nonconstallcont[allcontindex].bfestimindex = -1; + } + if constexpr (USE_LUT_PHOTOION || USE_LUT_BFHEATING) { int index_in_groundlevelcontestimator = 0; nonconstallcont[allcontindex].index_in_groundphixslist = @@ -1495,6 +1503,7 @@ static void setup_phixs_list() { } } } + printout("[info] BF estimators activated for %d photoionisation transitions\n", globals::bfestimcount); assert_always(allcontindex == globals::nbfcontinua); assert_always(globals::nbfcontinua >= 0); // was initialised as -1 before startup @@ -2189,4 +2198,4 @@ void write_timestep_file() { globals::timesteps[n].width / DAY); } fclose(timestepfile); -} \ No newline at end of file +} diff --git a/ltepop.cc b/ltepop.cc index 376207534..e811e5550 100644 --- a/ltepop.cc +++ b/ltepop.cc @@ -241,8 +241,8 @@ auto calculate_levelpop_lte(int modelgridindex, int element, int ion, int level) exp(-E_aboveground / KB / T_exc)); } -static auto calculate_levelpop_nominpop(int modelgridindex, int element, int ion, int level, bool *skipminpop) - -> double { +static auto calculate_levelpop_nominpop(int modelgridindex, int element, int ion, int level, + bool *skipminpop) -> double { assert_testmodeonly(modelgridindex < grid::get_npts_model()); assert_testmodeonly(element < get_nelements()); assert_testmodeonly(ion < get_nions(element)); @@ -424,8 +424,8 @@ auto get_nnion(int modelgridindex, int element, int ion) -> double grid::modelgrid[modelgridindex].composition[element].partfunct[ion] / stat_weight(element, ion, 0); } -static auto find_uppermost_ion(const int modelgridindex, const int element, const double nne_hi, const bool force_lte) - -> int { +static auto find_uppermost_ion(const int modelgridindex, const int element, const double nne_hi, + const bool force_lte) -> int { const int nions = get_nions(element); if (nions == 0) { return -1; diff --git a/ltepop.h b/ltepop.h index cddec48fa..d48f165d3 100644 --- a/ltepop.h +++ b/ltepop.h @@ -11,13 +11,13 @@ auto get_groundlevelpop(int modelgridindex, int element, int ion) -> double; auto calculate_levelpop(int modelgridindex, int element, int ion, int level) -> double; auto calculate_levelpop_lte(int modelgridindex, int element, int ion, int level) -> double; auto get_levelpop(int modelgridindex, int element, int ion, int level) -> double; -[[nodiscard]] auto calculate_sahafact(int element, int ion, int level, int upperionlevel, double T, double E_threshold) - -> double; +[[nodiscard]] auto calculate_sahafact(int element, int ion, int level, int upperionlevel, double T, + double E_threshold) -> double; [[nodiscard]] auto get_nnion(int modelgridindex, int element, int ion) -> double; void calculate_ion_balance_nne(int modelgridindex); void calculate_cellpartfuncts(int modelgridindex, int element); -[[nodiscard]] auto calculate_ionfractions(int element, int modelgridindex, double nne, bool use_phi_lte) - -> std::vector; +[[nodiscard]] auto calculate_ionfractions(int element, int modelgridindex, double nne, + bool use_phi_lte) -> std::vector; void set_groundlevelpops(int modelgridindex, int element, float nne, bool force_lte); #endif // LTEPOP_H diff --git a/macroatom.cc b/macroatom.cc index b0c52af7a..7cc647cc0 100644 --- a/macroatom.cc +++ b/macroatom.cc @@ -744,8 +744,8 @@ auto rad_deexcitation_ratecoeff(const int modelgridindex, const int element, con } auto rad_excitation_ratecoeff(const int modelgridindex, const int element, const int ion, const int lower, - const int uptransindex, const double epsilon_trans, int lineindex, const double t_current) - -> double + const int uptransindex, const double epsilon_trans, int lineindex, + const double t_current) -> double /// radiative excitation rate: paperII 3.5.2 // multiply by lower level population to get a rate per second { diff --git a/macroatom.h b/macroatom.h index b3df8b56b..f43b788ac 100644 --- a/macroatom.h +++ b/macroatom.h @@ -40,8 +40,8 @@ auto rad_excitation_ratecoeff(int modelgridindex, int element, int ion, int lowe double epsilon_trans, int lineindex, double t_current) -> double; auto rad_recombination_ratecoeff(float T_e, float nne, int element, int upperion, int upperionlevel, int lowerionlevel, int modelgridindex) -> double; -auto stim_recombination_ratecoeff(float nne, int element, int upperion, int upper, int lower, int modelgridindex) - -> double; +auto stim_recombination_ratecoeff(float nne, int element, int upperion, int upper, int lower, + int modelgridindex) -> double; auto col_recombination_ratecoeff(int modelgridindex, int element, int upperion, int upper, int lower, double epsilon_trans) -> double; diff --git a/nltepop.cc b/nltepop.cc index 76f2d038f..a4fd570c4 100644 --- a/nltepop.cc +++ b/nltepop.cc @@ -119,8 +119,8 @@ static void filter_nlte_matrix(const int element, gsl_matrix *rate_matrix, gsl_v } static auto get_total_rate(const int index_selected, const gsl_matrix *rate_matrix, const gsl_vector *popvec, - const bool into_level, const bool only_levels_below, const bool only_levels_above) - -> double { + const bool into_level, const bool only_levels_below, + const bool only_levels_above) -> double { double total_rate = 0.; assert_always(!only_levels_below || !only_levels_above); @@ -174,8 +174,8 @@ static auto get_total_rate_in(const int index_to, const gsl_matrix *rate_matrix, return get_total_rate(index_to, rate_matrix, popvec, true, false, false); } -static auto get_total_rate_out(const int index_from, const gsl_matrix *rate_matrix, const gsl_vector *popvec) - -> double { +static auto get_total_rate_out(const int index_from, const gsl_matrix *rate_matrix, + const gsl_vector *popvec) -> double { return get_total_rate(index_from, rate_matrix, popvec, false, false, false); } @@ -393,7 +393,7 @@ static void nltepop_reset_element(const int modelgridindex, const int element) { for (int ion = 0; ion < nions; ion++) { const int nlte_start = globals::elements[element].ions[ion].first_nlte; const int nlevels_nlte = get_nlevels_nlte(element, ion); - for (int level = 1; level < nlevels_nlte; level++) { + for (int level = 1; level <= nlevels_nlte; level++) { grid::modelgrid[modelgridindex].nlte_pops[nlte_start + level - 1] = -1.0; // flag to indicate no useful data } diff --git a/nonthermal.cc b/nonthermal.cc index 907f5b424..28b57c94f 100644 --- a/nonthermal.cc +++ b/nonthermal.cc @@ -1551,8 +1551,8 @@ auto nt_ionisation_maxupperion(const int element, const int lowerion) -> int { return maxupper; } -auto nt_random_upperion(const int modelgridindex, const int element, const int lowerion, const bool energyweighted) - -> int { +auto nt_random_upperion(const int modelgridindex, const int element, const int lowerion, + const bool energyweighted) -> int { assert_testmodeonly(lowerion < get_nions(element) - 1); if (NT_SOLVE_SPENCERFANO && NT_MAX_AUGER_ELECTRONS > 0) { while (true) { @@ -1610,8 +1610,8 @@ auto nt_ionization_ratecoeff(const int modelgridindex, const int element, const static auto calculate_nt_excitation_ratecoeff_perdeposition(const int modelgridindex, const int element, const int ion, const int lower, const int uptransindex, - const double statweight_lower, const double epsilon_trans) - -> double + const double statweight_lower, + const double epsilon_trans) -> double // Kozma & Fransson equation 9 divided by level population and epsilon_trans { if (nt_solution[modelgridindex].yfunc == nullptr) { @@ -1962,8 +1962,8 @@ static void analyse_sf_solution(const int modelgridindex, const int timestep, co (excitationindex)++; } } // NT_EXCITATION_ON - } // for t - } // for lower + } // for t + } // for lower printout(" frac_excitation: %g\n", frac_excitation_ion); if (frac_excitation_ion > 1. || !std::isfinite(frac_excitation_ion)) { diff --git a/radfield.cc b/radfield.cc index dfb5749bb..29f89eb90 100644 --- a/radfield.cc +++ b/radfield.cc @@ -7,9 +7,11 @@ #include #include +#include #include #include "atomic.h" +#include "globals.h" #include "grid.h" #include "sn3d.h" #include "vectors.h" @@ -354,7 +356,7 @@ void init(int my_rank, int ndo_nonempty) if (globals::rank_in_node == 0) { my_rank_cells += nonempty_npts_model - (my_rank_cells * globals::node_nprocs); } - auto size = static_cast(my_rank_cells * globals::nbfcontinua * sizeof(float)); + auto size = static_cast(my_rank_cells * globals::bfestimcount * sizeof(float)); int disp_unit = sizeof(float); MPI_Win_allocate_shared(size, disp_unit, MPI_INFO_NULL, globals::mpi_comm_node, &prev_bfrate_normed, &win_prev_bfrate_normed); @@ -362,16 +364,16 @@ void init(int my_rank, int ndo_nonempty) } #else { - prev_bfrate_normed = static_cast(malloc(nonempty_npts_model * globals::nbfcontinua * sizeof(float))); + prev_bfrate_normed = static_cast(malloc(nonempty_npts_model * globals::bfestimcount * sizeof(float))); } #endif printout("[info] mem_usage: detailed bf estimators for non-empty cells occupy %.3f MB (node shared memory)\n", - nonempty_npts_model * globals::nbfcontinua * sizeof(float) / 1024. / 1024.); + nonempty_npts_model * globals::bfestimcount * sizeof(float) / 1024. / 1024.); - bfrate_raw = static_cast(malloc(nonempty_npts_model * globals::nbfcontinua * sizeof(double))); + bfrate_raw = static_cast(malloc(nonempty_npts_model * globals::bfestimcount * sizeof(double))); printout("[info] mem_usage: detailed bf estimator acculumators for non-empty cells occupy %.3f MB\n", - nonempty_npts_model * globals::nbfcontinua * sizeof(double) / 1024. / 1024.); + nonempty_npts_model * globals::bfestimcount * sizeof(double) / 1024. / 1024.); } for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { @@ -380,9 +382,9 @@ void init(int my_rank, int ndo_nonempty) if constexpr (MULTIBIN_RADFIELD_MODEL_ON) { if (globals::rank_in_node == 0) { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; radfieldbin_solutions[mgibinindex].W = -1.; radfieldbin_solutions[mgibinindex].T_R = -1.; } @@ -473,23 +475,21 @@ auto get_Jb_lu_contribcount(const int modelgridindex, const int jblueindex) -> i static auto get_bin_J(int modelgridindex, int binindex) -> double // get the normalised J_nu { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); assert_testmodeonly(J_normfactor[nonemptymgi] > 0.0); assert_testmodeonly(modelgridindex < grid::get_npts_model()); assert_testmodeonly(binindex >= 0); assert_testmodeonly(binindex < RADFIELDBINCOUNT); - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; - return radfieldbins[mgibinindex].J_raw * J_normfactor[nonemptymgi]; + return radfieldbins[nonemptymgi * RADFIELDBINCOUNT + binindex].J_raw * J_normfactor[nonemptymgi]; } static auto get_bin_nuJ(int modelgridindex, int binindex) -> double { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); assert_testmodeonly(J_normfactor[nonemptymgi] > 0.0); assert_testmodeonly(modelgridindex < grid::get_npts_model()); assert_testmodeonly(binindex >= 0); assert_testmodeonly(binindex < RADFIELDBINCOUNT); - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; - return radfieldbins[mgibinindex].nuJ_raw * J_normfactor[nonemptymgi]; + return radfieldbins[nonemptymgi * RADFIELDBINCOUNT + binindex].nuJ_raw * J_normfactor[nonemptymgi]; } static inline auto get_bin_nu_bar(int modelgridindex, int binindex) -> double @@ -508,18 +508,18 @@ static inline auto get_bin_nu_lower(int binindex) -> double { } static inline auto get_bin_contribcount(int modelgridindex, int binindex) -> int { - const int mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; - return radfieldbins[mgibinindex].contribcount; + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + return radfieldbins[nonemptymgi * RADFIELDBINCOUNT + binindex].contribcount; } static inline auto get_bin_W(int modelgridindex, int binindex) -> float { - const int mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; - return radfieldbin_solutions[mgibinindex].W; + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + return radfieldbin_solutions[nonemptymgi * RADFIELDBINCOUNT + binindex].W; } static inline auto get_bin_T_R(int modelgridindex, int binindex) -> float { - const int mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; - return radfieldbin_solutions[mgibinindex].T_R; + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + return radfieldbin_solutions[nonemptymgi * RADFIELDBINCOUNT + binindex].T_R; } static inline auto select_bin(double nu) -> int { @@ -540,7 +540,7 @@ static inline auto select_bin(double nu) -> int { void write_to_file(int modelgridindex, int timestep) { assert_always(MULTIBIN_RADFIELD_MODEL_ON); - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); #ifdef _OPENMP #pragma omp critical(out_file) { @@ -650,13 +650,13 @@ void zero_estimators(int modelgridindex) return; } - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); if constexpr (DETAILED_BF_ESTIMATORS_ON) { assert_always(bfrate_raw != nullptr); if (grid::get_numassociatedcells(modelgridindex) > 0) { - for (int i = 0; i < globals::nbfcontinua; i++) { - bfrate_raw[nonemptymgi * globals::nbfcontinua + i] = 0.; + for (int i = 0; i < globals::bfestimcount; i++) { + bfrate_raw[nonemptymgi * globals::bfestimcount + i] = 0.; } } } @@ -678,7 +678,7 @@ void zero_estimators(int modelgridindex) assert_always(radfieldbins != nullptr); for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; radfieldbins[mgibinindex].J_raw = 0.0; radfieldbins[mgibinindex].nuJ_raw = 0.0; radfieldbins[mgibinindex].contribcount = 0; @@ -708,8 +708,13 @@ static void update_bfestimators(const int nonemptymgi, const double distance_e_c const double nu_max_phixs = nu_edge * last_phixs_nuovernuedge; // nu of the uppermost point in the phixs table if (nu_cmf >= nu_edge && nu_cmf <= nu_max_phixs) { - safeadd(bfrate_raw[nonemptymgi * nbfcontinua + allcontindex], - globals::phixslist[tid].gamma_contr[allcontindex] * distance_e_cmf_over_nu); + int bf_estimator = globals::allcont[allcontindex].bfestimindex; + if (bf_estimator >= 0) { + assert_testmodeonly(bf_estimator < globals::bfestimcount); + + safeadd(bfrate_raw[nonemptymgi * globals::bfestimcount + bf_estimator], + globals::phixslist[tid].gamma_contr[allcontindex] * distance_e_cmf_over_nu); + } } else if (nu_cmf < nu_edge) { // list is sorted by nu_edge, so all remaining will have nu_cmf < nu_edge @@ -731,7 +736,7 @@ void update_estimators(const int nonemptymgi, const double distance_e_cmf, const const int binindex = select_bin(nu_cmf); if (binindex >= 0) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; safeadd(radfieldbins[mgibinindex].J_raw, distance_e_cmf); safeadd(radfieldbins[mgibinindex].nuJ_raw, distance_e_cmf * nu_cmf); safeincrement(radfieldbins[mgibinindex].contribcount); @@ -769,7 +774,7 @@ auto radfield(double nu, int modelgridindex) -> double if (globals::timestep >= FIRST_NLTE_RADFIELD_TIMESTEP) { const int binindex = select_bin(nu); if (binindex >= 0) { - const int mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = grid::get_modelcell_nonemptymgi(modelgridindex) * RADFIELDBINCOUNT + binindex; const struct radfieldbin_solution *const bin = &radfieldbin_solutions[mgibinindex]; if (bin->W >= 0.) { const double J_nu = dbb(nu, bin->T_R, bin->W); @@ -1114,7 +1119,7 @@ void fit_parameters(int modelgridindex, int timestep) // T_R_bin = -1; // W_bin = -1; // } - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; radfieldbin_solutions[mgibinindex].T_R = T_R_bin; radfieldbin_solutions[mgibinindex].W = W_bin; } @@ -1159,11 +1164,11 @@ void normalise_J(const int modelgridindex, const double estimator_normfactor_ove void normalise_bf_estimators(const int modelgridindex, const double estimator_normfactor_over_H) { if constexpr (DETAILED_BF_ESTIMATORS_ON) { printout("normalise_bf_estimators for cell %d with factor %g\n", modelgridindex, estimator_normfactor_over_H); - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); assert_always(nonemptymgi >= 0); - for (int i = 0; i < globals::nbfcontinua; i++) { - const int mgibfindex = nonemptymgi * globals::nbfcontinua + i; - prev_bfrate_normed[mgibfindex] = bfrate_raw[mgibfindex] * estimator_normfactor_over_H; + for (int i = 0; i < globals::bfestimcount; i++) { + const auto detailed_mgibfindex = nonemptymgi * globals::bfestimcount + i; + prev_bfrate_normed[detailed_mgibfindex] = bfrate_raw[detailed_mgibfindex] * estimator_normfactor_over_H; } } } @@ -1190,10 +1195,11 @@ auto get_bfrate_estimator(const int element, const int lowerion, const int lower if constexpr (!DETAILED_BF_ESTIMATORS_ON) { return -1; } else { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); const int allcontindex = get_bfcontindex(element, lowerion, lower, phixstargetindex); if (allcontindex >= 0) { - return prev_bfrate_normed[nonemptymgi * globals::nbfcontinua + allcontindex]; + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const int bfestimindex = globals::allcont[allcontindex].bfestimindex; + return (bfestimindex >= 0) ? prev_bfrate_normed[nonemptymgi * globals::bfestimcount + bfestimindex] : 0.; } printout("no bf rate for element Z=%d ion_stage %d lower %d phixstargetindex %d\n", get_atomicnumber(element), @@ -1257,8 +1263,8 @@ void reduce_estimators() MPI_Allreduce(MPI_IN_PLACE, nuJ.data(), nonempty_npts_model, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); if constexpr (DETAILED_BF_ESTIMATORS_ON) { - MPI_Allreduce(MPI_IN_PLACE, bfrate_raw, nonempty_npts_model * globals::nbfcontinua, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, bfrate_raw, grid::get_nonempty_npts_model() * globals::bfestimcount, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); } if constexpr (MULTIBIN_RADFIELD_MODEL_ON) { @@ -1268,7 +1274,7 @@ void reduce_estimators() for (int nonemptymgi = 0; nonemptymgi < nonempty_npts_model; nonemptymgi++) { for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; MPI_Allreduce(MPI_IN_PLACE, &radfieldbins[mgibinindex].J_raw, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE, &radfieldbins[mgibinindex].nuJ_raw, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE, &radfieldbins[mgibinindex].contribcount, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); @@ -1306,12 +1312,12 @@ void do_MPI_Bcast(const int modelgridindex, const int root, int root_node_id) return; } - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); MPI_Bcast(&J_normfactor[nonemptymgi], 1, MPI_DOUBLE, root, MPI_COMM_WORLD); if constexpr (MULTIBIN_RADFIELD_MODEL_ON) { for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; if (globals::rank_in_node == 0) { MPI_Bcast(&radfieldbin_solutions[mgibinindex].W, 1, MPI_FLOAT, root_node_id, globals::mpi_comm_internode); MPI_Bcast(&radfieldbin_solutions[mgibinindex].T_R, 1, MPI_FLOAT, root_node_id, globals::mpi_comm_internode); @@ -1324,8 +1330,8 @@ void do_MPI_Bcast(const int modelgridindex, const int root, int root_node_id) if constexpr (DETAILED_BF_ESTIMATORS_ON) { if (globals::rank_in_node == 0) { - MPI_Bcast(&prev_bfrate_normed[nonemptymgi * globals::nbfcontinua], globals::nbfcontinua, MPI_FLOAT, root_node_id, - globals::mpi_comm_internode); + MPI_Bcast(&prev_bfrate_normed[nonemptymgi * globals::bfestimcount], globals::bfestimcount, MPI_FLOAT, + root_node_id, globals::mpi_comm_internode); } } @@ -1358,12 +1364,15 @@ void write_restart_data(FILE *gridsave_file) { const int nbfcontinua = globals::nbfcontinua; fprintf(gridsave_file, "%d\n", nbfcontinua); + const int bfestimcount = globals::bfestimcount; + fprintf(gridsave_file, "%d\n", bfestimcount); + for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { if (grid::get_numassociatedcells(modelgridindex) > 0) { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); fprintf(gridsave_file, "%d\n", modelgridindex); - for (int i = 0; i < nbfcontinua; i++) { - fprintf(gridsave_file, "%a ", prev_bfrate_normed[nonemptymgi * nbfcontinua + i]); + for (int i = 0; i < bfestimcount; i++) { + fprintf(gridsave_file, "%a ", prev_bfrate_normed[nonemptymgi * bfestimcount + i]); } } } @@ -1379,13 +1388,13 @@ void write_restart_data(FILE *gridsave_file) { for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { if (grid::get_numassociatedcells(modelgridindex) > 0) { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); assert_testmodeonly(nonemptymgi >= 0); fprintf(gridsave_file, "%d %la\n", modelgridindex, J_normfactor[nonemptymgi]); if constexpr (MULTIBIN_RADFIELD_MODEL_ON) { for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; fprintf(gridsave_file, "%la %la %a %a %d\n", radfieldbins[mgibinindex].J_raw, radfieldbins[mgibinindex].nuJ_raw, radfieldbin_solutions[mgibinindex].W, radfieldbin_solutions[mgibinindex].T_R, radfieldbins[mgibinindex].contribcount); @@ -1456,17 +1465,21 @@ void read_restart_data(FILE *gridsave_file) { assert_always(fscanf(gridsave_file, "%d\n", &gridsave_nbf_in) == 1); assert_always(gridsave_nbf_in == globals::nbfcontinua); + int gridsave_nbfestim_in = 0; + assert_always(fscanf(gridsave_file, "%d\n", &gridsave_nbfestim_in) == 1); + assert_always(gridsave_nbfestim_in == globals::bfestimcount); + for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) { if (grid::get_numassociatedcells(modelgridindex) > 0) { - const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); + const ptrdiff_t nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex); int mgi_in = 0; assert_always(fscanf(gridsave_file, "%d\n", &mgi_in) == 1); assert_always(mgi_in == modelgridindex); - for (int i = 0; i < globals::nbfcontinua; i++) { + for (int i = 0; i < globals::bfestimcount; i++) { float bfrate_normed = 0; assert_always(fscanf(gridsave_file, "%a ", &bfrate_normed) == 1); - const int mgibfindex = nonemptymgi * globals::nbfcontinua + i; + const auto mgibfindex = nonemptymgi * globals::bfestimcount + i; #ifdef MPI_ON if (globals::rank_in_node == 0) #endif @@ -1505,7 +1518,7 @@ void read_restart_data(FILE *gridsave_file) { if constexpr (MULTIBIN_RADFIELD_MODEL_ON) { for (int binindex = 0; binindex < RADFIELDBINCOUNT; binindex++) { - const int mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; + const ptrdiff_t mgibinindex = nonemptymgi * RADFIELDBINCOUNT + binindex; float W = 0; float T_R = 0; assert_always(fscanf(gridsave_file, "%la %la %a %a %d\n", &radfieldbins[mgibinindex].J_raw, diff --git a/ratecoeff.cc b/ratecoeff.cc index bbdf776d8..c7420387c 100644 --- a/ratecoeff.cc +++ b/ratecoeff.cc @@ -709,8 +709,8 @@ auto get_spontrecombcoeff(int element, int ion, int level, int phixstargetindex, auto calculate_ionrecombcoeff(const int modelgridindex, const float T_e, const int element, const int upperion, const bool assume_lte, const bool collisional_not_radiative, const bool printdebug, - const bool lower_superlevel_only, const bool per_groundmultipletpop, const bool stimonly) - -> double + const bool lower_superlevel_only, const bool per_groundmultipletpop, + const bool stimonly) -> double // multiply by upper ion population (or ground population if per_groundmultipletpop is true) and nne to get a rate { const int lowerion = upperion - 1; diff --git a/rpkt.h b/rpkt.h index 55f5ea3e3..e32dc23a9 100644 --- a/rpkt.h +++ b/rpkt.h @@ -33,8 +33,8 @@ constexpr auto get_linedistance(const double prop_time, const double nu_cmf, con return CLIGHT * prop_time * (nu_cmf / nu_trans - 1); } -[[nodiscard]] inline auto get_ionestimindex_nonemptymgi(const int nonemptymgi, const int element, const int ion) - -> int { +[[nodiscard]] inline auto get_ionestimindex_nonemptymgi(const int nonemptymgi, const int element, + const int ion) -> int { assert_testmodeonly(ion >= 0); assert_testmodeonly(ion < get_nions(element) - 1); return nonemptymgi * get_includedions() + get_uniqueionindex(element, ion); diff --git a/stats.cc b/stats.cc index f584e2ec8..8ba380327 100644 --- a/stats.cc +++ b/stats.cc @@ -110,8 +110,8 @@ void increment_ion_stats_contabsorption(const struct packet *const pkt_ptr, cons } } -auto get_ion_stats(const int modelgridindex, const int element, const int ion, enum ionstattypes ionstattype) - -> double { +auto get_ion_stats(const int modelgridindex, const int element, const int ion, + enum ionstattypes ionstattype) -> double { assert_always(ion < get_nions(element)); assert_always(ionstattype < ION_STAT_COUNT); const int uniqueionindex = get_uniqueionindex(element, ion); diff --git a/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_final.txt b/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_final.txt new file mode 100644 index 000000000..9a425378d --- /dev/null +++ b/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_final.txt @@ -0,0 +1,20 @@ +77a8cefec72da8c0bcebdf4e04d2fcf5 absorption.out +aff44ecad68d986f04fcd6110331ba8b bflist.out +9cfd74fab2244cafd4bce62dc6fe2ad7 deposition.out +986056dc9a141a00e7b62e653b2b4f4b emission.out +8607d1185f06e29637171e28955097e1 emissiontrue.out +9aacf70462a00d805820041d07deaa02 gamma_light_curve.out +035228757b6e30f01ff6a04a6a94d364 gamma_spec.out +057b226c371f3819cba5e04bfea3d114 gammalinelist.out +47dd2e31cc98d5d1dda67520fb399429 grid.out +2c51d13a376994eb179c6e97d44f1dfd light_curve.out +9773becefdfe4afffb041b703210c67c linestat.out +98195ccb920b6831d137861c8c6cbfa7 modelgridrankassignments.out +ce06681bb20c3a895bc724caa303cb54 packets00_0000.out +195ea92baa7650d5ae6cab37393a1483 packets00_0001.out +7ff320b87eb952c5e24295c2e12a6258 spec.out +8e7163982f1aa938dc1505478b8c60d1 timesteps.out +463e7d1a82b94f8c0fa7083c1f9b59d7 job1/estimators_0000.out +ebc4df506dc9890bde3a762f581bd5ad job1/nlte_0000.out +1c4b1ccb2ea3bd56f3ef976e4dce870e job1/nonthermalspec_0000.out +ef9c7843dfc4b7dc5e6908270a6bd780 job1/radfield_0000.out \ No newline at end of file diff --git a/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_job0.txt b/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_job0.txt new file mode 100644 index 000000000..fe98f1ce9 --- /dev/null +++ b/tests/nebularonezone_1d_3dgrid_limitbfest_inputfiles/results_md5_job0.txt @@ -0,0 +1,19 @@ +f5e5d8fc1f4eb7d1d2e7624e6929954c absorption.out +aff44ecad68d986f04fcd6110331ba8b bflist.out +b1c057449d45d3eafc42a30c2530cb92 deposition.out +91195eb22d4b77ed93bfc5edc96fbd79 emission.out +6bdc6267ae27feab57a34c8a91f78c37 emissiontrue.out +59e10adcf5e9d1c05f935e5583ad9e1e gamma_light_curve.out +057b226c371f3819cba5e04bfea3d114 gammalinelist.out +47dd2e31cc98d5d1dda67520fb399429 grid.out +e0d00b6a32be200c5cf0e9962089c0f3 light_curve.out +9773becefdfe4afffb041b703210c67c linestat.out +98195ccb920b6831d137861c8c6cbfa7 modelgridrankassignments.out +353ee551590fbe95446d10f2554adb9c packets00_0000.out +58ae1da9bf04c21767d40fd722bf737e packets00_0001.out +65539a0900e70d998038852567575617 spec.out +8e7163982f1aa938dc1505478b8c60d1 timesteps.out +2a274951a4fff75b7962281cf7bd043c job0/estimators_0000.out +7d86d86b36f9208ba00af17f31846025 job0/nlte_0000.out +1c4b1ccb2ea3bd56f3ef976e4dce870e job0/nonthermalspec_0000.out +fc10e7d1fd915e9372566c536d9d4652 job0/radfield_0000.out \ No newline at end of file diff --git a/tests/setup_nebularonezone_1d_3dgrid_limitbfest.sh b/tests/setup_nebularonezone_1d_3dgrid_limitbfest.sh new file mode 100755 index 000000000..a2b0c6083 --- /dev/null +++ b/tests/setup_nebularonezone_1d_3dgrid_limitbfest.sh @@ -0,0 +1,45 @@ +#!/usr/bin/env zsh + +set -x + +runfolder=nebularonezone_1d_3dgrid_limitbfest_testrun + +rsync -av nebularonezone_1d_3dgrid_inputfiles/ nebularonezone_1d_3dgrid_limitbfest_testrun/ + +rsync -av nebularonezone_1d_3dgrid_limitbfest_inputfiles/ nebularonezone_1d_3dgrid_limitbfest_testrun/ + +if [ ! -f atomicdata_feconi.tar.xz ]; then curl -O https://theory.gsi.de/~lshingle/artis_http_public/artis/atomicdata_feconi.tar.xz; fi + +tar -xf atomicdata_feconi.tar.xz --directory nebularonezone_1d_3dgrid_limitbfest_testrun/ + +cp ../data/* nebularonezone_1d_3dgrid_limitbfest_testrun/ + +cp ../artisoptions_nltenebular.h nebularonezone_1d_3dgrid_limitbfest_testrun/artisoptions.h + +cd nebularonezone_1d_3dgrid_limitbfest_testrun + +sed -i'' -e 's/constexpr int MPKTS.*/constexpr int MPKTS = 1000000;/g' artisoptions.h + +sed -i'' -e 's/constexpr int GRID_TYPE.*/constexpr int GRID_TYPE = GRID_CARTESIAN3D;/g' artisoptions.h + +sed -i'' -e 's/constexpr int CUBOID_NCOORDGRID_X.*/constexpr int CUBOID_NCOORDGRID_X = 50;/g' artisoptions.h +sed -i'' -e 's/constexpr int CUBOID_NCOORDGRID_Y.*/constexpr int CUBOID_NCOORDGRID_Y = 50;/g' artisoptions.h +sed -i'' -e 's/constexpr int CUBOID_NCOORDGRID_Z.*/constexpr int CUBOID_NCOORDGRID_Z = 50;/g' artisoptions.h + +sed -i'' -e 's/constexpr int TABLESIZE.*/constexpr int TABLESIZE = 20;/g' artisoptions.h +sed -i'' -e 's/constexpr double MINTEMP.*/constexpr double MINTEMP = 2000.;/g' artisoptions.h +sed -i'' -e 's/constexpr double MAXTEMP.*/constexpr double MAXTEMP = 10000.;/g' artisoptions.h + +sed -i'' -e 's/constexpr int FIRST_NLTE_RADFIELD_TIMESTEP.*/constexpr int FIRST_NLTE_RADFIELD_TIMESTEP = 7;/g' artisoptions.h + +sed -i'' -e 's/constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP.*/constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 7;/g' artisoptions.h + +sed -i'' -e 's/constexpr bool SF_AUGER_CONTRIBUTION_ON.*/constexpr bool SF_AUGER_CONTRIBUTION_ON = false;/g' artisoptions.h + +sed -i'' -e 's/constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC.*/constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = true;/g' artisoptions.h + +sed -i'' -e 's/ \/\/ return LEVEL_IS_NLTE.*/return LEVEL_IS_NLTE(element_z, ionstage, level);/g' artisoptions.h + +cd - + +set +x diff --git a/thermalbalance.cc b/thermalbalance.cc index 41518c48d..d733285fe 100644 --- a/thermalbalance.cc +++ b/thermalbalance.cc @@ -76,8 +76,8 @@ static auto integrand_bfheatingcoeff_custom_radfield(double nu, void *voidparas) return sigma_bf * (1 - nu_edge / nu) * radfield::radfield(nu, modelgridindex) * (1 - exp(-HOVERKB * nu / T_R)); } -static auto calculate_bfheatingcoeff(int element, int ion, int level, int phixstargetindex, int modelgridindex) - -> double { +static auto calculate_bfheatingcoeff(int element, int ion, int level, int phixstargetindex, + int modelgridindex) -> double { double error = 0.0; const double epsrel = 1e-3; const double epsrelwarning = 1e-1; diff --git a/vpkt.cc b/vpkt.cc index 0f5bb28fb..fedde4168 100644 --- a/vpkt.cc +++ b/vpkt.cc @@ -933,8 +933,8 @@ auto vpkt_call_estimators(struct packet *pkt_ptr, const enum packet_type realtyp } } -auto rot_angle(std::span n1, std::span n2, std::span ref1, std::span ref2) - -> double { +auto rot_angle(std::span n1, std::span n2, std::span ref1, + std::span ref2) -> double { // Rotation angle from the scattering plane // We need to rotate Stokes Parameters to (or from) the scattering plane from (or to) // the meridian frame such that Q=1 is in the scattering plane and along ref1