PLaSK library

Loading...
Searching...
No Matches
freecarrier.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__GAIN_FREECARRIER_FREECARRIER_H
15#define PLASK__SOLVER__GAIN_FREECARRIER_FREECARRIER_H
16
17#include <plask/plask.hpp>
18
19#include <boost/math/tools/roots.hpp>
20 using boost::math::tools::toms748_solve;
21
22 namespace plask { namespace gain { namespace freecarrier {
23
24 template <typename BaseT> struct GainSpectrum;
25
26 extern OmpNestedLock gain_omp_lock;
27
31 template <typename BaseT> struct PLASK_SOLVER_API FreeCarrierGainSolver : public BaseT {
33 enum WhichLevel : size_t { EL = EnergyLevels::ELECTRONS, HH = EnergyLevels::HEAVY_HOLES, LH = EnergyLevels::LIGHT_HOLES };
34
35 typedef typename BaseT::SpaceType GeometryType;
36
37 enum { DIM = GeometryType::DIM };
38
39 struct ActiveRegionParams;
40
42 struct Level {
43 double E;
44 Tensor2<double> M;
45 double thickness;
46
47 Level(double E, const Tensor2<double>& M) : E(E), M(M) {}
48
49 Level(double E, const Tensor2<double>& M, double t) : E(E), M(M), thickness(t) {}
50
51 // TODO Better effective thickness computation (based on wavefunction)?
52 Level(double E, const Tensor2<double>& M, WhichLevel which, const ActiveRegionParams& params);
53 };
54
56 struct ActiveRegionInfo {
57 shared_ptr<StackContainer<DIM>> layers;
58 Vec<DIM> origin;
59
60 ActiveRegionInfo(Vec<DIM> origin) : layers(plask::make_shared<StackContainer<DIM>>()), origin(origin) {}
61
63 size_t size() const { return layers->getChildrenCount(); }
64
66 shared_ptr<Material> getLayerMaterial(size_t n) const {
67 auto block = static_cast<Block<DIM>*>(static_cast<Translation<DIM>*>(layers->getChildNo(n).get())->getChild().get());
68 if (auto m = block->singleMaterial()) return m;
69 throw plask::Exception("freeCarrierGainSolver requires solid layers.");
70 }
71
73 typename Primitive<DIM>::Box getLayerBox(size_t n) const {
74 return static_cast<GeometryObjectD<DIM>*>(layers->getChildNo(n).get())->getBoundingBox() + origin;
75 }
76
78 bool isQW(size_t n) const { return static_cast<Translation<DIM>*>(layers->getChildNo(n).get())->getChild()->hasRole("QW"); }
79
81 typename Primitive<DIM>::Box getBoundingBox() const { return layers->getBoundingBox() + origin; }
82
84
86 bool contains(const Vec<DIM>& point) const { return getBoundingBox().contains(point); }
87
89 bool inQW(const Vec<DIM>& point) const {
90 if (!contains(point)) return false;
91 assert(layers->getChildForHeight(point.vert() - origin.vert()));
92 return layers->getChildForHeight(point.vert() - origin.vert())->getChild()->hasRole("QW");
93 }
94
95 double averageNr(double lam, double T, double conc = 0.) const {
96 double nr = 0.;
97 for (size_t i = 0; i != materials.size(); ++i)
98 if (isQW(i)) nr += thicknesses[i] * materials[i]->Nr(lam, T, conc).real();
99 return nr / totalqw;
100 }
101
102 std::vector<shared_ptr<Material>> materials;
103 std::vector<double> thicknesses;
104 std::vector<size_t> wells;
105
106 double total;
107 double totalqw;
108 double bottom;
109 double top;
110
111 enum ConsideredHoles : unsigned {
112 NO_HOLES = 0,
113 HEAVY_HOLES = 1,
114 LIGHT_HOLES = 2,
115 BOTH_HOLES = 3
116 } holes;
117
122 void summarize(const FreeCarrierGainSolver<BaseT>* solver);
123 };
124
126 struct PLASK_SOLVER_API ActiveRegionParams {
127 const ActiveRegionInfo& region;
128 std::vector<double> U[3];
129 std::vector<Tensor2<double>> M[3];
130 double Mt;
131
132 std::vector<Level> levels[3];
133 double Eg;
134 size_t nhh,
135 nlh;
136
137 ActiveRegionParams(const FreeCarrierGainSolver* solver,
138 const ActiveRegionInfo& region,
139 double T,
140 bool quiet = false,
141 double mt = 0.);
142
143 ActiveRegionParams(const FreeCarrierGainSolver* solver, const ActiveRegionInfo& region, bool quiet = false, double mt = 0.)
144 : ActiveRegionParams(solver, region, solver->T0, quiet, mt) {}
145
146 explicit ActiveRegionParams(const FreeCarrierGainSolver* solver, const ActiveRegionParams& ref, double T, bool quiet = true)
147 : ActiveRegionParams(solver, ref.region, T, quiet, ref.Mt) {
148 nhh = ref.nhh;
149 nlh = ref.nlh;
150 for (size_t which = 0; which < 3; ++which) {
151 double shift = delta(WhichLevel(which), ref);
152 levels[which].reserve(ref.levels[which].size());
153 for (Level level : ref.levels[which]) levels[which].emplace_back(level.E + shift, level.M, level.thickness);
154 }
155 }
156
157 double sideU(WhichLevel which) const { return 0.5 * (U[which][0] + U[which][U[which].size() - 1]); }
158
159 Tensor2<double> sideM(WhichLevel which) const { return 0.5 * (M[which][0] + M[which][M[which].size() - 1]); }
160
161 double delta(WhichLevel which, const ActiveRegionParams& ref) const {
162 assert(U[which].size() == ref.U[which].size());
163 double delta = 0;
164 for (size_t i = 0; i < U[which].size(); ++i) {
165 delta += U[which][i] - ref.U[which][i];
166 }
167 return delta / double(U[which].size());
168 }
169 };
170
172 std::vector<ActiveRegionInfo> regions;
173
175 ReceiverFor<Temperature, GeometryType> inTemperature;
176
178 ReceiverFor<CarriersConcentration, GeometryType> inCarriersConcentration;
179
181 ReceiverFor<BandEdges, GeometryType> inBandEdges;
182
184 ReceiverFor<FermiLevels, GeometryType> inFermiLevels;
185
187 typename ProviderFor<Gain, GeometryType>::Delegate outGain;
188
190 typename ProviderFor<EnergyLevels, GeometryType>::Delegate outEnergyLevels;
191
192 FreeCarrierGainSolver(const std::string& name = "");
193
194 virtual ~FreeCarrierGainSolver();
195
196 void loadConfiguration(plask::XMLReader& reader, plask::Manager& manager) override;
197
198 protected:
200 shared_ptr<Material> substrateMaterial;
201
203 bool explicitSubstrate = false;
204
209 virtual void detectActiveRegions() = 0;
210
212 double level(WhichLevel which, double E, const ActiveRegionParams& params, size_t start, size_t stop) const;
213
214 double level(WhichLevel which, double E, const ActiveRegionParams& params) const {
215 return level(which, E, params, 0, params.region.materials.size() - 1);
216 }
217
218 double level(WhichLevel which, double E, const ActiveRegionParams& params, size_t well) const {
219 return level(which, E, params, params.region.wells[well], params.region.wells[well + 1]);
220 }
221
223 void estimateWellLevels(WhichLevel which, ActiveRegionParams& params, size_t qw) const;
224
226 void estimateAboveLevels(WhichLevel which, ActiveRegionParams& params) const;
227
229 template <class F>
230 std::pair<double, double> fermi_bracket_and_solve(F f, double guess, double step, boost::uintmax_t& max_iter) const {
231 double a = guess - 0.5 * step, b = guess + 0.5 * step;
232 double fa = f(a), fb = f(b);
233 if (fa == 0.) return std::make_pair(a, a);
234 if (fb == 0.) return std::make_pair(b, b);
235 boost::uintmax_t count = max_iter - 1;
236 if ((fa < 0.) == (fa < fb)) {
237 while ((fa < 0.) == (fb < 0.)) {
238 if (count == 0) return std::make_pair(a, b);
239 a = b;
240 fa = fb;
241 b += step;
242 fb = f(b);
243 if (fb == 0.) return std::make_pair(b, b);
244 --count;
245 }
246 } else {
247 while ((fb < 0.) == (fa < 0.)) {
248 if (count == 0) return std::make_pair(a, b);
249 b = a;
250 fb = fa;
251 a -= step;
252 fa = f(a);
253 if (fa == 0.) return std::make_pair(a, a);
254 --count;
255 }
256 }
257 auto res = toms748_solve(
258 f, a, b, fa, fb, [this](double l, double r) { return r - l < levelsep; }, count);
259 max_iter += count;
260 return res;
261 }
262
263#ifndef NDEBUG
264 public:
265#endif
267 double getN(double F, double T, const ActiveRegionParams& params) const;
268
270 double getP(double F, double T, const ActiveRegionParams& params) const;
271
272 protected:
273 double lifetime;
274 double matrixelem;
275
276 double T0;
277
278 protected:
279 double levelsep;
280
281 bool strained;
282
284 void estimateLevels();
285
287 void onInitialize() override;
288
290 void onInvalidate() override;
291
293 void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) {
294 outGain.fireChanged(); // the input changed, so we inform the world that everybody should get the new gain
295 }
296
305 virtual const LazyData<Tensor2<double>> getGainData(Gain::EnumType what,
306 const shared_ptr<const MeshD<DIM>>& dst_mesh,
307 double wavelength,
308 InterpolationMethod interp = INTERPOLATION_DEFAULT) = 0;
309
314 virtual const LazyData<std::vector<double>> getEnergyLevels(EnergyLevels::EnumType which,
315 const shared_ptr<const MeshD<DIM>>& dst_mesh,
316 InterpolationMethod interp = INTERPOLATION_DEFAULT) = 0;
317
318 public:
319 std::vector<ActiveRegionParams> params0;
320
321 bool quick_levels;
322
323#ifndef NDEBUG
324 double detEl(double E, const ActiveRegionParams& params, size_t well = 0) {
325 if (well)
326 return level(EL, E, params, well - 1);
327 else
328 return level(EL, E, params);
329 }
330 double detHh(double E, const ActiveRegionParams& params, size_t well = 0) {
331 if (well)
332 return level(HH, E, params, well - 1);
333 else
334 return level(HH, E, params);
335 }
336 double detLh(double E, const ActiveRegionParams& params, size_t well = 0) {
337 if (well)
338 return level(LH, E, params, well - 1);
339 else
340 return level(LH, E, params);
341 }
342#endif
343
345 shared_ptr<Material> getSubstrate() const { return substrateMaterial; }
346
348 void setSubstrate(shared_ptr<Material> material) {
349 bool invalid = substrateMaterial != material;
350 substrateMaterial = material;
351 explicitSubstrate = bool(material);
352 if (invalid) this->invalidate();
353 }
354
356 void findFermiLevels(double& Fc, double& Fv, double n, double T, const ActiveRegionParams& params) const;
357
359 Tensor2<double> getGain0(double hw, double Fc, double Fv, double T, double nr, const ActiveRegionParams& params) const;
360
362 Tensor2<double> getGain(double hw, double Fc, double Fv, double T, double nr, const ActiveRegionParams& params) const;
363
364 double getT0() const { return T0; }
365 void setT0(double T) {
366 T0 = T;
367 this->invalidate();
368 }
369
370 double getLifeTime() const { return lifetime; }
371 void setLifeTime(double iLifeTime) { lifetime = iLifeTime; }
372
373 double getMatrixElem() const { return matrixelem; }
374 void setMatrixElem(double iMatrixElem) {
375 matrixelem = iMatrixElem;
376 this->invalidate();
377 }
378
379 bool getStrained() const { return strained; }
380 void setStrained(bool value) {
381 strained = value;
382 this->invalidate();
383 }
384
385 friend struct GainSpectrum<BaseT>;
386
387 shared_ptr<GainSpectrum<BaseT>> getGainSpectrum(const Vec<DIM>& point) {
388 this->initCalculation();
389 return make_shared<GainSpectrum<BaseT>>(this, point);
390 }
391
392 typedef GainSpectrum<BaseT> GainSpectrumType;
393};
394
398 template <typename BaseT> struct GainSpectrum {
399 typedef typename FreeCarrierGainSolver<BaseT>::ActiveRegionParams ActiveRegionParams;
400
401 FreeCarrierGainSolver<BaseT>* solver;
402 plask::optional<Vec<BaseT::SpaceType::DIM>> point;
403
404 size_t reg;
405 double temp;
406 double conc;
407 double Fc, Fv;
409 std::unique_ptr<ActiveRegionParams> params;
410
411 GainSpectrum(FreeCarrierGainSolver<BaseT>* solver, const Vec<BaseT::SpaceType::DIM> point) : solver(solver), point(point) {
412 for (size_t i = 0; i != solver->regions.size(); ++i) {
413 if (solver->regions[i].contains(point)) {
414 reg = i;
415 solver->inTemperature.changedConnectMethod(this, &GainSpectrum::onChange);
416 solver->inCarriersConcentration.changedConnectMethod(this, &GainSpectrum::onChange);
417 temp = solver->inTemperature(plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(point))[0];
418 conc = solver->inCarriersConcentration(CarriersConcentration::PAIRS,
419 plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(point))[0];
420 updateParams();
421 return;
422 };
423 }
424 throw BadInput(solver->getId(), "point {0} does not belong to any active region", point);
425 }
426
427 GainSpectrum(FreeCarrierGainSolver<BaseT>* solver, double T, double n, size_t reg)
428 : solver(solver), temp(T), conc(n), reg(reg) {
429 updateParams();
430 }
431
432 void updateParams() {
433 params.reset(new ActiveRegionParams(solver, solver->params0[reg], temp));
434 Fc = Fv = NAN;
435 solver->findFermiLevels(Fc, Fv, conc, temp, *params);
436 }
437
438 void onChange(ReceiverBase&, ReceiverBase::ChangeReason) {
439 temp = solver->inTemperature(plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(*point))[0];
440 conc = solver->inCarriersConcentration(CarriersConcentration::PAIRS,
441 plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(*point))[0];
442 updateParams();
443 }
444
445 ~GainSpectrum() {
446 solver->inTemperature.changedDisconnectMethod(this, &GainSpectrum::onChange);
447 solver->inCarriersConcentration.changedDisconnectMethod(this, &GainSpectrum::onChange);
448 }
449
455 Tensor2<double> getGain(double wavelength) const {
456 double nr = solver->regions[reg].averageNr(wavelength, temp, conc);
457 return solver->getGain(phys::h_eVc1e9 / wavelength, Fc, Fv, temp, nr, *params);
458 }
459};
460
461}}} // namespace plask::gain::freecarrier
462
463#endif // PLASK__SOLVER__GAIN_FREECARRIER_FREECARRIER_H