PLaSK library

Loading...
Searching...
No Matches
expansion.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_EXPANSION_H
15#define PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_H
16
17#include <plask/plask.hpp>
18
19#ifdef OPENMP_FOUND
20# include <omp.h>
21#endif
22
23#include "solverbase.hpp"
24#include "matrices.hpp"
25#include "meshadapter.hpp"
26#include "temp_matrix.hpp"
27
28 namespace plask { namespace optical { namespace modal {
29
30 struct PLASK_SOLVER_API Expansion {
31
33 enum Component {
34 E_UNSPECIFIED = 0,
35 E_TRAN = 1,
36 E_LONG = 2
37 };
38
39 WhichField which_field;
40 InterpolationMethod field_interpolation;
41
43 ModalBase* solver;
44
46 dcomplex k0;
47
49 double lam0;
50
52 LazyData<double> temperature;
53
55 bool gain_connected;
56
58 bool epsilon_connected;
59
61 LazyData<Tensor2<double>> gain;
62
64 LazyData<Tensor3<dcomplex>> epsilons;
65
67 LazyData<double> carriers;
68
69 Expansion(ModalBase* solver): solver(solver), k0(NAN), lam0(NAN) {}
70
71 virtual ~Expansion() {}
72
73 TempMatrix getTempMatrix() {
74 size_t N = matrixSize();
75 return temporary.get(N, N);
76 }
77
78 private:
79 dcomplex glambda;
80
81 protected:
82
83 TempMatrixPool temporary;
84
88 virtual void beforeLayersIntegrals(dcomplex PLASK_UNUSED(lam), dcomplex PLASK_UNUSED(glam)) {}
89
93 virtual void afterLayersIntegrals() {
94 temperature.reset();
95 gain.reset();
96 carriers.reset();
97 }
98
105 virtual void layerIntegrals(size_t layer, double lam, double glam) = 0;
106
107 public:
108
110 virtual void beforeGetEpsilon() {
111 computeIntegrals();
112 }
113
115 virtual void afterGetEpsilon() {}
116
118 double getLam0() const { return lam0; }
120 void setLam0(double lam) {
121 if (lam != lam0 && !(isnan(lam0) && isnan(lam))) {
122 lam0 = lam;
123 solver->recompute_integrals = true;
124 solver->clearFields();
125 }
126 }
128 void clearLam0() {
129 if (!isnan(lam0)) {
130 lam0 = NAN;
131 solver->recompute_integrals = true;
132 solver->clearFields();
133 }
134 }
135
137 dcomplex getK0() const { return k0; }
139 void setK0(dcomplex k) {
140 if (k != k0) {
141 k0 = k;
142 if (k0 == 0.) k0 = 1e-12;
143 if (isnan(lam0)) solver->recompute_integrals = true;
144 solver->clearFields();
145 }
146 }
147
149 void computeIntegrals() {
150 dcomplex lambda = 2e3*PI/k0;
151 if (solver->recompute_integrals) {
152 dcomplex lam;
153 if (!isnan(lam0)) {
154 lam = lam0;
155 glambda = (solver->always_recompute_gain)? lambda : lam;
156 } else{
157 lam = glambda = lambda;
158 }
159 size_t nlayers = solver->lcount;
160 std::exception_ptr error;
161 beforeLayersIntegrals(lam, glambda);
162 PLASK_OMP_PARALLEL_FOR
163 for (plask::openmp_size_t l = 0; l < nlayers; ++l) {
164 if (error) continue;
165 try {
166 layerIntegrals(l, real(lam), real(glambda));
167 } catch(...) {
168 #pragma omp critical
169 error = std::current_exception();
170 }
171 }
172 afterLayersIntegrals();
173 if (error) std::rethrow_exception(error);
174 solver->recompute_integrals = false;
175 solver->recompute_gain_integrals = false;
176 } else if (solver->recompute_gain_integrals ||
177 (solver->always_recompute_gain && !is_zero(lambda - glambda))) {
178 dcomplex lam = isnan(lam0)? lambda : solver->lam0;
179 glambda = solver->always_recompute_gain? lambda : lam;
180 std::vector<size_t> recomputed_layers;
181 size_t nlayers = solver->lcount;
182 recomputed_layers.reserve(nlayers);
183 for (size_t l = 0; l != nlayers; ++l) if (solver->lgained[l] || solver->lcomputed[l])
184 recomputed_layers.push_back(l);
185 std::exception_ptr error;
186 beforeLayersIntegrals(lam, glambda);
187 PLASK_OMP_PARALLEL_FOR
188 for (plask::openmp_size_t l = 0; l < recomputed_layers.size(); ++l) {
189 if (error) continue;
190 try {
191 layerIntegrals(recomputed_layers[l], real(lam), real(glambda));
192 } catch(...) {
193 #pragma omp critical
194 error = std::current_exception();
195 }
196 }
197 afterLayersIntegrals();
198 if (error) std::rethrow_exception(error);
199 solver->recompute_gain_integrals = false;
200 }
201 }
202
208 virtual bool diagonalQE(size_t PLASK_UNUSED(l)) const { return false; }
209
214 virtual size_t matrixSize() const = 0;
215
221 virtual void getMatrices(size_t layer, cmatrix& RE, cmatrix& RH) = 0;
222
230 virtual LazyData<Tensor3<dcomplex>> getMaterialEps(size_t lay,
231 const shared_ptr<const typename LevelsAdapter::Level> &level,
232 InterpolationMethod interp) = 0;
240 virtual void getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal& gamma);
241
247 void initField(WhichField which, InterpolationMethod method) {
248 which_field = which;
249 field_interpolation = method;
250 prepareField();
251 }
252
261 virtual double integrateField(WhichField field, size_t layer, const cmatrix& TE, const cmatrix& TH,
262 const std::function<std::pair<dcomplex,dcomplex>(size_t, size_t)>& vertical) = 0;
263
270 virtual double integratePoyntingVert(const cvector& E, const cvector& H) = 0;
271
279 inline double getModeFlux(size_t n, const cmatrix& TE, const cmatrix& TH) {
280 const cvector E(const_cast<dcomplex*>(TE.data()) + n*TE.rows(), TE.rows()),
281 H(const_cast<dcomplex*>(TH.data()) + n*TH.rows(), TH.rows());
282 return integratePoyntingVert(E, H);
283 }
284
285 protected:
286
290 virtual void prepareField() {}
291
292 public:
296 virtual void cleanupField() {}
297
306 virtual LazyData<Vec<3,dcomplex>> getField(size_t l,
307 const shared_ptr<const typename LevelsAdapter::Level>& level,
308 const cvector& E,
309 const cvector& H) = 0;
310};
311
312
313}}} // namespace plask
314
315#endif // PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_H