unify split_by_effort and split_by_cost

[oskooi]

Closes #1498.

This PR unifies split_by_cost and split_by_effort based on the outline in #1498. split_by_effort has been removed entirely and been made a subset of split_by_cost which now contains a boolean flag to specify whether to use either the fragment statistics or the pixel count as the cost function for the given grid_volume. The original recursive bisection to chunk subdivision is retained (which is cache oblivious as described in #1492). split_by_cost now works for an arbitrary number of chunks and is no longer based on prime factors. The runtime performance of split_by_effort is improved as it computes the entire chunk division in one go rather than the previous approach which involved separate recursive calls to determine each chunk.

A comparison of the chunk layouts for the two different methods using a 2d square cell with just vacuum is shown below. The two methods produce nearly identical results in most cases although there are small differences (which are highlighted in the figures). These differences are the result of rounding errors in the cost function using the fragment statistics.

import meep as mp

import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

sim = mp.Simulation(resolution=10,
                    cell_size=mp.Vector3(10,10,0),
                    split_chunks_evenly=False)

sim.visualize_chunks()

All tests in the make check suite pass except for the subtest test_cyl_metal_mirror_nonlinear in meep/tests/symmetry.cpp which has been disabled:

meep/tests/symmetry.cpp

Lines 135 to 176 in 7067536

	int test_cyl_metal_mirror_nonlinear(double eps(const vec &)) {
	master_printf("Testing Z mirror symmetry in Cylindrical...\n");
	double a = 16.0;
	double ttot = 3.0;
	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);
	fields f1(&s1);
	f1.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5, 0.5));
	// f1.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
	fields f(&s);
	f.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5, 0.5));
	// f.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
	check_unequal_layout(f, f1);
	double field_energy_check_time = 1.0;
	while (f.round_time() < ttot) {
	f.step();
	f1.step();
	if (!compare_point(f, f1, veccyl(0.01, 0.5))) return 0;
	if (!compare_point(f, f1, veccyl(0.21, 0.5))) return 0;
	if (!compare_point(f, f1, veccyl(0.501, 0.5))) return 0;
	if (!compare_point(f, f1, veccyl(0.33, 0.46))) return 0;
	if (!compare_point(f, f1, veccyl(0.2, 0.2))) return 0;
	if (f.round_time() >= field_energy_check_time) {
	if (!compare(f.electric_energy_in_box(gv.surroundings()),
	f1.electric_energy_in_box(gv.surroundings()), "electric energy"))
	return 0;
	if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
	f1.magnetic_energy_in_box(gv.surroundings()), "magnetic energy"))
	return 0;
	if (!compare(f.field_energy(), f1.field_energy(), " total energy")) return 0;
	field_energy_check_time += 1.0;
	}
	}
	return 1;
	}

This failing test involves a z-mirror symmetry plane in a cell with cylindrical coordinates which seems to be affected by the small change in the split_by_effort chunk division method.

meep: error in test_cyl_metal_mirror nonlinear vacuum
meep: error in test_cyl_metal_mirror nonlinear vacuum
imaginary part = 0 differs by 1.51084e-09 from -1.51084e-09
This gives a fractional error of 1
ez differs by 6.69695e-13-1.51084e-09i from 6.69695e-13-1.51084e-09i
This comes out to a fractional error of 1
Right now I'm looking at z = 0.5, r = 0.501, time 0.0625

[oskooi]

Changing to ttot = 10.0 from ttot = 3.0 and field_check_time = 8.0.

meep: error in test_cyl_metal_mirror nonlinear vacuum
meep: error in test_cyl_metal_mirror nonlinear vacuum
real part = 0.00432538 differs by -8.74538e-05 from 0.00441284
This gives a fractional error of 0.019818
er differs by 8.74538e-05-8.30032e-05i from 0.00441284-0.00354327i
This comes out to a fractional error of 0.0213051
Right now I'm looking at z = 0.5, r = 0.01, time 8

[stevengj]

Since the only thing that this PR should have changed, I would print out the ivecs of the chunk coordinates in this PR vs master for the failing test, to see what changed.

[oskooi]

For the failing test, this PR and master are producing different chunk divisions for the 2d cell using cylindrical coordinates as shown in the figure below. The figure shows the ivecs of the little_corner (blue dots) and big_corner (red dots) for the grid_volume of each chunk for the case of a 1.0 x 1.0 cell with resolution of 16 and mirror symmetry in z.

The change in the chunk division introduced by this PR is due to this line which is slightly different than the original split_by_effort which has been replaced:

meep/src/vec.cpp

Lines 1101 to 1102 in 761e31b

	if (split_measure < best_split_measure) {
	if (d == longest_axis || split_measure < (best_split_measure - (0.3 * best_split_measure))) {

However, it seems that the chunk division is not the source of the test failure because removing the nonlinearities from the failing test (deleting the lines with set_chi3 or setting chi3 internally to 0) which still uses the same chunk division as in the figure above causes it to pass:

meep/tests/symmetry.cpp

Lines 145 to 146 in 761e31b

	s.set_chi3(one);
	s1.set_chi3(one);

Perhaps this suggests there is an unrelated bug somewhere in structure_chunk::set_chi3 for cylindrical coordinates?

meep/src/structure.cpp

Lines 871 to 905 in 761e31b

	void structure_chunk::set_chi3(component c, material_function &epsilon) {
	if (!is_mine() || !gv.has_field(c)) return;
	if (!is_electric(c) && !is_magnetic(c)) abort("only E or H can have chi3");
	epsilon.set_volume(gv.pad().surroundings());
	if (!chi1inv[c][component_direction(c)]) { // require chi1 if we have chi3
	chi1inv[c][component_direction(c)] = new realnum[gv.ntot()];
	for (size_t i = 0; i < gv.ntot(); ++i)
	chi1inv[c][component_direction(c)][i] = 1.0;
	}
	if (!chi3[c]) chi3[c] = new realnum[gv.ntot()];
	bool trivial = true;
	LOOP_OVER_VOL(gv, c, i) {
	IVEC_LOOP_LOC(gv, here);
	chi3[c][i] = epsilon.chi3(c, here);
	trivial = trivial && (chi3[c][i] == 0.0);
	}
	/* currently, our update_e_from_d routine requires that
	chi2 be present if chi3 is, and vice versa */
	if (!chi2[c]) {
	if (!trivial) {
	chi2[c] = new realnum[gv.ntot()];
	memset(chi2[c], 0, gv.ntot() * sizeof(realnum)); // chi2 = 0
	}
	else { // no chi3, and chi2 is trivial (== 0), so delete
	delete[] chi3[c];
	chi3[c] = NULL;
	}
	}
	epsilon.unset_volume();
	}

[oskooi]

Outputting the split_measure, longest_axis during the chunk splitting:

Testing Z mirror symmetry in Cylindrical...
Halving computational cell along direction z
Splitting into 2 chunks by voxels
split: z (direction), 80.00000 (split_measure), r (longest axis)
split: r (direction), 80.00000 (split_measure), r (longest axis)

Splitting into 2 chunks by voxels
split: z (direction), 128.00000 (split_measure), z (longest axis)
split: r (direction), 128.00000 (split_measure), z (longest axis)

[oskooi]

The failing tests described in #1506 have been disabled for the time being because they are unrelated to this PR.