Incorrect bounds checking on conductive_aux.rho_*

[prmunro]

Describe the bug
When the multilayer variable in the input file is not set to [], there are two errors, one within iteratefdtd_matrix.m and one within the TDMS executeable. The first error (iteratefdtd_matrix) is due to a check to ensure the entries in multilayer increase monotomically, however, this test was specified incorrectly. The second error assumes that resulting data structures (conductive_aux.rho_x and conductive_aux.rho_y) should be vectors, which is incorrect.

To Reproduce

1. Run gen_input_file.m in the attached, path to TDMS will need updating.
2. This will generate an error in iteratefdtd_matrix.m, which may be fixed using the following change:

line 299 of iteratefdtd_matrix.m, currently has:

elseif all(sort(multilayer)==multilayer)

whereas is should be

elseif ~all(sort(multilayer)==multilayer)

3. With the above correction made, tdms can be run, which will generate this error:

terminate called after throwing an instance of 'std::runtime_error'
  what():  Incorrect dimension conductive_aux.rho_x. Required 1D

Inspecting the mat file produced by iteratefdtd_matrix reveals:

>> dat.conductive_aux

ans = 

  struct with fields:

    rho_x: [41x61 double]
    rho_y: [41x61 double]
    rho_z: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ]

meaning that rho_x and rho_y are not vectors, which is the expected behavior. Looking at the original version of the PSTD/FDTD code, we see the following error checking:

---

  /*Get conductive_aux */

  if(mxIsStruct(prhs[input_counter])){
    num_fields = mxGetNumberOfFields(prhs[input_counter]);
    if(num_fields != 3){
      sprintf(message_buffer, "conductive_aux should have 3 members, it has %d",num_fields);
      mexfprintf(message_buffer);
    }
    for(int i=0;i<3;i++){
      element = mxGetField( (mxArray *)prhs[input_counter], 0, conductive_aux_elements[i]);
      ndims = mxGetNumberOfDimensions(element);
      if( ndims == 2 ){
	dimptr_out = mxGetDimensions(element);
	if( !(dimptr_out[0] == 1 ||dimptr_out[0] == 0) ){ 
	  //sprintf(message_buffer, "Incorrect dimension on conductive_aux.%s",conductive_aux_elements[i]);
	  //mexfprintf(message_buffer);
	}
      }
      if(!strcmp(conductive_aux_elements[i],"rho_x")){
	rho_x = mxGetPr(element);
      }
      else if(!strcmp(conductive_aux_elements[i],"rho_y")){
	rho_y = mxGetPr(element);
      }
      else if(!strcmp(conductive_aux_elements[i],"rho_z")){
	rho_z = mxGetPr(element);
      }
      else{
	sprintf(message_buffer, "element conductive_aux.%s not handled",conductive_aux_elements[i]);
	mexfprintf(message_buffer);
      }
    }
  }
  else{
    sprintf(message_buffer, "Argument %d was expected to be a structure (conductive_aux)",input_counter);
    mexfprintf(message_buffer);
  }
  input_counter++;
  /*Get conductive_aux */
************************************************************** 

which does not require elements to be vectors. I am happy to not check the dimension of these variables since they should be set correctly in iteratefdtd_matrix

input_file_ml.zip

This bug is platform independent.

[willGraham01]

Hi Peter - Sam and I are taking a look at this today.

We've managed to get TDMS to run without an error (and without editing the C++/executable code) after applying your suggested change to iteratefdtd_matrix: the following output file is produced. @samcunliffe & I of course don't know if this is the output you'd expect, so please let us know if that output is wrong.

output-on-main-with-ifdtd-patch.zip

This output the result of running the "development version" of TDMS - it includes the partial switch from the MATLAB reader to the HDF5 reader. It seems that these changes (which aren't in version 1.0.1) at least prevent the segmentation fault (likely because we're no longer relying on MATLAB-allocated memory) - but as I mentioned we don't know if it actually fixes the bug or just prevents the seg-fault and produces an incorrect output!

