feat: flux and poynting properties for FieldProjectionCartesianData

[yaugenst-flex]

Fixes #1870

[momchil-flex]

So actually I think this is going to produce the same result as the previous version (could you check?). I think the issue is that we are providing the total power going through the surface in all directions, while the user's expectation for a cartesian monitor would be to get the power that's normal to the plane, like we do for other field monitors. In the field projection cartesian monitor, this normal direction would be the proj_axis.

@tomflexcompute fyi

[yaugenst-flex]

@momchil-flex yeah that's correct, they are the same. In particular, field projection power as implemented (before this) calculates $S_r$, and if $E_r = 0$ this corresponds exactly to the calculation in Cartesian too. Actually the calculation as it currently is makes sense and is also the convention? Probably the user should use the cartesian field projection monitor, that should already do the right thing?
Maybe there was also some confusion about obtaining total power, because you can't just sum up the fluxes but need to integrate over the area elements? @tomflexcompute

But yeah we can probably close this

[tomflexcompute]

@momchil-flex yeah that's correct, they are the same. In particular, field projection power as implemented (before this) calculates Sr, and if Er=0 this corresponds exactly to the calculation in Cartesian too. Actually the calculation as it currently is makes sense and is also the convention? Probably the user should use the cartesian field projection monitor, that should already do the right thing? Maybe there was also some confusion about obtaining total power, because you can't just sum up the fluxes but need to integrate over the area elements? @tomflexcompute

But yeah we can probably close this

Not only integrating over the area but taking the dot product of the S components with the surface normal, which is a bit inconvenient. If we don't change the definition of power, maybe add another method flux that directly returns the net flux through the projected surface?

[momchil-flex]

Yeah I think the point is that the cartesian projected data really looks like a regular FieldData, for which .flux returns the flux through the plane of the monitor, rather than the summed up poynting vector. I think @tomflexcompute 's suggestion of introducing a flux property then makes sense.

[yaugenst-flex]

So should the user not call .power at all for those monitors? I was thinking that for cartesian monitors, projected_fields.power could return the flux through the monitor plane, and for angular monitors, it should return the current value? A .flux property probably makes more sense, but should we warn the user if they call .power on field projection data from FieldProjectionAngleMonitors?

[tomflexcompute]

For angle monitor and k space monitor I think the original implementation of power is fine and useful. The issue is only for the cartesian monitor.

[yaugenst-flex]

I added flux and poynting properties to the cartesian field projection data now @momchil-flex

[momchil-flex]

Great, thanks. Could you try a real test where you place a FieldMonitor at some distance, and a cartesian field projection monitor to the same location, and compare the two .flux-es? There could be a small difference because of some of the colocation/integration details, but it should be almost the same.

[yaugenst-flex]

I set up a little simulation to test this:

This is the script:

import matplotlib.pyplot as plt
import numpy as np

import tidy3d as td
from tidy3d.web import run

is_2d = False

lcen = 1.0
min_steps_per_wvl = 20
monitor_buffer = 0.5
sim_size_z = 20.0
fcen = td.C_0 / lcen
fwidth = fcen / 10

box = td.Structure(
    geometry=td.Box(center=(0, 0, -sim_size_z / 2 + 1), size=(2.0, 2.0, 0.5)),
    medium=td.Medium(permittivity=2),
)

source = td.PointDipole(
    center=box.geometry.center,
    source_time=td.GaussianPulse(freq0=fcen, fwidth=fwidth),
    polarization="Ey",
)

sim_size_xy = np.array(box.geometry.size[:2]) + 4 * monitor_buffer

near_monitor = td.FieldMonitor(
    center=(0, 0, box.geometry.bounds[1][2] + monitor_buffer),
    size=(td.inf, td.inf, 0),
    freqs=(fcen,),
    colocate=False,
    name="near_fields",
)

flux_monitor = td.FluxMonitor(
    center=(0, 0, sim_size_z / 2 - monitor_buffer),
    size=(3, 3, 0),
    freqs=near_monitor.freqs,
    name="flux",
)

sim = td.Simulation(
    size=(sim_size_xy[0], 0 if is_2d else sim_size_xy[1], sim_size_z),
    center=(0, 0, 0),
    sources=(source,),
    monitors=(near_monitor, flux_monitor),
    structures=(box,),
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
    boundary_spec=td.BoundarySpec.pml(x=True, y=not is_2d, z=True),
    run_time=1e-12,
)

