Skip to content

File latticeforwardgradient.h

File List > algebra > spatialderivatives > latticeforwardgradient.h

Go to the documentation of this file

#ifndef TEMPLAT_LATTICE_ALGEBRA_SPATIALDERIVATIVES_FORWARDGRADIENTLOCAL_H
#define TEMPLAT_LATTICE_ALGEBRA_SPATIALDERIVATIVES_FORWARDGRADIENTLOCAL_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, Franz R. Sattler, Year: 2025

#include "TempLat/lattice/algebra/operators/unaryoperator.h"
#include "TempLat/lattice/algebra/helpers/getstring.h"
#include "TempLat/lattice/algebra/helpers/getgetreturntype.h"
#include "TempLat/lattice/algebra/helpers/getfloattype.h"
#include "TempLat/lattice/algebra/helpers/gettoolbox.h"
#include "TempLat/lattice/algebra/helpers/getvectorcomponent.h"
#include "TempLat/lattice/algebra/helpers/getderiv.h"
#include "TempLat/lattice/algebra/helpers/getkir.h"
#include "TempLat/lattice/algebra/helpers/getdx.h"
#include "TempLat/lattice/algebra/helpers/doeval.h"
#include "TempLat/lattice/algebra/listoperators/vectordotter.h"
#include "TempLat/util/tuple_tools.h"
#include "TempLat/lattice/algebra/spatialderivatives/forwdiff.h"
#include "TempLat/lattice/algebra/helpers/getndim.h"

namespace TempLat
{

  template <typename R> class LatticeForwardGradient : public UnaryOperator<R>
  {
  private:
    using UnaryOperator<R>::mR;

  public:
    // Put public methods here. These should change very little over time.
    using GetReturnType = typename GetGetReturnType<R>::type;
    using FloatType = typename GetFloatType<GetReturnType>::type;

    static constexpr size_t NDim = GetNDim::get<R>();

    LatticeForwardGradient(const R &pR) : UnaryOperator<R>(pR), dx(GetDx::getDx(mR)) {}

    template <typename... IDX>
      requires requires(std::decay_t<R> r, IDX... idx) {
        requires IsVariadicIndex<IDX...>;
        DoEval::eval(r, idx...);
      }
    DEVICE_INLINE_FUNCTION auto eval(const IDX &...idx) const
    {
      if constexpr (UnaryOperator<R>::getNDim() == 0)
        return ZeroType();
      else {
        device::array<GetReturnType, NDim> result{};
        constexpr_for<0, NDim>([&](const auto _i) {
          constexpr int i = decltype(_i)::value;

          result[i] = -DoEval::eval(mR, idx...);
          device::apply([&](const auto &...shifted_idx) { result[i] += DoEval::eval(mR, shifted_idx...); },
                        tuple_add_to_nth<i, 1>(device::tie(idx...)));
          result[i] /= dx;
        });
        return result;
      }
    }

    template <int N> auto vectorGet(Tag<N> t) const
    {
      static_assert(N > 0 && N <= NDim, "VectorGet: N must be between 1 and NDim for LatticeForwardGradient");
      return ForwDiff<N, R>(mR);
    }

    template <int N> auto operator()(Tag<N> t) const
    {
      static_assert(N > 0 && N <= NDim, "Operator(): N must be between 1 and NDim for LatticeForwardGradient");
      return vectorGet(t);
    }

    static constexpr size_t getVectorSize() { return NDim; }

    auto norm2() { return dot(*this, *this); }
    auto norm() { return pow(this->norm2(), 0.5); }

    std::string toString() const { return "Grad(" + GetString::get(mR) + ")"; }

    void doWeNeedGhosts() const { mR.confirmGhostsUpToDate(); }

    auto getDx() const { return dx; }
    auto getKIR() const { return GetKIR::getKIR(mR); }

    void confirmSpace(device::Idx i, const LayoutStruct<NDim> &newLayout, const SpaceStateType &spaceType) const
    {
      ConfirmSpace::apply(mR, i, newLayout, spaceType);
    }

    template <int N> device::Idx confirmGhostsUpToDate(Tag<N> i) const { return ConfirmGhosts::apply(mR, i); }

    inline device::memory::host_ptr<MemoryToolBox<NDim>> getToolBox() const { return GetToolBox::get(mR); }

    template <typename S> inline auto d(const S &other) { return LatForwardGrad(GetDeriv::get(mR, other)); }

  private:
    /* Put all member variables and private methods here. These may change arbitrarily. */
    const FloatType dx;
  };

  template <size_t NDim_ = 0, typename R> auto LatForwardGrad(R pR)
  {
    static_assert(NDim_ == 0 || NDim_ == GetNDim::get<R>(),
                  "Explicit NDim does not match the NDim deduced from expression type R.");
    return LatticeForwardGradient<R>(pR);
  }

} // namespace TempLat

#endif