diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index c9a7a314b..04bc34050 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -11,23 +11,34 @@ jobs: strategy: fail-fast: false matrix: - gpu: ["amd", "nvidia"] - precision: ["double", "single"] - runs-on: [self-hosted, "${{ matrix.gpu }}-gpu"] + device: [amd-gpu, nvidia-gpu] + precision: [double, single] + exclude: + - device: amd-gpu + precision: double + # my AMD GPUs doesn't support fp64 atomics : ( + runs-on: [self-hosted, "${{ matrix.device }}"] steps: - name: Checkout uses: actions/checkout@v3.3.0 - name: Configure run: | - if [ "${{ matrix.gpu }}" = "nvidia" ]; then - FLAGS="-D Kokkos_ENABLE_CUDA=ON -D Kokkos_ARCH_AMPERE86=ON" - elif [ "${{ matrix.gpu }}" = "amd" ]; then + if [ "${{ matrix.device }}" = "nvidia-gpu" ]; then + FLAGS="-D Kokkos_ENABLE_CUDA=ON" + if [[ ! -z $(nvidia-smi | grep "V100") ]]; then + FLAGS+=" -D Kokkos_ARCH_VOLTA70=ON" + elif [[ ! -z $(nvidia-smi | grep "A100") ]]; then + FLAGS+=" -D Kokkos_ARCH_AMPERE80=ON" + else + FLAGS+=" -D Kokkos_ARCH_AMPERE86=ON" + fi + elif [ "${{ matrix.device }}" = "amd-gpu" ]; then FLAGS="-D Kokkos_ENABLE_HIP=ON -D Kokkos_ARCH_AMD_GFX1100=ON" fi cmake -B build -D TESTS=ON -D output=ON -D precision=${{ matrix.precision }} $FLAGS - name: Compile run: | - cmake --build build -j $(exec nproc) + cmake --build build -j $(nproc) - name: Run tests run: | - ctest --test-dir build --output-on-failure --verbose \ No newline at end of file + ctest --test-dir build --output-on-failure --verbose diff --git a/.gitignore b/.gitignore index 76f69e364..20bfe33a3 100644 --- a/.gitignore +++ b/.gitignore @@ -53,3 +53,4 @@ Testing/ .schema.json *_old/ +action-token diff --git a/cmake/report.cmake b/cmake/report.cmake index e15074943..fe914baa8 100644 --- a/cmake/report.cmake +++ b/cmake/report.cmake @@ -150,6 +150,16 @@ PrintChoices("CUDA" 0 42 ) +PrintChoices("HIP" + "Kokkos_ENABLE_HIP" + "${ON_OFF_VALUES}" + ${Kokkos_ENABLE_HIP} + "OFF" + "${Green}" + HIP_REPORT + 0 + 42 +) PrintChoices("OpenMP" "Kokkos_ENABLE_OPENMP" "${ON_OFF_VALUES}" @@ -183,6 +193,28 @@ PrintChoices("C compiler" 42 ) +get_cmake_property(_variableNames VARIABLES) +foreach (_variableName ${_variableNames}) + string(REGEX MATCH "Kokkos_ARCH_*" _isMatched ${_variableName}) + if(_isMatched) + get_property(isSet CACHE ${_variableName} PROPERTY VALUE) + if(isSet STREQUAL "ON") + string(REGEX REPLACE "Kokkos_ARCH_" "" ARCH ${_variableName}) + break() + endif() + endif() +endforeach() +PrintChoices("Architecture" + "Kokkos_ARCH_*" + "${ARCH}" + "${ARCH}" + "N/A" + "${ColorReset}" + ARCH_REPORT + 0 + 42 +) + if(${Kokkos_ENABLE_CUDA}) if("${CMAKE_CUDA_COMPILER}" STREQUAL "") execute_process(COMMAND which nvcc OUTPUT_VARIABLE CUDACOMP) @@ -194,14 +226,13 @@ if(${Kokkos_ENABLE_CUDA}) message(STATUS "CUDA compiler: ${CUDACOMP}") execute_process(COMMAND bash -c "${CUDACOMP} --version | grep release | sed -e 's/.*release //' -e 's/,.*//'" - OUTPUT_VARIABLE CUDACOMP_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) PrintChoices("CUDA compiler" "CMAKE_CUDA_COMPILER" - "${CUDACOMP} v${CUDACOMP_VERSION}" - "${CUDACOMP} v${CUDACOMP_VERSION}" + "${CUDACOMP}" + "${CUDACOMP}" "N/A" "${ColorReset}" CUDA_COMPILER_REPORT @@ -210,6 +241,12 @@ if(${Kokkos_ENABLE_CUDA}) ) endif() +if (${Kokkos_ENABLE_HIP}) + execute_process(COMMAND bash -c "hipcc --version | grep HIP | cut -d ':' -f 2 | tr -d ' '" + OUTPUT_VARIABLE ROCM_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE) +endif() + set(DOT_SYMBOL "${ColorReset}.") set(DOTTED_LINE_SYMBOL "${ColorReset}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ") @@ -242,7 +279,9 @@ message(" ${OUTPUT_REPORT}") message("${DASHED_LINE_SYMBOL} Compile configurations") +message(" ${ARCH_REPORT}") message(" ${CUDA_REPORT}") +message(" ${HIP_REPORT}") message(" ${OPENMP_REPORT}") message(" ${C_COMPILER_REPORT}") @@ -259,6 +298,11 @@ message(" ${DEBUG_REPORT}") message("${DASHED_LINE_SYMBOL}\nDependencies") +if (NOT "${CUDACOMP_VERSION}" STREQUAL "") + message(" - CUDA:\tv${CUDACOMP_VERSION}") +elseif(NOT "${ROCM_VERSION}" STREQUAL "") + message(" - ROCm:\tv${ROCM_VERSION}") +endif() message(" - Kokkos:\tv${Kokkos_VERSION}") if(${output}) message(" - ADIOS2:\tv${adios2_VERSION}") diff --git a/dev/runners/Dockerfile.runner.nvidia b/dev/runners/Dockerfile.runner.cuda similarity index 93% rename from dev/runners/Dockerfile.runner.nvidia rename to dev/runners/Dockerfile.runner.cuda index d7925ec6a..4ff132990 100644 --- a/dev/runners/Dockerfile.runner.nvidia +++ b/dev/runners/Dockerfile.runner.cuda @@ -59,14 +59,14 @@ ARG HOME=/home/$USER WORKDIR $HOME # gh runner -ARG TOKEN RUN mkdir actions-runner WORKDIR $HOME/actions-runner -RUN curl -o actions-runner-linux-x64-2.317.0.tar.gz \ +RUN --mount=type=secret,id=ghtoken \ + curl -o actions-runner-linux-x64-2.317.0.tar.gz \ -L https://github.com/actions/runner/releases/download/v2.317.0/actions-runner-linux-x64-2.317.0.tar.gz && \ tar xzf ./actions-runner-linux-x64-2.317.0.tar.gz && \ sudo ./bin/installdependencies.sh && \ - ./config.sh --url https://github.com/entity-toolkit/entity --token $TOKEN --labels nvidia-gpu + ./config.sh --url https://github.com/entity-toolkit/entity --token "$(sudo cat /run/secrets/ghtoken)" --labels nvidia-gpu ENTRYPOINT ["./run.sh"] diff --git a/dev/runners/Dockerfile.runner.rocm b/dev/runners/Dockerfile.runner.rocm new file mode 100644 index 000000000..63579dcd5 --- /dev/null +++ b/dev/runners/Dockerfile.runner.rocm @@ -0,0 +1,88 @@ +FROM rocm/rocm-terminal:latest + +USER root +ENV PATH=/opt/rocm/bin:$PATH + +ARG DEBIAN_FRONTEND=noninteractive +ENV DISPLAY=host.docker.internal:0.0 + +# upgrade +RUN apt-get update && apt-get upgrade -y + +# cmake & build tools +RUN apt-get remove -y --purge cmake && \ + apt-get install -y sudo wget curl build-essential && \ + wget "https://github.com/Kitware/CMake/releases/download/v3.29.6/cmake-3.29.6-linux-x86_64.tar.gz" -P /opt && \ + tar xvf /opt/cmake-3.29.6-linux-x86_64.tar.gz -C /opt && \ + rm /opt/cmake-3.29.6-linux-x86_64.tar.gz +ENV PATH=/opt/cmake-3.29.6-linux-x86_64/bin:$PATH + +# adios2 +RUN apt-get update && apt-get install -y git libhdf5-dev && \ + git clone https://github.com/ornladios/ADIOS2.git /opt/adios2-src && \ + cd /opt/adios2-src && \ + cmake -B build \ + -D CMAKE_CXX_STANDARD=17 \ + -D CMAKE_CXX_EXTENSIONS=OFF \ + -D CMAKE_POSITION_INDEPENDENT_CODE=TRUE \ + -D BUILD_SHARED_LIBS=ON \ + -D ADIOS2_USE_HDF5=ON \ + -D ADIOS2_USE_Python=OFF \ + -D ADIOS2_USE_Fortran=OFF \ + -D ADIOS2_USE_ZeroMQ=OFF \ + -D BUILD_TESTING=OFF \ + -D ADIOS2_BUILD_EXAMPLES=OFF \ + -D ADIOS2_USE_MPI=OFF \ + -D ADIOS2_HAVE_HDF5_VOL=OFF \ + -D CMAKE_INSTALL_PREFIX=/opt/adios2 && \ + cmake --build build -j && \ + cmake --install build && \ + rm -rf /opt/adios2-src + +ENV HDF5_ROOT=/usr +ENV ADIOS2_DIR=/opt/adios2 +ENV PATH=/opt/adios2/bin:$PATH + +# additional ROCm packages +RUN git clone -b release/rocm-rel-6.1.1.1 https://github.com/ROCm/rocThrust.git /opt/rocthrust-src && \ + git clone -b release/rocm-rel-6.1.00.36 https://github.com/ROCm/rocPRIM.git /opt/rocprim-src && \ + cd /opt/rocthrust-src && ./install --install && \ + cd /opt/rocprim-src && ./install --install && \ + rm -rf /opt/rocthrust-src /opt/rocprim-src + +ENV CMAKE_PREFIX_PATH=/opt/rocm +ENV CC=hipcc +ENV CXX=hipcc + +# cleanup +RUN apt-get clean && \ + apt-get autoclean && \ + apt-get autoremove -y && \ + rm -rf /var/lib/cache/* && \ + rm -rf /var/lib/log/* && \ + rm -rf /var/lib/apt/lists/* + +ARG USER=runner +RUN useradd -ms /usr/bin/zsh $USER && \ + usermod -aG sudo $USER && \ + echo '%sudo ALL=(ALL) NOPASSWD:ALL' >> /etc/sudoers + +USER $USER +ARG HOME=/home/$USER +WORKDIR $HOME + +# gh runner +ARG RUNNER_VERSION=2.317.0 +RUN mkdir actions-runner +WORKDIR $HOME/actions-runner + +RUN curl -o actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz \ + -L https://github.com/actions/runner/releases/download/v${RUNNER_VERSION}/actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz && \ + tar xzf ./actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz && \ + sudo ./bin/installdependencies.sh + +ADD start.sh start.sh +RUN sudo chown $USER:$USER start.sh && \ + sudo chmod +x start.sh + +ENTRYPOINT ["./start.sh"] diff --git a/dev/runners/README.md b/dev/runners/README.md index 526068bea..08d0cd176 100644 --- a/dev/runners/README.md +++ b/dev/runners/README.md @@ -4,12 +4,18 @@ GitHub allows to listen to repository changes and run the so-called "actions" (e To do that, one needs to create an image with the corresponding `Dockerfile`, and then launch a Docker container which will run in the background, listening to commands and running any actions forwarded from the GitHub. +First, you will need to obtain a runner token from the Entity GitHub repo, by going to Settings -> Actions -> Runners -> New self-hosted runner. Copy the token to use it later. The images differ slightly from the type of runner. + ### NVIDIA GPUs ```sh -# 1. Create the image -docker build -t ghrunner:nvidia -f Dockerfile.runner.nvidia . -# 2. Run a container from the image with GPU support -# ... (see wiki for instructions on NVIDIA runtime) -docker run --runtime=nvidia --gpus=all -dt ghrunner:nvidia +docker build -t ghrunner:nvidia -f Dockerfile.runner.cuda . +docker run -e TOKEN= -e LABEL=nvidia-gpu --runtime=nvidia --gpus=all -dt ghrunner:nvidia +``` + +### AMD GPUs + +```sh +docker build -t ghrunner:amd -f Dockerfile.runner.rocm . +docker run -e TOKEN= -e LABEL=amd-gpu --device=/dev/kfd --device=/dev/dri --security-opt seccomp=unconfined --group-add video -dt ghrunner:amd ``` diff --git a/dev/runners/start.sh b/dev/runners/start.sh new file mode 100644 index 000000000..c44e26f16 --- /dev/null +++ b/dev/runners/start.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +./config.sh --unattended --url https://github.com/entity-toolkit/entity --token ${TOKEN} --labels ${LABEL} + +if [[ ${LABEL} == "amd-gpu" ]]; then + echo "AMD GPU runner detected" + export HSA_OVERRIDE_GFX_VERSION=11.0.0 HIP_VISIBLE_DEVICES=0 ROCR_VISIBLE_DEVICES=0 +else + echo "Non-AMD runner" +fi + +cleanup() { + echo "Removing runner..." + ./config.sh remove --unattended --token ${TOKEN} +} + +trap 'cleanup; exit 130' INT +trap 'cleanup; exit 143' TERM + +./run.sh & +wait $! diff --git a/dev/welcome.rocm b/dev/welcome.rocm index 35e5e15d9..780375684 100644 --- a/dev/welcome.rocm +++ b/dev/welcome.rocm @@ -68,7 +68,8 @@ else python_version=${error_msg} fi -str=$(cat <${Reset} ## Testing the code * To build the tests, remove the \`build\` directory, and preconfigure with: - ${FgGreen}$ cmake -B build -D TESTS=ON -D output=ON -D Kokkos_ENABLE_HIP=ON -D Kokkos_ARCH_***=ON${Reset} + ${FgGreen}$ cmake -B build -D TESTS=ON -D output=ON -D Kokkos_ENABLE_HIP=ON -D Kokkos_ARCH_AMD_***=ON${Reset} * Then run the tests with: ${FgGreen}$ ctest --test-dir build${Reset} @@ -117,5 +119,4 @@ ${Red}Enjoy!${Reset} 🚀 EOF ) - printf "$str\n" diff --git a/docker-compose.yml b/docker-compose.yml index 3de0516a3..cb0093f4b 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -30,7 +30,7 @@ services: entity-rocm: container_name: entity_rocm devices: - - "/dev/fdk" + - "/dev/kfd" - "/dev/dri" security_opt: - seccomp:unconfined diff --git a/extern/Kokkos b/extern/Kokkos index 486cc745c..eb11070f6 160000 --- a/extern/Kokkos +++ b/extern/Kokkos @@ -1 +1 @@ -Subproject commit 486cc745cb9a287f3915061455105a3ee588c616 +Subproject commit eb11070f67565b2e660659f5207f0363bdf3b882 diff --git a/extern/adios2 b/extern/adios2 index c52d755bb..b8761e2af 160000 --- a/extern/adios2 +++ b/extern/adios2 @@ -1 +1 @@ -Subproject commit c52d755bb2ef1730a4a93b9243e2414e5a2d8456 +Subproject commit b8761e2afab2cd05b89d09b2ee4da1cd7a834225 diff --git a/extern/toml11 b/extern/toml11 index 85faca9cb..12c0f379f 160000 --- a/extern/toml11 +++ b/extern/toml11 @@ -1 +1 @@ -Subproject commit 85faca9cbe8d76324ff38c1801be44c63e12d5be +Subproject commit 12c0f379f2e865b4ce984758d5ae004f9de07d69 diff --git a/input.example.toml b/input.example.toml index 2423ffff0..88589495c 100644 --- a/input.example.toml +++ b/input.example.toml @@ -336,6 +336,10 @@ # @note: For T, in cartesian can also use "x" "y" "z" instead of "1" "2" "3" # @note: By default, we accumulate moments from all massive species, one can specify only specific species: e.g., Ttt_1_2, Rho_1, Rho_3_4 quantities = "" + # Custom (user-defined) field quantities: + # @type: array of strings + # @default: [] + custom = "" # @NOT_IMPLEMENTED: Stride for the output of fields: # @type: unsigned short: > 1 # @default: 1 diff --git a/setups/wip/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp similarity index 65% rename from setups/wip/shock/pgen.hpp rename to setups/srpic/shock/pgen.hpp index 4b1cd392e..f07b99878 100644 --- a/setups/wip/shock/pgen.hpp +++ b/setups/srpic/shock/pgen.hpp @@ -14,23 +14,6 @@ namespace user { using namespace ntt; - template - struct DriftDist : public arch::EnergyDistribution { - DriftDist(const M& metric, real_t ux) - : arch::EnergyDistribution { metric } - , ux { ux } {} - - Inline void operator()(const coord_t&, - vec_t& v, - unsigned short) const override { - v[0] = -ux; - v[1] = 0.1 * ux; - } - - private: - const real_t ux; - }; - template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator @@ -68,13 +51,6 @@ namespace user { local_domain, injector, 1.0); - // const auto energy_dist = DriftDist(local_domain.mesh.metric, drift_ux); - // const auto injector = arch::UniformInjector(energy_dist, - // { 1, 2 }); - // arch::InjectUniform>(params, - // local_domain, - // injector, - // 1.0); } }; diff --git a/setups/srpic/shock/shock.py b/setups/srpic/shock/shock.py new file mode 100644 index 000000000..64224c728 --- /dev/null +++ b/setups/srpic/shock/shock.py @@ -0,0 +1,75 @@ +import nt2.read as nt2r +import matplotlib.pyplot as plt +import matplotlib as mpl + +data = nt2r.Data("shock-03.h5") + + +def frame(ti, f): + quantities = [ + { + "name": "density", + "compute": lambda f: f.N_2 + f.N_1, + "cmap": "inferno", + "norm": mpl.colors.Normalize(0, 5), + }, + { + "name": r"$E_x$", + "compute": lambda f: f.Ex, + "cmap": "RdBu_r", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + { + "name": r"$E_y$", + "compute": lambda f: f.Ey, + "cmap": "RdBu_r", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + { + "name": r"$E_z$", + "compute": lambda f: f.Ez, + "cmap": "RdBu_r", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + { + "name": r"$B_x$", + "compute": lambda f: f.Bx, + "cmap": "BrBG", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + { + "name": r"$B_y$", + "compute": lambda f: f.By, + "cmap": "BrBG", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + { + "name": r"$B_z$", + "compute": lambda f: f.Bz, + "cmap": "BrBG", + "norm": mpl.colors.Normalize(-0.05, 0.05), + }, + ] + fig = plt.figure(figsize=(12, 5.5), dpi=300) + gs = fig.add_gridspec(len(quantities), 1, hspace=0.02) + axs = [fig.add_subplot(gs[i]) for i in range(len(quantities))] + + for ax, q in zip(axs, quantities): + q["compute"](f).coarsen(x=2, y=2).mean().plot( + ax=ax, + cmap=q["cmap"], + norm=q["norm"], + cbar_kwargs={"label": q["name"], "shrink": 0.8, "aspect": 10, "pad": 0.005}, + ) + for i, ax in enumerate(axs): + ax.set(aspect=1) + if i != 0: + ax.set(title=None) + if i != len(axs) - 1: + ax.set( + xticks=[], + xticklabels=[], + xlabel=None, + title=ax.get_title().split(",")[0], + ) + return fig diff --git a/setups/wip/shock/shock.toml b/setups/srpic/shock/shock.toml similarity index 60% rename from setups/wip/shock/shock.toml rename to setups/srpic/shock/shock.toml index c67ee181a..f48edb2d6 100644 --- a/setups/wip/shock/shock.toml +++ b/setups/srpic/shock/shock.toml @@ -1,7 +1,7 @@ [simulation] name = "shock" engine = "srpic" - runtime = 100.0 + runtime = 50.0 [grid] resolution = [2048, 128] @@ -16,34 +16,36 @@ [scales] larmor0 = 1e-2 - skindepth0 = 1e-3 + skindepth0 = 1e-2 [algorithms] - current_filters = 4 + current_filters = 8 [algorithms.timestep] CFL = 0.5 [particles] - ppc0 = 4.0 + ppc0 = 16.0 [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e7 + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e8 [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e7 + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e8 [setup] - drift_ux = -0.1 - temperature = 0.01 + drift_ux = 0.1 + temperature = 1e-3 [output] - fields = ["N_1", "N_2", "E", "B", "J"] - format = "hdf5" interval_time = 0.1 + format = "hdf5" + + [output.fields] + quantities = ["N_1", "N_2", "E", "B", "T0i_1", "T0i_2", "J"] diff --git a/setups/tests/customout/customout.toml b/setups/tests/customout/customout.toml new file mode 100644 index 000000000..497b96dc2 --- /dev/null +++ b/setups/tests/customout/customout.toml @@ -0,0 +1,50 @@ +[simulation] + name = "customout" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [256, 256] + extent = [[-1.0, 1.0], [-1.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.01 + skindepth0 = 0.01 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 20.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e7 + pusher = "Boris" + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e7 + pusher = "Boris" + +[output] + format = "hdf5" + interval_time = 0.02 + + [output.fields] + quantities = ["E", "B", "J"] + custom = ["mybuff", "EdotB+1"] diff --git a/setups/tests/customout/pgen.hpp b/setups/tests/customout/pgen.hpp new file mode 100644 index 000000000..22c8f6564 --- /dev/null +++ b/setups/tests/customout/pgen.hpp @@ -0,0 +1,86 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "arch/traits.h" + +#include "archetypes/problem_generator.h" +#include "framework/domain/metadomain.h" + +namespace user { + using namespace ntt; + + template + struct PGen : public arch::ProblemGenerator { + // compatibility traits for the problem generator + static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto metrics { traits::compatible_with::value }; + static constexpr auto dimensions { traits::compatible_with::value }; + + // for easy access to variables in the child class + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + array_t cbuff; + + inline PGen(const SimulationParams& p, const Metadomain&) + : arch::ProblemGenerator(p) {} + + inline PGen() {} + + void CustomPostStep(std::size_t step, long double, Domain& domain) { + if (step == 0) { + // allocate the array at time = 0 + cbuff = array_t("cbuff", + domain.mesh.n_all(in::x1), + domain.mesh.n_all(in::x2)); + } + // fill with zeros + Kokkos::deep_copy(cbuff, ZERO); + // populate the array atomically (here it's not strictly necessary) + auto cbuff_sc = Kokkos::Experimental::create_scatter_view(cbuff); + Kokkos::parallel_for( + "FillCbuff", + domain.mesh.rangeActiveCells(), + Lambda(index_t i1, index_t i2) { + auto cbuff_acc = cbuff_sc.access(); + cbuff_acc(i1, i2) += static_cast(i1 + i2); + }); + Kokkos::Experimental::contribute(cbuff, cbuff_sc); + } + + void CustomFieldOutput(const std::string& name, + ndfield_t buffer, + std::size_t index, + const Domain& domain) { + printf("CustomFieldOutput: %s\n", name.c_str()); + // examples for 2D + if (name == "mybuff") { + // copy the custom buffer to the buffer output + Kokkos::deep_copy(Kokkos::subview(buffer, Kokkos::ALL, Kokkos::ALL, index), + cbuff); + } else if (name == "EdotB+1") { + // calculate the custom buffer from EM fields + const auto& EM = domain.fields.em; + Kokkos::parallel_for( + "EdotB+1", + domain.mesh.rangeActiveCells(), + Lambda(index_t i1, index_t i2) { + buffer(i1, i2, index) = EM(i1, i2, em::ex1) * EM(i1, i2, em::bx1) + + EM(i1, i2, em::ex2) * EM(i1, i2, em::bx2) + + EM(i1, i2, em::ex3) * EM(i1, i2, em::bx3) + + ONE; + }); + } else { + raise::Error("Custom output not provided", HERE); + } + } + }; + +} // namespace user + +#endif diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index dd1aa957c..4e75003bb 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -1,13 +1,17 @@ /** * @file archetypes/particle_injector.h - * @brief Particle injector routines - * ... + * @brief Particle injector routines and classes + * @implements + * - arch::UniformInjector<> + * - arch::NonUniformInjector<> + * - arch::AtmosphereInjector<> + * - arch::InjectUniform<> -> void + * - arch::InjectGlobally<> -> void + * - arch::InjectNonUniform<> -> void + * @namespaces: + * - arch:: */ -/* -------------------------------------------------------------------------- */ -/* This header file is still under construction */ -/* -------------------------------------------------------------------------- */ - #ifndef ARCHETYPES_PARTICLE_INJECTOR_H #define ARCHETYPES_PARTICLE_INJECTOR_H diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp index d94d00376..90dec3326 100644 --- a/src/engines/engine_printer.cpp +++ b/src/engines/engine_printer.cpp @@ -137,6 +137,7 @@ namespace ntt { const auto hash = std::string(ENTITY_GIT_HASH); const auto pgen = std::string(PGEN); const auto nspec = metadomain.species_params().size(); + const auto precision = (sizeof(real_t) == 4) ? "single" : "double"; #if defined(__clang__) const std::string ccx = "Clang/LLVM " __clang_version__; @@ -219,6 +220,7 @@ namespace ntt { add_param(report, 4, "HDF5", "%s", hdf5_version.c_str()); add_param(report, 4, "Kokkos", "%s", kokkos_version.c_str()); add_param(report, 4, "ADIOS2", "%s", adios2_version.c_str()); + add_param(report, 4, "Precision", "%s", precision); add_param(report, 4, "Debug", "%s", dbg.c_str()); report += "\n"; add_category(report, 4, "Configuration"); diff --git a/src/engines/engine_run.cpp b/src/engines/engine_run.cpp index f69b4301e..4485e0e40 100644 --- a/src/engines/engine_run.cpp +++ b/src/engines/engine_run.cpp @@ -9,6 +9,8 @@ #include "metrics/qspherical.h" #include "metrics/spherical.h" +#include "framework/domain/domain.h" + #include "engines/engine.hpp" namespace ntt { @@ -61,7 +63,21 @@ namespace ntt { auto print_output = false; #if defined(OUTPUT_ENABLED) timers.start("Output"); - print_output = m_metadomain.Write(m_params, step, time); + if constexpr ( + traits::has_method::value) { + auto lambda_custom_field_output = [&](const std::string& name, + ndfield_t& buff, + std::size_t idx, + const Domain& dom) { + m_pgen.CustomFieldOutput(name, buff, idx, dom); + }; + print_output = m_metadomain.Write(m_params, + step, + time, + lambda_custom_field_output); + } else { + print_output = m_metadomain.Write(m_params, step, time); + } timers.stop("Output"); #endif diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index ea01f591b..bddc557c9 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -208,6 +208,7 @@ namespace ntt { m_params.template get( "algorithms.timestep.correction") * dt; + auto range = range_with_axis_BCs(domain); if constexpr (M::CoordType == Coord::Cart) { // minkowski case const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); @@ -222,15 +223,9 @@ namespace ntt { Kokkos::parallel_for( "Ampere", - domain.mesh.rangeActiveCells(), + range, kernel::mink::Ampere_kernel(domain.fields.em, coeff1, coeff2)); } else { - range_t range {}; - if constexpr (M::Dim == Dim::_2D) { - range = CreateRangePolicy( - { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, - { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); - } const auto ni2 = domain.mesh.n_active(in::x2); Kokkos::parallel_for("Ampere", range, @@ -538,13 +533,8 @@ namespace ntt { coeff / V0, ONE / n0)); } else { - range_t range {}; - if constexpr (M::Dim == Dim::_2D) { - range = CreateRangePolicy( - { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, - { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); - } - const auto ni2 = domain.mesh.n_active(in::x2); + auto range = range_with_axis_BCs(domain); + const auto ni2 = domain.mesh.n_active(in::x2); Kokkos::parallel_for( "Ampere", range, @@ -560,34 +550,7 @@ namespace ntt { void CurrentsFilter(domain_t& domain) { logger::Checkpoint("Launching currents filtering kernels", HERE); - range_t range = domain.mesh.rangeActiveCells(); - if constexpr (M::CoordType != Coord::Cart) { - /** - * @brief taking one extra cell in the x2 direction if AXIS BCs - * . . . . . - * . ^= =^ . - * . |* *\*. - * . |* *\*. - * . ^- -^ . - * . . . . . - */ - if constexpr (M::Dim == Dim::_2D) { - if (domain.mesh.flds_bc_in({ +1, 0 }) == FldsBC::AXIS) { - range = CreateRangePolicy( - { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, - { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); - } - } else if constexpr (M::Dim == Dim::_3D) { - if (domain.mesh.flds_bc_in({ +1, 0, 0 }) == FldsBC::AXIS) { - range = CreateRangePolicy({ domain.mesh.i_min(in::x1), - domain.mesh.i_min(in::x2), - domain.mesh.i_min(in::x3) }, - { domain.mesh.i_max(in::x1), - domain.mesh.i_max(in::x2) + 1, - domain.mesh.i_max(in::x3) }); - } - } - } + auto range = range_with_axis_BCs(domain); const auto nfilter = m_params.template get( "algorithms.current_filters"); tuple_t size; @@ -1228,6 +1191,32 @@ namespace ntt { } return { sign, dim, xg_min, xg_max }; } + + auto range_with_axis_BCs(const domain_t& domain) -> range_t { + auto range = domain.mesh.rangeActiveCells(); + if constexpr (M::CoordType != Coord::Cart) { + /** + * @brief taking one extra cell in the x2 direction if AXIS BCs + */ + if constexpr (M::Dim == Dim::_2D) { + if (domain.mesh.flds_bc_in({ 0, +1 }) == FldsBC::AXIS) { + range = CreateRangePolicy( + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + } + } else if constexpr (M::Dim == Dim::_3D) { + if (domain.mesh.flds_bc_in({ 0, +1, 0 }) == FldsBC::AXIS) { + range = CreateRangePolicy({ domain.mesh.i_min(in::x1), + domain.mesh.i_min(in::x2), + domain.mesh.i_min(in::x3) }, + { domain.mesh.i_max(in::x1), + domain.mesh.i_max(in::x2) + 1, + domain.mesh.i_max(in::x3) }); + } + } + } + return range; + } }; } // namespace ntt diff --git a/src/framework/domain/comm_mpi.hpp b/src/framework/domain/comm_mpi.hpp index 2854599ef..63dd8271a 100644 --- a/src/framework/domain/comm_mpi.hpp +++ b/src/framework/domain/comm_mpi.hpp @@ -55,9 +55,6 @@ namespace comm { if ((send_idx == idx) and (recv_idx == idx)) { // trivial copy if sending to self and receiving from self - // raise::ErrorIf((recv_idx != idx) || (send_idx != idx), - // "Cannot send to self and receive from another domain", - // HERE); if (not additive) { // simply filling the ghost cells if constexpr (D == Dim::_1D) { @@ -446,214 +443,4 @@ namespace comm { } // namespace comm -/** - * attempt to receive directly into fields - */ -// template -// auto to_tuple_impl(const Array& a, std::index_sequence) { -// return std::make_tuple(a[I]...); -// } - -// template > auto to_tuple(const std::array& a) { -// return to_tuple_impl(a, Indices {}); -// } - -// template -// inline void CommunicateField(unsigned int idx, -// ndfield_t& fld, -// unsigned int send_idx, -// unsigned int recv_idx, -// int send_rank, -// int recv_rank, -// const std::vector& send_slice, -// const std::vector& recv_slice, -// const range_tuple_t& comps, -// bool additive) { -// raise::ErrorIf(send_rank < 0 && recv_rank < 0, -// "CommunicateField called with negative ranks", -// HERE); - -// int rank; -// MPI_Comm_rank(MPI_COMM_WORLD, &rank); - -// raise::ErrorIf( -// (send_rank == rank && send_idx != idx) || -// (recv_rank == rank && recv_idx != idx), -// "Multiple-domain single-rank communication not yet implemented", -// HERE); - -// // trivial copy if sending to self and receiving from self -// if ((send_idx == idx) || (recv_idx == idx)) { -// raise::ErrorIf((recv_idx != idx) || (send_idx != idx), -// "Cannot send to self and receive from another domain", -// HERE); -// // sending/recv to/from self -// if constexpr (D == Dim::_1D) { -// Kokkos::deep_copy(Kokkos::subview(fld, recv_slice[0], comps), -// Kokkos::subview(fld, send_slice[0], comps)); -// } else if constexpr (D == Dim::_2D) { -// Kokkos::deep_copy( -// Kokkos::subview(fld, recv_slice[0], recv_slice[1], comps), -// Kokkos::subview(fld, send_slice[0], send_slice[1], comps)); -// } else if constexpr (D == Dim::_3D) { -// Kokkos::deep_copy( -// Kokkos::subview(fld, recv_slice[0], recv_slice[1], recv_slice[2], comps), -// Kokkos::subview(fld, send_slice[0], send_slice[1], send_slice[2], comps)); -// } -// return; -// } - -// std::size_t nsend { comps.second - comps.first }, -// nrecv { comps.second - comps.first }; -// // ndarray_t(D) + 1> send_fld, recv_fld; - -// for (short d { 0 }; d < (short)D; ++d) { -// if (send_rank >= 0) { -// nsend *= (send_slice[d].second - send_slice[d].first); -// } -// if (recv_rank >= 0) { -// nrecv *= (recv_slice[d].second - recv_slice[d].first); -// } -// } - -// std::vector> send_init_list; -// std::vector> recv_init_list; - -// if (send_rank >= 0) { -// for (unsigned short d = 0; d < D; ++d) { -// send_init_list.push_back({ send_slice[d].first, send_slice[d].second }); -// recv_init_list.push_back({ recv_slice[d].first, recv_slice[d].second }); -// } -// send_init_list.push_back({ comps.first, comps.second }); -// recv_init_list.push_back({ comps.first, comps.second }); -// } - -// std::array, D + 1> send_array; -// std::copy(send_init_list.begin(), send_init_list.end(), send_array.begin()); - -// std::array, D + 1> recv_array; -// std::copy(recv_init_list.begin(), recv_init_list.end(), recv_array.begin()); - -// if (send_rank >= 0 && recv_rank >= 0) { -// auto send_subview = std::apply( -// [&](auto... args) { -// return Kokkos::subview(fld, args...); -// }, -// to_tuple(send_array)); -// auto recv_subview = std::apply( -// [&](auto... args) { -// return Kokkos::subview(fld, args...); -// }, -// to_tuple(recv_array)); - -// MPI_Sendrecv(send_subview.data(), -// nsend, -// mpi::get_type(), -// send_rank, -// 0, -// recv_subview.data(), -// nrecv, -// mpi::get_type(), -// recv_rank, -// 0, -// MPI_COMM_WORLD, -// MPI_STATUS_IGNORE); -// } else if (send_rank >= 0) { -// auto send_subview = std::apply( -// [&](auto... args) { -// return Kokkos::subview(fld, args...); -// }, -// to_tuple(send_array)); - -// MPI_Send(send_subview.data(), -// nsend, -// mpi::get_type(), -// send_rank, -// 0, -// MPI_COMM_WORLD); -// } else if (recv_rank >= 0) { -// auto recv_subview = std::apply( -// [&](auto... args) { -// return Kokkos::subview(fld, args...); -// }, -// to_tuple(recv_array)); - -// MPI_Recv(recv_subview.data(), -// nrecv, -// mpi::get_type(), -// recv_rank, -// 0, -// MPI_COMM_WORLD, -// MPI_STATUS_IGNORE); -// } else { -// raise::Error("CommunicateField called with negative ranks", HERE); -// } -// if (recv_rank >= 0) { -// // // !TODO: perhaps directly recv to the fld? -// // if (not additive) { -// // if constexpr (D == Dim::_1D) { -// // Kokkos::deep_copy(Kokkos::subview(fld, recv_slice[0], comps), recv_fld); -// // } else if constexpr (D == Dim::_2D) { -// // Kokkos::deep_copy( -// // Kokkos::subview(fld, recv_slice[0], recv_slice[1], comps), -// // recv_fld); -// // } else if constexpr (D == Dim::_3D) { -// // Kokkos::deep_copy( -// // Kokkos::subview(fld, recv_slice[0], recv_slice[1], recv_slice[2], comps), -// // recv_fld); -// // } -// // } else { -// // if constexpr (D == Dim::_1D) { -// // const auto offset_x1 = recv_slice[0].first; -// // const auto offset_c = comps.first; -// // Kokkos::parallel_for( -// // "CommunicateField-extract", -// // Kokkos::MDRangePolicy, AccelExeSpace>( -// // { recv_slice[0].first, comps.first }, -// // { recv_slice[0].second, comps.second }), -// // Lambda(index_t i1, index_t ci) { -// // fld(i1, ci) += recv_fld(i1 - offset_x1, ci - offset_c); -// // }); -// // } else if constexpr (D == Dim::_2D) { -// // const auto offset_x1 = recv_slice[0].first; -// // const auto offset_x2 = recv_slice[1].first; -// // const auto offset_c = comps.first; -// // Kokkos::parallel_for( -// // "CommunicateField-extract", -// // Kokkos::MDRangePolicy, AccelExeSpace>( -// // { recv_slice[0].first, recv_slice[1].first, comps.first }, -// // { recv_slice[0].second, recv_slice[1].second, comps.second }), -// // Lambda(index_t i1, index_t i2, index_t ci) { -// // fld(i1, i2, ci) += recv_fld(i1 - offset_x1, -// // i2 - offset_x2, -// // ci - offset_c); -// // }); -// // } else if constexpr (D == Dim::_3D) { -// // const auto offset_x1 = recv_slice[0].first; -// // const auto offset_x2 = recv_slice[1].first; -// // const auto offset_x3 = recv_slice[2].first; -// // const auto offset_c = comps.first; -// // Kokkos::parallel_for( -// // "CommunicateField-extract", -// // Kokkos::MDRangePolicy, AccelExeSpace>( -// // { recv_slice[0].first, -// // recv_slice[1].first, -// // recv_slice[2].first, -// // comps.first }, -// // { recv_slice[0].second, -// // recv_slice[1].second, -// // recv_slice[2].second, -// // comps.second }), -// // Lambda(index_t i1, index_t i2, index_t i3, index_t ci) { -// // fld(i1, i2, i3, ci) += recv_fld(i1 - offset_x1, -// // i2 - offset_x2, -// // i3 - offset_x3, -// // ci - offset_c); -// // }); -// // } -// // } -// } -// } - #endif // FRAMEWORK_DOMAIN_COMM_MPI_HPP diff --git a/src/framework/domain/grid.cpp b/src/framework/domain/grid.cpp index d189eaf5d..9302386e1 100644 --- a/src/framework/domain/grid.cpp +++ b/src/framework/domain/grid.cpp @@ -9,19 +9,19 @@ namespace ntt { template <> - auto Grid::rangeAllCells() -> range_t { + auto Grid::rangeAllCells() const -> range_t { box_region_t region { CellLayer::allLayer }; return rangeCells(region); } template <> - auto Grid::rangeAllCells() -> range_t { + auto Grid::rangeAllCells() const -> range_t { box_region_t region { CellLayer::allLayer, CellLayer::allLayer }; return rangeCells(region); } template <> - auto Grid::rangeAllCells() -> range_t { + auto Grid::rangeAllCells() const -> range_t { box_region_t region { CellLayer::allLayer, CellLayer::allLayer, CellLayer::allLayer }; @@ -29,19 +29,19 @@ namespace ntt { } template <> - auto Grid::rangeActiveCells() -> range_t { + auto Grid::rangeActiveCells() const -> range_t { box_region_t region { CellLayer::activeLayer }; return rangeCells(region); } template <> - auto Grid::rangeActiveCells() -> range_t { + auto Grid::rangeActiveCells() const -> range_t { box_region_t region { CellLayer::activeLayer, CellLayer::activeLayer }; return rangeCells(region); } template <> - auto Grid::rangeActiveCells() -> range_t { + auto Grid::rangeActiveCells() const -> range_t { box_region_t region { CellLayer::activeLayer, CellLayer::activeLayer, CellLayer::activeLayer }; @@ -49,7 +49,7 @@ namespace ntt { } template - auto Grid::rangeCells(const box_region_t& region) -> range_t { + auto Grid::rangeCells(const box_region_t& region) const -> range_t { tuple_t imin, imax; for (unsigned short i = 0; i < (unsigned short)D; i++) { switch (region[i]) { @@ -87,7 +87,8 @@ namespace ntt { // !TODO: too ugly, implement a better solution (combine with device) template - auto Grid::rangeCellsOnHost(const box_region_t& region) -> range_h_t { + auto Grid::rangeCellsOnHost(const box_region_t& region) const + -> range_h_t { tuple_t imin, imax; for (unsigned short i = 0; i < (unsigned short)D; i++) { switch (region[i]) { @@ -123,19 +124,19 @@ namespace ntt { } template <> - auto Grid::rangeAllCellsOnHost() -> range_h_t { + auto Grid::rangeAllCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::allLayer }; return rangeCellsOnHost(region); } template <> - auto Grid::rangeAllCellsOnHost() -> range_h_t { + auto Grid::rangeAllCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::allLayer, CellLayer::allLayer }; return rangeCellsOnHost(region); } template <> - auto Grid::rangeAllCellsOnHost() -> range_h_t { + auto Grid::rangeAllCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::allLayer, CellLayer::allLayer, CellLayer::allLayer }; @@ -143,19 +144,19 @@ namespace ntt { } template <> - auto Grid::rangeActiveCellsOnHost() -> range_h_t { + auto Grid::rangeActiveCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::activeLayer }; return rangeCellsOnHost(region); } template <> - auto Grid::rangeActiveCellsOnHost() -> range_h_t { + auto Grid::rangeActiveCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::activeLayer, CellLayer::activeLayer }; return rangeCellsOnHost(region); } template <> - auto Grid::rangeActiveCellsOnHost() -> range_h_t { + auto Grid::rangeActiveCellsOnHost() const -> range_h_t { box_region_t region { CellLayer::activeLayer, CellLayer::activeLayer, CellLayer::activeLayer }; @@ -163,7 +164,7 @@ namespace ntt { } template - auto Grid::rangeCells(const tuple_t, D>& ranges) + auto Grid::rangeCells(const tuple_t, D>& ranges) const -> range_t { tuple_t imin, imax; for (unsigned short i = 0; i < (unsigned short)D; i++) { @@ -182,4 +183,4 @@ namespace ntt { template struct Grid; template struct Grid; -} // namespace ntt \ No newline at end of file +} // namespace ntt diff --git a/src/framework/domain/grid.h b/src/framework/domain/grid.h index 176263189..97a939117 100644 --- a/src/framework/domain/grid.h +++ b/src/framework/domain/grid.h @@ -159,19 +159,19 @@ namespace ntt { * @brief Loop over all active cells (disregard ghost cells) * @returns Kokkos range policy with proper min/max indices and dimension */ - auto rangeActiveCells() -> range_t; + auto rangeActiveCells() const -> range_t; /** * @brief Loop over all cells * @returns Kokkos range policy with proper min/max indices and dimension */ - auto rangeAllCells() -> range_t; + auto rangeAllCells() const -> range_t; /** * @brief Pick a particular region of cells * @param box_region_t region of cells to pick: tuple of cellLayer objects * @returns Kokkos range policy with proper min/max indices and dimension */ - auto rangeCells(const box_region_t&) -> range_t; + auto rangeCells(const box_region_t&) const -> range_t; /** * @brief Pick a particular region of cells * @overload @@ -180,7 +180,7 @@ namespace ntt { * @example {{0, 0}, {0, 0}, {0, 0}} corresponds to allActiveLayer in all 3 dimensions * @returns Kokkos range policy with proper min/max indices and dimension */ - auto rangeCells(const tuple_t, D>&) -> range_t; + auto rangeCells(const tuple_t, D>&) const -> range_t; /* Ranges in the host execution space ----------------------------------- */ /** @@ -188,20 +188,20 @@ namespace ntt { * @returns Kokkos range policy in the host space with proper min/max * indices and dimension. */ - auto rangeActiveCellsOnHost() -> range_h_t; + auto rangeActiveCellsOnHost() const -> range_h_t; /** * @brief Loop over all cells * @returns Kokkos range policy in the host space with proper min/max * indices and dimension. */ - auto rangeAllCellsOnHost() -> range_h_t; + auto rangeAllCellsOnHost() const -> range_h_t; /** * @brief Pick a particular region of cells * @param box_region_t region of cells to pick: tuple of cellLayer objects * @returns Kokkos range policy in the host space with proper min/max * indices and dimension. */ - auto rangeCellsOnHost(const box_region_t&) -> range_h_t; + auto rangeCellsOnHost(const box_region_t&) const -> range_h_t; protected: std::vector m_resolution; diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index 74e5da2ce..fb81fcfca 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -18,6 +18,7 @@ #include "enums.h" #include "global.h" +#include "arch/kokkos_aliases.h" #include "utils/timer.h" #include "framework/containers/species.h" @@ -33,6 +34,7 @@ #include "output/writer.h" #endif +#include #include #include #include @@ -111,7 +113,13 @@ namespace ntt { #if defined(OUTPUT_ENABLED) void InitWriter(const SimulationParams&); - auto Write(const SimulationParams&, std::size_t, long double) -> bool; + auto Write(const SimulationParams&, + std::size_t, + long double, + std::function&, + std::size_t, + const Domain&)> = {}) -> bool; #endif Metadomain(const Metadomain&) = delete; diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 77c41ee25..cabb37093 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -14,6 +14,7 @@ #include "metrics/spherical.h" #include "framework/containers/particles.h" +#include "framework/domain/domain.h" #include "framework/domain/metadomain.h" #include "framework/parameters.h" @@ -25,6 +26,8 @@ #include #include +#include +#include #include namespace ntt { @@ -61,9 +64,17 @@ namespace ntt { M::CoordType); const auto fields_to_write = params.template get>( "output.fields.quantities"); + const auto custom_fields_to_write = params.template get>( + "output.fields.custom"); + std::vector all_fields_to_write; + std::merge(fields_to_write.begin(), + fields_to_write.end(), + custom_fields_to_write.begin(), + custom_fields_to_write.end(), + std::back_inserter(all_fields_to_write)); const auto species_to_write = params.template get>( "output.particles.species"); - g_writer.defineFieldOutputs(S, fields_to_write); + g_writer.defineFieldOutputs(S, all_fields_to_write); g_writer.defineParticleOutputs(M::PrtlDim, species_to_write); // spectra write all particle species std::vector spectra_species {}; @@ -153,9 +164,13 @@ namespace ntt { } template - auto Metadomain::Write(const SimulationParams& params, - std::size_t step, - long double time) -> bool { + auto Metadomain::Write( + const SimulationParams& params, + std::size_t step, + long double time, + std::function< + void(const std::string&, ndfield_t&, std::size_t, const Domain&)> + CustomFieldOutput) -> bool { raise::ErrorIf( local_subdomain_indices().size() != 1, "Output for now is only supported for one subdomain per rank", @@ -278,14 +293,24 @@ namespace ntt { } else { raise::Error("Wrong moment requested for output", HERE); } - SynchronizeFields(*local_domain, - Comm::Bckp, - { addresses.back(), addresses.back() + 1 }); + } else if (fld.is_custom()) { + if (CustomFieldOutput) { + CustomFieldOutput(fld.name().substr(1), + local_domain->fields.bckp, + addresses.back(), + *local_domain); + } else { + raise::Error("Custom output requested but no function provided", + HERE); + } } else { - raise::Error( - "Wrong # of components requested for non-moment output", - HERE); + raise::Error("Wrong # of components requested for " + "non-moment/non-custom output", + HERE); } + SynchronizeFields(*local_domain, + Comm::Bckp, + { addresses.back(), addresses.back() + 1 }); } else if (fld.comp.size() == 3) { // vector for (auto i = 0; i < 3; ++i) { names.push_back(fld.name(i)); diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 19f9a4683..6e581a6b7 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -383,15 +383,21 @@ namespace ntt { promiseToDefine("output.spectra.interval_time"); promiseToDefine("output.spectra.enable"); - const auto flds_out = toml::find_or(raw_data, + const auto flds_out = toml::find_or(raw_data, "output", "fields", "quantities", std::vector {}); + const auto custom_flds_out = toml::find_or(raw_data, + "output", + "fields", + "custom", + std::vector {}); if (flds_out.size() == 0) { raise::Warning("No fields output specified", HERE); } set("output.fields.quantities", flds_out); + set("output.fields.custom", custom_flds_out); set("output.fields.mom_smooth", toml::find_or(raw_data, "output", diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h index ea93a8571..e915bdf1a 100644 --- a/src/global/arch/traits.h +++ b/src/global/arch/traits.h @@ -13,6 +13,7 @@ * - traits::pgen::field_driver_t * - traits::pgen::init_prtls_t * - traits::pgen::custom_fields_t + * - traits::pgen::custom_field_output_t * - traits::pgen::custom_poststep_t * - traits::check_compatibility<> * - traits::compatibility<> @@ -100,6 +101,9 @@ namespace traits { template using custom_poststep_t = decltype(&T::CustomPostStep); + + template + using custom_field_output_t = decltype(&T::CustomFieldOutput); } // namespace pgen // for pgen extforce diff --git a/src/global/enums.h b/src/global/enums.h index 2dc492dac..57822dec4 100644 --- a/src/global/enums.h +++ b/src/global/enums.h @@ -13,7 +13,7 @@ * - enum ntt::PrtlPusher // boris, vay, photon, none * - enum ntt::Cooling // synchrotron, none * - enum ntt::FldsID // e, dive, d, divd, b, h, j, - * a, t, rho, charge, n, nppc + * a, t, rho, charge, n, nppc, custom * @namespaces: * - ntt:: * @note Enums of the same type can be compared with each other and with strings @@ -288,16 +288,17 @@ namespace ntt { Charge = 11, N = 12, Nppc = 13, + Custom = 14, }; constexpr FldsID(uint8_t c) : enums_hidden::BaseEnum { c } {} - static constexpr type variants[] = { E, divE, D, divD, B, H, J, - A, T, Rho, Charge, N, Nppc }; - static constexpr const char* lookup[] = { "e", "dive", "d", "divd", - "b", "h", "j", "a", - "t", "rho", "charge", "n", - "nppc" }; + static constexpr type variants[] = { E, divE, D, divD, B, H, J, + A, T, Rho, Charge, N, Nppc, Custom }; + static constexpr const char* lookup[] = { "e", "dive", "d", "divd", + "b", "h", "j", "a", + "t", "rho", "charge", "n", + "nppc", "custom" }; static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]); }; diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp index 81e63d57a..1fc57398f 100644 --- a/src/global/tests/enums.cpp +++ b/src/global/tests/enums.cpp @@ -66,8 +66,9 @@ auto main() -> int { enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" }; enum_str_t all_coolings = { "synchrotron", "none" }; - enum_str_t all_out_flds = { "e", "dive", "d", "divd", "b", "h", "j", - "a", "t", "rho", "charge", "n", "nppc" }; + enum_str_t all_out_flds = { "e", "dive", "d", "divd", "b", + "h", "j", "a", "t", "rho", + "charge", "n", "nppc", "custom" }; checkEnum(all_coords); checkEnum(all_metrics); diff --git a/src/kernels/ampere_sr.hpp b/src/kernels/ampere_sr.hpp index 6dde42ff0..e4faec6ce 100644 --- a/src/kernels/ampere_sr.hpp +++ b/src/kernels/ampere_sr.hpp @@ -122,7 +122,8 @@ namespace kernel::sr { */ template class CurrentsAmpere_kernel { - static constexpr auto D = M::Dim; + static constexpr auto D = M::Dim; + static constexpr std::size_t i2min = N_GHOSTS; ndfield_t E; ndfield_t J; @@ -161,9 +162,8 @@ namespace kernel::sr { Inline void operator()(index_t i1, index_t i2) const { if constexpr (D == Dim::_2D) { - constexpr std::size_t i2min { N_GHOSTS }; - const real_t i1_ { COORD(i1) }; - const real_t i2_ { COORD(i2) }; + const real_t i1_ { COORD(i1) }; + const real_t i2_ { COORD(i2) }; // convert the "curly" current, to contravariant, normalized to // `J0=n0*q0` then add "curly" current with the right coefficient diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index ae8440342..9d3fd7d81 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -512,37 +512,37 @@ namespace kernel { if (ppc == 0) { return; } - const auto index = Kokkos::atomic_fetch_add(&idx(), ppc); - auto rand_gen = random_pool.get_state(); + auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { - const auto dx1 = Random(rand_gen); + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); - i1s_1(index + offset1 - p) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1 - p) = dx1; - i1s_2(index + offset2 - p) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2 - p) = dx1; + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; vec_t v_T { ZERO }, v_XYZ { ZERO }; energy_dist(x_Ph, v_T, spidx1); metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_1(index + offset1 - p) = v_XYZ[0]; - ux2s_1(index + offset1 - p) = v_XYZ[1]; - ux3s_1(index + offset1 - p) = v_XYZ[2]; + ux1s_1(index + offset1) = v_XYZ[0]; + ux2s_1(index + offset1) = v_XYZ[1]; + ux3s_1(index + offset1) = v_XYZ[2]; energy_dist(x_Ph, v_T, spidx2); metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_2(index + offset2 - p) = v_XYZ[0]; - ux2s_2(index + offset2 - p) = v_XYZ[1]; - ux3s_2(index + offset2 - p) = v_XYZ[2]; + ux1s_2(index + offset2) = v_XYZ[0]; + ux2s_2(index + offset2) = v_XYZ[1]; + ux3s_2(index + offset2) = v_XYZ[2]; - tags_1(index + offset1 - p) = ParticleTag::alive; - tags_2(index + offset2 - p) = ParticleTag::alive; + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; if (M::CoordType == Coord::Cart) { - weights_1(index + offset1 - p) = ONE; - weights_2(index + offset2 - p) = ONE; + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; } else { - const auto wei = metric.sqrt_det_h({ i1_ + HALF }); - weights_1(index + offset1 - p) = wei; - weights_2(index + offset2 - p) = wei; + const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; } } random_pool.free_state(rand_gen); @@ -568,21 +568,21 @@ namespace kernel { if (ppc == 0) { return; } - const auto index = Kokkos::atomic_fetch_add(&idx(), ppc); - auto rand_gen = random_pool.get_state(); + auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); - i1s_1(index + offset1 - p) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1 - p) = dx1; - i1s_2(index + offset2 - p) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2 - p) = dx1; + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; - i2s_1(index + offset1 - p) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1 - p) = dx2; - i2s_2(index + offset2 - p) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2 - p) = dx2; + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; vec_t v_T { ZERO }, v_Cd { ZERO }; energy_dist(x_Ph, v_T, spidx1); @@ -591,28 +591,28 @@ namespace kernel { } else if constexpr (S == SimEngine::GRPIC) { metric.template transform(x_Cd_, v_T, v_Cd); } - ux1s_1(index + offset1 - p) = v_Cd[0]; - ux2s_1(index + offset1 - p) = v_Cd[1]; - ux3s_1(index + offset1 - p) = v_Cd[2]; + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; energy_dist(x_Ph, v_T, spidx2); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd_, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { metric.template transform(x_Cd_, v_T, v_Cd); } - ux1s_2(index + offset2 - p) = v_Cd[0]; - ux2s_2(index + offset2 - p) = v_Cd[1]; - ux3s_2(index + offset2 - p) = v_Cd[2]; + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; - tags_1(index + offset1 - p) = ParticleTag::alive; - tags_2(index + offset2 - p) = ParticleTag::alive; + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; if (M::CoordType == Coord::Cart) { - weights_1(index + offset1 - p) = ONE; - weights_2(index + offset2 - p) = ONE; + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; } else { const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0; - weights_1(index + offset1 - p) = wei; - weights_2(index + offset2 - p) = wei; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; } } random_pool.free_state(rand_gen); @@ -635,27 +635,27 @@ namespace kernel { if (ppc == 0) { return; } - const auto index = Kokkos::atomic_fetch_add(&idx(), ppc); - auto rand_gen = random_pool.get_state(); + auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - const auto dx3 = Random(rand_gen); - - i1s_1(index + offset1 - p) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1 - p) = dx1; - i1s_2(index + offset2 - p) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2 - p) = dx1; - - i2s_1(index + offset1 - p) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1 - p) = dx2; - i2s_2(index + offset2 - p) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2 - p) = dx2; - - i3s_1(index + offset1 - p) = static_cast(i3) - N_GHOSTS; - dx3s_1(index + offset1 - p) = dx3; - i3s_2(index + offset2 - p) = static_cast(i3) - N_GHOSTS; - dx3s_2(index + offset2 - p) = dx3; + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; + dx3s_1(index + offset1) = dx3; + i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; + dx3s_2(index + offset2) = dx3; vec_t v_T { ZERO }, v_Cd { ZERO }; energy_dist(x_Ph, v_T, spidx1); @@ -664,30 +664,30 @@ namespace kernel { } else if constexpr (S == SimEngine::GRPIC) { metric.template transform(x_Cd, v_T, v_Cd); } - ux1s_1(index + offset1 - p) = v_Cd[0]; - ux2s_1(index + offset1 - p) = v_Cd[1]; - ux3s_1(index + offset1 - p) = v_Cd[2]; + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; energy_dist(x_Ph, v_T, spidx2); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { metric.template transform(x_Cd, v_T, v_Cd); } - ux1s_2(index + offset2 - p) = v_Cd[0]; - ux2s_2(index + offset2 - p) = v_Cd[1]; - ux3s_2(index + offset2 - p) = v_Cd[2]; + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; - tags_1(index + offset1 - p) = ParticleTag::alive; - tags_2(index + offset2 - p) = ParticleTag::alive; + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; if (M::CoordType == Coord::Cart) { - weights_1(index + offset1 - p) = ONE; - weights_2(index + offset2 - p) = ONE; + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; } else { const auto wei = metric.sqrt_det_h( { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * inv_V0; - weights_1(index + offset1 - p) = wei; - weights_2(index + offset2 - p) = wei; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; } } random_pool.free_state(rand_gen); diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 5776fad38..8b668a036 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -57,7 +57,7 @@ namespace kernel { const float charge; const bool use_weights; const M metric; - const std::size_t ni2; + const int ni2; const real_t inv_n0; const unsigned short window; @@ -111,7 +111,7 @@ namespace kernel { , charge { charge } , use_weights { use_weights } , metric { metric } - , ni2 { ni2 } + , ni2 { static_cast(ni2) } , inv_n0 { inv_n0 } , window { window } , contrib { get_contrib(mass, charge) } diff --git a/src/kernels/tests/ampere_mink.cpp b/src/kernels/tests/ampere_mink.cpp index a7e2594a0..80af88fa5 100644 --- a/src/kernels/tests/ampere_mink.cpp +++ b/src/kernels/tests/ampere_mink.cpp @@ -23,8 +23,7 @@ void errorIf(bool condition, const std::string& message) { } } -Inline auto equal(const real_t& a, const real_t& b, const char* msg, const real_t acc) - -> bool { +Inline auto equal(real_t a, real_t b, const char* msg, real_t acc) -> bool { if (not(math::abs(a - b) < acc)) { printf("%.12e != %.12e [%.12e] %s\n", a, b, math::abs(a - b), msg); return false; @@ -368,4 +367,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/kernels/tests/deposit.cpp b/src/kernels/tests/deposit.cpp index 02292c643..9a8ae1cc6 100644 --- a/src/kernels/tests/deposit.cpp +++ b/src/kernels/tests/deposit.cpp @@ -29,10 +29,8 @@ void errorIf(bool condition, const std::string& message) { inline static constexpr auto epsilon = std::numeric_limits::epsilon(); -Inline auto equal(const real_t& a, - const real_t& b, - const char* msg = "", - const real_t acc = ONE) -> bool { +Inline auto equal(real_t a, real_t b, const char* msg = "", real_t acc = ONE) + -> bool { const auto eps = epsilon * acc; if (not cmp::AlmostEqual(a, b, eps)) { printf("%.12e != %.12e %s\n", a, b, msg); @@ -89,8 +87,8 @@ void testDeposit(const std::vector& res, const prtldx_t dxi = 0.53, dxf = 0.47; const prtldx_t dyi = 0.34, dyf = 0.52; - const real_t xi = (real_t)i0 + (real_t)dxi, xf = (real_t)i0 + (real_t)dxf; - const real_t yi = (real_t)j0 + (real_t)dyi, yf = (real_t)j0 + (real_t)dyf; + const real_t xi = (real_t)i0 + (real_t)dxi, xf = (real_t)i0 + (real_t)dxf; + const real_t yi = (real_t)j0 + (real_t)dyi, yf = (real_t)j0 + (real_t)dyf; const real_t xr = 0.5 * (xi + xf); const real_t yr = 0.5 * (yi + yf); diff --git a/src/kernels/tests/faraday_mink.cpp b/src/kernels/tests/faraday_mink.cpp index e11c8d1ff..74c2b9b1a 100644 --- a/src/kernels/tests/faraday_mink.cpp +++ b/src/kernels/tests/faraday_mink.cpp @@ -23,8 +23,7 @@ void errorIf(bool condition, const std::string& message) { } } -Inline auto equal(const real_t& a, const real_t& b, const char* msg, const real_t acc) - -> bool { +Inline auto equal(real_t a, real_t b, const char* msg, real_t acc) -> bool { if (not(math::abs(a - b) < acc)) { printf("%.12e != %.12e [%.12e] %s\n", a, b, math::abs(a - b), msg); return false; @@ -344,4 +343,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/kernels/tests/prtl_bc.cpp b/src/kernels/tests/prtl_bc.cpp index fa1103043..c8f9eae04 100644 --- a/src/kernels/tests/prtl_bc.cpp +++ b/src/kernels/tests/prtl_bc.cpp @@ -27,8 +27,7 @@ void errorIf(bool condition, const std::string& message = "") { } } -Inline auto equal(const real_t& a, const real_t& b, const std::string& msg) - -> bool { +Inline auto equal(real_t a, real_t b, const std::string& msg) -> bool { if (not(math::abs(a - b) < 1e-4)) { printf("%.12e != %.12e %s\n", a, b, msg.c_str()); return false; diff --git a/src/metrics/tests/coord_trans.cpp b/src/metrics/tests/coord_trans.cpp index 027b7677a..f3779a852 100644 --- a/src/metrics/tests/coord_trans.cpp +++ b/src/metrics/tests/coord_trans.cpp @@ -29,7 +29,7 @@ template Inline auto equal(const coord_t& a, const coord_t& b, const char* msg, - const real_t acc = ONE) -> bool { + real_t acc = ONE) -> bool { const auto eps = epsilon * acc; for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], eps)) { @@ -149,7 +149,7 @@ auto main(int argc, char* argv[]) -> int { testMetric>({ 128 }, ext1dcart); testMetric>(res2d, ext2dcart, 200); - testMetric>(res3d, ext3dcart, 200); + testMetric>(res3d, ext3dcart, 500); testMetric>(res2d, extsph, 10); testMetric>(res2d, extsph, 100, params); @@ -189,4 +189,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/metrics/tests/ks-qks.cpp b/src/metrics/tests/ks-qks.cpp index 195465701..bed051e16 100644 --- a/src/metrics/tests/ks-qks.cpp +++ b/src/metrics/tests/ks-qks.cpp @@ -23,11 +23,10 @@ template Inline auto equal(const vec_t& a, const vec_t& b, const char* msg, - const real_t acc = ONE) -> bool { + real_t acc = ONE) -> bool { const auto eps = epsilon * acc; for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], eps)) { - printf("%d : %.12e != %.12e %s\n", d, a[d], b[d], msg); return false; } } @@ -64,13 +63,12 @@ void testMetric(const std::vector& res, npts *= res[d]; } - unsigned long all_wrongs = 0; - const auto rg = metric.rg(); - const auto a = metric.spin(); - Kokkos::parallel_reduce( + const auto rg = metric.rg(); + const auto a = metric.spin(); + Kokkos::parallel_for( "h_ij/hij", npts, - Lambda(index_t n, unsigned long& wrongs) { + Lambda(index_t n) { tuple_t idx; unravel(n, idx, res_tup); coord_t x_Code { ZERO }; @@ -138,17 +136,35 @@ void testMetric(const std::vector& res, vec_t hij_expect { h11_expect, h22_expect, h33_expect }; vec_t h_ij_expect { h_11_expect, h_22_expect, h_33_expect }; - wrongs += not equal(h_ij_predict, h_ij_expect, "h_ij", acc); - wrongs += not equal(h_13_predict, { h_13_expect }, "h_13", acc); - wrongs += not equal(hij_predict, hij_expect, "hij", acc); - wrongs += not equal(h13_predict, { h13_expect }, "h13", acc); - }, - all_wrongs); - - errorIf(all_wrongs != 0, - "wrong h_ij/hij for " + std::to_string(M::Dim) + "D " + - std::string(metric.Label) + " with " + std::to_string(all_wrongs) + - " errors"); + if (not equal(h_ij_predict, h_ij_expect, "h_ij", acc)) { + printf("h_ij: %.12e %.12e %.12e : %.12e %.12e %.12e\n", + h_ij_predict[0], + h_ij_predict[1], + h_ij_predict[2], + h_ij_expect[0], + h_ij_expect[1], + h_ij_expect[2]); + Kokkos::abort("h_ij"); + } + if (not equal(h_13_predict, { h_13_expect }, "h_13", acc)) { + printf("h_13: %.12e : %.12e\n", h_13_predict[0], h_13_expect); + Kokkos::abort("h_13"); + } + if (not equal(hij_predict, hij_expect, "hij", acc)) { + printf("hij: %.12e %.12e %.12e : %.12e %.12e %.12e\n", + hij_predict[0], + hij_predict[1], + hij_predict[2], + hij_expect[0], + hij_expect[1], + hij_expect[2]); + Kokkos::abort("hij"); + } + if (not equal(h13_predict, { h13_expect }, "h13", acc)) { + printf("h13: %.12e : %.12e\n", h13_predict[0], h13_expect); + Kokkos::abort("h13"); + } + }); } auto main(int argc, char* argv[]) -> int { @@ -168,7 +184,7 @@ auto main(int argc, char* argv[]) -> int { testMetric>( { - 64, + 32, 42 }, { { 0.8, 10.0 }, { 0.0, constant::PI } }, @@ -182,4 +198,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/metrics/tests/minkowski.cpp b/src/metrics/tests/minkowski.cpp index 8fce450e9..1ef27b4fa 100644 --- a/src/metrics/tests/minkowski.cpp +++ b/src/metrics/tests/minkowski.cpp @@ -19,7 +19,7 @@ void errorIf(bool condition, const std::string& message) { inline static constexpr auto epsilon = std::numeric_limits::epsilon(); template -Inline auto equal(const coord_t& a, const coord_t& b, const real_t acc = ONE) +Inline auto equal(const coord_t& a, const coord_t& b, real_t acc = ONE) -> bool { for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], epsilon * acc)) { @@ -111,4 +111,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/metrics/tests/sph-qsph.cpp b/src/metrics/tests/sph-qsph.cpp index 13b1101a4..230a763e1 100644 --- a/src/metrics/tests/sph-qsph.cpp +++ b/src/metrics/tests/sph-qsph.cpp @@ -23,7 +23,7 @@ template Inline auto equal(const vec_t& a, const vec_t& b, const char* msg, - const real_t acc = ONE) -> bool { + real_t acc = ONE) -> bool { const auto eps = epsilon * acc; for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], eps)) { @@ -133,4 +133,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/metrics/tests/sr-cart-sph.cpp b/src/metrics/tests/sr-cart-sph.cpp index 126dea027..ec2f6ddc0 100644 --- a/src/metrics/tests/sr-cart-sph.cpp +++ b/src/metrics/tests/sr-cart-sph.cpp @@ -26,7 +26,7 @@ template Inline auto equal(const coord_t& a, const coord_t& b, const char* msg, - const real_t acc = ONE) -> bool { + real_t acc = ONE) -> bool { const auto eps = epsilon * acc; for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], eps)) { @@ -145,7 +145,7 @@ auto main(int argc, char* argv[]) -> int { testMetric>({ 128 }, ext1dcart); testMetric>(res2d, ext2dcart, 200); - testMetric>(res3d, ext3dcart, 200); + testMetric>(res3d, ext3dcart, 500); testMetric>(res2d, extsph, 10); testMetric>(res2d, extsph, 200, params); @@ -156,4 +156,4 @@ auto main(int argc, char* argv[]) -> int { } Kokkos::finalize(); return 0; -} \ No newline at end of file +} diff --git a/src/metrics/tests/vec_trans.cpp b/src/metrics/tests/vec_trans.cpp index a249c1d41..31015115c 100644 --- a/src/metrics/tests/vec_trans.cpp +++ b/src/metrics/tests/vec_trans.cpp @@ -29,7 +29,7 @@ template Inline auto equal(const vec_t& a, const vec_t& b, const char* msg, - const real_t acc = ONE) -> bool { + real_t acc = ONE) -> bool { const auto eps = epsilon * acc; for (unsigned short d = 0; d < D; ++d) { if (not cmp::AlmostEqual(a[d], b[d], eps)) { @@ -199,24 +199,24 @@ auto main(int argc, char* argv[]) -> int { { "h", (real_t)0.25} }; - testMetric>({ 128 }, ext1dcart); - testMetric>(res2d, ext2dcart, 200); - testMetric>(res3d, ext3dcart, 200); - testMetric>(res2d, extsph, 10); - testMetric>(res2d, extsph, 10, params); - - const auto resks = std::vector { 64, 54 }; - const auto extsgr = boundaries_t { - {0.8, 50.0}, - {0.0, constant::PI} - }; - const auto paramsks = std::map { - {"r0", -TWO}, - { "h", ZERO}, - { "a", (real_t)0.95} - }; - testMetric>(resks, extsgr, 150, paramsks); - + // testMetric>({ 128 }, ext1dcart); + // testMetric>(res2d, ext2dcart, 200); + // testMetric>(res3d, ext3dcart, 200); + // testMetric>(res2d, extsph, 10); + // testMetric>(res2d, extsph, 10, params); + + // const auto resks = std::vector { 64, 54 }; + // const auto extsgr = boundaries_t { + // {0.8, 50.0}, + // {0.0, constant::PI} + // }; + // const auto paramsks = std::map { + // {"r0", -TWO}, + // { "h", ZERO}, + // { "a", (real_t)0.95} + // }; + // testMetric>(resks, extsgr, 150, paramsks); + // const auto resqks = std::vector { 64, 42 }; const auto extqks = boundaries_t { {0.8, 10.0}, @@ -228,13 +228,13 @@ auto main(int argc, char* argv[]) -> int { { "a", (real_t)0.8} }; testMetric>(resqks, extqks, 500, paramsqks); - - const auto resks0 = std::vector { 64, 54 }; - const auto extks0 = boundaries_t { - {0.5, 20.0}, - {0.0, constant::PI} - }; - testMetric>(resks0, extks0, 10); + // + // const auto resks0 = std::vector { 64, 54 }; + // const auto extks0 = boundaries_t { + // {0.5, 20.0}, + // {0.0, constant::PI} + // }; + // testMetric>(resks0, extks0, 10); } catch (std::exception& e) { std::cerr << e.what() << std::endl; diff --git a/src/output/fields.cpp b/src/output/fields.cpp index e6a1f4e12..aa5a752d4 100644 --- a/src/output/fields.cpp +++ b/src/output/fields.cpp @@ -24,7 +24,11 @@ namespace out { const auto pos = name.find("_"); auto name_raw = (pos == std::string::npos) ? name : name.substr(0, pos); name_raw = name_raw.substr(0, name_raw.find_first_of("0123ijxyzt")); - m_id = FldsID::pick(fmt::toLower(name_raw).c_str()); + if (FldsID::contains(fmt::toLower(name_raw).c_str())) { + m_id = FldsID::pick(fmt::toLower(name_raw).c_str()); + } else { + m_id = FldsID::Custom; + } // determine the species and components to output if (is_moment()) { species = InterpretSpecies(name); @@ -41,11 +45,11 @@ namespace out { // energy-momentum tensor comp = InterpretComponents({ name.substr(1, 1), name.substr(2, 1) }); } else { - // scalar (Rho, divE, etc.) + // scalar (Rho, divE, Custom, etc.) comp = {}; } // data preparation flags - if (not is_moment()) { + if (not is_moment() && not is_custom()) { if (S == SimEngine::SRPIC) { prepare_flag = PrepareOutput::ConvertToHat; } else { @@ -58,7 +62,8 @@ namespace out { interp_flag = PrepareOutput::InterpToCellCenterFromEdges; } else if (is_field() || is_gr_aux_field()) { interp_flag = PrepareOutput::InterpToCellCenterFromFaces; - } else if (not(is_moment() || is_vpotential() || is_divergence())) { + } else if ( + not(is_moment() || is_vpotential() || is_divergence() || is_custom())) { raise::Error("Unrecognized field type for output", HERE); } } diff --git a/src/output/fields.h b/src/output/fields.h index d1399e981..a520a246d 100644 --- a/src/output/fields.h +++ b/src/output/fields.h @@ -77,27 +77,37 @@ namespace out { return (id() == FldsID::A); } + [[nodiscard]] + auto is_custom() const -> bool { + return (id() == FldsID::Custom); + } + [[nodiscard]] inline auto name() const -> std::string { // generate the name - auto tmp = std::string(id().to_string()); - if (id() == FldsID::T) { - tmp += m_name.substr(1, 2); - } else if (id() == FldsID::A) { - tmp += "3"; - } else if (is_field()) { - tmp += "i"; - } - if (species.size() > 0) { - tmp += "_"; - for (auto& s : species) { - tmp += std::to_string(s); + std::string tmp; + if (id() == FldsID::Custom) { + tmp = m_name; + } else { + tmp = std::string(id().to_string()); + if (id() == FldsID::T) { + tmp += m_name.substr(1, 2); + } else if (id() == FldsID::A) { + tmp += "3"; + } else if (is_field()) { + tmp += "i"; + } + if (species.size() > 0) { tmp += "_"; + for (auto& s : species) { + tmp += std::to_string(s); + tmp += "_"; + } + tmp.pop_back(); } - tmp.pop_back(); + // capitalize the first letter + tmp[0] = std::toupper(tmp[0]); } - // capitalize the first letter - tmp[0] = std::toupper(tmp[0]); return "f" + tmp; }