Problems with point source `python: step_db.cpp:40: void meep::fields::step_db(meep::field_type): Assertion 'changed_materials' failed.` #2787

niccolot · 2024-02-20T11:28:43Z

niccolot
Feb 20, 2024

Hello, i'm trying to modify the tutorials on adjoint optimization, in particular the 3rd one (https://nbviewer.org/github/NanoComp/meep/blob/master/python/examples/adjoint_optimization/03-Filtered_Waveguide_Bend.ipynb)

I've set the waveguide to be straight and i'm trying to maximize the coupling of the waveguide with a source placed at the center. The optimization variables will be constrained to be fixed in a small radius around the source.

The problem is that the code works with the source provided in the tutorial, i.e.

source_size = mp.Vector3(0, Sy, 0)
kpoint = mp.Vector3(1, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src,
        eig_band=1,
        direction=mp.NO_DIRECTION,
        eig_kpoint=kpoint,
        size=source_size,
        center=design_region.center,
    )
]

but when using a point dipole source

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [mp.Source(src, component=mp.Ez, center=design_region.center)]

or even with a small size like size=mp.Vector3(0.1,0.1,0.1)

the code throws the error

python: step_db.cpp:40: void meep::fields::step_db(meep::field_type): Assertion 'changed_materials' failed.

The design plot is

while the code is

import meep as mp
import meep.adjoint as mpa
from meep import Animate2D
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product
import nlopt
from matplotlib import pyplot as plt

mp.verbosity(1)

Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)

waveguide_width = 0.5
design_region_width = 2.5
design_region_height = 2.5
protected_radius = 0.2

waveguide_length = 0.5

pml_size = 1.0

resolution = 20

minimum_length = 0.09  # minimum length scale (microns)
eta_i = (
    0.5  # blueprint (or intermediate) design field thresholding point (between 0 and 1)
)
eta_e = 0.55  # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e  # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
design_region_resolution = int(1 * resolution)

Sx = 2 * pml_size + 2 * waveguide_length + design_region_width
Sy = 2 * pml_size + design_region_height + 0.5
cell_size = mp.Vector3(Sx, Sy)

pml_layers = [mp.PML(pml_size)]

fcen = 1 / 1.55
width = 0.1
fwidth = width * fcen
freqs = 1 / np.linspace(1.5, 1.6, 10)

Nx = int(design_region_resolution * design_region_width) + 1
Ny = int(design_region_resolution * design_region_height) + 1

design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), SiO2, Si, grid_type="U_MEAN", do_averaging=False)
design_region = mpa.DesignRegion(
    design_variables,
    volume=mp.Volume(
        center=mp.Vector3(),
        size=mp.Vector3(design_region_width, design_region_height, 0),
    ),
)

x_g = np.linspace(-design_region_width / 2, design_region_width / 2, Nx)
y_g = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g, Y_g = np.meshgrid(x_g, y_g, sparse=True, indexing="ij")

left_wg_mask = (X_g == -design_region_width / 2) & (np.abs(Y_g) <= waveguide_width / 2)
right_wg_mask = (X_g == design_region_width / 2) & (np.abs(Y_g) <= waveguide_width / 2)
radius_mask = (X_g**2 + Y_g**2 < protected_radius**2) # force to have material around the dipole
Si_mask = left_wg_mask | right_wg_mask | radius_mask

border_mask = (
    (X_g == -design_region_width / 2)
    | (X_g == design_region_width / 2)
    | (Y_g == -design_region_height / 2)
    | (Y_g == design_region_height / 2)
)

SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False

def mapping(x, eta, beta):
    x = npa.where(Si_mask.flatten(), 1, npa.where(SiO2_mask.flatten(), 0, x))

    # filter
    filtered_field = mpa.conic_filter(
        x,
        filter_radius,
        design_region_width,
        design_region_height,
        design_region_resolution,
    )

    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)
    
    projected_field = (npa.fliplr(projected_field) + projected_field) / 2  # up-down symmetry
    projected_field = (npa.flipud(projected_field) + projected_field) / 2  # left-right symmetry

    # interpolate to actual materials
    return projected_field.flatten()


