PLaSK library

Loading...
Searching...
No Matches
freecarrier2d.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 "freecarrier2d.hpp"
15#include "fd.hpp"
16#include "gauss_matrix.hpp"
17
18 namespace plask { namespace gain { namespace freecarrier {
19
20 constexpr double DIFF_STEP = 0.001;
21
22 template <typename GeometryT>
23 FreeCarrierGainSolver2D<GeometryT>::FreeCarrierGainSolver2D(const std::string& name)
24 : FreeCarrierGainSolver<SolverWithMesh<GeometryT, MeshAxis>>(name) {}
25
26 struct Region2D {
27 size_t left, right, bottom, top;
28 size_t rowl, rowr;
29 Region2D()
30 : left(0),
31 right(0),
32 bottom(std::numeric_limits<size_t>::max()),
33 top(std::numeric_limits<size_t>::max()),
34 rowl(std::numeric_limits<size_t>::max()),
35 rowr(0) {}
36};
37
38 template <typename GeometryT> void FreeCarrierGainSolver2D<GeometryT>::detectActiveRegions() {
39 shared_ptr<RectangularMesh<2>> mesh = makeGeometryGrid(this->geometry->getChild());
40 shared_ptr<RectangularMesh<2>> points = mesh->getElementMesh();
41
42 std::vector<Region2D> regions;
43
44 for (size_t r = 0; r < points->vert()->size(); ++r) {
45 size_t prev = 0;
46 shared_ptr<Material> material;
47 for (size_t c = 0; c < points->tran()->size(); ++c) { // In the (possible) active region
48 auto point = points->at(c, r);
49 size_t num = 0;
50
51 auto roles = this->geometry->getRolesAt(point);
52 for (const auto& role : roles) {
53 if (role.substr(0, 6) == "active") {
54 if (num != 0) throw BadInput(this->getId(), "multiple 'active' roles specified");
55 if (role.size() == 6) {
56 num = 1;
57 } else {
58 try {
59 num = boost::lexical_cast<size_t>(role.substr(6)) + 1;
60 } catch (boost::bad_lexical_cast&) {
61 throw BadInput(this->getId(), "bad active region number in role '{0}'", role);
62 }
63 }
64 } else if (role == "substrate") {
65 if (this->explicitSubstrate)
66 this->writelog(LOG_WARNING, "Explicit substrate layer specified, role 'substrate' ignored");
67 else {
68 if (!this->substrateMaterial)
69 this->substrateMaterial = this->geometry->getMaterial(point);
70 else if (*this->substrateMaterial != *this->geometry->getMaterial(point))
71 throw Exception("{0}: Non-uniform substrate layer.", this->getId());
72 }
73 }
74 }
75 if (num == 0 && roles.find("QW") != roles.end())
76 throw Exception("{0}: All marked quantum wells must belong to marked active region.", this->getId());
77
78 if (num) { // here we are inside the active region
79 if (regions.size() >= num) {
80 if (!material)
81 material = this->geometry->getMaterial(points->at(c, r));
82 else if (*material != *this->geometry->getMaterial(points->at(c, r))) {
83 throw Exception("{0}: Active region {1} is laterally non-uniform", this->getId(), num - 1);
84 }
85 }
86 regions.resize(max(regions.size(), num));
87 auto& reg = regions[num - 1];
88 if (prev != num) { // this region starts in the current row
89 if (reg.top < r) {
90 throw Exception("{0}: Active region {1} is disjoint", this->getId(), num - 1);
91 }
92 if (reg.bottom >= r)
93 reg.bottom = r; // first row
94 else if (reg.rowr <= c)
95 throw Exception("{0}: Active region {1} is disjoint", this->getId(), num - 1);
96 reg.top = r + 1;
97 reg.rowl = c;
98 if (reg.left > reg.rowl) reg.left = reg.rowl;
99 }
100 }
101 if (prev && prev != num) { // previous region ended
102 auto& reg = regions[prev - 1];
103 if (reg.bottom < r && reg.rowl >= c) throw Exception("{0}: Active region {1} is disjoint", this->getId(), prev - 1);
104 reg.rowr = c;
105 if (reg.right < reg.rowr) reg.right = reg.rowr;
106 }
107 prev = num;
108 }
109 if (prev) // junction reached the edge
110 regions[prev - 1].rowr = regions[prev - 1].right = points->tran()->size();
111 }
112
113 this->regions.clear();
114 size_t act = 0;
115 for (auto& reg : regions) {
116 if (reg.bottom == std::numeric_limits<size_t>::max()) {
117 ++act;
118 continue;
119 }
120 if (reg.bottom == 0) //
121 throw Exception("{0}: Active region cannot be located at the bottom of the structure.", this->getId());
122 if (reg.top == points->axis[1]->size())
123 throw Exception("{0}: Active region cannot be located at the top of the structure.", this->getId());
124 this->regions.emplace_back(mesh->at(reg.left, reg.bottom - 1));
125 auto region = &this->regions.back();
126 region->bottom = mesh->axis[1]->at(reg.bottom) - mesh->axis[1]->at(reg.bottom - 1);
127 region->top = mesh->axis[1]->at(reg.top + 1) - mesh->axis[1]->at(reg.top);
128 double width = mesh->axis[0]->at(reg.right) - mesh->axis[0]->at(reg.left);
129 for (size_t r = reg.bottom - 1, j = 0; r <= reg.top; ++r, ++j) {
130 bool layerQW = false;
131 for (size_t c = reg.left; c < reg.right; ++c) {
132 shared_ptr<Material> material;
133 auto point = points->at(c, r);
134 double height = mesh->axis[1]->at(r + 1) - mesh->axis[1]->at(r);
135 auto roles = this->geometry->getRolesAt(point);
136 bool QW = roles.find("QW") != roles.end() /* || roles.find("QD") != roles.end() */;
137 if (c == reg.left) {
138 auto material = this->geometry->getMaterial(point);
139 size_t n = region->layers->getChildrenCount();
140 shared_ptr<Block<2>> last;
141 if (n > 0) {
142 last = static_pointer_cast<Block<2>>(
143 static_pointer_cast<Translation<2>>(region->layers->getChildNo(n - 1))->getChild());
144 assert(!last || last->size.c0 == width);
145 if (last && material == last->getRepresentativeMaterial() && QW == region->isQW(region->size() - 1)) {
146 // TODO check if usage of getRepresentativeMaterial is fine here (was material)
147 last->setSize(width, last->size.c1 + height);
148 continue;
149 }
150 }
151 auto layer = plask::make_shared<Block<2>>(Vec<2>(width, height), material);
152 if (QW) layer->addRole("QW");
153 region->layers->push_back(layer);
154 layerQW = QW;
155 } else if (layerQW != QW)
156 throw Exception("{}: Quantum wells in active region {} are not consistent", this->getId(), act);
157 }
158 }
159 ++act;
160 }
161
162 if (this->strained && !this->substrateMaterial)
163 throw BadInput(this->getId(), "strained quantum wells requested but no layer with substrate role set");
164
165 this->writelog(LOG_DETAIL, "Found {0} active region{1}", this->regions.size(), (this->regions.size() == 1) ? "" : "s");
166 for (auto& region : this->regions) region.summarize(this);
167}
168
169 static const shared_ptr<OrderedAxis> zero_axis(new OrderedAxis({0.}));
170
172 template <typename GeometryT> template <typename DT> struct FreeCarrierGainSolver2D<GeometryT>::DataBase : public LazyDataImpl<DT> {
173 struct AveragedData {
174 shared_ptr<const RectangularMesh<2>> mesh;
175 LazyData<double> data;
176 double factor;
177 const FreeCarrierGainSolver2D<GeometryT>* solver;
178 const char* name;
179
180 AveragedData(const FreeCarrierGainSolver2D<GeometryT>* solver,
181 const char* name,
182 const shared_ptr<const MeshAxis>& haxis,
183 const ActiveRegionInfo& region)
184 : solver(solver), name(name) {
185 auto vaxis = plask::make_shared<OrderedAxis>();
186 OrderedAxis::WarningOff vaxiswoff(vaxis);
187 for (size_t n = 0; n != region.size(); ++n) {
188 if (region.isQW(n)) {
189 auto box = region.getLayerBox(n);
190 vaxis->addPoint(0.5 * (box.lower.c1 + box.upper.c1));
191 }
192 }
193 mesh = plask::make_shared<const RectangularMesh<2>>(const_pointer_cast<MeshAxis>(haxis), vaxis,
194 RectangularMesh<2>::ORDER_01);
195 factor = 1. / double(vaxis->size());
196 }
197
198 size_t size() const { return mesh->axis[0]->size(); }
199
200 double operator[](size_t i) const {
201 double val = 0.;
202 for (size_t j = 0; j != mesh->axis[1]->size(); ++j) {
203 double v = data[mesh->index(i, j)];
204 if (isnan(v)) throw ComputationError(solver->getId(), "wrong {0} ({1}) at {2}", name, v, mesh->at(i, j));
205 val += v;
206 }
207 return val * factor;
208 }
209 };
210
211 typedef typename FreeCarrierGainSolver2D<GeometryT>::ActiveRegionParams ActiveRegionParams;
212
213 FreeCarrierGainSolver2D<GeometryT>* solver;
214 std::vector<shared_ptr<MeshAxis>> regpoints;
215 shared_ptr<const MeshD<2>> dest_mesh;
216 InterpolationFlags interpolation_flags;
217
218 void setupFromAxis(const shared_ptr<MeshAxis>& axis) {
219 regpoints.reserve(solver->regions.size());
220 InterpolationFlags flags(solver->geometry);
221 for (size_t r = 0; r != solver->regions.size(); ++r) {
222 std::set<double> pts;
223 auto box = solver->regions[r].getBoundingBox();
224 double y = 0.5 * (box.lower.c1 + box.upper.c1);
225 for (double x : *axis) {
226 auto p = flags.wrap(vec(x, y));
227 if (solver->regions[r].contains(p)) pts.insert(p.c0);
228 }
229 auto msh = plask::make_shared<OrderedAxis>();
230 OrderedAxis::WarningOff mshw(msh);
231 ;
232 msh->addOrderedPoints(pts.begin(), pts.end(), pts.size());
233 regpoints.emplace_back(std::move(msh));
234 }
235 }
236
237 DataBase(FreeCarrierGainSolver2D<GeometryT>* solver, const shared_ptr<const MeshD<2>>& dst_mesh)
238 : solver(solver), dest_mesh(dst_mesh), interpolation_flags(solver->geometry) {
239 // Create horizontal points lists
240 if (solver->mesh) {
241 setupFromAxis(solver->mesh);
242 } else if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh)) {
243 setupFromAxis(rect_mesh->axis[0]);
244 } else {
245 regpoints.reserve(solver->regions.size());
246 InterpolationFlags flags(solver->geometry);
247 for (size_t r = 0; r != solver->regions.size(); ++r) {
248 std::set<double> pts;
249 for (auto point : *dest_mesh) {
250 auto p = flags.wrap(point);
251 if (solver->regions[r].contains(p)) pts.insert(p.c0);
252 }
253 auto msh = plask::make_shared<OrderedAxis>();
254 OrderedAxis::WarningOff mshw(msh);
255 msh->addOrderedPoints(pts.begin(), pts.end(), pts.size());
256 regpoints.emplace_back(std::move(msh));
257 }
258 }
259 }
260
261 size_t size() const override { return dest_mesh->size(); }
262};
263
264 template <typename GeometryT>
265 struct FreeCarrierGainSolver2D<GeometryT>::ComputedData : public FreeCarrierGainSolver2D<GeometryT>::DataBaseTensor2 {
266 using typename DataBaseTensor2::AveragedData;
267
269 std::vector<LazyData<Tensor2<double>>> data;
270
271 template <typename... Args> ComputedData(Args... args) : DataBaseTensor2(args...) {}
272
273 void compute(double wavelength, InterpolationMethod interp) {
274 // Compute gains on mesh for each active region
275 OmpLockGuard lock(gain_omp_lock);
276 this->data.resize(this->solver->regions.size());
277 for (size_t reg = 0; reg != this->solver->regions.size(); ++reg) {
278 if (this->regpoints[reg]->size() == 0) {
279 this->data[reg] = LazyData<Tensor2<double>>(this->dest_mesh->size(), Tensor2<double>(0., 0.));
280 continue;
281 }
282 AveragedData temps(this->solver, "temperature", this->regpoints[reg], this->solver->regions[reg]);
283 AveragedData concs(temps);
284 concs.name = "carriers concentration";
285 temps.data = this->solver->inTemperature(temps.mesh, interp);
286 concs.data = this->solver->inCarriersConcentration(CarriersConcentration::PAIRS, temps.mesh, interp);
287 this->data[reg] =
288 interpolate(plask::make_shared<RectangularMesh<2>>(this->regpoints[reg], zero_axis),
289 getValues(wavelength, interp, reg, concs, temps), this->dest_mesh, interp, this->interpolation_flags);
290 }
291 }
292
293 virtual DataVector<Tensor2<double>> getValues(double wavelength,
294 InterpolationMethod interp,
295 size_t reg,
296 const AveragedData& concs,
297 const AveragedData& temps) = 0;
298
299 Tensor2<double> at(size_t i) const override {
300 for (size_t reg = 0; reg != this->solver->regions.size(); ++reg)
301 if (this->solver->regions[reg].inQW(this->interpolation_flags.wrap(this->dest_mesh->at(i)))) return this->data[reg][i];
302 return Tensor2<double>(0., 0.);
303 }
304};
305
306 template <typename GeometryT>
307 struct FreeCarrierGainSolver2D<GeometryT>::GainData : public FreeCarrierGainSolver2D<GeometryT>::ComputedData {
308 using typename DataBaseTensor2::AveragedData;
309
310 template <typename... Args> GainData(Args... args) : ComputedData(args...) {}
311
312 DataVector<Tensor2<double>> getValues(double wavelength,
313 InterpolationMethod interp,
314 size_t reg,
315 const AveragedData& concs,
316 const AveragedData& temps) override {
317 double hw = phys::h_eVc1e9 / wavelength;
318 DataVector<Tensor2<double>> values(this->regpoints[reg]->size());
319 std::exception_ptr error;
320
321 if (this->solver->inFermiLevels.hasProvider()) {
322 AveragedData Fcs(temps);
323 Fcs.name = "quasi Fermi level for electrons";
324 AveragedData Fvs(temps);
325 Fvs.name = "quasi Fermi level for holes";
326 Fcs.data = this->solver->inFermiLevels(FermiLevels::ELECTRONS, temps.mesh, interp);
327 Fvs.data = this->solver->inFermiLevels(FermiLevels::HOLES, temps.mesh, interp);
328 plask::openmp_size_t end = this->regpoints[reg]->size();
329 PLASK_OMP_PARALLEL_FOR
330 for (plask::openmp_size_t i = 0; i < end; ++i) {
331 if (error) continue;
332 try {
333 double T = temps[i];
334 double conc = max(concs[i], 1e-6); // To avoid hangs
335 double nr = this->solver->regions[reg].averageNr(wavelength, T, conc);
336 ActiveRegionParams params(this->solver, this->solver->params0[reg], T, bool(i));
337 values[i] = this->solver->getGain(hw, Fcs[i], Fvs[i], T, nr, params);
338 } catch (...) {
339 #pragma omp critical
340 error = std::current_exception();
341 }
342 }
343 if (error) std::rethrow_exception(error);
344 } else {
345 plask::openmp_size_t end = this->regpoints[reg]->size();
346 PLASK_OMP_PARALLEL_FOR
347 for (plask::openmp_size_t i = 0; i < end; ++i) {
348 if (error) continue;
349 try {
350 double T = temps[i];
351 double conc = max(concs[i], 1e-6); // To avoid hangs
352 double nr = this->solver->regions[reg].averageNr(wavelength, T, conc);
353 ActiveRegionParams params(this->solver, this->solver->params0[reg], T, bool(i));
354 double Fc = NAN, Fv = NAN;
355 this->solver->findFermiLevels(Fc, Fv, conc, T, params);
356 values[i] = this->solver->getGain(hw, Fc, Fv, T, nr, params);
357 } catch (...) {
358 #pragma omp critical
359 error = std::current_exception();
360 }
361 }
362 if (error) std::rethrow_exception(error);
363 }
364 return values;
365 }
366};
367
368 template <typename GeometryT>
369 struct FreeCarrierGainSolver2D<GeometryT>::DgdnData : public FreeCarrierGainSolver2D<GeometryT>::ComputedData {
370 using typename DataBaseTensor2::AveragedData;
371
372 template <typename... Args> DgdnData(Args... args) : ComputedData(args...) {}
373
374 DataVector<Tensor2<double>> getValues(double wavelength,
375 InterpolationMethod /*interp*/,
376 size_t reg,
377 const AveragedData& concs,
378 const AveragedData& temps) override {
379 double hw = phys::h_eVc1e9 / wavelength;
380 const double h = 0.5 * DIFF_STEP;
381 DataVector<Tensor2<double>> values(this->regpoints[reg]->size());
382 std::exception_ptr error;
383 plask::openmp_size_t end = this->regpoints[reg]->size();
384 PLASK_OMP_PARALLEL_FOR
385 for (plask::openmp_size_t i = 0; i < end; ++i) {
386 if (error) continue;
387 try {
388 double T = temps[i];
389 double conc = max(concs[i], 1e-6); // To avoid hangs
390 double nr = this->solver->regions[reg].averageNr(wavelength, T, conc);
391 ActiveRegionParams params(this->solver, this->solver->params0[reg], T, bool(i));
392 double Fc = NAN, Fv = NAN;
393 this->solver->findFermiLevels(Fc, Fv, (1. - h) * conc, T, params);
394 Tensor2<double> gain1 = this->solver->getGain(hw, Fc, Fv, T, nr, params);
395 this->solver->findFermiLevels(Fc, Fv, (1. + h) * conc, T, params);
396 Tensor2<double> gain2 = this->solver->getGain(hw, Fc, Fv, T, nr, params);
397 values[i] = (gain2 - gain1) / (2. * h * conc);
398 } catch (...) {
399 #pragma omp critical
400 error = std::current_exception();
401 }
402 }
403 if (error) std::rethrow_exception(error);
404 return values;
405 }
406};
407
408 template <typename GeometryT>
409 const LazyData<Tensor2<double>> FreeCarrierGainSolver2D<GeometryT>::getGainData(Gain::EnumType what,
410 const shared_ptr<const MeshD<2>>& dst_mesh,
411 double wavelength,
412 InterpolationMethod interp) {
413 if (what == Gain::GAIN) {
414 this->initCalculation(); // This must be called before any calculation!
415 this->writelog(LOG_DETAIL, "Calculating gain");
416 GainData* data = new GainData(this, dst_mesh);
417 data->compute(wavelength, getInterpolationMethod<INTERPOLATION_SPLINE>(interp));
418 return LazyData<Tensor2<double>>(data);
419 } else if (what == Gain::DGDN) {
420 this->initCalculation(); // This must be called before any calculation!
421 this->writelog(LOG_DETAIL, "Calculating gain over carriers concentration derivative");
422 DgdnData* data = new DgdnData(this, dst_mesh);
423 data->compute(wavelength, getInterpolationMethod<INTERPOLATION_SPLINE>(interp));
424 return LazyData<Tensor2<double>>(data);
425 } else {
426 throw BadInput(this->getId(), "wrong gain type requested");
427 }
428}
429
430 template <typename GeometryT>
431 struct FreeCarrierGainSolver2D<GeometryT>::EnergyLevelsData : public FreeCarrierGainSolver2D<GeometryT>::DataBaseVector {
432 using typename DataBaseVector::AveragedData;
433
434 size_t which;
435 std::vector<LazyData<double>> temps;
436
437 EnergyLevelsData(EnergyLevels::EnumType which,
438 FreeCarrierGainSolver2D<GeometryT>* solver,
439 const shared_ptr<const MeshD<2>>& dst_mesh,
440 InterpolationMethod interp)
441 : DataBaseVector(solver, dst_mesh), which(size_t(which)) {
442 temps.reserve(solver->regions.size());
443 for (size_t reg = 0; reg != solver->regions.size(); ++reg) {
444 AveragedData temp(this->solver, "temperature", this->regpoints[reg], this->solver->regions[reg]);
445 temp.data = this->solver->inTemperature(temp.mesh, interp);
446 temps.emplace_back(
447 interpolate(temp.mesh, DataVector<const double>(temp.data), this->dest_mesh, interp, this->solver->geometry));
448 }
449 }
450
451 std::vector<double> at(size_t i) const override {
452 for (size_t reg = 0; reg != this->solver->regions.size(); ++reg)
453 if (this->solver->regions[reg].contains(this->interpolation_flags.wrap(this->dest_mesh->at(i)))) {
454 double T = temps[reg][i];
455 ActiveRegionParams params(this->solver, this->solver->params0[reg], T, bool(i));
456 std::vector<double> result;
457 result.reserve(params.levels[which].size());
458 for (const auto& level : params.levels[which]) result.push_back(level.E);
459 return result;
460 }
461 return std::vector<double>();
462 }
463};
464
465 template <typename GeometryT>
466 const LazyData<std::vector<double>> FreeCarrierGainSolver2D<GeometryT>::getEnergyLevels(EnergyLevels::EnumType which,
467 const shared_ptr<const MeshD<2>>& dst_mesh,
468 InterpolationMethod interp) {
469 this->initCalculation();
470 EnergyLevelsData* data = new EnergyLevelsData(which, this, dst_mesh, getInterpolationMethod<INTERPOLATION_LINEAR>(interp));
471 return LazyData<std::vector<double>>(data);
472}
473
474 template <> std::string FreeCarrierGainSolver2D<Geometry2DCartesian>::getClassName() const { return "gain.FreeCarrier2D"; }
475 template <> std::string FreeCarrierGainSolver2D<Geometry2DCylindrical>::getClassName() const { return "gain.FreeCarrierCyl"; }
476
477 template struct PLASK_SOLVER_API FreeCarrierGainSolver2D<Geometry2DCartesian>;
478 template struct PLASK_SOLVER_API FreeCarrierGainSolver2D<Geometry2DCylindrical>;
479
480}}} // namespace plask::gain::freecarrier