E- and H-field interpolation tests, frequency domain interpolation (#80)

* Extract checks for field extraction being within FDTD grid checker into its own function * Tidy up a couple more docstrings * Create test_interpolate.cpp placeholder * Update interpolateFieldCentralE/H_TE/M error checks to match docstrings * Update docstrings to match function performance * Add unit test to build; add unit test for exception handling near edges of computational grid * Add matrix.h header to interpolate.h so it compiles * Update comments and docstring for checkInterpolationPoints * Update interpolate header with new functions * Add test for determineInterpScheme * Tests for determineInterpScheme * Add tests for cubic interpolation functions * Account for FLOP error in cubic interpolation tests * Update test for FLOP accuracy * Added interpolation function for Hx component in 3D * Update tdms/src/interpolate.cpp Co-authored-by: Tom Young <[email protected]> * Avoid redeclarations of arrays * Update error in test criteria * Fix determineInterpScheme's determination of right_cell_buffer * Create enums for identifiers, update coding practices * Refactor exception checks for checkInterpolationPoints * Refactor H-field interp into new file * Refactor interpolation methods into separate files for readability * Refactor E-field interpolation methods * Complete refactor of E-field interpolation methods. Complex data yet to be coded up * Function reorder * Add extra aux functions for the H-field * Create H-field x-component extraction functions * Update test framework with reworked enums * Update interpolate header * Add interpolation_methods include to interpolate * Looking into header redefinition issue * Header def fix2 * Add interpolation_methods to CMakeLists sources * Test typo, <=4 -> <4 * Correct interpolation_methods.cpp logic flagged by test * Update CMakeLists to identify new functions * Update tdms/include/interpolation_methods.h Co-authored-by: Tom Young <[email protected]> * Make interpScheme class and update syntaxes * docstrings -> headers * Update test framework with new interpolation checks * Revert changes to checkInterpolationPoints due to errors * Revert changes to test_interpolate.cpp * Add overloaded call operator to interpScheme class * Update interpScheme class and make constant instances * Don't double-define default value * Update docstring, impliment BLi to cell 0 if we want to include this functionality in future * Update test_interpolate get_value -> get_priority * unit tests... * Update names to constant refs * Update bandlimited interpolation test * Switch -1 -> +1s error in E-field component schemes * Hx component interpolation written * Add interpolation methods for other H-field components * Update CMakelists from main * fix bracket mismatch * Update BLi coeff sum tests * breakout scheme determination tests into separate file * typo in filenames * Breakout interpolation test functions into separate cases * move abs() mismatch * remove double test * E-field component tests * E-field test * E-field test compares errors to MATLAB * Update tdms/tests/unit/test_interp_fns.cpp Co-authored-by: Mosè Giordano <[email protected]> * update comments prior to tests * Add Bli vs Cubic test, update misleading docstring * further benchmarking bli vs cubic * add all tests, some fail to compete with MATLAB benchmark * include field and benchmark tests - currently fail * Update H/E field test readability * UPDATE: Adjust interpolation logic and process for H-field * Add definition of interpolateTimeDomainHField * H-field MATLAB benchmark test running * Update tdms/include/interpolate_Hfield.h Co-authored-by: Tom Young <[email protected]> * Update tdms/include/interpolate_Hfield.h Co-authored-by: Tom Young <[email protected]> * Update I -> nI and simila * Pull-out general, bulky comment on interpolate_Hfield.cpp * Condense assertions into single error check for each component * add BLi vs cubic test * Add CB{Lst,Mid,Fst} to cubic test * Fix failure-test case * Field interpolation tests updated: to decide on rtol condition? * remove repeated test * remove test that was duplicated in merge with main * Update tdms/tests/unit/test_interp_fns.cpp Co-authored-by: Sam Cunliffe <[email protected]> * Update tdms/tests/unit/test_interp_fns.cpp Co-authored-by: Sam Cunliffe <[email protected]> * Update tdms/tests/unit/test_interp_determination.cpp Co-authored-by: Sam Cunliffe <[email protected]> * Update tdms/src/interpolate_Efield.cpp Co-authored-by: Sam Cunliffe <[email protected]> * Update tdms/include/interpolate_Hfield.h Co-authored-by: Sam Cunliffe <[email protected]> * Apply suggestions from code review Co-authored-by: Sam Cunliffe <[email protected]> * Quick review comment fixes * Change index -> number_of_datapoints_to_left * Actually merge CMakeLists properly * Fix MATLAB scripts for benchmarking * Make BLi n_trials a const so memory assignment works on MacOS build? * Use rigorous delete syntax * Add complex interpolation overloads * Add freq-domain interpolation functions * Add namespacing to TDMS globals * Update electric_field with new globals.h * Update iterator with new globals convention * Update global refs in other files * Revert one that was meant to be an I * macOS → latest * Add MATLAB benchmark scripts * Remove unused globals * Fix namespacing issues with globals * Change interpScheme.interpolate to templated function * Revert "Change interpScheme.interpolate to templated function" This reverts commit e6d5a0c. * Fix complex typing * Clean up file docstrings and includes * General formatting tidyup * Typesetting for codebase * Fix macOS build (#117) * macOS → latest * Remove DOMP_ROOT flag * Remove DCXX_ROOT flag * Use dev build 4127d37 Co-authored-by: willGraham01 <[email protected]> * make macos tests consistent with main merge fix * Lowercase namespaces, include tidy up * check test namespacing and iterator includes * Update test --verbose output * Change to floats in multiplications * turns out it helps if you initalise to 0 before adding things * round > ceil * One j that should've been an i * Move E-field interpolation to templated types * Change H-field interpolation functions to templates * Remove redundant source files now we have templates * Remove unneccessary include * Add split/unified field tests to catch2 * Add newlines to things * Sam's comments #1 * Sam's comments #2 * Fix namespacing issue affecting tests * Peter's comments #1 * Peter's function suggestions * interpScheme additions post-merge * Interpolation Functions are now Field (Class) Methods (#130) * Fix tensor3D allocation * Add split-E-field interpolation as class method * Update split-H-field class with interpolation method * Overhaul ElectricField to have interpolation method * H-field has interpolation methods * Remove redundant headers now class methods are established * Update tdms/include/arrays.h Co-authored-by: Sam Cunliffe <[email protected]> * Move interpolation functionality higher up class hierarchy * Add overrides for MacOS compilation Co-authored-by: Sam Cunliffe <[email protected]> * Update the matlab benchmarking scripts with Peter's function suggestions * Fix E-field test offsets * Pre index overhaul * Index issue fixed, but H-field tests need examination * Fix mag field mismatches Co-authored-by: Tom Young <[email protected]> Co-authored-by: Mosè Giordano <[email protected]> Co-authored-by: Sam Cunliffe <[email protected]>
UCL · Oct 17, 2022 · d98a564 · d98a564
1 parent 55e3d91
commit d98a564
Show file tree

Hide file tree

Showing 26 changed files with 2,024 additions and 635 deletions.
diff --git 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
@@ -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
@@ -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<double> 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<std::complex<double>>
   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<double> 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<double> The interpolated field value
+   */
+  virtual std::complex<double> 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<double> The interpolated component value
+   */
+  std::complex<double> 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<double> The interpolated component value
+   */
+  std::complex<double> 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
@@ -38,6 +38,12 @@ typedef struct complex_vector
 	std::complex<double> 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
diff --git a/tdms/include/interpolate_Hfield.h b/tdms/include/interpolate_Hfield.h
diff --git a/tdms/include/interpolation_methods.h b/tdms/include/interpolation_methods.h
@@ -1,12 +1,14 @@
 /**
- * @file interpolation_methods.cpp
+ * @file interpolation_methods.h
  * @author William Graham ([email protected])
  * @brief InterpScheme class methods and supporting functions
  * 
  * Non InterpScheme methods are required to preserve functionality whilst testing new schemes
  */
 #pragma once
 
+#include <complex>
+
 /*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<typename T>
+        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