KfField.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 "KfFieldRegion.h"
#include "KfFieldSlice.h"
 
#include <boost/serialization/access.hpp>
#include <boost/serialization/split_free.hpp>
 
#include <array>
#include <optional>
#include <set>
 
namespace cbm::algo::kf
{
  namespace detail
  {
    template<typename T, EFieldMode FldMode>
    class alignas(VcMemAlign) FieldBase {
    };
 
    template<typename T>
    class alignas(VcMemAlign) FieldBase<T, EFieldMode::Orig> {
      template<typename, EFieldMode>
      friend class FieldBase;
 
     public:
      FieldBase() = default;
 
      template<typename I>
      FieldBase(const FieldBase<I, EFieldMode::Orig>& other) : fFieldFn(other.fFieldFn)
      {
        fvFieldSliceZ.reserve(other.fvFieldSliceZ.size());
        for (const auto& slice : other.fvFieldSliceZ) {
          fvFieldSliceZ.emplace_back(utils::simd::Cast<I, T>(slice));
        }
      }
 
      FieldBase& operator=(const FieldBase& other)
      {
        if (this != &other) {
          fvFieldSliceZ = other.fvFieldSliceZ;
          fFieldFn      = other.fFieldFn;
        }
        return *this;
      }
 
      void AddFieldSlice(const T& fieldSliceZref) { fvFieldSliceZ.push_back(fieldSliceZref); }
 
      FieldValue<T> GetFieldValue(int sliceID, const T& x, const T& y) const
      {
        return GlobalField::GetFieldValue(fFieldFn, x, y, fvFieldSliceZ[sliceID]);
      }
 
      const FieldFn_t& GetFieldFunction() const { return fFieldFn; }
 
      int GetNofFieldSlices() const { return fvFieldSliceZ.size(); }
 
      void RemoveSlice(int iLayer) { fvFieldSliceZ.erase(fvFieldSliceZ.begin() + iLayer); }
 
      void SetFieldFunction(const FieldFn_t& fieldFn) { fFieldFn = fieldFn; }
 
      std::string ToString(int indentLevel, int verbose) const;
 
     private:
      std::vector<T> fvFieldSliceZ{};         
      FieldFn_t fFieldFn{defs::ZeroFieldFn};  
    };
 
 
    template<typename T>
    class alignas(VcMemAlign) FieldBase<T, EFieldMode::Intrpl> {
      template<typename, EFieldMode>
      friend class FieldBase;
 
      using SlicesContainer_t = std::vector<FieldSlice<T>>;
 
     public:
      FieldBase() = default;
 
      template<typename I>
      FieldBase(const FieldBase<I, EFieldMode::Intrpl>& other)
      {
        fvFieldSlices.reserve(other.fvFieldSlices.size());
        for (const auto& slice : other.fvFieldSlices) {
          fvFieldSlices.emplace_back(FieldSlice<T>(slice));
        }
      }
 
      FieldBase& operator=(const FieldBase& other)
      {
        if (this != &other) {
          fvFieldSlices = other.fvFieldSlices;
        }
        return *this;
      }
 
      void AddFieldSlice(const FieldSlice<T>& fieldSlice) { fvFieldSlices.push_back(fieldSlice); }
 
      const auto& GetFieldSlice(int sliceID) const { return fvFieldSlices[sliceID]; }
 
      FieldValue<T> GetFieldValue(int sliceID, const T& x, const T& y) const
      {
        return fvFieldSlices[sliceID].GetFieldValue(x, y);
      }
 
      int GetNofFieldSlices() const { return fvFieldSlices.size(); }
 
      void RemoveSlice(int iLayer) { fvFieldSlices.erase(fvFieldSlices.begin() + iLayer); }
 
      std::string ToString(int indentLevel, int verbose) const;
 
     protected:
      SlicesContainer_t fvFieldSlices{};  
 
