Wrong size of Block in simulation

[leSemaleon]

I want to obtain a set of spectra from a cross-shaped metasurface. I've encountered a problem: the geometry displayed by get_array doesn't match my expectations.

Crosses have three parameters: period, arm length, and arm width, P, L, W.

I'm constructing the geometry based on these values, but get_array shows that the crosses are always symmetrical about 0,0, but the point 0,0 itself isn't in the eps array. In my opinion, for a grid size that's a multiple of two, 00 shouldn't be included, but for a grid size that's not a multiple of two, 00 should be included.

Currently, all dimensions are padded so that 00 doesn't appear, and accordingly, I don't get the size I expect. Instead of 5 µm, the step becomes 10, and the 60 µm I specify becomes 70 µm on the get_array

In the figures below, you can see that increasing the period from 300 to 305, that is, by one grid step, does not change anything. The grid size remains at 62 pixels, which is equal to 310 µm.

A similar situation occurs with other parameters.

The code is hidden under the spoiler

Code

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from datetime import datetime
import matplotlib.pyplot as plt
import meep as mp
import numpy as np
import materials
from tqdm import tqdm
from matplotlib.ticker import MultipleLocator


NFREQ = 2000


def plot_epsdata(sim: mp.Simulation, pos, title):
    sim.init_sim()
    eps_data = sim.get_array(
        center=mp.Vector3(0, 0, pos),
        size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0),
        component=mp.Dielectric,
    )
    eps_data1 = sim.get_array(
        center=mp.Vector3(0, 0, 0),
        size=mp.Vector3(sim.cell_size.x, 0, sim.cell_size.z),
        component=mp.Dielectric,
    )
    eps_data2 = sim.get_array(
        center=mp.Vector3(0, 0, 0),
        size=mp.Vector3(0, sim.cell_size.y, sim.cell_size.z),
        component=mp.Dielectric,
    )
    # plt.ion()
    _, ax = plt.subplots()
    ax.pcolormesh(eps_data, edgecolor="gray", linewidth=0.5)
    ax.set_aspect("equal")
    ax.set_title(title)
   
    plt.savefig("eps_" + datetime.now().strftime("%H:%M:%S") + ".svg")
    plt.show()

    plt.pause(0.5)


def simulation(period, widthCross, widthShouder, iteration):
    f_min_thz = 0.21  # source min frequency THz
    f_max_thz = 1.2  # source max frequency THz
    f_min = f_min_thz / 300
    f_max = f_max_thz / 300
    f_cen = 0.5 * (f_min + f_max)
    df = f_max - f_min

    gridstep = 5

    dpml = int((1 / f_min) * 0.6 / gridstep + 0.5) * gridstep
    dAirIn = 500
    dMetal = 10
    dAirOut = 500

    layers = [
        dpml,
        dAirIn,
        dMetal,
        dAirOut,
        dpml,
    ]

    accuracityOfDecay = 1e-7

    sx = period
    sy = period
    sz = sum(layers)

    field_component = mp.Ey

    cell_size = mp.Vector3(sx, sy, sz)

    resolution = 1 / gridstep

    nfreq = NFREQ  # number of frequencies at which to compute flux

    sources = [
        mp.Source(
            mp.GaussianSource(f_cen, fwidth=df, is_integrated=True),
            component=field_component,
            center=mp.Vector3(0, 0, -0.5 * sz + dpml),
            size=mp.Vector3(sx, sy, 0),
            # amplitude=100,
        )
    ]

    boundary_layers = [
        mp.PML(
            thickness=dpml, direction=mp.Z, side=mp.Low, pml_profile=lambda u: u * u * u
        ),
        mp.PML(
            thickness=dpml,
            direction=mp.Z,
            side=mp.High,
            pml_profile=lambda u: u * u * u,
        ),
    ]

    # default_material = mp.Medium(epsilon=1, D_conductivity=1e-7)
    default_material = mp.air

    geometry = []
    rel_z = -0.5 * sz + dpml + dAirIn
    if iteration:
        meta = [
            mp.Block(
                size=mp.Vector3(sx, sy, dMetal),
                center=mp.Vector3(0, 0, rel_z + 0.5 * dMetal),
                # material=mp.metal,
                material=materials.Copper(),
            ),
            mp.Block(
                size=mp.Vector3(widthCross, widthShouder, dMetal),
                center=mp.Vector3(
                    0,
                    0,
                    rel_z + 0.5 * dMetal,
                ),
                material=default_material,
            ),
            mp.Block(
                size=mp.Vector3(widthShouder, widthCross, dMetal),
                center=mp.Vector3(
                    0,
                    0,
                    rel_z + 0.5 * dMetal,
                ),
                material=default_material,
            ),
        ]
        posMetal = rel_z + 0.5 * dMetal
        rel_z += dMetal
        geometry.extend(meta)

    sim = mp.Simulation(
        cell_size=cell_size,
        boundary_layers=boundary_layers,
        sources=sources,
        geometry=geometry,
        default_material=default_material,
        k_point=mp.Vector3(),
        resolution=resolution,
        eps_averaging=False,
    )

    if iteration and 1:
        plot_epsdata(
            sim, posMetal, title=f"P: {period}, L: {widthCross}, W: {widthShouder}"
        )


def main():
    P = [305]
    L = [130]
    W = [65]

    for i in tqdm(range(len(P)), ncols=50, disable=True):
        for j in tqdm(range(len(L)), ncols=50, disable=True):
            for k in tqdm(range(len(W)), ncols=50, disable=True):
                simulation(P[i], L[j], W[k], 1)
                continue


if __name__ == "__main__":
    main()