geometry = [
    mp.Block(
        center=mp.Vector3(x=-Sx / 4), material=Si, size=mp.Vector3(Sx / 2, 0.5, 0)
    ),  # left waveguide
    mp.Block(
        center=mp.Vector3(x=Sx / 4), material=Si, size=mp.Vector3(Sx / 2, 0.5, 0)
    ),  # right waveguide
    mp.Block(
        center=design_region.center, size=design_region.size, material=design_variables
    ),  # design region
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [mp.Source(src, component=mp.Ez, center=design_region.center)]

""" # commenting the previous source and uncommentig this one is ok
source_size = mp.Vector3(0, Sy, 0)
kpoint = mp.Vector3(1, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src,
        eig_band=1,
        direction=mp.NO_DIRECTION,
        eig_kpoint=kpoint,
        size=source_size,
        center=design_region.center,
    )
]
"""


sim = mp.Simulation(
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    default_material=SiO2,
    resolution=resolution,
)

mode = 1

TE_left = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
    ),
    mode,
)

TE_right = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
    ),
    mode,
)

ob_list = [TE_left, TE_right]

def J(left, right):
    return npa.mean(npa.abs(left) ** 2 + npa.abs(right) ** 2)

opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=J,
    objective_arguments=ob_list,
    design_regions=[design_region],
    frequencies=freqs
)

evaluation_history = []
cur_iter = [0]

def f(v, gradient, cur_beta):
    print("Current iteration: {}".format(cur_iter[0] + 1))

    f0, dJ_du = opt([mapping(v, eta_i, cur_beta)])  # compute objective and gradient

    if gradient.size > 0:
        gradient[:] = tensor_jacobian_product(mapping, 0)(
            v, eta_i, cur_beta, np.sum(dJ_du, axis=1)
        )  # backprop

    evaluation_history.append(np.real(f0))

    cur_iter[0] = cur_iter[0] + 1

    return np.real(f0)

algorithm = nlopt.LD_MMA

# Initial guess
x = np.ones((Nx,Ny)) * 0.5

x[Si_mask] = 1  # set the edges of waveguides to silicon
x[SiO2_mask] = 0  # set the other edges to SiO2

x = x.flatten()

# lower and upper bounds
lb = np.zeros((Nx*Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Ny,))
ub[SiO2_mask.flatten()] = 0



cur_beta = 4
beta_scale = 2
num_betas = 6
update_factor = 12
for iters in range(num_betas):
    solver = nlopt.opt(algorithm, Nx*Ny)
    solver.set_lower_bounds(lb)
    solver.set_upper_bounds(ub)
    solver.set_max_objective(lambda a, g: f(a, g, cur_beta))
    solver.set_maxeval(update_factor)
    x[:] = solver.optimize(x)
    cur_beta = cur_beta * beta_scale


plt.figure()
plt.plot(10 * np.log10(evaluation_history), "o-")
plt.grid(True)
plt.xlabel("Iteration")
plt.ylabel("Average Insertion Loss (dB)")
plt.savefig("straightWG_evaluation_history_dipole.png")


opt.update_design([mapping(x, eta_i, cur_beta)])

