KfFieldRegion.h

Go to the documentation of this file.

/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */

 
#pragma once
 
#include "KfDefs.h"
#include "KfFieldValue.h"
#include "KfMath.h"
#include "KfUtils.h"
#include "KfVisit.h"
 
#include <boost/serialization/access.hpp>
#include <boost/serialization/array.hpp>
#include <boost/serialization/optional.hpp>
#include <boost/serialization/split_free.hpp>
 
#include <array>
#include <optional>
#include <string>
#include <type_traits>
 
namespace cbm::algo::kf
{

  class GlobalField {
   public:
    ~GlobalField() = delete;

    static FieldFn_t fgOriginalField;

    static EFieldType fgOriginalFieldType;  

    static std::tuple<double, double, double> UndefField(double, double, double);

    static void SetFieldFunction(EFieldType fldType, const FieldFn_t& fn)
    {
      fgOriginalFieldType = fldType;
      fgOriginalField     = fn;
    }

    template<typename T>
    [[gnu::always_inline]] static FieldValue<T> GetFieldValue(const FieldFn_t& fieldFn, T x, T y, T z);

    static void ForceUseOfOriginalField(bool v = true) { fgForceUseOfOriginalField = v; }

    static bool IsUsingOriginalFieldForced() { return fgForceUseOfOriginalField; }
 
   private:
    static bool fgForceUseOfOriginalField;  
  };
 
 
  namespace detail
  {
    template<typename T, EFieldMode FldMode>
    class alignas(VcMemAlign) ConcreteFieldRegion {
    };
 

    template<typename T>
    class alignas(VcMemAlign) ConcreteFieldRegion<T, EFieldMode::Interpolated> {
      template<typename, EFieldMode>
      friend class ConcreteFieldRegion;
 
      using CoeffArray_t = std::array<std::array<T, 3>, 3>;  
 
     public:
      using FieldValue_t = FieldValue<T>;

      ConcreteFieldRegion() = default;

      ConcreteFieldRegion(const FieldValue_t& b0, const FieldValue_t& b1, const FieldValue_t& b2)
      {
        this->Set(b0, b1, b2);
      }

      template<typename I>
      ConcreteFieldRegion(const ConcreteFieldRegion<I, EFieldMode::Interpolated>& other)
        : fZfirst(utils::simd::Cast<I, T>(other.fZfirst))
        , fFieldType(other.fFieldType)
      {
        for (size_t iDim = 0; iDim < fCoeff.size(); ++iDim) {
          for (size_t iPow = 0; iPow < fCoeff[iDim].size(); ++iPow) {
            fCoeff[iDim][iPow] = utils::simd::Cast<I, T>(other.fCoeff[iDim][iPow]);
          }
        }
      }

      ConcreteFieldRegion(const ConcreteFieldRegion&) = default;

      ConcreteFieldRegion(ConcreteFieldRegion&&) = default;

      ~ConcreteFieldRegion() = default;

      ConcreteFieldRegion& operator=(const ConcreteFieldRegion&) = default;

      ConcreteFieldRegion& operator=(ConcreteFieldRegion&&) = default;

      FieldValue<T> Get(T x, T y, T z) const;

      const auto& GetCoefficients() const { return fCoeff; }

      std::tuple<T, T, T> GetDoubleIntegrals(T x1, T y1, T z1, T x2, T y2, T z2) const;

      EFieldType GetFieldType() const { return fFieldType; }

      const T& GetZfirst() const { return fZfirst; }

      void Set(const FieldValue<T>& b0, const FieldValue<T>& b1, const FieldValue<T>& b2);

      void Set(const CoeffArray_t& coeff, T z0)
      {
        fCoeff  = coeff;
        fZfirst = z0;
      }

      void SetFieldType(EFieldType type) { fFieldType = type; }

      void Shift(T z);

      void CheckConsistency() const;

      std::string ToString(int indentLevel = 0, int verbose = 1) const;

      static constexpr EFieldMode Mode() { return EFieldMode::Interpolated; }
 
     private:
      CoeffArray_t fCoeff{};                      
      T fZfirst{Literal_t<T>()};                  
      EFieldType fFieldType{EFieldType::Normal};  

      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive& ar, const unsigned int)
      {
        ar& fCoeff;
        ar& fZfirst;
        ar& fFieldType;
      }
    };
 

