Learning Forum

[longch]

We modified the modematch code and added phase control to the optimization objective to achieve the inverse design of a 2x3 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°. 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!

code of modematch：
""" Copyright chriskeraly
Copyright (c) 2019 Lumerical Inc. """

import sys
import numpy as np
import scipy as sp
import scipy.constants
import lumapi

from lumopt.utilities.wavelengths import Wavelengths

""" Helper function to determine if a variable can be converted to an integer """
def is_int(x):
try:
int(x)
return True
except ValueError:
return False

class ModeMatch(object):

""" Calculates the figure of merit from an overlap integral between the fields recorded by a field monitor and the slected mode.
A mode expansion monitor is added to the field monitor to calculate the overlap result, which appears as T_forward in the
list of mode expansion monitor results. The T_forward result is described in the following page:

https://kb.lumerical.com/ref_sim_obj_using_mode_expansion_monitors.html

This result is equivalent to equation (7) in the following paper:

C. Lalau-Keraly, S. Bhargava, O. Miller, and E. Yablonovitch, "Adjoint shape optimization applied to electromagnetic design,"
Opt. Express 21, 21693-21701 (2013). https://doi.org/10.1364/OE.21.021693

Parameters
----------
:param monitor_name: name of the field monitor that records the fields to be used in the mode overlap calculation.
:param mode_number: mode number in the list of modes generated by the mode expansion monitor.
:param direction: direction of propagation ('Forward' or 'Backward') of the mode injected by the source.
:param multi_freq_src: bool flag to enable / disable a multi-frequency mode calculation and injection for the adjoint source.
:param target_T_fwd: function describing the target T_forward vs wavelength (see documentation for mode expansion monitors).
:param norm_p: exponent of the p-norm used to generate the figure of merit; use to generate the FOM.
:param target_fom: A target value for the figure of merit. This allows to print/plot the distance of the current
design from a target value
"""

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):
self.monitor_name = str(monitor_name)
if not self.monitor_name:
raise UserWarning('empty monitor name.')
self.mode_expansion_monitor_name = monitor_name + '_mode_exp'
self.adjoint_source_name = monitor_name + '_mode_src'
self.mode_number = mode_number
self.target_fom = target_fom
# 新增参数...
self.target_phase = np.radians(target_phase) # 转换为弧度
self.power_weight = power_weight
self.phase_weight = phase_weight

if is_int(mode_number):
self.mode_number = int(mode_number)
if self.mode_number <= 0:
raise UserWarning('mode number should be positive.')
else:
self.mode_number = mode_number
self.direction = str(direction)
self.multi_freq_src = bool(multi_freq_src)
if self.direction != 'Forward' and self.direction != 'Backward':
raise UserWarning('invalid propagation direction.')
target_T_fwd_result = target_T_fwd(np.linspace(0.1e-6, 10.0e-6, 1000))
if target_T_fwd_result.size != 1000:
raise UserWarning('target transmission must return a flat vector with the requested number of wavelength samples.')
elif np.any(target_T_fwd_result.min() < 0.0) or np.any(target_T_fwd_result.max() > 1.0):
raise UserWarning('target transmission must always return numbers between zero and one.')
else:
self.target_T_fwd = target_T_fwd
self.norm_p = int(norm_p)
if self.norm_p < 1:
raise UserWarning('exponent p for norm must be positive.')

def initialize(self, sim):
self.check_monitor_alignment(sim)

ModeMatch.add_mode_expansion_monitor(sim, self.monitor_name, self.mode_expansion_monitor_name, self.mode_number)
adjoint_injection_direction = 'Backward' if self.direction == 'Forward' else 'Forward'
ModeMatch.add_mode_source(sim, self.monitor_name, self.adjoint_source_name, adjoint_injection_direction, self.mode_number, self.multi_freq_src)

def make_forward_sim(self, sim):
sim.fdtd.setnamed(self.adjoint_source_name, 'enabled', False)

def make_adjoint_sim(self, sim):
sim.fdtd.setnamed(self.adjoint_source_name, 'enabled', True)

def check_monitor_alignment(self, sim):

