diff --git a/.gitmodules b/.gitmodules index 577d08ea6..835fbe5b8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "extern/Kokkos"] path = extern/Kokkos url = https://github.com/kokkos/kokkos.git +[submodule "extern/entity-pgens"] + path = extern/entity-pgens + url = https://github.com/entity-toolkit/entity-pgens.git diff --git a/.vscode/.taplo.toml b/.vscode/.taplo.toml index 0bfa6bec9..c24ab0926 100644 --- a/.vscode/.taplo.toml +++ b/.vscode/.taplo.toml @@ -1,4 +1,4 @@ -include = ["input.example.toml", "setups/**/*.toml"] +include = ["input.example.toml", "pgens/**/*.toml"] exclude = [".taplo.toml"] [formatting] diff --git a/cmake/config.cmake b/cmake/config.cmake index 7214812cd..97ed658e3 100644 --- a/cmake/config.cmake +++ b/cmake/config.cmake @@ -18,25 +18,60 @@ endfunction() # ---------------------------- Problem generator --------------------------- # function(set_problem_generator pgen_name) - file(GLOB_RECURSE PGENS "${CMAKE_CURRENT_SOURCE_DIR}/setups/**/pgen.hpp" - "${CMAKE_CURRENT_SOURCE_DIR}/setups/pgen.hpp") + if(pgen_name STREQUAL ".") + message(FATAL_ERROR "Problem generator not specified") + endif() + + file(GLOB_RECURSE PGENS "${CMAKE_CURRENT_SOURCE_DIR}/pgens/**/pgen.hpp") + foreach(PGEN ${PGENS}) get_filename_component(PGEN_NAME ${PGEN} DIRECTORY) - string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/setups/" "" PGEN_NAME - ${PGEN_NAME}) - string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/setups" "" PGEN_NAME + string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/pgens/" "" PGEN_NAME ${PGEN_NAME}) list(APPEND PGEN_NAMES ${PGEN_NAME}) endforeach() + list(FIND PGEN_NAMES ${pgen_name} PGEN_FOUND) - if(NOT ${pgen_name} STREQUAL "." AND ${PGEN_FOUND} EQUAL -1) - message(FATAL_ERROR "Invalid problem generator: " - "${pgen_name}\nValid options are: ${PGEN_NAMES}") + + file(GLOB_RECURSE EXTRA_PGENS + "${CMAKE_CURRENT_SOURCE_DIR}/extern/entity-pgens/**/pgen.hpp") + foreach(EXTRA_PGEN ${EXTRA_PGENS}) + get_filename_component(EXTRA_PGEN_NAME ${EXTRA_PGEN} DIRECTORY) + string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/extern/entity-pgens/" "" + EXTRA_PGEN_NAME ${EXTRA_PGEN_NAME}) + list(APPEND PGEN_NAMES "pgens/${EXTRA_PGEN_NAME}") + endforeach() + + if(${PGEN_FOUND} EQUAL -1) + if(${pgen_name} MATCHES "^pgens/") + get_filename_component(pgen_name ${pgen_name} NAME) + set(pgen_path + "${CMAKE_CURRENT_SOURCE_DIR}/extern/entity-pgens/${pgen_name}") + set(pgen_name "pgens/${pgen_name}") + else() + set(pgen_path ${pgen_name}) + get_filename_component(pgen_path ${pgen_path} ABSOLUTE) + string(REGEX REPLACE ".*/" "" pgen_name ${pgen_name}) + list(APPEND PGEN_NAMES ${pgen_name}) + endif() + else() + set(pgen_path ${CMAKE_CURRENT_SOURCE_DIR}/pgens/${pgen_name}) + endif() + + file(GLOB_RECURSE PGEN_FILES "${pgen_path}/pgen.hpp") + if(NOT PGEN_FILES) + message(FATAL_ERROR "pgen.hpp file not found in ${pgen_path}") endif() + + add_library(ntt_pgen INTERFACE) + target_link_libraries(ntt_pgen INTERFACE ntt_global ntt_framework + ntt_archetypes ntt_kernels) + + target_include_directories(ntt_pgen INTERFACE ${pgen_path}) + set(PGEN ${pgen_name} PARENT_SCOPE) - include_directories(${CMAKE_CURRENT_SOURCE_DIR}/setups/${pgen_name}) set(PGEN_FOUND TRUE PARENT_SCOPE) diff --git a/cmake/report.cmake b/cmake/report.cmake index b0e299d87..7bd623943 100644 --- a/cmake/report.cmake +++ b/cmake/report.cmake @@ -4,7 +4,7 @@ if(${PGEN_FOUND}) "pgen" "${problem_generators}" ${PGEN} - ${default_pgen} + "" "${Blue}" PGEN_REPORT 0) diff --git a/extern/entity-pgens b/extern/entity-pgens new file mode 160000 index 000000000..386eefc80 --- /dev/null +++ b/extern/entity-pgens @@ -0,0 +1 @@ +Subproject commit 386eefc80e2f63d0e29168869f881dd0b288952d diff --git a/setups/CMakeLists.txt b/pgens/CMakeLists.txt similarity index 76% rename from setups/CMakeLists.txt rename to pgens/CMakeLists.txt index c92c1d345..e3d047a98 100644 --- a/setups/CMakeLists.txt +++ b/pgens/CMakeLists.txt @@ -3,17 +3,7 @@ # # @includes: # -# * ../ -# -# @depends: -# -# * ntt_pgen [required] -# -# @uses: -# -# * kokkos [required] -# * plog [required] -# * mpi [optional] +# * ../src/ # ------------------------------ add_library(ntt_pgen INTERFACE) @@ -22,4 +12,3 @@ target_link_libraries(ntt_pgen INTERFACE ntt_global ntt_framework target_include_directories(ntt_pgen INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/${PGEN}) - diff --git a/setups/srpic/magnetosphere/magnetosphere.py b/pgens/magnetosphere/magnetosphere.py similarity index 100% rename from setups/srpic/magnetosphere/magnetosphere.py rename to pgens/magnetosphere/magnetosphere.py diff --git a/setups/srpic/magnetosphere/magnetosphere.toml b/pgens/magnetosphere/magnetosphere.toml similarity index 100% rename from setups/srpic/magnetosphere/magnetosphere.toml rename to pgens/magnetosphere/magnetosphere.toml diff --git a/setups/srpic/magnetosphere/pgen.hpp b/pgens/magnetosphere/pgen.hpp similarity index 100% rename from setups/srpic/magnetosphere/pgen.hpp rename to pgens/magnetosphere/pgen.hpp diff --git a/setups/pgen.hpp b/pgens/pgen.hpp similarity index 100% rename from setups/pgen.hpp rename to pgens/pgen.hpp diff --git a/setups/srpic/reconnection/pgen.hpp b/pgens/reconnection/pgen.hpp similarity index 100% rename from setups/srpic/reconnection/pgen.hpp rename to pgens/reconnection/pgen.hpp diff --git a/setups/srpic/reconnection/reconnection.toml b/pgens/reconnection/reconnection.toml similarity index 100% rename from setups/srpic/reconnection/reconnection.toml rename to pgens/reconnection/reconnection.toml diff --git a/setups/srpic/shock/pgen.hpp b/pgens/shock/pgen.hpp similarity index 100% rename from setups/srpic/shock/pgen.hpp rename to pgens/shock/pgen.hpp diff --git a/setups/legacy/shocktest/shock.py b/pgens/shock/shock.py similarity index 100% rename from setups/legacy/shocktest/shock.py rename to pgens/shock/shock.py diff --git a/setups/srpic/shock/shock.toml b/pgens/shock/shock.toml similarity index 100% rename from setups/srpic/shock/shock.toml rename to pgens/shock/shock.toml diff --git a/setups/srpic/streaming/pgen.hpp b/pgens/streaming/pgen.hpp similarity index 100% rename from setups/srpic/streaming/pgen.hpp rename to pgens/streaming/pgen.hpp diff --git a/setups/srpic/streaming/twostream.toml b/pgens/streaming/twostream.toml similarity index 100% rename from setups/srpic/streaming/twostream.toml rename to pgens/streaming/twostream.toml diff --git a/setups/srpic/streaming/weibel.toml b/pgens/streaming/weibel.toml similarity index 100% rename from setups/srpic/streaming/weibel.toml rename to pgens/streaming/weibel.toml diff --git a/setups/srpic/turbulence/pgen.hpp b/pgens/turbulence/pgen.hpp similarity index 100% rename from setups/srpic/turbulence/pgen.hpp rename to pgens/turbulence/pgen.hpp diff --git a/setups/srpic/turbulence/turbulence.toml b/pgens/turbulence/turbulence.toml similarity index 100% rename from setups/srpic/turbulence/turbulence.toml rename to pgens/turbulence/turbulence.toml diff --git a/setups/grpic/pgen_grpic_example.hpp b/setups/grpic/pgen_grpic_example.hpp deleted file mode 100644 index f553ae849..000000000 --- a/setups/grpic/pgen_grpic_example.hpp +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/traits.h" - -#include "archetypes/problem_generator.h" - -namespace user { - using namespace ntt; - - template - struct PGen : public 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 ProblemGenerator::D; - using ProblemGenerator::C; - using ProblemGenerator::params; - using ProblemGenerator::domain; - - inline PGen(SimulationParams& p, const Metadomain&) : - ProblemGenerator(p) {} - - inline PGen() {} - }; - -} // namespace user - -#endif diff --git a/setups/legacy/em_vacuum/em_vacuum.py b/setups/legacy/em_vacuum/em_vacuum.py deleted file mode 100644 index 13a62ea3b..000000000 --- a/setups/legacy/em_vacuum/em_vacuum.py +++ /dev/null @@ -1,14 +0,0 @@ -import nt2.read as nt2r -import matplotlib.pyplot as plt - -data = nt2r.Data("em_vacuum.h5") - - -def plot(ti): - fig = plt.figure(figsize=(10, 5), dpi=150) - ax = fig.add_subplot(121) - data.Bz.isel(t=ti).plot(ax=ax, cmap="BrBG") - ax = fig.add_subplot(122) - data.Ey.isel(t=ti).plot(ax=ax, cmap="RdBu_r") - for ax in fig.axes[::2]: - ax.set_aspect("equal") diff --git a/setups/legacy/em_vacuum/em_vacuum.toml b/setups/legacy/em_vacuum/em_vacuum.toml deleted file mode 100644 index 23381b1c6..000000000 --- a/setups/legacy/em_vacuum/em_vacuum.toml +++ /dev/null @@ -1,43 +0,0 @@ -[simulation] - name = "em_vacuum" - engine = "srpic" - runtime = 2.0 - -[grid] - resolution = [256, 512] - extent = [[-1.0, 1.0], [-2.0, 2.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.01 - -[algorithms] - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 1.0 - -[setup] - amplitude = 1.0 - kx1 = 1 - kx2 = 1 - kx3 = 0 - -[output] - format = "hdf5" - interval_time = 0.1 - - [output.fields] - quantities = ["E", "B"] - -[diagnostics] - colored_stdout = true diff --git a/setups/legacy/em_vacuum/pgen.hpp b/setups/legacy/em_vacuum/pgen.hpp deleted file mode 100644 index 52368cbd8..000000000 --- a/setups/legacy/em_vacuum/pgen.hpp +++ /dev/null @@ -1,108 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/comparators.h" -#include "utils/numeric.h" - -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct InitFields { - InitFields(real_t a, real_t sx1, real_t sx2, real_t sx3, int k1, int k2, int k3) - : amplitude { a } - , kx1 { (sx1 > ZERO) ? (real_t)(constant::TWO_PI) * (real_t)k1 / sx1 : ZERO } - , kx2 { (sx2 > ZERO) ? (real_t)(constant::TWO_PI) * (real_t)k2 / sx2 : ZERO } - , kx3 { (sx3 > ZERO) ? (real_t)(constant::TWO_PI) * (real_t)k3 / sx3 : ZERO } - , kmag13 { math::sqrt(SQR(kx1) + SQR(kx3)) } - , kmag { math::sqrt(SQR(kx1) + SQR(kx2) + SQR(kx3)) } { - raise::ErrorIf(cmp::AlmostZero_host(kx1) and cmp::AlmostZero_host(kx3), - "kx1 and kx3 cannot be zero", - HERE); - } - - // B is in k x y - // E is in -k x B - - Inline auto arg(const coord_t& x_Ph) const -> real_t { - if constexpr (D == Dim::_1D) { - return kx1 * x_Ph[0]; - } else if constexpr (D == Dim::_2D) { - return kx1 * x_Ph[0] + kx2 * x_Ph[1]; - } else { - return kx1 * x_Ph[0] + kx2 * x_Ph[1] + kx3 * x_Ph[2]; - } - } - - Inline auto ex1(const coord_t& x_Ph) const -> real_t { - return -amplitude * kx1 * kx2 / (kmag13 * kmag) * math::sin(arg(x_Ph)); - } - - Inline auto ex2(const coord_t& x_Ph) const -> real_t { - return amplitude * (SQR(kx1) + SQR(kx3)) / (kmag13 * kmag) * - math::sin(arg(x_Ph)); - } - - Inline auto ex3(const coord_t& x_Ph) const -> real_t { - return -amplitude * kx3 * kx2 / (kmag13 * kmag) * math::sin(arg(x_Ph)); - } - - Inline auto bx1(const coord_t& x_Ph) const -> real_t { - return -amplitude * (kx3 / kmag13) * math::sin(arg(x_Ph)); - } - - // skipping bx2 - - Inline auto bx3(const coord_t& x_Ph) const -> real_t { - return amplitude * (kx1 / kmag13) * math::sin(arg(x_Ph)); - } - - private: - const real_t amplitude; - const real_t kx1, kx2, kx3, kmag13, kmag; - }; - - 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; - - const real_t amplitude; - const int kx1, kx2, kx3; - const real_t sx1, sx2, sx3; - InitFields init_flds; - - inline PGen(const SimulationParams& p, const Metadomain& global_domain) - : arch::ProblemGenerator { p } - , amplitude { params.template get("setup.amplitude", 1.0) } - , kx1 { params.template get("setup.kx1", 1) } - , kx2 { params.template get("setup.kx2", 0) } - , kx3 { params.template get("setup.kx3", 0) } - , sx1 { global_domain.mesh().extent(in::x1).second - - global_domain.mesh().extent(in::x1).first } - , sx2 { global_domain.mesh().extent(in::x2).second - - global_domain.mesh().extent(in::x2).first } - , sx3 { global_domain.mesh().extent(in::x3).second - - global_domain.mesh().extent(in::x3).first } - , init_flds { amplitude, sx1, sx2, sx3, kx1, kx2, kx3 } {} - }; - -} // namespace user - -#endif diff --git a/setups/legacy/example/pgen.hpp b/setups/legacy/example/pgen.hpp deleted file mode 100644 index cf6c12b7a..000000000 --- a/setups/legacy/example/pgen.hpp +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/numeric.h" - -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" - -#include - -namespace user { - using namespace ntt; - - template - struct InitFields { - - InitFields(real_t a, real_t sx2, int kx2) - : amplitude { a } - , sx2 { sx2 } - , kx2 { kx2 } {} - - // only set ex2 and bx3 - - Inline auto ex2(const coord_t& x_Ph) const -> real_t { - return amplitude * math::sin(constant::TWO_PI * (x_Ph[1] / sx2) * - static_cast(kx2)); - } - - Inline auto bx3(const coord_t& x_Ph) const -> real_t { - return -amplitude * math::cos(constant::TWO_PI * (x_Ph[1] / sx2) * - static_cast(kx2)); - } - - private: - const real_t amplitude; - const real_t sx2; - const int kx2; - }; - - template - struct ExtForce { - const std::vector species { 1, 2 }; - - ExtForce() = default; - - Inline auto fx1(const spidx_t& sp, - const simtime_t& time, - const coord_t& x_Ph) const -> real_t { - (void)sp; - (void)time; - (void)x_Ph; - return ZERO; - } - - Inline auto fx2(const spidx_t& sp, - const simtime_t& time, - const coord_t& x_Ph) const -> real_t { - (void)sp; - (void)time; - (void)x_Ph; - return ZERO; - } - - Inline auto fx3(const spidx_t& sp, - const simtime_t& time, - const coord_t& x_Ph) const -> real_t { - (void)sp; - (void)time; - (void)x_Ph; - return ZERO; - } - }; - - 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; - - InitFields init_flds; - ExtForce ext_force; - - inline PGen(const SimulationParams& p, const Metadomain& global_domain) - : arch::ProblemGenerator { p } - , init_flds { params.template get("setup.amplitude", 1.0), - global_domain.mesh().extent(in::x2).second - - global_domain.mesh().extent(in::x2).first, - params.template get("setup.kx2", 2) } - , ext_force {} {} - }; - -} // namespace user - -#endif diff --git a/setups/legacy/langmuir/langmuir.py b/setups/legacy/langmuir/langmuir.py deleted file mode 100644 index a880bc00b..000000000 --- a/setups/legacy/langmuir/langmuir.py +++ /dev/null @@ -1,17 +0,0 @@ -import nt2.read as nt2r -import matplotlib.pyplot as plt - -data = nt2r.Data("langmuir.h5") - - -def plot(ti, d): - # for 2D - fig = plt.figure(figsize=(10, 5), dpi=150) - ax = fig.add_subplot(211) - d.Rho.isel(t=ti).plot(ax=ax, cmap="inferno", vmin=0, vmax=4) - ax = fig.add_subplot(212) - d.Ex.isel(t=ti).plot(ax=ax, cmap="RdBu_r", vmin=-1, vmax=1) - for ax in fig.get_axes()[::2]: - ax.set_aspect("equal") - fig.get_axes()[0].set(xlabel="", xticks=[]) - fig.get_axes()[2].set(title=None) diff --git a/setups/legacy/langmuir/langmuir.toml b/setups/legacy/langmuir/langmuir.toml deleted file mode 100644 index b054a940d..000000000 --- a/setups/legacy/langmuir/langmuir.toml +++ /dev/null @@ -1,55 +0,0 @@ -[simulation] - name = "langmuir" - engine = "srpic" - runtime = 1.0 - -[grid] - resolution = [2048, 512] - extent = [[0.0, 1.0], [0.0, 0.25]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.01 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 14.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e7 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e7 - -[setup] - vmax = 0.1 - nx1 = 4 - nx2 = 2 - -[output] - format = "hdf5" - interval_time = 0.0025 - - [output.fields] - quantities = ["Rho", "E"] - -[diagnostics] - colored_stdout = true diff --git a/setups/legacy/langmuir/pgen.hpp b/setups/legacy/langmuir/pgen.hpp deleted file mode 100644 index 28dbd24c5..000000000 --- a/setups/legacy/langmuir/pgen.hpp +++ /dev/null @@ -1,124 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/numeric.h" - -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/domain.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct SinEDist : public arch::EnergyDistribution { - SinEDist(const M& metric, - real_t v_max, - const std::vector& n, - const std::vector& s) - : arch::EnergyDistribution { metric } - , v_max { v_max } - , kx1 { s.size() > 0 ? static_cast(constant::TWO_PI) * - static_cast(n[0]) / s[0] - : ZERO } - , kx2 { s.size() > 1 ? static_cast(constant::TWO_PI) * - static_cast(n[1]) / s[1] - : ZERO } - , kx3 { s.size() > 2 ? static_cast(constant::TWO_PI) * - static_cast(n[2]) / s[2] - : ZERO } {} - - Inline void operator()(const coord_t& x_Ph, - vec_t& v, - spidx_t sp) const override { - if (sp == 1) { - const auto k = math::sqrt(SQR(kx1) + SQR(kx2) + SQR(kx3)); - if constexpr (M::Dim == Dim::_1D) { - v[0] = v_max * math::sin(kx1 * x_Ph[0]); - } else if constexpr (M::Dim == Dim::_2D) { - v[0] = v_max * kx1 / k * math::sin(kx1 * x_Ph[0] + kx2 * x_Ph[1]); - v[1] = v_max * kx2 / k * math::sin(kx1 * x_Ph[0] + kx2 * x_Ph[1]); - } else { - v[0] = v_max * kx1 / k * - math::sin(kx1 * x_Ph[0] + kx2 * x_Ph[1] + kx3 * x_Ph[2]); - v[1] = v_max * kx2 / k * - math::sin(kx1 * x_Ph[0] + kx2 * x_Ph[1] + kx3 * x_Ph[2]); - v[2] = v_max * kx3 / k * - math::sin(kx1 * x_Ph[0] + kx2 * x_Ph[1] + kx3 * x_Ph[2]); - } - } else { - v[0] = ZERO; - v[1] = ZERO; - v[2] = ZERO; - } - } - - private: - const real_t v_max, kx1, kx2, kx3; - }; - - 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; - - const real_t sx1, sx2, sx3; - const real_t vmax; - const int nx1, nx2, nx3; - - std::vector svec; - std::vector nvec; - - inline PGen(const SimulationParams& p, const Metadomain& global_domain) - : arch::ProblemGenerator { p } - , sx1 { global_domain.mesh().extent(in::x1).second - - global_domain.mesh().extent(in::x1).first } - , sx2 { global_domain.mesh().extent(in::x2).second - - global_domain.mesh().extent(in::x2).first } - , sx3 { global_domain.mesh().extent(in::x3).second - - global_domain.mesh().extent(in::x3).first } - , vmax { p.get("setup.vmax", 0.01) } - , nx1 { p.get("setup.nx1", 10) } - , nx2 { p.get("setup.nx2", 10) } - , nx3 { p.get("setup.nx3", 10) } { - const auto sxs = std::vector { sx1, sx2, sx3 }; - const auto nxs = std::vector { nx1, nx2, nx3 }; - for (auto d = 0u; d < M::Dim; ++d) { - svec.push_back(sxs[d]); - nvec.push_back(nxs[d]); - } - } - - inline void InitPrtls(Domain& local_domain) { - const auto energy_dist = SinEDist(local_domain.mesh.metric, - vmax, - nvec, - svec); - const auto injector = arch::UniformInjector(energy_dist, - { 1, 2 }); - arch::InjectUniform>(params, - local_domain, - injector, - 1.0); - } - }; - -} // namespace user - -#endif diff --git a/setups/legacy/magnetar/magnetar.py b/setups/legacy/magnetar/magnetar.py deleted file mode 100644 index 0bbf790e5..000000000 --- a/setups/legacy/magnetar/magnetar.py +++ /dev/null @@ -1,9 +0,0 @@ -import nt2.read as nt2r -import matplotlib as mpl - -data = nt2r.Data("magnetar.h5") - -def plot (ti, data): - (data.Bph*(data.r*np.sin(data.th))).isel(t=ti).polar.pcolor( - norm=mpl.colors.Normalize(vmin=-0.075, vmax=0.075), - cmap="PuOr") \ No newline at end of file diff --git a/setups/legacy/magnetar/magnetar.toml b/setups/legacy/magnetar/magnetar.toml deleted file mode 100644 index fab2eb01c..000000000 --- a/setups/legacy/magnetar/magnetar.toml +++ /dev/null @@ -1,108 +0,0 @@ -[simulation] - name = "magnetar" - engine = "srpic" - runtime = 50.0 - -[grid] - resolution = [2048, 1024] - extent = [[1.0, 400.0]] - - [grid.metric] - metric = "qspherical" - - [grid.boundaries] - fields = [["ATMOSPHERE", "ABSORB"]] - particles = [["ATMOSPHERE", "ABSORB"]] - - [grid.boundaries.absorb] - ds = 1.0 - - [grid.boundaries.atmosphere] - temperature = 0.1 - density = 40.0 - height = 0.02 - species = [1, 2] - ds = 0.5 - -[scales] - larmor0 = 1e-5 - skindepth0 = 0.01 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - - [algorithms.gca] - e_ovr_b_max = 0.9 - larmor_max = 100.0 - -[particles] - ppc0 = 4.0 - use_weights = true - clear_interval = 100 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 5e7 - pusher = "Boris,GCA" - -[setup] - Bsurf = 1.0 - omega = 0.0125 - pp_thres = 10.0 - gamma_pairs = 1.75 - -[output] - format = "hdf5" - - [output.fields] - interval_time = 0.5 - quantities = ["N_1", "N_2", "N_3", "N_4", "N_5", "N_6", "B", "E", "J"] - - [output.particles] - enable = false - - [output.spectra] - enable = false - -[diagnostics] - interval = 1 diff --git a/setups/legacy/magnetar/pgen.hpp b/setups/legacy/magnetar/pgen.hpp deleted file mode 100644 index 10f98ea5d..000000000 --- a/setups/legacy/magnetar/pgen.hpp +++ /dev/null @@ -1,281 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/numeric.h" - -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct InitFields { - InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {} - - Inline auto bx1(const coord_t& x_Ph) const -> real_t { - return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); - } - - Inline auto bx2(const coord_t& x_Ph) const -> real_t { - return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); - } - - private: - const real_t Bsurf, Rstar; - }; - - template - struct DriveFields : public InitFields { - DriveFields(real_t time, real_t bsurf, real_t rstar, real_t omega) - : InitFields { bsurf, rstar } - , time { time } - , Omega { omega } {} - - using InitFields::bx1; - using InitFields::bx2; - - Inline auto bx3(const coord_t&) const -> real_t { - return ZERO; - } - - Inline auto ex1(const coord_t& x_Ph) const -> real_t { - auto sigma = (x_Ph[1] - HALF * constant::PI) / - (static_cast(0.2) * constant::PI); - return Omega * bx2(x_Ph) * x_Ph[0] * math::sin(x_Ph[1]) * sigma * - math::exp((ONE - SQR(SQR(sigma))) * INV_4); - } - - Inline auto ex2(const coord_t& x_Ph) const -> real_t { - auto sigma = (x_Ph[1] - HALF * constant::PI) / - (static_cast(0.2) * constant::PI); - return -Omega * bx1(x_Ph) * x_Ph[0] * math::sin(x_Ph[1]) * sigma * - math::exp((ONE - SQR(SQR(sigma))) * INV_4); - } - - Inline auto ex3(const coord_t&) const -> real_t { - return ZERO; - } - - private: - const real_t time, Omega; - }; - - 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; - - const Metadomain& global_domain; - - const real_t Bsurf, Rstar, Omega, gamma_pairs, pp_thres; - InitFields init_flds; - - inline PGen(const SimulationParams& p, const Metadomain& m) - : arch::ProblemGenerator(p) - , global_domain { m } - , Bsurf { p.template get("setup.Bsurf", ONE) } - , Rstar { m.mesh().extent(in::x1).first } - , Omega { p.template get("setup.omega") } - , pp_thres { p.template get("setup.pp_thres") } - , gamma_pairs { p.template get("setup.gamma_pairs") } - , init_flds { Bsurf, Rstar } {} - - inline PGen() {} - - auto AtmFields(real_t time) const -> DriveFields { - const real_t omega_t = - Omega * - ((ONE - math::tanh((static_cast(5.0) - time) * HALF)) * - (ONE + (-ONE + math::tanh((static_cast(45.0) - time) * HALF)) * - HALF)) * - HALF; - return DriveFields { time, Bsurf, Rstar, omega_t }; - } - - auto MatchFields(real_t) const -> InitFields { - return InitFields { Bsurf, Rstar }; - } - - void CustomPostStep(std::size_t, long double, Domain& domain) { - - // Ad-hoc PP kernel - { - - auto& species2_e = domain.species[2]; - auto& species2_p = domain.species[3]; - auto& species3_e = domain.species[4]; - auto& species3_p = domain.species[5]; - auto metric = domain.mesh.metric; - auto pp_thres_ = this->pp_thres; - auto gamma_pairs_ = this->gamma_pairs; - - for (std::size_t s { 0 }; s < 6; ++s) { - if (s == 1) { - continue; - } - - array_t elec_ind("elec_ind"); - array_t pos_ind("pos_ind"); - - auto offset_e = species3_e.npart(); - auto offset_p = species3_p.npart(); - - auto ux1_e = species3_e.ux1; - auto ux2_e = species3_e.ux2; - auto ux3_e = species3_e.ux3; - auto i1_e = species3_e.i1; - auto i2_e = species3_e.i2; - auto dx1_e = species3_e.dx1; - auto dx2_e = species3_e.dx2; - auto phi_e = species3_e.phi; - auto weight_e = species3_e.weight; - auto tag_e = species3_e.tag; - - auto ux1_p = species3_p.ux1; - auto ux2_p = species3_p.ux2; - auto ux3_p = species3_p.ux3; - auto i1_p = species3_p.i1; - auto i2_p = species3_p.i2; - auto dx1_p = species3_p.dx1; - auto dx2_p = species3_p.dx2; - auto phi_p = species3_p.phi; - auto weight_p = species3_p.weight; - auto tag_p = species3_p.tag; - - if (s == 0) { - - offset_e = species2_e.npart(); - offset_p = species2_p.npart(); - - ux1_e = species2_e.ux1; - ux2_e = species2_e.ux2; - ux3_e = species2_e.ux3; - i1_e = species2_e.i1; - i2_e = species2_e.i2; - dx1_e = species2_e.dx1; - dx2_e = species2_e.dx2; - phi_e = species2_e.phi; - weight_e = species2_e.weight; - tag_e = species2_e.tag; - - ux1_p = species2_p.ux1; - ux2_p = species2_p.ux2; - ux3_p = species2_p.ux3; - i1_p = species2_p.i1; - i2_p = species2_p.i2; - dx1_p = species2_p.dx1; - dx2_p = species2_p.dx2; - phi_p = species2_p.phi; - weight_p = species2_p.weight; - tag_p = species2_p.tag; - } - - auto& species = domain.species[s]; - auto ux1 = species.ux1; - auto ux2 = species.ux2; - auto ux3 = species.ux3; - auto i1 = species.i1; - auto i2 = species.i2; - auto dx1 = species.dx1; - auto dx2 = species.dx2; - auto phi = species.phi; - auto weight = species.weight; - auto tag = species.tag; - - Kokkos::parallel_for( - "InjectPairs", - species.rangeActiveParticles(), - Lambda(index_t p) { - if (tag(p) == ParticleTag::dead) { - return; - } - - auto px = ux1(p); - auto py = ux2(p); - auto pz = ux3(p); - auto gamma = math::sqrt(ONE + SQR(px) + SQR(py) + SQR(pz)); - - const coord_t xCd { static_cast(i1(p)) + dx1(p), - static_cast(i2(p)) + dx2(p) }; - - coord_t xPh { ZERO }; - metric.template convert(xCd, xPh); - - if ((gamma > pp_thres_) && (math::sin(xPh[1]) > 0.1)) { - - auto new_gamma = gamma - 2.0 * gamma_pairs_; - auto new_fac = math::sqrt(SQR(new_gamma) - 1.0) / - math::sqrt(SQR(gamma) - 1.0); - auto pair_fac = math::sqrt(SQR(gamma_pairs_) - 1.0) / - math::sqrt(SQR(gamma) - 1.0); - - auto elec_p = Kokkos::atomic_fetch_add(&elec_ind(), 1); - auto pos_p = Kokkos::atomic_fetch_add(&pos_ind(), 1); - - i1_e(elec_p + offset_e) = i1(p); - dx1_e(elec_p + offset_e) = dx1(p); - i2_e(elec_p + offset_e) = i2(p); - dx2_e(elec_p + offset_e) = dx2(p); - phi_e(elec_p + offset_e) = phi(p); - ux1_e(elec_p + offset_e) = px * pair_fac; - ux2_e(elec_p + offset_e) = py * pair_fac; - ux3_e(elec_p + offset_e) = pz * pair_fac; - weight_e(elec_p + offset_e) = weight(p); - tag_e(elec_p + offset_e) = ParticleTag::alive; - - i1_p(pos_p + offset_p) = i1(p); - dx1_p(pos_p + offset_p) = dx1(p); - i2_p(pos_p + offset_p) = i2(p); - dx2_p(pos_p + offset_p) = dx2(p); - phi_p(pos_p + offset_p) = phi(p); - ux1_p(pos_p + offset_p) = px * pair_fac; - ux2_p(pos_p + offset_p) = py * pair_fac; - ux3_p(pos_p + offset_p) = pz * pair_fac; - weight_p(pos_p + offset_p) = weight(p); - tag_p(pos_p + offset_p) = ParticleTag::alive; - - ux1(p) *= new_fac; - ux2(p) *= new_fac; - ux3(p) *= new_fac; - } - }); - - auto elec_ind_h = Kokkos::create_mirror(elec_ind); - Kokkos::deep_copy(elec_ind_h, elec_ind); - if (s == 0) { - species2_e.set_npart(offset_e + elec_ind_h()); - } else { - species3_e.set_npart(offset_e + elec_ind_h()); - } - - auto pos_ind_h = Kokkos::create_mirror(pos_ind); - Kokkos::deep_copy(pos_ind_h, pos_ind); - if (s == 0) { - species2_p.set_npart(offset_p + pos_ind_h()); - } else { - species3_p.set_npart(offset_p + pos_ind_h()); - } - } - } // Ad-hoc PP kernel - } - }; - -} // namespace user - -#endif diff --git a/setups/legacy/magpump/pgen.hpp b/setups/legacy/magpump/pgen.hpp deleted file mode 100644 index 045552aff..000000000 --- a/setups/legacy/magpump/pgen.hpp +++ /dev/null @@ -1,170 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/traits.h" - -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" - -#include - -namespace user { - using namespace ntt; - - template - struct InitFields { - InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {} - - Inline auto bx1(const coord_t& x_Ph) const -> real_t { - return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); - } - - Inline auto bx2(const coord_t& x_Ph) const -> real_t { - return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); - } - - private: - const real_t Bsurf, Rstar; - }; - - template - struct DriveFields : public InitFields { - DriveFields(real_t time, real_t bsurf, real_t rstar) - : InitFields { bsurf, rstar } - , time { time } {} - - using InitFields::bx1; - using InitFields::bx2; - - Inline auto bx3(const coord_t&) const -> real_t { - return ZERO; - } - - Inline auto ex1(const coord_t&) const -> real_t { - return ZERO; - } - - Inline auto ex2(const coord_t& x_Ph) const -> real_t { - return ZERO; - } - - Inline auto ex3(const coord_t&) const -> real_t { - return ZERO; - } - - private: - const real_t time; - }; - - template - struct Inflow : public arch::EnergyDistribution { - Inflow(const M& metric, real_t vin) - : arch::EnergyDistribution { metric } - , vin { vin } {} - - Inline void operator()(const coord_t&, - vec_t& v_Ph, - spidx_t) const override { - v_Ph[0] = -vin; - } - - private: - const real_t vin; - }; - - template - struct Sphere : public arch::SpatialDistribution { - Sphere(const M& metric, real_t r0, real_t dr) - : arch::SpatialDistribution { metric } - , r0 { r0 } - , dr { dr } {} - - Inline auto operator()(const coord_t& x_Ph) const -> real_t override { - return math::exp(-SQR((x_Ph[0] - r0) / dr)) * - (x_Ph[1] > 0.25 && x_Ph[1] < constant::PI - 0.25); - } - - private: - const real_t r0, dr; - }; - - template - struct PGen : public arch::ProblemGenerator { - static constexpr auto engines { traits::compatible_with::value }; - static constexpr auto metrics { - traits::compatible_with::value - }; - static constexpr auto dimensions { traits::compatible_with::value }; - - using arch::ProblemGenerator::D; - using arch::ProblemGenerator::C; - using arch::ProblemGenerator::params; - - const real_t Bsurf, pump_period, pump_ampl, pump_radius, Rstar; - const real_t vin, drinj; - InitFields init_flds; - - inline PGen(const SimulationParams& p, const Metadomain& m) - : arch::ProblemGenerator { p } - , Bsurf { p.template get("setup.Bsurf", ONE) } - , pump_period { p.template get("setup.pump_period") } - , pump_ampl { p.template get("setup.pump_ampl") } - , pump_radius { p.template get("setup.pump_radius") } - , Rstar { m.mesh().extent(in::x1).first } - , vin { p.template get("setup.vin") } - , drinj { p.template get("setup.drinj") } - , init_flds { Bsurf, Rstar } {} - - auto FieldDriver(real_t time) const -> DriveFields { - return DriveFields { time, Bsurf, Rstar }; - } - - void CustomPostStep(std::size_t, long double time, Domain& domain) { - const real_t radius = pump_radius + - pump_ampl * - math::sin(time * constant::TWO_PI / pump_period); - const real_t dr = 1.0; - const auto& metric = domain.mesh.metric; - auto EM = domain.fields.em; - Kokkos::parallel_for( - "outerBC", - domain.mesh.rangeActiveCells(), - Lambda(index_t i1, index_t i2) { - const auto i1_ = COORD(i1), i2_ = COORD(i2); - const auto r = metric.template convert<1, Crd::Cd, Crd::Ph>(i1_); - if (r > radius - 5 * dr) { - const auto smooth = HALF * (ONE - math::tanh((r - radius) / dr)); - EM(i1, i2, em::ex1) = smooth * EM(i1, i2, em::ex1); - EM(i1, i2, em::ex2) = smooth * EM(i1, i2, em::ex2); - EM(i1, i2, em::ex3) = smooth * EM(i1, i2, em::ex3); - EM(i1, i2, em::bx1) = smooth * EM(i1, i2, em::bx1); - EM(i1, i2, em::bx2) = smooth * EM(i1, i2, em::bx2); - EM(i1, i2, em::bx3) = smooth * EM(i1, i2, em::bx3); - } - }); - - if (time < pump_period * 0.25) { - const auto energy_dist = Inflow(domain.mesh.metric, vin); - const auto spatial_dist = Sphere(domain.mesh.metric, radius, drinj); - const auto injector = arch::NonUniformInjector( - energy_dist, - spatial_dist, - { 1, 2 }); - - arch::InjectNonUniform>( - params, - domain, - injector, - ONE, - true); - } - } - }; - -} // namespace user - -#endif diff --git a/setups/legacy/rec-gravity/pgen.hpp b/setups/legacy/rec-gravity/pgen.hpp deleted file mode 100644 index e8c461418..000000000 --- a/setups/legacy/rec-gravity/pgen.hpp +++ /dev/null @@ -1,174 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/directions.h" -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/numeric.h" - -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "archetypes/spatial_dist.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct Gravity { - const std::vector species { 1, 2 }; - - Gravity(real_t f, real_t tscale, real_t ymid) - : force { f } - , tscale { tscale } - , ymid { ymid } {} - - Inline auto fx1(const spidx_t&, const simtime_t&, const coord_t&) const - -> real_t { - return ZERO; - } - - Inline auto fx2(const spidx_t&, const simtime_t& t, const coord_t& x_Ph) const - -> real_t { - const auto sign { (x_Ph[1] < ymid) ? ONE : -ONE }; - const auto t_ { static_cast(t) }; - if (t_ > tscale) { - return sign * force; - } else { - return sign * force * (ONE - math::cos(constant::PI * t_ / tscale)) / TWO; - } - } - - Inline auto fx3(const spidx_t&, const simtime_t&, const coord_t&) const - -> real_t { - return ZERO; - } - - private: - const real_t force, tscale, ymid; - }; - - template - struct CurrentLayer : public arch::SpatialDistribution { - CurrentLayer(const M& metric, real_t width, real_t yi) - : arch::SpatialDistribution { metric } - , width { width } - , yi { yi } {} - - Inline auto operator()(const coord_t& x_Ph) const -> real_t override { - return ONE / SQR(math::cosh((x_Ph[1] - yi) / width)); - } - - private: - const real_t yi, width; - }; - - template - struct InitFields { - InitFields(real_t Bmag, real_t width, real_t y1, real_t y2) - : Bmag { Bmag } - , width { width } - , y1 { y1 } - , y2 { y2 } {} - - Inline auto bx3(const coord_t& x_Ph) const -> real_t { - return Bmag * (math::tanh((x_Ph[1] - y1) / width) - - math::tanh((x_Ph[1] - y2) / width) - 1); - } - - private: - const real_t Bmag, width, y1, y2; - }; - - 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; - - const real_t Bmag, width, overdensity, y1, y2, bg_temp; - InitFields init_flds; - - Gravity ext_force; - - inline PGen(const SimulationParams& p, const Metadomain& m) - : arch::ProblemGenerator(p) - , Bmag { p.template get("setup.Bmag", 1.0) } - , width { p.template get("setup.width") } - , overdensity { p.template get("setup.overdensity") } - , bg_temp { p.template get("setup.bg_temp") } - , y1 { m.mesh().extent(in::x2).first + - INV_4 * - (m.mesh().extent(in::x2).second - m.mesh().extent(in::x2).first) } - , y2 { m.mesh().extent(in::x2).first + - 3 * INV_4 * - (m.mesh().extent(in::x2).second - m.mesh().extent(in::x2).first) } - , init_flds { Bmag, width, y1, y2 } - , ext_force { - p.template get("setup.fmag", 0.1) * - p.template get("scales.omegaB0"), - (m.mesh().extent(in::x1).second - m.mesh().extent(in::x1).first), - INV_2 * (m.mesh().extent(in::x2).second + m.mesh().extent(in::x2).first) - } {} - - inline PGen() {} - - inline void InitPrtls(Domain& local_domain) { - // background - const auto energy_dist = arch::Maxwellian(local_domain.mesh.metric, - local_domain.random_pool, - bg_temp); - const auto injector = arch::UniformInjector( - energy_dist, - { 1, 2 }); - arch::InjectUniform(params, - local_domain, - injector, - ONE); - // current layers - const auto sigma = params.template get("scales.sigma0"); - const auto c_omp = params.template get("scales.skindepth0"); - const auto cs_drift_beta = math::sqrt(sigma) * c_omp / (width * overdensity); - const auto cs_drift_gamma = ONE / math::sqrt(ONE - SQR(cs_drift_beta)); - const auto cs_drift_u = cs_drift_beta * cs_drift_gamma; - const auto cs_temp = HALF * sigma / overdensity; - - for (auto i = 0; i < 2; ++i) { - const auto drift_vel = (i == 0) ? cs_drift_u : -cs_drift_u; - const auto y_cs = (i == 0) ? y1 : y2; - auto edist_cs = arch::Maxwellian(local_domain.mesh.metric, - local_domain.random_pool, - cs_temp, - drift_vel, - in::x1, - false); - const auto sdist_cs = CurrentLayer(local_domain.mesh.metric, - width, - y_cs); - const auto inj_cs = arch::NonUniformInjector( - edist_cs, - sdist_cs, - { 1, 2 }); - arch::InjectNonUniform(params, - local_domain, - inj_cs, - overdensity); - } - } // namespace user - }; - -} // namespace user - -#endif diff --git a/setups/legacy/rec-gravity/rec-gravity.toml b/setups/legacy/rec-gravity/rec-gravity.toml deleted file mode 100644 index f29b090e3..000000000 --- a/setups/legacy/rec-gravity/rec-gravity.toml +++ /dev/null @@ -1,54 +0,0 @@ -[simulation] - name = "rec-gravity" - engine = "srpic" - runtime = 20.0 - -[grid] - resolution = [2000, 4000] - extent = [[-0.5, 0.5], [-1.0, 1.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 3.1e-4 - skindepth0 = 1e-3 - -[algorithms] - current_filters = 8 - - [algorithms.timestep] - CFL = 0.45 - -[particles] - ppc0 = 8.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e8 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e8 - -[setup] - Bmag = 1.0 - width = 0.04 - bg_temp = 1e-4 - overdensity = 3.0 - angle = 0.0 - -[output] - format = "hdf5" - interval_time = 0.36 - - [output.fields] - quantities = ["N_1", "N_2", "E", "B", "J", "T00_1", "T00_2"] diff --git a/setups/legacy/shocktest/pgen.hpp b/setups/legacy/shocktest/pgen.hpp deleted file mode 100644 index a7bcf5c3d..000000000 --- a/setups/legacy/shocktest/pgen.hpp +++ /dev/null @@ -1,155 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/traits.h" -#include "utils/error.h" -#include "utils/numeric.h" - -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" - -#include -#include - -namespace user { - using namespace ntt; - - template - struct InitFields { - /* - Sets up magnetic and electric field components for the simulation. - Must satisfy E = -v x B for Lorentz Force to be zero. - - @param bmag: magnetic field scaling - @param btheta: magnetic field polar angle - @param bphi: magnetic field azimuthal angle - @param drift_ux: drift velocity in the x direction - */ - InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux) - : Bmag { bmag } - , Btheta { btheta * static_cast(convert::deg2rad) } - , Bphi { bphi * static_cast(convert::deg2rad) } - , Vx { drift_ux } {} - - // magnetic field components - Inline auto bx1(const coord_t& x_ph) const -> real_t { - return Bmag * math::cos(Btheta); - } - - Inline auto bx2(const coord_t& x_ph) const -> real_t { - return Bmag * math::sin(Btheta) * math::sin(Bphi); - } - - Inline auto bx3(const coord_t& x_ph) const -> real_t { - return Bmag * math::sin(Btheta) * math::cos(Bphi); - } - - // electric field components - Inline auto ex1(const coord_t& x_ph) const -> real_t { - return ZERO; - } - - Inline auto ex2(const coord_t& x_ph) const -> real_t { - return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi); - } - - Inline auto ex3(const coord_t& x_ph) const -> real_t { - return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi); - } - - private: - const real_t Btheta, Bphi, Vx, Bmag; - }; - - 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; - - const real_t drift_ux, temperature; - - const std::vector x1arr_e, x2arr_e, ux1arr_e, ux2arr_e, ux3arr_e; - const std::vector x1arr_i, x2arr_i, ux1arr_i, ux2arr_i, ux3arr_i; - - const real_t Btheta, Bphi, Bmag; - InitFields init_flds; - const Metadomain* metadomain; - - inline PGen(const SimulationParams& p, const Metadomain& m) - : arch::ProblemGenerator { p } - , drift_ux { p.template get("setup.drift_ux") } - , temperature { p.template get("setup.temperature") } - , x1arr_e { p.template get>("setup.x_e") } - , x2arr_e { p.template get>("setup.y_e") } - , ux1arr_e { p.template get>("setup.ux_e") } - , ux2arr_e { p.template get>("setup.uy_e") } - , ux3arr_e { p.template get>("setup.uz_e") } - , x1arr_i { p.template get>("setup.x_i") } - , x2arr_i { p.template get>("setup.y_i") } - , ux1arr_i { p.template get>("setup.ux_i") } - , ux2arr_i { p.template get>("setup.uy_i") } - , ux3arr_i { p.template get>("setup.uz_i") } - , Btheta { p.template get("setup.Btheta", ZERO) } - , Bmag { p.template get("setup.Bmag", ZERO) } - , Bphi { p.template get("setup.Bphi", ZERO) } - , init_flds { Bmag, Btheta, Bphi, drift_ux } - , metadomain { &m } {} - - inline PGen() {} - - auto FixFieldsConst(const bc_in&, const em& comp) const - -> std::pair { - if (comp == em::ex2) { - return { init_flds.ex2({ ZERO }), true }; - } else if (comp == em::ex3) { - return { init_flds.ex3({ ZERO }), true }; - } else { - return { ZERO, false }; - } - } - - auto MatchFields(real_t time) const -> InitFields { - return init_flds; - } - - inline void InitPrtls(Domain& domain) { - arch::InjectGlobally(*metadomain, - domain, - 1, - { - { "x1", x1arr_e }, - { "x2", x2arr_e }, - { "ux1", ux1arr_e }, - { "ux2", ux1arr_e }, - { "ux3", ux3arr_e } - }); - arch::InjectGlobally(*metadomain, - domain, - 2, - { - { "x1", x1arr_i }, - { "x2", x2arr_i }, - { "ux1", ux1arr_i }, - { "ux2", ux1arr_i }, - { "ux3", ux3arr_i } - }); - } - }; - -} // namespace user - -#endif diff --git a/setups/legacy/shocktest/shock.toml b/setups/legacy/shocktest/shock.toml deleted file mode 100644 index 6c0c9a3a0..000000000 --- a/setups/legacy/shocktest/shock.toml +++ /dev/null @@ -1,69 +0,0 @@ -[simulation] - name = "shock" - engine = "srpic" - runtime = 50.0 - -[grid] - resolution = [2048, 128] - extent = [[0.0, 10.0], [-0.3125, 0.3125]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["CONDUCTOR", "FIXED"], ["PERIODIC"]] - particles = [["REFLECT", "ABSORB"], ["PERIODIC"]] - -[scales] - larmor0 = 1e-2 - skindepth0 = 1e-2 - -[algorithms] - current_filters = 8 - fieldsolver = "false" - deposit = "false" - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 16.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e8 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e8 - -[setup] - drift_ux = 0.1 - temperature = 1e-3 - Bmag = 1.0 - Btheta = 0.0 - Bphi = 0.0 - x_e = [0.05] - y_e = [0.0] - ux_e = [-0.01] - uy_e = [0.01] - uz_e = [0.001] - x_i = [0.05] - y_i = [0.0] - ux_i = [0.01] - uy_i = [-0.01] - uz_i = [-0.001] - -[output] - interval = 1 - format = "hdf5" - - [output.fields] - quantities = ["N_1", "N_2", "E", "B", "T0i_1", "T0i_2", "J"] - - [output.debug] - ghosts = true diff --git a/setups/legacy/spider/pgen.hpp b/setups/legacy/spider/pgen.hpp deleted file mode 100644 index 27d9504b7..000000000 --- a/setups/legacy/spider/pgen.hpp +++ /dev/null @@ -1,38 +0,0 @@ -#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; - - inline PGen(const SimulationParams& p, const Metadomain& m) - : arch::ProblemGenerator(p) {} - - inline PGen() {} - }; - -} // namespace user - -#endif diff --git a/setups/srpic/pgen_srpic_example.hpp b/setups/srpic/pgen_srpic_example.hpp deleted file mode 100644 index de3bae408..000000000 --- a/setups/srpic/pgen_srpic_example.hpp +++ /dev/null @@ -1,39 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.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; - - inline PGen(const SimulationParams& p, const Metadomain&) - : arch::ProblemGenerator(p) {} - - inline PGen() {} - }; - -} // namespace user - -#endif diff --git a/setups/srpic/shock/shock.py b/setups/srpic/shock/shock.py deleted file mode 100644 index dc1565572..000000000 --- a/setups/srpic/shock/shock.py +++ /dev/null @@ -1,75 +0,0 @@ -import nt2.read as nt2r -import matplotlib.pyplot as plt -import matplotlib as mpl - -data = nt2r.Data("shock.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.isel(t=ti)).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/tests/blob/blob.py b/setups/tests/blob/blob.py deleted file mode 100644 index 77337d3b2..000000000 --- a/setups/tests/blob/blob.py +++ /dev/null @@ -1,62 +0,0 @@ -import h5py -import numpy as np -import matplotlib.pyplot as plt - -f = open("report", "r") -Lines = f.readlines() -f.close() - -em_new = [] -ep_new = [] -time_new = [] -for i in range (len(Lines)): - line = Lines[i] - line = line.strip() - arr = line.split() - - if (len(arr)>0 and arr[0]=='species'): - nparts = arr[2].split("..") - if (nparts[0]=="(e-_p)"): - em_new.append(float(nparts[-1])) - if (nparts[0]=="(e+_p)"): - ep_new.append(float(nparts[-1])) - - if (len(arr)>0 and arr[0]=='Time:'): - time_new.append(float(arr[1])) - -f = h5py.File('blob.h5', 'r') - -Nsteps = len(f.keys()) -print(list(f['Step0'].keys())) - -for i in range (Nsteps): - print (i) - fig = plt.figure(dpi=300, figsize=(8,8), facecolor='white') - - densMax = max(np.max(f['Step'+str(i)]['fN_1']),np.max(f['Step'+str(i)]['fN_2'])) - print(densMax) - ax1 = fig.add_axes([0.05,0.05,0.4,0.4]) - im1=ax1.pcolormesh(f['Step'+str(i)]['X1'],f['Step'+str(i)]['X2'],f['Step'+str(i)]['fN_1'],cmap='turbo',vmin=0,vmax=1.0) - ax1.set_title(r"$N_1$") - ax1.vlines(0,-10.0,10.0,color='white') - - ax1 = fig.add_axes([0.48,0.05,0.4,0.4]) - ax1.pcolormesh(f['Step'+str(i)]['X1'],f['Step'+str(i)]['X2'],f['Step'+str(i)]['fN_2'],cmap='turbo',vmin=0,vmax=1.0) - ax1.set_yticklabels([]) - ax1.set_title(r"$N_2$") - ax1.vlines(0,-10.0,10.0,color='white') - - ax4cb = fig.add_axes([0.89, 0.05, 0.01, 0.4]) - cbar4 = fig.colorbar(im1,cax=ax4cb) - - ax1= fig.add_axes([0.05,0.5,0.83,0.4]) - ax1.plot(time_new,em_new, color='blue', label=r'$e^-$, new') - ax1.plot(time_new,ep_new, color='red', label=r'$e^+$, new') - ax1.legend() - ax1.set_ylim(0,1.8e5) - ax1.set_xlim(0,100) - ax1.vlines(i, 0,1.8e5, color='green',linewidth=0.6) - - - fig.savefig("%05d"%i+".png",dpi=300,bbox_inches='tight') - plt.close() diff --git a/setups/tests/blob/blob.toml b/setups/tests/blob/blob.toml deleted file mode 100644 index 7a047f348..000000000 --- a/setups/tests/blob/blob.toml +++ /dev/null @@ -1,66 +0,0 @@ -[simulation] - name = "blob" - engine = "srpic" - runtime = 100.0 - - [simulation.domain] - decomposition = [2, 1, 1] - -[grid] - resolution = [1024, 1024] - extent = [[-10.0, 10.0], [-10.0, 10.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 1.0 - skindepth0 = 1.0 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 16.0 - - [[particles.species]] - label = "e-_p" - mass = 1.0 - charge = -1.0 - maxnpart = 1e7 - - [[particles.species]] - label = "e+_p" - mass = 1.0 - charge = 1.0 - maxnpart = 1e7 - -[setup] - temp_1 = 1e-4 - x1c = -5.0 - x2c = 0.0 - v_max = 50.0 - dr = 1.0 - -[output] - format = "hdf5" - interval_time = 1.0 - - [output.fields] - quantities = ["N_1", "N_2", "B", "E"] - - [output.particles] - enable = false - - [output.spectra] - enable = false - -[diagnostics] - colored_stdout = false diff --git a/setups/tests/blob/nparts.py b/setups/tests/blob/nparts.py deleted file mode 100644 index e759422c0..000000000 --- a/setups/tests/blob/nparts.py +++ /dev/null @@ -1,38 +0,0 @@ -import h5py -import numpy as np -import matplotlib.pyplot as plt - -f = open("report", "r") -Lines = f.readlines() -f.close() - -em_new = [] -ep_new = [] -time_new = [] -for i in range (len(Lines)): - line = Lines[i] - line = line.strip() - arr = line.split() - - if (len(arr)>0 and arr[0]=='species'): - nparts = arr[2].split("..") - if (nparts[0]=="(e-_p)"): - em_new.append(float(nparts[-1])) - if (nparts[0]=="(e+_p)"): - ep_new.append(float(nparts[-1])) - - if (len(arr)>0 and arr[0]=='Time:'): - time_new.append(float(arr[1])) - - -fig = plt.figure(dpi=300, figsize=(8,8), facecolor='white') - -ax1= fig.add_axes([0.05,0.5,0.83,0.4]) -ax1.plot(time_new,em_new, color='blue', label=r'$e^-$, new') -ax1.plot(time_new,ep_new, color='red', label=r'$e^+$, new') -ax1.legend() -ax1.set_ylim(0,1.8e5) -ax1.set_xlim(0,100) - -fig.savefig("nparts.png",dpi=300,bbox_inches='tight') -plt.close() diff --git a/setups/tests/blob/pgen.hpp b/setups/tests/blob/pgen.hpp deleted file mode 100644 index 9120e244f..000000000 --- a/setups/tests/blob/pgen.hpp +++ /dev/null @@ -1,103 +0,0 @@ -#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/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/domain.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct CounterstreamEnergyDist : public arch::EnergyDistribution { - CounterstreamEnergyDist(const M& metric, real_t v_max) - : arch::EnergyDistribution { metric } - , v_max { v_max } {} - - Inline void operator()(const coord_t& x_Ph, - vec_t& v, - spidx_t sp) const override { - v[0] = v_max; - } - - private: - const real_t v_max; - }; - - template - struct GaussianDist : public arch::SpatialDistribution { - GaussianDist(const M& metric, real_t x1c, real_t x2c, real_t dr) - : arch::SpatialDistribution { metric } - , x1c { x1c } - , x2c { x2c } - , dr { dr } {} - - // to properly scale the number density, the probability should be normalized to 1 - Inline auto operator()(const coord_t& x_Ph) const -> real_t override { - if (math::abs(x_Ph[0] - x1c) < dr && math::abs(x_Ph[1] - x2c) < dr) { - return 1.0; - } else { - return 0.0; - } - } - - private: - const real_t x1c, x2c, dr; - }; - - 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; - - const real_t temp_1, x1c, x2c, dr, v_max; - - inline PGen(const SimulationParams& p, const Metadomain& global_domain) - : arch::ProblemGenerator { p } - , temp_1 { p.template get("setup.temp_1") } - , x1c { p.template get("setup.x1c") } - , x2c { p.template get("setup.x2c") } - , v_max { p.template get("setup.v_max") } - , dr { p.template get("setup.dr") } {} - - inline void InitPrtls(Domain& local_domain) { - const auto energy_dist = CounterstreamEnergyDist(local_domain.mesh.metric, - v_max); - const auto spatial_dist = GaussianDist(local_domain.mesh.metric, - x1c, - x2c, - dr); - const auto injector = - arch::NonUniformInjector( - energy_dist, - spatial_dist, - { 1, 2 }); - - arch::InjectNonUniform>( - params, - local_domain, - injector, - 1.0); - } - }; - -} // namespace user - -#endif diff --git a/setups/tests/customout/customout.toml b/setups/tests/customout/customout.toml deleted file mode 100644 index 497b96dc2..000000000 --- a/setups/tests/customout/customout.toml +++ /dev/null @@ -1,50 +0,0 @@ -[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 deleted file mode 100644 index 22c8f6564..000000000 --- a/setups/tests/customout/pgen.hpp +++ /dev/null @@ -1,86 +0,0 @@ -#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/setups/tests/deposit/deposit-mink.toml b/setups/tests/deposit/deposit-mink.toml deleted file mode 100644 index 2dc953896..000000000 --- a/setups/tests/deposit/deposit-mink.toml +++ /dev/null @@ -1,71 +0,0 @@ -[simulation] - name = "deposit-test-mink" - engine = "srpic" - runtime = 5.0 - -[grid] - resolution = [32, 32] - extent = [[0.0, 1.0], [0.0, 1.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.1 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 10.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e2 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e2 - -[setup] - x1_e = [0.25] - x2_e = [0.85] - x3_e = [0.33] - ux1_e = [0.6] - ux2_e = [-0.3] - ux3_e = [-0.2] - - x1_i = [0.25] - x2_i = [0.85] - x3_i = [0.33] - ux1_i = [-0.6] - ux2_i = [0.3] - ux3_i = [0.2] - -[output] - format = "hdf5" - interval = 5 - - [output.fields] - quantities = ["N_1", "N_2", "E", "B", "J"] - - [output.particles] - enable = false - - [output.spectra] - enable = false - -[checkpoint] - keep = 0 diff --git a/setups/tests/deposit/deposit-sr.toml b/setups/tests/deposit/deposit-sr.toml deleted file mode 100644 index 0e1648d12..000000000 --- a/setups/tests/deposit/deposit-sr.toml +++ /dev/null @@ -1,71 +0,0 @@ -[simulation] - name = "deposit-sr" - engine = "srpic" - runtime = 10.0 - -[grid] - resolution = [64, 64] - extent = [[1.0, 5.0]] - - [grid.metric] - metric = "qspherical" - - [grid.boundaries] - fields = [["FIXED", "FIXED"]] - particles = [["REFLECT", "REFLECT"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.1 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 10.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e2 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e2 - -[setup] - x1_e = [2.25] - x2_e = [1.25] - phi_e = [0.0] - ux1_e = [0.6] - ux2_e = [-0.3] - ux3_e = [-0.2] - - x1_i = [2.25] - x2_i = [1.25] - phi_i = [0.0] - ux1_i = [-0.6] - ux2_i = [0.3] - ux3_i = [0.2] - -[output] - format = "hdf5" - interval = 5 - - [output.fields] - quantities = ["N_1", "N_2", "E", "B", "J"] - - [output.particles] - enable = false - - [output.spectra] - enable = false - -[checkpoint] - keep = 0 diff --git a/setups/tests/deposit/pgen.hpp b/setups/tests/deposit/pgen.hpp deleted file mode 100644 index 0080af8fd..000000000 --- a/setups/tests/deposit/pgen.hpp +++ /dev/null @@ -1,113 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/traits.h" - -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/domain.h" -#include "framework/domain/metadomain.h" - -#include - -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; - - const Metadomain& global_domain; - - inline PGen(const SimulationParams& p, const Metadomain& global_domain) - : arch::ProblemGenerator { p } - , global_domain { global_domain } {} - - inline void InitPrtls(Domain& local_domain) { - const auto empty = std::vector {}; - const auto x1_e = params.template get>("setup.x1_e", - empty); - const auto x2_e = params.template get>("setup.x2_e", - empty); - const auto x3_e = params.template get>("setup.x3_e", - empty); - const auto phi_e = params.template get>("setup.phi_e", - empty); - const auto ux1_e = params.template get>("setup.ux1_e", - empty); - const auto ux2_e = params.template get>("setup.ux2_e", - empty); - const auto ux3_e = params.template get>("setup.ux3_e", - empty); - - const auto x1_i = params.template get>("setup.x1_i", - empty); - const auto x2_i = params.template get>("setup.x2_i", - empty); - const auto x3_i = params.template get>("setup.x3_i", - empty); - const auto phi_i = params.template get>("setup.phi_i", - empty); - const auto ux1_i = params.template get>("setup.ux1_i", - empty); - const auto ux2_i = params.template get>("setup.ux2_i", - empty); - const auto ux3_i = params.template get>("setup.ux3_i", - empty); - std::map> data_e { - { "x1", x1_e }, - { "x2", x2_e }, - { "ux1", ux1_e }, - { "ux2", ux2_e }, - { "ux3", ux3_e } - }; - std::map> data_i { - { "x1", x1_i }, - { "x2", x2_i }, - { "ux1", ux1_i }, - { "ux2", ux2_i }, - { "ux3", ux3_i } - }; - if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { - data_e["x3"] = x3_e; - data_i["x3"] = x3_i; - } else if constexpr (D == Dim::_2D) { - data_e["phi"] = phi_e; - data_i["phi"] = phi_i; - } - - arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_e); - arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_i); - } - - auto FixFieldsConst(const bc_in&, const em&) const -> std::pair { - return { ZERO, false }; - } - }; - -} // namespace user - -#endif diff --git a/setups/tests/deposit/plot-mink.py b/setups/tests/deposit/plot-mink.py deleted file mode 100644 index 9d6760613..000000000 --- a/setups/tests/deposit/plot-mink.py +++ /dev/null @@ -1,62 +0,0 @@ -import nt2 -import matplotlib.pyplot as plt -import matplotlib as mpl - -datas = [] -cpus = [1, 2, 3, 4, 5, 6, 8] -for i in cpus: - datas.append(nt2.Data(path=f"mink-np{i}")) - - -def plot(ti): - fig = plt.figure(figsize=(16, 7), dpi=300) - gs = mpl.gridspec.GridSpec(3, 7, figure=fig) - - for p, quant in enumerate(["Jx", "Jy", "Jz"]): - axs = [fig.add_subplot(gs[p, i]) for i in range(7)] - (datas[0].fields[quant]).isel(t=ti).plot( - ax=axs[0], - cmap="seismic", - add_colorbar=False, - norm=mpl.colors.SymLogNorm( - linthresh=1e-5, - linscale=1, - vmin=-1e-2, - vmax=1e-2, - ), - ) - for i, (d, ax) in enumerate(zip(datas[1:], axs[1:])): - (d.fields[quant] - datas[0].fields[quant]).isel(t=ti).plot( - ax=ax, - cmap="seismic", - add_colorbar=False, - norm=mpl.colors.SymLogNorm( - linthresh=1e-10, - linscale=1, - vmin=-1e-7, - vmax=1e-7, - ), - ) - - for i, ax in enumerate(axs): - ax.set_aspect(1) - if i > 0: - if p == 0: - ax.set_title(f"np{cpus[i]} - np1") - else: - ax.set_title(None) - ax.set_yticklabels([]) - ax.set_ylabel(None) - else: - if p == 0: - ax.set_title(f"np1") - else: - ax.set_title(None) - - if p != 2: - ax.set_xticklabels([]) - ax.set_xlabel(None) - - -nt2.export.makeFrames(plot, datas[0].fields.s[::4], "mink-diff", num_cpus=4) -nt2.export.makeMovie(framerate=10, input="mink-diff/", number=5, output="mink-diff.mp4") diff --git a/setups/tests/deposit/plot-sr.py b/setups/tests/deposit/plot-sr.py deleted file mode 100644 index 0ba5dff24..000000000 --- a/setups/tests/deposit/plot-sr.py +++ /dev/null @@ -1,29 +0,0 @@ -import nt2 -import matplotlib.pyplot as plt -import matplotlib as mpl - -data = nt2.Data(path=f"sr-np8") - - -def plot(ti): - fig = plt.figure(figsize=(9, 6), dpi=300) - gs = mpl.gridspec.GridSpec(1, 3, figure=fig) - axs = [fig.add_subplot(gs[0, i]) for i in range(3)] - for i, (ax, j) in enumerate(zip(axs, ["Jr", "Jth", "Jph"])): - data.fields.isel(t=ti)[j].polar.pcolor( - ax=ax, - cbar_position="top", - cbar_size="2%", - norm=mpl.colors.SymLogNorm(linthresh=1e-8, vmin=-1e-4, vmax=1e-4), - cmap="seismic", - ) - ax.set_title(None) - ax.add_artist(mpl.patches.Circle((0, 0), 1, color="k", alpha=0.2)) - ax.add_artist(mpl.patches.Circle((0, 0), 5, edgecolor="k", facecolor="none")) - if i > 0: - ax.set_yticklabels([]) - ax.set_ylabel(None) - - -nt2.export.makeFrames(plot, data.fields.s, "sr-dep", num_cpus=4) -nt2.export.makeMovie(framerate=10, input="sr-dep/", number=5, output="sr-dep.mp4") diff --git a/setups/tests/deposit/run-mink.sh b/setups/tests/deposit/run-mink.sh deleted file mode 100755 index 4b52d6642..000000000 --- a/setups/tests/deposit/run-mink.sh +++ /dev/null @@ -1,14 +0,0 @@ -#!/usr/bin/env bash - -for i in {1..8}; do - if [ $i -eq 7 ]; then - continue - fi - run=$(echo "np${i}") - cp deposit-mink.toml deposit-mink-${run}.toml && \ - sed -i 's/name[[:space:]]*=[[:space:]]*".*\?"/name = "mink-'${run}'"/g' deposit-mink-${run}.toml && \ - mpiexec -np ${i} ./entity.xc -input deposit-mink-${run}.toml && \ - rm deposit-mink-${run}.toml -done - -rm *.info *.err *.log *.csv diff --git a/setups/tests/injector/injector.toml b/setups/tests/injector/injector.toml deleted file mode 100644 index 10fdaa251..000000000 --- a/setups/tests/injector/injector.toml +++ /dev/null @@ -1,62 +0,0 @@ -[simulation] - name = "injector-test" - engine = "srpic" - runtime = 2.0 - -[grid] - resolution = [512, 512] - extent = [[-1.0, 1.0], [-1.0, 1.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["ABSORB", "ABSORB"], ["ABSORB", "ABSORB"]] - particles = [["ABSORB", "ABSORB"], ["ABSORB", "ABSORB"]] - - [grid.boundaries.absorb] - ds = 0.15 - -[scales] - larmor0 = 0.1 - skindepth0 = 0.1 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 1.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e6 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e6 - -[setup] - period = 0.1 - vmax = 1.0 - x1c = 0.25 - x2c = -0.32 - dr = 1e-2 - rate = 0.1 - -[output] - format = "hdf5" - interval_time = 0.01 - - [output.fields] - quantities = ["N_1", "N_2", "E"] - -[diagnostics] - interval = 10 - colored_stdout = true diff --git a/setups/tests/injector/pgen.hpp b/setups/tests/injector/pgen.hpp deleted file mode 100644 index 0b4a34a07..000000000 --- a/setups/tests/injector/pgen.hpp +++ /dev/null @@ -1,103 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/numeric.h" - -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "archetypes/spatial_dist.h" -#include "framework/domain/domain.h" -#include "framework/domain/metadomain.h" - -namespace user { - using namespace ntt; - - template - struct Firehose : public arch::EnergyDistribution { - Firehose(const M& metric, real_t time, real_t period, real_t vmax) - : arch::EnergyDistribution { metric } - , phase { (real_t)(constant::TWO_PI)*time / period } - , vmax { vmax } {} - - Inline void operator()(const coord_t&, - vec_t& v_Ph, - spidx_t) const override { - v_Ph[0] = vmax * math::cos(phase); - v_Ph[1] = vmax * math::sin(phase); - } - - private: - const real_t phase, vmax; - }; - - template - struct PointDistribution : public arch::SpatialDistribution { - PointDistribution(const M& metric, real_t x1c, real_t x2c, real_t dr) - : arch::SpatialDistribution { metric } - , x1c { x1c } - , x2c { x2c } - , dr { dr } {} - - Inline auto operator()(const coord_t& x_Ph) const -> real_t override { - return math::exp(-(SQR(x_Ph[0] - x1c) + SQR(x_Ph[1] - x2c)) / SQR(dr)); - } - - private: - const real_t x1c, x2c, dr; - }; - - 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; - - const real_t period, vmax, x1c, x2c, dr, rate; - - inline PGen(const SimulationParams& p, const Metadomain&) - : arch::ProblemGenerator { p } - , period { params.template get("setup.period", 1.0) } - , vmax { params.template get("setup.vmax", 1.0) } - , x1c { params.template get("setup.x1c", 0.0) } - , x2c { params.template get("setup.x2c", 0.0) } - , dr { params.template get("setup.dr", 0.1) } - , rate { params.template get("setup.rate", 1.0) } {} - - void CustomPostStep(std::size_t, long double time, Domain& domain) { - const auto energy_dist = Firehose(domain.mesh.metric, - (real_t)time, - period, - vmax); - const auto spatial_dist = PointDistribution(domain.mesh.metric, - x1c, - x2c, - dr); - const auto injector = arch::NonUniformInjector( - energy_dist, - spatial_dist, - { 1, 2 }); - - arch::InjectNonUniform>( - params, - domain, - injector, - rate); - } - }; - -} // namespace user - -#endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 715183b64..a54a119b5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -37,8 +37,6 @@ if(${output}) add_subdirectory(${SRC_DIR}/checkpoint ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) endif() -add_subdirectory(${SRC_DIR}/../setups ${CMAKE_CURRENT_BINARY_DIR}/setups) - set(ENTITY ${PROJECT_NAME}.xc) set(SOURCES ${SRC_DIR}/entity.cpp) diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h index 67f1bdf8a..e9ef8b0b7 100644 --- a/src/archetypes/energy_dist.h +++ b/src/archetypes/energy_dist.h @@ -184,12 +184,13 @@ namespace arch { } template - Inline void SampleFromMaxwellian(vec_t& v, - const real_t& temperature, - const real_t& boost_velocity, - const in& boost_direction, - bool flip_velocity, - const random_number_pool_t& pool) { + Inline void SampleFromMaxwellian( + vec_t& v, + const random_number_pool_t& pool, + const real_t& temperature, + const real_t& boost_velocity = static_cast(0), + const in& boost_direction = in::x1, + bool flip_velocity = false) { if (cmp::AlmostZero(temperature)) { v[0] = ZERO; v[1] = ZERO; @@ -242,7 +243,7 @@ namespace arch { "Maxwellian: Temperature must be non-negative", HERE); raise::ErrorIf( - (not cmp::AlmostZero(boost_vel, ZERO)) && (M::CoordType != Coord::Cart), + (not cmp::AlmostZero_host(boost_vel, ZERO)) && (M::CoordType != Coord::Cart), "Maxwellian: Boosting is only supported in Cartesian coordinates", HERE); } @@ -251,12 +252,12 @@ namespace arch { vec_t& v, spidx_t sp = 0) const override { SampleFromMaxwellian(v, + pool, temperature, boost_velocity, boost_direction, not zero_current and - sp % 2 == 0, - pool); + sp % 2 == 0); if constexpr (S == SimEngine::GRPIC) { // convert from the tetrad basis to covariant vec_t v_Hat; @@ -312,11 +313,11 @@ namespace arch { spidx_t sp = 0) const override { SampleFromMaxwellian( v, + pool, (sp == sp_1) ? temperature_1 : temperature_2, boost_velocity, boost_direction, - not zero_current and sp == sp_1, - pool); + not zero_current and sp == sp_1); if constexpr (S == SimEngine::GRPIC) { // convert from the tetrad basis to covariant vec_t v_Hat;