Skip to content

Initialization of Boundaries

Currently, the simulation boundary can either be a Perfectly Matched Layer for absorbing all incoming light, or a periodic boundary which wraps around to the other side of the simulation.

fdtdx.objects.boundaries.BoundaryConfig

Bases: ExtendedTreeClass

Configuration class for boundary conditions.

This class stores parameters for boundary conditions in all six directions (min/max x/y/z). Supports both PML and periodic boundaries. For PML, the parameters control the absorption properties and physical size of the PML regions.

Attributes:

Name Type Description
boundary_type_minx str

Boundary type at minimum x ("pml" or "periodic"). Default "pml".
boundary_type_maxx str

Boundary type at maximum x ("pml" or "periodic"). Default "pml".
boundary_type_miny str

Boundary type at minimum y ("pml" or "periodic"). Default "pml".
boundary_type_maxy str

Boundary type at maximum y ("pml" or "periodic"). Default "pml".
boundary_type_minz str

Boundary type at minimum z ("pml" or "periodic"). Default "pml".
boundary_type_maxz str

Boundary type at maximum z ("pml" or "periodic"). Default "pml".
thickness_grid_minx int

Number of grid cells for PML at minimum x boundary. Default 10.
thickness_grid_maxx int

Number of grid cells for PML at maximum x boundary. Default 10.
thickness_grid_miny int

Number of grid cells for PML at minimum y boundary. Default 10.
thickness_grid_maxy int

Number of grid cells for PML at maximum y boundary. Default 10.
thickness_grid_minz int

Number of grid cells for PML at minimum z boundary. Default 10.
thickness_grid_maxz int

Number of grid cells for PML at maximum z boundary. Default 10.
kappa_start_minx float

Initial kappa value at min x boundary. Default 1.0.
kappa_end_minx float

Final kappa value at min x boundary. Default 1.5.
kappa_start_maxx float

Initial kappa value at max x boundary. Default 1.0.
kappa_end_maxx float

Final kappa value at max x boundary. Default 1.5.
kappa_start_miny float

Initial kappa value at min y boundary. Default 1.0.
kappa_end_miny float

Final kappa value at min y boundary. Default 1.5.
kappa_start_maxy float

Initial kappa value at max y boundary. Default 1.0.
kappa_end_maxy float

Final kappa value at max y boundary. Default 1.5.
kappa_start_minz float

Initial kappa value at min z boundary. Default 1.0.
kappa_end_minz float

Final kappa value at min z boundary. Default 1.5.
kappa_start_maxz float

Initial kappa value at max z boundary. Default 1.0.
kappa_end_maxz float