## Here, we check that the FOM_monitor is properly aligned with the mesh
if sim.fdtd.getnamednumber(self.monitor_name) != 1:
raise UserWarning('monitor could not be found or the specified name is not unique.')

# Get the orientation
monitor_type = sim.fdtd.getnamed(self.monitor_name, 'monitor type')

if (monitor_type == 'Linear X') or (monitor_type == '2D X-normal'):
orientation = 'x'
elif (monitor_type == 'Linear Y') or (monitor_type == '2D Y-normal'):
orientation = 'y'
elif (monitor_type == 'Linear Z') or (monitor_type == '2D Z-normal'):
orientation = 'z'
else:
raise UserWarning('monitor should be 2D or linear for a mode expansion to be meaningful.')

monitor_pos = sim.fdtd.getnamed(self.monitor_name, orientation)
if sim.fdtd.getnamednumber('FDTD') == 1:
grid = sim.fdtd.getresult('FDTD', orientation)
elif sim.fdtd.getnamednumber('varFDTD') == 1:
grid = sim.fdtd.getresult('varFDTD', orientation)
else:
raise UserWarning('no FDTD or varFDTD solver object could be found.')
## Check if this is exactly aligned with the simulation mesh. It exactly aligns if we find a point
## along the grid which is no more than 'tol' away from the position
tol = 1e-9
dist_from_nearest_mesh_point = min(abs(grid-monitor_pos))
if dist_from_nearest_mesh_point > tol:
print('WARNING: The monitor "{}" 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.'.format(self.monitor_name,dist_from_nearest_mesh_point))


@staticmethod
def add_mode_expansion_monitor(sim, monitor_name, mode_expansion_monitor_name, mode):
# modify existing DFT monitor
if sim.fdtd.getnamednumber(monitor_name) != 1:
raise UserWarning('monitor could not be found or the specified name is not unique.')
sim.fdtd.setnamed(monitor_name, 'override global monitor settings', False)
# append a mode expansion monitor to the existing DFT monitor
if sim.fdtd.getnamednumber(mode_expansion_monitor_name) == 0:
sim.fdtd.addmodeexpansion()
sim.fdtd.set('name', mode_expansion_monitor_name)
sim.fdtd.setexpansion(mode_expansion_monitor_name, monitor_name)

sim.fdtd.setnamed(mode_expansion_monitor_name, 'auto update before analysis', True)
sim.fdtd.setnamed(mode_expansion_monitor_name, 'override global monitor settings', False)
# properties that must be synchronized
props = ['monitor type']
monitor_type = sim.fdtd.getnamed(monitor_name, 'monitor type')
geo_props, normal = ModeMatch.cross_section_monitor_props(monitor_type)
props.extend(geo_props)
# synchronize properties
for prop_name in props:
prop_val = sim.fdtd.getnamed(monitor_name, prop_name)
sim.fdtd.setnamed(mode_expansion_monitor_name, prop_name, prop_val)

# select mode
sim.fdtd.select(mode_expansion_monitor_name)

## If "mode" is an integer, it means that we want a 'user select' mode. Otherwise we try using it as a string directly
if is_int(mode):
sim.fdtd.setnamed(mode_expansion_monitor_name, 'mode selection', 'user select')
sim.fdtd.updatemodes(mode)
else:
sim.fdtd.setnamed(mode_expansion_monitor_name, 'mode selection', mode)
sim.fdtd.updatemodes()
else:
raise UserWarning('there is already a expansion monitor with the same name.')

@staticmethod
def cross_section_monitor_props(monitor_type):
geometric_props = ['x', 'y', 'z']
normal = ''
if monitor_type == '2D X-normal':
geometric_props.extend(['y span','z span'])
normal = 'x'
elif monitor_type == '2D Y-normal':
geometric_props.extend(['x span','z span'])
normal = 'y'
elif monitor_type == '2D Z-normal':
geometric_props.extend(['x span','y span'])
normal = 'z'
elif monitor_type == 'Linear X':
geometric_props.append('x span')
normal = 'y'
elif monitor_type == 'Linear Y':
geometric_props.append('y span')
normal = 'x'
elif monitor_type == 'Linear Z':
geometric_props.append('z span')
else:
raise UserWarning('monitor should be 2D or linear for a mode expansion to be meaningful.')
return geometric_props, normal