f0, dJ_du = opt([mapping(x, eta_i, cur_beta // 2)], need_gradient=False)
frequencies = opt.frequencies
right_coef = opt.get_objective_arguments()

right_profile = np.abs(right_coef) ** 2

plt.figure()
plt.plot(1 / frequencies, right_profile * 100, "-o")
plt.grid(True)
plt.xlabel("Wavelength (microns)")
plt.ylabel("Transmission (%)")
plt.savefig("straightWG_transmission_dipole.png")


plt.figure()
plt.plot(1 / frequencies, 10 * np.log10(right_profile), "-o")
plt.grid(True)
plt.xlabel("Wavelength (microns)")
plt.ylabel("Insertion Loss (dB)")
plt.savefig("straightWG_insertion_loss_dipole.png")

Answered by smartalecH

Feb 29, 2024

This looks to be the same issue encountered in #2309. There's a simple hack, as described in my comment: set force_all_components = True in your Simulation constructor.

Hopefully we can fix the root someday...

View full answer

oskooi · 2024-02-21T06:20:24Z

oskooi
Feb 21, 2024
Collaborator

Placing a source within the design region requires some care. This is something that @mochen4 faced when setting up the cavity-design problem which involved an objective function based on the LDOS at the center of the cavity. The details are described in Section E of JOSA B, Vol. 41, pp. A161-A176 (2024).

Perhaps @mochen4 can provide some guidance for how to set this up.

1 reply

niccolot Feb 21, 2024
Author

The weird thing is that placing and EigenMode source with directionality works even when placed at the center, it seems the problem is with the normal mp.Source object, regardless of the point-size or finite size. Some guidance would be kindly accepted

mochen4 · 2024-02-21T15:10:11Z

mochen4
Feb 21, 2024
Collaborator

Actually, I don't think we need any extra treatment for sources inside the design region. I used a point dipole source and it worked fine. I can check your code and see what is the issue. In the meantime, can you try putting a small circle of air around the source on top of the design region? Can you also try running the simulation itself, with some random configuration of design variables, without the EigenMode monitors?

3 replies

niccolot Feb 21, 2024
Author

I tryed to put a small circle (mp.Cylinder with heigth=0 and radius=0.2 and tried with a small square also) of both Si and air on top of the design region (appending this at the end of the geometry list) but the same problem persist in both cases. Running the simulation (sim.run(until=10000) right after the sim object and commenting everything below) with some random design variables (design_region.update_parameters(np.random.rand(Nx,Ny)) right after defining the design_region object) is fine, it ends the simulation without a problem.

Thank you a lot, every help is higly appreciated

mochen4 Feb 21, 2024
Collaborator

In that case, I don't think the issue is about putting point source in the design region.
I assume the error occurs during the forward run not the adjoint run, right?
If you don't have the materialgrid/design region at all (equivalently, configure the design variables to be 0 everywhere), and just add the eigenmode monitor (sim.add_mode_monitor) yourself, do you still see the error?

(btw, sim.run(until=100) should be sufficient for debugging purposes.)

niccolot Feb 21, 2024
Author

Yes, the errors appear just after running the code (mpirun -np 4 python codename.py). I added the two monitors (left and right) at the same positions than before with

mon_lef = sim.add_mode_monitor(fcen, fwidth, len(freqs), mp.ModeRegion(center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),))

mon_right = sim.add_mode_monitor(fcen, fwidth, len(freqs), mp.ModeRegion(center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),))

just after the simulation obj definition and before the sim.run. The code runs fine with both mp.FLuxRegion and mp.ModeRegion

mochen4 · 2024-02-22T03:39:20Z

mochen4
Feb 22, 2024
Collaborator

Actually, I ran the code and found the error occurs during the adjoint run, so the issue is related to the mpa.EigenmodeCoefficient. Any ideas @smartalecH ?

2 replies

niccolot Feb 23, 2024
Author

I tried to put some mpa.Near2farFields in the left and right parts of the waveguide and the optimization gives sensible results with the point source inside the design region also, so the problem seems indeed to be related with the mpa.EigenmodeCoefficient monitors

farPoints = [mp.Vector3(x=-2),
             mp.Vector3(x=2)]

NearRegions = [
    mp.Near2FarRegion(
        center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
        weight=-1,
    ),
    mp.Near2FarRegion(
        center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
        weight=+1,
    ),
]

FarFields = mpa.Near2FarFields(sim, NearRegions, farPoints)
ob_list = [FarFields]

def J(FF):
    return npa.mean(npa.abs(FF[0, :, 2]) ** 2)


opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=[J],
    objective_arguments=ob_list,
    design_regions=[design_region],
    frequencies=freqs
)

niccolot Feb 23, 2024
Author

Also by using 2 separate design region with a waveguide in the middle i get the same error (with point source and mpa.EigenmodeCoefficient monitors)