    template<typename T>
    class alignas(VcMemAlign) ConcreteFieldRegion<T, EFieldMode::Original> {
      template<typename, EFieldMode>
      friend class ConcreteFieldRegion;
 
     public:
      ConcreteFieldRegion() = default;

      ConcreteFieldRegion(EFieldType type, const FieldFn_t& fieldFn) : fFieldFn(fieldFn), fFieldType(type) {}

      template<typename I>
      ConcreteFieldRegion(const ConcreteFieldRegion<I, EFieldMode::Original>& other)
        : fFieldFn(other.fFieldFn)
        , fFieldType(other.fFieldType)
      {
      }

      ConcreteFieldRegion(const ConcreteFieldRegion&) = default;

      ConcreteFieldRegion(ConcreteFieldRegion&&) = default;

      ~ConcreteFieldRegion() = default;

      ConcreteFieldRegion& operator=(const ConcreteFieldRegion&) = default;

      ConcreteFieldRegion& operator=(ConcreteFieldRegion&&) = default;

      FieldValue<T> Get(T x, T y, T z) const { return GlobalField::GetFieldValue(fFieldFn, x, y, z); }

      std::tuple<T, T, T> GetDoubleIntegrals(T /*x1*/, T /*y1*/, T /*z1*/, T /*x2*/, T /*y2*/, T /*z2*/) const
      {
        assert("cbm::algo::kf::ConcreteFieldRegion<T, EFieldMode::Original>::GetDoubleIntegrals() is not implemented");
        return {defs::Undef<T>, defs::Undef<T>, defs::Undef<T>};
      }

      EFieldType GetFieldType() const { return fFieldType; }

      void CheckConsistency() const;

      // TODO: remove?
      void SetFieldType(EFieldType type) { fFieldType = type; }

      std::string ToString(int indentLevel = 0, int verbose = 1) const;

      static constexpr EFieldMode Mode() { return EFieldMode::Original; }
 
