# Electromagnetics: Frequency Domain Maxwell’s Equation¶

## Introduction¶

This tutorial demonstrates how to use Modulus to do the electromagnetic (EM) simulation. Currently, Modulus offers the following features for frequency domain EM simulation:

1. Frequency domain Maxwell’s equation in scalar form. This is same to Helmholtz equation. (This is available for 1D, 2D, and 3D. Only real form is available now.)

2. Frequency domain Maxwell’s equation in vector form. (This is available for 3D case and only real form is available now.)

3. Perfect electronic conductor (PEC) boundary conditions for 2D and 3D cases.

4. Radiation boundary condition (or, absorbing boundary condition) for 3D.

5. 1D waveguide port solver for 2D waveguide source.

Two electromagnetic problems are solved in this tutorial. All the simulations are appropriately non-dimensionalized.

Note

This tutorial assumes that you have completed the tutorial Lid Driven Cavity Background and are familiar with Modulus APIs

All the scripts referred in this tutorial can be found in examples/waveguide/.

## Problem 1: 2D Waveguide Cavity¶

Consider a 2D domain $$\Omega=[0, 2]\times [0, 2]$$ as shown below. The whole domain is vacuum. Say, relative permittivity $$\epsilon_r = 1$$. The left boundary is a waveguide port while the right boundary is absorbing boundary (or ABC). The top and the bottom is PEC.

In this example, the waveguide problem is solved by transverse-magnetic ($$TM_z$$) mode, so that our unknown variable is $$E_z(x,y)$$. The governing equation in $$\Omega$$ is

(179)$\Delta E_z(x,y) + k^2E_z(x,y) = 0,$

where $$k$$ is the wavenumber. Notice in 2D scalar case, the PEC and ABC will be simplified in the following form, respectively:

(180)$E_z(x,y)=0\mbox{ on top and bottom boundaries, }\frac{\partial E_z}{\partial y}=0\mbox{ on right boundary.}$

### Case Setup¶

This subsection shows how to use Modulus to setup the EM solver. Similar to our previous tutorials, we first import the necessary libraries.

from sympy import Symbol, pi, sin, Number, Eq
from sympy.logic.boolalg import Or

import modulus
from modulus.hydra import to_absolute_path, to_yaml, instantiate_arch
from modulus.hydra.config import ModulusConfig
from modulus.csv_utils.csv_rw import csv_to_dict
from modulus.continuous.solvers.solver import Solver
from modulus.continuous.domain.domain import Domain
from modulus.geometry.csg.csg_2d import Rectangle
from modulus.continuous.constraints.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
)
from modulus.continuous.validator.validator import PointwiseValidator
from modulus.continuous.inferencer.inferencer import PointwiseInferencer
from modulus.tensorboard_utils.plotter import ValidatorPlotter, InferencerPlotter
from modulus.key import Key
from modulus.PDES.wave_equation import HelmholtzEquation
from modulus.PDES.navier_stokes import GradNormal


Then, define the variables for sympy symbolic calculation and parameters for geometry. Also, before you define the main classes for Modulus, you need to compute the eigenmode for waveguide solver. Since the material is uniform (vacuum), the closed form of the eigenmode is of the form $$\sin(\frac{k\pi y}{L})$$, where $$L$$ is the length of the port, and $$k = 1, 2,\cdots$$. Then define the waveguide port directly by using sympy function. The code for the geometry and computing eigenmode can be found below.

@modulus.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
print(to_yaml(cfg))

# params for domain
height = 2
width = 2

eigenmode = [2]
wave_number = 32.0  # wave_number = freq/c
waveguide_port = Number(0)
for k in eigenmode:
waveguide_port += sin(k * pi * y / height)

# define geometry
rec = Rectangle((0, 0), (width, height))


For wave simulation, since the result is always periodic, Fourier feature will be greatly helpful for the convergence and accuracy. The frequency of the Fourier feature can be implied by the wavenumber. This block of code shows the solver setup. Also, define the normal gradient for the boundary conditions. Finally, make the domain for training.

    # make list of nodes to unroll graph on
hm = HelmholtzEquation(u="u", k=wave_number, dim=2)
gn = GradNormal(T="u", dim=2, time=False)
wave_net = instantiate_arch(
input_keys=[Key("x"), Key("y")],
output_keys=[Key("u")],
frequencies=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number) * 2 + 1)],
),
frequencies_params=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number) * 2 + 1)],
),
cfg=cfg.arch.modified_fourier,
)
nodes = (
hm.make_nodes()
+ gn.make_nodes()
+ [wave_net.make_node(name="wave_network", jit=cfg.jit)]
)

