PLaSK library

Loading...
Searching...
No Matches
expansion3d.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_PW3D_H
15#define PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_PW3D_H
16
17#include <plask/plask.hpp>
18
19#include "../expansion.hpp"
20#include "fft.hpp"
21
22 namespace plask { namespace optical { namespace modal {
23
24 struct FourierSolver3D;
25
26 struct PLASK_SOLVER_API GradientFunctions: public MultiFieldProperty<double> {
27 enum EnumType {
28 COS2 = 0,
29 COSSIN = 1
30 };
31 static constexpr size_t NUM_VALS = 2;
32 static constexpr const char* NAME = "refractive index gradient functions cos² and cos·sin";
33 static constexpr const char* UNIT = "-";
34};
35
36 struct PLASK_SOLVER_API ExpansionPW3D: public Expansion {
37
38 dcomplex klong,
39 ktran;
40
41 size_t Nl,
42 Nt;
43 size_t nNl,
44 nNt;
45 size_t eNl,
46 eNt;
47
48 double left;
49 double right;
50 double back;
51 double front;
52 bool periodic_long,
53 periodic_tran;
54 bool initialized;
55
56 Component symmetry_long,
57 symmetry_tran;
58
59 size_t pil,
60 pir,
61 pif,
62 pib;
63
65 struct Coeff {
66 dcomplex c22, c00, ic00, c11, ic11, c01;
67 Coeff() {}
68 Coeff(const Coeff&) = default;
69 Coeff(const dcomplex& val): c22(val),
70 c00(val), ic00((val != 0.)? 1. / val : 0.),
71 c11(val), ic11((val != 0.)? 1. / val : 0.),
72 c01(0.) {}
73 Coeff(const Tensor3<dcomplex>& eps): c22(eps.c22), c00(eps.c00), c11(eps.c11), c01(eps.c01) {
74 dcomplex det = eps.c00 * eps.c11 - eps.c01 * conj(eps.c01);
75 ic00 = (det != 0.)? eps.c11 / det : 0.;
76 ic11 = (det != 0.)? eps.c00 / det : 0.;
77 }
78 Coeff& operator*=(dcomplex a) {
79 c22 *= a; c00 *= a; ic00 *= a; c11 *= a; ic11 *= a; c01 *= a;
80 return *this;
81 }
82 Coeff& operator=(const Tensor3<dcomplex>& eps) {
83 c22 = eps.c22;
84 c00 = eps.c00; ic00 = (eps.c00 != 0.)? 1. / eps.c00 : 0.;
85 c11 = eps.c11; ic11 = (eps.c11 != 0.)? 1. / eps.c11 : 0.;
86 c01 = eps.c01;
87 return *this;
88 }
89 bool differs(const Coeff& other) const {
90 return !(is_zero(other.c22-c22) && is_zero(other.c00-c00) && is_zero(other.c11-c11) && is_zero(other.c01-c01));
91 }
92 operator Tensor3<dcomplex>() const { return Tensor3<dcomplex>(c00, c11, c22, c01); }
93
94 Tensor3<dcomplex> toInverseTensor() const {
95 return Tensor3<dcomplex>(ic00, ic11, c22, c01);
96 // return Tensor3<dcomplex>(ic00, ic11, c22, -c01 / (c00 * c11 - c01 * conj(c01)));
97 }
98 };
99
101 std::vector<DataVector<Coeff>> coeffs;
102
104 struct Gradient {
105 struct Vertex;
106 dcomplex c2, cs;
107 Gradient(const Gradient&) = default;
108 Gradient(double c2, double cs): c2(c2), cs(cs) {}
109 Gradient(const Vec<2>& norm): c2(norm.c0 * norm.c0), cs(norm.c0 * norm.c1) {}
110 Gradient& operator=(const Gradient& norm) = default;
111 Gradient& operator=(const Vec<2>& norm) {
112 // double f = 1. / (norm.c0*norm.c0 + norm.c1*norm.c1);
113 c2 = /* f * */norm.c0 * norm.c0;
114 cs = /* f * */norm.c0 * norm.c1;
115 return *this;
116 }
117 Gradient operator*(double f) const { return Gradient(c2.real() * f, cs.real() * f); }
118 Gradient operator/(double f) const { return Gradient(c2.real() / f, cs.real() / f); }
119 Gradient& operator+=(const Gradient& norm) { c2 += norm.c2; cs += norm.cs; return *this; }
120 Gradient& operator*=(double f) { c2 *= f; cs *= f; return *this; }
121 Gradient& operator/=(double f) { c2 /= f; cs /= f; return *this; }
122 // Gradient operator/(size_t n) const { double f = 1. / double(n); return Gradient(c2.real() * f, cs.real() * f); }
123 bool isnan() const { return ::isnan(c2.real()); }
124 };
125
127 std::vector<DataVector<Gradient>> gradients;
128
130 std::vector<cmatrix> coeffs_ezz, coeffs_dexx, coeffs_deyy;
131
133 std::vector<bool> diagonals;
134
136 shared_ptr<RectangularMesh<3>> mesh;
137
142 ExpansionPW3D(FourierSolver3D* solver);
143
145 bool symmetric_long() const { return symmetry_long != E_UNSPECIFIED; }
146
148 bool symmetric_tran() const { return symmetry_tran != E_UNSPECIFIED; }
149
154 void init();
155
157 void reset();
158
159 bool diagonalQE(size_t l) const override {
160 return diagonals[l];
161 }
162
163 size_t matrixSize() const override { return 2*Nl*Nt; }
164
165 void getMatrices(size_t l, cmatrix& RE, cmatrix& RH) override;
166
167 void prepareField() override;
168
169 void cleanupField() override;
170
171 LazyData<Vec<3,dcomplex>> getField(size_t l,
172 const shared_ptr<const typename LevelsAdapter::Level>& level,
173 const cvector& E, const cvector& H) override;
174
175 LazyData<Tensor3<dcomplex>> getMaterialEps(size_t lay,
176 const shared_ptr<const typename LevelsAdapter::Level>& level,
177 InterpolationMethod interp) override;
178
179 LazyData<double> getGradients(GradientFunctions::EnumType what,
180 const shared_ptr<const typename LevelsAdapter::Level>& level,
181 InterpolationMethod interp);
182
183 double integrateField(WhichField field, size_t layer, const cmatrix& TE, const cmatrix& TH,
184 const std::function<std::pair<dcomplex,dcomplex>(size_t, size_t)>& vertical) override;
185
186 double integratePoyntingVert(const cvector& E, const cvector& H) override;
187
188 void getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal& gamma) override;
189
190 private:
191
192 DataVector<Vec<3,dcomplex>> field;
193 FFT::Backward2D fft_x, fft_y, fft_z;
194
195 void addToeplitzMatrix(cmatrix& work, int ordl, int ordt, size_t lay, int c, char syml, char symt, double a = 1.) {
196 for (int it = (symt ? 0 : -ordt); it <= ordt; ++it) {
197 size_t It = (it >= 0)? it : it + Nt;
198 for (int il = (syml ? 0 : -ordl); il <= ordl; ++il) {
199 size_t Il = (il >= 0)? il : il + Nl;
200 for (int jt = -ordt; jt <= ordt; ++jt) {
201 size_t Jt = (jt >= 0)? jt : (symt)? -jt : jt + Nt;
202 int ijt = it - jt; if (symt && ijt < 0) ijt = -ijt;
203 for (int jl = -ordl; jl <= ordl; ++jl) {
204 size_t Jl = (jl >= 0)? jl : (syml)? -jl : jl + Nl;
205 double f = 1.;
206 if (syml && jl < 0) { f *= syml; }
207 if (symt && jt < 0) { f *= symt; }
208 int ijl = il - jl; if (syml && ijl < 0) ijl = -ijl;
209 work(Nl * It + Il, Nl * Jt + Jl) += a * f * eps(lay, ijl, ijt, c);
210 }
211 }
212 }
213 }
214 }
215
216 void makeToeplitzMatrix(cmatrix& work, int ordl, int ordt, size_t lay, int c, char syml, char symt, double a = 1.) {
217 zero_matrix(work);
218 addToeplitzMatrix(work, ordl, ordt, lay, c, syml, symt, a);
219 }
220
221 void makeToeplitzMatrix(cmatrix& workc2, cmatrix& workcs, const DataVector<Gradient>& norms,
222 int ordl, int ordt, char syml, char symt) {
223 zero_matrix(workc2);
224 zero_matrix(workcs);
225 for (int it = (symt ? 0 : -ordt); it <= ordt; ++it) {
226 size_t It = (it >= 0)? it : it + Nt;
227 for (int il = (syml ? 0 : -ordl); il <= ordl; ++il) {
228 size_t Il = (il >= 0)? il : il + Nl;
229 size_t I = Nl * It + Il;
230 for (int jt = -ordt; jt <= ordt; ++jt) {
231 size_t Jt = (jt >= 0)? jt : (symt)? -jt : jt + Nt;
232 int ijt = it - jt; if (ijt < 0) ijt = symt? -ijt : ijt + nNt;
233 for (int jl = -ordl; jl <= ordl; ++jl) {
234 size_t Jl = (jl >= 0)? jl : (syml)? -jl : jl + Nl;
235 double fc = 1., fs = 1.;
236 if (syml && jl < 0) { fc *= syml; fs *= -syml; }
237 if (symt && jt < 0) { fc *= symt; fs *= -symt; }
238 int ijl = il - jl; if (ijl < 0) ijl = syml? -ijl : ijl + nNl;
239 size_t J = Nl * Jt + Jl, ij = nNl * ijt + ijl;
240 workc2(I,J) += fc * norms[ij].c2;
241 workcs(I,J) += fs * norms[ij].cs;
242 }
243 }
244 }
245 }
246 }
247
248 protected:
249
250 DataVector<Tensor2<dcomplex>> mag_long;
251 DataVector<Tensor2<dcomplex>> mag_tran;
252
253 FFT::Forward2D matFFT;
254 FFT::Forward2D cos2FFT, cssnFFT;
255
256 void beforeLayersIntegrals(dcomplex lam, dcomplex glam) override;
257
258 void layerIntegrals(size_t layer, double lam, double glam) override;
259
260 public:
261
262 dcomplex getKlong() const { return klong; }
263 void setKlong(dcomplex k) {
264 if (k != klong) {
265 klong = k;
266 solver->clearFields();
267 }
268 }
269
270 dcomplex getKtran() const { return ktran; }
271 void setKtran(dcomplex k) {
272 if (k != ktran) {
273 ktran = k;
274 solver->clearFields();
275 }
276 }
277
278 Component getSymmetryLong() const { return symmetry_long; }
279 void setSymmetryLong(Component sym) {
280 if (sym != symmetry_long) {
281 symmetry_long = sym;
282 solver->clearFields();
283 }
284 }
285
286 Component getSymmetryTran() const { return symmetry_tran; }
287 void setSymmetryTran(Component sym) {
288 if (sym != symmetry_tran) {
289 symmetry_tran = sym;
290 solver->clearFields();
291 }
292 }
293
295 size_t idx(int l, int t) {
296 if (l < 0) { if (symmetric_long()) l = -l; else l += int(Nl); }
297 if (t < 0) { if (symmetric_tran()) t = -t; else t += int(Nt); }
298 assert(0 <= l && std::size_t(l) < Nl);
299 assert(0 <= t && std::size_t(t) < Nt);
300 return Nl * t + l;
301 }
302
304 size_t idxe(int l, int t) {
305 if (l < 0) l += int(eNl);
306 if (t < 0) t += int(eNt);
307 assert(0 <= l && std::size_t(l) < eNl);
308 assert(0 <= t && std::size_t(t) < eNt);
309 return eNl * t + l;
310 }
311
313 dcomplex epsxx(size_t lay, int l, int t) {
314 if (l < 0) l += int(nNl);
315 if (t < 0) t += int(nNt);
316 return coeffs[lay][nNl * t + l].c00;
317 }
318
320 dcomplex epsyy(size_t lay, int l, int t) {
321 if (l < 0) l += int(nNl);
322 if (t < 0) t += int(nNt);
323 return coeffs[lay][nNl * t + l].c11;
324 }
325
327 dcomplex iepszz(size_t lay, int l, int t) {
328 if (l < 0) l += int(nNl);
329 if (t < 0) t += int(nNt);
330 return coeffs[lay][nNl * t + l].c22;
331 }
332
334 dcomplex iepszz(size_t lay, int il, int jl, int it, int jt) {
335 if (jl < 0) { if (symmetric_long()) return 0.; else jl += int(Nl); }
336 if (jt < 0) { if (symmetric_tran()) return 0.; else jt += int(Nt); }
337 if (il < 0) il += int(Nl);
338 if (it < 0) it += int(Nt);
339 return coeffs_ezz[lay](Nl * it + il, Nl * jt + jl);
340 }
341
343 dcomplex epsxy(size_t lay, int l, int t) {
344 if (l < 0) l += int(nNl);
345 if (t < 0) t += int(nNt);
346 return coeffs[lay][nNl * t + l].c01;
347 }
348
350 dcomplex eps(size_t lay, int l, int t, int c) {
351 if (l < 0) l += int(nNl);
352 if (t < 0) t += int(nNt);
353 return *(reinterpret_cast<dcomplex*>(coeffs[lay].data() + (nNl * t + l)) + c);
354 }
355
357 dcomplex epsyx(size_t lay, int l, int t) { return conj(epsxy(lay, l, t)); }
358
360 dcomplex muxx(size_t PLASK_UNUSED(lay), int l, int t) { return mag_long[(l>=0)?l:l+nNl].c11 * mag_tran[(t>=0)?t:t+nNt].c00; }
361
363 dcomplex muyy(size_t PLASK_UNUSED(lay), int l, int t) { return mag_long[(l>=0)?l:l+nNl].c00 * mag_tran[(t>=0)?t:t+nNt].c11; }
364
366 dcomplex imuzz(size_t PLASK_UNUSED(lay), int l, int t) { return mag_long[(l>=0)?l:l+nNl].c11 * mag_tran[(t>=0)?t:t+nNt].c11; }
367
369 size_t iEx(int l, int t) {
370 return 2 * idx(l, t);
371 }
372
374 size_t iEy(int l, int t) {
375 return 2 * idx(l, t) + 1;
376 }
377
379 size_t iHx(int l, int t) {
380 return 2 * idx(l, t) + 1;
381 }
382
384 size_t iHy(int l, int t) {
385 return 2 * idx(l, t);
386 }
387};
388
389 struct ExpansionPW3D::Gradient::Vertex {
390 int l, t;
391 Gradient val;
392 Vertex(int l, int t, const Gradient& src): l(l), t(t), val(src) {}
393};
394
395}}} // namespace plask
396
397#endif // PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_PW3D_H