PLaSK library
Loading...
Searching...
No Matches
solver.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 "solver.hpp"
15#include "../plask/math.hpp"
16#include "admittance.hpp"
17#include "brent.hpp"
18#include "broyden.hpp"
19#include "diagonalizer.hpp"
20#include "expansion.hpp"
21#include "impedance.hpp"
22#include "meshadapter.hpp"
23#include "muller.hpp"
24#include "reflection.hpp"
25
26namespace plask { namespace optical { namespace modal {
27
28void ModalBase::initTransfer(Expansion& expansion, bool reflection) {
29 switch (transfer_method) {
30 case Transfer::METHOD_REFLECTION_ADMITTANCE:
31 case Transfer::METHOD_REFLECTION_IMPEDANCE: reflection = true; break;
32 case Transfer::METHOD_ADMITTANCE:
33 case Transfer::METHOD_IMPEDANCE: reflection = false; break;
34 default: break;
35 }
36 if (reflection) {
37 ReflectionTransfer::Matching matching = (transfer_method == Transfer::METHOD_REFLECTION_IMPEDANCE)
38 ? ReflectionTransfer::MATCH_IMPEDANCE
39 : ReflectionTransfer::MATCH_ADMITTANCE;
40 if (!this->transfer) {
41 ReflectionTransfer* transfer = dynamic_cast<ReflectionTransfer*>(this->transfer.get());
42 if (!transfer || transfer->diagonalizer->source() != &expansion || transfer->matching != matching)
43 this->transfer.reset(new ReflectionTransfer(this, expansion, matching));
44 }
45 } else {
46 if (transfer_method == Transfer::METHOD_IMPEDANCE) {
47 if (!this->transfer || !dynamic_cast<ImpedanceTransfer*>(this->transfer.get()) ||
48 this->transfer->diagonalizer->source() != &expansion)
49 this->transfer.reset(new ImpedanceTransfer(this, expansion));
50 } else {
51 if (!this->transfer || !dynamic_cast<AdmittanceTransfer*>(this->transfer.get()) ||
52 this->transfer->diagonalizer->source() != &expansion)
53 this->transfer.reset(new AdmittanceTransfer(this, expansion));
54 }
55 }
56}
57
58template <typename BaseT>
59ModalSolver<BaseT>::ModalSolver(const std::string& name)
60 : BaseT(name),
61 smooth(0.),
62 outEpsilon(this, &ModalSolver<BaseT>::getEpsilonProfile),
63 outRefractiveIndex(this, &ModalSolver<BaseT>::getRefractiveIndex),
64 outWavelength(this, &ModalSolver<BaseT>::getWavelength, &ModalSolver<BaseT>::nummodes),
65 outLightMagnitude(this, &ModalSolver<BaseT>::getLightMagnitude, &ModalSolver<BaseT>::nummodes),
66 outLightE(this, &ModalSolver<BaseT>::getLightE<>, &ModalSolver<BaseT>::nummodes),
67 outLightH(this, &ModalSolver<BaseT>::getLightH<>, &ModalSolver<BaseT>::nummodes),
68 outUpwardsLightE(this, &ModalSolver<BaseT>::getLightE<PROPAGATION_UPWARDS>, &ModalSolver<BaseT>::nummodes),
69 outUpwardsLightH(this, &ModalSolver<BaseT>::getLightH<PROPAGATION_UPWARDS>, &ModalSolver<BaseT>::nummodes),
70 outDownwardsLightE(this, &ModalSolver<BaseT>::getLightE<PROPAGATION_DOWNWARDS>, &ModalSolver<BaseT>::nummodes),
71 outDownwardsLightH(this, &ModalSolver<BaseT>::getLightH<PROPAGATION_DOWNWARDS>, &ModalSolver<BaseT>::nummodes) {
72 inTemperature = 300.; // temperature receiver has some sensible value
73 this->inTemperature.changedConnectMethod(this, &ModalSolver<BaseT>::onInputChanged);
74 this->inGain.changedConnectMethod(this, &ModalSolver<BaseT>::onGainChanged);
75 this->inCarriersConcentration.changedConnectMethod(this, &ModalSolver<BaseT>::onInputChanged);
76}
77
78template <typename BaseT> ModalSolver<BaseT>::~ModalSolver() {
79 this->inTemperature.changedDisconnectMethod(this, &ModalSolver<BaseT>::onInputChanged);
80 this->inGain.changedDisconnectMethod(this, &ModalSolver<BaseT>::onGainChanged);
81 this->inCarriersConcentration.changedDisconnectMethod(this, &ModalSolver<BaseT>::onInputChanged);
82}
83
84std::unique_ptr<RootDigger> ModalBase::getRootDigger(const RootDigger::function_type& func, const char* name) {
85 typedef std::unique_ptr<RootDigger> Res;
86 if (root.method == RootDigger::ROOT_MULLER)
87 return Res(new RootMuller(*this, func, root, name));
88 else if (root.method == RootDigger::ROOT_BROYDEN)
89 return Res(new RootBroyden(*this, func, root, name));
90 else if (root.method == RootDigger::ROOT_BRENT)
91 return Res(new RootBrent(*this, func, root, name));
92 throw BadInput(getId(), "wrong root finding method");
93 return Res();
94}
95
96template <typename BaseT> struct LateralMeshAdapter {
97 shared_ptr<RectangularMesh<2>> mesh;
98
99 LateralMeshAdapter(const BaseT* solver) : mesh(makeGeometryGrid(solver->getGeometry())) {}
100
101 void resetMidpoints(const shared_ptr<MeshAxis>& vbounds) {
102 mesh = make_shared<RectangularMesh<2>>(mesh->axis[0]->getMidpointAxis(), vbounds, RectangularMesh<2>::ORDER_10);
103 }
104
105 void resetMidpoints(const shared_ptr<MeshAxis>& vbounds, double spacing) {
106 mesh = make_shared<RectangularMesh<2>>(refineAxis(mesh->axis[0], spacing)->getMidpointAxis(), vbounds,
107 RectangularMesh<2>::ORDER_10);
108 }
109
110 void reset(const shared_ptr<MeshAxis>& verts) {
111 mesh = make_shared<RectangularMesh<2>>(mesh->axis[0], verts, RectangularMesh<2>::ORDER_10);
112 }
113
114 shared_ptr<RectangularMesh<2>> makeMesh(const shared_ptr<MeshAxis>& verts) {
115 return make_shared<RectangularMesh<2>>(mesh->axis[0], verts, RectangularMesh<2>::ORDER_10);
116 }
117
118 shared_ptr<OrderedAxis> vert() { return dynamic_pointer_cast<OrderedAxis>(mesh->vert()); }
119
120 shared_ptr<RectangularMesh<2>> midmesh() const {
121 return make_shared<RectangularMesh<2>>(mesh->axis[0], mesh->axis[1]->getMidpointAxis());
122 }
123
124 size_t size() const { return mesh->axis[0]->size(); }
125
126 size_t idx(size_t i, size_t v) const { return mesh->index(i, v); }
127
128 Vec<2> at(size_t i, size_t v) const { return mesh->at(i, v); }
129
130 shared_ptr<RectangularMesh<2>> makeLine(size_t i, size_t v, double spacing) const {
131 shared_ptr<OrderedAxis> vaxis(new OrderedAxis({mesh->axis[1]->at(v - 1), mesh->axis[1]->at(v)}));
132 vaxis = refineAxis(vaxis, spacing);
133 return make_shared<RectangularMesh<2>>(make_shared<OnePointAxis>(mesh->axis[0]->at(i)), vaxis);
134 }
135};
136
137template <> struct LateralMeshAdapter<SolverOver<Geometry3D>> {
138 private:
139 size_t _size;
140
141 public:
142 shared_ptr<RectangularMesh<3>> mesh;
143
144 LateralMeshAdapter(const SolverOver<Geometry3D>* solver) : mesh(make_shared<RectangularMesh<3>>()) {
145 std::set<double> xys[2];
146 for (const auto& bbox : solver->getGeometry()->getLeafsBoundingBoxes()) {
147 xys[0].insert(bbox.lower[0]);
148 xys[0].insert(bbox.upper[0]);
149 xys[1].insert(bbox.lower[1]);
150 xys[1].insert(bbox.upper[1]);
151 }
152 // We divide each rectangle into three points to correctly consider circles and triangles
153 for (int ax = 0; ax != 2; ++ax) {
154 std::vector<double> verts;
155 if (!xys[ax].empty()) {
156 verts.reserve(3 * xys[ax].size() - 2);
157 auto it = xys[ax].begin();
158 double x = *it;
159 verts.push_back(x);
160 for (++it; it != xys[ax].end(); ++it) {
161 verts.push_back((2. * x + *it) / 3.);
162 verts.push_back((x + 2. * *it) / 3.);
163 verts.push_back(*it);
164 x = *it;
165 }
166 }
167 mesh->setAxis(ax, shared_ptr<OrderedAxis>(new OrderedAxis(verts)), false);
168 }
169 _size = mesh->axis[0]->size() * mesh->axis[1]->size();
170 mesh->setAxis(2, makeGeometryAxis(solver->getGeometry()->getChild(), Primitive<3>::DIRECTION_VERT), false);
171 }
172
173 void resetMidpoints(const shared_ptr<MeshAxis>& vbounds) {
174 mesh = make_shared<RectangularMesh<3>>(mesh->axis[0]->getMidpointAxis(), mesh->axis[1]->getMidpointAxis(), vbounds,
175 RectangularMesh<3>::ORDER_201);
176 _size = mesh->axis[0]->size() * mesh->axis[1]->size();
177 }
178
179 void resetMidpoints(const shared_ptr<MeshAxis>& vbounds, double spacing) {
180 mesh = make_shared<RectangularMesh<3>>(refineAxis(mesh->axis[0], spacing)->getMidpointAxis(),
181 refineAxis(mesh->axis[1], spacing)->getMidpointAxis(), vbounds,
182 RectangularMesh<3>::ORDER_201);
183 _size = mesh->axis[0]->size() * mesh->axis[1]->size();
184 }
185
186 void reset(const shared_ptr<MeshAxis>& verts) {
187 mesh = make_shared<RectangularMesh<3>>(mesh->axis[0], mesh->axis[1], verts, RectangularMesh<3>::ORDER_210);
188 }
189
190 shared_ptr<RectangularMesh<3>> makeMesh(const shared_ptr<MeshAxis>& verts) {
191 return make_shared<RectangularMesh<3>>(mesh->axis[0], mesh->axis[1], verts, RectangularMesh<3>::ORDER_210);
192 }
193
194 shared_ptr<OrderedAxis> vert() { return dynamic_pointer_cast<OrderedAxis>(mesh->vert()); }
195
196 shared_ptr<RectangularMesh<3>> midmesh() const {
197 return make_shared<RectangularMesh<3>>(mesh->axis[0], mesh->axis[1], mesh->axis[2]->getMidpointAxis());
198 }
199
200 size_t size() const { return _size; }
201
202 size_t idx(size_t i, size_t v) const { return _size * v + i; }
203
204 Vec<3> at(size_t i, size_t v) const { return mesh->RectilinearMesh3D::at(idx(i, v)); }
205
206 shared_ptr<RectangularMesh<3>> makeLine(size_t i, size_t v, double spacing) const {
207 shared_ptr<OrderedAxis> vaxis(new OrderedAxis({mesh->axis[2]->at(v - 1), mesh->axis[2]->at(v)}));
208 vaxis = refineAxis(vaxis, spacing);
209 return make_shared<RectangularMesh<3>>(make_shared<OnePointAxis>(mesh->axis[0]->at(mesh->index0(i))),
210 make_shared<OnePointAxis>(mesh->axis[1]->at(mesh->index1(i))), vaxis);
211 }
212};
213
214template <typename BaseT> void ModalSolver<BaseT>::parseCommonModalConfiguration(XMLReader& reader, Manager& manager) {
215 std::string param = reader.getNodeName();
216 if (param == "interface") {
217 if (reader.hasAttribute("index")) {
218 throw XMLException(reader, "Setting interface by layer index is not supported anymore (set it by object or position)");
219 } else if (reader.hasAttribute("position")) {
220 if (reader.hasAttribute("object")) throw XMLConflictingAttributesException(reader, "index", "object");
221 if (reader.hasAttribute("path")) throw XMLConflictingAttributesException(reader, "index", "path");
222 setInterfaceAt(reader.requireAttribute<double>("position"));
223 } else if (reader.hasAttribute("object")) {
224 auto object = manager.requireGeometryObject<GeometryObject>(reader.requireAttribute("object"));
225 PathHints path;
226 if (auto pathattr = reader.getAttribute("path")) path = manager.requirePathHints(*pathattr);
227 setInterfaceOn(object, path);
228 } else if (reader.hasAttribute("path")) {
229 throw XMLUnexpectedAttrException(reader, "path");
230 }
231 reader.requireTagEnd();
232 } else if (param == "vpml") {
233 vpml.factor = reader.getAttribute<dcomplex>("factor", vpml.factor);
234 vpml.size = reader.getAttribute<double>("size", vpml.size);
235 vpml.dist = reader.getAttribute<double>("dist", vpml.dist);
236 if (reader.hasAttribute("order")) { // TODO Remove in the future
237 writelog(LOG_WARNING, "XML line {:d} in <vpml>: Attribute 'order' is obsolete, use 'shape' instead",
238 reader.getLineNr());
239 vpml.order = reader.requireAttribute<double>("order");
240 }
241 vpml.order = reader.getAttribute<double>("shape", vpml.order);
242 reader.requireTagEnd();
243 } else if (param == "transfer") {
244 transfer_method = reader.enumAttribute<Transfer::Method>("method")
245 .value("auto", Transfer::METHOD_AUTO)
246 .value("reflection", Transfer::METHOD_REFLECTION_ADMITTANCE)
247 .value("reflection-admittance", Transfer::METHOD_REFLECTION_ADMITTANCE)
248 .value("reflection-impedance", Transfer::METHOD_REFLECTION_IMPEDANCE)
249 .value("admittance", Transfer::METHOD_ADMITTANCE)
250 .value("impedance", Transfer::METHOD_IMPEDANCE)
251 .get(transfer_method);
252 determinant_type = reader.enumAttribute<Transfer::Determinant>("determinant")
253 .value("eigen", Transfer::DETERMINANT_EIGENVALUE)
254 .value("eigenvalue", Transfer::DETERMINANT_EIGENVALUE)
255 .value("full", Transfer::DETERMINANT_FULL)
256 .get(determinant_type);
257 reader.requireTagEnd();
258 } else if (param == "root") {
259 readRootDiggerConfig(reader);
260 } else {
261 this->parseStandardConfiguration(reader, manager);
262 }
263}
264
265template <typename BaseT> void ModalSolver<BaseT>::setupLayers() {
266 if (!this->geometry) throw NoGeometryException(this->getId());
267
268 LateralMeshAdapter<BaseT> adapter(this);
269
270 vbounds = adapter.vert();
271 if (this->geometry->isSymmetric(Geometry::DIRECTION_VERT)) {
272 std::deque<double> zz;
273 for (double z : *vbounds) zz.push_front(-z);
274 OrderedAxis::WarningOff nowarn(vbounds);
275 vbounds->addOrderedPoints(zz.begin(), zz.end(), zz.size());
276 }
277
278 if (inTemperature.hasProvider() && !isnan(max_temp_diff) && !isinf(max_temp_diff) && !isnan(temp_dist) && !isinf(temp_dist) &&
279 !isnan(temp_layer) && !isinf(temp_layer))
280 adapter.resetMidpoints(vbounds, temp_dist);
281 else
282 adapter.resetMidpoints(vbounds);
283
284 // Divide layers with too large temperature gradient
285 if (inTemperature.hasProvider() && !isnan(max_temp_diff) && !isinf(max_temp_diff) && !isnan(temp_dist) && !isinf(temp_dist) &&
286 !isnan(temp_layer) && !isinf(temp_layer)) {
287 auto temp = inTemperature(adapter.mesh);
288 std::deque<double> refines;
289 for (size_t v = 1; v != vbounds->size(); ++v) {
290 double mdt = 0.;
291 size_t idt;
292 for (size_t i = 0; i != adapter.size(); ++i) {
293 double dt = abs(temp[adapter.idx(i, v)] - temp[adapter.idx(i, v - 1)]);
294 if (dt > mdt) {
295 mdt = dt;
296 idt = i;
297 }
298 }
299 if (mdt > max_temp_diff) {
300 // We need to divide the layer.
301 auto line_mesh = adapter.makeLine(idt, v, temp_layer);
302 auto tmp = inTemperature(line_mesh);
303 auto line = line_mesh->vert();
304 size_t li = 0;
305 for (size_t i = 2; i != line->size(); ++i) {
306 if (abs(tmp[i] - tmp[li]) > max_temp_diff) {
307 li = i - 1;
308 refines.push_back(line->at(li));
309 }
310 }
311 }
312 }
313 vbounds->addOrderedPoints(refines.begin(), refines.end(), refines.size());
314 }
315
316 adapter.reset(vbounds->getMidpointAxis());
317
318 // Add layers below bottom boundary and above top one
319 verts = dynamic_pointer_cast<OrderedAxis>(adapter.mesh->vert());
320 OrderedAxis::WarningOff nowarn(verts);
321 verts->addPoint(vbounds->at(0) - 2. * OrderedAxis::MIN_DISTANCE);
322 verts->addPoint(vbounds->at(vbounds->size() - 1) + 2. * OrderedAxis::MIN_DISTANCE);
323
324 lgained.clear();
325 lcomputed.clear();
326 stack.clear();
327 stack.reserve(verts->size());
328 lcount = 0;
329
330 struct LayerItem {
331 shared_ptr<Material> material;
332 std::set<std::string> roles;
333 bool operator==(const LayerItem& other) { return *material == *other.material && roles == other.roles; }
334 bool operator!=(const LayerItem& other) { return !(*this == other); }
335 };
336 std::vector<std::vector<LayerItem>> layers;
337
338 for (size_t v = 0; v != verts->size(); ++v) {
339 bool gain = false, computed = false;
340 bool unique = !group_layers || layers.size() == 0;
341
342 std::vector<LayerItem> layer(adapter.size());
343 for (size_t i = 0; i != adapter.size(); ++i) {
344 auto p(adapter.at(i, v));
345 layer[i].material = this->geometry->getMaterial(p);
346 for (const std::string& role : this->geometry->getRolesAt(p)) {
347 if (role.substr(0, 3) == "opt")
348 layer[i].roles.insert(role);
349 else if (role == "unique") {
350 layer[i].roles.insert(role);
351 unique = true;
352 } else if (role == "inEpsilon") {
353 layer[i].roles.insert(role);
354 computed = true;
355 } else if (role == "QW" || role == "QD" || role == "gain") {
356 layer[i].roles.insert(role);
357 gain = true;
358 }
359 }
360 }
361
362 if (!unique) {
363 for (size_t l = 0; l != layers.size(); ++l) {
364 unique = false;
365 for (size_t i = 0; i != layers[l].size(); ++i) {
366 if (layers[l][i] != layer[i]) {
367 unique = true;
368 break;
369 }
370 }
371 if (!unique) {
372 stack.push_back(l);
373 break;
374 }
375 }
376 }
377 if (unique) {
378 stack.push_back(lcount++);
379 layers.emplace_back(std::move(layer));
380 lgained.push_back(gain);
381 lcomputed.push_back(computed);
382 }
383 }
384 assert(vbounds->size() == stack.size() - 1);
385 assert(verts->size() == stack.size());
386 assert(layers.size() == lcount);
387
388 // Split groups with too large temperature gradient
389 // We are using CLINK naive algorithm for this purpose
390 if (group_layers && inTemperature.hasProvider() && !isnan(max_temp_diff) && !isinf(max_temp_diff) && !isnan(temp_dist) &&
391 !isinf(temp_dist) && !isnan(temp_layer) && !isinf(temp_layer)) {
392 auto temp = inTemperature(adapter.mesh);
393 size_t nl = lcount; // number of independent layers to consider (stays fixed)
394 for (size_t l = 0; l != nl; ++l) {
395 std::vector<size_t> indices;
396 std::list<std::list<size_t>> groups;
397 typedef std::list<std::list<size_t>>::iterator ListIterator;
398 typedef std::list<size_t>::const_iterator ItemIterator;
399 for (size_t i = 0, j = 0; i != stack.size(); ++i) {
400 if (stack[i] == l) {
401 indices.push_back(i);
402 groups.emplace_back(std::list<size_t>({j++}));
403 }
404 }
405 size_t n = indices.size();
406 // Make distance matrix
407 std::unique_ptr<double[]> dists(new double[n * n]);
408#define dists_at(a, b) dists[(a) * n + (b)] // TODO develop plask::UniquePtr2D<> and remove this macro
409 for (size_t i = 0; i != n; ++i) {
410 dists_at(i, i) = INFINITY; // the simplest way to avoid clustering with itself
411 for (size_t j = i + 1; j != n; ++j) {
412 double mdt = 0.;
413 for (size_t k = 0; k != adapter.size(); ++k) {
414 double dt = abs(temp[adapter.idx(k, indices[i])] - temp[adapter.idx(k, indices[j])]);
415 if (dt > mdt) mdt = dt;
416 }
417 dists_at(i, j) = dists_at(j, i) = mdt;
418 }
419 }
420 // Go and merge groups with the smallest distances
421 while (true) {
422 double mdt = INFINITY;
423 ListIterator mg1, mg2;
424 for (ListIterator g1 = groups.begin(); g1 != groups.end(); ++g1) {
425 ListIterator g2 = g1;
426 for (++g2; g2 != groups.end(); ++g2) {
427 double dt = 0.;
428 for (ItemIterator i1 = g1->begin(); i1 != g1->end(); ++i1)
429 for (ItemIterator i2 = g2->begin(); i2 != g2->end(); ++i2) dt = max(dists_at(*i1, *i2), dt);
430 if (dt < mdt) {
431 mg1 = g1;
432 mg2 = g2;
433 mdt = dt;
434 }
435 }
436 }
437 if (mdt > 0.66667 * max_temp_diff) break;
438 for (ItemIterator i2 = mg2->begin(); i2 != mg2->end(); ++i2) mg1->push_back(*i2);
439 groups.erase(mg2);
440 }
441 // Now update the stack
442 ListIterator g = groups.begin();
443 for (++g; g != groups.end(); ++g) {
444 for (ItemIterator i = g->begin(); i != g->end(); ++i) stack[indices[*i]] = lcount;
445 lgained.push_back(lgained[l]);
446 lcomputed.push_back(lcomputed[l]);
447 ++lcount;
448 }
449 }
450 }
451 assert(lgained.size() == lcount);
452 assert(lcomputed.size() == lcount);
453
454 if (!isnan(interface_position)) {
455 interface = std::lower_bound(vbounds->begin(), vbounds->end(), interface_position - 0.5 * OrderedAxis::MIN_DISTANCE) -
456 vbounds->begin() + 1;
457 // OrderedAxis::MIN_DISTANCE to compensate for truncation errors
458 if (std::size_t(interface) > vbounds->size()) interface = vbounds->size();
459 } else
460 interface = -1;
461
462 // Merge identical adjacent layers
463 std::ptrdiff_t i1 = interface - 1;
464 for (size_t i = 0; i < vbounds->size(); ++i) {
465 if (stack[i] == stack[i + 1] && i != i1) {
466 stack.erase(stack.begin() + i + 1);
467 vbounds->removePoint(i);
468 verts->removePoints(i, i + 2);
469 if (i == 0) {
470 OrderedAxis::WarningOff nowarn(verts);
471 verts->addPoint(vbounds->at(i) - 2. * OrderedAxis::MIN_DISTANCE);
472 } else if (i == vbounds->size()) {
473 OrderedAxis::WarningOff nowarn(verts);
474 verts->addPoint(vbounds->at(i - 1) + 2. * OrderedAxis::MIN_DISTANCE);
475 } else {
476 verts->addPoint(0.5 * (vbounds->at(i - 1) + vbounds->at(i)));
477 }
478 --i;
479 if (i < i1) {
480 --interface;
481 --i1;
482 }
483 assert(vbounds->size() == stack.size() - 1);
484 assert(verts->size() == stack.size());
485 if (vbounds->size() == 1) break; // We always have minimum two layers
486 }
487 }
488
489 Solver::writelog(LOG_DETAIL, "Detected {0} {1}layers", lcount, group_layers ? "distinct " : "");
490
491 // DEBUG
492 // for (size_t i = 0; i < stack.size(); ++i) {
493 // if (i != 0) std::cerr << "---- " << vbounds->at(i-1) << "\n";
494 // std::cerr << stack[i] << (lgained[stack[i]]? "+" : "") << (lgained[stack[i]]? "*\n" : "\n");
495 // }
496
497 if (interface >= 0) {
498 double pos = vbounds->at(interface - 1);
499 if (abs(pos) < OrderedAxis::MIN_DISTANCE) pos = 0.;
500 Solver::writelog(LOG_DEBUG, "Interface is at layer {:d} (exact position {:g}um)", interface, pos);
501 } else {
502 interface = -1;
503 }
504}
505
506template <typename BaseT>
507DataVector<const Tensor3<dcomplex>> ModalSolver<BaseT>::getEpsilonProfile(
508 const shared_ptr<const MeshD<BaseT::SpaceType::DIM>>& dst_mesh,
509 dcomplex lam,
510 InterpolationMethod interp) {
511 if (!isnan(real(lam))) throw BadInput(this->getId(), "wavelength cannot be specified for outEpsilon in this solver");
512 Solver::initCalculation();
513 Expansion& expansion = getExpansion();
514 setExpansionDefaults(false);
515 if (isnan(expansion.lam0) || always_recompute_gain || isnan(expansion.k0)) expansion.setK0(isnan(k0) ? 2e3 * PI / lam0 : k0);
516 // initTransfer(expansion, false);
517 expansion.beforeGetEpsilon();
518
519 // TODO maybe there is a more efficient way to implement this
520 DataVector<Tensor3<dcomplex>> result(dst_mesh->size());
521 auto levels = makeLevelsAdapter(dst_mesh);
522
523 while (auto level = levels->yield()) {
524 double h = level->vpos();
525 size_t n = getLayerFor(h);
526 size_t l = stack[n];
527 auto data = expansion.getMaterialEps(l, level, interp);
528 for (size_t i = 0; i != level->size(); ++i) result[level->index(i)] = data[i];
529 }
530
531 expansion.afterGetEpsilon();
532
533 return result;
534}
535
536template <typename BaseT>
537LazyData<dcomplex> ModalSolver<BaseT>::getRefractiveIndex(RefractiveIndex::EnumType component,
538 const shared_ptr<const MeshD<BaseT::SpaceType::DIM>>& dst_mesh,
539 dcomplex lam,
540 InterpolationMethod interp) {
541 if (!isnan(real(lam))) throw BadInput(this->getId(), "wavelength cannot be specified for outRefractiveIndex in this solver");
542 Solver::initCalculation();
543 DataVector<const Tensor3<dcomplex>> epsilon = getEpsilonProfile(dst_mesh, NAN, interp);
544 switch (component) {
545 case RefractiveIndex::COMPONENT_LONG:
546 return LazyData<dcomplex>(epsilon.size(), [epsilon](size_t i) { return sqrt(epsilon[i].c00); });
547 case RefractiveIndex::COMPONENT_TRAN:
548 return LazyData<dcomplex>(epsilon.size(), [epsilon](size_t i) { return sqrt(epsilon[i].c11); });
549 case RefractiveIndex::COMPONENT_VERT:
550 return LazyData<dcomplex>(epsilon.size(), [epsilon](size_t i) { return sqrt(epsilon[i].c22); });
551 }
552 throw BadInput(getId(), "wrong refractive index component");
553}
554
555#ifndef NDEBUG
556void ModalBase::getMatrices(size_t layer, cmatrix& RE, cmatrix& RH) {
557 initCalculation();
558 computeIntegrals();
559 size_t N = this->getExpansion().matrixSize();
560 if (std::size_t(RE.cols()) != N || std::size_t(RE.rows()) != N) RE.reset(N, N);
561 if (std::size_t(RH.cols()) != N || std::size_t(RH.rows()) != N) RH.reset(N, N);
562 this->getExpansion().getMatrices(layer, RE, RH);
563}
564#endif
565
566template <typename BaseT> size_t ModalSolver<BaseT>::initIncidence(Transfer::IncidentDirection side, dcomplex lam) {
567 Expansion& expansion = getExpansion();
568 bool changed = Solver::initCalculation() || setExpansionDefaults(isnan(lam));
569 if (!isnan(lam)) {
570 dcomplex k0 = 2e3 * M_PI / lam;
571 if (!is_zero(k0 - expansion.getK0())) {
572 expansion.setK0(k0);
573 changed = true;
574 }
575 }
576 size_t layer = stack[(side == Transfer::INCIDENCE_BOTTOM) ? 0 : stack.size() - 1];
577 if (!transfer) {
578 initTransfer(expansion, true);
579 changed = true;
580 }
581 if (changed) {
582 transfer->initDiagonalization();
583 transfer->diagonalizer->diagonalizeLayer(layer);
584 } else if (!transfer->diagonalizer->isDiagonalized(layer))
585 transfer->diagonalizer->diagonalizeLayer(layer);
586 return layer;
587}
588
589template <typename BaseT> cvector ModalSolver<BaseT>::incidentVector(Transfer::IncidentDirection side, size_t idx, dcomplex lam) {
590 size_t layer = initIncidence(side, lam);
591 if (idx >= transfer->diagonalizer->matrixSize()) throw BadInput(getId(), "wrong incident eignenmode index");
592 cvector incident(transfer->diagonalizer->matrixSize(), 0.);
593 incident[idx] = 1.;
594
595 scaleIncidentVector(incident, layer);
596 return incident;
597}
598
599template <typename BaseT>
600cvector ModalSolver<BaseT>::incidentVector(Transfer::IncidentDirection side, const cvector& incident, dcomplex lam) {
601 size_t layer = initIncidence(side, lam);
602 if (incident.size() != transfer->diagonalizer->matrixSize()) throw BadInput(getId(), "wrong incident vector size");
603 cvector result = incident.claim();
604 scaleIncidentVector(result, layer);
605 return result;
606}
607
608void ModalBase::scaleIncidentVector(cvector& incident, size_t layer, double size_factor) {
609 double norm2 = 0.;
610 size_t N = transfer->diagonalizer->matrixSize();
611 for (size_t i = 0; i != N; ++i) {
612 double P = real(incident[i] * conj(incident[i]));
613 if (P != 0.)
614 norm2 += P * getExpansion().getModeFlux(i, transfer->diagonalizer->TE(layer), transfer->diagonalizer->TH(layer));
615 }
616
617 double norm = size_factor / sqrt(abs(norm2));
618 for (size_t i = 0; i != N; ++i) incident[i] *= norm;
619}
620
621template <> void ModalSolver<SolverWithMesh<Geometry2DCartesian, MeshAxis>>::scaleIncidentVector(cvector& incident, size_t layer) {
622 ModalBase::scaleIncidentVector(incident, layer, 1e-3); // sqrt(µm -> m)
623}
624
625template <> void ModalSolver<SolverWithMesh<Geometry2DCylindrical, MeshAxis>>::scaleIncidentVector(cvector& incident, size_t layer) {
626 throw NotImplemented(getId(), "CylindicalSolver::incidentVector");
627}
628
629template <> void ModalSolver<SolverOver<Geometry3D>>::scaleIncidentVector(cvector& incident, size_t layer) {
630 ModalBase::scaleIncidentVector(incident, layer, 1e-6); // sqrt(µm² -> m²)
631}
632
633dvector ModalBase::getIncidentFluxes(const cvector& incident, Transfer::IncidentDirection side) {
634 initCalculation();
635 if (!transfer) initTransfer(getExpansion(), true);
636
637 dvector result(incident.size());
638
639 size_t n = (side == Transfer::INCIDENCE_BOTTOM) ? 0 : stack.size() - 1;
640 size_t l = stack[n];
641
642 size_t N = transfer->diagonalizer->matrixSize();
643
644 Expansion& expansion = getExpansion();
645
646 double input_flux = 0.;
647 for (size_t i = 0; i != N; ++i) {
648 double P = real(incident[i] * conj(incident[i]));
649 if (P != 0.) {
650 result[i] = P * expansion.getModeFlux(i, transfer->diagonalizer->TE(l), transfer->diagonalizer->TH(l));
651 input_flux += result[i];
652 } else
653 result[i] = 0.;
654 }
655
656 result /= input_flux;
657
658 return result;
659}
660
661dvector ModalBase::getReflectedFluxes(const cvector& incident, Transfer::IncidentDirection side) {
662 cvector reflected = getReflectedCoefficients(incident, side);
663 dvector result(reflected.size());
664
665 size_t n = (side == Transfer::INCIDENCE_BOTTOM) ? 0 : stack.size() - 1;
666 size_t l = stack[n];
667
668 size_t N = transfer->diagonalizer->matrixSize();
669
670 Expansion& expansion = getExpansion();
671
672 double input_flux = 0.;
673 for (size_t i = 0; i != N; ++i) {
674 double P = real(incident[i] * conj(incident[i]));
675 if (P != 0.) input_flux += P * expansion.getModeFlux(i, transfer->diagonalizer->TE(l), transfer->diagonalizer->TH(l));
676 }
677
678 for (size_t i = 0; i != N; ++i) {
679 double R = real(reflected[i] * conj(reflected[i]));
680 if (R != 0.)
681 result[i] = R * expansion.getModeFlux(i, transfer->diagonalizer->TE(l), transfer->diagonalizer->TH(l)) / input_flux;
682 else
683 result[i] = 0.;
684 }
685
686 return result;
687}
688
689dvector ModalBase::getTransmittedFluxes(const cvector& incident, Transfer::IncidentDirection side) {
690 cvector transmitted = getTransmittedCoefficients(incident, side);
691 dvector result(transmitted.size());
692
693 size_t ni = (side == Transfer::INCIDENCE_BOTTOM) ? 0 : stack.size() - 1;
694 size_t li = stack[ni];
695
696 size_t nt = stack.size() - 1 - ni; // opposite side than ni
697 size_t lt = stack[nt];
698
699 size_t N = transfer->diagonalizer->matrixSize();
700
701 Expansion& expansion = getExpansion();
702
703 double input_flux = 0.;
704 for (size_t i = 0; i != N; ++i) {
705 double P = real(incident[i] * conj(incident[i]));
706 if (P != 0.) input_flux += P * expansion.getModeFlux(i, transfer->diagonalizer->TE(li), transfer->diagonalizer->TH(li));
707 }
708
709 for (size_t i = 0; i != N; ++i) {
710 double T = real(transmitted[i] * conj(transmitted[i]));
711 if (T != 0.)
712 result[i] = T * expansion.getModeFlux(i, transfer->diagonalizer->TE(lt), transfer->diagonalizer->TH(lt)) / input_flux;
713 else
714 result[i] = 0.;
715 }
716
717 return result;
718}
719
720cvector ModalBase::getReflectedCoefficients(const cvector& incident, Transfer::IncidentDirection side) {
721 initCalculation();
722 if (!transfer) initTransfer(getExpansion(), true);
723
724 return transfer->getReflectionVector(incident, side);
725}
726
727cvector ModalBase::getTransmittedCoefficients(const cvector& incident, Transfer::IncidentDirection side) {
728 initCalculation();
729 if (!transfer) initTransfer(getExpansion(), true);
730
731 return transfer->getTransmissionVector(incident, side);
732}
733
734template <typename BaseT>
735template <PropagationDirection part>
736LazyData<Vec<3, dcomplex>> ModalSolver<BaseT>::getLightE(size_t num,
737 shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
738 InterpolationMethod method) {
739 assert(transfer);
740 double power = applyMode(num);
741 return transfer->getFieldE(power, dst_mesh, method, part);
742}
743
744template <typename BaseT>
745template <PropagationDirection part>
746LazyData<Vec<3, dcomplex>> ModalSolver<BaseT>::getLightH(size_t num,
747 shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
748 InterpolationMethod method) {
749 assert(transfer);
750 double power = applyMode(num);
751 return transfer->getFieldH(power, dst_mesh, method, part);
752}
753
754template <typename BaseT>
755LazyData<double> ModalSolver<BaseT>::getLightMagnitude(size_t num,
756 shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
757 InterpolationMethod method) {
758 assert(transfer);
759 double power = applyMode(num);
760 return transfer->getFieldMagnitude(power, dst_mesh, method);
761}
762
763// template class PLASK_SOLVER_API ModalSolver<SolverOver<Geometry2DCartesian>>;
764template class PLASK_SOLVER_API ModalSolver<SolverWithMesh<Geometry2DCartesian, MeshAxis>>;
765template class PLASK_SOLVER_API ModalSolver<SolverWithMesh<Geometry2DCylindrical, MeshAxis>>;
766template class PLASK_SOLVER_API ModalSolver<SolverOver<Geometry3D>>;
767
768// FiltersFactory::RegisterStandard<Epsilon> registerEpsilonFilters;
769
770}}} // namespace plask::optical::modal