waveguide_domain = Domain()


Now, define the constraints for PDEs and boundary conditions. The BCs are defined based on the explanations provided above. In the interior domain, the weights of the PDE is 1.0/wave_number**2. This is because when the wavenumber is large, the PDE loss will be very large in the beginning and will potentially break the training. Using this weighting method you can eliminate this phenomenon.

    PEC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"u": 0.0},
batch_size=cfg.batch_size.PEC,
lambda_weighting={"u": 100.0},
criteria=Or(Eq(y, 0), Eq(y, height)),
)

Waveguide_port = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"u": waveguide_port},
batch_size=cfg.batch_size.Waveguide_port,
lambda_weighting={"u": 100.0},
criteria=Eq(x, 0),
)

ABC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
batch_size=cfg.batch_size.ABC,
criteria=Eq(x, width),
)

Interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=rec,
outvar={"helmholtz": 0.0},
batch_size=cfg.batch_size.Interior,
bounds={x: (0, width), y: (0, height)},
lambda_weighting={
"helmholtz": 1.0 / wave_number ** 2,
},
)


To validate the result, we import the csv files for the validation domain below.

    # add validation data
mapping = {"x": "x", "y": "y", "u": "u"}
validation_var = csv_to_dict(
to_absolute_path("../validation/2Dwaveguide_32_2.csv"), mapping
)
validation_invar_numpy = {
key: value for key, value in validation_var.items() if key in ["x", "y"]
}
validation_outvar_numpy = {
key: value for key, value in validation_var.items() if key in ["u"]
}

csv_validator = PointwiseValidator(
validation_invar_numpy,
validation_outvar_numpy,
nodes,
batch_size=2048,
plotter=ValidatorPlotter(),
)


Inferencer has been implemented to plot the results.

    # add inferencer data
csv_inference = PointwiseInferencer(
validation_invar_numpy,
["u"],
nodes,
plotter=InferencerPlotter(),
batch_size=2048,
)


### Results¶

The full code of this example can be found in examples/waveguide/cavity_2D/waveguide2D_TMz.py. The simulation with wavenumber equals $$32$$. The solution from comercial solver, Modulus prediction, and their difference are shown below.

## Problem 2: 2D Dielectric slab waveguide¶

This section, demonstrates using the 2D waveguide simulation with a dielectric slab. The problem setup is almost same as before except there is a horizontal dielectric slab in the middle of the domain. The domain is shown below.

In the dielectric, set the relative permittivity $$\epsilon_r=2$$. That is,

(181)$\begin{split}\epsilon_r = \begin{cases} 2 & \mbox{ in dielectric slab,}\\ 1 & \mbox{ otherwise.} \end{cases}\end{split}$

All the other settings are kept the same to the previous example.

### Case setup¶

Here, for simplicity, only the parts of the code that are different than the previous example are shown. The main difference is the spatially dependent permittivity. First, compute the eigenfunctions on the left boundary.

# helper function for computing laplacian eigen values
def Laplacian_1D_eig(a, b, N, eps=lambda x: np.ones_like(x), k=3):
n = N - 2
h = (b - a) / (N - 1)

L = diags([1, -2, 1], [-1, 0, 1], shape=(n, n))
L = -L / h ** 2

x = np.linspace(a, b, num=N)
M = diags([eps(x[1:-1])], [0])

eigvals, eigvecs = eigsh(L, k=k, M=M, which="SM")
eigvecs = np.vstack((np.zeros((1, k)), eigvecs, np.zeros((1, k))))
norm_eigvecs = np.linalg.norm(eigvecs, axis=0)
eigvecs /= norm_eigvecs
return eigvals.astype(np.float32), eigvecs.astype(np.float32), x.astype(np.float32)


For the geometry part, we will need to define the slab and corresponding permittivity function. There is a square root in eps_sympy because in the HelmholtzEquation, the wavenumber will be squared. Next, based on the permittivity function, use the eigensolver to get the numerical waveguide port.

    # params for domain
height = 2
width = 2