[samcunliffe]

Hey @prmunro,

We travelled back in time to the very first version of the code you gave us in Jan 2022, and we think the above output file is consistent with what was being done back then.

If the above file is wrong, could you give us an example pair of input/output files that you are confident in the output's correctness, please?

[prmunro]

Hi @samcunliffe , sorry for the delay in responding. I'll have a look at the development version, however, the main problem with v 1.0.1, or at least the initial problem, is that it is generating an error when dat.conductive_aux.rho_x and dat.conductive_aux.rho_y are matrices. Are you able to see in v 1.0.1 that this is occuring?

[willGraham01]

Hi @prmunro,

is that it is generating an error when dat.conductive_aux.rho_x and dat.conductive_aux.rho_y are matrices. Are you able to see in v 1.0.1 that this is occuring?

Yes we are able to see this occurring in version 1.0.1 on our machines too. We also have an idea as to why this bug/error occurs, given the fact that the "development" version is working and 1.0.1 is not:

rho_x is read into a "flattened vector", but it looks like the TDMS executable is setup so as to expect the input to have already been flattened by MATLAB. EG If rho_x is 5-by-10, tdms is expecting to read a 50-by-1 input from the file (which incidentally is why rho_x & rho_y being row vectors has not thrown the error). This is assuming that the output from the development version (attached above) is correct - otherwise, it is likely a slightly more complex error in how the MATLAB memory is being read.

[prmunro]

Hi @willGraham01 , my understanding is that TDMS should not be expecting a 50-by-1 vector. I think in the first version of the code, rho_x could be, for exampl, 5-by-10, and TDMS was also expecting that.

Please tell me if I'm missing your point above :-)

[willGraham01]

I think in the first version of the code, rho_x could be, for exampl, 5-by-10, and TDMS was also expecting that.

Hi Peter - Sam & I came to this conclusion as well; apologies my explanation above wasn't super clear.

