PLaSK library

Loading...
Searching...
No Matches
solver3d.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 "solver3d.hpp"
15#include "expansion3d.hpp"
16#include "../diagonalizer.hpp"
17
18 namespace plask { namespace optical { namespace modal {
19
20 FourierSolver3D::FourierSolver3D(const std::string& name): ModalSolver<SolverOver<Geometry3D>>(name),
21 size_long(12), size_tran(12),
22 klong(0.), ktran(0.),
23 symmetry_long(Expansion::E_UNSPECIFIED), symmetry_tran(Expansion::E_UNSPECIFIED),
24 dct(2),
25 expansion_rule(RULE_DIRECT),
26 expansion(this),
27 refine_long(16), refine_tran(16),
28 grad_smooth(1e-3),
29 outGradients(this, &FourierSolver3D::getGradients)
30{
31 pml_tran.factor = {1., -2.};
32 pml_long.factor = {1., -2.};
33}
34
35 static inline void updatePML(PML& pml, XMLReader& reader) {
36 pml.factor = reader.getAttribute<dcomplex>("factor", pml.factor);
37 pml.size = reader.getAttribute<double>("size", pml.size);
38 pml.dist = reader.getAttribute<double>("dist", pml.dist);
39 pml.order = reader.getAttribute<double>("shape", pml.order);
40}
41
42 template <typename T>
43 static inline void readComaAttr(XMLReader& reader, const std::string& attr, T& long_field, T& tran_field, bool nocommon=false) {
44 if (reader.hasAttribute(attr)) {
45 std::string value = reader.requireAttribute<std::string>(attr);
46 if (value.find(',') == std::string::npos) {
47 if (nocommon) throw XMLBadAttrException(reader, attr, value);
48 long_field = tran_field = boost::lexical_cast<T>(value);
49 } else {
50 auto values = splitString2(value, ',');
51 long_field = boost::lexical_cast<T>(values.first);
52 tran_field = boost::lexical_cast<T>(values.second);
53 }
54 if (reader.hasAttribute(attr+"-long")) throw XMLConflictingAttributesException(reader, attr, attr+"-long");
55 if (reader.hasAttribute(attr+"-tran")) throw XMLConflictingAttributesException(reader, attr, attr+"-tran");
56 } else {
57 long_field = reader.getAttribute<T>(attr+"-long", long_field);
58 tran_field = reader.getAttribute<T>(attr+"-tran", tran_field);
59 }
60}
61
62 static inline Expansion::Component readSymmetry(const FourierSolver3D* solver, const XMLReader& reader, const std::string& repr) {
63 AxisNames* axes = nullptr;
64 if (solver->getGeometry()) axes = &solver->getGeometry()->axisNames;
65 if (repr == "none" || repr == "NONE" || repr == "None")
66 return Expansion::E_UNSPECIFIED;
67 else if (repr == "Etran" || repr == "Et" || (axes && repr == "E"+axes->getNameForTran()) ||
68 repr == "Hlong" || repr == "Hl" || (axes && repr == "H"+axes->getNameForLong()))
69 return Expansion::E_TRAN;
70 else if (repr == "Elong" || repr == "El" || (axes && repr == "E"+axes->getNameForLong()) ||
71 repr == "Htran" || repr == "Ht" || (axes && repr == "H"+axes->getNameForTran()))
72 return Expansion::E_LONG;
73 else
74 throw XMLBadAttrException(reader, "symmetry", repr, "symmetric field component name (maybe you need to specify the geometry first)");
75}
76
77 void FourierSolver3D::loadConfiguration(XMLReader& reader, Manager& manager)
78{
79 while (reader.requireTagOrEnd()) {
80 std::string param = reader.getNodeName();
81 if (param == "expansion") {
82 readComaAttr(reader, "size", size_long, size_tran);
83 readComaAttr(reader, "refine", refine_long, refine_tran);
84 if (reader.hasAttribute("oversampling")) {
85 reader.ignoreAttribute("oversampling");
86 writelog(LOG_WARNING, "obsolete 'oversampling' attribute in XML line {}", reader.getLineNr());
87 }
88 if (reader.hasAttribute("oversampling-long")) {
89 reader.ignoreAttribute("oversampling-long");
90 writelog(LOG_WARNING, "obsolete 'oversampling-long' attribute in XML line {}", reader.getLineNr());
91 }
92 if (reader.hasAttribute("oversampling-tran")) {
93 reader.ignoreAttribute("oversampling-tran");
94 writelog(LOG_WARNING, "obsolete 'oversampling-tran' attribute in XML line {}", reader.getLineNr());
95 }
96 smooth = reader.getAttribute<double>("smooth", smooth);
97 int dc = reader.getAttribute<int>("dct", dct);
98 if (dc != 1 && dc != 2)
99 throw XMLBadAttrException(reader, "dct", boost::lexical_cast<std::string>(dc), "\"1\" or \"2\"");
100 dct = dc;
101 group_layers = reader.getAttribute<bool>("group-layers", group_layers);
102 lam0 = reader.getAttribute<double>("lam0", NAN);
103 always_recompute_gain = reader.getAttribute<bool>("update-gain", always_recompute_gain);
104 max_temp_diff = reader.getAttribute<double>("temp-diff", max_temp_diff);
105 temp_dist = reader.getAttribute<double>("temp-dist", temp_dist);
106 temp_layer = reader.getAttribute<double>("temp-layer", temp_layer);
107 expansion_rule = reader.enumAttribute<ExpansionRule>("rule")
108 .value("direct", RULE_DIRECT)
109 .value("inverse", RULE_INVERSE)
110 .value("combined", RULE_COMBINED)
111 .value("old", RULE_OLD)
112 .get(expansion_rule);
113 grad_smooth = reader.getAttribute<double>("grad-smooth", grad_smooth);
114 reader.requireTagEnd();
115 } else if (param == "pmls") {
116 updatePML(pml_long, reader);
117 updatePML(pml_tran, reader);
118 while (reader.requireTagOrEnd()) {
119 std::string node = reader.getNodeName();
120 if (node == "long") {
121 updatePML(pml_long, reader);
122 reader.requireTagEnd();
123 } else if (node == "tran") {
124 updatePML(pml_tran, reader);
125 reader.requireTagEnd();
126 } else throw XMLUnexpectedElementException(reader, "<tran>, <long>, or </pmls>", node);
127 }
128 } else if (param == "mode") {
129 emission = reader.enumAttribute<Emission>("emission")
130 .value("undefined", EMISSION_UNSPECIFIED)
131 .value("top", EMISSION_TOP)
132 .value("bottom", EMISSION_BOTTOM)
133 .get(emission);
134 if (reader.hasAttribute("wavelength")) { //TODO Remove in the future
135 writelog(LOG_WARNING, "XML line {:d} in <mode>: Attribute 'wavelength' is obsolete, use 'lam' instead", reader.getLineNr());
136 if (reader.hasAttribute("lam")) throw XMLConflictingAttributesException(reader, "wavelength", "lam");
137 k0 = 2e3*PI / reader.requireAttribute<dcomplex>("wavelength");
138 }
139 if (reader.hasAttribute("lam")) k0 = 2e3*PI / reader.requireAttribute<dcomplex>("lam");
140 ktran = reader.getAttribute<dcomplex>("k-tran", ktran);
141 klong = reader.getAttribute<dcomplex>("k-long", klong);
142 std::string sym_tran, sym_long;
143 readComaAttr(reader, "symmetry", sym_long, sym_tran, true);
144 if (sym_long != "") setSymmetryLong(readSymmetry(this, reader, sym_long));
145 if (sym_tran != "") setSymmetryTran(readSymmetry(this, reader, sym_tran));
146 reader.requireTagEnd();
147 } else
148 parseCommonModalConfiguration(reader, manager);
149 }
150}
151
152
153 void FourierSolver3D::onInitialize()
154{
155 this->setupLayers();
156 if (this->interface == -1)
157 Solver::writelog(LOG_DETAIL, "Initializing Fourier3D solver ({0} layers in the stack)",
158 this->stack.size());
159 else
160 Solver::writelog(LOG_DETAIL, "Initializing Fourier3D solver ({0} layers in the stack, interface after {1} layer{2})",
161 this->stack.size(), this->interface, (this->interface==1)? "" : "s");
162 setExpansionDefaults();
163 expansion.init();
164 this->recompute_integrals = true;
165}
166
167
168 void FourierSolver3D::onInvalidate()
169{
170 modes.clear();
171 expansion.reset();
172 transfer.reset();
173}
174
175
176 size_t FourierSolver3D::findMode(FourierSolver3D::What what, dcomplex start)
177{
178 expansion.setSymmetryLong(symmetry_long);
179 expansion.setSymmetryTran(symmetry_tran);
180 expansion.setLam0(this->lam0);
181 Solver::initCalculation();
182 ensureInterface();
183 if (!transfer) initTransfer(expansion, false);
184 std::unique_ptr<RootDigger> root;
185 switch (what) {
186 case FourierSolver3D::WHAT_WAVELENGTH:
187 expansion.setKlong(klong);
188 expansion.setKtran(ktran);
189 root = getRootDigger([this](const dcomplex& x) {
190 if (isnan(x)) throw ComputationError(this->getId(), "'lam' converged to NaN");
191 expansion.setK0(2e3*PI/x); return transfer->determinant();
192 }, "lam");
193 break;
194 case FourierSolver3D::WHAT_K0:
195 expansion.setKlong(klong);
196 expansion.setKtran(ktran);
197 root = getRootDigger([this](const dcomplex& x) {
198 if (isnan(x)) throw ComputationError(this->getId(), "'k0' converged to NaN");
199 expansion.setK0(x); return transfer->determinant();
200 }, "k0");
201 break;
202 case FourierSolver3D::WHAT_KLONG:
203 if (expansion.symmetric_long())
204 throw Exception("{}: Cannot search for longitudinal wavevector with longitudinal symmetry", getId());
205 expansion.setK0(this->k0);
206 expansion.setKtran(ktran);
207 transfer->fields_determined = Transfer::DETERMINED_NOTHING;
208 root = getRootDigger([this](const dcomplex& x) {
209 if (isnan(x)) throw ComputationError(this->getId(), "'klong' converged to NaN");
210 expansion.klong = x; return transfer->determinant();
211 }, "klong");
212 break;
213 case FourierSolver3D::WHAT_KTRAN:
214 if (expansion.symmetric_tran())
215 throw Exception("{}: Cannot search for transverse wavevector with transverse symmetry", getId());
216 expansion.setK0(this->k0);
217 expansion.setKlong(klong);
218 transfer->fields_determined = Transfer::DETERMINED_NOTHING;
219 root = getRootDigger([this](const dcomplex& x) {
220 if (isnan(x)) throw ComputationError(this->getId(), "'ktran' converged to NaN");
221 expansion.ktran = x; return transfer->determinant();
222 }, "ktran");
223 break;
224 }
225 root->find(start);
226 return insertMode();
227}
228
229
230 double FourierSolver3D::getWavelength(size_t n) {
231 if (n >= modes.size()) throw NoValue(ModeWavelength::NAME);
232 return real(2e3*M_PI / modes[n].k0);
233}
234
235
236 size_t FourierSolver3D::initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam)
237{
238
239 if (polarization == ExpansionPW3D::E_UNSPECIFIED)
240 throw BadInput(getId(), "unspecified incident polarization for reflectivity computation");
241 if (expansion.symmetry_long == Expansion::Component(3-polarization))
242 throw BadInput(getId(), "current longitudinal symmetry is inconsistent with the specified incident polarization");
243 if (expansion.symmetry_tran == Expansion::Component(3-polarization))
244 throw BadInput(getId(), "current transverse symmetry is inconsistent with the specified incident polarization");
245 return ModalSolver<SolverOver<Geometry3D>>::initIncidence(side, lam);
246}
247
248 cvector FourierSolver3D::incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam)
249{
250 size_t layer = initIncidence(side, polarization, lam);
251
252 size_t idx = (polarization == ExpansionPW3D::E_LONG)? expansion.iEx(0,0) : expansion.iEy(0,0);
253 cvector physical(expansion.matrixSize(), 0.);
254 physical[idx] = 1.;
255
256 cvector incident = transfer->diagonalizer->invTE(layer) * physical;
257 scaleIncidentVector(incident, layer);
258 return incident;
259}
260
261
262 cvector FourierSolver3D::incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization,
263 double sigma_long, double sigma_tran, double center_long, double center_tran,
264 dcomplex lam)
265{
266 size_t layer = initIncidence(side, polarization, lam);
267
268 double bl = 2.*PI / (expansion.front-expansion.back) * (expansion.symmetric_long()? 0.5 : 1.0),
269 bt = 2.*PI / (expansion.right-expansion.left) * (expansion.symmetric_tran()? 0.5 : 1.0);
270 dcomplex dl = I * bl * (center_long - expansion.back), dt = I * bt * (center_tran - expansion.left);
271 double cl2 = - 0.5 * sigma_long*sigma_long * bl*bl, ct2 = - 0.5 * sigma_tran*sigma_tran * bt*bt;
272
273 cvector physical(expansion.matrixSize(), 0.);
274 for (int it = -int(size_tran); it <= int(size_tran); ++it) {
275 dcomplex vt = exp(ct2 * double(it*it) - dt*double(it));
276 for (int il = -int(size_long); il <= int(size_long); ++il) {
277 size_t idx = (polarization == Expansion::E_LONG)? expansion.iEx(il, it) : expansion.iEy(il, it);
278 physical[idx] = vt * exp(cl2 * double(il*il) - dl*double(il));
279 }
280 }
281
282 cvector incident = transfer->diagonalizer->invTE(layer) * physical;
283 scaleIncidentVector(incident, layer);
284 return incident;
285}
286
287 LazyData<double> FourierSolver3D::getGradients(GradientFunctions::EnumType what,
288 const shared_ptr<const MeshD<3>>& dst_mesh,
289 InterpolationMethod interp) {
290 initCalculation();
291 computeIntegrals();
292 DataVector<double> destination(dst_mesh->size());
293 auto levels = makeLevelsAdapter(dst_mesh);
294 while (auto level = levels->yield()) {
295 auto dest = expansion.getGradients(what, level, interp);
296 for (size_t i = 0; i != level->size(); ++i) destination[level->index(i)] = dest[i];
297 }
298 return destination;
299}
300
301}}} // namespace