WIP: allow add_srcdata to place sources on non-owned points

[stevengj]

For ongoing work in parallelizing the placement of adjoint sources for DFT fields using add_srcdata, we discovered that we need to place sources on non-owned points (since these are used for voxel-center averaging from the Yee grid). After the sources are added, this unfortunately requires a "cleanup" step, implemented in this PR, which transfers sources from non-owned points to the processors where they are actually located.

[mawc2019]

There are some compilation errors. In fix_boundary_sources.cpp, after adding

#ifdef HAVE_MPI
static MPI_Comm mycomm = MPI_COMM_WORLD;
#endif

below this line and adding #include <algorithm> below this line, the compilation errors are gone.

In order to test the code, I simply changed this line in simulation.py as

self.fields.add_srcdata(src.srcdata, src.src.swigobj, src.num_pts, src.amp_arr,needs_boundary_fix=False)

but NotImplementedError appeared when running my test example. After removing the default value of needs_boundary_fix at this line in meep.hpp, NotImplementedError is gone.

[smartalecH]

@mawc2019 any updates on this?

[mawc2019]

@smartalecH , Steven pointed out some revisions:

1. Add id = t.id; in src_time(const src_time &t) below this line of meep.hpp;
2. To allocate an id for the source time, register the source time in _create_time_profile of objective.py as

    def _create_time_profile(self, fwidth_frac=0.1):
        srct = mp.GaussianSource(
            np.mean(self._frequencies),
            fwidth=fwidth_frac * np.mean(self._frequencies))
        self.sim.fields.register_src_time(srct.swigobj)
        return srct

But the id registration of source time still fails in the adjoint run of my modified FourierFields, which places the adjoint sources in a way similar to that in Near2FarFields. I also tried to register the source time of the adjoint source by adding self.fields.register_src_time(self.sources.swigobj) in change_sources of simulation.py, but got 'NoneType' object has no attribute 'register_src_time'.

I plan to upload my modified FourierFields at first and then upload a revised version of this fix_boundary_sources branch, so that the problem may be clearer. However, I cannot upload my modified FourierFields at the moment, because the adjoint and finite-difference gradients are inconsistent when my modified FourierFields is based on the current master branch. In contrast, they were consistent when my modified FourierFields was based on the old master branch a few month ago. I am looking for the reason.

[smartalecH]

Thanks for the update. Maybe submit a PR for what you have, and then we can all look at it.

The adjoint recombination step was recently overhauled. In theory it shouldn't cause any problems, but you never know.

[mawc2019]

With the aforementioned corrections added and the new FourierFields merged, the revised version of this WIP is here.

When the FourierField adjoint solver needs the boundary fix (which may happen when yee_grid=False), the error message bug: fix_boundary_sources called for non-registered source is reported during the adjoint run. An example is as follows.

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
Si, SiO2 = mp.Medium(index=3.4), mp.Medium(index=1.44)

resolution,design_resolution = 20,20
Lx,Ly,pml_size = 2,4,1
cell_size = mp.Vector3(Lx,Ly+2*pml_size)
pml_layers = [mp.PML(pml_size, direction=mp.Y)]
design_width,design_height = Lx,0.4

fcen,fwidth = 1/1.55,0.01/1.55
source_center,source_size = mp.Vector3(0,-Ly/4),mp.Vector3(Lx,0)

src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,component=mp.Ez,center=source_center,size=source_size)]
Nx,Ny = int(round(design_resolution*design_width)),int(round(design_resolution*design_height))

design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_width, design_height)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)] # design region
sim = mp.Simulation(
    cell_size=cell_size,boundary_layers=pml_layers,geometry=geometry,
    sources=source,default_material=SiO2,resolution=resolution)

nf = mpa.FourierFields(sim,mp.Volume(center=mp.Vector3(0,design_height+0.2),size=mp.Vector3(Lx,0)),mp.Ez,yee_grid=False)
ob_list = [nf]

def J(near_data):
    return npa.sum(npa.abs(near_data)**2)

opt = mpa.OptimizationProblem(
    simulation = sim,objective_functions = J,objective_arguments = ob_list,design_regions = [design_region],
    fcen=fcen,df = 0,nf = 1,decay_by = 1e-5)

x0 = 0.5*np.ones(Nx*Ny)
opt.update_design([x0])
f0, dJ_deps = opt()

It seems that register_src_time in sources.cpp does not fulfill its function for the FourierFields adjoint source, which is placed in a way similar to that in Near2FarFields.

But the id registration of source time still fails in the adjoint run of my modified FourierFields, which places the adjoint sources in a way similar to that in Near2FarFields. I also tried to register the source time of the adjoint source by adding self.fields.register_src_time(self.sources.swigobj) in change_sources of simulation.py, but got 'NoneType' object has no attribute 'register_src_time'.

[stevengj]

Closed favor of #1959.