Final kappa value at max z boundary. Default 1.5.
Source code in src/fdtdx/objects/boundaries/initialization.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
@extended_autoinit
class BoundaryConfig(ExtendedTreeClass):
    """Configuration class for boundary conditions.

    This class stores parameters for boundary conditions in all six directions (min/max x/y/z).
    Supports both PML and periodic boundaries. For PML, the parameters control the absorption
    properties and physical size of the PML regions.

    Attributes:
        boundary_type_minx (str): Boundary type at minimum x ("pml" or "periodic"). Default "pml".
        boundary_type_maxx (str): Boundary type at maximum x ("pml" or "periodic"). Default "pml".
        boundary_type_miny (str): Boundary type at minimum y ("pml" or "periodic"). Default "pml".
        boundary_type_maxy (str): Boundary type at maximum y ("pml" or "periodic"). Default "pml".
        boundary_type_minz (str): Boundary type at minimum z ("pml" or "periodic"). Default "pml".
        boundary_type_maxz (str): Boundary type at maximum z ("pml" or "periodic"). Default "pml".
        thickness_grid_minx (int): Number of grid cells for PML at minimum x boundary. Default 10.
        thickness_grid_maxx (int): Number of grid cells for PML at maximum x boundary. Default 10.
        thickness_grid_miny (int): Number of grid cells for PML at minimum y boundary. Default 10.
        thickness_grid_maxy (int): Number of grid cells for PML at maximum y boundary. Default 10.
        thickness_grid_minz (int): Number of grid cells for PML at minimum z boundary. Default 10.
        thickness_grid_maxz (int): Number of grid cells for PML at maximum z boundary. Default 10.
        kappa_start_minx (float): Initial kappa value at min x boundary. Default 1.0.
        kappa_end_minx (float): Final kappa value at min x boundary. Default 1.5.
        kappa_start_maxx (float): Initial kappa value at max x boundary. Default 1.0.
        kappa_end_maxx (float): Final kappa value at max x boundary. Default 1.5.
        kappa_start_miny (float): Initial kappa value at min y boundary. Default 1.0.
        kappa_end_miny (float): Final kappa value at min y boundary. Default 1.5.
        kappa_start_maxy (float): Initial kappa value at max y boundary. Default 1.0.
        kappa_end_maxy (float): Final kappa value at max y boundary. Default 1.5.
        kappa_start_minz (float): Initial kappa value at min z boundary. Default 1.0.
        kappa_end_minz (float): Final kappa value at min z boundary. Default 1.5.
        kappa_start_maxz (float): Initial kappa value at max z boundary. Default 1.0.
        kappa_end_maxz (float): Final kappa value at max z boundary. Default 1.5.
    """

    boundary_type_minx: str = "pml"
    boundary_type_maxx: str = "pml"
    boundary_type_miny: str = "pml"
    boundary_type_maxy: str = "pml"
    boundary_type_minz: str = "pml"
    boundary_type_maxz: str = "pml"
    thickness_grid_minx: int = 10
    thickness_grid_maxx: int = 10
    thickness_grid_miny: int = 10
    thickness_grid_maxy: int = 10
    thickness_grid_minz: int = 10
    thickness_grid_maxz: int = 10
    kappa_start_minx: float = 1.0
    kappa_end_minx: float = 1.5
    kappa_start_maxx: float = 1.0
    kappa_end_maxx: float = 1.5
    kappa_start_miny: float = 1.0
    kappa_end_miny: float = 1.5
    kappa_start_maxy: float = 1.0
    kappa_end_maxy: float = 1.5
    kappa_start_minz: float = 1.0
    kappa_end_minz: float = 1.5
    kappa_start_maxz: float = 1.0
    kappa_end_maxz: float = 1.5

    def get_dict(self) -> dict[str, int]:
        """Gets a dictionary mapping boundary names to their grid thicknesses.

        Returns:
            dict[str, int]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
                mapping to their respective grid thickness values.
        """
        return {
            "min_x": self.thickness_grid_minx,
            "max_x": self.thickness_grid_maxx,
            "min_y": self.thickness_grid_miny,
            "max_y": self.thickness_grid_maxy,
            "min_z": self.thickness_grid_minz,
            "max_z": self.thickness_grid_maxz,
        }

    def get_type_dict(self) -> dict[str, str]:
        """Gets a dictionary mapping boundary names to their boundary types.

        Returns:
            dict[str, str]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
                mapping to their respective boundary types ("pml" or "periodic").
        """
        return {
            "min_x": self.boundary_type_minx,
            "max_x": self.boundary_type_maxx,
            "min_y": self.boundary_type_miny,
            "max_y": self.boundary_type_maxy,
            "min_z": self.boundary_type_minz,
            "max_z": self.boundary_type_maxz,
        }

    def get_kappa_dict(
        self,
        prop: Literal["kappa_start", "kappa_end"],
    ) -> dict[str, float]:
        """Gets a dictionary mapping boundary names to their kappa values.

        Args:
            prop: Which kappa property to get, either "kappa_start" or "kappa_end"

        Returns:
            dict[str, float]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
                mapping to their respective kappa values.

        Raises:
            Exception: If prop is not "kappa_start" or "kappa_end"
        """
        if prop == "kappa_start":
            return {
                "min_x": self.kappa_start_minx,
                "max_x": self.kappa_start_maxx,
                "min_y": self.kappa_start_miny,
                "max_y": self.kappa_start_maxy,
                "min_z": self.kappa_start_minz,
                "max_z": self.kappa_start_maxz,
            }
        elif prop == "kappa_end":
            return {
                "min_x": self.kappa_end_minx,
                "max_x": self.kappa_end_maxx,
                "min_y": self.kappa_end_miny,
                "max_y": self.kappa_end_maxy,
                "min_z": self.kappa_end_minz,
                "max_z": self.kappa_end_maxz,
            }
        else:
            raise Exception(f"Unknown: {prop=}")

    def get_inside_boundary_slice(self) -> tuple[slice, slice, slice]:
        """Gets slice objects for the non-PML interior region of the simulation volume.

        Returns:
            tuple[slice, slice, slice]: Three slice objects for indexing the x, y, z dimensions
                respectively, excluding the PML boundary regions.
        """
        return (
            slice(
                self.thickness_grid_minx + 1 if self.boundary_type_minx == "pml" else 0,
                -self.thickness_grid_maxx - 1 if self.boundary_type_maxx == "pml" else None,
            ),
            slice(
                self.thickness_grid_miny + 1 if self.boundary_type_miny == "pml" else 0,
                -self.thickness_grid_maxy - 1 if self.boundary_type_maxy == "pml" else None,
            ),
            slice(
                self.thickness_grid_minz + 1 if self.boundary_type_minz == "pml" else 0,
                -self.thickness_grid_maxz - 1 if self.boundary_type_maxz == "pml" else None,
            ),
        )

    @classmethod
    def from_uniform_bound(
        cls,
        thickness: int = 10,
        boundary_type: str = "pml",
        kappa_start: float = 1,
        kappa_end: float = 1.5,
    ) -> "BoundaryConfig":
        """Creates a BoundaryConfig with uniform parameters for all boundaries.

        Args:
            thickness: Grid thickness to use for all PML boundaries
            boundary_type: Type of boundary to use ("pml" or "periodic"). Defaults to "pml".
            kappa_start: Initial kappa value for all boundaries. Defaults to 1.0.
            kappa_end: Final kappa value for all boundaries. Defaults to 1.5.

        Returns:
            BoundaryConfig: New config object with uniform parameters
        """
        return cls(
            boundary_type_minx=boundary_type,
            boundary_type_maxx=boundary_type,
            boundary_type_miny=boundary_type,
            boundary_type_maxy=boundary_type,
            boundary_type_minz=boundary_type,
            boundary_type_maxz=boundary_type,
            thickness_grid_minx=thickness,
            thickness_grid_maxx=thickness,
            thickness_grid_miny=thickness,
            thickness_grid_maxy=thickness,
            thickness_grid_minz=thickness,
            thickness_grid_maxz=thickness,
            kappa_start_minx=kappa_start,
            kappa_end_minx=kappa_end,
            kappa_start_maxx=kappa_start,
            kappa_end_maxx=kappa_end,
            kappa_start_miny=kappa_start,
            kappa_end_miny=kappa_end,
            kappa_start_maxy=kappa_start,
            kappa_end_maxy=kappa_end,
            kappa_start_minz=kappa_start,
            kappa_end_minz=kappa_end,
            kappa_start_maxz=kappa_start,
            kappa_end_maxz=kappa_end,
        )