sim_data = run(sim, task_name="projection_test")

sx, sy = flux_monitor.size[:2]
px = np.linspace(-sx / 2, sx / 2, 11)
py = np.linspace(-sy / 2, sy / 2, 11)

monitor_far = td.FieldProjectionCartesianMonitor(
    center=near_monitor.center,
    size=near_monitor.size,
    freqs=near_monitor.freqs,
    x=px,
    y=py,
    proj_distance=flux_monitor.center[2] - near_monitor.center[2],
    proj_axis=2,
    far_field_approx=False,
    name="far_field",
)

projector = td.FieldProjector.from_near_field_monitors(
    sim_data=sim_data,
    near_monitors=[near_monitor],
    normal_dirs=["+"],
)

projected_fields = projector.project_fields(monitor_far)

print("real flux: ", sim_data["flux"].flux.item())
print("proj flux: ", projected_fields.flux.item())

And this outputs:

real flux:  0.6501566171646118
proj flux:  0.7093686717166056

This seems to be a pretty significant difference? Not sure what to make of that. What do you think? Note that it's not really related to the discretization of the projection monitor (px, py), I can 10x that and the flux changes only in the third digit.

[momchil-flex]

Note that the user-side field projector does not support far_field_approx=False. However, it seems like we don't explicitly error about this - could you make a small PR to fix that?

I modified your script to use server-side projection, and I also modified the monitor to be 3D surrounding the box, because there could be some effects from the PML in your setup. Then I got

real flux:  0.6501566171646118
proj flux:  0.6609804962995253

Increasing resolution to 30 yields

real flux:  0.6668662428855896
proj flux:  0.6695659347379561

So generally, this seems good. I think the difference with resolution may be due to a couple things, like 1. the fact that fdtd field propagation is not the same as the analytical field propagation, and 2. details of the numerical integration on the far field plane.

Here's the updated code:

import matplotlib.pyplot as plt
import numpy as np

import tidy3d as td
from tidy3d.web import run

is_2d = False

lcen = 1.0
min_steps_per_wvl = 30
monitor_buffer = 0.5
sim_size_z = 20.0
fcen = td.C_0 / lcen
fwidth = fcen / 10

box = td.Structure(
    geometry=td.Box(center=(0, 0, -sim_size_z / 2 + 1), size=(2.0, 2.0, 0.5)),
    medium=td.Medium(permittivity=2),
)

source = td.PointDipole(
    center=box.geometry.center,
    source_time=td.GaussianPulse(freq0=fcen, fwidth=fwidth),
    polarization="Ey",
)

sim_size_xy = np.array(box.geometry.size[:2]) + 4 * monitor_buffer

near_monitor = td.FieldMonitor(
    center=(0, 0, box.geometry.bounds[1][2] + monitor_buffer),
    size=(td.inf, td.inf, 0),
    freqs=(fcen,),
    colocate=False,
    name="near_fields",
)

flux_monitor = td.FluxMonitor(
    center=(0, 0, sim_size_z / 2 - monitor_buffer),
    size=(3, 3, 0),
    freqs=near_monitor.freqs,
    name="flux",
)

sx, sy = flux_monitor.size[:2]
px = np.linspace(-sx / 2, sx / 2, 11)
py = np.linspace(-sy / 2, sy / 2, 11)

monitor_far = td.FieldProjectionCartesianMonitor(
    center=box.geometry.center,
    size=(2.5, 2.5, 1),
    freqs=near_monitor.freqs,
    x=px,
    y=py,
    proj_distance=flux_monitor.center[2] - box.geometry.center[2],
    proj_axis=2,
    far_field_approx=False,
    name="far_field",
)

sim = td.Simulation(
    size=(sim_size_xy[0], 0 if is_2d else sim_size_xy[1], sim_size_z),
    center=(0, 0, 0),
    sources=(source,),
    monitors=(near_monitor, flux_monitor, monitor_far),
    structures=(box,),
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
    boundary_spec=td.BoundarySpec.pml(x=True, y=not is_2d, z=True),
    run_time=1e-12,
)

sim_data = run(sim, task_name="projection_test")

print("real flux: ", sim_data["flux"].flux.item())
print("proj flux: ", sim_data["far_field"].flux.item())