parallel_sim.py

import numpy as np
import meep as mp


def simulate_G(comp, src_pos, res, size, phi, fcen, eps, df=0.1, pml=0.7):
    """Simulate and calculate the comp-column of the Greens tensor.
     Args:
        size (mp.Vector3): size of the simulation box
        src_pos (mp.Vector3): position of the gaussian source current
        fcen (float): frequency of the source
        df (float): frequency width of the gaussian
        res (int): resolution of the simulation
        pml (float): size of perfectly matching layer
    Returns:
        np.array: 2d array of the column of the Greens Tensor.
    """
    geom = _phi_to_geom(phi, res, size, eps)
    gSrc = mp.GaussianSource(fcen, fwidth=df)
    pml = [mp.PML(pml)]
    src = [mp.Source(gSrc, component=comp, center=src_pos)]
    sim = mp.Simulation(cell_size=size,
                        sources=src, resolution=res, boundary_layers=pml,
                        force_complex_fields=True, force_all_components=True,
                        geometry=geom, eps_averaging=False)
    dft_obj = sim.add_dft_fields(
        [mp.Ex, mp.Ey, mp.Ez], fcen, fcen, 1,
        center=mp.Vector3(), size=size)
    sim.run(until_after_sources=100)
    jw = gSrc.fourier_transform(fcen)
    ex = sim.get_dft_array(dft_obj, mp.Ex, 0).transpose()
    ey = sim.get_dft_array(dft_obj, mp.Ey, 0).transpose()
    ez = sim.get_dft_array(dft_obj, mp.Ez, 0).transpose()
    ex = _interpolate_E(ex, mp.Ex, res)
    ey = _interpolate_E(ey, mp.Ey, res)
    ez = _interpolate_E(ez, mp.Ez, res)
    G_col = _calc_G(comp, ex, ey, ez, jw, fcen, res, size)
    return G_col


def _calc_G(comp, ex, ey, ez, jw, fcen, res, size):
    """Calculate the components of G generated by j_comp.
    Args:
        comp (mp.Const): mp.Ex, mp.Ey, mp.Ez
        jw (np.array): 2d array of the fourier transformed source current
    Returns:
        np.array: 2d array of the Green's Tensor.
    """
    G = np.zeros((3, int(res*size.y), int(res*size.x)), dtype=complex)
    G[0] = ex/(1j*2*np.pi*fcen*jw)
    G[1] = ey/(1j*2*np.pi*fcen*jw)
    G[2] = ez/(1j*2*np.pi*fcen*jw)
    return G


def r_cellmiddle(res, size):
    """Calculate the position values for the middle of the cells.
    In the format r[y,x] = (x,y,z)."""
    r = np.zeros((int(size.y)*res, int(size.x)*res, 3))
    for y in range(len(r)):
        for x in range(len(r[0])):
            r[y, x] = np.array([x/res-size.x/2+1/(2*res),
                                y/res-size.y/2+1/(2*res), 0])
    return r


def _interpolate_E(E, comp, res):
    """Linear interpolation of E to the middle of the cell.
    Args:
        comp (mp.Const): mp.Ex, mp.Ey, mp.Ez
    Returns:
        np.arrays: position r and the E-field
    """
    centE = np.zeros_like(E)
    if comp == mp.Ex:
        centE[0: -1, :] = (E[0: -1, :] + E[1:, :])/2
        m = (E[-1, :] - E[-2, :]) * res
        centE[-1, :] = m/(2*res)
    elif comp == mp.Ey:
        centE[:, 0: -1] = (E[:, 0: -1] + E[:, 1:])/2
        m = (E[:, -1] - E[:, -2])*res
        centE[:, -1] = m/(2*res)
    elif comp == mp.Ez:
        centE[: -1, : -1] = (E[: -1, : -1] + E[1:, 1:])/2
        E[-1, :] = E[-2, :]
        E[:, -1] = E[:, -2]
    return centE


def _phi_to_geom(phi, res, size, eps):
    """Create a geometry from the level set function phi (array).
    Where phi is negative, matter will be placed.
    Args:
        phi (np.array): 2d array representing the level set function
    Returns:
        list: list of the geometries
    """
    geom = []
    r = r_cellmiddle(res, size)
    size = mp.Vector3(1.0001/res, 1.0001/res, mp.inf)
    for y in range(len(phi)):
        for x in range(len(phi[0])):
            if phi[y, x] <= 0:
                cent = mp.Vector3(r[y, x][0], r[y, x][1])
                geom.append(mp.Block(size, center=cent,
                                     material=mp.Medium(epsilon_diag=eps)))
    return geom


def sim_to_file(direct, comp, src_pos, size, res, pml, phi, fcen, eps=mp.Vector3(12, 12, 12)):
    """Simulate continuous source and write to a hdf5 file.
    Args:
        filename (str): filename
        comp (mp.Const): component of the field to simulate
    Returns:
        None.
    """
    pml = [mp.PML(pml)]
    geom = _phi_to_geom(phi, res, size, eps)
    if direct[-1] == '/':
        direct = direct[:-1]
    src = [mp.Source(gSrc=mp.GaussianSource(fcen, fwidth=0.1),
                     component=comp,
                     center=src_pos)]
    sim = mp.Simulation(cell_size=size,
                        sources=src, resolution=res, boundary_layers=pml,
                        force_complex_fields=True, force_all_components=True,
                        geometry=geom, eps_averaging=False)

    A = '{:}/analyze-eps-000000.00.h5'.format(direct)
    a = 'yarg:0.5'
    out = mp.output_png(
        comp, '-Zc dkbluered -m 0.5 -M 0.5 -A {:} -a {:}'.format(A, a))
    sim.use_output_directory(dname=direct)
    sim.run(mp.at_beginning(mp.output_epsilon),
            mp.at_every(0.3, out),
            until=200)