len_slab = 0.6
eps0 = 1.0
eps1 = 2.0
eps_numpy = lambda y: np.where(
np.logical_and(y > (height - len_slab) / 2, y < (height + len_slab) / 2),
eps1,
eps0,
)
eps_sympy = sqrt(
eps0
+ (
Heaviside(y - (height - len_slab) / 2)
- Heaviside(y - (height + len_slab) / 2)
)
* (eps1 - eps0)
)
eigvals, eigvecs, yv = Laplacian_1D_eig(0, height, 1000, eps=eps_numpy, k=3)
yv = yv.reshape((-1, 1))
eigenmode = [1]
wave_number = 32.0  # wave_number = freq/c
waveguide_port_invar_numpy = {"x": np.zeros_like(yv), "y": yv}
waveguide_port_outvar_numpy = {"u": 10 * eigvecs[:, 0:1]}

# define geometry
rec = Rectangle((0, 0), (width, height))


In the definition of the PDEs and neural network’s structure, set the k in HelmholtzEquation as the product of wavenumber and permittivity function. Also, update the frequency for the Fourier features to suit the problem.

    # make list of nodes to unroll graph on
hm = HelmholtzEquation(u="u", k=wave_number * eps_sympy, dim=2)
gn = GradNormal(T="u", dim=2, time=False)
wave_net = instantiate_arch(
input_keys=[Key("x"), Key("y")],
output_keys=[Key("u")],
frequencies=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number * np.sqrt(eps1)) * 2 + 1)],
),
frequencies_params=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number * np.sqrt(eps1)) * 2 + 1)],
),
cfg=cfg.arch.modified_fourier,
)
nodes = (
hm.make_nodes()
+ gn.make_nodes()
+ [wave_net.make_node(name="wave_network", jit=cfg.jit)]
)


Now, define the constraints. The only difference here is the left boundary, which will be given by a numpy array. Only the modified BC is shown below:

    Waveguide_port = PointwiseConstraint.from_numpy(
nodes=nodes,
invar=waveguide_port_invar_numpy,
outvar=waveguide_port_outvar_numpy,
batch_size=cfg.batch_size.Waveguide_port,
lambda_weighting={"u": np.full_like(yv, 0.5)},
)


### Results¶

The full code of this example can be found in examples/waveguide/slab_2D/slab_2D.py. We do the simulation with wavenumber equals $$16$$ and $$32$$, respectively. The results are shown in figure below

## Problem 3: 3D waveguide cavity¶

This example, shows how to setup a 3D waveguide simulation in Modulus. Unlike the previous examples, the features in Modulus to define the boundary condition are used. The geometry is $$\Omega = [0,2]^3$$, as shown below.

### Problem setup¶

We will solve the 3D frequency domain Maxwell’s equation for electronic field $$\mathbf{E}=(E_x, E_y, E_z)$$:

(182)$\nabla\times \nabla\times \mathbf{E}+\epsilon_rk^2\mathbf{E} = 0,$

where $$\epsilon_r$$ is the permittivity, and the $$k$$ is the wavenumber. Note that, currently Modulus only support real permittivity and wavenumber. For the sake of simplicity, assume the permeability $$\mu_r=1$$. As before, waveguide port has been applied on the left. We apply absorbing boundary condition on the right side and PEC for the rest. In 3D, the absorbing boundary condition for real form reads

(183)$\mathbf{n}\times\nabla\times \mathbf{E} = 0,$

while the PEC is

(184)$\mathbf{n}\times \mathbf{E} = 0.$

### Case setup¶

This section shows how to use Modulus to setup the 3D frequency EM solver, especially for the boundary conditions.

First import the necessary libraries.

from sympy import Symbol, pi, sin, Number, Eq, And

import modulus
from modulus.hydra import to_yaml, instantiate_arch
from modulus.hydra.config import ModulusConfig
from modulus.continuous.solvers.solver import Solver
from modulus.continuous.domain.domain import Domain
from modulus.geometry.csg.csg_3d import Box
from modulus.continuous.constraints.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
)
from modulus.continuous.inferencer.inferencer import PointwiseInferencer
from modulus.tensorboard_utils.plotter import InferencerPlotter
from modulus.key import Key
from modulus.PDES.electromagnetic import PEC, SommerfeldBC, MaxwellFreqReal


Define sympy variables, geometry and waveguide function.

@modulus.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
print(to_yaml(cfg))

# params for domain
length = 2
height = 2
width = 2

eigenmode = [1]
wave_number = 16.0  # wave_number = freq/c
waveguide_port = Number(0)
for k in eigenmode:
waveguide_port += sin(k * pi * y / length) * sin(k * pi * z / height)

# define geometry
rec = Box((0, 0, 0), (width, length, height))


Define the PDE class and neural network structure.

    # make list of nodes to unroll graph on
