negative array indices in step_update_EDHB

[jianguan210]

When the test cases "convergence_cyl_waveguide" and "symmetry" were ran with ASan, it reports the heap-buffer-overflow meep/src/step_generic_stride1.cpp:539:29 and heap-buffer-overflow meep/src/step_generic_stride1.cpp:551:46

My study shows that the index values i - s1 and i - s2 will hit negative values in LOOP_OVER_VOL_OWNED. When I used if statements (i >=s1) and (i >= s2) to conditionally access the arrays, test was passed without warnings.
Since this loop macro is used in many places, could you please take a look and find a generic solution to fix this issue?

[oskooi]

To provide some additional information, this issue is related to a heap-buffer-overflow error due to accessing out-of-bound array elements (in particular negative array indices) in step_update_EDHB of step_generic_stride1.cpp which were traced to the following two sets of lines in step_generic.cpp.

1. linear, 2x2 off-diagonal ε (tests/convergence_cyl_waveguide.cpp):

meep/src/step_generic.cpp

Line 536 in 2dcfa9f

f[i] = (gs * us + OFFDIAG(u1, g1, s1));

In this case, OFFDIAG is a function macro defined as:

meep/src/step_generic.cpp

Lines 376 to 378 in 2dcfa9f

	// stable averaging of offdiagonal components
	#define OFFDIAG(u, g, sx) \
	(0.25 * ((g[i] + g[i - sx]) * u[i] + (g[i + s] + g[(i + s) - sx]) * u[i + s]))

In this case, the problem is with the g array index i - sx which is equivalent to i - s1.

2. nonlinear, diagonal ε (tests/symmetry.cpp):

meep/src/step_generic.cpp

Lines 547 to 548 in 2dcfa9f

	realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)];
	realnum g2s = g2[i] + g2[i + s] + g2[i - s2] + g2[i + (s - s2)];

The similar problem here are the g1 and g2 array indices i - s1 and i - s2, respectively.

The errors can be removed simply by enforcing the rule that these lines are invoked only for non-negative indices, that is when i >= s1 and i >= s2.

[stevengj]

Weird! At first glance, the field components that we average for off-diagοnal ε

should always be in the chunk volume if we are averaging to an "owned" point (via LOOP_OVER_VOL_OWNED):

[stevengj]

It might be helpful to print out the coordinates of the problematic point, i.e. add something like

if (i - s1 < 0) {
    IVEC_LOOP_ILOC(gv, iloc);
   // print out iloc
}

[oskooi]

Looking into this further, it seems the error is not related to negative array indices which must therefore mean it is due to the array indices being exceeded in the other (i.e., positive) direction.

Some more clues: the failing subtest in tests/symmetry.cpp (which according to git blame was added many years ago) also involves cylindrical coordinates (similar to the other failing test in tests/convergence_cyl_waveguide.cpp):

meep/tests/symmetry.cpp

Lines 984 to 985 in f2bae0c

	if ((count_processors() == 1) && (!test_cyl_metal_mirror_nonlinear(one)))
	abort("error in test_cyl_metal_mirror nonlinear vacuum\n");

The error is related to the chi3 nonlinearities since removing it (lines 145 and 146) from this subtest caused the error to go away:

meep/tests/symmetry.cpp

Lines 140 to 146 in f2bae0c

	const grid_volume gv = volcyl(1.0, 1.0, a);
	the_center = gv.center();
	const symmetry S = mirror(Z, gv);
	structure s(gv, eps, no_pml(), S);
	structure s1(gv, eps);
	s.set_chi3(one);
	s1.set_chi3(one);

(The error is not related to the z-mirror symmetry.)

[oskooi]

The failing subtest in tests/convergence_cyl_waveguide.cpp involves subpixel smoothing test_convergence_with_averaging which produces non-zero off-diagonal elements for the ε tensor (whereas the no smoothing case test_convergence_without_averaging runs fine):

meep/tests/convergence_cyl_waveguide.cpp

Lines 175 to 176 in f2bae0c

	test_convergence_without_averaging();
	test_convergence_with_averaging();

Following the suggestion of using an if-statement inside step_generic_stride1.cpp to print out the coordinates of the ivec at which the array index is negative reveals: i = 0, s1 = 2 (i.e., i - s1 < 0)at the location (r,z) = (0,0).

Note that this particular test setup uses a Bloch-periodic boundary condition k=(k_r,k_z)=(0.1,0) in cylindrical coordinates. Removing this boundary condition did not affect the negative array index error.