PLaSK library

Loading...
Searching...
No Matches
expansion2d.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_PW2D_H
15#define PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_PW2D_H
16
17#include <plask/plask.hpp>
18
19#include "../expansion.hpp"
20#include "../meshadapter.hpp"
21#include "fft.hpp"
22
23 namespace plask { namespace optical { namespace modal {
24
25 struct FourierSolver2D;
26
27 struct PLASK_SOLVER_API ExpansionPW2D : public Expansion {
28 dcomplex beta,
29 ktran;
30
31 size_t N;
32 size_t nN;
33 double left;
34 double right;
35 bool periodic;
36 bool initialized;
37
38 Component symmetry;
39 Component polarization;
40
41 size_t pil,
42 pir;
43
44 struct Coeffs {
45 DataVector<dcomplex> zz, rxx, yy, zx;
46 };
48 std::vector<Coeffs> coeffs;
49
50 struct CoeffMatrices {
51 cmatrix exx, reyy;
52 };
54 std::vector<CoeffMatrices> coeff_matrices;
55
57 cmatrix coeff_matrix_mxx, coeff_matrix_rmyy;
58
60 std::vector<bool> diagonals;
61
63 shared_ptr<RectangularMesh<2>> mesh;
64
66 shared_ptr<MeshAxis> original_mesh;
67
72 ExpansionPW2D(FourierSolver2D* solver);
73
75 bool symmetric() const { return symmetry != E_UNSPECIFIED; }
76
78 bool separated() const { return polarization != E_UNSPECIFIED; }
79
84 void init();
85
87 void reset();
88
89 bool diagonalQE(size_t l) const override { return diagonals[l]; }
90
91 size_t matrixSize() const override { return separated() ? N : 2 * N; }
92
93 void getMatrices(size_t l, cmatrix& RE, cmatrix& RH) override;
94
95 void prepareField() override;
96
97 void cleanupField() override;
98
99 LazyData<Vec<3, dcomplex>> getField(size_t l,
100 const shared_ptr<const typename LevelsAdapter::Level>& level,
101 const cvector& E,
102 const cvector& H) override;
103
104 LazyData<Tensor3<dcomplex>> getMaterialEps(size_t l,
105 const shared_ptr<const typename LevelsAdapter::Level>& level,
106 InterpolationMethod interp = INTERPOLATION_DEFAULT) override;
107
108 double integrateField(WhichField field,
109 size_t layer,
110 const cmatrix& TE,
111 const cmatrix& TH,
112 const std::function<std::pair<dcomplex, dcomplex>(size_t, size_t)>& vertical) override;
113
114 double integratePoyntingVert(const cvector& E, const cvector& H) override;
115
116 void getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal& gamma) override;
117
118 private:
119 DataVector<Vec<3, dcomplex>> field;
120 FFT::Backward1D fft_x, fft_yz;
121
122 void add_coeffs(int start, int end, double b, double l, double r, DataVector<dcomplex>& dst, dcomplex val) {
123 for (int k = start; k != end; ++k) {
124 size_t j = (k >= 0) ? k : k + nN;
125 dcomplex ff = (j) ? (dcomplex(0., 0.5 / PI / k) * (exp(dcomplex(0., -b * k * r)) - exp(dcomplex(0., -b * k * l))))
126 : ((r - l) * b * (0.5 / PI));
127 dst[j] += val * ff;
128 }
129 }
130
131 void make_permeability_matrices(cmatrix& work);
132
133 protected:
134 DataVector<dcomplex> mag;
135 DataVector<dcomplex> rmag;
136
137 FFT::Forward1D matFFT;
138
139 void beforeLayersIntegrals(dcomplex lam, dcomplex glam) override;
140
141 void layerIntegrals(size_t layer, double lam, double glam) override;
142
143 Tensor3<dcomplex> getEpsilon(const shared_ptr<GeometryD<2>>& geometry,
144 size_t layer,
145 double maty,
146 double lam,
147 double glam,
148 size_t j) {
149 Tensor3<dcomplex> eps;
150 std::set<std::string> roles;
151 if (epsilon_connected && solver->lcomputed[layer] || gain_connected && solver->lgained[layer])
152 roles = geometry->getRolesAt(vec(mesh->tran()->at(j), maty));
153 bool computed = epsilon_connected && solver->lcomputed[layer] && roles.find("inEpsilon") != roles.end();
154 if (computed) {
155 eps = Zero<Tensor3<dcomplex>>();
156 double W = 0.;
157 for (size_t k = 0, v = j * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
158 if (solver->stack[k] == layer) {
159 if (isnan(epsilons[v]))
160 throw BadInput(solver->getId(), "complex permittivity tensor got from inEpsilon is NaN at {}", mesh->at(v));
161 double w =
162 (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
163 eps += w * epsilons[v];
164 W += w;
165 }
166 }
167 eps /= W;
168 } else {
169 double T = 0., W = 0., C = 0.;
170 for (size_t k = 0, v = j * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
171 if (solver->stack[k] == layer) {
172 double w =
173 (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
174 T += w * temperature[v];
175 C += w * carriers[v];
176 W += w;
177 }
178 }
179 T /= W;
180 C /= W;
181 {
182 OmpLockGuard lock; // this must be declared before `material` to guard its destruction
183 auto material = geometry->getMaterial(vec(mesh->tran()->at(j), maty));
184 lock = material->lock();
185 eps = material->Eps(lam, T, C);
186 if (isnan(eps))
187 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) for {} is NaN at lam={}nm, T={}K, n={}/cm3",
188 material->name(), lam, T, C);
189 }
190 }
191 if (eps.c01 != 0. || eps.c10 != 0.) {
192 if (symmetric())
193 throw BadInput(solver->getId(), "symmetry not allowed for structure with non-diagonal permittivity tensor");
194 if (separated())
195 throw BadInput(solver->getId(),
196 "single polarization not allowed for structure with non-diagonal permittivity tensor");
197 }
198 if (eps.c02 != 0. || eps.c21 != 0. || eps.c12 != 0. || eps.c20 != 0.) {
199 throw BadInput(solver->getId(), "symmetry not allowed for structure with non-diagonal permittivity tensor");
200 }
201 if (!computed && gain_connected && solver->lgained[layer] && roles.find("QW") != roles.end() ||
202 roles.find("QD") != roles.end() || roles.find("gain") != roles.end()) {
203 Tensor2<double> g = 0.;
204 double W = 0.;
205 for (size_t k = 0, v = j * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
206 if (solver->stack[k] == layer) {
207 double w =
208 (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
209 g += w * gain[v];
210 W += w;
211 }
212 }
213 Tensor2<double> ni = glam * g / W * (0.25e-7 / PI);
214 double n00 = sqrt(eps.c00).real(), n11 = sqrt(eps.c11).real(), n22 = sqrt(eps.c22).real();
215 eps.c00 = dcomplex(n00 * n00 - ni.c00 * ni.c00, 2 * n00 * ni.c00);
216 eps.c11 = dcomplex(n11 * n11 - ni.c00 * ni.c00, 2 * n11 * ni.c00);
217 eps.c22 = dcomplex(n22 * n22 - ni.c11 * ni.c11, 2 * n22 * ni.c11);
218 eps.c01.imag(0.);
219 eps.c02.imag(0.);
220 eps.c10.imag(0.);
221 eps.c12.imag(0.);
222 eps.c20.imag(0.);
223 eps.c21.imag(0.);
224 }
225 return eps;
226 }
227
228 public:
229 dcomplex getBeta() const { return beta; }
230 void setBeta(dcomplex b) {
231 if (b != beta) {
232 beta = b;
233 solver->clearFields();
234 }
235 }
236
237 dcomplex getKtran() const { return ktran; }
238 void setKtran(dcomplex k) {
239 if (k != ktran) {
240 ktran = k;
241 solver->clearFields();
242 }
243 }
244
245 Component getSymmetry() const { return symmetry; }
246 void setSymmetry(Component sym) {
247 if (sym != symmetry) {
248 symmetry = sym;
249 solver->clearFields();
250 solver->recompute_integrals = true;
251 }
252 }
253
254 Component getPolarization() const { return polarization; }
255 void setPolarization(Component pol);
256
257 const DataVector<dcomplex>& epszz(size_t l) { return coeffs[l].zz; }
258 const DataVector<dcomplex>& epsyy(size_t l) { return coeffs[l].yy; }
259 const DataVector<dcomplex>& repsxx(size_t l) { return coeffs[l].rxx; }
260 const DataVector<dcomplex>& epszx(size_t l) { return coeffs[l].zx; }
261 const DataVector<dcomplex>& muzz() { return mag; }
262 const DataVector<dcomplex>& muxx() { return mag; }
263 const DataVector<dcomplex>& muyy() { return mag; }
264 const DataVector<dcomplex>& rmuzz() { return rmag; }
265 const DataVector<dcomplex>& rmuxx() { return rmag; }
266 const DataVector<dcomplex>& rmuyy() { return rmag; }
267
268 dcomplex epszz(size_t l, int i) { return coeffs[l].zz[(i >= 0) ? i : i + nN]; }
269 dcomplex epsyy(size_t l, int i) { return coeffs[l].yy[(i >= 0) ? i : i + nN]; }
270 dcomplex repsxx(size_t l, int i) {
271 return coeffs[l].rxx[(i >= 0) ? i : i + nN];
272 }
273 dcomplex epszx(size_t l, int i) { return coeffs[l].zx[(i >= 0) ? i : i + nN]; }
274 dcomplex epsxz(size_t l, int i) {
275 return conj(coeffs[l].zx[(i >= 0) ? i : i + nN]);
276 }
277 dcomplex muzz(int i) { return mag[(i >= 0) ? i : i + nN]; }
278 dcomplex muxx(int i) { return mag[(i >= 0) ? i : i + nN]; }
279 dcomplex muyy(int i) { return mag[(i >= 0) ? i : i + nN]; }
280 dcomplex rmuzz(int i) { return rmag[(i >= 0) ? i : i + nN]; }
281 dcomplex rmuxx(int i) { return rmag[(i >= 0) ? i : i + nN]; }
282
283 size_t iEx(int i) { return 2 * ((i >= 0) ? i : i + N); }
284 size_t iEz(int i) { return 2 * ((i >= 0) ? i : i + N) + 1; }
285 size_t iHx(int i) { return 2 * ((i >= 0) ? i : i + N) + 1; }
286 size_t iHz(int i) { return 2 * ((i >= 0) ? i : i + N); }
287 size_t iEH(int i) { return (i >= 0) ? i : i + N; }
288};
289
290}}} // namespace plask::optical::modal
291
292#endif // PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_PW2D_H