3D Maxwell Equation

Mathematical description

The PDEs is given by

$\begin{aligned} \nabla \cdot \bold{E} &= 0\\ \nabla \cdot \bold{H} &= 0\\ \nabla \times \bold{E} &= -\frac{\mu}{c} \frac{\partial \bold{H}}{\partial t} \\ \nabla \times \bold{H} &= \frac{\varepsilon}{c} \frac{\partial \bold{E}}{\partial t} \end{aligned}$

where $c＝299792458.0$ is the speed of light in vacuum, $\epsilon=10$ is the relative permittivity, and $\mu=1$ is the relative permeability. This formulation has omitted the permittivity and permeability in the free space by using $E\leftarrow \sqrt{{{\varepsilon }_{0}}}E,H\leftarrow \sqrt{{{\mu }_{0}}}H$ ，as described in the fdtd document.

The domain is defined as

$(\bold{x} \times t) \in [0, 1.6^{-5}]^3 \times [0, 1.91^{-13}]$

Dataset file information

File hierarchy:

/
|-- 0001 (Group)
|   |-- data (Dataset)
|   |-- global_maximums (Dataset)
|   |-- grid (Group)
|       |-- t (Dataset)
|       |-- x (Dataset)
|       |-- y (Dataset)
|       |-- z (Dataset)
|-- ...
|-- 0999 (Group)
|   |-- data (Dataset)
|   |-- global_maximums (Dataset)
|   |-- grid (Group)
|       |-- t (Dataset)
|       |-- x (Dataset)
|       |-- y (Dataset)
|       |-- z (Dataset)

Datasets information:

• data
• Shape: (8, 32, 32, 32, 6)
• Data Type: float
• Description: one solution sample of shape (Nt, Nx, Ny, Nz, Nc), with the channels corresponding to [Ex,Ey,Ez,Hx,Hy,Hz]
• global_maximums
• Shape: (6,)
• Data Type: float
• Description: the global maximums of each channels over all samples for normalization.
• t
• Shape: (8,)
• Data Type: float
• Description: temporal coordinates
• x
• Shape: (32,)
• Data Type: float
• Description: x-axis spatial coordinates
• y
• Shape: (32,)
• Data Type: float
• Description: y-axis spatial coordinates
• z
• Shape: (32,)
• Data Type: float
• Description: z-axis spatial coordinates

Data generation method

This dataset is from the following paper

@article{brandstetter2022clifford,
  title={Clifford neural layers for PDE modeling},
  author={Brandstetter, Johannes and Berg, Rianne van den and Welling, Max and Gupta, Jayesh K},
  journal={arXiv preprint arXiv:2209.04934},
  year={2022}
}

and the generation code is originated from https://microsoft.github.io/pdearena/data/

The data is generated by ftdt. A cubic space ($64\times64\times64$ grid) of the length 3.2e-5m was set up, where observations is made in the central subspace （a $32\times32\times32$ grid). A total of 18 emission sources with a wavelength of 1.0e-5m were randomly placed near the x-y, x-z, and y-z planes (6 for each) outside the observation subspace. The time step of the ftdt simulation is determined by the Courant number. Here, the snapshot is sampled every 25 time steps, with a total of eight snapshots taken, resulting in a temporal-spatial resolution of $8\times32\times32\times32$ .

Detailed data generation script is below:

import fdtd
from joblib import Parallel, delayed
import tqdm


class Maxwell3D(PDEConfig):
    wavelength: float = 1.0e-5
    sol: float = 299_792_458.0
    amplitude: float = 1.0
    permittivity: float = 10.0
    permeability: float = 1.0
    L: float = 3.2e-5
    n: int = 32
    n_large: int = 64
    nt: int = 12
    skip_nt: int = 250
    sample_rate: int = 15
    device: str = "cpu"

    def __repr__(self):
        return "Maxwell3D"

    @property
    def trajlen(self):
        return self.nt

    @property
    def grid_size(self):
        return (self.nt, self.n, self.n, self.n)

    @property
    def spatial_grid_size(self):
        return (self.n, self.n, self.n)

    @property
    def grid_spacing(self):
        return self.L / self.n_large

