{"id":427139,"date":"2025-03-28T00:32:59","date_gmt":"2025-03-28T00:32:59","guid":{"rendered":"https:\/\/innovationspace.ansys.com\/forum\/forums\/topic\/phase-optimization-of-23-mmi-using-lumopt-is-ineffective\/"},"modified":"2025-03-28T12:57:46","modified_gmt":"2025-03-28T12:57:46","slug":"phase-optimization-of-23-mmi-using-lumopt-is-ineffective","status":"closed","type":"topic","link":"https:\/\/innovationspace.ansys.com\/forum\/forums\/topic\/phase-optimization-of-23-mmi-using-lumopt-is-ineffective\/","title":{"rendered":"Phase optimization of 2\u00d73 MMI using Lumopt is ineffective."},"content":{"rendered":"<p>&lt;p&gt;We modified the modematch code and added phase control to the optimization objective to achieve the inverse design of a 2&#215;3 MMI device. When two mode light sources enter through two input ports respectively, we aim to obtain equal power at the three output ports with a phase difference of 120&deg;. However, we observed that the actual FOM curve does not trend towards a gradual reduction in insertion loss, and the phases of the monitors at each port are all 0 as the number of optimization iterations increases. Why is this the case? Please provide me with professional advice and assistance. Thank you!&lt;br&gt;&lt;br&gt;code of modematch\uff1a&lt;br&gt;&#8221;&#8221;&#8221; Copyright chriskeraly&lt;br&gt;Copyright (c) 2019 Lumerical Inc. &#8220;&#8221;&#8221;&lt;br&gt;&lt;br&gt;import sys&lt;br&gt;import numpy as np&lt;br&gt;import scipy as sp&lt;br&gt;import scipy.constants&lt;br&gt;import lumapi&lt;br&gt;&lt;br&gt;from lumopt.utilities.wavelengths import Wavelengths&lt;br&gt;&lt;br&gt;&#8221;&#8221;&#8221; Helper function to determine if a variable can be converted to an integer &#8220;&#8221;&#8221;&lt;br&gt;def is_int(x):&lt;br&gt;try:&lt;br&gt;int(x)&lt;br&gt;return True&lt;br&gt;except ValueError:&lt;br&gt;return False&lt;br&gt;&lt;br&gt;class ModeMatch(object):&lt;br&gt;&lt;br&gt;&#8221;&#8221;&#8221; Calculates the figure of merit from an overlap integral between the fields recorded by a field monitor and the slected mode.&lt;br&gt;A mode expansion monitor is added to the field monitor to calculate the overlap result, which appears as T_forward in the&lt;br&gt;list of mode expansion monitor results. The T_forward result is described in the following page:&lt;br&gt;&lt;br&gt;https:\/\/kb.lumerical.com\/ref_sim_obj_using_mode_expansion_monitors.html&lt;br&gt;&lt;br&gt;This result is equivalent to equation (7) in the following paper:&lt;br&gt;&lt;br&gt;C. Lalau-Keraly, S. Bhargava, O. Miller, and E. Yablonovitch, &#8220;Adjoint shape optimization applied to electromagnetic design,&#8221;&lt;br&gt;Opt. Express 21, 21693-21701 (2013). https:\/\/doi.org\/10.1364\/OE.21.021693&lt;br&gt;&lt;br&gt;Parameters&lt;br&gt;&#8212;&#8212;&#8212;-&lt;br&gt;:param monitor_name: name of the field monitor that records the fields to be used in the mode overlap calculation.&lt;br&gt;:param mode_number: mode number in the list of modes generated by the mode expansion monitor.&lt;br&gt;:param direction: direction of propagation (&#8216;Forward&#8217; or &#8216;Backward&#8217;) of the mode injected by the source.&lt;br&gt;:param multi_freq_src: bool flag to enable \/ disable a multi-frequency mode calculation and injection for the adjoint source.&lt;br&gt;:param target_T_fwd: function describing the target T_forward vs wavelength (see documentation for mode expansion monitors).&lt;br&gt;:param norm_p: exponent of the p-norm used to generate the figure of merit; use to generate the FOM.&lt;br&gt;:param target_fom: A target value for the figure of merit. This allows to print\/plot the distance of the current&lt;br&gt;design from a target value&lt;br&gt;&#8221;&#8221;&#8221;&lt;br&gt;&lt;br&gt;def __init__(self, monitor_name, mode_number, direction, target_phase=0, power_weight=1.0, phase_weight=1.0, multi_freq_src = False, target_T_fwd = lambda wl: np.ones(wl.size), norm_p = 1, target_fom = 0):&lt;br&gt;self.monitor_name = str(monitor_name)&lt;br&gt;if not self.monitor_name:&lt;br&gt;raise UserWarning(&#8217;empty monitor name.&#8217;)&lt;br&gt;self.mode_expansion_monitor_name = monitor_name + &#8216;_mode_exp'&lt;br&gt;self.adjoint_source_name = monitor_name + &#8216;_mode_src'&lt;br&gt;self.mode_number = mode_number&lt;br&gt;self.target_fom = target_fom&lt;br&gt;# \u65b0\u589e\u53c2\u6570&#8230;&lt;br&gt;self.target_phase = np.radians(target_phase) # \u8f6c\u6362\u4e3a\u5f27\u5ea6&lt;br&gt;self.power_weight = power_weight&lt;br&gt;self.phase_weight = phase_weight&lt;br&gt;&lt;br&gt;if is_int(mode_number):&lt;br&gt;self.mode_number = int(mode_number)&lt;br&gt;if self.mode_number &lt;= 0:&lt;br&gt;raise UserWarning(&#8216;mode number should be positive.&#8217;)&lt;br&gt;else:&lt;br&gt;self.mode_number = mode_number&lt;br&gt;self.direction = str(direction)&lt;br&gt;self.multi_freq_src = bool(multi_freq_src)&lt;br&gt;if self.direction != &#8216;Forward&#8217; and self.direction != &#8216;Backward&#8217;:&lt;br&gt;raise UserWarning(&#8216;invalid propagation direction.&#8217;)&lt;br&gt;target_T_fwd_result = target_T_fwd(np.linspace(0.1e-6, 10.0e-6, 1000))&lt;br&gt;if target_T_fwd_result.size != 1000:&lt;br&gt;raise UserWarning(&#8216;target transmission must return a flat vector with the requested number of wavelength samples.&#8217;)&lt;br&gt;elif np.any(target_T_fwd_result.min() &lt; 0.0) or np.any(target_T_fwd_result.max() &gt; 1.0):&lt;br&gt;raise UserWarning(&#8216;target transmission must always return numbers between zero and one.&#8217;)&lt;br&gt;else:&lt;br&gt;self.target_T_fwd = target_T_fwd&lt;br&gt;self.norm_p = int(norm_p)&lt;br&gt;if self.norm_p &lt; 1:&lt;br&gt;raise UserWarning(&#8216;exponent p for norm must be positive.&#8217;)&lt;br&gt;&lt;br&gt;def initialize(self, sim):&lt;br&gt;self.check_monitor_alignment(sim)&lt;br&gt;&lt;br&gt;ModeMatch.add_mode_expansion_monitor(sim, self.monitor_name, self.mode_expansion_monitor_name, self.mode_number)&lt;br&gt;adjoint_injection_direction = &#8216;Backward&#8217; if self.direction == &#8216;Forward&#8217; else &#8216;Forward'&lt;br&gt;ModeMatch.add_mode_source(sim, self.monitor_name, self.adjoint_source_name, adjoint_injection_direction, self.mode_number, self.multi_freq_src)&lt;br&gt;&lt;br&gt;def make_forward_sim(self, sim):&lt;br&gt;sim.fdtd.setnamed(self.adjoint_source_name, &#8216;enabled&#8217;, False)&lt;br&gt;&lt;br&gt;def make_adjoint_sim(self, sim):&lt;br&gt;sim.fdtd.setnamed(self.adjoint_source_name, &#8216;enabled&#8217;, True)&lt;br&gt;&lt;br&gt;def check_monitor_alignment(self, sim):&lt;br&gt;&lt;br&gt;## Here, we check that the FOM_monitor is properly aligned with the mesh&lt;br&gt;if sim.fdtd.getnamednumber(self.monitor_name) != 1:&lt;br&gt;raise UserWarning(&#8216;monitor could not be found or the specified name is not unique.&#8217;)&lt;br&gt;&lt;br&gt;# Get the orientation&lt;br&gt;monitor_type = sim.fdtd.getnamed(self.monitor_name, &#8216;monitor type&#8217;)&lt;br&gt;&lt;br&gt;if (monitor_type == &#8216;Linear X&#8217;) or (monitor_type == &#8216;2D X-normal&#8217;):&lt;br&gt;orientation = &#8216;x'&lt;br&gt;elif (monitor_type == &#8216;Linear Y&#8217;) or (monitor_type == &#8216;2D Y-normal&#8217;):&lt;br&gt;orientation = &#8216;y'&lt;br&gt;elif (monitor_type == &#8216;Linear Z&#8217;) or (monitor_type == &#8216;2D Z-normal&#8217;):&lt;br&gt;orientation = &#8216;z'&lt;br&gt;else:&lt;br&gt;raise UserWarning(&#8216;monitor should be 2D or linear for a mode expansion to be meaningful.&#8217;)&lt;br&gt;&lt;br&gt;monitor_pos = sim.fdtd.getnamed(self.monitor_name, orientation)&lt;br&gt;if sim.fdtd.getnamednumber(&#8216;FDTD&#8217;) == 1:&lt;br&gt;grid = sim.fdtd.getresult(&#8216;FDTD&#8217;, orientation)&lt;br&gt;elif sim.fdtd.getnamednumber(&#8216;varFDTD&#8217;) == 1:&lt;br&gt;grid = sim.fdtd.getresult(&#8216;varFDTD&#8217;, orientation)&lt;br&gt;else:&lt;br&gt;raise UserWarning(&#8216;no FDTD or varFDTD solver object could be found.&#8217;)&lt;br&gt;## Check if this is exactly aligned with the simulation mesh. It exactly aligns if we find a point&lt;br&gt;## along the grid which is no more than &#8216;tol&#8217; away from the position&lt;br&gt;tol = 1e-9&lt;br&gt;dist_from_nearest_mesh_point = min(abs(grid-monitor_pos))&lt;br&gt;if dist_from_nearest_mesh_point &gt; tol:&lt;br&gt;print(&#8216;WARNING: The monitor &#8220;{}&#8221; is not aligned with the grid. Its distance to the nearest mesh point is {}. This can introduce small phase errors which sometimes result in inaccurate gradients.&#8217;.format(self.monitor_name,dist_from_nearest_mesh_point))&lt;br&gt;&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def add_mode_expansion_monitor(sim, monitor_name, mode_expansion_monitor_name, mode):&lt;br&gt;# modify existing DFT monitor&lt;br&gt;if sim.fdtd.getnamednumber(monitor_name) != 1:&lt;br&gt;raise UserWarning(&#8216;monitor could not be found or the specified name is not unique.&#8217;)&lt;br&gt;sim.fdtd.setnamed(monitor_name, &#8216;override global monitor settings&#8217;, False)&lt;br&gt;# append a mode expansion monitor to the existing DFT monitor&lt;br&gt;if sim.fdtd.getnamednumber(mode_expansion_monitor_name) == 0:&lt;br&gt;sim.fdtd.addmodeexpansion()&lt;br&gt;sim.fdtd.set(&#8216;name&#8217;, mode_expansion_monitor_name)&lt;br&gt;sim.fdtd.setexpansion(mode_expansion_monitor_name, monitor_name)&lt;br&gt;&lt;br&gt;sim.fdtd.setnamed(mode_expansion_monitor_name, &#8216;auto update before analysis&#8217;, True)&lt;br&gt;sim.fdtd.setnamed(mode_expansion_monitor_name, &#8216;override global monitor settings&#8217;, False)&lt;br&gt;# properties that must be synchronized&lt;br&gt;props = [&#8216;monitor type&#8217;]&lt;br&gt;monitor_type = sim.fdtd.getnamed(monitor_name, &#8216;monitor type&#8217;)&lt;br&gt;geo_props, normal = ModeMatch.cross_section_monitor_props(monitor_type)&lt;br&gt;props.extend(geo_props)&lt;br&gt;# synchronize properties&lt;br&gt;for prop_name in props:&lt;br&gt;prop_val = sim.fdtd.getnamed(monitor_name, prop_name)&lt;br&gt;sim.fdtd.setnamed(mode_expansion_monitor_name, prop_name, prop_val)&lt;br&gt;&lt;br&gt;# select mode&lt;br&gt;sim.fdtd.select(mode_expansion_monitor_name)&lt;br&gt;&lt;br&gt;## If &#8220;mode&#8221; is an integer, it means that we want a &#8216;user select&#8217; mode. Otherwise we try using it as a string directly&lt;br&gt;if is_int(mode):&lt;br&gt;sim.fdtd.setnamed(mode_expansion_monitor_name, &#8216;mode selection&#8217;, &#8216;user select&#8217;)&lt;br&gt;sim.fdtd.updatemodes(mode)&lt;br&gt;else:&lt;br&gt;sim.fdtd.setnamed(mode_expansion_monitor_name, &#8216;mode selection&#8217;, mode)&lt;br&gt;sim.fdtd.updatemodes()&lt;br&gt;else:&lt;br&gt;raise UserWarning(&#8216;there is already a expansion monitor with the same name.&#8217;)&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def cross_section_monitor_props(monitor_type):&lt;br&gt;geometric_props = [&#8216;x&#8217;, &#8216;y&#8217;, &#8216;z&#8217;]&lt;br&gt;normal = &#8221;&lt;br&gt;if monitor_type == &#8216;2D X-normal&#8217;:&lt;br&gt;geometric_props.extend([&#8216;y span&#8217;,&#8217;z span&#8217;])&lt;br&gt;normal = &#8216;x'&lt;br&gt;elif monitor_type == &#8216;2D Y-normal&#8217;:&lt;br&gt;geometric_props.extend([&#8216;x span&#8217;,&#8217;z span&#8217;])&lt;br&gt;normal = &#8216;y'&lt;br&gt;elif monitor_type == &#8216;2D Z-normal&#8217;:&lt;br&gt;geometric_props.extend([&#8216;x span&#8217;,&#8217;y span&#8217;])&lt;br&gt;normal = &#8216;z'&lt;br&gt;elif monitor_type == &#8216;Linear X&#8217;:&lt;br&gt;geometric_props.append(&#8216;x span&#8217;)&lt;br&gt;normal = &#8216;y'&lt;br&gt;elif monitor_type == &#8216;Linear Y&#8217;:&lt;br&gt;geometric_props.append(&#8216;y span&#8217;)&lt;br&gt;normal = &#8216;x'&lt;br&gt;elif monitor_type == &#8216;Linear Z&#8217;:&lt;br&gt;geometric_props.append(&#8216;z span&#8217;)&lt;br&gt;else:&lt;br&gt;raise UserWarning(&#8216;monitor should be 2D or linear for a mode expansion to be meaningful.&#8217;)&lt;br&gt;return geometric_props, normal&lt;br&gt;&lt;br&gt;def get_fom(self, sim):&lt;br&gt;trans_coeff, self.wavelengths = ModeMatch.get_transmission_coefficient(sim, self.direction, self.monitor_name, self.mode_expansion_monitor_name)&lt;br&gt;assert np.allclose(self.wavelengths, ModeMatch.get_wavelengths(sim))&lt;br&gt;source_power = ModeMatch.get_source_power(sim, self.wavelengths)&lt;br&gt;&lt;br&gt;# Calculate the actual power |T|&sup2; \/ P_source (normalized)&lt;br&gt;power_term = np.real(trans_coeff * trans_coeff.conj() \/ source_power)&lt;br&gt;power_error = (power_term &#8211; self.target_fom **2 ) ** 2 # target_power=target_fom&lt;br&gt;# self.T_fwd_vs_wavelength = np.real(trans_coeff * trans_coeff.conj() \/ source_power)&lt;br&gt;# self.phase_prefactors = trans_coeff \/ 4.0 \/ source_power&lt;br&gt;&lt;br&gt;# Calculate the phase error term (&ang;T &#8211; target_phase)^2 (mapped to [-&pi;, &pi;])&lt;br&gt;phase = np.angle(trans_coeff)&lt;br&gt;phase_diff = phase &#8211; self.target_phase&lt;br&gt;phase_diff = np.mod(phase_diff + np.pi, 2 * np.pi) &#8211; np.pi&lt;br&gt;phase_error = phase_diff ** 2&lt;br&gt;&lt;br&gt;# Weighted combination of the figure of merit (FOM) (adjust the sign according to the optimization direction)&lt;br&gt;self.T_fwd_vs_wavelength = (&lt;br&gt;self.power_weight * power_error +&lt;br&gt;self.phase_weight * phase_error&lt;br&gt;)&lt;br&gt;fom = ModeMatch.fom_wavelength_integral(self.T_fwd_vs_wavelength, self.wavelengths, self.target_T_fwd, self.norm_p)&lt;br&gt;return fom&lt;br&gt;&lt;br&gt;def get_adjoint_field_scaling(self, sim):&lt;br&gt;# Obtain the transmission coefficient and wavelength from the original code&lt;br&gt;trans_coeff, self.wavelengths = ModeMatch.get_transmission_coefficient(sim, self.direction, self.monitor_name,&lt;br&gt;self.mode_expansion_monitor_name)&lt;br&gt;assert np.allclose(self.wavelengths, ModeMatch.get_wavelengths(sim))&lt;br&gt;source_power = ModeMatch.get_source_power(sim, self.wavelengths)&lt;br&gt;# Power gradient term&lt;br&gt;# grad_power = 2 * trans_coeff.conj() # when the FOM=max(power)&lt;br&gt;grad_power = 2 * (np.real(trans_coeff * trans_coeff.conj() \/ source_power) &#8211; self.target_fom) * 2 * trans_coeff.conj()&lt;br&gt;&lt;br&gt;# Calculate the gradient of the phase term d( (&ang;T &#8211; &theta;)^2 )\/dT = 2(&ang;T &#8211; &theta;) * (j\/T*) and normalize it within the range of [-&pi;, &pi;]&lt;br&gt;phase_diff = np.angle(trans_coeff) &#8211; self.target_phase&lt;br&gt;phase_diff = np.mod(phase_diff + np.pi, 2 * np.pi) &#8211; np.pi&lt;br&gt;grad_phase = 2 * phase_diff * (1j \/ trans_coeff.conj())&lt;br&gt;&lt;br&gt;# Combine gradient terms&lt;br&gt;total_grad = 1\/3 * (&lt;br&gt;self.power_weight * grad_power +&lt;br&gt;self.phase_weight * grad_phase&lt;br&gt;)&lt;br&gt;&lt;br&gt;omega = 2.0 * np.pi * sp.constants.speed_of_light \/ self.wavelengths&lt;br&gt;adjoint_source_power = ModeMatch.get_source_power(sim, self.wavelengths)&lt;br&gt;&lt;br&gt;# Final scaling factor (verify the sign)&lt;br&gt;scaling_factor = total_grad * omega * 1j \/ np.sqrt(adjoint_source_power)&lt;br&gt;return scaling_factor&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def get_wavelengths(sim):&lt;br&gt;return Wavelengths(sim.fdtd.getglobalsource(&#8216;wavelength start&#8217;),&lt;br&gt;sim.fdtd.getglobalsource(&#8216;wavelength stop&#8217;),&lt;br&gt;sim.fdtd.getglobalmonitor(&#8216;frequency points&#8217;)).asarray()&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def get_transmission_coefficient(sim, direction, monitor_name, mode_exp_monitor_name):&lt;br&gt;mode_exp_result_name = &#8216;expansion for &#8216; + mode_exp_monitor_name&lt;br&gt;if not sim.fdtd.haveresult(mode_exp_monitor_name, mode_exp_result_name):&lt;br&gt;raise UserWarning(&#8216;unable to calcualte mode expansion.&#8217;)&lt;br&gt;mode_exp_data_set = sim.fdtd.getresult(mode_exp_monitor_name, mode_exp_result_name)&lt;br&gt;mode_exp_wl = mode_exp_data_set[&#8216;lambda&#8217;]&lt;br&gt;fwd_trans_coeff = mode_exp_data_set[&#8216;a&#8217;] * np.sqrt(mode_exp_data_set[&#8216;N&#8217;].real)&lt;br&gt;back_trans_coeff = mode_exp_data_set[&#8216;b&#8217;] * np.sqrt(mode_exp_data_set[&#8216;N&#8217;].real)&lt;br&gt;if direction == &#8216;Backward&#8217;:&lt;br&gt;fwd_trans_coeff, back_trans_coeff = back_trans_coeff, fwd_trans_coeff&lt;br&gt;return fwd_trans_coeff.flatten(), mode_exp_wl.flatten()&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def get_source_power(sim, wavelengths):&lt;br&gt;frequency = sp.constants.speed_of_light \/ wavelengths&lt;br&gt;source_power = sim.fdtd.sourcepower(frequency)&lt;br&gt;return np.asarray(source_power).flatten()&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def fom_wavelength_integral(T_fwd_vs_wavelength, wavelengths, target_T_fwd, norm_p):&lt;br&gt;target_T_fwd_vs_wavelength = target_T_fwd(wavelengths).flatten()&lt;br&gt;if len(wavelengths) &gt; 1:&lt;br&gt;wavelength_range = wavelengths.max() &#8211; wavelengths.min()&lt;br&gt;assert wavelength_range &gt; 0.0, &#8220;wavelength range must be positive.&#8221;&lt;br&gt;T_fwd_integrand = np.power(np.abs(target_T_fwd_vs_wavelength), norm_p) \/ wavelength_range&lt;br&gt;const_term = np.power(np.trapz(y = T_fwd_integrand, x = wavelengths), 1.0 \/ norm_p)&lt;br&gt;T_fwd_error = np.abs(T_fwd_vs_wavelength.flatten() &#8211; target_T_fwd_vs_wavelength)&lt;br&gt;T_fwd_error_integrand = np.power(T_fwd_error, norm_p) \/ wavelength_range&lt;br&gt;error_term = np.power(np.trapz(y = T_fwd_error_integrand, x = wavelengths), 1.0 \/ norm_p)&lt;br&gt;fom = const_term &#8211; error_term&lt;br&gt;else:&lt;br&gt;fom = np.abs(target_T_fwd_vs_wavelength) &#8211; np.abs(T_fwd_vs_wavelength.flatten() &#8211; target_T_fwd_vs_wavelength)&lt;br&gt;return fom.real&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def add_mode_source(sim, monitor_name, source_name, direction, mode, multi_freq_src):&lt;br&gt;if sim.fdtd.getnamednumber(&#8216;FDTD&#8217;) == 1:&lt;br&gt;sim.fdtd.addmode()&lt;br&gt;elif sim.fdtd.getnamednumber(&#8216;varFDTD&#8217;) == 1:&lt;br&gt;sim.fdtd.addmodesource()&lt;br&gt;else:&lt;br&gt;raise UserWarning(&#8216;no FDTD or varFDTD solver object could be found.&#8217;)&lt;br&gt;sim.fdtd.set(&#8216;name&#8217;, source_name)&lt;br&gt;sim.fdtd.select(source_name)&lt;br&gt;monitor_type = sim.fdtd.getnamed(monitor_name, &#8216;monitor type&#8217;)&lt;br&gt;geo_props, normal = ModeMatch.cross_section_monitor_props(monitor_type)&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;injection axis&#8217;, normal.lower() + &#8216;-axis&#8217;)&lt;br&gt;if sim.fdtd.getnamednumber(&#8216;varFDTD&#8217;) == 1:&lt;br&gt;geo_props.remove(&#8216;z&#8217;)&lt;br&gt;for prop_name in geo_props:&lt;br&gt;prop_val = sim.fdtd.getnamed(monitor_name, prop_name)&lt;br&gt;sim.fdtd.setnamed(source_name, prop_name, prop_val)&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;override global source settings&#8217;, False)&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;direction&#8217;, direction)&lt;br&gt;if sim.fdtd.haveproperty(&#8216;multifrequency mode calculation&#8217;):&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;multifrequency mode calculation&#8217;, multi_freq_src)&lt;br&gt;if multi_freq_src:&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;frequency points&#8217;, sim.fdtd.getglobalmonitor(&#8216;frequency points&#8217;))&lt;br&gt;&lt;br&gt;if is_int(mode):&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;mode selection&#8217;, &#8216;user select&#8217;)&lt;br&gt;sim.fdtd.select(source_name)&lt;br&gt;sim.fdtd.updatesourcemode(int(mode))&lt;br&gt;else:&lt;br&gt;sim.fdtd.setnamed(source_name, &#8216;mode selection&#8217;, mode)&lt;br&gt;sim.fdtd.select(source_name)&lt;br&gt;sim.fdtd.updatesourcemode()&lt;br&gt;&lt;br&gt;def fom_gradient_wavelength_integral(self, T_fwd_partial_derivs_vs_wl, wl):&lt;br&gt;assert np.allclose(wl, self.wavelengths)&lt;br&gt;return ModeMatch.fom_gradient_wavelength_integral_impl(self.T_fwd_vs_wavelength, T_fwd_partial_derivs_vs_wl, self.target_T_fwd(wl).flatten(), self.wavelengths, self.norm_p)&lt;br&gt;&lt;br&gt;def fom_gradient_wavelength_integral_on_cad(self, sim, grad_var_name, wl):&lt;br&gt;assert np.allclose(wl, self.wavelengths)&lt;br&gt;&lt;br&gt;target_T_fwd_vs_wavelength = self.target_T_fwd(wl).flatten()&lt;br&gt;T_fwd_error = self.T_fwd_vs_wavelength &#8211; target_T_fwd_vs_wavelength&lt;br&gt;&lt;br&gt;if wl.size &gt; 1:&lt;br&gt;wavelength_range = wl.max() &#8211; wl.min()&lt;br&gt;T_fwd_error_integrand = np.power(np.abs(T_fwd_error), self.norm_p) \/ wavelength_range&lt;br&gt;const_factor = -1.0 * np.power(np.trapz(y = T_fwd_error_integrand, x = wl), 1.0 \/ self.norm_p &#8211; 1.0)&lt;br&gt;integral_kernel = np.power(np.abs(T_fwd_error), self.norm_p &#8211; 1) * np.sign(T_fwd_error) \/ wavelength_range&lt;br&gt;&lt;br&gt;d = np.diff(wl)&lt;br&gt;quad_weight = np.append(np.append(d[0], d[0:-1]+d[1:]),d[-1])\/2 #&lt; There is probably a more elegant way to do this&lt;br&gt;v = const_factor * integral_kernel * quad_weight&lt;br&gt;&lt;br&gt;lumapi.putMatrix(sim.fdtd.handle, &#8220;wl_scaled_integral_kernel&#8221;, v)&lt;br&gt;sim.fdtd.eval((&#8216;dF_dp_s=size({0});'&lt;br&gt;&#8217;dF_dp2 = reshape(permute({0},[3,2,1]),[dF_dp_s(3),dF_dp_s(2)*dF_dp_s(1)]);'&lt;br&gt;&#8217;T_fwd_partial_derivs=real(mult(transpose(wl_scaled_integral_kernel),dF_dp2));&#8217;).format(grad_var_name) )&lt;br&gt;T_fwd_partial_derivs_on_cad = sim.fdtd.getv(&#8220;T_fwd_partial_derivs&#8221;)&lt;br&gt;&lt;br&gt;else:&lt;br&gt;sim.fdtd.eval((&#8216;T_fwd_partial_derivs=real({0});&#8217;).format(grad_var_name) )&lt;br&gt;T_fwd_partial_derivs_on_cad = sim.fdtd.getv(&#8220;T_fwd_partial_derivs&#8221;)&lt;br&gt;T_fwd_partial_derivs_on_cad*= -1.0 * np.sign(T_fwd_error)&lt;br&gt;&lt;br&gt;return T_fwd_partial_derivs_on_cad.flatten()&lt;br&gt;&lt;br&gt;&lt;br&gt;@staticmethod&lt;br&gt;def fom_gradient_wavelength_integral_impl(T_fwd_vs_wavelength, T_fwd_partial_derivs_vs_wl, target_T_fwd_vs_wavelength, wl, norm_p):&lt;br&gt;&lt;br&gt;if wl.size &gt; 1:&lt;br&gt;assert T_fwd_partial_derivs_vs_wl.shape[1] == wl.size&lt;br&gt;&lt;br&gt;wavelength_range = wl.max() &#8211; wl.min()&lt;br&gt;T_fwd_error = T_fwd_vs_wavelength &#8211; target_T_fwd_vs_wavelength&lt;br&gt;T_fwd_error_integrand = np.power(np.abs(T_fwd_error), norm_p) \/ wavelength_range&lt;br&gt;const_factor = -1.0 * np.power(np.trapz(y = T_fwd_error_integrand, x = wl), 1.0 \/ norm_p &#8211; 1.0)&lt;br&gt;integral_kernel = np.power(np.abs(T_fwd_error), norm_p &#8211; 1) * np.sign(T_fwd_error) \/ wavelength_range&lt;br&gt;&lt;br&gt;## Implement the trapezoidal integration as a matrix-vector-product for performance reasons&lt;br&gt;d = np.diff(wl)&lt;br&gt;quad_weight = np.append(np.append(d[0], d[0:-1]+d[1:]),d[-1])\/2 #&lt; There is probably a more elegant way to do this&lt;br&gt;v = const_factor * integral_kernel * quad_weight&lt;br&gt;T_fwd_partial_derivs = T_fwd_partial_derivs_vs_wl.dot(v)&lt;br&gt;&lt;br&gt;## This is the much slower (but possibly more readable) code&lt;br&gt;# num_opt_param = T_fwd_partial_derivs_vs_wl.shape[0]&lt;br&gt;# T_fwd_partial_derivs = np.zeros(num_opt_param, dtype = &#8216;complex&#8217;)&lt;br&gt;# for i in range(num_opt_param):&lt;br&gt;# T_fwd_partial_deriv = np.take(T_fwd_partial_derivs_vs_wl.transpose(), indices = i, axis = 1)&lt;br&gt;# T_fwd_partial_derivs[i] = const_factor * np.trapz(y = integral_kernel * T_fwd_partial_deriv, x = wl)&lt;br&gt;else:&lt;br&gt;T_fwd_partial_derivs = -1.0 * np.sign(T_fwd_vs_wavelength &#8211; target_T_fwd_vs_wavelength) * T_fwd_partial_derivs_vs_wl.flatten()&lt;br&gt;&lt;br&gt;return T_fwd_partial_derivs.flatten().real&lt;br&gt;&lt;br&gt;&lt;br&gt;&lt;br&gt;The code for calling the modematch function is as follows:&lt;br&gt;geometry = TopologyOptimization2D(params=params, eps_min=eps_bg, eps_max=eps_wg, x=x_pos, y=y_pos, z=0, filter_R=filter_R, beta=beta)&lt;br&gt;&lt;br&gt;######## DEFINE FIGURE OF MERIT FOR EACH OUTPUT WAVEGUIDE ########&lt;br&gt;fom1 = ModeMatch(&lt;br&gt;monitor_name=&#8217;fom_1&#8242;, mode_number=&#8217;Fundamental TE mode&#8217;, direction=&#8217;Forward&#8217;,&lt;br&gt;target_phase=0, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1\/3 * np.ones(wl.size),&lt;br&gt;norm_p=2, target_fom=1\/3)&lt;br&gt;fom2 = ModeMatch(&lt;br&gt;monitor_name=&#8217;fom_2&#8242;, mode_number=&#8217;Fundamental TE mode&#8217;, direction=&#8217;Forward&#8217;,&lt;br&gt;target_phase=120, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1\/3 * np.ones(wl.size),&lt;br&gt;norm_p=2, target_fom=1\/3)&lt;br&gt;fom3 = ModeMatch(&lt;br&gt;monitor_name=&#8217;fom_3&#8242;, mode_number=&#8217;Fundamental TE mode&#8217;, direction=&#8217;Forward&#8217;,&lt;br&gt;target_phase=240, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1\/3 * np.ones(wl.size),&lt;br&gt;norm_p=2, target_fom=1\/3)&lt;br&gt;&lt;br&gt;######## DEFINE OPTIMIZATION ALGORITHM ########&lt;br&gt;optimizer = ScipyOptimizers(max_iter=400, method=&#8217;L-BFGS-B&#8217;, pgtol=1e-6, ftol=1e-5, scale_initial_gradient_to=0.25)&lt;br&gt;&lt;br&gt;######## DEFINE SETUP SCRIPT AND INDIVIDUAL OPTIMIZERS ########&lt;br&gt;script1 = load_from_lsf(&#8216;construct 2&#215;3 MMI 6&#215;6 rect input1.lsf&#8217;)&lt;br&gt;&lt;br&gt;&lt;br&gt;wavelengths1 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)&lt;br&gt;opt1 = Optimization(base_script=script1, wavelengths = wavelengths1, fom=fom1, geometry=geometry, optimizer=optimizer, use_deps=False, hide_fdtd_cad=False, plot_history=False, store_all_simulations=False, save_global_index=True)&lt;br&gt;wavelengths2 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)&lt;br&gt;opt2 = Optimization(base_script=script1, wavelengths = wavelengths2, fom=fom2, geometry=geometry, optimizer=optimizer, use_deps=False, hide_fdtd_cad=False, plot_history=False, store_all_simulations=False, save_global_index=False)&lt;br&gt;wavelengths3 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)&lt;br&gt;opt3 = Optimization(base_script=script1, wavelengths = wavelengths3, fom=fom3, geometry=geometry, optimizer=optimizer, use_deps=False, hide_fdtd_cad=False, plot_history=False, store_all_simulations=False, save_global_index=False)&lt;br&gt;wavelengths1 = Wavelengths(start=1550e-9, stop=1550e-9, points=1)&lt;br&gt;&lt;br&gt;&lt;br&gt;######## PUT EVERYTHING TOGETHER AND RUN ########&lt;br&gt;opt = opt1+opt2+opt3&lt;br&gt;opt.run(working_dir = working_dir)&lt;\/p&gt;<\/p>\n","protected":false},"template":"","class_list":["post-427139","topic","type-topic","status-closed","hentry","topic-tag-lumopt-1","topic-tag-mmi-2","topic-tag-phase-1"],"aioseo_notices":[],"acf":[],"custom_fields":[{"0":{"_bbp_forum_id":["27833"],"_bbp_topic_id":["427139"],"_bbp_author_ip":["2408:8452:2af:d311:7922:2b94:fe68:28be"],"_bbp_last_reply_id":["427204"],"_bbp_last_active_id":["427204"],"_bbp_last_active_time":["2025-03-28 12:57:35"],"_bbp_reply_count":["1"],"_bbp_reply_count_hidden":["0"],"_bbp_voice_count":["2"],"_bbp_engagement":["520400","473055"],"_btv_view_count":["519"],"_bbp_topic_status":["unanswered"],"_bbp_notification_enabled":["473055"],"_bbp_subscription":["473055"],"_bbp_status":["publish"]},"test":"longchdtu-dk"}],"_links":{"self":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/427139","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics"}],"about":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/types\/topic"}],"version-history":[{"count":0,"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/topics\/427139\/revisions"}],"wp:attachment":[{"href":"https:\/\/innovationspace.ansys.com\/forum\/wp-json\/wp\/v2\/media?parent=427139"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}