better source-file test

[stevengj]

Now that #388 is merged, we need a better test for this feature.

I would suggest doing:

1. an off-center source box. e.g. if the simulation size is 1x1, then a 0.3x0.2 source box centered at (0.1,0.2).

2. a source amplitude function that is linear (so that it is exactly reproduced by bilinear interpolation) and non-symmetrical (to catch bugs where the data is rotated/flipped). e.g. just the function amp(p) = p.x + 2*p.y

The corresponding HDF5 file can be created by just creating a NumPy array at any reasonably high resolution, say NxM=100x200 (non-square to test for bugs where dimensions get swapped) that stores the function amp(p). Some care is required to make sure that the coordinates in the array are mapped to points p: for a 0.3x0.2 box, you'll want the array coordinate a[i,j] (where i=0..N-1, j=0..M-1 in the usual Python indexing) to correspond (I think) to x = i/N - 0.15 and y = j/M - 0.1.

[ChristopherHogan]

I believe I've set up the test as described in this issue:

def amp_fun(p):
    return p.x + 2 * p.y

class TestAmpFileFunc(unittest.TestCase):
    fname = 'amp_func_file.h5'
    dset = 'amp_data'

    def create_h5data(self):
        N = 100
        M = 200
        self.amp_data = np.zeros((N, M, 1), dtype=np.complex128)

        for i in range(N):
            for j in range(M):
                v = mp.Vector3(i / N - 0.15, j / M - 0.1)
                self.amp_data[i, j] = amp_fun(v)

        with h5py.File(self.fname, 'w') as f:
            f[self.dset + '.re'] = np.real(self.amp_data)
            f[self.dset + '.im'] = np.imag(self.amp_data)

    def init_and_run(self, test_type):
        cell = mp.Vector3(1, 1)
        resolution = 10
        fcen = 0.8
        df = 0.02

        cen = mp.Vector3(0.1, 0.2)
        sz = mp.Vector3(0.3, 0.2)

        if test_type == 'file':
            sources = [mp.Source(mp.ContinuousSource(fcen, fwidth=df), component=mp.Ez, center=cen,
                                 size=sz, amp_func_file=self.fname, amp_func_dataset=self.dset)]
        elif test_type == 'func':
            sources = [mp.Source(mp.ContinuousSource(fcen, fwidth=df), component=mp.Ez, center=cen,
                                 size=sz, amp_func=amp_fun)]

        sim = mp.Simulation(cell_size=cell, resolution=resolution, sources=sources)
        sim.run(until=200)
        return sim.get_field_point(mp.Ez, mp.Vector3())

    def test_amp_file_func(self):
        self.create_h5data()
        field_point_amp_file = self.init_and_run(test_type='file')
        field_point_amp_func = self.init_and_run(test_type='func')
        self.assertAlmostEqual(field_point_amp_file, field_point_amp_func)

The test fails though

(-0.014293682535476492+0j) != (0.00037394548515602883+0j) within 7 places

Since i / N and j / M are between 0 and 1, shouldn't we multiply those by the source size? If I change the vector construction to

 v = mp.Vector3((i / N) * 0.3 - 0.15, (j / M) * 0.2 - 0.1)

the result is closer, but still not accurate to 7 places

(0.0004777429808667141+0j) != (0.00037394548515602883+0j) within 7 places

Am I setting up the test correctly, or does this look like a bug in the source file code?

[oskooi]

It doesn't seem as though the source amplitude values are the same for the file and function examples. The source amplitude function is called with a different range of (x,y) coordinates for the rectangular area (bottom-left to top-right corner) of the two cases: the file example has (-0.15,-0.1) to (0.84,0.895) while the function example has (-0.05,0.1) to (0.25,0.3). Note that the origin of the computational cell is (0,0) which ranges from (-0.5,-0.5) [bottom left] to (0.5,0.5) [top right]. The file example thus has values which defined outside of the computational cell.

[ChristopherHogan]

the function example has (-0.05,0.1) to (0.25,0.3)

If I print out all the points in the amplitude function, I get (-0.2, -0.1) to (0.2, 0.2). I guess those are relative to the center of the source, so adding the center (0.1, 0.2) gives (-0.1, 0.1) to (0.3, 0.4), which still doesn't match what you say the function example has. Am I missing something?

For the function example, amp_fun is evaluated at 20 points, 11 of which are outside of the source area. The size of the area evaluated is (0.4, 0.3) but the source size is only (0.3, 0.2). So what is the amplitude file supposed to be? I assume it's amp_fun evaluated at (100 x 200) points, but is that distributed over the size of the source (0.3, 0.2) or the size of the area evaluated by the function example (0.4, 0.3). Are the points relative to the center of the source or to the origin? Is file[0, 0] supposed to correspond to the bottom left?

[oskooi]

You are correct: current amplitude functions in Meep are defined relative to the center of the current source. It might be helpful to define the amplitude function amp_fun relative to a fixed point (i.e., the center of the source amplitude box) so that the "origin" is fixed, similar to the example in pw-source.py.