Total dft_chunks size does not agree with size allocated for output array

[zebihcetin]

Hello everyone,

I’m encountering the following error with Meep:

meep: Total dft_chunks size does not agree with size allocated for output array.
Abort(1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
.
.
.
meep: Total dft_chunks size does not agree with size allocated for output array.
Abort(1) on node 447 (rank 447 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 447

The issue arises in a Python script that first runs a simulation without a structure. After calling reset_meep() to reset the simulation environment, the script sets up a structure and runs the simulation again. The first part of the script works as expected. However, after the first part completes, Meep aborts with the error mentioned above.

The code runs without issues at low resolution, but the error occurs when the resolution is increased. The only parameter that changes is the resolution; everything else remains the same. I am running Meep on 4 nodes, each with 112 cores and 256 GB of RAM.

Has anyone encountered a similar issue or have any suggestions for resolving this

[stevengj]

That seems like a bug. Can you post a minimal script that triggers it? (i.e. delete as much as possible from your script)

[zebihcetin]

Hello @stevengj,

Here is the minimal script and the command line arguments that triggers the error.

mpirun -np 448 python3 scripy.py -resol 48 -PML -freq_ind

It works without any error with -resol 32.

(MPI version is 4.1 and MEEP version is 1.28.0)

import numpy as np
np.finfo(np.dtype("float32"))
np.finfo(np.dtype("float64"))
import meep as mp
import argparse
from meep.materials import Al

def main(args):
    dz=0.5
    resol = args.resol
    lat_const = 1
    run_time = 2
    num_layers = 1
    dpmlx = 1
    dpmly = 1
    dpmlz = 3
    epsa, epsb = 3.00, 1
    nx, ny, nz = 4, 4, 3
    dpadz = nz*lat_const

    # Computational Dim and cell_size-------------------------------------------------------------------------
    if args.PBC:
        sx, sy, sz = nx*lat_const, ny*lat_const, num_layers*lat_const + 2*dpadz + 2*dpmlz
    else:
        sx, sy, sz = nx*lat_const + 2*dpmlx, ny*lat_const + 2*dpmly, num_layers*lat_const + 2*dpadz + 2*dpmlz

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

    # Geometry---------------------------------------------------------------------------------------------------------
    geometry = []

    # Sources----------------------------------------------------------------------------------------------------------
    # fmin, fmax = Al.valid_freq_range[0], Al.valid_freq_range[1]
    lmin = 0.25
    lmax = 45.0
    fmin = 1.0/lmax
    fmax = 1.0/lmin
    fcen = 0.5*(fmin + fmax)
    df = fmax - fmin
    nfreqs = 9000
    print('fmax',fmin, 'fmin', fmax, 'fcen', fcen)
    courant = 0.4

    sources = [mp.Source(src=mp.GaussianSource(frequency=fcen, fwidth=df, cutoff=10, is_integrated=False if args.PBC else True),
                         center=mp.Vector3(0, 0, -0.5 * sz + dpmlz + 0.50),
                         size=mp.Vector3(sx, sy, 0),
                         component=mp.Ex,
                         amplitude=1)]

    # Boundary Layers, PML or Absorber--------------------------------------------------------------------------------------
    if args.PML:
        pml_layers = [mp.PML(dpmlx, direction=mp.X, pml_profile=lambda u: u * u * u),
                      mp.PML(dpmly, direction=mp.Y, pml_profile=lambda u: u * u * u),
                      mp.PML(dpmlz, direction=mp.Z, pml_profile=lambda u: u * u * u)]
    else:
        pml_layers = [mp.Absorber(dpmlx, direction=mp.X, pml_profile=lambda u: u * u * u),
                      mp.Absorber(dpmly, direction=mp.Y, pml_profile=lambda u: u * u * u),
                      mp.Absorber(dpmlz, direction=mp.Z, pml_profile=lambda u: u * u * u)]

    # Simulation Class-------------------------------------------------------------------------------------------------
    fname_prefix = "norm-run"
    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        k_point=mp.Vector3(x=0.0, y=0.0, z=0.0) if args.PBC else False,
                        geometry=geometry,
                        sources=sources,
                        resolution=resol,
                        filename_prefix=fname_prefix,
                        Courant=courant)
    # Added to set boundary metallic-----------------------------------------------------------------------------------
    # sim.set_boundary(mp.ALL, mp.ALL, mp.Metallic)

    # Reflected flux---------------------------------------------------------------------------------------------------
    refl_flux_region = mp.FluxRegion(center=mp.Vector3(0, 0, -0.5*sz+dpmlz+1.0), size=mp.Vector3(sx, sy, 0))
    refl_flux_monitor = sim.add_flux(fcen, df, nfreqs, refl_flux_region)

    # Output flux in Z direction---------------------------------------------------------------------------------------
    tran_flux_region = mp.FluxRegion(center=mp.Vector3(0, 0, 0.5 * sz - dpmlz - 0.5), size=mp.Vector3(sx, sy, 0), direction=mp.Z)
    tran_flux_monitor = sim.add_flux(fcen, df, nfreqs, tran_flux_region)

    # Run--------------------------------------------------------------------------------------------------------------
    decay_by = np.power(10, float(-run_time))  # decay_by = 1e-12
    dt = 20
    sim.run(mp.at_beginning(mp.output_epsilon),until_after_sources=mp.stop_when_energy_decayed(dt=dt, decay_by=decay_by))

    # get output flux-------------------------------------------------------------------------------------------------
    empt_tran_flux = mp.get_fluxes(tran_flux_monitor)

    # Reflected flux with empty simulation region---------------------------------------------------------------------
    empt_refl_flux = mp.get_fluxes(refl_flux_monitor)

    # Reflected data, will be used in the next run with geometry-------------------------------------------------------
    empt_refl_data = sim.get_flux_data(refl_flux_monitor)

    # print fluxes to terminal-----------------------------------------------------------------------------------------
    sim.display_fluxes(refl_flux_monitor, tran_flux_monitor)

    #
    sim.reset_meep()  # -------------------------------------------------------------------------------------------------
    #

    # Geometry---------------------------------------------------------------------------------------------------------
    mymat = mp.Medium(epsilon=epsa) if args.freq_ind else Al #mymaterial

    geometry = [mp.Block(size=mp.Vector3(mp.inf, mp.inf, dz),
                         material=mymat,
                         center=mp.Vector3(x=0, y=0, z=0))]

    # Simulation ------------------------------------------------------------------------------------------------------
    fname_prefix = "str-run"
    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        k_point=mp.Vector3(x=0.0, y=0.0, z=0.0) if args.PBC else False,
                        geometry=geometry,
                        sources=sources,
                        resolution=resol,
                        filename_prefix=fname_prefix,
                        Courant=courant)

    # Reflected flux---------------------------------------------------------------------------------------------------
    refl_flux_monitor = sim.add_flux(fcen, df, nfreqs, refl_flux_region)

    # Output flux in Z direction---------------------------------------------------------------------------------------
    tran_flux_monitor = sim.add_flux(fcen, df, nfreqs, tran_flux_region)

    # Load negatif flux data to refl_flux_monitor----------------------------------------------------------------------
    sim.load_minus_flux_data(refl_flux_monitor, empt_refl_data)

    # Run--------------------------------------------------------------------------------------------------------------
    sim.run(mp.at_beginning(mp.output_epsilon), until_after_sources=mp.stop_when_energy_decayed(dt=dt, decay_by=decay_by))

    # output flux------------------------------------------------------------------------------------------------------
    flux_freqs = mp.get_flux_freqs(tran_flux_monitor)

    # get tran flux
    geom_tran_flux = mp.get_fluxes(tran_flux_monitor)

    # get refl flux
    geom_refl_flux = mp.get_fluxes(refl_flux_monitor)

    # print fluxes to terminal-----------------------------------------------------------------------------------------
    sim.display_fluxes(refl_flux_monitor, tran_flux_monitor)

    mp.all_wait()
    # ------------------------------------------------------------------------------------------------------------------
    if mp.am_master():
        filename = "fluxes.dat"
        with open(filename, 'wb') as f:
            for i in range(nfreqs):
                data = np.column_stack((flux_freqs[i], -geom_refl_flux[i]/empt_tran_flux[i], geom_tran_flux[i]/empt_tran_flux[i]))
                np.savetxt(f, data, delimiter=', ', newline='\n')

# ---------------------------------------------------------------------------------------------------------------------
if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-resol',  type=int, default=64, help='resolution of simulation (default: 64)')

    parser.add_argument('-PBC',      action='store_true', default=False, help='True if empty (default: False)')
    parser.add_argument('-PML',      action='store_true', default=False, help='True if empty (default: False)')
    parser.add_argument('-freq_ind', action='store_true', default=False, help='True if empty (default: False)')
    args = parser.parse_args()
    main(args)