From ee7e823b9d0ba5e9fb96450f5b9b2eba88aa094f Mon Sep 17 00:00:00 2001 From: Ian Williamson Date: Tue, 28 Sep 2021 15:32:13 -0700 Subject: [PATCH] Original fix --- python/adjoint/objective.py | 30 +++++++++++++++++++----------- python/tests/test_adjoint_jax.py | 19 ++++++++++--------- src/mpb.cpp | 5 +---- 3 files changed, 30 insertions(+), 24 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index d653c10e0..10f37022e 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -92,6 +92,10 @@ def _adj_src_scale(self, include_resolution=True): else: # multi frequency simulations scale = dV * iomega / adj_src_phase + # compensate for the fact that real fields take the real part of the current, + # which halves the Fourier amplitude at the positive frequency (Re[J] = (J + J*)/2) + if self.sim.using_real_fields(): + scale *= 2 return scale def _create_time_profile(self, fwidth_frac=0.1): @@ -110,10 +114,6 @@ def _create_time_profile(self, fwidth_frac=0.1): ) -_KPOINT_POS_DIR = mp.Vector3(0.0, 0.0, 0.0) -_KPOINT_NEG_DIR = mp.Vector3(-0.0, -0.0, -0.0) - - class EigenmodeCoefficient(ObjectiveQuantity): """A frequency-dependent eigenmode coefficient. Attributes: @@ -134,6 +134,7 @@ def __init__(self, forward=True, kpoint_func=None, kpoint_func_overlap_idx=0, + decimation_factor=0, **kwargs): super().__init__(sim) if kpoint_func_overlap_idx not in [0, 1]: @@ -148,6 +149,7 @@ def __init__(self, self.eigenmode_kwargs = kwargs self._monitor = None self._cscale = None + self.decimation_factor = decimation_factor def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) @@ -155,6 +157,7 @@ def register_monitors(self, frequencies): frequencies, mp.ModeRegion(center=self.volume.center, size=self.volume.size), yee_grid=True, + decimation_factor=self.decimation_factor, ) return self._monitor @@ -169,7 +172,8 @@ def place_adjoint_source(self, dJ): if self.kpoint_func: eig_kpoint = -1 * self.kpoint_func(time_src.frequency, self.mode) else: - eig_kpoint = _KPOINT_NEG_DIR if self.forward else _KPOINT_POS_DIR + direction = mp.Vector3(*[float(v == 0) for v in self.volume.size]) + eig_kpoint = -1 * direction if self.forward else direction if self._frequencies.size == 1: amp = da_dE * dJ * scale @@ -201,7 +205,8 @@ def __call__(self): kpoint_func = self.kpoint_func overlap_idx = self.kpoint_func_overlap_idx else: - kpoint_func = lambda *not_used: _KPOINT_POS_DIR if self.forward else _KPOINT_NEG_DIR + direction = mp.Vector3(*[float(v == 0) for v in self.volume.size]) + kpoint_func = lambda *not_used: direction if self.forward else -1 * direction overlap_idx = 0 ob = self.sim.get_eigenmode_coefficients( self._monitor, @@ -218,10 +223,11 @@ def __call__(self): class FourierFields(ObjectiveQuantity): - def __init__(self, sim, volume, component): + def __init__(self, sim, volume, component, decimation_factor=0): super().__init__(sim) self.volume = volume self.component = component + self.decimation_factor = decimation_factor def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) @@ -230,6 +236,7 @@ def register_monitors(self, frequencies): self._frequencies, where=self.volume, yee_grid=False, + decimation_factor=self.decimation_factor, ) return self._monitor @@ -281,7 +288,7 @@ def place_adjoint_source(self, dJ): for yi in range(y_dim): for xi in range(x_dim): '''We only need to add a current source if the - jacobian is nonzero for all frequencies at + jacobian is nonzero for all frequencies at that particular point. Otherwise, the fitting algorithm is going to fail. ''' @@ -316,18 +323,19 @@ def __call__(self): class Near2FarFields(ObjectiveQuantity): - def __init__(self, sim, Near2FarRegions, far_pts): + def __init__(self, sim, Near2FarRegions, far_pts, decimation_factor=0): super().__init__(sim) self.Near2FarRegions = Near2FarRegions self.far_pts = far_pts #list of far pts self._nfar_pts = len(far_pts) + self.decimation_factor = decimation_factor def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) self._monitor = self.sim.add_near2far( self._frequencies, *self.Near2FarRegions, - yee_grid=True, + decimation_factor=self.decimation_factor, ) return self._monitor @@ -359,7 +367,7 @@ def place_adjoint_source(self, dJ): time_src.frequency, self._frequencies, scale, - dt, + self.sim.fields.dt, ) (num_basis, num_pts) = src.nodes.shape for basis_i in range(num_basis): diff --git a/python/tests/test_adjoint_jax.py b/python/tests/test_adjoint_jax.py index 8481956eb..f15e36478 100644 --- a/python/tests/test_adjoint_jax.py +++ b/python/tests/test_adjoint_jax.py @@ -1,7 +1,7 @@ import unittest import parameterized -from utils import VectorComparisonMixin +from utils import ApproxComparisonTestCase import jax import jax.numpy as jnp @@ -16,7 +16,7 @@ _FD_STEP = 1e-4 # The tolerance for the adjoint and finite difference gradient comparison -_TOL = 2e-2 +_TOL = 0.1 if mp.is_single_precision() else 0.02 mp.verbosity(0) @@ -125,7 +125,8 @@ def build_straight_wg_simulation( simulation, mp.Volume(center=center, size=monitor_size), mode=1, - forward=forward) for center in monitor_centers for forward in [True, False] + forward=True, + decimation_factor=1) for center in monitor_centers for forward in [True, False] ] return simulation, sources, monitors, design_regions, frequencies @@ -174,15 +175,15 @@ def test_design_region_monitor_helpers(self): self.assertEqual(value.dtype, onp.complex128) -class WrapperTest(VectorComparisonMixin, unittest.TestCase): +class WrapperTest(ApproxComparisonTestCase): @parameterized.parameterized.expand([ ('1500_1550bw_01relative_gaussian_port1', onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 1.0, 0), - # ('1550_1600bw_02relative_gaussian_port1', onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 0), - # ('1500_1600bw_03relative_gaussian_port1', onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 0), + ('1550_1600bw_02relative_gaussian_port1', onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 0), + ('1500_1600bw_03relative_gaussian_port1', onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 0), ('1500_1550bw_01relative_gaussian_port2', onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 1.0, 1), - # ('1550_1600bw_02relative_gaussian_port2', onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 1), - # ('1500_1600bw_03relative_gaussian_port2', onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 1), + ('1550_1600bw_02relative_gaussian_port2', onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 1), + ('1500_1600bw_03relative_gaussian_port2', onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 1), ]) def test_wrapper_gradients(self, _, frequencies, gaussian_rel_width, design_variable_fill_value, excite_port_idx): """Tests gradient from the JAX-Meep wrapper against finite differences.""" @@ -246,7 +247,7 @@ def loss_fn(x, excite_port_idx=0): fd_projection = onp.stack(fd_projection) # Check that dp . ∇T ~ T(p + dp) - T(p) - self.assertVectorsClose( + self.assertClose( projection, fd_projection, epsilon=_TOL, diff --git a/src/mpb.cpp b/src/mpb.cpp index a3b736c31..6bc865a18 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -484,7 +484,7 @@ void *fields::get_eigenmode(double frequency, direction d, const volume where, c // which we automatically pick if kmatch == 0. if (match_frequency && kmatch == 0) { vec cen = eig_vol.center(); - kmatch = copysign(frequency * sqrt(real(get_eps(cen, frequency)) * real(get_mu(cen, frequency))), _kpoint.in_direction(d)); + kmatch = frequency * sqrt(real(get_eps(cen, frequency)) * real(get_mu(cen, frequency))); if (d == NO_DIRECTION) { for (int i = 0; i < 3; ++i) k[i] = dot_product(R[i], kdir) * kmatch; // kdir*kmatch in reciprocal basis @@ -831,9 +831,6 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d // electric current K = nHat \times H */ // magnetic current N = -nHat \times E */ /*--------------------------------------------------------------*/ - if (global_eigenmode_data->group_velocity < 0) - master_printf("vg is negative!\n"); - amp *= -1; // equivalent to flipping the direction of nhat. if (is_D(c0)) c0 = direction_component(Ex, component_direction(c0)); if (is_B(c0)) c0 = direction_component(Hx, component_direction(c0)); component cE[3] = {Ex, Ey, Ez}, cH[3] = {Hx, Hy, Hz};