Added topology optimization and enhancements

Added a major new feature to perform topology optimization. This
supports 2D optimization and 3D optimization with an extruded quasi-2D
layer. The topology optimization is constrained to binary materials
and respects a minimum feature size which the user can control.

Added a number of performance enhancements to the shape based
optimization to support grating coupler devices with large numbers of
parameters and using the high presision volume averrage mesh refnement
strategy.

Example files have been added and updated to demonstrate the new
features included in this commit.

This commit corresponds to the version shipped with Lumerical FDTD 2019b
R2

QA/adaptive_gradient_descent_maximization_test.py

@@ -0,0 +1,105 @@
+""" Copyright chriskeraly
+    Copyright (c) 2019 Lumerical Inc. """
+
+
+import sys
+sys.path.append(".")
+import os
+import numpy as np
+
+from qatools import *
+
+from lumopt.optimizers.adaptive_gradient_descent import AdaptiveGradientDescent
+
+class AdaptiveGradientDescentMaximizationTest(TestCase):
+    """ 
+        Unit test for the AdaptiveGradientDescent class. It performs a quick sanity check that the optimizer can indeed maximize a figure of merit
+        with a known unique maximum within the provided bounds.
+    """
+
+    machine_eps = np.finfo(float).eps
+
+    def test_single_parameter_maximization(self):
+        optimizer = AdaptiveGradientDescent(max_dx = 0.5,
+                                            min_dx = 1.0e-12, 
+                                            max_iter = 7, 
+                                            dx_regrowth_factor = 2.0,
+                                            all_params_equal = True,
+                                            scaling_factor = 1.0)
+        fom = lambda param: -1.0/np.sin(np.pi*param)
+        jac = lambda param: np.pi*np.cos(np.pi*param) / np.power(np.sin(np.pi*param), 2)
+        start = np.array([1.0e-6])
+        bounds = np.array([(self.machine_eps, 1.0-self.machine_eps)])
+        def plot_fun(): pass
+        optimizer.initialize(start_params = start, callable_fom = fom, callable_jac = jac, bounds = bounds, plotting_function = plot_fun)
+        results = optimizer.run()
+        self.assertAlmostEqual(results['x'].size, 1)
+        self.assertAlmostEqual(results['x'][0], 0.5, 8)
+        self.assertAlmostEqual(results['fun'], -1.0, 8)
+        self.assertLessEqual(results['nit'], 7)
+
+    def test_single_parameter_maximization_with_scaling(self):
+        span = 1.0e-9
+        optimizer = AdaptiveGradientDescent(max_dx = 0.5*span,
+                                           min_dx = 1.0e-12*span,
+                                           max_iter = 7, 
+                                           dx_regrowth_factor = 2.0,
+                                           all_params_equal = True,
+                                           scaling_factor = 1.0/span)
+        fom = lambda param: -1.0/np.sin(np.pi/span*param[0])
+        jac = lambda param: np.pi*np.cos(np.pi/span*param[0])/np.power(np.sin(np.pi/span*param[0]),2)/span
+        start = np.array(1.0e-6)*span
+        bounds = np.array([(self.machine_eps, 1.0-self.machine_eps)])*span
+        def plot_fun(): pass
+        optimizer.initialize(start_params = start, callable_fom = fom, callable_jac = jac, bounds = bounds, plotting_function = plot_fun)
+        results = optimizer.run()
+        self.assertAlmostEqual(results['x'].size, 1)
+        self.assertAlmostEqual(results['x'][0], 0.5*span, 8)
+        self.assertAlmostEqual(results['fun'], -1.0, 8)
+        self.assertLessEqual(results['nit'], 7)
+
+    def test_two_parameter_maximization(self):
+        optimizer = AdaptiveGradientDescent(max_dx = 0.5*np.ones(2),
+                                            min_dx = 1.0e-12*np.ones(2),
+                                            max_iter = 8, 
+                                            dx_regrowth_factor = 2.0,
+                                            all_params_equal = False,
+                                            scaling_factor = np.ones(2))
+        fom = lambda params: np.sin(np.pi*params[0])*np.sin(np.pi*params[1])
+        jac = lambda params: np.pi*np.array([np.cos(np.pi*params[0])*np.sin(np.pi*params[1]),
+                                             np.sin(np.pi*params[0])*np.cos(np.pi*params[1])])
+        start = np.array([1.0e-5, 1.0-1e-5])
+        bounds = np.array([(self.machine_eps, 1.0-self.machine_eps), (self.machine_eps, 1.0-self.machine_eps)])
+        def plot_fun(): pass
+        optimizer.initialize(start_params = start, callable_fom = fom, callable_jac = jac, bounds = bounds, plotting_function = plot_fun)
+        results = optimizer.run()
+        self.assertAlmostEqual(results['x'].size, 2)
+        self.assertAlmostEqual(results['x'][0], 0.5, 8)
+        self.assertAlmostEqual(results['x'][1], 0.5, 8)
+        self.assertAlmostEqual(results['fun'], 1.0, 8)
+        self.assertLessEqual(results['nit'], 8)
+
+    def test_two_parameter_maximization(self):
+        span = 1.0e-6
+        optimizer = AdaptiveGradientDescent(max_dx = np.array([0.5*span, 0.5]),
+                                            min_dx = np.array([1.0e-12*span, 1.0e-12]),
+                                            max_iter = 8, 
+                                            dx_regrowth_factor = 2.0,
+                                            all_params_equal = False,
+                                            scaling_factor = np.array([1.0/span, 1.0]))
+        fom = lambda params: np.sin(np.pi/span*params[0])*np.sin(np.pi*params[1])
+        jac = lambda params: np.pi*np.array([np.cos(np.pi/span*params[0])*np.sin(np.pi*params[1])/span,
+                                             np.sin(np.pi/span*params[0])*np.cos(np.pi*params[1])])
+        start = np.array([1.0e-5*span, 1.0-1e-5])
+        bounds = np.array([(self.machine_eps*span, (1.0-self.machine_eps)*span), (self.machine_eps, 1.0-self.machine_eps)])
+        def plot_fun(): pass
+        optimizer.initialize(start_params = start, callable_fom = fom, callable_jac = jac, bounds = bounds, plotting_function = plot_fun)
+        results = optimizer.run()
+        self.assertAlmostEqual(results['x'].size, 2)
+        self.assertAlmostEqual(results['x'][0], 0.5*span, 8)
+        self.assertAlmostEqual(results['x'][1], 0.5, 8)
+        self.assertAlmostEqual(results['fun'], 1.0, 8)
+        self.assertLessEqual(results['nit'], 8)
+
+if __name__ == "__main__":
+    run([__file__])

