From ed3e6b412f3dc4fb4b72019ececb8549e33c146c Mon Sep 17 00:00:00 2001 From: Denis Langevin <68906630+Milloupe@users.noreply.github.com> Date: Wed, 12 Feb 2025 15:57:26 +0100 Subject: [PATCH 1/5] Create ceviche_loss_fn.py --- .../functions/photonics/ceviche_loss_fn.py | 117 ++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 nevergrad/functions/photonics/ceviche_loss_fn.py diff --git a/nevergrad/functions/photonics/ceviche_loss_fn.py b/nevergrad/functions/photonics/ceviche_loss_fn.py new file mode 100644 index 000000000..e7f86312d --- /dev/null +++ b/nevergrad/functions/photonics/ceviche_loss_fn.py @@ -0,0 +1,117 @@ +import numpy as np +import typing as tp +import mode_converter + + +first_time_ceviche = True + +def ceviche( + x: np.ndarray, benchmark_type: int = 0, discretize=False, wantgrad=False, wantfields=False +) -> tp.Any: + """ + x = 2d or 3d array of scalars + Inputs: + 1. benchmark_type = 0 1 2 or 3, depending on which benchmark you want + 2. discretize = True if we want to force x to be in {0,1} (just checks sign(x-0.5) ) + 3. wantgrad = True if we want to know the gradient + 4. wantfields = True if we want the fields of the simulation + Returns: + 1. the loss (to be minimized) + 2. the gradient or none (depending on wantgrad) + 3. the fields or none (depending on wantfields + """ + global first_time_ceviche + global model + import autograd # type: ignore + import autograd.numpy as npa # type: ignore + import ceviche_challenges # type: ignore + import autograd # type: ignore + + # ceviche_challenges.beam_splitter.prefabs + # ceviche_challenges.mode_converter.prefabs + # ceviche_challenges.waveguide_bend.prefabs + # ceviche_challenges.wdm.prefabs + + if first_time_ceviche or x is None: + if benchmark_type == 2: + print("Mode converter benchmark") + design_variable_shape = (345,345) + simulate = mode_converter.mode_converter + else: + print("Benchmark not defined yet") + return None + # elif benchmark_type == 0: + # spec = ceviche_challenges.waveguide_bend.prefabs.waveguide_bend_2umx2um_spec() + # params = ceviche_challenges.waveguide_bend.prefabs.waveguide_bend_sim_params() + # model = ceviche_challenges.waveguide_bend.model.WaveguideBendModel(params, spec) + # elif benchmark_type == 1: + # spec = ceviche_challenges.beam_splitter.prefabs.pico_splitter_spec() + # params = ceviche_challenges.beam_splitter.prefabs.pico_splitter_sim_params() + # model = ceviche_challenges.beam_splitter.model.BeamSplitterModel(params, spec) + # elif benchmark_type == 3: + # spec = ceviche_challenges.wdm.prefabs.wdm_spec() + # params = ceviche_challenges.wdm.prefabs.wdm_sim_params() + # model = ceviche_challenges.wdm.model.WdmModel(params, spec) + if discretize: + # x between 0 and 1 shifted to either 0 or 1 + x = np.round(x) + design = x + + # if isinstance(x, str) and x == "name": + # return {0: "waveguide-bend", 1: "beam-splitter", 2: "mode-converter", 3: "wdm"}[benchmark_type] + # elif x is None: + # return model.design_variable_shape + + assert x.shape == design_variable_shape, f"Expected shape: {design_variable_shape}" # type: ignore + + # The model class provides a convenience property, `design_variable_shape` + # which specifies the design shape it expects. + + # Construct a loss function, assuming the `model` and `design` from the code + # snippet above are instantiated. + + def loss_fn(x, fields_are_needed=False): + """A simple loss function""" + the_loss, field = simulate(x) + + if fields_are_needed: + return the_loss, field + return the_loss + + def loss_fn_nograd(x, fields_are_needed=False): + """A simple loss function""" + the_loss, field = simulate(x) + + if fields_are_needed: + return the_loss, field + return the_loss + + if wantgrad: + loss_value_npa, grad = autograd.value_and_grad(loss_fn)(design) # type: ignore + else: + loss_value = loss_fn_nograd(design) # type: ignore + # print("DIFF:", loss_value, loss_value_npa) + # loss_value, loss_grad = autograd.value_and_grad(loss_fn)(design) # type: ignore + first_time_ceviche = False + assert not (wantgrad and wantfields) + if wantgrad: + # if loss_value_npa< 0.0: + # print(x, "NP:", loss_fn_nograd(x), "NPA:", loss_fn(x)) + # assert False, f"superweird_{loss_fn_nograd(x)}_{loss_fn(x)}" + return loss_value_npa, grad + if wantfields: + loss_value, fields = loss_fn_nograd(design, fields_are_needed=True) + # if loss_value < 0.0: + # print(x, "NP:", loss_fn_nograd(x), "NPA:", loss_fn(x)) + # assert False, f"superweird_{loss_fn_nograd(x)}_{loss_fn(x)}" + return loss_value, fields + # if loss_value < 0.0: + # print(x, "NP:", loss_fn_nograd(x), "NPA:", loss_fn(x)) + # assert False, f"superweird_{loss_fn_nograd(x)}_{loss_fn(x)}" + return loss_value + +if __name__ == "__main__": + design_variable_shape = (345,345) + x = np.random.random(design_variable_shape) + loss = ceviche(x, benchmark_type=2, discretize=False) + print(loss) From 1d8add4d852fac978c245eef6dbff351b5f6c367 Mon Sep 17 00:00:00 2001 From: Denis Langevin <68906630+Milloupe@users.noreply.github.com> Date: Wed, 12 Feb 2025 15:57:51 +0100 Subject: [PATCH 2/5] Create mode_converter.py --- nevergrad/functions/mode_converter.py | 432 ++++++++++++++++++++++++++ 1 file changed, 432 insertions(+) create mode 100644 nevergrad/functions/mode_converter.py diff --git a/nevergrad/functions/mode_converter.py b/nevergrad/functions/mode_converter.py new file mode 100644 index 000000000..01fea33ce --- /dev/null +++ b/nevergrad/functions/mode_converter.py @@ -0,0 +1,432 @@ +import numpy as np +import scipy.sparse as sp +import scipy.sparse.linalg as lin +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +from functions import addupml2d, yeeder2d, block + +import time + +#%% CALCULATE S-PARAMETERS FOR WGs + + +c0 = 299792458 +epsi0 = 8.85418782e-12 +mu0 = 12.566370614e-7 +Z0 = np.sqrt(mu0/epsi0) + +# UNITS in micrometers +micrometers = 1 +nanometers = micrometers / 1000 + +######################################################### +## Using Ceviche mode converter dimentions +######################################################### + +def mode_converter(X, ev_out=-3.5**2/4, ev_in=-3.5**2, plot=False, show=False): + """ + Computes the conversion efficiency between the mode of index ev_in + in the input wg, into the mode of index ev_out in the output wg, + given a central block of dimensions: + centerWG_w = 2 * micrometers # Width + centerWG_L = 2 * micrometers # Length + described in X float array of size (centerWG_w x dx , centerWG_L x dy) + X is the permittivities + Default ev values compute fundamental mode to first order mode + """ + + # SOURCE + lam0 = 1.55 * micrometers + k0 = 2 * np.pi/lam0 + + MODE = 'E' # 'E' or 'H' + # TODO: change mode choice + # MODE_IN = 1 # MODE NUMBER INPUT PORT + # MODE_OUT = 1 # MODE NUMBER OUTPUT PORT + + # CENTER modifiable box PARAMETERS + centerWG_n = 3.5 # refractive index, if discretized + centerWG_w = 1.5 * micrometers # Width + centerWG_L = 1.5 * micrometers # Length + + # INPUT / OUTPUT WGs PARAMETERS + inputWG_n = 3.5 # refractive index + inputWG_L = 750 * nanometers # Width + inputWG_w = 200 * nanometers # Length + shift_in = 00 * nanometers # y-shift + + outputWG_n = 3.5 # refractive index + outputWG_L = 750 * nanometers # Width + outputWG_w = 400 * nanometers # Length + shift_out = -00 * nanometers # y-shift + + # FDFD PARAMETERS + NRES = 50 # GRID RESOLUTION + SPACER = lam0 * np.array([1, 1]) # Y SPACERS + NPML = [20, 20, 20, 20] # NB PML + nmax = max([inputWG_n, centerWG_n, outputWG_n]) + + + ### CALCULATE OPTIMIZED GRID + # GRID RESOLUTION + dx = lam0 / (nmax*NRES) + dy = lam0 / (nmax*NRES) + + # SNAP GRID TO CRITICAL DIMENSIONS + a = np.min([inputWG_w, centerWG_w, outputWG_w]) + nx = np.ceil(a/dx) + dx = a/nx + ny = np.ceil(a/dy) + dy = a/ny + + # GRID SIZE + grid_size_x = inputWG_L + centerWG_L + outputWG_L + Nx = int(NPML[0] + np.ceil(grid_size_x/dx) + NPML[1]) + grid_size_x = Nx * dx + + grid_size_y = SPACER[0] + centerWG_w + SPACER[1] + Ny = int(NPML[2] + np.ceil(grid_size_y/dy) + NPML[3]) + grid_size_y = Ny * dy + + # 2X GRID + Nx2 = 2 * Nx + dx2 = dx/2 + Ny2 = 2 * Ny + dy2 = dy/2 + + # CALCULATE AXIS VECTORS + xa = np.arange(1, Nx+1) * dx + ya = np.arange(1, Ny+1) * dy + xa2 = np.arange(1, Nx2+1) * dx2 + ya2 = np.arange(1, Ny2+1) * dy2 + # CENTER WINDOW (with (x,y) = (0,0) in the center) + xa = xa - np.mean(xa) + ya = ya - np.mean(ya) + xa2 = xa2 - np.mean(xa2) + ya2 = ya2 -np.mean(ya2) + + PML = [np.abs(xa2[2 * NPML[0] - 1]-xa2[0]), np.abs(xa2[Nx2-1]-xa2[Nx2-1-2 * NPML[1]])] + + + #%% BUILD OPTICAL INTEGRATED CIRCUIT + # INITIALIZE MATERIALS ARRAYS + ER2 = np.ones((Nx2,Ny2)) + UR2 = np.ones((Nx2,Ny2)) + + # INPUT WAVEGUIDE + pos_x = -grid_size_x/2 + pos_y = -inputWG_w/2 + Len_x = PML[0]+ inputWG_L + Len_y = inputWG_w + pos_y = pos_y + shift_in + ER2 = block(xa2, ya2, ER2, pos_x, pos_y, Len_x, Len_y, inputWG_n) + + + # OUTPUT WAVEGUIDE + pos_x = -grid_size_x/2 + PML[0] + inputWG_L+ centerWG_L + pos_y = - outputWG_w/2 + Len_x = outputWG_L + PML[1] + Len_y = outputWG_w + pos_y = pos_y + shift_out + ER2 = block(xa2, ya2, ER2, pos_x, pos_y, Len_x, Len_y, outputWG_n) + + + # ev_in = -inputWG_n**2 + # ev_out = - 2 * inputWG_n**2 # second mode ? + + nb_struct_x = int(centerWG_L / dx2) + nb_struct_y = int(centerWG_w / dy2) + assert (len(X) == nb_struct_x and len(X[0]) == nb_struct_y), f"nb_struct_x={nb_struct_x}, nb_struct_y={nb_struct_y}, np.shape(X)={np.shape(X)}" + + # CENTER WAVEGUIDE + X = 1 + X * (centerWG_n**2-1) # shift value range from 0-1 to 1-n + + center_x_beg = int(Nx - nb_struct_x/2) + center_y_beg = int(Ny - nb_struct_y/2) + + ER2[center_x_beg:center_x_beg+nb_struct_x, center_y_beg:center_y_beg+nb_struct_y] = X + + if plot: + plt.figure(0) + gs = gridspec.GridSpec(4,4) + plt.subplot(gs[2:, 1:3]) + plt.pcolor(xa2, ya2, np.sqrt(ER2.T)) + plt.colorbar() + plt.title('Refractive Index') + plt.xlabel("x ($\mu m$)") + plt.ylabel("y ($\mu m$)") + + # INCORPORATE PML + [inv_ERxx,inv_ERyy,ERzz,inv_URxx,inv_URyy,URzz] = addupml2d(ER2,UR2,NPML) + + ######################################################### + ## ANALYZE INPUT / OUTPUT WAVEGUIDES + ######################################################### + + # EXTRACT INPUT SLAB WAVEGUIDE FROM GRID + nx = 2 * NPML[0] + 1 + inv_erxx = sp.diags_array(1/ER2[nx, 0:Ny2:2], format="csc") + inv_eryy = sp.diags_array(1/ER2[nx, 1:Ny2:2], format="csc") + erzz = sp.diags_array(ER2[nx, 0:Ny2:2], format="csc") + inv_urxx = sp.diags_array(1/UR2[nx, 1:Ny2:2], format="csc") + inv_uryy = sp.diags_array(1/UR2[nx, 0:Ny2:2], format="csc") + urzz = sp.diags_array(UR2[nx, 1:Ny2:2], format="csc") + + # BUILD DERIVATIVE MATRICES + NS = [1, Ny] + RES = np.array([1, dy]) + BC = [0, 0] + [_,DEY,_,DHY] = yeeder2d(NS,k0 * RES,BC) + + # BUILD EIGEN-VALUE PROBLEM + if MODE == 'E': + A = -(DHY @ inv_urxx @ DEY + erzz) + B = inv_uryy + else: + A = -(DEY @ inv_erxx @ DHY + urzz) + B = inv_eryy + + + # SOLVE FULL EIGEN-VALUE PROBLEM + ev = ev_in + [D, input_mode] = lin.eigs(A, k=1, M=B, sigma=ev) # Find the eigen value closest to the -n**2 + D = np.sqrt(D) + NEFF = -1j * D + # GET SOURCE MODE + neff_in = NEFF[0] + normalization = np.sqrt(neff_in/(2*Z0) * (np.sum(np.abs(input_mode[:,0]) ** 2) * dy)) + input_mode = input_mode[:,0] / normalization # NORMALIZED MODE + + if plot: + plt.subplot(gs[:2, :2]) + plt.plot(np.real(input_mode), ya) + plt.xlabel('Real(U)'), + plt.ylabel('y ($\mu$m)') + plt.title(f"Input Mode, neff = {np.round(neff_in,3)}") + + # EXTRACT OUPUT SLAB WAVEGUIDE FROM GRID + nx = Nx2 - 2 * NPML[0] + 1 + erzz = sp.diags_array(ER2[nx, 0:Ny2:2], format="csc") + urxx = sp.diags_array(UR2[nx, 1:Ny2:2], format="csc") + inv_uryy = sp.diags_array(1/UR2[nx+1, 0:Ny2:2], format="csc") + + # BUILD DERIVATIVE MATRICES + NS = [1, Ny] + RES = np.array([1, dy]) + BC = [0, 0] + [_,DEY,_,DHY] = yeeder2d(NS,k0 * RES,BC) + + # BUILD EIGEN-VALUE PROBLEM + if MODE == 'E': + A = -(DHY @ inv_urxx @ DEY + erzz) + B = inv_uryy + else: + A = -(DEY @ inv_erxx @ DHY + urzz) + B = inv_eryy + + # np.savetxt("Debug", B.toarray(), fmt="%.3f") + + # SOLVE FULL EIGEN-VALUE PROBLEM + ev = ev_out + [D, output_mode] = lin.eigs(A, k=1, M=B, sigma=ev) # Find the eigen value closest to the -n**2 + D = np.sqrt(D) + NEFF = -1j * D + # GET SOURCE MODE + neff_out = NEFF[0] + normalization = np.sqrt(neff_out/(2*Z0) * (np.sum(np.abs(output_mode[:,0]) ** 2) * dy)) + output_mode = output_mode[:,0] /normalization + + if plot: + plt.subplot(gs[:2, 2:]) + plt.plot(np.real(output_mode), ya) + plt.xlabel('Real(U)') + plt.ylabel("y ($\mu$m)") + plt.title(f"Output Mode, neff = {np.round(neff_out,3)}") + plt.tight_layout() + plt.savefig("Testing_shapes.png") + + ######################################################### + ## PERFORM FDFD ANALYSIS + ######################################################### + + # BUILD DERIVATIVE MATRICES + NS = [Nx, Ny] + RES = np.array([dx, dy]) + BC = [0, 0] + [DEX,DEY,DHX,DHY] = yeeder2d(NS,k0 * RES,BC) + # print(DEX, DEY, DHX, DHY) + + + # BUILD WAVE MATRIX + if MODE == 'E': + A = DHX @ inv_URyy @ DEX + DHY @ inv_URxx @ DEY + ERzz + else: + A = DEX @ inv_ERyy @ DHX + DEY @ inv_ERxx @ DHY + URzz + + # np.savetxt("Debug", A.toarray()[:100,:100], fmt="%.3f") + # CALCULATE SOURCE FIELD + fsrc = np.zeros((Nx,Ny), dtype=complex) + for nx in range(Nx): + fsrc[nx,:] = input_mode * np.exp(-1j * k0 * neff_in * (nx+1) * dx) + + # CALCULATE SCATTERED-FIELD MASKING MATRIX + nx = NPML[0] + 2 + Q = np.zeros((Nx,Ny)) + Q[:nx,:] = 1 + Q = sp.diags_array(Q.flatten('F'), format="csc") + + + # CALCULATE SOURCE VECTOR + b = (Q @ A - A @ Q) @ fsrc.flatten('F') + # SOLVE FOR FIELD + f = lin.spsolve(A, b) + U = np.reshape(f,(Nx,Ny), order="F") + + ######################################################### + ## ANALYZE TRANSMITTED AND REFLECTED + ######################################################### + + # EXTRACT FIELDS + nx = NPML[0] + fref = U[nx,:] + nx = Nx - NPML[1]-2 + ftrn = U[nx,:] + + # CALCULATE MODE AMPLITUDES + #aref = input_mode\fref[:] + #atrn = output_mode\ftrn[:] + # CALCULATE REFLECTANCE AND TRANSMITTANCE + #R = abs(aref[0])^2 + #T = abs(atrn[0])^2 + + # CALCULATE S PARAMETERS + S11 = fref @ np.conj(input_mode) / (np.conj(input_mode).T @ input_mode) + S21 = ftrn @ np.conj(output_mode) / (np.conj(output_mode).T @ output_mode) + SR = abs(S11)**2 # REFLECTION + ST = abs(S21)**2 # TRANSMISSION + ######################################################### + ## DRAW SOLUTION + ######################################################### + if plot: + plt.figure(1) + plt.subplot(2, 1, 2) + plt.pcolormesh(xa, ya, np.real(U).T / np.max(np.abs(U))) + plt.contour(xa2, ya2, np.sqrt(ER2.T)) + plt.xlabel('x ($\\mu$m)') + plt.ylabel('y ($\\mu$m)') + plt.colorbar() + + plt.subplot(2,1,1) + plt.pie([SR, ST, 1-SR-ST]) + plt.legend(['S_R','S_T','Loss']) + plt.tight_layout() + plt.savefig("Testing_Ceviche_modeconverter.png") + if show: + plt.show() + return 1-ST, U + +#%% +if __name__ == "__main__": + + # SOURCE + lam0 = 1.55 * micrometers + k0 = 2 * np.pi/lam0 + + MODE = 'E' # 'E' or 'H' + # TODO: change mode choice + # MODE_IN = 1 # MODE NUMBER INPUT PORT + # MODE_OUT = 1 # MODE NUMBER OUTPUT PORT + + # CENTER modifiable box PARAMETERS + centerWG_n = 3.5 # refractive index, if discretized + centerWG_w = 1.5 * micrometers # Width + centerWG_L = 1.5 * micrometers # Length + + # INPUT / OUTPUT WGs PARAMETERS + inputWG_n = 3.5 # refractive index + inputWG_L = 750 * nanometers # Width + inputWG_w = 200 * nanometers # Length + shift_in = 00 * nanometers # y-shift + + outputWG_n = 3.5 # refractive index + outputWG_L = 750 * nanometers # Width + outputWG_w = 400 * nanometers # Length + shift_out = -00 * nanometers # y-shift + + # FDFD PARAMETERS + NRES = 50 # GRID RESOLUTION + SPACER = lam0 * np.array([1, 1]) # Y SPACERS + NPML = [20, 20, 20, 20] # NB PML + nmax = max([inputWG_n, centerWG_n, outputWG_n]) + + + ### CALCULATE OPTIMIZED GRID + # GRID RESOLUTION + dx = lam0 / (nmax*NRES) + dy = lam0 / (nmax*NRES) + + # SNAP GRID TO CRITICAL DIMENSIONS + a = np.min([inputWG_w, centerWG_w, outputWG_w]) + nx = np.ceil(a/dx) + dx = a/nx + ny = np.ceil(a/dy) + dy = a/ny + + # GRID SIZE + grid_size_x = inputWG_L + centerWG_L + outputWG_L + Nx = int(NPML[0] + np.ceil(grid_size_x/dx) + NPML[1]) + grid_size_x = Nx * dx + + grid_size_y = SPACER[0] + centerWG_w + SPACER[1] + Ny = int(NPML[2] + np.ceil(grid_size_y/dy) + NPML[3]) + grid_size_y = Ny * dy + + # 2X GRID + Nx2 = 2 * Nx + dx2 = dx/2 + Ny2 = 2 * Ny + dy2 = dy/2 + + # CALCULATE AXIS VECTORS + xa = np.arange(1, Nx+1) * dx + ya = np.arange(1, Ny+1) * dy + xa2 = np.arange(1, Nx2+1) * dx2 + ya2 = np.arange(1, Ny2+1) * dy2 + # CENTER WINDOW (with (x,y) = (0,0) in the center) + xa = xa - np.mean(xa) + ya = ya - np.mean(ya) + xa2 = xa2 - np.mean(xa2) + ya2 = ya2 -np.mean(ya2) + + PML = [np.abs(xa2[2 * NPML[0] - 1]-xa2[0]), np.abs(xa2[Nx2-1]-xa2[Nx2-1-2 * NPML[1]])] + + + #%% BUILD OPTICAL INTEGRATED CIRCUIT + # INITIALIZE MATERIALS ARRAYS + ER2 = np.ones((Nx2,Ny2)) + UR2 = np.ones((Nx2,Ny2)) + + # INPUT WAVEGUIDE + pos_x = -grid_size_x/2 + pos_y = -inputWG_w/2 + Len_x = PML[0]+ inputWG_L + Len_y = inputWG_w + pos_y = pos_y + shift_in + ER2 = block(xa2, ya2, ER2, pos_x, pos_y, Len_x, Len_y, inputWG_n) + + + # OUTPUT WAVEGUIDE + pos_x = -grid_size_x/2 + PML[0] + inputWG_L+ centerWG_L + pos_y = - outputWG_w/2 + Len_x = outputWG_L + PML[1] + Len_y = outputWG_w + pos_y = pos_y + shift_out + ER2 = block(xa2, ya2, ER2, pos_x, pos_y, Len_x, Len_y, outputWG_n) + + + nb_struct_x = int(centerWG_L / dx2) + nb_struct_y = int(centerWG_w / dy2) + X = X = np.ones((nb_struct_x, nb_struct_y))#(np.random.random((nb_struct_x, nb_struct_y))+1)**2 + print(nb_struct_x, nb_struct_y) + ST, field = mode_converter(X, ev_out=-inputWG_n**2/4, ev_in=-inputWG_n**2, plot=True, show=False) + print(ST) From 9d70ad5c8871f00fa0ce54f42d5e3f91e66d9575 Mon Sep 17 00:00:00 2001 From: Denis Langevin <68906630+Milloupe@users.noreply.github.com> Date: Wed, 12 Feb 2025 15:58:07 +0100 Subject: [PATCH 3/5] Create functions.py --- nevergrad/functions.py | 241 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 241 insertions(+) create mode 100644 nevergrad/functions.py diff --git a/nevergrad/functions.py b/nevergrad/functions.py new file mode 100644 index 000000000..9f03e2261 --- /dev/null +++ b/nevergrad/functions.py @@ -0,0 +1,241 @@ +import numpy as np +import scipy.sparse as sp + +def addupml2d(ER2,UR2,NPML): + """ + ADDUPML2D Add UPML to a 2D Yee Grid + + [ERxx,ERyy,ERzz,URxx,URyy,URzz] = addupml2d(ER2,UR2,NPML) + + INPUT ARGUMENTS + ================ + ER2 Relative Permittivity on 2x Grid + UR2 Relative Permeability on 2x Grid + NPML [NXLO NXHI NYLO NYHI] Size of UPML on 1x Grid + + OUTPUT ARGUMENTS + ================ + ERxx xx Tensor Element for Relative Permittivity + ERyy yy Tensor Element for Relative Permittivity + ERzz zz Tensor Element for Relative Permittivity + URxx xx Tensor Element for Relative Permeability + URyy yy Tensor Element for Relative Permeability + URzz zz Tensor Element for Relative Permeability + """ + ######################################################### + ## INITIALIZE FUNCTION + ######################################################### + + # DEFINE PML PARAMETERS + amax = 4 + cmax = 1 + p = 3 + + # EXTRACT GRID PARAMETERS + [Nx2,Ny2] = np.shape(ER2) + + # EXTRACT PML PARAMETERS + NXLO = 2 * NPML[0] + NXHI = 2 * NPML[1] + NYLO = 2 * NPML[2] + NYHI = 2 * NPML[3] + + ######################################################### + ## CALCULATE PML PARAMETERS + ######################################################### + + # INITIALIZE PML PARAMETERS TO PROBLEM SPACE + sx = np.ones((Nx2,Ny2), dtype=complex) + sy = np.ones((Nx2,Ny2), dtype=complex) + + # ADD XLO PML + for nx in range(1, NXLO+1): + ax = 1 + (amax - 1) * (nx/NXLO) ** p + cx = cmax * np.sin(0.5 * np.pi * nx/NXLO) ** 2 + sx[NXLO - nx,:] = ax * (1 - 1j * 60 * cx) + + + # ADD XHI PML + for nx in range(1, NXHI+1): + ax = 1 + (amax - 1) * (nx/NXHI) ** p + cx = cmax * np.sin(0.5 * np.pi * nx/NXHI) ** 2 + sx[Nx2 - NXHI + nx - 1,:] = ax * (1 - 1j * 60 * cx) + + + # ADD YLO PML + for ny in range(1, NYLO+1): + ay = 1 + (amax - 1) * (ny/NYLO) ** p + cy = cmax * np.sin(0.5 * np.pi * ny/NYLO) ** 2 + sy[:,NYLO - ny] = ay * (1 - 1j * 60 * cy) + + + # ADD YHI PML + for ny in range(1, NYHI+1): + ay = 1 + (amax - 1) * (ny/NYHI) ** p + cy = cmax * np.sin(0.5 * np.pi * ny/NYHI) ** 2 + sy[:,Ny2 - NYHI + ny - 1] = ay * (1 - 1j * 60 * cy) + + + ######################################################### + ## INCORPORATE PML + ######################################################### + + # CALCULATE TENSOR ELEMENTS WITH UPML + ERxx = ER2/sx * sy + ERyy = ER2 * sx/sy + ERzz = ER2 * sx * sy + + URxx = UR2/sx * sy + URyy = UR2 * sx/sy + URzz = UR2 * sx * sy + + # EXTRACT TENSOR ELEMENTS ON YEE GRID + ERxx = ERxx[1:Nx2:2, 0:Ny2:2] + ERyy = ERyy[0:Nx2:2, 1:Ny2:2] + ERzz = ERzz[0:Nx2:2, 0:Ny2:2] + + URxx = URxx[0:Nx2:2, 1:Ny2:2] + URyy = URyy[1:Nx2:2, 0:Ny2:2] + URzz = URzz[1:Nx2:2, 1:Ny2:2] + + # DIAGONALIZE MATERIAL TENSORS + inv_ERxx = sp.diags_array(1/ERxx.flatten('F'), format="csc") + inv_ERyy = sp.diags_array(1/ERyy.flatten('F'), format="csc") + ERzz = sp.diags_array(ERzz.flatten('F'), format="csc") + inv_URxx = sp.diags_array(1/URxx.flatten('F'), format="csc") + inv_URyy = sp.diags_array(1/URyy.flatten('F'), format="csc") + URzz = sp.diags_array(URzz.flatten('F'), format="csc") + + return [inv_ERxx,inv_ERyy,ERzz,inv_URxx,inv_URyy,URzz] + +def yeeder2d(NS,RES,BC,kinc=0): + """ + YEEDER2D Derivative Matrices on a 2D Yee Grid + + [DEX,DEY,DHX,DHY] = yeeder2d(NS,RES,BC,kinc) + + INPUT ARGUMENTS + ================= + NS [Nx Ny] Grid Size + RES [dx dy] Grid Resolution + BC [xbc ybc] Boundary Conditions + 0: Dirichlet boundary conditions + 1: Periodic boundary conditions + kinc [kx ky] Incident Wave Vector + This argument is only needed for PBCs. + + Note: For normalized grids use k0 * RES and kinc/k0 + + OUTPUT ARGUMENTS + ================= + DEX Derivative Matrix wrt x for Electric Fields + DEY Derivative Matrix wrt to y for Electric Fields + DHX Derivative Matrix wrt to x for Magnetic Fields + DHY Derivative Matrix wrt to y for Magnetic Fields + """ + + ######################################################### + ## HANDLE INPUT ARGUMENTS + ######################################################### + + # EXTRACT GRID PARAMETERS + Nx = NS[0] + dx = RES[0] + Ny = NS[1] + dy = RES[1] + + # DEFAULT KINC + if not(kinc): + kinc = [0, 0] + + + # DETERMINE MATRIX SIZE + M = Nx * Ny + + # ZERO MATRIX + Z = sp.csc_array((M,M)) + + ######################################################### + ## BUILD DEX + ######################################################### + + # HANDLE IF SIZE IS 1 CELL + if Nx==1: + DEX = -1j * kinc[0] * sp.eye_array(M) + + # HANDLE ALL OTHER CASES + else: + + # Center Diagonal + d0 = -np.ones(M) + + # Upper Diagonal + d1 = np.ones(M-1) + d1[Nx-1:M:Nx] = 0 + + # Build Derivative Matrix with Dirichlet BCs + DEX = Z + sp.diags_array(d0/dx) + DEX = DEX + sp.diags_array(d1/dx, offsets=1) + + # Incorporate Periodic Boundary Conditions + if BC[0]==1: + d1 = sp.csc_array(M) + d1[:M:Nx] = np.exp(-1j * kinc[0] * Nx * dx)/dx + DEX = DEX + sp.diags_array(d1, offsets=1-Nx) + + + + + ######################################################### + ## BUILD DEY + ######################################################### + + # HANDLE IF SIZE IS 1 CELL + if Ny==1: + DEY = -1j * kinc[1] * sp.eye_array(M) + + # HANDLE ALL OTHER CASES + else: + + # Center Diagonal + d0 = -np.ones((M)) + + # Upper Diagonal + d1 = np.ones((M-Nx)) + + # Build Derivative Matrix with Dirichlet BCs + DEY = Z + sp.diags_array(d0/dy) + DEY = DEY + sp.diags_array(d1/dy, offsets=Nx) + + # Incorporate Periodic Boundary Conditions + if BC[1]==1: + d1 = (np.exp(-1j * kinc[1] * Ny * dy)/dy) * np.ones((M,1)) + DEY = DEY + sp.diags_array(d1, offsets=Nx-M) + + + + + ######################################################### + ## BUILD DHX AND DHY + ######################################################### + + DHX = -np.conj(DEX.T) + DHY = -np.conj(DEY.T) + return [DEX,DEY,DHX,DHY] + + +def block(xa2,ya2,ER2,pos_x,pos_y,Len_x,Len_y,n): + dx2 = xa2[1]-xa2[0] + dy2 = ya2[1]-ya2[0] + # MAKE n-index block whose corner left-down is at pos_x and pos_y + # of size Len_x X Len_y + # refractive index n + + ind_X = np.argmin(np.abs(xa2-pos_x)) + ind_Y = np.argmin(np.abs(ya2-pos_y)) + + nx = int(ind_X + round(Len_x/dx2)) + ny = int(ind_Y + round(Len_y/dy2)) + + ER2 [ind_X:nx, ind_Y:ny] = n ** 2 + return ER2 From 86fd97634aa147757e7a2b3bf5f6f78af6612e7b Mon Sep 17 00:00:00 2001 From: Denis Langevin <68906630+Milloupe@users.noreply.github.com> Date: Wed, 12 Feb 2025 16:00:22 +0100 Subject: [PATCH 4/5] Rename nevergrad/functions.py to nevergrad/functions/photonics/functions.py --- nevergrad/{ => functions/photonics}/functions.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename nevergrad/{ => functions/photonics}/functions.py (100%) diff --git a/nevergrad/functions.py b/nevergrad/functions/photonics/functions.py similarity index 100% rename from nevergrad/functions.py rename to nevergrad/functions/photonics/functions.py From dc292177b2ef6d03a34accd67064d567888e212b Mon Sep 17 00:00:00 2001 From: Denis Langevin <68906630+Milloupe@users.noreply.github.com> Date: Wed, 12 Feb 2025 16:00:57 +0100 Subject: [PATCH 5/5] Rename nevergrad/functions/mode_converter.py to nevergrad/functions/photonics/mode_converter.py --- nevergrad/functions/{ => photonics}/mode_converter.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename nevergrad/functions/{ => photonics}/mode_converter.py (100%) diff --git a/nevergrad/functions/mode_converter.py b/nevergrad/functions/photonics/mode_converter.py similarity index 100% rename from nevergrad/functions/mode_converter.py rename to nevergrad/functions/photonics/mode_converter.py