def get_fom(self, sim):
trans_coeff, self.wavelengths = ModeMatch.get_transmission_coefficient(sim, self.direction, self.monitor_name, self.mode_expansion_monitor_name)
assert np.allclose(self.wavelengths, ModeMatch.get_wavelengths(sim))
source_power = ModeMatch.get_source_power(sim, self.wavelengths)

# Calculate the actual power |T|² / P_source (normalized)
power_term = np.real(trans_coeff * trans_coeff.conj() / source_power)
power_error = (power_term - self.target_fom **2 ) ** 2 # target_power=target_fom
# self.T_fwd_vs_wavelength = np.real(trans_coeff * trans_coeff.conj() / source_power)
# self.phase_prefactors = trans_coeff / 4.0 / source_power

# Calculate the phase error term (∠T - target_phase)^2 (mapped to [-π, π])
phase = np.angle(trans_coeff)
phase_diff = phase - self.target_phase
phase_diff = np.mod(phase_diff + np.pi, 2 * np.pi) - np.pi
phase_error = phase_diff ** 2

# Weighted combination of the figure of merit (FOM) (adjust the sign according to the optimization direction)
self.T_fwd_vs_wavelength = (
self.power_weight * power_error +
self.phase_weight * phase_error
)
fom = ModeMatch.fom_wavelength_integral(self.T_fwd_vs_wavelength, self.wavelengths, self.target_T_fwd, self.norm_p)
return fom

def get_adjoint_field_scaling(self, sim):
# Obtain the transmission coefficient and wavelength from the original code
trans_coeff, self.wavelengths = ModeMatch.get_transmission_coefficient(sim, self.direction, self.monitor_name,
self.mode_expansion_monitor_name)
assert np.allclose(self.wavelengths, ModeMatch.get_wavelengths(sim))
source_power = ModeMatch.get_source_power(sim, self.wavelengths)
# Power gradient term
# grad_power = 2 * trans_coeff.conj() # when the FOM=max(power)
grad_power = 2 * (np.real(trans_coeff * trans_coeff.conj() / source_power) - self.target_fom) * 2 * trans_coeff.conj()

# Calculate the gradient of the phase term d( (∠T - θ)^2 )/dT = 2(∠T - θ) * (j/T*) and normalize it within the range of [-π, π]
phase_diff = np.angle(trans_coeff) - self.target_phase
phase_diff = np.mod(phase_diff + np.pi, 2 * np.pi) - np.pi
grad_phase = 2 * phase_diff * (1j / trans_coeff.conj())

# Combine gradient terms
total_grad = 1/3 * (
self.power_weight * grad_power +
self.phase_weight * grad_phase
)

omega = 2.0 * np.pi * sp.constants.speed_of_light / self.wavelengths
adjoint_source_power = ModeMatch.get_source_power(sim, self.wavelengths)

# Final scaling factor (verify the sign)
scaling_factor = total_grad * omega * 1j / np.sqrt(adjoint_source_power)
return scaling_factor

@staticmethod
def get_wavelengths(sim):
return Wavelengths(sim.fdtd.getglobalsource('wavelength start'),
sim.fdtd.getglobalsource('wavelength stop'),
sim.fdtd.getglobalmonitor('frequency points')).asarray()

@staticmethod
def get_transmission_coefficient(sim, direction, monitor_name, mode_exp_monitor_name):
mode_exp_result_name = 'expansion for ' + mode_exp_monitor_name
if not sim.fdtd.haveresult(mode_exp_monitor_name, mode_exp_result_name):
raise UserWarning('unable to calcualte mode expansion.')
mode_exp_data_set = sim.fdtd.getresult(mode_exp_monitor_name, mode_exp_result_name)
mode_exp_wl = mode_exp_data_set['lambda']
fwd_trans_coeff = mode_exp_data_set['a'] * np.sqrt(mode_exp_data_set['N'].real)
back_trans_coeff = mode_exp_data_set['b'] * np.sqrt(mode_exp_data_set['N'].real)
if direction == 'Backward':
fwd_trans_coeff, back_trans_coeff = back_trans_coeff, fwd_trans_coeff
return fwd_trans_coeff.flatten(), mode_exp_wl.flatten()