pde = Maxwell3D()

def genfunc(idx, s):
        rng = np.random.RandomState(idx + s)
        # Initialize grid and light sources
        grid = fdtd.Grid(
            (pde.L, pde.L, pde.L),
            grid_spacing=pde.grid_spacing,
            permittivity=pde.permittivity,
            permeability=pde.permeability,
        )

        grid[0, :, :] = fdtd.PeriodicBoundary(name="xbounds")
        grid[:, 0, :] = fdtd.PeriodicBoundary(name="ybounds")
        grid[:, :, 0] = fdtd.PeriodicBoundary(name="zbounds")

        outer_area = (pde.n_large - pde.n) // 2

        # light source planes in the xy plane
        for i in range(6):
            lengthx = rng.randint(2, 6)
            startx = rng.randint(0, outer_area - lengthx)
            lengthy = rng.randint(2, 6)
            starty = rng.randint(0, 16 - lengthy)
            pointz = rng.randint(0, 16)
            ampl = rng.rand() * pde.amplitude
            ps = rng.uniform(low=0.0, high=2 * math.pi)
            p = rng.randint(0, 2)
            polar = ["x", "y"]
            period = pde.wavelength / pde.sol * rng.uniform(low=0.001, high=1e3)
            grid[startx : startx + lengthx, starty : starty + lengthy, pointz] = fdtd.PlaneSource(
                period=period,
                amplitude=ampl,
                name=f"planesourcexy{i}",
                phase_shift=ps,
                polarization=polar[p],
            )

        for i in range(6):
            lengthx = rng.randint(2, 6)
            startx = rng.randint(0, 16 - lengthx)
            pointy = rng.randint(0, 16)
            lengthz = rng.randint(2, 6)
            startz = rng.randint(0, 16 - lengthz)
            ampl = rng.rand() * pde.amplitude
            ps = rng.uniform(low=0.0, high=2 * math.pi)
            p = rng.randint(0, 2)
            polar = ["x", "z"]
            period = pde.wavelength / pde.sol * rng.uniform(low=0.001, high=1e3)
            grid[startx : startx + lengthx, pointy, startz : startz + lengthz] = fdtd.PlaneSource(
                period=period,
                amplitude=ampl,
                name=f"planesourcexz{i}",
                phase_shift=ps,
                polarization=polar[p],
            )

        for i in range(6):
            pointx = rng.randint(0, 16)
            lengthy = rng.randint(2, 6)
            starty = rng.randint(0, 16 - lengthy)
            lengthz = rng.randint(2, 6)
            startz = rng.randint(0, 16 - lengthz)
            ampl = rng.rand() * pde.amplitude
            ps = rng.uniform(low=0.0, high=2 * math.pi)
            p = rng.randint(0, 2)
            polar = ["y", "z"]
            period = pde.wavelength / pde.sol * rng.uniform(low=0.001, high=1e3)
            grid[pointx, starty : starty + lengthy, startz : startz + lengthz] = fdtd.PlaneSource(
                period=period,
                amplitude=ampl,
                name=f"planesourceyz{i}",
                phase_shift=ps,
                polarization=polar[p],
            )

        d_field_, h_field_ = [], []
        grid.run(pde.skip_nt, progress_bar=False)
        for i in range(0, pde.nt):
            grid.run(pde.sample_rate, progress_bar=False)
            e_field_.append(grid.E[outer_area:-outer_area, outer_area:-outer_area, outer_area:-outer_area, :].copy())
            h_field_.append(grid.H[outer_area:-outer_area, outer_area:-outer_area, outer_area:-outer_area, :].copy())

        return np.array(d_field_), np.array(h_field_)


rngs = np.random.randint(np.iinfo(np.int32).max, size=num_samples)
d_field_ls, h_field_ls = zip(
    *Parallel(n_jobs=n_parallel)(delayed(genfunc)(idx, rngs[idx]) for idx in tqdm(range(num_samples)))
)


for idx in range(num_samples):
    # Saving the trajectories
    e_field[idx : (idx + 1), ...] = e_field_ls[idx]
    h_field[idx : (idx + 1), ...] = h_field_ls[idx]