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:

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

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.

2Dwaveguide.png

Fig. 118 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!
            

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 from modulus.sym.eq.pdes.wave_equation import HelmholtzEquation from modulus.sym.eq.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 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!
            

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

Copy
Copied!
            

# 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")] ) 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.

Copy
Copied!
            

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, }, ) waveguide_domain.add_constraint(Interior, "Interior")

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

Copy
Copied!
            

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(), ) waveguide_domain.add_validator(csv_validator)

Inferencer has been implemented to plot the results.

Copy
Copied!
            

# add inferencer data csv_inference = PointwiseInferencer( nodes=nodes, invar=validation_invar_numpy, output_names=["u"], plotter=InferencerPlotter(), batch_size=2048, ) waveguide_domain.add_inferencer(csv_inference, "inf_data")

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.

2Dwaveguide_modulus.png

Fig. 119 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.

2Dslab_geo.png

Fig. 120 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!
            

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 waveguide_port_invar_numpy = {"x": np.zeros_like(yv), "y": yv} waveguide_port_outvar_numpy = {"u": 10 * eigvecs[:, 0:1]}

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!
            

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

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!
            

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")] )

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_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)}, ) waveguide_domain.add_constraint(Waveguide_port, "Waveguide_port")

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

2Dslab_16.png

Fig. 121 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

3Dwaveguide_geo.png

Fig. 122 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!
            

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 from modulus.sym.key import Key from modulus.sym.eq.pdes.electromagnetic import PEC, SommerfeldBC, MaxwellFreqReal

Define sympy variables, geometry and waveguide function.

Copy
Copied!
            

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

Copy
Copied!
            

# 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")] )

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, ) waveguide_domain.add_constraint(Interior, "Interior")

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!
            

# 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, ) 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 slices of the three components are shown in below.

3Dwaveguide_Ex.png

Fig. 123 3D waveguide, \(E_x\)

3Dwaveguide_Ey.png

Fig. 124 3D waveguide, \(E_y\)

3Dwaveguide_Ez.png

Fig. 125 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. 126 and Fig. 127 shows the whole geometry and an \(xz\) cross-section that shows the dielectric slab.

3Dslab_geo.png

Fig. 126 Geometry for 3D dielectric slab waveguide

3Dslab_geo_xz.png

Fig. 127 \(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!
            

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

Copy
Copied!
            

# 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")] )

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, ) waveguide_domain.add_constraint(Interior, "Interior")

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

Copy
Copied!
            

# 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"], ) waveguide_domain.add_inferencer(slab_inference, "Inf" + str(int(wave_number)))

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

3Dslab_16_Ex.png

Fig. 128 3D dielectric slab, \(E_x\)

3Dslab_16_Ey.png

Fig. 129 3D dielectric slab, \(E_y\)

3Dslab_16_Ez.png

Fig. 130 3D dielectric slab, \(E_z\)

Also the results of higher wavenumber \(32\) are shown below.

3Dslab_32_Ex.png

Fig. 131 3D dielectric slab, \(E_x\)

3Dslab_32_Ey.png

Fig. 132 3D dielectric slab, \(E_y\)

3Dslab_32_Ez.png

Fig. 133 3D dielectric slab, \(E_z\)

Previous Moving Time Window: Taylor Green Vortex Decay
Next Fully Developed Turbulent Channel Flow
© Copyright 2023, NVIDIA Modulus Team. Last updated on Jan 25, 2024.