The makefiles in iterater/ should be removed in favour of a more generic build system (likely CMake). A pair of makefiles suitable for building on an M1 mac are here.

I noticed that there was a small difference between calculated
fields between the cases where eyi is set to 0 and where eyi is set to
[]. Although this was probably less than numerical precision, this
should be checked.

Currently there is only one piece of example/demonstration code, examples/arc-01.
It would be beneficial for both documentation and testing purposes if we had additional examples to hand, and additionally emphasised calls to the TDMS executable rather than the supporting MATLAB scripts:

• Suppressing the MATLAB GUI from opening when setting up the input files, and using a shell script to avoid hiding the call to tdms within MATLAB's eval()?
• Save visualisations of the output to files, rather than opening them as figures in MATLAB

@prmunro very sensibly followed the instructions and:

• the build dir didn't get populated (why?)
• cmake .. -DCMAKE_BUILT_TYPE=Test didn't create the tdms_test file (nor did ctest work).

Very embarrassing. We should double-check the instructions vs. GitHub actions are doing the same thing.

I noticed that parallelism wasn't being properly exploited in the
case of 2D simulations as the for loop executed in parallel had only
one iteration. I have hard-coded a solution to this which needs to be
made more general. In particular, I moved the pragma omp for statement
for the Exz and Hyz update equations to a different for loop in the
nest of for loops.

There should be a better way of handling the case of when
fieldsample=[], at present I just set the number elements in the
locally stored versions of the members of fieldsample to 0.

Currently there is only system tests, it would be helpful for development if there were some fast unit tests using e.g. catch2. This will depend on some refactoring (splitting functions)

The miscellaneous matlab scripts in the root directory should be moved/removed and maybe a more conventional project layout be adopted.

$ tree -L 1 
├── CMakeLists.txt
├── CONTRIBUTING.md
├── LICENSE
├── README.md
├── doc
├── example
├── include
├── src
└── tests

Thoughts much appreciated on this one!

Current suggestion:

TDMS: Time Domain Maxwell Solver

Introduce an option to the input file to choose whether to write of the Ex_out, Ey_out, etc., volume fields. These are the non-interpolated versions of the fields and are generally not used. Should be able to choose which field components will be written.

fieldsample_i and similar should really be integer arrays, however,
they are double at the moment

There should be an error message when a user specifies a particular
efield and hfield file when they also pass an illfile

Hi @samcunliffe and @willGraham01, in the process of preparing a test case for the steady state option, I found an error in my version (and therefore the TDSM version) of the code. What’s the best was to make this change on the TDSM repository. It’s a relatively simple update

The code appears in the vicinity of line 5900 in the old version of the code as follows, though let me know if you want me to make this change directly in the code. This code is now in the vicinity of line 5762 in the TDMS repo, and has been changed. Although the changes in the old code will transfer easily to the new code.

