WIP: Adjoint poynting flux

[smartalecH]

Closes #1307 and #1299.

[smartalecH]

I'm slightly confused by the default behavior of the get_dft_array() routine. Note that the C implementation of get_dft_array() contains two basic steps: grab the DFT fields themselves (process_dft_component()) and format the data array with collapse_array().

meep/src/dft.cpp

Lines 1082 to 1090 in d29f211

	cdouble *fields::get_dft_array(dft_fields fdft, component c, int num_freq, int *rank,
	size_t dims[3]) {
	dft_chunk *chunklists[1];
	chunklists[0] = fdft.chunks;
	cdouble *array;
	direction dirs[3];
	process_dft_component(chunklists, 1, num_freq, c, 0, &array, rank, dims, dirs);
	return collapse_array(array, rank, dims, dirs, fdft.where);
	}

When the user specifies a volume other than a plane/volume (e.g. point or line) and the user is pulling the interpolated fields from the center of each voxel, there's some odd behavior.

For a simple point, process_dft_component() returns a 2x2 array. This array is then passed to collapse_array(), which seems to grab the first point and throw the rest away.

In the case of a line monitor across 10 pixels, process_dft_component() returns a 10x2 array. Again, collapse_array() grabs the first column and throws the rest away.

Is this the expected behavior?

This makes it tricky when trying to use the indices of the DFT array data values for a broadband adjoint source (for both Poynting fluxes and basic Fourier fields).

[smartalecH]

There are two options for placing the adjoint source in this case:

1. In the C++ code, uncollapse the jacobian passed from python, combine it with all the uncollapsed dft chunk data, and add the result is an index source.
2. In the python code, use the collapsed dft fields and the collapsed Jacobian to add a source using the amp_data functionality for any arbitrary source. The python source will automatically distribute the fields onto the correct yee grid points (assuming the interpolation is taken into account).

Option 1 is rather difficult -- I don't see a reliable way to uncollapse the Jacobian and map it to the dft fields (which are just stored as sequential indices in the C++ code). Option 2 requires some slight modifications to the amp_data routine, specifically with how the array is interpolated onto the grid.

I implemented option 2 for the broadband Poynting flux (10 frequency points) and got the following result:

The result isn't perfect (1%-10% error) but it's a good start. I've incorporated the interpolation weights into the adjoint scale factor (as we've discussed numerous times) but they only slightly improve the accuracy. Perhaps I still need to take into account the interpolation of the final source object and amp_array routine?