fix memory leaks in structure and fields load during checkpointing

[oskooi]

Fixes two memory leaks introduced in #1738 and reported by the address sanitizer.

cc @kkg4theweb

[oskooi]

While this PR with the latest commit (e2eae75) fixes the memory leak in the C++ unit test dump_load.cpp, it is causing the Python unit test test_dump_load.py to fail. The failing Python unit test involves time stepping the fields in a reference simulation until t=15, saving the fields at t=15, and running until t=20. The fields at some arbitrary monitor point are saved at every 1 time units and the results are compared to a second test simulation which involves loading the fields from the reference simulation at t=15. What is strange is that both simulations produce identical results for t=15 and t=16 but then the test simulation starts to diverge.

Using MPI version 3.1, 1 processes
Saving temp files to dir: /tmp/meepTwK9OZ
Running test_load_dump_fields_2d
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.000271299 s
Working in 2D dimensions.
Computational cell is 5 x 5 x 0 with resolution 50
     block, center = (0,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (10.24,10.24,10.24)
     block, center = (1,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (13,13,13)
time for set_epsilon = 0.0707238 s
-----------
run 0 finished at t = 15.0 (1500 timesteps)
creating epsilon from file "/tmp/meepTwK9OZ/test_load_dump_fields/structure.h5" (1)...
Dumped structure to file: /tmp/meepTwK9OZ/test_load_dump_fields/structure.h5 (True)
creating fields output file "/tmp/meepTwK9OZ/test_load_dump_fields/fields.h5" (1)...
Dumped fields to file: /tmp/meepTwK9OZ/test_load_dump_fields/fields.h5 (True)
run 1 finished at t = 20.0 (2000 timesteps)
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 9.4034e-05 s
Working in 2D dimensions.
Computational cell is 5 x 5 x 0 with resolution 50
time for set_epsilon = 0.0330667 s
-----------
Working in 2D dimensions.
Computational cell is 5 x 5 x 0 with resolution 50
time for set_epsilon = 0.0333062 s
-----------
reading epsilon from file "/tmp/meepTwK9OZ/test_load_dump_fields/structure.h5" (1)...
Loaded structure from file: /tmp/meepTwK9OZ/test_load_dump_fields/structure.h5 (True)
reading fields from file "/tmp/meepTwK9OZ/test_load_dump_fields/fields.h5" (1)...
Loaded fields from file: /tmp/meepTwK9OZ/test_load_dump_fields/fields.h5 (True)
on time step 1500 (time=15), 2.0487 s/step
run 0 finished at t = 20.0 (2000 timesteps)
F
======================================================================
FAIL: test_load_dump_fields_2d (__main__.TestLoadDump)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "test_dump_load.py", line 257, in test_load_dump_fields_2d
    self._load_dump_fields_2d()
  File "test_dump_load.py", line 251, in _load_dump_fields_2d
    self.assertAlmostEqual(ref_field_points[t], v)
AssertionError: 0.3447239547967911 != 0.34472334384918213 within 7 places

----------------------------------------------------------------------
Ran 1 test in 0.464s

FAILED (failures=1)

Here are the results from the two runs.

ref:, 15.0, 0.42263978719711304
ref:, 16.0, 0.3860078901052475
ref:, 17.0, 0.3447239547967911
ref:, 18.0, 0.30143244564533234
ref:, 19.0, 0.25082115083932877
ref:, 20.0, 0.19293657690286636

test:, 15.0, 0.42263978719711304
test:, 16.0, 0.3860078901052475
test:, 17.0, 0.34472334384918213 ***
test:, 18.0, 0.3033006340265274 ***
test:, 19.0, 0.2514386847615242 ***
test:, 20.0, 0.18865110725164413 ***

[oskooi]

I think I have tracked down the cause of the two failing tests. There seems to be a bug in saving/load the fields in the PML regions. This can be demonstrated by simply removing pml_layers from the Simulation object which causes both of the failing tests shown below to pass. It's unclear though why this is related to this PR.

meep/python/tests/test_dump_load.py

Lines 180 to 196 in e41a735

	# Tests dumping/loading of fields & structure in 2d.
	def _load_dump_fields_2d(self, single_parallel_file=True):
	resolution = 50
	cell = mp.Vector3(5, 5)
	sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.4), center=mp.Vector3(), component=mp.Ez)
	one_by_one = mp.Vector3(1, 1, mp.inf)
	geometry = [mp.Block(material=mp.Medium(index=3.2), center=mp.Vector3(), size=one_by_one),
	mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)]
	pml_layers = [mp.PML(0.5)]
	symmetries = [mp.Mirror(mp.Y)]
	sim1 = mp.Simulation(resolution=resolution,
	cell_size=cell,
	boundary_layers=pml_layers,
	geometry=geometry,
	symmetries=symmetries,
	sources=[sources])

