PLaSK library

Loading...
Searching...
No Matches
transfer.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#include "transfer.hpp"
15#include "diagonalizer.hpp"
16#include "expansion.hpp"
17#include "fortran.hpp"
18#include "solver.hpp"
19
20 namespace plask { namespace optical { namespace modal {
21
22 Transfer::Transfer(ModalBase* solver, Expansion& expansion)
23 : solver(solver),
24 diagonalizer(new SimpleDiagonalizer(&expansion)), // TODO add other diagonalizer types
25 fields_determined(DETERMINED_NOTHING) {
26 // Reserve space for matrix multiplications...
27 const std::size_t N0 = diagonalizer->source()->matrixSize();
28 const std::size_t N = diagonalizer->matrixSize();
29 M = cmatrix(N0, N0);
30 temp = cmatrix(N, N);
31
32 // ...and eigenvalues determination
33 evals = aligned_new_array<dcomplex>(N0);
34 rwrk = aligned_new_array<double>(2 * N0);
35 lwrk = max(std::size_t(2 * N0), N0 * N0);
36 wrk = aligned_new_array<dcomplex>(lwrk);
37
38 // Nothing found so far
39 fields_determined = DETERMINED_NOTHING;
40 interface_field = nullptr;
41}
42
43 Transfer::~Transfer() {
44 // no need for aligned_delete_array because array members have trivial destructor
45 aligned_free<dcomplex>(evals);
46 evals = nullptr;
47 aligned_free<double>(rwrk);
48 rwrk = nullptr;
49 aligned_free<dcomplex>(wrk);
50 wrk = nullptr;
51}
52
53 void Transfer::initDiagonalization() {
54 // Get new coefficients if needed
55 solver->computeIntegrals();
56 this->diagonalizer->initDiagonalization();
57}
58
59dcomplex Transfer::determinant() {
60 // We change the matrices M and A so we will have to find the new fields
61 fields_determined = DETERMINED_NOTHING;
62
63 initDiagonalization();
64
65 // Obtain admittance
66 getFinalMatrix();
67
68 const std::size_t N = M.rows();
69
70 // This is probably expensive but necessary check to avoid hangs
71 const std::size_t NN = N * N;
72 for (std::size_t i = 0; i < NN; i++) {
73 if (isnan(real(M[i])) || isnan(imag(M[i]))) throw ComputationError(solver->getId(), "NaN in discontinuity matrix");
74 }
75
76 // Find the eigenvalues of M using LAPACK
77 dcomplex result;
78
79 switch (solver->determinant_type) {
80 case DETERMINANT_EIGENVALUE: {
81 // TODO add some consideration for degenerate modes
82 int info;
83 zgeev('N', 'N', int(N), M.data(), int(N), evals, nullptr, 1, nullptr, 1, wrk, int(lwrk), rwrk, info);
84 if (info != 0) throw ComputationError(solver->getId(), "eigenvalue determination failed");
85 // Find the smallest eigenvalue
86 double min_mag = 1e32;
87 for (std::size_t i = 0; i < N; i++) {
88 dcomplex val = evals[i];
89 double mag = abs2(val);
90 if (mag < min_mag) {
91 min_mag = mag;
92 result = val;
93 }
94 }
95 } break;
96
97 case DETERMINANT_FULL:
98 result = det(M);
99 break;
100
101 default:
102 assert(false);
103 }
104
105 interface_field = nullptr;
106
107 return result;
108}
109
110 const_cvector Transfer::getInterfaceVector() {
111 const std::size_t N = M.rows();
112
113 // Check if the necessary memory is already allocated
114 if (interface_field_matrix.rows() != N) {
115 interface_field_matrix = cmatrix(N, N);
116 interface_field = nullptr;
117 }
118
119 // If the field already found, don't compute again
120 if (!interface_field) {
121 // We change the matrices M and A so we will have to find the new fields
122 fields_determined = DETERMINED_NOTHING;
123
124 initDiagonalization();
125
126 // Obtain admittance
127 getFinalMatrix();
128
129 // Find the eigenvalues of M using LAPACK
130 int info;
131 zgeev('N', 'V', int(N), M.data(), int(N), evals, nullptr, 1, interface_field_matrix.data(), int(N), wrk, int(lwrk), rwrk,
132 info);
133 if (info != 0) throw ComputationError(solver->getId(), "interface field: zgeev failed");
134
135 // Find the number of the smallest eigenvalue
136 double min_mag = 1e32;
137 std::size_t n;
138 for (std::size_t i = 0; i < N; i++) {
139 double mag = abs2(evals[i]);
140 if (mag < min_mag) {
141 min_mag = mag;
142 n = i;
143 }
144 }
145
146 // Error handling
147 if (min_mag > solver->root.tolf_max * solver->root.tolf_max)
148 throw BadInput(solver->getId(), "interface field: determinant not sufficiently close to 0 (det={})", str(evals[n]));
149
150 // Chose the eigenvector corresponding to the smallest eigenvalue
151 interface_field = interface_field_matrix.data() + n * N;
152 }
153
154 return DataVector<const dcomplex>(interface_field, N);
155}
156
157 LazyData<Vec<3, dcomplex>> Transfer::computeFieldE(double power,
158 const shared_ptr<const Mesh>& dst_mesh,
159 InterpolationMethod method,
160 bool reflected,
161 PropagationDirection part) {
162 double fact = sqrt(2e-3 * power);
163 double zlim = solver->vpml.dist + solver->vpml.size;
164 DataVector<Vec<3, dcomplex>> destination(dst_mesh->size());
165 auto levels = makeLevelsAdapter(dst_mesh);
166 diagonalizer->source()->initField(FIELD_E, method);
167 while (auto level = levels->yield()) {
168 double z = level->vpos();
169 const std::size_t n = solver->getLayerFor(z);
170 if (!reflected) {
171 if (n == 0 && z < -zlim)
172 z = -zlim;
173 else if (n == solver->stack.size() - 1 && z > zlim)
174 z = zlim;
175 }
176 cvector E = getFieldVectorE(z, n, part);
177 cvector H = getFieldVectorH(z, n, part);
178 if (std::ptrdiff_t(n) >= solver->interface)
179 for (auto& h : H) h = -h;
180 size_t layer = solver->stack[n];
181 auto dest = fact * diagonalizer->source()->getField(layer, level, E, H);
182 for (size_t i = 0; i != level->size(); ++i) destination[level->index(i)] = dest[i];
183 }
184 diagonalizer->source()->cleanupField();
185 return destination;
186}
187
188 LazyData<Vec<3, dcomplex>> Transfer::computeFieldH(double power,
189 const shared_ptr<const Mesh>& dst_mesh,
190 InterpolationMethod method,
191 bool reflected,
192 PropagationDirection part) {
193 double fact = 1. / Z0 * sqrt(2e-3 * power);
194 double zlim = solver->vpml.dist + solver->vpml.size;
195 DataVector<Vec<3, dcomplex>> destination(dst_mesh->size());
196 auto levels = makeLevelsAdapter(dst_mesh);
197 diagonalizer->source()->initField(FIELD_H, method);
198 while (auto level = levels->yield()) {
199 double z = level->vpos();
200 size_t n = solver->getLayerFor(z);
201 if (!reflected) {
202 if (n == 0 && z < -zlim)
203 z = -zlim;
204 else if (n == solver->stack.size() - 1 && z > zlim)
205 z = zlim;
206 }
207 cvector E = getFieldVectorE(z, n, part);
208 cvector H = getFieldVectorH(z, n, part);
209 if (std::ptrdiff_t(n) >= solver->interface)
210 for (auto& h : H) h = -h;
211 size_t layer = solver->stack[n];
212 auto dest = fact * diagonalizer->source()->getField(layer, level, E, H);
213 for (size_t i = 0; i != level->size(); ++i) destination[level->index(i)] = dest[i];
214 }
215 diagonalizer->source()->cleanupField();
216 return destination;
217}
218
219 cvector Transfer::getFieldVectorE(double z, PropagationDirection part) {
220 determineFields();
221 const std::size_t n = solver->getLayerFor(z);
222 return getFieldVectorE(z, n, part);
223}
224
225 cvector Transfer::getFieldVectorH(double z, PropagationDirection part) {
226 determineFields();
227 const std::size_t n = solver->getLayerFor(z);
228 cvector H = getFieldVectorH(z, n, part);
229 if (std::ptrdiff_t(n) >= solver->interface)
230 for (auto& h : H) h = -h;
231 return H;
232}
233
234 cvector Transfer::getScatteredFieldVectorE(const cvector& incident, IncidentDirection side, double z, PropagationDirection part) {
235 determineReflectedFields(incident, side);
236 return getFieldVectorE(z, solver->getLayerFor(z), part);
237}
238
239 cvector Transfer::getScatteredFieldVectorH(const cvector& incident, IncidentDirection side, double z, PropagationDirection part) {
240 determineReflectedFields(incident, side);
241 return getFieldVectorH(z, solver->getLayerFor(z), part);
242}
243
244 double Transfer::getFieldIntegral(WhichField field, double z1, double z2, double power) {
245 determineFields();
246 if (z1 > z2) std::swap(z1, z2);
247 size_t end = solver->getLayerFor(z2);
248 if (is_zero(z2) && end != 0) {
249 --end;
250 z2 = solver->vbounds->at(end) - (end ? solver->vbounds->at(end-1) : solver->vbounds->at(end));
251 }
252 double result = 0.;
253 for (size_t n = solver->getLayerFor(z1); n <= end; ++n) {
254 result += integrateField(field, n, z1, (n != end) ? n ? solver->vbounds->at(n) - solver->vbounds->at(n-1) : 0. : z2);
255 z1 = 0.;
256 }
257 return ((field == FIELD_E) ? 2e-3 : 2e-3 / (Z0 * Z0)) * power * result;
258}
259
260 double Transfer::getScatteredFieldIntegral(WhichField field,
261 const cvector& incident,
262 IncidentDirection side,
263 double z1,
264 double z2) {
265 determineReflectedFields(incident, side);
266 if (z1 > z2) std::swap(z1, z2);
267 size_t end = solver->getLayerFor(z2);
268 if (is_zero(z2) && end != 0) {
269 --end;
270 z2 = solver->vbounds->at(end) - (end ? solver->vbounds->at(end-1) : solver->vbounds->at(end));
271 }
272 double result = 0.;
273 for (size_t n = solver->getLayerFor(z1); n <= end; ++n) {
274 result += integrateField(field, n, z1, (n != end) ? n ? solver->vbounds->at(n) - solver->vbounds->at(n-1) : 0. : z2);
275 z1 = 0.;
276 }
277 return ((field == FIELD_E) ? 2. * Z0 : 2. / Z0) * result;
278}
279
280}}} // namespace plask::optical::modal