From 59d8df56f95fee2556776d088b4553bfe0f70251 Mon Sep 17 00:00:00 2001 From: Casey Wojcik Date: Mon, 1 Dec 2025 09:02:23 -0800 Subject: [PATCH] feat: add EMESimulationData.coeffs (FXC-4275) --- CHANGELOG.md | 1 + schemas/EMESimulation.json | 4 + tests/sims/full_fdtd.h5 | Bin 479632 -> 479632 bytes tests/test_components/test_eme.py | 138 ++++++++++++++++++++++++- tidy3d/__init__.py | 8 ++ tidy3d/components/data/data_array.py | 45 ++++++++ tidy3d/components/eme/data/dataset.py | 138 ++++++++++++++++++++++++- tidy3d/components/eme/data/sim_data.py | 8 +- tidy3d/components/eme/simulation.py | 75 +++++++++++++- 9 files changed, 405 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f51376c4f0..bf12238276 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added `symmetrize_mirror`, `symmetrize_rotation`, `symmetrize_diagonal` functions to the autograd plugin. They can be used for enforcing symmetries in topology optimization. - Klayout plugin automatically finds klayout installation path at common locations. - Added autograd support for `TriangleMesh`, allowing gradient computation with respect to mesh vertices for inverse design. +- Added `EMESimulationData.coeffs` to store coefficients from the EME solver, including mode overlaps, interface S matrices, and effective propagation indices. ### Changed - Removed validator that would warn if `PerturbationMedium` values could become numerically unstable, since an error will anyway be raised if this actually happens when the medium is converted using actual perturbation data. diff --git a/schemas/EMESimulation.json b/schemas/EMESimulation.json index 42e7a1ac19..df704ebeeb 100644 --- a/schemas/EMESimulation.json +++ b/schemas/EMESimulation.json @@ -12879,6 +12879,10 @@ }, "type": "array" }, + "store_coeffs": { + "default": true, + "type": "boolean" + }, "store_port_modes": { "default": true, "type": "boolean" diff --git a/tests/sims/full_fdtd.h5 b/tests/sims/full_fdtd.h5 index 4f383a0e7317580c664d6f56a8776f8696f1e3fe..ebf763524cb0e4fada8b04664bab47710658e6ef 100644 GIT binary patch delta 2296 zcmYL~dsvL?8pdm?N#{gS6!A?5Q6WjH=j&v5LK2-OL_)H(Y$dTGD(RpSDjlSQs94x# zLd;BRX$`53N>iePRn-A84pZF-q#8HS%B--&1E`pq1h_u=zV8Cd{(Q-d1!{gRhn&f zXgWW?vtv&FhE5~!xS$hW=_LY_ql8BBHX$rDSKWAfJqIe(|5#TRj)8jkrRDRKWMCB+ zO=};;09lWj%j{*a)QfHq?)OWld5fv)j1n=%-_`R$&0|J<1sy) zxlbub3^9Md)d@^y0YRa6#8||I8PEQ7bF$Clz(KZ-FKT4bLzYCi-=0zzDO{C*Mpf9Zs){;O{$HynT zvPAH5mbjaezz4O{*KA8l1n@M)CC58i2xF`ukBllI+-2TN$`AXI^Z)s2Cntpm7_p$K z?yFVYYYBJ>kGhFBiD1~`<)>5)5qK~jw!}NIpxs)dnY@n!Csnscgg%qOlu_fgo9^uC zeV+8$-?G|_gLKr+R&39K3s)8q3aT)$be6rUcE=#CfgjRzO9niv%%_DhQt=%g=%dW`s<_|D} zN9#E9CLn{5lf7z!33-m&jAP~cakG8jEO_@6G6eY3L&Ruk$J6n9lzM zb-PZWx)*HXeAE$?cv$rp`h09#7VSvpu542dzfhj%1 z%8o&K&kcy{f*cOiu zxl(v+sO3o`$|(*@OQr{7a9r*8Vc&OB82;ARqVr^DrH7e7^N?` zOrMECkuuZ0%3umEZ(aEK@ZUHjT1cIbrBA`TU1jb3Mj>ELjMH5KQdnZ+a@=1*`cp-c zqIT?+z5tjf&nNLBgphT4o9#yt4qcvWIP(s0phZD@_kb9KsP;dWcP8RM9FLkHhOht= zst=KzB?WRvbaY~W&5@`3aR{qrCk5H?q0Yd2K;@+bZg+9x^Q}36WXz&ZsPI7BKGoF4 z7YFjd0XOzu8Emo&j?&6vLH_#}b16F*inmWe>sT*tdKH6$-LyW+Y9aLO-&b_nKne!M z2@7uAV!aa)4nRc9nQgVF> zClw)~YWMNHSsW;roX~Fe7Q*dtPkr$gHmnX+>fb!T2m3j$y#oM;-3+p5X5vW+jB|;TTw|HeM|(cUR}}WHt^ukW8FvT1RiV%BKq6L$)Pef z%s5TOp*^d$RLMmQ%pmd=t05s+?k#^lqAdh#@59nZ*JW^%>~NKGhyxb(o!r7RLI`>B zP&2|#0yN{|LtiX;K&)E1C%}sfgBia?ns1d+Is*uK(Lhc+RlIkc4Ql(_3nn7iFszdL zifPY-jqzI)=a%z<5Y`<#Lr(y&A94y1`3d4$MsfN|F?f(X%T06DKe%W(sby5sgb+y% z9Q>Nig^eq12S26Dpa)G!3AGeKPQWL@%o-`gezIM&HUon)+uuWj^Kj5Itv_uRCj42C z$rCBJQPxMT8h>YC-``E$GXI_#ML0%HY5e$2eUd@-iI2UQxVfL&ue%|iR`ii-Q9)cd z--|}Y?FhSb2rsDx@_J||da_92o^Lsz7jiXE#^wPp@ZR&vItk=0c5au&OTZ6%Tc>wd z07jNoO`Tj8G~kO_9$qU*~T%9$HBpb)U5)vHeI1hJ2ej+tU+{B`6Zy^-@8{Ph1S>a$Af*x1Q+O<}W zoK9ccV((=LdOx3WPiY+qoj3pSmk&1)wBp!ya*X`Vh^EJ+l_KbYw-uAN0#BYS;C9fgZ z=FH7DvNrcoj+~hsC8ZfNWA5WOueR(T&%f{2^LgIS>v=MZ)iR6KYFWcz@a9SH^;;Vd zi!Q=P88g89j=~d<6BLl$PZ~>$WP$oXqx+1`2AoTs_CRcfknoaPHW=IqIB~>B@v-4HX{y!)TzfSR7zKCxAsxZA_Sf0B}W0v!S~f zpgwhh6Jt(+lZH(mhjTb!y0%cFe1fyZ*&ZC$Mzq|3Fp@XjNR*@9+P+YoECf`Z_76FG z`Cx-Qnoi1~0p@&VYRoSJz+3ChE;_P++i{2zxou;Mqcb+;B`YR`Qu&#ZV_YU&TFQ^w zRwIPoW6Igq@_eWxX{NXx6N19irH;OS9w=NKzyIEm2bVWy$uZ;fEzY5rA!|;mLb&7S z?N?;W2D9bg(&PCw@QR$U$h*yiT Nf(sQIs$KBoEi6#Wk61}OBZSt;xllPn$~O-3 zI&bBz(Iy574YVId4`7k{Cj+redo1$vC#|*H3mC-T#3{@C6$U9fuuQx_!6G4K!nEZY z7Lk^@_<1A@gWyMpl)nh0mw$d2e2@s)YcxIV}e5HJj=Iw6XLb#cTN=v zAz9bIq8vj3Bt>%ZRT3YNsy2(tLcy1Y3P;ga;h*$h9%jJg#NMN4N*T}}9X*6M7C=tY zQfkdbCSc_+`{{>p;cCBDMY{_Zj69YjuI#1*s-z?%^|0_;E?}h=g$9x7KtwwnGfduu zXCWWt$4rIL&tK0}yUc@ggOdB?feWXLoYZiWbkKMmyKViX$p6RNwk9|J@%83rbkHJ& z51RSS$IX`*@Fw43T%v&nuW1k07sI*GwsN<0h)RWkV*H1^AU4Pou`hb8>vL9fM}c#t6v$S1a2M_V$1kyzz>M8rQO^EMr^ z6oSdnxx+CBx4a_>^@#*4qj+N59mNKN9~X$_8g~gBPuGcxE<28h4U&jgNJC2h@n0dj z?9p0@jvgQ;5`qYoo@+#0XM3m7ji#SA1?HtrCHp>Bt8hWAtG(;D6>t=DY-KF2JoT8cUPlP` zt8nK%u*giBX=TrrN$W@gZiEMog@>q6n=t%=CLOnP=mN1Gz z5QlM#;XBwbTbTQ~OMC|89s8{n?HRE8*Me2+OB(dWq27Av(V+)zk5?h_;ZVAtXIPd1 zRCsZ=GF4Q#ar(lF-Zn0@EZvm6y2Nk3=OB?N4#(BWcHU?~2&5vDf_tPG+on_i z89PIU(>z7D$o%^?F*^oe#7Fe~258`X?d(CjBGJuAlNapE-T=G17A+Rx-zeo(VTU4l zG?>a69Fl*`0gm4Z@^Ty>=F8AN7XRCTAd`$n=>ayd?!}I#dkH}KL4(`waXw_LzRPlz zW#@_#JxAONBy6b3qI2gwD6=`w0C&p>Yp)q0WJJ7dtgYh#Zf@9&hvI_OmE6zivI3}H zGbJB+&jzQ-U$3QIX9J!ztkUe;FJ5~KhCk`1o=;`NqbSU2#t$q&X`@u1ykLQ*j_~l^ zL>{<2)GOJUzyb5e9rHdm0x<9}jPaMF0FU5aV4Fhwdp3yYTe&&AZVck(UO90+3X5pW zmLw(jVUbBh)#PM17C~K7#TD|g$o7*h6Td5B5mniS#t|F_dG#Dyq>zb0nll|F62w7Q zBP8=|5(DNo_@-eW1mNB{o$q;#0gS1g+&3Of@M?;^(rwFw0z5%mp+Phx>avcd<|d5B zu6XlOzUq#bW&NA|BbE=9E3sK0wu@5NYm?e^g^>R2fiJp|4myjbTE905z{IVU@HUzQ z4pRDxM;!UEn#MS=D~Ar0(OKPn1wsgCnm1=&;Q|SB_hDlXAF4YR%y$1l1zk*+Q8JDL zX(^sV+Gzsls(C(3?W9ARR6~`RB(HhIO(H80vWwCMnxa=#&C_!UsH1~#Wse5y86PrK z7W|TmsnE2B`t#rz54s0~X(Mqwa8t75Dj(&5{*4|TbUFR&(&-r}5v!*~CNn*pwf1lT z7smNlMwt-yJCHaeG!2rC@yr@|7UVI_3e1cta7|cC_dh2DoaDf5eBDO`!(ERH7NfJq+TuQ%*sz0fU4PV=u}Dj7pL694Eb=8w3d=-NR8pdn O7S*<9mWfQB>i+@UZi{yS diff --git a/tests/test_components/test_eme.py b/tests/test_components/test_eme.py index 058ae384ac..b1ba84b853 100644 --- a/tests/test_components/test_eme.py +++ b/tests/test_components/test_eme.py @@ -285,6 +285,9 @@ def test_eme_monitor(): def test_eme_simulation(): sim = make_eme_sim() + # no log except deprecated coeffs monitor + with AssertLogLevel(None): + _ = sim.updated_copy(monitors=[sim.monitors[0], *list(sim.monitors[2:])]) _ = sim.plot(x=0, ax=AX) _ = sim.plot(y=0, ax=AX) _ = sim.plot(z=0, ax=AX) @@ -311,16 +314,16 @@ def test_eme_simulation(): # test warning for not providing wavelength in autogrid grid_spec = td.GridSpec.auto(min_steps_per_wvl=20) + sim = sim.updated_copy(grid_spec=grid_spec) with AssertLogLevel("INFO", contains_str="wavelength"): - sim = sim.updated_copy(grid_spec=grid_spec) + _ = sim.updated_copy(monitors=[]) # multiple freqs are ok, but not for autogrid _ = sim.updated_copy( grid_spec=td.GridSpec.uniform(dl=0.2), freqs=[10000000000.0, *list(sim.freqs)] ) with AssertLogLevel("INFO", contains_str="wavelength"): _ = sim.updated_copy( - freqs=[*list(sim.freqs), 10000000000.0], - grid_spec=grid_spec, + freqs=[*list(sim.freqs), 10000000000.0], grid_spec=grid_spec, monitors=[] ) # test port offsets @@ -460,6 +463,17 @@ def test_eme_simulation(): ) with AssertLogLevel("WARNING", contains_str="estimated storage"): sim_bad.validate_pre_upload() + # coeffs warning + sim_bad = sim.updated_copy( + size=(10, 10, 10), + monitors=[], + store_port_modes=False, + freqs=list(1e14 * np.linspace(1, 2, 100)), + eme_grid_spec=td.EMEUniformGrid(mode_spec=td.EMEModeSpec(num_modes=100), num_cells=100), + grid_spec=sim.grid_spec.updated_copy(wavelength=1), + ) + with AssertLogLevel("WARNING", contains_str="store_coeffs"): + sim_bad.validate_pre_upload() sim_bad = sim.updated_copy( size=(10, 10, 10), monitors=[large_monitor], @@ -558,12 +572,14 @@ def test_eme_simulation(): constraint="passive", eme_grid_spec=td.EMEUniformGrid(num_cells=1, mode_spec=td.EMEModeSpec(num_modes=40)), grid_spec=sim.grid_spec.updated_copy(wavelength=1), + monitors=[], ) sim_good.validate_pre_upload() sim_good = sim.updated_copy( constraint=None, eme_grid_spec=td.EMEUniformGrid(num_cells=1, mode_spec=td.EMEModeSpec(num_modes=60)), grid_spec=sim.grid_spec.updated_copy(wavelength=1), + monitors=[], ) sim_good.validate_pre_upload() # warn about num modes with constraint @@ -731,6 +747,47 @@ def _get_eme_smatrix_data_array(num_modes_in=2, num_modes_out=3, num_freqs=2, nu return smatrix_entry +def _get_eme_interface_smatrix_data_array( + num_modes_in=2, num_modes_out=3, num_freqs=2, num_sweep=0 +): + if num_modes_in != 0: + mode_index_in = np.arange(num_modes_in) + else: + mode_index_in = [0] + if num_modes_out != 0: + mode_index_out = np.arange(num_modes_out) + else: + mode_index_out = [0] + if num_sweep != 0: + sweep_index = np.arange(num_sweep) + else: + sweep_index = [0] + eme_cell_index = np.arange(3) + + f = td.C_0 * np.linspace(1, 2, num_freqs) + + data = (1 + 1j) * np.random.random( + (len(f), len(sweep_index), len(eme_cell_index), len(mode_index_out), len(mode_index_in)) + ) + coords = { + "f": f, + "sweep_index": sweep_index, + "eme_cell_index": eme_cell_index, + "mode_index_out": mode_index_out, + "mode_index_in": mode_index_in, + } + smatrix_entry = td.EMEInterfaceSMatrixDataArray(data, coords=coords) + + if num_modes_in == 0: + smatrix_entry = smatrix_entry.drop_vars("mode_index_in") + if num_modes_out == 0: + smatrix_entry = smatrix_entry.drop_vars("mode_index_out") + if num_sweep == 0: + smatrix_entry = smatrix_entry.drop_vars("sweep_index") + + return smatrix_entry + + def _get_eme_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): S11 = _get_eme_smatrix_data_array( num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep @@ -747,6 +804,35 @@ def _get_eme_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): return td.EMESMatrixDataset(S11=S11, S12=S12, S21=S21, S22=S22) +def _get_eme_interface_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): + S11 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + S12 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + S21 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + S22 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + return td.EMEInterfaceSMatrixDataset(S11=S11, S12=S12, S21=S21, S22=S22) + + +def _get_eme_overlaps_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): + O11 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + O12 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + O21 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + return td.EMEOverlapDataset(O11=O11, O12=O12, O21=O21) + + def _get_eme_coeff_data_array(num_sweep=0): f = [2e14] mode_index_out = [0, 1] @@ -787,7 +873,26 @@ def _get_eme_coeff_data_array(num_sweep=0): def _get_eme_coeff_dataset(num_sweep=0): A = _get_eme_coeff_data_array(num_sweep=num_sweep) B = _get_eme_coeff_data_array(num_sweep=num_sweep) - return td.EMECoefficientDataset(A=A, B=B) + flux = _get_eme_flux_data_array(num_sweep=num_sweep) + n_complex = _get_eme_mode_index_data_array(num_sweep=num_sweep) + interface_smatrices = _get_eme_interface_smatrix_dataset(num_sweep=num_sweep) + overlaps = _get_eme_overlaps_dataset(num_sweep=num_sweep) + return td.EMECoefficientDataset( + A=A, + B=B, + flux=flux, + n_complex=n_complex, + interface_smatrices=interface_smatrices, + overlaps=overlaps, + ) + + +def test_eme_normalize_coeff_dataset(): + coeffs = _get_eme_coeff_dataset() + coeffs_normalized = coeffs.normalized_copy + assert coeffs_normalized.flux is None + with pytest.raises(ValidationError): + _ = coeffs_normalized.normalized_copy def test_eme_coeff_data_array(): @@ -819,6 +924,29 @@ def _get_eme_mode_index_data_array(num_sweep=0): return data +def _get_eme_flux_data_array(num_sweep=0): + f = [td.C_0, 3e14] + mode_index = np.arange(10) + eme_cell_index = np.arange(7) + if num_sweep != 0: + sweep_index = np.arange(num_sweep) + else: + sweep_index = [0] + coords = { + "f": f, + "sweep_index": sweep_index, + "eme_cell_index": eme_cell_index, + "mode_index": mode_index, + } + data = td.EMEFluxDataArray( + np.random.random((len(f), len(sweep_index), len(eme_cell_index), len(mode_index))), + coords=coords, + ) + if num_sweep == 0: + data = data.drop_vars("sweep_index") + return data + + def test_eme_mode_index_data_array(): _ = _get_eme_mode_index_data_array() @@ -1305,7 +1433,7 @@ def test_eme_periodicity(): # remove the field monitor, now it passes desired_cell_index_pairs = set([(i, i + 1) for i in range(6)] + [(5, 1)]) - with AssertLogLevel(None): + with AssertLogLevel("WARNING", contains_str="deprecated"): sim = sim.updated_copy( monitors=[m for m in sim.monitors if not isinstance(m, td.EMEFieldMonitor)] ) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 18f7d99047..13934377f7 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -185,6 +185,8 @@ ChargeDataArray, DiffractionDataArray, EMECoefficientDataArray, + EMEFluxDataArray, + EMEInterfaceSMatrixDataArray, EMEModeIndexDataArray, EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, @@ -243,7 +245,9 @@ from .components.eme.data.dataset import ( EMECoefficientDataset, EMEFieldDataset, + EMEInterfaceSMatrixDataset, EMEModeSolverDataset, + EMEOverlapDataset, EMESMatrixDataset, ) from .components.eme.data.monitor_data import EMECoefficientData, EMEFieldData, EMEModeSolverData @@ -605,8 +609,11 @@ def set_logging_level(level: str) -> None: "EMEFieldData", "EMEFieldDataset", "EMEFieldMonitor", + "EMEFluxDataArray", "EMEFreqSweep", "EMEGrid", + "EMEInterfaceSMatrixDataArray", + "EMEInterfaceSMatrixDataset", "EMELengthSweep", "EMEModeIndexDataArray", "EMEModeSolverData", @@ -615,6 +622,7 @@ def set_logging_level(level: str) -> None: "EMEModeSpec", "EMEModeSweep", "EMEMonitor", + "EMEOverlapDataset", "EMEPeriodicitySweep", "EMESMatrixDataArray", "EMESMatrixDataset", diff --git a/tidy3d/components/data/data_array.py b/tidy3d/components/data/data_array.py index 47f90e6035..7c7cbae61f 100644 --- a/tidy3d/components/data/data_array.py +++ b/tidy3d/components/data/data_array.py @@ -1185,6 +1185,31 @@ class EMESMatrixDataArray(DataArray): _data_attrs = {"long_name": "scattering matrix element"} +class EMEInterfaceSMatrixDataArray(DataArray): + """Scattering matrix elements at a single cell interface for a fixed pair of ports, + possibly with an extra sweep index. + Example + ------- + >>> mode_index_in = [0, 1] + >>> mode_index_out = [0, 1, 2] + >>> eme_cell_index = [2, 4] + >>> f = [2e14] + >>> sweep_index = np.arange(10) + >>> coords = dict( + ... f=f, + ... sweep_index=sweep_index, + ... eme_cell_index=eme_cell_index, + ... mode_index_out=mode_index_out, + ... mode_index_in=mode_index_in, + ... ) + >>> fd = EMEInterfaceSMatrixDataArray((1 + 1j) * np.random.random((1, 10, 2, 3, 2)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "sweep_index", "eme_cell_index", "mode_index_out", "mode_index_in") + _data_attrs = {"long_name": "scattering matrix element"} + + class EMEModeIndexDataArray(DataArray): """Complex-valued effective propagation index of an EME mode, also indexed by EME cell. @@ -1203,6 +1228,24 @@ class EMEModeIndexDataArray(DataArray): _data_attrs = {"long_name": "Propagation index"} +class EMEFluxDataArray(DataArray): + """Power flux of an EME mode, also indexed by EME cell. + + Example + ------- + >>> f = [2e14, 3e14] + >>> sweep_index = np.arange(2) + >>> eme_cell_index = np.arange(5) + >>> mode_index = np.arange(4) + >>> coords = dict(f=f, sweep_index=sweep_index, eme_cell_index=eme_cell_index, mode_index=mode_index) + >>> data = EMEFluxDataArray(np.random.random((2,2,5,4)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "sweep_index", "eme_cell_index", "mode_index") + _data_attrs = {"units": WATT, "long_name": "flux"} + + class ChargeDataArray(DataArray): """Charge data array. @@ -1592,8 +1635,10 @@ def _make_impedance_data_array(result: DataArray) -> ImpedanceResultType: EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, EMESMatrixDataArray, + EMEInterfaceSMatrixDataArray, EMECoefficientDataArray, EMEModeIndexDataArray, + EMEFluxDataArray, EMEFreqModeDataArray, ChargeDataArray, SteadyVoltageDataArray, diff --git a/tidy3d/components/eme/data/dataset.py b/tidy3d/components/eme/data/dataset.py index 33f45cef2c..3b3e7a3722 100644 --- a/tidy3d/components/eme/data/dataset.py +++ b/tidy3d/components/eme/data/dataset.py @@ -2,16 +2,23 @@ from __future__ import annotations +from typing import Optional + +import numpy as np import pydantic.v1 as pd +from tidy3d.components.base import cached_property from tidy3d.components.data.data_array import ( EMECoefficientDataArray, + EMEFluxDataArray, + EMEInterfaceSMatrixDataArray, EMEModeIndexDataArray, EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, EMESMatrixDataArray, ) from tidy3d.components.data.dataset import Dataset, ElectromagneticFieldDataset +from tidy3d.exceptions import ValidationError class EMESMatrixDataset(Dataset): @@ -39,22 +46,143 @@ class EMESMatrixDataset(Dataset): ) +class EMEInterfaceSMatrixDataset(Dataset): + """Dataset storing S matrices associated with EME cell interfaces.""" + + S11: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S11 matrix", + description="S matrix relating output modes at port 1 to input modes at port 1.", + ) + S12: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S12 matrix", + description="S matrix relating output modes at port 1 to input modes at port 2.", + ) + S21: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S21 matrix", + description="S matrix relating output modes at port 2 to input modes at port 1.", + ) + S22: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S22 matrix", + description="S matrix relating output modes at port 2 to input modes at port 2.", + ) + + +class EMEOverlapDataset(Dataset): + """Dataset storing overlaps between EME modes. + + ``Oij`` is the unconjugated overlap computed using the E field of cell ``i`` + and the H field of cell ``j``. + + For consistency with ``Sij``, ``mode_index_out`` refers to the mode index + in cell ``i``, and ``mode_index_in`` refers to the mode index in cell ``j``. + """ + + O11: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O11 matrix", + description="Overlap integral between E field and H field in the same cell.", + ) + O12: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O12 matrix", + description="Overlap integral between E field on side 1 and H field on side 2.", + ) + O21: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O21 matrix", + description="Overlap integral between E field on side 2 and H field on side 1.", + ) + + class EMECoefficientDataset(Dataset): - """Dataset storing expansion coefficients for the modes in a cell. + """Dataset storing various coefficients related to the EME simulation. + These coefficients can be used for debugging or optimization. + + The ``A`` and ``B`` fields store the expansion coefficients for the modes in a cell. These are defined at the cell centers. + + The ``n_complex`` and ``flux`` fields respectively store the complex-valued effective + propagation index and the power flux associated with the modes used in the + EME calculation. + + The ``interface_Sij`` fields store the S matrices associated with the interfaces + between EME cells. """ - A: EMECoefficientDataArray = pd.Field( - ..., + A: Optional[EMECoefficientDataArray] = pd.Field( + None, title="A coefficient", description="Coefficient for forward mode in this cell.", ) - B: EMECoefficientDataArray = pd.Field( - ..., + + B: Optional[EMECoefficientDataArray] = pd.Field( + None, title="B coefficient", description="Coefficient for backward mode in this cell.", ) + n_complex: Optional[EMEModeIndexDataArray] = pd.Field( + None, + title="Propagation Index", + description="Complex-valued effective propagation indices associated with the EME modes.", + ) + + flux: Optional[EMEFluxDataArray] = pd.Field( + None, + title="Flux", + description="Power flux of the EME modes.", + ) + + interface_smatrices: Optional[EMEInterfaceSMatrixDataset] = pd.Field( + None, + title="Interface S Matrices", + description="S matrices associated with the interfaces between EME cells.", + ) + + overlaps: Optional[EMEOverlapDataset] = pd.Field( + None, title="Overlaps", description="Overlaps between EME modes." + ) + + @cached_property + def normalized_copy(self) -> EMECoefficientDataset: + """Return a copy of the ``EMECoefficientDataset`` where + the ``A`` and ``B`` coefficients as well as the ``interface_smatrices`` + are normalized by flux.""" + if self.flux is None: + raise ValidationError( + "The 'flux' field of the 'EMECoefficientDataset' is 'None', " + "so normalization cannot be performed." + ) + fields = {"A": self.A, "B": self.B} + flux_out = self.flux.rename(mode_index="mode_index_out") + for key, field in fields.items(): + if field is not None: + fields[key] = field * np.sqrt(flux_out) + if self.interface_smatrices is not None: + num_cells = len(self.flux.eme_cell_index) + flux1 = self.flux.isel(eme_cell_index=np.arange(num_cells - 1)) + flux2 = self.flux.isel(eme_cell_index=np.arange(1, num_cells)) + flux2 = flux2.assign_coords(eme_cell_index=np.arange(num_cells - 1)) + interface_S12 = self.interface_smatrices.S12 + flux_out = flux1.rename(mode_index="mode_index_out") + flux_in = flux2.rename(mode_index="mode_index_in") + interface_S12 = interface_S12 * np.sqrt(flux_out / flux_in) + interface_S21 = self.interface_smatrices.S21 + flux_out = flux2.rename(mode_index="mode_index_out") + flux_in = flux1.rename(mode_index="mode_index_in") + interface_S21 = interface_S21 * np.sqrt(flux_out / flux_in) + + fields["interface_smatrices"] = self.interface_smatrices.updated_copy( + S12=interface_S12, S21=interface_S21 + ) + # for safety to prevent normalizing twice + fields["flux"] = None + return self.updated_copy(**fields) + class EMEFieldDataset(ElectromagneticFieldDataset): """Dataset storing scalar components of E and H fields as a function of freq, mode_index, and port_index.""" diff --git a/tidy3d/components/eme/data/sim_data.py b/tidy3d/components/eme/data/sim_data.py index 05c1e6e739..750c088f2a 100644 --- a/tidy3d/components/eme/data/sim_data.py +++ b/tidy3d/components/eme/data/sim_data.py @@ -17,7 +17,7 @@ from tidy3d.exceptions import SetupError from tidy3d.log import log -from .dataset import EMESMatrixDataset +from .dataset import EMECoefficientDataset, EMESMatrixDataset from .monitor_data import EMEFieldData, EMEModeSolverData, EMEMonitorDataType @@ -39,6 +39,12 @@ class EMESimulationData(AbstractYeeGridSimulationData): None, title="S Matrix", description="Scattering matrix of the EME simulation." ) + coeffs: Optional[EMECoefficientDataset] = pd.Field( + None, + title="Coefficients", + description="Coefficients from the EME simulation. Useful for debugging and optimization.", + ) + port_modes_raw: Optional[EMEModeSolverData] = pd.Field( None, title="Port Modes", diff --git a/tidy3d/components/eme/simulation.py b/tidy3d/components/eme/simulation.py index 9f9c0841c7..808e32bf32 100644 --- a/tidy3d/components/eme/simulation.py +++ b/tidy3d/components/eme/simulation.py @@ -48,6 +48,7 @@ MAX_SIMULATION_DATA_SIZE_GB = 50 WARN_MODE_NUM_CELLS = 1e5 MAX_MODE_NUM_CELLS = 5e6 +WARN_COEFF_DATA_SIZE_GB = 0.5 # eme specific simulation parameters @@ -233,6 +234,13 @@ class EMESimulation(AbstractYeeGridSimulation): "Required to find scattering matrix in basis besides the computational basis.", ) + store_coeffs: bool = pd.Field( + True, + title="Store Coefficients", + description="Whether to store the internal coefficients from the EME simulation. " + "The results are stored in 'EMESimulationData.coeffs'.", + ) + normalize: bool = pd.Field( True, title="Normalize Scattering Matrix", @@ -752,6 +760,12 @@ def _validate_sweep_spec(self) -> None: def _validate_monitor_setup(self) -> None: """Check monitor setup.""" for i, monitor in enumerate(self.monitors): + if isinstance(monitor, EMECoefficientMonitor): + log.warning( + "'EMECoefficientMonitor' is deprecated. " + "The full coefficient data is stored in " + "'EMESimulationData.coeffs'." + ) if isinstance(monitor, EMEMonitor): _ = self._monitor_eme_cell_indices(monitor=monitor) if ( @@ -854,11 +868,70 @@ def _validate_monitor_size(self) -> None: total_size_gb += monitor_size_gb + # coefficients + if self.store_coeffs: + coeffs_size_b = 0 + bytes_complex = 8 + num_freqs = len(self.freqs) + num_modes = self.max_num_modes + num_eme_cells = self.eme_grid.num_cells + num_sweep = self._num_sweep + # A and B coefficients + coeffs_size_b += ( + 4 * bytes_complex * num_freqs * num_modes * num_modes * num_eme_cells * num_sweep + ) + # interface smatrices + coeffs_size_b += ( + 4 + * bytes_complex + * num_freqs + * num_modes + * num_modes + * (num_eme_cells - 1) + * self._num_sweep_interfaces + ) + # n_complex and flux + coeffs_size_b += ( + 2 * bytes_complex * num_freqs * num_modes * num_eme_cells * self._num_sweep_modes + ) + # overlaps + coeffs_size_b += ( + 2 + * bytes_complex + * num_freqs + * num_modes + * num_modes + * (num_eme_cells - 1) + * self._num_sweep_modes + ) + # self-overlaps + coeffs_size_b += ( + bytes_complex + * num_freqs + * num_modes + * num_modes + * num_eme_cells + * self._num_sweep_modes + ) + + coeffs_size_gb = coeffs_size_b / 1e9 + if coeffs_size_gb > WARN_COEFF_DATA_SIZE_GB: + log.warning( + "Simulation 'coeffs' have estimated storage size " + f"{coeffs_size_gb:1.2f}GB. " + "Consider setting 'store_coeffs=False' " + "or reducing the number of frequencies, modes, " + "EME cells, or sweep indices." + ) + + total_size_gb += coeffs_size_gb + if total_size_gb > MAX_SIMULATION_DATA_SIZE_GB: raise SetupError( f"Simulation's monitors have {total_size_gb:.2f}GB of estimated storage, " f"a maximum of {MAX_SIMULATION_DATA_SIZE_GB:.2f}GB are allowed. Note that " - "this estimate includes the port modes if 'store_port_modes' is 'True'." + "this estimate includes the port modes if 'store_port_modes' is 'True' " + "and the 'coeffs' if 'store_coeffs' is 'True'." ) # Make sure that internal storage from mode solvers also does not exceed the limit.