Incorrect version

     if(sourcemode==sm_steadystate){
	extractPhasorENorm(&E_norm_an, fte, tind, *omega_an, *dt, Nsteps);
	extractPhasorHNorm(&H_norm_an, fth, tind, *omega_an, *dt, Nsteps);
	for(int ifx=0;ifx<N_f_ex_vec;ifx++){
	  extractPhasorENorm(&E_norm[ifx], fte, tind, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
	  extractPhasorHNorm(&H_norm[ifx], fth, tind, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
	}

Correct version

      if(sourcemode==sm_steadystate){
	if( (tind % Nsteps)==0 ){
	  E_norm_an = 0.;
	  H_norm_an = 0.;
	}
	extractPhasorENorm(&E_norm_an, fte, tind % Nsteps, *omega_an, *dt, Nsteps);
	extractPhasorHNorm(&H_norm_an, fth, tind % Nsteps, *omega_an, *dt, Nsteps);
	//fprintf(stderr,"Enorm = %e+i%e, fte=%e, tinf=%d, %d\n", real(E_norm_an), imag(E_norm_an), fte, tind % Nsteps, Nsteps);
	for(int ifx=0;ifx<N_f_ex_vec;ifx++){
	  extractPhasorENorm(&E_norm[ifx], fte, tind % Nsteps, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
	  extractPhasorHNorm(&H_norm[ifx], fth, tind % Nsteps, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
	}

Compiling with -Wall throws a warning

.../tdms_iterator.cpp:429:29: warning: variable 'mx_fieldsample_x' set but not used [-Wunused-but-set-variable]
  mxArray *mx_fieldsample, *mx_fieldsample_x, *mx_fieldsample_y, *mx_fieldsample_z;

as the place where it may be used is commented out

Can the commented out bock and associated allocations be removed, or would it be better as an input file configurable option?
@prmunro's thoughts would be much appreciated!

Range checking should be implemented on campssample.vertices and
campssample.components

I changed exi and eyi to have size Nx x Ny x Nt, I am not sure why
I had an additional dimension of size 1. Should check this. Should
also ensure some error checking about this.

Windowing the FFT in the detsensefun feature

[Placeholder issue, we might want to reject]

We are currently building windows and linux, worth to add macos? Mention this because @samcunliffe had some strange issues compiling. If this is a workflow we expect of future users/students, we probably want to test it.

On the other hand, perhaps this is a waste of time.

masterdoc.pdf was generated with latex rather than pdflatex. The tex file needs to be updated to be compilable, and added to the repository. It is also lacking information about some features in the code.

Add a OSS compliant licence, which is compatible with the FFTW dependency

Introducing a feature to extract complex amplitudes from specific
vertices, also allowing for particular field components to be
specified.

mexFunction in iterator.cpp is the best part of 7k lines. This should be broken out into smaller functions with structs/classes for the commonly grouped function arguments.

• Convert tabs -> spaces
• Smaller function
• Common classes
• Extract methods
• Add docstrings
• Create generic field classes
• #132

Update the README to include the number of omp threads that the code will use

Currently there are only PSTD tests, at least one FDTD test should also be included, and others.

• FDTD
• steadystate = true
• exdetIntegral = true
• exi is present

For consistency I need to implement a version of
extractPhasorsVertices that does not use interpolation

I should make the extractPhasorsVertices functionality available
if using steady state source condition

The MATLAB dependency from TDMS can be removed, which will (a) simplify the build, (b) allow running on machines without a MATLAB licence and (c) allow static linking such that no runtime paths need to be set and a compiled binary distributed.

Decisions

• Tensor library vs. nested std::vector vs. raw arrays
• #181

Dependencies

• #63

I picked up a bug introduced during recent refactoring work, whereby there are now two variables storing dt, in particular:

params.dt
dt

both are used. This becomes a problem since in some instances dt can be changed within the code. After refactoring, params.dt was not updated, leading to incorrect behaviour. I will add this as an issue. For example, dt and params.dt are assigned:

  if(mxIsDouble(prhs[input_counter])){
    dt = mxGetPr( (mxArray *)prhs[input_counter]);
    params.dt = *dt;
    input_counter++;
  }
  else{
    throw runtime_error("expected dt to be a double");
  }

Then, in steady state mode, dt is in general updated:

  if(sourcemode==sm_steadystate){
    dt_old = dt[0];
    Nsteps_tmp = ceil(2.*dcpi/omega_an[0]/dt[0]*3);
    dt[0] = 2.*dcpi/omega_an[0]*3/Nsteps_tmp;
    params.dt = dt[0];
  }

This is where the but was introduced. Later we see that dt and params.dt are still used, e.g.,

   if( exphasorssurface || exphasorsvolume || exdetintegral || (nvertices > 0) ){
      if(sourcemode==sm_steadystate){
	if( (tind % Nsteps)==0 ){
	  E.angular_norm = 0.0;
	  H.angular_norm = 0.0;
	}
	
        E.add_to_angular_norm(fte, tind % Nsteps, Nsteps, params);
        H.add_to_angular_norm(fth, tind % Nsteps, Nsteps, params);

        for(int ifx=0;ifx<N_f_ex_vec;ifx++){
          extractPhasorENorm(&E_norm[ifx], fte, tind % Nsteps, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
          extractPhasorHNorm(&H_norm[ifx], fth, tind % Nsteps, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);
        }
      }

In particular,

E.add_to_angular_norm(fte, tind % Nsteps, Nsteps, params);

uses dt stored within params, yet

extractPhasorENorm(&E_norm[ifx], fte, tind % Nsteps, f_ex_vec[ifx]*2*dcpi, *dt, Nsteps);

uses dt. It would be good to remove references to dt.

Odd behaviour in interpolate.cpp [Related to #44]

Several functions in interpolate.cpp seem to have redundant for loops, and mismatches between the function doc-string and what the function actually does.
Looking at interpolateFieldCentralE_TE for example:

• An error is thrown unless k_u $\leq 0\leq$ k_l, which is in contrast to what the docstring claims the acceptable values of k_l ( $\geq 2$ ) and k_u ( $\leq$ K $-2$) are.
• The variables i_u, j_u, k_u have intent[out], but are never assigned to, not pointers, and not returned.
• for(k=k_l; k<=k_u; k++) only executes when k_u $\leq 0\leq$ k_l due to the previous error checks, and so k can only take the value 0.

The same things occur in interpolateFieldCentralE_TM, interpolateFieldCentralH_TE, and interpolateFieldCentralH_TM.

If the intended behaviour is what's listed in the doc-strings; the code needs to be changed and we need to check whether the result of any demonstration code or tests depend on then also changes, due to dependence on these functions.
If what's currently in the code is the intended behaviour, it looks like it can be made more readable and the doc-strings updated to match.

Error checking where the contents of illfile are loaded, in
particular, i don't currently check that the loaded file has the
correct element names, only the number of elements

Overhaul interpolation from the split grid to the central grid.

1. Different scheme required for each field component.
2. Is implemented for the electric field, although interpolateTimeDomainFieldCentralEBandLimited was hacked to work for the 2D case, this should be checked.
3. More significant work required for the magnetic field
4. Enable interpolation to be performed right up to the edge of the computational domain.

Dependencies:

• Theoretical outline from @prmunro

Tasks specified in PSTD_interpolation_scheme.pdf:

• Write (time-domain) interpolation methods for E-field components. Should use BLi, and if not possible, fall back on cubic interpolation. Closed by #62.
• Write (time-domain) interpolation methods for H-field components using the 2D-method interpolation method. Interpolation in each dimension should again use BLi if possible, and failing that fall back on cubic. Closed by #62.
• #112
• Unit tests for cubic interpolation methods. Closed by #62.
• #111
• #110
• #109
• #113
• #114
• #115

Add an option to the input file, which requires implementation in the C code also, to control which field components are extracted when wanting to extract an entire volume (exphasorsvolume).

The demonstration script run_pstd_bscan.m throws a MATLAB error when run.

MATLAB attempts to save the generated input/output files (that are to be passed to tdms) to the sub-directories gridfiles/, in/, and out/, all of which don't exist after cloning the repository (and successfully installing tdms) which causes the save() function to error.

MATLAB needs to check whether these directories need to be created, then call to mkdir if necessary, prior to attempting to generate the input and output files.

Nice to have something like cppcheck (big, slow?), clang-tidy, or cppclean as a precommit hook.

Need to check freeCastMatlab4DArray and CastMatlab4DArray

The ../../../FDTD path parts in the following lines

https://github.com/UCL/isotropic-finite-difference-time-domain/blob/b6136cb6f16b270425b78b530545bf8e60df498e/iterater/MakefileOMP_PSTD_UCL2#L76

https://github.com/UCL/isotropic-finite-difference-time-domain/blob/b6136cb6f16b270425b78b530545bf8e60df498e/iterater/MakefileOMP_FDTD_UCL2#L19

https://github.com/UCL/isotropic-finite-difference-time-domain/blob/b6136cb6f16b270425b78b530545bf8e60df498e/iterater/MakefileOMP_FDTD_UCL2#L73-L74

should I think just be ..

The following lines also refer to directories outside of the project root but are inside a conditional which is not entered when running the file as is

https://github.com/UCL/isotropic-finite-difference-time-domain/blob/b6136cb6f16b270425b78b530545bf8e60df498e/iterater/MakefileOMP_PSTD_UCL2#L39-L40
https://github.com/UCL/isotropic-finite-difference-time-domain/blob/b6136cb6f16b270425b78b530545bf8e60df498e/iterater/MakefileOMP_FDTD_UCL2#L37-L38

For better future-proof development of this package, I think we should

1. clear up all warnings. In particular, we should
   1. i. build the package in CI with -Wall
   2. ii. once all warnings are solved, build in CI with -Werror, to make sure no more warnings are introduced. To be clear: I'm not suggesting to use -Wall -Werror by default for all users, that's a disruptive setup
2. use a linter, like clang-tidy in CI:
   1. i. I think we can start with a job which runs clang-tidy and reports annotations if there are warnings in the changed lines
   2. ii. Ideally, clang-tidy should always come back clean, and if it doesn't the job should fail. But this can be enabled at a later stage, when all problems are resolved

Note that point 1.ii depends on 1.i, and point 2.ii depends on 2.i, but point 2.i is independent from whole of point 1 (although warnings reported by clang-tidy will usually be a superset of those reported by -Wall).

https://github.com/cpp-linter/cpp-linter-action looks a decent linter action, for point 2. I don't think 2 (and in particular 2.ii) is a high priority (also because clang-tidy can be particularly over-zealous), but in the long term it may be helpful.

Adding a set of contributing guidelines (as CONTRIBUTING.md) in the root directory covering: main branch protection, PRs enforced, work on issues, max size of PRs etc.

[Low priority]

Now we have spdlog we should use it. (A few cerrs should be spdlog::error and even spdlog::info.)

There is currently no error checking at all on the structure
fieldsample

What kind of checks do there need to be?

I have introduced a way to export time domain fields. This is
purely experimental and needs to be formalised in every aspect.

Naively running Doxygen seems to work, so propose to put an action in that runs it and publishes to gh-pages.

apologies – don't really have time for a proper review today. but in general more tests = better so generally happy. only thought is to reduce the system size so we can be more comprehensive i.e. have both 3D and 2D examples, for instance(?)

Originally posted by @t-young31 in #98 (comment)

The {x_grid_labels_out, y_grid_labels_out, z_grid_labels_out} don't seem to have a consistent datatype

I wonder if this is correct? cc: @prmunro

Originally posted by @t-young31 in #85 (comment)

We build and run tests on all three major OS. Perhaps separate out the system tests. And perhaps test parallelisation.

The Cartesian coordinates of the vertices should be written to the
output file

TDMS has only been built on MacOS and linux thus far. We need to check this is also possible on Windows

Perhaps spdlog? Suggestions welcome!

Some utility functions in the matlab/ directory are currently unused. Unit tests should be added to test the functionality. Specifically,

• fdtdduration.m
• minsteps_fdtd.m
• minsteps_pstd.m
• fdtdts.m

Dependencies

• #37

Starting point

(from @prmunro)

%exports the constants used internally in the PSTD/FDTD
%code when using pulsed illumination
[to hwhm] = fdtdduration('pstd_input_file_2D_01.m');

%calculates the minimum number of time steps required for light to
%make one double pass of the computational domain. If using the
%FDTD algorithm, the script fdtdminsteps.m should be used
Nt = minsteps_pstd('pstd_input_file_2D_01.m');

%the maximum allowable timestep
dt_upper = fdtdts('pstd_input_file_2D_01.m');