bisection search for faster splitting

[stevengj]

This updates the split_by_cost function to do a bisection search along each axis to find a balanced splitting, rather than searching every possible splitting, which should speed things up for large simulations.

(In most cases the left and right costs should be monotonic functions of the splitting point, in which cases a bisection search should be equivalent to checking every possible point.)

split_into_three should have similar improvements, but that isn't touched by this PR.

[oskooi]

This is failing python/tests/ring-cyl.py because it contains split_chunks_evenly=False. Since the fragment statistics (#681) does not support cylindrical coordinates, the fix is simple, e.g.:

structure.cpp:187-189

change

  if (meep_geom::fragment_stats::resolution == 0 ||
      meep_geom::fragment_stats::has_non_medium_material() ||
      meep_geom::fragment_stats::split_chunks_evenly) {

to

  if (meep_geom::fragment_stats::resolution == 0 ||
      meep_geom::fragment_stats::has_non_medium_material() ||
      meep_geom::fragment_stats::split_chunks_evenly ||
      gv.dim == Dcyl) {

[oskooi]

python/tests/chunks.py (the only other test involving split_chunks_evenly=False) seems to hang indefinitely when run with 4, 6, or 8 processors whereas with 2 processors finishes in about a second:

$ mpirun -n 2 python/tests/chunks.py
Using MPI version 3.0, 2 processes
-----------
Initializing structure...
Splitting into 2 chunks by cost
time for choose_chunkdivision = 0.000123024 s

[stevengj]

The ring-cyl.py test was working before, so I don't see why it should have suddenly stopped working? It looks like @ChristopherHogan added cylindrical support in #697.

I noticed a problem where my bisection search might not converge and hopefully fixed it.

[oskooi]

Output for python/tests/chunks.py from sim.visualize_chunk() using 8 processors/chunks for this PR:

[oskooi]

Output from master for 8 chunks/processors:

[stevengj]

Good, so it is coming up with the same (or almost the same…could differ by ±1 pixel) chunking.

[stevengj]

Independent of this PR, it would be nice to understand why some of the corner chunks are a different process (color) than the surrounding chunks — to minimize communication costs, it seems like we would want to assign each process to a contiguous set of chunks where possible.

[oskooi]

The following are the chunk layouts for a tall, skinny, 2d cell with an expensive flux region at one end for the split_by_cost algorithm from master and this PR for a range of processors/chunks: 4, 6, 8, 10, and 12. The results are similar but not identical. In particular, note that for the case of 4 and 6 processors/chunks, the bisection search of this PR produces 4/6 equally-sized chunks (which is correct) but that the linear search of master produces unequal-sized chunks.

import meep as mp
import matplotlib.pyplot as plt

sx = 30
sy = 5
fcen = 1.0
df = 0.1
nfreq = 500

src = mp.Source(mp.GaussianSource(fcen,fwidth=df),
                component=mp.Ez,
                center=mp.Vector3())

sim = mp.Simulation(cell_size=mp.Vector3(sx,sy),
                    sources=[src],
                    resolution=20,
                    verbose=True,
                    k_point=mp.Vector3(0.53,0.14,0),
                    split_chunks_evenly=False)

flux = sim.add_flux(fcen,
                    df,
                    nfreq,
                    mp.FluxRegion(center=mp.Vector3(14),size=mp.Vector3(y=sy)))

sim.init_sim()
sim.visualize_chunks()
if mp.am_master():
    plt.savefig('imbalanced_2d_chunk_layout.png')

master

4 processors/chunks

6 processors/chunks

8 processors/chunks

10 processors/chunks

12 processors/chunks

this PR

4 processors/chunks

6 processors/chunks

8 processors/chunks

10 processors/chunks

12 processors/chunks