hm = MaxwellFreqReal(k=wave_number)
pec = PEC()
pml = SommerfeldBC()
wave_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("z")],
output_keys=[Key("ux"), Key("uy"), Key("uz")],
frequencies=("axis,diagonal", [i / 2.0 for i in range(int(wave_number) + 1)]),
frequencies_params=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number) + 1)],
),
cfg=cfg.arch.modified_fourier,
)
nodes = (
hm.make_nodes()
+ pec.make_nodes()
+ pml.make_nodes()
+ [wave_net.make_node(name="wave_network", jit=cfg.jit)]
)


Then, define the constraints for PDEs and boundary conditions, and all boundary conditions. The 3D Maxwell’s equations has been implemented in Maxwell_Freq_real_3D, PEC has been implemented in PEC_3D, and absorbing boundary condition has been implemented in SommerfeldBC_real_3D. We may use these features directly to apply the corresponding constraints.

    waveguide_domain = Domain()

wall_PEC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"PEC_x": 0.0, "PEC_y": 0.0, "PEC_z": 0.0},
batch_size=cfg.batch_size.PEC,
lambda_weighting={"PEC_x": 100.0, "PEC_y": 100.0, "PEC_z": 100.0},
criteria=And(~Eq(x, 0), ~Eq(x, width)),
)

Waveguide_port = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"uz": waveguide_port},
batch_size=cfg.batch_size.Waveguide_port,
lambda_weighting={"uz": 100.0},
criteria=Eq(x, 0),
)

ABC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={
"SommerfeldBC_real_x": 0.0,
"SommerfeldBC_real_y": 0.0,
"SommerfeldBC_real_z": 0.0,
},
batch_size=cfg.batch_size.ABC,
lambda_weighting={
"SommerfeldBC_real_x": 10.0,
"SommerfeldBC_real_y": 10.0,
"SommerfeldBC_real_z": 10.0,
},
criteria=Eq(x, width),
)

Interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=rec,
outvar={
"Maxwell_Freq_real_x": 0,
"Maxwell_Freq_real_y": 0.0,
"Maxwell_Freq_real_z": 0.0,
},
batch_size=cfg.batch_size.Interior,
bounds={x: (0, width), y: (0, length), z: (0, height)},
lambda_weighting={
"Maxwell_Freq_real_x": 1.0 / wave_number ** 2,
"Maxwell_Freq_real_y": 1.0 / wave_number ** 2,
"Maxwell_Freq_real_z": 1.0 / wave_number ** 2,
},
fixed_dataset=False,
)


Note that this is done in 3D, so the PDEs, PEC and absorbing boundaries have three output components.

Inferencer domain has been defined to check the result.

    # add inferencer data
interior_points = rec.sample_interior(
10000, bounds={x: (0, width), y: (0, length), z: (0, height)}
)
numpy_inference = PointwiseInferencer(
interior_points,
["ux", "uy", "uz"],
nodes,
plotter=InferencerPlotter(),
batch_size=2048,
)
waveguide_domain.add_inferencer(numpy_inference, "Inf" + str(wave_number).zfill(4))


### Results¶

The full code of this example can be found in examples/waveguide/waveguide3D.py. The wavenumber equals $$32$$ and use second eigenmode for $$y$$ and $$z$$. The silces of the three components are shown in below.

## Problem 4: 3D Dielectric slab waveguide¶

This example, shows a 3D dielectric slab waveguide. In this case, consider a unit cude $$[0,1]^3$$ with a dielectric slab centered in the middle along $$y$$ axis. The length of the slab is $$0.2$$. Fig. 102 and Fig. 103 shows the whole geometry and an $$xz$$ cross-section that shows the dielectric slab.

The permittivity is defined as follows

(185)$\begin{split}\epsilon_r = \begin{cases} 1.5 & \mbox{ in dielectric slab,}\\ 1 & \mbox{ otherwise.} \end{cases}\end{split}$

### Case setup¶

For the sake of simplicity, we only show the differences compared to the last example. The main difference for this simulation is that you need to import the waveguide result from a csv file and then use that as the waveguide port boundary condition.

First define the geometry and the sympy permittivity function. To define the piece-wise sympy functions, use Heaviside instead of Piecewise as the later cannot be complied in Modulus for the time being. The waveguide data can also be imported using the csv_to_dict() function.

@modulus.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
print(to_yaml(cfg))

# params for domain
length = 1
height = 1
width = 1

