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: | diff --git a/README.md b/README.md index ab7da9e..e45616f 100644 --- a/README.md +++ b/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) - 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 c5b6afd..7425682 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,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(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(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; + } + } + } + }; + + 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 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; 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();