PLaSK library
Loading...
Searching...
No Matches
expansioncyl-fini.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 "expansioncyl-fini.hpp"
15#include "solvercyl.hpp"
16#include "zeros-data.hpp"
17
18#include "../gauss_legendre.hpp"
19
20#include "besselj.hpp"
21
22#define SOLVER static_cast<BesselSolverCyl*>(solver)
23
24namespace plask { namespace optical { namespace modal {
25
26ExpansionBesselFini::ExpansionBesselFini(BesselSolverCyl* solver) : ExpansionBessel(solver) {}
27
28void ExpansionBesselFini::computeBesselZeros() {
29 size_t N = SOLVER->size;
30 size_t n = 0;
31 kpts.resize(N);
32 if (m < 5) {
33 n = min(N, size_t(100));
34 std::copy_n(bessel_zeros[m], n, kpts.begin());
35 }
36 if (n < N) {
37 SOLVER->writelog(LOG_DEBUG, "Computing Bessel function J_({:d}) zeros {:d} to {:d}", m, n + 1, N);
38 PLASK_NO_CONVERSION_WARNING_BEGIN
39 cyl_bessel_j_zero(double(m), n+1, N-n, kpts.begin()+n);
40 PLASK_NO_WARNING_END
41 }
42 // #ifndef NDEBUG
43 // for (size_t i = 0; i != N; ++i) {
44 // auto report = [i,m,this]()->bool{
45 // std::cerr << "J(" << m << ", " << kpts[i] << ") = " << cyl_bessel_j(m, kpts[i]) << "\n";
46 // return false;
47 // };
48 // assert(is_zero(cyl_bessel_j(m, kpts[i]), 1e-9) || report());
49 // }
50 // #endif
51}
52
53void ExpansionBesselFini::init2() {
54 SOLVER->writelog(LOG_DETAIL, "Preparing Bessel functions for m = {}", m);
55 computeBesselZeros();
56
57 init3();
58
59 // Compute integrals for permeability
60 auto raxis = mesh->tran();
61
62 size_t nr = raxis->size(), N = SOLVER->size;
63
64 double ib = 1. / rbounds[rbounds.size() - 1];
65
66 // Compute mu integrals
67 if (SOLVER->pml.size > 0. && SOLVER->pml.factor != 1.) {
68
69 SOLVER->writelog(LOG_DETAIL, "Computing permeability integrals with {} rule", SOLVER->ruleName());
70
71 size_t pmlseg = segments.size() - 1;
72
73 size_t pmli = raxis->size() - segments[pmlseg].weights.size();
74 double pmlr = rbounds[pmlseg];
75
76 aligned_unique_ptr<dcomplex> mu_data(aligned_malloc<dcomplex>(nr));
77 aligned_unique_ptr<dcomplex> imu_data(aligned_malloc<dcomplex>(nr));
78
79 for (size_t ri = 0, wi = 0, seg = 0, nw = segments[0].weights.size(); ri != nr; ++ri, ++wi) {
80 if (wi == nw) {
81 nw = segments[++seg].weights.size();
82 wi = 0;
83 }
84 double r = raxis->at(ri);
85 double w = segments[seg].weights[wi] * segments[seg].D;
86
87 dcomplex mu = 1., imu = 1.;
88 if (ri >= pmli) {
89 mu = 1. + (SOLVER->pml.factor - 1.) * pow((r - pmlr) / SOLVER->pml.size, SOLVER->pml.order);
90 imu = 1. / mu;
91 }
92
93 mu_data.get()[ri] = mu * w;
94 imu_data.get()[ri] = imu * w;
95 }
96
97 switch (SOLVER->rule) {
98 case BesselSolverCyl::RULE_COMBINED_1:
99 case BesselSolverCyl::RULE_COMBINED_2:
100 integrateParams(mu_integrals, mu_data.get(), mu_data.get(), mu_data.get(), 1., 1., 1.); break;
101 case BesselSolverCyl::RULE_DIRECT:
102 integrateParams(mu_integrals, mu_data.get(), imu_data.get(), mu_data.get(), 1., 1., 1.); break;
103 case BesselSolverCyl::RULE_OLD:
104 integrateParams(mu_integrals, mu_data.get(), imu_data.get(), imu_data.get(), 1., 1., 1.); break;
105 }
106
107 } else {
108
109 mu_integrals.reset(N);
110 zero_matrix(mu_integrals.V_k);
111 zero_matrix(mu_integrals.Tss);
112 zero_matrix(mu_integrals.Tsp);
113 zero_matrix(mu_integrals.Tps);
114 zero_matrix(mu_integrals.Tpp);
115 for (size_t i = 0; i < N; ++i) {
116 double k = kpts[i] * ib;
117 mu_integrals.V_k(i,i) = k;
118 mu_integrals.Tss(i,i) = mu_integrals.Tpp(i,i) = 2.;
119 }
120 }
121}
122
127
128void ExpansionBesselFini::getMatrices(size_t layer, cmatrix& RE, cmatrix& RH) {
129 assert(initialized);
130 if (isnan(k0)) throw BadInput(SOLVER->getId(), "wavelength or k0 not set");
131 if (isinf(k0.real())) throw BadInput(SOLVER->getId(), "wavelength must not be 0");
132
133 size_t N = SOLVER->size;
134 dcomplex ik0 = 1. / k0;
135 double ib = 1. / rbounds[rbounds.size() - 1];
136
137 const Integrals& eps = layers_integrals[layer];
138 #define mu mu_integrals
139
140 for (size_t j = 0; j != N; ++j) {
141 size_t js = idxs(j), jp = idxp(j);
142 // double k = kpts[j] * ib;
143 for (size_t i = 0; i != N; ++i) {
144 size_t is = idxs(i), ip = idxp(i);
145 double g = kpts[i] * ib;
146 dcomplex gVk = ik0 * g * eps.V_k(i,j);
147 RH(is,jp) = 0.5 * ( gVk - k0 * mu.Tsp(i,j) );
148 RH(is,js) = 0.5 * ( gVk - k0 * mu.Tss(i,j) );
149 RH(ip,jp) = 0.5 * ( -gVk + k0 * mu.Tpp(i,j) );
150 RH(ip,js) = 0.5 * ( -gVk + k0 * mu.Tps(i,j) );
151 }
152 }
153
154 for (size_t j = 0; j != N; ++j) {
155 size_t js = idxs(j), jp = idxp(j);
156 // double k = kpts[j] * ib;
157 for (size_t i = 0; i != N; ++i) {
158 size_t is = idxs(i), ip = idxp(i);
159 double g = kpts[i] * ib;
160 dcomplex gVk = ik0 * g * mu.V_k(i,j);
161 RE(ip,js) = 0.5 * ( -gVk + k0 * eps.Tps(i,j) );
162 RE(ip,jp) = 0.5 * ( -gVk + k0 * eps.Tpp(i,j) );
163 RE(is,js) = 0.5 * ( gVk - k0 * eps.Tss(i,j) );
164 RE(is,jp) = 0.5 * ( gVk - k0 * eps.Tsp(i,j) );
165 }
166 }
167 #undef mu
168}
169
170double ExpansionBesselFini::fieldFactor(size_t i) {
171 double eta = cyl_bessel_j(m + 1, kpts[i]) * rbounds[rbounds.size() - 1];
172 return 0.5 * eta * eta;
173}
174
175
176#ifndef NDEBUG
177cmatrix ExpansionBesselFini::muV_k() {
178 size_t N = SOLVER->size;
179 cmatrix result(N, N, 0.);
180 for (size_t i = 0; i != N; ++i)
181 for (size_t j = 0; j != N; ++j)
182 result(i,j) = mu_integrals.V_k(i,j);
183 return result;
184}
185cmatrix ExpansionBesselFini::muTss() {
186 size_t N = SOLVER->size;
187 cmatrix result(N, N, 0.);
188 for (size_t i = 0; i != N; ++i)
189 for (size_t j = 0; j != N; ++j)
190 result(i,j) = mu_integrals.Tss(i,j);
191 return result;
192}
193cmatrix ExpansionBesselFini::muTsp() {
194 size_t N = SOLVER->size;
195 cmatrix result(N, N, 0.);
196 for (size_t i = 0; i != N; ++i)
197 for (size_t j = 0; j != N; ++j)
198 result(i,j) = mu_integrals.Tsp(i,j);
199 return result;
200}
201cmatrix ExpansionBesselFini::muTps() {
202 size_t N = SOLVER->size;
203 cmatrix result(N, N, 0.);
204 for (size_t i = 0; i != N; ++i)
205 for (size_t j = 0; j != N; ++j)
206 result(i,j) = mu_integrals.Tps(i,j);
207 return result;
208}
209cmatrix ExpansionBesselFini::muTpp() {
210 size_t N = SOLVER->size;
211 cmatrix result(N, N, 0.);
212 for (size_t i = 0; i != N; ++i)
213 for (size_t j = 0; j != N; ++j)
214 result(i,j) = mu_integrals.Tpp(i,j);
215 return result;
216}
217#endif
218
219void ExpansionBesselFini::integrateParams(Integrals& integrals,
220 const dcomplex* datap, const dcomplex* datar, const dcomplex* dataz,
221 dcomplex datap0, dcomplex datar0, dcomplex dataz0) {
222 auto raxis = mesh->tran();
223
224 size_t nr = raxis->size(), N = SOLVER->size;
225 double R = rbounds[rbounds.size()-1];
226 double ib = 1. / R;
227 double* factors; // scale factors for making matrices orthonormal
228
229 integrals.reset(N);
230
231 TempMatrix temp = getTempMatrix();
232 aligned_unique_ptr<double> _tmp;
233
234 cmatrix JmJp, JpJm;
235 if (N < 2) {
236 _tmp.reset(aligned_malloc<double>(4*N));
237 factors = _tmp.get();
238 } else if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_1) {
239 factors = reinterpret_cast<double*>(integrals.V_k.data());
240 } else if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_2) {
241 _tmp.reset(aligned_malloc<double>(3*N + 4*N*N));
242 JmJp.reset(N, N, integrals.V_k.data());
243 JpJm.reset(N, N, reinterpret_cast<dcomplex*>(_tmp.get() + 3*N));
244 factors = _tmp.get();
245 } else {
246 factors = reinterpret_cast<double*>(temp.data());
247 }
248 for (size_t i = 0; i < N; ++i) {
249 double fact = R * cyl_bessel_j(m+1, kpts[i]);
250 factors[i] = 2. / (fact * fact);
251 }
252 double* Jm = factors + N;
253 double* Jp = factors + 2*N;
254
255 if (SOLVER->rule == BesselSolverCyl::RULE_OLD) {
256
257 double* J = factors + 3*N;
258
259 zero_matrix(integrals.V_k);
260 zero_matrix(integrals.Tss);
261 zero_matrix(integrals.Tsp);
262 zero_matrix(integrals.Tps);
263 zero_matrix(integrals.Tpp);
264
265 for (size_t ri = 0; ri != nr; ++ri) {
266 double r = raxis->at(ri);
267 dcomplex repst = r * (datar[ri] + datap[ri]), repsd = r * (datar[ri] - datap[ri]);
268 const dcomplex riepsz = r * dataz[ri];
269
270 for (size_t i = 0; i < N; ++i) {
271 double kr = kpts[i] * ib * r;
272 Jm[i] = cyl_bessel_j(m-1, kr);
273 J[i] = cyl_bessel_j(m, kr);
274 Jp[i] = cyl_bessel_j(m+1, kr);
275 }
276 for (size_t j = 0; j < N; ++j) {
277 double k = kpts[j] * ib;
278 double Jk = J[j], Jmk = Jm[j], Jpk = Jp[j];
279 for (size_t i = 0; i < N; ++i) {
280 double Jg = J[i], Jmg = Jm[i], Jpg = Jp[i];
281 integrals.V_k(i,j) += Jg * riepsz * Jk * k * factors[i];
282 integrals.Tss(i,j) += Jmg * repst * Jmk * factors[i];
283 integrals.Tsp(i,j) += Jmg * repsd * Jpk * factors[i];
284 integrals.Tps(i,j) += Jpg * repsd * Jmk * factors[i];
285 integrals.Tpp(i,j) += Jpg * repst * Jpk * factors[i];
286 }
287 }
288 }
289
290 } else {
291
292 if (SOLVER->rule == BesselSolverCyl::RULE_DIRECT) {
293
294 zero_matrix(integrals.Tss);
295 zero_matrix(integrals.Tsp);
296 zero_matrix(integrals.Tps);
297 zero_matrix(integrals.Tpp);
298
299 for (size_t ri = 0; ri != nr; ++ri) {
300 double r = raxis->at(ri);
301 dcomplex repst = r * (datar[ri] + datap[ri]), repsd = r * (datar[ri] - datap[ri]);
302 for (size_t i = 0; i < N; ++i) {
303 double kr = kpts[i] * ib * r;
304 Jm[i] = cyl_bessel_j(m-1, kr);
305 Jp[i] = cyl_bessel_j(m+1, kr);
306 }
307 for (size_t i = 0; i < N; ++i) {
308 dcomplex cs = factors[i] * repst, cd = factors[i] * repsd;
309 for (size_t j = 0; j < N; ++j) {
310 integrals.Tss(i,j) += cs * Jm[i] * Jm[j];
311 integrals.Tsp(i,j) += cd * Jm[i] * Jp[j];
312 integrals.Tps(i,j) += cd * Jp[i] * Jm[j];
313 integrals.Tpp(i,j) += cs * Jp[i] * Jp[j];
314 }
315 }
316 }
317
318 } else {
319
320 if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_1) {
321
322 cmatrix workess(N, N, temp.data()), workepp(N, N, temp.data()+N*N),
323 worksp(N, N, temp.data()+2*N*N), workps(N, N, temp.data()+3*N*N);
324
325 zero_matrix(workess);
326 zero_matrix(workepp);
327 zero_matrix(worksp);
328 zero_matrix(workps);
329
330 for (size_t ri = 0, wi = 0, seg = 0, nw = segments[0].weights.size(); ri != nr; ++ri, ++wi) {
331 if (wi == nw) {
332 nw = segments[++seg].weights.size();
333 wi = 0;
334 }
335 double r = raxis->at(ri);
336 double rw = r * segments[seg].weights[wi] * segments[seg].D;
337 dcomplex riepsr = r * datar[ri];
338 for (size_t i = 0; i < N; ++i) {
339 double kr = kpts[i] * ib * r;
340 Jm[i] = cyl_bessel_j(m-1, kr);
341 Jp[i] = cyl_bessel_j(m+1, kr);
342 }
343 for (size_t i = 0; i < N; ++i) {
344 double cw = factors[i] * rw;
345 dcomplex ce = factors[i] * riepsr;
346 for (size_t j = 0; j < N; ++j) {
347 workess(i,j) += ce * Jm[i] * Jm[j];
348 workepp(i,j) += ce * Jp[i] * Jp[j];
349 worksp(i,j) += cw * Jm[i] * Jp[j];
350 workps(i,j) += cw * Jp[i] * Jm[j];
351 }
352 }
353 }
354 make_unit_matrix(integrals.Tss);
355 make_unit_matrix(integrals.Tpp);
356
357 invmult(workess, integrals.Tss);
358 invmult(workepp, integrals.Tpp);
359
360 mult_matrix_by_matrix(integrals.Tss, worksp, integrals.Tsp);
361 mult_matrix_by_matrix(integrals.Tpp, workps, integrals.Tps);
362
363 if (N > 2) {
364 std::copy_n(factors, N, reinterpret_cast<double*>(temp.data()));
365 factors = reinterpret_cast<double*>(temp.data());
366 Jm = factors + N; Jp = factors + 2*N;
367 }
368
369 } else { // if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_2)
370
371 cmatrix work(temp);
372
373 zero_matrix(integrals.TT);
374 zero_matrix(work);
375 zero_matrix(JmJp);
376 zero_matrix(JpJm);
377
378 for (size_t ri = 0, wi = 0, seg = 0, nw = segments[0].weights.size(); ri != nr; ++ri, ++wi) {
379 if (wi == nw) {
380 nw = segments[++seg].weights.size();
381 wi = 0;
382 }
383 double r = raxis->at(ri);
384 double rw = r * segments[seg].weights[wi] * segments[seg].D;
385 dcomplex riepsr = r * datar[ri];
386 for (size_t i = 0; i < N; ++i) {
387 double kr = kpts[i] * ib * r;
388 Jm[i] = cyl_bessel_j(m-1, kr);
389 Jp[i] = cyl_bessel_j(m+1, kr);
390 }
391 for (size_t j = 0; j < N; ++j) {
392 for (size_t i = 0; i < N; ++i) {
393 dcomplex ce = riepsr * factors[i];
394 double cw = rw * factors[i];
395 integrals.TT(i,j) += ce * Jm[i] * Jm[j];
396 integrals.TT(i,j+N) += ce * Jm[i] * Jp[j];
397 integrals.TT(i+N,j) += ce * Jp[i] * Jm[j];
398 integrals.TT(i+N,j+N) += ce * Jp[i] * Jp[j];
399 double mp = cw * Jm[i] * Jp[j];
400 double pm = cw * Jp[i] * Jm[j];
401 work(i,j+N) += mp;
402 work(i+N,j) += pm;
403 JmJp(i,j) += mp;
404 JpJm(i,j) += pm;
405 }
406 }
407 }
408 for (size_t i = 0; i < N; ++i) work(i,i) = 1.;
409 for (size_t i = 0; i < N; ++i) work(i+N,i+N) = 1.;
410
411 invmult(integrals.TT, work);
412
413 for (size_t j = 0; j < N; ++j) {
414 for (size_t i = 0; i < N; ++i) integrals.Tss(i,j) = work(i,j);
415 for (size_t i = 0; i < N; ++i) integrals.Tps(i,j) = work(i+N,j);
416 }
417 for (size_t j = 0; j < N; ++j) {
418 for (size_t i = 0; i < N; ++i) integrals.Tsp(i,j) = work(i,j+N);
419 for (size_t i = 0; i < N; ++i) integrals.Tpp(i,j) = work(i+N,j+N);
420 }
421
422 zgemm('N', 'N', int(N), int(N), int(N), 1., JpJm.data(), int(N), work.data()+2*N*N, int(2*N), 1.,
423 integrals.Tpp.data(), int(N));
424 zgemm('N', 'N', int(N), int(N), int(N), 1., JmJp.data(), int(N), work.data()+N, int(2*N), 1.,
425 integrals.Tss.data(), int(N));
426 zgemm('N', 'N', int(N), int(N), int(N), 1., JpJm.data(), int(N), work.data(), int(2*N), 1.,
427 integrals.Tps.data(), int(N));
428 zgemm('N', 'N', int(N), int(N), int(N), 1., JmJp.data(), int(N), work.data()+2*N*N+N, int(2*N), 1.,
429 integrals.Tsp.data(), int(N));
430 }
431
432 for (size_t ri = 0; ri != nr; ++ri) {
433 double r = raxis->at(ri);
434 dcomplex repsp = r * datap[ri];
435 for (size_t i = 0; i < N; ++i) {
436 double kr = kpts[i] * ib * r;
437 Jm[i] = cyl_bessel_j(m-1, kr);
438 Jp[i] = cyl_bessel_j(m+1, kr);
439 }
440 for (size_t i = 0; i < N; ++i) {
441 dcomplex c = repsp * factors[i];
442 for (size_t j = 0; j < N; ++j) {
443 integrals.Tss(i,j) += c * Jm[i] * Jm[j];
444 integrals.Tsp(i,j) -= c * Jm[i] * Jp[j];
445 integrals.Tps(i,j) -= c * Jp[i] * Jm[j];
446 integrals.Tpp(i,j) += c * Jp[i] * Jp[j];
447 }
448 }
449 }
450 }
451
452 cmatrix work(N, N, temp.data()+N*N);
453
454 zero_matrix(work);
455 double* J = Jm;
456
457 for (size_t ri = 0; ri != nr; ++ri) {
458 double r = raxis->at(ri);
459 dcomplex repsz = r * dataz[ri];
460 for (size_t i = 0; i < N; ++i) {
461 double kr = kpts[i] * ib * r;
462 J[i] = cyl_bessel_j(m, kr);
463 }
464 for (size_t j = 0; j < N; ++j) {
465 for (size_t i = 0; i < N; ++i) {
466 work(i,j) += repsz * factors[i] * J[i] * J[j];
467 }
468 }
469 }
470
471 // make_unit_matrix(integrals.V_k);
472 zero_matrix(integrals.V_k);
473 for (int i = 0; i < N; ++i) {
474 double k = kpts[i] * ib;
475 integrals.V_k(i,i) = k;
476 }
477 invmult(work, integrals.V_k);
478 // for (size_t i = 0; i < N; ++i) {
479 // double g = kpts[i] * ib;
480 // for (size_t j = 0; j < N; ++j) integrals.V_k(i,j) *= g;
481 // }
482
483 }
484}
485
486}}} // namespace plask::optical::modal