QA/base_script_test.py → QA/base_fdtd_script_test.py

@@ -1,5 +1,4 @@
-""" Copyright chriskeraly
-    Copyright (c) 2019 Lumerical Inc. """
+""" Copyright (c) 2019 Lumerical Inc. """
 
 import sys
 sys.path.append(".")
@@ -10,18 +9,18 @@
 from lumopt.utilities.simulation import Simulation
 from lumopt.utilities.base_script import BaseScript
 
-class TestBaseScript(TestCase):
+class TestFDTDBaseScript(TestCase):
     """ 
         Unit test for BaseScript class. It verifies that the object is able to run an *.lsf script, a *.fsp project file or a plain script in a string.
     """
 
     file_dir = os.path.abspath(os.path.dirname(__file__))
 
     def setUp(self):
-        self.sim = Simulation(workingDir = self.file_dir, hide_fdtd_cad = True)
+        self.sim = Simulation(workingDir = self.file_dir, use_var_fdtd = False, hide_fdtd_cad = True)
 
     def test_eval_project_file(self):
-        my_project_file = os.path.join(self.file_dir,'base_script_test.fsp')
+        my_project_file = os.path.join(self.file_dir,'base_fdtd_script_test.fsp')
         base_script_obj = BaseScript(my_project_file)
         base_script_obj(self.sim.fdtd)
         self.assertTrue(self.sim.fdtd.getnamednumber('FDTD') == 1)
@@ -35,13 +34,13 @@ def test_eval_python_script(self):
         self.assertTrue(self.sim.fdtd.getnamednumber('FDTD') == 1)
 
     def test_eval_script_file(self):
-        my_script_file = os.path.join(self.file_dir,'base_script_test.lsf')
+        my_script_file = os.path.join(self.file_dir,'base_fdtd_script_test.lsf')
         base_script_obj = BaseScript(my_script_file)
         base_script_obj.eval(self.sim.fdtd)
         self.assertTrue(self.sim.fdtd.getnamednumber('FDTD') == 1)
 
     def test_eval_script(self):
-        my_script = "load('base_script_test.fsp');"
+        my_script = "load('base_fdtd_script_test.fsp');"
         base_script_obj = BaseScript(my_script)
         base_script_obj.eval(self.sim.fdtd)
         self.assertTrue(self.sim.fdtd.getnamednumber('FDTD') == 1)

QA/base_varfdtd_script_test.lsf

@@ -0,0 +1,10 @@
+            ################
+
+    deleteall;
+                                    addvarfdtd;
+
+            addpoly;   
+                           addmodesource;
+
+            ###########
+

QA/base_varfdtd_script_test.py

