reset_meep() does not fully reset shared C-level state, causing incorrect MPB ModeSolver results in subsequent calls

[giannischatzo]

When running MEEP's Simulation.init_sim() followed by reset_meep(), subsequent MPB ModeSolver.find_k() calls in the same process return incorrect effective index values. The neff drops from the expected guided mode value (~2.54 for a Si slab at 1550nm) to the cladding index (~1.45), suggesting MPB is no longer seeing the correct simulation setup.
This makes it impossible to reliably alternate between MEEP simulations and MPB mode solving in a loop within a single process — a common pattern in photonic simulation workflows (e.g. computing effective indices for 2D simulations, parameter sweeps, etc.).
Steps to Reproduce

1. Run MPB ModeSolver.find_k() for a 1D Si/SiO2 slab waveguide → neff ≈ 2.5445 (correct)
2. Create a 2D MEEP Simulation with MaterialGrid, non-zero geometry_center, and EigenModeSource; call init_sim() then reset_meep()
3. Run the same MPB ModeSolver.find_k() again with identical parameters → neff ≈ 1.4484 (incorrect, equals the cladding index)
   The corruption happens on the very first iteration after init_sim() and persists for all subsequent MPB calls.
   Suspected Cause
   We believe this is caused by MEEP and MPB sharing process-wide C-level libctl global variables (dimensions, default_material, geometry_center, etc.) declared in ctlgeom-types.h via the shared libctlgeom library. init_sim() writes MEEP's 2D simulation values into these globals, but reset_meep() only cleans up Python-level state (self.fields, self.structure, etc.) without restoring the C-level globals to their defaults. MPB's ModeSolver constructor resets some of these globals but apparently not all of them, so subsequent MPB runs inherit stale state from the MEEP simulation.
   The following MEEP simulation features appear necessary to trigger the issue:

• 2D cell (cell_size.z == 0)
• MaterialGrid inside a mp.Block
• Non-zero geometry_center
• EigenModeSource with eig_kpoint and eig_parity
  Minimal Reproducible Example
  Tested with pymeep 1.31.0 (conda-forge, nompi), Python 3.13.12, numpy 2.4.2. Only depends on meep and numpy.

import numpy as np
import meep as mp
import meep.mpb as mpb

# --- Material constants (Si/SiO2 at 1550nm) ---
N_SI = 3.4807
N_SIO2 = 1.4484
N_EFF_2D = 2.8528
CORE_HEIGHT = 0.22
BOX_HEIGHT = 2.0
TOX_HEIGHT = 0.7
WAVELENGTH = 1.55
CENTRAL_FREQ = 1 / WAVELENGTH

def run_mpb(core_height=0.15):
    """Run MPB find_k for a 1D slab waveguide (calc_effective_index mode)."""
    mp.verbosity(0)
    sc_z = core_height + BOX_HEIGHT + TOX_HEIGHT
    geometry_lattice = mp.Lattice(size=(0, 0, sc_z))
    core_material = mp.Medium(index=N_SI)
    bottom_material = mp.Medium(index=N_SIO2)
    top_material = mp.Medium(index=N_SIO2)
    geometry = [
        mp.Block(
            size=mp.Vector3(mp.inf, mp.inf, BOX_HEIGHT),
            center=mp.Vector3(z=-core_height / 2 - BOX_HEIGHT / 2),
            material=bottom_material,
        ),
        mp.Block(
            size=mp.Vector3(mp.inf, mp.inf, core_height),
            center=mp.Vector3(),
            material=core_material,
        ),
    ]
    ms = mpb.ModeSolver(
        geometry_lattice=geometry_lattice,
        geometry=geometry,
        default_material=top_material,
        resolution=100,
        num_bands=3,
        verbose=False,
    )
    f_mode = 1 / WAVELENGTH
    kz = ms.find_k(
        mp.NO_PARITY, f_mode, 1, 3, mp.Vector3(1),
        1e-7, f_mode * N_SI**2, f_mode * 0.1, f_mode * 4.0,
    )
    return [k / f_mode for k in kz]

def run_meep_init_sim():
    """Create and init a 2D MEEP simulation, then reset."""
    mp.verbosity(0)
    size_x, size_y = 25.0, 10.0
    cell_size = mp.Vector3(size_x, size_y, 0)
    cx, cy = 11.5, 0.0
    geometry_center = mp.Vector3(x=cx, y=cy)
    background_material = mp.Medium(index=N_SIO2)
    waveguide_material = mp.Medium(index=N_EFF_2D)
    grid_nx, grid_ny = 482, 191
    weights = np.zeros((grid_nx, grid_ny), dtype=np.float64)
    wg_width_pixels = max(1, int(0.45 / size_y * grid_ny))
    cy_px = grid_ny // 2
    weights[:, cy_px - wg_width_pixels // 2 : cy_px + wg_width_pixels // 2 + 1] = 1.0
    grid = mp.MaterialGrid(
        mp.Vector3(grid_nx, grid_ny), background_material, waveguide_material,
        weights=weights, grid_type="U_DEFAULT", do_averaging=False,
    )
    geometry = [mp.Block(center=mp.Vector3(x=cx, y=cy), size=mp.Vector3(size_x, size_y), material=grid)]
    source = mp.EigenModeSource(
        src=mp.GaussianSource(frequency=CENTRAL_FREQ, fwidth=0.2 * CENTRAL_FREQ),
        size=mp.Vector3(0, 3), center=mp.Vector3(x=cx - size_x / 2 + 1.0, y=cy),
        eig_band=1, eig_parity=mp.EVEN_Z, eig_kpoint=mp.Vector3(1, 0, 0),
        direction=mp.NO_DIRECTION,
    )
    pml_layers = [mp.PML(0.8, direction=mp.X), mp.PML(0.8, direction=mp.Y)]
    sim = mp.Simulation(
        cell_size=cell_size, boundary_layers=pml_layers, geometry=geometry,
        sources=[source], symmetries=[], default_material=background_material,
        resolution=19, split_chunks_evenly=True, geometry_center=geometry_center,
        eps_averaging=False,
    )
    sim.init_sim()
    sim.reset_meep()
    del sim

# --- Reproduce the bug ---
baseline = run_mpb()
print(f"Baseline neff: {baseline[0]:.4f}")  # Expected: ~2.5445
for i in range(3):
    run_meep_init_sim()
    result = run_mpb()
    neff = result[0] if result else float("nan")
    print(f"Iteration {i}: neff = {neff:.4f}")  # Actual: ~1.4484 (WRONG)

Output:
Baseline neff: 2.5445
Iteration 0: neff = 1.4484
Iteration 1: neff = 1.4484
Iteration 2: neff = 1.4484
Environment

• pymeep: 1.31.0 (conda-forge, nompi)
• mpb: 1.12.0
• libctl: 4.5.1
• Python: 3.13.12
• numpy: 2.4.2
• OS: Linux (Ubuntu)

[giannischatzo]

Running something like this before the second call to MPB seems to fix it

  _tmp = mp.Simulation(cell_size=mp.Vector3(1, 1, 1), resolution=2)
  mp.init_libctl(
      mp.Medium(),
      True,
      _tmp._create_grid_volume(False),
      mp.Vector3(),
      mp.Vector3(),
      [],
  )
  del _tmp