from_uniform_bound(thickness=10, boundary_type='pml', kappa_start=1, kappa_end=1.5) classmethod

Creates a BoundaryConfig with uniform parameters for all boundaries.

Parameters:

Name Type Description Default
thickness int

Grid thickness to use for all PML boundaries

 10
boundary_type str

Type of boundary to use ("pml" or "periodic"). Defaults to "pml".

 'pml'
kappa_start float

Initial kappa value for all boundaries. Defaults to 1.0.

 1
kappa_end float

Final kappa value for all boundaries. Defaults to 1.5.

 1.5

Returns:

Name Type Description
BoundaryConfig BoundaryConfig

New config object with uniform parameters
Source code in src/fdtdx/objects/boundaries/initialization.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
@classmethod
def from_uniform_bound(
    cls,
    thickness: int = 10,
    boundary_type: str = "pml",
    kappa_start: float = 1,
    kappa_end: float = 1.5,
) -> "BoundaryConfig":
    """Creates a BoundaryConfig with uniform parameters for all boundaries.

    Args:
        thickness: Grid thickness to use for all PML boundaries
        boundary_type: Type of boundary to use ("pml" or "periodic"). Defaults to "pml".
        kappa_start: Initial kappa value for all boundaries. Defaults to 1.0.
        kappa_end: Final kappa value for all boundaries. Defaults to 1.5.

    Returns:
        BoundaryConfig: New config object with uniform parameters
    """
    return cls(
        boundary_type_minx=boundary_type,
        boundary_type_maxx=boundary_type,
        boundary_type_miny=boundary_type,
        boundary_type_maxy=boundary_type,
        boundary_type_minz=boundary_type,
        boundary_type_maxz=boundary_type,
        thickness_grid_minx=thickness,
        thickness_grid_maxx=thickness,
        thickness_grid_miny=thickness,
        thickness_grid_maxy=thickness,
        thickness_grid_minz=thickness,
        thickness_grid_maxz=thickness,
        kappa_start_minx=kappa_start,
        kappa_end_minx=kappa_end,
        kappa_start_maxx=kappa_start,
        kappa_end_maxx=kappa_end,
        kappa_start_miny=kappa_start,
        kappa_end_miny=kappa_end,
        kappa_start_maxy=kappa_start,
        kappa_end_maxy=kappa_end,
        kappa_start_minz=kappa_start,
        kappa_end_minz=kappa_end,
        kappa_start_maxz=kappa_start,
        kappa_end_maxz=kappa_end,
    )

get_dict()

Gets a dictionary mapping boundary names to their grid thicknesses.

Returns:

Type Description
dict[str, int]

dict[str, int]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z' mapping to their respective grid thickness values.
Source code in src/fdtdx/objects/boundaries/initialization.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def get_dict(self) -> dict[str, int]:
    """Gets a dictionary mapping boundary names to their grid thicknesses.

    Returns:
        dict[str, int]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
            mapping to their respective grid thickness values.
    """
    return {
        "min_x": self.thickness_grid_minx,
        "max_x": self.thickness_grid_maxx,
        "min_y": self.thickness_grid_miny,
        "max_y": self.thickness_grid_maxy,
        "min_z": self.thickness_grid_minz,
        "max_z": self.thickness_grid_maxz,
    }

get_inside_boundary_slice()

Gets slice objects for the non-PML interior region of the simulation volume.

Returns:

Type Description
tuple[slice, slice, slice]

tuple[slice, slice, slice]: Three slice objects for indexing the x, y, z dimensions respectively, excluding the PML boundary regions.
Source code in src/fdtdx/objects/boundaries/initialization.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
def get_inside_boundary_slice(self) -> tuple[slice, slice, slice]:
    """Gets slice objects for the non-PML interior region of the simulation volume.

    Returns:
        tuple[slice, slice, slice]: Three slice objects for indexing the x, y, z dimensions
            respectively, excluding the PML boundary regions.
    """
    return (
        slice(
            self.thickness_grid_minx + 1 if self.boundary_type_minx == "pml" else 0,
            -self.thickness_grid_maxx - 1 if self.boundary_type_maxx == "pml" else None,
        ),
        slice(
            self.thickness_grid_miny + 1 if self.boundary_type_miny == "pml" else 0,
            -self.thickness_grid_maxy - 1 if self.boundary_type_maxy == "pml" else None,
        ),
        slice(
            self.thickness_grid_minz + 1 if self.boundary_type_minz == "pml" else 0,
            -self.thickness_grid_maxz - 1 if self.boundary_type_maxz == "pml" else None,
        ),
    )

