PLaSK library

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