#216 iteratefdtd_matrix: Updates for the main loop

tdms/include/field.h

@@ -61,9 +61,7 @@ class SplitFieldComponent : public Tensor3D<double> {
   fftw_plan *plan_b = nullptr;// Backward fftw plan
 
   double **operator[](int value) const { return tensor[value]; };
-  double operator[](CellCoordinate cell) const {
-    return tensor[cell.k][cell.j][cell.i];
-  }
+  double &operator[](ijk cell) const { return tensor[cell.k][cell.j][cell.i]; }
 
   void initialise_from_matlab(double ***tensor, Dimensions &dims);
 

tdms/include/simulation_manager/objects_from_infile.h

@@ -75,10 +75,10 @@ class IndependentObjectsFromInfile {
 
   /* DERIVED VARIABLES FROM INDEPENDENT INPUTS */
 
-  IJKDimensions IJK_tot;//< total number of Yee cells in the x,y,z directions
-                        // respectively
-  int Nsteps;           //< Number of dfts to perform before checking for phasor
-                        // convergence
+  IJKDimensions IJK_tot;//!< total number of Yee cells in the x,y,z directions
+                        //!< respectively
+  int Nsteps;//!< Number of dfts to perform before checking for phasor
+             //!< convergence
 
   IndependentObjectsFromInfile(
           InputMatrices matrices_from_input_file,
@@ -140,5 +140,11 @@ class ObjectsFromInfile : public IndependentObjectsFromInfile {
                     PreferredInterpolationMethods _pim =
                             PreferredInterpolationMethods::BandLimited);
 
+  /** @brief Determine whether the {IJK}source terms are empty (true) or not
+   * (false) */
+  bool all_sources_are_empty() {
+    return Isource.is_empty() && Jsource.is_empty() && Ksource.is_empty();
+  }
+
   ~ObjectsFromInfile();
 };

tdms/include/simulation_manager/simulation_manager.h

@@ -31,22 +31,26 @@
  */
 class SimulationManager {
 private:
-  ObjectsFromInfile
-          inputs;    //< The input objects that are generated from an input file
-  LoopTimers timers; //< Timers for tracking the execution of the simulation
-  PSTDVariables PSTD;//< PSTD-solver-specific variables
-  FDTDBootstrapper FDTD; //< FDTD bootstrapping variables
-  OutputMatrices outputs;//< Output object that will contain the results of this
-                         // simulation, given the input file
-
-  PreferredInterpolationMethods
-          pim;//< The interpolation methods to use in this simulation
-  SolverMethod solver_method;//< The solver method to use in this simulation
+  /*! The input objects that are generated from an input file */
+  ObjectsFromInfile inputs;
+
+  LoopTimers timers;    //!< Timers for tracking the execution of the simulation
+  PSTDVariables PSTD;   //!< PSTD-solver-specific variables
+  FDTDBootstrapper FDTD;//!< FDTD bootstrapping variables
+
+  /*! Output object that will contain the results of this simulation, given the
+   * input file */
+  OutputMatrices outputs;
+
+  /*! The solver method to use in this simulation */
+  SolverMethod solver_method;
+  /*! The interpolation methods to use in this simulation */
+  PreferredInterpolationMethods pim;
 
   EHVec eh_vec;//!< TODO
 
-  double ramp_width = 4.;//< Width of the ramp when introducing the waveform in
-                         // steady state mode
+  /*! Width of the ramp when introducing the waveform in steady state mode */
+  double ramp_width = 4.;
   /**
    * @brief Evaluates the linear ramp function r(t)
    *
@@ -65,10 +69,10 @@ class SimulationManager {
     return std::min(1., t / (ramp_width * period));
   }
 
-  std::vector<std::complex<double>>
-          E_norm;//< Holds the E-field phasors norm at each extraction frequency
-  std::vector<std::complex<double>>
-          H_norm;//< Holds the H-field phasors norm at each extraction frequency
+  /*! Holds the {E,H}-field phasors norm at each extraction frequency */
+  std::vector<std::complex<double>> E_norm;
+  /*! @copydoc E_norm */
+  std::vector<std::complex<double>> H_norm;
   /**
    * @brief Extracts the phasors norms at the given frequency (index)
    *
@@ -88,6 +92,122 @@ class SimulationManager {
    */
   void prepare_output(const mxArray *fieldsample, const mxArray *campssample);
 
+  /**
+   * @brief Update electric-split field components and current densities AT A
+   * PARTICULAR CELL in steady-state after an E-field timestep has been
+   * performed.
+   *
+   * @param time_H The time the magnetic field is currently sitting at.
+   * @param parallel The axis to which the plane of the Source term is parallel.
+   * EG X = Isource, Y = Jsource, etc
+   * @param C_axis Whether we are updating terms along the C-axis (true) or
+   * B-axis (false) of the Source plane. See Source doc for axis information.
+   * @param zero_plane Whether we are updating terms on the 0-plane (true) or
+   * the 1-plane (false). EG I0 would have this input as true, whereas J1 as
+   * false.
+   * @param is_conductive Whether the medium is conductive (so J_c needs to be
+   * updated)
+   * @param cell_b,cell_c The coordinates (cell_a, cell_b, cell_c) of the Yee
+   * cell in which we are updating the field. See the Source doc for notation
+   * information.
+   * @param array_ind The index of the various material property arrays that
+   * correspond to this particular Yee cell, and thus update equation.
+   * @param J_c The current density in the conductive medium
+   * @param J_s The current density in the dispersive medium
+   */
+  void E_source_update_steadystate(double time_H, AxialDirection parallel,
+                                   bool C_axis, bool zero_plane,
+                                   bool is_conductive, int cell_b, int cell_c,
+                                   int array_ind, CurrentDensitySplitField &J_c,
+                                   CurrentDensitySplitField &J_s);
+  /**
+   * @brief Update magnetic-split field components AT A
+   * PARTICULAR CELL in steady-state after an H-field timestep has been
+   * performed.
+   *
+   * @param time_E The time the electric field is currently sitting at
+   * @param parallel The axis to which the plane of the Source term is parallel.
+   * EG X = Isource, Y = Jsource, etc
+   * @param C_axis Whether we are updating terms along the C-axis (true) or
+   * B-axis (false) of the Source plane. See Source doc for axis information.
+   * @param zero_plane Whether we are updating terms on the 0-plane (true) or
+   * the 1-plane (false). EG I0 would have this input as true, whereas J1 as
+   * false.
+   * @param array_ind The index of the various material property arrays that
+   * correspond to this particular Yee cell, and thus update equation.
+   * @param cell_b,cell_c The coordinates (cell_a, cell_b, cell_c) of the Yee
+   * cell in which we are updating the field. See the Source doc for notation
+   * information.
+   */
+  void H_source_update_steadystate(double time_E, AxialDirection parallel,
+                                   bool zero_plane, bool C_axis, int array_ind,
+                                   int cell_b, int cell_c);
+  /**
+   * @brief [E-FIELD UPDATES] Performs updates to the electric-split field
+   * components and current density fields after an E-field timestep has been
+   * performed, in accordance with the I,J, and K-source terms.
+   *
+   * @param time_H The time the magnetic field is currently sitting at.
+   * @param is_conductive Whether the medium is conductive (so J_c needs to be
+   * updated)
+   * @param J_c The current density in the conductive medium
+   * @param J_s The current density in the dispersive medium
+   */
+  void E_source_update_all_steadystate(double time_H, bool is_conductive,
+                                       CurrentDensitySplitField &J_c,
+                                       CurrentDensitySplitField &J_s);
+  /**
+   * @brief [H-FIELD UPDATES] Performs updates to the magnetic-split field
+   * components after an H-field timestep has been performed, in accordance with
+   * the I,J, and K-source terms.
+   *
+   * @param time_E The time the electric field is currently sitting at.
+   */
+  void H_source_update_all_steadystate(double time_E);
+
+  /*! @copydoc E_source_update_all_steadystate */
+  void E_Isource_update_steadystate(double time_H, bool is_conductive,
+                                    CurrentDensitySplitField &J_c,
+                                    CurrentDensitySplitField &J_s);
+  /*! @copydoc E_source_update_all_steadystate */
+  void E_Jsource_update_steadystate(double time_H, bool is_conductive,
+                                    CurrentDensitySplitField &J_c,
+                                    CurrentDensitySplitField &J_s);
+  /*! @copydoc E_source_update_all_steadystate */
+  void E_Ksource_update_steadystate(double time_H, bool is_conductive,
+                                    CurrentDensitySplitField &J_c,
+                                    CurrentDensitySplitField &J_s);
+  /*! @copydoc H_source_update_all_steadystate */
+  void H_Isource_update_steadystate(double time_E);
+  /*! @copydoc H_source_update_all_steadystate */
+  void H_Jsource_update_steadystate(double time_E);
+  /*! @copydoc H_source_update_all_steadystate */
+  void H_Ksource_update_steadystate(double time_E);
+
+  /**
+   * @brief [E-FIELD UPDATES] Performs the updates to the electric-split field
+   * components nad current density fields after an E-field timestep has been
+   * performed, in accordance with the source terms.
+   *
+   * @param time_H The time the magnetic field is currently sitting at.
+   * @param is_conductive Whether the medium is conductive (so J_c needs to be
+   * updated)
+   * @param J_c The current density in the conductive medium
+   * @param J_s The current density in the dispersive medium
+   */
+  void update_source_terms_pulsed(double time_H, bool is_conductive,
+                                  CurrentDensitySplitField &J_c,
+                                  CurrentDensitySplitField &J_s);
+  /**
+   * @brief [H-FIELD UPDATES] Performs the updates to the magnetic-split field
+   * components and current density fields after an H-field timestep has been
+   * performed, in accordance with the source terms.
+   *
+   * @param time_E The time the electric field is currently sitting at.
+   * @param tind The current iteration number
+   */
+  void update_source_terms_pulsed(double time_E, int tind);
+
 public:
   SimulationManager(InputMatrices in_matrices, SolverMethod _solver_method,
                     PreferredInterpolationMethods _pim);

tdms/include/source.h

@@ -3,14 +3,65 @@
  */
 #pragma once
 
+#include <complex>
 #include <string>
 
+#include "cell_coordinate.h"
 #include "mat_io.h"
 
+/**
+ * @brief Typedef for providing indices to the Source class.
+ *
+ * Assume the notation as in the docstring of the Source class. Recall that data
+ * in a Source instance is indexed via
+ * {real,imag}[cell_c][cell_b][split_field_ID].
+ *
+ * SourceIndex is a container for indexing the Source.{real,imag} objects in an
+ * efficient and notationally-consistent manner.
+ */
+struct SourceIndex {
+  int split_field_ID = 0,//!< Index of the split field
+          cell_b = 0,    //!< Index in the minor (B-)axis
+          cell_c = 0;    //!< Index in the major (C-)axis
+};
+
+/**
+ * @brief The Source class stores values of the Source field across a particular
+ * plane.
+ *
+ * Let A (= {i,j,k}) be the axial direction that the plane the given Source
+ * instance is storing data for. Let B, C be the remaining axial directions,
+ * with C being the axial direction with the slower-varying index. That is: A =
+ * i : B = j : C = k, A = j : B = i : C = k, A = k : B = i : C = j.
+ *
+ * The Source data (real and imag) is indexed by 3 indices, accessed via
+ * {real,imag}[cell_c][cell_b][split_field_ID].
+ *
+ * split_field_ID ranges between 0-7 inclusive.
+ * TODO: Indices <-> sources need to be deduced from input file generator
+ * functions.
+ *
+ * Let cell_a be the A-index of the plane that the instance is storing data
+ * on. Then (cell_A, cell_B, cell_C) is the Yee-cell index whose (source) data
+ * we are accessing with this call.
+ */
 class Source {
+private:
+  /*! Flags if the array is empty to avoid pointer preservation */
+  bool no_data_stored = true;
+
 public:
-  double ***real = nullptr;
-  double ***imag = nullptr;
+  double ***real = nullptr;//!< Real data for the source term
+  double ***imag = nullptr;//!< Imag data for the source term
 
   Source(const mxArray *ptr, int dim1, int dim2, const std::string &name);
+
+  /** @brief Check if the source term is empty (true) or not (false) */
+  bool is_empty() { return no_data_stored; }
+
+  std::complex<double> operator[](SourceIndex index) {
+    return std::complex<double>(
+            real[index.cell_c][index.cell_b][index.split_field_ID],
+            imag[index.cell_c][index.cell_b][index.split_field_ID]);
+  }
 };