Electromagnetics: Frequency Domain Maxwell’s Equation
This tutorial demonstrates how to use Modulus Sym to do the electromagnetic (EM) simulation. Currently, Modulus Sym offers the following features for frequency domain EM simulation:
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.)
Frequency domain Maxwell’s equation in vector form. (This is available for 3D case and only real form is available now.)
Perfect electronic conductor (PEC) boundary conditions for 2D and 3D cases.
Radiation boundary condition (or, absorbing boundary condition) for 3D.
1D waveguide port solver for 2D waveguide source.
Two electromagnetic problems are solved in this tutorial. All the simulations are appropriately nondimensionalized.
This tutorial assumes that you have completed the tutorial Introductory Example and are familiar with Modulus Sym APIs
All the scripts referred in this tutorial can be found in examples/waveguide/
.
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.
Fig. 120 Domain of 2D waveguide
In this example, the waveguide problem is solved by transverse magnetic (\(TM_z\)) mode, so that the unknown variable is \(E_z(x,y)\). The governing equation in \(\Omega\) is
(189)\[\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:
(190)\[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 Sym to setup the EM solver. Similar to the previous tutorials, you will first import the necessary libraries.
# limitations under the License.
import os
import warnings
from sympy import Symbol, pi, sin, Number, Eq
from sympy.logic.boolalg import Or
import modulus.sym
from modulus.sym.hydra import to_absolute_path, instantiate_arch, ModulusConfig
from modulus.sym.utils.io import csv_to_dict
from modulus.sym.solver import Solver
from modulus.sym.domain import Domain
from modulus.sym.geometry.primitives_2d import Rectangle
from modulus.sym.domain.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
)
from modulus.sym.domain.validator import PointwiseValidator
from modulus.sym.domain.inferencer import PointwiseInferencer
from modulus.sym.utils.io.plotter import ValidatorPlotter, InferencerPlotter
from modulus.sym.key import Key
Then, define the variables for sympy
symbolic calculation and parameters for geometry.
Also, before you define the main classes for Modulus Sym, 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.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# 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)
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.
# define geometry
rec = Rectangle((0, 0), (width, height))
# 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")]
)
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.
waveguide_domain = Domain()
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_domain.add_constraint(PEC, "PEC")
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),
)
waveguide_domain.add_constraint(Waveguide_port, "Waveguide_port")
ABC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"normal_gradient_u": 0.0},
batch_size=cfg.batch_size.ABC,
lambda_weighting={"normal_gradient_u": 10.0},
criteria=Eq(x, width),
)
waveguide_domain.add_constraint(ABC, "ABC")
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, you can import the csv
files for the validation
domain below.
file_path = "../validation/2Dwaveguide_32_2.csv"
if os.path.exists(to_absolute_path(file_path)):
mapping = {"x": "x", "y": "y", "u": "u"}
validation_var = csv_to_dict(to_absolute_path(file_path), 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(
nodes=nodes,
invar=validation_invar_numpy,
true_outvar=validation_outvar_numpy,
batch_size=2048,
plotter=ValidatorPlotter(),
Inferencer has been implemented to plot the results.
waveguide_domain.add_validator(csv_validator)
# add inferencer data
csv_inference = PointwiseInferencer(
nodes=nodes,
invar=validation_invar_numpy,
output_names=["u"],
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 Sym prediction, and their difference
are shown below.
Fig. 121 Modulus Sym, wavenumber=\(32\)
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.
Fig. 122 Domain of 2D Dielectric slab waveguide
In the dielectric, set the relative permittivity \(\epsilon_r=2\). That is,
(191)\[\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.
* (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 = 16.0 # wave_number = freq/c
For the geometry part, you 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.
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# 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)
)
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.
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 = (
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_domain.add_constraint(PEC, "PEC")
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
Fig. 123 Modulus Sym, wavenumber=\(16\)
This example, shows how to setup a 3D waveguide simulation in Modulus Sym. Unlike the previous examples, the features in Modulus Sym to define the boundary condition are used. The geometry is \(\Omega = [0,2]^3\), as shown below.
Problem setup
Fig. 124 3D waveguide geometry
We will solve the 3D frequency domain Maxwell’s equation for electronic field \(\mathbf{E}=(E_x, E_y, E_z)\):
(192)\[\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 Sym 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
(193)\[\mathbf{n}\times\nabla\times \mathbf{E} = 0,\]
while the PEC is
(194)\[\mathbf{n}\times \mathbf{E} = 0.\]
Case setup
This section shows how to use Modulus Sym to setup the 3D frequency EM solver, especially for the boundary conditions.
First import the necessary libraries.
# limitations under the License.
import os
import warnings
from sympy import Symbol, pi, sin, Number, Eq, And
import modulus.sym
from modulus.sym.hydra import instantiate_arch, ModulusConfig
from modulus.sym.solver import Solver
from modulus.sym.domain import Domain
from modulus.sym.geometry.primitives_3d import Box
from modulus.sym.domain.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
)
from modulus.sym.domain.inferencer import PointwiseInferencer
from modulus.sym.utils.io.plotter import InferencerPlotter
Define sympy
variables, geometry and waveguide
function.
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# 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 the PDE class and neural network structure.
# define geometry
rec = Box((0, 0, 0), (width, length, height))
# 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()
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_domain.add_constraint(wall_PEC, "PEC")
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),
)
waveguide_domain.add_constraint(Waveguide_port, "Waveguide_port")
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),
)
waveguide_domain.add_constraint(ABC, "ABC")
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.
waveguide_domain.add_constraint(Interior, "Interior")
# add inferencer data
interior_points = rec.sample_interior(
10000, bounds={x: (0, width), y: (0, length), z: (0, height)}
)
numpy_inference = PointwiseInferencer(
nodes=nodes,
invar=interior_points,
output_names=["ux", "uy", "uz"],
plotter=InferencerPlotter(),
batch_size=2048,
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
slices of the three components are shown in below.
Fig. 125 3D waveguide, \(E_x\)
Fig. 126 3D waveguide, \(E_y\)
Fig. 127 3D waveguide, \(E_z\)
This example, shows a 3D dielectric slab waveguide. In this case, consider a unit cube \([0,1]^3\) with a dielectric slab centered in the middle along \(y\) axis. The length of the slab is \(0.2\). Fig. 128 and Fig. 129 shows the whole geometry and an \(xz\) cross-section that shows the dielectric slab.
Fig. 128 Geometry for 3D dielectric slab waveguide
Fig. 129 \(xz\) cross-sectional view
The permittivity is defined as follows
(195)\[\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, only the differences compared to the
last example are covered here. 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 piecewise sympy
functions, use
Heaviside
instead of Piecewise
as the later cannot be complied
in Modulus Sym for the time being. The waveguide data can also be imported using the csv_to_dict()
function.
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# 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 = 16.0 # wave_number = freq/c
file_path = "../validation/2Dwaveguideport.csv"
if not os.path.exists(to_absolute_path(file_path)):
warnings.warn(
f"Directory{file_path}does not exist. Cannot continue. Please download the additional files from NGC https://catalog.ngc.nvidia.com/orgs/nvidia/teams/modulus/resources/modulus_sym_examples_supplemental_materials"
)
sys.exit()
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"]) - 0.5,
"y": data_var["x"],
"z": data_var["y"],
}
waveguide_port_outvar_numpy = {"uz": data_var["u0"]}
# define geometry
rec = Box(
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.
(-width / 2, -length / 2, -height / 2), (width / 2, length / 2, height / 2)
)
# 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()
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_domain.add_constraint(wall_PEC, "PEC")
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)},
)
waveguide_domain.add_constraint(Waveguide_port, "Waveguide_port")
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,
)
waveguide_domain.add_constraint(ABC, "ABC")
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,
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.
waveguide_domain.add_constraint(Interior, "Interior")
# add inferencer data
slab_inference = VoxelInferencer(
bounds=[[-0.5, 0.5], [-0.5, 0.5], [-0.5, 0.5]],
npoints=[128, 128, 128],
nodes=nodes,
output_names=["ux", "uy", "uz"],
Results
The full code of this example can be found in
examples/waveguide/slab_3D/slab_3D.py
. Here the
simulation for different wavenumbers is done. Below figures show the result for wavenumber equals \(16\).
Fig. 130 3D dielectric slab, \(E_x\)
Fig. 131 3D dielectric slab, \(E_y\)
Fig. 132 3D dielectric slab, \(E_z\)
Also the results of higher wavenumber \(32\) are shown below.
Fig. 133 3D dielectric slab, \(E_x\)
Fig. 134 3D dielectric slab, \(E_y\)
Fig. 135 3D dielectric slab, \(E_z\)