Code:

import meep as mp
import meep.adjoint as mpa
from meep import Animate2D
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product
import nlopt
from matplotlib import pyplot as plt

mp.verbosity(1)

Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)

waveguide_width = 0.5
design_region_width = 2.5
design_region_height = 2.5
protected_radius = 0.2

waveguide_length = 0.25

pml_size = 1.0

resolution = 20

minimum_length = 0.09  # minimum length scale (microns)
eta_i = (
    0.5  # blueprint (or intermediate) design field thresholding point (between 0 and 1)
)
eta_e = 0.55  # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e  # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
design_region_resolution = int(1 * resolution)

Sx = 2 * pml_size + 3 * waveguide_length + 2 * design_region_width
Sy = 2 * pml_size + design_region_height + 0.5
cell_size = mp.Vector3(Sx, Sy)

pml_layers = [mp.PML(pml_size)]

fcen = 1 / 1.55
width = 0.1
fwidth = width * fcen
freqs = 1 / np.linspace(1.5, 1.6, 10)

Nx = int(design_region_resolution * design_region_width) + 1
Ny = int(design_region_resolution * design_region_height) + 1

design_variables1 = mp.MaterialGrid([Nx,Ny,0], SiO2, Si, grid_type="U_MEAN", do_averaging=False)
design_variables2 = mp.MaterialGrid([Nx,Ny,0], SiO2, Si, grid_type="U_MEAN", do_averaging=False)

design_region1 = mpa.DesignRegion(
    design_variables1,
    volume=mp.Volume(
        center=mp.Vector3(x=(waveguide_length/2) + (design_region_width/2)),
        size=mp.Vector3(design_region_width, design_region_height, 0),
    ),
)

design_region2 = mpa.DesignRegion(
    design_variables2,
    volume=mp.Volume(
        center=mp.Vector3(x=-(waveguide_length/2) - (design_region_width/2)),
        size=mp.Vector3(design_region_width, design_region_height, 0),
    ),
)

design1_left = - waveguide_length/2 - design_region_width
design1_right = -waveguide_length/2
design2_left  = waveguide_length/2
design2_right = waveguide_length/2 + design_region_width

leftWGcenter = -waveguide_length/2 - design_region_width - waveguide_length / 2
rightWGcenter = +waveguide_length/2 + design_region_width + waveguide_length / 2

x_g1 = np.linspace(design1_left, design1_right, Nx)
y_g1 = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g1, Y_g1 = np.meshgrid(x_g1, y_g1, sparse=True, indexing="ij")

x_g2 = np.linspace(design2_left, design2_right, Nx)
y_g2 = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g2, Y_g2 = np.meshgrid(x_g2, y_g2, sparse=True, indexing="ij")

left_wg_mask1 = (X_g1 == design1_left) & (np.abs(Y_g1) <= waveguide_width / 2)
right_wg_mask1 = (X_g1 == design1_right) & (np.abs(Y_g1) <= waveguide_width / 2)

left_wg_mask2 = (X_g2 == design2_left) & (np.abs(Y_g2) <= waveguide_width / 2)
right_wg_mask2 = (X_g2 == design2_right) & (np.abs(Y_g2) <= waveguide_width / 2)

Si_mask1 = left_wg_mask1 | right_wg_mask1
Si_mask2 = left_wg_mask2 | right_wg_mask2 

Si_mask1 = np.array(Si_mask1)
Si_mask2 = np.array(Si_mask2)

Si_mask = np.stack((Si_mask1, Si_mask2))

border_mask1 = (
    (X_g1 == design1_left)
    | (X_g1 == design1_right)
    | (Y_g1 == -design_region_height / 2)
    | (Y_g1 == design_region_height / 2)
)

border_mask2 = (
    (X_g2 == design2_left)
    | (X_g2 == design2_right)
    | (Y_g2 == -design_region_height / 2)
    | (Y_g2 == design_region_height / 2)
)

border_mask1 = np.array(border_mask1)
border_mask2 = np.array(border_mask2)

