PLaSK library

Loading...
Searching...
No Matches
diffusion3d.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 "diffusion3d.hpp"
15
16#define DEFAULT_MESH_SPACING 0.01 // µm
17
18 namespace plask {
19
20 namespace electrical { namespace diffusion {
21
22 constexpr double inv_hc = 1.0e-13 / (phys::c * phys::h_J);
23 using phys::Z0;
24
25 Diffusion3DSolver::Diffusion3DSolver(const std::string& name)
26 : FemSolverWithMesh<Geometry3D, RectangularMesh<3>>(name),
27 loopno(0),
28 maxerr(0.05),
29 outCarriersConcentration(this, &Diffusion3DSolver::getConcentration) {
30 this->algorithm = ALGORITHM_ITERATIVE;
31 inTemperature = 300.;
32}
33
34 void Diffusion3DSolver::loadConfiguration(XMLReader& source, Manager& manager) {
35 while (source.requireTagOrEnd()) parseConfiguration(source, manager);
36}
37
38 void Diffusion3DSolver::parseConfiguration(XMLReader& source, Manager& manager) {
39 std::string param = source.getNodeName();
40
41 if (param == "loop") {
42 maxerr = source.getAttribute<double>("maxerr", maxerr);
43 source.requireTagEnd();
44 }
45
46 else if (!this->parseFemConfiguration(source, manager)) {
47 this->parseStandardConfiguration(source, manager);
48 }
49}
50
51 Diffusion3DSolver::~Diffusion3DSolver() {}
52
53 void Diffusion3DSolver::summarizeActiveRegion(std::map<size_t, ActiveRegion3D::Region>& regions,
54 size_t num,
55 size_t bottom,
56 size_t top,
57 size_t lon,
58 size_t tra,
59 const std::vector<bool>& isQW,
60 const shared_ptr<RectangularMesh3D::ElementMesh>& points) {
61 if (!num) return;
62
63 auto found = regions.find(num);
64 ActiveRegion3D::Region& region = (found == regions.end()) ? regions
65 .emplace(std::piecewise_construct, std::forward_as_tuple(num),
66 std::forward_as_tuple(bottom, top, lon, tra, isQW))
67 .first->second
68 : found->second;
69
70 if (bottom != region.bottom || top != region.top)
71 throw Exception("{0}: Active region {1} does not have top and bottom edges at constant heights", this->getId(), num - 1);
72
73 shared_ptr<Material> material;
74 for (size_t i = region.bottom; i < region.top; ++i) {
75 bool QW = isQW[i];
76 if (QW != region.isQW[i])
77 throw Exception("{0}: Active region {1} does not have QWs at constant heights", this->getId(), num - 1);
78 if (QW) {
79 auto point = points->at(lon, tra, i);
80 if (!material)
81 material = this->geometry->getMaterial(point);
82 else if (*material != *this->geometry->getMaterial(point)) {
83 throw Exception("{}: Quantum wells in active region {} are not identical", this->getId(), num - 1);
84 }
85 }
86 }
87}
88
89 void Diffusion3DSolver::setupActiveRegions() {
90 if (!this->geometry || !this->mesh) return;
91
92 auto points = this->mesh->getElementMesh();
93
94 std::map<size_t, ActiveRegion3D::Region> regions;
95 std::vector<bool> isQW(points->axis[2]->size());
96
97 for (size_t lon = 0; lon < points->axis[0]->size(); ++lon) {
98 for (size_t tra = 0; tra < points->axis[1]->size(); ++tra) {
99 std::fill(isQW.begin(), isQW.end(), false);
100 size_t num = 0;
101 size_t start = 0;
102 for (size_t ver = 0; ver < points->axis[2]->size(); ++ver) {
103 auto point = points->at(lon, tra, ver);
104 auto roles = this->geometry->getRolesAt(point);
105 size_t cur = 0;
106 for (auto role : roles) { // find the active region, point belongs to
107 size_t l = 0;
108 if (role.substr(0, 6) == "active")
109 l = 6;
110 else if (role.substr(0, 8) == "junction")
111 l = 8;
112 else
113 continue;
114 if (cur != 0) throw BadInput(this->getId(), "multiple 'active'/'junction' roles specified");
115 if (role.size() == l)
116 cur = 1;
117 else {
118 try {
119 cur = boost::lexical_cast<size_t>(role.substr(l)) + 1;
120 } catch (boost::bad_lexical_cast&) {
121 throw BadInput(this->getId(), "bad active region number in role '{0}'", role);
122 }
123 }
124 }
125 if (cur != num) {
126 summarizeActiveRegion(regions, num, start, ver, lon, tra, isQW, points);
127 num = cur;
128 start = ver;
129 }
130 if (cur) {
131 isQW[ver] =
132 roles.find("QW") != roles.end() || roles.find("QD") != roles.end() || roles.find("carriers") != roles.end();
133 auto found = regions.find(cur);
134 if (found != regions.end()) {
135 ActiveRegion3D::Region& region = found->second;
136 if (region.warn && lon != region.lon && tra != region.tra &&
137 *this->geometry->getMaterial(points->at(lon, tra, ver)) !=
138 *this->geometry->getMaterial(points->at(region.lon, region.tra, ver))) {
139 writelog(LOG_WARNING, "Active region {} is laterally non-uniform", num - 1);
140 region.warn = false;
141 }
142 }
143 }
144 }
145 // Summarize the last region
146 summarizeActiveRegion(regions, num, start, points->axis[2]->size(), lon, tra, isQW, points);
147 }
148 }
149
150 active.clear();
151 for (auto& ireg : regions) {
152 size_t act = ireg.first - 1;
153 ActiveRegion3D::Region& reg = ireg.second;
154 // Detect quantum wells in the active region
155 std::vector<double> QWz;
156 std::vector<std::pair<size_t, size_t>> QWbt;
157 double QWheight = 0.;
158 for (size_t i = reg.bottom; i < reg.top; ++i) {
159 if (reg.isQW[i]) {
160 if (QWbt.empty() || QWbt.back().second != i)
161 QWbt.emplace_back(i, i + 1);
162 else
163 QWbt.back().second = i + 1;
164 QWz.push_back(points->axis[2]->at(i));
165 QWheight += this->mesh->vert()->at(i + 1) - this->mesh->vert()->at(i);
166 }
167 }
168 if (QWz.empty()) {
169 throw Exception("{}: Active region {} does not contain quantum wells", this->getId(), act);
170 }
171
172 active.emplace(std::piecewise_construct, std::forward_as_tuple(act),
173 std::forward_as_tuple(this, reg.bottom, reg.top, QWheight, std::move(QWz), std::move(QWbt)));
174
175 this->writelog(LOG_DETAIL, "Total QWs thickness in active region {}: {}nm", act, 1e3 * QWheight);
176 this->writelog(LOG_DEBUG, "Active region {0} vertical span: {1}-{2}", act, reg.bottom, reg.top);
177 }
178}
179
180 void Diffusion3DSolver::onInitialize() {
181 if (!this->geometry) throw NoGeometryException(this->getId());
182 if (!this->mesh) {
183 auto mesh1 = makeGeometryGrid(this->geometry);
184 this->mesh = make_shared<RectangularMesh<3>>(refineAxis(mesh1->lon(), DEFAULT_MESH_SPACING),
185 refineAxis(mesh1->tran(), DEFAULT_MESH_SPACING), mesh1->vert());
186 writelog(LOG_DETAIL, "{}: Setting up default mesh [{}]", this->getId(), this->mesh->tran()->size());
187 }
188 setupActiveRegions();
189 loopno = 0;
190}
191
192 void Diffusion3DSolver::onInvalidate() { active.clear(); }
193
194// clang-format off
195 inline void Diffusion3DSolver::setLocalMatrix(FemMatrix& K, DataVector<double>& F, const ElementParams3D e,
196 const double A, const double B, const double C, const double D, const double* U, const double* J) {
197# include "diffusion3d-eval.ipp"
198}
199
200 inline void Diffusion3DSolver::addLocalBurningMatrix(FemMatrix& K, DataVector<double>& F, const ElementParams3D e,
201 const Tensor2<double> G, const Tensor2<double> dG, const double Ug, const Tensor2<double>* P) {
202# include "diffusion3d-eval-shb.ipp"
203}
204// clang-format on
205
206 double Diffusion3DSolver::compute(unsigned loops, bool shb, size_t act) {
207 this->initCalculation();
208
209 auto found = this->active.find(act);
210 if (found == this->active.end()) throw Exception("{}: Active region {} does not exist", this->getId(), act);
211 auto& active = found->second;
212
213 double z = active.vert();
214
215 size_t nn = active.mesh2->size(), ne = active.emesh2->size();
216 size_t N = 3 * nn;
217 size_t nm = active.mesh2->lateralMesh->fullMesh.minorAxis()->size();
218
219 size_t nmodes = 0;
220
221 if (!active.U) active.U.reset(N, 0.);
222
223 DataVector<double> A(ne), B(ne), C(ne), D(ne);
224
225 auto temperature = active.verticallyAverage(inTemperature, active.emesh3, InterpolationMethod::INTERPOLATION_SPLINE);
226 for (size_t i = 0; i != ne; ++i) {
227 auto material = this->geometry->getMaterial(active.emesh2->at(i));
228 double T = temperature[i];
229 A[i] = material->A(T);
230 B[i] = material->B(T);
231 C[i] = material->C(T);
232 D[i] = 1e8 * material->D(T); // cm²/s -> µm²/s
233 }
234
235 DataVector<double> J(nn);
236 double js = 1e7 / (phys::qe * active.QWheight);
237 size_t i = 0;
238 for (auto j : inCurrentDensity(active.mesh2, InterpolationMethod::INTERPOLATION_SPLINE)) {
239 J[i] = abs(js * j.c2);
240 ++i;
241 }
242
243 std::vector<DataVector<Tensor2<double>>> Ps;
244 std::vector<DataVector<double>> nrs;
245
246 this->writelog(LOG_INFO, "Running diffusion calculations");
247
248 if (shb) {
249 nmodes = inWavelength.size();
250
251 if (inLightE.size() != nmodes)
252 throw BadInput(this->getId(), "number of modes in inWavelength ({}) and inLightE ({}) differ", inWavelength.size(),
253 inLightE.size());
254
255 active.modesP.assign(inWavelength.size(), 0.);
256
257 Ps.reserve(nmodes);
258 nrs.reserve(nmodes);
259 for (size_t i = 0; i != nmodes; ++i) {
260 Ps.emplace_back(nn);
261 nrs.emplace_back(ne);
262 }
263
264 for (size_t i = 0; i != ne; ++i) {
265 auto material = this->geometry->getMaterial(active.emesh2->at(i));
266 for (size_t n = 0; n != nmodes; ++n) nrs[n][i] = material->Nr(real(inWavelength(n)), temperature[i]).real();
267 }
268
269 for (size_t n = 0; n != nmodes; ++n) {
270 double wavelength = real(inWavelength(n));
271 writelog(LOG_DEBUG, "Mode {} wavelength: {} nm", n, wavelength);
272 auto P = active.verticallyAverage(inLightE, active.mesh3);
273 for (size_t i = 0; i != nn; ++i) {
274 Ps[n][i].c00 = (0.5 / Z0) * real(P[i].c0 * conj(P[i].c0) + P[i].c1 * conj(P[i].c1));
275 Ps[n][i].c11 = (0.5 / Z0) * real(P[i].c2 * conj(P[i].c2));
276 }
277 }
278 }
279
280 unsigned loop = 0;
281
282 std::unique_ptr<FemMatrix> K;
283
284 toterr = 0.;
285
286 DataVector<double> F(N);
287 DataVector<double> resid(N);
288
289 switch (this->algorithm) {
290 case ALGORITHM_CHOLESKY: K.reset(new DpbMatrix(this, N, 3 * nm + 5)); break;
291 case ALGORITHM_GAUSS: K.reset(new DgbMatrix(this, N, 3 * nm + 5)); break;
292 case ALGORITHM_ITERATIVE: K.reset(new SparseFreeMatrix(this, N, 78 * ne)); break;
293 }
294
295 while (true) {
296 // Set stiffness matrix and load vector
297 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", K->describe());
298 K->clear();
299 F.fill(0.);
300
301 for (size_t ie = 0; ie < ne; ++ie)
302 setLocalMatrix(*K, F, ElementParams3D(active, ie), A[ie], B[ie], C[ie], D[ie], active.U.data(), J.data());
303
304 write_debug("{}: Iteration {}", this->getId(), loop);
305
306 // Add SHB
307 if (shb) {
308 std::fill(active.modesP.begin(), active.modesP.end(), 0.);
309 for (size_t n = 0; n != nmodes; ++n) {
310 double wavelength = real(inWavelength(n));
311 double factor = inv_hc * wavelength;
312 auto gain = inGain(active.emesh2, wavelength, InterpolationMethod::INTERPOLATION_SPLINE);
313 auto dgdn = inGain(Gain::DGDN, active.emesh2, wavelength, InterpolationMethod::INTERPOLATION_SPLINE);
314 const Tensor2<double>* Pdata = Ps[n].data();
315 for (size_t ie = 0; ie < ne; ++ie) {
316 ElementParams3D el(active, ie);
317
318 Tensor2<double> g = nrs[n][ie] * gain[ie];
319 Tensor2<double> dg = nrs[n][ie] * dgdn[ie];
320 Tensor2<double> p = integrateBilinear(el.X, el.Y, Pdata + ie);
321 active.modesP[n] += p.c00 * g.c00 + p.c11 * g.c11;
322 g *= factor;
323 dg *= factor;
324 double ug =
325 0.25 * (active.U[el.i00] + active.U[el.i02] + active.U[el.i20] + active.U[el.i22] +
326 0.25 * (el.X * (active.U[el.i01] - active.U[el.i03] + active.U[el.i21] - active.U[el.i23]) +
327 el.Y * (active.U[el.i10] + active.U[el.i12] - active.U[el.i30] - active.U[el.i32])));
328 addLocalBurningMatrix(*K, F, el, g, dg, ug, Pdata);
329 }
330 active.modesP[n] *= 1e-13 * active.QWheight;
331 // 10⁻¹³ from µm to cm conversion and conversion to mW (r dr), (...) — photon energy
332 writelog(LOG_DEBUG, "{}: Mode {} burned power: {} mW", this->getId(), n, active.modesP[n]);
333 }
334 }
335
336 // // Set derivatives to 0 at the edges
337 // size_t nl = active.emesh2->lon()->size(), nt = active.emesh2->tran()->size();
338 // for (size_t i = 0; i <= nt; ++i) {
339 // K->setBC(F, active.mesh2->index(0, i, 0) + 1, 0.);
340 // K->setBC(F, active.mesh2->index(nl, i, 0) + 1, 0.);
341 // }
342 // for (size_t i = 0; i <= nl; ++i) {
343 // K->setBC(F, active.mesh2->index(i, 0, 0) + 2, 0.);
344 // K->setBC(F, active.mesh2->index(i, nt, 0) + 2, 0.);
345 // }
346
347#ifndef NDEBUG
348 double* kend = K->data + K->size;
349 for (double* pk = K->data; pk != kend; ++pk) {
350 if (isnan(*pk) || isinf(*pk))
351 throw ComputationError(this->getId(), "error in stiffness matrix at position {0} ({1})", pk - K->data,
352 isnan(*pk) ? "nan" : "inf");
353 }
354 for (auto f = F.begin(); f != F.end(); ++f) {
355 if (isnan(*f) || isinf(*f))
356 throw ComputationError(this->getId(), "error in load vector at position {0} ({1})", f - F.begin(),
357 isnan(*f) ? "nan" : "inf");
358 }
359#endif
360
361 // Compute current error
362 for (auto f = F.begin(), r = resid.begin(); f != F.end(); ++f, ++r) *r = -*f;
363 K->addmult(active.U, resid);
364
365 double err = 0.;
366 for (auto r = resid.begin(); r != resid.end(); ++r) err += *r * *r;
367 double denorm = 0.;
368 for (auto f = F.begin(); f != F.end(); ++f) denorm += *f * *f;
369 err = 100. * sqrt(err / denorm);
370
371 // Do next calculation step
372 if (loop != 0) this->writelog(LOG_RESULT, "Loop {:d}({:d}) @ active region {}: error = {:g}%", loop, loopno, act, err);
373 ++loopno;
374 ++loop;
375 if (err < maxerr || ((loops != 0 && loop >= loops))) break;
376
377 // TODO add linear mixing with the previous solution
378 K->solve(F, active.U);
379 }
380
381 outCarriersConcentration.fireChanged();
382
383 return toterr;
384}
385
386 double Diffusion3DSolver::get_burning_integral_for_mode(size_t mode) const {
387 if (mode >= inLightE.size()) throw BadInput(this->getId(), "mode index out of range");
388 size_t i = 0;
389 double res = 0.;
390 for (const auto& iactive : this->active) {
391 const auto& active = iactive.second;
392 if (mode >= active.modesP.size()) throw Exception("{}: SHB not computed for active region {}", this->getId(), i);
393 res += active.modesP[mode];
394 ++i;
395 }
396 return res;
397}
398
399 double Diffusion3DSolver::get_burning_integral() const {
400 double res = 0.;
401 for (size_t i = 0; i != inLightE.size(); ++i) res += get_burning_integral_for_mode(i);
402 return res;
403}
404
405 Diffusion3DSolver::ConcentrationDataImpl::ConcentrationDataImpl(const Diffusion3DSolver* solver,
406 shared_ptr<const plask::MeshD<3>> dest_mesh,
407 InterpolationMethod interp)
408 : solver(solver), destination_mesh(dest_mesh), interpolationFlags(InterpolationFlags(solver->geometry)) {
409 concentrations.reserve(solver->active.size());
410
411 if (interp == InterpolationMethod::INTERPOLATION_DEFAULT || interp == InterpolationMethod::INTERPOLATION_SPLINE) {
412 for (const auto& iactive : solver->active) {
413 const auto& active = iactive.second;
414 if (!active.U) throw NoValue("carriers concentration");
415 concentrations.emplace_back(LazyData<double>(dest_mesh->size(), [this, active](size_t i) -> double {
416 auto point = destination_mesh->at(i);
417
418 Vec<2> wrapped_point;
419 size_t index0_lo, index0_hi, index1_lo, index1_hi;
420
421 if (!active.mesh2->lateralMesh->prepareInterpolation(vec(point.c0, point.c1), wrapped_point, index0_lo, index0_hi,
422 index1_lo, index1_hi, interpolationFlags))
423 return 0.; // point is outside the active region
424
425 double x = wrapped_point.c0 - active.mesh2->lateralMesh->fullMesh.getAxis0()->at(index0_lo);
426 double y = wrapped_point.c1 - active.mesh2->lateralMesh->fullMesh.getAxis1()->at(index1_lo);
427
428 ElementParams3D e(active, active.emesh2->lateralMesh->index(index0_lo, index1_lo));
429
430 const double x2 = x * x, y2 = y * y;
431 const double x3 = x2 * x, y3 = y2 * y;
432 const double X2 = e.X * e.X, Y2 = e.Y * e.Y;
433 const double X3 = X2 * e.X, Y3 = Y2 * e.Y;
434
435 return (e.X * x *
436 (-x * y2 * (e.X - x) * (3 * e.Y - 2 * y) * active.U[e.i32] -
437 x * (e.X - x) * (Y3 - 3 * e.Y * y2 + 2 * y3) * active.U[e.i30] +
438 y2 * (3 * e.Y - 2 * y) * (X2 - 2 * e.X * x + x2) * active.U[e.i12] +
439 (X2 - 2 * e.X * x + x2) * (Y3 - 3 * e.Y * y2 + 2 * y3) *
440 active.U[e.i10]) +
441 e.Y * y *
442 (-x2 * y * (3 * e.X - 2 * x) * (e.Y - y) * active.U[e.i23] +
443 x2 * (3 * e.X - 2 * x) * (Y2 - 2 * e.Y * y + y2) * active.U[e.i21] -
444 y * (e.Y - y) * (X3 - 3 * e.X * x2 + 2 * x3) * active.U[e.i03] +
445 (X3 - 3 * e.X * x2 + 2 * x3) * (Y2 - 2 * e.Y * y + y2) *
446 active.U[e.i01]) +
447 x2 * y2 * (3 * e.X - 2 * x) * (3 * e.Y - 2 * y) * active.U[e.i22] +
448 x2 * (3 * e.X - 2 * x) * (Y3 - 3 * e.Y * y2 + 2 * y3) * active.U[e.i20] +
449 y2 * (3 * e.Y - 2 * y) * (X3 - 3 * e.X * x2 + 2 * x3) * active.U[e.i02] +
450 (X3 - 3 * e.X * x2 + 2 * x3) *
451 (Y3 - 3 * e.Y * y2 + 2 * y3) * active.U[e.i00]) /
452 (X3 * Y3);
453 }));
454 }
455 } else {
456 for (const auto& iactive : solver->active) {
457 const auto& active = iactive.second;
458 if (!active.U) throw NoValue("carriers concentration");
459 DataVector<double> conc(active.U.size() / 3);
460 DataVector<double>::iterator c = conc.begin();
461 for (auto u = active.U.begin(); u < active.U.end(); u += 3, ++c) *c = *u;
462 assert(active.mesh2->size() == conc.size());
463 LazyData<double> interpolation = interpolate(active.mesh2, conc, dest_mesh, interp, interpolationFlags);
464 concentrations.emplace_back(LazyData<double>(dest_mesh->size(), [interpolation](size_t i) -> double {
465 auto res = interpolation.at(i);
466 if (isnan(res)) return 0.;
467 return res;
468 }));
469 }
470 }
471}
472
473 double Diffusion3DSolver::ConcentrationDataImpl::at(size_t i) const {
474 auto point = interpolationFlags.wrap(destination_mesh->at(i));
475 bool found = false;
476 size_t an;
477 for (const auto& iactive : solver->active) {
478 const auto& active = iactive.second;
479 if (solver->mesh->vert()->at(active.bottom) <= point.c2 && point.c2 <= solver->mesh->vert()->at(active.top)) {
480 // Make sure we have concentration only in the quantum wells
481 // TODO maybe more optimal approach would be reasonable?
482 for (auto qw : active.QWs)
483 if (qw.first <= point.c2 && point.c2 < qw.second) {
484 found = true;
485 an = iactive.first;
486 break;
487 }
488 break;
489 }
490 }
491 if (!found) return 0.;
492 return concentrations[an][i];
493}
494
495 const LazyData<double> Diffusion3DSolver::getConcentration(CarriersConcentration::EnumType what,
496 shared_ptr<const plask::MeshD<3>> dest_mesh,
497 InterpolationMethod interpolation) const {
498 if (what != CarriersConcentration::MAJORITY && what != CarriersConcentration::PAIRS) {
499 return LazyData<double>(dest_mesh->size(), NAN);
500 }
501 return LazyData<double>(new Diffusion3DSolver::ConcentrationDataImpl(this, dest_mesh, interpolation));
502}
503
504}} // namespace electrical::diffusion
505} // namespace plask