@staticmethod
def get_source_power(sim, wavelengths):
frequency = sp.constants.speed_of_light / wavelengths
source_power = sim.fdtd.sourcepower(frequency)
return np.asarray(source_power).flatten()

@staticmethod
def fom_wavelength_integral(T_fwd_vs_wavelength, wavelengths, target_T_fwd, norm_p):
target_T_fwd_vs_wavelength = target_T_fwd(wavelengths).flatten()
if len(wavelengths) > 1:
wavelength_range = wavelengths.max() - wavelengths.min()
assert wavelength_range > 0.0, "wavelength range must be positive."
T_fwd_integrand = np.power(np.abs(target_T_fwd_vs_wavelength), norm_p) / wavelength_range
const_term = np.power(np.trapz(y = T_fwd_integrand, x = wavelengths), 1.0 / norm_p)
T_fwd_error = np.abs(T_fwd_vs_wavelength.flatten() - target_T_fwd_vs_wavelength)
T_fwd_error_integrand = np.power(T_fwd_error, norm_p) / wavelength_range
error_term = np.power(np.trapz(y = T_fwd_error_integrand, x = wavelengths), 1.0 / norm_p)
fom = const_term - error_term
else:
fom = np.abs(target_T_fwd_vs_wavelength) - np.abs(T_fwd_vs_wavelength.flatten() - target_T_fwd_vs_wavelength)
return fom.real

@staticmethod
def add_mode_source(sim, monitor_name, source_name, direction, mode, multi_freq_src):
if sim.fdtd.getnamednumber('FDTD') == 1:
sim.fdtd.addmode()
elif sim.fdtd.getnamednumber('varFDTD') == 1:
sim.fdtd.addmodesource()
else:
raise UserWarning('no FDTD or varFDTD solver object could be found.')
sim.fdtd.set('name', source_name)
sim.fdtd.select(source_name)
monitor_type = sim.fdtd.getnamed(monitor_name, 'monitor type')
geo_props, normal = ModeMatch.cross_section_monitor_props(monitor_type)
sim.fdtd.setnamed(source_name, 'injection axis', normal.lower() + '-axis')
if sim.fdtd.getnamednumber('varFDTD') == 1:
geo_props.remove('z')
for prop_name in geo_props:
prop_val = sim.fdtd.getnamed(monitor_name, prop_name)
sim.fdtd.setnamed(source_name, prop_name, prop_val)
sim.fdtd.setnamed(source_name, 'override global source settings', False)
sim.fdtd.setnamed(source_name, 'direction', direction)
if sim.fdtd.haveproperty('multifrequency mode calculation'):
sim.fdtd.setnamed(source_name, 'multifrequency mode calculation', multi_freq_src)
if multi_freq_src:
sim.fdtd.setnamed(source_name, 'frequency points', sim.fdtd.getglobalmonitor('frequency points'))

if is_int(mode):
sim.fdtd.setnamed(source_name, 'mode selection', 'user select')
sim.fdtd.select(source_name)
sim.fdtd.updatesourcemode(int(mode))
else:
sim.fdtd.setnamed(source_name, 'mode selection', mode)
sim.fdtd.select(source_name)
sim.fdtd.updatesourcemode()

def fom_gradient_wavelength_integral(self, T_fwd_partial_derivs_vs_wl, wl):
assert np.allclose(wl, self.wavelengths)
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)

def fom_gradient_wavelength_integral_on_cad(self, sim, grad_var_name, wl):
assert np.allclose(wl, self.wavelengths)

target_T_fwd_vs_wavelength = self.target_T_fwd(wl).flatten()
T_fwd_error = self.T_fwd_vs_wavelength - target_T_fwd_vs_wavelength

if wl.size > 1:
wavelength_range = wl.max() - wl.min()
T_fwd_error_integrand = np.power(np.abs(T_fwd_error), self.norm_p) / wavelength_range
const_factor = -1.0 * np.power(np.trapz(y = T_fwd_error_integrand, x = wl), 1.0 / self.norm_p - 1.0)
integral_kernel = np.power(np.abs(T_fwd_error), self.norm_p - 1) * np.sign(T_fwd_error) / wavelength_range