     private:
      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive& ar, const unsigned int)
      {
        ar& fvFieldSlices;
      }
    };
  }  // namespace detail
 
 
  template<typename T>
  class alignas(VcMemAlign) Field {
    template<typename>
    friend class Field;
    friend class FieldFactory;
    friend class boost::serialization::access;
 
   public:
    Field(EFieldMode fldMode, EFieldType fldType)
      : foFldIntrpl(fldMode == EFieldMode::Intrpl ? std::make_optional(detail::FieldBase<T, EFieldMode::Intrpl>())
                                                  : std::nullopt)
      , foFldOrig(fldMode == EFieldMode::Orig ? std::make_optional(detail::FieldBase<T, EFieldMode::Orig>())
                                              : std::nullopt)
      , fPrimVertexField(FieldRegion<T>(fldMode, fldType))
      , fFieldType(fldType)
      , fFieldMode(fldMode)
    {
    }
 
    template<typename I>
    Field(const Field<I>& other)
      : foFldIntrpl(other.foFldIntrpl.has_value()
                      ? std::make_optional(detail::FieldBase<T, EFieldMode::Intrpl>(*other.foFldIntrpl))
                      : std::nullopt)
      , foFldOrig(other.foFldOrig.has_value()
                    ? std::make_optional(detail::FieldBase<T, EFieldMode::Orig>(*other.foFldOrig))
                    : std::nullopt)
      , fPrimVertexField(FieldRegion<I>(other.fPrimVertexField))
      , fFieldType(other.fFieldType)
      , fFieldMode(other.fFieldMode)
    {
    }
 
    ~Field() = default;
 
    Field& operator=(const Field& other)
    {
      if (this != &other) {
        fPrimVertexField = other.fPrimVertexField;
        foFldIntrpl      = other.foFldIntrpl;
        foFldOrig        = other.foFldOrig;
        fFieldType       = other.fFieldType;
        fFieldMode       = other.fFieldMode;
      }
      return *this;
    }
 
    FieldValue<T> GetFieldValue(int sliceID, const T& x, const T& y) const
    {
      return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->GetFieldValue(sliceID, x, y)
                                              : foFldOrig->GetFieldValue(sliceID, x, y);
    }
 
    EFieldType GetFieldType() const { return this->fFieldType; }
 
    const FieldRegion<T>& GetPrimVertexField() const { return fPrimVertexField; }
 
    FieldRegion<T> GetFieldRegion(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
                                  const FieldValue<T>& b2, const T& z2) const
    {
      return (fFieldMode == EFieldMode::Intrpl ? FieldRegion<T>(b0, z0, b1, z1, b2, z2)
                                               : FieldRegion<T>(fFieldType, foFldOrig->GetFieldFunction()));
    }
 
    void RemoveSlice(int iLayer);
 
    std::string ToString(int indentLevel, int verbose) const;
 
   private:
    Field() = default;
 
    template<class Archive>
    void load(Archive& ar, const unsigned int /*version*/)
    {
      auto field = detail::FieldBase<T, EFieldMode::Intrpl>{};
      ar >> field;
      foFldIntrpl = std::make_optional(field);
      ar >> fPrimVertexField;
      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 << fPrimVertexField;
        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::FieldBase<T, EFieldMode::Intrpl>> foFldIntrpl{std::nullopt};  
    std::optional<detail::FieldBase<T, EFieldMode::Orig>> foFldOrig{std::nullopt};      
    FieldRegion<T> fPrimVertexField;            
    EFieldType fFieldType{EFieldType::Normal};  
    EFieldMode fFieldMode;                      
  };
 
 
  class FieldFactory {
    struct SliceRef {
      double fHalfSizeX{defs::Undef<double>};  
      double fHalfSizeY{defs::Undef<double>};  
      double fRefZ{defs::Undef<double>};       
 
      SliceRef(double halfX, double halfY, double refZ) : fHalfSizeX(halfX), fHalfSizeY(halfY), fRefZ(refZ) {}
 
      bool operator<(const SliceRef& r) const { return fRefZ < r.fRefZ; }
    };
 
   public:
    FieldFactory() = default;
 
    FieldFactory(const FieldFactory&) = default;
 
    ~FieldFactory() = default;
 
    FieldFactory& operator=(const FieldFactory&) = default;
 
    void AddSliceReference(double halfSizeX, double halfSizeY, double refZ);
 
    EFieldMode GetFieldMode() const { return fFieldMode; }
 
    EFieldType GetFieldType() const { return fFieldType; }
 
    template<typename T>
    Field<T> MakeField() const;
 
    void Reset() { *this = FieldFactory(); }
 
    void ResetSliceReferences() { fSliceReferences.clear(); }
 
    void SetFieldFunction(const FieldFn_t& fieldFn, EFieldType fldType)
    {
      fFieldFn   = fieldFn;
      fFieldType = fldType;
    }
 
    void SetFieldMode(EFieldMode fldMode) { fFieldMode = fldMode; }
 
    void SetStep(double step = 2.5) { fTargetStep = step; }
 
    void SetTarget(double x, double y, double z) { fTarget = {x, y, z}; }
 
   private:
    std::set<SliceRef> fSliceReferences;      
    FieldFn_t fFieldFn{defs::ZeroFieldFn};    
    double fTargetStep{2.5};                  
    EFieldType fFieldType{EFieldType::Null};  
    EFieldMode fFieldMode{EFieldMode::Intrpl};  
 
    std::array<double, 3> fTarget{{defs::Undef<double>, defs::Undef<double>, defs::Undef<double>}};
  };
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<typename T>
  Field<T> FieldFactory::MakeField() const
  {
    auto field = Field<T>(fFieldMode, fFieldType);
 
    // Check initialization
    if (std::any_of(fTarget.begin(), fTarget.end(), [](double x) { return !utils::IsFinite(x); })) {
      throw std::logic_error("FieldFactory::MakeField: target is undefined");
    }
    if (!fSliceReferences.size()) {  // TODO: Remove requirement of slice references
      throw std::logic_error("FieldFactory::MakeField: no slice references were provided");
    }
    if (!fFieldFn) {
      throw std::logic_error("FieldFactory::CreateField: no field function is provided");
    }
 
    // Initialize the Field object
    if (fFieldMode == EFieldMode::Orig) {
      field.foFldOrig->SetFieldFunction(fFieldFn);
      field.fPrimVertexField = FieldRegion<T>(fFieldType, fFieldFn);
      for (const auto& sliceRef : fSliceReferences) {
        field.foFldOrig->AddFieldSlice(sliceRef.fRefZ);
      }
    }
    else if (fFieldMode == EFieldMode::Intrpl) {  // A version with the interpolated field
      for (const auto& sliceRef : fSliceReferences) {
        field.foFldIntrpl->AddFieldSlice(
          FieldSlice<T>(fFieldFn, sliceRef.fHalfSizeX, sliceRef.fHalfSizeY, sliceRef.fRefZ));
      }
      {
        // PV field initialization
        constexpr size_t nNodes = 3;
        double z                = fTarget[2];
        std::array<FieldValue<T>, nNodes> fld;
        std::array<T, nNodes> pos;
        for (size_t iNode = 0; iNode < nNodes; ++iNode) {
          pos[iNode] = utils::simd::Cast<double, T>(z);
          fld[iNode] = GlobalField::GetFieldValue<double>(fFieldFn, fTarget[0], fTarget[1], z);
          z += fTargetStep;
        }
        field.fPrimVertexField = FieldRegion<T>(fld[0], pos[0], fld[1], pos[1], fld[2], pos[2]);
        // TODO: Provide alternative method for Orig
      }
    }
    return field;
  }
}  // namespace cbm::algo::kf