get_kappa_dict(prop)

Gets a dictionary mapping boundary names to their kappa values.

Parameters:

Name Type Description Default
prop Literal['kappa_start', 'kappa_end']

Which kappa property to get, either "kappa_start" or "kappa_end"

 required

Returns:

Type Description
dict[str, float]

dict[str, float]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z' mapping to their respective kappa values.

Raises:

Type Description
Exception

If prop is not "kappa_start" or "kappa_end"
Source code in src/fdtdx/objects/boundaries/initialization.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def get_kappa_dict(
    self,
    prop: Literal["kappa_start", "kappa_end"],
) -> dict[str, float]:
    """Gets a dictionary mapping boundary names to their kappa values.

    Args:
        prop: Which kappa property to get, either "kappa_start" or "kappa_end"

    Returns:
        dict[str, float]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
            mapping to their respective kappa values.

    Raises:
        Exception: If prop is not "kappa_start" or "kappa_end"
    """
    if prop == "kappa_start":
        return {
            "min_x": self.kappa_start_minx,
            "max_x": self.kappa_start_maxx,
            "min_y": self.kappa_start_miny,
            "max_y": self.kappa_start_maxy,
            "min_z": self.kappa_start_minz,
            "max_z": self.kappa_start_maxz,
        }
    elif prop == "kappa_end":
        return {
            "min_x": self.kappa_end_minx,
            "max_x": self.kappa_end_maxx,
            "min_y": self.kappa_end_miny,
            "max_y": self.kappa_end_maxy,
            "min_z": self.kappa_end_minz,
            "max_z": self.kappa_end_maxz,
        }
    else:
        raise Exception(f"Unknown: {prop=}")

get_type_dict()

Gets a dictionary mapping boundary names to their boundary types.

Returns:

Type Description
dict[str, str]

dict[str, str]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z' mapping to their respective boundary types ("pml" or "periodic").
Source code in src/fdtdx/objects/boundaries/initialization.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def get_type_dict(self) -> dict[str, str]:
    """Gets a dictionary mapping boundary names to their boundary types.

    Returns:
        dict[str, str]: Dictionary with keys 'min_x', 'max_x', 'min_y', 'max_y', 'min_z', 'max_z'
            mapping to their respective boundary types ("pml" or "periodic").
    """
    return {
        "min_x": self.boundary_type_minx,
        "max_x": self.boundary_type_maxx,
        "min_y": self.boundary_type_miny,
        "max_y": self.boundary_type_maxy,
        "min_z": self.boundary_type_minz,
        "max_z": self.boundary_type_maxz,
    }

Configuration object for specifying at which side of the simulation wich type of boundary should be used. Also allows specification of the PML thickness and other parameters if used.

fdtdx.objects.boundaries.boundary_objects_from_config(config, volume)

Creates boundary objects from a boundary configuration.

Creates PerfectlyMatchedLayer or PeriodicBoundary objects for all six boundaries (min/max x/y/z) based on the provided configuration. Also generates position constraints to properly place the boundary objects relative to the simulation volume.

Parameters:

Name Type Description Default
config BoundaryConfig

Configuration object containing boundary parameters

 required
volume SimulationObject

The main simulation volume object that the boundaries will surround

 required

Returns:

Type Description
tuple[dict[str, Union[PerfectlyMatchedLayer, PeriodicBoundary]], list[PositionConstraint]]

tuple containing: - dict mapping boundary names ('min_x', 'max_x', etc) to boundary objects - list of PositionConstraint objects for placing the boundaries
Source code in src/fdtdx/objects/boundaries/initialization.py 
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def boundary_objects_from_config(
    config: BoundaryConfig,
    volume: SimulationObject,
) -> tuple[dict[str, Union[PerfectlyMatchedLayer, PeriodicBoundary]], list[PositionConstraint]]:
    """Creates boundary objects from a boundary configuration.

    Creates PerfectlyMatchedLayer or PeriodicBoundary objects for all six boundaries
    (min/max x/y/z) based on the provided configuration. Also generates position
    constraints to properly place the boundary objects relative to the simulation volume.

    Args:
        config: Configuration object containing boundary parameters
        volume: The main simulation volume object that the boundaries will surround

    Returns:
        tuple containing:
            - dict mapping boundary names ('min_x', 'max_x', etc) to boundary objects
            - list of PositionConstraint objects for placing the boundaries
    """
    boundaries, constraints = {}, []
    thickness_dict = config.get_dict()
    type_dict = config.get_type_dict()
    kappa_start_dict = config.get_kappa_dict("kappa_start")
    kappa_end_dict = config.get_kappa_dict("kappa_end")

    for kind, thickness in thickness_dict.items():
        axis, direction = axis_direction_from_kind(kind)
        boundary_type = type_dict[kind]
        kappa_start, kappa_end = kappa_start_dict[kind], kappa_end_dict[kind]

        grid_shape_list: list[int | None] = [None, None, None]
        grid_shape_list[axis] = thickness if boundary_type == "pml" else 1
        grid_shape: PartialGridShape3D = tuple(grid_shape_list)  # type: ignore

        other_axes = [0, 1, 2]
        del other_axes[axis]

        if boundary_type == "pml":
            cur_boundary = PerfectlyMatchedLayer(
                axis=axis,
                partial_grid_shape=grid_shape,
                kappa_start=kappa_start,
                kappa_end=kappa_end,
                direction=direction,
            )
        else:  # periodic
            cur_boundary = PeriodicBoundary(
                axis=axis,
                partial_grid_shape=grid_shape,
                direction=direction,
            )

        direction_int = -1 if direction == "-" else 1
        pos_constraint = cur_boundary.place_relative_to(
            volume,
            axes=(axis, other_axes[0], other_axes[1]),
            own_positions=(direction_int, 0, 0),
            other_positions=(direction_int, 0, 0),
        )

        boundaries[kind] = cur_boundary
        constraints.append(pos_constraint)

    return boundaries, constraints

