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.

Note

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/.

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.

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.

Copy
Copied!

            
            # 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.

Copy
Copied!

            
            @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.

Copy
Copied!

            
            # 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.

Copy
Copied!

            
            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.

Copy
Copied!

            
            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.

Copy
Copied!

            
            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\)

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.

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.

Copy
Copied!

            
                    * (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.

Copy
Copied!

            
            @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.

Copy
Copied!

            
            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:

Copy
Copied!

            
            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\)

Problem 3: 3D waveguide cavity

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.

Copy
Copied!

            
            # 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.

Copy
Copied!

            
            @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.

Copy
Copied!

            
            # 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.

Copy
Copied!

            
            )

    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.

Copy
Copied!

            
            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\)

Problem 4: 3D Dielectric slab waveguide

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.

Copy
Copied!

            
            @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.

Copy
Copied!

            
            (-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.

Copy
Copied!

            
            )

    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.

Copy
Copied!

            
            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\).