load_minus_flux_data() crashing 2D simulations with 3D geometry?

[FrockSaints]

I'm trying to simulate nanoparticles in a coating.
3D seems to be working, although it's surprisingly slow (months/years when it would take a few hours in 2D), but 2D crashes when using load_minus_flux_data(). I've tried commenting out that line and it works, but I don't understand why it would cause a crash.

More precisely, an initial run is done just with air, for normalization purposes, then another one with the coating plus the nanoparticles.
According to the documentation, "In a 2d calculation, only the intersections of the objects with the xy plane are considered.". This would be ideal, but for some reason it's crashing with anything other than spheres centered at z=0.
Spheres are easy to fix, as a sphere offset from z=0 has, for the 2D simulation, an easy equivalent centered at z=0. I've already tested the 2D simulation works fine in this case. But implementing equivalents for the other objects, with arbitrary offsets and rotations, is a bit of a hassle.

Two easier alternatives would be to specify the 3D geometry, grab the z=0 slice with get_array() and initialize the 2D simulation with epsilon_input_file, or to simply normalize it by hand (probably ideal, since it's just air, but haven't looked into it for a generic source yet), but now I'm curious regarding the crash.

Advice on how to proceed and insight into this problem are much appreciated.
Apologies if the exposition was long-winded and thank you for your time.

[smartalecH]

Can you post some (minimal) code? Have you tried debugging the crash (e.g. using gdb)?

[FrockSaints]

Here's a trimmed version that still crashes the jupyter kernel:

import meep as mp
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rd

def sim_init(sim_params, geometry=[]):
    
    eps_coat, resolution, sx, thicc, sh, dpml, sz, time, nfreq = sim_params
    
    pml_layers = [mp.PML(thickness=dpml, direction=mp.Y)]
    sy = thicc+sh+2*dpml  # size of cell in Y direction
    cell = mp.Vector3(sx,sy,sz)

    wvl_min = 0.3           # min wavelength
    wvl_max = 2.5           # max wavelength
    fmin = 1/wvl_max        # min frequency
    fmax = 1/wvl_min        # max frequency
    fcen = 0.5*(fmin+fmax)  # center frequency
    df = fmax-fmin          # frequency width

    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                         component=-mp.Ey,
                         center=mp.Vector3(0,-sy/2+dpml,0),
                         size=mp.Vector3(sx,0,sz))]
    
    sim = mp.Simulation(cell_size=cell,
                        dimensions=3,
                        boundary_layers=pml_layers,
                        geometry=geometry,
                        sources=sources,
                        resolution=resolution,
                        k_point=mp.Vector3())
    
    return sy, fcen, df, sources, sim

def air_flux(sim_params):
    
    sy, fcen, df, sources, sim = sim_init(sim_params)
    eps_coat, resolution, sx, thicc, sh, dpml, sz, time, nfreq = sim_params
    
    # incident flux
    inc_fr = mp.FluxRegion(center=mp.Vector3(0,-sy/2+dpml+sh,0), size=mp.Vector3(sx,0,sz))
    inc = sim.add_flux(fcen, df, nfreq, inc_fr)

    sim.run(until=time)
    
    inc_data = sim.get_flux_data(inc)
    inc_flux = mp.get_fluxes(inc)
    
    return inc_data, inc_flux
    
def fluxes(sim_params, part_params, crash):
    
    eps_coat, resolution, sx, thicc, sh, dpml, sz, time, nfreq = sim_params
    
    inc_data, inc_flux = air_flux(sim_params)
    
    geometry = [mp.Block(size=mp.Vector3(mp.inf,thicc,mp.inf),
                         center=mp.Vector3(0,sh/2,0),
                         material=mp.Medium(epsilon=eps_coat))]
    
    len1, eps, offset = part_params
    center = (rd.random()-0.5)*sx, sh/2 + (rd.random()-0.5)*thicc, crash*offset
    geometry.append(mp.Sphere(center=mp.Vector3(*center), radius=len1/2, material=mp.Medium(epsilon=eps)))
    
    sy, fcen, df, sources, sim = sim_init(sim_params, geometry)
   
    # reflected flux
    refl_fr = mp.FluxRegion(center=mp.Vector3(0,-sy/2+dpml+sh,0), size=mp.Vector3(sx,0,sz))
    
    refl = sim.add_flux(fcen, df, nfreq, refl_fr)

    # transmitted flux
    tran_fr = mp.FluxRegion(center=mp.Vector3(0,sy/2-dpml,0), size=mp.Vector3(sx,0,sz))
    tran = sim.add_flux(fcen, df, nfreq, tran_fr)
    
    # for normal run, load negated fields to subtract incident from refl. fields
    sim.load_minus_flux_data(refl, inc_data)

    sim.run(until=time)

    flux_freqs = np.array(mp.get_flux_freqs(refl))
    refl_flux = np.array(mp.get_fluxes(refl))
    tran_flux = np.array(mp.get_fluxes(tran))
    
    return inc_flux, flux_freqs, refl_flux, tran_flux

def spectrum(inc_flux, flux_freqs, refl_flux, tran_flux, plot=True):

    # The reflectance is the reflected flux divided by the incident flux.         
    wl = 1/flux_freqs
    Rs = -refl_flux/inc_flux
    Ts= tran_flux/inc_flux
    
    if plot==True:
        plt.figure()
        plt.title('spectrum')
        plt.plot(wl,Rs,'bo-',label='reflectance')
        plt.plot(wl,Ts,'ro-',label='transmittance')

        plt.axis([0.3, 2.5, -0.1, 1.1])
        plt.xlabel("wavelength (μm)")
        plt.legend(loc="upper right")
    
    return wl, Rs

def main(sim_params, part_params, crash):
    
    # eps_coat, resolution, sx, thicc, sh, dpml, sz, time, nfreq = sim_params
    # len1, eps, offset = part_params
        
    inc_flux, flux_freqs, refl_flux, tran_flux = fluxes(sim_params, part_params, crash)
    
    wl, Rs = spectrum(inc_flux, flux_freqs, refl_flux, tran_flux)

main(sim_params=[2.5, 20, 5, 5, 1, 5, 0, 2000, 201], part_params=[0.2, 9, 0.09], crash=True)

I haven't tried debugging with gdb (didn't even know about it). My knowledge of C is limited at best but I'll definitely take a look when I can. Please let me know if you notice some crass mistake.

[FrockSaints]

Apologies if the code wasn't minimal enough - definitely could have simplified it further.
It's been a while, as I had some issues getting gdb to work at the time and then had to move on, but I've managed to get it to work this time, and got this back:
Thread 1 "python" received signal SIGSEGV, Segmentation fault. 0x00007f62ded7d1b0 in meep::timing_scope::exit() () from /home/frock/conda/envs/mp/lib/python3.10/site-packages/meep/../../../libmeep.so.28
Why this is happening is beyond me - it seems weird to segfault on exit(), but after a quick search, there are apparently many reasons for it. Found an answer on stackoverflow recommending Valgrind to debug it. Might try that in the future (and update the findings here) but that's another tool I'm not familiar with and I've found a suitable alternative for my use case in the meantime, so I'll leave it at this for now.

Any insight is greatly appreciated nonetheless.