Dear FenicsX Community,
I am currently working on a simulation of a Vertical-Cavity Surface-Emitting Laser (VCSEL) using FenicsX, and I have encountered an issue that I hope the community can help me with. I’ve attached my Python code below. The code is designed to generate a 3D mesh for the VCSEL structure, tag different regions within the mesh, and assign material properties to each tagged region.
from dolfinx import fem, io, mesh
import ufl
import numpy as np
from mpi4py import MPI
from dolfinx.mesh import locate_entities, meshtags, locate_entities_boundary
import dolfinx
# Define the geometric parameters of the VCSEL
d = 4e-6 # Aperture diameter, 4 micrometers
# Introduce a new parameter D representing the lateral dimension of the entire VCSEL structure, assuming d is 2/3 of D
D = d * 3 / 2
lambda_0 = 850e-9 # Central wavelength, 850 nm
bottom_dbr_pairs = 29.5 # Number of bottom DBR pairs
top_dbr_pairs = 24 # Number of top DBR pairs
qw_thickness = 5e-9 # Quantum well thickness
oxide_gaas_thickness = 69.49e-9 # Thickness of the GaAs layer in the oxide window
oxide_algaas_top_thickness = 63.71e-9 # Thickness of the upper AlGaAs layer in the oxide window
oxide_complex_thickness = 15.93e-9 # Thickness of the composite layer in the oxide window
active_gaas_thickness = 136.49e-9 # Thickness of the GaAs layer in the active region
ni = 0.01 # Assumed ni value, can be modified according to the actual situation
# Refractive indices of materials
n_GaAs = 3.53
n_AlGaAs = 3.08
n_AlAs = 2.95
n_AlOx = 1.60
# Given thickness of each DBR layer in the original text
dbr_GaAs_thickness = 69.49e-9
dbr_AlGaAs_thickness = 79.63e-9
# Assume x is the 5 varying positions of the oxide aperture
x_positions = np.linspace(0, 63.71e-9, 5)
# Process each position in a loop
for i, x in enumerate(x_positions):
# Calculate the lengths of each part
top_dbr_length = top_dbr_pairs * (dbr_GaAs_thickness + dbr_AlGaAs_thickness)
oxide_window_length = oxide_gaas_thickness + (oxide_algaas_top_thickness - x) + oxide_complex_thickness + x
active_layer_length = 2 * active_gaas_thickness + qw_thickness
bottom_dbr_length = bottom_dbr_pairs * (dbr_GaAs_thickness + dbr_AlGaAs_thickness)
total_length = top_dbr_length + oxide_window_length + active_layer_length + bottom_dbr_length
# Generate a 3D mesh, increase the number of meshes in the longitudinal direction for better resolution of each layer
nx = 20
ny = 20
nz = 1000
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([-D / 2, -D / 2, 0]), np.array([D / 2, D / 2, total_length])],
(nx, ny, nz), cell_type=mesh.CellType.hexahedron)
# Define the tagging functions for different regions
def region_func(z_start, z_end, r_cond=lambda r: True):
return lambda x_coords: np.logical_and(
np.logical_and(x_coords[2] >= z_start, x_coords[2] < z_end),
r_cond(np.sqrt(x_coords[0] ** 2 + x_coords[1] ** 2))
)
regions = {
'top_dbr': region_func(0, top_dbr_length),
'oxide_gaas': region_func(top_dbr_length, top_dbr_length + oxide_gaas_thickness),
'oxide_algaas_top': region_func(
top_dbr_length + oxide_gaas_thickness,
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x)
),
'oxide_complex_center': region_func(
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x),
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x) + oxide_complex_thickness,
lambda r: r < d / 2
),
'oxide_complex_outer': region_func(
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x),
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x) + oxide_complex_thickness,
lambda r: r >= d / 2
),
'oxide_algaas_bottom': region_func(
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x) + oxide_complex_thickness,
top_dbr_length + oxide_gaas_thickness + (oxide_algaas_top_thickness - x) + oxide_complex_thickness + x
),
'active_gaas_top': region_func(
top_dbr_length + oxide_window_length,
top_dbr_length + oxide_window_length + active_gaas_thickness
),
'quantum_well_center': region_func(
top_dbr_length + oxide_window_length + active_gaas_thickness,
top_dbr_length + oxide_window_length + active_gaas_thickness + qw_thickness,
lambda r: r < d / 2
),
'quantum_well_outer': region_func(
top_dbr_length + oxide_window_length + active_gaas_thickness,
top_dbr_length + oxide_window_length + active_gaas_thickness + qw_thickness,
lambda r: r >= d / 2
),
'active_gaas_bottom': region_func(
top_dbr_length + oxide_window_length + active_gaas_thickness + qw_thickness,
top_dbr_length + oxide_window_length + 2 * active_gaas_thickness + qw_thickness
),
'bottom_dbr': region_func(
top_dbr_length + oxide_window_length + active_layer_length,
top_dbr_length + oxide_window_length + active_layer_length + bottom_dbr_length
),
'substrate': region_func(
top_dbr_length + oxide_window_length + active_layer_length + bottom_dbr_length,
total_length
)
}
# Initialize the tagging array (the same size as the number of mesh cells)
tdim = domain.topology.dim
index_map = domain.topology.index_map(tdim)
num_local_cells = index_map.size_local
cell_indices = np.arange(num_local_cells, dtype=np.int32)
global_cell_indices = index_map.local_to_global(cell_indices)
values = np.zeros(num_local_cells, dtype=np.int32)
# Define the tagging values
TAGS = {
'top_dbr_GaAs': 1, 'top_dbr_AlGaAs': 2,
'oxide_gaas': 3, 'oxide_algaas_top': 4, 'oxide_complex_center': 5, 'oxide_complex_outer': 6, 'oxide_algaas_bottom': 7,
'active_gaas_top': 8, 'quantum_well_center': 9, 'quantum_well_outer': 10, 'active_gaas_bottom': 11,
'bottom_dbr_GaAs': 12, 'bottom_dbr_AlGaAs': 13,
'substrate': 14
}
# Tag the top DBR layer
z = top_dbr_length
for j in range(int(top_dbr_pairs)):
z -= dbr_AlGaAs_thickness + dbr_GaAs_thickness
gaas_region = region_func(z, z + dbr_GaAs_thickness)
gaas_cells = locate_entities(domain, tdim, gaas_region)
local_gaas_cells = np.where(np.isin(global_cell_indices, gaas_cells))[0]
values[local_gaas_cells] = TAGS['top_dbr_GaAs']
z += dbr_GaAs_thickness
algaas_region = region_func(z, z + dbr_AlGaAs_thickness)
algaas_cells = locate_entities(domain, tdim, algaas_region)
local_algaas_cells = np.where(np.isin(global_cell_indices, algaas_cells))[0]
values[local_algaas_cells] = TAGS['top_dbr_AlGaAs']
z += dbr_AlGaAs_thickness
# Tag the oxide window layer
for region_name, func in [
('oxide_gaas', regions['oxide_gaas']),
('oxide_algaas_top', regions['oxide_algaas_top']),
('oxide_algaas_bottom', regions['oxide_algaas_bottom'])
]:
cells = locate_entities(domain, tdim, func)
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS[region_name]
# Tag the oxide composite layer, distinguish between inner and outer parts
cells = locate_entities(domain, tdim, regions['oxide_complex_center'])
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS['oxide_complex_center']
cells = locate_entities(domain, tdim, regions['oxide_complex_outer'])
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS['oxide_complex_outer']
# Tag the active region
for region_name, func in [
('active_gaas_top', regions['active_gaas_top']),
('active_gaas_bottom', regions['active_gaas_bottom'])
]:
cells = locate_entities(domain, tdim, func)
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS[region_name]
# Tag the quantum well layer, distinguish between inner and outer parts
cells = locate_entities(domain, tdim, regions['quantum_well_center'])
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS['quantum_well_center']
cells = locate_entities(domain, tdim, regions['quantum_well_outer'])
local_cells = np.where(np.isin(global_cell_indices, cells))[0]
values[local_cells] = TAGS['quantum_well_outer']
# Tag the bottom DBR layer
z = top_dbr_length + oxide_window_length + active_layer_length
for j in range(int(bottom_dbr_pairs)):
z += dbr_AlGaAs_thickness + dbr_GaAs_thickness
algaas_region = region_func(z - dbr_AlGaAs_thickness, z)
algaas_cells = locate_entities(domain, tdim, algaas_region)
local_algaas_cells = np.where(np.isin(global_cell_indices, algaas_cells))[0]
values[local_algaas_cells] = TAGS['bottom_dbr_AlGaAs']
z -= dbr_AlGaAs_thickness
gaas_region = region_func(z - dbr_GaAs_thickness, z)
gaas_cells = locate_entities(domain, tdim, gaas_region)
local_gaas_cells = np.where(np.isin(global_cell_indices, gaas_cells))[0]
values[local_gaas_cells] = TAGS['bottom_dbr_GaAs']
z -= dbr_GaAs_thickness
# Tag the substrate
substrate_cells = locate_entities(domain, tdim, regions['substrate'])
local_substrate_cells = np.where(np.isin(global_cell_indices, substrate_cells))[0]
values[local_substrate_cells] = TAGS['substrate']
# Correctly initialize MeshTags
marker = meshtags(domain, tdim, global_cell_indices, values)
# # Save the tagging information (optional)
# with open(f'vcsel_mesh_tags_x_{i}.txt', 'w') as f:
# f.write("Cell index, Tag\n")
# for cell_idx, tag in zip(global_cell_indices, values):
# f.write(f"{cell_idx}, {tag}\n")
# Define the material properties
materials = {
1: n_GaAs,
2: n_AlGaAs,
3: n_GaAs,
4: n_AlGaAs,
5: n_AlOx,
6: n_AlAs,
7: n_AlGaAs,
8: n_GaAs,
9: n_GaAs + 1j * ni,
10: n_GaAs,
11: n_GaAs,
12: n_GaAs,
13: n_AlGaAs,
14: n_GaAs
}
material_values = np.array([materials[tag] for tag in values], dtype=np.complex128)
material = fem.Function(domain, ufl.FunctionSpace(domain, ("DG", 0)))
material.x.array[:] = material_values
When I run this code, I receive the following error:
Traceback (most recent call last):
File "/home/xxx_77/program/xxx.py", line 223, in <module>
material_values = np.array([materials[tag] for tag in values], dtype=np.complex128)
~~~~~~~~~^^^^^
KeyError: 0
I’m noob at python，please help me！ I would greatly appreciate any help, advice, or insights from the community.
Thank you in advance!