The short version is: we (incorrectly) added a check to rho_x/y that they should be vectors, as we were reorganising the code. We didn't catch this because there is no test data on CI that tries to use these values as matrices. My hypothesis: removing this "check" (which is what the development version effectively does) once again lets TDMS run as it did before. Removing the check in version 1.0.1's code I also think would do the same thing (but haven't tested it).

The long version:

TDMS originally read the (say 5-by-10) matrix into a pointer (float* to be precise), and then manually worked out the appropriate indices to respect the "shape" that's lost when using a single-index pointer. So in the code, rho_x is only ever accessed as rho_x[array_index], where array_index is something like K_tot * k + j mapping to the original (k,j) index in the 5-by-10 structure.

Flash forward to our first steps refactoring / reorganisation of the code. rho_x is still read into a flattened vector, but it now contains a hard check for the input data to be 1D (see here for the offending function and its "error check"). This has slipped under the radar because none of the tests give rho_x as a 2D array. However, rho_x is still indexed in the same way (array_ind = K_tot * k + j, etc) - so I suspect "removing" this error check will result in the original behaviour coming back (though haven't tried this myself). This is why v1.0.1 throws the error.

Now on the "development version" of TDMS, this bug is not seen because we use the HDF5 reader functions. In particular, we replaced the individual "read vector/matrix/tensor" etc methods with a general "read input data and convert to a flattened array". As such, whatever shape rho_x comes in as (vector or matrix) no longer matters from the point of view of the reader, and thus TDMS runs and produces an output. Again, rho_x is still being indexed in the same way, it's just the reader that has changed but the output should be consistent with the original.

[prmunro]

Thanks for the above, sorry it's taking me ages to get to this.

Can you please remind me how to access the development version of TDMS?

[willGraham01]

Hi Peter,

The development version of TDMS is the latest version on main. If you open a terminal and run the following;

git clone git@github.com:UCL/TDMS.git  # Don't need these two commands if you've already 
cd TDMS                                # got the TDMS repository cloned from GitHub, just
                                       # cd into the appropriate folder
git checkout main
git pull

you'll have the development version of the code, which you can then compile with the usual instructions on the website.

You can also download the code for the development version as a zip file from github, but the clone method is recommended if you've used the GitHub clone before.

[prmunro]

Hi @willGraham01 , sorry it's taking me a long time to get to this. I just noticed that test arc_15 tests this functionality, so I'm confused as to how the test was passing. Could you please check that this test was indeed passing on the stable version of TDMS? This test can be used to check the correctness of the development version too.

By the way, do we test for the creation of the .mat input file? Just asking as when I first noticed this problem, iteratefdtd_matrix.m would generate an error when multilayer was non-empty. So test arc_15 should also generate an error...unless I'm missing something!

Peter

[willGraham01]

Hi @prmunro,

I just noticed that test arc_15 tests this functionality, so I'm confused as to how the test was passing.

It looks like arc_15 (and arc_14 for that matter, if it exists) aren't tracked in the GitHub repository, so aren't being run when new versions are uploaded. See tests run on development version and tests run on stable version. This test data is also not on Zenodo. This would explain why we didn't catch this sooner.

By the way, do we test for the creation of the .mat input file?

We do when running the system tests (those linked to above). The first stage of these tests is to take the pstd/fdtd_input_file_XX.m and run it through the run_bscan.m script, to generate the input .mat file. Then that file gets passed to TDMS and the output is checked against the Zenodo data (expected output).

If you have arc_15 to hand (the input file, the particular version of run_bscan that it gets passed through, the expected .mat file that iteratefdtd_matrix should produce, and the TDMS output file) then I can try out running these on the stable and development versions of TDMS.

[prmunro]

Hi @willGraham01 , thanks for this. Here is a onedrive link to arc_14 and arc_15. I hope you can follow the logic of run_pstd_bscan.m. It generates the input .mat files (if recalc is set to 1) and then performs a test against analytically calculated fields. This test code uses an input .m file that is not compatible with TDMS, i.e., it's using my original file format. I know you have a way to deal with this, but let me know if that's a problem.

https://liveuclac-my.sharepoint.com/:f:/g/personal/rmaprtm_ucl_ac_uk/EoNTGsRrw4xHmYPA2KS7G-ABcnZzDM-RyCBplq5DTwESrA?email=william.graham%40ucl.ac.uk&e=0pZa7m

Cheers - Peter

[willGraham01]

Hi @prmunro -

I've just run these examples through the development version of TDMS. A mixture of good news and bad news with these tests as detailed below. (Applies to both arc_14 and arc_15).

For both tdms version 1.0.1 & the development version, the only way to generate the in/pstd_fs.mat files correctly is to set use_pstd = 0 and compactsource = 0. Doing so, then running the bscan script will produce an input file that is identical to the one provided in the zip files above. However passing this input into tdms results in use of the fdtd solver (since use_pstd = 0) which produces an incorrect out/pstd_fs.mat. Interestingly enough though, if I generate in/pstd_fs.mat as described above, then manually change use_pstd to 1 in the generated in/pstd_fs.mat, then run tdms, the resulting out/pstd_fs.mat seems to be within good (numerical) agreement of the output provided.

Attempting to generate in/pstd_fs.mat as I would expect (with use_pstd = 1, which requires compactsource = 1) results in in/pstd_fs.mat being incorrect (again compared to the provided file in the zip folder). The Ksource term in particular is generated incorrectly in this case, which is of course the one affected by the compactsource condition in this case.

So a couple of things here:

• Are arc_14 and arc_15 meant to be PSTD solver simulations?
  • If so, then it looks like the development version of tdms is producing the correct result given the input file and reference outputs you've provided.
  • I haven't checked this for the 1.0.1 version, but I imagine they will fail for the reason you described in the original bug.
• If this is meant to be a PSTD solver simulation, we have to be using the compact source condition. Is it possible that - because these input files are from an old version of tdms - that the zip file data is affected by the dz/2 offset bug (link)?