Skip to content

Commit

Permalink
Added a validator to enable FieldProjectionCartesianMonitor for 2D si…
Browse files Browse the repository at this point in the history
…mulations
  • Loading branch information
QimingFlex authored and momchil-flex committed Aug 5, 2024
1 parent cf4dfe7 commit 3de1e67
Show file tree
Hide file tree
Showing 5 changed files with 248 additions and 33 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Error if field projection monitors found in 2D simulations, except `FieldProjectionAngleMonitor` with `far_field_approx = True`. Support for other monitors and for exact field projection will be coming in a subsequent Tidy3D version.
- Mode solver now always operates on a reduced simulation copy.
- Moved `EMESimulation` size limit validators to preupload.
- Error if field projection monitors found in 2D simulations, except `FieldProjectionAngleMonitor` or `FieldProjectionCartesianMonitor` with `far_field_approx = True`. Support for other monitors and for exact field projection will be coming in a subsequent Tidy3D version.

### Fixed
- Error when loading a previously run `Batch` or `ComponentModeler` containing custom data.
Expand Down
189 changes: 189 additions & 0 deletions tests/test_components/test_field_projection.py
Original file line number Diff line number Diff line change
Expand Up @@ -370,3 +370,192 @@ def test_proj_clientside():
val.sel(f=f0)
with pytest.raises(DataError):
exact_fields_cartesian.renormalize_fields(proj_distance=5e6)


def make_2d_proj_monitors(center, size, freqs, plane):
"""Helper function to make near-to-far monitors in 2D simulations."""

if plane == "xy":
thetas = [np.pi / 2]
phis = np.linspace(0, 2 * np.pi, 100)
far_size = 10 * WAVELENGTH
Ns = 40
xs = np.linspace(-far_size, far_size, Ns)
ys = [0]
projection_axis = 0
elif plane == "yz":
thetas = np.linspace(0, np.pi, 1)
phis = [np.pi / 2]
far_size = 10 * WAVELENGTH
Ns = 40
xs = [0]
ys = np.linspace(-far_size, far_size, Ns)
projection_axis = 1
elif plane == "xz":
thetas = np.linspace(0, np.pi, 100)
phis = [0]
far_size = 10 * WAVELENGTH
Ns = 40
xs = [0]
ys = np.linspace(-far_size, far_size, Ns)
projection_axis = 0
else:
raise ValueError("Invalid plane. Use 'xy', 'yz', or 'xz'.")

n2f_angle_monitor_2d = td.FieldProjectionAngleMonitor(
center=center,
size=size,
freqs=freqs,
name="far_field_angle",
phi=list(phis),
theta=list(thetas),
proj_distance=R_FAR,
far_field_approx=True, # Fields are far enough for geometric far field approximations
)

n2f_car_monitor_2d = td.FieldProjectionCartesianMonitor(
center=center,
size=size,
freqs=freqs,
name="far_field_cartesian",
x=list(xs),
y=list(ys),
proj_axis=projection_axis,
proj_distance=R_FAR,
far_field_approx=True, # Fields are far enough for geometric far field approximations
)

return (n2f_angle_monitor_2d, n2f_car_monitor_2d)


def make_2d_proj(plane):
center = (0, 0, 0)
f0 = 1e13

if plane == "xy":
sim_size = (5, 5, 0)
monitor_size = (0, 2, td.inf)
# boundary conditions
boundary_conds = td.BoundarySpec(
x=td.Boundary.pml(),
y=td.Boundary.pml(),
z=td.Boundary.periodic(),
)
# data coordinates
x = np.array([0.0])
y = np.linspace(-1, 1, 10)
z = np.array([0.0])
coords = dict(x=x, y=y, z=z, f=[f0])
scalar_field = td.ScalarFieldDataArray(
(1 + 1j) * np.random.random((1, 10, 1, 1)), coords=coords
)
elif plane == "yz":
sim_size = (0, 5, 5)
monitor_size = (td.inf, 0, 2)
# boundary conditions
boundary_conds = td.BoundarySpec(
x=td.Boundary.periodic(),
y=td.Boundary.pml(),
z=td.Boundary.pml(),
)
# data coordinates
x = np.array([0.0])
y = np.array([0.0])
z = np.linspace(-1, 1, 10)
coords = dict(x=x, y=y, z=z, f=[f0])
scalar_field = td.ScalarFieldDataArray(
(1 + 1j) * np.random.random((1, 1, 10, 1)), coords=coords
)
elif plane == "xz":
sim_size = (5, 0, 5)
monitor_size = (0, td.inf, 2)
# boundary conditions
boundary_conds = td.BoundarySpec(
x=td.Boundary.pml(),
y=td.Boundary.periodic(),
z=td.Boundary.pml(),
)
# data coordinates
x = np.array([0.0])
y = np.array([0.0])
z = np.linspace(-1, 1, 10)
coords = dict(x=x, y=y, z=z, f=[f0])
scalar_field = td.ScalarFieldDataArray(
(1 + 1j) * np.random.random((1, 1, 10, 1)), coords=coords
)
else:
raise ValueError("Invalid plane. Use 'xy', 'yz', or 'xz'.")

