Gyrotropic Fails on Multiple materials (process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

[mojv]

hi, i've been trying to use the new gyrotropic feature, but it seems not to work when used for 2 or more materials. the following example is just a simple modification of the faraday_rotation found the meep repository:

import meep as mp
from meep import materials as mt

b0 = 0.15 # magnitude of bias vector

## Set up and run the Meep simulation:
tmax = 100
L = 20.0
cell = mp.Vector3(0, 0, L)
fsrc, src_z = 0.8, -8.5
pml_layers = [mp.PML(thickness=1.0, direction=mp.Z)]

######## With this block works good ##########
BK7 = gyrotropic_conversion(mt.BK7, bias=mp.Vector3(0, 0, b0))
geometry = [
    mp.Block(material=BK7, size=mp.Vector3(mp.inf, mp.inf, L), center=mp.Vector3(0, 0, 0)),
]
#############################################

sources = [mp.Source(mp.ContinuousSource(frequency=fsrc),
                     component=mp.Ex, center=mp.Vector3(0, 0, src_z))]

sim = mp.Simulation(cell_size=cell, geometry=geometry, sources=sources,
                    boundary_layers=pml_layers, resolution=50)
sim.run(until=tmax)

## Plot results:
import numpy as np
import matplotlib.pyplot as plt

ex_data = sim.get_efield_x().real
ey_data = sim.get_efield_y().real

z = np.linspace(-L/2, L/2, len(ex_data))
plt.figure(1)
plt.plot(z, ex_data, label='Ex')
plt.plot(z, ey_data, label='Ey')
plt.xlim(-L/2, L/2); plt.xlabel('z')
plt.legend()


plt.figure(2)
plt.subplot(2,1,1)
plt.plot(z, ex_data, label='Ex (MEEP)')
plt.xlim(-L/2, L/2); plt.xlabel('z')
plt.legend(loc='lower right')

plt.subplot(2,1,2)
plt.plot(z, ey_data, label='Ey (MEEP)')
plt.xlim(-L/2, L/2); plt.xlabel('z')
plt.legend(loc='lower right')
plt.tight_layout()
plt.show()

Ths previous code runs good and the only different thing is the function "gyrotropic_conversion" i created for convert existing meep.materials in gyrotropic.

def gyrotropic_conversion (material, bias):
    material_suc = []
    for susceptibility in material.E_susceptibilities:
        if type(susceptibility).__name__ == 'DrudeSusceptibility':
            material_suc.append(
                mp.GyrotropicDrudeSusceptibility(frequency=susceptibility.frequency, gamma=susceptibility.gamma,
                                                 sigma=susceptibility.sigma_diag.x, bias=bias))
        elif type(susceptibility).__name__ == 'LorentzianSusceptibility':
            material_suc.append(
                mp.GyrotropicLorentzianSusceptibility(frequency=susceptibility.frequency, gamma=susceptibility.gamma,
                                                 sigma=susceptibility.sigma_diag.x, bias=bias))

    return mp.Medium(epsilon=material.epsilon_diag.x, E_susceptibilities=material_suc, valid_freq_range=material.valid_freq_range)

Now if I add an extra material to the main code like this:

######## With this block doesnt work ##########
BK7 = gy.gyrotropic_conversion(mt.BK7, bias=mp.Vector3(0, 0, b0))
AU = gy.gyrotropic_conversion(mt.Au, bias=mp.Vector3(0, 0, b0))
geometry = [
    mp.Block(material=BK7, size=mp.Vector3(mp.inf, mp.inf, L), center=mp.Vector3(0, 0, 0)),
    mp.Block(material=AU, size=mp.Vector3(mp.inf, mp.inf, 2), center=mp.Vector3(0, 0, 0)),
]
#############################################

The execution stops with the error

process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

or even if we just want to pretend an empty space with a small layer of any material (in this case BK7) like in the following, it won't also compile.

######## With this block doesnt work ##########
BK7 = gy.gyrotropic_conversion(mt.BK7, bias=mp.Vector3(0, 0, b0))
geometry = [
    mp.Block(material=BK7, size=mp.Vector3(mp.inf, mp.inf, 2), center=mp.Vector3(0, 0, 0)),
]
#############################################

best regards

[mojv]

The above error was thrown by the PyCharm console. but if I run the file in the Ubuntu shell I get the following:

Segmentation fault (core dumped)

And with gdb I get the following

Thread 1 "python" received signal SIGSEGV, Segmentation fault.
0x00007ffff5ff3667 in meep::fields::step_boundaries(meep::field_type) () from /home/ubuntu/miniconda3/envs/mp/lib/python3.7/site-packages/meep/../../../libmeep.so.16

[stevengj]

@seewhydee, do you have time to take a look at this?

[mojv]

@seewhydee @stevengj, Any news on this, it will be of great help and greatly appreciated.

[seewhydee]

Sorry for the delay, looking into it now.

[seewhydee]

OK, this seems to be related to the allocation of polarization (P) data in gyrotropic_susceptibility::needs_P. There are some comments in susceptibility.cpp (line 66) talking about the subtleties of when to allocate the P data, which I'm afraid I don't understand.

As an experiment, if we over-relax the rules on needs_P for gyrotropic media (so that the P data allocation triggers more easily), and add a couple of extra checks to update_P, then the segfault no longer seems to happen:

seewhydee@6f0928f

However, I am pretty sure this is not exactly the right thing to do. If anyone has some insight about the P allocation, it would be really helpful.

[seewhydee]

@stevengj, any insight?

[stevengj]

In general, needs_P needs to return true for any polarization component that is needed by the update equations.

Since the gyro_tensor couples every Cartesian component of P to every other component of P, and will crash if another component of P is NULL, it seems like this code needs to be something like:

  FOR_DIRECTIONS(d) {
    if (!trivial_sigma[direction_component(c, d)][d]) return true;
  }

That is, if you have a nonzero sigma diagonal for any direction of the current component (E or H), then it should return true. (This is regardless of whether the corresponding field W is allocated — your update_P code has a special case for when a W component is NULL, but still requires the corresponding P component.)

[stevengj]

Note that the trivial_sigma array is synchronized among all chunks/processes, so that it tells you if sigma ≠ 0 in any chunk. This is necessary in general because off-diagonal polarization terms may come from different Yee grid points that need to be communicated from other processors, and so we store P even on chunks that have sigma = 0 in case neighboring chunks need them. (In principle, we could eliminate this storage by just realizing that the neighboring P is zero in this case, but I haven't gotten around to implementing that logic.)

Looking at your update_P code, one thing I noticed is that you aren't doing any OFFDIAG averaging when computing the influence of other P components from off-diagonal gyro_tensor entries. Doesn't this give only first-order accuracy?

[seewhydee]

Actually, using the default susceptibility::needs_P works with some minor modifications to gyrotropic_susceptibility::update_P. I have opened pull request #1219 to discuss this.

Regarding off-diag averaging, averaging is done to extract the other components of the E fields, which enter into the precession/forcing terms of the equations of motion for P. It doesn't seem like it should be necessary to average over P. Maybe I'm misunderstanding?

[stevengj]

Looking back at the code, I see that we are fine for the reason you discussed in this comment:

To implement gyrotropic susceptibilities, we track three polarization components (e.g. Px, Py, Pz) on EACH of the Yee cell's three driving field positions (e.g., Ex, Ey, and Ez), i.e. 9 numbers per cell. This takes 3x the memory and runtime compared to Lorentzian susceptibility. The advantage is that during update_P, we can directly access the value of P at each update point without averaging.

If we weren't paying this factor of 3 in memory cost, then we would need a more complicated averaging scheme for off-diagonal P terms, but as it is you are right.

[mojv]

@stevengj @seewhydee, thanks a lot for the help, we manage to recreate a Au/Co/Au nanodisk with MO response and got results pretty similar to what was measured experimentally here