diff --git a/input.example.toml b/input.example.toml index 4bf005b0e..33f094ca8 100644 --- a/input.example.toml +++ b/input.example.toml @@ -90,11 +90,11 @@ # Boundary conditions for fields: # @required # @type: 1/2/3-size array of string tuples, each of size 1 or 2 - # @valid: "PERIODIC", "ABSORB", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON" - # @example: [["CUSTOM", "ABSORB"]] (for 2D spherical [[rmin, rmax]]) + # @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON" + # @example: [["CUSTOM", "MATCH"]] (for 2D spherical [[rmin, rmax]]) # @note: When periodic in any of the directions, you should only set one value: [..., ["PERIODIC"], ...] - # @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "ABSORB"]] - # @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]] + # @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "match"]] + # @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["match"]] fields = "" # Boundary conditions for fields: # @required @@ -106,7 +106,7 @@ # @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]] particles = "" - [grid.boundaries.absorb] + [grid.boundaries.match] # Size of the absorption layer in physical (code) units: # @type: float # @default: 1% of the domain size (in shortest dimension) @@ -118,6 +118,14 @@ # @default: 1.0 coeff = "" + [grid.boundaries.absorb] + # Size of the absorption layer for particles in physical (code) units: + # @type: float + # @default: 1% of the domain size (in shortest dimension) + # @note: In spherical, this is the size of the layer in r from the outer wall + # @note: In cartesian, this is the same for all dimensions where applicable + ds = "" + [grid.boundaries.atmosphere] # @required: if ATMOSPHERE is one of the boundaries # Temperature of the atmosphere in units of m0 c^2 diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp index 55dffb5d9..b8f169521 100644 --- a/setups/srpic/shock/pgen.hpp +++ b/setups/srpic/shock/pgen.hpp @@ -13,6 +13,8 @@ #include "archetypes/problem_generator.h" #include "framework/domain/metadomain.h" +#include + namespace user { using namespace ntt; @@ -93,20 +95,21 @@ namespace user { inline PGen() {} - auto FixField(const em& comp) const -> real_t { + auto FixFieldsConst(const bc_in&, const em& comp) const + -> std::pair { if (comp == em::ex2) { - return init_flds.ex2({ ZERO }); + return { init_flds.ex2({ ZERO }), true }; } else if (comp == em::ex3) { - return init_flds.ex3({ ZERO }); - } else if (comp == em::bx1) { - return init_flds.bx1({ ZERO }); + return { init_flds.ex3({ ZERO }), true }; } else { - raise::Error("Other components should not be requested when BC is in X", - HERE); - return ZERO; + return { ZERO, false }; } } + auto MatchFields(real_t time) const -> InitFields { + return init_flds; + } + inline void InitPrtls(Domain& local_domain) { const auto energy_dist = arch::Maxwellian(local_domain.mesh.metric, local_domain.random_pool, diff --git a/setups/wip/magpump/pgen.hpp b/setups/wip/magpump/pgen.hpp new file mode 100644 index 000000000..21d4c8882 --- /dev/null +++ b/setups/wip/magpump/pgen.hpp @@ -0,0 +1,170 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.h" + +#include "archetypes/particle_injector.h" +#include "archetypes/problem_generator.h" +#include "framework/domain/metadomain.h" + +#include + +namespace user { + using namespace ntt; + + template + struct InitFields { + InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {} + + Inline auto bx1(const coord_t& x_Ph) const -> real_t { + return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + } + + Inline auto bx2(const coord_t& x_Ph) const -> real_t { + return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + } + + private: + const real_t Bsurf, Rstar; + }; + + template + struct DriveFields : public InitFields { + DriveFields(real_t time, real_t bsurf, real_t rstar) + : InitFields { bsurf, rstar } + , time { time } {} + + using InitFields::bx1; + using InitFields::bx2; + + Inline auto bx3(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex1(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex2(const coord_t& x_Ph) const -> real_t { + return ZERO; + } + + Inline auto ex3(const coord_t&) const -> real_t { + return ZERO; + } + + private: + const real_t time; + }; + + template + struct Inflow : public arch::EnergyDistribution { + Inflow(const M& metric, real_t vin) + : arch::EnergyDistribution { metric } + , vin { vin } {} + + Inline void operator()(const coord_t&, + vec_t& v_Ph, + unsigned short) const override { + v_Ph[0] = -vin; + } + + private: + const real_t vin; + }; + + template + struct Sphere : public arch::SpatialDistribution { + Sphere(const M& metric, real_t r0, real_t dr) + : arch::SpatialDistribution { metric } + , r0 { r0 } + , dr { dr } {} + + Inline auto operator()(const coord_t& x_Ph) const -> real_t override { + return math::exp(-SQR((x_Ph[0] - r0) / dr)) * + (x_Ph[1] > 0.25 && x_Ph[1] < constant::PI - 0.25); + } + + private: + const real_t r0, dr; + }; + + template + struct PGen : public arch::ProblemGenerator { + static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto metrics { + traits::compatible_with::value + }; + static constexpr auto dimensions { traits::compatible_with::value }; + + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + const real_t Bsurf, pump_period, pump_ampl, pump_radius, Rstar; + const real_t vin, drinj; + InitFields init_flds; + + inline PGen(const SimulationParams& p, const Metadomain& m) + : arch::ProblemGenerator { p } + , Bsurf { p.template get("setup.Bsurf", ONE) } + , pump_period { p.template get("setup.pump_period") } + , pump_ampl { p.template get("setup.pump_ampl") } + , pump_radius { p.template get("setup.pump_radius") } + , Rstar { m.mesh().extent(in::x1).first } + , vin { p.template get("setup.vin") } + , drinj { p.template get("setup.drinj") } + , init_flds { Bsurf, Rstar } {} + + auto FieldDriver(real_t time) const -> DriveFields { + return DriveFields { time, Bsurf, Rstar }; + } + + void CustomPostStep(std::size_t, long double time, Domain& domain) { + const real_t radius = pump_radius + + pump_ampl * + math::sin(time * constant::TWO_PI / pump_period); + const real_t dr = 1.0; + const auto& metric = domain.mesh.metric; + auto EM = domain.fields.em; + Kokkos::parallel_for( + "outerBC", + domain.mesh.rangeActiveCells(), + Lambda(index_t i1, index_t i2) { + const auto i1_ = COORD(i1), i2_ = COORD(i2); + const auto r = metric.template convert<1, Crd::Cd, Crd::Ph>(i1_); + if (r > radius - 5 * dr) { + const auto smooth = HALF * (ONE - math::tanh((r - radius) / dr)); + EM(i1, i2, em::ex1) = smooth * EM(i1, i2, em::ex1); + EM(i1, i2, em::ex2) = smooth * EM(i1, i2, em::ex2); + EM(i1, i2, em::ex3) = smooth * EM(i1, i2, em::ex3); + EM(i1, i2, em::bx1) = smooth * EM(i1, i2, em::bx1); + EM(i1, i2, em::bx2) = smooth * EM(i1, i2, em::bx2); + EM(i1, i2, em::bx3) = smooth * EM(i1, i2, em::bx3); + } + }); + + if (time < pump_period * 0.25) { + const auto energy_dist = Inflow(domain.mesh.metric, vin); + const auto spatial_dist = Sphere(domain.mesh.metric, radius, drinj); + const auto injector = arch::NonUniformInjector( + energy_dist, + spatial_dist, + { 1, 2 }); + + arch::InjectNonUniform>( + params, + domain, + injector, + ONE, + true); + } + } + }; + +} // namespace user + +#endif diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index 244a5f863..98396062a 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -579,8 +579,8 @@ namespace ntt { void FieldBoundaries(domain_t& domain, BCTags tags) { for (auto& direction : dir::Directions::orth) { - if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::ABSORB) { - AbsorbFieldsIn(direction, domain, tags); + if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::MATCH) { + MatchFieldsIn(direction, domain, tags); } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::AXIS) { if (domain.mesh.flds_bc_in(direction) == FldsBC::AXIS) { AxisFieldsIn(direction, domain, tags); @@ -601,14 +601,13 @@ namespace ntt { } // loop over directions } - void AbsorbFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { + void MatchFieldsIn(dir::direction_t direction, + domain_t& domain, + BCTags tags) { /** - * absorbing boundaries + * matching boundaries */ - const auto ds = m_params.template get( - "grid.boundaries.absorb.ds"); + const auto ds = 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(); @@ -647,40 +646,49 @@ namespace ntt { range_min[d] = intersect_range[d].first; range_max[d] = intersect_range[d].second; } - if (dim == in::x1) { - Kokkos::parallel_for( - "AbsorbFields", - CreateRangePolicy(range_min, range_max), - kernel::AbsorbBoundaries_kernel(domain.fields.em, - domain.mesh.metric, - xg_edge, - ds, - tags)); - } else if (dim == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - Kokkos::parallel_for( - "AbsorbFields", - CreateRangePolicy(range_min, range_max), - kernel::AbsorbBoundaries_kernel(domain.fields.em, - domain.mesh.metric, - xg_edge, - ds, - tags)); - } else { - raise::Error("Invalid dimension", HERE); - } - } else if (dim == in::x3) { - if constexpr (M::Dim == Dim::_3D) { + if constexpr (traits::has_member::value) { + auto match_fields = m_pgen.MatchFields(time); + if (dim == in::x1) { Kokkos::parallel_for( - "AbsorbFields", + "MatchFields", CreateRangePolicy(range_min, range_max), - kernel::AbsorbBoundaries_kernel(domain.fields.em, - domain.mesh.metric, - xg_edge, - ds, - tags)); - } else { - raise::Error("Invalid dimension", HERE); + kernel::MatchBoundaries_kernel( + domain.fields.em, + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags)); + } else if (dim == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "MatchFields", + CreateRangePolicy(range_min, range_max), + kernel::MatchBoundaries_kernel( + domain.fields.em, + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags)); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dim == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "MatchFields", + CreateRangePolicy(range_min, range_max), + kernel::MatchBoundaries_kernel( + domain.fields.em, + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags)); + } else { + raise::Error("Invalid dimension", HERE); + } } } } @@ -769,40 +777,51 @@ namespace ntt { if (tags & BC::B) { comps.push_back(normal_b_comp); } - if constexpr (traits::has_member::value) { - raise::Error("Field driver for fixed fields not implemented", HERE); - } else { - // if field driver not present, set fields to fixed values + if constexpr (traits::has_member::value) { + raise::Error("Non-const fixed fields not implemented", HERE); + } else if constexpr ( + traits::has_member::value) { for (const auto& comp : comps) { - auto value = ZERO; + auto value = ZERO; + bool shouldset = false; if constexpr ( - traits::has_member::value) { + traits::has_member::value) { // if fix field function present, read from it - value = m_pgen.FixField((em)comp); + const auto newset = m_pgen.FixFieldsConst( + (bc_in)(sign * ((short)dim + 1)), + (em)comp); + value = newset.first; + shouldset = newset.second; } - if constexpr (M::Dim == Dim::_1D) { - Kokkos::deep_copy(Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - comp), - value); - } else if constexpr (M::Dim == Dim::_2D) { - Kokkos::deep_copy(Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - std::make_pair(xi_min[1], xi_max[1]), - comp), - value); - } else if constexpr (M::Dim == Dim::_3D) { - Kokkos::deep_copy( - Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - std::make_pair(xi_min[1], xi_max[1]), - std::make_pair(xi_min[2], xi_max[2]), - comp), - value); - } else { - raise::Error("Invalid dimension", HERE); + if (shouldset) { + if constexpr (M::Dim == Dim::_1D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + comp), + value); + } else if constexpr (M::Dim == Dim::_2D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + std::make_pair(xi_min[1], xi_max[1]), + comp), + value); + } else if constexpr (M::Dim == Dim::_3D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + std::make_pair(xi_min[1], xi_max[1]), + std::make_pair(xi_min[2], xi_max[2]), + comp), + value); + } else { + raise::Error("Invalid dimension", HERE); + } } } + } else { + raise::Error("Fixed fields not present (both const and non-const)", HERE); } } @@ -812,7 +831,7 @@ namespace ntt { /** * atmosphere field boundaries */ - if constexpr (traits::has_member::value) { + if constexpr (traits::has_member::value) { const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction); const auto dd = static_cast(dim); boundaries_t box; @@ -841,7 +860,7 @@ namespace ntt { range_min[d] = intersect_range[d].first; range_max[d] = intersect_range[d].second; } - auto field_driver = m_pgen.FieldDriver(time); + auto atm_fields = m_pgen.AtmFields(time); std::size_t il_edge; if (sign > 0) { il_edge = range_min[dd] - N_GHOSTS; @@ -854,9 +873,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -864,9 +883,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -877,9 +896,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -887,9 +906,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -903,9 +922,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -913,9 +932,9 @@ namespace ntt { Kokkos::parallel_for( "AtmosphereBCFields", range, - kernel::EnforcedBoundaries_kernel( + kernel::EnforcedBoundaries_kernel( domain.fields.em, - field_driver, + atm_fields, domain.mesh.metric, il_edge, tags)); @@ -927,8 +946,7 @@ namespace ntt { raise::Error("Invalid dimension", HERE); } } else { - raise::Error("Field driver not implemented in PGEN for atmosphere BCs", - HERE); + raise::Error("Atm fields not implemented in PGEN for atmosphere BCs", HERE); } } diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index b667b5ac9..713a8a984 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -351,9 +351,9 @@ namespace ntt { auto atm_defined = false; for (const auto& bcs : flds_bc) { for (const auto& bc : bcs) { - if (fmt::toLower(bc) == "absorb") { - promiseToDefine("grid.boundaries.absorb.ds"); - promiseToDefine("grid.boundaries.absorb.coeff"); + if (fmt::toLower(bc) == "match") { + promiseToDefine("grid.boundaries.match.ds"); + promiseToDefine("grid.boundaries.match.coeff"); } if (fmt::toLower(bc) == "atmosphere") { raise::ErrorIf(atm_defined, @@ -386,7 +386,6 @@ namespace ntt { for (const auto& bc : bcs) { if (fmt::toLower(bc) == "absorb") { promiseToDefine("grid.boundaries.absorb.ds"); - promiseToDefine("grid.boundaries.absorb.coeff"); } if (fmt::toLower(bc) == "atmosphere") { raise::ErrorIf(atm_defined, @@ -730,6 +729,38 @@ namespace ntt { set("grid.boundaries.fields", flds_bc_pairwise); set("grid.boundaries.particles", prtl_bc_pairwise); + if (isPromised("grid.boundaries.match.ds")) { + if (coord_enum == Coord::Cart) { + auto min_extent = std::numeric_limits::max(); + for (const auto& e : extent_pairwise) { + min_extent = std::min(min_extent, e.second - e.first); + } + set("grid.boundaries.match.ds", + toml::find_or(toml_data, + "grid", + "boundaries", + "match", + "ds", + min_extent * defaults::bc::match::ds_frac)); + } else { + auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first; + set("grid.boundaries.match.ds", + toml::find_or(toml_data, + "grid", + "boundaries", + "match", + "ds", + r_extent * defaults::bc::match::ds_frac)); + } + set("grid.boundaries.match.coeff", + toml::find_or(toml_data, + "grid", + "boundaries", + "match", + "coeff", + defaults::bc::match::coeff)); + } + if (isPromised("grid.boundaries.absorb.ds")) { if (coord_enum == Coord::Cart) { auto min_extent = std::numeric_limits::max(); @@ -753,13 +784,6 @@ namespace ntt { "ds", r_extent * defaults::bc::absorb::ds_frac)); } - set("grid.boundaries.absorb.coeff", - toml::find_or(toml_data, - "grid", - "boundaries", - "absorb", - "coeff", - defaults::bc::absorb::coeff)); } if (isPromised("grid.boundaries.atmosphere.temperature")) { diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h index 9fd40e201..4cde4fca5 100644 --- a/src/global/arch/traits.h +++ b/src/global/arch/traits.h @@ -10,7 +10,11 @@ * - traits::run_t, traits::to_string_t * - traits::pgen::init_flds_t * - traits::pgen::ext_force_t - * - traits::pgen::field_driver_t + * - traits::pgen::atm_fields_t + * - traits::pgen::match_fields_const_t + * - traits::pgen::match_fields_t + * - traits::pgen::fix_fields_const_t + * - traits::pgen::fix_fields_t * - traits::pgen::init_prtls_t * - traits::pgen::custom_fields_t * - traits::pgen::custom_field_output_t @@ -94,10 +98,19 @@ namespace traits { using ext_force_t = decltype(&T::ext_force); template - using field_driver_t = decltype(&T::FieldDriver); + using atm_fields_t = decltype(&T::AtmFields); template - using fix_field_t = decltype(&T::FixField); + using match_fields_t = decltype(&T::MatchFields); + + template + using match_fields_const_t = decltype(&T::MatchFieldsConst); + + template + using fix_fields_t = decltype(&T::FixFields); + + template + using fix_fields_const_t = decltype(&T::FixFieldsConst); template using custom_fields_t = decltype(&T::CustomFields); diff --git a/src/global/defaults.h b/src/global/defaults.h index be92acbf9..7b2ee89ab 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -41,9 +41,13 @@ namespace ntt::defaults { } // namespace gr namespace bc { - namespace absorb { + namespace match { const real_t ds_frac = 0.01; const real_t coeff = 1.0; + } // namespace match + + namespace absorb { + const real_t ds_frac = 0.01; } // namespace absorb } // namespace bc diff --git a/src/global/enums.h b/src/global/enums.h index 283cb456d..8f2495c13 100644 --- a/src/global/enums.h +++ b/src/global/enums.h @@ -8,7 +8,7 @@ * - enum ntt::SimEngine // SRPIC, GRPIC * - enum ntt::PrtlBC // periodic, absorb, atmosphere, custom, * reflect, horizon, axis, sync - * - enum ntt::FldsBC // periodic, absorb, fixed, atmosphere, + * - enum ntt::FldsBC // periodic, match, fixed, atmosphere, * custom, horizon, axis, sync * - enum ntt::PrtlPusher // boris, vay, photon, none * - enum ntt::Cooling // synchrotron, none @@ -215,7 +215,7 @@ namespace ntt { enum type : uint8_t { INVALID = 0, PERIODIC = 1, - ABSORB = 2, + MATCH = 2, FIXED = 3, ATMOSPHERE = 4, CUSTOM = 5, @@ -226,9 +226,9 @@ namespace ntt { constexpr FldsBC(uint8_t c) : enums_hidden::BaseEnum { c } {} - static constexpr type variants[] = { PERIODIC, ABSORB, FIXED, ATMOSPHERE, + static constexpr type variants[] = { PERIODIC, MATCH, FIXED, ATMOSPHERE, CUSTOM, HORIZON, AXIS, SYNC }; - static constexpr const char* lookup[] = { "periodic", "absorb", "fixed", + static constexpr const char* lookup[] = { "periodic", "match", "fixed", "atmosphere", "custom", "horizon", "axis", "sync" }; static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]); diff --git a/src/global/global.h b/src/global/global.h index ad524fb0e..5fcbd9d06 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -13,6 +13,10 @@ * - enum PrepareOutput * - enum CellLayer // allLayer, activeLayer, minGhostLayer, * minActiveLayer, maxActiveLayer, maxGhostLayer + * - enum Idx // U, D, T, XYZ, Sph, PU, PD + * - enum Crd // Cd, Ph, XYZ, Sph + * - enum in // x1, x2, x3 + * - enum bc_in // Px1, Mx1, Px2, Mx2, Px3, Mx3 * - type box_region_t * - files::LogFile, files::ErrFile, files::InfoFile * - type prtldx_t @@ -184,6 +188,15 @@ enum class in : unsigned short { x3 = 2, }; +enum class bc_in : short { + Mx1 = -1, + Px1 = 1, + Mx2 = -2, + Px2 = 2, + Mx3 = -3, + Px3 = 3, +}; + template using box_region_t = CellLayer[D]; diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp index 4d678e85e..7785ec1a3 100644 --- a/src/global/tests/enums.cpp +++ b/src/global/tests/enums.cpp @@ -61,7 +61,7 @@ auto main() -> int { enum_str_t all_simulation_engines = { "srpic", "grpic" }; enum_str_t all_particle_bcs = { "periodic", "absorb", "atmosphere", "custom", "reflect", "horizon", "axis", "sync" }; - enum_str_t all_fields_bcs = { "periodic", "absorb", "fixed", "atmosphere", + enum_str_t all_fields_bcs = { "periodic", "match", "fixed", "atmosphere", "custom", "horizon", "axis", "sync" }; enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" }; enum_str_t all_coolings = { "synchrotron", "none" }; diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp index 2f2a458bb..4d52dd207 100644 --- a/src/kernels/fields_bcs.hpp +++ b/src/kernels/fields_bcs.hpp @@ -15,61 +15,133 @@ namespace kernel { using namespace ntt; - template - struct AbsorbBoundaries_kernel { + template + struct MatchBoundaries_kernel { static_assert(M::is_metric, "M must be a metric class"); - static_assert(i <= static_cast(M::Dim), + static_assert(static_cast(o) < + static_cast(M::Dim), "Invalid component index"); + static constexpr idx_t i = static_cast(o) + 1u; + static constexpr bool defines_dx1 = traits::has_method::value; + static constexpr bool defines_dx2 = traits::has_method::value; + static constexpr bool defines_dx3 = traits::has_method::value; + static constexpr bool defines_ex1 = traits::has_method::value; + static constexpr bool defines_ex2 = traits::has_method::value; + static constexpr bool defines_ex3 = traits::has_method::value; + static constexpr bool defines_bx1 = traits::has_method::value; + static constexpr bool defines_bx2 = traits::has_method::value; + static constexpr bool defines_bx3 = traits::has_method::value; ndfield_t Fld; + const I fset; const M metric; const real_t xg_edge; const real_t dx_abs; const BCTags tags; - AbsorbBoundaries_kernel(ndfield_t Fld, - const M& metric, - real_t xg_edge, - real_t dx_abs, - BCTags tags) + MatchBoundaries_kernel(ndfield_t Fld, + const I& fset, + const M& metric, + real_t xg_edge, + real_t dx_abs, + BCTags tags) : Fld { Fld } + , fset { fset } , metric { metric } , xg_edge { xg_edge } , dx_abs { dx_abs } , tags { tags } {} + Inline auto shape(const real_t& dx) const -> real_t { + return math::tanh(dx * FOUR / dx_abs); + } + Inline void operator()(index_t i1) const { if constexpr (M::Dim == Dim::_1D) { const auto i1_ = COORD(i1); - for (const auto comp : - { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) { - if ((comp == em::ex1) and not(tags & BC::Ex1)) { - continue; - } else if ((comp == em::ex2) and not(tags & BC::Ex2)) { - continue; - } else if ((comp == em::ex3) and not(tags & BC::Ex3)) { - continue; - } else if ((comp == em::bx1) and not(tags & BC::Bx1)) { - continue; - } else if ((comp == em::bx2) and not(tags & BC::Bx2)) { - continue; - } else if ((comp == em::bx3) and not(tags & BC::Bx3)) { - continue; - } - coord_t x_Cd { ZERO }; - if (comp == em::ex1 or comp == em::bx2 or comp == em::bx3) { - x_Cd[0] = i1_ + HALF; - } else if (comp == em::ex2 or comp == em::bx1 or comp == em::ex3) { - x_Cd[0] = i1_; - } - const auto dx = math::abs( - metric.template convert(x_Cd[i - 1]) - xg_edge); - Fld(i1, comp) *= math::tanh(dx / (INV_4 * dx_abs)); + + coord_t x_Ph_0 { ZERO }; + coord_t x_Ph_H { ZERO }; + metric.template convert({ i1_ }, x_Ph_0); + metric.template convert({ i1_ + HALF }, x_Ph_H); + + if constexpr (S == SimEngine::SRPIC) { + // 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)); + } + } + + { + const auto dx = math::abs( + metric.template convert(i1_ + HALF) - xg_edge); + const auto s = shape(dx); + if (tags & BC::E) { + // ex1 + Fld(i1, em::ex1) = s * Fld(i1, em::ex1) + (ONE - s) * ex1_U; + } + if (tags & BC::B) { + // bx2 + Fld(i1, em::bx2) = s * Fld(i1, em::bx2) + (ONE - s) * bx2_U; + // bx3 + Fld(i1, em::bx3) = s * Fld(i1, em::bx3) + (ONE - s) * bx3_U; + } + } + { + const auto dx = math::abs( + metric.template convert(i1_) - xg_edge); + const auto s = shape(dx); + if (tags & BC::B) { + // bx1 + Fld(i1, em::bx1) = s * Fld(i1, em::bx1) + (ONE - s) * bx1_U; + } + if (tags & BC::E) { + // ex2 + Fld(i1, em::ex2) = s * Fld(i1, em::ex2) + (ONE - s) * ex2_U; + // ex3 + Fld(i1, em::ex3) = s * Fld(i1, em::ex3) + (ONE - s) * ex3_U; + } + } + } else { + // GRPIC + raise::KernelError(HERE, "1D GRPIC not implemented"); } } else { raise::KernelError( HERE, - "AbsorbFields_kernel: 1D implementation called for D != 1"); + "MatchBoundaries_kernel: 1D implementation called for D != 1"); } } @@ -77,43 +149,130 @@ namespace kernel { if constexpr (M::Dim == Dim::_2D) { const auto i1_ = COORD(i1); const auto i2_ = COORD(i2); - for (const auto comp : - { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) { - if ((comp == em::ex1) and not(tags & BC::Ex1)) { - continue; - } else if ((comp == em::ex2) and not(tags & BC::Ex2)) { - continue; - } else if ((comp == em::ex3) and not(tags & BC::Ex3)) { - continue; - } else if ((comp == em::bx1) and not(tags & BC::Bx1)) { - continue; - } else if ((comp == em::bx2) and not(tags & BC::Bx2)) { - continue; - } else if ((comp == em::bx3) and not(tags & BC::Bx3)) { - continue; - } - coord_t x_Cd { ZERO }; - if (comp == em::ex1 or comp == em::bx2) { - x_Cd[0] = i1_ + HALF; - x_Cd[1] = i2_; - } else if (comp == em::ex2 or comp == em::bx1) { - x_Cd[0] = i1_; - x_Cd[1] = i2_ + HALF; - } else if (comp == em::ex3) { - x_Cd[0] = i1_; - x_Cd[1] = i2_; - } else if (comp == em::bx3) { - x_Cd[0] = i1_ + HALF; - x_Cd[1] = i2_ + HALF; - } - const auto dx = math::abs( - metric.template convert(x_Cd[i - 1]) - xg_edge); - Fld(i1, i2, comp) *= math::tanh(dx / (INV_4 * dx_abs)); + + if constexpr (S == SimEngine::SRPIC) { + // SRPIC + { + 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_; + } + const auto dx = math::abs( + metric.template convert(xi_Cd) - xg_edge); + const auto s = shape(dx); + if (tags & BC::E) { + auto ex1_U { ZERO }; + if constexpr (defines_ex1) { + ex1_U = metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.ex1(x_Ph_H0)); + } + // ex1 + Fld(i1, i2, em::ex1) = s * Fld(i1, i2, em::ex1) + (ONE - s) * ex1_U; + } + if (tags & BC::B) { + auto bx2_U { ZERO }; + if constexpr (defines_bx2) { + bx2_U = metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.bx2(x_Ph_H0)); + } + // bx2 + Fld(i1, i2, em::bx2) = s * Fld(i1, i2, em::bx2) + (ONE - s) * bx2_U; + } + } + { + coord_t x_Ph_0H { ZERO }; + metric.template convert({ i1_, i2_ + HALF }, x_Ph_0H); + // i1, i2 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_; + } else { + xi_Cd = i2_ + HALF; + } + const auto dx = math::abs( + metric.template convert(xi_Cd) - xg_edge); + const auto s = shape(dx); + if (tags & BC::E) { + auto ex2_U { ZERO }; + if constexpr (defines_ex2) { + ex2_U = metric.template transform<2, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.ex2(x_Ph_0H)); + } + // ex2 + Fld(i1, i2, em::ex2) = s * Fld(i1, i2, em::ex2) + (ONE - s) * ex2_U; + } + if (tags & BC::B) { + auto bx1_U { ZERO }; + if constexpr (defines_bx1) { + bx1_U = metric.template transform<1, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.bx1(x_Ph_0H)); + } + // bx1 + Fld(i1, i2, em::bx1) = s * Fld(i1, i2, em::bx1) + (ONE - s) * bx1_U; + } + } + if (tags & BC::E) { + auto ex3_U { ZERO }; + if constexpr (defines_ex3) { + 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); + // ex3 + Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3) + (ONE - s) * ex3_U; + } + if (tags & BC::B) { + auto bx3_U { ZERO }; + if constexpr (defines_bx3) { + 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; + } + } else { + // GRPIC + raise::KernelError(HERE, "GRPIC not implemented"); } } else { raise::KernelError( HERE, - "AbsorbFields_kernel: 2D implementation called for D != 2"); + "MatchBoundaries_kernel: 2D implementation called for D != 2"); } } @@ -122,55 +281,180 @@ namespace kernel { const auto i1_ = COORD(i1); const auto i2_ = COORD(i2); const auto i3_ = COORD(i3); - for (const auto comp : - { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) { - if ((comp == em::ex1) and not(tags & BC::Ex1)) { - continue; - } else if ((comp == em::ex2) and not(tags & BC::Ex2)) { - continue; - } else if ((comp == em::ex3) and not(tags & BC::Ex3)) { - continue; - } else if ((comp == em::bx1) and not(tags & BC::Bx1)) { - continue; - } else if ((comp == em::bx2) and not(tags & BC::Bx2)) { - continue; - } else if ((comp == em::bx3) and not(tags & BC::Bx3)) { - continue; - } - coord_t x_Cd { ZERO }; - if (comp == em::ex1) { - x_Cd[0] = i1_ + HALF; - x_Cd[1] = i2_; - x_Cd[2] = i3_; - } else if (comp == em::ex2) { - x_Cd[0] = i1_; - x_Cd[1] = i2_ + HALF; - x_Cd[2] = i3_; - } else if (comp == em::ex3) { - x_Cd[0] = i1_; - x_Cd[1] = i2_; - x_Cd[2] = i3_ + HALF; - } else if (comp == em::bx1) { - x_Cd[0] = i1_; - x_Cd[1] = i2_ + HALF; - x_Cd[2] = i3_ + HALF; - } else if (comp == em::bx2) { - x_Cd[0] = i1_ + HALF; - x_Cd[1] = i2_; - x_Cd[2] = i3_ + HALF; - } else if (comp == em::bx3) { - x_Cd[0] = i1_ + HALF; - x_Cd[1] = i2_ + HALF; - x_Cd[2] = i3_; - } - const auto dx = math::abs( - metric.template convert(x_Cd[i - 1]) - xg_edge); - Fld(i1, i2, i3, comp) *= math::tanh(dx / (INV_4 * dx_abs)); + + if constexpr (S == SimEngine::SRPIC) { + // SRPIC + if (tags & BC::E) { + { + // i1 + 1/2, i2, i3 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_ + HALF; + } else if constexpr (o == in::x2) { + xi_Cd = i2_; + } 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 }; + if constexpr (defines_ex1) { + 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)); + } + // ex1 + Fld(i1, i2, i3, em::ex1) = s * Fld(i1, i2, i3, em::ex1) + + (ONE - s) * ex1_U; + } + { + // i1, i2 + 1/2, i3 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_; + } else if constexpr (o == in::x2) { + xi_Cd = i2_ + HALF; + } 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 }; + if constexpr (defines_ex2) { + 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)); + } + // ex2 + Fld(i1, i2, i3, em::ex2) = s * Fld(i1, i2, i3, em::ex2) + + (ONE - s) * ex2_U; + } + { + // i1, i2, i3 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_; + } else if constexpr (o == in::x2) { + xi_Cd = i2_; + } 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 }; + if constexpr (defines_ex3) { + 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)); + } + // ex3 + Fld(i1, i2, i3, em::ex3) = s * Fld(i1, i2, i3, em::ex3) + + (ONE - s) * ex3_U; + } + } + if (tags & BC::B) { + { + // i1, i2 + 1/2, i3 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_; + } else if constexpr (o == in::x2) { + xi_Cd = i2_ + HALF; + } 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; + } + { + // i1 + 1/2, i2, i3 + 1/2 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_ + HALF; + } else if constexpr (o == in::x2) { + xi_Cd = i2_; + } 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 }; + if constexpr (defines_bx2) { + 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)); + } + // bx2 + Fld(i1, i2, i3, em::bx2) = s * Fld(i1, i2, i3, em::bx2) + + (ONE - s) * bx2_U; + } + { + // i1 + 1/2, i2 + 1/2, i3 + real_t xi_Cd; + if constexpr (o == in::x1) { + xi_Cd = i1_ + HALF; + } else if constexpr (o == in::x2) { + xi_Cd = i2_ + HALF; + } 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 }; + if constexpr (defines_bx3) { + 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)); + } + // bx3 + Fld(i1, i2, i3, em::bx3) = s * Fld(i1, i2, i3, em::bx3) + + (ONE - s) * bx3_U; + } + } + } else { + // GRPIC + raise::KernelError(HERE, "GRPIC not implemented"); } } else { raise::KernelError( HERE, - "AbsorbFields_kernel: 3D implementation called for D != 3"); + "MatchBoundaries_kernel: 3D implementation called for D != 3"); } } }; @@ -226,31 +510,28 @@ namespace kernel { static constexpr bool defines_bx2 = traits::has_method::value; static constexpr bool defines_bx3 = traits::has_method::value; - static_assert(defines_ex1 and defines_ex2 and defines_ex3 and - defines_bx1 and defines_bx2 and defines_bx3, - "not all components of E or B are specified in PGEN"); + static_assert(defines_ex1 or defines_ex2 or defines_ex3 or defines_bx1 or + defines_bx2 or defines_bx3, + "none of the components of E or B are specified in PGEN"); static_assert(M::is_metric, "M must be a metric class"); static_assert(static_cast(O) < static_cast(M::Dim), "Invalid Orientation"); ndfield_t Fld; - const I finit; + const I fset; const M metric; const std::size_t i_edge; - const bool setE, setB; EnforcedBoundaries_kernel(ndfield_t& Fld, - const I& finit, + const I& fset, const M& metric, std::size_t i_edge, BCTags tags) : Fld { Fld } - , finit { finit } + , fset { fset } , metric { metric } - , i_edge { i_edge + N_GHOSTS } - , setE { tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3 } - , setB { tags & BC::Bx1 or tags & BC::Bx2 or tags & BC::Bx3 } {} + , i_edge { i_edge + N_GHOSTS } {} Inline void operator()(index_t i1) const { if constexpr (D == Dim::_1D) { @@ -259,8 +540,8 @@ namespace kernel { coord_t x_Ph_H { ZERO }; metric.template convert({ i1_ }, x_Ph_0); metric.template convert({ i1_ + HALF }, x_Ph_H); - bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB, - setBx2 = setB, setBx3 = setB; + bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3, + setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3; if constexpr (O == in::x1) { // x1 -- normal // x2,x3 -- tangential @@ -276,35 +557,47 @@ namespace kernel { } else { raise::KernelError(HERE, "Invalid Orientation"); } - if (setEx1) { - Fld(i1, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF }, - finit.ex1(x_Ph_H)); + if constexpr (defines_ex1) { + if (setEx1) { + Fld(i1, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.ex1(x_Ph_H)); + } } - if (setEx2) { - Fld(i1, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_ }, - finit.ex2(x_Ph_0)); + if constexpr (defines_ex2) { + if (setEx2) { + Fld(i1, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_ }, + fset.ex2(x_Ph_0)); + } } - if (setEx3) { - Fld(i1, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_ }, - finit.ex3(x_Ph_0)); + if constexpr (defines_ex3) { + if (setEx3) { + Fld(i1, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_ }, + fset.ex3(x_Ph_0)); + } } - if (setBx1) { - Fld(i1, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_ }, - finit.bx1(x_Ph_0)); + if constexpr (defines_bx1) { + if (setBx1) { + Fld(i1, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_ }, + fset.bx1(x_Ph_0)); + } } - if (setBx2) { - Fld(i1, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF }, - finit.bx2(x_Ph_H)); + if constexpr (defines_bx2) { + if (setBx2) { + Fld(i1, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.bx2(x_Ph_H)); + } } - if (setBx3) { - Fld(i1, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF }, - finit.bx3(x_Ph_H)); + if constexpr (defines_bx3) { + if (setBx3) { + Fld(i1, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF }, + fset.bx3(x_Ph_H)); + } } } else { raise::KernelError(HERE, "Invalid Dimension"); @@ -324,8 +617,8 @@ namespace kernel { metric.template convert({ i1_ + HALF, i2_ }, x_Ph_H0); metric.template convert({ i1_ + HALF, i2_ + HALF }, x_Ph_HH); - bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB, - setBx2 = setB, setBx3 = setB; + bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3, + setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3; if constexpr (O == in::x1) { // x1 -- normal // x2,x3 -- tangential @@ -353,35 +646,47 @@ namespace kernel { } else { raise::KernelError(HERE, "Invalid Orientation"); } - if (setEx1) { - Fld(i1, i2, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF, i2_ }, - finit.ex1(x_Ph_H0)); + if constexpr (defines_ex1) { + if (setEx1) { + Fld(i1, i2, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.ex1(x_Ph_H0)); + } } - if (setEx2) { - Fld(i1, i2, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_, i2_ + HALF }, - finit.ex2(x_Ph_0H)); + if constexpr (defines_ex2) { + if (setEx2) { + Fld(i1, i2, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.ex2(x_Ph_0H)); + } } - if (setEx3) { - Fld(i1, i2, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_ }, - finit.ex3(x_Ph_00)); + if constexpr (defines_ex3) { + if (setEx3) { + Fld(i1, i2, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_ }, + fset.ex3(x_Ph_00)); + } } - if (setBx1) { - Fld(i1, i2, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_, i2_ + HALF }, - finit.bx1(x_Ph_0H)); + if constexpr (defines_bx1) { + if (setBx1) { + Fld(i1, i2, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_, i2_ + HALF }, + fset.bx1(x_Ph_0H)); + } } - if (setBx2) { - Fld(i1, i2, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF, i2_ }, - finit.bx2(x_Ph_H0)); + if constexpr (defines_bx2) { + if (setBx2) { + Fld(i1, i2, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF, i2_ }, + fset.bx2(x_Ph_H0)); + } } - if (setBx3) { - Fld(i1, i2, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF, i2_ + HALF }, - finit.bx3(x_Ph_HH)); + if constexpr (defines_bx3) { + if (setBx3) { + Fld(i1, i2, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF, i2_ + HALF }, + fset.bx3(x_Ph_HH)); + } } } else { raise::KernelError(HERE, "Invalid Dimension"); @@ -412,8 +717,8 @@ namespace kernel { x_Ph_H0H); metric.template convert({ i1_, i2_ + HALF, i3_ + HALF }, x_Ph_0HH); - bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB, - setBx2 = setB, setBx3 = setB; + bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3, + setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3; if constexpr (O == in::x1) { // x1 -- normal // x2,x3 -- tangential @@ -453,35 +758,47 @@ namespace kernel { } else { raise::KernelError(HERE, "Invalid Orientation"); } - if (setEx1) { - Fld(i1, i2, i3, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_ + HALF, i2_, i3_ }, - finit.ex1(x_Ph_H00)); + if constexpr (defines_ex1) { + if (setEx1) { + Fld(i1, i2, i3, em::ex1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_ + HALF, i2_, i3_ }, + fset.ex1(x_Ph_H00)); + } } - if (setEx2) { - Fld(i1, i2, i3, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_, i2_ + HALF, i3_ }, - finit.ex2(x_Ph_0H0)); + if constexpr (defines_ex2) { + if (setEx2) { + Fld(i1, i2, i3, em::ex2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_, i2_ + HALF, i3_ }, + fset.ex2(x_Ph_0H0)); + } } - if (setEx3) { - Fld(i1, i2, i3, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_, i3_ + HALF }, - finit.ex3(x_Ph_00H)); + if constexpr (defines_ex3) { + if (setEx3) { + Fld(i1, i2, i3, em::ex3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_, i3_ + HALF }, + fset.ex3(x_Ph_00H)); + } } - if (setBx1) { - Fld(i1, i2, i3, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( - { i1_, i2_ + HALF, i3_ + HALF }, - finit.bx1(x_Ph_0HH)); + if constexpr (defines_bx1) { + if (setBx1) { + Fld(i1, i2, i3, em::bx1) = metric.template transform<1, Idx::T, Idx::U>( + { i1_, i2_ + HALF, i3_ + HALF }, + fset.bx1(x_Ph_0HH)); + } } - if (setBx2) { - Fld(i1, i2, i3, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( - { i1_ + HALF, i2_, i3_ + HALF }, - finit.bx2(x_Ph_H0H)); + if constexpr (defines_bx2) { + if (setBx2) { + Fld(i1, i2, i3, em::bx2) = metric.template transform<2, Idx::T, Idx::U>( + { i1_ + HALF, i2_, i3_ + HALF }, + fset.bx2(x_Ph_H0H)); + } } - if (setBx3) { - Fld(i1, i2, i3, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( - { i1_ + HALF, i2_ + HALF, i3_ }, - finit.bx3(x_Ph_HH0)); + if constexpr (defines_bx3) { + if (setBx3) { + Fld(i1, i2, i3, em::bx3) = metric.template transform<3, Idx::T, Idx::U>( + { i1_ + HALF, i2_ + HALF, i3_ }, + fset.bx3(x_Ph_HH0)); + } } } else { raise::KernelError(HERE, "Invalid Dimension");