Add CIC interpolation sample mode (OpenMP sample)

.github/workflows/main.yml

@@ -45,10 +45,10 @@ jobs:
           find ${PWD}/build -name kokkos_sample
           ${PWD}/bin/kokkos_sample 512 25
 
-      - name: Run sample
+      - name: Run CIC interpolation sample
         shell: bash
         run: |
-          ${PWD}/bin/sample 512 25
+          FDTD_SAMPLE_MODE=interpolation ${PWD}/bin/sample 512 25
 
       - name: Install Intel OneAPI (ifort)
         run: |

README.md

@@ -34,7 +34,15 @@ This repository contains a C++ project with the main implementation of the metho
   ./kokkos_sample
   ```
 
+The OpenMP `sample` also includes an optional cloud-in-cell (CIC) interpolation
+example. It reads a Yee-grid field at an arbitrary physical position with
+component-specific spatial offsets and trilinear weights from eight neighbouring
+grid values. Enable it with:
+
+```
+FDTD_SAMPLE_MODE=interpolation ./sample
+```
+
 # Visualization
 
 ![](https://github.com/Amazingkivas/FDTD_Method/blob/main/animations/animation_Ez.gif)
-

include/FDTD/FDTD.h

@@ -36,6 +36,9 @@ class FDTD {
     FDTD(Parameters _parameters, double _dt);
 
     Field& get_field(Component this_field);
+    // Returns a field component at physical coordinates using trilinear CIC
+    // interpolation. Component-specific Yee-grid spatial offsets are applied.
+    FP get_field_CIC(Component this_field, FP x, FP y, FP z) const;
     virtual void update_fields();
     void zeroed_currents();
 };

perf-tests/sample/sample.cpp

@@ -4,6 +4,8 @@
 #include <filesystem>
 #include <cmath>
 #include <chrono>
+#include <cstdlib>
+#include <string>
 
 #include "test_FDTD.h"
 #include "FDTD_PML.h"
@@ -147,18 +149,79 @@ void spherical_wave(int n, int it, std::string base_path = "") {
 #endif //__PML_TEST__
 }
 
+// Demonstrates trilinear CIC interpolation of a Yee-grid field at an
+// arbitrary physical position. Each component is sampled with its own spatial
+// offset and therefore uses eight neighbouring grid values.
+void interpolated_field_example(int n, int) {
+    if (n < 2) {
+        throw std::invalid_argument("CIC interpolation requires at least two cells per axis");
+    }
+
+    double d = FDTD_const::C;
+    double boundary = static_cast<double>(n) / 2.0 * d;
+    Parameters params {
+        n, n, n,
+        -boundary, boundary,
+        -boundary, boundary,
+        -boundary, boundary,
+        d, d, d
+    };
+
+    FDTD_openmp::FDTD method(params, 0.2);
+    const double source_cell = static_cast<double>(n - 1) / 2.0;
+    const double x = params.ax + (source_cell + 0.15) * params.dx;
+    const double y = params.ay + (source_cell + 0.20) * params.dy;
+    const double z = params.az + (source_cell + 0.25) * params.dz;
+    const double expected = x + 2.0 * y + 3.0 * z;
+
+    const auto fill_linear_field = [&params](Field& field, double sx, double sy, double sz) {
+        for (int k = 0; k < params.Nk; ++k) {
+            for (int j = 0; j < params.Nj; ++j) {
+                for (int i = 0; i < params.Ni; ++i) {
+                    const int index = i + j * params.Ni + k * params.Ni * params.Nj;
+                    const double field_x = params.ax + (i + sx) * params.dx;
+                    const double field_y = params.ay + (j + sy) * params.dy;
+                    const double field_z = params.az + (k + sz) * params.dz;
+                    field[index] = field_x + 2.0 * field_y + 3.0 * field_z;
+                }
+            }
+        }
+    };
+
+    fill_linear_field(method.get_field(Component::EX), 0.0, 0.5, 0.5);
+    fill_linear_field(method.get_field(Component::BX), 0.5, 0.0, 0.0);
+
+    std::cout << "CIC trilinear interpolation at (" << x << ", " << y << ", " << z
+              << "):\n  Ex = " << method.get_field_CIC(Component::EX, x, y, z)
+              << " (expected " << expected << ")\n  Bx = "
+              << method.get_field_CIC(Component::BX, x, y, z)
+              << " (expected " << expected << ")" << std::endl;
+}
+
 int main(int argc, char* argv[]) {
     std::ifstream source_fin;
     std::vector<char*> arguments(argv, argv + argc);
+    const char* sample_mode = std::getenv("FDTD_SAMPLE_MODE");
+    bool use_interpolation = sample_mode != nullptr &&
+        std::string(sample_mode) == "interpolation";
+
     if (argc == 1) {
         int N = 32;
         int Iterations = 100;
-        spherical_wave(N, Iterations, "../../");
+        if (use_interpolation) {
+            interpolated_field_example(N, Iterations);
+        } else {
+            spherical_wave(N, Iterations, "../../");
+        }
     }
     else if (argc == 3) {
         int N = std::atoi(arguments[1]);
         int Iterations = std::atoi(arguments[2]);
-        spherical_wave(N, Iterations);
+        if (use_interpolation) {
+            interpolated_field_example(N, Iterations);
+        } else {
+            spherical_wave(N, Iterations);
+        }
     }
     else {
         std::cout << "ERROR: Incorrect number of parameters" << std::endl;

src/FDTD/FDTD.cpp

@@ -1,5 +1,37 @@
 #include "FDTD.h"
 
+namespace {
+
+struct SpatialShift {
+    FP x;
+    FP y;
+    FP z;
+};
+
+SpatialShift yee_shift(FDTD_enums::Component component) {
+    using FDTD_enums::Component;
+
+    switch (component) {
+        case Component::EX:
+        case Component::JX: return {0.0, 0.5, 0.5};
+        case Component::EY:
+        case Component::JY: return {0.5, 0.0, 0.5};
+        case Component::EZ:
+        case Component::JZ: return {0.5, 0.5, 0.0};
+        case Component::BX: return {0.5, 0.0, 0.0};
+        case Component::BY: return {0.0, 0.5, 0.0};
+        case Component::BZ: return {0.0, 0.0, 0.5};
+        default: throw std::logic_error("ERROR: Invalid field component");
+    }
+}
+
+int periodic_index(int index, int size) {
+    int wrapped = index % size;
+    return wrapped < 0 ? wrapped + size : wrapped;
+}
+
+}
+
 FDTD_openmp::FDTD::FDTD(Parameters _parameters, FP _dt)
     : parameters(_parameters), dt(_dt) {
     if (parameters.Ni <= 0 || parameters.Nj <= 0 || parameters.Nk <= 0 || dt <= 0) {
@@ -150,6 +182,50 @@ FDTD_openmp::Field& FDTD_openmp::FDTD::get_field(Component this_field) {
     }
 }
 
+FP FDTD_openmp::FDTD::get_field_CIC(Component this_field, FP x, FP y, FP z) const {
+    const Field* field = nullptr;
+    switch (this_field) {
+        case Component::JX: field = &Jx; break;
+        case Component::JY: field = &Jy; break;
+        case Component::JZ: field = &Jz; break;
+        case Component::EX: field = &Ex; break;
+        case Component::EY: field = &Ey; break;
+        case Component::EZ: field = &Ez; break;
+        case Component::BX: field = &Bx; break;
+        case Component::BY: field = &By; break;
+        case Component::BZ: field = &Bz; break;
+        default: throw std::logic_error("ERROR: Invalid field component");
+    }
+
+    const SpatialShift shift = yee_shift(this_field);
+    const FP grid_x = (x - parameters.ax) / dx - shift.x;
+    const FP grid_y = (y - parameters.ay) / dy - shift.y;
+    const FP grid_z = (z - parameters.az) / dz - shift.z;
+    const int i0 = static_cast<int>(std::floor(grid_x));
+    const int j0 = static_cast<int>(std::floor(grid_y));
+    const int k0 = static_cast<int>(std::floor(grid_z));
+    const FP wx1 = grid_x - static_cast<FP>(i0);
+    const FP wy1 = grid_y - static_cast<FP>(j0);
+    const FP wz1 = grid_z - static_cast<FP>(k0);
+    const FP wx[2] = {1.0 - wx1, wx1};
+    const FP wy[2] = {1.0 - wy1, wy1};
+    const FP wz[2] = {1.0 - wz1, wz1};
+
+    FP value = 0.0;
+    for (int dk = 0; dk < 2; ++dk) {
+        const int k = periodic_index(k0 + dk, Nk);
+        for (int dj = 0; dj < 2; ++dj) {
+            const int j = periodic_index(j0 + dj, Nj);
+            for (int di = 0; di < 2; ++di) {
+                const int i = periodic_index(i0 + di, Ni);
+                const int index = i + j * Ni + k * Ni * Nj;
+                value += wx[di] * wy[dj] * wz[dk] * (*field)[index];
+            }
+        }
+    }
+    return value;
+}
+
 void FDTD_openmp::FDTD::update_fields() {
     update_B();
     update_E();