@@ -0,0 +1,49 @@
+""" Copyright (c) 2019 Lumerical Inc. """
+
+import sys
+sys.path.append(".")
+import os
+
+from qatools import *
+
+from lumopt.utilities.simulation import Simulation
+from lumopt.utilities.base_script import BaseScript
+
+class TestVarFDTDBaseScript(TestCase):
+    """ 
+        Unit test for BaseScript class. It verifies that the object is able to run an *.lsf script, a *.lms project file or a plain script in a string.
+    """
+
+    file_dir = os.path.abspath(os.path.dirname(__file__))
+
+    def setUp(self):
+        self.sim = Simulation(workingDir = self.file_dir, use_var_fdtd = True, hide_fdtd_cad = True)
+
+    def test_eval_project_file(self):
+        my_project_file = os.path.join(self.file_dir,'base_varfdtd_script_test.lms')
+        base_script_obj = BaseScript(my_project_file)
+        base_script_obj(self.sim.fdtd)
+        self.assertTrue(self.sim.fdtd.getnamednumber('varFDTD') == 1)
+        self.assertTrue(self.sim.fdtd.getnamednumber('polygon') == 1)
+
+    def test_eval_python_script(self):
+        my_fun = lambda fdtd_handle: fdtd_handle.addvarfdtd()
+        base_script_obj = BaseScript(my_fun)
+        base_script_obj.eval(self.sim.fdtd)
+        self.assertTrue(self.sim.fdtd.getnamednumber('varFDTD') == 1)
+
+    def test_eval_script_file(self):
+        my_script_file = os.path.join(self.file_dir,'base_varfdtd_script_test.lsf')
+        base_script_obj = BaseScript(my_script_file)
+        base_script_obj.eval(self.sim.fdtd)
+        self.assertTrue(self.sim.fdtd.getnamednumber('varFDTD') == 1)
+
+    def test_eval_script(self):
+        my_script = "load('base_varfdtd_script_test.fsp');"
+        base_script_obj = BaseScript(my_script)
+        base_script_obj.eval(self.sim.fdtd)
+        self.assertTrue(self.sim.fdtd.getnamednumber('varFDTD') == 1)
+        self.assertTrue(self.sim.fdtd.getnamednumber('polygon') == 1)
+
+if __name__ == "__main__":
+    run([__file__])

QA/co_optimization_parallel_plate_waveguide_test.py

