Merge branch 'NanoComp:master' into sharded-dump

.github/workflows/build-ci.yml

@@ -140,7 +140,7 @@ jobs:
         echo "MEEP_VERSION=${MEEP_VERSION}" >> $GITHUB_ENV
 
     - name: Run configure
-      if: ${{ !(matrix.enable-mpi == false && matrix.python-version == 3.6) }}
+      if: ${{ !(matrix.enable-mpi == false && matrix.python-version == 3.6) && !(matrix.enable-mpi == true && matrix.python-version == 3.9) }}
       run: |
         mkdir -p build &&
         pushd build &&
@@ -149,7 +149,15 @@ jobs:
 
     - name: Run configure
       if: ${{ matrix.enable-mpi == false && matrix.python-version == 3.6 }}
-      run: ./configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF}
+      run: ./configure --enable-maintainer-mode --with-coverage --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF}
+
+    - name: Run configure with single-precision floating point
+      if: ${{ matrix.enable-mpi == true && matrix.python-version == 3.9 }}
+      run: |
+        mkdir -p build &&
+        pushd build &&
+        ../configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF} --enable-single &&
+        popd
 
     - name: Run make
       if: ${{ matrix.enable-mpi == false && matrix.python-version == 3.6 }}

configure.ac

@@ -514,7 +514,7 @@ fi
 
 ##############################################################################
 # Miscellaneous function and header checks
-AC_CHECK_HEADERS([sys/time.h])
+AC_CHECK_HEADERS([sys/time.h immintrin.h])
 AC_CHECK_FUNCS([BSDgettimeofday gettimeofday cblas_ddot cblas_daxpy jn])
 
 ##############################################################################
