From 49f4f4c0ea72afd9ee081f734501f50152aa5399 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Tue, 31 Jan 2023 16:09:15 +0000 Subject: [PATCH 01/25] Cherry-pick Sam's linting --- .pre-commit-config.yaml | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b2725b857..23fbfa058 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ repos: # Sort order of Python imports - repo: https://github.com/pycqa/isort - rev: 5.10.1 + rev: 5.12.0 hooks: - id: isort # Python code formatting @@ -21,3 +21,18 @@ repos: hooks: - id: trailing-whitespace - id: end-of-file-fixer + # Add C++ linting + - repo: https://github.com/pocc/pre-commit-hooks + rev: v1.3.5 + hooks: + - id: clang-format + args: ["--style=file"] + - id: clang-tidy + args: ["-p=tdms/build/compile_commands.json"] + # - id: cppcheck + # args: ["--project=tdms/build/compile_commands.json"] + - id: cpplint + args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] + # suppress whitespace because it conflicts with clang-tidy and doxygen + # suppress copyright in file headers because we don't require them in the developer docs + # run quietly so we don't see you doing this if everything is fine From c16f9faee17d3a8d88557a7dcace28ec1f78e1fd Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Tue, 31 Jan 2023 16:23:36 +0000 Subject: [PATCH 02/25] Add flag to Source class to indicate empty array --- .pre-commit-config.yaml | 6 +- .../simulation_manager/objects_from_infile.h | 59 ++++++++++--------- tdms/include/source.h | 11 +++- tdms/src/source.cpp | 29 +++++---- 4 files changed, 59 insertions(+), 46 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 23fbfa058..c78c25b90 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -26,13 +26,13 @@ repos: rev: v1.3.5 hooks: - id: clang-format - args: ["--style=file"] + args: ["--style=file", "-i"] - id: clang-tidy args: ["-p=tdms/build/compile_commands.json"] # - id: cppcheck # args: ["--project=tdms/build/compile_commands.json"] - - id: cpplint - args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] + # - id: cpplint + # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] # suppress whitespace because it conflicts with clang-tidy and doxygen # suppress copyright in file headers because we don't require them in the developer docs # run quietly so we don't see you doing this if everything is fine diff --git a/tdms/include/simulation_manager/objects_from_infile.h b/tdms/include/simulation_manager/objects_from_infile.h index c1aabd2c0..e93a1aa2a 100644 --- a/tdms/include/simulation_manager/objects_from_infile.h +++ b/tdms/include/simulation_manager/objects_from_infile.h @@ -9,9 +9,9 @@ #include "matrix.h" #include "arrays.h" +#include "cell_coordinate.h" #include "field.h" #include "fieldsample.h" -#include "cell_coordinate.h" #include "globals.h" #include "grid_labels.h" #include "input_matrices.h" @@ -32,20 +32,20 @@ class IndependentObjectsFromInfile { public: SolverMethod solver_method;//< Either PSTD (default) or FDTD, the solver method - int skip_tdf; //< Either 1 if we are using PSTD, or 6 if using FDTD + int skip_tdf; //< Either 1 if we are using PSTD, or 6 if using FDTD PreferredInterpolationMethods interpolation_methods = - BandLimited; //< Either band_limited or cubic, the preferred interpolation methods + BandLimited;//< Either band_limited or cubic, the preferred interpolation methods - SimulationParameters params; //< The parameters for this simulation + SimulationParameters params;//< The parameters for this simulation - ElectricSplitField E_s; //< The split electric-field values - MagneticSplitField H_s; //< The split magnetic-field values + ElectricSplitField E_s;//< The split electric-field values + MagneticSplitField H_s;//< The split magnetic-field values - uint8_t ***materials; //< TODO - CMaterial Cmaterial; //< TODO - DMaterial Dmaterial; //< TODO - CCollection C; //< TODO - DCollection D; //< TODO + uint8_t ***materials;//< TODO + CMaterial Cmaterial; //< TODO + DMaterial Dmaterial; //< TODO + CCollection C; //< TODO + DCollection D; //< TODO double *freespace_Cbx; //< freespace constants double *alpha, *beta, *gamma;//< dispersion parameters @@ -53,22 +53,22 @@ class IndependentObjectsFromInfile { InterfaceComponent I0, I1, J0, J1, K0, K1;//< user-defined interface components Cuboid cuboid; //< user-defined surface to extract phasors over - XYZVectors rho_cond; //< conductive aux - DispersiveMultiLayer matched_layer; //< dispersive aux + XYZVectors rho_cond; //< conductive aux + DispersiveMultiLayer matched_layer;//< dispersive aux - IncidentField Ei; //< time-domain field + IncidentField Ei;//< time-domain field - FrequencyVectors f_vec; //< frequency vector - Pupil pupil; //< TODO - DTilde D_tilde; //< TODO - TDFieldExporter2D ex_td_field_exporter; //< two-dimensional field exporter + FrequencyVectors f_vec; //< frequency vector + Pupil pupil; //< TODO + DTilde D_tilde; //< TODO + TDFieldExporter2D ex_td_field_exporter;//< two-dimensional field exporter GridLabels input_grid_labels;//< cartesian labels of the Yee cells /* 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 + int Nsteps; //< Number of dfts to perform before checking for phasor convergence IndependentObjectsFromInfile( InputMatrices matrices_from_input_file, @@ -113,14 +113,19 @@ class IndependentObjectsFromInfile { */ class ObjectsFromInfile : public IndependentObjectsFromInfile { public: - Source Isource, Jsource, Ksource;//< TODO - GratingStructure structure;//< TODO - FrequencyExtractVector f_ex_vec;//< Vector of frequencies to extract field & phasors at + Source Isource, Jsource, Ksource;//< TODO + GratingStructure structure; //< TODO + FrequencyExtractVector f_ex_vec; //< Vector of frequencies to extract field & phasors at - ObjectsFromInfile( - InputMatrices matrices_from_input_file, - SolverMethod _solver_method = SolverMethod::PseudoSpectral, - PreferredInterpolationMethods _pim = PreferredInterpolationMethods::BandLimited); + ObjectsFromInfile( + InputMatrices matrices_from_input_file, + SolverMethod _solver_method = SolverMethod::PseudoSpectral, + 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(); + ~ObjectsFromInfile(); }; diff --git a/tdms/include/source.h b/tdms/include/source.h index 8d4c4afbe..bc440c593 100644 --- a/tdms/include/source.h +++ b/tdms/include/source.h @@ -7,10 +7,15 @@ #include "mat_io.h" -class Source{ +class Source { +private: + bool no_data_stored = true;//!< Flags if the array is empty to avoid pointer preservation 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; } }; diff --git a/tdms/src/source.cpp b/tdms/src/source.cpp index b80317bce..b20881779 100644 --- a/tdms/src/source.cpp +++ b/tdms/src/source.cpp @@ -1,36 +1,39 @@ #include "source.h" -#include #include #include +#include -#include "matlabio.h" #include "dimensions.h" +#include "matlabio.h" using namespace std; -Source::Source(const mxArray *ptr, int dim1, int dim2, const std::string &name){ +Source::Source(const mxArray *ptr, int dim1, int dim2, const std::string &name) { if (mxIsEmpty(ptr)) { spdlog::info("{} is empty", name); } else { + // fetch dimensions of the input array, and check they are what we expect auto dims = Dimensions(ptr); - if (dims.are_1d()){ - throw runtime_error(name+" should be 3- or 2-dimensional"); - } - if (dims.are_2d()){ - dim2 = 0; - } - if (!(dims[0] == 8 && dims[1] == dim1 && dims[2] == dim2)){ + if (dims.are_1d()) { throw runtime_error(name + " should be 3- or 2-dimensional"); } + if (dims.are_2d()) { dim2 = 0; } + if (!(dims[0] == 8 && dims[1] == dim1 && dims[2] == dim2)) { cerr << name << " has incorrect size" << endl; } - if (!mxIsComplex(ptr)){ - throw runtime_error(name+" should be complex, use a call of " - "complex(real(Isource),imag(Isource)) in matlab if necessary"); + + // check that complex data has been passed + if (!mxIsComplex(ptr)) { + throw runtime_error(name + " should be complex, use a call of " + "complex(real(Isource),imag(Isource)) in matlab if necessary"); } + + // cast MATLAB arrays to C++ datatypes real = cast_matlab_3D_array(mxGetPr(ptr), dims[0], dims[1], dims[2]); imag = cast_matlab_3D_array(mxGetPi(ptr), dims[0], dims[1], dims[2]); + // flag Source as non-empty + no_data_stored = false; } } From d11dc2e5aaba50b78927b83eb600b2fdaa5cc8a7 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Tue, 31 Jan 2023 16:24:14 +0000 Subject: [PATCH 03/25] Add flag to Source class to indicate empty array --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c78c25b90..37acc617b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -31,8 +31,8 @@ repos: args: ["-p=tdms/build/compile_commands.json"] # - id: cppcheck # args: ["--project=tdms/build/compile_commands.json"] - # - id: cpplint - # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] + - id: cpplint + args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] # suppress whitespace because it conflicts with clang-tidy and doxygen # suppress copyright in file headers because we don't require them in the developer docs # run quietly so we don't see you doing this if everything is fine From 0a9b35af8f56456fdee91f216ec29d1ed73c4f4b Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Tue, 31 Jan 2023 16:25:57 +0000 Subject: [PATCH 04/25] Indentation --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 37acc617b..c78c25b90 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -31,8 +31,8 @@ repos: args: ["-p=tdms/build/compile_commands.json"] # - id: cppcheck # args: ["--project=tdms/build/compile_commands.json"] - - id: cpplint - args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] + # - id: cpplint + # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] # suppress whitespace because it conflicts with clang-tidy and doxygen # suppress copyright in file headers because we don't require them in the developer docs # run quietly so we don't see you doing this if everything is fine From fb5c0b5ebd049ccb8792984b58bbe79b571a4256 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Tue, 31 Jan 2023 16:30:45 +0000 Subject: [PATCH 05/25] Initialise values rather than assign --- .../simulation_manager/simulation_manager.cpp | 48 +++++++++---------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/tdms/src/simulation_manager/simulation_manager.cpp b/tdms/src/simulation_manager/simulation_manager.cpp index bab5df0db..edd78ebd8 100644 --- a/tdms/src/simulation_manager/simulation_manager.cpp +++ b/tdms/src/simulation_manager/simulation_manager.cpp @@ -9,11 +9,8 @@ using namespace tdms_math_constants; SimulationManager::SimulationManager(InputMatrices in_matrices, SolverMethod _solver_method, PreferredInterpolationMethods _pim) - : inputs(in_matrices, _solver_method, _pim), - FDTD(n_Yee_cells()) { - solver_method = _solver_method; - pim = _pim; - + : inputs(in_matrices, _solver_method, _pim), FDTD(n_Yee_cells()), solver_method(_solver_method), + pim(_pim) { // read number of Yee cells IJKDimensions IJK_tot = n_Yee_cells(); @@ -70,27 +67,30 @@ void SimulationManager::prepare_output(const mxArray *fieldsample, const mxArray } void SimulationManager::post_loop_processing() { - // normalise the phasors in the volume, if we extracted there - if (inputs.params.run_mode == RunMode::complete && inputs.params.exphasorsvolume) { - outputs.E.normalise_volume(); - outputs.H.normalise_volume(); - } + // we only normalise if we were provided with at least one non-empty source term + if (!inputs.all_sources_are_empty()) { + // normalise the phasors in the volume, if we extracted there + if (inputs.params.run_mode == RunMode::complete && inputs.params.exphasorsvolume) { + outputs.E.normalise_volume(); + outputs.H.normalise_volume(); + } - // normalise the phasors on the surface, if we extracted there - if (inputs.params.run_mode == RunMode::complete && inputs.params.exphasorssurface) { - spdlog::info("Surface phasors"); - for (int ifx = 0; ifx < inputs.f_ex_vec.size(); ifx++) { - outputs.surface_phasors.normalise_surface(ifx, E_norm[ifx], H_norm[ifx]); - spdlog::info("\tE_norm[{0:d}]: {1:.5e} {2:.5e}", ifx, real(E_norm[ifx]), imag(E_norm[ifx])); + // normalise the phasors on the surface, if we extracted there + if (inputs.params.run_mode == RunMode::complete && inputs.params.exphasorssurface) { + spdlog::info("Surface phasors"); + for (int ifx = 0; ifx < inputs.f_ex_vec.size(); ifx++) { + outputs.surface_phasors.normalise_surface(ifx, E_norm[ifx], H_norm[ifx]); + spdlog::info("\tE_norm[{0:d}]: {1:.5e} {2:.5e}", ifx, real(E_norm[ifx]), imag(E_norm[ifx])); + } } - } - // normalise the phasors on the vertices, if we extracted at any - if (inputs.params.run_mode == RunMode::complete && - outputs.vertex_phasors.there_are_vertices_to_extract_at()) { - spdlog::info("Vertex phasors"); - for (int ifx = 0; ifx < inputs.f_ex_vec.size(); ifx++) { - outputs.vertex_phasors.normalise_vertices(ifx, E_norm[ifx], H_norm[ifx]); - spdlog::info("\tE_norm[{0:d}]: {1:.5e} {2:.5e}", ifx, real(E_norm[ifx]), imag(E_norm[ifx])); + // normalise the phasors on the vertices, if we extracted at any + if (inputs.params.run_mode == RunMode::complete && + outputs.vertex_phasors.there_are_vertices_to_extract_at()) { + spdlog::info("Vertex phasors"); + for (int ifx = 0; ifx < inputs.f_ex_vec.size(); ifx++) { + outputs.vertex_phasors.normalise_vertices(ifx, E_norm[ifx], H_norm[ifx]); + spdlog::info("\tE_norm[{0:d}]: {1:.5e} {2:.5e}", ifx, real(E_norm[ifx]), imag(E_norm[ifx])); + } } } From 61651cb77a47f09962f07dee6c235e28f52a4e7c Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 09:46:29 +0000 Subject: [PATCH 06/25] Refactored source updates - flagged mystery indexing --- .../simulation_manager/simulation_manager.h | 38 +- .../simulation_manager/execute_simulation.cpp | 558 ++++-------------- .../source_term_updates.cpp | 451 ++++++++++++++ 3 files changed, 608 insertions(+), 439 deletions(-) create mode 100644 tdms/src/simulation_manager/source_term_updates.cpp diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index cda31ff74..f98b33126 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -5,8 +5,8 @@ #pragma once #include -#include #include +#include #include "arrays.h" #include "cell_coordinate.h" @@ -26,14 +26,15 @@ */ 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 + 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 + PreferredInterpolationMethods pim;//< The interpolation methods to use in this simulation + SolverMethod solver_method; //< The solver method to use in this simulation EHVec eh_vec;//< TODO @@ -56,8 +57,10 @@ class SimulationManager { return std::min(1., t / (ramp_width * period)); } - std::vector> E_norm;//< Holds the E-field phasors norm at each extraction frequency - std::vector> H_norm;//< Holds the H-field phasors norm at each extraction frequency + std::vector> + E_norm;//< Holds the E-field phasors norm at each extraction frequency + std::vector> + H_norm;//< Holds the H-field phasors norm at each extraction frequency /** * @brief Extracts the phasors norms at the given frequency (index) * @@ -73,10 +76,21 @@ class SimulationManager { * @param fieldsample The fieldsample input from the input file * @param campssample The complex amplitude sample input from the input file */ - void prepare_output(const mxArray *fieldsample, const mxArray* campssample); + void prepare_output(const mxArray *fieldsample, const mxArray *campssample); + + void update_Isource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s); + void update_Jsource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s); + void update_Ksource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s); public: - SimulationManager(InputMatrices in_matrices, SolverMethod _solver_method, PreferredInterpolationMethods _pim); + SimulationManager(InputMatrices in_matrices, SolverMethod _solver_method, + PreferredInterpolationMethods _pim); /** @brief Fetch the number of Yee cells in each dimension */ IJKDimensions n_Yee_cells() { return inputs.IJK_tot; } diff --git a/tdms/src/simulation_manager/execute_simulation.cpp b/tdms/src/simulation_manager/execute_simulation.cpp index 1e33965e3..05a6e6bd9 100644 --- a/tdms/src/simulation_manager/execute_simulation.cpp +++ b/tdms/src/simulation_manager/execute_simulation.cpp @@ -7,8 +7,8 @@ #include #include -#include "simulation_manager/loop_variables.h" #include "numerical_derivative.h" +#include "simulation_manager/loop_variables.h" using namespace std; using namespace tdms_math_constants; @@ -31,9 +31,10 @@ void SimulationManager::execute() { // DECLARE VARIABLES SCOPED TO THIS FUNCTION ONLY double rho; - double alpha_l, beta_l, gamma_l;//< alpha, beta, gamma parameters of the layer the local thread is examining + double alpha_l, beta_l, + gamma_l;//< alpha, beta, gamma parameters of the layer the local thread is examining double kappa_l, sigma_l;//< kappa, sigma parameters of the layer the local thread is examining - double t0;//< (Real) time since the last iteration log was written to the screen + double t0; //< (Real) time since the last iteration log was written to the screen double Ca, Cb, Cc;//used by interpolation scheme //the C and D vars for free space and pml @@ -46,10 +47,11 @@ void SimulationManager::execute() { double lambda_an_t;//< Wavelength of light in free space at the current frequency int i, j, k;//< Loop variables - int n;//< The thread number of the local OMP thread - int k_loc;//< Local thread copy of the variable k + int n; //< The thread number of the local OMP thread + int k_loc; //< Local thread copy of the variable k - int dft_counter = 0;//< Number of DFTs we have performed since last checking for phasor convergence + int dft_counter = + 0;//< Number of DFTs we have performed since last checking for phasor convergence complex Idxt, Idyt, kprop; @@ -315,8 +317,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -443,8 +447,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -588,8 +594,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -710,8 +718,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -846,8 +856,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) { - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -973,8 +985,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) { - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1117,8 +1131,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1238,8 +1254,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1379,8 +1397,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1505,8 +1525,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1648,8 +1670,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1743,8 +1767,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -1872,8 +1898,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2017,8 +2045,10 @@ void SimulationManager::execute() { rho = 0.; k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2105,373 +2135,19 @@ void SimulationManager::execute() { /********************/ //update terms for self consistency across scattered/total interface - E updates## - if (inputs.params.source_mode == SourceMode::steadystate) {//steadystate + + if (inputs.params.source_mode == SourceMode::steadystate) { + // Steady-state source term updates + update_Isource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, + loop_variables.J_s); + update_Jsource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, + loop_variables.J_s); + update_Ksource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, + loop_variables.J_s); + complex commonPhase = exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); double commonAmplitude = linear_ramp(time_H); - for (k = (inputs.K0.index); k <= (inputs.K1.index); k++) - for (j = (inputs.J0.index); j <= (inputs.J1.index); j++) { - if (inputs.I0.apply) {//Perform across I0 - - if (!inputs.params.is_multilayer) array_ind = inputs.I0.index; - else - array_ind = (I_tot + 1) * k + inputs.I0.index; - - if (k < (inputs.K1.index) || - inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zx[k][j][inputs.I0.index] = - inputs.E_s.zx[k][j][inputs.I0.index] - - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); - if (loop_variables.is_conductive) - loop_variables.J_c.zx[k][j][inputs.I0.index] += - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.zx[k][j][inputs.I0.index] += - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); - } - if (j < (inputs.J1.index)) { - inputs.E_s.yx[k][j][inputs.I0.index] = - inputs.E_s.yx[k][j][inputs.I0.index] + - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); - if (loop_variables.is_conductive) - loop_variables.J_c.yx[k][j][inputs.I0.index] -= - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.yx[k][j][inputs.I0.index] -= - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); - } - } - if (inputs.I1.apply) {//Perform across I1 - - if (!inputs.params.is_multilayer) array_ind = inputs.I1.index; - else - array_ind = (I_tot + 1) * k + inputs.I1.index; - - if (k < (inputs.K1.index) || - inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zx[k][j][inputs.I1.index] = - inputs.E_s.zx[k][j][inputs.I1.index] + - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); - if (loop_variables.is_conductive) - loop_variables.J_c.zx[k][j][inputs.I1.index] -= - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.zx[k][j][inputs.I1.index] -= - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); - } - if (j < (inputs.J1.index)) { - inputs.E_s.yx[k][j][inputs.I1.index] = - inputs.E_s.yx[k][j][inputs.I1.index] - - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); - if (loop_variables.is_conductive) - loop_variables.J_c.yx[k][j][inputs.I1.index] += - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.yx[k][j][inputs.I1.index] += - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); - } - } - } - - for (k = (inputs.K0.index); k <= (inputs.K1.index); k++) - for (i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - if (inputs.J0.apply) {//Perform across J0 - if (k < (inputs.K1.index) || - inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - - if (!inputs.params.is_multilayer) array_ind = inputs.J0.index; - else - array_ind = (J_tot + 1) * k + inputs.J0.index; - - inputs.E_s.zy[k][(inputs.J0.index)][i] = - inputs.E_s.zy[k][(inputs.J0.index)][i] + - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); - if (loop_variables.is_conductive) - loop_variables.J_c.zy[k][(inputs.J0.index)][i] -= - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.zy[k][(inputs.J0.index)][i] -= - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); - } - if (i < (inputs.I1.index)) { - inputs.E_s.xy[k][(inputs.J0.index)][i] = - inputs.E_s.xy[k][(inputs.J0.index)][i] - - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); - if (loop_variables.is_conductive) - loop_variables.J_c.xy[k][(inputs.J0.index)][i] += - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.xy[k][(inputs.J0.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); - } - } - if (inputs.J1.apply) {//Perform across J1 - - if (!inputs.params.is_multilayer) array_ind = inputs.J1.index; - else - array_ind = (J_tot + 1) * k + inputs.J1.index; - - if (k < (inputs.K1.index) || - inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zy[k][(inputs.J1.index)][i] = - inputs.E_s.zy[k][(inputs.J1.index)][i] - - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); - if (loop_variables.is_conductive) - loop_variables.J_c.zy[k][(inputs.J1.index)][i] += - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.zy[k][(inputs.J1.index)][i] -= - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); - } - if (i < (inputs.I1.index)) { - inputs.E_s.xy[k][(inputs.J1.index)][i] = - inputs.E_s.xy[k][(inputs.J1.index)][i] + - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); - if (loop_variables.is_conductive) - loop_variables.J_c.xy[k][(inputs.J1.index)][i] -= - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.xy[k][(inputs.J1.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); - } - } - } - - for (j = (inputs.J0.index); j <= (inputs.J1.index); j++) - for (i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - if (inputs.K0.apply) {//Perform across K0 - if (j < (inputs.J1.index)) { - inputs.E_s.yz[(inputs.K0.index)][j][i] = - inputs.E_s.yz[(inputs.K0.index)][j][i] - - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); - if (loop_variables.is_conductive) - loop_variables.J_c.yz[(inputs.K0.index)][j][i] += - inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.yz[(inputs.K0.index)][j][i] -= - inputs.matched_layer.kappa.z[(inputs.K0.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); - } - if (i < (inputs.I1.index)) { - inputs.E_s.xz[(inputs.K0.index)][j][i] = - inputs.E_s.xz[(inputs.K0.index)][j][i] + - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); - if (loop_variables.is_conductive) - loop_variables.J_c.xz[(inputs.K0.index)][j][i] -= - inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.xz[(inputs.K0.index)][j][i] += - inputs.matched_layer.kappa.z[(inputs.K0.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); - } - } - if (inputs.K1.apply) {//Perform across K1 - if (j < (inputs.J1.index)) { - inputs.E_s.yz[(inputs.K1.index)][j][i] = - inputs.E_s.yz[(inputs.K1.index)][j][i] + - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); - if (loop_variables.is_conductive) - loop_variables.J_c.yz[(inputs.K1.index)][j][i] -= - inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.yz[(inputs.K1.index)][j][i] += - inputs.matched_layer.kappa.z[(inputs.K1.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); - } - if (i < (inputs.I1.index)) { - inputs.E_s.xz[(inputs.K1.index)][j][i] = - inputs.E_s.xz[(inputs.K1.index)][j][i] - - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); - if (loop_variables.is_conductive) - loop_variables.J_c.xz[(inputs.K1.index)][j][i] += - inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); - if (inputs.params.is_disp_ml) - loop_variables.J_s.xz[(inputs.K1.index)][j][i] -= - inputs.matched_layer.kappa.z[(inputs.K1.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); - } - } - } outputs.H.ft = real(commonAmplitude * commonPhase); } else if (inputs.params.source_mode == SourceMode::pulsed) {//pulsed @@ -2665,8 +2341,10 @@ void SimulationManager::execute() { for (i = 0; i < (I_tot + 1); i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2699,8 +2377,10 @@ void SimulationManager::execute() { for (k = 0; k < K_tot; k++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2748,8 +2428,10 @@ void SimulationManager::execute() { for (i = 0; i < (I_tot + 1); i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2784,8 +2466,10 @@ void SimulationManager::execute() { for (j = 0; j < J_tot; j++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2834,8 +2518,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2871,8 +2557,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2921,8 +2609,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -2957,8 +2647,10 @@ void SimulationManager::execute() { for (k = 0; k < K_tot; k++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3011,8 +2703,10 @@ void SimulationManager::execute() { for (i = 0; i < (I_tot + 1); i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3046,8 +2740,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3095,8 +2791,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3131,8 +2829,10 @@ void SimulationManager::execute() { for (j = 0; j < J_tot; j++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3182,8 +2882,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3218,8 +2920,10 @@ void SimulationManager::execute() { for (i = 0; i < I_tot; i++) { k_loc = k; if (inputs.params.is_structure) - if (k > inputs.params.pml.Dzl && k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { - if ((k - inputs.structure[i][1]) < (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && + if (k > inputs.params.pml.Dzl && + k < (inputs.params.pml.Dzl + loop_variables.n_non_pml_cells_in_K)) { + if ((k - inputs.structure[i][1]) < + (loop_variables.n_non_pml_cells_in_K + inputs.params.pml.Dzl) && (k - inputs.structure[i][1]) > inputs.params.pml.Dzl) k_loc = k - inputs.structure[i][1]; else if ((k - inputs.structure[i][1]) >= @@ -3590,7 +3294,7 @@ void SimulationManager::execute() { outputs.H.add_to_angular_norm(tind, inputs.Nsteps, inputs.params); for (int ifx = 0; ifx < inputs.f_ex_vec.size(); ifx++) { - extract_phasor_norms(ifx, tind, inputs.Nsteps); + extract_phasor_norms(ifx, tind, inputs.Nsteps); } } else { if ((tind - inputs.params.start_tind) % inputs.params.Np == 0) { diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_term_updates.cpp new file mode 100644 index 000000000..8c2432d21 --- /dev/null +++ b/tdms/src/simulation_manager/source_term_updates.cpp @@ -0,0 +1,451 @@ +#include "simulation_manager/simulation_manager.h" + +using namespace std; +using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; + +void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { + // Only run update equations is source data was provided + if (inputs.Isource.is_empty()) { return; } + + int array_ind; + int I_tot = n_Yee_cells().i; + + complex commonPhase = + exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + double commonAmplitude = linear_ramp(time_H); + + // Update across I0, provided a source term is here + if (inputs.I0.apply) { + for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + if (!inputs.params.is_multilayer) { + array_ind = inputs.I0.index; + } else { + array_ind = (I_tot + 1) * k + inputs.I0.index; + } + + if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + inputs.E_s.zx[k][j][inputs.I0.index] = + inputs.E_s.zx[k][j][inputs.I0.index] - + inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource + .real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2])); + if (is_conductive) { + J_c.zx[k][j][inputs.I0.index] += + inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2])); + } + if (inputs.params.is_disp_ml) { + J_s.zx[k][j][inputs.I0.index] += + inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2])); + } + } + + if (j < (inputs.J1.index)) { + inputs.E_s.yx[k][j][inputs.I0.index] = + inputs.E_s.yx[k][j][inputs.I0.index] + + inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource + .real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3])); + if (is_conductive) { + J_c.yx[k][j][inputs.I0.index] -= + inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3])); + } + if (inputs.params.is_disp_ml) { + J_s.yx[k][j][inputs.I0.index] -= + inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3])); + } + } + } + } + } + + // Update across I1, provided a source term is here + if (inputs.I1.apply) { + for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + if (!inputs.params.is_multilayer) { + array_ind = inputs.I1.index; + } else { + array_ind = (I_tot + 1) * k + inputs.I1.index; + } + + if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + inputs.E_s.zx[k][j][inputs.I1.index] = + inputs.E_s.zx[k][j][inputs.I1.index] + + inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource + .real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6])); + if (is_conductive) { + J_c.zx[k][j][inputs.I1.index] -= + inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6])); + } + if (inputs.params.is_disp_ml) { + J_s.zx[k][j][inputs.I1.index] -= + inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6])); + } + } + + if (j < (inputs.J1.index)) { + inputs.E_s.yx[k][j][inputs.I1.index] = + inputs.E_s.yx[k][j][inputs.I1.index] - + inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource + .real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7])); + if (is_conductive) { + J_c.yx[k][j][inputs.I1.index] += + inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7])); + } + if (inputs.params.is_disp_ml) { + J_s.yx[k][j][inputs.I1.index] += + inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7])); + } + } + } + } + } +} + +void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { + // Only run update equations is source data was provided + if (inputs.Jsource.is_empty()) { return; } + + int array_ind; + int J_tot = n_Yee_cells().j; + + complex commonPhase = + exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + double commonAmplitude = linear_ramp(time_H); + + // Update across J0, provided a source term is there + if (inputs.J0.apply) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int i = inputs.I0.index; i < inputs.I1.index; i++) { + if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + if (!inputs.params.is_multilayer) { + array_ind = inputs.J0.index; + } else { + array_ind = (J_tot + 1) * k + inputs.J0.index; + } + + inputs.E_s.zy[k][(inputs.J0.index)][i] = + inputs.E_s.zy[k][(inputs.J0.index)][i] + + inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource + .real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2])); + if (is_conductive) { + J_c.zy[k][(inputs.J0.index)][i] -= + inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2])); + } + if (inputs.params.is_disp_ml) { + J_s.zy[k][(inputs.J0.index)][i] -= + inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2])); + } + } + if (i < (inputs.I1.index)) { + inputs.E_s.xy[k][(inputs.J0.index)][i] = + inputs.E_s.xy[k][(inputs.J0.index)][i] - + inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource + .real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3])); + if (is_conductive) { + J_c.xy[k][(inputs.J0.index)][i] += + inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3])); + } + if (inputs.params.is_disp_ml) { + J_s.xy[k][(inputs.J0.index)][i] += + inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3])); + } + } + } + } + } + + // Update across J1, provided a source term is there + if (inputs.J1.apply) { + for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { + for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { + if (!inputs.params.is_multilayer) { + array_ind = inputs.J1.index; + } else { + array_ind = (J_tot + 1) * k + inputs.J1.index; + } + + if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + inputs.E_s.zy[k][(inputs.J1.index)][i] = + inputs.E_s.zy[k][(inputs.J1.index)][i] - + inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource + .real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6])); + if (is_conductive) { + J_c.zy[k][(inputs.J1.index)][i] += + inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6])); + } + if (inputs.params.is_disp_ml) { + J_s.zy[k][(inputs.J1.index)][i] -= + inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6])); + } + } + if (i < (inputs.I1.index)) { + inputs.E_s.xy[k][(inputs.J1.index)][i] = + inputs.E_s.xy[k][(inputs.J1.index)][i] + + inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource + .real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7])); + if (is_conductive) { + J_c.xy[k][(inputs.J1.index)][i] -= + inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7])); + } + if (inputs.params.is_disp_ml) { + J_s.xy[k][(inputs.J1.index)][i] += + inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / + (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + real(commonAmplitude * commonPhase * + (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7])); + } + } + } + } + } +} + +void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { + // Only run update equations is source data was provided + if (inputs.Ksource.is_empty()) { return; } + + // TODO: Line of doom hangover! k is effectively this at all times + int k = inputs.K1.index; + + int array_ind; + int K_tot = n_Yee_cells().k; + + complex commonPhase = + exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + double commonAmplitude = linear_ramp(time_H); + + // Update across K0, provided a source term is there + if (inputs.K0.apply) { + for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { + for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { + if (j < (inputs.J1.index)) { + inputs.E_s.yz[(inputs.K0.index)][j][i] = + inputs.E_s.yz[(inputs.K0.index)][j][i] - + inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource + .real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2])); + if (is_conductive) { + J_c.yz[(inputs.K0.index)][j][i] += + inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2])); + } + if (inputs.params.is_disp_ml) { + J_s.yz[(inputs.K0.index)][j][i] -= + inputs.matched_layer.kappa.z[(inputs.K0.index)] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2])); + } + } + if (i < (inputs.I1.index)) { + inputs.E_s.xz[(inputs.K0.index)][j][i] = + inputs.E_s.xz[(inputs.K0.index)][j][i] + + inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource + .real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3])); + if (is_conductive) { + J_c.xz[(inputs.K0.index)][j][i] -= + inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3])); + } + if (inputs.params.is_disp_ml) { + J_s.xz[(inputs.K0.index)][j][i] += + inputs.matched_layer.kappa.z[(inputs.K0.index)] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.z[inputs.K0.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3])); + } + } + } + } + } + + // Update across K1, provided a source term is there + if (inputs.K1.apply) { + for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { + for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { + if (j < (inputs.J1.index)) { + inputs.E_s.yz[(inputs.K1.index)][j][i] = + inputs.E_s.yz[(inputs.K1.index)][j][i] + + inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource + .real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6])); + if (is_conductive) { + J_c.yz[(inputs.K1.index)][j][i] -= + inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6])); + } + if (inputs.params.is_disp_ml) { + J_s.yz[(inputs.K1.index)][j][i] += + inputs.matched_layer.kappa.z[(inputs.K1.index)] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6])); + } + } + if (i < (inputs.I1.index)) { + inputs.E_s.xz[(inputs.K1.index)][j][i] = + inputs.E_s.xz[(inputs.K1.index)][j][i] - + inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource + .real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7])); + if (is_conductive) { + J_c.xz[(inputs.K1.index)][j][i] += + inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7])); + } + if (inputs.params.is_disp_ml) { + J_s.xz[(inputs.K1.index)][j][i] -= + inputs.matched_layer.kappa.z[(inputs.K1.index)] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.z[inputs.K1.index] * + real(commonAmplitude * commonPhase * + (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7])); + } + } + } + } + } +} From dc18a0459e744b18f866529b3ad04f432ce28f80 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 10:01:27 +0000 Subject: [PATCH 07/25] Pre-commit consistency with #223 --- .pre-commit-config.yaml | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c78c25b90..07ab04c7d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -26,13 +26,11 @@ repos: rev: v1.3.5 hooks: - id: clang-format - args: ["--style=file", "-i"] - - id: clang-tidy - args: ["-p=tdms/build/compile_commands.json"] - # - id: cppcheck - # args: ["--project=tdms/build/compile_commands.json"] - # - id: cpplint - # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] - # suppress whitespace because it conflicts with clang-tidy and doxygen - # suppress copyright in file headers because we don't require them in the developer docs - # run quietly so we don't see you doing this if everything is fine + args: ["--style=file"] + # - id: cppcheck + # args: ["--project=tdms/build/compile_commands.json"] + # - id: cpplint + # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] + # suppress whitespace because it conflicts with clang-tidy and doxygen + # suppress copyright in file headers because we don't require them in the developer docs + # run quietly so we don't see you doing this if everything is fine From 738628899995f6945440692d0d728ede92295a63 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 10:15:52 +0000 Subject: [PATCH 08/25] Update clang-format to #222 settings --- .clang-format | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.clang-format b/.clang-format index 6ffcd105b..1d87759b0 100644 --- a/.clang-format +++ b/.clang-format @@ -33,7 +33,7 @@ BreakBeforeBinaryOperators: None BreakBeforeTernaryOperators: true BreakConstructorInitializers: BeforeColon BreakInheritanceList: BeforeColon -ColumnLimit: 100 +ColumnLimit: 80 CompactNamespaces: false ContinuationIndentWidth: 8 IndentCaseLabels: true @@ -41,11 +41,11 @@ IndentPPDirectives: None IndentWidth: 2 KeepEmptyLinesAtTheStartOfBlocks: true MaxEmptyLinesToKeep: 2 -NamespaceIndentation: All +NamespaceIndentation: None ObjCSpaceAfterProperty: false ObjCSpaceBeforeProtocolList: true PointerAlignment: Right -ReflowComments: false +ReflowComments: true SpaceAfterCStyleCast: true SpaceAfterLogicalNot: false SpaceAfterTemplateKeyword: false From a3c015a1e4d225d5e2ea0a9f2aed8485e5919353 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 10:16:12 +0000 Subject: [PATCH 09/25] Fix out-of-order initilisation warning --- .../simulation_manager/simulation_manager.h | 64 ++++++++++++------- 1 file changed, 41 insertions(+), 23 deletions(-) diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index f98b33126..8b691c266 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -1,6 +1,8 @@ /** * @file simulation_manager.h - * @brief Class that runs the physics of the TDMS executable, and produces an output. Implementation is split between simulation_manager.cpp and execute_simulation.cpp. + * @brief Class that runs the physics of the TDMS executable, and produces an + * output. Implementation is split between simulation_manager.cpp and + * execute_simulation.cpp. */ #pragma once @@ -20,25 +22,35 @@ /** * @brief Manages the physics of TDMS and the simulation loop itself. * - * Takes in the contents of the input file and produces the output in the outputs member, so that it can be written to the specified output file. + * Takes in the contents of the input file and produces the output in the + * outputs member, so that it can be written to the specified output file. * - * Attributes are private to avoid external alterations to their content and potential reassignment of MATLAB pointers, resulting in leaking memory. The outputs can be accessed through a fetch-method within main.cpp if necessary. + * Attributes are private to avoid external alterations to their content and + * potential reassignment of MATLAB pointers, resulting in leaking memory. The + * outputs can be accessed through a fetch-method within main.cpp if necessary. */ 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 + /*! The input objects that are generated from an input file */ + ObjectsFromInfile inputs; - PreferredInterpolationMethods pim;//< The interpolation methods to use in this simulation - SolverMethod solver_method; //< The solver method to use in this simulation + LoopTimers timers; //!< Timers for tracking the execution of the simulation + PSTDVariables PSTD; //!< PSTD-solver-specific variables + FDTDBootstrapper FDTD;//!< FDTD bootstrapping variables - EHVec eh_vec;//< TODO + /*! Output object that will contain the results of this simulation, given the + * input file */ + OutputMatrices outputs; - double ramp_width = 4.;//< Width of the ramp when introducing the waveform in steady state mode + /*! 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 + + /*! Width of the ramp when introducing the waveform in steady state mode */ + double ramp_width = 4.; /** * @brief Evaluates the linear ramp function r(t) * @@ -57,21 +69,23 @@ class SimulationManager { return std::min(1., t / (ramp_width * period)); } - std::vector> - E_norm;//< Holds the E-field phasors norm at each extraction frequency - std::vector> - 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> E_norm; + /*! @copydoc E_norm */ + std::vector> H_norm; /** * @brief Extracts the phasors norms at the given frequency (index) * - * @param frequency_index The index of the frequency (in inputs.f_ex_vec) to extract at + * @param frequency_index The index of the frequency (in inputs.f_ex_vec) to + * extract at * @param tind The current timestep * @param Nt The number of timesteps in a sinusoidal period */ void extract_phasor_norms(int frequency_index, int tind, int Nt); /** - * @brief Creates MATLAB memory blocks that will be iteratively updated in the execute() method, and assigns array sizes for writing purposes later. + * @brief Creates MATLAB memory blocks that will be iteratively updated in the + * execute() method, and assigns array sizes for writing purposes later. * * @param fieldsample The fieldsample input from the input file * @param campssample The complex amplitude sample input from the input file @@ -99,19 +113,23 @@ class SimulationManager { void execute(); /** - * @brief Perform additional (mostly conditional) computations on the outputs, depending on the simulation inputs and run specifications. + * @brief Perform additional (mostly conditional) computations on the outputs, + * depending on the simulation inputs and run specifications. * * This should only be run AFTER a successful run of the execute() method. */ void post_loop_processing(); /** - * @brief Write the outputs to the file provided. Wrapper for outputs.save_outputs + * @brief Write the outputs to the file provided. Wrapper for + * outputs.save_outputs * * @param output_file The filename to write the outputs to - * @param compressed_format If true, write compressed output (do not write facets and vertices) + * @param compressed_format If true, write compressed output (do not write + * facets and vertices) */ - void write_outputs_to_file(std::string output_file, bool compressed_format = false) { + void write_outputs_to_file(std::string output_file, + bool compressed_format = false) { outputs.save_outputs(output_file, compressed_format); } }; From c20f1dbaacec5827e7fbdbc824119be2e54975cb Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 10:21:28 +0000 Subject: [PATCH 10/25] Revert pre-commit to suppress GitHub failures until #223 is ready --- .pre-commit-config.yaml | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 07ab04c7d..82ff6e4f3 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -21,16 +21,3 @@ repos: hooks: - id: trailing-whitespace - id: end-of-file-fixer - # Add C++ linting - - repo: https://github.com/pocc/pre-commit-hooks - rev: v1.3.5 - hooks: - - id: clang-format - args: ["--style=file"] - # - id: cppcheck - # args: ["--project=tdms/build/compile_commands.json"] - # - id: cpplint - # args: ["--filter=-whitespace/comments,-legal/copyright", "--linelength=120", "--quiet"] - # suppress whitespace because it conflicts with clang-tidy and doxygen - # suppress copyright in file headers because we don't require them in the developer docs - # run quietly so we don't see you doing this if everything is fine From 04a4353d884faa86fa8aa834fb2779240c19a285 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 10:42:42 +0000 Subject: [PATCH 11/25] Remove unused variables in updating Ksource terms --- .../source_term_updates.cpp | 445 +++++++++++------- 1 file changed, 283 insertions(+), 162 deletions(-) diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_term_updates.cpp index 8c2432d21..1df6a0ad6 100644 --- a/tdms/src/simulation_manager/source_term_updates.cpp +++ b/tdms/src/simulation_manager/source_term_updates.cpp @@ -3,17 +3,17 @@ using namespace std; using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; -void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s) { +void SimulationManager::update_Isource_terms_steadystate( + double time_H, bool is_conductive, CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided if (inputs.Isource.is_empty()) { return; } int array_ind; int I_tot = n_Yee_cells().i; - complex commonPhase = - exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + complex commonPhase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); double commonAmplitude = linear_ramp(time_H); // Update across I0, provided a source term is here @@ -26,31 +26,42 @@ void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_ array_ind = (I_tot + 1) * k + inputs.I0.index; } - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { inputs.E_s.zx[k][j][inputs.I0.index] = inputs.E_s.zx[k][j][inputs.I0.index] - inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)] + [2])); if (is_conductive) { J_c.zx[k][j][inputs.I0.index] += inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2])); } if (inputs.params.is_disp_ml) { J_s.zx[k][j][inputs.I0.index] += - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + inputs.matched_layer.kappa.x[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][2])); } } @@ -59,26 +70,36 @@ void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_ inputs.E_s.yx[k][j][inputs.I0.index] + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)] + [3])); if (is_conductive) { J_c.yx[k][j][inputs.I0.index] -= inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3])); } if (inputs.params.is_disp_ml) { J_s.yx[k][j][inputs.I0.index] -= - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + inputs.matched_layer.kappa.x[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][3])); } } } @@ -95,31 +116,42 @@ void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_ array_ind = (I_tot + 1) * k + inputs.I1.index; } - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { inputs.E_s.zx[k][j][inputs.I1.index] = inputs.E_s.zx[k][j][inputs.I1.index] + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)] + [6])); if (is_conductive) { J_c.zx[k][j][inputs.I1.index] -= inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6])); } if (inputs.params.is_disp_ml) { J_s.zx[k][j][inputs.I1.index] -= - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + inputs.matched_layer.kappa.x[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][6])); } } @@ -128,26 +160,36 @@ void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_ inputs.E_s.yx[k][j][inputs.I1.index] - inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)] + [7])); if (is_conductive) { J_c.yx[k][j][inputs.I1.index] += inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7])); } if (inputs.params.is_disp_ml) { J_s.yx[k][j][inputs.I1.index] += - inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.x[array_ind] * + inputs.matched_layer.kappa.x[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.x[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)][j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); + (inputs.Isource.real[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7] + + IMAGINARY_UNIT * + inputs.Isource + .imag[k - (inputs.K0.index)] + [j - (inputs.J0.index)][7])); } } } @@ -155,24 +197,25 @@ void SimulationManager::update_Isource_terms_steadystate(double time_H, bool is_ } } -void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s) { +void SimulationManager::update_Jsource_terms_steadystate( + double time_H, bool is_conductive, CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided if (inputs.Jsource.is_empty()) { return; } int array_ind; int J_tot = n_Yee_cells().j; - complex commonPhase = - exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + complex commonPhase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); double commonAmplitude = linear_ramp(time_H); // Update across J0, provided a source term is there if (inputs.J0.apply) { for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { for (int i = inputs.I0.index; i < inputs.I1.index; i++) { - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { if (!inputs.params.is_multilayer) { array_ind = inputs.J0.index; } else { @@ -183,26 +226,36 @@ void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_ inputs.E_s.zy[k][(inputs.J0.index)][i] + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)] + [2])); if (is_conductive) { J_c.zy[k][(inputs.J0.index)][i] -= inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2])); } if (inputs.params.is_disp_ml) { J_s.zy[k][(inputs.J0.index)][i] -= - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + inputs.matched_layer.kappa.y[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][2])); } } if (i < (inputs.I1.index)) { @@ -210,26 +263,36 @@ void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_ inputs.E_s.xy[k][(inputs.J0.index)][i] - inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)] + [3])); if (is_conductive) { J_c.xy[k][(inputs.J0.index)][i] += inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3])); } if (inputs.params.is_disp_ml) { J_s.xy[k][(inputs.J0.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + inputs.matched_layer.kappa.y[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][3])); } } } @@ -246,31 +309,42 @@ void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_ array_ind = (J_tot + 1) * k + inputs.J1.index; } - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { inputs.E_s.zy[k][(inputs.J1.index)][i] = inputs.E_s.zy[k][(inputs.J1.index)][i] - inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)] + [6])); if (is_conductive) { J_c.zy[k][(inputs.J1.index)][i] += inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6])); } if (inputs.params.is_disp_ml) { J_s.zy[k][(inputs.J1.index)][i] -= - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + inputs.matched_layer.kappa.y[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][6])); } } if (i < (inputs.I1.index)) { @@ -278,26 +352,36 @@ void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_ inputs.E_s.xy[k][(inputs.J1.index)][i] + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)] + [7])); if (is_conductive) { J_c.xy[k][(inputs.J1.index)][i] -= inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7])); } if (inputs.params.is_disp_ml) { J_s.xy[k][(inputs.J1.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / - (2. * inputs.params.dt) * inputs.C.b.y[array_ind] * + inputs.matched_layer.kappa.y[array_ind] * + inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * + inputs.C.b.y[array_ind] * real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Jsource.real[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Jsource + .imag[k - (inputs.K0.index)] + [i - (inputs.I0.index)][7])); } } } @@ -305,20 +389,17 @@ void SimulationManager::update_Jsource_terms_steadystate(double time_H, bool is_ } } -void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s) { +void SimulationManager::update_Ksource_terms_steadystate( + double time_H, bool is_conductive, CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided if (inputs.Ksource.is_empty()) { return; } // TODO: Line of doom hangover! k is effectively this at all times int k = inputs.K1.index; - int array_ind; - int K_tot = n_Yee_cells().k; - - complex commonPhase = - exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + complex commonPhase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); double commonAmplitude = linear_ramp(time_H); // Update across K0, provided a source term is there @@ -330,17 +411,24 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.E_s.yz[(inputs.K0.index)][j][i] - inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)] + [2])); if (is_conductive) { J_c.yz[(inputs.K0.index)][j][i] += - inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * + inputs.rho_cond.z[(inputs.K0.index)] * + inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2])); } if (inputs.params.is_disp_ml) { J_s.yz[(inputs.K0.index)][j][i] -= @@ -348,9 +436,12 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][2])); } } if (i < (inputs.I1.index)) { @@ -358,17 +449,24 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.E_s.xz[(inputs.K0.index)][j][i] + inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)] + [3])); if (is_conductive) { J_c.xz[(inputs.K0.index)][j][i] -= - inputs.rho_cond.z[(inputs.K0.index)] * inputs.C.b.z[inputs.K0.index] * + inputs.rho_cond.z[(inputs.K0.index)] * + inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3])); } if (inputs.params.is_disp_ml) { J_s.xz[(inputs.K0.index)][j][i] += @@ -376,9 +474,12 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * inputs.C.b.z[inputs.K0.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][3])); } } } @@ -394,17 +495,24 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.E_s.yz[(inputs.K1.index)][j][i] + inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)] + [6])); if (is_conductive) { J_c.yz[(inputs.K1.index)][j][i] -= - inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * + inputs.rho_cond.z[(inputs.K1.index)] * + inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6])); } if (inputs.params.is_disp_ml) { J_s.yz[(inputs.K1.index)][j][i] += @@ -412,9 +520,12 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][6])); } } if (i < (inputs.I1.index)) { @@ -422,17 +533,24 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.E_s.xz[(inputs.K1.index)][j][i] - inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)] + [7])); if (is_conductive) { J_c.xz[(inputs.K1.index)][j][i] += - inputs.rho_cond.z[(inputs.K1.index)] * inputs.C.b.z[inputs.K1.index] * + inputs.rho_cond.z[(inputs.K1.index)] * + inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7])); } if (inputs.params.is_disp_ml) { J_s.xz[(inputs.K1.index)][j][i] -= @@ -440,9 +558,12 @@ void SimulationManager::update_Ksource_terms_steadystate(double time_H, bool is_ inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * inputs.C.b.z[inputs.K1.index] * real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); + (inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][7])); } } } From 92e64e55629987ce9ca9c92cd340b9f24d07a58f Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 12:05:20 +0000 Subject: [PATCH 12/25] Start refactoring the refactored functions. Flag more inconsistencies --- tdms/include/cell_coordinate.h | 14 +- tdms/include/source.h | 45 +- .../source_term_updates.cpp | 432 +++++------------- tdms/src/source.cpp | 13 +- 4 files changed, 172 insertions(+), 332 deletions(-) diff --git a/tdms/include/cell_coordinate.h b/tdms/include/cell_coordinate.h index 1fa1370d8..7a953a47c 100644 --- a/tdms/include/cell_coordinate.h +++ b/tdms/include/cell_coordinate.h @@ -8,9 +8,12 @@ #include /** - * @brief A structure for holding three values, which typically pertain to the same quantity but for each of the axial directions. + * @brief A structure for holding three values, which typically pertain to the + * same quantity but for each of the axial directions. * - * Effectively stores the 3-vector (i,j,k). This is typically used to represent Yee cell indices, or 3D-array dimensions, or the maximum number of Yee cells in each coordinate direction, for example. + * Effectively stores the 3-vector (i,j,k). This is typically used to represent + * Yee cell indices, or 3D-array dimensions, or the maximum number of Yee cells + * in each coordinate direction, for example. */ struct ijk { // The values (by dimension) that the object contains @@ -22,7 +25,12 @@ struct ijk { /** @brief Print the (i,j,k) values */ void print() const { spdlog::info("ijk: ({},{},{})", i, j, k); } - ijk &operator+=(int n) { this->i += n; this->j += n; this->k += n; return *this; } + ijk &operator+=(int n) { + this->i += n; + this->j += n; + this->k += n; + return *this; + } }; /* Synonyms for code readability */ diff --git a/tdms/include/source.h b/tdms/include/source.h index bc440c593..277695f61 100644 --- a/tdms/include/source.h +++ b/tdms/include/source.h @@ -3,13 +3,51 @@ */ #pragma once +#include #include +#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 overlays the ijk struct for code re-use, via the association: + * ijk.k => SourceIndex.cell_c + * ijk.j => SourceIndex.cell_b + * ijk.i => SourceIndex.split_field_ID + */ +typedef ijk SourceIndex; + +/** + * @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]. + * + * source_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: - bool no_data_stored = true;//!< Flags if the array is empty to avoid pointer preservation + /*! Flags if the array is empty to avoid pointer preservation */ + bool no_data_stored = true; + public: double ***real = nullptr;//!< Real data for the source term double ***imag = nullptr;//!< Imag data for the source term @@ -18,4 +56,9 @@ class Source { /** @brief Check if the source term is empty (true) or not (false) */ bool is_empty() { return no_data_stored; } + + std::complex operator[](SourceIndex index) { + return complex(real[index.k][index.j][index.i], + imag[index.k][index.j][index.i]); + } }; diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_term_updates.cpp index 1df6a0ad6..d31cf219f 100644 --- a/tdms/src/simulation_manager/source_term_updates.cpp +++ b/tdms/src/simulation_manager/source_term_updates.cpp @@ -1,5 +1,7 @@ #include "simulation_manager/simulation_manager.h" +#include "cell_coordinate.h" + using namespace std; using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; @@ -26,80 +28,46 @@ void SimulationManager::update_Isource_terms_steadystate( array_ind = (I_tot + 1) * k + inputs.I0.index; } + // Update zx if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zx[k][j][inputs.I0.index] = - inputs.E_s.zx[k][j][inputs.I0.index] - + // Compute common prefactor + SourceIndex s_index{2, j - inputs.J0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)] - [2])); + real(commonAmplitude * commonPhase * inputs.Isource[s_index]); + + inputs.E_s.zx[k][j][inputs.I0.index] -= common_factor; if (is_conductive) { J_c.zx[k][j][inputs.I0.index] += - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); + inputs.rho_cond.x[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.zx[k][j][inputs.I0.index] += inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][2])); + common_factor; } } + // Update yx if (j < (inputs.J1.index)) { - inputs.E_s.yx[k][j][inputs.I0.index] = - inputs.E_s.yx[k][j][inputs.I0.index] + + // Compute common prefactor + SourceIndex s_index{3, j - inputs.J0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)] - [3])); + real(commonAmplitude * commonPhase * inputs.Isource[s_index]); + + inputs.E_s.yx[k][j][inputs.I0.index] += common_factor; if (is_conductive) { J_c.yx[k][j][inputs.I0.index] -= - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); + inputs.rho_cond.x[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.yx[k][j][inputs.I0.index] -= inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][3])); + common_factor; } } } @@ -116,80 +84,45 @@ void SimulationManager::update_Isource_terms_steadystate( array_ind = (I_tot + 1) * k + inputs.I1.index; } + // Update zx if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zx[k][j][inputs.I1.index] = - inputs.E_s.zx[k][j][inputs.I1.index] + + // Common prefactors + SourceIndex s_index{6, j - inputs.J0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)] - [6])); + real(commonAmplitude * commonPhase * inputs.Isource[s_index]); + + inputs.E_s.zx[k][j][inputs.I1.index] += common_factor; if (is_conductive) { J_c.zx[k][j][inputs.I1.index] -= - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); + inputs.rho_cond.x[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.zx[k][j][inputs.I1.index] -= inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][6])); + common_factor; } } if (j < (inputs.J1.index)) { - inputs.E_s.yx[k][j][inputs.I1.index] = - inputs.E_s.yx[k][j][inputs.I1.index] - + // Common prefactors + SourceIndex s_index{7, j - inputs.J0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)] - [7])); + real(commonAmplitude * commonPhase * inputs.Isource[s_index]); + + inputs.E_s.yx[k][j][inputs.I1.index] -= common_factor; if (is_conductive) { J_c.yx[k][j][inputs.I1.index] += - inputs.rho_cond.x[array_ind] * inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); + inputs.rho_cond.x[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.yx[k][j][inputs.I1.index] += inputs.matched_layer.kappa.x[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource.real[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7] + - IMAGINARY_UNIT * - inputs.Isource - .imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][7])); + common_factor; } } } @@ -216,83 +149,48 @@ void SimulationManager::update_Jsource_terms_steadystate( for (int i = inputs.I0.index; i < inputs.I1.index; i++) { if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + // shouldn't this be one-level up to be consistent with the other + // source update steps?!?!?!?! if (!inputs.params.is_multilayer) { array_ind = inputs.J0.index; } else { array_ind = (J_tot + 1) * k + inputs.J0.index; } - inputs.E_s.zy[k][(inputs.J0.index)][i] = - inputs.E_s.zy[k][(inputs.J0.index)][i] + + SourceIndex s_index{2, i - inputs.I0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)] - [2])); + real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); + + inputs.E_s.zy[k][inputs.J0.index][i] += common_factor; if (is_conductive) { - J_c.zy[k][(inputs.J0.index)][i] -= - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); + J_c.zy[k][inputs.J0.index][i] -= + inputs.rho_cond.y[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { - J_s.zy[k][(inputs.J0.index)][i] -= + J_s.zy[k][inputs.J0.index][i] -= inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][2])); + common_factor; } } - if (i < (inputs.I1.index)) { - inputs.E_s.xy[k][(inputs.J0.index)][i] = - inputs.E_s.xy[k][(inputs.J0.index)][i] - + + if (i < inputs.I1.index) { + SourceIndex s_index{3, i - inputs.I0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)] - [3])); + real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); + + inputs.E_s.xy[k][inputs.J0.index][i] -= common_factor; if (is_conductive) { - J_c.xy[k][(inputs.J0.index)][i] += - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); + J_c.xy[k][inputs.J0.index][i] += + inputs.rho_cond.y[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.xy[k][(inputs.J0.index)][i] += inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][3])); + common_factor; } } } @@ -301,8 +199,8 @@ void SimulationManager::update_Jsource_terms_steadystate( // Update across J1, provided a source term is there if (inputs.J1.apply) { - for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { - for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (!inputs.params.is_multilayer) { array_ind = inputs.J1.index; } else { @@ -311,77 +209,40 @@ void SimulationManager::update_Jsource_terms_steadystate( if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - inputs.E_s.zy[k][(inputs.J1.index)][i] = - inputs.E_s.zy[k][(inputs.J1.index)][i] - + SourceIndex s_index{6, i - inputs.I0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)] - [6])); + real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); + + inputs.E_s.zy[k][(inputs.J1.index)][i] -= common_factor; if (is_conductive) { J_c.zy[k][(inputs.J1.index)][i] += - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); + inputs.rho_cond.y[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.zy[k][(inputs.J1.index)][i] -= inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][6])); + common_factor; } } + if (i < (inputs.I1.index)) { - inputs.E_s.xy[k][(inputs.J1.index)][i] = - inputs.E_s.xy[k][(inputs.J1.index)][i] + + SourceIndex s_index{7, i - inputs.I0.index, k - inputs.K0.index}; + double common_factor = inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)] - [7])); + real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); + + inputs.E_s.xy[k][(inputs.J1.index)][i] += common_factor; if (is_conductive) { J_c.xy[k][(inputs.J1.index)][i] -= - inputs.rho_cond.y[array_ind] * inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); + inputs.rho_cond.y[array_ind] * common_factor; } if (inputs.params.is_disp_ml) { J_s.xy[k][(inputs.J1.index)][i] += inputs.matched_layer.kappa.y[array_ind] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource.real[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Jsource - .imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][7])); + common_factor; } } } @@ -395,7 +256,8 @@ void SimulationManager::update_Ksource_terms_steadystate( // Only run update equations is source data was provided if (inputs.Ksource.is_empty()) { return; } - // TODO: Line of doom hangover! k is effectively this at all times + // TODO: Line of doom hangover! k is effectively this at all times, due to + // previous scoping errors int k = inputs.K1.index; complex commonPhase = exp( @@ -407,79 +269,40 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - inputs.E_s.yz[(inputs.K0.index)][j][i] = - inputs.E_s.yz[(inputs.K0.index)][j][i] - + SourceIndex s_index{2, i - inputs.I0.index, j - inputs.J0.index}; + double common_factor = inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)] - [2])); + real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); + + inputs.E_s.yz[(inputs.K0.index)][j][i] -= common_factor; if (is_conductive) { J_c.yz[(inputs.K0.index)][j][i] += - inputs.rho_cond.z[(inputs.K0.index)] * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); + inputs.rho_cond.z[(inputs.K0.index)] * common_factor; } if (inputs.params.is_disp_ml) { J_s.yz[(inputs.K0.index)][j][i] -= inputs.matched_layer.kappa.z[(inputs.K0.index)] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2])); + common_factor; } } + if (i < (inputs.I1.index)) { - inputs.E_s.xz[(inputs.K0.index)][j][i] = - inputs.E_s.xz[(inputs.K0.index)][j][i] + + SourceIndex s_index{3, i - inputs.I0.index, j - inputs.J0.index}; + double common_factor = inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)] - [3])); + real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); + + inputs.E_s.xz[(inputs.K0.index)][j][i] += common_factor; if (is_conductive) { J_c.xz[(inputs.K0.index)][j][i] -= - inputs.rho_cond.z[(inputs.K0.index)] * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); + inputs.rho_cond.z[(inputs.K0.index)] * common_factor; } if (inputs.params.is_disp_ml) { J_s.xz[(inputs.K0.index)][j][i] += inputs.matched_layer.kappa.z[(inputs.K0.index)] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3])); + common_factor; } } } @@ -491,79 +314,40 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - inputs.E_s.yz[(inputs.K1.index)][j][i] = - inputs.E_s.yz[(inputs.K1.index)][j][i] + + SourceIndex s_index{6, i - inputs.I0.index, j - inputs.J0.index}; + double common_factor = inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)] - [6])); + real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); + + inputs.E_s.yz[(inputs.K1.index)][j][i] += common_factor; if (is_conductive) { - J_c.yz[(inputs.K1.index)][j][i] -= - inputs.rho_cond.z[(inputs.K1.index)] * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); + J_c.yz[inputs.K1.index][j][i] -= + inputs.rho_cond.z[(inputs.K1.index)] * common_factor; } if (inputs.params.is_disp_ml) { J_s.yz[(inputs.K1.index)][j][i] += inputs.matched_layer.kappa.z[(inputs.K1.index)] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][6])); + common_factor; } } + if (i < (inputs.I1.index)) { - inputs.E_s.xz[(inputs.K1.index)][j][i] = - inputs.E_s.xz[(inputs.K1.index)][j][i] - + SourceIndex s_index{7, i - inputs.I0.index, j - inputs.J0.index}; + double common_factor = inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)] - [7])); + real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); + + inputs.E_s.xz[(inputs.K1.index)][j][i] -= common_factor; if (is_conductive) { J_c.xz[(inputs.K1.index)][j][i] += - inputs.rho_cond.z[(inputs.K1.index)] * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); + inputs.rho_cond.z[(inputs.K1.index)] * common_factor; } if (inputs.params.is_disp_ml) { J_s.xz[(inputs.K1.index)][j][i] -= inputs.matched_layer.kappa.z[(inputs.K1.index)] * inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * - (inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][7])); + common_factor; } } } diff --git a/tdms/src/source.cpp b/tdms/src/source.cpp index b20881779..a49d44807 100644 --- a/tdms/src/source.cpp +++ b/tdms/src/source.cpp @@ -10,7 +10,8 @@ using namespace std; -Source::Source(const mxArray *ptr, int dim1, int dim2, const std::string &name) { +Source::Source(const mxArray *ptr, int dim1, int dim2, + const std::string &name) { if (mxIsEmpty(ptr)) { spdlog::info("{} is empty", name); @@ -18,7 +19,9 @@ Source::Source(const mxArray *ptr, int dim1, int dim2, const std::string &name) // fetch dimensions of the input array, and check they are what we expect auto dims = Dimensions(ptr); - if (dims.are_1d()) { throw runtime_error(name + " should be 3- or 2-dimensional"); } + if (dims.are_1d()) { + throw runtime_error(name + " should be 3- or 2-dimensional"); + } if (dims.are_2d()) { dim2 = 0; } if (!(dims[0] == 8 && dims[1] == dim1 && dims[2] == dim2)) { cerr << name << " has incorrect size" << endl; @@ -26,8 +29,10 @@ Source::Source(const mxArray *ptr, int dim1, int dim2, const std::string &name) // check that complex data has been passed if (!mxIsComplex(ptr)) { - throw runtime_error(name + " should be complex, use a call of " - "complex(real(Isource),imag(Isource)) in matlab if necessary"); + throw runtime_error( + name + + " should be complex, use a call of " + "complex(real(Isource),imag(Isource)) in matlab if necessary"); } // cast MATLAB arrays to C++ datatypes From 10a41b17c6e993f4265c6929366d4f54261f0edb Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 12:06:50 +0000 Subject: [PATCH 13/25] Missing std:: --- tdms/include/source.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tdms/include/source.h b/tdms/include/source.h index 277695f61..fc695e2bc 100644 --- a/tdms/include/source.h +++ b/tdms/include/source.h @@ -58,7 +58,7 @@ class Source { bool is_empty() { return no_data_stored; } std::complex operator[](SourceIndex index) { - return complex(real[index.k][index.j][index.i], - imag[index.k][index.j][index.i]); + return std::complex(real[index.k][index.j][index.i], + imag[index.k][index.j][index.i]); } }; From d3e838de889c2f71a8797b43b5623972e8fccefd Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 16:51:46 +0000 Subject: [PATCH 14/25] Merge E_s steady-state updates into one function (revert files to previous commit if undesired) --- tdms/include/field.h | 2 +- .../simulation_manager/simulation_manager.h | 40 ++ .../source_term_updates.cpp | 411 ++++++++---------- 3 files changed, 222 insertions(+), 231 deletions(-) diff --git a/tdms/include/field.h b/tdms/include/field.h index 29a426e89..d1452a322 100644 --- a/tdms/include/field.h +++ b/tdms/include/field.h @@ -55,7 +55,7 @@ class SplitFieldComponent: public Tensor3D{ 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); diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index 8b691c266..b1500bca7 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -92,12 +92,52 @@ class SimulationManager { */ void prepare_output(const mxArray *fieldsample, const mxArray *campssample); + /** + * @brief Update electric-split field components and current densities 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 update_source_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 Performs updates to the electric-split field components and current + * density fields after an E-field timestep has been performed, in accordance + * with the {IJK} 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_Isource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); + /*! @copydoc update_Isource_terms_steadystate */ void update_Jsource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); + /*! @copydoc update_Isource_terms_steadystate */ void update_Ksource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_term_updates.cpp index d31cf219f..657d72fb5 100644 --- a/tdms/src/simulation_manager/source_term_updates.cpp +++ b/tdms/src/simulation_manager/source_term_updates.cpp @@ -12,11 +12,6 @@ void SimulationManager::update_Isource_terms_steadystate( if (inputs.Isource.is_empty()) { return; } int array_ind; - int I_tot = n_Yee_cells().i; - - complex commonPhase = exp( - -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); - double commonAmplitude = linear_ramp(time_H); // Update across I0, provided a source term is here if (inputs.I0.apply) { @@ -25,50 +20,16 @@ void SimulationManager::update_Isource_terms_steadystate( if (!inputs.params.is_multilayer) { array_ind = inputs.I0.index; } else { - array_ind = (I_tot + 1) * k + inputs.I0.index; + array_ind = (n_Yee_cells().i + 1) * k + inputs.I0.index; } - - // Update zx if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - // Compute common prefactor - SourceIndex s_index{2, j - inputs.J0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * inputs.Isource[s_index]); - - inputs.E_s.zx[k][j][inputs.I0.index] -= common_factor; - if (is_conductive) { - J_c.zx[k][j][inputs.I0.index] += - inputs.rho_cond.x[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.zx[k][j][inputs.I0.index] += - inputs.matched_layer.kappa.x[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::X, true, true, + is_conductive, j, k, array_ind, J_c, J_s); } - - // Update yx if (j < (inputs.J1.index)) { - // Compute common prefactor - SourceIndex s_index{3, j - inputs.J0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * inputs.Isource[s_index]); - - inputs.E_s.yx[k][j][inputs.I0.index] += common_factor; - if (is_conductive) { - J_c.yx[k][j][inputs.I0.index] -= - inputs.rho_cond.x[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.yx[k][j][inputs.I0.index] -= - inputs.matched_layer.kappa.x[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::X, false, true, + is_conductive, j, k, array_ind, J_c, J_s); } } } @@ -81,49 +42,16 @@ void SimulationManager::update_Isource_terms_steadystate( if (!inputs.params.is_multilayer) { array_ind = inputs.I1.index; } else { - array_ind = (I_tot + 1) * k + inputs.I1.index; + array_ind = (n_Yee_cells().i + 1) * k + inputs.I1.index; } - - // Update zx if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - // Common prefactors - SourceIndex s_index{6, j - inputs.J0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * inputs.Isource[s_index]); - - inputs.E_s.zx[k][j][inputs.I1.index] += common_factor; - if (is_conductive) { - J_c.zx[k][j][inputs.I1.index] -= - inputs.rho_cond.x[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.zx[k][j][inputs.I1.index] -= - inputs.matched_layer.kappa.x[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::X, true, false, + is_conductive, j, k, array_ind, J_c, J_s); } - if (j < (inputs.J1.index)) { - // Common prefactors - SourceIndex s_index{7, j - inputs.J0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.x[array_ind] * - real(commonAmplitude * commonPhase * inputs.Isource[s_index]); - - inputs.E_s.yx[k][j][inputs.I1.index] -= common_factor; - if (is_conductive) { - J_c.yx[k][j][inputs.I1.index] += - inputs.rho_cond.x[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.yx[k][j][inputs.I1.index] += - inputs.matched_layer.kappa.x[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::X, false, false, + is_conductive, j, k, array_ind, J_c, J_s); } } } @@ -137,11 +65,6 @@ void SimulationManager::update_Jsource_terms_steadystate( if (inputs.Jsource.is_empty()) { return; } int array_ind; - int J_tot = n_Yee_cells().j; - - complex commonPhase = exp( - -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); - double commonAmplitude = linear_ramp(time_H); // Update across J0, provided a source term is there if (inputs.J0.apply) { @@ -154,44 +77,15 @@ void SimulationManager::update_Jsource_terms_steadystate( if (!inputs.params.is_multilayer) { array_ind = inputs.J0.index; } else { - array_ind = (J_tot + 1) * k + inputs.J0.index; + array_ind = (n_Yee_cells().j + 1) * k + inputs.J0.index; } - SourceIndex s_index{2, i - inputs.I0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); - - inputs.E_s.zy[k][inputs.J0.index][i] += common_factor; - if (is_conductive) { - J_c.zy[k][inputs.J0.index][i] -= - inputs.rho_cond.y[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.zy[k][inputs.J0.index][i] -= - inputs.matched_layer.kappa.y[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, true, true, + is_conductive, i, k, array_ind, J_c, J_s); } - if (i < inputs.I1.index) { - SourceIndex s_index{3, i - inputs.I0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); - - inputs.E_s.xy[k][inputs.J0.index][i] -= common_factor; - if (is_conductive) { - J_c.xy[k][inputs.J0.index][i] += - inputs.rho_cond.y[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.xy[k][(inputs.J0.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, false, true, + is_conductive, i, k, array_ind, J_c, J_s); } } } @@ -204,46 +98,16 @@ void SimulationManager::update_Jsource_terms_steadystate( if (!inputs.params.is_multilayer) { array_ind = inputs.J1.index; } else { - array_ind = (J_tot + 1) * k + inputs.J1.index; + array_ind = (n_Yee_cells().j + 1) * k + inputs.J1.index; } - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - SourceIndex s_index{6, i - inputs.I0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); - - inputs.E_s.zy[k][(inputs.J1.index)][i] -= common_factor; - if (is_conductive) { - J_c.zy[k][(inputs.J1.index)][i] += - inputs.rho_cond.y[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.zy[k][(inputs.J1.index)][i] -= - inputs.matched_layer.kappa.y[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, true, false, + is_conductive, i, k, array_ind, J_c, J_s); } - if (i < (inputs.I1.index)) { - SourceIndex s_index{7, i - inputs.I0.index, k - inputs.K0.index}; - double common_factor = - inputs.C.b.y[array_ind] * - real(commonAmplitude * commonPhase * inputs.Jsource[s_index]); - - inputs.E_s.xy[k][(inputs.J1.index)][i] += common_factor; - if (is_conductive) { - J_c.xy[k][(inputs.J1.index)][i] -= - inputs.rho_cond.y[array_ind] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.xy[k][(inputs.J1.index)][i] += - inputs.matched_layer.kappa.y[array_ind] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, false, false, + is_conductive, i, k, array_ind, J_c, J_s); } } } @@ -256,54 +120,19 @@ void SimulationManager::update_Ksource_terms_steadystate( // Only run update equations is source data was provided if (inputs.Ksource.is_empty()) { return; } - // TODO: Line of doom hangover! k is effectively this at all times, due to - // previous scoping errors - int k = inputs.K1.index; - - complex commonPhase = exp( - -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); - double commonAmplitude = linear_ramp(time_H); - // Update across K0, provided a source term is there if (inputs.K0.apply) { for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - SourceIndex s_index{2, i - inputs.I0.index, j - inputs.J0.index}; - double common_factor = - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); - - inputs.E_s.yz[(inputs.K0.index)][j][i] -= common_factor; - if (is_conductive) { - J_c.yz[(inputs.K0.index)][j][i] += - inputs.rho_cond.z[(inputs.K0.index)] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.yz[(inputs.K0.index)][j][i] -= - inputs.matched_layer.kappa.z[(inputs.K0.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, true, true, + is_conductive, i, j, inputs.K0.index, J_c, + J_s); } - if (i < (inputs.I1.index)) { - SourceIndex s_index{3, i - inputs.I0.index, j - inputs.J0.index}; - double common_factor = - inputs.C.b.z[inputs.K0.index] * - real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); - - inputs.E_s.xz[(inputs.K0.index)][j][i] += common_factor; - if (is_conductive) { - J_c.xz[(inputs.K0.index)][j][i] -= - inputs.rho_cond.z[(inputs.K0.index)] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.xz[(inputs.K0.index)][j][i] += - inputs.matched_layer.kappa.z[(inputs.K0.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Y, false, true, + is_conductive, i, j, inputs.K0.index, J_c, + J_s); } } } @@ -314,43 +143,165 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - SourceIndex s_index{6, i - inputs.I0.index, j - inputs.J0.index}; - double common_factor = - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); - - inputs.E_s.yz[(inputs.K1.index)][j][i] += common_factor; - if (is_conductive) { - J_c.yz[inputs.K1.index][j][i] -= - inputs.rho_cond.z[(inputs.K1.index)] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.yz[(inputs.K1.index)][j][i] += - inputs.matched_layer.kappa.z[(inputs.K1.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Z, true, false, + is_conductive, i, j, inputs.K1.index, J_c, + J_s); } - if (i < (inputs.I1.index)) { - SourceIndex s_index{7, i - inputs.I0.index, j - inputs.J0.index}; - double common_factor = - inputs.C.b.z[inputs.K1.index] * - real(commonAmplitude * commonPhase * inputs.Ksource[s_index]); - - inputs.E_s.xz[(inputs.K1.index)][j][i] -= common_factor; - if (is_conductive) { - J_c.xz[(inputs.K1.index)][j][i] += - inputs.rho_cond.z[(inputs.K1.index)] * common_factor; - } - if (inputs.params.is_disp_ml) { - J_s.xz[(inputs.K1.index)][j][i] -= - inputs.matched_layer.kappa.z[(inputs.K1.index)] * - inputs.matched_layer.gamma[k] / (2. * inputs.params.dt) * - common_factor; - } + update_source_steadystate(time_H, AxialDirection::Z, false, false, + is_conductive, i, j, inputs.K1.index, J_c, + J_s); } } } } } + +void SimulationManager::update_source_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) { + // Common phase and amplitude terms to apply in all computations + complex common_phase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + double common_amplitude = linear_ramp(time_H); + + double c_constant; //< C.b.{x,y,z} constant that appears in the update + // equations + double conductive_aux;//< Conductivity, only present in a conductive material. + //= inputs.rho_cond.{x,y,z} + + /*! Dispersive parameter. = inputs.matched_layer.kappa.{x,y,z} */ + double kappa; + /*! NOTE: if dealing with Ksource, then the index used to fetch the gamma + * parameter is k = inputs.K1.index, rather than cell_c. This inconsistency & + * potential error is flagged in #221. Have a feeling that, in the K-case, + * this should be inputs.K{0,1}.index depending on the plane we're in */ + int gamma_ind = (parallel == AxialDirection::Z) ? inputs.K1.index : cell_c; + /*! Dispersive parameter. = inputs.matched_layer.gamma[k]. */ + double gamma = inputs.matched_layer.gamma[gamma_ind]; + + /*! Value of the source term in Yee cell (cell_a, cell_b, cell_c) */ + complex source_value; + /*! The split-field cell index to update */ + CellCoordinate cell_to_update; + /*! Field updates are signed based on the component, plane, and axis. [0] = + * E_s sign, [1] = J_c sign, [2] = J_s sign */ + double update_sign[3] = {1., 1., 1.}; + + SplitFieldComponent + *E_s_component,//< The E_s component that should be updated + *J_c_component,//< The J_c component that should be updated + *J_s_component;//< The J_s component that should be updated + + // If we are currently updating along the C-axis, the split_field_ID is one + // less than if we are updating along the B-axis + int split_field_ID = (C_axis) ? 2 : 3; + // If we are updating in the {IJK}1, rather than the {IJK}0 plane, the + // split_field_ID is 4 greater than usual (uses indices 6/7 as opposed to 2/3) + if (!zero_plane) { split_field_ID += 4; } + /*! This is the index that is to be passed to the (relevent) Source instance, + * to recover the source value at the relevant Yee cell. + + cell_b and cell_c need to be offset by the corresponding {IJK}0.index before + we cna retrieve the source value. */ + SourceIndex s_index = {split_field_ID, cell_b, cell_c}; + + // setup axis-dependent parameters + switch (parallel) { + case AxialDirection::X: + // Dealing with Isource. j = cell_b, k = cell_c + c_constant = inputs.C.b.x[array_ind]; + + s_index.j -= inputs.J0.index; + s_index.k -= inputs.K0.index; + source_value = inputs.Isource[s_index]; + + cell_to_update = {inputs.I0.index, cell_b, cell_c}; + if (!zero_plane) { cell_to_update.i = inputs.I1.index; } + + conductive_aux = inputs.rho_cond.x[array_ind]; + kappa = inputs.matched_layer.kappa.x[array_ind]; + + if (C_axis == zero_plane) { + update_sign[0] = -1.; + } else { + update_sign[1] = -1.; + update_sign[2] = -1.; + } + + E_s_component = (C_axis) ? &inputs.E_s.zx : &inputs.E_s.yx; + J_c_component = (C_axis) ? &J_c.zx : &J_c.yx; + J_s_component = (C_axis) ? &J_s.zx : &J_s.yx; + break; + case AxialDirection::Y: + // Dealing with Jsource. i = cell_b, k = cell_c + c_constant = inputs.C.b.y[array_ind]; + + s_index.j -= inputs.I0.index; + s_index.k -= inputs.K0.index; + source_value = inputs.Jsource[s_index]; + + cell_to_update = {cell_b, inputs.J0.index, cell_c}; + if (!zero_plane) { cell_to_update.j = inputs.J1.index; } + + conductive_aux = inputs.rho_cond.y[array_ind]; + kappa = inputs.matched_layer.kappa.y[array_ind]; + + // again, odd - the others have a symmetric whereas here we have an + // asymmetry + if (C_axis != zero_plane) { + update_sign[0] = -1.; + } else { + update_sign[1] = -1.; + } + if (C_axis) { update_sign[2] = -1.; } + + E_s_component = (C_axis) ? &inputs.E_s.zy : &inputs.E_s.xy; + J_c_component = (C_axis) ? &J_c.zy : &J_c.xy; + J_s_component = (C_axis) ? &J_s.zy : &J_s.xy; + break; + case AxialDirection::Z: + // Dealing with Ksource. i = cell_b, j = cell_c + c_constant = inputs.C.b.z[array_ind]; + + s_index.j -= inputs.I0.index; + s_index.k -= inputs.I1.index; + source_value = inputs.Ksource[s_index]; + + cell_to_update = {cell_b, cell_c, inputs.K0.index}; + if (!zero_plane) { cell_to_update.k = inputs.K1.index; } + + conductive_aux = inputs.rho_cond.z[array_ind]; + kappa = inputs.matched_layer.kappa.z[array_ind]; + + if (C_axis != zero_plane) { + update_sign[1] = -1.; + } else { + update_sign[0] = -1.; + update_sign[2] = -1.; + } + + E_s_component = (C_axis) ? &inputs.E_s.yz : &inputs.E_s.xz; + J_c_component = (C_axis) ? &J_c.yz : &J_c.xz; + J_s_component = (C_axis) ? &J_s.yz : &J_s.xz; + break; + } + + // we can now compute the field update values + double E_split_update = + c_constant * real(common_amplitude * common_phase * source_value); + + // update the relevant split-field component + (*E_s_component)[cell_to_update] += update_sign[0] * E_split_update; + // update the current density in a conductive medium + if (is_conductive) { + (*J_c_component)[cell_to_update] += + update_sign[1] * conductive_aux * E_split_update; + } + // update the current density in a dispersive medium + if (inputs.params.is_disp_ml) { + (*J_s_component)[cell_to_update] += update_sign[2] * kappa * gamma / + (2 * inputs.params.dt) * E_split_update; + } +} From aaf61e6d989dfe6513bc3951bc4e4926df939802 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 2 Feb 2023 17:05:36 +0000 Subject: [PATCH 15/25] Fix bugs in unification (IE Will can't type) --- tdms/src/simulation_manager/source_term_updates.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_term_updates.cpp index 657d72fb5..59b9b382d 100644 --- a/tdms/src/simulation_manager/source_term_updates.cpp +++ b/tdms/src/simulation_manager/source_term_updates.cpp @@ -125,12 +125,12 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - update_source_steadystate(time_H, AxialDirection::Y, true, true, + update_source_steadystate(time_H, AxialDirection::Z, true, true, is_conductive, i, j, inputs.K0.index, J_c, J_s); } if (i < (inputs.I1.index)) { - update_source_steadystate(time_H, AxialDirection::Y, false, true, + update_source_steadystate(time_H, AxialDirection::Z, false, true, is_conductive, i, j, inputs.K0.index, J_c, J_s); } @@ -266,7 +266,7 @@ void SimulationManager::update_source_steadystate( c_constant = inputs.C.b.z[array_ind]; s_index.j -= inputs.I0.index; - s_index.k -= inputs.I1.index; + s_index.k -= inputs.J0.index; source_value = inputs.Ksource[s_index]; cell_to_update = {cell_b, cell_c, inputs.K0.index}; From c146fa40fbf95fc754ce525c00e1b40e987d6230 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 11:17:12 +0000 Subject: [PATCH 16/25] Pull out E_s source updates, and wrap to avoid empty-source updates --- .../simulation_manager/simulation_manager.h | 25 ++- .../simulation_manager/execute_simulation.cpp | 211 +++--------------- .../source_updates_pulsed.cpp | 99 ++++++++ ...tes.cpp => source_updates_steadystate.cpp} | 9 + 4 files changed, 161 insertions(+), 183 deletions(-) create mode 100644 tdms/src/simulation_manager/source_updates_pulsed.cpp rename tdms/src/simulation_manager/{source_term_updates.cpp => source_updates_steadystate.cpp} (96%) diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index b1500bca7..fc7548439 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -122,7 +122,7 @@ class SimulationManager { /** * @brief Performs updates to the electric-split field components and current * density fields after an E-field timestep has been performed, in accordance - * with the {IJK} source terms. + * 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 @@ -130,18 +130,37 @@ class SimulationManager { * @param J_c The current density in the conductive medium * @param J_s The current density in the dispersive medium */ + void update_source_terms_steadystate(double time_H, bool is_conductive, + CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s); + /*! @copydoc update_source_terms_steadystate */ void update_Isource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); - /*! @copydoc update_Isource_terms_steadystate */ + /*! @copydoc update_source_terms_steadystate */ void update_Jsource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); - /*! @copydoc update_Isource_terms_steadystate */ + /*! @copydoc update_source_terms_steadystate */ void update_Ksource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); + /** + * @brief 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); + public: SimulationManager(InputMatrices in_matrices, SolverMethod _solver_method, PreferredInterpolationMethods _pim); diff --git a/tdms/src/simulation_manager/execute_simulation.cpp b/tdms/src/simulation_manager/execute_simulation.cpp index 05a6e6bd9..178eec73c 100644 --- a/tdms/src/simulation_manager/execute_simulation.cpp +++ b/tdms/src/simulation_manager/execute_simulation.cpp @@ -2134,191 +2134,42 @@ void SimulationManager::execute() { if (TIME_EXEC) { timers.click_timer(TimersTrackingLoop::INTERNAL); } /********************/ - //update terms for self consistency across scattered/total interface - E updates## + /* Update source terms for self consistency across scattered/total + * interface. E_s updates use time_H simulation time. */ if (inputs.params.source_mode == SourceMode::steadystate) { // Steady-state source term updates - update_Isource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, - loop_variables.J_s); - update_Jsource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, - loop_variables.J_s); - update_Ksource_terms_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, - loop_variables.J_s); - - complex commonPhase = - exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_H, 2. * DCPI)); - double commonAmplitude = linear_ramp(time_H); - outputs.H.ft = real(commonAmplitude * commonPhase); - } else if (inputs.params.source_mode == SourceMode::pulsed) {//pulsed - - if (J_tot == 0) { - j = 0; - for (i = 0; i < (I_tot + 1); i++) { - inputs.E_s.yz[inputs.K0.index][j][i] = - inputs.E_s.yz[inputs.K0.index][j][i] - - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[0][i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + - inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //E_s.yz[(int)K0[0]][j][i] = E_s.yz[(int)K0[0]][j][i] - C.b.z[(int)K0[0]]*real((Ksource.real[0][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[0][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - if (loop_variables.is_conductive) - loop_variables.J_c.yz[inputs.K0.index][j][i] += - inputs.rho_cond.z[inputs.K0.index] * inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[0][i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_c.yz[(int)K0[0]][j][i] += rho_cond.z[(int)K0[0]]*C.b.z[(int)K0[0]]*real((Ksource.real[0][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[0][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - if (inputs.params.is_disp_ml) { - loop_variables.J_s.yz[inputs.K0.index][j][i] -= - inputs.matched_layer.kappa.z[inputs.K0.index] * - inputs.matched_layer.gamma[inputs.K0.index] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[0][i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_s.yz[(int)K0[0]][j][i] -= matched_layer.kappa.z[(int)K0[0]]*matched_layer.gamma[(int)K0[0]]/(2.*params.dt)*C.b.z[(int)K0[0]]*real((Ksource.real[0][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[0][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - } - } - } else - for (j = 0; j < J_tot; j++) - for (i = 0; i < (I_tot + 1); i++) { - inputs.E_s.yz[inputs.K0.index][j][i] = - inputs.E_s.yz[inputs.K0.index][j][i] - - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + - inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //E_s.yz[(int)K0[0]][j][i] = E_s.yz[(int)K0[0]][j][i] - C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - if (loop_variables.is_conductive) - loop_variables.J_c.yz[inputs.K0.index][j][i] += - inputs.rho_cond.z[inputs.K0.index] * inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + - inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_c.yz[(int)K0[0]][j][i] += rho_cond.z[(int)K0[0]]*C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - if (inputs.params.is_disp_ml) { - loop_variables.J_s.yz[inputs.K0.index][j][i] -= - inputs.matched_layer.kappa.z[inputs.K0.index] * - inputs.matched_layer.gamma[inputs.K0.index] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][2] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][2]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2. * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + - inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_s.yz[(int)K0[0]][j][i] -= matched_layer.kappa.z[(int)K0[0]]*matched_layer.gamma[(int)K0[0]]/(2.*params.dt)*C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][2] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][2])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - } - } - for (j = 0; j < (J_tot + 1); j++) - for (i = 0; i < I_tot; i++) { - inputs.E_s.xz[inputs.K0.index][j][i] = - inputs.E_s.xz[inputs.K0.index][j][i] + - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2 * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + - inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //E_s.xz[(int)K0[0]][j][i] = E_s.xz[(int)K0[0]][j][i] + C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][3] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][3])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2 )); - if (loop_variables.is_conductive) - loop_variables.J_c.xz[inputs.K0.index][j][i] -= - inputs.rho_cond.z[inputs.K0.index] * inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2 * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_c.xz[(int)K0[0]][j][i] -= rho_cond.z[(int)K0[0]]*C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][3] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][3])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2 )); - if (inputs.params.is_disp_ml) - loop_variables.J_s.xz[inputs.K0.index][j][i] += - inputs.matched_layer.kappa.z[inputs.K0.index] * - inputs.matched_layer.gamma[inputs.K0.index] / (2. * inputs.params.dt) * - inputs.C.b.z[inputs.K0.index] * - real((inputs.Ksource.real[j - (inputs.J0.index)][i - (inputs.I0.index)][3] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][3]) * - (-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), - 2 * DCPI))) * - exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), - 2)); - //J_s.xz[(int)K0[0]][j][i] += matched_layer.kappa.z[(int)K0[0]]*matched_layer.gamma[(int)K0[0]]/(2.*params.dt)*C.b.z[(int)K0[0]]*real((Ksource.real[j-((int)J0[0])][i-((int)I0[0])][3] + IMAGINARY_UNIT*Ksource.imag[j-((int)J0[0])][i-((int)I0[0])][3])*(-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2 )); - } - //fth = real((-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); - outputs.H.ft = - real((-1.0 * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), 2. * DCPI))) * + update_source_terms_steadystate(time_H, loop_variables.is_conductive, + loop_variables.J_c, loop_variables.J_s); + + // Common phase term in update equations + complex common_phase = + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * time_H, 2. * DCPI)); + // Common amplitude factor in update equations + double common_amplitude = linear_ramp(time_H); + // Update output H-field + outputs.H.ft = real(common_amplitude * common_phase); + } else if (inputs.params.source_mode == SourceMode::pulsed) { + // Pulsed source term updates + update_source_terms_pulsed(time_H, loop_variables.is_conductive, + loop_variables.J_c, loop_variables.J_s); + + // Common amplitude factor in update equations + double common_amplitude = exp(-1.0 * DCPI * - pow((time_H - inputs.params.to_l + inputs.params.delta.dz / LIGHT_V / 2.) / - (inputs.params.hwhm), + pow((time_H - inputs.params.to_l + + inputs.params.delta.dz / LIGHT_V / 2.) / + inputs.params.hwhm, 2)); - //fth = real((-1.0*IMAGINARY_UNIT)*exp(-IMAGINARY_UNIT*fmod(params.omega_an*(time_H - params.to_l),2.*DCPI)))*exp( -1.0*DCPI*pow((time_H - params.to_l)/(params.hwhm),2)); + // Common phase term in update equations + complex common_phase = + -1.0 * IMAGINARY_UNIT * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), + 2. * DCPI)); + // Update output H-field + outputs.H.ft = real(common_phase) * common_amplitude; } //end of source terms diff --git a/tdms/src/simulation_manager/source_updates_pulsed.cpp b/tdms/src/simulation_manager/source_updates_pulsed.cpp new file mode 100644 index 000000000..89cd4268f --- /dev/null +++ b/tdms/src/simulation_manager/source_updates_pulsed.cpp @@ -0,0 +1,99 @@ +#include "simulation_manager/simulation_manager.h" + +#include "cell_coordinate.h" + +using namespace std; +using namespace tdms_phys_constants; +using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; + +void SimulationManager::update_source_terms_pulsed( + double time_H, bool is_conductive, CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { + //! Exit now if Ksource is empty, to avoid seg-faults + if (inputs.Ksource.is_empty()) { return; } + + //! Simulation dimensions + int I_tot = n_Yee_cells().i, J_tot = n_Yee_cells().j; + //! The material constant that appears in the update equation + double c_constant = inputs.C.b.z[inputs.K0.index]; + //! Conductive aux of the cell being updated + double conductive_aux = inputs.rho_cond.z[inputs.K0.index]; + //! Dispersion prefactor in update equations + double dispersion_factor = inputs.matched_layer.kappa.z[inputs.K0.index] * + inputs.matched_layer.gamma[inputs.K0.index] / + (2. * inputs.params.dt); + //! Common amplitude factor in update equations + double common_amplitude = + c_constant * exp(-1.0 * DCPI * + pow((time_H - inputs.params.to_l + + inputs.params.delta.dz / LIGHT_V / 2.) / + inputs.params.hwhm, + 2)); + //! Common phase term in update equations + complex common_phase = + -1.0 * IMAGINARY_UNIT * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * (time_H - inputs.params.to_l), + 2. * DCPI)); + //! Source term index to read from + SourceIndex s_index; + //! split-field cell index to update + CellCoordinate cell_to_update; + //! The update to apply to the split-field + double split_field_update; + + if (J_tot == 0) { + int j = 0; + for (int i = 0; i < I_tot + 1; i++) { + s_index = {2, i - inputs.I0.index, j}; + cell_to_update = {i, j, inputs.K0.index}; + + split_field_update = + common_amplitude * real(inputs.Ksource[s_index] * common_phase); + + inputs.E_s.yz[cell_to_update] -= split_field_update; + if (is_conductive) { + J_c.yz[cell_to_update] += conductive_aux * split_field_update; + } + if (inputs.params.is_disp_ml) { + J_s.yz[cell_to_update] -= dispersion_factor * split_field_update; + } + } + } else { + for (int j = 0; j < J_tot; j++) { + for (int i = 0; i < I_tot + 1; i++) { + s_index = {2, i - inputs.I0.index, j - inputs.J0.index}; + cell_to_update = {i, j, inputs.K0.index}; + + split_field_update = + common_amplitude * real(inputs.Ksource[s_index] * common_phase); + + inputs.E_s.yz[cell_to_update] -= split_field_update; + if (is_conductive) { + J_c.yz[cell_to_update] += conductive_aux * split_field_update; + } + if (inputs.params.is_disp_ml) { + J_s.yz[cell_to_update] -= dispersion_factor * split_field_update; + } + } + } + } + + for (int j = 0; j < J_tot + 1; j++) { + for (int i = 0; i < I_tot; i++) { + s_index = {3, i - inputs.I0.index, j - inputs.J0.index}; + cell_to_update = {i, j, inputs.K0.index}; + + split_field_update = + common_amplitude * real(inputs.Ksource[s_index] * common_phase); + + inputs.E_s.xz[cell_to_update] += split_field_update; + if (is_conductive) { + J_c.xz[cell_to_update] -= conductive_aux * split_field_update; + } + if (inputs.params.is_disp_ml) { + J_s.xz[cell_to_update] += dispersion_factor * split_field_update; + } + } + } +} diff --git a/tdms/src/simulation_manager/source_term_updates.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp similarity index 96% rename from tdms/src/simulation_manager/source_term_updates.cpp rename to tdms/src/simulation_manager/source_updates_steadystate.cpp index 59b9b382d..0f5cfd857 100644 --- a/tdms/src/simulation_manager/source_term_updates.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -3,8 +3,17 @@ #include "cell_coordinate.h" using namespace std; +using namespace tdms_phys_constants; using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; +void SimulationManager::update_source_terms_steadystate( + double time_H, bool is_conductive, CurrentDensitySplitField &J_c, + CurrentDensitySplitField &J_s) { + update_Isource_terms_steadystate(time_H, is_conductive, J_c, J_s); + update_Jsource_terms_steadystate(time_H, is_conductive, J_c, J_s); + update_Ksource_terms_steadystate(time_H, is_conductive, J_c, J_s); +} + void SimulationManager::update_Isource_terms_steadystate( double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s) { From 962c4547ca0e11ddf8c80f7f73c1e5ef7886a54c Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 14:21:42 +0000 Subject: [PATCH 17/25] Pull out H-source updates, don't add empty array checks yet --- .../simulation_manager/simulation_manager.h | 30 +- .../simulation_manager/execute_simulation.cpp | 286 ++---------------- .../source_updates_pulsed.cpp | 81 +++++ .../source_updates_steadystate.cpp | 204 ++++++++++++- 4 files changed, 337 insertions(+), 264 deletions(-) diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index fc7548439..df966cfb2 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -120,9 +120,9 @@ class SimulationManager { int array_ind, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); /** - * @brief 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. + * @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 @@ -133,6 +133,15 @@ class SimulationManager { void update_source_terms_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 update_source_terms_steadystate(double time_E); + /*! @copydoc update_source_terms_steadystate */ void update_Isource_terms_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, @@ -147,9 +156,9 @@ class SimulationManager { CurrentDensitySplitField &J_s); /** - * @brief 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. + * @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 @@ -160,6 +169,15 @@ class SimulationManager { 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, diff --git a/tdms/src/simulation_manager/execute_simulation.cpp b/tdms/src/simulation_manager/execute_simulation.cpp index 178eec73c..5a11ee413 100644 --- a/tdms/src/simulation_manager/execute_simulation.cpp +++ b/tdms/src/simulation_manager/execute_simulation.cpp @@ -2818,265 +2818,37 @@ void SimulationManager::execute() { } //end parallel if (TIME_EXEC) { timers.click_timer(TimersTrackingLoop::INTERNAL); } - //update terms for self consistency across scattered/total interface - E updates - if (inputs.params.source_mode == SourceMode::steadystate) {//steadystate - complex commonPhase = - exp(-IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_E, 2. * DCPI)); - double commonAmplitude = linear_ramp(time_E); - for (k = (inputs.K0.index); k <= (inputs.K1.index); k++) - for (j = (inputs.J0.index); j <= (inputs.J1.index); j++) { - if (inputs.I0.apply) {//Perform across I0 - - if (!inputs.params.is_multilayer) array_ind = inputs.I0.index - 1; - else - array_ind = (I_tot + 1) * k + inputs.I0.index - 1; - - if (j < (inputs.J1.index)) - inputs.H_s.zx[k][j][(inputs.I0.index) - 1] = - inputs.H_s.zx[k][j][(inputs.I0.index) - 1] + - inputs.D.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][0] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][0])); - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) - inputs.H_s.yx[k][j][(inputs.I0.index) - 1] = - inputs.H_s.yx[k][j][(inputs.I0.index) - 1] - - inputs.D.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][1] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][1])); - } - if (inputs.I1.apply) {//Perform across I1 - - if (!inputs.params.is_multilayer) array_ind = inputs.I1.index; - else - array_ind = (I_tot + 1) * k + inputs.I1.index; - - if (j < (inputs.J1.index)) - inputs.H_s.zx[k][j][(inputs.I1.index)] = - inputs.H_s.zx[k][j][(inputs.I1.index)] - - inputs.D.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][4] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][4])); - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) - inputs.H_s.yx[k][j][(inputs.I1.index)] = - inputs.H_s.yx[k][j][(inputs.I1.index)] + - inputs.D.b.x[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Isource - .real[k - (inputs.K0.index)][j - (inputs.J0.index)][5] + - IMAGINARY_UNIT * - inputs.Isource.imag[k - (inputs.K0.index)] - [j - (inputs.J0.index)][5])); - } - } + /* Update source terms for self consistency across scattered/total interface + * - H updates (use time_E) */ - for (k = (inputs.K0.index); k <= (inputs.K1.index); k++) - for (i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - if (inputs.J0.apply) {//Perform across J0 - - if (!inputs.params.is_multilayer) array_ind = inputs.J0.index; - else - array_ind = (J_tot + 1) * k + inputs.J0.index; - - if (i < (inputs.I1.index)) - inputs.H_s.zy[k][(inputs.J0.index) - 1][i] = - inputs.H_s.zy[k][(inputs.J0.index) - 1][i] - - inputs.D.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][0])); - - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) - inputs.H_s.xy[k][(inputs.J0.index) - 1][i] = - inputs.H_s.xy[k][(inputs.J0.index) - 1][i] + - inputs.D.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][1])); - } - if (inputs.J1.apply) {//Perform across J1 - - if (!inputs.params.is_multilayer) array_ind = inputs.J1.index; - else - array_ind = (J_tot + 1) * k + inputs.J1.index; - - if (i < (inputs.I1.index)) - inputs.H_s.zy[k][(inputs.J1.index)][i] = - inputs.H_s.zy[k][(inputs.J1.index)][i] + - inputs.D.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][4] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][4])); - if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) - inputs.H_s.xy[k][(inputs.J1.index)][i] = - inputs.H_s.xy[k][(inputs.J1.index)][i] - - inputs.D.b.y[array_ind] * - real(commonAmplitude * commonPhase * - (inputs.Jsource - .real[k - (inputs.K0.index)][i - (inputs.I0.index)][5] + - IMAGINARY_UNIT * - inputs.Jsource.imag[k - (inputs.K0.index)] - [i - (inputs.I0.index)][5])); - } - } + if (inputs.params.source_mode == SourceMode::steadystate) { + // Steady-state source term updates + update_source_terms_steadystate(time_E); - for (j = (inputs.J0.index); j <= (inputs.J1.index); j++) - for (i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - if (inputs.K0.apply) {//Perform across K0 - if (i < (inputs.I1.index)) - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][0])); - if (j < (inputs.J1.index)) - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][1])); - } - if (inputs.K1.apply) {//Perform across K1 - if (i < (inputs.I1.index)) - inputs.H_s.yz[(inputs.K1.index)][j][i] = - inputs.H_s.yz[(inputs.K1.index)][j][i] - - inputs.D.b.z[(inputs.K1.index)] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][4] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][4])); - if (j < (inputs.J1.index)) - inputs.H_s.xz[(inputs.K1.index)][j][i] = - inputs.H_s.xz[(inputs.K1.index)][j][i] + - inputs.D.b.z[(inputs.K1.index)] * - real(commonAmplitude * commonPhase * - (inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][5] + - IMAGINARY_UNIT * - inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][5])); - } - } - outputs.E.ft = real(commonAmplitude * commonPhase); - } else if (inputs.params.source_mode == 1) {//pulsed - if (J_tot == 0) { - j = 0; - for (i = 0; i < (I_tot + 1); i++) { - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * inputs.Ksource.imag[0][i - (inputs.I0.index)][1]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); - //broadband source term - if (inputs.params.eyi_present) - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; - } - for (i = 0; i < I_tot; i++) { - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * inputs.Ksource.imag[0][i - (inputs.I0.index)][0]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); - //broadband source term - if (inputs.params.exi_present) - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; - //if(i==511) - } - } else { - for (j = 0; j < J_tot; j++) - for (i = 0; i < (I_tot + 1); i++) { - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][1]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); - //broadband source term - if (inputs.params.eyi_present) - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; - } - for (j = 0; j < (J_tot + 1); j++) - for (i = 0; i < I_tot; i++) { - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource - .real[j - (inputs.J0.index)][i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * inputs.Ksource.imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][0]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); - //broadband source term - if (inputs.params.exi_present) - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; - } - } - outputs.E.ft = - real((-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), 2 * DCPI))) * - exp(-1. * DCPI * pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); + //! Common phase term in update equations + complex common_phase = + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * time_E, 2. * DCPI)); + //! Common amplitude term in update equations + double common_amplitude = linear_ramp(time_E); + // Update output field + outputs.E.ft = real(common_amplitude * common_phase); + } else if (inputs.params.source_mode == SourceMode::pulsed) { + // Pulsed source term updates + update_source_terms_pulsed(time_E, tind); + + //! Common amplitude factor in update equations + double common_amplitude = + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / inputs.params.hwhm, 2.)); + //! Common phase term in update equations + complex common_phase = + -1. * IMAGINARY_UNIT * + exp(-1. * IMAGINARY_UNIT * + fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), + 2 * DCPI)); + // Update output field + outputs.E.ft = common_amplitude * real(common_phase); } if (TIME_EXEC) { timers.click_timer(TimersTrackingLoop::INTERNAL); } diff --git a/tdms/src/simulation_manager/source_updates_pulsed.cpp b/tdms/src/simulation_manager/source_updates_pulsed.cpp index 89cd4268f..42760621f 100644 --- a/tdms/src/simulation_manager/source_updates_pulsed.cpp +++ b/tdms/src/simulation_manager/source_updates_pulsed.cpp @@ -97,3 +97,84 @@ void SimulationManager::update_source_terms_pulsed( } } } + +void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { + // Exit if there are no sources to update + if (inputs.Ksource.is_empty()) { return; } + + //! Simulation dimensions + int I_tot = n_Yee_cells().i, J_tot = n_Yee_cells().j; + //! The material constant that appears in the update equation + double d_constant = inputs.D.b.z[inputs.K0.index - 1]; + //! Common amplitude factor in update equations + double common_amplitude = + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / inputs.params.hwhm, 2.)); + //! Common phase term in update equations + complex common_phase = + -1. * IMAGINARY_UNIT * + exp(-1. * IMAGINARY_UNIT * + fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), + 2 * DCPI)); + //! Source term index to read from + SourceIndex s_index; + //! split-field cell index to update + CellCoordinate cell_to_update; + + if (J_tot == 0) { + int j = 0; + for (int i = 0; i < I_tot + 1; i++) { + s_index = {1, i - inputs.I0.index, j}; + cell_to_update = {i, j, inputs.K0.index - 1}; + + // Update magnetic split field + inputs.H_s.xz[cell_to_update] -= + d_constant * common_amplitude * + real(inputs.Ksource[s_index] * common_phase); + // Update broadband source term + if (inputs.params.eyi_present) { + inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y[tind][j][i]; + } + } + for (int i = 0; i < I_tot; i++) { + s_index = {0, i - inputs.I0.index, j}; + cell_to_update = {i, j, inputs.K0.index - 1}; + + // Update magnetic split field + inputs.H_s.yz[cell_to_update] += + d_constant * common_amplitude * + real(inputs.Ksource[s_index] * common_phase); + // Update broadband source term + if (inputs.params.exi_present) { + inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x[tind][j][i]; + } + } + } else { + for (int j = 0; j < J_tot; j++) { + for (int i = 0; i < I_tot + 1; i++) { + s_index = {1, i - inputs.I0.index, j - inputs.J0.index}; + cell_to_update = {i, j, inputs.K0.index - 1}; + + inputs.H_s.xz[cell_to_update] -= + d_constant * common_amplitude * + real(inputs.Ksource[s_index] * common_phase); + if (inputs.params.eyi_present) { + inputs.H_s.xz[cell_to_update] -= d_constant * inputs.Ei.y[tind][j][i]; + } + } + } + for (int j = 0; j < J_tot + 1; j++) { + for (int i = 0; i < I_tot; i++) { + s_index = {0, i - inputs.I0.index, j - inputs.J0.index}; + cell_to_update = {i, j, inputs.K0.index}; + + inputs.H_s.yz[cell_to_update] += + d_constant * common_amplitude * + real(inputs.Ksource[s_index] * common_phase); + if (inputs.params.exi_present) { + inputs.H_s.yz[cell_to_update] += d_constant * inputs.Ei.x[tind][j][i]; + } + } + } + } +} diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index 0f5cfd857..1111e4455 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -213,7 +213,7 @@ void SimulationManager::update_source_steadystate( * to recover the source value at the relevant Yee cell. cell_b and cell_c need to be offset by the corresponding {IJK}0.index before - we cna retrieve the source value. */ + we can retrieve the source value. */ SourceIndex s_index = {split_field_ID, cell_b, cell_c}; // setup axis-dependent parameters @@ -314,3 +314,205 @@ void SimulationManager::update_source_steadystate( (2 * inputs.params.dt) * E_split_update; } } + +void SimulationManager::update_source_terms_steadystate(double time_E) { + //! Simulation dimensions + int I_tot = n_Yee_cells().i, J_tot = n_Yee_cells().j; + //! The index of the various material property arrays that correspond to this + //! particular Yee cell, and thus update equation. + int array_ind; + //! Common phase factor that appears in update equations + complex common_phase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_E, 2. * DCPI)); + //! Common amplitude factor that appears in update equations + double common_amplitude = linear_ramp(time_E); + + /*! Value of the source term in Yee cell (cell_a, cell_b, cell_c) */ + complex source_value; + /*! The split-field cell index to update */ + CellCoordinate cell_to_update; + /*! This is the index that is to be passed to the (relevent) Source instance, + * to recover the source value at the relevant Yee cell. + + cell_b and cell_c need to be offset by the corresponding {IJK}0.index before + we can retrieve the source value. */ + SourceIndex s_index; + + // TODO: Pull if statements outside loops. Refactor akin to e-field: have + // central function with switch to setup, then call that rather than repeating + // yourself over and over again + + // Update Isource terms + for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + // Update across I0 + if (inputs.I0.apply) { + cell_to_update = {inputs.I0.index - 1, j, k}; + + if (!inputs.params.is_multilayer) { + array_ind = inputs.I0.index - 1; + } else { + array_ind = (I_tot + 1) * k + inputs.I0.index - 1; + } + + if (j < inputs.J1.index) { + s_index = {0, j - inputs.J0.index, k - inputs.K0.index}; + + inputs.H_s.zx[cell_to_update] += + inputs.D.b.x[array_ind] * + real(common_amplitude * common_phase * + inputs.Isource[s_index]); + } + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + s_index = {1, j - inputs.J0.index, k - inputs.K0.index}; + + inputs.H_s.yx[cell_to_update] -= + inputs.D.b.x[array_ind] * + real(common_amplitude * common_phase * + inputs.Isource[s_index]); + } + } + + // Update across I1 + if (inputs.I1.apply) { + cell_to_update = {inputs.I1.index, j, k}; + + if (!inputs.params.is_multilayer) { + array_ind = inputs.I1.index; + } else { + array_ind = (I_tot + 1) * k + inputs.I1.index; + } + + if (j < inputs.J1.index) { + s_index = {4, j - inputs.J0.index, k - inputs.K0.index}; + + inputs.H_s.zx[cell_to_update] -= + inputs.D.b.x[array_ind] * + real(common_amplitude * common_phase * + inputs.Isource[s_index]); + } + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + s_index = {5, j - inputs.J0.index, k - inputs.K0.index}; + + inputs.H_s.yx[cell_to_update] += + inputs.D.b.x[array_ind] * + real(common_amplitude * common_phase * + inputs.Isource[s_index]); + } + } + } + } + + // Update Jsource terms + for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { + for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { + // Perform across J0 + if (inputs.J0.apply) { + cell_to_update = {i, inputs.J0.index - 1, k}; + + if (!inputs.params.is_multilayer) { + array_ind = inputs.J0.index; + } else { + array_ind = (J_tot + 1) * k + inputs.J0.index; + } + + if (i < inputs.I1.index) { + s_index = {0, inputs.I0.index, inputs.K0.index}; + + inputs.H_s.zy[s_index] -= inputs.D.b.y[array_ind] * + real(common_amplitude * common_phase * + inputs.Jsource[s_index]); + } + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + s_index = {1, inputs.I0.index, inputs.K0.index}; + + inputs.H_s.xy[cell_to_update] += + inputs.D.b.y[array_ind] * + real(common_amplitude * common_phase * + inputs.Jsource[s_index]); + } + } + + // Perform across J1 + if (inputs.J1.apply) { + cell_to_update = {i, inputs.J1.index, k}; + + if (!inputs.params.is_multilayer) { + array_ind = inputs.J1.index; + } else { + array_ind = (J_tot + 1) * k + inputs.J1.index; + } + + if (i < (inputs.I1.index)) { + s_index = {4, i - inputs.I0.index, k - inputs.K0.index}; + + inputs.H_s.zy[cell_to_update] += + inputs.D.b.y[array_ind] * + real(common_amplitude * common_phase * + inputs.Jsource[s_index]); + } + if (k < (inputs.K1.index) || + inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { + s_index = {5, i - inputs.I0.index, k - inputs.K0.index}; + + inputs.H_s.xy[cell_to_update] -= + inputs.D.b.y[array_ind] * + real(common_amplitude * common_phase * + inputs.Jsource[s_index]); + } + } + } + } + + // Update Ksource terms + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { + // Perform across K0 + if (inputs.K0.apply) { + cell_to_update = {i, j, inputs.K0.index - 1}; + + if (i < inputs.I1.index) { + s_index = {0, i - inputs.I0.index, j - inputs.J0.index}; + + inputs.H_s.yz[cell_to_update] += + inputs.D.b.z[inputs.K0.index - 1] * + real(common_amplitude * common_phase * + inputs.Ksource[s_index]); + } + if (j < inputs.J1.index) { + s_index = {1, i - inputs.I0.index, j - inputs.J0.index}; + + inputs.H_s.xz[cell_to_update] -= + inputs.D.b.z[inputs.K0.index - 1] * + real(common_amplitude * common_phase * + inputs.Ksource[s_index]); + } + } + + // Perform across K1 + if (inputs.K1.apply) { + cell_to_update = {i, j, inputs.K1.index}; + + if (i < inputs.I1.index) { + s_index = {4, i - inputs.I0.index, j - inputs.J0.index}; + + inputs.H_s.yz[cell_to_update] -= + inputs.D.b.z[inputs.K1.index] * + real(common_amplitude * common_phase * + inputs.Ksource[s_index]); + } + if (j < (inputs.J1.index)) { + s_index = {5, i - inputs.I0.index, j - inputs.J0.index}; + + inputs.H_s.xz[cell_to_update] += + inputs.D.b.z[inputs.K1.index] * + real(common_amplitude * common_phase * + inputs.Ksource[s_index]); + } + } + } + } +} \ No newline at end of file From 553f75acfd0b61ffb7b3b9f9584438fef987a734 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 15:49:20 +0000 Subject: [PATCH 18/25] Fully refactor H-field, now to find the bug... --- .../simulation_manager/simulation_manager.h | 73 +++- .../simulation_manager/execute_simulation.cpp | 4 +- .../source_updates_steadystate.cpp | 386 +++++++++--------- 3 files changed, 257 insertions(+), 206 deletions(-) diff --git a/tdms/include/simulation_manager/simulation_manager.h b/tdms/include/simulation_manager/simulation_manager.h index df966cfb2..9407fd6d7 100644 --- a/tdms/include/simulation_manager/simulation_manager.h +++ b/tdms/include/simulation_manager/simulation_manager.h @@ -93,8 +93,9 @@ class SimulationManager { void prepare_output(const mxArray *fieldsample, const mxArray *campssample); /** - * @brief Update electric-split field components and current densities in - * steady-state after an E-field timestep has been performed. + * @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. @@ -114,11 +115,33 @@ class SimulationManager { * @param J_c The current density in the conductive medium * @param J_s The current density in the dispersive medium */ - void update_source_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); + 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 @@ -130,7 +153,7 @@ class SimulationManager { * @param J_c The current density in the conductive medium * @param J_s The current density in the dispersive medium */ - void update_source_terms_steadystate(double time_H, bool is_conductive, + void E_source_update_all_steadystate(double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s); /** @@ -140,20 +163,26 @@ class SimulationManager { * * @param time_E The time the electric field is currently sitting at. */ - void update_source_terms_steadystate(double time_E); - - /*! @copydoc update_source_terms_steadystate */ - void update_Isource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s); - /*! @copydoc update_source_terms_steadystate */ - void update_Jsource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s); - /*! @copydoc update_source_terms_steadystate */ - void update_Ksource_terms_steadystate(double time_H, bool is_conductive, - CurrentDensitySplitField &J_c, - CurrentDensitySplitField &J_s); + 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 diff --git a/tdms/src/simulation_manager/execute_simulation.cpp b/tdms/src/simulation_manager/execute_simulation.cpp index 5a11ee413..dd3575bf2 100644 --- a/tdms/src/simulation_manager/execute_simulation.cpp +++ b/tdms/src/simulation_manager/execute_simulation.cpp @@ -2139,7 +2139,7 @@ void SimulationManager::execute() { if (inputs.params.source_mode == SourceMode::steadystate) { // Steady-state source term updates - update_source_terms_steadystate(time_H, loop_variables.is_conductive, + E_source_update_all_steadystate(time_H, loop_variables.is_conductive, loop_variables.J_c, loop_variables.J_s); // Common phase term in update equations @@ -2823,7 +2823,7 @@ void SimulationManager::execute() { if (inputs.params.source_mode == SourceMode::steadystate) { // Steady-state source term updates - update_source_terms_steadystate(time_E); + H_source_update_all_steadystate(time_E); //! Common phase term in update equations complex common_phase = diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index 1111e4455..33754c73b 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -6,15 +6,15 @@ using namespace std; using namespace tdms_phys_constants; using tdms_math_constants::DCPI, tdms_math_constants::IMAGINARY_UNIT; -void SimulationManager::update_source_terms_steadystate( +void SimulationManager::E_source_update_all_steadystate( double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s) { - update_Isource_terms_steadystate(time_H, is_conductive, J_c, J_s); - update_Jsource_terms_steadystate(time_H, is_conductive, J_c, J_s); - update_Ksource_terms_steadystate(time_H, is_conductive, J_c, J_s); + E_Isource_update_steadystate(time_H, is_conductive, J_c, J_s); + E_Jsource_update_steadystate(time_H, is_conductive, J_c, J_s); + E_Ksource_update_steadystate(time_H, is_conductive, J_c, J_s); } -void SimulationManager::update_Isource_terms_steadystate( +void SimulationManager::E_Isource_update_steadystate( double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided @@ -33,12 +33,12 @@ void SimulationManager::update_Isource_terms_steadystate( } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - update_source_steadystate(time_H, AxialDirection::X, true, true, - is_conductive, j, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::X, true, true, + is_conductive, j, k, array_ind, J_c, J_s); } if (j < (inputs.J1.index)) { - update_source_steadystate(time_H, AxialDirection::X, false, true, - is_conductive, j, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::X, false, true, + is_conductive, j, k, array_ind, J_c, J_s); } } } @@ -55,19 +55,19 @@ void SimulationManager::update_Isource_terms_steadystate( } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - update_source_steadystate(time_H, AxialDirection::X, true, false, - is_conductive, j, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::X, true, false, + is_conductive, j, k, array_ind, J_c, J_s); } if (j < (inputs.J1.index)) { - update_source_steadystate(time_H, AxialDirection::X, false, false, - is_conductive, j, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::X, false, false, + is_conductive, j, k, array_ind, J_c, J_s); } } } } } -void SimulationManager::update_Jsource_terms_steadystate( +void SimulationManager::E_Jsource_update_steadystate( double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided @@ -89,12 +89,12 @@ void SimulationManager::update_Jsource_terms_steadystate( array_ind = (n_Yee_cells().j + 1) * k + inputs.J0.index; } - update_source_steadystate(time_H, AxialDirection::Y, true, true, - is_conductive, i, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::Y, true, true, + is_conductive, i, k, array_ind, J_c, J_s); } if (i < inputs.I1.index) { - update_source_steadystate(time_H, AxialDirection::Y, false, true, - is_conductive, i, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::Y, false, true, + is_conductive, i, k, array_ind, J_c, J_s); } } } @@ -111,19 +111,19 @@ void SimulationManager::update_Jsource_terms_steadystate( } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - update_source_steadystate(time_H, AxialDirection::Y, true, false, - is_conductive, i, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::Y, true, false, + is_conductive, i, k, array_ind, J_c, J_s); } if (i < (inputs.I1.index)) { - update_source_steadystate(time_H, AxialDirection::Y, false, false, - is_conductive, i, k, array_ind, J_c, J_s); + E_source_update_steadystate(time_H, AxialDirection::Y, false, false, + is_conductive, i, k, array_ind, J_c, J_s); } } } } } -void SimulationManager::update_Ksource_terms_steadystate( +void SimulationManager::E_Ksource_update_steadystate( double time_H, bool is_conductive, CurrentDensitySplitField &J_c, CurrentDensitySplitField &J_s) { // Only run update equations is source data was provided @@ -134,14 +134,14 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - update_source_steadystate(time_H, AxialDirection::Z, true, true, - is_conductive, i, j, inputs.K0.index, J_c, - J_s); + E_source_update_steadystate(time_H, AxialDirection::Z, true, true, + is_conductive, i, j, inputs.K0.index, J_c, + J_s); } if (i < (inputs.I1.index)) { - update_source_steadystate(time_H, AxialDirection::Z, false, true, - is_conductive, i, j, inputs.K0.index, J_c, - J_s); + E_source_update_steadystate(time_H, AxialDirection::Z, false, true, + is_conductive, i, j, inputs.K0.index, J_c, + J_s); } } } @@ -152,21 +152,21 @@ void SimulationManager::update_Ksource_terms_steadystate( for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { if (j < (inputs.J1.index)) { - update_source_steadystate(time_H, AxialDirection::Z, true, false, - is_conductive, i, j, inputs.K1.index, J_c, - J_s); + E_source_update_steadystate(time_H, AxialDirection::Z, true, false, + is_conductive, i, j, inputs.K1.index, J_c, + J_s); } if (i < (inputs.I1.index)) { - update_source_steadystate(time_H, AxialDirection::Z, false, false, - is_conductive, i, j, inputs.K1.index, J_c, - J_s); + E_source_update_steadystate(time_H, AxialDirection::Z, false, false, + is_conductive, i, j, inputs.K1.index, J_c, + J_s); } } } } } -void SimulationManager::update_source_steadystate( +void SimulationManager::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) { @@ -315,204 +315,226 @@ void SimulationManager::update_source_steadystate( } } -void SimulationManager::update_source_terms_steadystate(double time_E) { - //! Simulation dimensions - int I_tot = n_Yee_cells().i, J_tot = n_Yee_cells().j; - //! The index of the various material property arrays that correspond to this - //! particular Yee cell, and thus update equation. - int array_ind; - //! Common phase factor that appears in update equations - complex common_phase = exp( - -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_E, 2. * DCPI)); - //! Common amplitude factor that appears in update equations - double common_amplitude = linear_ramp(time_E); - - /*! Value of the source term in Yee cell (cell_a, cell_b, cell_c) */ - complex source_value; - /*! The split-field cell index to update */ - CellCoordinate cell_to_update; - /*! This is the index that is to be passed to the (relevent) Source instance, - * to recover the source value at the relevant Yee cell. - - cell_b and cell_c need to be offset by the corresponding {IJK}0.index before - we can retrieve the source value. */ - SourceIndex s_index; - - // TODO: Pull if statements outside loops. Refactor akin to e-field: have - // central function with switch to setup, then call that rather than repeating - // yourself over and over again +void SimulationManager::H_source_update_all_steadystate(double time_E) { + H_Isource_update_steadystate(time_E); + H_Jsource_update_steadystate(time_E); + H_Ksource_update_steadystate(time_E); +} - // Update Isource terms - for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { - for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { - // Update across I0 - if (inputs.I0.apply) { - cell_to_update = {inputs.I0.index - 1, j, k}; +void SimulationManager::H_Isource_update_steadystate(double time_E) { + if (inputs.Isource.is_empty()) { return; } + int array_ind; - if (!inputs.params.is_multilayer) { - array_ind = inputs.I0.index - 1; - } else { - array_ind = (I_tot + 1) * k + inputs.I0.index - 1; + // Update across I0 + array_ind = inputs.I0.index - 1; + if (inputs.I0.apply) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { + if (inputs.params.is_multilayer) { + array_ind = (n_Yee_cells().i + 1) * k + inputs.I0.index - 1; } if (j < inputs.J1.index) { - s_index = {0, j - inputs.J0.index, k - inputs.K0.index}; - - inputs.H_s.zx[cell_to_update] += - inputs.D.b.x[array_ind] * - real(common_amplitude * common_phase * - inputs.Isource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::X, false, true, + array_ind, j, k); } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - s_index = {1, j - inputs.J0.index, k - inputs.K0.index}; - - inputs.H_s.yx[cell_to_update] -= - inputs.D.b.x[array_ind] * - real(common_amplitude * common_phase * - inputs.Isource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::X, true, true, + array_ind, j, k); } } - - // Update across I1 - if (inputs.I1.apply) { - cell_to_update = {inputs.I1.index, j, k}; - - if (!inputs.params.is_multilayer) { - array_ind = inputs.I1.index; - } else { - array_ind = (I_tot + 1) * k + inputs.I1.index; + } + } + // Update across I1 + array_ind = inputs.I1.index; + if (inputs.I1.apply) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int j = inputs.J0.index; j <= inputs.J1.index; j++) { + if (inputs.params.is_multilayer) { + array_ind = (n_Yee_cells().i + 1) * k + inputs.I1.index; } if (j < inputs.J1.index) { - s_index = {4, j - inputs.J0.index, k - inputs.K0.index}; - - inputs.H_s.zx[cell_to_update] -= - inputs.D.b.x[array_ind] * - real(common_amplitude * common_phase * - inputs.Isource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::X, false, false, + array_ind, j, k); } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - s_index = {5, j - inputs.J0.index, k - inputs.K0.index}; - - inputs.H_s.yx[cell_to_update] += - inputs.D.b.x[array_ind] * - real(common_amplitude * common_phase * - inputs.Isource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::X, true, false, + array_ind, j, k); } } } } +} - // Update Jsource terms - for (int k = (inputs.K0.index); k <= (inputs.K1.index); k++) { - for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - // Perform across J0 - if (inputs.J0.apply) { - cell_to_update = {i, inputs.J0.index - 1, k}; +void SimulationManager::H_Jsource_update_steadystate(double time_E) { + if (inputs.Jsource.is_empty()) { return; } + int array_ind; - if (!inputs.params.is_multilayer) { - array_ind = inputs.J0.index; - } else { - array_ind = (J_tot + 1) * k + inputs.J0.index; + // Perform across J0 + array_ind = inputs.J0.index; + if (inputs.J0.apply) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { + if (inputs.params.is_multilayer) { + array_ind = (n_Yee_cells().j + 1) * k + inputs.J0.index; } if (i < inputs.I1.index) { - s_index = {0, inputs.I0.index, inputs.K0.index}; - - inputs.H_s.zy[s_index] -= inputs.D.b.y[array_ind] * - real(common_amplitude * common_phase * - inputs.Jsource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Y, false, true, + array_ind, i, k); } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - s_index = {1, inputs.I0.index, inputs.K0.index}; - - inputs.H_s.xy[cell_to_update] += - inputs.D.b.y[array_ind] * - real(common_amplitude * common_phase * - inputs.Jsource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Y, true, true, + array_ind, i, k); } } - - // Perform across J1 - if (inputs.J1.apply) { - cell_to_update = {i, inputs.J1.index, k}; - - if (!inputs.params.is_multilayer) { - array_ind = inputs.J1.index; - } else { - array_ind = (J_tot + 1) * k + inputs.J1.index; + } + } + // Perform across J1 + array_ind = inputs.J1.index; + if (inputs.J1.apply) { + for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { + for (int i = inputs.I0.index; i <= inputs.I1.index; i++) { + if (inputs.params.is_multilayer) { + array_ind = (n_Yee_cells().j + 1) * k + inputs.J1.index; } if (i < (inputs.I1.index)) { - s_index = {4, i - inputs.I0.index, k - inputs.K0.index}; - - inputs.H_s.zy[cell_to_update] += - inputs.D.b.y[array_ind] * - real(common_amplitude * common_phase * - inputs.Jsource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Y, false, false, + array_ind, i, k); } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - s_index = {5, i - inputs.I0.index, k - inputs.K0.index}; - - inputs.H_s.xy[cell_to_update] -= - inputs.D.b.y[array_ind] * - real(common_amplitude * common_phase * - inputs.Jsource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Y, true, false, + array_ind, i, k); } } } } +} - // Update Ksource terms - for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { - for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { - // Perform across K0 - if (inputs.K0.apply) { - cell_to_update = {i, j, inputs.K0.index - 1}; +void SimulationManager::H_Ksource_update_steadystate(double time_E) { + if (inputs.Ksource.is_empty()) { return; } + // Perform across K0 + if (inputs.K0.apply) { + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { + // Perform across K0 if (i < inputs.I1.index) { - s_index = {0, i - inputs.I0.index, j - inputs.J0.index}; - - inputs.H_s.yz[cell_to_update] += - inputs.D.b.z[inputs.K0.index - 1] * - real(common_amplitude * common_phase * - inputs.Ksource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Z, false, true, + inputs.K0.index - 1, i, j); } if (j < inputs.J1.index) { - s_index = {1, i - inputs.I0.index, j - inputs.J0.index}; - - inputs.H_s.xz[cell_to_update] -= - inputs.D.b.z[inputs.K0.index - 1] * - real(common_amplitude * common_phase * - inputs.Ksource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Z, true, true, + inputs.K0.index - 1, i, j); } } - - // Perform across K1 - if (inputs.K1.apply) { - cell_to_update = {i, j, inputs.K1.index}; - + } + } + // Perform across K1 + if (inputs.K1.apply) { + for (int j = (inputs.J0.index); j <= (inputs.J1.index); j++) { + for (int i = (inputs.I0.index); i <= (inputs.I1.index); i++) { if (i < inputs.I1.index) { - s_index = {4, i - inputs.I0.index, j - inputs.J0.index}; - - inputs.H_s.yz[cell_to_update] -= - inputs.D.b.z[inputs.K1.index] * - real(common_amplitude * common_phase * - inputs.Ksource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Z, false, false, + inputs.K1.index, i, j); } if (j < (inputs.J1.index)) { - s_index = {5, i - inputs.I0.index, j - inputs.J0.index}; - - inputs.H_s.xz[cell_to_update] += - inputs.D.b.z[inputs.K1.index] * - real(common_amplitude * common_phase * - inputs.Ksource[s_index]); + H_source_update_steadystate(time_E, AxialDirection::Z, true, false, + inputs.K1.index, i, j); } } } } +} + +void SimulationManager::H_source_update_steadystate( + double time_E, AxialDirection parallel, bool C_axis, bool zero_plane, + int array_ind, int cell_b, int cell_c) { + //! Common phase factor that appears in update equations + complex common_phase = exp( + -IMAGINARY_UNIT * fmod(inputs.params.omega_an * time_E, 2. * DCPI)); + //! Common amplitude factor that appears in update equations + double common_amplitude = linear_ramp(time_E); + + /*! D.b.{x,y,z} constant that appears in the update equations */ + double d_constant; + /*! Value of the source term in Yee cell (cell_a, cell_b, cell_c) */ + complex source_value; + /*! The split-field cell index to update */ + CellCoordinate cell_to_update; + /*! This is the index that is to be passed to the (relevent) Source instance, + * to recover the source value at the relevant Yee cell. + + cell_b and cell_c need to be offset by the corresponding {IJK}0.index before + we can retrieve the source value. */ + SourceIndex s_index = {0, cell_b, cell_c}; + // If working along the C-axis, split_field_ID is one greater than the B-axis + if (C_axis) { s_index.i += 1; } + // If working along the {IJK}1 plane, split_field_Id is 4 greater than the + // 0-plane + if (!zero_plane) { s_index.i += 4; } + + /*! Field updates are signed based on the component, plane, and axis */ + double update_sign = 1.; + + SplitFieldComponent + *H_s_component;//< The H_s component that should be updated + + switch (parallel) { + case AxialDirection::X: + // Dealing with Isource + d_constant = inputs.D.b.x[array_ind]; + + s_index.j -= inputs.J0.index; + s_index.k -= inputs.K0.index; + source_value = inputs.Isource[s_index]; + + cell_to_update.i = (zero_plane) ? inputs.I0.index - 1 : inputs.I1.index; + cell_to_update.j = cell_b; + cell_to_update.k = cell_c; + + H_s_component = (C_axis) ? &inputs.H_s.yx : &inputs.H_s.zx; + if (C_axis == zero_plane) { update_sign = -1.; } + break; + case AxialDirection::Y: + // Dealing with Jsource + d_constant = inputs.D.b.y[array_ind]; + + s_index.j -= inputs.I0.index; + s_index.k -= inputs.K0.index; + source_value = inputs.Jsource[s_index]; + + cell_to_update.i = cell_b; + cell_to_update.j = (zero_plane) ? inputs.J0.index - 1 : inputs.J1.index; + cell_to_update.k = cell_c; + + H_s_component = (C_axis) ? &inputs.H_s.xy : &inputs.H_s.zy; + if (C_axis != zero_plane) { update_sign = -1.; } + break; + case AxialDirection::Z: + // Dealing wtih Ksource + d_constant = inputs.D.b.z[array_ind]; + + s_index.j -= inputs.I0.index; + s_index.k -= inputs.J0.index; + source_value = inputs.Ksource[s_index]; + + cell_to_update.i = cell_b; + cell_to_update.j = cell_c; + cell_to_update.k = (zero_plane) ? inputs.K0.index - 1 : inputs.K1.index; + + H_s_component = (C_axis) ? &inputs.H_s.xz : &inputs.H_s.yz; + if (C_axis == zero_plane) { update_sign = -1.; } + break; + } + + (*H_s_component)[cell_to_update] += + update_sign * d_constant * + real(common_amplitude * common_phase * source_value); } \ No newline at end of file From b79d266ad774e3920bb74ee82b848174de85bff0 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 15:59:08 +0000 Subject: [PATCH 19/25] Add that newline to make pre-commit happy --- tdms/src/simulation_manager/source_updates_steadystate.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index 33754c73b..e55513db7 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -537,4 +537,4 @@ void SimulationManager::H_source_update_steadystate( (*H_s_component)[cell_to_update] += update_sign * d_constant * real(common_amplitude * common_phase * source_value); -} \ No newline at end of file +} From 9bf6d11a014c1ed5084b27134018f7c89f5b0fb6 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 16:30:02 +0000 Subject: [PATCH 20/25] Missed a -1 in ONE place... does this fix the bug? --- .../source_updates_pulsed.cpp | 120 +++++++++++++++++- 1 file changed, 119 insertions(+), 1 deletion(-) diff --git a/tdms/src/simulation_manager/source_updates_pulsed.cpp b/tdms/src/simulation_manager/source_updates_pulsed.cpp index 42760621f..4f51974b5 100644 --- a/tdms/src/simulation_manager/source_updates_pulsed.cpp +++ b/tdms/src/simulation_manager/source_updates_pulsed.cpp @@ -166,7 +166,7 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { for (int j = 0; j < J_tot + 1; j++) { for (int i = 0; i < I_tot; i++) { s_index = {0, i - inputs.I0.index, j - inputs.J0.index}; - cell_to_update = {i, j, inputs.K0.index}; + cell_to_update = {i, j, inputs.K0.index - 1}; inputs.H_s.yz[cell_to_update] += d_constant * common_amplitude * @@ -178,3 +178,121 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { } } } + +/* +void original_code() { + if (J_tot == 0) { + j = 0; + for (i = 0; i < (I_tot + 1); i++) { + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - + inputs.D.b.z[(inputs.K0.index) - 1] * + real((inputs.Ksource.real[0][i - (inputs.I0.index)][1] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[0][i - (inputs.I0.index)] + [1]) * + (-1. * IMAGINARY_UNIT) * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * + (time_E - inputs.params.to_l), + 2 * DCPI))) * + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / + (inputs.params.hwhm), + 2.)); + // broadband source term + if (inputs.params.eyi_present) + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - + inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; + } + for (i = 0; i < I_tot; i++) { + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + + inputs.D.b.z[(inputs.K0.index) - 1] * + real((inputs.Ksource.real[0][i - (inputs.I0.index)][0] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[0][i - (inputs.I0.index)] + [0]) * + (-1. * IMAGINARY_UNIT) * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * + (time_E - inputs.params.to_l), + 2 * DCPI))) * + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / + (inputs.params.hwhm), + 2.)); + // broadband source term + if (inputs.params.exi_present) + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + + inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; + // if(i==511) + } + } else { + for (j = 0; j < J_tot; j++) + for (i = 0; i < (I_tot + 1); i++) { + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - + inputs.D.b.z[(inputs.K0.index) - 1] * + real((inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][1] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][1]) * + (-1. * IMAGINARY_UNIT) * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * + (time_E - inputs.params.to_l), + 2 * DCPI))) * + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / + (inputs.params.hwhm), + 2.)); + // broadband source term + if (inputs.params.eyi_present) + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - + inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; + } + for (j = 0; j < (J_tot + 1); j++) + for (i = 0; i < I_tot; i++) { + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + + inputs.D.b.z[(inputs.K0.index) - 1] * + real((inputs.Ksource.real[j - (inputs.J0.index)] + [i - (inputs.I0.index)][0] + + IMAGINARY_UNIT * + inputs.Ksource + .imag[j - (inputs.J0.index)] + [i - (inputs.I0.index)][0]) * + (-1. * IMAGINARY_UNIT) * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * + (time_E - inputs.params.to_l), + 2 * DCPI))) * + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / + (inputs.params.hwhm), + 2.)); + // broadband source term + if (inputs.params.exi_present) + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = + inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + + inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; + } + } + outputs.E.ft = + real((-1. * IMAGINARY_UNIT) * + exp(-IMAGINARY_UNIT * + fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), + 2 * DCPI))) * + exp(-1. * DCPI * + pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); +} + +*/ \ No newline at end of file From 15de02fc9cf687536b2c9e9271b2b4f987fc495a Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Fri, 3 Feb 2023 16:49:32 +0000 Subject: [PATCH 21/25] Remove useless comment --- .../source_updates_pulsed.cpp | 118 ------------------ 1 file changed, 118 deletions(-) diff --git a/tdms/src/simulation_manager/source_updates_pulsed.cpp b/tdms/src/simulation_manager/source_updates_pulsed.cpp index 4f51974b5..3529eae32 100644 --- a/tdms/src/simulation_manager/source_updates_pulsed.cpp +++ b/tdms/src/simulation_manager/source_updates_pulsed.cpp @@ -178,121 +178,3 @@ void SimulationManager::update_source_terms_pulsed(double time_E, int tind) { } } } - -/* -void original_code() { - if (J_tot == 0) { - j = 0; - for (i = 0; i < (I_tot + 1); i++) { - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[0][i - (inputs.I0.index)] - [1]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * - (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / - (inputs.params.hwhm), - 2.)); - // broadband source term - if (inputs.params.eyi_present) - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; - } - for (i = 0; i < I_tot; i++) { - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[0][i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[0][i - (inputs.I0.index)] - [0]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * - (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / - (inputs.params.hwhm), - 2.)); - // broadband source term - if (inputs.params.exi_present) - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; - // if(i==511) - } - } else { - for (j = 0; j < J_tot; j++) - for (i = 0; i < (I_tot + 1); i++) { - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][1] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][1]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * - (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / - (inputs.params.hwhm), - 2.)); - // broadband source term - if (inputs.params.eyi_present) - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.xz[(inputs.K0.index) - 1][j][i] - - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.y[tind][j][i]; - } - for (j = 0; j < (J_tot + 1); j++) - for (i = 0; i < I_tot; i++) { - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * - real((inputs.Ksource.real[j - (inputs.J0.index)] - [i - (inputs.I0.index)][0] + - IMAGINARY_UNIT * - inputs.Ksource - .imag[j - (inputs.J0.index)] - [i - (inputs.I0.index)][0]) * - (-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * - (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / - (inputs.params.hwhm), - 2.)); - // broadband source term - if (inputs.params.exi_present) - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] = - inputs.H_s.yz[(inputs.K0.index) - 1][j][i] + - inputs.D.b.z[(inputs.K0.index) - 1] * inputs.Ei.x[tind][j][i]; - } - } - outputs.E.ft = - real((-1. * IMAGINARY_UNIT) * - exp(-IMAGINARY_UNIT * - fmod(inputs.params.omega_an * (time_E - inputs.params.to_l), - 2 * DCPI))) * - exp(-1. * DCPI * - pow((time_E - inputs.params.to_l) / (inputs.params.hwhm), 2.)); -} - -*/ \ No newline at end of file From 577cdfc45bb91c954012fcb188e11624b933449c Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Wed, 8 Feb 2023 09:22:39 +0000 Subject: [PATCH 22/25] Apply Peter's suggestion on #226 --- .../source_updates_steadystate.cpp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index e55513db7..decb1c0e1 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -79,16 +79,13 @@ void SimulationManager::E_Jsource_update_steadystate( if (inputs.J0.apply) { for (int k = inputs.K0.index; k <= inputs.K1.index; k++) { for (int i = inputs.I0.index; i < inputs.I1.index; i++) { + if (!inputs.params.is_multilayer) { + array_ind = inputs.J0.index; + } else { + array_ind = (n_Yee_cells().j + 1) * k + inputs.J0.index; + } if (k < (inputs.K1.index) || inputs.params.dimension == Dimension::TRANSVERSE_MAGNETIC) { - // shouldn't this be one-level up to be consistent with the other - // source update steps?!?!?!?! - if (!inputs.params.is_multilayer) { - array_ind = inputs.J0.index; - } else { - array_ind = (n_Yee_cells().j + 1) * k + inputs.J0.index; - } - E_source_update_steadystate(time_H, AxialDirection::Y, true, true, is_conductive, i, k, array_ind, J_c, J_s); } From 8b18880d3bb537c80cde596b6086c0ae3fc3da70 Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Wed, 8 Feb 2023 09:39:21 +0000 Subject: [PATCH 23/25] Add Peter's resolution to #225 --- .../source_updates_steadystate.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index decb1c0e1..b93ac2635 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -180,11 +180,18 @@ void SimulationManager::E_source_update_steadystate( /*! Dispersive parameter. = inputs.matched_layer.kappa.{x,y,z} */ double kappa; /*! NOTE: if dealing with Ksource, then the index used to fetch the gamma - * parameter is k = inputs.K1.index, rather than cell_c. This inconsistency & - * potential error is flagged in #221. Have a feeling that, in the K-case, - * this should be inputs.K{0,1}.index depending on the plane we're in */ - int gamma_ind = (parallel == AxialDirection::Z) ? inputs.K1.index : cell_c; - /*! Dispersive parameter. = inputs.matched_layer.gamma[k]. */ + * parameter is k = inputs.K1.index or inputs.K0.index, depending on whether + * we are updating in the zero_plnae or not. If dealing with Isource or + * Jsource, cell_c (which = the k index) is the index for the gamma parameter. + * See https://github.com/UCL/TDMS/issues/225 + */ + int gamma_ind; + if (parallel == AxialDirection::Z) { + gamma_ind = (zero_plane) ? inputs.K0.index : inputs.K1.index; + } else { + gamma_ind = cell_c; + } + /*! Dispersive parameter. Pulls values from inputs.matched_layer.gamma */ double gamma = inputs.matched_layer.gamma[gamma_ind]; /*! Value of the source term in Yee cell (cell_a, cell_b, cell_c) */ From 7ece7e926c509d5f0207a1d0dd104273ca56f39f Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 23 Feb 2023 09:37:38 +0000 Subject: [PATCH 24/25] Apply Sam's suggestions for SourceIndex struct --- tdms/include/source.h | 17 ++++++----- .../source_updates_steadystate.cpp | 29 ++++++++++--------- 2 files changed, 25 insertions(+), 21 deletions(-) diff --git a/tdms/include/source.h b/tdms/include/source.h index fc695e2bc..a86c07bb9 100644 --- a/tdms/include/source.h +++ b/tdms/include/source.h @@ -16,12 +16,14 @@ * in a Source instance is indexed via * {real,imag}[cell_c][cell_b][split_field_ID]. * - * SourceIndex overlays the ijk struct for code re-use, via the association: - * ijk.k => SourceIndex.cell_c - * ijk.j => SourceIndex.cell_b - * ijk.i => SourceIndex.split_field_ID + * SourceIndex is a container for indexing the Source.{real,imag} objects in an + * efficient and notationally-consistent manner. */ -typedef ijk SourceIndex; +struct SourceIndex { + int split_field_ID = 0,//!< Value in the i/x direction + cell_b = 0, //!< Value in the j/y direction + cell_c = 0; //!< Value in the k/z direction +}; /** * @brief The Source class stores values of the Source field across a particular @@ -58,7 +60,8 @@ class Source { bool is_empty() { return no_data_stored; } std::complex operator[](SourceIndex index) { - return std::complex(real[index.k][index.j][index.i], - imag[index.k][index.j][index.i]); + return std::complex( + real[index.cell_c][index.cell_b][index.split_field_ID], + imag[index.cell_c][index.cell_b][index.split_field_ID]); } }; diff --git a/tdms/src/simulation_manager/source_updates_steadystate.cpp b/tdms/src/simulation_manager/source_updates_steadystate.cpp index b93ac2635..9c880ac32 100644 --- a/tdms/src/simulation_manager/source_updates_steadystate.cpp +++ b/tdms/src/simulation_manager/source_updates_steadystate.cpp @@ -1,6 +1,7 @@ #include "simulation_manager/simulation_manager.h" #include "cell_coordinate.h" +#include "source.h" using namespace std; using namespace tdms_phys_constants; @@ -226,8 +227,8 @@ void SimulationManager::E_source_update_steadystate( // Dealing with Isource. j = cell_b, k = cell_c c_constant = inputs.C.b.x[array_ind]; - s_index.j -= inputs.J0.index; - s_index.k -= inputs.K0.index; + s_index.cell_b -= inputs.J0.index; + s_index.cell_c -= inputs.K0.index; source_value = inputs.Isource[s_index]; cell_to_update = {inputs.I0.index, cell_b, cell_c}; @@ -251,8 +252,8 @@ void SimulationManager::E_source_update_steadystate( // Dealing with Jsource. i = cell_b, k = cell_c c_constant = inputs.C.b.y[array_ind]; - s_index.j -= inputs.I0.index; - s_index.k -= inputs.K0.index; + s_index.cell_b -= inputs.I0.index; + s_index.cell_c -= inputs.K0.index; source_value = inputs.Jsource[s_index]; cell_to_update = {cell_b, inputs.J0.index, cell_c}; @@ -278,8 +279,8 @@ void SimulationManager::E_source_update_steadystate( // Dealing with Ksource. i = cell_b, j = cell_c c_constant = inputs.C.b.z[array_ind]; - s_index.j -= inputs.I0.index; - s_index.k -= inputs.J0.index; + s_index.cell_b -= inputs.I0.index; + s_index.cell_c -= inputs.J0.index; source_value = inputs.Ksource[s_index]; cell_to_update = {cell_b, cell_c, inputs.K0.index}; @@ -479,10 +480,10 @@ void SimulationManager::H_source_update_steadystate( we can retrieve the source value. */ SourceIndex s_index = {0, cell_b, cell_c}; // If working along the C-axis, split_field_ID is one greater than the B-axis - if (C_axis) { s_index.i += 1; } + if (C_axis) { s_index.split_field_ID += 1; } // If working along the {IJK}1 plane, split_field_Id is 4 greater than the // 0-plane - if (!zero_plane) { s_index.i += 4; } + if (!zero_plane) { s_index.split_field_ID += 4; } /*! Field updates are signed based on the component, plane, and axis */ double update_sign = 1.; @@ -495,8 +496,8 @@ void SimulationManager::H_source_update_steadystate( // Dealing with Isource d_constant = inputs.D.b.x[array_ind]; - s_index.j -= inputs.J0.index; - s_index.k -= inputs.K0.index; + s_index.cell_b -= inputs.J0.index; + s_index.cell_c -= inputs.K0.index; source_value = inputs.Isource[s_index]; cell_to_update.i = (zero_plane) ? inputs.I0.index - 1 : inputs.I1.index; @@ -510,8 +511,8 @@ void SimulationManager::H_source_update_steadystate( // Dealing with Jsource d_constant = inputs.D.b.y[array_ind]; - s_index.j -= inputs.I0.index; - s_index.k -= inputs.K0.index; + s_index.cell_b -= inputs.I0.index; + s_index.cell_c -= inputs.K0.index; source_value = inputs.Jsource[s_index]; cell_to_update.i = cell_b; @@ -525,8 +526,8 @@ void SimulationManager::H_source_update_steadystate( // Dealing wtih Ksource d_constant = inputs.D.b.z[array_ind]; - s_index.j -= inputs.I0.index; - s_index.k -= inputs.J0.index; + s_index.cell_b -= inputs.I0.index; + s_index.cell_c -= inputs.J0.index; source_value = inputs.Ksource[s_index]; cell_to_update.i = cell_b; From 147d8bb06a6f5f636d6abdb8a353b4f1026999bf Mon Sep 17 00:00:00 2001 From: willGraham01 <1willgraham@gmail.com> Date: Thu, 23 Feb 2023 10:05:11 +0000 Subject: [PATCH 25/25] Small doc changes that didn't get caught --- tdms/include/source.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tdms/include/source.h b/tdms/include/source.h index a86c07bb9..b554ae2d2 100644 --- a/tdms/include/source.h +++ b/tdms/include/source.h @@ -20,9 +20,9 @@ * efficient and notationally-consistent manner. */ struct SourceIndex { - int split_field_ID = 0,//!< Value in the i/x direction - cell_b = 0, //!< Value in the j/y direction - cell_c = 0; //!< Value in the k/z direction + 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 }; /** @@ -37,7 +37,7 @@ struct SourceIndex { * The Source data (real and imag) is indexed by 3 indices, accessed via * {real,imag}[cell_c][cell_b][split_field_ID]. * - * source_field_ID ranges between 0-7 inclusive. + * split_field_ID ranges between 0-7 inclusive. * TODO: Indices <-> sources need to be deduced from input file generator * functions. *