d = np.diff(wl)
quad_weight = np.append(np.append(d[0], d[0:-1]+d[1:]),d[-1])/2 #< There is probably a more elegant way to do this
v = const_factor * integral_kernel * quad_weight

lumapi.putMatrix(sim.fdtd.handle, "wl_scaled_integral_kernel", v)
sim.fdtd.eval(('dF_dp_s=size({0});'
'dF_dp2 = reshape(permute({0},[3,2,1]),[dF_dp_s(3),dF_dp_s(2)*dF_dp_s(1)]);'
'T_fwd_partial_derivs=real(mult(transpose(wl_scaled_integral_kernel),dF_dp2));').format(grad_var_name) )
T_fwd_partial_derivs_on_cad = sim.fdtd.getv("T_fwd_partial_derivs")

else:
sim.fdtd.eval(('T_fwd_partial_derivs=real({0});').format(grad_var_name) )
T_fwd_partial_derivs_on_cad = sim.fdtd.getv("T_fwd_partial_derivs")
T_fwd_partial_derivs_on_cad*= -1.0 * np.sign(T_fwd_error)

return T_fwd_partial_derivs_on_cad.flatten()


@staticmethod
def fom_gradient_wavelength_integral_impl(T_fwd_vs_wavelength, T_fwd_partial_derivs_vs_wl, target_T_fwd_vs_wavelength, wl, norm_p):

if wl.size > 1:
assert T_fwd_partial_derivs_vs_wl.shape[1] == wl.size

wavelength_range = wl.max() - wl.min()
T_fwd_error = T_fwd_vs_wavelength - target_T_fwd_vs_wavelength
T_fwd_error_integrand = np.power(np.abs(T_fwd_error), norm_p) / wavelength_range
const_factor = -1.0 * np.power(np.trapz(y = T_fwd_error_integrand, x = wl), 1.0 / norm_p - 1.0)
integral_kernel = np.power(np.abs(T_fwd_error), norm_p - 1) * np.sign(T_fwd_error) / wavelength_range

## Implement the trapezoidal integration as a matrix-vector-product for performance reasons
d = np.diff(wl)
quad_weight = np.append(np.append(d[0], d[0:-1]+d[1:]),d[-1])/2 #< There is probably a more elegant way to do this
v = const_factor * integral_kernel * quad_weight
T_fwd_partial_derivs = T_fwd_partial_derivs_vs_wl.dot(v)

## This is the much slower (but possibly more readable) code
# num_opt_param = T_fwd_partial_derivs_vs_wl.shape[0]
# T_fwd_partial_derivs = np.zeros(num_opt_param, dtype = 'complex')
# for i in range(num_opt_param):
# T_fwd_partial_deriv = np.take(T_fwd_partial_derivs_vs_wl.transpose(), indices = i, axis = 1)
# T_fwd_partial_derivs[i] = const_factor * np.trapz(y = integral_kernel * T_fwd_partial_deriv, x = wl)
else:
T_fwd_partial_derivs = -1.0 * np.sign(T_fwd_vs_wavelength - target_T_fwd_vs_wavelength) * T_fwd_partial_derivs_vs_wl.flatten()

return T_fwd_partial_derivs.flatten().real



The code for calling the modematch function is as follows:
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)

######## DEFINE FIGURE OF MERIT FOR EACH OUTPUT WAVEGUIDE ########
fom1 = ModeMatch(
monitor_name='fom_1', mode_number='Fundamental TE mode', direction='Forward',
target_phase=0, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1/3 * np.ones(wl.size),
norm_p=2, target_fom=1/3)
fom2 = ModeMatch(
monitor_name='fom_2', mode_number='Fundamental TE mode', direction='Forward',
target_phase=120, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1/3 * np.ones(wl.size),
norm_p=2, target_fom=1/3)
fom3 = ModeMatch(
monitor_name='fom_3', mode_number='Fundamental TE mode', direction='Forward',
target_phase=240, power_weight=0.5, phase_weight=0.5, target_T_fwd = lambda wl: 1/3 * np.ones(wl.size),
norm_p=2, target_fom=1/3)