@@ -538,6 +538,9 @@ if test "$USE_MAINTAINER_MODE" = yes && (test "x$with_python" = "xyes" || test "
   AC_MSG_ERROR([SWIG not found; configure --without-python --without-scheme or use a release tarball])
 fi
 
+AC_ARG_WITH(coverage, [AS_HELP_STRING([--with-coverage],[enable Python coverage tests])],
+            with_coverage=$withval, with_coverage=no)
+
 if test "x$with_python" = xno; then
   have_python=no
 else
@@ -572,14 +575,15 @@ else
           AC_MSG_WARN([disabling Python wrappers])
           have_python=no],[#include <Python.h>])
 
-        AC_MSG_CHECKING([for coverage module])
-        $PYTHON -c 'import coverage' 2>/dev/null
-        if test $? = 0; then
-          AC_MSG_RESULT([yes])
-          have_coverage=yes
-        else
-          AC_MSG_RESULT([no])
-          have_coverage=no
+        if test x"$with_coverage" = xyes; then
+          AC_MSG_CHECKING([for coverage module])
+          $PYTHON -c 'import coverage' 2>/dev/null
+          if test $? = 0; then
+            AC_MSG_RESULT([yes])
+          else
+            AC_MSG_RESULT([no])
+            with_coverage=no
+          fi
         fi
       fi
 
@@ -591,7 +595,7 @@ fi # with_python
 
 AC_SUBST(PYTHON_INCLUDES)
 AM_CONDITIONAL(WITH_PYTHON, test x"$have_python" = "xyes")
-AM_CONDITIONAL(WITH_COVERAGE, test x"$have_coverage" = "xyes")
+AM_CONDITIONAL(WITH_COVERAGE, test x"$with_coverage" = "xyes")
 
 if test "x$with_scheme" = xyes; then
    # Copy/symlink casimir.scm and materials.scm to builddir for out-of-tree builds

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 not self.sim.force_complex_fields:
+            scale *= 2
         return scale
 
     def _create_time_profile(self, fwidth_frac=0.1):
@@ -262,7 +266,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.
                         '''
@@ -340,7 +344,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):

python/geom.py

@@ -517,7 +517,7 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct
         # Account for conductivity term (only multiply if nonzero to avoid unnecessary complex numbers)
         conductivity = np.expand_dims(Matrix(diag=conductivity_diag,offdiag=conductivity_offdiag),axis=0)
         if np.count_nonzero(conductivity) > 0:
-            epsmu = (1 + 1j/freqs * conductivity) * epsmu
+            epsmu = (1 + 1j/(2*np.pi*freqs) * conductivity) * epsmu
 
         # Convert list matrix to 3D numpy array size [freqs,3,3]
         return np.squeeze(epsmu)

python/simulation.py

@@ -40,6 +40,7 @@
 
 verbosity = Verbosity(mp.cvar, 'meep', 1)
 
+mp.set_zero_subnormals(True)
 
 # Send output from Meep, ctlgeom, and MPB to Python's stdout
 mp.set_meep_printf_callback(mp.py_master_printf_wrap)

python/tests/test_3rd_harm_1d.py

@@ -2,10 +2,9 @@
 
 import unittest
 import meep as mp
-import numpy as np
+from utils import ApproxComparisonTestCase
 
-
-class Test3rdHarm1d(unittest.TestCase):
+class Test3rdHarm1d(ApproxComparisonTestCase):
 
     def setUp(self):
         self.sz = 100
@@ -53,7 +52,8 @@ def test_3rd_harm_1d(self):
 
         harmonics = [self.k, self.amp, mp.get_fluxes(self.trans1)[0], mp.get_fluxes(self.trans3)[0]]
 
-        np.testing.assert_allclose(expected_harmonics, harmonics)
+        tol = 3e-5 if mp.is_single_precision() else 1e-7
+        self.assertClose(expected_harmonics, harmonics, epsilon=tol)
 
 
 if __name__ == '__main__':

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)
 
@@ -166,7 +166,7 @@ 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', onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 1.0),
@@ -237,7 +237,7 @@ def loss_fn(x):
     fd_projection = onp.stack(fd_projection)
 
     # Check that dp . ∇T ~ T(p + dp) - T(p)
-    self.assertVectorsClose(
+    self.assertClose(
       projection,
       fd_projection,
       epsilon=_TOL,

python/tests/test_adjoint_solver.py

@@ -8,6 +8,7 @@
 from autograd import tensor_jacobian_product
 import unittest
 from enum import Enum
+from utils import ApproxComparisonTestCase
 
 MonitorObject = Enum('MonitorObject', 'EIGENMODE DFT')
 
@@ -52,13 +53,13 @@
                               eig_parity=eig_parity)]
 
 
-def forward_simulation(design_params,mon_type,frequencies=None):
+def forward_simulation(design_params,mon_type,frequencies=None, use_complex=False):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
                               silicon,
                               weights=design_params.reshape(Nx,Ny),
                               grid_type='U_MEAN')
-            
+
     matgrid_geometry = [mp.Block(center=mp.Vector3(),
                                  size=mp.Vector3(design_shape.x,design_shape.y,0),
                                  material=matgrid)]
@@ -69,7 +70,8 @@ def forward_simulation(design_params,mon_type,frequencies=None):
                         cell_size=cell_size,
                         boundary_layers=boundary_layers,
                         sources=sources,
-                        geometry=geometry)
+                        geometry=geometry,
+                        force_complex_fields=use_complex)
     if not frequencies:
         frequencies = [fcen]
 
@@ -106,7 +108,7 @@ def forward_simulation(design_params,mon_type,frequencies=None):
         return Ez2
 
 
-def adjoint_solver(design_params, mon_type, frequencies=None):
+def adjoint_solver(design_params, mon_type, frequencies=None, use_complex=False):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
                               silicon,
@@ -126,7 +128,8 @@ def adjoint_solver(design_params, mon_type, frequencies=None):
                         cell_size=cell_size,
                         boundary_layers=boundary_layers,
                         sources=sources,
-                        geometry=geometry)
+                        geometry=geometry,
+                        force_complex_fields=use_complex)
     if not frequencies:
         frequencies = [fcen]
 
@@ -170,7 +173,7 @@ def mapping(x,filter_radius,eta,beta):
     return projected_field.flatten()
 
 
-class TestAdjointSolver(unittest.TestCase):
+class TestAdjointSolver(ApproxComparisonTestCase):
 
     def test_adjoint_solver_DFT_fields(self):
         print("*** TESTING DFT ADJOINT FEATURES ***")
@@ -181,47 +184,50 @@ def test_adjoint_solver_DFT_fields(self):
 
             ## compute unperturbed S12
             S12_unperturbed = forward_simulation(p, MonitorObject.DFT, frequencies)
-            
+
             ## compare objective results
             print("|Ez|^2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+            self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(p+dp, MonitorObject.DFT, frequencies)
-            
+
             ## compare gradients
             if adjsol_grad.ndim < 2:
                 adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
             adj_scale = (dp[None,:]@adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+            tol = 0.3 if mp.is_single_precision() else 0.01
+            self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
 
     def test_adjoint_solver_eigenmode(self):
         print("*** TESTING EIGENMODE ADJOINT FEATURES ***")
         ## test the single frequency and multi frequency cases
-        for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]:
-            ## compute gradient using adjoint solver
-            adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies)
+        for use_complex in [True, False]:
+            for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]:
+                ## compute gradient using adjoint solver
+                adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies, use_complex)
 
-            ## compute unperturbed S12
-            S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies)
-
-            ## compare objective results
-            print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+                ## compute unperturbed S12
+                S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, use_complex)
 
-            ## compute perturbed S12
-            S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies)
-
-            ## compare gradients
-            if adjsol_grad.ndim < 2:
-                adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
-            adj_scale = (dp[None,:]@adjsol_grad).flatten()
-            fd_grad = S12_perturbed-S12_unperturbed
-            print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+                ## compare objective results
+                print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
+                self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
+
+                ## compute perturbed S12
+                S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, use_complex)
+
+                ## compare gradients
+                if adjsol_grad.ndim < 2:
+                    adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
+                adj_scale = (dp[None,:]@adjsol_grad).flatten()
+                fd_grad = S12_perturbed-S12_unperturbed
+                print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
+                tol = 0.1 if mp.is_single_precision() else 0.03
+                self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
 
     def test_gradient_backpropagation(self):
@@ -250,7 +256,8 @@ def test_gradient_backpropagation(self):
 
             ## compare objective results
             print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+            tol = 1e-2 if mp.is_single_precision() else 1e-3
+            self.assertClose(adjsol_obj,S12_unperturbed,epsilon=tol)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(mapping(p+dp,filter_radius,eta,beta), MonitorObject.EIGENMODE,frequencies)
@@ -260,7 +267,8 @@ def test_gradient_backpropagation(self):
             adj_scale = (dp[None,:]@bp_adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+            tol = 0.04 if mp.is_single_precision() else 0.03
+            self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
 
 if __name__ == '__main__':

python/tests/test_array_metadata.py

@@ -1,8 +1,9 @@
-import unittest
 import meep as mp
+import unittest
 import numpy as np
+from utils import ApproxComparisonTestCase
 
-class TestArrayMetadata(unittest.TestCase):
+class TestArrayMetadata(ApproxComparisonTestCase):
 
     def test_array_metadata(self):
         resolution = 25
@@ -80,7 +81,8 @@ def vec_func(r):
         vec_func_sum = np.sum(W*(xm**2 + 2*ym**2))
         pulse_modal_volume = np.sum(W*EpsE2)/np.max(EpsE2) * vec_func_sum
 
-        self.assertAlmostEqual(cw_modal_volume/pulse_modal_volume, 1.00, places=2)
+        tol = 5e-2 if mp.is_single_precision() else 1e-2
+        self.assertClose(cw_modal_volume/pulse_modal_volume, 1.0, epsilon=tol)
 
 if __name__ == '__main__':
     unittest.main()