PLaSK library

Loading...
Searching...
No Matches
solverbase.hpp
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#ifndef PLASK__SOLVER_OPTICAL_MODAL_SOLVERBASE_H
15#define PLASK__SOLVER_OPTICAL_MODAL_SOLVERBASE_H
16
17#include <plask/plask.hpp>
18#include "rootdigger.hpp"
19#include "transfer.hpp"
20
21#undef interface
22
23 namespace plask { namespace optical { namespace modal {
24
26 struct PML {
27 dcomplex factor;
28 double size;
29 double dist;
30 double order;
31 PML() : factor(1., 0.), size(1.), dist(0.5), order(1) {}
32 PML(dcomplex factor, double size, double dist, double order) : factor(factor), size(size), dist(dist), order(order) {}
33};
34
38 struct PLASK_SOLVER_API ModalBase {
40 enum Emission {
41 EMISSION_UNSPECIFIED = 0,
42 EMISSION_TOP,
43 EMISSION_BOTTOM,
44 EMISSION_FRONT,
45 EMISSION_BACK
46 };
47
49 Emission emission;
50
52 Transfer::Determinant determinant_type;
53
54 protected:
56 std::unique_ptr<RootDigger> getRootDigger(const RootDigger::function_type& func, const char* name);
57
59 Transfer::Method transfer_method;
60
65 void readRootDiggerConfig(XMLReader& reader) {
66 root.tolx = reader.getAttribute<double>("tolx", root.tolx);
67 root.tolf_min = reader.getAttribute<double>("tolf-min", root.tolf_min);
68 root.tolf_max = reader.getAttribute<double>("tolf-max", root.tolf_max);
69 root.maxstep = reader.getAttribute<double>("maxstep", root.maxstep);
70 root.maxiter = reader.getAttribute<int>("maxiter", root.maxiter);
71 root.alpha = reader.getAttribute<double>("alpha", root.alpha);
72 root.lambda_min = reader.getAttribute<double>("lambd", root.lambda_min);
73 root.initial_dist = reader.getAttribute<dcomplex>("initial-range", root.initial_dist);
74 root.method = reader.enumAttribute<RootDigger::Method>("method")
75 .value("brent", RootDigger::ROOT_BRENT)
76 .value("broyden", RootDigger::ROOT_BROYDEN)
77 .value("muller", RootDigger::ROOT_MULLER)
78 .get(root.method);
79 reader.requireTagEnd();
80 }
81
82 public:
84 std::unique_ptr<Transfer> transfer;
85
87 void initTransfer(Expansion& expansion, bool reflection);
88
90 shared_ptr<OrderedAxis> vbounds;
91
93 shared_ptr<OrderedAxis> verts;
94
96 size_t lcount;
97
99 std::vector<bool> lgained;
100
102 std::vector<bool> lcomputed;
103
105 std::vector<std::size_t> stack;
106
108 std::ptrdiff_t interface;
109
111 double interface_position;
112
114 double lam0;
115
117 dcomplex k0;
118
120 PML vpml;
121
123 RootDigger::Params root;
124
126 bool recompute_integrals;
127
129 bool recompute_gain_integrals;
130
132 bool always_recompute_gain;
133
134 protected:
136 bool group_layers;
137
139 double max_temp_diff;
140
142 double temp_dist;
143
145 double temp_layer;
146
148 void scaleIncidentVector(cvector& incident, size_t layer, double size_factor);
149
150 public:
151 ModalBase()
152 : emission(EMISSION_UNSPECIFIED),
153 determinant_type(Transfer::DETERMINANT_EIGENVALUE),
154 transfer_method(Transfer::METHOD_AUTO),
155 interface(-1),
156 interface_position(NAN),
157 lam0(NAN),
158 k0(NAN),
159 vpml(dcomplex(1., -2.), 2.0, 10., 0),
160 recompute_integrals(true),
161 always_recompute_gain(true),
162 group_layers(true),
163 max_temp_diff(NAN),
164 temp_dist(0.5),
165 temp_layer(0.05) {}
166
167 virtual ~ModalBase() {}
168
170 double getLam0() const { return lam0; }
172 void setLam0(double lam) {
173 lam0 = lam;
174 if (!isnan(lam) && (isnan(real(k0)) || isnan(imag(k0)))) k0 = 2e3 * PI / lam;
175 }
177 void clearLam0() { lam0 = NAN; }
178
180 dcomplex getK0() const { return k0; }
181
183 void setK0(dcomplex k) {
184 k0 = k;
185 if (k0 == 0.) k0 = 1e-12;
186 }
187
189 dcomplex getLam() const { return 2e3 * PI / k0; }
190
192 void setLam(dcomplex lambda) { k0 = 2e3 * PI / lambda; }
193
195 void clearFields() {
196 if (transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
197 }
198
205 size_t getLayerFor(double& h) const {
206 size_t n = vbounds->findUpIndex(h + 1e-15);
207 if (n == 0)
208 h -= vbounds->at(0);
209 else
210 h -= vbounds->at(n - 1);
211 return n;
212 }
213
215 virtual std::string getId() const = 0;
216
218 virtual bool initCalculation() = 0;
219
222 const std::vector<std::size_t>& getStack() const { return stack; }
223
225 virtual void computeIntegrals() = 0;
226
228 virtual void clearModes() = 0;
229
234 virtual bool setExpansionDefaults(bool with_k0 = true) = 0;
235
237 void ensureInterface() {
238 if (interface == -1) throw BadInput(this->getId(), "no interface position set");
239 if (interface == 0 || interface >= std::ptrdiff_t(stack.size()))
240 throw BadInput(this->getId(), "wrong interface position {0} (min: 1, max: {1})", interface, stack.size() - 1);
241 }
242
244 virtual Expansion& getExpansion() = 0;
245
251 dvector getIncidentFluxes(const cvector& incident, Transfer::IncidentDirection side);
252
258 dvector getReflectedFluxes(const cvector& incident, Transfer::IncidentDirection side);
259
265 dvector getTransmittedFluxes(const cvector& incident, Transfer::IncidentDirection side);
266
272 cvector getReflectedCoefficients(const cvector& incident, Transfer::IncidentDirection side);
273
279 cvector getTransmittedCoefficients(const cvector& incident, Transfer::IncidentDirection side);
280
286 double getReflection(const cvector& incident, Transfer::IncidentDirection side) {
287 double R = 0.;
288 for (double r : getReflectedFluxes(incident, side)) R += r;
289 return R;
290 }
291
297 double getTransmission(const cvector& incident, Transfer::IncidentDirection side) {
298 double T = 0.;
299 for (double t : getTransmittedFluxes(incident, side)) T += t;
300 return T;
301 }
302
303#ifndef NDEBUG
304 void getMatrices(size_t layer, cmatrix& RE, cmatrix& RH);
305#endif
306};
307
308}}} // namespace plask::optical::modal
309
310#endif // PLASK__SOLVER_OPTICAL_MODAL_SOLVERBASE_H