File su2liealgebrafield.h

File List > algebra > su2algebra > su2liealgebrafield.h

Go to the documentation of this file

#ifndef COSMOINTERFACE_SU2ALGEBRA_SU2LIEALGEBRAFIELD_H
#define COSMOINTERFACE_SU2ALGEBRA_SU2LIEALGEBRAFIELD_H

/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
   Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Adrien Florio, Year: 2021

#include "TempLat/lattice/algebra/su2algebra/su2field.h"

namespace TempLat
{
  template <typename T, size_t _NDim = 0> class SU2LieAlgebraField
  {
  public:
    // Put public methods here. These should change very little over time.
    static_assert(_NDim != 0, "NDim template parameter is required. Use e.g. SU2LieAlgebraField<double, 3>.");

    static constexpr size_t NDim = _NDim;

    SU2LieAlgebraField(Field<T, NDim> f1, Field<T, NDim> f2, Field<T, NDim> f3)
        : fs{{f1, f2, f3}}, mName("NoName"), mLayout(fs[0].getToolBox()->mLayouts.getConfigSpaceLayout())
    {
    }

    SU2LieAlgebraField(std::string name, device::memory::host_ptr<MemoryToolBox<NDim>> toolBox,
                       LatticeParameters<T> pLatPar)
        : fs{{
              Field<T, NDim>(name + "_1", toolBox, pLatPar), //
              Field<T, NDim>(name + "_2", toolBox, pLatPar), //
              Field<T, NDim>(name + "_3", toolBox, pLatPar)  //
          }},
          mName(name), mLayout(toolBox->mLayouts.getConfigSpaceLayout())
    {
    }

    template <int N> auto SU2Get(Tag<N> t) const
    {
      static_assert(N >= 0 && N <= 3, "SU2Get: N must be between 0 and 3 for SU2LieAlgebraField");
      return operator()(t);
    }

    template <int N> auto SU2LieAlgebraGet(Tag<N> t) const { return 2 * SU2Get(t); }

    ZeroType operator()(Tag<0> t) const { return {}; }

    template <int M> auto &operator()(Tag<M> t)
    {
      static_assert(M >= 1 && M <= 3, "Operator(): M must be between 0 and 3 for SU2LieAlgebraField");
      return fs[M - 1];
    }
    template <int M> const auto &operator()(Tag<M> t) const
    {
      static_assert(M >= 1 && M <= 3, "Operator(): M must be between 0 and 3 for SU2LieAlgebraField");
      return fs[M - 1];
    }

    template <typename R> void operator=(R &&r)
    {
      // Confirm config-space / ghost requirements by walking the WHOLE SU(2) expression r (linear),
      // instead of building r.SU2Get(k) per component (the exponential per-component expansion built
      // only for confirmation, then discarded - values use the fused DoEval::eval below). One call per
      // target component is kept (onBeforeAssignment also confirms the target field fs[k]).
      fs[0].onBeforeAssignment(r);
      fs[1].onBeforeAssignment(r);
      fs[2].onBeforeAssignment(r);

      PreGet::apply(r);

      const auto view1 = fs[0].getView();
      const auto view2 = fs[1].getView();
      const auto view3 = fs[2].getView();

      auto functor = DEVICE_CLASS_LAMBDA(const device::IdxArray<NDim> &idx)
      {
        device::apply(
            [&](const auto &...args) {
              auto result = DoEval::eval(r, args...);
              view1(args...) = result[1];
              view2(args...) = result[2];
              view3(args...) = result[3];
            },
            idx);
      };
      device::iteration::foreach ("SU2AlgebraConfigViewAssign", mLayout, functor);

      PostGet::apply(r);

      fs[0].setGhostsAreStale();
      fs[1].setGhostsAreStale();
      fs[2].setGhostsAreStale();
    }

    template <typename... IDX>
      requires requires(Field<T, NDim> f, IDX... idx) {
        requires IsVariadicIndex<IDX...>;
        DoEval::eval(f, idx...);
      }
    DEVICE_INLINE_FUNCTION auto eval(const IDX &...idx) const
    {
      device::array<T, 4> result;
      result[0] = T(0);
      result[1] = fs[0].eval(idx...);
      result[2] = fs[1].eval(idx...);
      result[3] = fs[2].eval(idx...);
      return result;
    }

    std::string toString() const { return *mName; }

    auto getDx() const { return GetDx::getDx(fs[0]); }
    auto getKIR() const { return GetKIR::getKIR(fs[0]); }

    inline auto getToolBox() const { return GetToolBox::get(fs[0]); }

    // Space/ghost confirmation forwarded to the 3 algebra component fields (see operator=).
    void doWeNeedGhosts() const
    {
      GhostsHunter::apply(fs[0]);
      GhostsHunter::apply(fs[1]);
      GhostsHunter::apply(fs[2]);
    }
    device::Idx confirmGhostsUpToDate() const
    {
      MemoryManager<T, NDim> *mgrs[] = {fs[0].getMemoryManager().get(), fs[1].getMemoryManager().get(),
                                        fs[2].getMemoryManager().get()};
      return MemoryManager<T, NDim>::confirmGhostsUpToDateBatch(mgrs);
    }
    void confirmSpace(const LayoutStruct<NDim> &newLayout, const SpaceStateType &spaceType) const
    {
      ConfirmSpace::apply(fs[0], newLayout, spaceType);
      ConfirmSpace::apply(fs[1], newLayout, spaceType);
      ConfirmSpace::apply(fs[2], newLayout, spaceType);
    }

    inline void updateGhosts()
    {
      MemoryManager<T, NDim> *mgrs[] = {fs[0].getMemoryManager().get(), fs[1].getMemoryManager().get(),
                                        fs[2].getMemoryManager().get()};
      MemoryManager<T, NDim>::updateGhostsBatch(mgrs);
    }

    using Getter = SU2Getter;
    static constexpr size_t SHIFTIND = 0;
    static constexpr size_t size = 4;
    static constexpr size_t numberToSkipAsTuple = 1;

  private:
    device::array<Field<T, NDim>, 3> fs;
    const device::memory::host_string mName;
    LayoutStruct<NDim> mLayout;
  };
} // namespace TempLat

#endif