     private:
      FieldFn_t fFieldFn{defs::ZeroFieldFn};      
      EFieldType fFieldType{EFieldType::Normal};  
    };
  }  // namespace detail
 

  template<typename T>
  class alignas(VcMemAlign) FieldRegion {
    template<typename>
    friend class FieldRegion;
 
   public:
    template<EFieldMode Mode>
    using ConcreteFieldRegion_t = typename detail::ConcreteFieldRegion<T, Mode>;
    using Interpolated_t        = ConcreteFieldRegion_t<EFieldMode::Interpolated>;
    using Original_t            = ConcreteFieldRegion_t<EFieldMode::Original>;
    using FieldRegionVariant_t  = typename std::variant<Interpolated_t, Original_t>;

    FieldRegion() = default;

    FieldRegion(const FieldRegion&) = default;

    FieldRegion(FieldRegion&&) = default;

    FieldRegion(const FieldValue<T>& b0, const FieldValue<T>& b1, const FieldValue<T>& b2)
      : fFrVariant(FieldRegionVariant_t(std::in_place_type<Interpolated_t>, b0, b1, b2))
    {
    }

    FieldRegion(EFieldType type, const FieldFn_t& fieldFn)
      : fFrVariant(FieldRegionVariant_t(std::in_place_type<Original_t>, type, fieldFn))
    {
    }


    ~FieldRegion() = default;

    FieldRegion& operator=(const FieldRegion&) = default;

    FieldRegion& operator=(FieldRegion&&) = default;

    template<typename I>
    static FieldRegion CopyConvert(const FieldRegion<I>& other)
    {
      return kf::visit(
        [&other](const auto& c) -> FieldRegion { return FieldRegion(other, Tag<std::decay_t<decltype(c)>::Mode()>{}); },
        other.fFrVariant);
    }

    FieldValue<T> Get(T x, T y, T z) const
    {
      return kf::visit([&x, &y, &z](const auto& c) -> FieldValue<T> { return c.Get(x, y, z); }, fFrVariant);
    }

    template<EFieldMode FieldMode>
    ConcreteFieldRegion_t<FieldMode>* GetConcrete()
    {
      return std::get_if<ConcreteFieldRegion_t<FieldMode>>(&fFrVariant);
    }

    template<EFieldMode FieldMode>
    const ConcreteFieldRegion_t<FieldMode>* GetConcrete() const
    {
      return std::get_if<ConcreteFieldRegion_t<FieldMode>>(&fFrVariant);
    }

    std::tuple<T, T, T> GetDoubleIntegrals(T x1, T y1, T z1, T x2, T y2, T z2) const
    {
      return kf::visit([&x1, &y1, &z1, &x2, &y2, &z2](
                         const auto& c) -> std::tuple<T, T, T> { return c.GetDoubleIntegrals(x1, y1, z1, x2, y2, z2); },
                       fFrVariant);
    }

    EFieldMode GetFieldMode() const
    {
      return kf::visit([](const auto& c) { return std::decay_t<decltype(c)>::Mode(); }, fFrVariant);
    }

    EFieldType GetFieldType() const
    {
      return kf::visit([](const auto& c) -> EFieldType { return c.GetFieldType(); }, fFrVariant);
    }

    template<EFieldMode FieldMode, typename... Args>
    static constexpr FieldRegion Create(Args... args)
    {
      return FieldRegion(args..., Tag<FieldMode>{});
    }

    template<typename... Args>
    static constexpr FieldRegion Create(EFieldMode mode, Args... args)
    {
      switch (mode) {
        case EFieldMode::Interpolated: return FieldRegion(args..., Tag<EFieldMode::Interpolated>{});
        case EFieldMode::Original:
        default: return FieldRegion(args..., Tag<EFieldMode::Original>{});
      }
      __builtin_unreachable();  // No return is reachable
    }

    void Set(const FieldValue<T>& b0, const FieldValue<T>& b1, const FieldValue<T>& b2)
    {
      kf::visit(
        [&b0, &b1, &b2](auto&& c) -> void {
          if constexpr (std::decay_t<decltype(c)>::Mode() == EFieldMode::Interpolated) {
            c.Set(b0, b1, b2);
          }
        },
        fFrVariant);
    }

    void SetFieldType(EFieldType type)
    {
      kf::visit([type](auto&& c) { c.SetFieldType(type); }, fFrVariant);
    }

    void Shift(T z)
    {
      kf::visit(
        [&z](auto&& c) -> void {
          if constexpr (std::decay_t<decltype(c)>::Mode() == EFieldMode::Interpolated) {
            c.Shift(z);
          }
        },
        fFrVariant);
    }

    void CheckConsistency() const
    {
      kf::visit([](const auto& c) -> void { return c.CheckConsistency(); }, fFrVariant);
    }

    std::string ToString(int indentLevel = 0, int verbose = 1) const
    {
      return kf::visit(
        [indentLevel, verbose](const auto& c) -> std::string { return c.ToString(indentLevel, verbose); }, fFrVariant);
    }
 
   private:
    template<typename I, EFieldMode FieldMode>
    FieldRegion(const FieldRegion<I>& other, Tag<FieldMode>)
      : fFrVariant(FieldRegionVariant_t(
        std::in_place_type<ConcreteFieldRegion_t<FieldMode>>,
        ConcreteFieldRegion_t<FieldMode>(std::get<detail::ConcreteFieldRegion<I, FieldMode>>(other.fFrVariant))))
    {
    }

    template<class Archive>
    void load(Archive& ar, const unsigned int /*version*/)
    {
      auto fieldRegion = Interpolated_t{};
      ar >> fieldRegion;
      fFrVariant.template emplace<Interpolated_t>(std::move(fieldRegion));
    }

    template<class Archive>
    void save(Archive& ar, const unsigned int /*version*/) const
    {
      if (std::holds_alternative<Interpolated_t>(fFrVariant)) {
        ar << std::get<Interpolated_t>(fFrVariant);
      }
      else {
        throw std::logic_error("Attempt to serialize a kf::FieldRegion object with fFieldMode == EFieldMode::Original");
      }
    }

    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int version)
    {
      boost::serialization::split_member(ar, *this, version);
    }
 
    FieldRegionVariant_t fFrVariant;  
  };
 
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<typename T>
  inline FieldValue<T> GlobalField::GetFieldValue(const FieldFn_t& fieldFn, T x, T y, T z)
  {
    FieldValue<T> B;
    using utils::simd::Cast;
    for (size_t i = 0; i < utils::simd::Size<T>(); ++i) {
      auto [bx, by, bz] = fieldFn(Cast<T, double>(x, i), Cast<T, double>(y, i), Cast<T, double>(z, i));
      B.SetSimdEntry(bx, by, bz, Cast<T, double>(z, i), i);
    }
    return B;
  }
 
}  // namespace cbm::algo::kf