border_mask = np.stack((border_mask1, border_mask2))

SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False

def mapping1(x, eta, beta):
    x = x.reshape((Nx,Ny))
    x = npa.where(Si_mask[0], 1, npa.where(SiO2_mask[0], 0, x))

    # filter
    filtered_field = mpa.conic_filter(
        x,
        filter_radius,
        design_region_width,
        design_region_height,
        design_region_resolution,
    )

    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)
    
    projected_field = (npa.fliplr(projected_field) + projected_field) / 2  # up-down symmetry
    projected_field = (npa.flipud(projected_field) + projected_field) / 2  # left-right symmetry

    # interpolate to actual materials
    return projected_field.flatten()

def mapping2(x, eta, beta):
    x = x.reshape((Nx,Ny))
    x = npa.where(Si_mask[1], 1, npa.where(SiO2_mask[1], 0, x))

    # filter
    filtered_field = mpa.conic_filter(
        x,
        filter_radius,
        design_region_width,
        design_region_height,
        design_region_resolution,
    )

    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)
    
    projected_field = (npa.fliplr(projected_field) + projected_field) / 2  # up-down symmetry
    projected_field = (npa.flipud(projected_field) + projected_field) / 2  # left-right symmetry

    # interpolate to actual materials
    return projected_field.flatten()


geometry = [
    mp.Block(
        center=mp.Vector3(x=leftWGcenter), material=Si, size=mp.Vector3(waveguide_length, 0.5, 0)
    ),  # left waveguide
    mp.Block(
        center=mp.Vector3(), material=Si, size=mp.Vector3(waveguide_length, 0.5, 0)
    ),  # center waveguide
    mp.Block(
        center=mp.Vector3(x=rightWGcenter), material=Si, size=mp.Vector3(waveguide_length, 0.5, 0)
    ),  # right waveguide
    mp.Block(center=design_region1.center, size=design_region1.size, material=design_variables1), # left design region
    mp.Block(center=design_region2.center, size=design_region2.size, material=design_variables2), # right design region
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [mp.Source(src, component=mp.Ez, center=mp.Vector3())]

sim = mp.Simulation(
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    default_material=SiO2,
    resolution=resolution,
)

mode = 1

TE_left = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
    ),
    mode,
)

TE_right = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
    ),
    mode,
)

ob_list = [TE_left, TE_right]

def J(left, right):
    return npa.mean(npa.abs(left) ** 2 + npa.abs(right) ** 2)

opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=J,
    objective_arguments=ob_list,
    design_regions=[design_region1, design_region2],
    frequencies=freqs
)

evaluation_history = []
cur_iter = [0]

def f(v, gradient, cur_beta):
    print("Current iteration: {}".format(cur_iter[0] + 1))
    v = v.reshape((2,Nx,Ny))
    f0, dJ_du = opt([mapping1(v[0], eta_i, cur_beta), mapping2(v[1], eta_i, cur_beta)])  # compute objective and gradient

    if gradient.size > 0:
        gradient[:] = tensor_jacobian_product(mapping1, mapping2, 0)(
            v, eta_i, cur_beta, np.sum(dJ_du, axis=1)
        )  # backprop

    evaluation_history.append(np.real(f0))

    cur_iter[0] = cur_iter[0] + 1

    return np.real(f0)

algorithm = nlopt.LD_MMA

# Initial guess
x1 = np.ones((Nx,Ny)) * 0.5
x2 = np.ones((Nx,Ny)) * 0.5

x1[Si_mask[0]] = 1  # set the edges of waveguides to silicon
x1[SiO2_mask[0]] = 0  # set the other edges to SiO2
x2[Si_mask[1]] = 1  # set the edges of waveguides to silicon
x2[SiO2_mask[1]] = 0  # set the other edges to SiO2
x = np.stack((x1,x2)).flatten()


