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