diff --git a/tdms/CMakeLists.txt b/tdms/CMakeLists.txt index 38a946162..c9ac7b582 100644 --- a/tdms/CMakeLists.txt +++ b/tdms/CMakeLists.txt @@ -154,6 +154,8 @@ if (BUILD_TESTING) add_executable(tdms_tests tests/unit/test_fields.cpp + tests/unit/test_BLi_vs_cubic_interpolation.cpp + tests/unit/test_field_interpolation.cpp tests/unit/test_interpolation_determination.cpp tests/unit/test_interpolation_functions.cpp tests/unit/test_numerical_derivative.cpp diff --git a/tdms/include/arrays.h b/tdms/include/arrays.h index ff28a6284..33b37428b 100644 --- a/tdms/include/arrays.h +++ b/tdms/include/arrays.h @@ -25,6 +25,13 @@ class XYZTensor3D { default: throw std::runtime_error("Have no element " + to_string(c)); } } + + /** + * @brief Allocates x, y, and z as (K_total+1) * (J_total+1) * (I_total+1) arrays + * + * @param I_total,J_total,K_total Dimensions of the tensor size to set + */ + void allocate(int I_total, int J_total, int K_total); }; class XYZVectors { @@ -222,7 +229,7 @@ class Tensor3D{ void allocate(int nK, int nJ, int nI){ n_layers = nK, n_cols = nJ, n_rows = nI; - tensor = (T ***)malloc(n_layers * sizeof(T *)); + tensor = (T ***)malloc(n_layers * sizeof(T **)); for(int k=0; k < n_layers; k++){ tensor[k] = (T **)malloc(n_cols * sizeof(T *)); @@ -235,6 +242,13 @@ class Tensor3D{ } }; + /** + * @brief Computes the Frobenius norm of the tensor + * + * fro_norm = \f$\sqrt{ \sum_{i=0}^{I_tot}\sum_{j=0}^{J_tot}\sum_{k=0}^{K_tot} |t[k][j][i]|^2 }\f$ + */ + double frobenius(); + ~Tensor3D(){ if (tensor == nullptr) return; if (is_matlab_initialised){ diff --git a/tdms/include/field.h b/tdms/include/field.h index 772d8e487..2cad7d710 100644 --- a/tdms/include/field.h +++ b/tdms/include/field.h @@ -10,6 +10,7 @@ #include "mat_io.h" #include "simulation_parameters.h" #include "utils.h" +#include "globals.h" /** @@ -22,6 +23,9 @@ * (0, 0, 0) (1, 0, 0) (2, 0, 0) * * has I_tot = 2, J_tot = 1, K_tot = 0. + * + * NOTE: For storage purposes, this means that field values associated to cells are stored _to the left_. + * That is, Grid(0,0,0) is associated to the cell (-1,-1,-1). This is contrary to the way values are associated to cells, where cell (0,0,0) is associated to the field values (0,0,0). */ class Grid{ @@ -111,6 +115,15 @@ class SplitField : public Grid{ * @param eh_vec // TODO */ void initialise_fftw_plan(int n_threads, EHVec &eh_vec); + + /** + * @brief Interpolates a SplitField component to the centre of a Yee cell + * + * @param d SplitField component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return double The interpolated field value + */ + virtual double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0; }; class ElectricSplitField: public SplitField{ @@ -126,6 +139,15 @@ class ElectricSplitField: public SplitField{ */ ElectricSplitField(int I_total, int J_total, int K_total) : SplitField(I_total, J_total, K_total){}; + + /** + * @brief Interpolates a split E-field component to the centre of a Yee cell + * + * @param d Field component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return double The interpolated component value + */ + double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; }; class MagneticSplitField: public SplitField{ @@ -141,6 +163,15 @@ class MagneticSplitField: public SplitField{ */ MagneticSplitField(int I_total, int J_total, int K_total) : SplitField(I_total, J_total, K_total){}; + + /** + * @brief Interpolates a split E-field component to the centre of a Yee cell + * + * @param d Field component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return double The interpolated component value + */ + double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; }; class CurrentDensitySplitField: public SplitField{ @@ -156,6 +187,8 @@ class CurrentDensitySplitField: public SplitField{ */ CurrentDensitySplitField(int I_total, int J_total, int K_total) : SplitField(I_total, J_total, K_total){}; + + double interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override { return 0.; }; }; /** @@ -169,7 +202,7 @@ class Field : public Grid{ std::complex angular_norm = 0.; - // TODO: this is likely better as a set of complex arrays + // TODO: this is likely better as a set of complex arrays - use Tensor3d> XYZTensor3D real; XYZTensor3D imag; @@ -184,10 +217,35 @@ class Field : public Grid{ void normalise_volume(); /** - * Zero all components of the real and imaginary parts of the field + * Default no arguments constructor + */ + Field() = default; + + /** + * Constructor of the field with a defined size in the x, y, z Cartesian + * dimensions + */ + Field(int I_total, int J_total, int K_total); + + /** + * Allocate the memory appropriate for all the 3D tensors associated with + * this split field + */ + void allocate(); + + /** + * Set all the values of all components of the field to zero */ void zero(); + /** + * Allocate and set to zero all components of the field + */ + void allocate_and_zero() { + allocate(); + zero(); + } + /** * Set the phasors for this field, given a split field. Result gives field according to the * exp(-iwt) convention @@ -206,6 +264,14 @@ class Field : public Grid{ std::complex phasor_norm(double f, int n, double omega, double dt, int Nt); virtual double phase(int n, double omega, double dt) = 0; + /** + * @brief Interpolates a Field component to the centre of a Yee cell + * + * @param d Field component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return std::complex The interpolated field value + */ + virtual std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) = 0; /** * Set the values of all components in this field from another, equally sized field @@ -219,12 +285,38 @@ class ElectricField: public Field{ private: double phase(int n, double omega, double dt) override; + +public: + ElectricField() = default; + ElectricField(int I_total, int J_total, int K_total) : Field(I_total, J_total, K_total){}; + + /** + * @brief Interpolates an E-field component to the centre of a Yee cell + * + * @param d Field component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return std::complex The interpolated component value + */ + std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; }; class MagneticField: public Field{ private: double phase(int n, double omega, double dt) override; + +public: + MagneticField() = default; + MagneticField(int I_total, int J_total, int K_total) : Field(I_total, J_total, K_total){}; + + /** + * @brief Interpolates an H-field component to the centre of a Yee cell + * + * @param d Field component to interpolate + * @param i,j,k Index (i,j,k) of the Yee cell to interpolate to the centre of + * @return std::complex The interpolated component value + */ + std::complex interpolate_to_centre_of(AxialDirection d, int i, int j, int k) override; }; /** diff --git a/tdms/include/globals.h b/tdms/include/globals.h index 8158fb421..087572e58 100644 --- a/tdms/include/globals.h +++ b/tdms/include/globals.h @@ -38,6 +38,12 @@ typedef struct complex_vector std::complex Z; } complex_vector; +enum AxialDirection +{ + X = 'x', + Y = 'y', + Z = 'z' +}; // ********************** // Enumerated constants diff --git a/tdms/include/interpolate_Efield.h b/tdms/include/interpolate_Efield.h deleted file mode 100644 index 901011794..000000000 --- a/tdms/include/interpolate_Efield.h +++ /dev/null @@ -1,48 +0,0 @@ -/** - * @file interpolate_Efield.h - * @brief Functions for interpolating the Electric field. - */ -#pragma once - -/** - * @brief Interpolate the Ex field component to the centre of a Yee cell - * - * @param[in] Exy,Exz split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nI Number of Yee cells in the x-dimension - * @param[out] Ex Interpolated value of the Ex field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainEx(double ***Exy, double ***Exz, int i, int j, int k, int nI, double *Ex); - -/** - * @brief Interpolate the Ey field component to the centre of a Yee cell - * - * @param[in] Eyx,Eyz split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nJ Number of Yee cells in the y-dimension - * @param[out] Ey Interpolated value of the Ey field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainEy(double ***Eyx, double ***Eyz, int i, int j, int k, int nJ, double *Ey); - -/** - * @brief Interpolate the Ez field component to the centre of a Yee cell - * - * @param[in] Ezx,Ezy split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nK Number of Yee cells in the z-dimension - * @param[out] Ez Interpolated value of the Ez field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainEz(double ***Ezx, double ***Ezy, int i, int j, int k, int nK, double *Ez); - -/** - * @brief Interpolate the E-field to the centre of Yee cell i,j,k - * - * @param[in] Exy,Exz,Eyx,Eyz,Ezx,Ezy Split components of the Yee cell - * @param[in] i,j.k Index of the Yee cell to interpolate to the centre of - * @param[in] nI,nJ,nK Number of Yee cells in the i,j,k directions (respectively) - * @param[out] Ex,Ey,Ez Interpolated values of the x,y,z (respectively) field component - */ -void interpolateTimeDomainEField(double ***Exy, double ***Exz, double ***Eyx, - double ***Eyz, double ***Ezx, double ***Ezy, - int i, int j, int k, int nI, int nJ, int nK, - double *Ex, double *Ey, double *Ez); diff --git a/tdms/include/interpolate_Hfield.h b/tdms/include/interpolate_Hfield.h deleted file mode 100644 index f2a7c9316..000000000 --- a/tdms/include/interpolate_Hfield.h +++ /dev/null @@ -1,61 +0,0 @@ -/** - * @file interpolate_Hfield.h - * @brief Functions for interpolating the Electric field. - */ -#pragma once - -/** - * @brief Interpolate the x-component of the magnetic field to the centre of Yee cell i,j,k - * - * Interpolation for the magnetic field must be done across two dimensions. - * The stored Hx-component at position (i, j+Dy, k+Dz) is associated to the Yee cell (i,j,k). - * To interpolate Hx to the centre of cell (i,j,k), we must therefore interpolate in the y-direction, then the z-direction, or vice-versa. - * - * @param[in] Hxy,Hxz Split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nJ,nK Number of Yee cells in the j,k directions respectively - * @param[out] Hx Interpolated value of the Hx field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainHx(double ***Hxy, double ***Hxz, int i, int j, int k, int nJ, int nK, double *Hx); - -/** - * @brief Interpolate the y-component of the magnetic field to the centre of Yee cell i,j,k - * - * Interpolation for the magnetic field must be done across two dimensions. - * The stored Hy-component at position (i+Dx, j, k+Dz) is associated to the Yee cell (i,j,k). - * To interpolate Hy to the centre of cell (i,j,k), we must therefore interpolate in the x-direction, then the z-direction, or vice-versa. - * - * @param[in] Hyx,Hyz Split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nI,nK Number of Yee cells in the i,k directions respectively - * @param[out] Hy Interpolated value of the Hy field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainHy(double ***Hxy, double ***Hxz, int i, int j, int k, int nI, int nK, double *Hy); - -/** - * @brief Interpolate the z-component of the magnetic field to the centre of Yee cell i,j,k - * - * Interpolation for the magnetic field must be done across two dimensions. - * The stored Hz-component at position (i+Dx, j+Dy, k) is associated to the Yee cell (i,j,k). - * To interpolate Hz to the centre of cell (i,j,k), we must therefore interpolate in the x-direction, then the y-direction, or vice-versa. - * - * @param[in] Hzx,Hzy Split components of the Yee cell - * @param[in] i,j,k Yee cell index - * @param[in] nI,nJ Number of Yee cells in the i,j directions respectively - * @param[out] Hz Interpolated value of the Hz field at centre of Yee cell i,j,k - */ -void interpolateTimeDomainHz(double ***Hxy, double ***Hxz, int i, int j, int k, int nI, int nJ, double *Hz); - -/** - * @brief Interpolate the H-field to the centre of Yee cell i,j,k - * - * @param[in] Hxy,Hxz,Hyx,Hyz,Hzx,Hzy Split components of the Yee cell - * @param[in] i,j.k Index of the Yee cell to interpolate to the centre of - * @param[in] nI,nJ,nK Number of Yee cells in the i,j,k directions (respectively) - * @param[out] Hx,Hy,Hz Interpolated values of the x,y,z (respectively) field component - */ -void interpolateTimeDomainHField(double ***Hxy, double ***Hxz, double ***Hyx, - double ***Hyz, double ***Hzx, double ***Hzy, - int i, int j, int k, int nI, int nJ, int nK, - double *Hx, double *Hy, double *Hz); - diff --git a/tdms/include/interpolation_methods.h b/tdms/include/interpolation_methods.h index fa8b2d224..3a32c599f 100644 --- a/tdms/include/interpolation_methods.h +++ b/tdms/include/interpolation_methods.h @@ -1,5 +1,5 @@ /** - * @file interpolation_methods.cpp + * @file interpolation_methods.h * @author William Graham (ccaegra@ucl.ac.uk) * @brief InterpScheme class methods and supporting functions * @@ -7,6 +7,8 @@ */ #pragma once +#include + /*Use cubic interpolation to interpolate between the middle 2 of 4 points * v0 v1 v2 v3 * o o x o o @@ -78,7 +80,7 @@ enum scheme_value BAND_LIMITED_5 = 5, // use bandlimited interpolation w/ interp position = 5 BAND_LIMITED_6 = 3, // use bandlimited interpolation w/ interp position = 6 BAND_LIMITED_7 = -1, // use bandlimited interpolation w/ interp position = 7 [Only applicable if we want to extend beyond the final Yee cell, current code functionality is to throw an error in the case where this would be used.] - BAND_LIMITED_CELL_ZERO = -2, // use bandlimited interpolation to interpolate to the centre of Yee cell 0 [implemented, but current code functionality is to throw an error here] + BAND_LIMITED_CELL_ZERO = -2, // use bandlimited interpolation to interpolate to the centre of Yee cell <0 [implemented, but current code functionality is to throw an error here] CUBIC_INTERP_MIDDLE = 2, // cubic interpolation to middle 2 of 4 points (interp1) CUBIC_INTERP_FIRST = 1, // cubic interpolation to first 2 of 4 points (interp2) CUBIC_INTERP_LAST = 0 // cubic interpolation to last 2 of 4 points (interp3) @@ -129,14 +131,24 @@ class InterpolationScheme { * @brief Executes the interpolation scheme on the data provided * * The interpolation schemes are all of the form - * interpolated_value = \sum_{i=first_nonzero_coeff}^{last_nonzero_coeff} scheme_coeffs[i] * v[i], + * interpolated_value = \sum_{i=0}^{7} scheme_coeffs[i] * v[i], * so provided that the coefficients have been set correctly in construction (and the data gathered appropriately), we can run the same for loop for each interpolation scheme. * - * @param v Sample datapoints to use in interpolation + * For slight speedup, the actual sum performed loops over those i such that + * 0 <= first_nonzero_coeff <= i <= last_nonzero_coeff <= 7. + * + * @param v Sample datapoints to use in interpolation; v[0] should be the first of 8 values * @param offset [Default 0] Read buffer from v[offset] rather than v[0] * @return double Interpolated value */ - double interpolate(const double *v, const int offset = 0) const; + template + T interpolate(const T *v, const int offset = 0) const { + T interp_value = 0.; + for (int ind = first_nonzero_coeff; ind <= last_nonzero_coeff; ind++) { + interp_value += scheme_coeffs[ind] * v[ind + offset]; + } + return interp_value; + }; /** * @brief Determines whether another interpScheme has greater value than this one diff --git a/tdms/src/arrays.cpp b/tdms/src/arrays.cpp index 58aa5bfa9..7085c54c5 100644 --- a/tdms/src/arrays.cpp +++ b/tdms/src/arrays.cpp @@ -7,6 +7,22 @@ using namespace std; using namespace tdms_math_constants; +void XYZTensor3D::allocate(int I_total, int J_total, int K_total) { + x = (double ***) malloc(K_total * sizeof(double **)); + y = (double ***) malloc(K_total * sizeof(double **)); + z = (double ***) malloc(K_total * sizeof(double **)); + for (int k = 0; k < K_total; k++) { + x[k] = (double **) malloc(J_total * sizeof(double *)); + y[k] = (double **) malloc(J_total * sizeof(double *)); + z[k] = (double **) malloc(J_total * sizeof(double *)); + for (int j = 0; j < J_total; j++) { + x[k][j] = (double *) malloc(I_total * sizeof(double)); + y[k][j] = (double *) malloc(I_total * sizeof(double)); + z[k][j] = (double *) malloc(I_total * sizeof(double)); + } + } +} + void XYZVectors::set_ptr(const char c, double* ptr){ switch (c) { case 'x': {x = ptr; break;} @@ -216,6 +232,27 @@ void Tensor3D::zero() { } } +template<> +double Tensor3D::frobenius() { + double norm_val = 0.; + for (int i1 = 0; i1 < n_layers; i1++) { + for (int i2 = 0; i2 < n_cols; i2++) { + for (int i3 = 0; i3 < n_rows; i3++) { norm_val += tensor[i1][i2][i3] * tensor[i1][i2][i3]; } + } + } + return sqrt(norm_val); +} +template<> +double Tensor3D>::frobenius() { + double norm_val = 0.; + for (int i1 = 0; i1 < n_layers; i1++) { + for (int i2 = 0; i2 < n_cols; i2++) { + for (int i3 = 0; i3 < n_rows; i3++) { norm_val += norm(tensor[i1][i2][i3]); } + } + } + return sqrt(norm_val); +} + void DTilde::set_component(Tensor3D> &tensor, const mxArray *ptr, const string &name, int n_rows, int n_cols){ diff --git a/tdms/src/electric_field.cpp b/tdms/src/electric_field.cpp index ed8959a2d..e00568e1d 100644 --- a/tdms/src/electric_field.cpp +++ b/tdms/src/electric_field.cpp @@ -1,11 +1,11 @@ #include -#include "globals.h" + #include "field.h" +#include "globals.h" using namespace std; using namespace tdms_math_constants; - void ElectricField::add_to_angular_norm(double f, int n, int Nt, SimulationParameters ¶ms) { angular_norm += phasor_norm(f, n, params.omega_an, params.dt, Nt); } diff --git a/tdms/src/fields/base.cpp b/tdms/src/fields/base.cpp index 9a8098c98..cc082c5af 100644 --- a/tdms/src/fields/base.cpp +++ b/tdms/src/fields/base.cpp @@ -27,8 +27,19 @@ void Field::normalise_volume() { } } -void Field::zero() { +Field::Field(int I_total, int J_total, int K_total) { + I_tot = I_total; + J_tot = J_total; + K_tot = K_total; +} +void Field::allocate() { + for (auto arr : {&real, &imag}) { + arr->allocate(I_tot+1, J_tot+1, K_tot+1); + } +}; + +void Field::zero() { for (auto &arr : {real, imag}) for (char c : {'x', 'y', 'z'}) for (int k = 0; k < K_tot; k++) diff --git a/tdms/src/fields/electric.cpp b/tdms/src/fields/electric.cpp index 582a6555f..2341778ab 100644 --- a/tdms/src/fields/electric.cpp +++ b/tdms/src/fields/electric.cpp @@ -1,8 +1,108 @@ #include "field.h" +#include "globals.h" +#include "interpolation_methods.h" using namespace std; - +using namespace tdms_math_constants; double ElectricField::phase(int n, double omega, double dt){ return omega * ((double) n + 1) * dt; } + +/* +You may notice several +1's appearing when the interpolation data (interp_data) arrays are constructed from the field values. +This is because E.x[k][j][i] is, by construction of the Grid class, not the value associated to Yee cell (i,j,k) but rather cell (i-1,j,k). Similarly Ey[k][j][i] is the value of cell (i,j-1,k) and Ez[k][j][i] of cell (i,j,k-1). +*/ + +complex ElectricField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { + const InterpolationScheme *scheme; + // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind + complex interp_data[8]; + switch (d) { + case X: + // determine the interpolation scheme to use + scheme = &(best_scheme(I_tot, i)); + // now fill the interpolation data + // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = + real.x[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind] + + IMAGINARY_UNIT * imag.x[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind]; + } + break; + case Y: + // determine the interpolation scheme to use + scheme = &(best_scheme(J_tot, j)); + + // now fill the interpolation data + // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = + real.y[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i] + + IMAGINARY_UNIT * imag.y[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i]; + } + break; + case Z: + // determine the interpolation scheme to use + scheme = &(best_scheme(K_tot, k)); + + // now fill the interpolation data + // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = + real.z[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i] + + IMAGINARY_UNIT * imag.z[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i]; + } + break; + default: + throw runtime_error("Invalid axial direction selected for interpolation!\n"); + break; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); +} + +double ElectricSplitField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { + const InterpolationScheme *scheme; + // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind + double interp_data[8]; + + switch (d) { + case X: + scheme = &(best_scheme(I_tot, i)); + // now fill the interpolation data + // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = xy[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind] + + xz[k][j][i + 1 - scheme->number_of_datapoints_to_left + ind]; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); + break; + case Y: + scheme = &(best_scheme(J_tot, j)); + // now fill the interpolation data + // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = yx[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i] + + yz[k][j + 1 - scheme->number_of_datapoints_to_left + ind][i]; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); + break; + case Z: + scheme = &(best_scheme(K_tot, k)); + // now fill the interpolation data + // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = scheme->first_nonzero_coeff; ind <= scheme->last_nonzero_coeff; ind++) { + interp_data[ind] = zx[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i] + + zy[k + 1 - scheme->number_of_datapoints_to_left + ind][j][i]; + } + // now run the interpolation scheme and place the result into the output + return scheme->interpolate(interp_data); + break; + default: + throw runtime_error("Invalid axial direction selected for interpolation!\n"); + break; + } +} diff --git a/tdms/src/fields/magnetic.cpp b/tdms/src/fields/magnetic.cpp index 66d691b06..e9a367e6a 100644 --- a/tdms/src/fields/magnetic.cpp +++ b/tdms/src/fields/magnetic.cpp @@ -1,8 +1,374 @@ #include "field.h" +#include "interpolation_methods.h" +#include "globals.h" using namespace std; - +using namespace tdms_math_constants; double MagneticField::phase(int n, double omega, double dt){ return omega * ((double) n + 0.5) * dt; // 0.5 added because it's known half a time step after E } + +/* 2D INTERPOLATION SCHEMES (FOR THE MAGNETIC FIELD IN 3D SIMULATIONS) + +Unlike the E-field, the H-field components associated with Yee cell i,j,k are _not_ aligned with the centre of the Yee cells. +Instead, the position of the field components (relative to the Yee cell centre is): + +Hx | (0.0, 0.5, 0.5) .* (Dx, Dy, Dx) +Hy | (0.5, 0.0, 0.5) .* (Dx, Dy, Dz) +Hz | (0.5, 0.5, 0.0) .* (Dx, Dy, Dz) + +where Dx, Dy, Dz are the dimensions of the Yee cell. +This requires us to interpolate twice to recover (any of the) field components at the centre of Yee cell i,j,k, when running a 3D simulation. + +Henceforth, let {a,b,c} = {x,y,z} be some 1-to-1 assignment of the co-ordinate axes to the labels a,b,c. +The values Da, Db, Dc are the corresponding permutation of Dx, Dy, Dz. +Suppose that we wish to interpolate the Ha field component to the centre a particular Yee cell. + +We will use the notation (a_i,b_j,c_k) for the Yee cell indices, although bear in mind that this does not reflect the order the indices appear in the code. +For example; if a = y, b = x, c = z, then the value of Ha at cell (a_i,b_j,c_k) is accessed via Ha[c_k][a_i][b_j], due to the interchanging of the x and y directions. +Similarly; we will write Ha[a_i, b_j, c_k] to refer to the value of Ha associated to cell (a_i,b_j,c_k). + +Suppose now that we have selected cell a_i,b_j,c_k to interpolate to the centre of. +We must interpolate in the b-direction, then c-direction, or vice-versa. +For optimal accuracy we must determine the best interpolation scheme we can use in each direction at cell (a_i,b_j,c_k), and use the WORSE scheme second. + +Let us assume WLOG that the b-direction interpolation scheme (b_scheme) is inferior to that of the c-direction (c_scheme). +We now need to interpolate Ha in the c-direction to obtain the value of Ha at the spatial positions (a_i, cell_b + Db, c_k) where +b_j - b_scheme.number_of_datapoints_to_left + b_scheme.first_nonzero_coeff <= cell_b <= b_j - b_scheme.number_of_datapoints_to_left + b_scheme.last_nonzero_coeff. + + Let cell_b be one particular index in this range. + To apply c_scheme to obtain the value of Ha at (a_i, cell_b + Db, c_k), we require the values Ha[a_i, cell_b, cell_c] where + c_k - c_scheme.number_of_datapoints_to_left + c_scheme.first_nonzero_coeff <= cell_c <= c_k - c_scheme.number_of_datapoints_to_left + c_scheme.last_nonzero_coeff. + These values are fed to c_scheme.interpolate, which provides us with Ha at the spatial position (a_i, cell_b + Db, c_k). + +Now with approximations of Ha at each (a_i, cell_b + Db, c_k), we can pass this information to b_scheme.interpolate to recover the value of Ha at the centre of Yee cell (a_i, b_j, c_k). + +In a 2D simulation, J_tot = 0. This means we cannot interpolate in two directions for the Hx and Hz fields, since there is no y-direction to interpolate in. As a result, we only interpolate in the z-direction for Hx and x-direction for Hz when running a two-dimensional simulation. + +You may notice several +1's appearing when the interpolation data (interp_data) arrays are constructed from the field values. +This is because H.x[k][j][i] is, by construction of the Grid class, not the value associated to Yee cell (i,j,k) but rather cell (i,j-1,k-1). Similarly Ey[k][j][i] is the value of cell (i-1,j,k) and Ez[k][j][i] of cell (i-1,j-1,k). +*/ + +complex MagneticField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { + // this data will be passed to the second interpolation scheme + complex data_for_second_scheme[8]; + // this data will hold values for the interpolation in the first interpolation scheme + complex data_for_first_scheme[8]; + // the interpolation schemes that are to be used + const InterpolationScheme *b_scheme, *c_scheme; + + switch (d) { + case X: + // Associations: a = x, b = y, c = z + c_scheme = &(best_scheme(K_tot, k)); + b_scheme = &(best_scheme(J_tot, j)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the z-direction first, then in y + for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { + // this is the j-index of the cell we are looking at + int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // the k-index of the current cell we are looking at (readability, don't need to define this here) + int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + // gather the data for interpolating in the z dimension + data_for_first_scheme[kk] = + real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; + } + // interpolate in z to obtain a value for the Hx field at position (i, cell_j+Dy, k) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[jj] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the y-direction first, then in z + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // this is the k-index of the cell we are looking at + int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme + for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + // gather the data for interpolating in the y dimension + data_for_first_scheme[jj] = + real.x[cell_k][cell_j][i] + IMAGINARY_UNIT * imag.x[cell_k][cell_j][i]; + } + // interpolate in y to obtain a value for the Hx field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[kk] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the z-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } + break; + case Y: + // Associations: a = y, b = z, c = x + c_scheme = &(best_scheme(I_tot, i)); + b_scheme = &(best_scheme(K_tot, k)); + + if (b_scheme->is_better_than(*c_scheme)) { + // we will be interpolating in the z-direction first, then in x + for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { + // this is the i-index of the cell we are looking at + int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + // determine the Hy values of cells (cell_i, j, k-z_scheme.index-1) through (cell_i, j, k-z_scheme.index-1+7), and interpolate them via z_scheme + for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { + // the k-index of the current cell we are looking at (readability, don't need to define this here) + int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + // gather the data for interpolating in the z dimension + data_for_first_scheme[kk] = + real.y[cell_k][j][cell_i] + IMAGINARY_UNIT * imag.y[cell_k][j][cell_i]; + } + // interpolate in z to obtain a value for the Hy field at position (cell_i+Dx, j, k) + // place this into the appropriate index in the data being passed to the x_scheme + data_for_second_scheme[ii] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the x-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the x-direction first, then in z + for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { + // this is the k-index of the cell we are looking at + int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + // determine the Hy values of cells (i - x_scheme.index-1, j, cell_k) through (i- x_scheme.index-1+7, j, cell_k), and interpolate them via x_scheme + for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { + // the i-index of the current cell we are looking at (readability, don't need to define this here) + int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + // gather the data for interpolating in the x dimension + data_for_first_scheme[ii] = + real.y[cell_k][j][cell_i] + IMAGINARY_UNIT * imag.y[cell_k][j][cell_i]; + } + // interpolate in x to obtain a value for the Hy field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[kk] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the z-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } + break; + case Z: + // Associations: a = z, b = x, c = y + b_scheme = &(best_scheme(I_tot, i)); + c_scheme = &(best_scheme(J_tot, j)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the y-direction first, then in x + for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { + // this is the i-index of the cell we are looking at + int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + // gather the data for interpolating in the y dimension + data_for_first_scheme[jj] = + real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; + } + // interpolate in y to obtain a value for the Hz field at position (cell_i+Dx, j, k) + // place this into the appropriate index in the data being passed to the x_scheme + data_for_second_scheme[ii] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the x-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the x-direction first, then in y + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // this is the j-index of the cell we are looking at + int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme + for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { + // the i-index of the current cell we are looking at (readability, don't need to define this here) + int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + // gather the data for interpolating in the x dimension + data_for_first_scheme[ii] = + real.z[k][cell_j][cell_i] + IMAGINARY_UNIT * imag.z[k][cell_j][cell_i]; + } + // interpolate in x to obtain a value for the Hz field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[jj] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } + break; + default: + throw runtime_error("Error: invalid axial direction selected for interpolation\n"); + break; + } +} + +double MagneticSplitField::interpolate_to_centre_of(AxialDirection d, int i, int j, int k) { + const InterpolationScheme *b_scheme, *c_scheme; + // this data will be passed to the second interpolation scheme + double data_for_second_scheme[8]; + // this data will hold values for the interpolation in the first interpolation scheme + double data_for_first_scheme[8]; + + switch (d) { + case X: + if (J_tot <= 0) { + // in a 2D simulation, we must interpolate in the z-direction to recover Hx, due to the magnetic-field offsets from the centre + b_scheme = &(best_scheme(K_tot, k)); + // now fill the interpolation data + // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { + data_for_first_scheme[ind] = + xy[k + 1 - b_scheme->number_of_datapoints_to_left + ind][j][i] + + xz[k + 1 - b_scheme->number_of_datapoints_to_left + ind][j][i]; + } + + // now run the interpolation scheme and place the result into the output + return b_scheme->interpolate(data_for_first_scheme); + } else { + // Associations: a = x, b = y, c = z + c_scheme = &(best_scheme(K_tot, k)); + b_scheme = &(best_scheme(J_tot, j)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the z-direction first, then in y + for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { + // this is the j-index of the cell we are looking at + int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // the k-index of the current cell we are looking at (readability, don't need to define this here) + int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + // gather the data for interpolating in the z dimension + data_for_first_scheme[kk] = xy[cell_k][cell_j][i] + xz[cell_k][cell_j][i]; + } + // interpolate in z to obtain a value for the Hx field at position (i, cell_j+Dy, k) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[jj] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the y-direction first, then in z + for (int kk = c_scheme->first_nonzero_coeff; kk <= c_scheme->last_nonzero_coeff; kk++) { + // this is the k-index of the cell we are looking at + int cell_k = k + 1 - c_scheme->number_of_datapoints_to_left + kk; + // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme + for (int jj = b_scheme->first_nonzero_coeff; jj <= b_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + int cell_j = j + 1 - b_scheme->number_of_datapoints_to_left + jj; + // gather the data for interpolating in the y dimension + data_for_first_scheme[jj] = xy[cell_k][cell_j][i] + xz[cell_k][cell_j][i]; + } + // interpolate in y to obtain a value for the Hx field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[kk] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the z-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } + } + break; + case Y: + // we can always interpolate in two directions for Hy, since even in a 2D simulation the x- and z-directions still exist + // Associations: a = y, b = z, c = x + c_scheme = &(best_scheme(I_tot, i)); + b_scheme = &(best_scheme(K_tot, k)); + + if (b_scheme->is_better_than(*c_scheme)) { + // we will be interpolating in the z-direction first, then in x + for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { + // this is the i-index of the cell we are looking at + int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + // determine the Hy values of cells (cell_i, j, k-z_scheme.index-1) through (cell_i, j, k-z_scheme.index-1+7), and interpolate them via z_scheme + for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { + // the k-index of the current cell we are looking at (readability, don't need to define this here) + int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + // gather the data for interpolating in the z dimension + data_for_first_scheme[kk] = yx[cell_k][j][cell_i] + yz[cell_k][j][cell_i]; + } + // interpolate in z to obtain a value for the Hy field at position (cell_i+Dx, j, k) + // place this into the appropriate index in the data being passed to the x_scheme + data_for_second_scheme[ii] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the x-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the x-direction first, then in z + for (int kk = b_scheme->first_nonzero_coeff; kk <= b_scheme->last_nonzero_coeff; kk++) { + // this is the k-index of the cell we are looking at + int cell_k = k + 1 - b_scheme->number_of_datapoints_to_left + kk; + // determine the Hy values of cells (i - x_scheme.index-1, j, cell_k) through (i- x_scheme.index-1+7, j, cell_k), and interpolate them via x_scheme + for (int ii = c_scheme->first_nonzero_coeff; ii <= c_scheme->last_nonzero_coeff; ii++) { + // the i-index of the current cell we are looking at (readability, don't need to define this here) + int cell_i = i + 1 - c_scheme->number_of_datapoints_to_left + ii; + // gather the data for interpolating in the x dimension + data_for_first_scheme[ii] = yx[cell_k][j][cell_i] + yz[cell_k][j][cell_i]; + } + // interpolate in x to obtain a value for the Hy field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[kk] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the z-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } + break; + case Z: + if (J_tot <= 0) { + // in a 2D simulation, we must interpolate in the x-direction to recover Hz, due to the magnetic-field offsets from the centre + b_scheme = &(best_scheme(I_tot, i)); + // now fill the interpolation data + // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation + for (int ind = b_scheme->first_nonzero_coeff; ind <= b_scheme->last_nonzero_coeff; ind++) { + data_for_first_scheme[ind] = + zx[k][j][i + 1 - b_scheme->number_of_datapoints_to_left + ind] + + zy[k][j][i + 1 - b_scheme->number_of_datapoints_to_left + ind]; + } + + // now run the interpolation scheme and place the result into the output + return b_scheme->interpolate(data_for_first_scheme); + } else { + // Associations: a = z, b = x, c = y + b_scheme = &(best_scheme(I_tot, i)); + c_scheme = &(best_scheme(J_tot, j)); + + if (c_scheme->is_better_than(*b_scheme)) { + // we will be interpolating in the y-direction first, then in x + for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { + // this is the i-index of the cell we are looking at + int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // the j-index of the current cell we are looking at (readability, don't need to define this here) + int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + // gather the data for interpolating in the y dimension + data_for_first_scheme[jj] = zx[k][cell_j][cell_i] + zy[k][cell_j][cell_i]; + } + // interpolate in y to obtain a value for the Hz field at position (cell_i+Dx, j, k) + // place this into the appropriate index in the data being passed to the x_scheme + data_for_second_scheme[ii] = c_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the x-direction to the centre of Yee cell (i,j,k) + return b_scheme->interpolate(data_for_second_scheme); + } else { + // we will be interpolating in the x-direction first, then in y + for (int jj = c_scheme->first_nonzero_coeff; jj <= c_scheme->last_nonzero_coeff; jj++) { + // this is the j-index of the cell we are looking at + int cell_j = j + 1 - c_scheme->number_of_datapoints_to_left + jj; + // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme + for (int ii = b_scheme->first_nonzero_coeff; ii <= b_scheme->last_nonzero_coeff; ii++) { + // the i-index of the current cell we are looking at (readability, don't need to define this here) + int cell_i = i + 1 - b_scheme->number_of_datapoints_to_left + ii; + // gather the data for interpolating in the x dimension + data_for_first_scheme[ii] = zx[k][cell_j][cell_i] + zy[k][cell_j][cell_i]; + } + // interpolate in x to obtain a value for the Hz field at position (i, j, cell_k+Dz) + // place this into the appropriate index in the data being passed to the y_scheme + data_for_second_scheme[jj] = b_scheme->interpolate(data_for_first_scheme); + } + // now interpolate in the y-direction to the centre of Yee cell (i,j,k) + return c_scheme->interpolate(data_for_second_scheme); + } + } + break; + default: + throw runtime_error("Error: invalid axial direction selected for interpolation\n"); + } +} diff --git a/tdms/src/fields/split.cpp b/tdms/src/fields/split.cpp index 2a8569e4d..5f1633506 100644 --- a/tdms/src/fields/split.cpp +++ b/tdms/src/fields/split.cpp @@ -1,5 +1,6 @@ #include #include "field.h" +#include "interpolation_methods.h" using namespace std; diff --git a/tdms/src/interpolate_Efield.cpp b/tdms/src/interpolate_Efield.cpp deleted file mode 100644 index 38f08fc5e..000000000 --- a/tdms/src/interpolate_Efield.cpp +++ /dev/null @@ -1,66 +0,0 @@ -#include "interpolate_Efield.h" -#include "interpolation_methods.h" - - -void interpolateTimeDomainEx(double ***Exy, double ***Exz, int i, int j, int k, int nI, double *Ex) { - - // determine the interpolation scheme to use - const InterpolationScheme &scheme = best_scheme(nI, i); - // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind - double interp_data[8]; - - // now fill the interpolation data - // i - (scheme.number_of_datapoints_to_left) is the index of the Yee cell that plays the role of v0 in the interpolation - for(int ind=scheme.first_nonzero_coeff; ind<=scheme.last_nonzero_coeff; ind++) { - interp_data[ind] = Exy[k][j][i - scheme.number_of_datapoints_to_left + ind] + Exz[k][j][i - scheme.number_of_datapoints_to_left + ind]; - } - - // now run the interpolation scheme and place the result into the output - *Ex = scheme.interpolate(interp_data); -} - -void interpolateTimeDomainEy(double ***Eyx, double ***Eyz, int i, int j, int k, int nJ, double *Ey) -{ - // determine the interpolation scheme to use - const InterpolationScheme &scheme = best_scheme(nJ, j); - // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind - double interp_data[8]; - - // now fill the interpolation data - // j - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation - for (int ind = scheme.first_nonzero_coeff; ind <= scheme.last_nonzero_coeff; ind++) - { - interp_data[ind] = Eyx[k][j - scheme.number_of_datapoints_to_left + ind][i] + Eyz[k][j - scheme.number_of_datapoints_to_left + ind][j]; - } - - // now run the interpolation scheme and place the result into the output - *Ey = scheme.interpolate(interp_data); -} - -void interpolateTimeDomainEz(double ***Ezx, double ***Ezy, int i, int j, int k, int nK, double *Ez) { - - // determine the interpolation scheme to use - const InterpolationScheme &scheme = best_scheme(nK, k); - // prepare input data - if using a cubic scheme we have reserved more memory than necessary but nevermind - double interp_data[8]; - - // now fill the interpolation data - // k - scheme.number_of_datapoints_to_left is the index of the Yee cell that plays the role of v0 in the interpolation - for (int ind = scheme.first_nonzero_coeff; ind <= scheme.last_nonzero_coeff; ind++) - { - interp_data[ind] = Ezx[k - scheme.number_of_datapoints_to_left + ind][j][i] + Ezy[k - scheme.number_of_datapoints_to_left + ind][j][i]; - } - - // now run the interpolation scheme and place the result into the output - *Ez = scheme.interpolate(interp_data); -} - -void interpolateTimeDomainEField(double ***Exy, double ***Exz, double ***Eyx, - double ***Eyz, double ***Ezx, double ***Ezy, - int i, int j, int k, int nI, int nJ, int nK, - double *Ex, double *Ey, double *Ez) -{ - interpolateTimeDomainEx(Exy, Exz, i, j, k, nI, Ex); - interpolateTimeDomainEy(Eyx, Eyz, i, j, k, nJ, Ey); - interpolateTimeDomainEz(Ezx, Ezy, i, j, k, nK, Ez); -} diff --git a/tdms/src/interpolate_Hfield.cpp b/tdms/src/interpolate_Hfield.cpp deleted file mode 100644 index c411de3f2..000000000 --- a/tdms/src/interpolate_Hfield.cpp +++ /dev/null @@ -1,232 +0,0 @@ -#include -#include "interpolate_Hfield.h" -#include "interpolation_methods.h" - -/* INTERPOLATION SCHEMES FOR THE MAGNETIC FIELD - -Unlike the E-field, the H-field components associated with Yee cell i,j,k are _not_ aligned with the centre of the Yee cells. -Instead, the position of the field components (relative to the Yee cell centre is): - -Hx | (0.0, 0.5, 0.5) .* (Dx, Dy, Dx) -Hy | (0.5, 0.0, 0.5) .* (Dx, Dy, Dz) -Hz | (0.5, 0.5, 0.0) .* (Dx, Dy, Dz) - -where Dx, Dy, Dz are the dimensions of the Yee cell. -This requires us to interpolate twice to recover (any of the) field components at the centre of Yee cell i,j,k. - -Henceforth, let {a,b,c} = {x,y,z} be some 1-to-1 assignment of the co-ordinate axes to the labels a,b,c. -The values Da, Db, Dc are the corresponding permutation of Dx, Dy, Dz. -Suppose that we wish to interpolate the Ha field component to the centre a particular Yee cell. - -We will use the notation (a_i,b_j,c_k) for the Yee cell indices, although bear in mind that this does not reflect the order the indices appear in the code. -For example; if a = y, b = x, c = z, then the value of Ha at cell (a_i,b_j,c_k) is accessed via Ha[c_k][a_i][b_j], due to the interchanging of the x and y directions. -Similarly; we will write Ha[a_i, b_j, c_k] to refer to the value of Ha associated to cell (a_i,b_j,c_k). - -Suppose now that we have selected cell a_i,b_j,c_k to interpolate to the centre of. -We must interpolate in the b-direction, then c-direction, or vice-versa. -For optimal accuracy we must determine the best interpolation scheme we can use in each direction at cell (a_i,b_j,c_k), and use the WORSE scheme second. - -Let us assume WLOG that the b-direction interpolation scheme (b_scheme) is inferior to that of the c-direction (c_scheme). -We now need to interpolate Ha in the c-direction to obtain the value of Ha at the spatial positions (a_i, cell_b + Db, c_k) where -b_j - b_scheme.number_of_datapoints_to_left + b_scheme.first_nonzero_coeff <= cell_b <= b_j - b_scheme.number_of_datapoints_to_left + b_scheme.last_nonzero_coeff. - - Let cell_b be one particular index in this range. - To apply c_scheme to obtain the value of Ha at (a_i, cell_b + Db, c_k), we require the values Ha[a_i, cell_b, cell_c] where - c_k - c_scheme.number_of_datapoints_to_left + c_scheme.first_nonzero_coeff <= cell_c <= c_k - c_scheme.number_of_datapoints_to_left + c_scheme.last_nonzero_coeff. - These values are fed to c_scheme.interpolate, which provides us with Ha at the spatial position (a_i, cell_b + Db, c_k). - -Now with approximations of Ha at each (a_i, cell_b + Db, c_k), we can pass this information to b_scheme.interpolate to recover the value of Ha at the centre of Yee cell (a_i, b_j, c_k). -*/ - -void interpolateTimeDomainHx(double ***Hxy, double ***Hxz, int i, int j, int k, int nJ, int nK, double *Hx) -{ - // Associations: a = x, b = y, c = z - - // determine the z-direction scheme - const InterpolationScheme &z_scheme = best_scheme(nK, k); - // determine the y-direction scheme - const InterpolationScheme &y_scheme = best_scheme(nJ, j); - - // this data will be passed to the second interpolation scheme - double data_for_second_scheme[8]; - // this data will hold values for the interpolation in the first interpolation scheme - double data_for_first_scheme[8]; - - // which of z_scheme and y_scheme is WORSE? We will interpolate in this direction SECOND - if (z_scheme.is_better_than(y_scheme)) - { - // we will be interpolating in the z-direction first, then in y - for (int jj = y_scheme.first_nonzero_coeff; jj <= y_scheme.last_nonzero_coeff; jj++) - { - // this is the j-index of the cell we are looking at - int cell_j = j - y_scheme.number_of_datapoints_to_left + jj; - // determine the Hx values of cells (i, cell_j, k-z_scheme.index-1) through (i, cell_j, k-z_scheme.index-1+7), and interpolate them via z_scheme - for (int kk = z_scheme.first_nonzero_coeff; kk <= z_scheme.last_nonzero_coeff; kk++) - { - // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k - z_scheme.number_of_datapoints_to_left + kk; - // gather the data for interpolating in the z dimension - data_for_first_scheme[kk] = Hxy[cell_k][cell_j][i] + Hxz[cell_k][cell_j][i]; - } - // interpolate in z to obtain a value for the Hx field at position (i, cell_j+Dy, k) - // place this into the appropriate index in the data being passed to the y_scheme - data_for_second_scheme[jj] = z_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the y-direction to the centre of Yee cell (i,j,k) - *Hx = y_scheme.interpolate(data_for_second_scheme); - } - else - { - // we will be interpolating in the y-direction first, then in z - for (int kk = z_scheme.first_nonzero_coeff; kk <= z_scheme.last_nonzero_coeff; kk++) - { - // this is the k-index of the cell we are looking at - int cell_k = k - z_scheme.number_of_datapoints_to_left + kk; - // determine the Hx values of cells (i, j - y_scheme.index-1, cell_k) through (i, j - y_scheme.index-1+7, cell_k), and interpolate them via y_scheme - for (int jj = y_scheme.first_nonzero_coeff; jj <= y_scheme.last_nonzero_coeff; jj++) - { - // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j - y_scheme.number_of_datapoints_to_left + jj; - // gather the data for interpolating in the y dimension - data_for_first_scheme[jj] = Hxy[cell_k][cell_j][i] + Hxz[cell_k][cell_j][i]; - } - // interpolate in y to obtain a value for the Hx field at position (i, j, cell_k+Dz) - // place this into the appropriate index in the data being passed to the y_scheme - data_for_second_scheme[kk] = y_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the z-direction to the centre of Yee cell (i,j,k) - *Hx = z_scheme.interpolate(data_for_second_scheme); - } -} - -void interpolateTimeDomainHy(double ***Hyx, double ***Hyz, int i, int j, int k, int nI, int nK, double *Hy) -{ - // Associations: a = y, b = z, c = x - - // determine the x-direction scheme - const InterpolationScheme &x_scheme = best_scheme(nI, i); - // determine the z-direction scheme - const InterpolationScheme &z_scheme = best_scheme(nK, k); - - // this data will be passed to the second interpolation scheme - double data_for_second_scheme[8]; - // this data will hold values for the interpolation in the first interpolation scheme - double data_for_first_scheme[8]; - - // which of x_scheme and z_scheme is WORSE? We will interpolate in this direction SECOND - if (z_scheme.is_better_than(x_scheme)) - { - // we will be interpolating in the z-direction first, then in x - for (int ii = x_scheme.first_nonzero_coeff; ii <= x_scheme.last_nonzero_coeff; ii++) - { - // this is the i-index of the cell we are looking at - int cell_i = i - x_scheme.number_of_datapoints_to_left + ii; - // determine the Hy values of cells (cell_i, j, k-z_scheme.index-1) through (cell_i, j, k-z_scheme.index-1+7), and interpolate them via z_scheme - for (int kk = z_scheme.first_nonzero_coeff; kk <= z_scheme.last_nonzero_coeff; kk++) - { - // the k-index of the current cell we are looking at (readability, don't need to define this here) - int cell_k = k - z_scheme.number_of_datapoints_to_left + kk; - // gather the data for interpolating in the z dimension - data_for_first_scheme[kk] = Hyx[cell_k][j][cell_i] + Hyz[cell_k][j][cell_i]; - } - // interpolate in z to obtain a value for the Hy field at position (cell_i+Dx, j, k) - // place this into the appropriate index in the data being passed to the x_scheme - data_for_second_scheme[ii] = z_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the x-direction to the centre of Yee cell (i,j,k) - *Hy = x_scheme.interpolate(data_for_second_scheme); - } - else - { - // we will be interpolating in the x-direction first, then in z - for (int kk = z_scheme.first_nonzero_coeff; kk <= z_scheme.last_nonzero_coeff; kk++) - { - // this is the k-index of the cell we are looking at - int cell_k = k - z_scheme.number_of_datapoints_to_left + kk; - // determine the Hy values of cells (i - x_scheme.index-1, j, cell_k) through (i- x_scheme.index-1+7, j, cell_k), and interpolate them via x_scheme - for (int ii = x_scheme.first_nonzero_coeff; ii <= x_scheme.last_nonzero_coeff; ii++) - { - // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i - x_scheme.number_of_datapoints_to_left + ii; - // gather the data for interpolating in the x dimension - data_for_first_scheme[ii] = Hyx[cell_k][j][cell_i] + Hyz[cell_k][j][cell_i]; - } - // interpolate in x to obtain a value for the Hy field at position (i, j, cell_k+Dz) - // place this into the appropriate index in the data being passed to the y_scheme - data_for_second_scheme[kk] = x_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the z-direction to the centre of Yee cell (i,j,k) - *Hy = z_scheme.interpolate(data_for_second_scheme); - } -} - -void interpolateTimeDomainHz(double ***Hzx, double ***Hzy, int i, int j, int k, int nI, int nJ, double *Hz) -{ - // Associations: a = z, b = x, c = y - - // determine the x-direction scheme - const InterpolationScheme &x_scheme = best_scheme(nI, i); - // determine the y-direction scheme - const InterpolationScheme &y_scheme = best_scheme(nJ, j); - - // this data will be passed to the second interpolation scheme - double data_for_second_scheme[8]; - // this data will hold values for the interpolation in the first interpolation scheme - double data_for_first_scheme[8]; - - // which of x_scheme and y_scheme is WORSE? We will interpolate in this direction SECOND - if (y_scheme.is_better_than(x_scheme)) - { - // we will be interpolating in the y-direction first, then in x - for (int ii = x_scheme.first_nonzero_coeff; ii <= x_scheme.last_nonzero_coeff; ii++) - { - // this is the i-index of the cell we are looking at - int cell_i = i - x_scheme.number_of_datapoints_to_left + ii; - // determine the Hz values of cells (cell_i, j-y_scheme.index-1, k) through (cell_i, j-y_scheme.index-1+7, k), and interpolate them via y_scheme - for (int jj = y_scheme.first_nonzero_coeff; jj <= y_scheme.last_nonzero_coeff; jj++) - { - // the j-index of the current cell we are looking at (readability, don't need to define this here) - int cell_j = j - y_scheme.number_of_datapoints_to_left + jj; - // gather the data for interpolating in the y dimension - data_for_first_scheme[jj] = Hzx[k][cell_j][cell_i] + Hzy[k][cell_j][cell_i]; - } - // interpolate in y to obtain a value for the Hz field at position (cell_i+Dx, j, k) - // place this into the appropriate index in the data being passed to the x_scheme - data_for_second_scheme[ii] = y_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the x-direction to the centre of Yee cell (i,j,k) - *Hz = x_scheme.interpolate(data_for_second_scheme); - } - else - { - // we will be interpolating in the x-direction first, then in y - for (int jj = y_scheme.first_nonzero_coeff; jj <= y_scheme.last_nonzero_coeff; jj++) - { - // this is the j-index of the cell we are looking at - int cell_j = j - y_scheme.number_of_datapoints_to_left + jj; - // determine the Hz values of cells (i - x_scheme.index-1, cell_j, k) through (i- x_scheme.index-1+7, cell_j, k), and interpolate them via x_scheme - for (int ii = x_scheme.first_nonzero_coeff; ii <= x_scheme.last_nonzero_coeff; ii++) - { - // the i-index of the current cell we are looking at (readability, don't need to define this here) - int cell_i = i - x_scheme.number_of_datapoints_to_left + ii; - // gather the data for interpolating in the x dimension - data_for_first_scheme[ii] = Hzx[k][cell_j][cell_i] + Hzy[k][cell_j][cell_i]; - } - // interpolate in x to obtain a value for the Hz field at position (i, j, cell_k+Dz) - // place this into the appropriate index in the data being passed to the y_scheme - data_for_second_scheme[jj] = x_scheme.interpolate(data_for_first_scheme); - } - // now interpolate in the y-direction to the centre of Yee cell (i,j,k) - *Hz = y_scheme.interpolate(data_for_second_scheme); - } -} - -void interpolateTimeDomainHField(double ***Hxy, double ***Hxz, double ***Hyx, - double ***Hyz, double ***Hzx, double ***Hzy, - int i, int j, int k, int nI, int nJ, int nK, - double *Hx, double *Hy, double *Hz) -{ - interpolateTimeDomainHx(Hxy, Hxz, i, j, k, nJ, nK, Hx); - interpolateTimeDomainHy(Hyx, Hyz, i, j, k, nI, nK, Hy); - interpolateTimeDomainHz(Hzx, Hzy, i, j, k, nI, nJ, Hz); -} diff --git a/tdms/src/interpolation_methods.cpp b/tdms/src/interpolation_methods.cpp index db8d9e3e8..e8048b4a1 100644 --- a/tdms/src/interpolation_methods.cpp +++ b/tdms/src/interpolation_methods.cpp @@ -237,41 +237,32 @@ int InterpolationScheme::num_nonzero_coeffs() const { return last_nonzero_coeff - first_nonzero_coeff + 1; } -double InterpolationScheme::interpolate(const double *v, const int offset) const { - - double interp_value = 0.; - for(int ind=first_nonzero_coeff; ind<=last_nonzero_coeff; ind++) { - interp_value += scheme_coeffs[ind] * v[ind+offset]; - } - return interp_value; -} - bool InterpolationScheme::is_better_than(const InterpolationScheme s) const { return (priority > s.get_priority()); } const InterpolationScheme &best_scheme(int cells_in_direction, int cell_id) { - // interpolation is impossible with fewer than 4 cells in a dimension - if (cells_in_direction < 4) { - throw out_of_range("Error: computational domain has fewer than 4 cells in at least 1 dimension, cubic and bandlimited interpolation impossible.\n"); + // interpolation is impossible with fewer than 4 datapoints in a dimension + // 4 datapoints => 3 cells in a given direction + if (cells_in_direction < 3) { + throw out_of_range("Error: domain axis has < 3 cells (< 4 datapoints), cubic and bandlimited interpolation impossible.\n"); } - // Yee cell with index <0 does exist (starts from 0) - // Current code functionality is to error if we attempt to interpolate to cell 0, because the cubic interpolation found this impossible - // However, BL_TO_CELL_0 can handle this, provided cells_in_direction >= 8. - else if (cell_id <= 0) { - throw out_of_range("Error: Yee cell index <=0 requested (must be >=1).\n"); + // Yee cell with index <0 doesn't allow for interpolation (this is the PML) + else if (cell_id < 0) { + throw out_of_range("Error: Yee cell index < 0 requested.\n"); } // Yee cell with index >= cells_in_direction doesn't exist // Again, we can in theory determine the value "here" using BL7, however cubic interpolation cannot do this, and so current code functionality is to throw an error here. else if (cell_id >= cells_in_direction) { throw out_of_range("Error: Yee cell index beyond maximum number of Yee cells requested.\n"); } - else if (cells_in_direction < 8) { + else if (cells_in_direction < 7) { // we do not have enough cells to use bandlimited interpolation, but can use cubic - // by definition, cell_id = 1 requires us to use CBFst, cell_id = cells_in_direction-1 CBLast, - // and everything else CBMid - if (cell_id == 1) { + // cell_id = 0 requires us to use CBFst, + // cell_id = 1,2,...,cells_in_direction-2 requires us to use CBMid, + // cell_id = cells_in_direction-1 requires us to use CBLast + if (cell_id == 0) { return CBFst; } else if (cell_id == cells_in_direction-1) { @@ -283,17 +274,17 @@ const InterpolationScheme &best_scheme(int cells_in_direction, int cell_id) { } else { // we can apply bandlimited interpolation. - // unless we are <=3 cells away from either boundary, we will want to use BL3 - if ((cell_id >= 4) && (cell_id <= cells_in_direction - 4)) { + // if there are 4 datapoints either side, we want to use BL3 + if ((cell_id > 2) && (cell_id < cells_in_direction - 3)) { return BL3; } - else if (cell_id == 1) { + else if (cell_id == 0) { return BL0; } - else if (cell_id == 2) { + else if (cell_id == 1) { return BL1; } - else if (cell_id == 3) { + else if (cell_id == 2) { return BL2; } else if (cell_id == cells_in_direction - 3) { diff --git a/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_BLi_vs_cubic_validation.m b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_BLi_vs_cubic_validation.m new file mode 100644 index 000000000..28e03bc48 --- /dev/null +++ b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_BLi_vs_cubic_validation.m @@ -0,0 +1,90 @@ +%{ +BENCHMARK ERROR VALUES FOR BAND-LIMITED INTERPOLATION AGAINST CUBIC +INTERPOLATION + +Determines the error in interpolated values when using MATLAB's interp +function and cubic interpolation methods, to be used as the benchmark for +test_BLi_vs_cubic_interpolation.cpp. + +The function we use is a Cauchy distribution; f(x) = 1/(10x^2+1), over the +range x\in[-2,2] + +Neither cubic nor BLi should be the superior interpolation method for all +datapoint-separation distances. + +Results are printed out for each cell size to stdout. +%} +clear; +close all; + +% interp inputs +r = 2; +N = 4; +% the spacing between datapoints, mimicing the Yee cell dimensions +cell_Sizes = [0.25, 0.1, 0.05, 0.01]; + +% for each trial compute BLi and cubic errors +for trial=1:4 + cellSize = cell_Sizes(trial); + n_datapts = ceil(4/cellSize); % number of Yee cells + x_lower = -2.; % lower value of the range over which we shall test + + cell_centres = zeros(1,n_datapts-1); % interpolation positions (Yee cell centres) + field_pos = zeros(1,n_datapts); % data sample positions (field component associated to Yee cells) + for i=1:n_datapts + if i~=n_datapts + cell_centres(i) = x_lower + (i-0.5)*cellSize; + end + field_pos(i) = x_lower + (i-1)*cellSize; + end + field_samples = BLi_vs_cubic_fn(field_pos); % compute data samples + exact_vals = BLi_vs_cubic_fn(cell_centres); % compute exact values + + BLi_interp_full = interp(field_samples, r, N); % perform BLi interpolation + BLi_interp = BLi_interp_full(2:2:end-1); % final point is beyond last Yee cell + + % perform cubic interpolation + cubic_interp = zeros(1,n_datapts-1); + cubic_interp(1) = interp_cubic(field_samples(1:4),-1); + for i=2:n_datapts-2 + cubic_interp(i) = interp_cubic(field_samples(i-1:i+2),0); + end + cubic_interp(n_datapts-1) = interp_cubic(field_samples(end-3:end),1); + + % compute errors + BLi_err = abs( BLi_interp - exact_vals ); + BLi_err_max = max(BLi_err); + BLi_err_norm = norm( BLi_interp - exact_vals ); + cub_err = abs( cubic_interp - exact_vals ); + cub_err_max = max(cub_err); + cub_err_norm = norm( cubic_interp - exact_vals ); + + BLi_was_better = ones(size(BLi_err)); + BLi_was_better( cub_err < BLi_err ) = 0.; + + fprintf(" ====== \n Cellsize %.3f : \n", cellSize); + fprintf("MAX Errors: \n"); + fprintf("Cubic: \t %.8e\n", cub_err_max); + fprintf("BLi: \t %.8e\n", BLi_err_max); + fprintf("NORM Errors: \n"); + fprintf("Cubic: \t %.8e\n", cub_err_norm); + fprintf("BLi: \t %.8e\n", BLi_err_norm); + fprintf("Pointwise, BLi was better at %d points.\n", sum(BLi_was_better)); +end + +function [out] = interp_cubic(v, pos) + if pos==0 + % midpoint interpolate + out = -(1/16) * v(1) + (9/16)*v(2) + (9/16)*v(3) - (1/16)*v(4); + elseif pos==-1 + % left interp + out = (5/16)*v(1) + (15/16)*v(2) - (5/16)*v(3) + (1/16)*v(4); + elseif pos==1 + % right interp + out = (1/16)*v(1) - (5/16)*v(2) + (15/16)*v(3) + (5/16)*v(4); + end +end + +function [out] = BLi_vs_cubic_fn(x) + out = 1 ./ (10*x.^2 + 1); +end diff --git a/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_E.m b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_E.m new file mode 100644 index 000000000..54f2c6792 --- /dev/null +++ b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_E.m @@ -0,0 +1,170 @@ +%{ + BENCHMARK FOR E-FIELD INTERPOLATION METHODS + +Tests the interpolation accuracy for the E-field components, when +interpolated in each direction to the centre of a Yee cell - benchmarking +test_field_interpolation::TEST_CASE("E-field interpolation check"). + +The E-field used will have identical field components, each of which are +given by the expression: + +E_t(tt) = sin(2\pi tt) exp(-tt^2), + +and will be tested over the range $[-2,4]^3$. + +Errors are computed as the Frobenius norm of the pointwise difference +(between the exact and interpolated values) matrix, and as the maximum +error in the interpolated values along each axial direction (slice-norm +metric). + +Results are displayed to stdout for each component; Ex, Ey, Ez. +%} +clear; +close all; + +% interp parameters +r = 2; +N = 4; +% bad practice but saves us passing dimensions to subfunctions +% with clear and close all, this should be safe +global cellDims; + +% (1,2,3) = (Dx, Dy, Dz) +cellDims = [0.25, 0.1, 0.05]; +% domain range is [-2,4] +x_lower = -2.; extent_x = 4; +y_lower = -2.; extent_y = 4; +z_lower = -2.; extent_z = 4; +% function handle with these _lower values +GCC = @(x,y,z) GetCellCentre(x,y,z,x_lower,y_lower,z_lower); +FP = @(x,y,z,c) GetFieldPos(x,y,z,x_lower,y_lower,z_lower,c); + +% number of Yee cells in each dimension +nX = round(extent_x/cellDims(1)); +nY = round(extent_y/cellDims(2)); +nZ = round(extent_z/cellDims(3)); + +% components of the field at the Yee cell centres, +% and at the field-sample positions +Ex_exact = zeros(nX, nY, nZ); Ex_sample = zeros(nX, nY, nZ); +Ey_exact = zeros(nX, nY, nZ); Ey_sample = zeros(nX, nY, nZ); +Ez_exact = zeros(nX, nY, nZ); Ez_sample = zeros(nX, nY, nZ); +% populate exact values and sample datapoints +for i=1:nX + for j=1:nY + for k=1:nZ + centre = GCC(i,j,k); + fp_ex = FP(i,j,k,"Ex"); + fp_ey = FP(i,j,k,"Ey"); + fp_ez = FP(i,j,k,"Ez"); + Ex_exact(i,j,k) = Ex_field(centre(1), centre(2), centre(3)); + Ey_exact(i,j,k) = Ey_field(centre(1), centre(2), centre(3)); + Ez_exact(i,j,k) = Ez_field(centre(1), centre(2), centre(3)); + Ex_sample(i,j,k) = Ex_field(fp_ex(1), fp_ex(2), fp_ex(3)); + Ey_sample(i,j,k) = Ey_field(fp_ey(1), fp_ey(2), fp_ey(3)); + Ez_sample(i,j,k) = Ez_field(fp_ez(1), fp_ez(2), fp_ez(3)); + end + end +end +% NOTE: we don't use the field values at cells with index 0 +% so we'll only use E{x,y,z}_exact(2:end,2:end,2:end) from here on out + +% interpolate Ex in the x direction to find the values at the cell centres +Ex_interp = zeros(nX-1, nY, nZ); +slice_norms = zeros(nY, nZ); +for k=1:nZ + for j=1:nY + temp_x = interp(reshape(Ex_sample(:,j,k),1,nX),r,N); + % final point is beyond final Yee cell centre + Ex_interp(:,j,k) = temp_x(2:2:end-1); + + % manual check along this axis of the square-norm error + slice_norms(j,k) = norm(Ex_interp(:,j,k) - Ex_exact(1:end-1,j,k)); + end +end +Ex_errs = Ex_interp - Ex_exact(2:end,:,:); +Ex_max_slice_norm = max( slice_norms, [], "all" ); +Ex_norm_err = norm( Ex_errs, "fro" ); + +% interpolate Ey in the y direction to find the values at the cell centres +Ey_interp = zeros(nX, nY-1, nZ); +slice_norms = zeros(nX, nZ); +for k=1:nZ + for i=1:nX + temp_y = interp(reshape(Ey_sample(i,:,k),1,nY),r,N); + % final point is beyond final Yee cell centre + Ey_interp(i,:,k) = temp_y(2:2:end-1); + + % manual check along this axis of the square-norm error + slice_norms(i,k) = norm(Ey_interp(i,:,k) - Ey_exact(i,1:end-1,k)); + end +end +Ey_errs = Ey_interp - Ey_exact(:,2:end,:); +Ey_max_slice_norm = max( slice_norms, [], "all" ); +Ey_norm_err = norm( Ey_errs, "fro" ); + +% interpolate Ez in the z direction to find the values at the cell centres +Ez_interp = zeros(nX, nY, nZ-1); +slice_norms = zeros(nX, nY); +for i=1:nX + for j=1:nY + temp_z = interp(reshape(Ez_sample(i,j,:),1,nZ),r,N); + % final point is beyond final Yee cell centre + Ez_interp(i,j,:) = temp_z(2:2:end-1); + + % manual check along this axis of the square-norm error + slice_norms(i,j) = norm(reshape(Ez_interp(i,j,:) - Ez_exact(i,j,1:end-1), 1, nZ-1)); + end +end +Ez_errs = Ez_interp - Ez_exact(:,:,2:end); +Ez_max_slice_norm = max( slice_norms, [], "all" ); +Ez_norm_err = norm( Ez_errs, "fro" ); + +fprintf("Ex Frobenius norm error: \t %.16e \n", Ex_norm_err); +fprintf("Ex max-slice norm error: \t %.16e \n", Ex_max_slice_norm); +fprintf("Ey Frobenius norm error: \t %.16e \n", Ey_norm_err); +fprintf("Ey max-slice norm error: \t %.16e \n", Ey_max_slice_norm); +fprintf("Ez Frobenius norm error: \t %.16e \n", Ez_norm_err); +fprintf("Ez max-slice norm error: \t %.16e \n", Ez_max_slice_norm); + +%% E-field component functions +% Field components E_t(tt) = sin(2\pi tt) exp(-tt^2) +function [value] = Ex_field(~,y,~) + value = sin(2.*pi*y) * exp(-y.^2); +end +function [value] = Ey_field(~,~,z) + value = sin(2.*pi*z) * exp(-z.^2); +end +function [value] = Ez_field(x,~,~) + value = sin(2.*pi*x) * exp(-x.^2); +end + +%% Computes the coordinates of the centre of a Yee cell +function [centre] = GetCellCentre(i,j,k,x_lower,y_lower,z_lower) + global cellDims; + centre = [x_lower + (i-0.5)*cellDims(1), y_lower + (j-0.5)*cellDims(2), z_lower + (k-0.5)*cellDims(3)]; +end + +%% Computes the position of the field component associated to a Yee cell +function [fieldPos] = GetFieldPos(i,j,k,x_lower,y_lower,z_lower,component) + global cellDims; + fieldPos = GetCellCentre(i,j,k,x_lower,y_lower,z_lower); + if component=="Ex" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + elseif component=="Ey" + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + elseif component=="Ez" + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hx" + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hy" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hz" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + else + error("Not a recognised component") + end +end diff --git a/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_H.m b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_H.m new file mode 100644 index 000000000..77c1aa18a --- /dev/null +++ b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_field_interpolation_H.m @@ -0,0 +1,199 @@ +%{ + BENCHMARK FOR H-FIELD INTERPOLATION METHODS + +Tests the interpolation accuracy for the H-field components, when +interpolated in each direction to the centre of a Yee cell - benchmarking +test_field_interpolation::TEST_CASE("H-field interpolation check"). + +The H-field used will have identical field components, each of which are +given by the expression: + +H_t(tt) = sin(2\pi tt) exp(-tt^2), + +and will be tested over the range $[-2,4]^3$. + +Errors are computed as the Frobenius norm of the pointwise difference +(between the exact and interpolated values) matrix. Since we need to slice +in two dimensions for the H-field components, we do not test the slice-norm +metric. + +Results are displayed to stdout for each component; Hx, Hy, Hz. +%} +clear; +close all; + +% interp parameters +r = 2; +N = 4; +% bad practice but saves us passing dimensions to subfunctions +% with clear and close all, this should be safe +global cellDims; + +% (1,2,3) = (Dx, Dy, Dz) +cellDims = [0.1, 0.05, 0.025]; +x_lower = -2.; extent_x = 4; +y_lower = -2.; extent_y = 4; +z_lower = -2.; extent_z = 4; +% function handle with these _lower values +GCC = @(x,y,z) GetCellCentre(x,y,z,x_lower,y_lower,z_lower); +FP = @(x,y,z,c) GetFieldPos(x,y,z,x_lower,y_lower,z_lower,c); + +% number of Yee cells in each dimension +nX = ceil(extent_x/cellDims(1)); +nY = ceil(extent_y/cellDims(2)); +nZ = ceil(extent_z/cellDims(3)); + +% components of the field at the Yee cell centres, +% and at the field-sample positions +Hx_exact = zeros(nX, nY, nZ); Hx_sample = zeros(nX, nY, nZ); +Hy_exact = zeros(nX, nY, nZ); Hy_sample = zeros(nX, nY, nZ); +Hz_exact = zeros(nX, nY, nZ); Hz_sample = zeros(nX, nY, nZ); +% populate exact values and sample datapoints +for i=1:nX + for j=1:nY + for k=1:nZ + centre = GCC(i,j,k); + fp_ex = FP(i,j,k,"Hx"); + fp_ey = FP(i,j,k,"Hy"); + fp_ez = FP(i,j,k,"Hz"); + Hx_exact(i,j,k) = Hx_field(centre(1), centre(2), centre(3)); + Hy_exact(i,j,k) = Hy_field(centre(1), centre(2), centre(3)); + Hz_exact(i,j,k) = Hz_field(centre(1), centre(2), centre(3)); + Hx_sample(i,j,k) = Hx_field(fp_ex(1), fp_ex(2), fp_ex(3)); + Hy_sample(i,j,k) = Hy_field(fp_ey(1), fp_ey(2), fp_ey(3)); + Hz_sample(i,j,k) = Hz_field(fp_ez(1), fp_ez(2), fp_ez(3)); + end + end +end +% NOTE: we don't use the field values at cells with index 0 +% so we'll only use H{x,y,z}_exact(2:end,2:end,2:end) from here on out + +% interpolate Hx... do so in y then z, then z then y, and take the worse +% result? +Hx_ypullback = zeros(nX, nY-1, nZ); +Hx_zpullback = zeros(nX, nY, nZ-1); +Hx_interp_yz = zeros(nX, nY-1, nZ-1); +Hx_interp_zy = zeros(nX, nY-1, nZ-1); +for i=1:nX + for j=1:nY + temp_z = interp(reshape( Hx_sample(i,j,:),1,nZ),r,N); + Hx_zpullback(i,j,:) = temp_z(2:2:end-1); + end + for k=1:nZ + temp_y = interp(reshape( Hx_sample(i,:,k),1,nY),r,N); + Hx_ypullback(i,:,k) = temp_y(2:2:end-1); + end + for j=1:nY-1 + temp_z = interp(reshape( Hx_ypullback(i,j,:),1,nZ),r,N); + Hx_interp_yz(i,j,:) = temp_z(2:2:end-1); + end + for k=1:nZ-1 + temp_y = interp(reshape( Hx_zpullback(i,:,k),1,nY),r,N); + Hx_interp_zy(i,:,k) = temp_y(2:2:end-1); + end +end +% take worst-case scenario, in the C++ scheme chooses to use one dimension +% over another arbitrarily +Hx_errors = max(abs(Hx_interp_zy - Hx_exact(:,1:end-1,1:end-1)), abs(Hx_interp_yz - Hx_exact(:,1:end-1,1:end-1))); +Hx_fro_err = norm( Hx_errors, "fro" ); +fprintf("Hx Frobenius norm error: \t %.16e \n", Hx_fro_err); + +% interpolate Hy... do so in x then z, then z then x, and take the worse +% result? +Hy_xpullback = zeros(nX-1, nY, nZ); +Hy_zpullback = zeros(nX, nY, nZ-1); +Hy_interp_xz = zeros(nX-1, nY, nZ-1); +Hy_interp_zx = zeros(nX-1, nY, nZ-1); +for j=1:nY + for i=1:nX + temp_z = interp(reshape( Hy_sample(i,j,:),1,nZ),r,N); + Hy_zpullback(i,j,:) = temp_z(2:2:end-1); + end + for k=1:nZ + temp_x = interp(reshape( Hy_sample(:,j,k),1,nX),r,N); + Hy_xpullback(:,j,k) = temp_x(2:2:end-1); + end + for i=1:nX-1 + temp_z = interp(reshape( Hy_xpullback(i,j,:),1,nZ),r,N); + Hy_interp_xz(i,j,:) = temp_z(2:2:end-1); + end + for k=1:nZ-1 + temp_x = interp(reshape( Hy_zpullback(:,j,k),1,nX),r,N); + Hy_interp_zx(:,j,k) = temp_x(2:2:end-1); + end +end +% take worst-case scenario, in the C++ scheme chooses to use one dimension +% over another arbitrarily +Hy_errors = max(abs(Hy_interp_zx - Hy_exact(1:end-1,:,1:end-1)), abs(Hy_interp_xz - Hy_exact(1:end-1,:,1:end-1))); +Hy_fro_err = norm( Hy_errors, "fro" ); +fprintf("Hy Frobenius norm error: \t %.16e \n", Hy_fro_err); + +% interpolate Hz... do so in x then y, then y then x, and take the worse +% result? +Hz_xpullback = zeros(nX-1, nY, nZ); +Hz_ypullback = zeros(nX, nY-1, nZ); +Hz_interp_xy = zeros(nX-1, nY-1, nZ); +Hz_interp_yx = zeros(nX-1, nY-1, nZ); +for k=1:nZ + for j=1:nY + temp_x = interp(reshape( Hz_sample(:,j,k),1,nX),r,N); + Hz_xpullback(:,j,k) = temp_x(2:2:end-1); + end + for i=1:nX + temp_y = interp(reshape( Hz_sample(i,:,k),1,nY),r,N); + Hz_ypullback(i,:,k) = temp_y(2:2:end-1); + end + for j=1:nY-1 + temp_x = interp(reshape( Hz_ypullback(:,j,k),1,nX),r,N); + Hz_interp_yx(:,j,k) = temp_x(2:2:end-1); + for i=1:nX-1 + temp_y = interp(reshape( Hz_xpullback(i,:,k),1,nY),r,N); + Hz_interp_xy(i,:,k) = temp_y(2:2:end-1); + end + end +end +Hz_errors = max(abs(Hz_interp_yx - Hz_exact(1:end-1,1:end-1,:)), abs(Hz_interp_xy - Hz_exact(1:end-1,1:end-1,:))); +Hz_fro_err = norm( Hz_errors, "fro" ); +fprintf("Hz Frobenius norm error: \t %.16e \n", Hz_fro_err); + +%% H-field component function +% Field components H_t(tt) = sin(2\pi tt) exp(-tt^2) +function [value] = Hx_field(~,y,~) + value = sin(2. * pi * y) * exp(-y.^2); +end +function [value] = Hy_field(~,~,z) + value = sin(2. * pi * z) * exp(-z.^2); +end +function [value] = Hz_field(x,~,~) + value = sin(2. * pi * x) * exp(-x.^2); +end + +%% Computes the coordinates of the centre of a Yee cell +function [centre] = GetCellCentre(i,j,k,x_lower,y_lower,z_lower) + global cellDims; + centre = [x_lower + (i-0.5)*cellDims(1), y_lower + (j-0.5)*cellDims(2), z_lower + (k-0.5)*cellDims(3)]; +end + +%% Computes the position of the field component associated to a Yee cell +function [fieldPos] = GetFieldPos(i,j,k,x_lower,y_lower,z_lower,component) + global cellDims; + fieldPos = GetCellCentre(i,j,k,x_lower,y_lower,z_lower); + if component=="Ex" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + elseif component=="Ey" + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + elseif component=="Ez" + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hx" + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hy" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + fieldPos(3) = fieldPos(3) - 0.5*cellDims(3); + elseif component=="Hz" + fieldPos(1) = fieldPos(1) - 0.5*cellDims(1); + fieldPos(2) = fieldPos(2) - 0.5*cellDims(2); + else + error("Not a recognised component") + end +end diff --git a/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_interpolation_functions.m b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_interpolation_functions.m new file mode 100644 index 000000000..f31afd439 --- /dev/null +++ b/tdms/tests/unit/matlab_benchmark_scripts/benchmark_test_interpolation_functions.m @@ -0,0 +1,72 @@ +%{ + BENCHMARK ERROR VALUES FOR BAND-LIMITED INTERPOLATION + +Determines the error in MATLAB's BLi scheme for the test functions used in +the BLi unit tests in test_interpolation_functions.cpp: +- The constant function +- sin(2\pi x) +- A pulse function +- A complex-valued function of the form sin(2\pi x) + pulse(x)*i + +These errors are printed to the screen upon execution of the script. +%} + +clear; +close all; + +nSamples = 100; % number of sample points to use +r = 2; % interpolate to midpoints of samples +N = 4; % use 2*4=8 sample points in interpolation +x = linspace(0,1,r*nSamples); % sample points + +% Constant function interpolation +const_fn_data = const_fn(x); +const_fn_interp = interp(const_fn_data(1:r:end),r,N); +const_fn_err = max(abs( const_fn_interp(2:r:end-1) - const_fn_data(2:r:end-1) )); + +% sin(2\pi x) interpolation +s2pi_data = s2pi(x); +s2pi_interp = interp(s2pi_data(1:r:end),r,N); +s2pi_err = max(abs( s2pi_interp(2:r:end-1) - s2pi_data(2:r:end-1) )); + +% pulse function interpolation +pulse_data = pulse_fn(x); +pulse_interp = interp(pulse_data(1:r:end),r,N); +pulse_err = max(abs( pulse_interp(2:r:end-1) - pulse_data(2:r:end-1) )); + +% complex: s2pi(x) + i*pulse_fn(x) interpolation +complex_data = complex_test(x); +complex_interp = interp(complex_data(1:r:end),r,N); +complex_err = max(abs( complex_interp(2:r:end-1) - complex_data(2:r:end-1) )); + +fprintf("Constant function error: %.8e \n", const_fn_err); +fprintf("sin function error: %.8e \n", s2pi_err); +fprintf("pulse function error: %.8e \n", pulse_err); +fprintf("complex fn error: %.8e \n", complex_err); + +%% constant function +function [out] = const_fn(x) + out = 1. + x*0.; +end + +%% sin(2\pi x) +function [out] = s2pi(x) + out = sin(2*pi*x); +end + +%% compact pulse function +% f(x) = exp( -1 / (1-3|2x-1|)^2 ) when 3|2x-1|<1, +% f(x) = 0 otherwise. +function [out] = pulse_fn(x) + xhat = 3 * (2*x - 1); + out = zeros(size(xhat)); + out(abs(xhat)>=1) = 0.; + out(abs(xhat)<1) = exp( -1 ./ ( 1 - abs( xhat(abs(xhat)<1) ).^2 )); +end + +%% complex-valued function +% Real part is s2pi(x) +% Imag part is pulse_fn(x) +function [out] = complex_test(x) + out = s2pi(x) + pulse_fn(x)*1.i; +end diff --git a/tdms/tests/unit/matlab_benchmark_scripts/readme.md b/tdms/tests/unit/matlab_benchmark_scripts/readme.md new file mode 100644 index 000000000..d4801d3d4 --- /dev/null +++ b/tdms/tests/unit/matlab_benchmark_scripts/readme.md @@ -0,0 +1,8 @@ +# **`MATLAB` Benchmarking Scripts** + +The directory `matlab_benchmark_scripts` contains scripts which perform band-limited interpolation (BLi) using `MATLAB`'s `interp` function. +`TDMS`'s interpolation schemes are based off this `MATLAB` function (specficially, in the coefficients the scheme uses to interpolate). + +In order to test that the interpolation is correctly implimented in the `TDMS` source, we provide unit tests that benchmark against `MATLAB`'s implimentations. These scripts are provided here for developer use and documentation. + +Currently; the unit tests have the output error values from `MATLAB` hard-coded into the source, and are required to achieve a comparable level of accuracy for the tests to pass. Each script contains a note referencing which unit test it provides a benchmark for. \ No newline at end of file diff --git a/tdms/tests/unit/test_BLi_vs_cubic_interpolation.cpp b/tdms/tests/unit/test_BLi_vs_cubic_interpolation.cpp new file mode 100644 index 000000000..0d6c1987c --- /dev/null +++ b/tdms/tests/unit/test_BLi_vs_cubic_interpolation.cpp @@ -0,0 +1,143 @@ +/** + * @file test_BLi_vs_cubic_interpolation.cpp + * @author William Graham (ccaegra@ucl.ac.uk) + * @brief Tests the performance of band-limited interpolation against cubic interpolation + */ +#include +#include +#include +#include + +#include "interpolation_methods.h" + +using namespace std; + +// Computes the 2-norm of the vector v from buffer start to buffer end +inline double norm(double *v, int end, int start = 0) { + double norm_val = 0.; + for (int i = start; i <= end; i++) { + norm_val += v[i] * v[i]; + } + return sqrt(norm_val); +} + +// function to interpolate +inline double f_BLi_vs_Cubic(double x) { + return 1. / ((10. * x * x) + 1.); +} + +// computes order of magnitude of a double +int orderOfMagnitude(double x) { + return floor(log10(x)); +} + +/** + * @brief We will check that BLi gives a better approximation than cubic interpolation. + * + * The metric we will use is the norm of the error-vector. + * As the cell-size becomes larger (across the same domain), BLi should begin to perform better than cubic interpolation. + * + * For the following cell sizes, over the range -2 to 2, and for the function f(x) = 1/(10x^2+1), MATLAB predicts the following errors: + * Cell size | BLi err | Cubic err + * 0.25 | 1.25953429e-02 | 3.28105674e-02 + * 0.1 | 7.10587711e-04 | 5.49601748e-03 + * 0.05 | 7.80325615e-04 | 5.53430587e-04 + * 0.01 | 1.97562511e-03 | 2.06917814e-06 + * Values computed via the benchmark_test_BLi_vs_cubic.m script. + * + * We will test the following: + * - The norm-error of BLi is lesser/greater than cubic in each case + * - The order of the magnitude (of BOTH) errors is the same as that obtained from MATLAB (for validation reasons) + */ +TEST_CASE("Benchmark: BLi is better than cubic interpolation") { + SPDLOG_INFO("===== Benchmarking BLi against cubic interpolation ====="); + // number of cell sizes to test at + const int n_trials = 4; + // Yee cell sizes to test across + double cell_sizes[n_trials] = {0.25, 0.1, 0.05, 0.01}; + // Flags for whether or not we expect BLi to perform better + bool BLi_is_better[n_trials] = {true, true, false, false}; + // MATLAB errors for BLi + double MATLAB_BLi_errs[n_trials] = {1.25953429e-02, 7.10587711e-04, 7.80325615e-04, + 1.97562511e-03}; + // MATLAB errors for cubic interp + double MATLAB_cub_errs[n_trials] = {3.28105674e-02, 5.49601748e-03, 5.53430587e-04, + 2.06917814e-06}; + // coordinate of the LHS of the first Yee cell + double x_lower = -2.; + // spatial extent of the domain will be [x_lower, x_lower+extent] + double extent = 4.; + // norm errors recorded across runs + double BLi_norm_errs[n_trials], cub_norm_errs[n_trials]; + + for (int trial = 0; trial < n_trials; trial++) { + // Yee cell "size" to use + double cellSize = cell_sizes[trial]; + // the number of datapoints that we have, the number of "cells" we have is one less than this + int n_datapts = round(extent / cellSize); + CHECK(n_datapts > 8); // BLi cannot run otherwise + + // the centres of the Yee cells + double cell_centres[n_datapts - 1]; + // the exact field values at the Yee cell centres + double exact_field_values[n_datapts - 1]; + // the spatial coordinates of the samples of our field + double field_positions[n_datapts]; + // the field values at the sample positions + double field_samples[n_datapts]; + for (int i = 0; i < n_datapts; i++) { + if (i != n_datapts - 1) { cell_centres[i] = x_lower + (((double) i) + 0.5) * cellSize; } + field_positions[i] = x_lower + ((double) i) * cellSize; + + exact_field_values[i] = f_BLi_vs_Cubic(cell_centres[i]); + field_samples[i] = f_BLi_vs_Cubic(field_positions[i]); + } + + // the interpolated field values via BLi + double BLi_interp[n_datapts - 1]; + // the interpolated field values via cubic interp + double cub_interp[n_datapts - 1]; + + // pointwise-interpolation error in BLi + double BLi_err[n_datapts - 1]; + // pointwise-interpolation error in cubic + double cub_err[n_datapts - 1]; + + // perform interpolation - this is manual since best_interp_scheme will always want to do BLi + // cubic interpolation only changes at the first and last cells + cub_interp[0] = CBFst.interpolate(field_samples); + for (int i = 1; i < n_datapts - 2; i++) { + cub_interp[i] = CBMid.interpolate(field_samples, i + 1 - CBMid.number_of_datapoints_to_left); + } + cub_interp[n_datapts - 2] = CBLst.interpolate(field_samples, n_datapts - 1 - CBLst.number_of_datapoints_to_left); + + // BLi interpolation + for (int i = 0; i < n_datapts - 1; i++) { + // we checked earlier that BLi should always be available, so this is always a BLi scheme + InterpolationScheme use_scheme = best_scheme(n_datapts - 1, i); + // to be on the safe side, check that we have a BLi scheme. BLi schemes are always strictly better than CUBIC_INTERP_MIDDLE, and not equal or worse. + CHECK(use_scheme.is_better_than(CUBIC_INTERP_MIDDLE)); + // interpolate using the cell data we have provided + // i + 1 - use_scheme.number_of_datapoints_to_left is the index of field_samples to read as the point v[0] in the schemes + BLi_interp[i] = use_scheme.interpolate(field_samples, i + 1 - use_scheme.number_of_datapoints_to_left); + } + + // compare to exact values - again ignore cell 0 for the time being + for (int i = 0; i < n_datapts - 1; i++) { + BLi_err[i] = abs(BLi_interp[i] - exact_field_values[i]); + cub_err[i] = abs(cub_interp[i] - exact_field_values[i]); + } + // compute square-norm error + BLi_norm_errs[trial] = norm(BLi_err, n_datapts - 2); + cub_norm_errs[trial] = norm(cub_err, n_datapts - 2); + + // value-testing commences: the better method should be superior, and the worse method should be worse. + // assert this is the case for each cellSize we used. + CHECK((BLi_norm_errs[trial] <= cub_norm_errs[trial]) == BLi_is_better[trial]); + SPDLOG_INFO("BLi err: {0:.8e} | Cubic err : {1:.8e} | Expected BLi to be better: {2:d}", BLi_norm_errs[trial], cub_norm_errs[trial], BLi_is_better[trial]); + + // in all cases, assert that the order of magnitude of error is what we expect, if we had used MATLAB + CHECK(orderOfMagnitude(BLi_norm_errs[trial]) == orderOfMagnitude(MATLAB_BLi_errs[trial])); + CHECK(orderOfMagnitude(cub_norm_errs[trial]) == orderOfMagnitude(MATLAB_cub_errs[trial])); + } +} diff --git a/tdms/tests/unit/test_field_interpolation.cpp b/tdms/tests/unit/test_field_interpolation.cpp new file mode 100644 index 000000000..0ccc574c3 --- /dev/null +++ b/tdms/tests/unit/test_field_interpolation.cpp @@ -0,0 +1,405 @@ +/** + * @file test_field_interpolation.cpp + * @author William Graham (ccaegra@ucl.ac.uk) + * @brief Tests interpolation of E- and H-fields and compares the errors against MATLAB benchmarks + */ +#include +#include +#include + +#include "globals.h" +#include "field.h" + +using namespace tdms_math_constants; + +/* Overview of Tests + +We will test the performance of BLi at interpolating a known field to the centre of each Yee cell, for both the E- and H- fields. +The benchmark for success will be superior or equivalent error (Frobenius and 1D-dimension-wise norm) to that produced by MATLAB performing the same functionality. +Frobenious error is the Frobenius norm of the 3D matrix whose (i,j,k)-th entry is the error in the interpolated value at Yee cell centre (i,j,k). +The slice error, for a fixed j,k, is the norm-error of the 1D array of interpolated values along the axis (:,j,k). +The maximum of this is then the max-slice error: keeping track of this ensures us that the behaviour of BLi is consistent, and does not dramaically over-compensate in some areas and under-compensate in others. + +All tests will be performed with cell sizes Dx = 0.25, Dy = 0.1, Dz = 0.05, over the range [-2,2]. +*/ + +// computes the Euclidean norm of a 1d-array (or pointer thereto) +inline double euclidean(double *v, int end, int start = 0) { + double norm_val = 0.; + for (int i = start; i < end; i++) { + norm_val += v[i] * v[i]; + } + return sqrt(norm_val); +} + +// functional form for the {E,H}-field components +inline double field_component(double t) { + return sin(2. * DCPI * t) * exp(-t * t); +} + +/** + * @brief Test the interpolation of the E-field components to the centre of the Yee cells + * + * Each component of the E-field will take the form + * E_t(tt) = sin(2\pi tt) * exp(-tt^2) + * + * We test both the Fro- and slice-norm metrics, since interpolation only happens along one axis + */ +TEST_CASE("E-field interpolation check") { + SPDLOG_INFO("===== Testing E-field BLi ====="); + // error tolerance, based on MATLAB performance + // script: benchmark_test_field_interpolation_H.m + double Ex_fro_tol = 2.4535659128911650e-02, Ex_ms_tol = 1.0294466335592440e-03; + double Ey_fro_tol = 2.5021563893394754e-02, Ey_ms_tol = 1.6590813667643383e-03; + double Ez_fro_tol = 2.5181324617226587e-02, Ez_ms_tol = 1.7507103927894884e-03; + + // additional tolerance to allow for floating-point rounding imprecisions, etc + double acc_tol = 1e-12; + + // fake domain setup + double x_lower = -2., y_lower = -2., z_lower = -2.; + double extent_x = 4., extent_y = 4., extent_z = 4.; + double cellDims[3] = {0.25, 0.1, 0.05}; + // The number of cells in each direction is then 16 = 4/0.25, 40 = 4/0.1, 80 = 4/0.05. + // Note that due to the possibility that Nx/cellDims[0] computing something that is not quite an integer, we need to use round() to get an int safely + int Nx = round(extent_x / cellDims[0]), Ny = round(extent_y / cellDims[1]), + Nz = round(extent_z / cellDims[2]); + SPDLOG_INFO("(Nx, Ny, Nz) = ({},{},{})", Nx, Ny, Nz); + + // setup the "split" E-field components + ElectricSplitField E_split(Nx, Ny, Nz); + E_split.allocate(); + // setup for non-split field components + ElectricField E(Nx, Ny, Nz); + E.allocate(); + // storage for pointwise errors + Tensor3D Ex_error, Ex_split_error, Ey_error, Ey_split_error, Ez_error, Ez_split_error; + Ex_error.allocate(Nz, Ny, Nx - 1); + Ex_split_error.allocate(Nz, Ny, Nx - 1); + Ey_error.allocate(Nz, Ny - 1, Nx); + Ey_split_error.allocate(Nz, Ny - 1, Nx); + Ez_error.allocate(Nz - 1, Ny, Nx); + Ez_split_error.allocate(Nz - 1, Ny, Nx); + + // compute the exact field and the "split field" components + // indices range from (0,0,0) to (Nx, Ny, Nz) INCLUSIVE! + for (int ii = 0; ii <= Nx; ii++) { + for (int jj = 0; jj <= Ny; jj++) { + for (int kk = 0; kk <= Nz; kk++) { + /* The coordinates of the grids are different, which throws a spanner in the works. + Recall that the "x-grid" is at the position of the (i,j,k)-th Ex field sample. + x-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + (k-0.5)*cellDims[2]) + y-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + j*cellDims[1], z_lower + (k-0.5)*cellDims[2]) + z-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + k*cellDims[2]) + The components only depend on one of the coordinate directions, so we don't need to compute all 9 values. + */ + double x_comp_value = + field_component(y_lower + (((double) jj) + 0.5) * cellDims[1]);// Ex depends on y + double y_comp_value = + field_component(z_lower + (((double) kk) + 0.5) * cellDims[2]);// Ey depends on z + double z_comp_value = + field_component(x_lower + (((double) ii) + 0.5) * cellDims[0]);// Ez depends on x + // assign to "freq domain" ElectricField + E.real.x[kk][jj][ii] = x_comp_value; + E.imag.x[kk][jj][ii] = 0.; + E.real.y[kk][jj][ii] = y_comp_value; + E.imag.y[kk][jj][ii] = 0.; + E.real.z[kk][jj][ii] = z_comp_value; + E.imag.z[kk][jj][ii] = 0.; + // assign to "time domain" ElectricSplitField - use weighting that sums to 1 to check addition is behaving as planned + E_split.xy[kk][jj][ii] = x_comp_value; + E_split.xz[kk][jj][ii] = 0.; + E_split.yx[kk][jj][ii] = y_comp_value * .5; + E_split.yz[kk][jj][ii] = y_comp_value * .5; + E_split.zx[kk][jj][ii] = z_comp_value * .25; + E_split.zy[kk][jj][ii] = z_comp_value * .75; + } + } + } + // the fields are now assigned, we next need to test the interpolation itself + // we perform Nx-1, Ny-1, Nz-1 interpolations here (for Ex, Ey, Ez respectively) so with some logic checks we can do everything in one for loop + // however, cell "1" has centre 0.5*cellDims, yet we are storing the true field value and interpolated values at index "0", so there are some index-offsets here + for (int ii = 0; ii < Nx; ii++) { + for (int jj = 0; jj < Ny; jj++) { + for (int kk = 0; kk < Nz; kk++) { + // coordinates of the cell that we are interested in + double cell_centre[3]; + cell_centre[0] = x_lower + ((double) ii + 0.5) * cellDims[0]; + cell_centre[1] = y_lower + ((double) jj + 0.5) * cellDims[1]; + cell_centre[2] = z_lower + ((double) kk + 0.5) * cellDims[2]; + // Ex interpolation only goes up to Nx-1 + if (ii != Nx - 1) { + // exact field value + double Ex_exact = field_component(cell_centre[1]); + // split field interpolation + double Ex_split_interp = + E_split.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk); + // freq domain field interpolation - take real part since using entirely real field + double Ex_interp = E.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk).real(); + // compute the errors + Ex_error[kk][jj][ii] = Ex_interp - Ex_exact; + Ex_split_error[kk][jj][ii] = Ex_split_interp - Ex_exact; + } + // Ey interpolation only goes up to Ny-1 + if (jj != Ny - 1) { + double Ey_exact = field_component(cell_centre[2]); + double Ey_split_interp = + E_split.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk); + double Ey_interp = E.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk).real(); + Ey_error[kk][jj][ii] = Ey_interp - Ey_exact; + Ey_split_error[kk][jj][ii] = Ey_split_interp - Ey_exact; + } + // Ez interpolation only goes up to Nz-1 + if (kk != Nz - 1) { + double Ez_exact = field_component(cell_centre[0]); + double Ez_split_interp = + E_split.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk); + double Ez_interp = E.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk).real(); + Ez_error[kk][jj][ii] = Ez_interp - Ez_exact; + Ez_split_error[kk][jj][ii] = Ez_split_interp - Ez_exact; + } + } + } + } + + // compute error-matrix Frobenius norms + double Ex_fro_err = Ex_error.frobenius(), Ey_fro_err = Ey_error.frobenius(), + Ez_fro_err = Ez_error.frobenius(), Ex_split_fro_err = Ex_split_error.frobenius(), + Ey_split_fro_err = Ey_split_error.frobenius(), + Ez_split_fro_err = Ez_split_error.frobenius(); + + // compute max-slice errors + double Ex_ms_err = 0., Ey_ms_err = 0., Ez_ms_err = 0., Ex_split_ms_err = 0., Ey_split_ms_err = 0., + Ez_split_ms_err = 0.; + // Ex-slices + for (int jj = 0; jj < Ny; jj++) { + for (int kk = 0; kk < Nz; kk++) { + // "slices" might not constitute sequential memory + // as such, make a new array for safety + double jk_errors[Nx - 1], jk_split_errors[Nx - 1]; + for (int ii = 0; ii < Nx - 1; ii++) { + jk_errors[ii] = Ex_error[kk][jj][ii]; + jk_split_errors[ii] = Ex_split_error[kk][jj][ii]; + } + // compute norm-error of this slice + double jk_slice_error = euclidean(jk_errors, Nx - 1), + jk_split_slice_error = euclidean(jk_split_errors, Nx - 1); + // if this exceeds the current recorded maximum error, record this + if (jk_slice_error > Ex_ms_err) { Ex_ms_err = jk_slice_error; } + if (jk_split_slice_error > Ex_split_ms_err) { Ex_split_ms_err = jk_split_slice_error; } + } + } + // Ey-slices + for (int ii = 0; ii < Nx; ii++) { + for (int kk = 0; kk < Nz; kk++) { + double ik_errors[Ny - 1], ik_split_errors[Ny - 1]; + for (int jj = 0; jj < Ny - 1; jj++) { + ik_errors[jj] = Ey_error[kk][jj][ii]; + ik_split_errors[jj] = Ey_split_error[kk][jj][ii]; + } + double ik_slice_error = euclidean(ik_errors, Ny - 1), + ik_split_slice_error = euclidean(ik_split_errors, Ny - 1); + if (ik_slice_error > Ey_ms_err) { Ey_ms_err = ik_slice_error; } + if (ik_split_slice_error > Ey_split_ms_err) { Ey_split_ms_err = ik_split_slice_error; } + } + } + // Ez-slices + for (int ii = 0; ii < Nx; ii++) { + for (int jj = 0; jj < Ny; jj++) { + double ij_errors[Nz - 1], ij_split_errors[Nz - 1]; + for (int kk = 0; kk < Nz - 1; kk++) { + ij_errors[kk] = Ez_error[kk][jj][ii]; + ij_split_errors[kk] = Ez_split_error[kk][jj][ii]; + } + double ij_slice_error = euclidean(ij_errors, Nz - 1), + ij_split_slice_error = euclidean(ij_split_errors, Nz - 1); + if (ij_slice_error > Ez_ms_err) { Ez_ms_err = ij_slice_error; } + if (ij_split_slice_error > Ez_split_ms_err) { Ez_split_ms_err = ij_split_slice_error; } + } + } + + // check Frobenius errors are acceptable + CHECK(Ex_fro_err <= Ex_fro_tol + acc_tol); + CHECK(Ey_fro_err <= Ey_fro_tol + acc_tol); + CHECK(Ez_fro_err <= Ez_fro_tol + acc_tol); + CHECK(Ex_split_fro_err <= Ex_fro_tol + acc_tol); + CHECK(Ey_split_fro_err <= Ey_fro_tol + acc_tol); + CHECK(Ez_split_fro_err <= Ez_fro_tol + acc_tol); + + // check max-slice errors are acceptable + CHECK(Ex_ms_err <= Ex_ms_tol + acc_tol); + CHECK(Ey_ms_err <= Ey_ms_tol + acc_tol); + CHECK(Ez_ms_err <= Ez_ms_tol + acc_tol); + CHECK(Ex_split_ms_err <= Ex_ms_tol + acc_tol); + CHECK(Ey_split_ms_err <= Ey_ms_tol + acc_tol); + CHECK(Ez_split_ms_err <= Ez_ms_tol + acc_tol); + + // print information to the debugger/log + SPDLOG_INFO(" Component | Frobenius err. : ( benchmark ) | Max-slice err. : ( benchmark )"); + SPDLOG_INFO(" x | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ex_fro_err, Ex_fro_tol, + Ex_ms_err, Ex_ms_tol); + SPDLOG_INFO(" y | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ey_fro_err, Ey_fro_tol, + Ey_ms_err, Ey_ms_tol); + SPDLOG_INFO(" z | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ez_fro_err, Ez_fro_tol, + Ez_ms_err, Ez_ms_tol); + SPDLOG_INFO(" [Split] Component | Frobenius err. : ( benchmark ) | Max-slice err. : ( " + "benchmark )"); + SPDLOG_INFO(" x | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ex_split_fro_err, + Ex_fro_tol, Ex_split_ms_err, Ex_ms_tol); + SPDLOG_INFO(" y | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ey_split_fro_err, + Ey_fro_tol, Ey_split_ms_err, Ey_ms_tol); + SPDLOG_INFO(" z | {0:.8e} : ({1:.8e}) | {2:.8e} : ({3:.8e})", Ez_split_fro_err, + Ez_fro_tol, Ez_split_ms_err, Ez_ms_tol); +} + +/** + * @brief Test the interpolation of the H-field components to the centre of the Yee cells + * + * Each component of the H-field will take the form + * H_t(tt) = sin(2\pi tt) * exp(-tt^2) + * The decision to make this function wholey imaginary is just to test whether the imaginary part of the field is correctly picked up and worked with. + * + * We only test Fro-norm error metrics, since interpolation must occur along two axes for each component. + */ +TEST_CASE("H-field interpolation check") { + SPDLOG_INFO("===== Testing H-field BLi ====="); + // error tolerance, based on MATLAB performance + // script: benchmark_test_field_interpolation_H.m + double Hx_fro_tol = 6.1860732269207769e-02; + double Hy_fro_tol = 1.2622903101481411e-01; + double Hz_fro_tol = 7.3136417157073502e-02; + + // additional tolerance to allow for floating-point rounding imprecisions, etc + double acc_tol = 1e-12; + + // fake domain setup + double x_lower = -2., y_lower = -2., z_lower = -2.; + double extent_x = 4., extent_y = 4., extent_z = 4.; + double cellDims[3] = {0.1, 0.05, 0.025}; + // The number of cells in each direction is then 40, 80, 160 respectively + // Note that due to the possibility that Nx/cellDims[0] computing something that is not quite an integer, we need to use round() to get an int safely + int Nx = round(extent_x / cellDims[0]), Ny = round(extent_y / cellDims[1]), + Nz = round(extent_z / cellDims[2]); + SPDLOG_INFO("(Nx, Ny, Nz) = ({},{},{})", Nx, Ny, Nz); + + // setup the "split" H-field components + MagneticSplitField H_split(Nx, Ny, Nz); + H_split.allocate(); + // setup the non-split field components + MagneticField H(Nx, Ny, Nz); + H.allocate(); + // storage for pointwise errors (-1 since we don't interpolate to cell 0's centre) + Tensor3D Hx_error, Hx_split_error, Hy_error, Hy_split_error, Hz_error, Hz_split_error; + Hx_error.allocate(Nz - 1, Ny - 1, Nx); + Hx_split_error.allocate(Nz - 1, Ny - 1, Nx); + Hy_error.allocate(Nz - 1, Ny, Nx - 1); + Hy_split_error.allocate(Nz - 1, Ny, Nx - 1); + Hz_error.allocate(Nz, Ny - 1, Nx - 1); + Hz_split_error.allocate(Nz, Ny - 1, Nx - 1); + + // compute the exact field and the "split field" components + // indices range from (0,0,0) to (Nx, Ny, Nz) INCLUSIVE! + for (int ii = 0; ii <= Nx; ii++) { + for (int jj = 0; jj <= Ny; jj++) { + for (int kk = 0; kk <= Nz; kk++) { + /* The coordinates of the grids are different, which throws a spanner in the works. + Recall that the "x-grid" is at the position of the (i,j,k)-th Hx field sample. + x-grid coordinates: (i,j,k) = (x_lower + (i-0.5)*cellDims[0], y_lower + j*cellDims[1], z_lower + k*cellDims[2]) + y-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + (j-0.5)*cellDims[1], z_lower + k*cellDims[2]) + z-grid coordinates: (i,j,k) = (x_lower + i*cellDims[0], y_lower + j*cellDims[1], z_lower + (k-0.5)*cellDims[2]) + The components only depend on one of the coordinate directions, so we don't need to compute all 9 values. + */ + double x_comp_value = + field_component(y_lower + ((double) jj) * cellDims[1]);// Hx depends on y + double y_comp_value = + field_component(z_lower + ((double) kk) * cellDims[2]);// Hy depends on z + double z_comp_value = + field_component(x_lower + ((double) ii) * cellDims[0]);// Hz depends on x + // assign to "freq domain" ElectricField + H.imag.x[kk][jj][ii] = x_comp_value; + H.real.x[kk][jj][ii] = 0.; + H.imag.y[kk][jj][ii] = y_comp_value; + H.real.y[kk][jj][ii] = 0.; + H.imag.z[kk][jj][ii] = z_comp_value; + H.real.z[kk][jj][ii] = 0.; + // assign to "time domain" ElectricSplitField - use weighting that sums to 1 to check addition is behaving as planned + H_split.xy[kk][jj][ii] = x_comp_value; + H_split.xz[kk][jj][ii] = 0.; + H_split.yx[kk][jj][ii] = y_comp_value * .25; + H_split.yz[kk][jj][ii] = y_comp_value * .75; + H_split.zx[kk][jj][ii] = z_comp_value * .125; + H_split.zy[kk][jj][ii] = z_comp_value * .875; + } + } + } + // the fields are now assigned, we next need to test the interpolation itself + for (int ii = 0; ii < Nx; ii++) { + for (int jj = 0; jj < Ny; jj++) { + for (int kk = 0; kk < Nz; kk++) { + // coordinates of the cell that we are interested in + double cell_centre[3]; + cell_centre[0] = x_lower + ((double) ii + 0.5) * cellDims[0]; + cell_centre[1] = y_lower + ((double) jj + 0.5) * cellDims[1]; + cell_centre[2] = z_lower + ((double) kk + 0.5) * cellDims[2]; + // Hx interpolation only goes up to Ny-1, Nz-1 + if ((jj != Ny - 1) && (kk != Nz - 1)) { + // exact field value + double Hx_exact = field_component(cell_centre[1]); + // split field interpolation + double Hx_split_interp = + H_split.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk); + // freq domain field interpolation - take real part since using entirely real field + double Hx_interp = + H.interpolate_to_centre_of(AxialDirection::X, ii, jj, kk).imag(); + // compute the errors + Hx_error[kk][jj][ii] = Hx_interp - Hx_exact; + Hx_split_error[kk][jj][ii] = Hx_split_interp - Hx_exact; + } + // Hy interpolation only goes up to Nx-1,Nz-1 + if ((ii != Nx - 1) && (kk != Nz - 1)) { + double Hy_exact = field_component(cell_centre[2]); + double Hy_split_interp = + H_split.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk); + double Hy_interp = + H.interpolate_to_centre_of(AxialDirection::Y, ii, jj, kk).imag(); + Hy_error[kk][jj][ii] = Hy_interp - Hy_exact; + Hy_split_error[kk][jj][ii] = Hy_split_interp - Hy_exact; + } + // Hz interpolation only goes up to Nx-1,Ny-1 + if ((ii != Nx - 1) && (jj != Ny - 1)) { + double Hz_exact = field_component(cell_centre[0]); + double Hz_split_interp = + H_split.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk); + double Hz_interp = H.interpolate_to_centre_of(AxialDirection::Z, ii, jj, kk).imag(); + Hz_error[kk][jj][ii] = Hz_interp - Hz_exact; + Hz_split_error[kk][jj][ii] = Hz_split_interp - Hz_exact; + } + } + } + } + // compute Frobenius norms + double Hx_fro_err = Hx_error.frobenius(), + Hy_fro_err = Hy_error.frobenius(), + Hz_fro_err = Hz_error.frobenius(), + Hx_split_fro_err = Hx_split_error.frobenius(), + Hy_split_fro_err = Hy_split_error.frobenius(), + Hz_split_fro_err = Hz_split_error.frobenius(); + + // check Frobenius errors are acceptable + CHECK(Hx_fro_err <= Hx_fro_tol + acc_tol); + CHECK(Hy_fro_err <= Hy_fro_tol + acc_tol); + CHECK(Hz_fro_err <= Hz_fro_tol + acc_tol); + CHECK(Hx_split_fro_err <= Hx_fro_tol + acc_tol); + CHECK(Hy_split_fro_err <= Hy_fro_tol + acc_tol); + CHECK(Hz_split_fro_err <= Hz_fro_tol + acc_tol); + + // print information to the debugger/log + SPDLOG_INFO(" Component | Frobenius err. : ( benchmark )"); + SPDLOG_INFO(" x | {0:.8e} : ({1:.8e})", Hx_fro_err, Hx_fro_tol); + SPDLOG_INFO(" y | {0:.8e} : ({1:.8e})", Hy_fro_err, Hy_fro_tol); + SPDLOG_INFO(" z | {0:.8e} : ({1:.8e})", Hz_fro_err, Hz_fro_tol); + SPDLOG_INFO(" [Split] Component | Frobenius err. : ( benchmark )"); + SPDLOG_INFO(" x | {0:.8e} : ({1:.8e})", Hx_split_fro_err, Hx_fro_tol); + SPDLOG_INFO(" y | {0:.8e} : ({1:.8e})", Hy_split_fro_err, Hy_fro_tol); + SPDLOG_INFO(" z | {0:.8e} : ({1:.8e})", Hz_split_fro_err, Hz_fro_tol); +} diff --git a/tdms/tests/unit/test_fields.cpp b/tdms/tests/unit/test_fields.cpp index da461393f..61dbe8fdf 100644 --- a/tdms/tests/unit/test_fields.cpp +++ b/tdms/tests/unit/test_fields.cpp @@ -4,10 +4,12 @@ */ #include #include + #include "field.h" #include "globals.h" using namespace std; +using namespace tdms_math_constants; template inline bool is_close(T a, T b){ @@ -27,7 +29,6 @@ TEST_CASE("Test electric field angular norm addition") { double DT = 0.1; int N = 3; int N_T = 5; - auto I = tdms_math_constants::IMAGINARY_UNIT; auto E = ElectricField(); E.ft = 1.0; @@ -36,7 +37,7 @@ TEST_CASE("Test electric field angular norm addition") { params.dt = DT; // z = e^(iω(n+1)dt) / N_t - auto expected = exp(OMEGA * ((double)N + 1) * DT * I) / ((double) N_T); + auto expected = exp(OMEGA * ((double)N + 1) * DT * IMAGINARY_UNIT) / ((double) N_T); E.add_to_angular_norm(N, N_T, params); REQUIRE(is_close(E.angular_norm, expected)); @@ -48,17 +49,15 @@ TEST_CASE("Test magnetic field angular norm addition") { double DT = 0.2; int N = 7; int N_T = 5; - auto I = tdms_math_constants::IMAGINARY_UNIT; auto H = MagneticField(); H.ft = 1.0; - auto params = SimulationParameters(); params.omega_an = OMEGA; params.dt = DT; // z = e^(iω(n+1/2)dt) / N_t - auto expected = exp(OMEGA * ((double)N + 0.5) * DT * I) / ((double) N_T); + auto expected = exp(OMEGA * ((double)N + 0.5) * DT * IMAGINARY_UNIT) / ((double) N_T); H.add_to_angular_norm(N, N_T, params); REQUIRE(is_close(H.angular_norm, expected)); diff --git a/tdms/tests/unit/test_interpolation_determination.cpp b/tdms/tests/unit/test_interpolation_determination.cpp index 0c3ad828a..d2794e916 100644 --- a/tdms/tests/unit/test_interpolation_determination.cpp +++ b/tdms/tests/unit/test_interpolation_determination.cpp @@ -1,10 +1,12 @@ /** * @file test_interpolation_determination.cpp - * @brief Tests of the interpolation functions. + * @author William Graham (ccaegra@ucl.ac.uk) + * @brief Tests the logic that determines which interpolation schemes are appropriate. */ #include #include +#include using namespace std; @@ -15,9 +17,7 @@ using namespace std; * * THIS FUNCTION WILL BE DEPRECIATED UPON SWITCHING TO THE BLI FRAMEWORK */ -TEST_CASE("checkInterpolationPoints: exceptions thrown") -{ - +TEST_CASE("checkInterpolationPoints: exceptions thrown") { // setup some fake field dimensions int I = 6, J = 7, K = 8; @@ -49,9 +49,7 @@ TEST_CASE("checkInterpolationPoints: exceptions thrown") * * THIS FUNCTION WILL BE DEPRECIATED UPON SWITCHING TO THE BLI FRAMEWORK */ -TEST_CASE("checkInterpolationPoints: check valid inputs") -{ - +TEST_CASE("checkInterpolationPoints: check valid inputs") { // setup some fake field dimensions int I = 6, J = 7, K = 8; @@ -76,14 +74,15 @@ TEST_CASE("checkInterpolationPoints: check valid inputs") * @brief Test whether best_scheme correctly determines the appropriate interpolation scheme to use, given the number of Yee cells either side of cell (i,j,k) * */ -TEST_CASE("best_scheme: correct interpolation chosen") -{ - +TEST_CASE("best_interp_scheme: correct interpolation chosen") { + SPDLOG_INFO("===== Testing interpolation scheme selection logic ====="); int N = 10; - // should throw out_of_range exception if interpolation is impossible (<4 Yee cells in direction) - REQUIRE_THROWS_AS(best_scheme(3, 2), out_of_range); - // should throw out_of_range exception if Yee cell of invalid index is requested - REQUIRE_THROWS_AS(best_scheme(N, 0), out_of_range); + // should throw out_of_range exception if interpolation is impossible (<3 Yee cells in direction) + REQUIRE_THROWS_AS(best_scheme(1, 0), out_of_range); + REQUIRE_THROWS_AS(best_scheme(2, 0), out_of_range); + REQUIRE_THROWS_AS(best_scheme(2, 1), out_of_range); + // should throw out_of_range exception if Yee cell if invalid index is requested + REQUIRE_THROWS_AS(best_scheme(N, -1), out_of_range); REQUIRE_THROWS_AS(best_scheme(N, N), out_of_range); /* Suppose we have N >= 8 Yee cells in a dimension. The program should determine: @@ -93,10 +92,10 @@ TEST_CASE("best_scheme: correct interpolation chosen") - cell_id == N-3,N-2,N-1 : Use BAND_LIMITED_(4,5,6) scheme respectively */ N = 10; - CHECK(best_scheme(N, 1).get_priority() == BAND_LIMITED_0); - CHECK(best_scheme(N, 2).get_priority() == BAND_LIMITED_1); - CHECK(best_scheme(N, 3).get_priority() == BAND_LIMITED_2); - for (int i = 4; i <= N - 4; i++) { + CHECK(best_scheme(N, 0).get_priority() == BAND_LIMITED_0); + CHECK(best_scheme(N, 1).get_priority() == BAND_LIMITED_1); + CHECK(best_scheme(N, 2).get_priority() == BAND_LIMITED_2); + for (int i = 3; i <= N - 4; i++) { CHECK(best_scheme(N, i).get_priority() == BAND_LIMITED_3); } CHECK(best_scheme(N, N - 3).get_priority() == BAND_LIMITED_4); @@ -109,9 +108,8 @@ TEST_CASE("best_scheme: correct interpolation chosen") - cell_id == 2,...,N-2 : Use CUBIC_MIDDLE - cell_id == N-1 : Use CUBIC_LAST */ - N = 7; - CHECK(best_scheme(N, 1).get_priority() == CUBIC_INTERP_FIRST); - for(int i=2; i<=N-2; i++) {CHECK(best_scheme(N, i).get_priority() == CUBIC_INTERP_MIDDLE);} + N = 6; + CHECK(best_scheme(N, 0).get_priority() == CUBIC_INTERP_FIRST); + for(int i=1; i<=N-2; i++) {CHECK(best_scheme(N, i).get_priority() == CUBIC_INTERP_MIDDLE);} CHECK(best_scheme(N, N-1).get_priority() == CUBIC_INTERP_LAST); } - diff --git a/tdms/tests/unit/test_interpolation_functions.cpp b/tdms/tests/unit/test_interpolation_functions.cpp index 71073b229..e6d4c5a63 100644 --- a/tdms/tests/unit/test_interpolation_functions.cpp +++ b/tdms/tests/unit/test_interpolation_functions.cpp @@ -4,12 +4,17 @@ * @brief Tests the performance of the interpolation functions, using 1D data mimicing a coordinate axes */ #include +#include #include #include -#include +#include +#include +#include "globals.h" #include "interpolation_methods.h" +using namespace tdms_math_constants; +using Catch::Approx; using namespace std; /** @@ -17,10 +22,10 @@ using namespace std; * polynomial fields up to cubic order are interpolated exactly (to within * machine error) * + * Checks are run on both the old interp{1,2,3} functions and newer const Interp_scheme instances. Old cubic methods will be redundant upon integration of BLi into the codebase. */ -TEST_CASE("test_interpolation_functions: testing that cubic interpolation is exact") -{ - +TEST_CASE("test_interpolation_functions: testing that cubic interpolation is exact") { + SPDLOG_INFO("===== Testing exact cubic interpolation ====="); // equidistant points double x[] = {0., 1., 2., 3.}; // test acceptence tolerance. Allow for FLOP imprecision and rounding errors @@ -28,55 +33,90 @@ TEST_CASE("test_interpolation_functions: testing that cubic interpolation is exa double tol = 4e1 * __DBL_EPSILON__; // constant field - double c0 = 3.1415; + double c0 = M_PI; + // array that will be passed to Interp_scheme::interpolate + double interp_data[4] = {c0, c0, c0, c0}; + double v1 = c0, v2 = c0, v3 = c0, v4 = c0; double v12 = c0, v23 = c0, v34 = c0; - CHECK(abs(v12 - interp2(v1, v2, v3, v4)) <= tol); - CHECK(abs(v23 - interp1(v1, v2, v3, v4)) <= tol); - CHECK(abs(v34 - interp3(v1, v2, v3, v4)) <= tol); + // check old interp methods + CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); + // check Interp_scheme class method + CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); + CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); + CHECK(v34 == Approx(CBLst.interpolate(interp_data)).epsilon(tol)); // linear double c1 = -2.7182818; v1 += c1 * x[0]; + interp_data[0] += c1 * x[0]; v2 += c1 * x[1]; + interp_data[1] += c1 * x[1]; v3 += c1 * x[2]; + interp_data[2] += c1 * x[2]; v4 += c1 * x[3]; + interp_data[3] += c1 * x[3]; v12 += c1 * (x[1] + x[0]) / 2.; v23 += c1 * (x[2] + x[1]) / 2.; v34 += c1 * (x[3] + x[2]) / 2.; - CHECK(abs(v12 - interp2(v1, v2, v3, v4)) <= tol); - CHECK(abs(v23 - interp1(v1, v2, v3, v4)) <= tol); - CHECK(abs(v34 - interp3(v1, v2, v3, v4)) <= tol); + // check old interp methods + CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); + // check Interp_scheme class method + CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); + CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); + CHECK(v34 == Approx(CBLst.interpolate(interp_data)).epsilon(tol)); // quadratic double c2 = 9.81; v1 += c2 * x[0] * x[0]; + interp_data[0] += c2 * x[0] * x[0]; v2 += c2 * x[1] * x[1]; + interp_data[1] += c2 * x[1] * x[1]; v3 += c2 * x[2] * x[2]; + interp_data[2] += c2 * x[2] * x[2]; v4 += c2 * x[3] * x[3]; + interp_data[3] += c2 * x[3] * x[3]; v12 += c2 * (x[1] + x[0]) * (x[1] + x[0]) / 4.; v23 += c2 * (x[2] + x[1]) * (x[2] + x[1]) / 4.; v34 += c2 * (x[3] + x[2]) * (x[3] + x[2]) / 4.; - CHECK(abs(v12 - interp2(v1, v2, v3, v4)) <= tol); - CHECK(abs(v23 - interp1(v1, v2, v3, v4)) <= tol); - CHECK(abs(v34 - interp3(v1, v2, v3, v4)) <= tol); + // check old interp methods + CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); + // check Interp_scheme class method + CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); + CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); + CHECK(v34 == Approx(CBLst.interpolate(interp_data)).epsilon(tol)); // cubic double c3 = 4.2; v1 += c3 * x[0] * x[0] * x[0]; + interp_data[0] += c3 * x[0] * x[0] * x[0]; v2 += c3 * x[1] * x[1] * x[1]; + interp_data[1] += c3 * x[1] * x[1] * x[1]; v3 += c3 * x[2] * x[2] * x[2]; + interp_data[2] += c3 * x[2] * x[2] * x[2]; v4 += c3 * x[3] * x[3] * x[3]; + interp_data[3] += c3 * x[3] * x[3] * x[3]; v12 += c3 * (x[1] + x[0]) * (x[1] + x[0]) * (x[1] + x[0]) / 8.; v23 += c3 * (x[2] + x[1]) * (x[2] + x[1]) * (x[2] + x[1]) / 8.; v34 += c3 * (x[3] + x[2]) * (x[3] + x[2]) * (x[3] + x[2]) / 8.; - CHECK(abs(v12 - interp2(v1, v2, v3, v4)) <= tol); - CHECK(abs(v23 - interp1(v1, v2, v3, v4)) <= tol); - CHECK(abs(v34 - interp3(v1, v2, v3, v4)) <= tol); + // check old interp methods + CHECK(v12 == Approx(interp2(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v23 == Approx(interp1(v1, v2, v3, v4)).epsilon(tol)); + CHECK(v34 == Approx(interp3(v1, v2, v3, v4)).epsilon(tol)); + // check Interp_scheme class method + CHECK(v12 == Approx(CBFst.interpolate(interp_data)).epsilon(tol)); + CHECK(v23 == Approx(CBMid.interpolate(interp_data)).epsilon(tol)); + CHECK(v34 == Approx(CBLst.interpolate(interp_data)).epsilon(tol)); } /** @@ -85,13 +125,12 @@ TEST_CASE("test_interpolation_functions: testing that cubic interpolation is exa * * Note - the coefficients are not required to sum to unity! */ -TEST_CASE("bandlimited_interpolation: check that the interpolation constant values all sum to the same value") -{ - +TEST_CASE("bandlimited_interpolation: check that the interpolation constant values all sum to the same value") { + SPDLOG_INFO("===== Testing BLi coefficient sum ====="); /* Tolerance to accept imprecision to max. 16 FLOPs implies max discrepency of 8*_DBL_EPSILON__ error (8 additions, error in each number max __DBL_EPSILON__). - tol should be \approx 10 * __DBL_EPSILON__ , + tol should be \approx 10 * __DBL_EPSILON__ , after accounting for funny-business multiplying everything by 1 */ double tol = 1e1 * __DBL_EPSILON__; @@ -109,192 +148,233 @@ TEST_CASE("bandlimited_interpolation: check that the interpolation constant valu coeff_sums[7] = BL7.interpolate(a); // now check that the entries of coeff_sums are the same int n = 8; - while (--n > 0 && abs(coeff_sums[n] - coeff_sums[0]) < tol); + while (--n > 0 && abs(coeff_sums[n] - coeff_sums[0]) < tol) + ; // we only reach the end of the while loop, IE get to n==0, when all elements in the array are the same REQUIRE(n == 0); } /** - * @brief We will check that BLi interpolation gives comparible error to the equivalent functions in MATLAB + * @brief We will check that BLi interpolation over real-valued data gives comparible error to the equivalent functions in MATLAB * * For 100 sample points, we will use BLi to interpolate the following functions with 100 sample points: * - The constant function 1 : range 0,1 : max. element-wise error (MATLAB) 2.82944733e-04 * - sin(2\pi x) : range 0,1 : max. element-wise error (MATLAB) 2.63468327e-04 - * - e^(-1/(2|x-0.5|-1))^2) : range 0,1 : max. element-wise error (MATLAB) 1.57403846e-03 + * - pulse function : range 0,1 : max. element-wise error (MATLAB) 4.87599933e-04 + * Error values produced from benchmark_test_interpolation_functions.m * * We will then compare the maximum of the pointwise error between the interpolated values and exact values to the same quantity computed via MATLAB's interp function, and determine whether the order of magnitude of the errors is the same. * * For readability, we break this out into separate test cases for each function. - * + * */ -TEST_CASE("bandlimited_interpolation: order of error, constant function") -{ - - int nSamples = 100; - double field_positions[nSamples], Yee_cell_centres[nSamples]; - - // constant function variables - double const_fn_vals[nSamples], const_fn_interp[nSamples], const_fn_exact[nSamples], const_fn_errors[nSamples-1], const_fn_max_error, const_fn_MATLAB_error; - - // input the MATLAB error - const_fn_MATLAB_error = 2.82944733e-04; - - // setup the sample points, Yee cell centres, and function values (for sampling and exactness) - for (int i = 0; i < nSamples; i++) - { - field_positions[i] = ((double)i / (double)nSamples); - Yee_cell_centres[i] = field_positions[i] - 1./(2.*(double)nSamples); - - const_fn_vals[i] = 1.; const_fn_exact[i] = 1.; - } - // Yee cell 0 has no value "to the left" - this will change with BL_TO_CELL_0 being included. - // also recall that best_scheme(nSamples, i) returns the scheme that interpolates to the centre of cell i, IE, to position Yee_cell_centres[i]. +/** + * @brief Test BLi performance on the constant function. + * - The constant function 1 : range 0,1 : max. element-wise error (MATLAB) 2.82944733e-04 + */ +TEST_CASE("(real-valued) Band limited interpolation: constant function") { + SPDLOG_INFO("===== (real valued) BLi: constant function ====="); + int nSamples = 100; // number of datapoints in this dimension + double const_fn_data[nSamples]; // function data (only need 8 data points, but simulate 100 cells) + fill_n(const_fn_data, nSamples, 1); // initalise to f(x)=1 + double const_fn_interp[nSamples-1]; // interpolated values + double const_fn_errors[nSamples-1]; // error in interpolated values + double const_fn_max_error; // maximum error + + // hardcoding the MATLAB error (benchmark_test_interpolation_functions.m) + double const_fn_MATLAB_error = 2.82944733e-04; // constant function interpolation - const_fn_interp[1] = best_scheme(nSamples, 1).interpolate(const_fn_vals); - const_fn_interp[2] = best_scheme(nSamples, 2).interpolate(const_fn_vals); - const_fn_interp[3] = best_scheme(nSamples, 3).interpolate(const_fn_vals); - for (int i=4; i= 1, + * e^{-1/(1-|x|^2)} & when |x| < 1 + * \end{cases} + * This function is compactly supported over the interval [-1,1]. + * By evaluating pulse(x) = \phi(3[2x-1]), we obtain a smooth function with support in [1/3,2/3] \subset [0,1] + * + * @param x Point of evaluation + * @return double Evaluted value + */ +inline double pulse(double x) { + double absxhat = abs(3. * (2.*x - 1.)); + if (absxhat >= 1) { + return 0.; + } + else { + return exp( -1. / (1 - absxhat*absxhat) ); + } +} - // input the MATLAB error - pulse_MATLAB_error = 1.57403846e-03; +/** + * @brief Test BLi performance on the compact pulse. + * - pulse function : range 0,1 : max. element-wise error (MATLAB) 4.87599933e-04 + */ +TEST_CASE("(real-valued) Band limited interpolation: compact pulse") { + SPDLOG_INFO("===== (real valued) BLi: compact pulse ====="); + int nSamples = 100; // number of datapooints in this dimension + double spacing = 1. / (double) (nSamples - 1);// spacing between datapoints + double xi[nSamples]; // coordinates of the samples + double xi5[nSamples]; // coordinates of the interpolation points + double f_data[nSamples]; // function data at xi + double f_exact[nSamples - 1]; // function exact values at xi5 + double f_interp[nSamples - 1]; // interpolated values at xi5 + double f_errors[nSamples - 1]; // error at xi5 + double max_error; // maximum error across xi5 points + + // hardcoding the MATLAB error (benchmark_test_interpolation_functions.m) + double pulse_MATLAB_error = 4.87599933e-04; // setup the sample points, Yee cell centres, and function values (for sampling and exactness) - for (int i = 0; i < nSamples; i++) - { - field_positions[i] = ((double)i / (double)nSamples); - Yee_cell_centres[i] = field_positions[i] - 1. / (2. * (double)nSamples); - - if (2. * abs(field_positions[i] - 0.5) >= 1) - { - pulse_vals[i] = 0.; - } - else - { - double x = 2. * abs(field_positions[i] - 0.5); - pulse_vals[i] = exp(-1. / ((x - 1.) * (x - 1.))); - } - if (2. * abs(Yee_cell_centres[i] - 0.5) >= 1) - { - pulse_exact[i] = 0.; - } - else - { - double x = 2. * abs(Yee_cell_centres[i] - 0.5); - pulse_exact[i] = exp(-1. / ((x - 1.) * (x - 1.))); - } + for (int i = 0; i < nSamples; i++) { + xi[i] = ((double) i) * spacing; + f_data[i] = pulse(xi[i]); + if (i != nSamples - 1) { + xi5[i] = xi[i] + spacing / 2.; + f_exact[i] = pulse(xi5[i]); + } } - // Yee cell 0 has no value "to the left" - this will change with BL_TO_CELL_0 being included. - // Also recall that best_scheme(nSamples, i) returns the scheme that interpolates to the centre of cell i, IE, to position Yee_cell_centres[i]. + // sin function interpolation + for (int i = 0; i < nSamples - 1; i++) { + InterpolationScheme scheme = best_scheme(nSamples - 1, i); + f_interp[i] = scheme.interpolate(f_data, i + 1 - scheme.number_of_datapoints_to_left); + // Compare interpolated values to the true values + f_errors[i] = abs(f_exact[i] - f_interp[i]); + } - // pulse function interpolation - pulse_interp[1] = best_scheme(nSamples, 1).interpolate(pulse_vals); - pulse_interp[2] = best_scheme(nSamples, 2).interpolate(pulse_vals); - pulse_interp[3] = best_scheme(nSamples, 3).interpolate(pulse_vals); - for (int i = 4; i < nSamples - 4; i++) - { - // need to offset now so that the correct sample points are provided - pulse_interp[i] = best_scheme(nSamples, i).interpolate(pulse_vals, i - 4); + // get maximum error + max_error = *max_element(f_errors, f_errors + nSamples - 1); + // compare O.o.Mag of error - fail if we are orders of magnitude out + REQUIRE(floor(log10(max_error)) <= floor(log10(pulse_MATLAB_error))); + // compare absolute error - flag (and fail, but less harshly) if we are doing worse than we expect (but are close) + CHECK(max_error < pulse_MATLAB_error); + // report test results + SPDLOG_INFO("Error: {0:.8e} | Benchmark: {1:.8e}", max_error, pulse_MATLAB_error); +} + +/** + * @brief We will check that BLi interpolation over complex-valued data gives comparible error to the equivalent functions in MATLAB + * + * For 100 sample points, we will use BLi to interpolate the following complex-valued function with 100 sample points: + * real part of sin(2\pi x) + * imag part of pulse function. + * + * Interoplation will then be tested against over the range [0,1], the max element-wise error (by absolute value) will be determined. We will then check that this is of the same order of magnitude as the error produced by MATLAB, 5.35317432e-04. + */ +TEST_CASE("(complex-valued) Band limited interpolation") { + SPDLOG_INFO("===== (complex valued) BLi: complex function test case ====="); + int nSamples = 100; // number of "Yee cells" in this dimension + double spacing = 1. / (double)(nSamples - 1); // spacing between Yee cell centres + double xi[nSamples]; // positions of the "field components" + double xi5[nSamples]; // positions of the "Yee cell" centres, xi5[i] = centre of cell i + complex f_data[nSamples]; // function data at xi + complex f_exact[nSamples - 1]; // function exact values at xi5 + complex f_interp[nSamples - 1]; // interpolated values at xi5 + double f_abs_errors[nSamples - 1]; // error at xi5 + + double max_error; // maximum (abs) error across xi5 points + // hardcoding the MATLAB error (benchmark_test_interpolation_functions.m) + double MATLAB_error = 5.35317432e-04; + + // setup the sample points, Yee cell centres, and function values (for sampling and exactness) + for (int i = 0; i < nSamples; i++) { + xi[i] = ((double) i) * spacing; + f_data[i] = s2pi(xi[i]) + IMAGINARY_UNIT*pulse(xi[i]); + if (i != nSamples - 1) { + xi5[i] = xi[i] + spacing / 2.; + f_exact[i] = s2pi(xi5[i]) + IMAGINARY_UNIT*pulse(xi5[i]); + } } - pulse_interp[nSamples - 4] = best_scheme(nSamples, nSamples - 4).interpolate(pulse_vals, nSamples - 8); - pulse_interp[nSamples - 3] = best_scheme(nSamples, nSamples - 3).interpolate(pulse_vals, nSamples - 8); - pulse_interp[nSamples - 2] = best_scheme(nSamples, nSamples - 2).interpolate(pulse_vals, nSamples - 8); - pulse_interp[nSamples - 1] = best_scheme(nSamples, nSamples - 1).interpolate(pulse_vals, nSamples - 8); - // compare interpolated values to the true values. NOTE: index 0 is invalid as we currently don't interpolate to here - for (int i = 0; i < nSamples - 1; i++) - { - pulse_errors[i] = abs(pulse_exact[i + 1] - pulse_interp[i + 1]); + // sin function interpolation + for (int i = 0; i < nSamples - 1; i++) { + InterpolationScheme scheme = best_scheme(nSamples - 1, i); + f_interp[i] = scheme.interpolate(f_data, i + 1 - scheme.number_of_datapoints_to_left); + // Compare interpolated values to the true values + f_abs_errors[i] = abs(f_exact[i] - f_interp[i]); } // get maximum error - pulse_max_error = *max_element(pulse_errors, pulse_errors + nSamples - 1); - - // compare O.o.Mag of error - CHECK(floor(log10(pulse_max_error)) <= floor(log10(pulse_MATLAB_error))); - cout << "pulse \t | " + to_string(pulse_max_error) + "\t\t | " + to_string(pulse_MATLAB_error) + "\n"; + max_error = *max_element(f_abs_errors, f_abs_errors + nSamples - 1); + // compare O.o.Mag of error - fail if we are orders of magnitude out + REQUIRE(floor(log10(max_error)) <= floor(log10(MATLAB_error))); + // compare absolute error - flag (and fail, but less harshly) if we are doing worse than we expect (but are close) + CHECK(max_error < MATLAB_error); + // report test results + SPDLOG_INFO("Error: {0:.8e} | Benchmark {1:.8e}", max_error, MATLAB_error); }