PLaSK library

Loading...
Searching...
No Matches
solvercyl.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 <memory>
15
16#include "solvercyl.hpp"
17
18 namespace plask { namespace optical { namespace modal {
19
20 BesselSolverCyl::BesselSolverCyl(const std::string& name)
21 : ModalSolver<SolverWithMesh<Geometry2DCylindrical, MeshAxis>>(name),
22 domain(DOMAIN_INFINITE),
23 m(1),
24 size(12),
25 rule(RULE_DIRECT),
26 kscale(1.),
27 kmax(5.),
28 kmethod(WAVEVECTORS_NONUNIFORM),
29 integral_error(1e-6),
30 max_integration_points(1000),
31 outLoss(this, &BesselSolverCyl::getModalLoss, &BesselSolverCyl::nummodes) {
32 pml.dist = 20.;
33 pml.size = 0.;
34 this->writelog(LOG_WARNING, "This is an EXPERIMENTAL solver! Calculation results may not be reliable!");
35}
36
37 void BesselSolverCyl::loadConfiguration(XMLReader& reader, Manager& manager) {
38 while (reader.requireTagOrEnd()) {
39 std::string param = reader.getNodeName();
40 if (param == "expansion") {
41 domain = reader.enumAttribute<BesselDomain>("domain")
42 .value("finite", DOMAIN_FINITE)
43 .value("infinite", DOMAIN_INFINITE)
44 .get(domain);
45 size = reader.getAttribute<size_t>("size", size);
46 group_layers = reader.getAttribute<bool>("group-layers", group_layers);
47 lam0 = reader.getAttribute<double>("lam0", NAN);
48 always_recompute_gain = reader.getAttribute<bool>("update-gain", always_recompute_gain);
49 max_temp_diff = reader.getAttribute<double>("temp-diff", max_temp_diff);
50 temp_dist = reader.getAttribute<double>("temp-dist", temp_dist);
51 temp_layer = reader.getAttribute<double>("temp-layer", temp_layer);
52 integral_error = reader.getAttribute<double>("integrals-error", integral_error);
53 max_integration_points = reader.getAttribute<size_t>("integrals-points", max_integration_points);
54 kscale = reader.getAttribute<double>("k-scale", kscale);
55 kmax = reader.getAttribute<double>("k-max", kmax);
56 kmethod = reader.enumAttribute<InfiniteWavevectors>("k-method")
57 .value("uniform", WAVEVECTORS_UNIFORM)
58 .value("nonuniform", WAVEVECTORS_NONUNIFORM)
59 .value("non-uniform", WAVEVECTORS_NONUNIFORM)
60 .value("laguerre", WAVEVECTORS_LAGUERRE)
61 .value("manual", WAVEVECTORS_MANUAL)
62 .get(kmethod);
63 if (reader.hasAttribute("k-list")) {
64 klist.clear();
65 for (auto val : boost::tokenizer<boost::char_separator<char>>(reader.requireAttribute("k-list"),
66 boost::char_separator<char>(" ,;\t\n"))) {
67 try {
68 double val = boost::lexical_cast<double>(val);
69 klist.push_back(val);
70 } catch (boost::bad_lexical_cast&) {
71 throw XMLException(reader, format("value '{0}' cannot be converted to float", val));
72 }
73 }
74 }
75 if (reader.hasAttribute("k-weights")) {
76 std::string value = reader.requireAttribute("k-weights");
77 if (value.empty() || value == "auto")
78 kweights.reset();
79 else {
80 kweights.reset(std::vector<double>());
81 kweights->clear();
82 for (auto val : boost::tokenizer<boost::char_separator<char>>(value, boost::char_separator<char>(" ,;\t\n"))) {
83 try {
84 double val = boost::lexical_cast<double>(val);
85 kweights->push_back(val);
86 } catch (boost::bad_lexical_cast&) {
87 throw XMLException(reader, format("value '{0}' cannot be converted to float", val));
88 }
89 }
90 }
91 }
92 if (reader.hasAttribute("rule")) {
93 rule = reader.enumAttribute<Rule>("rule")
94 .value("direct", RULE_DIRECT)
95 .value("combined1", RULE_COMBINED_1)
96 .value("combined2", RULE_COMBINED_2)
97 .value("old", RULE_OLD)
98 .require();
99 }
100 reader.requireTagEnd();
101 } else if (param == "mode") {
102 emission = reader.enumAttribute<Emission>("emission")
103 .value("undefined", EMISSION_UNSPECIFIED)
104 .value("top", EMISSION_TOP)
105 .value("bottom", EMISSION_BOTTOM)
106 .get(emission);
107 if (reader.hasAttribute("wavelength")) { // TODO Remove in the future
108 writelog(LOG_WARNING, "XML line {:d} in <mode>: Attribute 'wavelength' is obsolete, use 'lam' instead",
109 reader.getLineNr());
110 if (reader.hasAttribute("lam")) throw XMLConflictingAttributesException(reader, "wavelength", "lam");
111 k0 = 2e3 * PI / reader.requireAttribute<dcomplex>("wavelength");
112 }
113 if (reader.hasAttribute("lam")) k0 = 2e3 * PI / reader.requireAttribute<dcomplex>("lam");
114 reader.requireTagEnd();
115 } else if (param == "pml") {
116 pml.factor = reader.getAttribute<dcomplex>("factor", pml.factor);
117 pml.size = reader.getAttribute<double>("size", pml.size);
118 pml.dist = reader.getAttribute<double>("dist", pml.dist);
119 if (reader.hasAttribute("order")) { // TODO Remove in the future
120 writelog(LOG_WARNING, "XML line {:d} in <pml>: Attribute 'order' is obsolete, use 'shape' instead",
121 reader.getLineNr());
122 pml.order = reader.requireAttribute<double>("order");
123 }
124 pml.order = reader.getAttribute<double>("shape", pml.order);
125 reader.requireTagEnd();
126 } else
127 parseCommonModalConfiguration(reader, manager);
128 }
129}
130
131 void BesselSolverCyl::onInitialize() {
132 if (size == 0)
133 throw BadInput(getId(), "bessel solver size cannot be 0");
134 this->setupLayers();
135 std::string dom;
136 switch (domain) {
137 case DOMAIN_FINITE: dom = "finite"; break;
138 case DOMAIN_INFINITE: dom = "infinite"; break;
139 default: assert(0);
140 }
141
142 if (this->interface == -1)
143 Solver::writelog(LOG_DETAIL, "Initializing BesselCyl solver in {} domain ({} layers in the stack)", dom,
144 this->stack.size());
145 else
146 Solver::writelog(LOG_DETAIL,
147 "Initializing BesselCyl solver in {} domain ({} layers in the stack, interface after {} layer{})", dom,
148 this->stack.size(), this->interface, (this->interface == 1) ? "" : "s");
149 switch (domain) {
150 case DOMAIN_FINITE: expansion.reset(new ExpansionBesselFini(this)); break;
151 case DOMAIN_INFINITE: expansion.reset(new ExpansionBesselInfini(this)); break;
152 default: assert(0);
153 }
154 setExpansionDefaults();
155 expansion->init1();
156 this->recompute_integrals = true;
157}
158
159 void BesselSolverCyl::onInvalidate() {
160 modes.clear();
161 expansion->reset();
162 transfer.reset();
163}
164
165 size_t BesselSolverCyl::findMode(dcomplex start, int m) {
166 Solver::initCalculation();
167 ensureInterface();
168 expansion->setLam0(this->lam0);
169 expansion->setM(m);
170 initTransfer(*expansion, false);
171 std::unique_ptr<RootDigger> root = getRootDigger(
172 [this](const dcomplex& x) {
173 if (isnan(x)) throw ComputationError(this->getId(), "'lam' converged to NaN");
174 expansion->setK0(2e3 * PI / x);
175 return transfer->determinant();
176 },
177 "lam");
178 root->find(start);
179 return insertMode();
180}
181
182
183 double BesselSolverCyl::getWavelength(size_t n) {
184 if (n >= modes.size()) throw NoValue(ModeWavelength::NAME);
185 return real(2e3 * M_PI / modes[n].k0);
186}
187
188#ifndef NDEBUG
189 cmatrix BesselSolverCyl::epsV_k(size_t layer) {
190 Solver::initCalculation();
191 computeIntegrals();
192 return expansion->epsV_k(layer);
193}
194 cmatrix BesselSolverCyl::epsTss(size_t layer) {
195 Solver::initCalculation();
196 computeIntegrals();
197 return expansion->epsTss(layer);
198}
199 cmatrix BesselSolverCyl::epsTpp(size_t layer) {
200 Solver::initCalculation();
201 computeIntegrals();
202 return expansion->epsTpp(layer);
203}
204 cmatrix BesselSolverCyl::epsTsp(size_t layer) {
205 Solver::initCalculation();
206 computeIntegrals();
207 return expansion->epsTsp(layer);
208}
209 cmatrix BesselSolverCyl::epsTps(size_t layer) {
210 Solver::initCalculation();
211 computeIntegrals();
212 return expansion->epsTps(layer);
213}
214
215 cmatrix BesselSolverCyl::muV_k() {
216 Solver::initCalculation();
217 if (auto finite_expansion = dynamic_cast<ExpansionBesselFini*>(expansion.get())) {
218 computeIntegrals();
219 return finite_expansion->muV_k();
220 } else {
221 return cmatrix();
222 }
223}
224 cmatrix BesselSolverCyl::muTss() {
225 Solver::initCalculation();
226 if (auto finite_expansion = dynamic_cast<ExpansionBesselFini*>(expansion.get())) {
227 computeIntegrals();
228 return finite_expansion->muTss();
229 } else {
230 return cmatrix();
231 }
232}
233 cmatrix BesselSolverCyl::muTsp() {
234 Solver::initCalculation();
235 if (auto finite_expansion = dynamic_cast<ExpansionBesselFini*>(expansion.get())) {
236 computeIntegrals();
237 return finite_expansion->muTsp();
238 } else {
239 return cmatrix();
240 }
241}
242 cmatrix BesselSolverCyl::muTps() {
243 Solver::initCalculation();
244 if (auto finite_expansion = dynamic_cast<ExpansionBesselFini*>(expansion.get())) {
245 computeIntegrals();
246 return finite_expansion->muTps();
247 } else {
248 return cmatrix();
249 }
250}
251 cmatrix BesselSolverCyl::muTpp() {
252 Solver::initCalculation();
253 if (auto finite_expansion = dynamic_cast<ExpansionBesselFini*>(expansion.get())) {
254 computeIntegrals();
255 return finite_expansion->muTpp();
256 } else {
257 return cmatrix();
258 }
259}
260#endif
261
262}}} // namespace plask::optical::modal