diff --git a/extern/Kokkos b/extern/Kokkos index 5fc08a9a7..6f69d94a4 160000 --- a/extern/Kokkos +++ b/extern/Kokkos @@ -1 +1 @@ -Subproject commit 5fc08a9a7da14d8530f8c7035d008ef63ddb4e5c +Subproject commit 6f69d94a4abb153a22848d897a816f1a2fe4e3e3 diff --git a/extern/adios2 b/extern/adios2 index a6e8314cc..f80ad829d 160000 --- a/extern/adios2 +++ b/extern/adios2 @@ -1 +1 @@ -Subproject commit a6e8314cc3c0b28d496b44dcd4f15685013b887b +Subproject commit f80ad829d751241140c40923503e1888e27e22e1 diff --git a/extern/plog b/extern/plog index 85a871b13..96637a6e5 160000 --- a/extern/plog +++ b/extern/plog @@ -1 +1 @@ -Subproject commit 85a871b13be0bd1a9e0110744fa60cc9bd1e8380 +Subproject commit 96637a6e5e53f54e4e56d667d312c564d979ec0e diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index d082d617f..f90d90f56 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -4,7 +4,7 @@ * @implements * - ntt::GRPICEngine<> : ntt::Engine<> * @cpp: - * - srpic.cpp + * - grpic.cpp * @namespaces: * - ntt:: */ @@ -43,4 +43,8 @@ namespace ntt { } // namespace ntt + + + + #endif // ENGINES_GRPIC_GRPIC_H diff --git a/src/engines/grpic_1d.hpp b/src/engines/grpic_1d.hpp new file mode 100644 index 000000000..04706aa12 --- /dev/null +++ b/src/engines/grpic_1d.hpp @@ -0,0 +1,245 @@ +/** + * @file engines/grpic_1d.hpp + * @brief Simulation engien class which specialized on 1D GRPIC + * @implements + * - ntt::GRPICEngine_1D<> : ntt::Engine<> + * @cpp: + * - grpic_1d.cpp + * @namespaces: + * - ntt:: + */ + +#ifndef ENGINES_GRPIC_GRPIC_1D_H +#define ENGINES_GRPIC_GRPIC_1D_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "arch/traits.h" +#include "utils/log.h" +#include "utils/numeric.h" +#include "utils/timer.h" + +#include "archetypes/particle_injector.h" +#include "framework/domain/domain.h" +#include "framework/parameters.h" + +#include "metrics/boyer_lindq_tp.h" + +#include "engines/engine.hpp" +#include "kernels/ampere_gr.hpp" +#include "kernels/currents_deposit.hpp" +#include "kernels/digital_filter.hpp" +#include "kernels/fields_bcs.hpp" +#include "kernels/particle_moments.hpp" +#include "kernels/particle_pusher_1D_gr.hpp" +#include "pgen.hpp" + +#include +#include + +#include +#include + +namespace ntt { + + template + class GRPICEngine_1D : public Engine { + + using base_t = Engine; + using pgen_t = user::PGen; + using domain_t = Domain; + // constexprs + using base_t::pgen_is_ok; + // contents + using base_t::m_metadomain; + using base_t::m_params; + using base_t::m_pgen; + // methods + using base_t::init; + // variables + using base_t::dt; + using base_t::max_steps; + using base_t::runtime; + using base_t::step; + using base_t::time; + + auto constexpr D { M::Dim }; + + public: + static constexpr auto S { SimEngine::GRPIC }; + + GRPICEngine_1D(SimulationParams& params) + : Engine { params } { + raise::ErrorIf(D != Dim::_1D, "GRPICEngine_1D only works in 1D", HERE); + raise::ErrorIf(M.Lable != "boyer_lindq_tp", "GRPICEngine_1D only works with BoyerLindqTP metric", HERE); + } + + ~GRPICEngine_1D() = default; + + void step_forward(timer::Timers&, Domain&) override { + const auto sort_interval = m_params.template get( + "particles.sort_interval"); + + if (step == 0) { + m_metadomain.CommunicateFields(dom, Comm::D); + /** + * !CommunicateFields, ParticleInjector: Special version for 1D GRPIC needed + */ + ParticleInjector(dom); + } + + { + timers.start("ParticlePusher"); + ParticlePush(dom); + timers.stop("ParticlePusher"); + + timers.start("CurrentDeposit"); + Kokkos::deep_copy(dom.fields.cur, ZERO); + CurrentsDeposit(dom); + timers.stop("CurrentDeposit"); + + timers.start("Communications"); + m_metadomain.SynchronizeFields(dom, Comm::J); + /** + * !SynchronizeFields: Special version for 1D GRPIC needed + */ + m_metadomain.CommunicateFields(dom, Comm::J); + timers.stop("Communications"); + + timers.start("CurrentFiltering"); + CurrentsFilter(dom); + /** + * !CurrentsFilter: Special version for 1D GRPIC needed + */ + timers.stop("CurrentFiltering"); + + timers.start("Communications"); + if ((sort_interval > 0) and (step % sort_interval == 0)) { + m_metadomain.CommunicateParticles(dom, &timers); + } + timers.stop("Communications"); + } + + { + timers.start("FieldSolver"); + CurrentsAmpere(dom); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::D | Comm::J); + timers.stop("Communications"); + + } + + { + timers.start("Injector"); + DischargeInjector(dom); + timers.stop("Injector"); + } + } + + void ParticlePush(domain_t& domain) { + for (auto& species : domain.species) { + species.set_unsorted(); + logger::Checkpoint( + fmt::format("Launching particle pusher kernel for %d [%s] : %lu", + species.index(), + species.label().c_str(), + species.npart()), + HERE); + if (species.npart() == 0) { + continue; + } + const auto q_ovr_m = species.mass() > ZERO + ? species.charge() / species.mass() + : ZERO; + // coeff = q / m dt omegaB0 + const auto coeff = q_ovr_m * dt * + m_params.template get("scales.omegaB0"); + PrtlPusher::type pusher; + if (species.pusher() == PrtlPusher::FORCEFREE) { + pusher = PrtlPusher::FORCEFREE; + } else { + raise::Fatal("Invalid particle pusher", HERE); + } + // clang-format off + Kokkos::parallel_for( + "ParticlePusher", + species.rangeActiveParticles(), + kernel::gr::Pusher_kernel>( + domain.fields.em, + species.i1, + species.i1_prev, + species.dx1, + species.dx1_prev, + species.px1, + species.tag, + domain.mesh.metric, + coeff, dt, + domain.mesh.n_active(in::x1), + 1e-2, + 30, + domain.mesh.n_active(in::x1))); + } + } + + void CurrentsDeposit(domain_t& domain) { + auto scatter_cur = Kokkos::Experimental::create_scatter_view( + domain.fields.cur); + for (auto& species : domain.species) { + logger::Checkpoint( + fmt::format("Launching currents deposit kernel for %d [%s] : %lu %f", + species.index(), + species.label().c_str(), + species.npart(), + (double)species.charge()), + HERE); + if (species.npart() == 0 || cmp::AlmostZero(species.charge())) { + continue; + } + // clang-format off + /** + * !DepositCurrents_kernel: Special version for 1D GRPIC needed + */ + Kokkos::parallel_for("CurrentsDeposit", + species.rangeActiveParticles(), + kernel::DepositCurrents_kernel( + scatter_cur, + species.i1, + species.i1_prev, + species.dx1, + species.dx1_prev, + species.px1, + species.weight, + species.tag, + domain.mesh.metric, + (real_t)(species.charge()), + dt)); + } + Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); + } + + void CurrentsAmpere(domain_t& domain) { + logger::Checkpoint("Launching Ampere kernel for adding currents", HERE); + const auto q0 = m_params.template get("scales.q0"); + const auto n0 = m_params.template get("scales.n0"); + const auto B0 = m_params.template get("scales.B0"); + const auto coeff = -dt * q0 * n0 / B0; + // clang-format off + Kokkos::parallel_for( + "Ampere", + domain.mesh.rangeActiveCells(), + kernel::gr::CurrentsAmpere_kernel_1D(domain.fields.em, + domain.fields.cur, + domain.mesh.metric, + coeff)); + } + + +}; // class GRPICEngine_1D + +} // namespace ntt + +#endif // ENGINES_GRPIC_GRPIC_1D_H diff --git a/src/global/enums.h b/src/global/enums.h index 57822dec4..ff8a92a3f 100644 --- a/src/global/enums.h +++ b/src/global/enums.h @@ -2,9 +2,9 @@ * @file enums.h * @brief Special enum variables describing the simulation * @implements - * - enum ntt::Coord // Cart, Sph, Qsph + * - enum ntt::Coord // Cart, Sph, Qsph, Fs * - enum ntt::Metric // Minkowski, Spherical, QSpherical, - * Kerr_Schild, QKerr_Schild, Kerr_Schild_0 + * Kerr_Schild, QKerr_Schild, Kerr_Schild_0, BoyerLindqTP * - enum ntt::SimEngine // SRPIC, GRPIC * - enum ntt::PrtlBC // periodic, absorb, atmosphere, custom, * reflect, horizon, axis, sync @@ -134,12 +134,13 @@ namespace ntt { Cart = 1, Sph = 2, Qsph = 3, + Bltp = 4, }; constexpr Coord(uint8_t c) : enums_hidden::BaseEnum { c } {} - static constexpr type variants[] = { Cart, Sph, Qsph }; - static constexpr const char* lookup[] = { "cart", "sph", "qsph" }; + static constexpr type variants[] = { Cart, Sph, Qsph, Bltp }; + static constexpr const char* lookup[] = { "cart", "sph", "qsph", "bltp" }; static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]); }; @@ -154,16 +155,19 @@ namespace ntt { Kerr_Schild = 4, QKerr_Schild = 5, Kerr_Schild_0 = 6, + BoyerLindqTP = 7, }; constexpr Metric(uint8_t c) : enums_hidden::BaseEnum { c } {} static constexpr type variants[] = { Minkowski, Spherical, QSpherical, Kerr_Schild, - QKerr_Schild, Kerr_Schild_0 }; - static constexpr const char* lookup[] = { "minkowski", "spherical", - "qspherical", "kerr_schild", - "qkerr_schild", "kerr_schild_0" }; + QKerr_Schild, Kerr_Schild_0, + BoyerLindqTP }; + static constexpr const char* lookup[] = { "minkowski", "spherical", + "qspherical", "kerr_schild", + "qkerr_schild", "kerr_schild_0", + "boyer_lindq_tp" }; static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]); }; @@ -239,18 +243,23 @@ namespace ntt { static constexpr const char* label = "prtl_pusher"; enum type : uint8_t { - INVALID = 0, - BORIS = 1, - VAY = 2, - PHOTON = 3, - NONE = 4, + INVALID = 0, + BORIS = 1, + VAY = 2, + PHOTON = 3, + FORCEFREE = 4, + NONE = 5, }; constexpr PrtlPusher(uint8_t c) : enums_hidden::BaseEnum { c } {} - static constexpr type variants[] = { BORIS, VAY, PHOTON, NONE }; - static constexpr const char* lookup[] = { "boris", "vay", "photon", "none" }; + static constexpr type variants[] = { BORIS, VAY, PHOTON, FORCEFREE, NONE }; + static constexpr const char* lookup[] = { "boris", + "vay", + "photon", + "forcefree", + "none" }; static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]); }; diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp index 1fc57398f..941c15793 100644 --- a/src/global/tests/enums.cpp +++ b/src/global/tests/enums.cpp @@ -55,16 +55,17 @@ auto main() -> int { using enum_str_t = const std::vector; - enum_str_t all_coords = { "cart", "sph", "qsph" }; - enum_str_t all_metrics = { "minkowski", "spherical", "qspherical", - "kerr_schild", "qkerr_schild", "kerr_schild_0" }; + enum_str_t all_coords = { "cart", "sph", "qsph", "bltp" }; + enum_str_t all_metrics = { "minkowski", "spherical", "qspherical", + "kerr_schild", "qkerr_schild", "kerr_schild_0", + "boyer_lindq_tp" }; enum_str_t all_simulation_engines = { "srpic", "grpic" }; enum_str_t all_particle_bcs = { "periodic", "absorb", "atmosphere", "custom", "reflect", "horizon", "axis", "sync" }; enum_str_t all_fields_bcs = { "periodic", "absorb", "atmosphere", "custom", "horizon", "conductor", "axis", "sync" }; - enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" }; - enum_str_t all_coolings = { "synchrotron", "none" }; + enum_str_t all_particle_pushers = { "boris", "vay", "photon", "forcefree", "none" }; + enum_str_t all_coolings = { "synchrotron", "none" }; enum_str_t all_out_flds = { "e", "dive", "d", "divd", "b", "h", "j", "a", "t", "rho", diff --git a/src/kernels/ampere_gr.hpp b/src/kernels/ampere_gr.hpp index 5af0fa4ef..0dd7bf906 100644 --- a/src/kernels/ampere_gr.hpp +++ b/src/kernels/ampere_gr.hpp @@ -185,3 +185,235 @@ namespace kernel::gr { } // namespace kernel::gr #endif // KERNELS_AMPERE_GR_HPP + +/** + * @file kernels/ampere_gr.hpp + * @brief Algorithms for Ampere's law in GR + * @implements + * - kernel::gr::Ampere_kernel<> + * - kernel::gr::CurrentsAmpere_kernel<> + * @namespaces: + * - kernel::gr:: + * !TODO: + * - 3D implementation + */ + +#ifndef KERNELS_AMPERE_GR_HPP +#define KERNELS_AMPERE_GR_HPP + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/error.h" +#include "utils/numeric.h" + +namespace kernel::gr { + using namespace ntt; + + /** + * @brief Algorithms for Ampere's law: + * @brief `d(Din)^i / dt = curl H_j`, `Dout += dt * d(Din)/dt`. + * @tparam M Metric. + */ + template + class Ampere_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static constexpr auto D = M::Dim; + + const ndfield_t Din; + ndfield_t Dout; + const ndfield_t H; + const M metric; + const std::size_t i2max; + const real_t coeff; + bool is_axis_i2min { false }, is_axis_i2max { false }; + + public: + Ampere_kernel(const ndfield_t& Din, + const ndfield_t& Dout, + const ndfield_t& H, + const M& metric, + real_t coeff, + std::size_t ni2, + const boundaries_t& boundaries) + : Din { Din } + , Dout { Dout } + , H { H } + , metric { metric } + , i2max { ni2 + N_GHOSTS } + , coeff { coeff } { + if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { + raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE); + raise::ErrorIf(boundaries[1].size() < 2, + "boundaries defined incorrectly", + HERE); + is_axis_i2min = (boundaries[1].first == FldsBC::AXIS); + is_axis_i2max = (boundaries[1].second == FldsBC::AXIS); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + constexpr std::size_t i2min { N_GHOSTS }; + const real_t i1_ { COORD(i1) }; + const real_t i2_ { COORD(i2) }; + + const real_t inv_sqrt_detH_0pH { ONE / + metric.sqrt_det_h({ i1_, i2_ + HALF }) }; + + if ((i2 == i2min) && is_axis_i2min) { + // theta = 0 + const real_t inv_polar_area_pH { ONE / metric.polar_area(i1_ + HALF) }; + Dout(i1, i2, em::dx1) = Din(i1, i2, em::dx1) + + inv_polar_area_pH * coeff * H(i1, i2, em::hx3); + Dout(i1, i2, em::dx2) = Din(i1, i2, em::dx2) + + coeff * inv_sqrt_detH_0pH * + (H(i1 - 1, i2, em::hx3) - H(i1, i2, em::hx3)); + } else if ((i2 == i2max) && is_axis_i2max) { + // theta = pi + const real_t inv_polar_area_pH { ONE / metric.polar_area(i1_ + HALF) }; + Dout(i1, i2, em::dx1) = Din(i1, i2, em::dx1) - inv_polar_area_pH * coeff * + H(i1, i2 - 1, em::hx3); + } else { + const real_t inv_sqrt_detH_00 { ONE / metric.sqrt_det_h({ i1_, i2_ }) }; + const real_t inv_sqrt_detH_pH0 { ONE / metric.sqrt_det_h( + { i1_ + HALF, i2_ }) }; + Dout(i1, i2, em::dx1) = Din(i1, i2, em::dx1) + + coeff * inv_sqrt_detH_pH0 * + (H(i1, i2, em::hx3) - H(i1, i2 - 1, em::hx3)); + Dout(i1, i2, em::dx2) = Din(i1, i2, em::dx2) + + coeff * inv_sqrt_detH_0pH * + (H(i1 - 1, i2, em::hx3) - H(i1, i2, em::hx3)); + Dout(i1, i2, em::dx3) = Din(i1, i2, em::dx3) + + coeff * inv_sqrt_detH_00 * + ((H(i1, i2 - 1, em::hx1) - H(i1, i2, em::hx1)) + + (H(i1, i2, em::hx2) - H(i1 - 1, i2, em::hx2))); + } + } else { + raise::KernelError(HERE, "Ampere_kernel: 2D implementation called for D != 2"); + } + } + }; + + /** + * @brief Add the currents to the D field with the appropriate conversion. + */ + template + class CurrentsAmpere_kernel { + static constexpr auto D = M::Dim; + + ndfield_t Df; + const ndfield_t J; + const M metric; + const std::size_t i2max; + const real_t coeff; + bool is_axis_i2min { false }; + bool is_axis_i2max { false }; + + public: + /** + * @brief Constructor. + * @param mblock Meshblock. + */ + CurrentsAmpere_kernel(const ndfield_t& Df, + const ndfield_t& J, + const M& metric, + real_t coeff, + std::size_t ni2, + const boundaries_t& boundaries) + : Df { Df } + , J { J } + , metric { metric } + , i2max { ni2 + N_GHOSTS } + , coeff { coeff } { + if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { + raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE); + is_axis_i2min = (boundaries[1].first == FldsBC::AXIS); + is_axis_i2max = (boundaries[1].second == FldsBC::AXIS); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + constexpr std::size_t i2min { N_GHOSTS }; + const real_t i1_ { COORD(i1) }; + const real_t i2_ { COORD(i2) }; + + const real_t inv_sqrt_detH_0pH { ONE / + metric.sqrt_det_h({ i1_, i2_ + HALF }) }; + + if ((i2 == i2min) && is_axis_i2min) { + // theta = 0 (first active cell) + Df(i1, i2, em::dx1) += J(i1, i2, cur::jx1) * HALF * coeff / + metric.polar_area(i1_ + HALF); + Df(i1, i2, em::dx2) += J(i1, i2, cur::jx2) * coeff * inv_sqrt_detH_0pH; + } else if ((i2 == i2max) && is_axis_i2max) { + // theta = pi (first ghost cell from end) + Df(i1, i2, em::dx1) += J(i1, i2, cur::jx1) * HALF * coeff / + metric.polar_area(i1_ + HALF); + } else { + // 0 < theta < pi + const real_t inv_sqrt_detH_00 { ONE / metric.sqrt_det_h({ i1_, i2_ }) }; + const real_t inv_sqrt_detH_pH0 { ONE / metric.sqrt_det_h( + { i1_ + HALF, i2_ }) }; + + Df(i1, i2, em::dx1) += J(i1, i2, cur::jx1) * coeff * inv_sqrt_detH_pH0; + Df(i1, i2, em::dx2) += J(i1, i2, cur::jx2) * coeff * inv_sqrt_detH_0pH; + Df(i1, i2, em::dx3) += J(i1, i2, cur::jx3) * coeff * inv_sqrt_detH_00; + } + } else { + raise::KernelError( + HERE, + "CurrentsAmpere_kernel: 2D implementation called for D != 2"); + } + } + }; + + /** + * @brief Ampere's law for only 1D-GRPIC. + */ + template + class Ampere_kernel_1D { + static constexpr auto D = M::Dim; + + ndfield_t Df; + const ndfield_t J; + const M metric; + const real_t coeff; + + + public: + /** + * @brief Constructor. + * @param mblock Meshblock. + */ + CurrentsAmpere_kernel_1D(ndfield_t & Df, + const ndfield_t& J, + const M& metric, + const real_t coeff) + : Df { Df } + , J { J } + , metric { metric } + , coeff { coeff } { + } + + Inline void operator()(index_t i1) const { + if constexpr (D == Dim::_1D) { + const real_t i1_ { COORD(i1) }; + + const real_t inv_sqrt_detH { ONE / metric.sqrt_det_h(i1_ + HALF) }; + + Df(i1, em::dx1) += (J(i1, cur::jx1) * inv_sqrt_detH - metric.J_ff()) * coeff ; + } else { + raise::KernelError( + HERE, + "CurrentsAmpere_kernel: 1D implementation called for D != 1"); + } + } + }; + + +} // namespace kernel::gr + +#endif // KERNELS_AMPERE_GR_HPP diff --git a/src/kernels/particle_pusher_1D_gr.hpp b/src/kernels/particle_pusher_1D_gr.hpp new file mode 100644 index 000000000..ceb1c9fb3 --- /dev/null +++ b/src/kernels/particle_pusher_1D_gr.hpp @@ -0,0 +1,307 @@ +/** + * @file kernels/particle_pusher_1D_gr.h + * @brief Implementation of the particle pusher for 1D GR + * @implements + * - kernel::gr::Pusher_kernel_1D<> + * @namespaces: + * - kernel::gr:: + * @macros: + * - MPI_ENABLED +*/ + +#ifndef KERNELS_PARTICLE_PUSHER_1D_GR_HPP +#define KERNELS_PARTICLE_PUSHER_1D_GR_HPP + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/error.h" +#include "utils/numeric.h" + +#if defined(MPI_ENABLED) + #include "arch/mpi_tags.h" +#endif + +/* -------------------------------------------------------------------------- */ +/* Local macros */ +/* -------------------------------------------------------------------------- */ +#define from_Xi_to_i(XI, I) \ + { I = static_cast((XI)); } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + +#define DERIVATIVE(func, x) \ + ((func({ x + epsilon }) - func({ x - epsilon })) / \ + (TWO * epsilon)) + + +/* -------------------------------------------------------------------------- */ + +namespace kernel::gr { + using namespace ntt; + + struct Massive_t {}; + + struct Massless_t {}; + + /** + * @brief Algorithm for the Particle pusher + * @tparam M Metric + */ + template + class Pusher_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static constexpr auto D = M::Dim; + static_assert(D == Dim::_1D, "Only 1d available for this pusher"); + + const ndfield_t Dfield; + array_t i1; + array_t i1_prev; + array_t dx1; + array_t dx1_prev; + array_t px1; + array_t tag; + const M metric; + + const real_t coeff, dt; + const int ni1; + const real_t epsilon; + const int niter; + const int i1_absorb; + + + bool is_absorb_i1min { false }, is_absorb_i1max { false }; + + public: + Pusher_kernel(const ndfield_t& Dfield, + const array_t& i1, + const array_t& i1_prev, + const array_t& dx1, + const array_t& dx1_prev, + const array_t& px1, + const array_t& tag, + const M& metric, + const real_t& coeff, + const real_t& dt, + const int& ni1, + const real_t& epsilon, + const int& niter, + const boundaries_t& boundaries) + : Dfield { Dfield } + , i1 { i1 } + , i1_prev { i1_prev } + , dx1 { dx1 } + , dx1_prev { dx1_prev } + , px1 { px1 } + , tag { tag } + , metric { metric } + , coeff { coeff } + , dt { dt } + , ni1 { ni1 } + , epsilon { epsilon } + , niter { niter } + , i1_absorb { 2 } { + + raise::ErrorIf(boundaries.size() < 1, "boundaries defined incorrectly", HERE); + is_absorb_i1min = (boundaries[0].first == PrtlBC::ABSORB) || + (boundaries[0].first == PrtlBC::HORIZON); + is_absorb_i1max = (boundaries[0].second == PrtlBC::ABSORB) || + (boundaries[0].second == PrtlBC::HORIZON); + + } + + /** + * @brief Main pusher subroutine for photon particles. + */ + //Inline void operator()(Massless_t, index_t) const; + + /** + * @brief Main pusher subroutine for massive particles. + */ + Inline void operator()(index_t) const; + + /** + * @brief Iterative geodesic pusher substep for momentum only. + * @param xp particle coordinate. + * @param vp particle velocity. + * @param vp_upd updated particle velocity [return]. + */ + Inline void ForceFreePush( const coord_t& xp, + const real_t& pp, + const real_t& ex, + real_t& pp_upd) const; + + /** + * @param xp particle coordinate. + * @param vp particle velocity. + * @param xp_upd updated particle coordinate [return]. + */ + Inline void CoordinatePush(const coord_t& xp, + const real_t& vp, + coord_t& xp_upd) const; + + // Helper functions + + /** + * @brief First order Yee mesh field interpolation to particle position. + * @param p index of the particle. + * @param e interpolated e-field + */ + Inline void interpolateFields(index_t& p, + real_t& e) const { + const int i { i1(p) + static_cast(N_GHOSTS) }; + const auto dx1_ { static_cast(dx1(p)) }; + real_t c1 = HALF * (Dfield(i, em::dx1) + Dfield(i - 1, em::dx1)); + real_t c2 = HALF * (Dfield(i, em::dx1) + Dfield(i + 1, em::dx1)); + e = c1 * (ONE - dx1_) + c2 * dx1_; + } + + /** + * @brief Compute controvariant component u^0 for massive particles. + */ + + Inline auto compute_u0(const real_t& pxi, + const coord_t& xi) const -> real_t{ + return math::sqrt((SQR(pxi) + metric.f2(xi)) / + (metric.f2(xi) * (SQR(metric.alpha(xi)) - metric.f0(xi)) + SQR(metric.f1(xi)))); + } + + + // Extra + Inline void boundaryConditions(index_t& p) const{ + if (i1(p) < i1_absorb && is_absorb_i1min) { + tag(p) = ParticleTag::dead; + } else if (i1(p) >= ni1 && is_absorb_i1max) { + tag(p) = ParticleTag::dead; + } + +#if defined(MPI_ENABLED) + tag(p) = mpi::SendTag(tag(p), i1(p) < 0, i1(p) >= ni1); +#endif + } + + }; + + + + /* -------------------------------------------------------------------------- */ + /* 1D Pushers */ + /* -------------------------------------------------------------------------- */ + + /** + * massive particle momentum pusher under force-free constraint + */ + template + Inline void Pusher_kernel::ForceFreePush( const coord_t& xp, + const real_t& pp, + const real_t& ex, + real_t& pp_upd) const { + + + real_t pp_mid { pp }; + pp_upd = pp; + + + //printf("Iteration starts.\n"); + for (auto i { 0 }; i < niter; ++i) { + // find midpoint values + pp_mid = 0.5 * (pp + pp_upd); + + // find updated momentum + pp_upd = pp + + dt * (coeff * ex + + HALF * compute_u0(pp_mid, xp) * + (-TWO * metric.alpha(xp) * DERIVATIVE(metric.alpha, xp[0]) + + DERIVATIVE(metric.f2, xp[0]) * SQR((pp_mid / compute_u0(pp_mid, xp) - metric.f1(xp)) / metric.f2(xp)) + + TWO * DERIVATIVE(metric.f1, xp[0]) * (pp_mid / compute_u0(pp_mid, xp) - metric.f1(xp)) / metric.f2(xp) + + DERIVATIVE(metric.f0, xp[0]) + ) + ); + } + } + + /** + * coordinate pusher + */ + template + Inline void Pusher_kernel::CoordinatePush(const coord_t& xp, + const real_t& vp, + coord_t& xp_upd) const { + xp_upd[0] = dt * vp + xp[0]; + } + + /* ------------------------- Massive particle pusher ------------------------ */ +template + Inline void Pusher_kernel::operator()(index_t p) const { + if (tag(p) != ParticleTag::alive) { + if (tag(p) != ParticleTag::dead) { + raise::KernelError(HERE, "Invalid particle tag in pusher"); + } + return; + } + // record previous coordinate + i1_prev(p) = i1(p); + dx1_prev(p) = dx1(p); + + coord_t xp { ZERO }; + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + + real_t Dp_cntrv { ZERO }; + interpolateFields(p, Dp_cntrv); + + real_t ep_cntrv { metric.alpha(xp) * metric.template transform<1, Idx::U, Idx::D>(xp, Dp_cntrv) }; + + real_t pp { px1(p) }; + + + /* -------------------------------- Leapfrog -------------------------------- */ + /* u_i(n - 1/2) -> u_i(n + 1/2) */ + real_t pp_upd { ZERO }; + + ForceFreePush(xp, pp, ep_cntrv, pp_upd); + + + /* x^i(n) -> x^i(n + 1) */ + coord_t xp_upd { ZERO }; + CoordinatePush(xp, (pp_upd / compute_u0(pp_upd, xp) - metric.f1(xp)) / metric.f2(xp), xp_upd); + + /* update coordinate */ + int i1_; + prtldx_t dx1_; + from_Xi_to_i_di(xp_upd[0], i1_, dx1_); + i1(p) = i1_; + dx1(p) = dx1_; + + /* update velocity */ + px1(p) = pp_upd; + + + boundaryConditions(p); + } + + + + + + +} // namespace kernel::gr + + + + + +#undef DERIVATIVE + +#undef i_di_to_Xi +#undef from_Xi_to_i_di +#undef from_Xi_to_i + +#endif // KERNELS_PARTICLE_PUSHER_GR_1D_HPP + diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt index 10e8bb944..70146f06c 100644 --- a/src/kernels/tests/CMakeLists.txt +++ b/src/kernels/tests/CMakeLists.txt @@ -31,3 +31,4 @@ gen_test(fields_to_phys) gen_test(prtls_to_phys) gen_test(gca_pusher) gen_test(prtl_bc) +gen_test(ff_pusher) diff --git a/src/kernels/tests/ff_pusher.cpp b/src/kernels/tests/ff_pusher.cpp new file mode 100644 index 000000000..185e2b0fd --- /dev/null +++ b/src/kernels/tests/ff_pusher.cpp @@ -0,0 +1,245 @@ +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/comparators.h" +#include "utils/numeric.h" + +#include "metrics/boyer_lindq_tp.h" + +#include "kernels/particle_pusher_1D_gr.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + +#define from_Xi_to_i(XI, I) \ + { \ + I = static_cast((XI)); \ + } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define DERIVATIVE(func, x) \ + ((func({ x + epsilon }) - func({ x - epsilon })) / (TWO * epsilon)) + +inline static constexpr auto epsilon = std::numeric_limits::epsilon(); + +Inline auto equal(const real_t& a, + const real_t& b, + const char* msg, + const real_t acc = ONE) -> bool { + const auto eps = math::sqrt(epsilon) * acc; + if (not cmp::AlmostEqual(a, b, eps)) { + printf("%.12e != %.12e %s\n", a, b, msg); + return false; + } + return true; +} + +void errorIf(bool condition, const std::string& message) { + if (condition) { + throw std::runtime_error(message); + } +} + +template +void put_value(array_t& arr, T v, index_t p) { + auto h = Kokkos::create_mirror_view(arr); + Kokkos::deep_copy(h, arr); + h(p) = v; + Kokkos::deep_copy(arr, h); +} + +template +void testFFPusher(const std::vector& res, + const boundaries_t& ext, + const real_t acc, + const std::map& params = {}) { + static_assert(M::Dim == 1, "Only 1D is supported"); + errorIf(res.size() != M::Dim, "res.size() != M::Dim"); + + boundaries_t extent; + extent = ext; + + M metric { res, extent, params }; + + static constexpr auto D = M::Dim; + + const int nx1 = res[0]; + + const real_t d_eta_inv { nx1 / (metric.r2eta(extent[0].second) - + metric.r2eta(extent[0].first)) }; + + const auto coeff { d_eta_inv }; + const auto dt = real_t { 0.1 }; + + const auto range_ext = CreateRangePolicy({ 0 }, + { res[0] + 2 * N_GHOSTS }); + + auto efield = ndfield_t { "efield", res[0] + 2 * N_GHOSTS }; + + array_t i1 { "i1", 30 }; + array_t i1_prev { "i1_prev", 30 }; + array_t dx1 { "dx1", 30 }; + array_t dx1_prev { "dx1_prev", 30 }; + array_t px1 { "px1", 30 }; + array_t tag { "tag", 30 }; + + const real_t sep = 0.1 * res[0]; + real_t x1i { ZERO }; + int ii { 0 }; + prtldx_t dx1i { ZERO }; + + // Particle boundaries + auto boundaries = boundaries_t {}; + boundaries = { + { PrtlBC::PERIODIC, PrtlBC::PERIODIC } + }; + + // No efield + + Kokkos::parallel_for( + "init efield", + range_ext, + Lambda(index_t i1) { efield(i1, em::ex1) = ZERO; }); + + for (size_t n { 0 }; n < 30; n++) { + x1i = ZERO + (n % 10) * sep; + from_Xi_to_i_di(x1i, ii, dx1i) put_value(i1, ii, n); + put_value(dx1, dx1i, n); + put_value(px1, ONE, n); + put_value(tag, ParticleTag::alive, n); + } + + // Without electric field + + Kokkos::parallel_for( + "pusher", + CreateRangePolicy({ 0 }, { 10 }), + // clang-format off + kernel::gr::Pusher_kernel>(efield, + i1, i1_prev, dx1, dx1_prev, px1, + tag, metric, -coeff, dt, nx1, + 1e-2, 30, boundaries)); + // clang-format on + + // With electric field + // positive charge + + Kokkos::parallel_for( + "init efield", + range_ext, + Lambda(index_t i1) { efield(i1, em::ex1) = 99.0; }); + + Kokkos::parallel_for( + "pusher", + CreateRangePolicy({ 10 }, { 20 }), + // clang-format off + kernel::gr::Pusher_kernel>(efield, + i1, i1_prev, dx1, dx1_prev, px1, + tag, metric, coeff, dt, nx1, + 1e-2, 30, boundaries)); + // clang-format on + // negative charge + + Kokkos::parallel_for( + "pusher", + CreateRangePolicy({ 20 }, { 30 }), + // clang-format off + kernel::gr::Pusher_kernel>(efield, + i1, i1_prev, dx1, dx1_prev, px1, + tag, metric, -coeff, dt, nx1, + 1e-2, 30, boundaries)); + // clang-format on + + auto i1_ = Kokkos::create_mirror_view(i1); + Kokkos::deep_copy(i1_, i1); + auto i1_prev_ = Kokkos::create_mirror_view(i1_prev); + Kokkos::deep_copy(i1_prev_, i1_prev); + auto dx1_ = Kokkos::create_mirror_view(dx1); + Kokkos::deep_copy(dx1_, dx1); + auto dx1_prev_ = Kokkos::create_mirror_view(dx1_prev); + Kokkos::deep_copy(dx1_prev_, dx1_prev); + auto px1_ = Kokkos::create_mirror_view(px1); + Kokkos::deep_copy(px1_, px1); + + real_t pp_exp[10] { 1.00552757803207, 1.00905004166760, 1.01264731982748, + 1.01622119133582, 1.01961904791731, 1.02260412478936, + 1.02480734477237, 1.02564651052812, 1.02418169991879, + 1.01916798198126 }; + real_t pp_e_p_exp[10] { 1.12420507263443, 1.18040523823437, 1.27055279565053, + 1.42520265982451, 1.71347897729332, 2.31132787710123, + 3.74396948272974, 7.98290677490210, 25.6627981723171, + 170.710945349036 }; + real_t pp_e_n_exp[10] { 0.886850044532529, 0.837694736928196, + 0.754741504814654, 0.607238475100313, + 0.325753270778111, -0.266166995807237, + -1.67869862572720, -5.82318047383146, + -23.1877218243589, -166.709733472388 }; + + unsigned int wrongs { 0 }; + + for (size_t n { 0 }; n < 10; ++n) { + wrongs += !equal(px1_(n), pp_exp[n], "no efield", acc); + } + + for (size_t n { 0 }; n < 10; ++n) { + wrongs += !equal(px1_(n + 10), pp_e_p_exp[n], "positive charge", acc); + } + + for (size_t n { 0 }; n < 10; ++n) { + wrongs += !equal(px1_(n + 20), pp_e_n_exp[n], "negative charge", acc); + } + + if (wrongs) { + throw std::runtime_error("ff_usher failed."); + } + + // with electric field, positive charge +} + +auto main(int argc, char* argv[]) -> int { + Kokkos::initialize(argc, argv); + + try { + using namespace ntt; + + testFFPusher>( + { + 128 + }, + { { 2.0, 20.0 } }, + 100, + { { "a", (real_t)0.95 }, + { "psi0", (real_t)1.0 }, + { "theta0", (real_t)1.0 }, + { "Omega", (real_t)0.5 } }); + + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + Kokkos::finalize(); + return 1; + } + Kokkos::finalize(); + return 0; +} + +#undef DERIVATIVE +#undef i_di_to_Xi +#undef from_Xi_to_i_di +#undef from_Xi_to_i diff --git a/src/metrics/boyer_lindq_tp.h b/src/metrics/boyer_lindq_tp.h new file mode 100644 index 000000000..78c4880ee --- /dev/null +++ b/src/metrics/boyer_lindq_tp.h @@ -0,0 +1,362 @@ +/** + * @file metrics/boyer_lindq_tp.h + * @brief metric along constant flux surfaces in force-free fields, based on totoised Boyer-Lindquist-psi coordinates + * @implements + * - metric::BoyerLindqTP<> : metric::MetricBase<> + * @namespaces: + * - metric:: + * !TODO + * None radial surfaces needs to be implemented (dpsi_dr != 0). + */ + +#ifndef METRICS_BOYER_LINDQ_TP_H +#define METRICS_BOYER_LINDQ_TP_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/numeric.h" + +#include "metrics/metric_base.h" + +#include +#include +#include + + +namespace metric { + + template + class BoyerLindqTP : public MetricBase { + static_assert(D == Dim::_1D, "Only 1D boyer_lindq_tp is available"); + + private: + // Spin parameter, in [0,1[ + // and horizon size in units of rg + // all physical extents are in units of rg + const real_t a, rg_, psi0, th0; + const real_t rh, rh_m, psi, bt, Omega, dpsi_dth, dbt_dth; + const real_t eta_max, eta_min; + const real_t d_eta, d_eta_inv; + + Inline auto Delta(const real_t& r) const -> real_t { + return SQR(r) - TWO * r + SQR(a); + } + + Inline auto Sigma(const real_t& r) const -> real_t { + return SQR(r) + SQR(a) * SQR(math::cos(th0)); + } + + Inline auto A(const real_t& r) const -> real_t { + return SQR(SQR(r) + SQR(a)) - SQR(a) * Delta(r) * SQR(math::sin(th0)); + } + + Inline auto omega(const real_t& r) const -> real_t { + return TWO * a * r / A(r); + } + + + + + public: + static constexpr const char* Label { "boyer_lindq_tp" }; + static constexpr Dimension PrtlDim { D }; + static constexpr ntt::Coord::type CoordType { ntt::Coord::Bltp }; + static constexpr ntt::Metric::type MetricType { ntt::Metric::BoyerLindqTP }; + using MetricBase::x1_min; + using MetricBase::x1_max; + using MetricBase::nx1; + using MetricBase::set_dxMin; + + BoyerLindqTP(std::vector res, + boundaries_t ext, + const std::map& params) + : MetricBase { res, ext } + , a { params.at("a") } + , psi0 { params.at("psi0") } + , th0 { params.at("theta0") } + , rg_ { ONE } + , rh { ONE + math::sqrt(ONE - SQR(a)) } + , rh_m { ONE - math::sqrt(ONE - SQR(a)) } + , psi { psi0 * (1 - math::cos(th0)) } + , bt { -HALF * psi0 * a * math::sin(th0) * math::cos(th0) / Sigma(rh) } + , Omega { params.at("Omega") * a / (SQR(a) + SQR(rh)) } + , dpsi_dth { psi0 * math::sin(th0) } + , dbt_dth { -HALF * psi0 * a * (SQR(a * math::cos(th0)) + SQR(rh) * math::cos(TWO * th0)) / SQR(Sigma(rh)) } + , eta_min { r2eta(x1_min) } + , eta_max { r2eta(x1_max) } + , d_eta { (eta_max - eta_min) / nx1 } + , d_eta_inv { ONE / d_eta }{ + set_dxMin(find_dxMin()); + } + + ~BoyerLindqTP() = default; + + [[nodiscard]] + Inline auto spin() const -> real_t { + return a; + } + + [[nodiscard]] + Inline auto rhorizon() const -> real_t { + return rh; + } + + [[nodiscard]] + Inline auto rhorizon_minus() const -> real_t { + return rh_m; + } + + [[nodiscard]] + Inline auto rg() const -> real_t { + return rg_; + } + + [[nodiscard]] + Inline auto OmegaF() const -> real_t { + return Omega; + } + + + /** + * lapse function + * @param x coordinate array in code units + */ + Inline auto alpha(const coord_t& xi) const -> real_t { + const real_t r_ { eta2r(xi[0] * d_eta + eta_min) }; + return math::sqrt(Sigma(r_) * Delta(r_) / A(r_)); + } + + /** + * shift vector, only third covariant component is non-zero + * @param x coordinate array in code units + */ + + Inline auto beta(const coord_t& xi) const -> real_t { + return math::sqrt(h<3, 3>(xi)) * omega(eta2r(xi[0] * d_eta + eta_min)); + } + + Inline auto beta3(const coord_t& xi) const -> real_t { + return -omega(eta2r(xi[0] * d_eta + eta_min)); + } + + /** + * @brief Compute helper functions f0,f1,f2 in 1D-GRPIC + */ + Inline auto f2(const coord_t& xi) const -> real_t { + const real_t r_ { eta2r(xi[0] * d_eta + eta_min) }; + return SQR(d_eta) * Sigma(r_) * (Delta(r_) + + A(r_) * SQR(bt / dpsi_dth ) + ); + } + + Inline auto f1(const coord_t& xi) const -> real_t { + const real_t r_ { eta2r(xi[0] * d_eta + eta_min) }; + return d_eta * A(r_) * bt * (Omega + beta3(xi)) / psi0; + } + + Inline auto f0(const coord_t& xi) const -> real_t { + return h_<3, 3>(xi) * SQR(Omega + beta3(xi)); + } + + /** + * @brief force-free charge densities and currents + */ + + Inline auto rho_ff(const coord_t& xi) const -> real_t { + const real_t r_ { eta2r(xi[0] * d_eta + eta_min) }; + return psi0 / sqrt_det_h(xi) * math::sin(TWO * th0) * + ((omega(r_) - Omega) * (ONE + SQR(a * math::sin(th0)) / Sigma(r_)) / SQR(alpha(xi)) + + SQR(a * math::sin(th0)) * Omega / Sigma(r_)); + } + + Inline auto J_ff() const -> real_t { + return -d_eta_inv * HALF * psi0 * a * + (ONE - TWO * SQR(rh * math::sin(th0)) / Sigma(rh)) + / Sigma(rh); + } + + + + /** + * minimum effective cell size for a given metric (in physical units) + */ + [[nodiscard]] + auto find_dxMin() const -> real_t override { + real_t min_dx { -ONE }; + for (int i { 0 }; i < nx1; ++i){ + real_t i_ { static_cast(i) + HALF }; + coord_t xi { i_ }; + real_t dx = ONE / math::sqrt(h<1, 1>(xi)); + if ((min_dx > dx) || (min_dx < 0.0)) { + min_dx = dx; + } + } + return min_dx; + } + + /** + * metric component with lower indices: h_ij + * @param x coordinate array in code units + */ + template + Inline auto h_(const coord_t& x) const -> real_t { + static_assert(i > 0 && i <= 3, "Invalid index i"); + static_assert(j > 0 && j <= 3, "Invalid index j"); + + const real_t r_ { eta2r(x[0] * d_eta + eta_min) }; + if constexpr (i == 1 && j == 1) { + // h_11 + return SQR(d_eta) * Sigma(r_) * Delta(r_); + } else if constexpr (i == 2 && j == 2) { + // h_22 + return Sigma(r_) / SQR(dpsi_dth) ; + } else if constexpr (i == 3 && j == 3) { + // h_33 + return A(r_) * SQR(math::sin(th0)) / Sigma(r_); + }else { + return ZERO; + } + } + + /** + * metric component with upper indices: h^ij + * @param x coordinate array in code units + */ + template + Inline auto h(const coord_t& x) const -> real_t { + static_assert(i > 0 && i <= 3, "Invalid index i"); + static_assert(j > 0 && j <= 3, "Invalid index j"); + + if constexpr (i == j) { + return ONE / h_(x); + }else { + return ZERO; + } + } + + /** + * sqrt(det(h_ij)) + * @param x coordinate array in code units + */ + Inline auto sqrt_det_h(const coord_t& x) const -> real_t { + return math::sqrt(h_<1, 1>(x) * h_<2, 2>(x) * h_<3, 3>(x)); + } + + /** + * coordinate conversions + */ + template + Inline void convert(const coord_t& x_in, coord_t& x_out) const { + static_assert(in != out, "Invalid coordinate conversion"); + static_assert((in != Crd::XYZ) && + (out != Crd::XYZ), + "Invalid coordinate conversion: XYZ not allowed in GR"); + if constexpr (in == Crd::Cd && (out == Crd::Sph || out == Crd::Ph)){ + // code -> sph/phys + x_out[0] = eta2r(x_in[0] * d_eta + eta_min); + }else if constexpr ((in == Crd::Sph || in == Crd::Ph) && out == Crd::Cd){ + //sph/phys -> code + x_out[0] = (r2eta(x_in[0]) - eta_min) / d_eta; + } + } + + + /** + * component wise vector transformations + */ + template + Inline auto transform(const coord_t& x, + const real_t& v_in) const -> real_t { + static_assert(in != out, "Invalid vector transformation"); + static_assert(in != Idx::XYZ && out != Idx::XYZ, + "Invalid vector transformation: XYZ not allowed in GR"); + static_assert(i > 0 && i <= 3, "Invalid index i"); + if constexpr ((in == Idx::T && out == Idx::Sph) || + (in == Idx::Sph && out == Idx::T)) { + // tetrad <-> sph + return v_in; + } else if constexpr (((in == Idx::T || in == Idx::Sph) && out == Idx::U) || + (in == Idx::D && (out == Idx::T || out == Idx::Sph))){ + // tetrad/sph -> cntrv ; cov -> tetrad/sph + if constexpr (i == 1){ + return v_in / math::sqrt(h_(x)); + }else if constexpr (i == 2){ + return v_in / math::sqrt(h_(x)); + }else{ + return v_in / math::sqrt(h_(x)); + } + } else if constexpr ((in == Idx::U && (out == Idx::T || out == Idx::Sph)) || + ((in == Idx::T || in == Idx::Sph) && out == Idx::D)){ + // cntrv -> tetrad/sph ; tetrad/sph -> cov + if constexpr (i == 1){ + return v_in * math::sqrt(h_(x)); + }else if constexpr (i == 2){ + return v_in * math::sqrt(h_(x)); + }else{ + return v_in * math::sqrt(h_(x)); + } + } else if constexpr (in == Idx::U && out == Idx::D) { + // cntrv -> cov + return v_in * h_(x); + } else if constexpr (in == Idx::D && out == Idx::U) { + // cov -> cntrv + return v_in * h(x); + } else if constexpr ((in == Idx::U && out == Idx::PU) || + (in == Idx::PD && out == Idx::D)) { + // cntrv -> phys cntrv || phys cov -> cov + if constexpr (i == 1){ + return v_in * Delta(eta2r(x[0] * d_eta + eta_min)) * d_eta ; + }else if constexpr (i == 2){ + return v_in / dpsi_dth; + }else{ + return v_in; + } + } else if constexpr ((in == Idx::PU && out == Idx::U) || + (in == Idx::D && out == Idx::PD)) { + // phys cntrv -> cntrv || cov -> phys cov + if constexpr (i == 1){ + return v_in * d_eta_inv / Delta(eta2r(x[0] * d_eta + eta_min)); + }else if constexpr (i == 2){ + return v_in * dpsi_dth; + }else{ + return v_in; + } + } else { + raise::KernelError(HERE, "Invalid transformation"); + } + } + + /** + * full vector transformations + */ + template + Inline void transform(const coord_t& xi, + const vec_t& v_in, + vec_t& v_out) const { + static_assert(in != out, "Invalid vector transformation"); + if constexpr (in != Idx::XYZ && out != Idx::XYZ) { + v_out[0] = transform<1, in, out>(xi, v_in[0]); + v_out[1] = transform<2, in, out>(xi, v_in[1]); + v_out[2] = transform<3, in, out>(xi, v_in[2]); + } else { + raise::KernelError(HERE, "Invalid vector transformation"); + } + } + + Inline auto r2eta(const real_t& r) const -> real_t{ + return math::log((r - rh) / (r - rh_m)) / (rh - rh_m); + } + + Inline auto eta2r(const real_t& eta) const -> real_t{ + return rh_m - TWO * math::sqrt(ONE - SQR(a)) / math::expm1(eta * TWO * math::sqrt(ONE - SQR(a))); + } + }; + + + +} + + +#endif // METRICS_BOYER_LINDQ_TP_H diff --git a/src/metrics/tests/CMakeLists.txt b/src/metrics/tests/CMakeLists.txt index c997ab079..410497e25 100644 --- a/src/metrics/tests/CMakeLists.txt +++ b/src/metrics/tests/CMakeLists.txt @@ -28,4 +28,4 @@ gen_test(coord_trans) gen_test(sph-qsph) gen_test(ks-qks) gen_test(sr-cart-sph) - +gen_test(bltp) diff --git a/src/metrics/tests/bltp.cpp b/src/metrics/tests/bltp.cpp new file mode 100644 index 000000000..b5d4326cd --- /dev/null +++ b/src/metrics/tests/bltp.cpp @@ -0,0 +1,199 @@ +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/comparators.h" + +#include "metrics/boyer_lindq_tp.h" + +#include +#include +#include +#include + +void errorIf(bool condition, const std::string& message) { + if (condition) { + throw std::runtime_error(message); + } +} + +inline static constexpr auto epsilon = std::numeric_limits::epsilon(); + +template +Inline auto equal(const vec_t& a, + const vec_t& b, + const char* msg, + const real_t acc = ONE) -> bool { + const auto eps = epsilon * acc; + for (unsigned short d = 0; d < D; ++d) { + if (not cmp::AlmostEqual(a[d], b[d], eps)) { + printf("%d : %.12e != %.12e %s\n", d, a[d], b[d], msg); + return false; + } + } + return true; +} + +Inline auto equal(const real_t& a, + const real_t& b, + const char* msg, + const real_t acc = ONE) -> bool { + const auto eps = epsilon * acc; + if (not cmp::AlmostEqual(a, b, eps)) { + printf("%.12e != %.12e %s\n", a, b, msg); + return false; + } + return true; +} + +template +Inline void unravel(std::size_t idx, + tuple_t& ijk, + const tuple_t& res) { + for (unsigned short d = 0; d < D; ++d) { + ijk[d] = idx % res[d]; + idx /= res[d]; + } +} + +template +void testMetric(const std::vector& res, + const boundaries_t& ext, + const real_t acc = ONE, + const std::map& params = {}) { + static_assert(M::Dim == 1, "Dim != 1"); + errorIf(res.size() != (std::size_t)(M::Dim), "res.size() != M.dim"); + errorIf(ext.size() != (std::size_t)(M::Dim), "ext.size() != M.dim"); + for (const auto& e : ext) { + errorIf(e.first >= e.second, "e.first >= e.second"); + } + + M metric(res, ext, params); + tuple_t res_tup; + std::size_t npts { 1 }; + for (auto d = 0; d < M::Dim; ++d) { + res_tup[d] = res[d]; + npts *= res[d]; + } + + const auto x_min { ext[0].first }; + const auto x_max { ext[0].second }; + + + unsigned long all_wrongs { 0 }; + const auto rg { metric.rg() }; + const auto rh { metric.rhorizon() }; + const auto a { metric.spin() }; + const auto th0 { params.at("theta0") }; + const auto psi0 { params.at("psi0") }; + const auto bt { -HALF * psi0 * a * math::sin(th0) * math::cos(th0) / + (SQR(rh) + SQR(a * math::cos(th0))) }; + const auto Omega { params.at("Omega") * a / (SQR(a) + SQR(rh)) }; + const auto dpsi_dth { -psi0 * math::sin(th0) }; + const auto dbt_dth { -HALF * psi0 * a * (SQR(a * math::cos(th0)) + SQR(rh) * math::cos(TWO * th0)) + / SQR(SQR(rh) + SQR(a * math::cos(th0))) }; + const auto rh_m { rh - TWO * math::sqrt(ONE - SQR(a)) }; + const auto eta_min { math::log((x_min - rh) / (x_min - rh_m)) / (rh - rh_m) }; + const auto eta_max { math::log((x_max - rh) / (x_max - rh_m)) / (rh - rh_m) }; + const auto d_eta { (eta_max - eta_min) / ((real_t)res[0]) }; + + Kokkos::parallel_reduce( + "h_ij/hij", + npts, + Lambda(index_t n, unsigned long& wrongs) { + tuple_t idx; + unravel(n, idx, res_tup); + coord_t x_Code { ZERO }; + coord_t x_Phys { ZERO }; + + for (unsigned short d = 0; d < M::Dim; ++d) { + x_Code[d] = (real_t)(idx[d]) + HALF; + } + + const auto h11 { metric.template h<1, 1>(x_Code) }; + const auto h22 { metric.template h<2, 2>(x_Code) }; + const auto h33 { metric.template h<3, 3>(x_Code) }; + const auto h_11 { metric.template h_<1, 1>(x_Code) }; + const auto h_22 { metric.template h_<2, 2>(x_Code) }; + const auto h_33 { metric.template h_<3, 3>(x_Code) }; + + metric.template convert(x_Code, x_Phys); + const auto r { x_Phys[0] }; + + + const auto Sigma { SQR(r) + SQR(a * math::cos(th0)) }; + const auto Delta { SQR(r) - TWO * rg * r + SQR(a) }; + const auto A { SQR(SQR(r) + SQR(a)) - SQR(a * math::sin(th0)) * Delta }; + const auto omega { TWO * a * r / A }; + + const auto h_11_expect { Sigma / Delta }; + const auto h_22_expect { Sigma }; + const auto h_33_expect { A * SQR(math::sin(th0)) / Sigma }; + const auto h11_expect { ONE / h_11_expect }; + const auto h22_expect { ONE / h_22_expect }; + const auto h33_expect { ONE / h_33_expect }; + + const auto f0_expect { h_33_expect * SQR(Omega - omega) }; + const auto f1_expect { d_eta * A * bt * (Omega - omega) / psi0 }; + const auto f2_expect { SQR(d_eta) * Sigma * (Delta + A * SQR(bt / dpsi_dth)) }; + + const auto f0_predict { metric.f0(x_Code) }; + const auto f1_predict { metric.f1(x_Code) }; + const auto f2_predict { metric.f2(x_Code) }; + + const auto rho_ff_expect { psi0 * math::sin(TWO * th0) / (d_eta * math::sqrt(h_11_expect * h_22_expect * h_33_expect)) + * ( A / (Sigma * Delta) * (omega - Omega) * (ONE + SQR(a * math::sin(th0)) / Sigma) + + SQR(a * math::sin(th0)) / Sigma * Omega) }; + const auto rho_ff_predict { metric.rho_ff(x_Code) }; + + + vec_t hij_temp { ZERO }, hij_predict { ZERO }; + metric.template transform(x_Code, { h11, h22, h33 }, hij_temp); + metric.template transform(x_Code, hij_temp, hij_predict); + + + + vec_t hij_expect { h11_expect, h22_expect, h33_expect }; + + wrongs += not equal(hij_expect, hij_predict, "hij", acc); + wrongs += not equal(f0_expect, f0_predict, "f0", acc); + wrongs += not equal(f1_expect, f1_predict, "f1", acc); + wrongs += not equal(f2_expect, f2_predict, "f2", acc); + + }, + all_wrongs); + + const auto J_ff_expect { -HALF / d_eta * psi0 * a / (SQR(rh) + SQR(a * math::cos(th0))) + * (ONE - TWO * SQR(rh * math::sin(th0)) / (SQR(rh) + SQR(a * math::cos(th0))) )}; + + const auto J_ff_predict { metric.J_ff() }; + + all_wrongs += not equal(J_ff_expect, J_ff_predict, "J_ff", acc); + + errorIf(all_wrongs != 0, + "wrong h_ij/hij for " + std::to_string(M::Dim) + "D " + + std::string(metric.Label) + " with " + std::to_string(all_wrongs) + + " errors"); +} + +auto main(int argc, char* argv[]) -> int { + Kokkos::initialize(argc, argv); + + try { + using namespace ntt; + using namespace metric; + testMetric>( + { 128 }, + { { 2.0, 20.0 } }, + 30, + { { "a", (real_t)0.95 } , + { "psi0", (real_t)1.0 } , + { "theta0", (real_t)1.0 } , + { "Omega", (real_t)0.5 }}); + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + Kokkos::finalize(); + return 1; + } + Kokkos::finalize(); + return 0; +} diff --git a/src/metrics/tests/coord_trans.cpp b/src/metrics/tests/coord_trans.cpp index f3779a852..f507063c1 100644 --- a/src/metrics/tests/coord_trans.cpp +++ b/src/metrics/tests/coord_trans.cpp @@ -9,6 +9,7 @@ #include "metrics/qkerr_schild.h" #include "metrics/qspherical.h" #include "metrics/spherical.h" +#include "metrics/boyer_lindq_tp.h" #include #include @@ -182,6 +183,15 @@ auto main(int argc, char* argv[]) -> int { }; testMetric>(resks0, extks0, 150); + const auto psramsbltp = std::map { + {"a" , (real_t)0.95}, + {"psi0" , ONE}, + {"theta0", ONE}, + {"Omega" , HALF} + }; + + testMetric>({ 128 }, ext1dcart, 200, psramsbltp); + } catch (std::exception& e) { std::cerr << e.what() << std::endl; Kokkos::finalize(); diff --git a/src/metrics/tests/vec_trans.cpp b/src/metrics/tests/vec_trans.cpp index 31015115c..3561b731e 100644 --- a/src/metrics/tests/vec_trans.cpp +++ b/src/metrics/tests/vec_trans.cpp @@ -9,6 +9,7 @@ #include "metrics/qkerr_schild.h" #include "metrics/qspherical.h" #include "metrics/spherical.h" +#include "metrics/boyer_lindq_tp.h" #include #include @@ -236,6 +237,15 @@ auto main(int argc, char* argv[]) -> int { // }; // testMetric>(resks0, extks0, 10); + const auto psramsbltp = std::map { + {"a" , (real_t)0.95}, + {"psi0" , ONE}, + {"theta0", ONE}, + {"Omega" , HALF} + }; + + testMetric>({ 128 }, ext1dcart, 200, psramsbltp); + } catch (std::exception& e) { std::cerr << e.what() << std::endl; Kokkos::finalize();