PLaSK library

Loading...
Searching...
No Matches
solver3d.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_FOURIER_REFLECTION_3D_H
15#define PLASK__SOLVER_OPTICAL_MODAL_FOURIER_REFLECTION_3D_H
16
17#include <plask/plask.hpp>
18
19#include "../solver.hpp"
20#include "../reflection.hpp"
21#include "expansion3d.hpp"
22
23#ifdef minor
24# undef minor
25#endif
26
27 namespace plask { namespace optical { namespace modal {
28
32 struct PLASK_SOLVER_API FourierSolver3D: public ModalSolver<SolverOver<Geometry3D>> {
33
34 typedef ModalSolver<SolverOver<Geometry3D>> BaseType;
35
36 friend struct ExpansionPW3D;
37
38 std::string getClassName() const override { return "optical.Fourier3D"; }
39
41 enum What {
42 WHAT_WAVELENGTH,
43 WHAT_K0,
44 WHAT_KLONG,
45 WHAT_KTRAN
46 };
47
49 enum ExpansionRule {
50 RULE_OLD,
51 RULE_DIRECT,
52 RULE_INVERSE,
53 RULE_COMBINED
54 };
55
56 struct Mode {
57 Expansion::Component symmetry_long;
58 Expansion::Component symmetry_tran;
59 double lam0;
60 dcomplex k0;
61 dcomplex klong;
62 dcomplex ktran;
63 double power;
64 double tolx;
65
66 Mode(const ExpansionPW3D& expansion, double tolx):
67 symmetry_long(expansion.symmetry_long),
68 symmetry_tran(expansion.symmetry_tran),
69 lam0(expansion.lam0),
70 k0(expansion.k0),
71 klong(expansion.klong),
72 ktran(expansion.ktran),
73 power(1.),
74 tolx(tolx) {}
75
76 bool operator==(const Mode& other) const {
77 return is_equal(k0, other.k0) && is_equal(klong, other.klong) && is_equal(ktran, other.ktran)
78 && symmetry_long == other.symmetry_long && symmetry_tran == other.symmetry_tran &&
79 ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
80 ;
81 }
82
83 bool operator==(const ExpansionPW3D& other) const {
84 return is_equal(k0, other.k0) && is_equal(klong, other.klong) && is_equal(ktran, other.ktran)
85 && symmetry_long == other.symmetry_long && symmetry_tran == other.symmetry_tran &&
86 ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
87 ;
88 }
89
90 template <typename T>
91 bool operator!=(const T& other) const {
92 return !(*this == other);
93 }
94
95 private:
96
98 template <typename T>
99 bool is_equal(T a, T b) const {
100 return abs(a-b) <= tolx;
101 }
102 };
103
105 size_t size_long;
107 size_t size_tran;
108
109 protected:
110
111 dcomplex klong,
112 ktran;
113
114 Expansion::Component symmetry_long,
115 symmetry_tran;
116
117 void onInitialize() override;
118
119 void onInvalidate() override;
120
121 void computeIntegrals() override {
122 expansion.computeIntegrals();
123 }
124
126 int dct;
127
129 ExpansionRule expansion_rule;
130
131 size_t initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam);
132
133 public:
134
136 ExpansionPW3D expansion;
137
139 std::vector<Mode> modes;
140
141 void clearModes() override {
142 modes.clear();
143 }
144
145 bool setExpansionDefaults(bool with_k0=true) override {
146 bool changed = false;
147 if (expansion.getLam0() != getLam0()) { changed = true; expansion.setLam0(getLam0()); }
148 if (with_k0) {
149 if (expansion.getK0() != getK0()) { changed = true; expansion.setK0(getK0()); }
150 }
151 if (expansion.getKlong() != getKlong()) { changed = true; expansion.setKlong(getKlong()); }
152 if (expansion.getKtran() != getKtran()) { changed = true; expansion.setKtran(getKtran()); }
153 if (expansion.getSymmetryLong() != getSymmetryLong()) { changed = true; expansion.setSymmetryLong(getSymmetryLong()); }
154 if (expansion.getSymmetryTran() != getSymmetryTran()) { changed = true; expansion.setSymmetryTran(getSymmetryTran()); }
155 return changed;
156 }
157
159 size_t refine_long;
161 size_t refine_tran;
162
164 double grad_smooth;
165
167 PML pml_long;
169 PML pml_tran;
170
172 ProviderFor<GradientFunctions, Geometry3D>::Delegate outGradients;
173
174 FourierSolver3D(const std::string& name="");
175
176 void loadConfiguration(XMLReader& reader, Manager& manager) override;
177
185 size_t findMode(What what, dcomplex start);
186
188 size_t getLongSize() const { return size_long; }
189
191 size_t getTranSize() const { return size_tran; }
192
194 void setLongSize(size_t n) {
195 size_long = n;
196 invalidate();
197 }
199 void setTranSize(size_t n) {
200 size_tran = n;
201 invalidate();
202 }
203
205 void setSizes(size_t nl, size_t nt) {
206 size_long = nl;
207 size_tran = nt;
208 invalidate();
209 }
210
212 Expansion::Component getSymmetryLong() const { return symmetry_long; }
213
215 void setSymmetryLong(Expansion::Component symmetry) {
216 if (symmetry != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry3D::DIRECTION_LONG))
217 throw BadInput(getId(), "longitudinal symmetry not allowed for asymmetric structure");
218 if ((symmetry_long == Expansion::E_UNSPECIFIED) != (symmetry == Expansion::E_UNSPECIFIED))
219 invalidate();
220 if (klong != 0. && symmetry != Expansion::E_UNSPECIFIED) {
221 Solver::writelog(LOG_WARNING, "Resetting klong to 0.");
222 klong = 0.;
223 expansion.setKlong(0.);
224 }
225 symmetry_long = symmetry;
226 }
227
229 Expansion::Component getSymmetryTran() const { return symmetry_tran; }
230
232 void setSymmetryTran(Expansion::Component symmetry) {
233 if (symmetry != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry3D::DIRECTION_TRAN))
234 throw BadInput(getId(), "transverse symmetry not allowed for asymmetric structure");
235 if ((symmetry_tran == Expansion::E_UNSPECIFIED) != (symmetry == Expansion::E_UNSPECIFIED))
236 invalidate();
237 if (ktran != 0. && symmetry != Expansion::E_UNSPECIFIED) {
238 Solver::writelog(LOG_WARNING, "Resetting ktran to 0.");
239 ktran = 0.;
240 expansion.setKtran(0.);
241 }
242 symmetry_tran = symmetry;
243 }
244
246 bool symmetricLong() const { return expansion.symmetric_long(); }
247
249 bool symmetricTran() const { return expansion.symmetric_tran(); }
250
252 dcomplex getKlong() const { return klong; }
253
255 void setKlong(dcomplex k) {
256 if (k != 0. && (expansion.symmetric_long() || symmetry_long != Expansion::E_UNSPECIFIED)) {
257 Solver::writelog(LOG_WARNING, "Resetting longitudinal mode symmetry");
258 symmetry_long = Expansion::E_UNSPECIFIED;
259 invalidate();
260 }
261 klong = k;
262 }
263
265 dcomplex getKtran() const { return ktran; }
266
268 void setKtran(dcomplex k) {
269 if (k != 0. && (expansion.symmetric_tran() || symmetry_tran != Expansion::E_UNSPECIFIED)) {
270 Solver::writelog(LOG_WARNING, "Resetting transverse mode symmetry");
271 symmetry_tran = Expansion::E_UNSPECIFIED;
272 invalidate();
273 }
274 ktran = k;
275 }
276
278 double getGradSmooth() const { return grad_smooth; }
279
281 void setGradSmooth(double value) {
282 bool changed = grad_smooth != value;
283 grad_smooth = value;
284 if (changed) this->invalidate();
285 }
286
288 int getDCT() const { return dct; }
290 void setDCT(int n) {
291 if (n != 1 && n != 2)
292 throw BadInput(getId(), "bad DCT type (can be only 1 or 2)");
293 if (dct != n) {
294 dct = n;
295 if (expansion.symmetric_long() || expansion.symmetric_tran()) invalidate();
296 }
297 }
299 bool dct2() const { return dct == 2; }
300
302 ExpansionRule getRule() const { return expansion_rule; }
304 void setRule(ExpansionRule rule) {
305 if (rule != expansion_rule) {
306 expansion_rule = rule;
307 invalidate();
308 }
309 }
310
311 // /// Get mesh at which material parameters are sampled along longitudinal axis
312 // RegularAxis getLongMesh() const { return expansion.mesh->lon(); }
313 //
314 // /// Get mesh at which material parameters are sampled along transverse axis
315 // RegularAxis getTranMesh() const { return expansion.mesh->tran(); }
316
317 Expansion& getExpansion() override { return expansion; }
318
320 size_t minor() const { return expansion.Nl; }
321
329 cvector incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
330
331 using BaseType::incidentVector;
332
342 cvector incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization,
343 double sigma_long, double sigma_tran, double center_long=0., double center_tran=0., dcomplex lam=NAN);
344
353 LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
354 Transfer::IncidentDirection side,
355 const shared_ptr<const MeshD<3>>& dst_mesh,
356 InterpolationMethod method,
357 PropagationDirection part = PROPAGATION_TOTAL) {
358 if (!Solver::initCalculation()) setExpansionDefaults(false);
359 if (!transfer) initTransfer(expansion, true);
360 return transfer->getScatteredFieldE(incident, side, dst_mesh, method, part);
361 }
362
371 LazyData<Vec<3,dcomplex>> getScatteredFieldH(const cvector& incident,
372 Transfer::IncidentDirection side,
373 const shared_ptr<const MeshD<3>>& dst_mesh,
374 InterpolationMethod method,
375 PropagationDirection part = PROPAGATION_TOTAL) {
376 if (!Solver::initCalculation()) setExpansionDefaults(false);
377 if (!transfer) initTransfer(expansion, true);
378 return transfer->getScatteredFieldH(incident, side, dst_mesh, method, part);
379 }
380
388 LazyData<double> getScatteredFieldMagnitude(const cvector& incident,
389 Transfer::IncidentDirection side,
390 const shared_ptr<const MeshD<3>>& dst_mesh,
391 InterpolationMethod method) {
392 if (!Solver::initCalculation()) setExpansionDefaults(false);
393 if (!transfer) initTransfer(expansion, true);
394 return transfer->getScatteredFieldMagnitude(incident, side, dst_mesh, method);
395 }
396
403 cvector getFieldVectorE(size_t num, double z) {
404 applyMode(modes[num]);
405 return transfer->getFieldVectorE(z);
406 }
407
414 cvector getFieldVectorH(size_t num, double z) {
415 applyMode(modes[num]);
416 return transfer->getFieldVectorH(z);
417 }
418
426 double getIntegralEE(size_t num, double z1, double z2) {
427 applyMode(modes[num]);
428 return transfer->getFieldIntegral(FIELD_E, z1, z2, modes[num].power);
429 }
430
438 double getIntegralHH(size_t num, double z1, double z2) {
439 applyMode(modes[num]);
440 return transfer->getFieldIntegral(FIELD_H, z1, z2, modes[num].power);
441 }
442
450 cvector getScatteredFieldVectorE(const cvector& incident, Transfer::IncidentDirection side, double z) {
451 if (!Solver::initCalculation()) setExpansionDefaults(false);
452 if (!transfer) initTransfer(expansion, true);
453 return transfer->getScatteredFieldVectorE(incident, side, z, PROPAGATION_TOTAL);
454 }
455
463 cvector getScatteredFieldVectorH(const cvector& incident, Transfer::IncidentDirection side, double z) {
464 if (!Solver::initCalculation()) setExpansionDefaults(false);
465 if (!transfer) initTransfer(expansion, true);
466 return transfer->getScatteredFieldVectorH(incident, side, z, PROPAGATION_TOTAL);
467 }
468
477 double getScatteredIntegralEE(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
478 if (!Solver::initCalculation()) setExpansionDefaults(false);
479 if (!transfer) initTransfer(expansion, true);
480 return transfer->getScatteredFieldIntegral(FIELD_E, incident, side, z1, z2);
481 }
482
491 double getScatteredIntegralHH(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
492 if (!Solver::initCalculation()) setExpansionDefaults(false);
493 if (!transfer) initTransfer(expansion, true);
494 return transfer->getScatteredFieldIntegral(FIELD_H, incident, side, z1, z2);
495 }
496
498 size_t setMode() {
499 if (abs2(this->getDeterminant()) > root.tolf_max*root.tolf_max)
500 throw BadInput(this->getId(), "cannot set the mode, determinant too large");
501 return insertMode();
502 }
503
504 protected:
505
507 size_t insertMode() {
508 static bool warn = true;
509 if (warn && emission != EMISSION_TOP && emission != EMISSION_BOTTOM) {
510 writelog(LOG_WARNING, "Mode fields are not normalized unless emission is set to 'top' or 'bottom'");
511 warn = false;
512 }
513 Mode mode(expansion, root.tolx);
514 for (size_t i = 0; i != modes.size(); ++i)
515 if (modes[i] == mode) return i;
516 modes.push_back(mode);
517 outLightMagnitude.fireChanged();
518 outLightE.fireChanged();
519 outLightH.fireChanged();
520 return modes.size()-1;
521 }
522
523 size_t nummodes() const override { return modes.size(); }
524
525 double applyMode(size_t n) override {
526 if (n >= modes.size()) throw BadInput(this->getId(), "mode {0} has not been computed", n);
527 applyMode(modes[n]);
528 return modes[n].power;
529 }
530
535 dcomplex getEffectiveIndex(size_t n) {
536 if (n >= modes.size()) throw NoValue(ModeEffectiveIndex::NAME);
537 return modes[n].klong / modes[n].k0;
538 }
539
540 void applyMode(const Mode& mode) {
541 writelog(LOG_DEBUG, "Current mode <lam: {}nm, klong: {}/um, ktran: {}/um, symmetry: ({},{})>",
542 str(2e3*PI/mode.k0, "({:.3f}{:+.3g}j)", "{:.3f}"),
543 str(mode.klong, "({:.3f}{:+.3g}j)", "{:.3f}"),
544 str(mode.ktran, "({:.3f}{:+.3g}j)", "{:.3f}"),
545 (mode.symmetry_long == Expansion::E_LONG)? "El" : (mode.symmetry_long == Expansion::E_TRAN)? "Et" : "none",
546 (mode.symmetry_tran == Expansion::E_LONG)? "El" : (mode.symmetry_tran == Expansion::E_TRAN)? "Et" : "none"
547 );
548 if (mode != expansion) {
549 expansion.setLam0(mode.lam0);
550 expansion.setK0(mode.k0);
551 expansion.klong = mode.klong;
552 expansion.ktran = mode.ktran;
553 expansion.symmetry_long = mode.symmetry_long;
554 expansion.symmetry_tran = mode.symmetry_tran;
555 clearFields();
556 }
557 }
558
559 double getWavelength(size_t n) override;
560
561 LazyData<double> getGradients(GradientFunctions::EnumType what,
562 const shared_ptr<const MeshD<3>>& dst_mesh,
563 InterpolationMethod interp);
564
565};
566
567
568}}} // namespace
569
570#endif