meep/python/tests/test_dump_load.py

Lines 261 to 277 in e41a735

	# Tests dumping/loading of fields & structure in 3d.
	def _load_dump_fields_3d(self, single_parallel_file=True):
	resolution = 15
	cell = mp.Vector3(2.6, 2.2, 2.3)
	sources = mp.Source(src=mp.GaussianSource(1, fwidth=0.4), center=mp.Vector3(), component=mp.Hx)
	one_by_one_by_one = mp.Vector3(1, 1, 1)
	geometry = [mp.Block(material=mp.Medium(index=2.4), center=mp.Vector3(), size=one_by_one_by_one),
	mp.Block(material=mp.Medium(epsilon=7.9), center=mp.Vector3(1), size=one_by_one_by_one)]
	pml_layers = [mp.PML(0.5)]
	symmetries = [mp.Mirror(mp.Y,phase=-1),mp.Mirror(mp.Z,phase=-1)]
	sim1 = mp.Simulation(resolution=resolution,
	cell_size=cell,
	boundary_layers=pml_layers,
	geometry=geometry,
	symmetries=symmetries,
	sources=[sources])

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 50.28%. Comparing base (e28f391) to head (6349c9b).
Report is 313 commits behind head on master.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@             Coverage Diff             @@
##           master    #1872       +/-   ##
===========================================
- Coverage   73.40%   50.28%   -23.13%     
===========================================
  Files          17       17               
  Lines        4896     4914       +18     
===========================================
- Hits         3594     2471     -1123     
- Misses       1302     2443     +1141     

see 16 files with indirect coverage changes

[oskooi]

This is finally working. The bug was a bit subtle to track down and has to do with the fact that during the fields initialization of fields_chunk::alloc_f H = B in the PML region:

meep/src/fields.cpp

Lines 475 to 483 in e41a735

	/* initially, we just set H == B ... later on, we lazily allocate
	H fields if needed (if mu != 1 or in PML) in update_eh */
	component bc = direction_component(Bx, component_direction(c));
	if (!f[bc][cmp]) {
	f[bc][cmp] = new realnum[gv.ntot()];
	for (size_t i = 0; i < gv.ntot(); i++)
	f[bc][cmp][i] = 0.0;
	}
	f[c][cmp] = f[bc][cmp];

H is lazily allocated during the first timestep in fields_chunk::update_eh:

meep/src/update_eh.cpp

Lines 166 to 172 in e41a735

	// lazily allocate any E/H fields that are needed (H==B initially)
	if (i == 0 && f[ec][cmp] == f[dc][cmp] &&
	(s->chi1inv[ec][d_ec] || have_f_minus_p || dsigw != NO_DIRECTION)) {
	f[ec][cmp] = new realnum[gv.ntot()];
	memcpy(f[ec][cmp], f[dc][cmp], gv.ntot() * sizeof(realnum));
	allocated_eh = true;
	}

The problem here is that when the saved fields are being loaded from the HDF5 file, H = B (because fields::step() and thus fields_chunk::update_eh has not yet been invoked) and thus the saved H fields were never actually being loaded into a separate H fields array. The H fields in the PML region were being assigned to the saved B fields in the PML region which is why the results were slightly different from the reference simulation. The fix involved manually allocating the fields array for H in the PML region within fields::load_fields_chunk_field (rather than in fields_chunk::update_eh) so that the saved H fields from could be then be correctly loaded from the file into the H fields array.

I added a comment to explain all this in fields::load_fields_chunk_field.