Initializes the corresponding boundary objects based on the config object.

Boundary Objects

fdtdx.objects.boundaries.PerfectlyMatchedLayer

Bases: BaseBoundary

Implements a Convolutional Perfectly Matched Layer (CPML) boundary condition.

The CPML absorbs outgoing electromagnetic waves with minimal reflection by using a complex coordinate stretching approach. This implementation supports arbitrary axis orientation and both positive/negative directions.

Attributes:

Name Type Description
axis int

Principal axis for PML (0=x, 1=y, 2=z)
direction Literal['+', '-']

Direction along axis ("+" or "-")
alpha float

Loss parameter for complex frequency shifting
kappa_start float

Initial kappa stretching coefficient
kappa_end float

Final kappa stretching coefficient
color tuple[float, float, float]

RGB color tuple for visualization
Source code in src/fdtdx/objects/boundaries/perfectly_matched_layer.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
@extended_autoinit
class PerfectlyMatchedLayer(BaseBoundary):
    """Implements a Convolutional Perfectly Matched Layer (CPML) boundary condition.

    The CPML absorbs outgoing electromagnetic waves with minimal reflection by using
    a complex coordinate stretching approach. This implementation supports arbitrary
    axis orientation and both positive/negative directions.

    Attributes:
        axis: Principal axis for PML (0=x, 1=y, 2=z)
        direction: Direction along axis ("+" or "-")
        alpha: Loss parameter for complex frequency shifting
        kappa_start: Initial kappa stretching coefficient
        kappa_end: Final kappa stretching coefficient
        color: RGB color tuple for visualization
    """

    axis: int = field(kind="KW_ONLY")  # type: ignore
    direction: Literal["+", "-"] = frozen_field(kind="KW_ONLY")  # type: ignore
    alpha: float = 1.0e-8
    kappa_start: float = 1.0
    kappa_end: float = 1.5
    color: tuple[float, float, float] = DARK_GREY

    @property
    def descriptive_name(self) -> str:
        """Gets a human-readable name describing this PML boundary's location.

        Returns:
            str: Description like "min_x" or "max_z" indicating position
        """
        axis_str = "x" if self.axis == 0 else "y" if self.axis == 1 else "z"
        direction_str = "min" if self.direction == "-" else "max"
        return f"{direction_str}_{axis_str}"

    @property
    def thickness(self) -> int:
        """Gets the thickness of the PML layer in grid points.

        Returns:
            int: Number of grid points in the PML along its axis
        """
        return self.grid_shape[self.axis]

    def init_state(
        self,
    ) -> BoundaryState:
        dtype = self._config.dtype
        sigma_E, sigma_H = standard_sigma_from_direction_axis(
            thickness=self.thickness,
            direction=self.direction,
            axis=self.axis,
            dtype=dtype,
        )

        kappa = kappa_from_direction_axis(
            kappa_start=self.kappa_start,
            kappa_end=self.kappa_end,
            thickness=self.thickness,
            direction=self.direction,
            axis=self.axis,
            dtype=dtype,
        )

        bE = jnp.exp(-self._config.courant_number * (sigma_E / kappa + self.alpha))
        bH = jnp.exp(-self._config.courant_number * (sigma_H / kappa + self.alpha))

        cE = (bE - 1) * sigma_E / (sigma_E * kappa + kappa**2 * self.alpha)
        cH = (bH - 1) * sigma_H / (sigma_H * kappa + kappa**2 * self.alpha)

        ext_shape = (3,) + self.grid_shape

        boundary_state = BoundaryState(
            psi_Ex=jnp.zeros(shape=ext_shape, dtype=dtype),
            psi_Ey=jnp.zeros(shape=ext_shape, dtype=dtype),
            psi_Ez=jnp.zeros(shape=ext_shape, dtype=dtype),
            psi_Hx=jnp.zeros(shape=ext_shape, dtype=dtype),
            psi_Hy=jnp.zeros(shape=ext_shape, dtype=dtype),
            psi_Hz=jnp.zeros(shape=ext_shape, dtype=dtype),
            bE=bE.astype(dtype),
            bH=bH.astype(dtype),
            cE=cE.astype(dtype),
            cH=cH.astype(dtype),
            kappa=kappa.astype(dtype),
        )
        return boundary_state

    def reset_state(self, state: BoundaryState) -> BoundaryState:
        dtype = self._config.dtype
        sigma_E, sigma_H = standard_sigma_from_direction_axis(
            thickness=self.thickness,
            direction=self.direction,
            axis=self.axis,
            dtype=dtype,
        )

        kappa = kappa_from_direction_axis(
            kappa_start=self.kappa_start,
            kappa_end=self.kappa_end,
            thickness=self.thickness,
            direction=self.direction,
            axis=self.axis,
            dtype=dtype,
        )

        bE = jnp.exp(-self._config.courant_number * (sigma_E / kappa + self.alpha))
        bH = jnp.exp(-self._config.courant_number * (sigma_H / kappa + self.alpha))

        cE = (bE - 1) * sigma_E / (sigma_E * kappa + kappa**2 * self.alpha)
        cH = (bH - 1) * sigma_H / (sigma_H * kappa + kappa**2 * self.alpha)

        new_state = BoundaryState(
            psi_Ex=state.psi_Ex * 0,
            psi_Ey=state.psi_Ey * 0,
            psi_Ez=state.psi_Ez * 0,
            psi_Hx=state.psi_Hx * 0,
            psi_Hy=state.psi_Hy * 0,
            psi_Hz=state.psi_Hz * 0,
            bE=bE.astype(dtype),
            bH=bH.astype(dtype),
            cE=cE.astype(dtype),
            cH=cH.astype(dtype),
            kappa=kappa.astype(dtype),
        )
        return new_state

    def boundary_interface_grid_shape(self) -> GridShape3D:
        if self.axis == 0:
            return 1, self.grid_shape[1], self.grid_shape[2]
        elif self.axis == 1:
            return self.grid_shape[0], 1, self.grid_shape[2]
        elif self.axis == 2:
            return self.grid_shape[0], self.grid_shape[1], 1
        raise Exception(f"Invalid axis: {self.axis=}")

    def boundary_interface_slice_tuple(self) -> SliceTuple3D:
        slice_list = [*self._grid_slice_tuple]
        if self.direction == "+":
            slice_list[self.axis] = (self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1)
        elif self.direction == "-":
            slice_list[self.axis] = (self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1])
        return slice_list[0], slice_list[1], slice_list[2]

    def boundary_interface_slice(self) -> Slice3D:
        slice_list = [*self.grid_slice]
        if self.direction == "+":
            slice_list[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )
        elif self.direction == "-":
            slice_list[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        return slice_list[0], slice_list[1], slice_list[2]

    def update_E_boundary_state(
        self,
        boundary_state: BoundaryState,
        H: jax.Array,
    ) -> BoundaryState:
        Hx = H[0, *self.grid_slice]
        Hy = H[1, *self.grid_slice]
        Hz = H[2, *self.grid_slice]

        psi_Ex = boundary_state.psi_Ex * boundary_state.bE
        psi_Ey = boundary_state.psi_Ey * boundary_state.bE
        psi_Ez = boundary_state.psi_Ez * boundary_state.bE

        psi_Ex = psi_Ex.at[1, :, 1:, :].add(
            (Hz[:, 1:, :] - Hz[:, :-1, :])
            * (boundary_state.cE[1, :, 1:, :] if self.axis == 1 else boundary_state.cE[1])
        )
        psi_Ex = psi_Ex.at[2, :, :, 1:].add(
            (Hy[:, :, 1:] - Hy[:, :, :-1])
            * (boundary_state.cE[2, :, :, 1:] if self.axis == 2 else boundary_state.cE[2])
        )

        psi_Ey = psi_Ey.at[2, :, :, 1:].add(
            (Hx[:, :, 1:] - Hx[:, :, :-1])
            * (boundary_state.cE[2, :, :, 1:] if self.axis == 2 else boundary_state.cE[2])
        )
        psi_Ey = psi_Ey.at[0, 1:, :, :].add(
            (Hz[1:, :, :] - Hz[:-1, :, :])
            * (boundary_state.cE[0, 1:, :, :] if self.axis == 0 else boundary_state.cE[0])
        )

        psi_Ez = psi_Ez.at[0, 1:, :, :].add(
            (Hy[1:, :, :] - Hy[:-1, :, :])
            * (boundary_state.cE[0, 1:, :, :] if self.axis == 0 else boundary_state.cE[0])
        )
        psi_Ez = psi_Ez.at[1, :, 1:, :].add(
            (Hx[:, 1:, :] - Hx[:, :-1, :])
            * (boundary_state.cE[1, :, 1:, :] if self.axis == 1 else boundary_state.cE[1])
        )

        boundary_state = boundary_state.at["psi_Ex"].set(psi_Ex)
        boundary_state = boundary_state.at["psi_Ey"].set(psi_Ey)
        boundary_state = boundary_state.at["psi_Ez"].set(psi_Ez)

        return boundary_state

    def update_H_boundary_state(
        self,
        boundary_state: BoundaryState,
        E: jax.Array,
    ) -> BoundaryState:
        Ex = E[0, *self.grid_slice]
        Ey = E[1, *self.grid_slice]
        Ez = E[2, *self.grid_slice]

        psi_Hx = boundary_state.psi_Hx * boundary_state.bH
        psi_Hy = boundary_state.psi_Hy * boundary_state.bH
        psi_Hz = boundary_state.psi_Hz * boundary_state.bH

        psi_Hx = psi_Hx.at[1, :, :-1, :].add(
            (Ez[:, 1:, :] - Ez[:, :-1, :])
            * (boundary_state.cH[1, :, :-1, :] if self.axis == 1 else boundary_state.cH[1])
        )
        psi_Hx = psi_Hx.at[2, :, :, :-1].add(
            (Ey[:, :, 1:] - Ey[:, :, :-1])
            * (boundary_state.cH[2, :, :, :-1] if self.axis == 2 else boundary_state.cH[2])
        )

        psi_Hy = psi_Hy.at[2, :, :, :-1].add(
            (Ex[:, :, 1:] - Ex[:, :, :-1])
            * (boundary_state.cH[2, :, :, :-1] if self.axis == 2 else boundary_state.cH[2])
        )
        psi_Hy = psi_Hy.at[0, :-1, :, :].add(
            (Ez[1:, :, :] - Ez[:-1, :, :])
            * (boundary_state.cH[0, :-1, :, :] if self.axis == 0 else boundary_state.cH[0])
        )

        psi_Hz = psi_Hz.at[0, :-1, :, :].add(
            (Ey[1:, :, :] - Ey[:-1, :, :])
            * (boundary_state.cH[0, :-1, :, :] if self.axis == 0 else boundary_state.cH[0])
        )
        psi_Hz = psi_Hz.at[1, :, :-1, :].add(
            (Ex[:, 1:, :] - Ex[:, :-1, :])
            * (boundary_state.cH[1, :, :-1, :] if self.axis == 1 else boundary_state.cH[1])
        )

        boundary_state = boundary_state.at["psi_Hx"].set(psi_Hx)
        boundary_state = boundary_state.at["psi_Hy"].set(psi_Hy)
        boundary_state = boundary_state.at["psi_Hz"].set(psi_Hz)

        return boundary_state

    def update_E(
        self,
        E: jax.Array,
        boundary_state: BoundaryState,
        inverse_permittivity: jax.Array,
    ) -> jax.Array:
        phi_Ex = boundary_state.psi_Ex[1] - boundary_state.psi_Ex[2]
        phi_Ey = boundary_state.psi_Ey[2] - boundary_state.psi_Ey[0]
        phi_Ez = boundary_state.psi_Ez[0] - boundary_state.psi_Ez[1]
        phi_E = jnp.stack((phi_Ex, phi_Ey, phi_Ez), axis=0)

        E = E.at[:, *self.grid_slice].divide(boundary_state.kappa)
        inv_perm_slice = inverse_permittivity[self.grid_slice]
        update = self._config.courant_number * inv_perm_slice * phi_E
        E = E.at[:, *self.grid_slice].add(update)
        return E

    def update_H(
        self,
        H: jax.Array,
        boundary_state: BoundaryState,
        inverse_permeability: jax.Array | float,
    ) -> jax.Array:
        phi_Hx = boundary_state.psi_Hx[1] - boundary_state.psi_Hx[2]
        phi_Hy = boundary_state.psi_Hy[2] - boundary_state.psi_Hy[0]
        phi_Hz = boundary_state.psi_Hz[0] - boundary_state.psi_Hz[1]
        phi_H = jnp.stack((phi_Hx, phi_Hy, phi_Hz), axis=0)

        H = H.at[:, *self.grid_slice].divide(boundary_state.kappa)
        if isinstance(inverse_permeability, jax.Array) and inverse_permeability.ndim > 0:
            inverse_permeability = inverse_permeability[self.grid_slice]
        update = -self._config.courant_number * inverse_permeability * phi_H
        H = H.at[:, *self.grid_slice].add(update)
        return H

descriptive_name: str property

Gets a human-readable name describing this PML boundary's location.

Returns:

Name Type Description
str str

Description like "min_x" or "max_z" indicating position

thickness: int property

Gets the thickness of the PML layer in grid points.

Returns:

Name Type Description
int int

Number of grid points in the PML along its axis

fdtdx.objects.boundaries.PeriodicBoundary

Bases: BaseBoundary

Implements periodic boundary conditions.

The periodic boundary connects opposite sides of the simulation domain, making waves that exit one side reenter from the opposite side.

Attributes:

Name Type Description
axis int

Principal axis for periodicity (0=x, 1=y, 2=z)
direction Literal['+', '-']

Direction along axis ("+" or "-")
color tuple[float, float, float]

RGB color tuple for visualization
Source code in src/fdtdx/objects/boundaries/periodic.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
@extended_autoinit
class PeriodicBoundary(BaseBoundary):
    """Implements periodic boundary conditions.

    The periodic boundary connects opposite sides of the simulation domain,
    making waves that exit one side reenter from the opposite side.

    Attributes:
        axis: Principal axis for periodicity (0=x, 1=y, 2=z)
        direction: Direction along axis ("+" or "-")
        color: RGB color tuple for visualization
    """

    axis: int = field(kind="KW_ONLY")  # type: ignore
    direction: Literal["+", "-"] = frozen_field(kind="KW_ONLY")  # type: ignore
    color: tuple[float, float, float] = LIGHT_BLUE

    @property
    def descriptive_name(self) -> str:
        """Gets a human-readable name describing this periodic boundary's location.

        Returns:
            str: Description like "min_x" or "max_z" indicating position
        """
        axis_str = "x" if self.axis == 0 else "y" if self.axis == 1 else "z"
        direction_str = "min" if self.direction == "-" else "max"
        return f"{direction_str}_{axis_str}"

    @property
    def thickness(self) -> int:
        """Gets the thickness of the periodic boundary layer in grid points.

        Returns:
            int: Number of grid points in the boundary layer (always 1 for periodic)
        """
        return 1

    def init_state(
        self,
    ) -> PeriodicBoundaryState:
        dtype = self._config.dtype
        ext_shape = (3,) + self.grid_shape

        boundary_state = PeriodicBoundaryState(
            E_opposite=jnp.zeros(shape=ext_shape, dtype=dtype),
            H_opposite=jnp.zeros(shape=ext_shape, dtype=dtype),
        )
        return boundary_state

    def reset_state(self, state: PeriodicBoundaryState) -> PeriodicBoundaryState:
        new_state = PeriodicBoundaryState(
            E_opposite=state.E_opposite * 0,
            H_opposite=state.H_opposite * 0,
        )
        return new_state

    def boundary_interface_grid_shape(self) -> GridShape3D:
        if self.axis == 0:
            return 1, self.grid_shape[1], self.grid_shape[2]
        elif self.axis == 1:
            return self.grid_shape[0], 1, self.grid_shape[2]
        elif self.axis == 2:
            return self.grid_shape[0], self.grid_shape[1], 1
        raise Exception(f"Invalid axis: {self.axis=}")

    def boundary_interface_slice_tuple(self) -> SliceTuple3D:
        slice_list = [*self._grid_slice_tuple]
        if self.direction == "+":
            slice_list[self.axis] = (self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1)
        elif self.direction == "-":
            slice_list[self.axis] = (self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1])
        return slice_list[0], slice_list[1], slice_list[2]

    def boundary_interface_slice(self) -> Slice3D:
        slice_list = [*self.grid_slice]
        if self.direction == "+":
            slice_list[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )
        elif self.direction == "-":
            slice_list[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        return slice_list[0], slice_list[1], slice_list[2]

    def update_E_boundary_state(
        self,
        boundary_state: PeriodicBoundaryState,
        H: jax.Array,
    ) -> PeriodicBoundaryState:
        # Get field values from opposite boundary
        opposite_slice = list(self.grid_slice)
        if self.direction == "+":
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        else:
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )

        # Store H field values from opposite boundary
        H_opposite = jnp.array(H[..., opposite_slice[0], opposite_slice[1], opposite_slice[2]])

        return PeriodicBoundaryState(
            E_opposite=boundary_state.E_opposite,  # Keep existing E values
            H_opposite=H_opposite,  # Update H values
        )

    def update_H_boundary_state(
        self,
        boundary_state: PeriodicBoundaryState,
        E: jax.Array,
    ) -> PeriodicBoundaryState:
        # Get field values from opposite boundary
        opposite_slice = list(self.grid_slice)
        if self.direction == "+":
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        else:
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )

        # Store E field values from opposite boundary
        E_opposite = jnp.array(E[..., opposite_slice[0], opposite_slice[1], opposite_slice[2]])

        return PeriodicBoundaryState(
            E_opposite=E_opposite,  # Update E values
            H_opposite=boundary_state.H_opposite,  # Keep existing H values
        )

    def update_E(
        self,
        E: jax.Array,
        boundary_state: PeriodicBoundaryState,
        inverse_permittivity: jax.Array,
    ) -> jax.Array:
        del boundary_state, inverse_permittivity
        # Get the boundary slice
        boundary_slice = list(self.grid_slice)
        if self.direction == "+":
            boundary_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )
            # Copy from opposite boundary (last slice)
            opposite_slice = list(self.grid_slice)
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        else:
            boundary_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
            # Copy from opposite boundary (first slice)
            opposite_slice = list(self.grid_slice)
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )

        # Copy field values from opposite boundary
        E = E.at[..., boundary_slice[0], boundary_slice[1], boundary_slice[2]].set(
            E[..., opposite_slice[0], opposite_slice[1], opposite_slice[2]]
        )

        return E

    def update_H(
        self,
        H: jax.Array,
        boundary_state: PeriodicBoundaryState,
        inverse_permeability: jax.Array | float,
    ) -> jax.Array:
        del boundary_state, inverse_permeability
        # Get the boundary slice
        boundary_slice = list(self.grid_slice)
        if self.direction == "+":
            boundary_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )
            # Copy from opposite boundary (last slice)
            opposite_slice = list(self.grid_slice)
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
        else:
            boundary_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][1] - 1, self._grid_slice_tuple[self.axis][1]
            )
            # Copy from opposite boundary (first slice)
            opposite_slice = list(self.grid_slice)
            opposite_slice[self.axis] = slice(
                self._grid_slice_tuple[self.axis][0], self._grid_slice_tuple[self.axis][0] + 1
            )

        # Copy field values from opposite boundary
        H = H.at[..., boundary_slice[0], boundary_slice[1], boundary_slice[2]].set(
            H[..., opposite_slice[0], opposite_slice[1], opposite_slice[2]]
        )

        return H

descriptive_name: str property

Gets a human-readable name describing this periodic boundary's location.

Returns:

Name Type Description
str str

Description like "min_x" or "max_z" indicating position

thickness: int property

Gets the thickness of the periodic boundary layer in grid points.

Returns:

Name Type Description
int int

Number of grid points in the boundary layer (always 1 for periodic)