From ecc8b586c5431bbdfc36377b90252ad232adcbeb Mon Sep 17 00:00:00 2001 From: Saveliy Borisov Date: Wed, 23 Sep 2026 23:00:40 +0300 Subject: [PATCH 1/3] Add CIC interpolation sample mode --- README.md | 9 +++- perf-tests/sample/sample.cpp | 91 +++++++++++++++++++++++++++++++++++- 2 files changed, 97 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index ab7da9e..4cd8411 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,14 @@ 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 deposits a point current between grid nodes using triangular weights +from its four neighbouring cells in the x-y plane. Enable it with: + +``` +FDTD_SAMPLE_MODE=interpolation ./sample +``` + # Visualization ![](https://github.com/Amazingkivas/FDTD_Method/blob/main/animations/animation_Ez.gif) - diff --git a/perf-tests/sample/sample.cpp b/perf-tests/sample/sample.cpp index c5b6afd..03afcab 100644 --- a/perf-tests/sample/sample.cpp +++ b/perf-tests/sample/sample.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include "test_FDTD.h" #include "FDTD_PML.h" @@ -147,18 +149,103 @@ void spherical_wave(int n, int it, std::string base_path = "") { #endif //__PML_TEST__ } +// Deposits a point current located between grid nodes with cloud-in-cell (CIC) +// interpolation. The triangular CIC weights distribute the source over the +// four adjacent cells in the x-y plane and sum to one. +void interpolated_point_source(int n, int it) { + if (n < 2) { + throw std::invalid_argument("CIC interpolation requires at least two cells per axis"); + } + + CurrentParameters cur_param { + 8, + 4, + 0.2, + }; + + double d = FDTD_const::C; + double boundary = static_cast(n) / 2.0 * d; + Parameters params { + n, n, n, + -boundary, boundary, + -boundary, boundary, + -boundary, boundary, + d, d, d + }; + + FDTD_openmp::FDTD method(params, cur_param.dt); + int source_time = std::min( + static_cast(static_cast(cur_param.period) / cur_param.dt), it); + + // An off-node position makes all four CIC weights non-zero. + double source_x = params.ax + + (static_cast(n - 2) / 2.0 + 0.35) * params.dx; + double source_y = params.ay + + (static_cast(n - 2) / 2.0 + 0.65) * params.dy; + double grid_x = (source_x - params.ax) / params.dx; + double grid_y = (source_y - params.ay) / params.dy; + int i0 = static_cast(std::floor(grid_x)); + int j0 = static_cast(std::floor(grid_y)); + double wx1 = grid_x - static_cast(i0); + double wy1 = grid_y - static_cast(j0); + double wx0 = 1.0 - wx1; + double wy0 = 1.0 - wy1; + int k = params.Nk / 2; + + auto index = [¶ms, k](int i, int j) { + return i + j * params.Ni + k * params.Ni * params.Nj; + }; + + std::cout << "Interpolation mode: CIC point source with weights " + << wx0 * wy0 << ", " << wx1 * wy0 << ", " + << wx0 * wy1 << ", " << wx1 * wy1 << std::endl; + + auto start = std::chrono::high_resolution_clock::now(); + for (int t = 0; t < it; ++t) { + method.zeroed_currents(); + if (t < source_time) { + double value = std::sin(2.0 * FDTD_const::PI * + static_cast(t + 1) * cur_param.dt / + static_cast(cur_param.period)); + Field& current = method.get_field(Component::JX); + current[index(i0, j0)] += value * wx0 * wy0; + current[index(i0 + 1, j0)] += value * wx1 * wy0; + current[index(i0, j0 + 1)] += value * wx0 * wy1; + current[index(i0 + 1, j0 + 1)] += value * wx1 * wy1; + } + method.update_fields(); + } + auto end = std::chrono::high_resolution_clock::now(); + + std::chrono::duration elapsed = end - start; + std::cout << "Execution time (CIC interpolation): " << elapsed.count() << " s" + << std::endl; +} + int main(int argc, char* argv[]) { std::ifstream source_fin; std::vector 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_point_source(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_point_source(N, Iterations); + } else { + spherical_wave(N, Iterations); + } } else { std::cout << "ERROR: Incorrect number of parameters" << std::endl; From 784c47b5975ec193b49df3664585e9505b4707a1 Mon Sep 17 00:00:00 2001 From: Saveliy Borisov Date: Fri, 25 Sep 2026 01:08:26 +0300 Subject: [PATCH 2/3] Implement Yee-grid CIC field interpolation --- README.md | 5 +- include/FDTD/FDTD.h | 3 ++ perf-tests/sample/sample.cpp | 88 +++++++++++++----------------------- src/FDTD/FDTD.cpp | 76 +++++++++++++++++++++++++++++++ 4 files changed, 114 insertions(+), 58 deletions(-) diff --git a/README.md b/README.md index 4cd8411..e45616f 100644 --- a/README.md +++ b/README.md @@ -35,8 +35,9 @@ This repository contains a C++ project with the main implementation of the metho ``` The OpenMP `sample` also includes an optional cloud-in-cell (CIC) interpolation -example. It deposits a point current between grid nodes using triangular weights -from its four neighbouring cells in the x-y plane. Enable it with: +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 diff --git a/include/FDTD/FDTD.h b/include/FDTD/FDTD.h index cda84dd..70423f6 100644 --- a/include/FDTD/FDTD.h +++ b/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(); }; diff --git a/perf-tests/sample/sample.cpp b/perf-tests/sample/sample.cpp index 03afcab..7425682 100644 --- a/perf-tests/sample/sample.cpp +++ b/perf-tests/sample/sample.cpp @@ -149,20 +149,14 @@ void spherical_wave(int n, int it, std::string base_path = "") { #endif //__PML_TEST__ } -// Deposits a point current located between grid nodes with cloud-in-cell (CIC) -// interpolation. The triangular CIC weights distribute the source over the -// four adjacent cells in the x-y plane and sum to one. -void interpolated_point_source(int n, int it) { +// 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"); } - CurrentParameters cur_param { - 8, - 4, - 0.2, - }; - double d = FDTD_const::C; double boundary = static_cast(n) / 2.0 * d; Parameters params { @@ -173,53 +167,35 @@ void interpolated_point_source(int n, int it) { d, d, d }; - FDTD_openmp::FDTD method(params, cur_param.dt); - int source_time = std::min( - static_cast(static_cast(cur_param.period) / cur_param.dt), it); - - // An off-node position makes all four CIC weights non-zero. - double source_x = params.ax + - (static_cast(n - 2) / 2.0 + 0.35) * params.dx; - double source_y = params.ay + - (static_cast(n - 2) / 2.0 + 0.65) * params.dy; - double grid_x = (source_x - params.ax) / params.dx; - double grid_y = (source_y - params.ay) / params.dy; - int i0 = static_cast(std::floor(grid_x)); - int j0 = static_cast(std::floor(grid_y)); - double wx1 = grid_x - static_cast(i0); - double wy1 = grid_y - static_cast(j0); - double wx0 = 1.0 - wx1; - double wy0 = 1.0 - wy1; - int k = params.Nk / 2; - - auto index = [¶ms, k](int i, int j) { - return i + j * params.Ni + k * params.Ni * params.Nj; + FDTD_openmp::FDTD method(params, 0.2); + const double source_cell = static_cast(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 = [¶ms](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; + } + } + } }; - std::cout << "Interpolation mode: CIC point source with weights " - << wx0 * wy0 << ", " << wx1 * wy0 << ", " - << wx0 * wy1 << ", " << wx1 * wy1 << std::endl; - - auto start = std::chrono::high_resolution_clock::now(); - for (int t = 0; t < it; ++t) { - method.zeroed_currents(); - if (t < source_time) { - double value = std::sin(2.0 * FDTD_const::PI * - static_cast(t + 1) * cur_param.dt / - static_cast(cur_param.period)); - Field& current = method.get_field(Component::JX); - current[index(i0, j0)] += value * wx0 * wy0; - current[index(i0 + 1, j0)] += value * wx1 * wy0; - current[index(i0, j0 + 1)] += value * wx0 * wy1; - current[index(i0 + 1, j0 + 1)] += value * wx1 * wy1; - } - method.update_fields(); - } - auto end = std::chrono::high_resolution_clock::now(); + 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::chrono::duration elapsed = end - start; - std::cout << "Execution time (CIC interpolation): " << elapsed.count() << " s" - << std::endl; + 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[]) { @@ -233,7 +209,7 @@ int main(int argc, char* argv[]) { int N = 32; int Iterations = 100; if (use_interpolation) { - interpolated_point_source(N, Iterations); + interpolated_field_example(N, Iterations); } else { spherical_wave(N, Iterations, "../../"); } @@ -242,7 +218,7 @@ int main(int argc, char* argv[]) { int N = std::atoi(arguments[1]); int Iterations = std::atoi(arguments[2]); if (use_interpolation) { - interpolated_point_source(N, Iterations); + interpolated_field_example(N, Iterations); } else { spherical_wave(N, Iterations); } diff --git a/src/FDTD/FDTD.cpp b/src/FDTD/FDTD.cpp index 61dde61..8c2bedb 100644 --- a/src/FDTD/FDTD.cpp +++ b/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(std::floor(grid_x)); + const int j0 = static_cast(std::floor(grid_y)); + const int k0 = static_cast(std::floor(grid_z)); + const FP wx1 = grid_x - static_cast(i0); + const FP wy1 = grid_y - static_cast(j0); + const FP wz1 = grid_z - static_cast(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(); From 735c09d0206ea46ef14db62f6e61f0fc56fce47d Mon Sep 17 00:00:00 2001 From: Saveliy Borisov Date: Fri, 25 Sep 2026 01:09:14 +0300 Subject: [PATCH 3/3] Run CIC interpolation sample in CI --- .github/workflows/main.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index e00cfcb..1384471 100644 --- a/.github/workflows/main.yml +++ b/.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: |