len_slab = 0.2
eps0 = 1.0
eps1 = 1.5
eps_sympy = sqrt(
eps0
+ (Heaviside(y + len_slab / 2) - Heaviside(y - len_slab / 2))
* (Heaviside(z + len_slab / 2) - Heaviside(z - len_slab / 2))
* (eps1 - eps0)
)

wave_number = 32.0  # wave_number = freq/c
mapping = {"x": "x", "y": "y", **{"u" + str(k): "u" + str(k) for k in range(6)}}
data_var = csv_to_dict(
to_absolute_path("../validation/2Dwaveguideport.csv"), mapping
)
waveguide_port_invar_numpy = {
"x": np.zeros_like(data_var["x"]),
"y": data_var["x"],
"z": data_var["y"],
}
waveguide_port_outvar_numpy = {"uz": data_var["u0"]}

# define geometry
rec = Box(
(-width / 2, -length / 2, -height / 2), (width / 2, length / 2, height / 2)
)


In validation/2Dwaveguideport.csv, there are six eigenmodes. You may try different modes to explore more interesting results.

Next, define the PDEs classes and neural network structure.

    # make list of nodes to unroll graph on
hm = MaxwellFreqReal(k=wave_number * eps_sympy)
pec = PEC()
pml = SommerfeldBC()
wave_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("z")],
output_keys=[Key("ux"), Key("uy"), Key("uz")],
frequencies=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number * np.sqrt(eps1)) * 2 + 1)],
),
frequencies_params=(
"axis,diagonal",
[i / 2.0 for i in range(int(wave_number * np.sqrt(eps1)) * 2 + 1)],
),
cfg=cfg.arch.modified_fourier,
)
nodes = (
hm.make_nodes()
+ pec.make_nodes()
+ pml.make_nodes()
+ [wave_net.make_node(name="wave_network", jit=cfg.jit)]
)


Then, define the constraints. Here imported data is used as the waveguide port boundary condition.

    waveguide_domain = Domain()

wall_PEC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"PEC_x": 0.0, "PEC_y": 0.0, "PEC_z": 0.0},
batch_size=cfg.batch_size.PEC,
lambda_weighting={"PEC_x": 100.0, "PEC_y": 100.0, "PEC_z": 100.0},
criteria=And(~Eq(x, -width / 2), ~Eq(x, width / 2)),
fixed_dataset=False,
)

Waveguide_port = PointwiseConstraint.from_numpy(
nodes=nodes,
invar=waveguide_port_invar_numpy,
outvar=waveguide_port_outvar_numpy,
batch_size=cfg.batch_size.Waveguide_port,
lambda_weighting={"uz": np.full_like(waveguide_port_invar_numpy["x"], 0.5)},
)

ABC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={
"SommerfeldBC_real_x": 0.0,
"SommerfeldBC_real_y": 0.0,
"SommerfeldBC_real_z": 0.0,
},
batch_size=cfg.batch_size.ABC,
lambda_weighting={
"SommerfeldBC_real_x": 10.0,
"SommerfeldBC_real_y": 10.0,
"SommerfeldBC_real_z": 10.0,
},
criteria=Eq(x, width / 2),
fixed_dataset=False,
)

Interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=rec,
outvar={
"Maxwell_Freq_real_x": 0,
"Maxwell_Freq_real_y": 0.0,
"Maxwell_Freq_real_z": 0.0,
},
batch_size=cfg.batch_size.Interior,
bounds={
x: (-width / 2, width / 2),
y: (-length / 2, length / 2),
z: (-height / 2, height / 2),
},
lambda_weighting={
"Maxwell_Freq_real_x": 1.0 / wave_number ** 2,
"Maxwell_Freq_real_y": 1.0 / wave_number ** 2,
"Maxwell_Freq_real_z": 1.0 / wave_number ** 2,
},
fixed_dataset=False,
)


Finally, define the inferencer. These are same as the previous example except the bounds for the domain.

    # add inferencer data
interior_points = rec.sample_interior(
10000,
bounds={
x: (-width / 2, width / 2),
y: (-length / 2, length / 2),
z: (-height / 2, height / 2),
},
)
numpy_inference = PointwiseInferencer(
interior_points,
["ux", "uy", "uz"],
nodes,
batch_size=1024,
plotter=InferencerPlotter(),
)
waveguide_domain.add_inferencer(numpy_inference, "Inf" + str(wave_number).zfill(4))


### Results¶

The full code of this example can be found in examples/waveguide/slab_3D/slab_3D.py. Here the simulation for different wavennumbers is done. Below figures show the result for wavenumber equals $$16$$.

Also the results of higher wavenumber $$32$$ are shown below.