From ce35970ceb0400d3e0a238e9f3c045ddc2114e80 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 6 Mar 2023 16:49:39 -0500 Subject: [PATCH 1/5] =?UTF-8?q?=E2=9C=A8=20nppc=20as=20output?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/definitions.h | 5 ++++- src/framework/io/output.cpp | 12 +++++++++- src/framework/meshblock/meshblock.cpp | 32 +++++++++++++++++---------- 3 files changed, 35 insertions(+), 14 deletions(-) diff --git a/src/definitions.h b/src/definitions.h index b5b89fffb..78c5c6134 100644 --- a/src/definitions.h +++ b/src/definitions.h @@ -150,7 +150,8 @@ namespace ntt { J, // Current density T, // Particle distribution moments Rho, // Particle density - N // Particle number density + N, // Particle number density + Nppc // Raw number of particles per each cell }; inline auto StringizeFieldID(const FieldID& id) -> std::string { @@ -167,6 +168,8 @@ namespace ntt { return "Rho"; case FieldID::N: return "N"; + case FieldID::Nppc: + return "Nppc"; default: return "UNKNOWN"; } diff --git a/src/framework/io/output.cpp b/src/framework/io/output.cpp index 86b42fc40..e59835269 100644 --- a/src/framework/io/output.cpp +++ b/src/framework/io/output.cpp @@ -118,6 +118,10 @@ namespace ntt { const auto id = FieldID::Rho; auto species = InterpretInputField_getspecies(fld); return InterpretInputField_helper(id, { {} }, species); + } else if (fld.find("Nppc") == 0) { + const auto id = FieldID::Nppc; + auto species = InterpretInputField_getspecies(fld); + return InterpretInputField_helper(id, { {} }, species); } else if (fld.find("N") == 0) { const auto id = FieldID::N; auto species = InterpretInputField_getspecies(fld); @@ -214,7 +218,13 @@ namespace ntt { } } else { for (std::size_t i { 0 }; i < comp.size(); ++i) { - mblock.ComputeMoments(params, m_id, comp[i], species, i % 3, params.outputMomSmooth()); + // no smoothing for FieldID::Nppc + mblock.ComputeMoments(params, + m_id, + comp[i], + species, + i % 3, + m_id == FieldID::Nppc ? 0 : params.outputMomSmooth()); PutField(io, writer, name(i), mblock.buff, i % 3); } } diff --git a/src/framework/meshblock/meshblock.cpp b/src/framework/meshblock/meshblock.cpp index 97115fd21..bfaac3531 100644 --- a/src/framework/meshblock/meshblock.cpp +++ b/src/framework/meshblock/meshblock.cpp @@ -82,7 +82,10 @@ namespace ntt { std::vector A = { this->buff_content[buff_ind] }; AssertEmptyContent(A); std::size_t ni1 = this->Ni1(), ni2 = this->Ni2(), ni3 = this->Ni3(); - real_t weight = (1.0 / params.ppc0()) / math::pow(2.0 * smooth + 1.0, static_cast(D)); + real_t weight = ONE / math::pow(2.0 * smooth + 1.0, static_cast(D)); + if (field != FieldID::Nppc) { + weight /= params.ppc0(); + } if constexpr (D == Dim1) { Kokkos::deep_copy(Kokkos::subview(this->buff, Kokkos::ALL(), (int)buff_ind), ZERO); } else if constexpr (D == Dim2) { @@ -131,7 +134,7 @@ namespace ntt { real_t contrib { ZERO }; if (field == FieldID::Rho) { contrib = ((mass == ZERO) ? ONE : mass); - } else if (field == FieldID::N) { + } else if ((field == FieldID::N) || (field == FieldID::Nppc)) { contrib = ONE; } else if (field == FieldID::T) { real_t energy { ((mass == ZERO) ? get_photon_Energy_SR(species, p) @@ -149,9 +152,11 @@ namespace ntt { } } } + if (fieldID != FieldID::Nppc) { + contrib *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1 }); + } for (int i1_ = i1_min; i1_ <= i1_max; ++i1_) { - buff_access(i1_, buff_ind) += contrib * weight * this_metric.min_cell_volume() - / this_metric.sqrt_det_h({ x1 }); + buff_access(i1_, buff_ind) += contrib * weight; } } }); @@ -171,7 +176,7 @@ namespace ntt { real_t contrib { ZERO }; if (field == FieldID::Rho) { contrib = ((mass == ZERO) ? ONE : mass); - } else if (field == FieldID::N) { + } else if ((field == FieldID::N) || (field == FieldID::Nppc)) { contrib = ONE; } else if (field == FieldID::T) { real_t energy { ((mass == ZERO) ? get_photon_Energy_SR(species, p) @@ -204,11 +209,12 @@ namespace ntt { } #endif } + if (fieldID != FieldID::Nppc) { + contrib *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1, x2 }); + } for (int i2_ = i2_min; i2_ <= i2_max; ++i2_) { for (int i1_ = i1_min; i1_ <= i1_max; ++i1_) { - buff_access(i1_, i2_, buff_ind) += contrib * weight - * this_metric.min_cell_volume() - / this_metric.sqrt_det_h({ x1, x2 }); + buff_access(i1_, i2_, buff_ind) += contrib * weight; } } } @@ -233,7 +239,7 @@ namespace ntt { real_t contrib { ZERO }; if (field == FieldID::Rho) { contrib = ((mass == ZERO) ? ONE : mass); - } else if (field == FieldID::N) { + } else if ((field == FieldID::N) || (field == FieldID::Nppc)) { contrib = ONE; } else if (field == FieldID::T) { real_t energy { ((mass == ZERO) ? get_photon_Energy_SR(species, p) @@ -264,12 +270,14 @@ namespace ntt { } #endif } + if (fieldID != FieldID::Nppc) { + contrib + *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1, x2, x3 }); + } for (int i3_ = i3_min; i3_ <= i3_max; ++i3_) { for (int i2_ = i2_min; i2_ <= i2_max; ++i2_) { for (int i1_ = i1_min; i1_ <= i1_max; ++i1_) { - buff_access(i1_, i2_, i3_, buff_ind) - += contrib * weight * this_metric.min_cell_volume() - / this_metric.sqrt_det_h({ x1, x2, x3 }); + buff_access(i1_, i2_, i3_, buff_ind) += contrib * weight; } } } From 820fa16af8f2483ffad8836455f4649b9d3aa9dd Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 6 Mar 2023 18:02:58 -0500 Subject: [PATCH 2/5] =?UTF-8?q?=F0=9F=9A=91=20random=20template=20added?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/wrapper/kokkos.h.in | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/src/wrapper/kokkos.h.in b/src/wrapper/kokkos.h.in index 306191fa2..f28727270 100644 --- a/src/wrapper/kokkos.h.in +++ b/src/wrapper/kokkos.h.in @@ -103,8 +103,22 @@ namespace ntt { Kokkos::MDRangePolicy, HostExeSpace>, std::nullptr_t>::type>::type>::type; - // Random number pool alias + // Random number pool/generator type alias using RandomNumberPool_t = Kokkos::Random_XorShift1024_Pool; + using RandomGenerator_t = typename RandomNumberPool_t::generator_type; + + template + Inline auto Random(RandomGenerator_t&) -> T; + + template <> + Inline auto Random(RandomGenerator_t& gen) -> float { + return gen.frand(); + } + + template <> + Inline auto Random(RandomGenerator_t& gen) -> double { + return gen.drand(); + } /** * @brief Function template for generating ND Kokkos range policy. From 6b93c373226ae0a46e6e17fe841c4a4a201b8e90 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 6 Mar 2023 19:24:49 -0500 Subject: [PATCH 3/5] =?UTF-8?q?=E2=9C=A8=20weights=20support=20added?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- input.toml.example | 1 + src/definitions.h | 2 + src/framework/meshblock/meshblock.cpp | 16 +++- src/framework/meshblock/particles.cpp | 2 +- src/framework/meshblock/particles.h | 4 +- src/framework/particle_macros.h | 133 +++++++++++++------------- src/framework/sim_params.cpp | 1 + src/framework/sim_params.h | 30 +++--- src/pic/currents/deposit_currents.cpp | 4 +- src/pic/currents/deposit_currents.hpp | 61 +++++++----- 10 files changed, 147 insertions(+), 107 deletions(-) diff --git a/input.toml.example b/input.toml.example index 19d624dbf..9209a34da 100644 --- a/input.toml.example +++ b/input.toml.example @@ -27,6 +27,7 @@ skindepth0 = n_species = shuffle_step = max_dead_frac = +use_weights = [species_] label = diff --git a/src/definitions.h b/src/definitions.h index 78c5c6134..c5f425e5f 100644 --- a/src/definitions.h +++ b/src/definitions.h @@ -135,6 +135,8 @@ namespace ntt { const unsigned short current_filters = 0; + const bool use_weights = false; + const int shuffle_interval = 0; const double max_dead_frac = 0.0; diff --git a/src/framework/meshblock/meshblock.cpp b/src/framework/meshblock/meshblock.cpp index bfaac3531..a4c63d15c 100644 --- a/src/framework/meshblock/meshblock.cpp +++ b/src/framework/meshblock/meshblock.cpp @@ -117,6 +117,7 @@ namespace ntt { } auto this_metric = this->metric; + auto use_weights = params.useWeights(); auto scatter_buff = Kokkos::Experimental::create_scatter_view(this->buff); for (auto& sp : out_species) { @@ -152,8 +153,11 @@ namespace ntt { } } } - if (fieldID != FieldID::Nppc) { + if (field != FieldID::Nppc) { contrib *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1 }); + if (use_weights) { + contrib *= species.weight(p); + } } for (int i1_ = i1_min; i1_ <= i1_max; ++i1_) { buff_access(i1_, buff_ind) += contrib * weight; @@ -209,8 +213,11 @@ namespace ntt { } #endif } - if (fieldID != FieldID::Nppc) { + if (field != FieldID::Nppc) { contrib *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1, x2 }); + if (use_weights) { + contrib *= species.weight(p); + } } for (int i2_ = i2_min; i2_ <= i2_max; ++i2_) { for (int i1_ = i1_min; i1_ <= i1_max; ++i1_) { @@ -270,9 +277,12 @@ namespace ntt { } #endif } - if (fieldID != FieldID::Nppc) { + if (field != FieldID::Nppc) { contrib *= this_metric.min_cell_volume() / this_metric.sqrt_det_h({ x1, x2, x3 }); + if (use_weights) { + contrib *= species.weight(p); + } } for (int i3_ = i3_min; i3_ <= i3_max; ++i3_) { for (int i2_ = i2_min; i2_ <= i2_max; ++i2_) { diff --git a/src/framework/meshblock/particles.cpp b/src/framework/meshblock/particles.cpp index 18bbcd324..7c9ada639 100644 --- a/src/framework/meshblock/particles.cpp +++ b/src/framework/meshblock/particles.cpp @@ -259,7 +259,7 @@ namespace ntt { Sorter.sort(Kokkos::subview(dx3_prev, slice)); } } - // !TODO: sort weights + Sorter.sort(Kokkos::subview(weight, slice)); } } // namespace ntt \ No newline at end of file diff --git a/src/framework/meshblock/particles.h b/src/framework/meshblock/particles.h index 4675308f8..816dcecef 100644 --- a/src/framework/meshblock/particles.h +++ b/src/framework/meshblock/particles.h @@ -33,11 +33,11 @@ namespace ntt { // Three spatial components of the covariant 4-velocity (physical units). array_t ux1, ux2, ux3; // Particle weights. - array_t weight; + array_t weight; // Additional variables (specific to different cases). // previous coordinates (GR specific) array_t i1_prev, i2_prev, i3_prev; - array_t dx1_prev, dx2_prev, dx3_prev; + array_t dx1_prev, dx2_prev, dx3_prev; // phi coordinate (for axisymmetry) array_t phi; // Array to tag the particles. diff --git a/src/framework/particle_macros.h b/src/framework/particle_macros.h index 706b62a0f..dbd606463 100644 --- a/src/framework/particle_macros.h +++ b/src/framework/particle_macros.h @@ -28,105 +28,108 @@ #define get_photon_Energy_SR(PARTICLES, P) (math::sqrt(get_prtl_Usqr_SR((PARTICLES), (P)))) +#define init_prtl_1d_i_di(SPECIES, INDEX, I1, DI1, U1, U2, U3, WEIGHT) \ + { \ + (SPECIES).i1((INDEX)) = I1; \ + (SPECIES).dx1((INDEX)) = DI1; \ + (SPECIES).ux1((INDEX)) = U1; \ + (SPECIES).ux2((INDEX)) = U2; \ + (SPECIES).ux3((INDEX)) = U3; \ + (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + (SPECIES).weight((INDEX)) = WEIGHT; \ + } + +#define init_prtl_2d_i_di(SPECIES, INDEX, I1, I2, DI1, DI2, U1, U2, U3, WEIGHT) \ + { \ + (SPECIES).i1((INDEX)) = I1; \ + (SPECIES).dx1((INDEX)) = DI1; \ + (SPECIES).i2((INDEX)) = I2; \ + (SPECIES).dx2((INDEX)) = DI2; \ + (SPECIES).ux1((INDEX)) = U1; \ + (SPECIES).ux2((INDEX)) = U2; \ + (SPECIES).ux3((INDEX)) = U3; \ + (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + (SPECIES).weight((INDEX)) = WEIGHT; \ + } + +#define init_prtl_3d_i_di(SPECIES, INDEX, I1, I2, I3, DI1, DI2, DI3, U1, U2, U3, WEIGHT) \ + { \ + (SPECIES).i1((INDEX)) = I1; \ + (SPECIES).dx1((INDEX)) = DI1; \ + (SPECIES).i2((INDEX)) = I2; \ + (SPECIES).dx2((INDEX)) = DI2; \ + (SPECIES).i3((INDEX)) = I3; \ + (SPECIES).dx3((INDEX)) = DI3; \ + (SPECIES).ux1((INDEX)) = U1; \ + (SPECIES).ux2((INDEX)) = U2; \ + (SPECIES).ux3((INDEX)) = U3; \ + (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + (SPECIES).weight((INDEX)) = WEIGHT; \ + } + #ifdef MINKOWSKI_METRIC -# define init_prtl_1d(MBLOCK, SPECIES, INDEX, X1, U1, U2, U3) \ +# define init_prtl_1d(MBLOCK, SPECIES, INDEX, X1, U1, U2, U3, WEIGHT) \ { \ coord_t X_CU; \ int I; \ float DX; \ ((MBLOCK).metric).x_Cart2Code({ (X1) }, X_CU); \ from_Xi_to_i_di(X_CU[0], I, DX); \ - (SPECIES).i1((INDEX)) = I; \ - (SPECIES).dx1((INDEX)) = DX; \ - (SPECIES).ux1((INDEX)) = U1; \ - (SPECIES).ux2((INDEX)) = U2; \ - (SPECIES).ux3((INDEX)) = U3; \ - (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + init_prtl_1d_i_di(SPECIES, INDEX, I, DX, U1, U2, U3, WEIGHT); \ } -# define init_prtl_2d(MBLOCK, SPECIES, INDEX, X1, X2, U1, U2, U3) \ +# define init_prtl_2d(MBLOCK, SPECIES, INDEX, X1, X2, U1, U2, U3, WEIGHT) \ { \ coord_t X_CU; \ - int I; \ - float DX; \ + int I1, I2; \ + float DX1, DX2; \ ((MBLOCK).metric).x_Cart2Code({ (X1), (X2) }, X_CU); \ - from_Xi_to_i_di(X_CU[0], I, DX); \ - (SPECIES).i1((INDEX)) = I; \ - (SPECIES).dx1((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[1], I, DX); \ - (SPECIES).i2((INDEX)) = I; \ - (SPECIES).dx2((INDEX)) = DX; \ - (SPECIES).ux1((INDEX)) = U1; \ - (SPECIES).ux2((INDEX)) = U2; \ - (SPECIES).ux3((INDEX)) = U3; \ - (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + from_Xi_to_i_di(X_CU[0], I1, DX1); \ + from_Xi_to_i_di(X_CU[1], I2, DX2); \ + init_prtl_2d_i_di(SPECIES, INDEX, I1, I2, DX1, DX2, U1, U2, U3, WEIGHT); \ } -# define init_prtl_3d(MBLOCK, SPECIES, INDEX, X1, X2, X3, U1, U2, U3) \ +# define init_prtl_3d(MBLOCK, SPECIES, INDEX, X1, X2, X3, U1, U2, U3, WEIGHT) \ { \ coord_t X_CU; \ - int I; \ - float DX; \ + int I1, I2, I3; \ + float DX1, DX2, DX3; \ ((MBLOCK).metric).x_Cart2Code({ (X1), (X2), (X3) }, X_CU); \ - from_Xi_to_i_di(X_CU[0], I, DX); \ - (SPECIES).i1((INDEX)) = I; \ - (SPECIES).dx1((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[1], I, DX); \ - (SPECIES).i2((INDEX)) = I; \ - (SPECIES).dx2((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[2], I, DX); \ - (SPECIES).i3((INDEX)) = I; \ - (SPECIES).dx3((INDEX)) = DX; \ - (SPECIES).ux1((INDEX)) = U1; \ - (SPECIES).ux2((INDEX)) = U2; \ - (SPECIES).ux3((INDEX)) = U3; \ - (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + from_Xi_to_i_di(X_CU[0], I1, DX1); \ + from_Xi_to_i_di(X_CU[1], I2, DX2); \ + from_Xi_to_i_di(X_CU[2], I3, DX3); \ + init_prtl_3d_i_di(SPECIES, INDEX, I1, I2, I3, DX1, DX2, DX3, U1, U2, U3, WEIGHT); \ } #else -# define init_prtl_2d(MBLOCK, SPECIES, INDEX, X1, X2, U1, U2, U3) \ +# define init_prtl_2d(MBLOCK, SPECIES, INDEX, X1, X2, U1, U2, U3, WEIGHT) \ { \ coord_t X_CU; \ vec_t U_C { ZERO, ZERO, ZERO }; \ - int I; \ - float DX; \ + int I1, I2; \ + float DX1, DX2; \ ((MBLOCK).metric).x_Sph2Code({ (X1), (X2) }, X_CU); \ ((MBLOCK).metric).v_Hat2Cart({ X_CU[0], X_CU[1], ZERO }, { U1, U2, U3 }, U_C); \ - from_Xi_to_i_di(X_CU[0], I, DX); \ - (SPECIES).i1((INDEX)) = I; \ - (SPECIES).dx1((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[1], I, DX); \ - (SPECIES).i2((INDEX)) = I; \ - (SPECIES).dx2((INDEX)) = DX; \ - (SPECIES).ux1((INDEX)) = U_C[0]; \ - (SPECIES).ux2((INDEX)) = U_C[1]; \ - (SPECIES).ux3((INDEX)) = U_C[2]; \ - (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + from_Xi_to_i_di(X_CU[0], I1, DX1); \ + from_Xi_to_i_di(X_CU[1], I2, DX2); \ + init_prtl_2d_i_di(SPECIES, INDEX, I1, I2, DX1, DX2, U_C[0], U_C[1], U_C[2], WEIGHT); \ } -# define init_prtl_3d(MBLOCK, SPECIES, INDEX, X1, X2, X3, U1, U2, U3) \ +# define init_prtl_3d(MBLOCK, SPECIES, INDEX, X1, X2, X3, U1, U2, U3, WEIGHT) \ { \ coord_t X_CU; \ vec_t U_C { ZERO, ZERO, ZERO }; \ - int I; \ - float DX; \ + int I1, I2, I3; \ + float DX1, DX2, DX3; \ ((MBLOCK).metric).x_Sph2Code({ (X1), (X2), (X3) }, X_CU); \ ((MBLOCK).metric).v_Hat2Cart(X_CU, { U1, U2, U3 }, U_C); \ - from_Xi_to_i_di(X_CU[0], I, DX); \ - (SPECIES).i1((INDEX)) = I; \ - (SPECIES).dx1((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[1], I, DX); \ - (SPECIES).i2((INDEX)) = I; \ - (SPECIES).dx2((INDEX)) = DX; \ - from_Xi_to_i_di(X_CU[2], I, DX); \ - (SPECIES).i3((INDEX)) = I; \ - (SPECIES).dx3((INDEX)) = DX; \ - (SPECIES).ux1((INDEX)) = U_C[0]; \ - (SPECIES).ux2((INDEX)) = U_C[1]; \ - (SPECIES).ux3((INDEX)) = U_C[2]; \ - (SPECIES).tag((INDEX)) = static_cast(ParticleTag::alive); \ + from_Xi_to_i_di(X_CU[0], I1, DX1); \ + from_Xi_to_i_di(X_CU[1], I2, DX2); \ + from_Xi_to_i_di(X_CU[2], I3, DX3); \ + init_prtl_3d_i_di( \ + SPECIES, INDEX, I1, I2, I3, DX1, DX2, DX3, U_C[0], U_C[1], U_C[2], WEIGHT); \ } #endif diff --git a/src/framework/sim_params.cpp b/src/framework/sim_params.cpp index 4eb734ca1..76d0c2dc7 100644 --- a/src/framework/sim_params.cpp +++ b/src/framework/sim_params.cpp @@ -55,6 +55,7 @@ namespace ntt { } m_shuffle_interval = get("particles", "shuffle_step", defaults::shuffle_interval); m_max_dead_frac = get("particles", "max_dead_frac", defaults::max_dead_frac); + m_use_weights = get("particles", "use_weights", defaults::use_weights); #ifdef MINKOWSKI_METRIC m_metric = "minkowski"; diff --git a/src/framework/sim_params.h b/src/framework/sim_params.h index 65520eccd..1ce407381 100644 --- a/src/framework/sim_params.h +++ b/src/framework/sim_params.h @@ -31,6 +31,9 @@ namespace ntt { // Vector of user-defined species parameters. std::vector m_species; + // Use particle weights + bool m_use_weights; + // Particle shuffle interval. int m_shuffle_interval; double m_max_dead_frac; @@ -130,6 +133,21 @@ namespace ntt { [[nodiscard]] auto species() const -> const std::vector& { return m_species; } + /** + * @brief Get the particle shuffle interval. + */ + [[nodiscard]] auto shuffleInterval() const -> const int& { + return m_shuffle_interval; + } + /** + * @brief Get maximum number of dead particles (as a fraction of current active particles). + */ + [[nodiscard]] auto maxDeadFraction() const -> const double& { + return m_max_dead_frac; + } + [[nodiscard]] auto useWeights() const -> const bool& { + return m_use_weights; + } /** * @brief Get the extent of the simulation box. */ @@ -208,18 +226,6 @@ namespace ntt { [[nodiscard]] auto outputMomSmooth() const -> const int& { return m_output_mom_smooth; } - /** - * @brief Get the particle shuffle interval. - */ - [[nodiscard]] auto shuffleInterval() const -> const int& { - return m_shuffle_interval; - } - /** - * @brief Get maximum number of dead particles (as a fraction of current active particles). - */ - [[nodiscard]] auto maxDeadFraction() const -> const double& { - return m_max_dead_frac; - } /** * @brief Get parameters read from the input (with default value if not found) diff --git a/src/pic/currents/deposit_currents.cpp b/src/pic/currents/deposit_currents.cpp index 7a4f165d4..cc5b3f80d 100644 --- a/src/pic/currents/deposit_currents.cpp +++ b/src/pic/currents/deposit_currents.cpp @@ -22,6 +22,7 @@ namespace ntt { template void PIC::CurrentsDeposit() { auto& mblock = this->meshblock; + auto params = *(this->params()); AssertEmptyContent(mblock.cur_content); @@ -30,7 +31,8 @@ namespace ntt { if (species.charge() != 0.0) { const real_t dt { mblock.timestep() }; const real_t charge { species.charge() }; - DepositCurrents_kernel deposit(mblock, species, scatter_cur, charge, dt); + DepositCurrents_kernel deposit( + mblock, species, scatter_cur, charge, params.useWeights(), dt); Kokkos::parallel_for("deposit", species.rangeActiveParticles(), deposit); } } diff --git a/src/pic/currents/deposit_currents.hpp b/src/pic/currents/deposit_currents.hpp index 7feefee55..9c546c42b 100644 --- a/src/pic/currents/deposit_currents.hpp +++ b/src/pic/currents/deposit_currents.hpp @@ -26,6 +26,7 @@ namespace ntt { scatter_ndfield_t m_scatter_cur; const real_t m_charge, m_dt; const real_t m_xi2max; + const bool m_use_weights; public: /** @@ -40,11 +41,13 @@ namespace ntt { const Particles& particles, const scatter_ndfield_t& scatter_cur, const real_t& charge, + const bool& use_weights, const real_t& dt) : m_mblock(mblock), m_particles(particles), m_scatter_cur(scatter_cur), m_charge(charge), + m_use_weights { use_weights }, m_dt(dt), m_xi2max((real_t)(m_mblock.i2_max()) - (real_t)(N_GHOSTS)) {} @@ -61,12 +64,20 @@ namespace ntt { // get [i, di]_init and [i, di]_final (per dimension) getDepositInterval(p, vp, Ip_f, Ip_i, xp_f, xp_i, xp_r); - depositCurrentsFromParticle(vp, Ip_f, Ip_i, xp_f, xp_i, xp_r); + depositCurrentsFromParticle(m_use_weights ? static_cast(m_particles.weight(p)) + : ONE, + vp, + Ip_f, + Ip_i, + xp_f, + xp_i, + xp_r); } } /** * @brief Deposit currents from a single particle. + * @param[in] weight Particle weight or 1 if no weights used. * @param[in] vp Particle 3-velocity. * @param[in] Ip_f Final position of the particle (cell index). * @param[in] Ip_i Initial position of the particle (cell index). @@ -74,7 +85,8 @@ namespace ntt { * @param[in] xp_i Previous step position. * @param[in] xp_r Intermediate point used in zig-zag deposit. */ - Inline void depositCurrentsFromParticle(const vec_t& vp, + Inline void depositCurrentsFromParticle(const real_t& weight, + const vec_t& vp, const tuple_t& Ip_f, const tuple_t& Ip_i, const coord_t& xp_f, @@ -171,8 +183,12 @@ namespace ntt { } }; + /** + * !TODO: fix the conversion to I+di + */ template <> Inline void DepositCurrents_kernel::depositCurrentsFromParticle( + const real_t& weight, const vec_t& vp, const tuple_t& Ip_f, const tuple_t& Ip_i, @@ -181,14 +197,14 @@ namespace ntt { const coord_t& xp_r) const { real_t Wx1_1 { HALF * (xp_i[0] + xp_r[0]) - static_cast(Ip_i[0]) }; real_t Wx1_2 { HALF * (xp_f[0] + xp_r[0]) - static_cast(Ip_f[0]) }; - real_t Fx1_1 { (xp_r[0] - xp_i[0]) * m_charge / m_dt }; - real_t Fx1_2 { (xp_f[0] - xp_r[0]) * m_charge / m_dt }; + real_t Fx1_1 { (xp_r[0] - xp_i[0]) * weight * m_charge / m_dt }; + real_t Fx1_2 { (xp_f[0] - xp_r[0]) * weight * m_charge / m_dt }; - real_t Fx2_1 { HALF * vp[1] * m_charge }; - real_t Fx2_2 { HALF * vp[1] * m_charge }; + real_t Fx2_1 { HALF * vp[1] * weight * m_charge }; + real_t Fx2_2 { HALF * vp[1] * weight * m_charge }; - real_t Fx3_1 { HALF * vp[2] * m_charge }; - real_t Fx3_2 { HALF * vp[2] * m_charge }; + real_t Fx3_1 { HALF * vp[2] * weight * m_charge }; + real_t Fx3_2 { HALF * vp[2] * weight * m_charge }; auto cur_access = m_scatter_cur.access(); ATOMIC_JX1(Ip_i[0]) += Fx1_1; @@ -205,11 +221,9 @@ namespace ntt { ATOMIC_JX3(Ip_f[0] + 1) += Fx3_2 * Wx1_2; } - /** - * !TODO: fix the conversion to I+di - */ template <> Inline void DepositCurrents_kernel::depositCurrentsFromParticle( + const real_t& weight, const vec_t& vp, const tuple_t& Ip_f, const tuple_t& Ip_i, @@ -218,16 +232,16 @@ namespace ntt { const coord_t& xp_r) const { real_t Wx1_1 { HALF * (xp_i[0] + xp_r[0]) - static_cast(Ip_i[0]) }; real_t Wx1_2 { HALF * (xp_f[0] + xp_r[0]) - static_cast(Ip_f[0]) }; - real_t Fx1_1 { (xp_r[0] - xp_i[0]) * m_charge / m_dt }; - real_t Fx1_2 { (xp_f[0] - xp_r[0]) * m_charge / m_dt }; + real_t Fx1_1 { (xp_r[0] - xp_i[0]) * weight * m_charge / m_dt }; + real_t Fx1_2 { (xp_f[0] - xp_r[0]) * weight * m_charge / m_dt }; real_t Wx2_1 { HALF * (xp_i[1] + xp_r[1]) - static_cast(Ip_i[1]) }; real_t Wx2_2 { HALF * (xp_f[1] + xp_r[1]) - static_cast(Ip_f[1]) }; - real_t Fx2_1 { (xp_r[1] - xp_i[1]) * m_charge / m_dt }; - real_t Fx2_2 { (xp_f[1] - xp_r[1]) * m_charge / m_dt }; + real_t Fx2_1 { (xp_r[1] - xp_i[1]) * weight * m_charge / m_dt }; + real_t Fx2_2 { (xp_f[1] - xp_r[1]) * weight * m_charge / m_dt }; - real_t Fx3_1 { HALF * vp[2] * m_charge }; - real_t Fx3_2 { HALF * vp[2] * m_charge }; + real_t Fx3_1 { HALF * vp[2] * weight * m_charge }; + real_t Fx3_2 { HALF * vp[2] * weight * m_charge }; auto cur_access = m_scatter_cur.access(); ATOMIC_JX1(Ip_i[0], Ip_i[1]) += Fx1_1 * (ONE - Wx2_1); @@ -253,6 +267,7 @@ namespace ntt { template <> Inline void DepositCurrents_kernel::depositCurrentsFromParticle( + const real_t& weight, const vec_t&, const tuple_t& Ip_f, const tuple_t& Ip_i, @@ -261,18 +276,18 @@ namespace ntt { const coord_t& xp_r) const { real_t Wx1_1 { HALF * (xp_i[0] + xp_r[0]) - static_cast(Ip_i[0]) }; real_t Wx1_2 { HALF * (xp_f[0] + xp_r[0]) - static_cast(Ip_f[0]) }; - real_t Fx1_1 { (xp_r[0] - xp_i[0]) * m_charge / m_dt }; - real_t Fx1_2 { (xp_f[0] - xp_r[0]) * m_charge / m_dt }; + real_t Fx1_1 { (xp_r[0] - xp_i[0]) * weight * m_charge / m_dt }; + real_t Fx1_2 { (xp_f[0] - xp_r[0]) * weight * m_charge / m_dt }; real_t Wx2_1 { HALF * (xp_i[1] + xp_r[1]) - static_cast(Ip_i[1]) }; real_t Wx2_2 { HALF * (xp_f[1] + xp_r[1]) - static_cast(Ip_f[1]) }; - real_t Fx2_1 { (xp_r[1] - xp_i[1]) * m_charge / m_dt }; - real_t Fx2_2 { (xp_f[1] - xp_r[1]) * m_charge / m_dt }; + real_t Fx2_1 { (xp_r[1] - xp_i[1]) * weight * m_charge / m_dt }; + real_t Fx2_2 { (xp_f[1] - xp_r[1]) * weight * m_charge / m_dt }; real_t Wx3_1 { HALF * (xp_i[2] + xp_r[2]) - static_cast(Ip_i[2]) }; real_t Wx3_2 { HALF * (xp_f[2] + xp_r[2]) - static_cast(Ip_f[2]) }; - real_t Fx3_1 { (xp_r[2] - xp_i[2]) * m_charge / m_dt }; - real_t Fx3_2 { (xp_f[2] - xp_r[2]) * m_charge / m_dt }; + real_t Fx3_1 { (xp_r[2] - xp_i[2]) * weight * m_charge / m_dt }; + real_t Fx3_2 { (xp_f[2] - xp_r[2]) * weight * m_charge / m_dt }; auto cur_access = m_scatter_cur.access(); ATOMIC_JX1(Ip_i[0], Ip_i[1], Ip_i[2]) += Fx1_1 * (ONE - Wx2_1) * (ONE - Wx3_1); From 51f5f81a1225fae4936982eb552368bc1f7ad94b Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 6 Mar 2023 19:25:27 -0500 Subject: [PATCH 4/5] =?UTF-8?q?=E2=9C=A8=20weights=20in=20pgen=20(explicit?= =?UTF-8?q?)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/pic/pgen/debug/deposit.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pic/pgen/debug/deposit.hpp b/src/pic/pgen/debug/deposit.hpp index b95578137..4de8fabfc 100644 --- a/src/pic/pgen/debug/deposit.hpp +++ b/src/pic/pgen/debug/deposit.hpp @@ -29,10 +29,10 @@ namespace ntt { positrons.setNpart(2); Kokkos::parallel_for( "init", CreateRangePolicy({ 0 }, { 1 }), Lambda(index_t) { - init_prtl_2d(mblock, electrons, 0, 1.2, 0.1, 0.0, -0.2, 0.0); - init_prtl_2d(mblock, positrons, 0, 1.2, 0.1, 0.0, 0.2, 0.0); - init_prtl_2d(mblock, electrons, 1, 1.2, constant::PI - 0.1, 0.0, 0.2, 0.0); - init_prtl_2d(mblock, positrons, 1, 1.2, constant::PI - 0.1, 0.0, -0.2, 0.0); + init_prtl_2d(mblock, electrons, 0, 1.2, 0.1, 0.0, -0.2, 0.0, 1.0); + init_prtl_2d(mblock, positrons, 0, 1.2, 0.1, 0.0, 0.2, 0.0, 1.0); + init_prtl_2d(mblock, electrons, 1, 1.2, constant::PI - 0.1, 0.0, 0.2, 0.0, 1.0); + init_prtl_2d(mblock, positrons, 1, 1.2, constant::PI - 0.1, 0.0, -0.2, 0.0, 1.0); }); } From a77429082e1b6bc0e5f9e44fb1e449b715399d64 Mon Sep 17 00:00:00 2001 From: hayk Date: Tue, 7 Mar 2023 01:53:44 -0500 Subject: [PATCH 5/5] =?UTF-8?q?=E2=9C=A8=20updated=20injector=20to=20use?= =?UTF-8?q?=20weights?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 2 +- src/framework/injector.hpp | 351 +++++++++++++++++---------------- src/pic/pgen/magnetosphere.hpp | 7 +- 3 files changed, 185 insertions(+), 175 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f5d8d900..9cecb4a3c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.16) set(PROJECT_NAME entity) project(${PROJECT_NAME} - VERSION 0.8.3 + VERSION 0.8.4.1 LANGUAGES CXX C ) diff --git a/src/framework/injector.hpp b/src/framework/injector.hpp index 3012466c8..d4ef2edac 100644 --- a/src/framework/injector.hpp +++ b/src/framework/injector.hpp @@ -47,9 +47,9 @@ namespace ntt { vec_t v { ZERO }; x[0] = rand_gen.frand(region[0], region[1]); energy_dist(x, v, species_index1); - init_prtl_1d(mblock, species1, p + offset1, x[0], v[0], v[1], v[2]); + init_prtl_1d(mblock, species1, p + offset1, x[0], v[0], v[1], v[2], ONE); energy_dist(x, v, species_index2); - init_prtl_1d(mblock, species2, p + offset2, x[0], v[0], v[1], v[2]); + init_prtl_1d(mblock, species2, p + offset2, x[0], v[0], v[1], v[2], ONE); pool.free_state(rand_gen); } @@ -94,9 +94,9 @@ namespace ntt { x[0] = rand_gen.frand(region[0], region[1]); x[1] = rand_gen.frand(region[2], region[3]); energy_dist(x, v, species_index1); - init_prtl_2d(mblock, species1, p + offset1, x[0], x[1], v[0], v[1], v[2]); + init_prtl_2d(mblock, species1, p + offset1, x[0], x[1], v[0], v[1], v[2], ONE); energy_dist(x, v, species_index2); - init_prtl_2d(mblock, species2, p + offset2, x[0], x[1], v[0], v[1], v[2]); + init_prtl_2d(mblock, species2, p + offset2, x[0], x[1], v[0], v[1], v[2], ONE); pool.free_state(rand_gen); } @@ -142,9 +142,9 @@ namespace ntt { x[1] = rand_gen.frand(region[2], region[3]); x[2] = rand_gen.frand(region[4], region[5]); energy_dist(x, v, species_index1); - init_prtl_3d(mblock, species1, p + offset1, x[0], x[1], x[2], v[0], v[1], v[2]); + init_prtl_3d(mblock, species1, p + offset1, x[0], x[1], x[2], v[0], v[1], v[2], ONE); energy_dist(x, v, species_index2); - init_prtl_3d(mblock, species2, p + offset2, x[0], x[1], x[2], v[0], v[1], v[2]); + init_prtl_3d(mblock, species2, p + offset2, x[0], x[1], x[2], v[0], v[1], v[2], ONE); pool.free_state(rand_gen); } @@ -248,7 +248,6 @@ namespace ntt { /* -------------------------------------------------------------------------- */ /* Volume injection kernels and routines */ /* -------------------------------------------------------------------------- */ - template class EnDist, @@ -274,70 +273,62 @@ namespace ntt { offset2 { sp2.npart() }, index { ind }, nppc { ppc }, + use_weights { params.useWeights() }, energy_dist { params, mblock }, spatial_dist { params, mblock }, inj_criterion { params, mblock }, pool { *(mblock.random_pool_ptr) } {} Inline void operator()(index_t i1) const { // cell node - coord_t xi { static_cast(static_cast(i1) - N_GHOSTS) }; - // cell center - coord_t xc { xi[0] + HALF }; - // physical coordinate - coord_t xph { ZERO }; - -#ifdef MINKOWSKI_METRIC - mblock.metric.x_Code2Cart(xc, xph); -#else - mblock.metric.x_Code2Sph(xc, xph); -#endif - - if (inj_criterion(xph)) { - typename RandomNumberPool_t::generator_type rand_gen = pool.get_state(); - real_t ninject = nppc * spatial_dist(xph); - while (ninject > ZERO) { - real_t random = rand_gen.frand(); - if (random < ninject) { - real_t dx1 = rand_gen.frand(); - real_t dx2 = rand_gen.frand(); - - auto p = Kokkos::atomic_fetch_add(&index(), 1); - - vec_t v { ZERO }, v_cart { ZERO }; - energy_dist(xph, v, 0); -#ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; -#else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); -#endif - species1.i1(offset1 + p) = static_cast(i1) - N_GHOSTS; - species1.dx1(offset1 + p) = dx1; - species1.ux1(offset1 + p) = v[0]; - species1.ux2(offset1 + p) = v[1]; - species1.ux3(offset1 + p) = v[2]; - species1.tag(offset1 + p) = static_cast(ParticleTag::alive); - - energy_dist(xph, v, 1); -#ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; -#else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); -#endif - species2.i1(offset2 + p) = static_cast(i1) - N_GHOSTS; - species2.dx1(offset2 + p) = dx1; - species2.ux1(offset2 + p) = v[0]; - species2.ux2(offset2 + p) = v[1]; - species2.ux3(offset2 + p) = v[2]; - species2.tag(offset2 + p) = static_cast(ParticleTag::alive); - } - ninject -= ONE; + coord_t xi { static_cast(static_cast(i1) - N_GHOSTS) }; + RandomGenerator_t rand_gen { pool.get_state() }; + real_t n_inject { nppc }; + coord_t xc { ZERO }; + coord_t xph { ZERO }; + float dx1, dx2; + vec_t v { ZERO }, v_cart { ZERO }; + real_t cell_vol; + + while (n_inject > ZERO) { + dx1 = Random(rand_gen); + xc[0] = xi[0] + dx1; + mblock.metric.x_Code2Cart(xc, xph); + if ((Random(rand_gen) < n_inject) && // # of prtls + inj_criterion(xph) && // injection criterion + (Random(rand_gen) < spatial_dist(xph)) // spatial distribution + ) { + auto p { Kokkos::atomic_fetch_add(&index(), 1) }; + cell_vol = mblock.metric.sqrt_det_h(xc) / mblock.metric.min_cell_volume(); + + energy_dist(xph, v, species_index1); + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; + init_prtl_1d_i_di(species1, + offset1 + p, + static_cast(i1) - N_GHOSTS, + dx1, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); + + energy_dist(xph, v, species_index2); + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; + init_prtl_1d_i_di(species2, + offset2 + p, + static_cast(i1) - N_GHOSTS, + dx1, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); } - pool.free_state(rand_gen); + n_inject -= ONE; } + pool.free_state(rand_gen); } private: @@ -348,6 +339,7 @@ namespace ntt { const std::size_t offset1, offset2; array_t index; const real_t nppc; + const bool use_weights; EnDist energy_dist; SpDist spatial_dist; InjCrit inj_criterion; @@ -379,75 +371,81 @@ namespace ntt { offset2 { sp2.npart() }, index { ind }, nppc { ppc }, + use_weights { params.useWeights() }, energy_dist { params, mblock }, spatial_dist { params, mblock }, inj_criterion { params, mblock }, pool { *(mblock.random_pool_ptr) } {} Inline void operator()(index_t i1, index_t i2) const { // cell node - coord_t xi { static_cast(static_cast(i1) - N_GHOSTS), + coord_t xi { static_cast(static_cast(i1) - N_GHOSTS), static_cast(static_cast(i2) - N_GHOSTS) }; - // cell center - coord_t xc { xi[0] + HALF, xi[1] + HALF }; - // physical coordinate - coord_t xph { ZERO }; - + RandomGenerator_t rand_gen { pool.get_state() }; + real_t n_inject { nppc }; + coord_t xc { ZERO }; + coord_t xph { ZERO }; + float dx1, dx2; + vec_t v { ZERO }, v_cart { ZERO }; + real_t cell_vol; + + while (n_inject > ZERO) { + dx1 = Random(rand_gen); + dx2 = Random(rand_gen); + xc[0] = xi[0] + dx1; + xc[1] = xi[1] + dx2; #ifdef MINKOWSKI_METRIC - mblock.metric.x_Code2Cart(xc, xph); + mblock.metric.x_Code2Cart(xc, xph); #else - mblock.metric.x_Code2Sph(xc, xph); + mblock.metric.x_Code2Sph(xc, xph); #endif - - if (inj_criterion(xph)) { - typename RandomNumberPool_t::generator_type rand_gen = pool.get_state(); - real_t ninject = nppc * spatial_dist(xph); - while (ninject > ZERO) { - real_t random = rand_gen.frand(); - if (random < ninject) { - real_t dx1 = rand_gen.frand(); - real_t dx2 = rand_gen.frand(); - - auto p = Kokkos::atomic_fetch_add(&index(), 1); - - vec_t v { ZERO }, v_cart { ZERO }; - energy_dist(xph, v, 0); + if ((Random(rand_gen) < n_inject) && // # of prtls + inj_criterion(xph) && // injection criterion + (Random(rand_gen) < spatial_dist(xph)) // spatial distribution + ) { + auto p { Kokkos::atomic_fetch_add(&index(), 1) }; + cell_vol = mblock.metric.sqrt_det_h(xc) / mblock.metric.min_cell_volume(); + + energy_dist(xph, v, species_index1); #ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; #else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); + mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); #endif - species1.i1(offset1 + p) = static_cast(i1) - N_GHOSTS; - species1.dx1(offset1 + p) = dx1; - species1.i2(offset1 + p) = static_cast(i2) - N_GHOSTS; - species1.dx2(offset1 + p) = dx2; - species1.ux1(offset1 + p) = v_cart[0]; - species1.ux2(offset1 + p) = v_cart[1]; - species1.ux3(offset1 + p) = v_cart[2]; - species1.tag(offset1 + p) = static_cast(ParticleTag::alive); - - energy_dist(xph, v, 1); + init_prtl_2d_i_di(species1, + offset1 + p, + static_cast(i1) - N_GHOSTS, + static_cast(i2) - N_GHOSTS, + dx1, + dx2, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); + + energy_dist(xph, v, species_index2); #ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; #else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); + mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); #endif - species2.i1(offset2 + p) = static_cast(i1) - N_GHOSTS; - species2.dx1(offset2 + p) = dx1; - species2.i2(offset2 + p) = static_cast(i2) - N_GHOSTS; - species2.dx2(offset2 + p) = dx2; - species2.ux1(offset2 + p) = v_cart[0]; - species2.ux2(offset2 + p) = v_cart[1]; - species2.ux3(offset2 + p) = v_cart[2]; - species2.tag(offset2 + p) = static_cast(ParticleTag::alive); - } - ninject -= ONE; + init_prtl_2d_i_di(species2, + offset2 + p, + static_cast(i1) - N_GHOSTS, + static_cast(i2) - N_GHOSTS, + dx1, + dx2, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); } - pool.free_state(rand_gen); + n_inject -= ONE; } + pool.free_state(rand_gen); } private: @@ -458,6 +456,7 @@ namespace ntt { const std::size_t offset1, offset2; array_t index; const real_t nppc; + const bool use_weights; EnDist energy_dist; SpDist spatial_dist; InjCrit inj_criterion; @@ -483,95 +482,105 @@ namespace ntt { mblock { mb }, species1 { sp1 }, species2 { sp2 }, + species_index1 { sp1.index() }, + species_index2 { sp2.index() }, offset1 { sp1.npart() }, offset2 { sp2.npart() }, index { ind }, nppc { ppc }, + use_weights { params.useWeights() }, energy_dist { params, mblock }, spatial_dist { params, mblock }, inj_criterion { params, mblock }, pool { *(mblock.random_pool_ptr) } {} Inline void operator()(index_t i1, index_t i2, index_t i3) const { // cell node - coord_t xi { static_cast(static_cast(i1) - N_GHOSTS), + coord_t xi { static_cast(static_cast(i1) - N_GHOSTS), static_cast(static_cast(i2) - N_GHOSTS), static_cast(static_cast(i3) - N_GHOSTS) }; - // cell center - coord_t xc { xi[0] + HALF, xi[1] + HALF, xi[2] + HALF }; - // physical coordinate - coord_t xph { ZERO }; - + RandomGenerator_t rand_gen { pool.get_state() }; + real_t n_inject { nppc }; + coord_t xc { ZERO }; + coord_t xph { ZERO }; + float dx1, dx2, dx3; + vec_t v { ZERO }, v_cart { ZERO }; + real_t cell_vol; + + while (n_inject > ZERO) { + dx1 = Random(rand_gen); + dx2 = Random(rand_gen); + dx3 = Random(rand_gen); + xc[0] = xi[0] + dx1; + xc[1] = xi[1] + dx2; + xc[2] = xi[2] + dx3; #ifdef MINKOWSKI_METRIC - mblock.metric.x_Code2Cart(xc, xph); + mblock.metric.x_Code2Cart(xc, xph); #else - mblock.metric.x_Code2Sph(xc, xph); + mblock.metric.x_Code2Sph(xc, xph); #endif - - if (inj_criterion(xph)) { - typename RandomNumberPool_t::generator_type rand_gen = pool.get_state(); - - real_t ninject = nppc * spatial_dist(xph); - while (ninject > ZERO) { - real_t random = rand_gen.frand(); - if (random < ninject) { - real_t dx1 = rand_gen.frand(); - real_t dx2 = rand_gen.frand(); - real_t dx3 = rand_gen.frand(); - - auto p = Kokkos::atomic_fetch_add(&index(), 1); - - vec_t v { ZERO }, v_cart { ZERO }; - energy_dist(xph, v, 0); + if ((Random(rand_gen) < n_inject) && // # of prtls + inj_criterion(xph) && // injection criterion + (Random(rand_gen) < spatial_dist(xph)) // spatial distribution + ) { + auto p { Kokkos::atomic_fetch_add(&index(), 1) }; + cell_vol = mblock.metric.sqrt_det_h(xc) / mblock.metric.min_cell_volume(); + + energy_dist(xph, v, species_index1); #ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; #else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); + mblock.metric.v_Hat2Cart({ xc[0], xc[1], xc[2] }, v, v_cart); #endif - species1.i1(offset1 + p) = static_cast(i1) - N_GHOSTS; - species1.dx1(offset1 + p) = dx1; - species1.i2(offset1 + p) = static_cast(i2) - N_GHOSTS; - species1.dx2(offset1 + p) = dx2; - species1.i3(offset1 + p) = static_cast(i3) - N_GHOSTS; - species1.dx3(offset1 + p) = dx3; - species1.ux1(offset1 + p) = v_cart[0]; - species1.ux2(offset1 + p) = v_cart[1]; - species1.ux3(offset1 + p) = v_cart[2]; - species1.tag(offset1 + p) = static_cast(ParticleTag::alive); - - energy_dist(xph, v, 1); + init_prtl_3d_i_di(species1, + offset1 + p, + static_cast(i1) - N_GHOSTS, + static_cast(i2) - N_GHOSTS, + static_cast(i3) - N_GHOSTS, + dx1, + dx2, + dx3, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); + + energy_dist(xph, v, species_index2); #ifdef MINKOWSKI_METRIC - v_cart[0] = v[0]; - v_cart[1] = v[1]; - v_cart[2] = v[2]; + v_cart[0] = v[0]; + v_cart[1] = v[1]; + v_cart[2] = v[2]; #else - mblock.metric.v_Hat2Cart({ xc[0], xc[1], ZERO }, v, v_cart); + mblock.metric.v_Hat2Cart({ xc[0], xc[1], xc[2] }, v, v_cart); #endif - species2.i1(offset2 + p) = static_cast(i1) - N_GHOSTS; - species2.dx1(offset2 + p) = dx1; - species2.i2(offset2 + p) = static_cast(i2) - N_GHOSTS; - species2.dx2(offset2 + p) = dx2; - species2.i3(offset2 + p) = static_cast(i3) - N_GHOSTS; - species2.dx3(offset2 + p) = dx3; - species2.ux1(offset2 + p) = v_cart[0]; - species2.ux2(offset2 + p) = v_cart[1]; - species2.ux3(offset2 + p) = v_cart[2]; - species2.tag(offset2 + p) = static_cast(ParticleTag::alive); - } - ninject -= ONE; + init_prtl_3d_i_di(species2, + offset2 + p, + static_cast(i1) - N_GHOSTS, + static_cast(i2) - N_GHOSTS, + static_cast(i3) - N_GHOSTS, + dx1, + dx2, + dx3, + v_cart[0], + v_cart[1], + v_cart[2], + use_weights ? cell_vol : ONE); } - pool.free_state(rand_gen); + n_inject -= ONE; } + pool.free_state(rand_gen); } private: SimulationParams params; Meshblock mblock; Particles species1, species2; + const int species_index1, species_index2; const std::size_t offset1, offset2; array_t index; const real_t nppc; + const bool use_weights; EnDist energy_dist; SpDist spatial_dist; InjCrit inj_criterion; diff --git a/src/pic/pgen/magnetosphere.hpp b/src/pic/pgen/magnetosphere.hpp index d9dd8c620..16b9bfb0b 100644 --- a/src/pic/pgen/magnetosphere.hpp +++ b/src/pic/pgen/magnetosphere.hpp @@ -153,8 +153,8 @@ namespace ntt { template struct InjectionShell : public SpatialDistribution { explicit InjectionShell(const SimulationParams& params, Meshblock& mblock) - : SpatialDistribution(params, mblock) { - inj_rmax = params.get("problem", "inj_rmax", 1.5 * mblock.metric.x1_min); + : SpatialDistribution(params, mblock), + inj_rmax { params.get("problem", "inj_rmax", 1.5 * mblock.metric.x1_min) } { const int buff_cells = 10; coord_t xcu { ZERO }, xph { ZERO }; xcu[0] = (real_t)buff_cells; @@ -166,7 +166,8 @@ namespace ntt { } private: - real_t inj_rmax, inj_rmin; + const real_t inj_rmax; + real_t inj_rmin; }; template