######## DEFINE OPTIMIZATION ALGORITHM ########
optimizer = ScipyOptimizers(max_iter=400, method='L-BFGS-B', pgtol=1e-6, ftol=1e-5, scale_initial_gradient_to=0.25)

######## DEFINE SETUP SCRIPT AND INDIVIDUAL OPTIMIZERS ########
script1 = load_from_lsf('construct 2x3 MMI 6x6 rect input1.lsf')


wavelengths1 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)
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)
wavelengths2 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)
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)
wavelengths3 = Wavelengths(start = 1550e-9, stop = 1550e-9, points = 1)
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)
wavelengths1 = Wavelengths(start=1550e-9, stop=1550e-9, points=1)


######## PUT EVERYTHING TOGETHER AND RUN ########
opt = opt1+opt2+opt3
opt.run(working_dir = working_dir)

[Kirill]

This topic has been moved to the Photonics forum. Thank you!

[longch]

Hello, Engineer Kirill. According to your suggestion, I reposted this post under the photonics topic, but you closed it. However, this post hasn't received any replies yet. What exactly should I do to get a response?

[George Patsamanis]

Hello,

thank you for reaching out to us.

For the 2x3 MMI you can find some feedback here:
Abnormal gradients in multi – port power optimization of 2×3 MMI using Lumopt

Regarding the phase optimization, currently we do not have a figure of merit for optimizing a geometry based on an output target phase. See "Set-up Classes" sections in here:
Getting Started with lumopt - Python API – Ansys Optics

However, once you have your geometry optimized based on power of modes, you can tune the phase of the output waveguides in separate simulations.

By any chance, have you used a working template that optimizes the final geometry based on ouput phases? Or did you setup a new parameter for the target phase? 

Regards,
George

[longch]

Thank you very much for your reply, Engineer George. However, I still want to consult you on some more questions.

 

1. Based on your reply, can I understand it in this way: I can first optimize the power equalization to obtain the device structure, and then adjust the phases of each output port in the generated structure? If so, could you provide specific adjustment methods or cases?
2. As shown in the code I provided, in order to achieve phase optimization, I added parameters such as target_phase, phase_weight, and power_weight, and modified modematch, incorporating phase optimization and gradient calculation. However, the gradient is abnormal. Are there still some places that haven't been modified properly? Or is there something wrong with my settings? Should I optimize the phase when two input light sources are incident simultaneously? Does this situation require a complete modification of the framework of the lumopt optimization algorithm? Do you have any better suggestions?

[George Patsamanis]

Hello,

thank you for the feedback.

1) Yes, you can first optimize the power splitting using the inverse design method. The phase adjustment can be performed in later step if you wish to do it passively. For example, you can add small S-bent waveguide structures along one arm, that will act as a delay line (zig-zag). However, in practice it is good to take into consideration fabrication errors that may occur. Optical phase may not be identical to the simulations, since the waveguide structures are simulated ideal most of the times. Usually, deployed active phase shifters ensure that the phase is corrected and properly adjusted.

2) It is most possible that these parameters regarding phase cannot be identified, since they are not defined in lumopt. This can be one explanation why the gradient is abnormal.

Optimizing the structure with 2 input sources at the same time will result in different geometry compared to when optimizing using only 1 input source. This is because the 2 modes will interfere along their way to the outputs. That said, if you optimize the structure for 2 input modes simultaneously, and then test the final structure with only one input mode, the power distribution will not be the same at the outputs because the optimization conditions were different. Varying the initial phase difference between the 2 input mode sources will also result in a different optimized geometry, since the field profile will be different due to interference.
Designing the device with 1 or 2 input modes at the same time depends on your application.

My suggestion is to start from scratch, having this example as template.
Inverse Design of a Splitter Using Topology Optimization – Ansys Optics

After getting familiar with the scripts, you can do your modifications one by one step by step to create your geometry and assign the foms and opts correctly. Make sure to remove the "y min bc" symmetry, as you don't need it (delete line 76 in the "splitter_base_2D_TE_topology.lsf).
(at the moment, there is an ongoing repair regarding the GDS export in step 3, but the geometry optimization is performed ok)

I hope this is helpful for your studies.

Regards,
George