diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index 8dfe1f38..e12b9d57 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -29,13 +29,15 @@ let getArch = _: if gpu != "NONE" && arch == "NATIVE" then - throw "Please specify an architecture when the GPU support is enabled. Available architectures: https://kokkos.org/kokkos-core-wiki/keywords.html#architectures" + throw "Please specify an architecture when the GPU support is enabled. Available architectures: https://kokkos.org/kokkos-core-wiki/get-started/configuration-guide.html#gpu-architectures" else arch; cmakeExtraFlags = { "HIP" = [ "-D Kokkos_ENABLE_HIP=ON" "-D Kokkos_ARCH_${getArch { }}=ON" + # remove leading AMD_ + "-D AMDGPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}" "-D CMAKE_C_COMPILER=hipcc" "-D CMAKE_CXX_COMPILER=hipcc" ]; diff --git a/pgens/accretion/accretion.toml b/pgens/accretion/accretion.toml new file mode 100644 index 00000000..1ec64143 --- /dev/null +++ b/pgens/accretion/accretion.toml @@ -0,0 +1,85 @@ +[simulation] + name = "wald" + engine = "grpic" + runtime = 500.0 + +[grid] + resolution = [256, 256] + extent = [[1.0, 6.0]] + + [grid.metric] + metric = "qkerr_schild" + qsph_r0 = 0.0 + qsph_h = 0.0 + ks_a = 0.95 + + [grid.boundaries] + fields = [["MATCH"]] + particles = [["ABSORB"]] + + [grid.boundaries.absorb] + ds = 1.0 + +[scales] + larmor0 = 0.025 + skindepth0 = 0.5 + +[algorithms] + current_filters = 4 + + [algorithms.gr] + pusher_niter = 10 + pusher_eps = 1e-2 + + [algorithms.timestep] + CFL = 0.5 + correction = 1.0 + + [algorithms.toggles] + deposit = true + fieldsolver = true + +[particles] + ppc0 = 4.0 + use_weights = true + clear_interval = 100 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 2e8 + pusher = "Boris" + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 2e8 + pusher = "Boris" + +[setup] + multiplicity = 1.0 + sigma_max = 1000.0 + temperature = 0.01 + xi_min = [1.5, 0.0] + xi_max = [4.0, 3.14159265] + m_eps = 1.0 + +[output] + format = "hdf5" + + [output.fields] + interval_time = 1.0 + quantities = ["D", "B", "N_1", "N_2", "A"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[diagnostics] + interval = 2 + colored_stdout = true + blocking_timers = true diff --git a/pgens/accretion/pgen.hpp b/pgens/accretion/pgen.hpp new file mode 100644 index 00000000..864140df --- /dev/null +++ b/pgens/accretion/pgen.hpp @@ -0,0 +1,281 @@ +#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/metadomain.h" + +#include "kernels/particle_moments.hpp" + +namespace user { + using namespace ntt; + + template + struct InitFields { + InitFields(M metric_, real_t m_eps) : metric { metric_ }, m_eps { m_eps } {} + + Inline auto A_3(const coord_t& x_Cd) const -> real_t { + return HALF * (metric.template h_<3, 3>(x_Cd) + + TWO * metric.spin() * metric.template h_<1, 3>(x_Cd) * + metric.beta1(x_Cd)); + } + + Inline auto A_1(const coord_t& x_Cd) const -> real_t { + return HALF * (metric.template h_<1, 3>(x_Cd) + + TWO * metric.spin() * metric.template h_<1, 1>(x_Cd) * + metric.beta1(x_Cd)); + } + + Inline auto A_0(const coord_t& x_Cd) const -> real_t { + real_t g_00 { -metric.alpha(x_Cd) * metric.alpha(x_Cd) + + metric.template h_<1, 1>(x_Cd) * metric.beta1(x_Cd) * + metric.beta1(x_Cd) }; + return HALF * (metric.template h_<1, 3>(x_Cd) * metric.beta1(x_Cd) + + TWO * metric.spin() * g_00); + } + + Inline auto bx1(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + x0m[0] = xi[0]; + x0m[1] = xi[1] - HALF * m_eps; + x0p[0] = xi[0]; + x0p[1] = xi[1] + HALF * m_eps; + + real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + + if (cmp::AlmostZero(x_Ph[1])) { + return ONE; + } else { + return (A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP / m_eps; + } + } + + Inline auto bx2(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + x0m[0] = xi[0] - HALF * m_eps; + x0m[1] = xi[1]; + x0p[0] = xi[0] + HALF * m_eps; + x0p[1] = xi[1]; + + real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + if (cmp::AlmostZero(x_Ph[1])) { + return ZERO; + } else { + return -(A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP / m_eps; + } + } + + Inline auto bx3(const coord_t& x_Ph) const -> real_t { + return ZERO; + } + + Inline auto dx1(const coord_t& x_Ph) const -> real_t { + return ZERO; + } + + Inline auto dx2(const coord_t& x_Ph) const -> real_t { + return ZERO; + } + + Inline auto dx3(const coord_t& x_Ph) const -> real_t { + return ZERO; + } + + private: + const M metric; + const real_t m_eps; + }; + + template + struct PointDistribution : public arch::SpatialDistribution { + PointDistribution(const std::vector& xi_min, + const std::vector& xi_max, + const real_t sigma_thr, + const real_t dens_thr, + const SimulationParams& params, + Domain* domain_ptr) + : arch::SpatialDistribution { domain_ptr->mesh.metric } + , metric { domain_ptr->mesh.metric } + , EM { domain_ptr->fields.em } + , density { domain_ptr->fields.buff } + , sigma_thr { sigma_thr } + , inv_n0 { ONE / params.template get("scales.n0") } + , dens_thr { dens_thr } { + std::copy(xi_min.begin(), xi_min.end(), x_min); + std::copy(xi_max.begin(), xi_max.end(), x_max); + + std::vector specs {}; + for (auto& sp : domain_ptr->species) { + if (sp.mass() > 0) { + specs.push_back(sp.index()); + } + } + + Kokkos::deep_copy(density, ZERO); + auto scatter_buff = Kokkos::Experimental::create_scatter_view(density); + // some parameters + auto& mesh = domain_ptr->mesh; + const auto use_weights = params.template get( + "particles.use_weights"); + const auto ni2 = mesh.n_active(in::x2); + + for (const auto& sp : specs) { + auto& prtl_spec = domain_ptr->species[sp - 1]; + // clang-format off + Kokkos::parallel_for( + "ComputeMoments", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMoments_kernel({}, scatter_buff, 0u, + prtl_spec.i1, prtl_spec.i2, prtl_spec.i3, + prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3, + prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3, + prtl_spec.phi, prtl_spec.weight, prtl_spec.tag, + prtl_spec.mass(), prtl_spec.charge(), + use_weights, + metric, mesh.flds_bc(), + ni2, inv_n0, ZERO)); + // clang-format on + } + Kokkos::Experimental::contribute(density, scatter_buff); + } + + Inline auto sigma_crit(const coord_t& x_Ph) const -> bool { + coord_t xi { ZERO }; + if constexpr (M::Dim == Dim::_2D) { + metric.template convert(x_Ph, xi); + const auto i1 = static_cast(xi[0]) + static_cast(N_GHOSTS); + const auto i2 = static_cast(xi[1]) + static_cast(N_GHOSTS); + const vec_t B_cntrv { EM(i1, i2, em::bx1), + EM(i1, i2, em::bx2), + EM(i1, i2, em::bx3) }; + const vec_t D_cntrv { EM(i1, i2, em::dx1), + EM(i1, i2, em::dx2), + EM(i1, i2, em::dx3) }; + vec_t B_cov { ZERO }; + metric.template transform(xi, B_cntrv, B_cov); + const auto bsqr = + DOT(B_cntrv[0], B_cntrv[1], B_cntrv[2], B_cov[0], B_cov[1], B_cov[2]); + const auto dens = density(i1, i2, 0); + return (bsqr > sigma_thr * dens) || (dens < dens_thr); + } + return false; + } + + Inline auto operator()(const coord_t& x_Ph) const -> real_t override { + auto fill = true; + for (auto d = 0u; d < M::Dim; ++d) { + fill &= x_Ph[d] > x_min[d] and x_Ph[d] < x_max[d] and sigma_crit(x_Ph); + } + return fill ? ONE : ZERO; + } + + private: + tuple_t x_min; + tuple_t x_max; + const real_t sigma_thr; + const real_t dens_thr; + const real_t inv_n0; + Domain* domain_ptr; + ndfield_t density; + ndfield_t EM; + const M metric; + }; + + 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 std::vector xi_min; + const std::vector xi_max; + const real_t sigma0, sigma_max, multiplicity, nGJ, temperature, m_eps; + + InitFields init_flds; + const Metadomain* metadomain; + + inline PGen(SimulationParams& p, const Metadomain& m) + : arch::ProblemGenerator(p) + , xi_min { p.template get>("setup.xi_min") } + , xi_max { p.template get>("setup.xi_max") } + , sigma_max { p.template get("setup.sigma_max") } + , sigma0 { p.template get("scales.sigma0") } + , multiplicity { p.template get("setup.multiplicity") } + , nGJ { p.template get("scales.B0") * + SQR(p.template get("scales.skindepth0")) } + , temperature { p.template get("setup.temperature") } + , m_eps { p.template get("setup.m_eps") } + , init_flds { m.mesh().metric, m_eps } + , metadomain { &m } {} + + inline void InitPrtls(Domain& local_domain) { + const auto energy_dist = arch::Maxwellian(local_domain.mesh.metric, + local_domain.random_pool, + temperature); + const auto spatial_dist = PointDistribution(xi_min, + xi_max, + sigma_max / sigma0, + multiplicity * nGJ, + params, + &local_domain); + + const auto injector = + arch::NonUniformInjector( + energy_dist, + spatial_dist, + { 1, 2 }); + arch::InjectNonUniform(params, + local_domain, + injector, + 1.0, + true); + } + + void CustomPostStep(std::size_t, long double time, Domain& local_domain) { + const auto energy_dist = arch::Maxwellian(local_domain.mesh.metric, + local_domain.random_pool, + temperature); + const auto spatial_dist = PointDistribution(xi_min, + xi_max, + sigma_max / sigma0, + multiplicity * nGJ, + params, + &local_domain); + + const auto injector = + arch::NonUniformInjector( + energy_dist, + spatial_dist, + { 1, 2 }); + arch::InjectNonUniform(params, + local_domain, + injector, + 1.0, + true); + } + }; + +} // namespace user + +#endif diff --git a/pgens/wald/pgen.hpp b/pgens/wald/pgen.hpp new file mode 100644 index 00000000..b2c2aeb0 --- /dev/null +++ b/pgens/wald/pgen.hpp @@ -0,0 +1,209 @@ +#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/formatting.h" +#include "utils/log.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" + +#include + +namespace user { + using namespace ntt; + + template + struct InitFields { + InitFields(M metric_) : metric { metric_ } {} + + Inline auto A_3(const coord_t& x_Cd) const -> real_t { + return HALF * (metric.template h_<3, 3>(x_Cd) + + TWO * metric.spin() * metric.template h_<1, 3>(x_Cd) * + metric.beta1(x_Cd)); + } + + Inline auto A_1(const coord_t& x_Cd) const -> real_t { + return HALF * (metric.template h_<1, 3>(x_Cd) + + TWO * metric.spin() * metric.template h_<1, 1>(x_Cd) * + metric.beta1(x_Cd)); + } + + Inline auto A_0(const coord_t& x_Cd) const -> real_t { + real_t g_00 { -metric.alpha(x_Cd) * metric.alpha(x_Cd) + + metric.template h_<1, 1>(x_Cd) * metric.beta1(x_Cd) * + metric.beta1(x_Cd) }; + return HALF * (metric.template h_<1, 3>(x_Cd) * metric.beta1(x_Cd) + + TWO * metric.spin() * g_00); + } + + Inline auto bx1(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + x0m[0] = xi[0]; + x0m[1] = xi[1] - HALF; + x0p[0] = xi[0]; + x0p[1] = xi[1] + HALF; + + real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + + if (cmp::AlmostZero(x_Ph[1])) { + return ONE; + } else { + return (A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP; + } + } + + Inline auto bx2(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + x0m[0] = xi[0] - HALF; + x0m[1] = xi[1]; + x0p[0] = xi[0] + HALF; + x0p[1] = xi[1]; + + real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + if (cmp::AlmostZero(x_Ph[1])) { + return ZERO; + } else { + return -(A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP; + } + } + + Inline auto bx3( + const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j + HALF ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + x0m[0] = xi[0]; + x0m[1] = xi[1] - HALF; + x0p[0] = xi[0]; + x0p[1] = xi[1] + HALF; + + real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP; + } + + Inline auto dx1(const coord_t& x_Ph) const -> real_t { // at ( i + HALF , j ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + + real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) }; + real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) }; + real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) }; + real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) }; + real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) }; + + // D1 at ( i + HALF , j ) + x0m[0] = xi[0] - HALF; + x0m[1] = xi[1]; + x0p[0] = xi[0] + HALF; + x0p[1] = xi[1]; + real_t E1d { (A_0(x0p) - A_0(x0m)) }; + real_t D1d { E1d / alpha_iPj }; + + // D3 at ( i , j ) + x0m[0] = xi[0] - HALF - HALF; + x0m[1] = xi[1]; + x0p[0] = xi[0] - HALF + HALF; + x0p[1] = xi[1]; + real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij }; + + real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d + + metric.template h<1, 3>({ xi[0], xi[1] }) * D3d }; + + return D1u; + } + + Inline auto dx2(const coord_t& x_Ph) const -> real_t { // at ( i , j + HALF ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + x0m[0] = xi[0]; + x0m[1] = xi[1] - HALF; + x0p[0] = xi[0]; + x0p[1] = xi[1] + HALF; + real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) }; + real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) }; + real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) }; + + real_t E2d { (A_0(x0p) - A_0(x0m)) }; + real_t D2d { E2d / alpha_ijP - (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP }; + real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d }; + + return D2u; + } + + Inline auto dx3(const coord_t& x_Ph) const -> real_t { // at ( i , j ) + coord_t xi { ZERO }, x0m { ZERO }, x0p { ZERO }; + metric.template convert(x_Ph, xi); + real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) }; + real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0], xi[1] }) }; + real_t beta_ij { metric.beta1({ xi[0], xi[1] }) }; + real_t alpha_ij { metric.alpha({ xi[0], xi[1] }) }; + real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) }; + + // D3 at ( i , j ) + x0m[0] = xi[0] - HALF; + x0m[1] = xi[1]; + x0p[0] = xi[0] + HALF; + x0p[1] = xi[1]; + real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij }; + + // D1 at ( i + HALF , j ) + x0m[0] = xi[0] + HALF - HALF; + x0m[1] = xi[1]; + x0p[0] = xi[0] + HALF + HALF; + x0p[1] = xi[1]; + real_t E1d { (A_0(x0p) - A_0(x0m)) }; + real_t D1d { E1d / alpha_iPj }; + + if (cmp::AlmostZero(x_Ph[1])) { + return metric.template h<1, 3>({ xi[0], xi[1] }) * D1d; + } else { + return metric.template h<3, 3>({ xi[0], xi[1] }) * D3d + + metric.template h<1, 3>({ xi[0], xi[1] }) * D1d; + } + } + + private: + const M metric; + }; + + 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; + const Metadomain& global_domain; + + inline PGen(const SimulationParams& p, const Metadomain& m) + : arch::ProblemGenerator { p } + , global_domain { m } + , init_flds { m.mesh().metric } {} + }; + +} // namespace user + +#endif diff --git a/pgens/wald/wald.toml b/pgens/wald/wald.toml new file mode 100644 index 00000000..1104f7a0 --- /dev/null +++ b/pgens/wald/wald.toml @@ -0,0 +1,59 @@ +[simulation] + name = "vacuum" + engine = "grpic" + runtime = 100.0 + +[grid] + resolution = [128, 128] + extent = [[1.0, 10.0]] + + [grid.metric] + metric = "kerr_schild" + qsph_r0 = 0.0 + qsph_h = 0.0 + ks_a = 0.95 + + [grid.boundaries] + fields = [["MATCH"]] + particles = [["ABSORB"]] + + [grid.boundaries.absorb] + ds = 1.0 + +[scales] + larmor0 = 0.0025 + skindepth0 = 0.05 + +[algorithms] + current_filters = 0 + + [algorithms.timestep] + CFL = 0.5 + + [algorithms.toggles] + deposit = false + fieldsolver = true + +[particles] + ppc0 = 2.0 + +[setup] + +[output] + format = "hdf5" + separate_files = false + + [output.fields] + interval_time = 1.0 + quantities = ["D", "H", "B", "A"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[diagnostics] + interval = 2 + colored_stdout = true + blocking_timers = true diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h index 901dbce3..24a3731d 100644 --- a/src/archetypes/energy_dist.h +++ b/src/archetypes/energy_dist.h @@ -83,8 +83,8 @@ namespace arch { , pl_ind { pl_ind } , pool { pool } {} - Inline void operator()(const coord_t& x_Code, - vec_t& v, + Inline void operator()(const coord_t&, + vec_t& v, spidx_t = 0) const override { auto rand_gen = pool.get_state(); auto rand_X1 = Random(rand_gen); @@ -108,15 +108,6 @@ namespace arch { v[1] = v[2] * math::cos(constant::TWO_PI * rand_X3); v[2] = v[2] * math::sin(constant::TWO_PI * rand_X3); - if constexpr (S == SimEngine::GRPIC) { - // convert from the tetrad basis to covariant - vec_t v_Hat; - v_Hat[0] = v[0]; - v_Hat[1] = v[1]; - v_Hat[2] = v[2]; - metric.template transform(x_Code, v_Hat, v); - } - pool.free_state(rand_gen); } @@ -248,9 +239,9 @@ namespace arch { HERE); } - Inline void operator()(const coord_t& x_Code, - vec_t& v, - spidx_t sp = 0) const override { + Inline void operator()(const coord_t&, + vec_t& v, + spidx_t sp = 0) const override { SampleFromMaxwellian(v, pool, temperature, @@ -258,14 +249,6 @@ namespace arch { boost_direction, not zero_current and sp % 2 == 0); - if constexpr (S == SimEngine::GRPIC) { - // convert from the tetrad basis to covariant - vec_t v_Hat; - v_Hat[0] = v[0]; - v_Hat[1] = v[1]; - v_Hat[2] = v[2]; - metric.template transform(x_Code, v_Hat, v); - } } private: @@ -308,9 +291,9 @@ namespace arch { HERE); } - Inline void operator()(const coord_t& x_Code, - vec_t& v, - spidx_t sp = 0) const override { + Inline void operator()(const coord_t&, + vec_t& v, + spidx_t sp = 0) const override { SampleFromMaxwellian( v, pool, @@ -318,14 +301,6 @@ namespace arch { boost_velocity, boost_direction, not zero_current and sp == sp_1); - if constexpr (S == SimEngine::GRPIC) { - // convert from the tetrad basis to covariant - vec_t v_Hat; - v_Hat[0] = v[0]; - v_Hat[1] = v[1]; - v_Hat[2] = v[2]; - metric.template transform(x_Code, v_Hat, v); - } } private: @@ -443,13 +418,6 @@ namespace arch { (drift_dir_x1 + ONE); } } - } else if constexpr (S == SimEngine::GRPIC) { - // convert from the tetrad basis to covariant - vec_t v_Hat; - v_Hat[0] = v[0]; - v_Hat[1] = v[1]; - v_Hat[2] = v[2]; - metric.template transform(x_Code, v_Hat, v); } } diff --git a/src/archetypes/field_setter.h b/src/archetypes/field_setter.h index 281c28df..5c5c4dbe 100644 --- a/src/archetypes/field_setter.h +++ b/src/archetypes/field_setter.h @@ -170,70 +170,30 @@ namespace arch { const real_t x2_H { metric.template convert<2, Crd::Cd, Crd::Ph>( i2_ + HALF) }; { // dx1 - vec_t d_PU { finit.dx1({ x1_H, x2_0 }), - finit.dx2({ x1_H, x2_0 }), - finit.dx3({ x1_H, x2_0 }) }; - vec_t d_U { ZERO }; - metric.template transform({ i1_ + HALF, i2_ }, - d_PU, - d_U); - EM(i1, i2, em::dx1) = d_U[0]; + EM(i1, i2, em::dx1) = finit.dx1({ x1_H, x2_0 }); } { // dx2 - vec_t d_PU { finit.dx1({ x1_0, x2_H }), - finit.dx2({ x1_0, x2_H }), - finit.dx3({ x1_0, x2_H }) }; - vec_t d_U { ZERO }; - metric.template transform({ i1_, i2_ + HALF }, - d_PU, - d_U); - EM(i1, i2, em::dx2) = d_U[1]; + EM(i1, i2, em::dx2) = finit.dx2({ x1_0, x2_H }); } { // dx3 - vec_t d_PU { finit.dx1({ x1_0, x2_0 }), - finit.dx2({ x1_0, x2_0 }), - finit.dx3({ x1_0, x2_0 }) }; - vec_t d_U { ZERO }; - metric.template transform({ i1_, i2_ }, d_PU, d_U); - EM(i1, i2, em::dx3) = d_U[2]; + EM(i1, i2, em::dx3) = finit.dx3({ x1_0, x2_0 }); } } if constexpr (defines_bx1 && defines_bx2 && defines_bx3) { const real_t x1_0 { metric.template convert<1, Crd::Cd, Crd::Ph>(i1_) }; const real_t x1_H { metric.template convert<1, Crd::Cd, Crd::Ph>( i1_ + HALF) }; - const real_t x2_0 { metric.template convert<1, Crd::Cd, Crd::Ph>(i2_) }; - const real_t x2_H { metric.template convert<1, Crd::Cd, Crd::Ph>( + const real_t x2_0 { metric.template convert<2, Crd::Cd, Crd::Ph>(i2_) }; + const real_t x2_H { metric.template convert<2, Crd::Cd, Crd::Ph>( i2_ + HALF) }; { // bx1 - vec_t b_PU { finit.dx1({ x1_0, x2_H }), - finit.dx2({ x1_0, x2_H }), - finit.dx3({ x1_0, x2_H }) }; - vec_t b_U { ZERO }; - metric.template transform({ i1_, i2_ + HALF }, - b_PU, - b_U); - EM(i1, i2, em::bx1) = b_U[0]; + EM(i1, i2, em::bx1) = finit.bx1({ x1_0, x2_H }); } { // bx2 - vec_t b_PU { finit.dx1({ x1_H, x2_0 }), - finit.dx2({ x1_H, x2_0 }), - finit.dx3({ x1_H, x2_0 }) }; - vec_t b_U { ZERO }; - metric.template transform({ i1_ + HALF, i2_ }, - b_PU, - b_U); - EM(i1, i2, em::bx2) = b_U[1]; + EM(i1, i2, em::bx2) = finit.bx2({ x1_H, x2_0 }); } { // bx3 - vec_t b_PU { finit.dx1({ x1_H, x2_H }), - finit.dx2({ x1_H, x2_H }), - finit.dx3({ x1_H, x2_H }) }; - vec_t b_U { ZERO }; - metric.template transform({ i1_ + HALF, i2_ + HALF }, - b_PU, - b_U); - EM(i1, i2, em::bx3) = b_U[2]; + EM(i1, i2, em::bx3) = finit.bx3({ x1_H, x2_H }); } } } else { diff --git a/src/archetypes/tests/powerlaw.cpp b/src/archetypes/tests/powerlaw.cpp index 58df1f4c..3cb76763 100644 --- a/src/archetypes/tests/powerlaw.cpp +++ b/src/archetypes/tests/powerlaw.cpp @@ -45,10 +45,12 @@ struct Caller { raise::KernelError(HERE, "Gamma out of bounds"); } } else { - vec_t vup { ZERO }; - metric.template transform(xp, vp, vup); + vec_t vd { ZERO }; + vec_t vu { ZERO }; + metric.template transform(xp, vp, vd); + metric.template transform(xp, vp, vu); const auto gamma = math::sqrt( - ONE + vup[0] * vp[0] + vup[1] * vp[1] + vup[2] * vp[2]); + ONE + vu[0] * vd[0] + vu[1] * vd[1] + vu[2] * vd[2]); if (gamma < 10 or gamma > 1000) { raise::KernelError(HERE, "Gamma out of bounds"); } diff --git a/src/checkpoint/tests/checkpoint-mpi.cpp b/src/checkpoint/tests/checkpoint-mpi.cpp index bc6d6038..2372d81b 100644 --- a/src/checkpoint/tests/checkpoint-mpi.cpp +++ b/src/checkpoint/tests/checkpoint-mpi.cpp @@ -20,7 +20,7 @@ using namespace checkpoint; void cleanup() { namespace fs = std::filesystem; - fs::path temp_path { "checkpoints" }; + fs::path temp_path { "chck" }; fs::remove_all(temp_path); } @@ -126,11 +126,12 @@ auto main(int argc, char* argv[]) -> int { } adios2::ADIOS adios; + const path_t checkpoint_path { "chck" }; { // write checkpoint Writer writer; - writer.init(&adios, 0, 0.0, 1); + writer.init(&adios, checkpoint_path, 0, 0.0, 1); writer.defineFieldVariables(SimEngine::GRPIC, { g_nx1_gh, g_nx2_gh }, @@ -190,7 +191,7 @@ auto main(int argc, char* argv[]) -> int { array_t plds1_read { "plds_1", npart1, 3 }; adios2::IO io = adios.DeclareIO("checkpointRead"); - adios2::Engine reader = io.Open("checkpoints/step-00000000.bp", + adios2::Engine reader = io.Open(checkpoint_path / "step-00000000.bp", adios2::Mode::Read); reader.BeginStep(); diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index d082d617..78966a25 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:: */ @@ -13,7 +13,11 @@ #define ENGINES_GRPIC_GRPIC_H #include "enums.h" +#include "global.h" +#include "arch/kokkos_aliases.h" +#include "utils/log.h" +#include "utils/numeric.h" #include "utils/timer.h" #include "utils/toml.h" @@ -21,26 +25,1163 @@ #include "framework/parameters.h" #include "engines/engine.hpp" +#include "kernels/ampere_gr.hpp" +#include "kernels/aux_fields_gr.hpp" +#include "kernels/currents_deposit.hpp" +#include "kernels/digital_filter.hpp" +#include "kernels/faraday_gr.hpp" +#include "kernels/fields_bcs.hpp" +#include "kernels/particle_pusher_gr.hpp" +#include "pgen.hpp" + +#include +#include + +#include +#include namespace ntt { + enum class gr_getE { + D0_B, + D_B0 + }; + enum class gr_getH { + D_B0, + D0_B0 + }; + enum class gr_faraday { + aux, + main + }; + enum class gr_ampere { + init, + aux, + main + }; + enum class gr_bc { + main, + aux, + curr + }; + template class GRPICEngine : public Engine { - using base_t = Engine; - - using Engine::m_params; - using Engine::m_metadomain; + 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; public: - static constexpr auto S { SimEngine::SRPIC }; + static constexpr auto S { SimEngine::GRPIC }; - GRPICEngine(const SimulationParams& params) : base_t { params } {} + GRPICEngine(SimulationParams& params) : base_t { params } {} ~GRPICEngine() = default; - void step_forward(timer::Timers&, Domain&) override {} - }; + void step_forward(timer::Timers& timers, domain_t& dom) override { + const auto fieldsolver_enabled = m_params.template get( + "algorithms.toggles.fieldsolver"); + const auto deposit_enabled = m_params.template get( + "algorithms.toggles.deposit"); + const auto clear_interval = m_params.template get( + "particles.clear_interval"); + + if (step == 0) { + if (fieldsolver_enabled) { + // communicate fields and apply BCs on the first timestep + /** + * Initially: em0::B -- + * em0::D -- + * em::B at -1/2 + * em::D at -1/2 + * + * cur0::J -- + * cur::J -- + * + * aux::E -- + * aux::H -- + * + * x_prtl at -1/2 + * u_prtl at -1/2 + */ + + /** + * em0::D, em::D, em0::B, em::B <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, + Comm::B | Comm::B0 | Comm::D | Comm::D0); + FieldBoundaries(dom, BC::B | BC::D, gr_bc::main); + + /** + * em0::B <- em::B + * em0::D <- em::D + * + * Now: em0::B & em0::D at -1/2 + */ + CopyFields(dom); + + /** + * aux::E <- alpha * em::D + beta x em0::B + * aux::H <- alpha * em::B0 - beta x em::D + * + * Now: aux::E & aux::H at -1/2 + */ + ComputeAuxE(dom, gr_getE::D_B0); + ComputeAuxH(dom, gr_getH::D_B0); + + /** + * aux::E, aux::H <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::H | Comm::E); + FieldBoundaries(dom, BC::H | BC::E, gr_bc::aux); + + /** + * em0::B <- (em0::B) <- -curl aux::E + * + * Now: em0::B at 0 + */ + Faraday(dom, gr_faraday::aux, HALF); + + /** + * em0::B, em::B <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::B | Comm::B0); + FieldBoundaries(dom, BC::B, gr_bc::main); + + /** + * em::D <- (em0::D) <- curl aux::H + * + * Now: em::D at 0 + */ + Ampere(dom, gr_ampere::init, HALF); + + /** + * em0::D, em::D <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::D | Comm::D0); + FieldBoundaries(dom, BC::D, gr_bc::main); + + /** + * aux::E <- alpha * em::D + beta x em0::B + * aux::H <- alpha * em0::B - beta x em::D + * + * Now: aux::E & aux::H at 0 + */ + ComputeAuxE(dom, gr_getE::D_B0); + ComputeAuxH(dom, gr_getH::D_B0); + + /** + * aux::E, aux::H <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::H | Comm::E); + FieldBoundaries(dom, BC::H | BC::E, gr_bc::aux); + + // !ADD: GR -- particles? + + /** + * em0::B <- (em::B) <- -curl aux::E + * + * Now: em0::B at 1/2 + */ + Faraday(dom, gr_faraday::main, ONE); + /** + * em0::B, em::B <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::B | Comm::B0); + FieldBoundaries(dom, BC::B, gr_bc::main); + + /** + * em0::D <- (em0::D) <- curl aux::H + * + * Now: em0::D at 1/2 + */ + Ampere(dom, gr_ampere::aux, ONE); + /** + * em0::D, em::D <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::D | Comm::D0); + FieldBoundaries(dom, BC::D, gr_bc::main); + + /** + * aux::H <- alpha * em0::B - beta x em0::D + * + * Now: aux::H at 1/2 + */ + ComputeAuxH(dom, gr_getH::D0_B0); + /** + * aux::H <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::H); + FieldBoundaries(dom, BC::H, gr_bc::aux); + + /** + * em0::D <- (em::D) <- curl aux::H + * + * Now: em0::D at 1 + * em::D at 0 + */ + Ampere(dom, gr_ampere::main, ONE); + /** + * em0::D, em::D <- boundary conditions + */ + m_metadomain.CommunicateFields(dom, Comm::D | Comm::D0); + FieldBoundaries(dom, BC::D, gr_bc::main); + + /** + * em::D <-> em0::D + * em::B <-> em0::B + * em::J <-> em0::J + */ + SwapFields(dom); + /** + * Finally: em0::B at -1/2 + * em0::D at 0 + * em::B at 1/2 + * em::D at 1 + * + * cur0::J -- + * cur::J -- + * + * aux::E -- + * aux::H -- + * + * x_prtl at 1 + * u_prtl at 1/2 + */ + } else { + /** + * em0::B <- em::B + * em0::D <- em::D + * + * Now: em0::B & em0::D at -1/2 + */ + CopyFields(dom); + } + } + + /** + * Initially: em0::B at n-3/2 + * em0::D at n-1 + * em::B at n-1/2 + * em::D at n + * + * cur0::J -- + * cur::J at n-1/2 + * + * aux::E -- + * aux::H -- + * + * x_prtl at n + * u_prtl at n-1/2 + */ + + if (fieldsolver_enabled) { + timers.start("FieldSolver"); + /** + * em0::D <- (em0::D + em::D) / 2 + * em0::B <- (em0::B + em::B) / 2 + * + * Now: em0::D at n-1/2 + * em0::B at n-1 + */ + TimeAverageDB(dom); + /** + * aux::E <- alpha * em0::D + beta x em::B + * + * Now: aux::E at n-1/2 + */ + ComputeAuxE(dom, gr_getE::D0_B); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::E); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + /** + * aux::E <- boundary conditions + */ + FieldBoundaries(dom, BC::E, gr_bc::aux); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * em0::B <- (em0::B) <- -curl aux::E + * + * Now: em0::B at n + */ + Faraday(dom, gr_faraday::aux, ONE); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::B | Comm::B0); + timers.stop("Communications"); + /** + * em0::B, em::B <- boundary conditions + */ + timers.start("FieldBoundaries"); + FieldBoundaries(dom, BC::B, gr_bc::main); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * aux::H <- alpha * em0::B - beta x em::D + * + * Now: aux::H at n + */ + ComputeAuxH(dom, gr_getH::D_B0); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::H); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + /** + * aux::H <- boundary conditions + */ + FieldBoundaries(dom, BC::H, gr_bc::aux); + timers.stop("FieldBoundaries"); + } + + { + /** + * x_prtl, u_prtl <- em::D, em0::B + * + * Now: x_prtl at n + 1, u_prtl at n + 1/2 + */ + timers.start("ParticlePusher"); + ParticlePush(dom); + timers.stop("ParticlePusher"); + + /** + * cur0::J <- current deposition + * + * Now: cur0::J at n+1/2 + */ + if (deposit_enabled) { + timers.start("CurrentDeposit"); + Kokkos::deep_copy(dom.fields.cur0, ZERO); + CurrentsDeposit(dom); + timers.stop("CurrentDeposit"); + + timers.start("Communications"); + m_metadomain.SynchronizeFields(dom, Comm::J); + m_metadomain.CommunicateFields(dom, Comm::J); + timers.stop("Communications"); + + timers.start("FieldBoundaries"); + FieldBoundaries(dom, BC::J, gr_bc::curr); + timers.stop("FieldBoundaries"); + + timers.start("CurrentFiltering"); + CurrentsFilter(dom); + timers.stop("CurrentFiltering"); + } + + timers.start("Communications"); + m_metadomain.CommunicateParticles(dom); + timers.stop("Communications"); + } + + if (fieldsolver_enabled) { + timers.start("FieldSolver"); + if (deposit_enabled) { + /** + * cur::J <- (cur0::J + cur::J) / 2 + * + * Now: cur::J at n + */ + TimeAverageJ(dom); + } + /** + * aux::Е <- alpha * em::D + beta x em0::B + * + * Now: aux::Е at n + */ + ComputeAuxE(dom, gr_getE::D_B0); + timers.stop("FieldSolver"); + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::E); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + /** + * aux::Е <- boundary conditions + */ + FieldBoundaries(dom, BC::E, gr_bc::aux); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * em0::B <- (em::B) <- -curl aux::E + * + * Now: em0::B at n+1/2 + * em::B at n-1/2 + */ + Faraday(dom, gr_faraday::main, ONE); + timers.stop("FieldSolver"); + + /** + * em0::B, em::B <- boundary conditions + */ + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::B | Comm::B0); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + FieldBoundaries(dom, BC::B, gr_bc::main); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * em0::D <- (em0::D) <- curl aux::H + * + * Now: em0::D at n+1/2 + */ + Ampere(dom, gr_ampere::aux, ONE); + timers.stop("FieldSolver"); + + if (deposit_enabled) { + timers.start("FieldSolver"); + /** + * em0::D <- (em0::D) <- cur::J + * + * Now: em0::D at n+1/2 + */ + AmpereCurrents(dom, gr_ampere::aux); + timers.stop("FieldSolver"); + } + + /** + * em0::D, em::D <- boundary conditions + */ + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::D | Comm::D0); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + FieldBoundaries(dom, BC::D, gr_bc::main); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * aux::H <- alpha * em0::B - beta x em0::D + * + * Now: aux::H at n+1/2 + */ + ComputeAuxH(dom, gr_getH::D0_B0); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::H); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + /** + * aux::H <- boundary conditions + */ + FieldBoundaries(dom, BC::B, gr_bc::aux); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + /** + * em0::D <- (em::D) <- curl aux::H + * + * Now: em0::D at n+1 + * em::D at n + */ + Ampere(dom, gr_ampere::main, ONE); + timers.stop("FieldSolver"); + + if (deposit_enabled) { + timers.start("FieldSolver"); + /** + * em0::D <- (em0::D) <- cur0::J + * + * Now: em0::D at n+1 + */ + AmpereCurrents(dom, gr_ampere::main); + timers.stop("FieldSolver"); + } + timers.start("FieldSolver"); + /** + * em::D <-> em0::D + * em::B <-> em0::B + * cur::J <-> cur0::J + */ + SwapFields(dom); + timers.stop("FieldSolver"); + + /** + * em0::D, em::D <- boundary conditions + */ + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::D | Comm::D0); + timers.stop("Communications"); + timers.start("FieldBoundaries"); + FieldBoundaries(dom, BC::D, gr_bc::main); + timers.stop("FieldBoundaries"); + } + + if (clear_interval > 0 and step % clear_interval == 0 and step > 0) { + timers.start("PrtlClear"); + m_metadomain.RemoveDeadParticles(dom); + timers.stop("PrtlClear"); + } + + /** + * Finally: em0::B at n-1/2 + * em0::D at n + * em::B at n+1/2 + * em::D at n+1 + * + * cur0::J (at n) + * cur::J at n+1/2 + * + * aux::E (at n+1/2) + * aux::H (at n) + * + * x_prtl at n+1 + * u_prtl at n+1/2 + */ + } + + /* algorithm substeps --------------------------------------------------- */ + void FieldBoundaries(domain_t& domain, BCTags tags, const gr_bc& g) { + if (g == gr_bc::main) { + for (auto& direction : dir::Directions::orth) { + if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::MATCH) { + MatchFieldsIn(direction, domain, tags, g); + } else if (domain.mesh.flds_bc_in(direction) == FldsBC::AXIS) { + AxisFieldsIn(direction, domain, tags); + } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CUSTOM) { + CustomFieldsIn(direction, domain, tags, g); + } else if (domain.mesh.flds_bc_in(direction) == FldsBC::HORIZON) { + HorizonFieldsIn(direction, domain, tags, g); + } + } // loop over directions + } else if (g == gr_bc::aux) { + for (auto& direction : dir::Directions::orth) { + if (domain.mesh.flds_bc_in(direction) == FldsBC::HORIZON) { + HorizonFieldsIn(direction, domain, tags, g); + } + } + } else if (g == gr_bc::curr) { + for (auto& direction : dir::Directions::orth) { + if (domain.mesh.prtl_bc_in(direction) == PrtlBC::ABSORB) { + MatchFieldsIn(direction, domain, tags, g); + } + } + } + } + + void MatchFieldsIn(dir::direction_t direction, + domain_t& domain, + BCTags tags, + const gr_bc& g) { + /** + * match boundaries + */ + const auto ds_array = m_params.template get>( + "grid.boundaries.match.ds"); + const auto dim = direction.get_dim(); + real_t xg_min, xg_max, xg_edge; + auto sign = direction.get_sign(); + + raise::ErrorIf(((dim != in::x1) or (sign < 0)) and (g == gr_bc::curr), + "Absorption of currents only possible in +x1 (+r)", + HERE); + + real_t ds; + if (sign > 0) { // + direction + ds = ds_array[(short)dim].second; + xg_max = m_metadomain.mesh().extent(dim).second; + xg_min = xg_max - ds; + xg_edge = xg_max; + } else { // - direction + ds = ds_array[(short)dim].first; + xg_min = m_metadomain.mesh().extent(dim).first; + xg_max = xg_min + ds; + xg_edge = xg_min; + } + boundaries_t box; + boundaries_t incl_ghosts; + for (unsigned short d { 0 }; d < M::Dim; ++d) { + if (d == static_cast(dim)) { + box.push_back({ xg_min, xg_max }); + incl_ghosts.push_back({ false, true }); + } else { + box.push_back(Range::All); + incl_ghosts.push_back({ true, true }); + } + } + if (not domain.mesh.Intersects(box)) { + return; + } + const auto intersect_range = domain.mesh.ExtentToRange(box, incl_ghosts); + tuple_t range_min { 0 }; + tuple_t range_max { 0 }; + + for (unsigned short d { 0 }; d < M::Dim; ++d) { + range_min[d] = intersect_range[d].first; + range_max[d] = intersect_range[d].second; + } + if (dim == in::x1) { + if (g != gr_bc::curr) { + Kokkos::parallel_for( + "MatchBoundaries", + CreateRangePolicy(range_min, range_max), + kernel::bc::MatchBoundaries_kernel( + domain.fields.em, + m_pgen.init_flds, + domain.mesh.metric, + xg_edge, + ds, + tags)); + Kokkos::parallel_for( + "MatchBoundaries", + CreateRangePolicy(range_min, range_max), + kernel::bc::MatchBoundaries_kernel( + domain.fields.em0, + m_pgen.init_flds, + domain.mesh.metric, + xg_edge, + ds, + tags)); + } else { + Kokkos::parallel_for( + "AbsorbCurrents", + CreateRangePolicy(range_min, range_max), + kernel::bc::gr::AbsorbCurrents_kernel(domain.fields.cur0, + domain.mesh.metric, + xg_edge, + ds)); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } + + // void AbsorbCurrentsIn(dir::direction_t direction, + // domain_t& domain, + // BCTags tags) { + // /** + // * absorbing currents at the boundaries + // */ + // // @TODO: also in MatchFieldsIn (need a better way) + // const auto ds_array = m_params.template get>( + // "grid.boundaries.match.ds"); + // const auto dim = in::x1; + // real_t xg_min, xg_max, xg_edge; + // real_t ds; + // ds = ds_array[(short)dim].second; + // xg_max = m_metadomain.mesh().extent(dim).second; + // xg_min = xg_max - ds; + // xg_edge = xg_max; + // + // boundaries_t box; + // boundaries_t incl_ghosts; + // for (unsigned short d { 0 }; d < M::Dim; ++d) { + // if (d == static_cast(dim)) { + // box.push_back({ xg_min, xg_max }); + // incl_ghosts.push_back({ false, true }); + // } else { + // box.push_back(Range::All); + // incl_ghosts.push_back({ true, true }); + // } + // } + // if (not domain.mesh.Intersects(box)) { + // return; + // } + // const auto intersect_range = domain.mesh.ExtentToRange(box, incl_ghosts); + // tuple_t range_min { 0 }; + // tuple_t range_max { 0 }; + // + // for (unsigned short d { 0 }; d < M::Dim; ++d) { + // range_min[d] = intersect_range[d].first; + // range_max[d] = intersect_range[d].second; + // } + // Kokkos::parallel_for( + // "AbsorbCurrentsGR", + // CreateRangePolicy(range_min, range_max), + // kernel::bc::AbsorbCurrentsGR_kernel(domain.fields.cur0, + // domain.mesh.metric, + // xg_edge, + // ds)); + // } + + void HorizonFieldsIn(dir::direction_t direction, + domain_t& domain, + BCTags tags, + const gr_bc& g) { + /** + * open boundaries + */ + raise::ErrorIf(M::CoordType == Coord::Cart, + "Invalid coordinate type for horizon BCs", + HERE); + raise::ErrorIf(direction.get_dim() != in::x1, + "Invalid horizon direction, should be x1", + HERE); + const auto i1_min = domain.mesh.i_min(in::x1); + auto range = CreateRangePolicy({ domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x2) + 1 }); + const auto nfilter = m_params.template get( + "algorithms.current_filters"); + if (g == gr_bc::main) { + Kokkos::parallel_for( + "OpenBCFields", + range, + kernel::bc::gr::HorizonBoundaries_kernel(domain.fields.em, + i1_min, + tags, + nfilter)); + Kokkos::parallel_for( + "OpenBCFields", + range, + kernel::bc::gr::HorizonBoundaries_kernel(domain.fields.em0, + i1_min, + tags, + nfilter)); + } + } + + void AxisFieldsIn(dir::direction_t direction, + domain_t& domain, + BCTags tags) { + /** + * axis boundaries + */ + raise::ErrorIf(M::CoordType == Coord::Cart, + "Invalid coordinate type for axis BCs", + HERE); + raise::ErrorIf(direction.get_dim() != in::x2, + "Invalid axis direction, should be x2", + HERE); + const auto i2_min = domain.mesh.i_min(in::x2); + const auto i2_max = domain.mesh.i_max(in::x2); + if (direction.get_sign() < 0) { + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em, + i2_min, + tags)); + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em0, + i2_min, + tags)); + } else { + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em, + i2_max, + tags)); + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em0, + i2_max, + tags)); + } + } + + void CustomFieldsIn(dir::direction_t direction, + domain_t& domain, + BCTags tags, + const gr_bc& g) { + (void)direction; + (void)domain; + (void)tags; + (void)g; + raise::Error("Custom boundaries not implemented", HERE); + // if constexpr ( + // traits::has_member::value) { + // const auto [box, custom_fields] = m_pgen.CustomFields(time); + // if (domain.mesh.Intersects(box)) { + // } + // + // } else { + // raise::Error("Custom boundaries not implemented", HERE); + // } + } + + /** + * @brief Swaps em and em0 fields, cur and cur0 currents. + */ + void SwapFields(domain_t& domain) { + std::swap(domain.fields.em, domain.fields.em0); + std::swap(domain.fields.cur, domain.fields.cur0); + } + + /** + * @brief Copies em fields into em0 + */ + void CopyFields(domain_t& domain) { + Kokkos::deep_copy(domain.fields.em0, domain.fields.em); + } + + void ComputeAuxE(domain_t& domain, const gr_getE& g) { + auto range = range_with_axis_BCs(domain); + if (g == gr_getE::D0_B) { + Kokkos::parallel_for( + "ComputeAuxE", + range, + kernel::gr::ComputeAuxE_kernel(domain.fields.em0, // D + domain.fields.em, // B + domain.fields.aux, // E + domain.mesh.metric)); + } else if (g == gr_getE::D_B0) { + Kokkos::parallel_for("ComputeAuxE", + range, + kernel::gr::ComputeAuxE_kernel(domain.fields.em, + domain.fields.em0, + domain.fields.aux, + domain.mesh.metric)); + } else { + raise::Error("Wrong option for `g`", HERE); + } + } + + void ComputeAuxH(domain_t& domain, const gr_getH& g) { + auto range = range_with_axis_BCs(domain); + if (g == gr_getH::D_B0) { + Kokkos::parallel_for( + "ComputeAuxH", + range, + kernel::gr::ComputeAuxH_kernel(domain.fields.em, // D + domain.fields.em0, // B + domain.fields.aux, // H + domain.mesh.metric)); + } else if (g == gr_getH::D0_B0) { + Kokkos::parallel_for("ComputeAuxH", + range, + kernel::gr::ComputeAuxH_kernel(domain.fields.em0, + domain.fields.em0, + domain.fields.aux, + domain.mesh.metric)); + } else { + raise::Error("Wrong option for `g`", HERE); + } + } + + auto range_with_axis_BCs(const domain_t& domain) -> range_t { + auto range = domain.mesh.rangeActiveCells(); + /** + * @brief taking one extra cell in the x1 and x2 directions if AXIS BCs + */ + if constexpr (M::Dim == Dim::_2D) { + if (domain.mesh.flds_bc_in({ 0, +1 }) == FldsBC::AXIS) { + range = CreateRangePolicy( + { domain.mesh.i_min(in::x1) - 1, domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + } + } else if constexpr (M::Dim == Dim::_3D) { + raise::Error("Invalid dimension", HERE); + } + return range; + } + + void Faraday(domain_t& domain, const gr_faraday& g, real_t fraction = ONE) { + logger::Checkpoint("Launching Faraday kernel", HERE); + const auto dT = fraction * + m_params.template get( + "algorithms.timestep.correction") * + dt; + if (g == gr_faraday::aux) { + Kokkos::parallel_for( + "Faraday", + domain.mesh.rangeActiveCells(), + kernel::gr::Faraday_kernel(domain.fields.em0, // Bin + domain.fields.em0, // Bout + domain.fields.aux, // E + domain.mesh.metric, + dT, + domain.mesh.n_active(in::x2), + domain.mesh.flds_bc())); + } else if (g == gr_faraday::main) { + Kokkos::parallel_for( + "Faraday", + domain.mesh.rangeActiveCells(), + kernel::gr::Faraday_kernel(domain.fields.em, + domain.fields.em0, + domain.fields.aux, + domain.mesh.metric, + dT, + domain.mesh.n_active(in::x2), + domain.mesh.flds_bc())); + + } else { + raise::Error("Wrong option for `g`", HERE); + } + } + + void Ampere(domain_t& domain, const gr_ampere& g, real_t fraction = ONE) { + logger::Checkpoint("Launching Ampere kernel", HERE); + const auto dT = fraction * + m_params.template get( + "algorithms.timestep.correction") * + dt; + auto range = CreateRangePolicy( + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + const auto ni2 = domain.mesh.n_active(in::x2); + + if (g == gr_ampere::aux) { + // First push, updates D0 with J. + Kokkos::parallel_for("Ampere-1", + range, + kernel::gr::Ampere_kernel(domain.fields.em0, // Din + domain.fields.em0, // Dout + domain.fields.aux, + domain.mesh.metric, + dT, + ni2, + domain.mesh.flds_bc())); + } else if (g == gr_ampere::main) { + // Second push, updates D with J0 but assigns it to D0. + Kokkos::parallel_for("Ampere-2", + range, + kernel::gr::Ampere_kernel(domain.fields.em, + domain.fields.em0, + domain.fields.aux, + domain.mesh.metric, + dT, + ni2, + domain.mesh.flds_bc())); + } else if (g == gr_ampere::init) { + // Second push, updates D with J0 and assigns it to D. + Kokkos::parallel_for("Ampere-3", + range, + kernel::gr::Ampere_kernel(domain.fields.em, + domain.fields.em, + domain.fields.aux, + domain.mesh.metric, + dT, + ni2, + domain.mesh.flds_bc())); + } else { + raise::Error("Wrong option for `g`", HERE); + } + } + + void AmpereCurrents(domain_t& domain, const gr_ampere& g) { + logger::Checkpoint("Launching Ampere kernel for adding currents", HERE); + const auto q0 = m_params.template get("scales.q0"); + const auto B0 = m_params.template get("scales.B0"); + const auto coeff = -dt * q0 / B0; + auto range = CreateRangePolicy( + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + const auto ni2 = domain.mesh.n_active(in::x2); + + if (g == gr_ampere::aux) { + // Updates D0 with J: D0(n-1/2) -> (J(n)) -> D0(n+1/2) + Kokkos::parallel_for( + "AmpereCurrentsAux", + range, + kernel::gr::CurrentsAmpere_kernel(domain.fields.em0, + domain.fields.cur, + domain.mesh.metric, + coeff, + ni2, + domain.mesh.flds_bc())); + } else if (g == gr_ampere::main) { + // Updates D0 with J0: D0(n) -> (J0(n+1/2)) -> D0(n+1) + Kokkos::parallel_for( + "AmpereCurrentsMain", + range, + kernel::gr::CurrentsAmpere_kernel(domain.fields.em0, + domain.fields.cur0, + domain.mesh.metric, + coeff, + ni2, + domain.mesh.flds_bc())); + } else { + raise::Error("Wrong option for `g`", HERE); + } + } + + void TimeAverageDB(domain_t& domain) { + Kokkos::parallel_for("TimeAverageDB", + domain.mesh.rangeActiveCells(), + kernel::gr::TimeAverageDB_kernel(domain.fields.em, + domain.fields.em0, + domain.mesh.metric)); + } + + void TimeAverageJ(domain_t& domain) { + Kokkos::parallel_for("TimeAverageJ", + domain.mesh.rangeActiveCells(), + kernel::gr::TimeAverageJ_kernel(domain.fields.cur, + domain.fields.cur0, + domain.mesh.metric)); + } + + void CurrentsDeposit(domain_t& domain) { + auto scatter_cur0 = Kokkos::Experimental::create_scatter_view( + domain.fields.cur0); + 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; + } + Kokkos::parallel_for("CurrentsDeposit", + species.rangeActiveParticles(), + kernel::DepositCurrents_kernel( + scatter_cur0, + species.i1, + species.i2, + species.i3, + species.i1_prev, + species.i2_prev, + species.i3_prev, + species.dx1, + species.dx2, + species.dx3, + species.dx1_prev, + species.dx2_prev, + species.dx3_prev, + species.ux1, + species.ux2, + species.ux3, + species.phi, + species.weight, + species.tag, + domain.mesh.metric, + (real_t)(species.charge()), + dt)); + } + Kokkos::Experimental::contribute(domain.fields.cur0, scatter_cur0); + } + + void CurrentsFilter(domain_t& domain) { + logger::Checkpoint("Launching currents filtering kernels", HERE); + auto range = CreateRangePolicy( + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + const auto nfilter = m_params.template get( + "algorithms.current_filters"); + tuple_t size; + size[0] = domain.mesh.n_active(in::x1); + size[1] = domain.mesh.n_active(in::x2); + + // !TODO: this needs to be done more efficiently + for (unsigned short i = 0; i < nfilter; ++i) { + Kokkos::deep_copy(domain.fields.buff, domain.fields.cur0); + Kokkos::parallel_for("CurrentsFilter", + range, + kernel::DigitalFilter_kernel( + domain.fields.cur0, + domain.fields.buff, + size, + domain.mesh.flds_bc())); + m_metadomain.CommunicateFields(domain, Comm::J); // J0 + } + } + + 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 / 2) omegaB0 + const auto coeff = q_ovr_m * HALF * dt * + m_params.template get( + "algorithms.timestep.correction") * + m_params.template get("scales.omegaB0"); + const auto eps = m_params.template get( + "algorithms.gr.pusher_eps"); + const auto niter = m_params.template get( + "algorithms.gr.pusher_niter"); + // clang-format off + if (species.pusher() == PrtlPusher::PHOTON) { + auto range_policy = Kokkos::RangePolicy( + 0, + species.npart()); + + Kokkos::parallel_for( + "ParticlePusher", + range_policy, + kernel::gr::Pusher_kernel( + domain.fields.em, + domain.fields.em0, + species.i1, species.i2, species.i3, + species.i1_prev, species.i2_prev, species.i3_prev, + species.dx1, species.dx2, species.dx3, + species.dx1_prev, species.dx2_prev, species.dx3_prev, + species.ux1, species.ux2, species.ux3, + species.phi, species.tag, + domain.mesh.metric, + coeff, dt, + domain.mesh.n_active(in::x1), + domain.mesh.n_active(in::x2), + domain.mesh.n_active(in::x3), + eps, niter, + domain.mesh.prtl_bc() + )); + } else if (species.pusher() == PrtlPusher::BORIS) { + auto range_policy = Kokkos::RangePolicy( + 0, + species.npart()); + Kokkos::parallel_for( + "ParticlePusher", + range_policy, + kernel::gr::Pusher_kernel( + domain.fields.em, + domain.fields.em0, + species.i1, species.i2, species.i3, + species.i1_prev, species.i2_prev, species.i3_prev, + species.dx1, species.dx2, species.dx3, + species.dx1_prev, species.dx2_prev, species.dx3_prev, + species.ux1, species.ux2, species.ux3, + species.phi, species.tag, + domain.mesh.metric, + coeff, dt, + domain.mesh.n_active(in::x1), + domain.mesh.n_active(in::x2), + domain.mesh.n_active(in::x3), + eps, niter, + domain.mesh.prtl_bc() + )); + } else if (species.pusher() == PrtlPusher::NONE) { + // do nothing + } else { + raise::Error("not implemented", HERE); + } + // clang-format on + } + } + }; } // namespace ntt #endif // ENGINES_GRPIC_GRPIC_H diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp index a538f853..507316b5 100644 --- a/src/framework/domain/communications.cpp +++ b/src/framework/domain/communications.cpp @@ -206,13 +206,22 @@ namespace ntt { template void Metadomain::CommunicateFields(Domain& domain, CommTags tags) { - const auto comm_fields = (tags & Comm::E) || (tags & Comm::B) || - (tags & Comm::J) || (tags & Comm::D) || - (tags & Comm::D0) || (tags & Comm::B0); - const bool comm_em = (tags & Comm::E) || (tags & Comm::B) || (tags & Comm::D); - const bool comm_em0 = (tags & Comm::B0) || (tags & Comm::D0); + // const auto comm_fields = (tags & Comm::E) or (tags & Comm::B) or + // (tags & Comm::J) or (tags & Comm::D) or + // (tags & Comm::D0) or (tags & Comm::B0) or + // (tags & Comm::H); + const auto comm_em = ((S == SimEngine::SRPIC) and + ((tags & Comm::E) or (tags & Comm::B))) or + ((S == SimEngine::GRPIC) and + ((tags & Comm::D) or (tags & Comm::B))); + const bool comm_em0 = (S == SimEngine::GRPIC) and + ((tags & Comm::B0) or (tags & Comm::D0)); const bool comm_j = (tags & Comm::J); - raise::ErrorIf(not comm_fields, "CommunicateFields called with no task", HERE); + const bool comm_aux = (S == SimEngine::GRPIC) and + ((tags & Comm::E) or (tags & Comm::H)); + raise::ErrorIf(not(comm_em or comm_em0 or comm_j or comm_aux), + "CommunicateFields called with no task", + HERE); std::string comms = ""; if (tags & Comm::E) { @@ -227,6 +236,9 @@ namespace ntt { if (tags & Comm::D) { comms += "D "; } + if (tags & Comm::H) { + comms += "H "; + } if (tags & Comm::D0) { comms += "D0 "; } @@ -243,16 +255,17 @@ namespace ntt { auto comp_range_fld = range_tuple_t {}; auto comp_range_cur = range_tuple_t {}; if constexpr (S == SimEngine::GRPIC) { - if (((tags & Comm::D) && (tags & Comm::B)) || - ((tags & Comm::D0) && (tags & Comm::B0))) { + if (((tags & Comm::D) and (tags & Comm::B)) or + ((tags & Comm::D0) and (tags & Comm::B0)) or + ((tags & Comm::E) and (tags & Comm::H))) { comp_range_fld = range_tuple_t(em::dx1, em::bx3 + 1); - } else if ((tags & Comm::D) || (tags & Comm::D0)) { + } else if ((tags & Comm::D) or (tags & Comm::D0) or (tags & Comm::E)) { comp_range_fld = range_tuple_t(em::dx1, em::dx3 + 1); - } else if ((tags & Comm::B) || (tags & Comm::B0)) { + } else if ((tags & Comm::B) or (tags & Comm::B0) or (tags & Comm::H)) { comp_range_fld = range_tuple_t(em::bx1, em::bx3 + 1); } } else if constexpr (S == SimEngine::SRPIC) { - if ((tags & Comm::E) && (tags & Comm::B)) { + if ((tags & Comm::E) and (tags & Comm::B)) { comp_range_fld = range_tuple_t(em::ex1, em::bx3 + 1); } else if (tags & Comm::E) { comp_range_fld = range_tuple_t(em::ex1, em::ex3 + 1); @@ -290,6 +303,19 @@ namespace ntt { false); } if constexpr (S == SimEngine::GRPIC) { + if (comm_aux) { + comm::CommunicateField(domain.index(), + domain.fields.aux, + domain.fields.aux, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_fld, + false); + } if (comm_em0) { comm::CommunicateField(domain.index(), domain.fields.em0, @@ -302,20 +328,59 @@ namespace ntt { recv_slice, comp_range_fld, false); + // @HACK_GR_1.2.0 -- this has to be done carefully + // comm::CommunicateField(domain.index(), + // domain.fields.aux, + // domain.fields.aux, + // send_ind, + // recv_ind, + // send_rank, + // recv_rank, + // send_slice, + // recv_slice, + // comp_range_fld, + // false); + } + if (comm_j) { + comm::CommunicateField(domain.index(), + domain.fields.cur0, + domain.fields.cur0, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_cur, + false); + } + } else { + if (comm_em) { + comm::CommunicateField(domain.index(), + domain.fields.em, + domain.fields.em, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_fld, + false); + } + if (comm_j) { + comm::CommunicateField(domain.index(), + domain.fields.cur, + domain.fields.cur, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_cur, + false); } - } - if (comm_j) { - comm::CommunicateField(domain.index(), - domain.fields.cur, - domain.fields.cur, - send_ind, - recv_ind, - send_rank, - recv_rank, - send_slice, - recv_slice, - comp_range_cur, - false); } } } @@ -435,17 +500,31 @@ namespace ntt { continue; } if (comm_j) { - comm::CommunicateField(domain.index(), - domain.fields.cur, - domain.fields.buff, - send_ind, - recv_ind, - send_rank, - recv_rank, - send_slice, - recv_slice, - comp_range_cur, - synchronize); + if constexpr (S == SimEngine::GRPIC) { + comm::CommunicateField(domain.index(), + domain.fields.cur0, + domain.fields.buff, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_cur, + synchronize); + } else { + comm::CommunicateField(domain.index(), + domain.fields.cur, + domain.fields.buff, + send_ind, + recv_ind, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_range_cur, + synchronize); + } } if (comm_bckp) { comm::CommunicateField(domain.index(), @@ -475,10 +554,17 @@ namespace ntt { } } if (comm_j) { - AddBufferedFields(domain.fields.cur, - domain.fields.buff, - domain.mesh.rangeActiveCells(), - comp_range_cur); + if constexpr (S == SimEngine::GRPIC) { + AddBufferedFields(domain.fields.cur0, + domain.fields.buff, + domain.mesh.rangeActiveCells(), + comp_range_cur); + } else { + AddBufferedFields(domain.fields.cur, + domain.fields.buff, + domain.mesh.rangeActiveCells(), + comp_range_cur); + } } if (comm_bckp) { AddBufferedFields(domain.fields.bckp, diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp index a0326913..32275448 100644 --- a/src/framework/domain/metadomain.cpp +++ b/src/framework/domain/metadomain.cpp @@ -377,7 +377,7 @@ namespace ntt { dx_min_from_domains = std::min(dx_min_from_domains, current_dx_min); } raise::ErrorIf( - not cmp::AlmostEqual(dx_min / dx_min_from_domains, ONE, epsilon), + not cmp::AlmostEqual_host(dx_min / dx_min_from_domains, ONE, epsilon), "dx_min is not the same across all domains: " + std::to_string(dx_min) + " " + std::to_string(dx_min_from_domains), HERE); @@ -392,10 +392,7 @@ namespace ntt { mpi::get_type(), MPI_COMM_WORLD); for (const auto& dx : dx_mins) { - raise::ErrorIf(!cmp::AlmostEqual(dx, - dx_min, - std::numeric_limits::epsilon() * - static_cast(10.0)), + raise::ErrorIf(not cmp::AlmostEqual_host(dx / dx_min, ONE, epsilon), "dx_min is not the same across all MPI ranks", HERE); } diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 63f4b6e7..e458de0e 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -88,7 +88,13 @@ namespace ntt { const auto species_to_write = params.template get>( "output.particles.species"); g_writer.defineFieldOutputs(S, all_fields_to_write); - g_writer.defineParticleOutputs(M::PrtlDim, species_to_write); + + Dimension dim = M::PrtlDim; + if constexpr (M::CoordType != Coord::Cart) { + dim = Dim::_3D; + } + g_writer.defineParticleOutputs(dim, species_to_write); + // spectra write all particle species std::vector spectra_species {}; for (const auto& sp : species_params()) { @@ -181,6 +187,39 @@ namespace ntt { } } + template + void ComputeVectorPotential(ndfield_t& buffer, + ndfield_t& EM, + unsigned short buff_idx, + const Mesh mesh) { + if constexpr (M::Dim == Dim::_2D) { + const auto i2_min = mesh.i_min(in::x2); + // !TODO: this is quite slow + Kokkos::parallel_for( + "ComputeVectorPotential", + mesh.rangeActiveCells(), + Lambda(index_t i, index_t j) { + const real_t i_ { static_cast(static_cast(i) - (N_GHOSTS)) }; + const auto k_min = (i2_min - (N_GHOSTS)) + 1; + const auto k_max = (j - (N_GHOSTS)); + real_t A3 = ZERO; + for (auto k { k_min }; k <= k_max; ++k) { + real_t k_ = static_cast(k); + real_t sqrt_detH_ij1 { mesh.metric.sqrt_det_h({ i_, k_ - HALF }) }; + real_t sqrt_detH_ij2 { mesh.metric.sqrt_det_h({ i_, k_ + HALF }) }; + auto k1 { k + N_GHOSTS }; + A3 += HALF * (sqrt_detH_ij1 * EM(i, k1 - 1, em::bx1) + + sqrt_detH_ij2 * EM(i, k1, em::bx1)); + } + buffer(i, j, buff_idx) = A3; + }); + } else { + raise::KernelError( + HERE, + "ComputeVectorPotential: 2D implementation called for D != 2"); + } + } + template auto Metadomain::Write( const SimulationParams& params, @@ -387,6 +426,17 @@ namespace ntt { raise::Error("Custom output requested but no function provided", HERE); } + } else if (fld.is_vpotential()) { + if (S == SimEngine::GRPIC && M::Dim == Dim::_2D) { + const auto c = static_cast(addresses.back()); + ComputeVectorPotential(local_domain->fields.bckp, + local_domain->fields.em, + c, + local_domain->mesh); + } else { + raise::Error("Vector potential can only be computed for GRPIC in 2D", + HERE); + } } else { raise::Error("Wrong # of components requested for " "non-moment/non-custom output", diff --git a/src/global/global.h b/src/global/global.h index 04552db7..adffcf6e 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -288,9 +288,17 @@ namespace BC { Dx1 = 1 << 0, Dx2 = 1 << 1, Dx3 = 1 << 2, + Hx1 = 1 << 3, + Hx2 = 1 << 4, + Hx3 = 1 << 5, + Jx1 = 1 << 6, + Jx2 = 1 << 7, + Jx3 = 1 << 8, B = Bx1 | Bx2 | Bx3, E = Ex1 | Ex2 | Ex3, D = Dx1 | Dx2 | Dx3, + H = Hx1 | Hx2 | Hx3, + J = Jx1 | Jx2 | Jx3, }; } // namespace BC diff --git a/src/kernels/ampere_gr.hpp b/src/kernels/ampere_gr.hpp index 782dece0..56f72290 100644 --- a/src/kernels/ampere_gr.hpp +++ b/src/kernels/ampere_gr.hpp @@ -57,9 +57,6 @@ namespace kernel::gr { , 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); } @@ -77,6 +74,8 @@ namespace kernel::gr { if ((i2 == i2min) && is_axis_i2min) { // theta = 0 const real_t inv_polar_area_pH { ONE / metric.polar_area(i1_ + HALF) }; + const real_t inv_sqrt_detH_0pH { ONE / + metric.sqrt_det_h({ 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) + diff --git a/src/kernels/aux_fields_gr.hpp b/src/kernels/aux_fields_gr.hpp index 5744c309..0983b6b7 100644 --- a/src/kernels/aux_fields_gr.hpp +++ b/src/kernels/aux_fields_gr.hpp @@ -230,6 +230,77 @@ namespace kernel::gr { } } }; + + /** + * @brief Kernel for computing time average of B and D + * @tparam M Metric + */ + template + class TimeAverageDB_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static constexpr auto D = M::Dim; + + const ndfield_t BDf; + ndfield_t BDf0; + const M metric; + + public: + TimeAverageDB_kernel(const ndfield_t& BDf, + const ndfield_t& BDf0, + const M& metric) + : BDf { BDf } + , BDf0 { BDf0 } + , metric { metric } {} + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + BDf0(i1, i2, em::bx1) = HALF * (BDf0(i1, i2, em::bx1) + BDf(i1, i2, em::bx1)); + BDf0(i1, i2, em::bx2) = HALF * (BDf0(i1, i2, em::bx2) + BDf(i1, i2, em::bx2)); + BDf0(i1, i2, em::bx3) = HALF * (BDf0(i1, i2, em::bx3) + BDf(i1, i2, em::bx3)); + BDf0(i1, i2, em::ex1) = HALF * (BDf0(i1, i2, em::ex1) + BDf(i1, i2, em::ex1)); + BDf0(i1, i2, em::ex2) = HALF * (BDf0(i1, i2, em::ex2) + BDf(i1, i2, em::ex2)); + BDf0(i1, i2, em::ex3) = HALF * (BDf0(i1, i2, em::ex3) + BDf(i1, i2, em::ex3)); + } else { + raise::KernelError( + HERE, + "ComputeAuxH_kernel: 2D implementation called for D != 2"); + } + } + }; + + /** + * @brief Kernel for computing time average of J + * @tparam M Metric + */ + template + class TimeAverageJ_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static constexpr auto D = M::Dim; + + ndfield_t Jf; + const ndfield_t Jf0; + const M metric; + + public: + TimeAverageJ_kernel(const ndfield_t& Jf, + const ndfield_t& Jf0, + const M& metric) + : Jf { Jf } + , Jf0 { Jf0 } + , metric { metric } {} + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + Jf(i1, i2, cur::jx1) = HALF * (Jf0(i1, i2, cur::jx1) + Jf(i1, i2, cur::jx1)); + Jf(i1, i2, cur::jx2) = HALF * (Jf0(i1, i2, cur::jx2) + Jf(i1, i2, cur::jx2)); + Jf(i1, i2, cur::jx3) = HALF * (Jf0(i1, i2, cur::jx3) + Jf(i1, i2, cur::jx3)); + } else { + raise::KernelError( + HERE, + "ComputeAuxH_kernel: 2D implementation called for D != 2"); + } + } + }; } // namespace kernel::gr #endif // KERNELS_AUX_FIELDS_GR_HPP diff --git a/src/kernels/faraday_gr.hpp b/src/kernels/faraday_gr.hpp index 7fb327e6..b79afa46 100644 --- a/src/kernels/faraday_gr.hpp +++ b/src/kernels/faraday_gr.hpp @@ -73,7 +73,9 @@ namespace kernel::gr { Bout(i1, i2, em::bx1) = Bin(i1, i2, em::bx1) + coeff * inv_sqrt_detH_0pH * (E(i1, i2, em::ex3) - E(i1, i2 + 1, em::ex3)); - if ((i2 != i2min) || !is_axis_i2min) { + if ((i2 == i2min) && is_axis_i2min) { + Bout(i1, i2, em::bx2) = ZERO; + } else { const real_t inv_sqrt_detH_pH0 { ONE / metric.sqrt_det_h( { i1_ + HALF, i2_ }) }; Bout(i1, i2, em::bx2) = Bin(i1, i2, em::bx2) + diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp index a520616b..40dd3d13 100644 --- a/src/kernels/fields_bcs.hpp +++ b/src/kernels/fields_bcs.hpp @@ -4,7 +4,10 @@ * @implements * - kernel::bc::MatchBoundaries_kernel<> * - kernel::bc::AxisBoundaries_kernel<> + * - kernel::bc::AxisBoundariesGR_kernel<> + * - kernel::bc::AbsorbCurrentsGR_kernel<> * - kernel::bc::EnforcedBoundaries_kernel<> + * - kernel::bc::HorizonBoundaries_kernel<> * - kernel::bc::ConductorBoundaries_kernel<> * @namespaces: * - kernel::bc:: @@ -13,6 +16,7 @@ #ifndef KERNELS_FIELDS_BCS_HPP #define KERNELS_FIELDS_BCS_HPP +#include "enums.h" #include "global.h" #include "arch/kokkos_aliases.h" @@ -89,80 +93,64 @@ namespace kernel::bc { metric.template convert({ i1_ }, x_Ph_0); metric.template convert({ i1_ + HALF }, x_Ph_H); - // SRPIC - auto ex1_U { ZERO }, ex2_U { ZERO }, ex3_U { ZERO }, bx1_U { ZERO }, - bx2_U { ZERO }, bx3_U { ZERO }; - if (tags & BC::E) { - if constexpr (defines_ex1) { - ex1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF }, - fset.ex1(x_Ph_H)); - } - if constexpr (defines_ex2) { - ex2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_ }, - fset.ex2(x_Ph_0)); - } - if constexpr (defines_ex3) { - ex3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_ }, - fset.ex3(x_Ph_0)); - } - } - if (tags & BC::B) { - if constexpr (defines_bx1) { - bx1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_ }, - fset.bx1(x_Ph_0)); - } - if constexpr (defines_bx2) { - bx2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF }, - fset.bx2(x_Ph_H)); - } - if constexpr (defines_bx3) { - bx3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF }, - fset.bx3(x_Ph_H)); - } - } - if constexpr (defines_ex1 or defines_bx2 or defines_bx3) { - const auto dx = math::abs( - metric.template convert(i1_ + HALF) - xg_edge); - const auto s = shape(dx); + const auto s = shape(math::abs( + metric.template convert(i1_ + HALF) - xg_edge)); if constexpr (defines_ex1) { if (tags & BC::E) { - Fld(i1, em::ex1) = s * Fld(i1, em::ex1) + (ONE - s) * ex1_U; + Fld(i1, em::ex1) = s * Fld(i1, em::ex1) + + (ONE - s) * + metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.ex1(x_Ph_H)); } } if constexpr (defines_bx2 or defines_bx3) { if (tags & BC::B) { if constexpr (defines_bx2) { - Fld(i1, em::bx2) = s * Fld(i1, em::bx2) + (ONE - s) * bx2_U; + Fld(i1, em::bx2) = s * Fld(i1, em::bx2) + + (ONE - s) * + metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.bx2(x_Ph_H)); } if constexpr (defines_bx3) { - Fld(i1, em::bx3) = s * Fld(i1, em::bx3) + (ONE - s) * bx3_U; + Fld(i1, em::bx3) = s * Fld(i1, em::bx3) + + (ONE - s) * + metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.bx3(x_Ph_H)); } } } } if constexpr (defines_bx1 or defines_ex2 or defines_ex3) { - const auto dx = math::abs( - metric.template convert(i1_) - xg_edge); - const auto s = shape(dx); + const auto s = shape(math::abs( + metric.template convert(i1_) - xg_edge)); if constexpr (defines_bx1) { if (tags & BC::B) { - Fld(i1, em::bx1) = s * Fld(i1, em::bx1) + (ONE - s) * bx1_U; + Fld(i1, em::bx1) = s * Fld(i1, em::bx1) + + (ONE - s) * + metric.template transform<1, Idx::T, Idx::U>( + { i1_ }, + fset.bx1(x_Ph_0)); } } if constexpr (defines_ex2 or defines_ex3) { if (tags & BC::E) { if constexpr (defines_ex2) { - Fld(i1, em::ex2) = s * Fld(i1, em::ex2) + (ONE - s) * ex2_U; + Fld(i1, em::ex2) = s * Fld(i1, em::ex2) + + (ONE - s) * + metric.template transform<2, Idx::T, Idx::U>( + { i1_ }, + fset.ex2(x_Ph_0)); } if constexpr (defines_ex3) { - Fld(i1, em::ex3) = s * Fld(i1, em::ex3) + (ONE - s) * ex3_U; + Fld(i1, em::ex3) = s * Fld(i1, em::ex3) + + (ONE - s) * + metric.template transform<3, Idx::T, Idx::U>( + { i1_ }, + fset.ex3(x_Ph_0)); } } } @@ -183,123 +171,156 @@ namespace kernel::bc { const auto i1_ = COORD(i1); const auto i2_ = COORD(i2); - if constexpr (S == SimEngine::SRPIC) { - // SRPIC - if constexpr (defines_ex1 or defines_bx2) { - coord_t x_Ph_H0 { ZERO }; - metric.template convert({ i1_ + HALF, i2_ }, x_Ph_H0); - // i1 + 1/2, i2 - real_t xi_Cd; - if constexpr (o == in::x1) { - xi_Cd = i1_ + HALF; - } else { - xi_Cd = i2_; + // SRPIC + if constexpr (defines_ex1 or defines_dx1 or defines_bx2) { + // i1 + 1/2, i2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_ + HALF; + } else { + xi_Cd = i2_; + } + + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + + coord_t x_Ph_H0 { ZERO }; + metric.template convert({ i1_ + HALF, i2_ }, x_Ph_H0); + + if constexpr (defines_ex1 or defines_dx1) { + if ((tags & BC::E) or (tags & BC::D)) { + if constexpr (defines_ex1 and S == SimEngine::SRPIC) { + Fld(i1, i2, em::ex1) = s * Fld(i1, i2, em::ex1) + + (ONE - s) * + metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.ex1(x_Ph_H0)); + } else if constexpr (defines_dx1 and S == SimEngine::GRPIC) { + Fld(i1, i2, em::dx1) = s * Fld(i1, i2, em::dx1) + + (ONE - s) * fset.dx1(x_Ph_H0); + } + } + } + + if constexpr (defines_bx2) { + if (tags & BC::B) { + if constexpr (S == SimEngine::SRPIC) { + Fld(i1, i2, em::bx2) = s * Fld(i1, i2, em::bx2) + + (ONE - s) * + metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.bx2(x_Ph_H0)); + } else if constexpr (S == SimEngine::GRPIC) { + Fld(i1, i2, em::bx2) = s * Fld(i1, i2, em::bx2) + + (ONE - s) * fset.bx2(x_Ph_H0); + } } + } + } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); + if constexpr (defines_ex2 or defines_dx2 or defines_bx1) { + // i1, i2 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_; + } else { + xi_Cd = i2_ + HALF; + } - if constexpr (defines_ex1) { - if (tags & BC::E) { - const auto ex1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF, i2_ }, - fset.ex1(x_Ph_H0)); - Fld(i1, i2, em::ex1) = s * Fld(i1, i2, em::ex1) + (ONE - s) * ex1_U; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + + coord_t x_Ph_0H { ZERO }; + metric.template convert({ i1_, i2_ + HALF }, x_Ph_0H); + + if constexpr (defines_ex2 or defines_dx2) { + if ((tags & BC::E) or (tags & BC::D)) { + if constexpr (defines_ex2 and S == SimEngine::SRPIC) { + Fld(i1, i2, em::ex2) = s * Fld(i1, i2, em::ex2) + + (ONE - s) * + metric.template transform<2, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.ex2(x_Ph_0H)); + } else if constexpr (defines_dx2 and S == SimEngine::GRPIC) { + Fld(i1, i2, em::dx2) = s * Fld(i1, i2, em::dx2) + + (ONE - s) * fset.dx2(x_Ph_0H); } } - if constexpr (defines_bx2) { - if (tags & BC::B) { - const auto bx2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF, i2_ }, - fset.bx2(x_Ph_H0)); - Fld(i1, i2, em::bx2) = s * Fld(i1, i2, em::bx2) + (ONE - s) * bx2_U; + } + + if constexpr (defines_bx1) { + if (tags & BC::B) { + if constexpr (S == SimEngine::SRPIC) { + Fld(i1, i2, em::bx1) = s * Fld(i1, i2, em::bx1) + + (ONE - s) * + metric.template transform<1, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.bx1(x_Ph_0H)); + } else if constexpr (S == SimEngine::GRPIC) { + Fld(i1, i2, em::bx1) = s * Fld(i1, i2, em::bx1) + + (ONE - s) * fset.bx1(x_Ph_0H); } } } + } - if constexpr (defines_ex2 or defines_bx1) { - coord_t x_Ph_0H { ZERO }; - metric.template convert({ i1_, i2_ + HALF }, x_Ph_0H); - // i1, i2 + 1/2 + if constexpr (defines_ex3 or defines_dx3) { + if (tags & BC::E) { + // i1, i2 real_t xi_Cd; if constexpr (o == in::x1) { xi_Cd = i1_; } else { - xi_Cd = i2_ + HALF; + xi_Cd = i2_; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - if constexpr (defines_ex2) { - if (tags & BC::E) { - auto ex2_U { ZERO }; - ex2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_, i2_ + HALF }, - fset.ex2(x_Ph_0H)); - Fld(i1, i2, em::ex2) = s * Fld(i1, i2, em::ex2) + (ONE - s) * ex2_U; - } - } - if constexpr (defines_bx1) { - if (tags & BC::B) { - auto bx1_U { ZERO }; - bx1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_, i2_ + HALF }, - fset.bx1(x_Ph_0H)); - Fld(i1, i2, em::bx1) = s * Fld(i1, i2, em::bx1) + (ONE - s) * bx1_U; - } + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + + coord_t x_Ph_00 { ZERO }; + metric.template convert({ i1_, i2_ }, x_Ph_00); + + if constexpr (defines_ex3 and S == SimEngine::SRPIC) { + Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3) + + (ONE - s) * + metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_ }, + fset.ex3(x_Ph_00)); + } else if constexpr (defines_dx3 and S == SimEngine::GRPIC) { + Fld(i1, i2, em::dx3) = s * Fld(i1, i2, em::dx3) + + (ONE - s) * fset.dx3(x_Ph_00); } } + } - if constexpr (defines_ex3) { - if (tags & BC::E) { - auto ex3_U { ZERO }; - coord_t x_Ph_00 { ZERO }; - metric.template convert({ i1_, i2_ }, x_Ph_00); - ex3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_ }, - fset.ex3(x_Ph_00)); - // i1, i2 - real_t xi_Cd; - if constexpr (o == in::x1) { - xi_Cd = i1_; - } else { - xi_Cd = i2_; - } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3) + (ONE - s) * ex3_U; + if constexpr (defines_bx3) { + if (tags & BC::B) { + // i1 + 1/2, i2 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_ + HALF; + } else { + xi_Cd = i2_ + HALF; } - } - if constexpr (defines_bx3) { - if (tags & BC::B) { - auto bx3_U { ZERO }; - coord_t x_Ph_HH { ZERO }; - metric.template convert({ i1_ + HALF, i2_ + HALF }, - x_Ph_HH); - bx3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF, i2_ + HALF }, - fset.bx3(x_Ph_HH)); - // i1 + 1/2, i2 + 1/2 - real_t xi_Cd; - if constexpr (o == in::x1) { - xi_Cd = i1_ + HALF; - } else { - xi_Cd = i2_ + HALF; - } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - // bx3 - Fld(i1, i2, em::bx3) = s * Fld(i1, i2, em::bx3) + (ONE - s) * bx3_U; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + + coord_t x_Ph_HH { ZERO }; + metric.template convert({ i1_ + HALF, i2_ + HALF }, + x_Ph_HH); + + if constexpr (S == SimEngine::SRPIC) { + Fld(i1, i2, em::bx3) = s * Fld(i1, i2, em::bx3) + + (ONE - s) * + metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF, i2_ + HALF }, + fset.bx3(x_Ph_HH)); + } else if constexpr (S == SimEngine::GRPIC) { + Fld(i1, i2, em::bx3) = s * Fld(i1, i2, em::bx3) + + (ONE - s) * fset.bx3(x_Ph_HH); } } - } else { - // GRPIC - raise::KernelError(HERE, "GRPIC not implemented"); } } else { raise::KernelError( @@ -328,18 +349,18 @@ namespace kernel::bc { } else { xi_Cd = i3_; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto ex1_U { ZERO }; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + coord_t x_Ph_H00 { ZERO }; metric.template convert({ i1_ + HALF, i2_, i3_ }, x_Ph_H00); - ex1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF, i2_, i3_ }, - fset.ex1(x_Ph_H00)); - Fld(i1, i2, i3, em::ex1) = s * Fld(i1, i2, i3, em::ex1) + - (ONE - s) * ex1_U; + + Fld(i1, i2, i3, em::ex1) = + s * Fld(i1, i2, i3, em::ex1) + + (ONE - s) * metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF, i2_, i3_ }, + fset.ex1(x_Ph_H00)); } if constexpr (defines_ex2) { @@ -352,18 +373,18 @@ namespace kernel::bc { } else { xi_Cd = i3_; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto ex2_U { ZERO }; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + coord_t x_Ph_0H0 { ZERO }; metric.template convert({ i1_, i2_ + HALF, i3_ }, x_Ph_0H0); - ex2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_, i2_ + HALF, i3_ }, - fset.ex2(x_Ph_0H0)); - Fld(i1, i2, i3, em::ex2) = s * Fld(i1, i2, i3, em::ex2) + - (ONE - s) * ex2_U; + + Fld(i1, i2, i3, em::ex2) = + s * Fld(i1, i2, i3, em::ex2) + + (ONE - s) * metric.template transform<2, Idx::T, Idx::U>( + { i1_, i2_ + HALF, i3_ }, + fset.ex2(x_Ph_0H0)); } if constexpr (defines_ex3) { @@ -376,18 +397,17 @@ namespace kernel::bc { } else { xi_Cd = i3_ + HALF; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto ex3_U { ZERO }; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + coord_t x_Ph_00H { ZERO }; metric.template convert({ i1_, i2_, i3_ + HALF }, x_Ph_00H); - ex3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_, i3_ + HALF }, - fset.ex3(x_Ph_00H)); - Fld(i1, i2, i3, em::ex3) = s * Fld(i1, i2, i3, em::ex3) + - (ONE - s) * ex3_U; + Fld(i1, i2, i3, em::ex3) = + s * Fld(i1, i2, i3, em::ex3) + + (ONE - s) * metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_, i3_ + HALF }, + fset.ex3(x_Ph_00H)); } } } @@ -404,22 +424,19 @@ namespace kernel::bc { } else { xi_Cd = i3_ + HALF; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto bx1_U { ZERO }; - if constexpr (defines_bx1) { - coord_t x_Ph_0HH { ZERO }; - metric.template convert( - { i1_, i2_ + HALF, i3_ + HALF }, - x_Ph_0HH); - bx1_U = metric.template transform<1, Idx::T, Idx::U>( - { i1_, i2_ + HALF, i3_ + HALF }, - fset.bx1(x_Ph_0HH)); - } - // bx1 - Fld(i1, i2, i3, em::bx1) = s * Fld(i1, i2, i3, em::bx1) + - (ONE - s) * bx1_U; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + + coord_t x_Ph_0HH { ZERO }; + metric.template convert( + { i1_, i2_ + HALF, i3_ + HALF }, + x_Ph_0HH); + + Fld(i1, i2, i3, em::bx1) = + s * Fld(i1, i2, i3, em::bx1) + + (ONE - s) * metric.template transform<1, Idx::T, Idx::U>( + { i1_, i2_ + HALF, i3_ + HALF }, + fset.bx1(x_Ph_0HH)); } if constexpr (defines_bx2) { @@ -432,19 +449,19 @@ namespace kernel::bc { } else { xi_Cd = i3_ + HALF; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto bx2_U { ZERO }; + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + coord_t x_Ph_H0H { ZERO }; metric.template convert( { i1_ + HALF, i2_, i3_ + HALF }, x_Ph_H0H); - bx2_U = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF, i2_, i3_ + HALF }, - fset.bx2(x_Ph_H0H)); - Fld(i1, i2, i3, em::bx2) = s * Fld(i1, i2, i3, em::bx2) + - (ONE - s) * bx2_U; + + Fld(i1, i2, i3, em::bx2) = + s * Fld(i1, i2, i3, em::bx2) + + (ONE - s) * metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF, i2_, i3_ + HALF }, + fset.bx2(x_Ph_H0H)); } if constexpr (defines_bx3) { @@ -457,19 +474,20 @@ namespace kernel::bc { } else { xi_Cd = i3_; } - const auto dx = math::abs( - metric.template convert(xi_Cd) - xg_edge); - const auto s = shape(dx); - auto bx3_U { ZERO }; + + const auto s = shape(math::abs( + metric.template convert(xi_Cd) - xg_edge)); + coord_t x_Ph_HH0 { ZERO }; metric.template convert( { i1_ + HALF, i2_ + HALF, i3_ }, x_Ph_HH0); - bx3_U = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF, i2_ + HALF, i3_ }, - fset.bx3(x_Ph_HH0)); - Fld(i1, i2, i3, em::bx3) = s * Fld(i1, i2, i3, em::bx3) + - (ONE - s) * bx3_U; + + Fld(i1, i2, i3, em::bx3) = + s * Fld(i1, i2, i3, em::bx3) + + (ONE - s) * metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF, i2_ + HALF, i3_ }, + fset.bx3(x_Ph_HH0)); } } } @@ -827,14 +845,58 @@ namespace kernel::bc { } }; - /* - * @tparam I: Field Setter class - * @tparam M: Metric - * @tparam P: Positive/Negative direction - * @tparam O: Orientation - * - * @brief Applies enforced boundary conditions (fixed value) - */ + // /* + // * @tparam I: Field Setter class + // * @tparam M: Metric + // * @tparam P: Positive/Negative direction + // * @tparam O: Orientation + // * + // * @brief Applies enforced boundary conditions (fixed value) + // */ + // template + // struct AxisBoundariesGR_kernel { + // ndfield_t Fld; + // const std::size_t i_edge; + // const bool setE, setB; + // + // AxisBoundariesGR_kernel(ndfield_t Fld, std::size_t i_edge, BCTags tags) + // : Fld { Fld } // , i_edge { i_edge } + // , i_edge { P ? (i_edge + 1) : i_edge } + // , setE { tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3 } + // , setB { tags & BC::Bx1 or tags & BC::Bx2 or tags & BC::Bx3 } {} + // + // Inline void operator()(index_t i1) const { + // if constexpr (D == Dim::_2D) { + // // if (setB) { + // // Fld(i1, i_edge, em::bx2) = ZERO; + // // } + // if constexpr (not P) { + // if (setE) { + // Fld(i1, i_edge - 1, em::ex2) = -Fld(i1, i_edge, em::ex2); + // Fld(i1, i_edge, em::ex3) = ZERO; + // } + // if (setB) { + // Fld(i1, i_edge - 1, em::bx1) = Fld(i1, i_edge, em::bx1); + // Fld(i1, i_edge, em::bx2) = ZERO; + // Fld(i1, i_edge - 1, em::bx3) = Fld(i1, i_edge, em::bx3); + // } + // } else { + // if (setE) { + // Fld(i1, i_edge + 1, em::ex2) = -Fld(i1, i_edge, em::ex2); + // Fld(i1, i_edge + 1, em::ex3) = ZERO; + // } + // if (setB) { + // Fld(i1, i_edge + 1, em::bx1) = Fld(i1, i_edge, em::bx1); + // Fld(i1, i_edge + 1, em::bx2) = ZERO; + // Fld(i1, i_edge + 1, em::bx3) = Fld(i1, i_edge, em::bx3); + // } + // } + // } else { + // raise::KernelError(HERE, "AxisBoundariesGR_kernel: D != 2"); + // } + // } + // }; + template struct EnforcedBoundaries_kernel { static constexpr Dimension D = M::Dim; @@ -1154,6 +1216,105 @@ namespace kernel::bc { } }; + namespace gr { + + template + struct HorizonBoundaries_kernel { + ndfield_t Fld; + const std::size_t i1_min; + const bool setE, setB; + const std::size_t nfilter; + + HorizonBoundaries_kernel(ndfield_t Fld, + std::size_t i1_min, + BCTags tags, + std::size_t nfilter) + : Fld { Fld } + , i1_min { i1_min } + , setE { (tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3) or + (tags & BC::Dx1 or tags & BC::Dx2 or tags & BC::Dx3) } + , setB { (tags & BC::Bx1 or tags & BC::Bx2 or tags & BC::Bx3) or + (tags & BC::Hx1 or tags & BC::Hx2 or tags & BC::Hx3) } + , nfilter { nfilter } {} + + Inline void operator()(index_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + if (setE) { + for (unsigned short i = 0; i <= 2 + nfilter; ++i) { + Fld(i1_min - N_GHOSTS + i, + i2, + em::dx1) = Fld(i1_min + 1 + nfilter, i2, em::dx1); + Fld(i1_min - N_GHOSTS + i, + i2, + em::dx2) = Fld(i1_min + 1 + nfilter, i2, em::dx2); + Fld(i1_min - N_GHOSTS + i, + i2, + em::dx3) = Fld(i1_min + 1 + nfilter, i2, em::dx3); + } + } + if (setB) { + for (unsigned short i = 0; i <= 2 + nfilter; ++i) { + Fld(i1_min - N_GHOSTS + i, + i2, + em::bx1) = Fld(i1_min + 1 + nfilter, i2, em::bx1); + Fld(i1_min - N_GHOSTS + i, + i2, + em::bx2) = Fld(i1_min + 1 + nfilter, i2, em::bx2); + Fld(i1_min - N_GHOSTS + i, + i2, + em::bx3) = Fld(i1_min + 1 + nfilter, i2, em::bx3); + } + } + } else { + raise::KernelError( + HERE, + "HorizonBoundaries_kernel: 2D implementation called for D != 2"); + } + } + }; + + template + struct AbsorbCurrents_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static_assert(i <= static_cast(M::Dim), + "Invalid component index"); + + ndfield_t J; + const M metric; + const real_t xg_edge; + const real_t dx_abs; + + AbsorbCurrents_kernel(ndfield_t J, + const M& metric, + real_t xg_edge, + real_t dx_abs) + : J { J } + , metric { metric } + , xg_edge { xg_edge } + , dx_abs { dx_abs } {} + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + coord_t x_Cd { ZERO }; + x_Cd[0] = i1_; + x_Cd[1] = i2_; + const auto dx = math::abs( + metric.template convert(x_Cd[i - 1]) - xg_edge); + J(i1, i2, 0) *= math::tanh(dx / (INV_4 * dx_abs)); + J(i1, i2, 1) *= math::tanh(dx / (INV_4 * dx_abs)); + J(i1, i2, 2) *= math::tanh(dx / (INV_4 * dx_abs)); + + } else { + raise::KernelError( + HERE, + "gr::AbsorbCurrents_kernel: 2D implementation called for D != 2"); + } + } + }; + } // namespace gr + } // namespace kernel::bc #endif // KERNELS_FIELDS_BCS_HPP diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 8406fc4e..bb4235dc 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -120,25 +120,23 @@ namespace kernel { if constexpr (M::CoordType == Coord::Cart) { vec_t v_Ph { ZERO }; energy_dist(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Ph, v_Ph, v1); energy_dist(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Ph, v_Ph, v2); } else if constexpr (S == SimEngine::SRPIC) { - coord_t x_Ph_ { ZERO }; - x_Ph_[0] = x_Ph[0]; - x_Ph_[1] = x_Ph[1]; - x_Ph_[2] = ZERO; // phi = 0 + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + x_Cd_[2] = ZERO; // phi = 0 vec_t v_Ph { ZERO }; energy_dist(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Ph_, v_Ph, v1); + metric.template transform_xyz(x_Cd_, v_Ph, v1); energy_dist(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Ph_, v_Ph, v2); + metric.template transform_xyz(x_Cd_, v_Ph, v2); } else if constexpr (S == SimEngine::GRPIC) { vec_t v_Ph { ZERO }; energy_dist(x_Ph, v_Ph, spidx1); - metric.template transform(x_Ph, v_Ph, v1); + metric.template transform(x_Cd, v_Ph, v1); energy_dist(x_Ph, v_Ph, spidx2); - metric.template transform(x_Ph, v_Ph, v2); + metric.template transform(x_Cd, v_Ph, v2); } else { raise::KernelError(HERE, "Unknown simulation engine"); } @@ -291,25 +289,23 @@ namespace kernel { if constexpr (M::CoordType == Coord::Cart) { vec_t v_Ph { ZERO }; energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Ph, v_Ph, v1); energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Ph, v_Ph, v2); } else if constexpr (S == SimEngine::SRPIC) { - coord_t x_Ph_ { ZERO }; - x_Ph_[0] = x_Ph[0]; - x_Ph_[1] = x_Ph[1]; - x_Ph_[2] = ZERO; // phi = 0 + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + x_Cd_[2] = ZERO; // phi = 0 vec_t v_Ph { ZERO }; energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Ph_, v_Ph, v1); + metric.template transform_xyz(x_Cd_, v_Ph, v1); energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Ph_, v_Ph, v2); + metric.template transform_xyz(x_Cd_, v_Ph, v2); } else if constexpr (S == SimEngine::GRPIC) { vec_t v_Ph { ZERO }; energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform(x_Ph, v_Ph, v1); + metric.template transform(x_Cd, v_Ph, v1); energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform(x_Ph, v_Ph, v2); + metric.template transform(x_Cd, v_Ph, v2); } else { raise::KernelError(HERE, "Unknown simulation engine"); } @@ -508,7 +504,7 @@ namespace kernel { if constexpr (S == SimEngine::SRPIC) { global_metric.template transform_xyz(x_Cd_, u_Ph, u_Cd); } else if constexpr (S == SimEngine::GRPIC) { - global_metric.template transform(x_Cd, u_Ph, u_Cd); + global_metric.template transform(x_Cd, u_Ph, u_Cd); } else { raise::KernelError(HERE, "Unknown simulation engine"); } @@ -553,7 +549,7 @@ namespace kernel { if constexpr (S == SimEngine::SRPIC) { global_metric.template transform_xyz(x_Cd, u_Ph, u_Cd); } else if constexpr (S == SimEngine::GRPIC) { - global_metric.template transform(x_Cd, u_Ph, u_Cd); + global_metric.template transform(x_Cd, u_Ph, u_Cd); } else { raise::KernelError(HERE, "Unknown simulation engine"); } diff --git a/src/kernels/particle_pusher_gr.hpp b/src/kernels/particle_pusher_gr.hpp index 547463fa..2318d543 100644 --- a/src/kernels/particle_pusher_gr.hpp +++ b/src/kernels/particle_pusher_gr.hpp @@ -29,7 +29,7 @@ /* Local macros */ /* -------------------------------------------------------------------------- */ #define from_Xi_to_i(XI, I) \ - { I = static_cast((XI)); } + { I = static_cast((XI + 1)) - 1; } #define from_Xi_to_i_di(XI, I, DI) \ { \ @@ -65,6 +65,7 @@ namespace kernel::gr { static_assert(M::is_metric, "M must be a metric class"); static constexpr auto D = M::Dim; + private: const randacc_ndfield_t DB; const randacc_ndfield_t DB0; array_t i1, i2, i3; @@ -76,44 +77,43 @@ namespace kernel::gr { array_t tag; const M metric; - const real_t coeff, dt; - const int ni1, ni2, ni3; - const real_t epsilon; - const int niter; - const int i1_absorb; + const real_t coeff, dt; + const int ni1, ni2, ni3; + const real_t epsilon; + const unsigned short niter; bool is_axis_i2min { false }, is_axis_i2max { false }; bool is_absorb_i1min { false }, is_absorb_i1max { false }; public: - Pusher_kernel(const ndfield_t& DB, - const ndfield_t& DB0, - const array_t& i1, - const array_t& i2, - const array_t& i3, - const array_t& i1_prev, - const array_t& i2_prev, - const array_t& i3_prev, - const array_t& dx1, - const array_t& dx2, - const array_t& dx3, - const array_t& dx1_prev, - const array_t& dx2_prev, - const array_t& dx3_prev, - const array_t& ux1, - const array_t& ux2, - const array_t& ux3, - const array_t& phi, - const array_t& tag, - const M& metric, - const real_t& coeff, - const real_t& dt, - const int& ni1, - const int& ni2, - const int& ni3, - const real_t& epsilon, - const int& niter, - const boundaries_t& boundaries) + Pusher_kernel(const ndfield_t& DB, + const ndfield_t& DB0, + array_t& i1, + array_t& i2, + array_t& i3, + array_t& i1_prev, + array_t& i2_prev, + array_t& i3_prev, + array_t& dx1, + array_t& dx2, + array_t& dx3, + array_t& dx1_prev, + array_t& dx2_prev, + array_t& dx3_prev, + array_t& ux1, + array_t& ux2, + array_t& ux3, + array_t& phi, + array_t& tag, + const M& metric, + real_t coeff, + real_t dt, + int ni1, + int ni2, + int ni3, + const real_t& epsilon, + const unsigned short& niter, + const boundaries_t& boundaries) : DB { DB } , DB0 { DB0 } , i1 { i1 } @@ -140,10 +140,7 @@ namespace kernel::gr { , ni2 { ni2 } , ni3 { ni3 } , epsilon { epsilon } - , niter { niter } - , i1_absorb { static_cast(metric.template convert<1, Crd::Ph, Crd::Cd>( - metric.rhorizon())) - - 5 } { + , niter { niter } { raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE); is_absorb_i1min = (boundaries[0].first == PrtlBC::ABSORB) || @@ -333,31 +330,53 @@ namespace kernel::gr { // find contravariant midpoint velocity metric.template transform(xp, vp_mid, vp_mid_cntrv); - // find Gamma / alpha at midpoint + // find Gamma / alpha at midpointы real_t u0 { computeGamma(T {}, vp_mid, vp_mid_cntrv) / metric.alpha(xp) }; // find updated velocity + // vp_upd[0] = + // vp[0] + + // dt * + // (-metric.alpha(xp) * u0 * DERIVATIVE_IN_R(metric.alpha, xp) + + // vp_mid[0] * DERIVATIVE_IN_R(metric.beta1, xp) - + // (HALF / u0) * + // (DERIVATIVE_IN_R((metric.template h<1, 1>), xp) * SQR(vp_mid[0]) + + // DERIVATIVE_IN_R((metric.template h<2, 2>), xp) * SQR(vp_mid[1]) + + // DERIVATIVE_IN_R((metric.template h<3, 3>), xp) * SQR(vp_mid[2]) + + // TWO * DERIVATIVE_IN_R((metric.template h<1, 3>), xp) * + // vp_mid[0] * vp_mid[2])); + // vp_upd[1] = + // vp[1] + + // dt * + // (-metric.alpha(xp) * u0 * DERIVATIVE_IN_TH(metric.alpha, xp) + + // vp_mid[0] * DERIVATIVE_IN_TH(metric.beta1, xp) - + // (HALF / u0) * + // (DERIVATIVE_IN_TH((metric.template h<1, 1>), xp) * SQR(vp_mid[0]) + + // DERIVATIVE_IN_TH((metric.template h<2, 2>), xp) * SQR(vp_mid[1]) + + // DERIVATIVE_IN_TH((metric.template h<3, 3>), xp) * SQR(vp_mid[2]) + + // TWO * DERIVATIVE_IN_TH((metric.template h<1, 3>), xp) * + // vp_mid[0] * vp_mid[2])); vp_upd[0] = vp[0] + dt * - (-metric.alpha(xp) * u0 * DERIVATIVE_IN_R(metric.alpha, xp) + - vp_mid[0] * DERIVATIVE_IN_R(metric.beta1, xp) - + (-metric.alpha(xp) * u0 * metric.dr_alpha(xp) + + vp_mid[0] * metric.dr_beta1(xp) - (HALF / u0) * - (DERIVATIVE_IN_R((metric.template h<1, 1>), xp) * SQR(vp_mid[0]) + - DERIVATIVE_IN_R((metric.template h<2, 2>), xp) * SQR(vp_mid[1]) + - DERIVATIVE_IN_R((metric.template h<3, 3>), xp) * SQR(vp_mid[2]) + - TWO * DERIVATIVE_IN_R((metric.template h<1, 3>), xp) * + (metric.dr_h11(xp) * SQR(vp_mid[0]) + + metric.dr_h22(xp) * SQR(vp_mid[1]) + + metric.dr_h33(xp) * SQR(vp_mid[2]) + + TWO * metric.dr_h13(xp) * vp_mid[0] * vp_mid[2])); vp_upd[1] = vp[1] + dt * - (-metric.alpha(xp) * u0 * DERIVATIVE_IN_TH(alpha, xp) + - vp_mid[1] * DERIVATIVE_IN_TH(beta1, xp) - + (-metric.alpha(xp) * u0 * metric.dt_alpha(xp) + + vp_mid[0] * metric.dt_beta1(xp) - (HALF / u0) * - (DERIVATIVE_IN_TH((metric.template h<1, 1>), xp) * SQR(vp_mid[0]) + - DERIVATIVE_IN_TH((metric.template h<2, 2>), xp) * SQR(vp_mid[1]) + - DERIVATIVE_IN_TH((metric.template h<3, 3>), xp) * SQR(vp_mid[2]) + - TWO * DERIVATIVE_IN_TH((metric.template h<1, 3>), xp) * + (metric.dt_h11(xp) * SQR(vp_mid[0]) + + metric.dt_h22(xp) * SQR(vp_mid[1]) + + metric.dt_h33(xp) * SQR(vp_mid[2]) + + TWO * metric.dt_h13(xp) * vp_mid[0] * vp_mid[2])); } } else if constexpr (D == Dim::_3D) { @@ -455,7 +474,7 @@ namespace kernel::gr { dt * (-metric.alpha(xp_mid) * u0 * DERIVATIVE_IN_TH(metric.alpha, xp_mid) + - vp_mid[1] * DERIVATIVE_IN_TH(metric.beta1, xp_mid) - + vp_mid[0] * DERIVATIVE_IN_TH(metric.beta1, xp_mid) - (HALF / u0) * (DERIVATIVE_IN_TH((metric.template h<1, 1>), xp_mid) * SQR(vp_mid[0]) + @@ -691,7 +710,7 @@ namespace kernel::gr { { (xp[0] + xp_upd[0]) * HALF, (xp[1] + xp_upd[1]) * HALF }, vp_upd, phi(p)); - + // update coordinate int i1_, i2_; prtldx_t dx1_, dx2_; @@ -718,19 +737,23 @@ namespace kernel::gr { template Inline void Pusher_kernel::boundaryConditions(index_t& p) const { if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { - if (i1(p) < i1_absorb && is_absorb_i1min) { + if (i1(p) < 0 && is_absorb_i1min) { tag(p) = ParticleTag::dead; } else if (i1(p) >= ni1 && is_absorb_i1max) { tag(p) = ParticleTag::dead; } } if constexpr (D == Dim::_2D || D == Dim::_3D) { - if (i2(p) < 1) { + if (i2(p) < 0) { if (is_axis_i2min) { + i2(p) = 0; + dx2(p) = ONE - dx2(p); ux2(p) = -ux2(p); } - } else if (i2(p) >= ni2 - 1) { - if (is_axis_i2min) { + } else if (i2(p) >= ni2) { + if (is_axis_i2max) { + i2(p) = ni2 - 1; + dx2(p) = ONE - dx2(p); ux2(p) = -ux2(p); } } diff --git a/src/kernels/tests/deposit.cpp b/src/kernels/tests/deposit.cpp index 1570937b..c9acafae 100644 --- a/src/kernels/tests/deposit.cpp +++ b/src/kernels/tests/deposit.cpp @@ -27,13 +27,16 @@ void errorIf(bool condition, const std::string& message) { } } -inline static constexpr auto epsilon = std::numeric_limits::epsilon(); +const real_t eps = std::is_same_v ? (real_t)(1e-6) + : (real_t)(1e-3); -Inline auto equal(real_t a, real_t b, const char* msg = "", real_t acc = ONE) - -> bool { - const auto eps = epsilon * acc; - if (not cmp::AlmostEqual(a, b, eps)) { +Inline auto equal(real_t a, real_t b, const char* msg, real_t eps) -> bool { + if ((a - b) >= eps * math::max(math::fabs(a), math::fabs(b))) { printf("%.12e != %.12e %s\n", a, b, msg); + printf("%.12e >= %.12e %s\n", + a - b, + eps * math::max(math::fabs(a), math::fabs(b)), + msg); return false; } return true; @@ -49,8 +52,8 @@ void put_value(array_t arr, T value, int i) { template void testDeposit(const std::vector& res, const boundaries_t& ext, - const std::map& params = {}, - const real_t acc = ONE) { + const std::map& params, + const real_t eps) { static_assert(M::Dim == 2); errorIf(res.size() != M::Dim, "res.size() != M::Dim"); using namespace ntt; @@ -158,13 +161,13 @@ void testDeposit(const std::vector& res, if (not cmp::AlmostZero(SumDivJ)) { throw std::logic_error("DepositCurrents_kernel::SumDivJ != 0"); } - errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + N_GHOSTS, cur::jx1), Jx1, "", acc), + errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + N_GHOSTS, cur::jx1), Jx1, "", eps), "DepositCurrents_kernel::Jx1 is incorrect"); - errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + 1 + N_GHOSTS, cur::jx1), Jx2, "", acc), + errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + 1 + N_GHOSTS, cur::jx1), Jx2, "", eps), "DepositCurrents_kernel::Jx2 is incorrect"); - errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + N_GHOSTS, cur::jx2), Jy1, "", acc), + errorIf(not equal(J_h(i0 + N_GHOSTS, j0 + N_GHOSTS, cur::jx2), Jy1, "", eps), "DepositCurrents_kernel::Jy1 is incorrect"); - errorIf(not equal(J_h(i0 + 1 + N_GHOSTS, j0 + N_GHOSTS, cur::jx2), Jy2, "", acc), + errorIf(not equal(J_h(i0 + 1 + N_GHOSTS, j0 + N_GHOSTS, cur::jx2), Jy2, "", eps), "DepositCurrents_kernel::Jy2 is incorrect"); } @@ -183,41 +186,18 @@ auto main(int argc, char* argv[]) -> int { { 0.0, 55.0 }, { 0.0, 55.0 } }; + const std::map params { + { "r0", 0.0 }, + { "h", 0.25 }, + { "a", 0.9 } + }; - testDeposit, SimEngine::SRPIC>(res, xy_extent, {}, 500); - - testDeposit, SimEngine::SRPIC>(res, r_extent, {}, 500); - - testDeposit, SimEngine::SRPIC>(res, - r_extent, - { - { "r0", 0.0 }, - { "h", 0.25 } - }, - 500); - - testDeposit, SimEngine::GRPIC>(res, - r_extent, - { - { "a", 0.9 } - }, - 500); - - testDeposit, SimEngine::GRPIC>(res, - r_extent, - { - { "r0", 0.0 }, - { "h", 0.25 }, - { "a", 0.9 } - }, - 500); - - testDeposit, SimEngine::GRPIC>(res, - r_extent, - { - { "a", 0.9 } - }, - 500); + testDeposit, SimEngine::SRPIC>(res, xy_extent, {}, eps); + testDeposit, SimEngine::SRPIC>(res, r_extent, {}, eps); + testDeposit, SimEngine::SRPIC>(res, r_extent, params, eps); + testDeposit, SimEngine::GRPIC>(res, r_extent, params, eps); + testDeposit, SimEngine::GRPIC>(res, r_extent, params, eps); + testDeposit, SimEngine::GRPIC>(res, r_extent, params, eps); } catch (std::exception& e) { std::cerr << e.what() << std::endl; diff --git a/src/kernels/tests/reduced_stats.cpp b/src/kernels/tests/reduced_stats.cpp index 499c24aa..ee036395 100644 --- a/src/kernels/tests/reduced_stats.cpp +++ b/src/kernels/tests/reduced_stats.cpp @@ -4,13 +4,11 @@ #include "global.h" #include "arch/kokkos_aliases.h" -#include "utils/comparators.h" #include "utils/error.h" #include "metrics/minkowski.h" #include -#include #include using namespace ntt; @@ -18,15 +16,17 @@ using namespace metric; template class Fill_kernel { - ndfield_t& arr; - real_t v; - unsigned short c; + ndfield_t arr; + real_t v; + unsigned short c; public: Fill_kernel(ndfield_t& arr_, real_t v_, unsigned short c_) : arr { arr_ } , v { v_ } - , c { c_ } {} + , c { c_ } { + raise::ErrorIf(c_ >= N, "c > N", HERE); + } Inline void operator()(index_t i1) const { arr(i1, c) = v; diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 74f5f1e4..186f160a 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -17,6 +17,7 @@ #include "arch/kokkos_aliases.h" #include "utils/numeric.h" +#include "utils/comparators.h" #include "metrics/metric_base.h" @@ -225,6 +226,28 @@ namespace metric { return ONE / math::sqrt(ONE + z(r, theta)); } + /** + * dr derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dr_alpha(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + + return - (dr * Sigma(r, theta) - r * dr_Sigma) * CUBE(alpha(x)) / SQR(Sigma(r, theta)); + } + + /** + * dtheta derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dt_alpha(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return CUBE(alpha(x)) * r * dt_Sigma(theta) / SQR(Sigma(r, theta)); + } + /** * radial component of shift vector * @param x coordinate array in code units @@ -234,6 +257,146 @@ namespace metric { return dr_inv * z_ / (ONE + z_); } + /** + * dr derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dr_beta1(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + + return dr_inv * TWO * (dr * Sigma(r, theta) - r * dr_Sigma) / SQR(Sigma(r, theta) + TWO * r); + } + + /** + * dtheta derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dt_beta1(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - dr_inv * TWO * r * dt_Sigma(theta) / SQR(Sigma(r, theta) * (ONE + z(r, theta))); + } + + /** + * dr derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dr_h11(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + const real_t dr_Delta {TWO * dr * (r - ONE)}; + const real_t dr_A {FOUR * r * dr * (SQR(r) + SQR(a)) - SQR(a) * SQR(math::sin(theta)) * dr_Delta}; + + return (Sigma(r, theta) * (Sigma(r, theta) + TWO * r) * dr_A + - TWO * A(r, theta) * (r * dr_Sigma + Sigma(r, theta) * (dr_Sigma + dr))) + / (SQR(Sigma(r, theta) * (Sigma(r, theta) + TWO * r))) * SQR(dr_inv); + } + + /** + * dr derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dr_h22(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + + return - dr_Sigma / SQR(Sigma(r, theta)) * SQR(dtheta_inv); + } + + /** + * dr derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dr_h33(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + + return - dr_Sigma / SQR(Sigma(r, theta)) / SQR(math::sin(theta)); + } + + /** + * dr derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dr_h13(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + const real_t dr_Sigma {TWO * r * dr}; + + return - a * dr_Sigma / SQR(Sigma(r, theta)) * dr_inv; + } + + /** + * dtheta derivative of Sigma + * @param x coordinate array in code units + */ + Inline auto dt_Sigma(const real_t& theta) const -> real_t { + const real_t dt_Sigma {- TWO * SQR(a) * math::sin(theta) * math::cos(theta) * dtheta}; + if (cmp::AlmostZero(dt_Sigma)) + return ZERO; + else + return dt_Sigma; + } + + /** + * dtheta derivative of A + * @param x coordinate array in code units + */ + Inline auto dt_A(const real_t& r, const real_t& theta) const -> real_t { + const real_t dt_A {- TWO * SQR(a) * math::sin(theta) * math::cos(theta) * Delta(r) * dtheta}; + if (cmp::AlmostZero(dt_A)) + return ZERO; + else + return dt_A; + } + + /** + * dtheta derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dt_h11(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return (Sigma(r, theta) * (Sigma(r, theta) + TWO * r) * dt_A(r, theta) + - TWO * A(r, theta) * dt_Sigma(theta) * (r + Sigma(r, theta))) + / (SQR(Sigma(r, theta) * (Sigma(r, theta) + TWO * r))) * SQR(dr_inv); + } + + /** + * dtheta derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dt_h22(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - dt_Sigma(theta) / SQR(Sigma(r, theta)) * SQR(dtheta_inv); + } + + /** + * dtheta derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dt_h33(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - TWO * dtheta * math::cos(theta) * (Sigma(r, theta) - SQR(a) * SQR(math::sin(theta))) / CUBE(math::sin(theta)) / SQR(Sigma(r, theta)); + } + + /** + * dtheta derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dt_h13(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - a * dt_Sigma(theta) / SQR(Sigma(r, theta)) * dr_inv; + } + /** * sqrt(det(h_ij)) * @param x coordinate array in code units diff --git a/src/metrics/kerr_schild_0.h b/src/metrics/kerr_schild_0.h index e7cfc825..3621b3e5 100644 --- a/src/metrics/kerr_schild_0.h +++ b/src/metrics/kerr_schild_0.h @@ -36,6 +36,7 @@ namespace metric { private: const real_t dr, dtheta, dphi; const real_t dr_inv, dtheta_inv, dphi_inv; + const real_t a, rg_, rh_; public: static constexpr const char* Label { "kerr_schild_0" }; @@ -57,6 +58,9 @@ namespace metric { const boundaries_t& ext, const std::map& = {}) : MetricBase { res, ext } + , a { ZERO } + , rg_ { ONE } + , rh_ { TWO } , dr { (x1_max - x1_min) / nx1 } , dtheta { (x2_max - x2_min) / nx2 } , dphi { (x3_max - x3_min) / nx3 } @@ -70,17 +74,17 @@ namespace metric { [[nodiscard]] Inline auto spin() const -> const real_t& { - return ZERO; + return a; } [[nodiscard]] Inline auto rhorizon() const -> const real_t& { - return ZERO; + return rh_; } [[nodiscard]] Inline auto rg() const -> const real_t& { - return ZERO; + return rg_; } /** @@ -181,6 +185,22 @@ namespace metric { return ONE; } + /** + * dr derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dr_alpha(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dtheta derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dt_alpha(const coord_t& x) const -> real_t { + return ZERO; + } + /** * radial component of shift vector * @param x coordinate array in code units @@ -189,6 +209,92 @@ namespace metric { return ZERO; } + /** + * dr derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dr_beta1(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dtheta derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dt_beta1(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dr derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dr_h11(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dr derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dr_h22(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - TWO / CUBE(r) * SQR(dtheta_inv) * dr; + } + + /** + * dr derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dr_h33(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - TWO / CUBE(r) / SQR(math::sin(theta)) * dr; + } + + /** + * dr derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dr_h13(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dtheta derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dt_h11(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dtheta derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dt_h22(const coord_t& x) const -> real_t { + return ZERO; + } + + /** + * dtheta derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dt_h33(const coord_t& x) const -> real_t { + const real_t r {x[0] * dr + x1_min}; + const real_t theta {x[1] * dtheta + x2_min}; + return - TWO * math::cos(theta) / SQR(r) / CUBE(math::sin(theta)) * dtheta; + } + + /** + * dtheta derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dt_h13(const coord_t& x) const -> real_t { + return ZERO; + } + /** * sqrt(det(h_ij)) * @param x coordinate array in code units diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index d39ddd72..efe44611 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -243,6 +243,33 @@ namespace metric { return ONE / math::sqrt(ONE + z(r, theta)); } + /** + * dr derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dr_alpha(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + + return - (dx_r * Sigma(r, theta) - r * dr_Sigma) * CUBE(alpha(x)) / SQR(Sigma(r, theta)); + } + + /** + * dtheta derivative of lapse function + * @param x coordinate array in code units + */ + Inline auto dt_alpha(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + const real_t dx_dt {deta * (ONE + TWO * h0 * static_cast(constant::INV_PI_SQR) * (TWO * THREE * SQR(eta) - TWO * THREE * static_cast(constant::PI) * eta + static_cast(constant::PI_SQR))) }; + const real_t dt_Sigma {- TWO * SQR(a) * math::sin(theta) * math::cos(theta) * dx_dt}; + + return r * dt_Sigma * CUBE(alpha(x)) / SQR(Sigma(r, theta)); + } + /** * radial component of shift vector * @param x coordinate array in code units @@ -255,6 +282,167 @@ namespace metric { return math::exp(-chi) * dchi_inv * z_ / (ONE + z_); } + /** + * dr derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dr_beta1(const coord_t& x) const -> real_t { + const real_t chi { x[0] * dchi + chi_min }; + const real_t r { r0 + math::exp(chi) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + const real_t z_ { z(r, theta) }; + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + + return math::exp(-chi) * dchi_inv * TWO * (dx_r * Sigma(r, theta) - r * dr_Sigma) / SQR(Sigma(r, theta) + TWO * r) + - dchi * math::exp(-chi) * dchi_inv * z_ / (ONE + z_); + } + + /** + * dr derivative of radial component of shift vector + * @param x coordinate array in code units + */ + Inline auto dt_beta1(const coord_t& x) const -> real_t { + const real_t chi { x[0] * dchi + chi_min }; + const real_t r { r0 + math::exp(chi) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + return - math::exp(-chi) * dchi_inv * TWO * r * dt_Sigma(eta) / SQR(Sigma(r, theta) * (ONE + z(r, theta))); + } + + /** + * dr derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dr_h11(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + const real_t dr_Delta {TWO * dx_r * (r - ONE)}; + const real_t dr_A {FOUR * r * dx_r * (SQR(r) + SQR(a)) - SQR(a) * SQR(math::sin(theta)) * dr_Delta}; + + return (math::exp(-TWO * (x[0] * dchi + chi_min)) / SQR(dchi) + * (Sigma(r, theta) * (Sigma(r, theta) + TWO * r) * dr_A + - TWO * A(r, theta) * (r * dr_Sigma + Sigma(r, theta) * (dr_Sigma + dx_r))) + / (SQR(Sigma(r, theta) * (Sigma(r, theta) + TWO * r))) ) + -TWO * dchi * math::exp(-TWO * (x[0] * dchi + chi_min)) / SQR(dchi) * A(r, theta) / (Sigma(r, theta) * (Sigma(r, theta) + TWO * r)); + } + + /** + * dr derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dr_h22(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + + return - dr_Sigma / SQR(Sigma(r, theta)) / SQR(deta); + } + + /** + * dr derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dr_h33(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + + return - dr_Sigma / SQR(Sigma(r, theta)) / SQR(math::sin(theta)); + } + + /** + * dr derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dr_h13(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t theta { eta2theta(x[1] * deta + eta_min) }; + const real_t dx_r {dchi * math::exp(x[0] * dchi + chi_min)}; + const real_t dr_Sigma {TWO * r * dx_r}; + + return - a * dr_Sigma / SQR(Sigma(r, theta)) * (math::exp(-(x[0] * dchi + chi_min)) * dchi_inv) + - dchi * (math::exp(-(x[0] * dchi + chi_min)) * dchi_inv) * a / Sigma(r, theta); + } + + /** + * dtheta derivative of Sigma + * @param x coordinate array in code units + */ + Inline auto dt_Sigma(const real_t& eta) const -> real_t { + const real_t theta { eta2theta(eta) }; + const real_t dt_Sigma {- TWO * SQR(a) * math::sin(theta) * math::cos(theta) * dx_dt(eta)}; + if (cmp::AlmostZero(dt_Sigma)) + return ZERO; + else + return dt_Sigma; + } + + /** + * dtheta derivative of A + * @param x coordinate array in code units + */ + Inline auto dt_A(const real_t& r, const real_t& eta) const -> real_t { + const real_t theta { eta2theta(eta) }; + const real_t dt_A {- TWO * SQR(a) * math::sin(theta) * math::cos(theta) * Delta(r) * dx_dt(eta)}; + if (cmp::AlmostZero(dt_A)) + return ZERO; + else + return dt_A; + } + + /** + * dtheta derivative of h^11 + * @param x coordinate array in code units + */ + Inline auto dt_h11(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + return math::exp(-TWO * (x[0] * dchi + chi_min)) / SQR(dchi) + * (Sigma(r, theta) * (Sigma(r, theta) + TWO * r) * dt_A(r, eta) + - TWO * A(r, theta) * dt_Sigma(eta) * (r + Sigma(r, theta))) + / (SQR(Sigma(r, theta) * (Sigma(r, theta) + TWO * r))); + } + + /** + * dtheta derivative of h^22 + * @param x coordinate array in code units + */ + Inline auto dt_h22(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + return - dt_Sigma(eta) / SQR(Sigma(r, theta)) / SQR(deta); + } + + /** + * dtheta derivative of h^33 + * @param x coordinate array in code units + */ + Inline auto dt_h33(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + return - (dt_Sigma(eta) + TWO * math::cos(theta) / math::sin(theta) * Sigma(r, theta) * dx_dt(eta)) / SQR(Sigma(r, theta) * math::sin(theta)); + } + + /** + * dtheta derivative of h^13 + * @param x coordinate array in code units + */ + Inline auto dt_h13(const coord_t& x) const -> real_t { + const real_t r { r0 + math::exp(x[0] * dchi + chi_min) }; + const real_t eta {x[1] * deta + eta_min }; + const real_t theta { eta2theta(eta) }; + return - a * dt_Sigma(eta) / SQR(Sigma(r, theta)) * (math::exp(-(x[0] * dchi + chi_min)) * dchi_inv); + } + /** * sqrt(det(h_ij)) * @param x coordinate array in code units @@ -274,7 +462,7 @@ namespace metric { } /** - * sqrt(det(h_ij)) + * sqrt(det(h_ij)) divided by sin(theta). * @param x coordinate array in code units */ Inline auto sqrt_det_h_tilde(const coord_t& x) const -> real_t { @@ -301,7 +489,7 @@ namespace metric { math::sqrt( ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) / (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) * - (ONE - math::cos(eta2theta(HALF * deta + eta_min))); + (ONE - math::cos(eta2theta(HALF * deta))); } /** @@ -465,6 +653,19 @@ namespace metric { } } + /** + * @brief quasi-spherical eta to spherical theta + */ + Inline auto dx_dt(const real_t& eta) const -> real_t { + if (cmp::AlmostZero(h0)) { + return deta; + } else { + return deta * (ONE + + TWO * h0 * constant::INV_PI_SQR * + (TWO * THREE * SQR(eta) - TWO * THREE * constant::PI * eta + constant::PI_SQR)); + } + } + /** * @brief spherical theta to quasi-spherical eta */