# lower and upper bounds
lb = np.zeros((2*Nx*Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((2*Nx*Ny,))
ub[SiO2_mask.flatten()] = 0


plt.figure()
rho_vector1 = np.random.rand(Nx,Ny)
rho_vector2 = np.random.rand(Nx,Ny)
opt.update_design([mapping1(rho_vector1, eta_i, 1e3), mapping2(rho_vector2, eta_i, 1e3)])
opt.plot2D(True, output_plane=mp.Volume(center=(0, 0, 0), size=(Sx, Sy, 0)))
plt.savefig("doubleRegion_designRegion.png")

cur_beta = 4
beta_scale = 2
num_betas = 6
update_factor = 12

for iters in range(num_betas):
    solver = nlopt.opt(algorithm, 2*Nx*Ny)
    solver.set_lower_bounds(lb)
    solver.set_upper_bounds(ub)
    solver.set_max_objective(lambda a, g: f(a, g, cur_beta))
    solver.set_maxeval(update_factor)
    x[:] = solver.optimize(x)
    cur_beta = cur_beta * beta_scale

    plt.figure()
plt.plot(10 * np.log10(evaluation_history), "o-")
plt.grid(True)
plt.xlabel("Iteration")
plt.ylabel("Average Insertion Loss (dB)")
plt.savefig("doubleDesign_evaluation_history_dipole.png")


opt.update_design([mapping1(x1, eta_i, cur_beta), mapping2(x2, eta_i, cur_beta)])

f0, dJ_du = opt([mapping1(x1, eta_i, cur_beta // 2), mapping2(x2, eta_i, cur_beta // 2)], need_gradient=False)
frequencies = opt.frequencies
right_coef = opt.get_objective_arguments()

right_profile = np.abs(right_coef) ** 2

plt.figure()
plt.plot(1 / frequencies, right_profile * 100, "-o")
plt.grid(True)
plt.xlabel("Wavelength (microns)")
plt.ylabel("Transmission (%)")
plt.savefig("doubleDesign_transmission_dipole.png")


plt.figure()
plt.plot(1 / frequencies, 10 * np.log10(right_profile), "-o")
plt.grid(True)
plt.xlabel("Wavelength (microns)")
plt.ylabel("Insertion Loss (dB)")
plt.savefig("doubleDesign_insertion_loss_dipole.png")

Design setup:

smartalecH · 2024-02-29T20:22:11Z

smartalecH
Feb 29, 2024
Collaborator

This looks to be the same issue encountered in #2309. There's a simple hack, as described in my comment: set force_all_components = True in your Simulation constructor.

Hopefully we can fix the root someday...

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Problems with point source `python: step_db.cpp:40: void meep::fields::step_db(meep::field_type): Assertion 'changed_materials' failed.` #2787

{{title}}

Replies: 4 comments 6 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Problems with point source python: step_db.cpp:40: void meep::fields::step_db(meep::field_type): Assertion 'changed_materials' failed. #2787

niccolot Feb 20, 2024

Replies: 4 comments · 6 replies

oskooi Feb 21, 2024 Collaborator

niccolot Feb 21, 2024 Author

mochen4 Feb 21, 2024 Collaborator

niccolot Feb 21, 2024 Author

mochen4 Feb 21, 2024 Collaborator

niccolot Feb 21, 2024 Author

mochen4 Feb 22, 2024 Collaborator

niccolot Feb 23, 2024 Author

niccolot Feb 23, 2024 Author

smartalecH Feb 29, 2024 Collaborator

Problems with point source `python: step_db.cpp:40: void meep::fields::step_db(meep::field_type): Assertion 'changed_materials' failed.` #2787

niccolot
Feb 20, 2024

Replies: 4 comments 6 replies

oskooi
Feb 21, 2024
Collaborator

niccolot Feb 21, 2024
Author

mochen4
Feb 21, 2024
Collaborator

niccolot Feb 21, 2024
Author

mochen4 Feb 21, 2024
Collaborator

niccolot Feb 21, 2024
Author

mochen4
Feb 22, 2024
Collaborator

niccolot Feb 23, 2024
Author

niccolot Feb 23, 2024
Author

smartalecH
Feb 29, 2024
Collaborator