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 <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:
    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, const T& x, const T& y,
                                                              const 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) FieldRegionBase {
    };
 
    template<typename T>
    class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Orig> {
      template<typename, EFieldMode>
      friend class FieldRegionBase;
 
     public:
      FieldRegionBase() = default;
 
      FieldRegionBase(const FieldFn_t& fieldFn) : fFieldFn(fieldFn) {}
 
      template<typename I>
      FieldRegionBase(const FieldRegionBase<I, EFieldMode::Orig>& other) : fFieldFn(other.fFieldFn)
      {
      }
 
      ~FieldRegionBase() = default;
 
      FieldRegionBase& operator=(const FieldRegionBase& other) = default;
 
      FieldValue<T> Get(const T& x, const T& y, const T& z) const
      {
        return GlobalField::GetFieldValue(fFieldFn, x, y, z);
      }
 
      std::tuple<T, T, T> GetDoubleIntegrals(const T& /*x1*/, const T& /*y1*/, const T& /*z1*/,  //
                                             const T& /*x2*/, const T& /*y2*/, const T& /*z2*/) const
      {
        assert("cbm::algo::kf::FieldRegionBase<T, EFieldMode::Orig>::GetDoubleIntegrals() is not implemented");
        return {defs::Undef<T>, defs::Undef<T>, defs::Undef<T>};
      }
 
      void CheckConsistency() const;
 
      std::string ToString(int indentLevel = 0, int verbose = 1) const;
 
     protected:
      FieldFn_t fFieldFn{defs::ZeroFieldFn};  
    };
 
 
    template<typename T>
    class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Intrpl> {
      template<typename, EFieldMode>
      friend class FieldRegionBase;
 
      using CoeffArray_t = std::array<std::array<T, 3>, 3>;  
 
     public:
      using FieldValue_t = FieldValue<T>;
 
      FieldRegionBase() = default;
 
      FieldRegionBase(const FieldValue_t& b0, const T& z0, const FieldValue_t& b1, const T& z1, const FieldValue_t& b2,
                      const T& z2)
      {
        this->Set(b0, z0, b1, z1, b2, z2);
      }
 
      template<typename I>
      FieldRegionBase(const FieldRegionBase<I, EFieldMode::Intrpl>& other)
        : fZfirst(utils::simd::Cast<I, T>(other.fZfirst))
      {
        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]);
          }
        }
      }
 
      ~FieldRegionBase() = default;
 
      FieldRegionBase& operator=(const FieldRegionBase& other) = default;
 
      FieldValue<T> Get(const T& x, const T& y, const T& z) const;
 
      std::tuple<T, T, T> GetDoubleIntegrals(const T& x1, const T& y1, const T& z1,  //
                                             const T& x2, const T& y2, const T& z2) const;
 
      const T& GetZfirst() const { return fZfirst; }
 
      const auto& GetCoefficients() const { return fCoeff; }
 
      void Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
               const T& z2);
 
      void Set(const CoeffArray_t& coeff, const T& z0)
      {
        fCoeff  = coeff;
        fZfirst = z0;
      }
 
      void Shift(const T& z);
 
      void CheckConsistency() const;
 
      std::string ToString(int indentLevel = 0, int verbose = 1) const;
 
     protected:
      CoeffArray_t fCoeff{{T(0.)}};  
      T fZfirst{Literal_t<T>()};     
 
     private:
      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive& ar, const unsigned int)
      {
        ar& fCoeff;
        ar& fZfirst;
      }
    };
  }  // namespace detail
 
 
  template<typename T>
  class alignas(VcMemAlign) FieldRegion {
    template<typename>
    friend class FieldRegion;
 
   public:
    FieldRegion(EFieldMode fldMode = EFieldMode::Intrpl, EFieldType fldType = EFieldType::Normal)
      : foFldIntrpl(fldMode == EFieldMode::Intrpl ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>()
                                                  : std::nullopt)
      , foFldOrig(fldMode == EFieldMode::Orig ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Orig>>()
                                              : std::nullopt)
      , fFieldType(fldType)
      , fFieldMode(fldMode)
    {
    }
 
    template<typename I>
    FieldRegion(const FieldRegion<I>& other)
      : foFldIntrpl(other.foFldIntrpl.has_value()
                      ? std::make_optional(detail::FieldRegionBase<T, EFieldMode::Intrpl>(*other.foFldIntrpl))
                      : std::nullopt)
      , foFldOrig(other.foFldOrig.has_value()
                    ? std::make_optional(detail::FieldRegionBase<T, EFieldMode::Orig>(*other.foFldOrig))
                    : std::nullopt)
      , fFieldType(other.fFieldType)
      , fFieldMode(other.fFieldMode)
    {
    }
 
    FieldRegion(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
                const T& z2)
      : foFldIntrpl(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>())
      , foFldOrig(std::nullopt)
      , fFieldType(EFieldType::Normal)
      , fFieldMode(EFieldMode::Intrpl)
    {
      // Set the field type automatically
      Set(b0, z0, b1, z1, b2, z2);
    }
 
    FieldRegion(EFieldType fieldType, const detail::FieldRegionBase<T, EFieldMode::Intrpl>& fld)
      : foFldIntrpl(std::make_optional(fld))
      , foFldOrig(std::nullopt)
      , fFieldType(fieldType)
      , fFieldMode(EFieldMode::Intrpl)
    {
    }
 
    FieldRegion(EFieldType fieldType, const FieldFn_t& fieldFn)
      : foFldIntrpl(std::nullopt)
      , foFldOrig(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Orig>>(fieldFn))
      , fFieldType(fieldType)
      , fFieldMode(EFieldMode::Orig)
    {
    }
 
    ~FieldRegion() = default;
 
    FieldRegion& operator=(const FieldRegion& other) = default;
 
    FieldValue<T> Get(const T& x, const T& y, const T& z) const
    {
      return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->Get(x, y, z) : foFldOrig->Get(x, y, z);
    }
 
    std::tuple<T, T, T> GetDoubleIntegrals(const T& x1, const T& y1, const T& z1, const T& x2, const T& y2,
                                           const T& z2) const
    {
      return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->GetDoubleIntegrals(x1, y1, z1, x2, y2, z2)
                                              : foFldOrig->GetDoubleIntegrals(x1, y1, z1, x2, y2, z2);
    }
 
    EFieldMode GetFieldMode() const { return fFieldMode; }
 
    EFieldType GetFieldType() const { return fFieldType; }
 
    void Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
             const T& z2)
    {
      fFieldMode = EFieldMode::Intrpl;
      fFieldType = EFieldType::Normal;
      // Set the field type automatically
      if (b0.IsZero() && b1.IsZero() && b2.IsZero()) {
        fFieldType = EFieldType::Null;
      }
      foFldIntrpl->Set(b0, z0, b1, z1, b2, z2);
    }
 
    void SetFieldType(EFieldType fldType) { fFieldType = fldType; }
 
    void Shift(const T& z);
 
    void CheckConsistency() const;
 
    std::string ToString(int indentLevel = 0, int verbose = 1) const;
 
    const auto& GetIntrplField() const { return foFldIntrpl; }
 
   private:
    friend class boost::serialization::access;
 
    template<class Archive>
    void load(Archive& ar, const unsigned int /*version*/)
    {
      auto field = detail::FieldRegionBase<T, EFieldMode::Intrpl>{};
      ar >> field;
      foFldIntrpl = std::make_optional(field);
      ar >> fFieldType;
      ar >> fFieldMode;
      foFldOrig = std::nullopt;  // Note: original field cannot be serialized
    }
 
    template<class Archive>
    void save(Archive& ar, const unsigned int /*version*/) const
    {
      if (fFieldMode == EFieldMode::Intrpl) {
        ar << (*foFldIntrpl);
        ar << fFieldType;
        ar << fFieldMode;
      }
      else if (fFieldMode == EFieldMode::Orig) {
        throw std::logic_error("Attempt to serialize a kf::FieldRegion object with fFieldMode == EFieldMode::Orig");
      }
    }
 
    template<class Archive>
    void serialize(Archive& ar, const unsigned int version)
    {
      boost::serialization::split_member(ar, *this, version);
    }
 
    std::optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>> foFldIntrpl{std::nullopt};
    std::optional<detail::FieldRegionBase<T, EFieldMode::Orig>> foFldOrig{std::nullopt};
    EFieldType fFieldType{EFieldType::Null};  
    EFieldMode fFieldMode;                    
  };
 
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<typename T>
  inline FieldValue<T> GlobalField::GetFieldValue(const FieldFn_t& fieldFn, const T& x, const T& y, const 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, i);
    }
    return B;
  }
 
}  // namespace cbm::algo::kf