monitor = td.FieldMonitor(
center=center, size=monitor_size, freqs=[f0], name="near_field", colocate=False
)

# Set up the simulation
sim = td.Simulation(
size=sim_size,
grid_spec=td.GridSpec.auto(wavelength=td.C_0 / f0),
boundary_spec=boundary_conds,
monitors=[monitor],
run_time=1e-12,
)

data = td.FieldData(
monitor=monitor,
Ex=scalar_field,
Ey=scalar_field,
Ez=scalar_field,
Hx=scalar_field,
Hy=scalar_field,
Hz=scalar_field,
symmetry=sim.symmetry,
symmetry_center=sim.center,
grid_expanded=sim.discretize_monitor(monitor),
)

sim_data = td.SimulationData(simulation=sim, data=(data,))

proj = td.FieldProjector.from_near_field_monitors(
sim_data=sim_data,
near_monitors=[monitor],
normal_dirs=["+"],
)

# make near-to-far monitors
(
n2f_angle_monitor_2d,
n2f_cart_monitor_2d,
) = make_2d_proj_monitors(center, monitor_size, [f0], plane)

far_fields_angular_2d = proj.project_fields(n2f_angle_monitor_2d)
far_fields_cartesian_2d = proj.project_fields(n2f_cart_monitor_2d)

# compute far field quantities
far_fields_angular_2d.r
far_fields_angular_2d.theta
far_fields_angular_2d.phi
far_fields_angular_2d.fields_spherical
far_fields_angular_2d.fields_cartesian
far_fields_angular_2d.radar_cross_section
far_fields_angular_2d.power
for val in far_fields_angular_2d.field_components.values():
val.sel(f=f0)
far_fields_angular_2d.renormalize_fields(proj_distance=5e6)

far_fields_cartesian_2d.x
far_fields_cartesian_2d.y
far_fields_cartesian_2d.z
far_fields_cartesian_2d.fields_spherical
far_fields_cartesian_2d.fields_cartesian
far_fields_cartesian_2d.radar_cross_section
far_fields_cartesian_2d.power
for val in far_fields_cartesian_2d.field_components.values():
val.sel(f=f0)
far_fields_cartesian_2d.renormalize_fields(proj_distance=5e6)


def test_2d_proj_clientside():
# Run simulations and tests for all three planes
planes = ["xy", "yz", "xz"]

