PLaSK library
Loading...
Searching...
No Matches
solver2d.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 "solver2d.hpp"
15#include "expansion2d.hpp"
16#include "../diagonalizer.hpp"
17
18namespace plask { namespace optical { namespace modal {
19
20FourierSolver2D::FourierSolver2D(const std::string& name):
21 ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>>(name),
22 beta(0.), ktran(0.),
23 symmetry(Expansion::E_UNSPECIFIED),
24 polarization(Expansion::E_UNSPECIFIED),
25 size(12),
26 dct(2),
27 ftt(FOURIER_DISCRETE),
28 expansion(this),
29 refine(32),
30 outNeff(this, &FourierSolver2D::getEffectiveIndex, &FourierSolver2D::nummodes)
31{}
32
33
34void FourierSolver2D::loadConfiguration(XMLReader& reader, Manager& manager)
35{
36 while (reader.requireTagOrEnd()) {
37 std::string param = reader.getNodeName();
38 if (param == "expansion") {
39 size = reader.getAttribute<size_t>("size", size);
40 refine = reader.getAttribute<size_t>("refine", refine);
41 smooth = reader.getAttribute<double>("smooth", smooth);
42 if (reader.hasAttribute("oversampling")) {
43 reader.ignoreAttribute("oversampling");
44 writelog(LOG_WARNING, "obsolete 'oversampling' attribute in XML line {}", reader.getLineNr());
45 }
46 ftt = reader.enumAttribute<FourierType>("ft")
47 .value("discrete", FOURIER_DISCRETE)
48 .value("analytic", FOURIER_ANALYTIC)
49 .get(ftt);
50 int dc = reader.getAttribute<int>("dct", dct);
51 if (dc != 1 && dc != 2)
52 throw XMLBadAttrException(reader, "dct", boost::lexical_cast<std::string>(dc), "\"1\" or \"2\"");
53 dct = dc;
54 group_layers = reader.getAttribute<bool>("group-layers", group_layers);
55 lam0 = reader.getAttribute<double>("lam0", NAN);
56 always_recompute_gain = reader.getAttribute<bool>("update-gain", always_recompute_gain);
57 max_temp_diff = reader.getAttribute<double>("temp-diff", max_temp_diff);
58 temp_dist = reader.getAttribute<double>("temp-dist", temp_dist);
59 temp_layer = reader.getAttribute<double>("temp-layer", temp_layer);
60 reader.requireTagEnd();
61 } else if (param == "pml") {
62 pml.factor = reader.getAttribute<dcomplex>("factor", pml.factor);
63 pml.size = reader.getAttribute<double>("size", pml.size);
64 pml.dist = reader.getAttribute<double>("dist", pml.dist);
65 if (reader.hasAttribute("order")) { //TODO Remove in the future
66 writelog(LOG_WARNING, "XML line {:d} in <pml>: Attribute 'order' is obsolete, use 'shape' instead", reader.getLineNr());
67 if (reader.hasAttribute("shape")) throw XMLConflictingAttributesException(reader, "order", "shape");
68 pml.order = reader.requireAttribute<double>("order");
69 }
70 pml.order = reader.getAttribute<double>("shape", pml.order);
71 reader.requireTagEnd();
72 } else if (param == "mode") {
73 emission = reader.enumAttribute<Emission>("emission")
74 .value("undefined", EMISSION_UNSPECIFIED)
75 .value("top", EMISSION_TOP)
76 .value("bottom", EMISSION_BOTTOM)
77 .get(emission);
78 if (reader.hasAttribute("wavelength")) { //TODO Remove in the future
79 writelog(LOG_WARNING, "XML line {:d} in <mode>: Attribute 'wavelength' is obsolete, use 'lam' instead", reader.getLineNr());
80 if (reader.hasAttribute("lam")) throw XMLConflictingAttributesException(reader, "wavelength", "lam");
81 k0 = 2e3*PI / reader.requireAttribute<dcomplex>("wavelength");
82 }
83 if (reader.hasAttribute("lam")) k0 = 2e3*PI / reader.requireAttribute<dcomplex>("lam");
84 ktran = reader.getAttribute<dcomplex>("k-tran", ktran);
85 if (reader.hasAttribute("beta")) {
86 if (reader.hasAttribute("k-long")) throw XMLConflictingAttributesException(reader, "beta", "k-long");
87 beta = reader.requireAttribute<dcomplex>("beta");
88 } else
89 beta = reader.getAttribute<dcomplex>("k-long", beta);
90 if (reader.hasAttribute("symmetry")) {
91 std::string repr = reader.requireAttribute("symmetry");
92 Expansion::Component val;
93 AxisNames* axes = nullptr;
94 if (geometry) axes = &geometry->axisNames;
95 if (repr == "none" || repr == "NONE" || repr == "None")
96 val = Expansion::E_UNSPECIFIED;
97 else if (repr == "Etran" || repr == "Et" || (axes && repr == "E"+axes->getNameForTran()) ||
98 repr == "Hlong" || repr == "Hl" || (axes && repr == "H"+axes->getNameForLong()))
99 val = Expansion::E_TRAN;
100 else if (repr == "Elong" || repr == "El" || (axes && repr == "E"+axes->getNameForLong()) ||
101 repr == "Htran" || repr == "Ht" || (axes && repr == "H"+axes->getNameForTran()))
102 val = Expansion::E_LONG;
103 else
104 throw XMLBadAttrException(reader, "symmetry", repr, "symmetric field component name (maybe you need to specify the geometry first)");
105 setSymmetry(val);
106 }
107 if (reader.hasAttribute("polarization")) {
108 std::string repr = reader.requireAttribute("polarization");
109 Expansion::Component val;
110 AxisNames* axes = nullptr;
111 if (geometry) axes = &geometry->axisNames;
112 if (repr == "none" || repr == "NONE" || repr == "None")
113 val = Expansion::E_UNSPECIFIED;
114 else if (repr == "Etran" || repr == "Et" || (axes && repr == "E"+axes->getNameForTran()) ||
115 repr == "Hlong" || repr == "Hl" || (axes && repr == "H"+axes->getNameForLong()) || repr == "TM")
116 val = Expansion::E_TRAN;
117 else if (repr == "Elong" || repr == "El" || (axes && repr == "E"+axes->getNameForLong()) ||
118 repr == "Htran" || repr == "Ht" || (axes && repr == "H"+axes->getNameForTran()) || repr == "TE")
119 val = Expansion::E_LONG;
120 else
121 throw XMLBadAttrException(reader, "polarization", repr, "existing field component name (maybe you need to specify the geometry first)");
122 setPolarization(val);
123 }
124 reader.requireTagEnd();
125 } else if (param == "mirrors") {
126 double R1 = reader.requireAttribute<double>("R1");
127 double R2 = reader.requireAttribute<double>("R2");
128 mirrors.reset(std::make_pair(R1,R2));
129 reader.requireTagEnd();
130 } else
131 parseCommonModalConfiguration(reader, manager);
132 }
133}
134
135
136void FourierSolver2D::onInitialize()
137{
138 this->setupLayers();
139 if (this->interface == -1)
140 Solver::writelog(LOG_DETAIL, "Initializing Fourier2D solver ({0} layers in the stack)",
141 this->stack.size());
142 else
143 Solver::writelog(LOG_DETAIL, "Initializing Fourier2D solver ({0} layers in the stack, interface after {1} layer{2})",
144 this->stack.size(), this->interface, (this->interface==1)? "" : "s");
145 setExpansionDefaults();
146 expansion.init();
147 this->recompute_integrals = true;
148}
149
150
151void FourierSolver2D::onInvalidate()
152{
153 modes.clear();
154 expansion.reset();
155 transfer.reset();
156}
157
158
159size_t FourierSolver2D::findMode(FourierSolver2D::What what, dcomplex start)
160{
161 expansion.setSymmetry(symmetry);
162 expansion.setPolarization(polarization);
163 expansion.setLam0(this->lam0);
164 Solver::initCalculation();
165 ensureInterface();
166 if (!transfer) initTransfer(expansion, false);
167 std::unique_ptr<RootDigger> root;
168 switch (what) {
169 case FourierSolver2D::WHAT_WAVELENGTH:
170 expansion.setBeta(beta);
171 expansion.setKtran(ktran);
172 root = getRootDigger([this](const dcomplex& x) {
173 if (isnan(x)) throw ComputationError(this->getId(), "'lam' converged to NaN");
174 expansion.setK0(2e3*PI/x); return transfer->determinant();
175 }, "lam");
176 break;
177 case FourierSolver2D::WHAT_K0:
178 expansion.setBeta(beta);
179 expansion.setKtran(ktran);
180 root = getRootDigger([this](const dcomplex& x) {
181 if (isnan(x)) throw ComputationError(this->getId(), "'k0' converged to NaN");
182 expansion.setK0(x); return transfer->determinant();
183 }, "k0");
184 break;
185 case FourierSolver2D::WHAT_NEFF:
186 if (expansion.separated())
187 throw Exception("{0}: Cannot search for effective index with polarization separation", getId());
188 expansion.setK0(k0);
189 expansion.setKtran(ktran);
190 clearFields();
191 root = getRootDigger([this](const dcomplex& x) {
192 if (isnan(x)) throw ComputationError(this->getId(), "'neff' converged to NaN");
193 expansion.setBeta(x * expansion.k0); return transfer->determinant();
194 }, "neff");
195 break;
196 case FourierSolver2D::WHAT_KTRAN:
197 if (expansion.symmetric())
198 throw Exception("{0}: Cannot search for transverse wavevector with symmetry", getId());
199 expansion.setK0(k0);
200 expansion.setBeta(beta);
201 root = getRootDigger([this](const dcomplex& x) {
202 if (isnan(x)) throw ComputationError(this->getId(), "'ktran' converged to NaN");
203 expansion.setKtran(x); return transfer->determinant();
204 }, "ktran");
205 break;
206 case FourierSolver2D::WHAT_BETA:
207 if (expansion.separated())
208 throw Exception("{0}: Cannot search for longitudinal wavevector with polarization separation", getId());
209 expansion.setK0(k0);
210 expansion.setKtran(ktran);
211 clearFields();
212 root = getRootDigger([this](const dcomplex& x) {
213 if (isnan(x)) throw ComputationError(this->getId(), "'beta' converged to NaN");
214 expansion.setBeta(x); return transfer->determinant();
215 }, "beta");
216 break;
217 }
218 root->find(start);
219 return insertMode();
220}
221
222
223double FourierSolver2D::getWavelength(size_t n) {
224 if (n >= modes.size()) throw NoValue(ModeWavelength::NAME);
225 return real(2e3*M_PI / modes[n].k0);
226}
227
228
229size_t FourierSolver2D::initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam) {
230 if (!isinf(geometry->getExtrusion()->getLength()))
231 throw Exception("{}: Reflectivity computation for 2D geometries possible only if the extrusion length is infinite",
232 getId());
233 if (polarization == Expansion::E_UNSPECIFIED)
234 throw BadInput(getId(), "unspecified incident polarization for reflectivity computation");
235 if (expansion.symmetric() && expansion.symmetry != polarization)
236 throw BadInput(getId(), "current solver symmetry is inconsistent with the specified incident polarization");
237 if (expansion.separated() && expansion.polarization != polarization)
238 throw BadInput(getId(), "current solver polarization is inconsistent with the specified incident polarization");
239 return ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>>::initIncidence(side, lam);
240}
241
242cvector FourierSolver2D::incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam)
243{
244 size_t layer = initIncidence(side, polarization, lam);
245
246 size_t idx;
247 if (expansion.separated()) idx = expansion.iEH(0);
248 else idx = (polarization == Expansion::E_TRAN)? expansion.iEx(0) : expansion.iEz(0);
249 cvector physical(expansion.matrixSize(), 0.);
250 physical[idx] = (polarization == Expansion::E_TRAN)? 1. : -1.;
251
252 cvector incident = transfer->diagonalizer->invTE(layer) * physical;
253 scaleIncidentVector(incident, layer);
254 return incident;
255}
256
257
258cvector FourierSolver2D::incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization, double sigma, double center, dcomplex lam)
259{
260 size_t layer = initIncidence(side, polarization, lam);
261
262 double b = 2.*PI / (expansion.right-expansion.left) * (expansion.symmetric()? 0.5 : 1.0);
263 dcomplex d = I * b * (center - expansion.left);
264 double c2 = - 0.5 * sigma*sigma * b*b;
265
266 cvector physical(expansion.matrixSize(), 0.);
267 for (int i = -int(size); i <= int(size); ++i) {
268 size_t idx;
269 if (expansion.separated()) idx = expansion.iEH(i);
270 else idx = (polarization == Expansion::E_TRAN)? expansion.iEx(i) : expansion.iEz(i);
271 dcomplex val = exp(c2 * double(i*i) - d*double(i));
272 physical[idx] = (polarization == Expansion::E_TRAN)? val : -val;
273 }
274
275 cvector incident = transfer->diagonalizer->invTE(layer) * physical;
276 scaleIncidentVector(incident, layer);
277 return incident;
278}
279
280
281}}} // namespace