Skip to content

File latticelaplacian.h

File List > algebra > spatialderivatives > latticelaplacian.h

Go to the documentation of this file

#ifndef TEMPLAT_LATTICE_ALGEBRA_SPATIALDERIVATIVES_LAPLACIANLOCAL_H
#define TEMPLAT_LATTICE_ALGEBRA_SPATIALDERIVATIVES_LAPLACIANLOCAL_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/constants/zerotype.h"
#include "TempLat/util/constexpr_for.h"
#include "TempLat/lattice/algebra/operators/unaryoperator.h"
#include "TempLat/lattice/algebra/helpers/getderiv.h"
#include "TempLat/lattice/algebra/helpers/getgetreturntype.h"
#include "TempLat/lattice/algebra/helpers/getfloattype.h"
#include "TempLat/lattice/algebra/operators/operators.h"
#include "TempLat/util/tuple_tools.h"
#include "TempLat/lattice/algebra/helpers/haseval.h"
#include "TempLat/lattice/algebra/helpers/getndim.h"

#include "TempLat/parallel/device.h"

namespace TempLat
{
  template <typename R> class LatticeLaplacian : public UnaryOperator<R>
  {
  private:
    using UnaryOperator<R>::mR;

  public:
    static constexpr size_t NDim = GetNDim::get<R>();
    using GetReturnType = typename GetGetReturnType<R>::type;
    using FloatType = typename GetFloatType<GetReturnType>::type;

    LatticeLaplacian(R pR) : UnaryOperator<R>(pR), dx2(pow(GetDx::getDx(pR), 2)) {}

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

    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 {
        auto result = (-2 * static_cast<device::Idx>(NDim) * DoEval::eval(mR, idx...));
        constexpr_for<0, NDim>([&](const auto _d) {
          constexpr size_t d = decltype(_d)::value;
          device::apply([&](const auto &...shifted_idx) { result += DoEval::eval(mR, shifted_idx...); },
                        tuple_add_to_nth<d, 1>(device::tie(idx...)));
          device::apply([&](const auto &...shifted_idx) { result += DoEval::eval(mR, shifted_idx...); },
                        tuple_add_to_nth<d, -1>(device::tie(idx...)));
        });
        return result / dx2;
      }
    }

    virtual std::string operatorString() const override { return "Laplacian"; }

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

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

  template <size_t NDim_ = 0, typename R>
    requires(HasEvalMethod<R> && GetNDim::get<std::decay_t<R>>() > 0)
  auto LatLapl(R pR)
  {
    static_assert(NDim_ == 0 || NDim_ == GetNDim::get<R>(),
                  "Explicit NDim does not match the NDim deduced from expression type R.");
    return LatticeLaplacian<R>(pR);
  }

  template <size_t NDim_ = 0, typename R>
    requires(!HasEvalMethod<R> || GetNDim::get<std::decay_t<R>>() == 0)
  constexpr auto LatLapl(R)
  {
    return ZeroType{};
  }
} // namespace TempLat

#endif