for plane in planes:
make_2d_proj(plane)
2 changes: 1 addition & 1 deletion tidy3d/components/data/monitor_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -2161,7 +2161,7 @@ def propagation_factor(dist: Union[float, None], k: complex, is_2d_simulation: b
return 1.0

if is_2d_simulation:
return -np.exp(1j * k * dist) * np.sqrt(1j * k / (8 * np.pi * dist))
return np.exp(1j * k * dist) * np.sqrt(-1j * k / (8 * np.pi * dist))

return -1j * k * np.exp(1j * k * dist) / (4 * np.pi * dist)

Expand Down
10 changes: 7 additions & 3 deletions tidy3d/components/field_projection.py
Original file line number Diff line number Diff line change
Expand Up @@ -413,10 +413,15 @@ def _far_fields_for_surface(
_, source_names = surface.monitor.pop_axis(("x", "y", "z"), axis=surface.axis)

# integration dimension for 2d far field projection
zero_dim = (dim for dim, size in enumerate(self.sim_data.simulation.size) if size == 0)
zero_dim = [dim for dim, size in enumerate(self.sim_data.simulation.size) if size == 0]
if self.is_2d_simulation:
# Ensure zero_dim has a single element since {zero_dim} expects a value
if len(zero_dim) != 1:
raise ValueError("Expected exactly one dimension with size 0 for 2D simulation")

zero_dim = zero_dim[0]
integration_axis = {0, 1, 2} - {zero_dim, surface.axis}
idx_int_1d = integration_axis.pop() # Get the remaining axis as an integer
idx_int_1d = integration_axis.pop()

idx_u, idx_v = idx_uv
cmp_1, cmp_2 = source_names
Expand Down Expand Up @@ -451,7 +456,6 @@ def integrate_for_one_theta(i_th: int):
J[idx_u, i_th, j_ph] = self.integrate_1d(
currents_f[f"E{cmp_1}"].values, phase_ij, pts[idx_int_1d]
)

J[idx_v, i_th, j_ph] = self.integrate_1d(
currents_f[f"E{cmp_2}"].values, phase_ij, pts[idx_int_1d]
)
Expand Down
79 changes: 50 additions & 29 deletions tidy3d/components/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -2882,11 +2882,11 @@ def _projection_monitors_distance(cls, val, values):
def _projection_monitors_2d(cls, val, values):
"""
Validate if the field projection monitor is set up for a 2D simulation and
ensure the observation angle is configured correctly.
ensure the observation parameters are configured correctly.
- For a 2D simulation in the x-y plane, 'theta' should be set to 'pi/2'.
- For a 2D simulation in the y-z plane, 'phi' should be set to 'pi/2'.
- For a 2D simulation in the x-z plane, 'phi' should be set to '0'.
- For a 2D simulation in the y-z plane, 'phi' should be set to 'pi/2' or '2*pi/3'.
- For a 2D simulation in the x-z plane, 'phi' should be set to '0' or 'pi'.
Note: Exact far field projection is not available yet. Currently, only
'far_field_approx = True' is supported.
Expand All @@ -2903,18 +2903,13 @@ def _projection_monitors_2d(cls, val, values):
return val

plane = None
phi_value = None
theta_value = None

if sim_size[0] == 0:
plane = "y-z"
phi_value = np.pi / 2
elif sim_size[1] == 0:
plane = "x-z"
phi_value = 0
elif sim_size[2] == 0:
plane = "x-y"
theta_value = np.pi / 2

for monitor in val:
if isinstance(monitor, AbstractFieldProjectionMonitor):
Expand All @@ -2923,37 +2918,63 @@ def _projection_monitors_2d(cls, val, values):
f"Monitor '{monitor.name}' is not supported in 1D simulations."
)

if isinstance(
monitor, (FieldProjectionCartesianMonitor, FieldProjectionKSpaceMonitor)
):
if isinstance(monitor, (FieldProjectionKSpaceMonitor)):
raise SetupError(
f"Monitor '{monitor.name}' in 2D simulations is coming soon. "
"Please use 'FieldProjectionAngleMonitor' instead."
"Please use 'FieldProjectionAngleMonitor' or 'FieldProjectionCartesianMonitor' instead."
)

if isinstance(monitor, FieldProjectionAngleMonitor):
if not monitor.far_field_approx:
raise SetupError(
"Exact far field projection in 2D simulations is coming soon."
"Please set 'far_field_approx = True'."
)
if plane == "y-z" and (len(monitor.phi) != 1 or monitor.phi[0] != phi_value):
if isinstance(monitor, (FieldProjectionCartesianMonitor)):
config = {
"y-z": {"valid_proj_axes": [1, 2], "coord": ["x", "x"]},
"x-z": {"valid_proj_axes": [0, 2], "coord": ["x", "y"]},
"x-y": {"valid_proj_axes": [0, 1], "coord": ["y", "y"]},
}[plane]

valid_proj_axes = config["valid_proj_axes"]
invalid_proj_axis = [i for i in range(3) if i not in valid_proj_axes]

if monitor.proj_axis in invalid_proj_axis:
raise SetupError(
"For a 2D simulation in the y-z plane, the observation angle 'phi' "
f"of monitor '{monitor.name}' should be set to 'np.pi/2'."
f"For a 2D simulation in the {plane} plane, the 'proj_axis' of "
f"monitor '{monitor.name}' should be set to one of {valid_proj_axes}."
)
elif plane == "x-z" and (len(monitor.phi) != 1 or monitor.phi[0] != phi_value):
raise SetupError(
"For a 2D simulation in the x-z plane, the observation angle 'phi' "
f"of monitor '{monitor.name}' should be set to '0'."

for idx, axis in enumerate(valid_proj_axes):
coord = getattr(monitor, config["coord"][idx])
if monitor.proj_axis == axis and not all(value in [0] for value in coord):
raise SetupError(
f"For a 2D simulation in the {plane} plane with "
f"'proj_axis = {monitor.proj_axis}', '{config['coord'][idx]}' of monitor "
f"'{monitor.name}' should be set to '[0]'."
)

if isinstance(monitor, FieldProjectionAngleMonitor):
config = {
"y-z": {"valid_value": [np.pi / 2, 3 * np.pi / 2], "coord": "phi"},
"x-z": {"valid_value": [0, np.pi], "coord": "phi"},
"x-y": {"valid_value": [np.pi / 2], "coord": "theta"},
}[plane]

coord = getattr(monitor, config["coord"])
if not all(value in config["valid_value"] for value in coord):
replacements = {
np.pi: "np.pi",
np.pi / 2: "np.pi/2",
3 * np.pi / 2: "3*np.pi/2",
0: "0",
}
valid_values_str = ", ".join(
replacements.get(val) for val in config["valid_value"]
)
elif plane == "x-y" and (
len(monitor.theta) != 1 or monitor.theta[0] != theta_value
):
raise SetupError(
"For a 2D simulation in the x-y plane, the observation angle 'theta' "
f"of monitor '{monitor.name}' should be set to 'np.pi/2'."
f"For a 2D simulation in the {plane} plane, the observation "
f"angle '{config['coord']}' of monitor "
f"'{monitor.name}' should be set to "
f"'{valid_values_str}'"
)

return val

@pydantic.validator("monitors", always=True)
Expand Down

0 comments on commit 3de1e67

Please sign in to comment.