operations.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd
// SPDX-License-Identifier: AGPL-3.0-or-later
 
#pragma once
 
#include <type_traits>
 
#include "openpfc/core/world.hpp"
#include "openpfc/fft.hpp"
#include "openpfc/model.hpp"
#include "openpfc/types.hpp"
// Local iteration implemented inline to work with HeFFTe inbox type
 
namespace pfc {
namespace field {
 
template <typename Fn>
inline void apply(RealField &field, const World &world, const FFT &fft, Fn &&fn) {
  const auto inbox = pfc::fft::get_inbox(fft);
  // Safety: ensure field size matches inbox voxel count
  const auto nx = inbox.size[0];
  const auto ny = inbox.size[1];
  const auto nz = inbox.size[2];
  const auto expected = static_cast<size_t>(nx) * ny * nz;
  if (field.size() != expected) {
    throw std::invalid_argument(
        "field::apply: field size does not match FFT inbox size");
  }
 
  size_t linear_idx = 0;
  for (int k = inbox.low[2]; k <= inbox.high[2]; ++k) {
    for (int j = inbox.low[1]; j <= inbox.high[1]; ++j) {
      for (int i = inbox.low[0]; i <= inbox.high[0]; ++i) {
        const pfc::Int3 idx{i, j, k};
        const auto x = pfc::world::to_coords(world, idx);
        field[linear_idx++] = static_cast<double>(fn(x));
      }
    }
  }
}
 
template <typename Fn>
inline void apply_with_time(RealField &field, const World &world, const FFT &fft,
                            double t, Fn &&fn) {
  const auto inbox = pfc::fft::get_inbox(fft);
  const auto nx = inbox.size[0];
  const auto ny = inbox.size[1];
  const auto nz = inbox.size[2];
  const auto expected = static_cast<size_t>(nx) * ny * nz;
  if (field.size() != expected) {
    throw std::invalid_argument(
        "field::apply_with_time: field size does not match FFT inbox size");
  }
 
  size_t linear_idx = 0;
  for (int k = inbox.low[2]; k <= inbox.high[2]; ++k) {
    for (int j = inbox.low[1]; j <= inbox.high[1]; ++j) {
      for (int i = inbox.low[0]; i <= inbox.high[0]; ++i) {
        const pfc::Int3 idx{i, j, k};
        const auto x = pfc::world::to_coords(world, idx);
        field[linear_idx++] = static_cast<double>(fn(x, t));
      }
    }
  }
}
 
template <typename Fn>
inline void apply_inplace(RealField &field, const World &world, const FFT &fft,
                          Fn &&fn) {
  const auto inbox = pfc::fft::get_inbox(fft);
  const auto nx = inbox.size[0];
  const auto ny = inbox.size[1];
  const auto nz = inbox.size[2];
  const auto expected = static_cast<size_t>(nx) * ny * nz;
  if (field.size() != expected) {
    throw std::invalid_argument(
        "field::apply_inplace: field size does not match FFT inbox size");
  }
 
  size_t linear_idx = 0;
  for (int k = inbox.low[2]; k <= inbox.high[2]; ++k) {
    for (int j = inbox.low[1]; j <= inbox.high[1]; ++j) {
      for (int i = inbox.low[0]; i <= inbox.high[0]; ++i) {
        const pfc::Int3 idx{i, j, k};
        const auto x = pfc::world::to_coords(world, idx);
        field[linear_idx] = static_cast<double>(fn(x, field[linear_idx]));
        ++linear_idx;
      }
    }
  }
}
 
template <typename Fn>
inline void apply_inplace_with_time(RealField &field, const World &world,
                                    const FFT &fft, double t, Fn &&fn) {
  const auto inbox = pfc::fft::get_inbox(fft);
  const auto nx = inbox.size[0];
  const auto ny = inbox.size[1];
  const auto nz = inbox.size[2];
  const auto expected = static_cast<size_t>(nx) * ny * nz;
  if (field.size() != expected) {
    throw std::invalid_argument(
        "field::apply_inplace_with_time: field size does not match FFT inbox size");
  }
 
  size_t linear_idx = 0;
  for (int k = inbox.low[2]; k <= inbox.high[2]; ++k) {
    for (int j = inbox.low[1]; j <= inbox.high[1]; ++j) {
      for (int i = inbox.low[0]; i <= inbox.high[0]; ++i) {
        const pfc::Int3 idx{i, j, k};
        const auto x = pfc::world::to_coords(world, idx);
        field[linear_idx] = static_cast<double>(fn(x, field[linear_idx], t));
        ++linear_idx;
      }
    }
  }
}
 
template <typename Fn>
inline void apply_inplace(Model &model, std::string_view field_name, Fn &&fn) {
  auto &f = model.get_real_field(field_name);
  apply_inplace(f, model.get_world(), model.get_fft(), std::forward<Fn>(fn));
}
 
template <typename Fn>
inline void apply_inplace_with_time(Model &model, std::string_view field_name,
                                    double t, Fn &&fn) {
  auto &f = model.get_real_field(field_name);
  apply_inplace_with_time(f, model.get_world(), model.get_fft(), t,
                          std::forward<Fn>(fn));
}
 
template <typename Fn>
inline void apply(Model &model, std::string_view field_name, Fn &&fn) {
  auto &f = model.get_real_field(field_name);
  apply(f, model.get_world(), model.get_fft(), std::forward<Fn>(fn));
}
 
template <typename Fn>
inline void apply_with_time(Model &model, std::string_view field_name, double t,
                            Fn &&fn) {
  auto &f = model.get_real_field(field_name);
  apply_with_time(f, model.get_world(), model.get_fft(), t, std::forward<Fn>(fn));
}
 
} // namespace field
} // namespace pfc