@@ -0,0 +1,104 @@
+""" Copyright chriskeraly
+    Copyright (c) 2019 Lumerical Inc. """
+
+import sys
+sys.path.append(".")
+import os
+import numpy as np
+
+from qatools import *
+
+from lumopt.utilities.wavelengths import Wavelengths
+from lumopt.utilities.materials import Material
+from lumopt.geometries.polygon import FunctionDefinedPolygon
+from lumopt.figures_of_merit.modematch import ModeMatch
+from lumopt.optimizers.generic_optimizers import ScipyOptimizers
+from lumopt.optimization import Optimization
+
+class TestCoOptimizationParallelPlateWaveguide(TestCase):
+    """ 
+        Unit test for the Optimization class. It performs a co-optimization using a parallel plate waveguide
+        filled by a dielectric excited with two different polarizations (TE and TM). The waveguide has a gap
+        that must be filled all the way to maximize transmission.
+
+    """
+
+    file_dir = os.path.abspath(os.path.dirname(__file__))
+
+    def setUp(self):
+        # Base simulation project files
+        self.base_TE_sim = os.path.join(self.file_dir, 'co_optimization_parallel_plate_waveguide_TE_base.fsp')
+        self.base_TM_sim = os.path.join(self.file_dir, 'co_optimization_parallel_plate_waveguide_TM_base.fsp')
+        # Simulation bandwidth        
+        self.wavelengths = Wavelengths(start = 1500e-9,
+                                       stop = 1600e-9,
+                                       points = 12)
+        # Polygon defining a rectangle that can grow or shrink along the y-axis to fill the gap
+        self.mesh_del = 10.0e-9; # must be kept in sych with self.base_script
+        initial_points_y = np.array([1.75 * self.mesh_del, 0.01 * self.mesh_del])
+        def wall(param = initial_points_y):
+            assert param.size == 2, "walls defined by two points."
+            self.wg_gap = 10.0 * self.mesh_del # must be kept in sych
+            points_x = 0.5 * np.array([-self.wg_gap, self.wg_gap, self.wg_gap, -self.wg_gap])
+            points_y = np.array([-param[0], -param[1], param[1], param[0]])
+            polygon_points = [(x, y) for x, y in zip(points_x, points_y)]
+            return np.array(polygon_points)
+        self.wg_width = 50.0 * self.mesh_del # must be kept in synch
+        bounds = [(0.0, self.wg_width / 2.0)] * initial_points_y.size 
+        self.geometry = FunctionDefinedPolygon(func = wall, 
+                                               initial_params = initial_points_y, 
+                                               bounds = bounds,
+                                               z = 0.0, # must be kept in sych
+                                               depth = self.wg_width,
+                                               eps_out = Material(base_epsilon = 1.0 ** 2, name = '<Object defined dielectric>', mesh_order = 2), # must be kept in synch
+                                               eps_in = Material(base_epsilon = 4.0 ** 2, name = '<Object defined dielectric>', mesh_order = 1), # must be kept in sych
+                                               edge_precision = 50,
+                                               dx = 1.0e-10)
+        # Figure of merit
+        self.fom = ModeMatch(monitor_name = 'fom', # must be kept in sych
+                             mode_number = 1, # must be kept in sych
+                             direction = 'Forward',
+                             multi_freq_src = True,
+                             target_T_fwd = lambda wl: np.ones(wl.size),
+                             norm_p = 1)
+        # Scipy optimizer
+        self.optimizer = ScipyOptimizers(max_iter = 5, 
+                                         method = 'L-BFGS-B',
+                                         scaling_factor = 1.0e7,
+                                         pgtol = 1.0e-5,
+                                         ftol = 1.0e-12,
+                                         target_fom = 0.0,
+                                         scale_initial_gradient_to = None)
+
+    def test_co_optimization_in_2D(self):
+        print("2D TE-TM co-optimization (use_deps = True): ")
+        #base_script, wavelengths, fom, geometry, optimizer, use_var_fdtd = False, hide_fdtd_cad = False, use_deps = True, plot_history = True, store_all_simulations = True
+        optTE = Optimization(base_script = self.base_TE_sim, 
+                             wavelengths = self.wavelengths,
+                             fom = self.fom,
+                             geometry = self.geometry,
+                             optimizer = self.optimizer,
+                             use_var_fdtd = False,
+                             hide_fdtd_cad = True,
+                             use_deps = True,
+                             plot_history = False,
+                             store_all_simulations = False)
+        optTM = Optimization(base_script = self.base_TM_sim, 
+                             wavelengths = self.wavelengths,
+                             fom = self.fom,
+                             geometry = self.geometry,
+                             optimizer = self.optimizer,
+                             use_var_fdtd = False,
+                             hide_fdtd_cad = True,
+                             use_deps = True,
+                             plot_history = False,
+                             store_all_simulations = False)
+        opt = optTE + optTM
+        fom, params = opt.run()
+        self.assertGreaterEqual(fom, 1.99990)
+        reference_value = self.wg_width / 2.0 * self.optimizer.scaling_factor[0]
+        self.assertAlmostEqual(params[0], reference_value)
+        self.assertAlmostEqual(params[1], reference_value)
+
+if __name__ == "__main__":
+    run([__file__])

QA/fixed_step_gradient_maximization_test.py

@@ -0,0 +1,40 @@
+""" Copyright chriskeraly
+    Copyright (c) 2019 Lumerical Inc. """
+
+
+import sys
+sys.path.append(".")
+import os
+import numpy as np
+
+from qatools import *
+
+from lumopt.optimizers.fixed_step_gradient_descent import FixedStepGradientDescent
+
+class FixedStepGradientDescentMaximizationTest(TestCase):
+    """ 
+        Unit test for the AdaptiveGradientDescent class. 
+    """
+
+    machine_eps = np.finfo(float).eps
+
+    def test_single_parameter_maximization(self):
+        optimizer = FixedStepGradientDescent(max_dx = 0.005,
+                                             max_iter = 99, 
+                                             all_params_equal = False,
+                                             noise_magnitude = 0.0,
+                                             scaling_factor = 1.0)
+        fom = lambda param: -1.0/np.sin(np.pi*param[0])
+        jac = lambda param: np.pi*np.cos(np.pi*param[0]) / np.power(np.sin(np.pi*param[0]), 2)
+        start = np.array(0.005)
+        bounds = np.array([(self.machine_eps, 1.0-self.machine_eps)])
+        def plot_fun(): pass
+        optimizer.initialize(start_params = start, callable_fom = fom, callable_jac = jac, bounds = bounds, plotting_function = plot_fun)
+        results = optimizer.run()
+        self.assertAlmostEqual(results['x'].size, 1)
+        self.assertAlmostEqual(results['x'][0], 0.5, 8)
+        self.assertAlmostEqual(results['fun'], -1.0, 3)
+        self.assertLessEqual(results['nit'], 99)
+
+if __name__ == "__main__":
+    run([__file__])