PLaSK library

Loading...
Searching...
No Matches
impedance.cpp
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/*#if defined(_WIN32) || defined(__WIN32__) || defined(WIN32)
15# define BOOST_USE_WINDOWS_H
16#endif*/
17
18#include "impedance.hpp"
19#include "solver.hpp"
20#include "diagonalizer.hpp"
21#include "expansion.hpp"
22#include "fortran.hpp"
23#include "meshadapter.hpp"
24
25 namespace plask { namespace optical { namespace modal {
26
27 ImpedanceTransfer::ImpedanceTransfer(ModalBase* solver, Expansion& expansion): XanceTransfer(solver, expansion)
28{
29 writelog(LOG_DETAIL, "{}: Initializing Impedance Transfer", solver->getId());
30}
31
32
33 void ImpedanceTransfer::getFinalMatrix()
34{
35 int N = int(diagonalizer->matrixSize()); // for using with LAPACK it is better to have int instead of std::size_t
36 int N0 = int(diagonalizer->source()->matrixSize());
37 size_t count = solver->stack.size();
38
39 // M = TE(interface) * Y(interface-1) * invTH(interface);
40 findImpedance(count-1, solver->interface-1);
41 zgemm('n','n', N, N0, N, 1., Y.data(), N, diagonalizer->invTH(solver->stack[solver->interface]).data(), N, 0., wrk, N);
42 zgemm('n','n', N0, N0, N, 1., diagonalizer->TE(solver->stack[solver->interface]).data(), N0, wrk, N, 0., M.data(), N0);
43
44 // Find the(diagonalized field) admittance matrix and store it for the future reference
45 findImpedance(0, solver->interface);
46 // M += TE(interface-1) * Y(interface) * invTH(interface-1);
47 zgemm('n','n', N, N0, N, 1., Y.data(), N, diagonalizer->invTH(solver->stack[solver->interface-1]).data(), N, 0., wrk, N);
48 zgemm('n','n', N0, N0, N, 1., diagonalizer->TE(solver->stack[solver->interface-1]).data(), N0, wrk, N, 1., M.data(), N0);
49}
50
51
52 void ImpedanceTransfer::findImpedance(std::ptrdiff_t start, std::ptrdiff_t end)
53{
54 const std::ptrdiff_t inc = (start < end) ? 1 : -1;
55
56 const std::size_t N = diagonalizer->matrixSize();
57 const std::size_t NN = N*N;
58
59 // Some temporary variables
60 cdiagonal gamma, y1(N), y2(N);
61
62 std::exception_ptr error;
63
64 PLASK_OMP_PARALLEL_FOR_1
65 for (int l = 0; l < int(diagonalizer->lcount); ++l) {
66 try {
67 if (!error) diagonalizer->diagonalizeLayer(l);
68 } catch(...) {
69 error = std::current_exception();
70 }
71 }
72 if (error) std::rethrow_exception(error);
73
74 // Now iteratively we find matrices Y[i]
75
76 // PML layer
77 #ifdef OPENMP_FOUND
78 write_debug("{}: Entering into single region of admittance search", solver->getId());
79 #endif
80 gamma = diagonalizer->Gamma(solver->stack[start]);
81 std::fill_n(y2.data(), N, dcomplex(1.)); // we use y2 for tracking sign changes
82 for (std::size_t i = 0; i < N; i++) {
83 y1[i] = gamma[i] * solver->vpml.factor;
84 if (real(y1[i]) < -SMALL) { y1[i] = -y1[i]; y2[i] = -y2[i]; }
85 if (imag(y1[i]) > SMALL) { y1[i] = -y1[i]; y2[i] = -y2[i]; }
86 }
87 get_y1(y1, solver->vpml.size, y1);
88 std::fill_n(Y.data(), NN, dcomplex(0.));
89 for (std::size_t i = 0; i < N; i++) Y(i,i) = - y2[i] / y1[i];
90
91 // First layer
92 double h = solver->vpml.dist;
93 gamma = diagonalizer->Gamma(solver->stack[start]);
94 get_y1(gamma, h, y1);
95 get_y2(gamma, h, y2);
96 // off-diagonal elements of Y are 0
97 for (std::size_t i = 0; i < N; i++) Y(i,i) = y2[i] * y2[i] / (y1[i] - Y(i,i)) - y1[i]; // Y = y2 * inv(y1-Y) * y2 - y1
98
99 // save the Y matrix for 1-st layer
100 storeY(start);
101
102 if (start == end) return;
103
104 // Declare temporary matrixH) on 'wrk' array
105 cmatrix work(N, N, wrk);
106
107 for (std::ptrdiff_t n = start+inc; n != end; n += inc)
108 {
109 gamma = diagonalizer->Gamma(solver->stack[n]);
110
111 h = solver->vbounds->at(n) - solver->vbounds->at(n-1);
112 get_y1(gamma, h, y1);
113 get_y2(gamma, h, y2);
114
115 // The main equation
116 // Y[n] = y2 * tH * inv(y1*tH - tE*Y[n-1]) * y2 - y1
117
118 mult_matrix_by_matrix(diagonalizer->TE(solver->stack[n-inc]), Y, temp); // work = tE * Y[n-1]
119 mult_matrix_by_matrix(diagonalizer->invTE(solver->stack[n]), temp, work); // ...
120
121 mult_matrix_by_matrix(diagonalizer->invTH(solver->stack[n]), diagonalizer->TH(solver->stack[n-inc]), temp); // compute tH
122
123 for (std::size_t j = 0; j < N; j++)
124 for (std::size_t i = 0; i < N; i++) Y(i,j) = y1[i]*temp(i,j) - work(i,j); // Y[n] = y1 * tH - work
125
126 for (std::size_t i = 0; i < NN; i++) work[i] = 0.;
127 for (std::size_t j = 0, i = 0; j < N; j++, i += N+1) work[i] = y2[j]; // work = y2
128
129 invmult(Y, work); // work = inv(Y[n]) * (work = y2)
130 mult_matrix_by_matrix(temp, work, Y); // Y[n] = tH * work
131
132 for (std::size_t j = 0; j < N; j++)
133 for (std::size_t i = 0; i < N; i++) Y(i,j) *= y2[i]; // Y[n] = y2 * Y[n]
134
135 for (std::size_t j = 0, i = 0; j < N; j++, i += N+1) Y[i] -= y1[j]; // Y[n] = Y[n] - y1
136
137 // Save the Y matrix for n-th layer
138 storeY(n);
139 }
140}
141
142
143
144 void ImpedanceTransfer::determineFields()
145{
146 if (fields_determined == DETERMINED_RESONANT) return;
147
148 writelog(LOG_DETAIL, solver->getId() + ": Determining optical fields");
149
150 const std::size_t N = diagonalizer->matrixSize();
151 const std::size_t N0 = diagonalizer->source()->matrixSize();
152 size_t count = solver->stack.size();
153
154 const std::size_t NN = N*N;
155
156 // Assign all the required space
157 cdiagonal gamma, y1(N), y2(N);
158
159 // Assign the space for the field vectors
160 fields.resize(count);
161
162 // Temporary vector for storing fields in the real domain
163 cvector tv(N0);
164
165 // Obtain the physical fields at the last layer
166 needAllY = true;
167 interface_field = nullptr;
168 auto H = getInterfaceVector();
169
170 // Declare temporary matrix on 'wrk' array
171 cmatrix work(N, N, wrk);
172
173 for (int pass = 0; pass < 1 || (pass < 2 && solver->interface != std::ptrdiff_t(count)); pass++)
174 {
175 // each pass for below and above the interface
176
177 std::ptrdiff_t start, end;
178 int inc;
179 switch (pass) {
180 case 0: start = solver->interface-1; end = -1; inc = 1; break;
181 case 1: start = solver->interface; end = count; inc = -1; break;
182 }
183
184 // Hd[start] = invTH[start] H
185 fields[start].Hd = cvector(N);
186 mult_matrix_by_vector(diagonalizer->invTH(solver->stack[start]), H, fields[start].Hd);
187
188 fields[start].Hd *= double(inc);
189
190 for (std::ptrdiff_t n = start; n != end; n -= inc)
191 {
192 const std::size_t curr = solver->stack[n];
193
194 double h = (n == 0 || n == std::ptrdiff_t(count)-1)? solver->vpml.dist : solver->vbounds->at(n) - solver->vbounds->at(n-1);
195 gamma = diagonalizer->Gamma(curr);
196 get_y1(gamma, h, y1);
197 get_y2(gamma, h, y2);
198
199 // work = Y[n] + y1
200 cmatrix Y = getY(n);
201 for (std::size_t i = 0; i < NN; i++) work[i] = Y[i];
202 for (std::size_t i = 0; i < N; i++) work (i,i) += y1[i];
203
204 // H0[n] = work * Hd[n]
205 fields[n].H0 = cvector(N);
206 mult_matrix_by_vector(work, fields[n].Hd, fields[n].H0);
207
208 // H0[n] = - inv(y2) * H0[0]
209 for (size_t i = 0; i < N; i++) {
210 if (abs(y2[i]) < SMALL) // Actually we cannot really compute H0 in this case.
211 fields[n].H0[i] = 0.; // So let's cheat a little, as the field cannot
212 else // increase to the boundaries.
213 fields[n].H0[i] /= - y2[i];
214 }
215
216 if (n != end+inc) { // not the last layer
217 const std::size_t prev = solver->stack[n-inc];
218 // Hd[n-inc] = invTH[n-inc] * TH[n] * H0[n]
219 fields[n-inc].Hd = cvector(N);
220 mult_matrix_by_vector(diagonalizer->TH(curr), fields[n].H0, tv);
221 mult_matrix_by_vector(diagonalizer->invTH(prev), tv, fields[n-inc].Hd);
222 }
223
224 // Now compute the electric fields
225
226 // Ed[n] = Y[n] * Hd[n]
227 fields[n].Ed = cvector(N);
228 mult_matrix_by_vector(Y, fields[n].Hd, fields[n].Ed);
229
230 if (n != start) {
231 std::size_t next = solver->stack[n+inc];
232 // E0[n+inc] = invTE[n+inc] * TE[n] * Ed[n]
233 fields[n+inc].E0 = cvector(N);
234 mult_matrix_by_vector(diagonalizer->TE(curr), fields[n].Ed, tv);
235 mult_matrix_by_vector(diagonalizer->invTE(next), tv, fields[n+inc].E0);
236 }
237
238 // An alternative method is to find the E0 from the following equation:
239 // E0 = y1 * H0 + y2 * Hd
240 // for (int i = 0; i < N; i++)
241 // fields[n].E0[i] = y1[i] * fields[n].H0[i] + y2[i] * fields[n].Hd[i];
242 // However in some cases this can make the electric field discontinuous.
243 }
244 if (start != end) {
245 // Zero electric field at the end
246 std::ptrdiff_t n = end + inc;
247 fields[n].E0 = cvector(N, 0.);
248 }
249 }
250
251 // Now fill the Y matrix with the one from the interface (necessary for interfaceField*)
252 memcpy(Y.data(), getY(solver->interface-1).data(), NN*sizeof(dcomplex));
253
254 needAllY = false;
255 fields_determined = DETERMINED_RESONANT;
256
257 // Finally normalize fields
258 if (solver->emission == ModalBase::EMISSION_BOTTOM || solver->emission == ModalBase::EMISSION_TOP) {
259 const std::size_t n = (solver->emission == ModalBase::EMISSION_BOTTOM)? 0 : count-1;
260 const std::size_t l = solver->stack[n];
261
262 cvector hv(N0);
263 mult_matrix_by_vector(diagonalizer->TE(l), fields[n].Ed, tv);
264 mult_matrix_by_vector(diagonalizer->TH(l), fields[n].Hd, hv);
265
266 double P = 1./Z0 * abs(diagonalizer->source()->integratePoyntingVert(tv, hv));
267
268 if (P < SMALL) {
269 writelog(LOG_WARNING, "Device is not emitting to the {} side: skipping normalization",
270 (solver->emission == ModalBase::EMISSION_TOP)? "top" : "bottom");
271 } else {
272 P = 1. / sqrt(P);
273 for (size_t i = 0; i < count; ++i) {
274 fields[i].E0 *= P;
275 fields[i].H0 *= P;
276 fields[i].Ed *= P;
277 fields[i].Hd *= P;
278 }
279 }
280 }
281}
282
283
284 cvector ImpedanceTransfer::getReflectionVector(const cvector& incident, IncidentDirection side)
285{
286 throw NotImplemented("reflection with impedance transfer");
287}
288
289 void ImpedanceTransfer::determineReflectedFields(const cvector& incident, IncidentDirection side)
290{
291 throw NotImplemented("reflection with impedance transfer");
292}
293
294
295}}} // namespace plask::optical::modal