PLaSK library

Loading...
Searching...
No Matches
admittance.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 "admittance.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 AdmittanceTransfer::AdmittanceTransfer(ModalBase* solver, Expansion& expansion): XanceTransfer(solver, expansion)
28{
29 writelog(LOG_DETAIL, "{}: Initializing Admittance Transfer", solver->getId());
30}
31
32
33 void AdmittanceTransfer::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 = TH(interface) * Y(interface-1) * invTE(interface);
40 findAdmittance(count-1, solver->interface-1);
41 zgemm('n','n', N, N0, N, 1., Y.data(), N, diagonalizer->invTE(solver->stack[solver->interface]).data(), N, 0., wrk, N);
42 zgemm('n','n', N0, N0, N, 1., diagonalizer->TH(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 findAdmittance(0, solver->interface);
46 // M += TH(interface-1) * Y(interface) * invTE(interface-1);
47 zgemm('n','n', N, N0, N, 1., Y.data(), N, diagonalizer->invTE(solver->stack[solver->interface-1]).data(), N, 0., wrk, N);
48 zgemm('n','n', N0, N0, N, 1., diagonalizer->TH(solver->stack[solver->interface-1]).data(), N0, wrk, N, 1., M.data(), N0);
49}
50
51
52 void AdmittanceTransfer::findAdmittance(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) = - y1[i] * y2[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 * tE * inv(y1*tE - tH*Y[n-1]) * y2 - y1
117
118 mult_matrix_by_matrix(diagonalizer->TH(solver->stack[n-inc]), Y, temp); // work = tH * Y[n-1]
119 mult_matrix_by_matrix(diagonalizer->invTH(solver->stack[n]), temp, work); // ...
120
121 mult_matrix_by_matrix(diagonalizer->invTE(solver->stack[n]), diagonalizer->TE(solver->stack[n-inc]), temp); // compute tE
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 * tE - 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] = tE * 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 AdmittanceTransfer::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 E = 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 // Ed[start] = invTE[start] E
185 fields[start].Ed = cvector(N);
186 mult_matrix_by_vector(diagonalizer->invTE(solver->stack[start]), E, fields[start].Ed);
187
188 for (std::ptrdiff_t n = start; n != end; n -= inc)
189 {
190 const std::size_t curr = solver->stack[n];
191
192 double h = (n == 0 || n == std::ptrdiff_t(count)-1)? solver->vpml.dist : solver->vbounds->at(n) - solver->vbounds->at(n-1);
193 gamma = diagonalizer->Gamma(curr);
194 get_y1(gamma, h, y1);
195 get_y2(gamma, h, y2);
196
197 // work = Y[n] + y1
198 cmatrix Y = getY(n);
199 for (std::size_t i = 0; i < NN; i++) work[i] = Y[i];
200 for (std::size_t i = 0; i < N; i++) work (i,i) += y1[i];
201
202 // E0[n] = work * Ed[n]
203 fields[n].E0 = cvector(N);
204 mult_matrix_by_vector(work, fields[n].Ed, fields[n].E0);
205
206 // E0[n] = - inv(y2) * E0[0]
207 for (size_t i = 0; i < N; i++) {
208 if (abs(y2[i]) < SMALL) // Actually we cannot really compute E0 in this case.
209 fields[n].E0[i] = 0.; // So let's cheat a little, as the field cannot
210 else // increase to the boundaries.
211 fields[n].E0[i] /= - y2[i];
212 }
213
214 if (n != end+inc) { // not the last layer
215 const std::size_t prev = solver->stack[n-inc];
216 // Ed[n-inc] = invTE[n-inc] * TE[n] * E0[n]
217 fields[n-inc].Ed = cvector(N);
218 mult_matrix_by_vector(diagonalizer->TE(curr), fields[n].E0, tv);
219 mult_matrix_by_vector(diagonalizer->invTE(prev), tv, fields[n-inc].Ed);
220 }
221
222 // Now compute the magnetic fields
223
224 // Hd[n] = Y[n] * Ed[n]
225 fields[n].Hd = cvector(N);
226 mult_matrix_by_vector(Y, fields[n].Ed, fields[n].Hd);
227
228 if (n != start) {
229 std::size_t next = solver->stack[n+inc];
230 // H0[n+inc] = invTH[n+inc] * TH[n] * Hd[n]
231 fields[n+inc].H0 = cvector(N);
232 mult_matrix_by_vector(diagonalizer->TH(curr), fields[n].Hd, tv);
233 mult_matrix_by_vector(diagonalizer->invTH(next), tv, fields[n+inc].H0);
234 }
235
236 // An alternative method is to find the H0 from the following equation:
237 // H0 = y1 * E0 + y2 * Ed
238 // for (int i = 0; i < N; i++)
239 // fields[n].H0[i] = y1[i] * fields[n].E0[i] + y2[i] * fields[n].Ed[i];
240 // However in some cases this can make the magnetic field discontinuous.
241 }
242 if (start != end) {
243 // anyway, we must do it in the last layer
244 // (y1 and y2 are already computed in the above loop)
245 std::ptrdiff_t n = end + inc;
246 fields[n].H0 = cvector(N);
247 for (std::size_t i = 0; i < N; i++)
248 fields[n].H0[i] = y1[i] * fields[n].E0[i] + y2[i] * fields[n].Ed[i];
249 }
250 }
251
252 // Now fill the Y matrix with the one from the interface (necessary for interfaceField*)
253 memcpy(Y.data(), getY(solver->interface-1).data(), NN*sizeof(dcomplex));
254
255 needAllY = false;
256 fields_determined = DETERMINED_RESONANT;
257
258 // Finally normalize fields
259 if (solver->emission == ModalBase::EMISSION_BOTTOM || solver->emission == ModalBase::EMISSION_TOP) {
260 const std::size_t n = (solver->emission == ModalBase::EMISSION_BOTTOM)? 0 : count-1;
261 const std::size_t l = solver->stack[n];
262
263 cvector hv(N0);
264 mult_matrix_by_vector(diagonalizer->TE(l), fields[n].Ed, tv);
265 mult_matrix_by_vector(diagonalizer->TH(l), fields[n].Hd, hv);
266
267 double P = 1./Z0 * abs(diagonalizer->source()->integratePoyntingVert(tv, hv));
268
269 if (P < SMALL) {
270 writelog(LOG_WARNING, "Device is not emitting to the {} side: skipping normalization",
271 (solver->emission == ModalBase::EMISSION_TOP)? "top" : "bottom");
272 } else {
273 P = 1. / sqrt(P);
274 for (size_t i = 0; i < count; ++i) {
275 fields[i].E0 *= P;
276 fields[i].H0 *= P;
277 fields[i].Ed *= P;
278 fields[i].Hd *= P;
279 }
280 }
281 }
282}
283
284
285 cvector AdmittanceTransfer::getReflectionVector(const cvector& incident, IncidentDirection side)
286{
287 std::size_t curr, prev;
288
289 initDiagonalization();
290
291 switch (side) {
292 case INCIDENCE_TOP:
293 findAdmittance(0, solver->stack.size()-1);
294 curr = solver->stack[solver->stack.size()-1];
295 prev = solver->stack[solver->stack.size()-2];
296 break;
297 case INCIDENCE_BOTTOM:
298 findAdmittance(solver->stack.size()-1, 0);
299 curr = solver->stack[0];
300 prev = solver->stack[1];
301 break;
302 }
303
304 std::size_t N = diagonalizer->matrixSize();
305 std::size_t NN = N * N;
306 cmatrix work(N, N, wrk); // we have Y, temp and work
307
308 // Transfer to the outermost layer:
309 if (prev != curr) {
310 mult_matrix_by_matrix(diagonalizer->invTE(prev), diagonalizer->TE(curr), work); // work = tE¯¹
311 mult_matrix_by_matrix(Y, work, temp); // temp = Y tE¯¹
312 mult_matrix_by_matrix(diagonalizer->invTH(curr), diagonalizer->TH(prev), work); // work = tH
313 mult_matrix_by_matrix(work, temp, Y); // Y = tH Y tE¯¹
314 }
315
316 std::copy_n(Y.data(), NN, temp.data()); // temp = Y
317 // Use Jacobi preconditioner for temp
318 for (std::size_t i = 0; i != N; ++i) {
319 temp(i,i) -= 1.; // temp = Y - I
320 dcomplex f = 1. / temp(i,i);
321 wrk[i] = f;
322 for (size_t j = 0; j != N; ++j)
323 temp(i,j) *= f;
324 }
325 cvector reflected(N);
326 for (std::size_t i = 0; i != N; ++i) reflected[i] = wrk[i] * incident[i];
327 invmult(temp, reflected);
328 for (std::size_t i = 0; i != N; ++i)
329 reflected[i] = - 2. * reflected[i] - incident[i]; // R = - 2 [Y-I]¯¹ P - P
330 return reflected;
331}
332
333 void AdmittanceTransfer::determineReflectedFields(const cvector& incident, IncidentDirection side)
334{
335 if (fields_determined == DETERMINED_REFLECTED && incident == incident_vector && side == incident_side) return;
336 incident_vector = incident.copy();
337 incident_side = side;
338
339 writelog(LOG_DETAIL, solver->getId() + ": Determining reflected optical fields");
340
341 const std::size_t N = diagonalizer->matrixSize();
342 const std::size_t N0 = diagonalizer->source()->matrixSize();
343 size_t count = solver->stack.size();
344
345 const std::size_t NN = N*N;
346
347 // Assign all the required space
348 cdiagonal gamma, y1(N), iy2(N);
349
350 // Assign the space for the field vectors
351 fields.resize(count);
352
353 // Temporary vector for storing fields in the real domain
354 cvector tv(N0);
355
356 ptrdiff_t start, end, inc;
357 switch (side) {
358 case INCIDENCE_TOP: start = count-1; end = -1; inc = 1; break;
359 case INCIDENCE_BOTTOM: start = 0; end = count; inc = -1; break;
360 }
361
362 // Obtain the physical fields at the incident-side layer
363 needAllY = true;
364 fields[start].E0 = getReflectionVector(incident, side);
365 fields[start].E0 += incident;
366 fields[start].H0 = Y * fields[start].E0;
367
368 // E(z) = cosh(iΓz) A + sinh(iΓz) B
369 // H(z) = -sinh(iΓz) A - cosh(iΓz) B
370 // A = E(0) B = -H(0)
371
372 gamma = diagonalizer->Gamma(solver->stack[start]);
373 get_y1(gamma, solver->vpml.dist, y1);
374
375 for (std::size_t i = 0; i != N; ++i) Y(i,i) -= y1[i];
376 fields[start].Ed = Y * fields[start].E0;
377 fields[start].Hd = cvector(N);
378 for (std::size_t i = 0; i < N; i++) {
379 dcomplex s = - sinh(I*gamma[i]*solver->vpml.dist);
380 double as = abs(s);
381 if (isinf(real(s)) || isinf(imag(s)) || as > 1./SMALL) {
382 fields[start].Ed[i] = 0.;
383 fields[start].Hd[i] = 0.;
384 } else {
385 if (as < SMALL) {
386 writelog(LOG_WARNING, "{}: Cannot compute fields at the structure input side (try changing vpml.dist a bit)", solver->getId());
387 }
388 fields[start].Ed[i] *= s;
389 fields[start].Hd[i] = - fields[start].E0[i] / s - y1[i] * fields[start].Ed[i];
390 }
391 }
392
393 // Declare temporary matrix on 'wrk' array
394 cmatrix work(N, N, wrk);
395
396 mult_matrix_by_matrix(diagonalizer->invTE(solver->stack[start-inc]), diagonalizer->TE(solver->stack[start]), work);
397 fields[start-inc].Ed = work * fields[start].E0;
398
399 for (std::ptrdiff_t n = start-inc; n != end; n -= inc)
400 {
401 const std::size_t curr = solver->stack[n];
402
403 double h = (n == 0 || n == int(count)-1)? solver->vpml.dist : solver->vbounds->at(n) - solver->vbounds->at(n-1);
404 gamma = diagonalizer->Gamma(curr);
405 get_y1(gamma, h, y1);
406
407 // work = Y[n] + y1
408 cmatrix Y = getY(n);
409 std::copy_n(Y.data(), NN, work.data());
410 for (std::size_t i = 0; i < N; i++) work (i,i) += y1[i];
411
412 // E0[n] = work * Ed[n]
413 fields[n].E0 = cvector(N);
414 mult_matrix_by_vector(work, fields[n].Ed, fields[n].E0);
415
416 // E0[n] = - inv(y2) * E0[0]
417 for (std::size_t i = 0; i < N; i++) {
418 iy2[i] = sinh(I*gamma[i]*h);
419 if (isinf(real(iy2[i])) || isinf(imag(iy2[i])) || abs(iy2[i]) > 1./SMALL) {
420 fields[n].E0[i] = 0.; // Actually we cannot really compute E0 in this case.
421 } else { // So let's cheat a little, as the field cannot
422 fields[n].E0[i] *= iy2[i]; // increase to the boundaries.
423 }
424 }
425
426 if (n != end+inc) { // not the last layer
427 const std::size_t prev = solver->stack[n-inc];
428 // Ed[n-inc] = invTE[n-inc] * TE[n] * E0[n]
429 mult_matrix_by_vector(diagonalizer->TE(curr), fields[n].E0, tv);
430 fields[n-inc].Ed = diagonalizer->invTE(prev) * tv;
431 } else {
432 fields[n].H0 = cvector(N);
433 for (std::size_t i = 0; i < N; i++)
434 //fields[end+inc].H0[i] = y2[i] * fields[end+inc].Ed[i];
435 fields[end+inc].H0[i] = double(inc) *
436 (y1[i] * fields[end+inc].E0[i] - fields[end+inc].Ed[i]) / iy2[i];
437 }
438
439 // Now compute the magnetic fields
440
441 // Hd[n] = Y[n] * Ed[n]
442 fields[n].Hd = Y * fields[n].Ed;
443
444 if (std::ptrdiff_t(n) != start-inc) {
445 std::size_t next = solver->stack[n+inc];
446 // H0[n+inc] = invTH[n+inc] * TH[n] * Hd[n]
447 fields[n+inc].H0 = cvector(N);
448 mult_matrix_by_vector(diagonalizer->TH(curr), fields[n].Hd, tv);
449 mult_matrix_by_vector(diagonalizer->invTH(next), tv, fields[n+inc].H0);
450 }
451
452 // An alternative method is to find the H0 from the following equation:
453 // H0 = y1 * E0 + y2 * Ed
454 // for (int i = 0; i < N; i++)
455 // fields[n].H0[i] = y1[i] * fields[n].E0[i] + y2[i] * fields[n].Ed[i];
456 // However in some cases this can make the magnetic field discontinuous
457 }
458
459 // Replace F and B at one side of the interface for consistency in getFieldVectorE and getFieldVectorH
460 size_t interface = size_t(max(solver->interface, ptrdiff_t(0)));
461 switch (side)
462 {
463 case INCIDENCE_TOP: start = interface; end = count; break;
464 case INCIDENCE_BOTTOM: start = 0; end = min(interface, count); break;
465 }
466 // start = size_t(max(solver->interface, ptrdiff_t(0))); end = count;
467 for (ptrdiff_t n = start; n < end; ++n) {
468 std::swap(fields[n].E0, fields[n].Ed);
469 std::swap(fields[n].H0, fields[n].Hd);
470 // TODO should I revert H?
471 }
472
473 needAllY = false;
474 fields_determined = DETERMINED_REFLECTED;
475}
476
477
478}}} // namespace plask::optical::modal