PLaSK library

Loading...
Searching...
No Matches
ddm2d.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 <boost/version.hpp>
15
16#if (BOOST_VERSION >= 105000)
17 #include <boost/algorithm/clamp.hpp>
18 using boost::algorithm::clamp;
19#else
20 template <typename T>
21 const T& clamp(const T& v, const T& min, const T& max) {
22 if (v < min) return min;
23 if (v > max) return max;
24 return v;
25 }
26#endif
27
28#include "ddm2d.hpp"
29
30 namespace plask { namespace electrical { namespace drift_diffusion {
31
36 static inline double Neff(Tensor2<double> M, double T) {
37 constexpr double fact = phys::me * phys::kB_eV / (2.*plask::PI * phys::hb_eV * phys::hb_J);
38 double m = pow(M.c00 * M.c00 * M.c11, 0.3333333333333333);
39 return 2e-6 * pow(fact * m * T, 1.5);
40}
41
46 static inline double Ni(double Nc, double Nv, double Eg, double T) {
47 return sqrt(Nc*Nv) * exp(-Eg/(2*phys::kB_eV*T));
48}
49
50 template <typename Geometry2DType>
51 DriftDiffusionModel2DSolver<Geometry2DType>::DriftDiffusionModel2DSolver(const std::string& name)
52 : FemSolverWithMesh<Geometry2DType, RectangularMesh<2>>(name),
53 mTx(300.),
54 mEx(phys::kB_eV*mTx),
55 mNx(1e18),
56 mEpsRx(12.9),
57 mXx(sqrt((phys::epsilon0*phys::kB_J*mTx*mEpsRx)/(phys::qe*phys::qe*mNx))*1e3),
58 //mKx(100.),
59 mMix(1000.),
60 mRx(((phys::kB_J*mTx*mMix*mNx)/(phys::qe*mXx*mXx))*1e8),
61 mJx(((phys::kB_J*mNx)*mTx*mMix/mXx)*10.),
62 //mtx(mNx/mRx),
63 mAx(mRx/mNx),
64 mBx(mRx/(mNx*mNx)),
65 mCx(mRx/(mNx*mNx*mNx)),
66 //mHx(((mKx*mTx)/(mXx*mXx))*1e12),
67 mPx((mXx*phys::qe*mNx)*1e-4), // polarization (C/m^2)
68 dU(0.002),
69 maxDelPsi0(2.),
70 maxDelPsi(0.1*dU),
71 maxDelFn(1e20),
72 maxDelFp(1e20),
73 stat(STAT_MB),
74 conttype(OHMIC),
75 needPsi0(true),
76 //loopno(0),
77 //maxerr(0.05),
78 outPotential(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getPotentials),
79 outFermiLevels(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getFermiLevels),
80 outBandEdges(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getBandEdges),
81 outCurrentDensityForElectrons(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getCurrentDensitiesForElectrons),
82 outCurrentDensityForHoles(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getCurrentDensitiesForHoles),
83 outCarriersConcentration(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getCarriersConcentration),
84 outHeat(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getHeatDensities),
85 //outConductivity(this, &DriftDiffusionModel2DSolver<Geometry2DType>::getConductivity),
86 mRsrh(false),
87 mRrad(false),
88 mRaug(false),
89 mPol(false),
90 mFullIon(true),
91 mSchottkyP(0.),
92 mSchottkyN(0.),
93 maxerrPsiI(1e-6),
94 maxerrPsi0(1e-6),
95 maxerrPsi(1e-6),
96 maxerrFn(1e-4),
97 maxerrFp(1e-4),
98 loopsPsiI(10000),
99 loopsPsi0(200),
100 loopsPsi(3),
101 loopsFn(3),
102 loopsFp(3),
103 T(300.), // TODO T=300 for kp method tests
104 T0(300.)
105{
106 onInvalidate();
107 inTemperature = 300.;
108 inTemperature.changedConnectMethod(this, &DriftDiffusionModel2DSolver<Geometry2DType>::onInputChange);
109}
110
111 template <typename Geometry2DType>
112 void DriftDiffusionModel2DSolver<Geometry2DType>::loadConfiguration(XMLReader &source, Manager &manager)
113{
114 while (source.requireTagOrEnd())
115 {
116 std::string param = source.getNodeName();
117
118 if (param == "voltage")
119 this->readBoundaryConditions(manager, source, voltage_boundary);
120 else if (param == "loop") {
121 stat = source.enumAttribute<Stat>("stat")
122 .value("MB", STAT_MB)
123 .value("FD", STAT_FD)
124 .value("Maxwell-Boltzmann", STAT_MB)
125 .value("Fermi-Dirac", STAT_FD)
126 .get(stat);
127 conttype = source.enumAttribute<ContType>("conttype")
128 .value("ohmic", OHMIC)
129 .value("Schottky", SCHOTTKY)
130 .get(conttype);
131 mSchottkyP = source.getAttribute<double>("SchottkyP", mSchottkyP);
132 mSchottkyN = source.getAttribute<double>("SchottkyN", mSchottkyN);
133 mRsrh = source.getAttribute<bool>("Rsrh", mRsrh);
134 mRrad = source.getAttribute<bool>("Rrad", mRrad);
135 mPol = source.getAttribute<bool>("Pol", mPol);
136 mRaug = source.getAttribute<bool>("Raug", mRaug);
137 mFullIon = source.getAttribute<bool>("FullIon", mFullIon);
138 maxerrPsiI = source.getAttribute<double>("maxerrVi", maxerrPsiI);
139 maxerrPsi0 = source.getAttribute<double>("maxerrV0", maxerrPsi0);
140 maxerrPsi = source.getAttribute<double>("maxerrV", maxerrPsi);
141 maxerrFn = source.getAttribute<double>("maxerrFn", maxerrFn);
142 maxerrFp = source.getAttribute<double>("maxerrFp", maxerrFp);
143 loopsPsiI = source.getAttribute<size_t>("loopsVi", loopsPsiI);
144 loopsPsi0 = source.getAttribute<size_t>("loopsV0", loopsPsi0);
145 loopsPsi = source.getAttribute<size_t>("loopsV", loopsPsi);
146 loopsFn = source.getAttribute<size_t>("loopsFn", loopsFn);
147 loopsFp = source.getAttribute<size_t>("loopsFp", loopsFp);
148 source.requireTagEnd();
149 } else if (param == "config") {
150 T0 = source.getAttribute<double>("T0", T0);
151 strained = source.getAttribute<bool>("strained", strained);
152 // quick_levels = reader.getAttribute<bool>("quick-levels", quick_levels);
153 source.requireTagEnd();
154 } else if (!this->parseFemConfiguration(source, manager)) {
155 this->parseStandardConfiguration(source, manager);
156 }
157 }
158}
159
160
161 template <typename Geometry2DType>
162 DriftDiffusionModel2DSolver<Geometry2DType>::~DriftDiffusionModel2DSolver() {
163}
164
165
166 template<typename Geometry2DType>
167 size_t DriftDiffusionModel2DSolver<Geometry2DType>::getActiveRegionMeshIndex(size_t actnum) const
168{
169 if (!this->geometry) throw NoGeometryException(this->getId());
170 if (!this->mesh) throw NoMeshException(this->getId());
171
172 size_t actlo, acthi, lon = 0, hin = 0;
173
174 shared_ptr<RectangularMesh<2>> points = this->mesh->getElementMesh();
175 size_t ileft = 0, iright = points->axis[0]->size();
176 bool in_active = false;
177 for (size_t r = 0; r < points->axis[1]->size(); ++r) {
178 bool had_active = false;
179 for (size_t c = 0; c < points->axis[0]->size(); ++c) { // In the (possible) active region
180 auto point = points->at(c,r);
181 bool active = isActive(point);
182 if (c >= ileft && c <= iright) {
183 // Here we are inside potential active region
184 if (active) {
185 if (!had_active) {
186 if (!in_active) { // active region is starting set-up new region info
187 ileft = c;
188 actlo = r;
189 lon++;
190 }
191 }
192 } else if (had_active) {
193 if (!in_active) iright = c;
194 else throw Exception("{}: Right edge of the active region not aligned.", this->getId());
195 }
196 had_active |= active;
197 }
198 }
199 in_active = had_active;
200 // Test if the active region has finished
201 if (!in_active && lon != hin) {
202 acthi = r;
203 if(hin++ == actnum) return (actlo + acthi) / 2;
204 }
205 }
206 // Test if the active region has finished
207 if (lon != hin) {
208 acthi = points->axis[1]->size();
209 if(hin++ == actnum) return (actlo + acthi) / 2;
210 }
211 throw BadInput(this->getId(), "wrong active region number {}", actnum);
212}
213
214
215
216 template <typename Geometry2DType>
217 void DriftDiffusionModel2DSolver<Geometry2DType>::onInitialize()
218{
219 if (!this->geometry) throw NoGeometryException(this->getId());
220 if (!this->mesh) throw NoMeshException(this->getId());
221 detectActiveRegions();
222
223 size = this->mesh->size();
224
225 dvnPsi0.reset(size);
226 dvnFnEta.reset(size, 1.);
227 dvnFpKsi.reset(size, 1.);
228
229 dvePsi.reset(this->mesh->getElementsCount());
230 dveFnEta.reset(this->mesh->getElementsCount(), 1.);
231 dveFpKsi.reset(this->mesh->getElementsCount(), 1.);
232 dveN.reset(this->mesh->getElementsCount());
233 dveP.reset(this->mesh->getElementsCount());
234
235 currentsN.reset(this->mesh->getElementsCount());
236 currentsP.reset(this->mesh->getElementsCount());
237
238 needPsi0 = true;
239}
240
241
242 template <typename Geometry2DType>
243 void DriftDiffusionModel2DSolver<Geometry2DType>::onInvalidate() {
244 dvnPsi0.reset();
245 dvnPsi.reset();
246 dvnFnEta.reset();
247 dvnFpKsi.reset();
248 dvePsi.reset();
249 dveFnEta.reset();
250 dveFpKsi.reset();
251 dveN.reset();
252 dveP.reset();
253 currentsN.reset();
254 currentsP.reset();
255 heats.reset();
256 regions.clear();
257 substrateMaterial.reset();
258}
259
260
261 template <>
262 inline void DriftDiffusionModel2DSolver<Geometry2DCartesian>::addCurvature(double&, double&, double&, double&,
263 double&, double&, double&, double&, double&, double&,
264 double, double, const Vec<2,double>&)
265{
266}
267
268 template <> // TODO czy to bedzie OK?
269 inline void DriftDiffusionModel2DSolver<Geometry2DCylindrical>::addCurvature(double& k44, double& k33, double& k22, double& k11,
270 double& k43, double& k21, double& k42, double& k31, double& k32, double& k41,
271 double, double, const Vec<2,double>& midpoint)
272{
273 double r = midpoint.rad_r();
274 k44 = r * k44;
275 k33 = r * k33;
276 k22 = r * k22;
277 k11 = r * k11;
278 k43 = r * k43;
279 k21 = r * k21;
280 k42 = r * k42;
281 k31 = r * k31;
282 k32 = r * k32;
283 k41 = r * k41;
284}
285
286 template <typename Geometry2DType>
287 template <CalcType calctype>
288 void DriftDiffusionModel2DSolver<Geometry2DType>::setMatrix(FemMatrix& A, DataVector<double>& B,
289 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double> &bvoltage)
290{
291 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
292
293 //auto iMesh = (this->mesh)->getElementMesh();
294 //auto temperatures = inTemperature(iMesh);
295 auto iMeshN = this->mesh;
296 auto temperaturesN = inTemperature(iMeshN);
297
298//TODO 2e-6*pow((Me(T,e,point).c00*plask::phys::me*plask::phys::kB_eV*300.)/(2.*PI*plask::phys::hb_eV*plask::phys::hb_J),1.5);
299
300 A.clear();
301 B.fill(0.);
302
303 // Set stiffness matrix and load vector
304 for (auto e: this->mesh->elements()) {
305
306 size_t i = e.getIndex();
307
308 // nodes numbers for the current element
309 size_t loleftno = e.getLoLoIndex();
310 size_t lorghtno = e.getUpLoIndex();
311 size_t upleftno = e.getLoUpIndex();
312 size_t uprghtno = e.getUpUpIndex();
313
314 // element size
315 double hx = (e.getUpper0() - e.getLower0()) / mXx; // normalised element width
316 double hy = (e.getUpper1() - e.getLower1()) / mXx; // normalised element height
317
318 Vec <2,double> midpoint = e.getMidpoint();
319 auto material = this->geometry->getMaterial(midpoint);
320
321 double T;//(300.); //TODO
322 // average temperature on the element
323 T = 0.25 * (temperaturesN[loleftno] + temperaturesN[lorghtno] + temperaturesN[upleftno] + temperaturesN[uprghtno]); // in (K)
324 double normT(T/mTx); // normalised temperature
325
326 double n, p;
327 if (calctype == CALC_PSI0) {
328 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
329 n = 0.;
330 p = 0.;
331 }
332 else {
333 double normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
334 double normEc0 = material->CB(T, 0., '*') / mEx;
335 double normNv = Neff(material->Mh(T, 0.), T) / mNx;
336 double normEv0 = material->VB(T, 0., '*') / mEx;
337 double normT = T / mTx;
338 double ePsi = 0.25 * (dvnPsi0[loleftno] + dvnPsi0[lorghtno] + dvnPsi0[upleftno] + dvnPsi0[uprghtno]);
339 n = calcN(normNc, 1., ePsi, normEc0, normT);
340 p = calcP(normNv, 1., ePsi, normEv0, normT);
341 }
342 }
343 else {
344 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
345 n = 0.;
346 p = 0.;
347 }
348 else { // earlier only this
349 n = dveN[i];
350 p = dveP[i];
351 }
352 }
353
354 double kk, kx, ky, gg, ff;
355
356 if (calctype == CALC_FN) {
357 double normEc0(0.), normNc(0.), normNv(0.), normNe(0.), normNi(0.), normMobN(0.), yn(0.);
358
359 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
360 yn = 1.; // ?
361 normMobN = 1e-3; // ?
362 normNe = 1e-20; // ?
363 }
364 else {
365 normEc0 = material->CB(T, 0., '*') / mEx;
366 normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
367 normNv = Neff(material->Mh(T, 0.), T) / mNx;
368 normNe = normNc * exp(dvePsi[i]-normEc0);
369 normNi = Ni(normNc,normNv,material->Eg(T, 0., '*'),T) / mNx;
370 normMobN = 0.5*(material->mobe(T).c00+material->mobe(T).c11) / mMix; // TODO
371
372 switch (stat) {
373 case STAT_MB: yn = 1.; break;
374 //case STAT_FD: yn = fermiDiracHalf(log(dveFnEta[i])+dvePsi[i]-normEc0)/(dveFnEta[i]*exp(dvePsi[i]-normEc0)); break;
375 case STAT_FD: yn = fermiDiracHalf((log(dveFnEta[i])+dvePsi[i]-normEc0)/normT) / (pow(dveFnEta[i],1./normT)*exp((dvePsi[i]-normEc0)/normT)); break;
376 }
377 }
378
379 kk = 1. / (3.*(hx*0.5)*(hy*0.5));
380 kx = normMobN * normNe * yn * (hy*0.5) * (hy*0.5);
381 ky = normMobN * normNe * yn * (hx*0.5) * (hx*0.5);
382 ff = gg = 0.;
383
384 if (material->kind() != Material::OXIDE && material->kind() != Material::DIELECTRIC && material->kind() != Material::EMPTY ) /*if (ttE->getL()->getID() == "QW")*/ { // TODO (only in active?)
385 if (mRsrh) {
386 //this->writelog(LOG_DATA, "Recombination SRH");
387 double normte = material->taue(T) * mAx * 1e-9; // 1e-9: ns -> s
388 double normth = material->tauh(T) * mAx * 1e-9;
389 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normNe * yn * (p + normNi) * (normNi * normth + p * normte)
390 / pow((n + normNi) * normth + (p + normNi) * normte, 2.));
391 ff += ((hx*0.5) * (hy*0.5) * (n * p - normNi * normNi) / ((n + normNi) * normth + (p + normNi) * normte));
392 }
393 if (mRrad) {
394 //this->writelog(LOG_DATA, "Recombination RAD");
395 double normB = material->B(T) / mBx;
396 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normB * normNe * yn * p);
397 ff += ((hx*0.5) * (hy*0.5) * normB * (n * p - normNi * normNi));
398 }
399 if (mRaug) {
400 //this->writelog(LOG_DATA, "Recombination AUG");
401 double normCe = material->Ce(T) / mCx;
402 double normCh = material->Ch(T) / mCx;
403 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normNe * yn * ((normCe * (2. * n * p - normNi * normNi) + normCh * p * p)));
404 ff += ((hx*0.5) * (hy*0.5) * (normCe * n + normCh * p) * (n * p - normNi * normNi));
405 }
406 }
407 }
408 else if (calctype == CALC_FP) {
409 double normEv0(0.), normNc(0.), normNv(0.), normNh(0.), normNi(0.), normMobP(0.), yp(0.);
410
411 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
412 yp = 1.; // ?
413 normMobP = 1e-3; // ?
414 normNh = 1e-20; // ?
415 }
416 else {
417 normEv0 = material->VB(T, 0., '*') / mEx;
418 normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
419 normNv = Neff(material->Mh(T, 0.), T) / mNx;
420 normNh = normNv * exp(-dvePsi[i]+normEv0);
421 normNi = Ni(normNc,normNv,material->Eg(T, 0., '*'),T) / mNx;
422 normMobP = 0.5*(material->mobh(T).c00+material->mobh(T).c11) / mMix; // TODO
423
424 switch (stat) {
425 case STAT_MB: yp = 1.; break;
426 //case STAT_FD: yp = fermiDiracHalf(log(dveFpKsi[i])-dvePsi[i]+normEv0)/(dveFpKsi[i]*exp(-dvePsi[i]+normEv0)); break;
427 case STAT_FD: yp = fermiDiracHalf((log(dveFpKsi[i])-dvePsi[i]+normEv0)/normT) / (pow(dveFpKsi[i],1./normT)*exp((-dvePsi[i]+normEv0)/normT)); break;
428 }
429 }
430
431 kk = 1. / (3.*(hx*0.5)*(hy*0.5));
432 kx = normMobP * normNh * yp * (hy*0.5) * (hy*0.5);
433 ky = normMobP * normNh * yp * (hx*0.5) * (hx*0.5);
434 ff = gg = 0.;
435
436 if (material->kind() != Material::OXIDE && material->kind() != Material::DIELECTRIC && material->kind() != Material::EMPTY ) /*if (ttE->getL()->getID() == "QW")*/ { // TODO (only in active?)
437 if (mRsrh) {
438 //this->writelog(LOG_DATA, "Recombination SRH");
439 double normte = material->taue(T) * mAx * 1e-9;
440 double normth = material->tauh(T) * mAx * 1e-9;
441 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normNh * yp * (n + normNi) * (normNi * normte + n * normth)
442 / pow((n + normNi) * normth + (p + normNi) * normte, 2.));
443 ff += ((hx*0.5) * (hy*0.5) * (n * p - normNi * normNi) / ((n + normNi) * normth + (p + normNi) * normte));
444 }
445 if (mRrad) {
446 //this->writelog(LOG_DATA, "Recombination RAD");
447 double normB = material->B(T) / mBx;
448 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normB * normNh * yp * n);
449 ff += ((hx*0.5) * (hy*0.5) * normB * (n * p - normNi * normNi));
450 }
451 if (mRaug) {
452 //this->writelog(LOG_DATA, "Recombination AUG");
453 double normCe = material->Ce(T) / mCx;
454 double normCh = material->Ch(T) / mCx;
455 gg += ((1./9.) * (hx*0.5) * (hy*0.5) * normNh * yp * ((normCh * (2. * n * p - normNi * normNi) + normCe * n * n)));
456 ff += ((hx*0.5) * (hy*0.5) * (normCe * n + normCh * p) * (n * p - normNi * normNi));
457 }
458 }
459 }
460 else { // CALC_PSI
461 double normEps = material->eps(T) / mEpsRx;
462
463 kk = 1. / (3.*(hx*0.5)*(hy*0.5));
464 kx = normT * normEps * (hy*0.5) * (hy*0.5);
465 ky = normT * normEps * (hx*0.5) * (hx*0.5);
466
467 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) /*if (ttE->getL()->getID() == "QW")*/ { // TODO (only in active?)
468 gg = 0.;
469 ff = 0.;
470 }
471 else {
472 gg = (1./9.) * (p + n) * (hx*0.5) * (hy*0.5);
473
474 double normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
475 double normNv = Neff(material->Mh(T, 0.), T) / mNx;
476 //double Ni = material->Ni(T) / mNx;
477 double normNd = material->Nd() / mNx;
478 double normNa = material->Na() / mNx;
479 double normNdIon = normNd;
480 double normNaIon = normNa;
481 if (!mFullIon) {
482 //this->writelog(LOG_RESULT, "Full ionization false");
483 double gD(2.), gA(4.);
484 double normEd = material->EactD(T) / mEx;
485 double normEa = material->EactA(T) / mEx;
486 double normNdTmp = (normNc/gD)*exp(-normEd);
487 double normNaTmp = (normNv/gA)*exp(-normEa);
488 normNdIon = normNd * (normNdTmp/(normNdTmp+n));
489 normNaIon = normNa * (normNaTmp/(normNaTmp+p));
490 }
491 ff = - (hx*0.5) * (hy*0.5) * (p - n + normNdIon - normNaIon);
492 if (mPol) {
493 double eII = (3.188 - material->lattC(T,'a')) / material->lattC(T,'a'); // TODO wstawic stala podloza
494 double eL = -2. * eII * material->c13(T) / material->c33(T); // TODO uzaleznic od kata teta
495 double Ppz = material->e33(T) * eL + 2. * material->e13(T) * eII; // TODO sprawdzic czy OK
496 double Ptot = material->Psp(T) + Ppz;
497 double normPtot = Ptot / mPx;
498 ff += normPtot;
499 }
500 }
501 }
502
503 // set symmetric matrix components
504 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
505 double g44, g33, g22, g11, g43, g21, g42, g31, g32, g41;
506 double v1, v2, v3, v4;
507
508 // local K
509 k44 = k33 = k22 = k11 = (kx+ky)*kk;
510 k43 = k21 = 0.5*(-2.*kx+ky)*kk;
511 k42 = k31 = 0.5*(-kx-ky)*kk;
512 k32 = k41 = 0.5*(kx-2.*ky)*kk;
513
514 // local G
515 g44 = g33 = g22 = g11 = 4.*gg;
516 g21 = g41 = g32 = g43 = 2.*gg;
517 g31 = g42 = gg;
518
519 // set stiffness matrix
520 addCurvature(k44, k33, k22, k11, k43, k21, k42, k31, k32, k41, ky, hx, midpoint); // TODO uncomment and correct after takng cylindrical structures into account
521
522 A(loleftno, loleftno) += k11 + g11;
523 A(lorghtno, lorghtno) += k22 + g22;
524 A(uprghtno, uprghtno) += k33 + g33;
525 A(upleftno, upleftno) += k44 + g44;
526
527 A(lorghtno, loleftno) += k21 + g21;
528 A(uprghtno, loleftno) += k31 + g31;
529 A(upleftno, loleftno) += k41 + g41;
530 A(uprghtno, lorghtno) += k32 + g32;
531 A(upleftno, lorghtno) += k42 + g42;
532 A(upleftno, uprghtno) += k43 + g43;
533
534 switch (calctype) {
535 case CALC_PSI0:
536 v1 = dvnPsi0[loleftno];
537 v2 = dvnPsi0[lorghtno];
538 v3 = dvnPsi0[uprghtno];
539 v4 = dvnPsi0[upleftno];
540 break;
541 case CALC_PSI:
542 v1 = dvnPsi[loleftno];
543 v2 = dvnPsi[lorghtno];
544 v3 = dvnPsi[uprghtno];
545 v4 = dvnPsi[upleftno];
546 break;
547 case CALC_FN:
548 v1 = dvnFnEta[loleftno];
549 v2 = dvnFnEta[lorghtno];
550 v3 = dvnFnEta[uprghtno];
551 v4 = dvnFnEta[upleftno];
552 break;
553 case CALC_FP:
554 v1 = dvnFpKsi[loleftno];
555 v2 = dvnFpKsi[lorghtno];
556 v3 = dvnFpKsi[uprghtno];
557 v4 = dvnFpKsi[upleftno];
558 }
559
560 B[loleftno] -= k11*v1 + k21*v2 + k31*v3 + k41*v4 + ff;
561 B[lorghtno] -= k21*v1 + k22*v2 + k32*v3 + k42*v4 + ff;
562 B[uprghtno] -= k31*v1 + k32*v2 + k33*v3 + k43*v4 + ff;
563 B[upleftno] -= k41*v1 + k42*v2 + k43*v3 + k44*v4 + ff;
564 }
565
566 // boundary conditions of the first kind
567 A.applyBC(bvoltage, B);
568
569#ifndef NDEBUG
570 double* aend = A.data + A.size;
571 for (double* pa = A.data; pa != aend; ++pa) {
572 if (isnan(*pa) || isinf(*pa))
573 throw ComputationError(this->getId(), "error in stiffness matrix at position {0} ({1})", pa-A.data, isnan(*pa)?"nan":"inf");
574 }
575#endif
576
577}
578
579
580 template <typename Geometry2DType>
581 void DriftDiffusionModel2DSolver<Geometry2DType>::savePsi()
582{
583 for (auto el: this->mesh->elements()) {
584 size_t i = el.getIndex();
585 size_t loleftno = el.getLoLoIndex();
586 size_t lorghtno = el.getUpLoIndex();
587 size_t upleftno = el.getLoUpIndex();
588 size_t uprghtno = el.getUpUpIndex();
589
590 dvePsi[i] = 0.25 * (dvnPsi[loleftno] + dvnPsi[lorghtno] + dvnPsi[upleftno] + dvnPsi[uprghtno]);
591 }
592}
593
594
595 template <typename Geometry2DType>
596 void DriftDiffusionModel2DSolver<Geometry2DType>::saveFnEta()
597{
598 for (auto el: this->mesh->elements()) {
599 size_t i = el.getIndex();
600 size_t loleftno = el.getLoLoIndex();
601 size_t lorghtno = el.getUpLoIndex();
602 size_t upleftno = el.getLoUpIndex();
603 size_t uprghtno = el.getUpUpIndex();
604
605 dveFnEta[i] = 0.25 * (dvnFnEta[loleftno] + dvnFnEta[lorghtno] + dvnFnEta[upleftno] + dvnFnEta[uprghtno]);
606 }
607}
608
609
610 template <typename Geometry2DType>
611 void DriftDiffusionModel2DSolver<Geometry2DType>::saveFpKsi()
612{
613 for (auto el: this->mesh->elements()) {
614 size_t i = el.getIndex();
615 size_t loleftno = el.getLoLoIndex();
616 size_t lorghtno = el.getUpLoIndex();
617 size_t upleftno = el.getLoUpIndex();
618 size_t uprghtno = el.getUpUpIndex();
619
620 dveFpKsi[i] = 0.25 * (dvnFpKsi[loleftno] + dvnFpKsi[lorghtno] + dvnFpKsi[upleftno] + dvnFpKsi[uprghtno]);
621 }
622}
623
624
625 template <typename Geometry2DType>
626 void DriftDiffusionModel2DSolver<Geometry2DType>::saveN()
627{
628 //this->writelog(LOG_DETAIL, "Saving electron concentration");
629
630 //auto iMesh = (this->mesh)->getElementMesh();
631 //auto temperatures = inTemperature(iMesh);
632 auto iMeshE = (this->mesh)->getElementMesh();
633 auto temperaturesE = inTemperature(iMeshE);
634
635 for (auto e: this->mesh->elements())
636 {
637 size_t i = e.getIndex();
638 Vec < 2,double> midpoint = e.getMidpoint();
639 auto material = this->geometry->getMaterial(midpoint);
640
641 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
642 dveN[i] = 0.;
643 continue;
644 }
645 //double T(300.); // TODO
646 double normNc = Neff(material->Me(temperaturesE[i], 0., '*'), temperaturesE[i]) / mNx;
647 double normEc0 = material->CB(temperaturesE[i], 0., '*') / mEx;
648 double normT = temperaturesE[i] / mTx;
649
650 dveN[i] = calcN(normNc, dveFnEta[i], dvePsi[i], normEc0, normT);
651 }
652}
653
654
655 template <typename Geometry2DType>
656 void DriftDiffusionModel2DSolver<Geometry2DType>::saveP()
657{
658 //this->writelog(LOG_DETAIL, "Saving hole concentration");
659
660 //auto iMesh = (this->mesh)->getElementMesh();
661 //auto temperatures = inTemperature(iMesh);
662 auto iMeshE = (this->mesh)->getElementMesh();
663 auto temperaturesE = inTemperature(iMeshE);
664
665 for (auto e: this->mesh->elements())
666 {
667 size_t i = e.getIndex();
668 Vec<2,double> midpoint = e.getMidpoint();
669 auto material = this->geometry->getMaterial(midpoint);
670
671 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
672 dveP[i] = 0.;
673 continue;
674 }
675
676 //double T(300.); // TODO
677 double normNv = Neff(material->Mh(temperaturesE[i], 0.), temperaturesE[i]) / mNx;
678 double normEv0 = material->VB(temperaturesE[i], 0., '*') / mEx;
679 double normT = temperaturesE[i] / mTx;
680
681 dveP[i] = calcP(normNv, dveFpKsi[i], dvePsi[i], normEv0, normT);
682 }
683}
684
685
686 template <typename Geometry2DType>
687 template <CalcType calctype>
688 double DriftDiffusionModel2DSolver<Geometry2DType>::addCorr(DataVector<double>& corr, const BoundaryConditionsWithMesh <RectangularMesh<2>::Boundary,double>& vconst)
689{
690 //this->writelog(LOG_DEBUG, "Adding corrections");
691
692 double err;
693
694 //double tMaxRelUpd = 0.; // update/old_value = this will be the result
695
696 for (auto cond: vconst)
697 for (auto i: cond.place)
698 corr[i] = 0.;
699
700 if (calctype == CALC_PSI0) {
701 err = 0.;
702 double normDel = maxDelPsi0 / mEx;
703 for (std::size_t i = 0; i < this->mesh->size(); ++i) {
704 corr[i] = clamp(corr[i], -normDel, normDel);
705 err = std::max(err, std::abs(corr[i]));
706 dvnPsi0[i] += corr[i];
707 }
708 this->writelog(LOG_DETAIL, "Maximum update for the built-in potential: {:g} V", err*mEx);
709 }
710 else if (calctype == CALC_PSI) {
711 err = 0.;
712 double normDel = maxDelPsi / mEx;
713 for (std::size_t i = 0; i < this->mesh->size(); ++i) {
714 corr[i] = clamp(corr[i], -normDel, normDel);
715 err = std::max(err, std::abs(corr[i]));
716 dvnPsi[i] += corr[i];
717 }
718 this->writelog(LOG_DETAIL, "Maximum update for the potential: {:g} V", err*mEx);
719 }
720 else if (calctype == CALC_FN) {
721 err = 0.;
722 //double normDel = maxDelFn / mEx;
723 for (std::size_t i = 0; i < this->mesh->size(); ++i) {
724 dvnFnEta[i] += corr[i];
725 err = std::max(err, std::abs(corr[i]/dvnFnEta[i]));
726 }
727 this->writelog(LOG_DETAIL, "Maximum relative update for the quasi-Fermi energy level for electrons: {0}.", err);
728 }
729 else if (calctype == CALC_FP) {
730 err = 0.;
731 //double normDel = maxDelFp / mEx;
732 for (std::size_t i = 0; i < this->mesh->size(); ++i) {
733 dvnFpKsi[i] += corr[i];
734 err = std::max(err, std::abs(corr[i]/dvnFpKsi[i]));
735 }
736 this->writelog(LOG_DETAIL, "Maximum relative update for the quasi-Fermi energy level for holes: {0}.", err);
737 }
738 return err; // for Psi -> normalised (max. delPsi)
739
740 /*double maxRelUpd(0.);
741 double mcNorm;
742 if (calctype == CALC_PSI0) {
743 mcNorm = maxDelPsi0/mEx;
744 for (int i = 0; i < this->mesh->size(); ++i) {
745 if (dvnDeltaPsi[i] > mcNorm) dvnDeltaPsi[i] = mcNorm;
746 else if (dvnDeltaPsi[i] < -mcNorm) dvnDeltaPsi[i] = -mcNorm;
747 if (std::abs(dvnDeltaPsi[i]/dvnPsi[i]) > maxRelUpd) maxRelUpd = std::abs(dvnDeltaPsi[i]/dvnPsi[i]);
748 dvnPsi[i] = dvnPsi[i] + dvnDeltaPsi[i];
749 }
750 }
751 else if (calctype == CALC_PSI) {
752 mcNorm = maxDelPsi/mEx;
753 for (int i = 0; i < this->mesh->size(); ++i) {
754 if (dvnDeltaPsi[i] > mcNorm) dvnDeltaPsi[i] = mcNorm;
755 else if (dvnDeltaPsi[i] < -mcNorm) dvnDeltaPsi[i] = -mcNorm;
756 if (std::abs(dvnDeltaPsi[i]/dvnPsi[i]) > maxRelUpd) maxRelUpd = std::abs(dvnDeltaPsi[i]/dvnPsi[i]);
757 dvnPsi[i] = dvnPsi[i] + dvnDeltaPsi[i];
758 }
759 }
760 return maxRelUpd;*/
761}
762
763 template <typename Geometry2DType>
764 void DriftDiffusionModel2DSolver<Geometry2DType>::computePsiI() {
765
766 this->writelog(LOG_INFO, "Calculating built-in potential");
767
768 typedef std::pair<const Material*, unsigned> KeyT;
769 std::map<KeyT, double> cache;
770
771 dvnPsi0.reset(size, 0.);
772
773 //auto iMesh = (this->mesh)->getElementMesh();
774 //auto temperatures = inTemperature(iMesh);
775 auto iMeshE = (this->mesh)->getElementMesh();
776 auto temperaturesE = inTemperature(iMeshE);
777
778 for (auto el: this->mesh->elements()) {
779 size_t i = el.getIndex();
780 // point and material in the middle of the element
781 Vec < 2,double> midpoint = el.getMidpoint();
782 auto material = this->geometry->getMaterial(midpoint);
783
784 // average temperature on the element
785 // double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]); // LP_09.2015
786 //double T(300.); // Temperature in the current element
787 double T = temperaturesE[i]; // Temperature in the current element
788
789 KeyT key = std::make_pair(material.get(), unsigned(0.5+T*100.)); // temperature precision 0.01 K
790 auto found = cache.find(key);
791
792 double epsi;
793 if (found != cache.end()) {
794 epsi = found->second;
795 }
796 else {
797 if (material->kind() == Material::OXIDE || material->kind() == Material::DIELECTRIC || material->kind() == Material::EMPTY ) { // 26.01.2016
798 cache[key] = epsi = 0.;
799 continue;
800 }
801 // normalise material parameters and temperature
802 double normEc0 = material->CB(T, 0., '*') / mEx;
803 double normEv0 = material->VB(T, 0., '*', 'h') / mEx;
804 double normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
805 double normNv = Neff(material->Mh(T, 0), T) / mNx;
806 double normNd = material->Nd() / mNx;
807 double normNa = material->Na() / mNx;
808 double normEd = material->EactD(T) / mEx;
809 double normEa = material->EactA(T) / mEx;
810 double normT = T / mTx;
811 std::size_t loop = 0;
812 cache[key] = epsi = findPsiI(normEc0, normEv0, normNc, normNv, normNd, normNa, normEd, normEa, 1., 1., normT, loop);
813 }
814
815 size_t loleftno = el.getLoLoIndex();
816 size_t lorghtno = el.getUpLoIndex();
817 size_t upleftno = el.getLoUpIndex();
818 size_t uprghtno = el.getUpUpIndex();
819 dvnPsi0[loleftno] += epsi;
820 dvnPsi0[lorghtno] += epsi;
821 dvnPsi0[upleftno] += epsi;
822 dvnPsi0[uprghtno] += epsi;
823 }
824 divideByElements(dvnPsi0);
825
826 if (conttype == SCHOTTKY) {
827 // Store boundary conditions for current mesh
828 auto vconst = voltage_boundary(this->mesh, this->geometry);
829 for (auto cond: vconst) {
830 for (auto i: cond.place) {
831 if (cond.value == 0)
832 dvnPsi0[i] += mSchottkyN/mEx;
833 else if (cond.value != 0)
834 dvnPsi0[i] += mSchottkyP/mEx;
835 }
836 }
837 }
838}
839
840 template <typename Geometry2DType>
841 double DriftDiffusionModel2DSolver<Geometry2DType>::findPsiI(double iEc0, double iEv0, double iNc, double iNv, double iNd, double iNa, double iEd, double iEa, double iFnEta, double iFpKsi, double iT, std::size_t &loop) const
842{
843 double tPsi0(0.), // calculated normalized initial potential
844 tPsi0a = (-15.) / mEx, // normalized edge of the initial range
845 tPsi0b = (15.) / mEx, // normalized edge of the initial range
846 tPsi0h = (0.1) / mEx, // normalized step in the initial range calculations
847 tN = 0., tP = 0., // normalized electron/hole concentrations
848 tNtot, tNtota = (-1e30) / mNx, tNtotb = (1e30) / mNx; // normalized carrier concentration and its initial values for potentials at range edges
849
850 // Initial calculations
851
852 int tPsi0n = static_cast<int>(round((tPsi0b-tPsi0a)/tPsi0h)) + 1 ; // number of initial normalized potential values
853
854 std::vector < double> tPsi0v(tPsi0n); // normalized potential values to check
855 for (int i = 0; i < tPsi0n; ++i)
856 tPsi0v[i] = tPsi0a + i*tPsi0h;
857
858 for (int i = 0; i < tPsi0n; ++i) {
859 tN = calcN(iNc, iFnEta, tPsi0v[i], iEc0, iT);
860 tP = calcP(iNv, iFpKsi, tPsi0v[i], iEv0, iT);
861
862 double iNdIon = iNd;
863 double iNaIon = iNa;
864
865 if (!mFullIon)
866 {
867 //this->writelog(LOG_RESULT, "Full ionization false");
868 double gD(2.), gA(4.);
869 double iNdTmp = (iNc/gD)*exp(-iEd);
870 double iNaTmp = (iNv/gA)*exp(-iEa);
871 iNdIon = iNd * (iNdTmp/(iNdTmp+tN));
872 iNaIon = iNa * (iNaTmp/(iNaTmp+tP));
873 }
874
875 tNtot = tP - tN + iNdIon - iNaIon; // total normalized carrier concentration
876
877 if (tNtot < 0.) {
878 if (tNtot > tNtota) {
879 tNtota = tNtot;
880 tPsi0b = tPsi0v[i];
881 }
882 }
883 else if (tNtot > 0.) {
884 if (tNtot < tNtotb) {
885 tNtotb = tNtot;
886 tPsi0a = tPsi0v[i];
887 }
888 }
889 else // found initial normalised potential
890 return tPsi0v[i];
891 }
892
893 // Precise calculations
894
895 double tPsiUpd = 1e30, // normalised potential update
896 tTmpA, tTmpB; // temporary data
897
898 std::size_t tL = 0; // loop counter
899 while ((std::abs(tPsiUpd) > (maxerrPsiI)/mEx) && (tL < loopsPsiI)) {
900 tTmpA = (tNtotb-tNtota) / (tPsi0b-tPsi0a);
901 tTmpB = tNtota - tTmpA*tPsi0a;
902 tPsi0 = - tTmpB/tTmpA; //Psi Check Value
903 tN = calcN(iNc, iFnEta, tPsi0, iEc0, iT);
904 tP = calcP(iNv, iFpKsi, tPsi0, iEv0, iT);
905
906 double iNdIon = iNd;
907 double iNaIon = iNa;
908
909 if (!mFullIon) {
910 double gD(2.), gA(4.);
911 double iNdTmp = (iNc/gD) * exp(-iEd);
912 double iNaTmp = (iNv/gA) * exp(-iEa);
913 iNdIon = iNd * (iNdTmp / (iNdTmp+tN));
914 iNaIon = iNa * (iNaTmp / (iNaTmp+tP));
915 }
916
917 tNtot = tP - tN + iNdIon - iNaIon; // total normalized carrier concentration
918
919 if (tNtot < 0.) {
920 tNtota = tNtot;
921 tPsi0b = tPsi0;
922 }
923 else if (tNtot > 0.) {
924 tNtotb = tNtot;
925 tPsi0a = tPsi0;
926 }
927 else { // found initial normalized potential
928 loop = tL;
929 //this->writelog(LOG_DEBUG, "{0} loops done. Calculated energy level corresponding to the initial potential: {1} eV", tL, (tPsi0)*mEx); // TEST
930 return tPsi0;
931 }
932
933 tPsiUpd = tPsi0b - tPsi0a;
934 #ifndef NDEBUG
935 if (!tL)
936 this->writelog(LOG_DEBUG, "Initial potential correction: {0} eV", (tPsiUpd)*mEx); // TEST
937 else
938 this->writelog(LOG_DEBUG, " {0} eV", (tPsiUpd)*mEx); // TEST
939 #endif
940 ++tL;
941 }
942
943 loop = tL;
944 //this->writelog(LOG_INFO, "{0} loops done. Calculated energy level corresponding to the initial potential: {1} eV", tL, (tPsi0)*mEx); // TEST
945
946 return tPsi0;
947}
948
949 template <typename Geometry2DType>
950 double DriftDiffusionModel2DSolver<Geometry2DType>::findEnergyLevels() {
951 this->writelog(LOG_INFO, "Finding energy levels..");
952
953 hh2m = 0.5 * phys::hb_eV * phys::hb_J * 1e9 * 1e9 / phys::me;
954 Eshift = 20.;
955
956 double kx = 0., ky = 0.;
957
958 bool potentialWell_el = false;
959 bool potentialWell_hh = false;
960 bool potentialWell_lh = false;
961
962 potentialWell_el = checkWell("el");
963
964 if (potentialWell_el)
965 {
966 double dzdz1 = 1. / (dz*dz*1e6);
967 double dz1 = 1. / (dz*1e3);
968
969 this->writelog(LOG_DETAIL, "Creating matrix for electrons..\n");
970
971 int K = 1;
972 int N = nn * K;
973
974 this->writelog(LOG_DETAIL, "\tsize of the matrix for CB: {0} x {1}", N, N);
975
976 Eigen::MatrixXcd Hc(N, N);
978 std::complex<double> Hc_zero(0., 0.);
979 for (size_t i=1; i<=N; ++i)
980 for (size_t j=1; j <= N; ++j)
981 Hc(i-1, j-1) = Hc_zero;
982
983 for (size_t i=1; i<=nn; ++i)
984 {
985 Vec<2, double> point_LE;
986 point_LE[0] = meshActMid->axis[0]->at(0); // TODO not only for 0
987 point_LE[1] = meshActMid->axis[1]->at(i-1);
988
989 Vec<2, double> point_RI;
990 point_RI[0] = meshActMid->axis[0]->at(0); // TODO not only for 0
991 point_RI[1] = meshActMid->axis[1]->at(i);
992
993 shared_ptr<Material> m_LE = this->geometry->getMaterial(point_LE);
994 shared_ptr<Material> m_RI = this->geometry->getMaterial(point_RI);
995
996 double CBel = 0.5 * (m_LE->CB(T, 0.,'*') + m_RI->CB(T, 0., '*')) - (getPotentials(meshAct, INTERPOLATION_LINEAR)).at(i);
997 this->writelog(LOG_DETAIL, "\tCBel for kp, node: {0}, CBel: {1}", i, CBel);
998
999 //this->writelog(LOG_DETAIL, "creating Hc matrix, central node: {0}, materials: {1} and {2}", i, m_le->name(), m_ri->name());
1000
1001 std::complex<double> Hc_11_LE(0., 0.), Hc_11_CE(0., 0.), Hc_11_RI(0., 0.);
1002 {
1003 double y0_LE = 1. / m_LE->Me(T).c00; // TODO or sth else than c00?
1004 double y0_RI = 1. / m_RI->Me(T).c00; // TODO or sth else than c00?
1005 double y0_CE = 0.5 * (y0_LE + y0_RI);
1006
1008 double beta_11_LE = - y0_LE * hh2m;
1009 double beta_11_RI = - y0_RI * hh2m;
1010 double alpha_11_CE = y0_CE * hh2m * (kx * kx + ky * ky);
1011 double Ec0_11_CE = CBel + Eshift; // TODO CBel[i] must be CB from elem + Psi fo elem and then calc average for node
1012 double v_11_CE = Ec0_11_CE;
1013 double s_11_CE = 0.;
1014 /*if (strain) // TODO
1015 {
1016 double s_11_LE = 2. * vE[e_LE].ac() * (1. - vE[e_LE].c12() / vE[e_LE].c11()) * vE[e_LE].gEps(); /// [001]
1017 double s_11_RI = 2. * vE[e_RI].ac() * (1. - vE[e_RI].c12() / vE[e_RI].c11()) * vE[e_RI].gEps(); /// [001]
1018 s_11_CE = 0.5 * (s_11_LE + s_11_RI);
1019 }*/
1020
1021 if (i>1)
1022 {
1023 Hc_11_LE.real(beta_11_LE * dzdz1);
1024 Hc(i - 1, i - 2) = Hc_11_LE;
1025 //this->writelog(LOG_DETAIL, "Hc_11_LE: {0}", Hc_11_LE());
1026 }
1027 {
1028 Hc_11_CE.real(alpha_11_CE - beta_11_LE * dzdz1 - beta_11_RI * dzdz1 + v_11_CE + s_11_CE);
1029 Hc(i - 1, i - 1) = Hc_11_CE;
1030 //this->writelog(LOG_DETAIL, "Hc_11_CE: {0}", Hc_11_CE());
1031 }
1032 if (i<nn)
1033 {
1034 Hc_11_RI.real(beta_11_RI * dzdz1);
1035 Hc(i - 1, i) = Hc_11_RI;
1036 //this->writelog(LOG_DETAIL, "Hc_11_RI: {0}", Hc_11_RI());
1037 }
1038 }
1039 }
1040 this->writelog(LOG_INFO, "Done.");
1041
1042 this->writelog(LOG_INFO, "Finding energy levels and wave functions for electrons..");
1043 Eigen::ComplexEigenSolver<Eigen::MatrixXcd> ces;
1044 ces.compute(Hc);
1045 int nEigVal = ces.eigenvalues().rows();
1046 this->writelog(LOG_INFO, "number of eigenvalues (Hc): {0}", nEigVal);
1047 if (nEigVal<1)
1048 return 1;
1049 this->writelog(LOG_INFO, "Done.");
1050
1051 {
1052 this->writelog(LOG_INFO, "Analyzing the solutions..");
1053 lev_el.clear();
1054 n_lev_el = 0;
1055 for (size_t i = 1; i < nEigVal; ++i)
1056 {
1057 if ((ces.eigenvalues()[i].real() - Eshift > CBelMin) && (ces.eigenvalues()[i].real() - Eshift < CBelMax))
1058 {
1059 std::vector<double> sum(K, 0.);
1060 for (size_t j = 1; j < nn*1; j += 1)
1061 {
1062 sum[0] += (pow(ces.eigenvectors().col(i)[j].real(), 2.) + pow(ces.eigenvectors().col(i)[j].imag(), 2.));
1063 }
1064 std::string carrier;
1065 if (sum[0])
1066 {
1067 carrier = "el";
1068 lev_el.push_back(ces.eigenvalues()[i].real() - Eshift);
1069 n_lev_el++;
1070 }
1071 }
1072 }
1073 this->writelog(LOG_INFO, "Done.");
1074
1075 this->writelog(LOG_INFO, "Sorting electron energy levels..");
1076 std::sort(lev_el.begin(), lev_el.end());
1077 for (size_t i = 0; i < n_lev_el; ++i)
1078 {
1079 this->writelog(LOG_INFO, "energy level for electron {0}: {1} eV", i, lev_el[i]);
1080 }
1081 this->writelog(LOG_INFO, "Done.");
1082 }
1083
1084 Hc.resize(0, 0);
1085
1086 // TODO: hh and lh
1087 }
1088
1089 this->writelog(LOG_INFO, "Done.");
1090
1091 return 0;
1092}
1093
1094 template <typename Geometry2DType>
1095 bool DriftDiffusionModel2DSolver<Geometry2DType>::checkWell(std::string _carrier) {
1096 if (_carrier == "el")
1097 {
1098 std::vector<double> CBel;
1099 this->writelog(LOG_DETAIL, "Checking the confinement for electrons..");
1100 CBel.clear();
1101 for (size_t i = 0; i < ne+2; ++i)
1102 {
1103 //double z_avg = 0.5*(z[i] + z[i+1]);
1104 Vec<2, double> point = meshActMid->at(0,i);
1105 //point[0] = meshActMid->axis[0]->at(0); // TODO tu musi byc jakis element haxis dla danego obszaru
1106 //point[1] = meshActMid->axis[1]->at(i);
1107
1108 //this->writelog(LOG_INFO, "position of element {0}: {1} um, {2} um", i, r_at_0, z_avg);
1109
1110 shared_ptr<Material> material = this->geometry->getMaterial(point);
1111 //this->writelog(LOG_DETAIL, "material found");
1112 //this->writelog(LOG_DETAIL, "element {0}: {1}", i, material->name());
1113
1114 //CBel.push_back(material->CB(T, 0., '*') - getPotentials(meshActMid->at(0,i)) * mEx);
1115 }
1117 //for (size_t i = 0; i < nn + 2; ++i)
1118 //{
1119 // CBel.push_back(material->CB(T, 0., '*') - dvnPsi[..] * mEx);
1120 //}
1121
1122
1123
1124 for (size_t i = 0; i < nn+2; ++i)
1125 CBel.push_back(5.0);
1126 for (size_t i = 60; i < 140; ++i)
1127 CBel[i] = 4.5;
1129 CBelMin = 1e6;
1130 CBelMax = -1e6;
1131 for (size_t i = 0; i < nn+2; ++i)
1132 {
1133 if (CBel[i] < CBelMin)
1134 CBelMin = CBel[i];
1135 if (CBel[i] > CBelMax)
1136 CBelMax = CBel[i];
1137 }
1139 CBel[0] = CBelMax;
1140 CBel[nn+1] = CBelMax;
1141 //for (size_t i = 0; i < nz+2; ++i) /// TEST
1142 // this->writelog(LOG_DETAIL, "node {0}: CBel = {1} eV", i, CBel[i]);
1143
1144 this->writelog(LOG_INFO, "Done.");
1145 }
1146
1147 return true;
1148}
1149
1150 template <typename Geometry2DType>
1151 double DriftDiffusionModel2DSolver<Geometry2DType>::compute(unsigned loops)
1152{
1153 bool was_initialized = !this->initCalculation();
1154 needPsi0 |= !was_initialized;
1155
1156 //heats.reset(); // LP_09.2015
1157
1158 auto iMesh = (this->mesh)->getElementMesh();
1159 auto temperatures = inTemperature(iMesh);
1160
1161 // Store boundary conditions for current mesh
1162 auto vconst = voltage_boundary(this->mesh, this->geometry);
1163
1164 this->writelog(LOG_INFO, "Running drift-diffusion calculations for a single voltage");
1165
1166 std::unique_ptr<FemMatrix> pA(this->getMatrix());
1167 FemMatrix& A = *pA.get();
1168
1169 DataVector<double> B(size);
1170
1171 double errorPsi0 = 0.;
1172
1173 if (needPsi0) {
1174 computePsiI();
1175 errorPsi0 = 2. * maxerrPsi0;
1176 unsigned iter = 0;
1177 while (errorPsi0 > maxerrPsi0 && iter < loopsPsi0) {
1178 setMatrix<CALC_PSI0>(A, B, vconst);
1179 A.solve(B);
1180 errorPsi0 = addCorr<CALC_PSI0>(B, vconst); // max. update
1181 this->writelog(LOG_DEBUG, "Initial potential maximum update: {0}", errorPsi0*mEx); // czy dla Fn i Fp tez bedzie mEx?
1182 iter += 1;
1183 }
1184 if (!dvnPsi) {
1185 dvnPsi = dvnPsi0.copy();
1186 savePsi();
1187 }
1188 saveN();
1189 saveP();
1190 needPsi0 = false;
1191 }
1192
1193 assert(dvnPsi);
1194
1195 // Apply boundary conditions of the first kind
1196 bool novoltage = true;
1197 for (auto cond: vconst) {
1198 for (auto i: cond.place) {
1199 double dU = cond.value / mEx;
1200 novoltage = novoltage && dU == 0.;
1201 dvnPsi[i] = dvnPsi0[i] + dU;
1202 dvnFnEta[i] = exp(-dU);
1203 dvnFpKsi[i] = exp(+dU);
1204 }
1205 }
1206 if (novoltage) {
1207 if (was_initialized) {
1208 std::copy(dvnPsi0.begin(), dvnPsi0.end(), dvnPsi.begin());
1209 dvnFnEta.fill(1.);
1210 dvnFpKsi.fill(1.);
1211 }
1212 return errorPsi0;
1213 }
1214
1215 savePsi();
1216 saveFnEta();
1217 saveFpKsi();
1218 saveN();
1219 saveP();
1220
1221 if (loops == 0) loops = std::numeric_limits<unsigned>::max();
1222 unsigned loopno = 0;
1223 double errorPsi = 2.*maxerrPsi, errorFn = 2.*maxerrFn, errorFp = 2.*maxerrFp, err;
1224
1225 while ((errorPsi > maxerrPsi || errorFn > maxerrFn || errorFp > maxerrFp) && loopno < loops) {
1226 this->writelog(LOG_DETAIL, "Calculating potential");
1227 unsigned itersPsi = 0;
1228 errorPsi = 0.;
1229 err = 2.*maxerrPsi;
1230 while(err > maxerrPsi && itersPsi < loopsPsi) {
1231 setMatrix<CALC_PSI>(A, B, vconst);
1232 A.solve(B);
1233 err = addCorr<CALC_PSI>(B, vconst); // max. update
1234 if (err > errorPsi) errorPsi = err;
1235 this->writelog(LOG_DETAIL, "Maximum potential update: {0}", err*mEx); // czy dla Fn i Fp tez bedzie mEx?
1236 savePsi();
1237 saveN();
1238 saveP();
1239 itersPsi += 1;
1240 }
1241 this->writelog(LOG_DETAIL, "Calculating quasi-Fermi level for electrons");
1242 unsigned itersFn = 0;
1243 errorFn = 0.;
1244 err = 2.*maxerrFn;
1245 while(err > maxerrFn && itersFn < loopsFn) {
1246 setMatrix<CALC_FN>(A, B, vconst);
1247 A.solve(B);
1248 err = addCorr<CALC_FN>(B, vconst); // max. update
1249 if (err > errorFn) errorFn = err;
1250 this->writelog(LOG_DETAIL, "Maximum electrons quasi-Fermi level update: {0}", err*mEx); // czy dla Fn i Fp tez bedzie mEx?
1251 saveFnEta();
1252 saveN();
1253 itersFn += 1;
1254 }
1255 this->writelog(LOG_DETAIL, "Calculating quasi-Fermi energy level for holes");
1256 unsigned itersFp = 0;
1257 errorFp = 0.;
1258 err = 2.*maxerrFp;
1259 while(err > maxerrFp && itersFp < loopsFp) {
1260 setMatrix<CALC_FP>(A, B, vconst);
1261 A.solve(B);
1262 err = addCorr<CALC_FP>(B, vconst); // max. update
1263 if (err > errorFp) errorFp = err;
1264 this->writelog(LOG_DETAIL, "Maximum holes quasi-Fermi level update: {0}", err*mEx); // czy dla Fn i Fp tez bedzie mEx?
1265 saveFpKsi();
1266 saveP();
1267 itersFp += 1;
1268 }
1269 }
1270
1271 // calculate electron and hole currents (jn and jp)
1272 for (auto el: this->mesh->elements()) { // PROBLEM
1273 size_t i = el.getIndex();
1274 size_t loleftno = el.getLoLoIndex();
1275 size_t lorghtno = el.getUpLoIndex();
1276 size_t upleftno = el.getLoUpIndex();
1277 size_t uprghtno = el.getUpUpIndex();
1278 double dFnx = 0.5 * (- log(dvnFnEta[loleftno]) + log(dvnFnEta[lorghtno]) - log(dvnFnEta[upleftno]) + log(dvnFnEta[uprghtno])) // + before 0.5 due to ln(FnEta)=Fn relation
1279 / ((el.getUpper0() - el.getLower0())/mXx); // normalised [dFn/dx]
1280 double dFny = 0.5 * (- log(dvnFnEta[loleftno]) - log(dvnFnEta[lorghtno]) + log(dvnFnEta[upleftno]) + log(dvnFnEta[uprghtno])) // + before 0.5 due to ln(FnEta)=Fn relation
1281 / ((el.getUpper1() - el.getLower1())/mXx); // normalised [dFn/dy]
1282 double dFpx = - 0.5 * (- log(dvnFpKsi[loleftno]) + log(dvnFpKsi[lorghtno]) - log(dvnFpKsi[upleftno]) + log(dvnFpKsi[uprghtno])) // - before 0.5 due to -ln(FpKsi)=Fp relation
1283 / ((el.getUpper0() - el.getLower0())/mXx); // normalised [dFp/dx]
1284 double dFpy = - 0.5 * (- log(dvnFpKsi[loleftno]) - log(dvnFpKsi[lorghtno]) + log(dvnFpKsi[upleftno]) + log(dvnFpKsi[uprghtno])) // - before 0.5 due to -ln(FpKsi)=Fp relation
1285 / ((el.getUpper1() - el.getLower1())/mXx); // normalised [dFp/dy]
1286
1287 double T = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]); // in (K)
1288 //double normT(T/mTx); // normalised temperature
1289
1290 /*double n, p;
1291 if (calctype == CALC_PSI0) {
1292 double normNc = Neff(material->Me(T, 0., '*'), T) / mNx;
1293 double normEc0 = material->CB(T, 0., '*') / mEx;
1294 double normNv = Neff(material->Mh(T, 0.), T) / mNx;
1295 double normEv0 = material->VB(T, 0., '*') / mEx;
1296 double normT = T / mTx;
1297 double ePsi = 0.25 * (dvnPsi0[loleftno] + dvnPsi0[lorghtno] + dvnPsi0[upleftno] + dvnPsi0[uprghtno]);
1298 n = calcN(normNc, 1., ePsi, normEc0, normT);
1299 p = calcP(normNv, 1., ePsi, normEv0, normT);
1300 } else {
1301 n = dveN[i];
1302 p = dveP[i];
1303 }*/
1304
1305 auto midpoint = el.getMidpoint();
1306 auto material = this->geometry->getMaterial(midpoint);
1307
1308 double normMobN = 0.5*(material->mobe(T).c00+material->mobe(T).c11) / mMix;
1309 auto curN = vec(normMobN * dveN[i] * dFnx * mJx, normMobN * dveN[i] * dFny * mJx);
1310
1311 double normMobP = 0.5*(material->mobh(T).c00+material->mobh(T).c11) / mMix;
1312 auto curP = vec(normMobP * dveP[i] * dFpx * mJx, normMobP * dveP[i] * dFpy * mJx);
1313
1314 currentsN[i] = curN; // in (kA/cm^2)
1315 currentsP[i] = curP; // in (kA/cm^2)
1316 }
1317
1318 outPotential.fireChanged();
1319 outFermiLevels.fireChanged();
1320 outBandEdges.fireChanged();
1321 outCurrentDensityForElectrons.fireChanged();
1322 outCurrentDensityForHoles.fireChanged();
1323 outCarriersConcentration.fireChanged();
1324 outHeat.fireChanged();
1325
1326 return errorPsi + errorFn + errorFp;
1327}
1328
1329 template <typename Geometry2DType>
1330 void DriftDiffusionModel2DSolver<Geometry2DType>::saveHeatDensities()
1331{
1332 this->writelog(LOG_DETAIL, "Computing heat densities");
1333
1334 heats.reset(this->mesh->getElementsCount());
1335
1336 auto iMesh = (this->mesh)->getElementMesh();
1337 auto temperatures = inTemperature(iMesh);
1338
1339 /*if (heatmet == HEAT_JOULES)*/ {
1340 for (auto e: this->mesh->elements()) {
1341 size_t i = e.getIndex();
1342 size_t loleftno = e.getLoLoIndex();
1343 size_t lorghtno = e.getUpLoIndex();
1344 size_t upleftno = e.getLoUpIndex();
1345 size_t uprghtno = e.getUpUpIndex();
1346 auto midpoint = e.getMidpoint();
1347 auto material = this->geometry->getMaterial(midpoint);
1348 if (material->kind() == Material::EMPTY || this->geometry->hasRoleAt("noheat", midpoint))
1349 heats[i] = 0.;
1350 else {
1351 double T = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]); // in (K)
1352
1353 double normMobN = 0.5*(material->mobe(T).c00+material->mobe(T).c11) / mMix;
1354 double normMobP = 0.5*(material->mobh(T).c00+material->mobh(T).c11) / mMix;
1355
1356 heats[i] = ((currentsN[i].c0*currentsN[i].c0+currentsN[i].c1*currentsN[i].c1) / (normMobN*dveN[i]) + (currentsP[i].c0*currentsP[i].c0+currentsP[i].c1*currentsP[i].c1) / (normMobP*dveP[i])) * (1e12 / phys::qe);
1357
1358 /*double dvx = 0.;//0.5e6 * (- potentials[loleftno] + potentials[lorghtno] - potentials[upleftno] + potentials[uprghtno]) // LP_09.2015
1359 // / (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
1360 double dvy = 0.; //0.5e6 * (- potentials[loleftno] - potentials[lorghtno] + potentials[upleftno] + potentials[uprghtno]) // LP_09.2015
1361 // / (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
1362 heats[i] = conds[i].c00 * dvx*dvx + conds[i].c11 * dvy*dvy;*/
1363 }
1364 }
1365 }
1366}
1367
1368
1369 template <> double DriftDiffusionModel2DSolver<Geometry2DCartesian>::integrateCurrent(size_t vindex, bool onlyactive)
1370{
1371 if (!dvnPsi) throw NoValue("current densities");
1372 this->writelog(LOG_DETAIL, "Computing total current");
1373 double result = 0.;
1374 for (size_t i = 0; i < mesh->axis[0]->size()-1; ++i) {
1375 auto element = mesh->element(i, vindex);
1376 if (!onlyactive || isActive(element.getMidpoint()))
1377 result += currentsN[element.getIndex()].c1 * element.getSize0() + currentsP[element.getIndex()].c1 * element.getSize0();
1378 }
1379 if (this->getGeometry()->isSymmetric(Geometry::DIRECTION_TRAN)) result *= 2.;
1380 return result * geometry->getExtrusion()->getLength() * 0.01; // kA/cm² µm² --> mA;
1381}
1382
1383
1384 template <> double DriftDiffusionModel2DSolver<Geometry2DCylindrical>::integrateCurrent(size_t vindex, bool onlyactive)
1385{
1386 if (!dvnPsi) throw NoValue("current densities");
1387 this->writelog(LOG_DETAIL, "Computing total current");
1388 double result = 0.;
1389 for (size_t i = 0; i < mesh->axis[0]->size()-1; ++i) {
1390 auto element = mesh->element(i, vindex);
1391 if (!onlyactive || isActive(element.getMidpoint())) {
1392 double rin = element.getLower0(), rout = element.getUpper0();
1393 result += currentsN[element.getIndex()].c1 * (rout*rout - rin*rin) + currentsP[element.getIndex()].c1 * (rout*rout - rin*rin);
1394 }
1395 }
1396 return result * plask::PI * 0.01; // kA/cm² µm² --> mA
1397}
1398
1399
1400 template <typename Geometry2DType>
1401 double DriftDiffusionModel2DSolver<Geometry2DType>::getTotalCurrent(size_t nact)
1402{
1403 size_t level = getActiveRegionMeshIndex(nact);
1404 return integrateCurrent(level, true);
1405}
1406
1407
1408 template <typename Geometry2DType>
1409 const LazyData < double> DriftDiffusionModel2DSolver<Geometry2DType>::getPotentials(shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method) const
1410{
1411 if (!dvnPsi) throw NoValue("potential");
1412 this->writelog(LOG_DEBUG, "Getting potentials");
1413 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1414 return interpolate(this->mesh, dvnPsi*mEx, dst_mesh, method, this->geometry); // here the potential is rescalled (*mEx)
1415}
1416
1417
1418 template <typename Geometry2DType>
1419 const LazyData <double> DriftDiffusionModel2DSolver<Geometry2DType>::getFermiLevels(
1420 FermiLevels::EnumType what, shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method) const
1421{
1422 if (what == FermiLevels::ELECTRONS) {
1423 if (!dvnFnEta) throw NoValue("quasi-Fermi electron level");
1424 this->writelog(LOG_DEBUG, "Getting quasi-Fermi electron level");
1425
1426 DataVector<double> dvnFn(size);
1427 for (size_t i = 0; i != dvnFnEta.size(); ++i) {
1428 if (dvnFnEta[i] > 0.) dvnFn[i] = log(dvnFnEta[i]) * mEx;
1429 else dvnFn[i] = 0.;
1430 }
1431
1432 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1433 return interpolate(this->mesh, dvnFn, dst_mesh, method, this->geometry); // here the quasi-Fermi electron level is rescalled (*mEx)
1434 } else if (what == FermiLevels::HOLES) {
1435 if (!dvnFpKsi) throw NoValue("quasi-Fermi hole level");
1436 this->writelog(LOG_DEBUG, "Getting quasi-Fermi hole level");
1437
1438 DataVector<double> dvnFp(size);
1439 for (size_t i = 0; i != dvnFpKsi.size(); ++i) {
1440 if (dvnFpKsi[i] > 0.) dvnFp[i] = - log(dvnFpKsi[i]) * mEx;
1441 else dvnFp[i] = 0.;
1442 }
1443
1444 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1445 return interpolate(this->mesh, dvnFp, dst_mesh, method, this->geometry); // here the quasi-Fermi hole level is rescalled (*mEx)
1446 }
1447 assert(0);
1448 std::abort(); // to silent warning in gcc/clang release build
1449}
1450
1451
1452 template <typename Geometry2DType>
1453 const LazyData < double> DriftDiffusionModel2DSolver<Geometry2DType>::getBandEdges(
1454 BandEdges::EnumType what, shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method)
1455{
1456 if (what == BandEdges::CONDUCTION) {
1457 if (!dvnPsi) throw NoValue("conduction band edge");
1458 this->writelog(LOG_DEBUG, "Getting conduction band edge");
1459
1460 DataVector<double> dvnEc(size, 0.);
1461
1462 auto iMeshE = (this->mesh)->getElementMesh();
1463 auto temperaturesE = inTemperature(iMeshE);
1464
1465 //double T(300.); // TODO
1466 double T;
1467
1468 for (auto e: this->mesh->elements()) {
1469 size_t i = e.getIndex();
1470 size_t loleftno = e.getLoLoIndex();
1471 size_t lorghtno = e.getUpLoIndex();
1472 size_t upleftno = e.getLoUpIndex();
1473 size_t uprghtno = e.getUpUpIndex();
1474
1475 Vec <2,double> midpoint = e.getMidpoint();
1476 auto material = this->geometry->getMaterial(midpoint);
1477
1478 T = temperaturesE[i]; // Temperature in the current element
1479
1480 dvnEc[loleftno] += material->CB(T, 0., '*') - dvnPsi[loleftno] * mEx;
1481 dvnEc[lorghtno] += material->CB(T, 0., '*') - dvnPsi[lorghtno] * mEx;
1482 dvnEc[upleftno] += material->CB(T, 0., '*') - dvnPsi[upleftno] * mEx;
1483 dvnEc[uprghtno] += material->CB(T, 0., '*') - dvnPsi[uprghtno] * mEx;
1484 }
1485 divideByElements(dvnEc);
1486
1487 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1488 return interpolate(this->mesh, dvnEc, dst_mesh, method, this->geometry); // here the conduction band edge is rescalled (*mEx)
1489 } else if (what == BandEdges::VALENCE_LIGHT || what == BandEdges::VALENCE_HEAVY) {
1490 if (!dvnPsi) throw NoValue("valence band edge");
1491 this->writelog(LOG_DEBUG, "Getting valence band edge");
1492
1493 DataVector<double> dvnEv(size, 0.);
1494
1495 auto iMeshE = (this->mesh)->getElementMesh();
1496 auto temperaturesE = inTemperature(iMeshE);
1497
1498 //double T(300.); // TODO
1499 double T;
1500
1501 for (auto e: this->mesh->elements()) {
1502 size_t i = e.getIndex();
1503 size_t loleftno = e.getLoLoIndex();
1504 size_t lorghtno = e.getUpLoIndex();
1505 size_t upleftno = e.getLoUpIndex();
1506 size_t uprghtno = e.getUpUpIndex();
1507
1508 Vec<2,double> midpoint = e.getMidpoint();
1509 auto material = this->geometry->getMaterial(midpoint);
1510
1511 T = temperaturesE[i]; // Temperature in the current element
1512
1513 dvnEv[loleftno] += material->VB(T, 0., '*') - dvnPsi[loleftno] * mEx;
1514 dvnEv[lorghtno] += material->VB(T, 0., '*') - dvnPsi[lorghtno] * mEx;
1515 dvnEv[upleftno] += material->VB(T, 0., '*') - dvnPsi[upleftno] * mEx;
1516 dvnEv[uprghtno] += material->VB(T, 0., '*') - dvnPsi[uprghtno] * mEx;
1517 }
1518 divideByElements(dvnEv);
1519
1520 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1521 return interpolate(this->mesh, dvnEv, dst_mesh, method, this->geometry); // here the valence band edge is rescalled (*mEx)
1522 }
1523 assert(0);
1524 std::abort(); // to silent warning in gcc/clang release build
1525}
1526
1527
1528 template <typename Geometry2DType>
1529 const LazyData < Vec<2>> DriftDiffusionModel2DSolver<Geometry2DType>::getCurrentDensitiesForElectrons(shared_ptr<const MeshD < 2> > dest_mesh, InterpolationMethod method)
1530{
1531 if (!dvnFnEta) throw NoValue("current density");
1532 this->writelog(LOG_DEBUG, "Getting current densities");
1533 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1534 InterpolationFlags flags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN);
1535 auto result = interpolate(this->mesh->getElementMesh(), currentsN, dest_mesh, method, flags);
1536 return LazyData < Vec<2>>(result.size(),
1537 [this, dest_mesh, result, flags](size_t i) {
1538 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i)))? result[i] : Vec<2>(0.,0.);
1539 }
1540 );
1541}
1542
1543
1544 template <typename Geometry2DType>
1545 const LazyData < Vec<2>> DriftDiffusionModel2DSolver<Geometry2DType>::getCurrentDensitiesForHoles(shared_ptr<const MeshD < 2> > dest_mesh, InterpolationMethod method)
1546{
1547 if (!dvnFpKsi) throw NoValue("current density");
1548 this->writelog(LOG_DEBUG, "Getting current densities");
1549 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1550 InterpolationFlags flags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN);
1551 auto result = interpolate(this->mesh->getElementMesh(), currentsP, dest_mesh, method, flags);
1552 return LazyData < Vec<2>>(result.size(),
1553 [this, dest_mesh, result, flags](size_t i) {
1554 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i)))? result[i] : Vec<2>(0.,0.);
1555 }
1556 );
1557}
1558
1559
1560 template <typename Geometry2DType>
1561 const LazyData < double> DriftDiffusionModel2DSolver<Geometry2DType>::getCarriersConcentration(
1562 CarriersConcentration::EnumType what, shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method)
1563{
1564 DataVector<double> dvn(size, 0.);
1565
1566 //double T(300.); // TODO
1567
1568 switch (what) {
1569 case CarriersConcentration::ELECTRONS:
1570 if (!dveN) throw NoValue("electron concentration");
1571 this->writelog(LOG_DEBUG, "Getting electron concentration");
1572
1573 for (auto e: this->mesh->elements()) {
1574 size_t i = e.getIndex();
1575 size_t loleftno = e.getLoLoIndex();
1576 size_t lorghtno = e.getUpLoIndex();
1577 size_t upleftno = e.getLoUpIndex();
1578 size_t uprghtno = e.getUpUpIndex();
1579
1580 //Vec <2,double> midpoint = e.getMidpoint();
1581 //auto material = this->geometry->getMaterial(midpoint);
1582
1583 dvn[loleftno] += dveN[i] * mNx;
1584 dvn[lorghtno] += dveN[i] * mNx;
1585 dvn[upleftno] += dveN[i] * mNx;
1586 dvn[uprghtno] += dveN[i] * mNx;
1587 }
1588 divideByElements(dvn);
1589
1590 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1591 return interpolate(this->mesh, dvn, dst_mesh, method, this->geometry); // here the electron concentration is rescalled (*mNx)*/
1592
1593 case CarriersConcentration::HOLES:
1594 if (!dveP) throw NoValue("hole concentration");
1595 this->writelog(LOG_DEBUG, "Getting hole concentration");
1596
1597 for (auto e: this->mesh->elements()) {
1598 size_t i = e.getIndex();
1599 size_t loleftno = e.getLoLoIndex();
1600 size_t lorghtno = e.getUpLoIndex();
1601 size_t upleftno = e.getLoUpIndex();
1602 size_t uprghtno = e.getUpUpIndex();
1603
1604 //Vec <2,double> midpoint = e.getMidpoint();
1605 //auto material = this->geometry->getMaterial(midpoint);
1606
1607 dvn[loleftno] += dveP[i] * mNx;
1608 dvn[lorghtno] += dveP[i] * mNx;
1609 dvn[upleftno] += dveP[i] * mNx;
1610 dvn[uprghtno] += dveP[i] * mNx;
1611 }
1612 divideByElements(dvn);
1613
1614 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1615 return interpolate(this->mesh, dvn, dst_mesh, method, this->geometry); // here the hole concentration is rescalled (*mNx)*/
1616
1617 default:
1618 throw NotImplemented("{}: Carriers concentration of this type", this->getId());
1619 }
1620}
1621
1622
1623 template <typename Geometry2DType>
1624 const LazyData < double> DriftDiffusionModel2DSolver<Geometry2DType>::getHeatDensities(shared_ptr<const MeshD < 2> > dest_mesh, InterpolationMethod method)
1625{
1626 if ((!dvnFnEta)||(!dvnFpKsi)) throw NoValue("heat density");
1627 this->writelog(LOG_DEBUG, "Getting heat density");
1628 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
1629 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
1630 InterpolationFlags flags(this->geometry);
1631 auto result = interpolate(this->mesh->getElementMesh(), heats, dest_mesh, method, flags);
1632 return LazyData<double>(result.size(),
1633 [this, dest_mesh, result, flags](size_t i) {
1634 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i)))? result[i] : 0.;
1635 }
1636 );
1637}
1638
1639
1640 template <typename Geometry2DType>
1641 void DriftDiffusionModel2DSolver<Geometry2DType>::detectActiveRegions()
1642{
1643 regions.clear();
1644
1645 shared_ptr<RectangularMesh<2>> mesh = makeGeometryGrid(this->geometry->getChild());
1646 shared_ptr<RectangularMesh<2>> points = mesh->getElementMesh();
1647
1648 size_t ileft = 0, iright = points->axis[0]->size();
1649 bool in_active = false;
1650
1651 bool added_bottom_cladding = false;
1652 bool added_top_cladding = false;
1653
1654 for (size_t r = 0; r < points->axis[1]->size(); ++r) {
1655 bool had_active = false; // indicates if we had active region in this layer
1656 shared_ptr<Material> layer_material;
1657 bool layer_QW = false;
1658
1659 for (size_t c = 0; c < points->axis[0]->size(); ++c)
1660 { // In the (possible) active region
1661 auto point = points->at(c, r);
1662 auto tags = this->geometry->getRolesAt(point);
1663 bool active = false; for (const auto& tag : tags) if (tag.substr(0, 6) == "active") { active = true; break; }
1664 bool QW = tags.find("QW") != tags.end()/* || tags.find("QD") != tags.end()*/;
1665 bool substrate = tags.find("substrate") != tags.end();
1666
1667 if (substrate) {
1668 if (!substrateMaterial)
1669 substrateMaterial = this->geometry->getMaterial(point);
1670 else if (*substrateMaterial != *this->geometry->getMaterial(point))
1671 throw Exception("{0}: Non-uniform substrate layer.", this->getId());
1672 }
1673
1674 if (QW && !active)
1675 throw Exception("{0}: All marked quantum wells must belong to marked active region.", this->getId());
1676
1677 if (c < ileft) {
1678 if (active)
1679 throw Exception("{0}: Left edge of the active region not aligned.", this->getId());
1680 }
1681 else if (c >= iright) {
1682 if (active)
1683 throw Exception("{0}: Right edge of the active region not aligned.", this->getId());
1684 }
1685 else {
1686 // Here we are inside potential active region
1687 if (active) {
1688 if (!had_active) {
1689 if (!in_active)
1690 { // active region is starting set-up new region info
1691 regions.emplace_back(mesh->at(c, r));
1692 ileft = c;
1693 }
1694 layer_material = this->geometry->getMaterial(point);
1695 layer_QW = QW;
1696 }
1697 else {
1698 if (*layer_material != *this->geometry->getMaterial(point))
1699 throw Exception("{0}: Non-uniform active region layer.", this->getId());
1700 if (layer_QW != QW)
1701 throw Exception("{0}: Quantum-well role of the active region layer not consistent.", this->getId());
1702 }
1703 }
1704 else if (had_active) {
1705 if (!in_active) {
1706 iright = c;
1707
1708 // add layer below active region (cladding) LUKASZ
1709 /*auto bottom_material = this->geometry->getMaterial(points->at(ileft,r-1));
1710 for (size_t cc = ileft; cc < iright; ++cc)
1711 if (*this->geometry->getMaterial(points->at(cc,r-1)) != *bottom_material)
1712 throw Exception("{0}: Material below quantum well not uniform.", this->getId());
1713 auto& region = regions.back();
1714 double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1715 double h = mesh->axis[1]->at(r) - mesh->axis[1]->at(r-1);
1716 region.origin += Vec<2>(0., -h);
1717 this->writelog(LOG_DETAIL, "Adding bottom cladding; h = {0}",h);
1718 region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), bottom_material));*/
1719 }
1720 else
1721 throw Exception("{0}: Right edge of the active region not aligned.", this->getId());
1722 }
1723 had_active |= active;
1724 }
1725 }
1726 in_active = had_active;
1727
1728 // Now fill-in the layer info
1729 ActiveRegionInfo* region = regions.empty() ? nullptr : &regions.back();
1730 if (region) {
1731 if (!added_bottom_cladding) {
1732 if (r == 0)
1733 throw Exception("{0}: Active region cannot start from the edge of the structure.", this->getId());
1734 // add layer below active region (cladding) LUKASZ
1735 auto bottom_material = this->geometry->getMaterial(points->at(ileft, r - 1));
1736 for (size_t cc = ileft; cc < iright; ++cc)
1737 if (*this->geometry->getMaterial(points->at(cc, r - 1)) != *bottom_material)
1738 throw Exception("{0}: Material below active region not uniform.", this->getId());
1739 auto& region = regions.back();
1740 double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1741 double h = mesh->axis[1]->at(r) - mesh->axis[1]->at(r - 1);
1742 region.origin += Vec<2>(0., -h);
1743 //this->writelog(LOG_DETAIL, "Adding bottom cladding; h = {0}",h);
1744 region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), bottom_material));
1745 region.bottom = h;
1746 added_bottom_cladding = true;
1747 }
1748
1749 double h = mesh->axis[1]->at(r + 1) - mesh->axis[1]->at(r);
1750 double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1751 if (in_active) {
1752 size_t n = region->layers->getChildrenCount();
1753 shared_ptr<Block<2>> last;
1754 if (n > 0) last = static_pointer_cast<Block<2>>(static_pointer_cast<Translation<2>>(region->layers->getChildNo(n - 1))->getChild());
1755 assert(!last || last->size.c0 == w);
1756 if (last && layer_material == last->getRepresentativeMaterial() && layer_QW == region->isQW(region->size() - 1)) {
1757 //TODO check if usage of getRepresentativeMaterial is fine here (was material)
1758 last->setSize(w, last->size.c1 + h);
1759 }
1760 else {
1761 auto layer = plask::make_shared<Block<2>>(Vec<2>(w, h), layer_material);
1762 if (layer_QW) layer->addRole("QW");
1763 region->layers->push_back(layer);
1764 }
1765 }
1766 else {
1767 if (!added_top_cladding) {
1768
1769 // add layer above active region (top cladding)
1770 auto top_material = this->geometry->getMaterial(points->at(ileft, r));
1771 for (size_t cc = ileft; cc < iright; ++cc)
1772 if (*this->geometry->getMaterial(points->at(cc, r)) != *top_material)
1773 throw Exception("{0}: Material above quantum well not uniform.", this->getId());
1774 region->layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), top_material));
1775 //this->writelog(LOG_DETAIL, "Adding top cladding; h = {0}",h);
1776
1777 ileft = 0;
1778 iright = points->axis[0]->size();
1779 region->top = h;
1780 added_top_cladding = true;
1781 }
1782 }
1783 }
1784 }
1785 if (!regions.empty() && regions.back().isQW(regions.back().size() - 1))
1786 throw Exception("{0}: Quantum-well cannot be located at the edge of the structure.", this->getId());
1787
1788 if (strained && !substrateMaterial)
1789 throw BadInput(this->getId(), "strained quantum wells requested but no layer with substrate role set");
1790
1791 this->writelog(LOG_DETAIL, "Found {0} active region{1}", regions.size(), (regions.size() == 1) ? "" : "s");
1792 for (auto& region : regions) region.summarize(this);
1793}
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806//template <typename GeometryType>
1807//void DriftDiffusionModel2DSolver<GeometryType>::detectActiveRegions()
1808//{
1809// this->writelog(LOG_INFO, "Detecting active regions..");
1810//
1811// regions.clear();
1812//
1813// shared_ptr< MeshAxis > axis_vert = this->mesh->vert(); /// for the whole structure
1814// shared_ptr< MeshAxis > axis_tran = this->mesh->tran(); /// for the whole structure
1815// double r_at_0 = 0.5 * (axis_tran->at(0) + axis_tran->at(1)); // TODO
1816//
1817// shared_ptr<RectangularMesh<2>> meshMid = this->mesh->getElementMesh();
1818// bool found_substrate = false;
1819// bool found_active = false;
1820//
1821// double z1(-1.), z2(-1.); /// min. and max. z for kp method [um]
1822// double zSub(-1.); /// z for substrate [um]
1823//
1824// //for (std::size_t i = 0; i < axis_vert->size(); ++i) // TEST
1825// //{
1826// // this->writelog(LOG_INFO, "axis_vert[{0}]: {1}", i, axis_vert->at(i));
1827// //}
1828// for (std::size_t i = 0; i < axis_vert->size()-1; ++i)
1829// {
1830// double zA = axis_vert->at(i); /// z bottom
1831// double zB = axis_vert->at(i+1); /// z top
1832// double z_avg = 0.5*(zA + zB);
1833// Vec<2, double> point;
1834// point[0] = r_at_0;
1835// point[1] = z_avg;
1836//
1837// std::string role_ = "-";
1838// auto tags = this->geometry->getRolesAt(point);
1839//
1840// if (tags.find("substrate") != tags.end())
1841// {
1842// //this->writelog(LOG_INFO, "axis_vert_mid[{0}] {1}, substrate", i, z_avg);
1843// if (!found_substrate) /// first time in substrate
1844// {
1845// zSub = z_avg;
1846// }
1847// found_substrate = true;
1848// }
1849// if (tags.find("active") != tags.end())
1850// {
1851// //in_active = true;
1852// //this->writelog(LOG_INFO, "axis_vert_mid[{0}] {1}, active", i, z_avg);
1853// if (!found_active) /// first time in active
1854// {
1855// z1 = zA;
1856// z2 = zB;
1857// }
1858// else
1859// {
1860// if (z1 > zA)
1861// z1 = zA;
1862// if (z2 < zB)
1863// z2 = zB;
1864// }
1865// found_active = true;
1866// }
1867// }
1868//
1869// this->writelog(LOG_INFO, "active region is from z = {0} um to z = {1} um", z1, z2); /// in [um]
1870// this->writelog(LOG_INFO, "active region thickness: {0} nm", (z2-z1)*1e3); /// [um] -> (nm)
1871//
1872// z.clear();
1873// //z.push_back(0.); /// so here z[0] = 0, but later the z.size will be stored here
1874// dz = 0.1e-3; /// in [um], default value: 0.1 nm
1875// nn = static_cast<int> ((z2 - z1 + 1e-6) / dz) + 1;
1876// this->writelog(LOG_INFO, "no of nodes for kp method: {0}", nn);
1877// z.push_back(z1 - dz); /// bottom cladding, left edge of the element is set here: z[0]
1878// for (std::size_t i = 0; i < nn; ++i)
1879// z.push_back(z1 + i*dz);
1880// z.push_back(z2 + dz); /// top cladding
1881// ne = nn - 1;
1882//
1883// auto vaxis = plask::make_shared<OrderedAxis>();
1884// auto haxis = plask::make_shared<OrderedAxis>();
1885// this->writelog(LOG_INFO, "vaxis size: {0}", vaxis->size());
1886// this->writelog(LOG_INFO, "haxis size: {0}", haxis->size());
1887// OrderedAxis::WarningOff vaxiswoff(vaxis);
1888// OrderedAxis::WarningOff haxiswoff(haxis);
1889// for (std::size_t i = 0; i != nn+2; ++i)
1890// {
1891// vaxis->addPoint(z[i]);
1892// }
1893// haxis->addPoint(1e-4); // TODO - tu musza byc wartosci z haxis dla danego obszaru czynnego
1894// haxis->addPoint(2e-4);
1895// haxis->addPoint(3e-4);
1896//
1897// this->writelog(LOG_INFO, "vaxis size: {0}", vaxis->size());
1898// this->writelog(LOG_INFO, "haxis size: {0}", haxis->size());
1899//
1900// meshAct = plask::make_shared<RectangularMesh<2>>(haxis, vaxis, RectangularMesh<2>::ORDER_01);
1901// meshActMid = meshAct->getElementMesh(); // LUKI TODO
1902//
1903// this->writelog(LOG_INFO, "MeshAct 0 1: {0} {1} {2}", meshAct->at(0,0), meshAct->at(0,1), meshAct->at(0,2));
1904// this->writelog(LOG_INFO, "MeshAct 0 1: {0} {1} {2}", meshAct->axis[1]->at(0), meshAct->axis[1]->at(1), meshAct->axis[1]->at(2));
1905//
1906// this->writelog(LOG_INFO, "Done.");
1907//
1908// /*shared_ptr<RectangularMesh<2>> points = mesh->getElementMesh();
1909//
1910// for (size_t r = 0; r < points->axis[1]->size(); ++r) {
1911// bool had_active = false; // indicates if we had active region in this layer
1912// shared_ptr<Material> layer_material;
1913// bool layer_QW = false;
1914//
1915// for (size_t c = 0; c < points->axis[0]->size(); ++c)
1916// { // In the (possible) active region
1917//
1918//
1919// else {
1920// // Here we are inside potential active region
1921// if (active) {
1922// if (!had_active) {
1923// if (!in_active)
1924// { // active region is starting set-up new region info
1925// regions.emplace_back(mesh->at(c,r));
1926// ileft = c;
1927// }
1928// layer_material = this->geometry->getMaterial(point);
1929// layer_QW = QW;
1930// } else {
1931// if (*layer_material != *this->geometry->getMaterial(point))
1932// throw Exception("{0}: Non-uniform active region layer.", this->getId());
1933// if (layer_QW != QW)
1934// throw Exception("{0}: Quantum-well role of the active region layer not consistent.", this->getId());
1935// }
1936// } else if (had_active) {
1937// if (!in_active) {
1938// iright = c;
1939//
1940// // add layer below active region (cladding) LUKASZ
1941// //auto bottom_material = this->geometry->getMaterial(points->at(ileft,r-1));
1942// //for (size_t cc = ileft; cc < iright; ++cc)
1943// // if (*this->geometry->getMaterial(points->at(cc,r-1)) != *bottom_material)
1944// // throw Exception("{0}: Material below quantum well not uniform.", this->getId());
1945// //auto& region = regions.back();
1946// //double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1947// //double h = mesh->axis[1]->at(r) - mesh->axis[1]->at(r-1);
1948// //region.origin += Vec<2>(0., -h);
1949// //this->writelog(LOG_DETAIL, "Adding bottom cladding; h = {0}",h);
1950// //region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), bottom_material));
1951// } else
1952// throw Exception("{0}: Right edge of the active region not aligned.", this->getId());
1953// }
1954// had_active |= active;
1955// }
1956// }
1957// in_active = had_active;
1958//
1959// // Now fill-in the layer info
1960// ActiveRegionInfo* region = regions.empty()? nullptr : &regions.back();
1961// if (region) {
1962// if (!added_bottom_cladding) {
1963// if (r == 0)
1964// throw Exception("{0}: Active region cannot start from the edge of the structure.", this->getId());
1965// // add layer below active region (cladding) LUKASZ
1966// auto bottom_material = this->geometry->getMaterial(points->at(ileft,r-1));
1967// for (size_t cc = ileft; cc < iright; ++cc)
1968// if (*this->geometry->getMaterial(points->at(cc,r-1)) != *bottom_material)
1969// throw Exception("{0}: Material below active region not uniform.", this->getId());
1970// auto& region = regions.back();
1971// double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1972// double h = mesh->axis[1]->at(r) - mesh->axis[1]->at(r-1);
1973// region.origin += Vec<2>(0., -h);
1974// //this->writelog(LOG_DETAIL, "Adding bottom cladding; h = {0}",h);
1975// region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), bottom_material));
1976// region.bottom = h;
1977// added_bottom_cladding = true;
1978// }
1979//
1980// double h = mesh->axis[1]->at(r+1) - mesh->axis[1]->at(r);
1981// double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
1982// if (in_active) {
1983// size_t n = region->layers->getChildrenCount();
1984// shared_ptr<Block<2>> last;
1985// if (n > 0) last = static_pointer_cast<Block<2>>(static_pointer_cast<Translation<2>>(region->layers->getChildNo(n-1))->getChild());
1986// assert(!last || last->size.c0 == w);
1987// if (last && layer_material == last->getRepresentativeMaterial() && layer_QW == region->isQW(region->size()-1)) {
1988// //TODO check if usage of getRepresentativeMaterial is fine here (was material)
1989// last->setSize(w, last->size.c1 + h);
1990// } else {
1991// auto layer = plask::make_shared<Block<2>>(Vec<2>(w,h), layer_material);
1992// if (layer_QW) layer->addRole("QW");
1993// region->layers->push_back(layer);
1994// }
1995// } else {
1996// if (!added_top_cladding) {
1997//
1998// // add layer above active region (top cladding)
1999// auto top_material = this->geometry->getMaterial(points->at(ileft,r));
2000// for (size_t cc = ileft; cc < iright; ++cc)
2001// if (*this->geometry->getMaterial(points->at(cc,r)) != *top_material)
2002// throw Exception("{0}: Material above quantum well not uniform.", this->getId());
2003// region->layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w,h), top_material));
2004// //this->writelog(LOG_DETAIL, "Adding top cladding; h = {0}",h);
2005//
2006// ileft = 0;
2007// iright = points->axis[0]->size();
2008// region->top = h;
2009// added_top_cladding = true;
2010// }
2011// }
2012// }
2013// }
2014// if (!regions.empty() && regions.back().isQW(regions.back().size()-1))
2015// throw Exception("{0}: Quantum-well cannot be located at the edge of the structure.", this->getId());
2016//
2017// if (strained && !substrateMaterial)
2018// throw BadInput(this->getId(), "strained quantum wells requested but no layer with substrate role set");
2019// */
2020// this->writelog(LOG_DETAIL, "Found {0} active region{1}", regions.size(), (regions.size()==1)?"":"s");
2021// for (auto& region: regions) region.summarize(this);
2022//}
2023
2024
2025 template <typename Geometry2DType>
2026 void DriftDiffusionModel2DSolver<Geometry2DType>::ActiveRegionInfo::summarize(const DriftDiffusionModel2DSolver<Geometry2DType>* solver)
2027{
2028 holes = BOTH_HOLES;
2029 auto bbox = layers->getBoundingBox();
2030 total = bbox.upper[1] - bbox.lower[1] - bottom - top;
2031 solver->writelog(LOG_DETAIL, "coordinates | bbox.upper: {0} um, bbox.lower: {1} um, bottom: {2} um, top: {3} um, total: {4} um", bbox.upper[1], bbox.lower[1], bottom, top, total);
2032 materials.clear(); materials.reserve(layers->children.size());
2033 thicknesses.clear(); thicknesses.reserve(layers->children.size());
2034 for (const auto& layer : layers->children) {
2035 auto block = static_cast<Block<2>*>(static_cast<Translation<2>*>(layer.get())->getChild().get());
2036 auto material = block->singleMaterial();
2037 if (!material) throw plask::Exception("{}: Active region can consist only of solid layers", solver->getId());
2038 auto bbox = static_cast<GeometryObjectD<2>*>(layer.get())->getBoundingBox();
2039 double thck = bbox.upper[1] - bbox.lower[1];
2040 solver->writelog(LOG_DETAIL, "layer | material: {0}, thickness: {1} um", material->name(), thck);
2041 materials.push_back(material);
2042 thicknesses.push_back(thck);
2043 }
2044 /*double substra = solver->strained ? solver->substrateMaterial->lattC(solver->T0, 'a') : 0.; // TODO moze cos z tego wziac???
2045 if (materials.size() > 2) {
2046 Material* material = materials[0].get();
2047 double e;
2048 if (solver->strained) { double latt = material->lattC(solver->T0, 'a'); e = (substra - latt) / latt; }
2049 else e = 0.;
2050 double el0 = material->CB(solver->T0, e, 'G'),
2051 hh0 = material->VB(solver->T0, e, 'G', 'H'),
2052 lh0 = material->VB(solver->T0, e, 'G', 'L');
2053 material = materials[1].get();
2054 if (solver->strained) { double latt = material->lattC(solver->T0, 'a'); e = (substra - latt) / latt; }
2055 else e = 0.;
2056 double el1 = material->CB(solver->T0, e, 'G'),
2057 hh1 = material->VB(solver->T0, e, 'G', 'H'),
2058 lh1 = material->VB(solver->T0, e, 'G', 'L');
2059 for (size_t i = 2; i < materials.size(); ++i) {
2060 material = materials[i].get();
2061 if (solver->strained) { double latt = material->lattC(solver->T0, 'a'); e = (substra - latt) / latt; }
2062 else e = 0.;
2063 double el2 = material->CB(solver->T0, e, 'G');
2064 double hh2 = material->VB(solver->T0, e, 'G', 'H');
2065 double lh2 = material->VB(solver->T0, e, 'G', 'L');
2066 if ((el0 < el1 && el1 > el2) || (hh0 > hh1 && hh1 < hh2) || (lh0 > lh1 && lh1 < lh2)) {
2067 if (i != 2 && i != materials.size() - 1) {
2068 bool eb = (el0 < el1 && el1 > el2);
2069 if (eb != (hh0 > hh1 && hh1 < hh2)) holes = ConsideredHoles(holes & ~HEAVY_HOLES);
2070 if (eb != (lh0 > lh1 && lh1 < lh2)) holes = ConsideredHoles(holes & ~LIGHT_HOLES);
2071 }
2072 if (holes == NO_HOLES)
2073 throw Exception("{0}: Quantum wells in conduction band do not coincide with wells is valence band", solver->getId());
2074 if ((el0 < el1 && el1 > el2) || (hh0 > hh1 && hh1 < hh2 && holes & HEAVY_HOLES) || (lh0 > lh1 && lh1 < lh2 && holes & LIGHT_HOLES))
2075 wells.push_back(i - 1);
2076 }
2077 else if (i == 2) wells.push_back(0);
2078 if (el2 != el1) { el0 = el1; el1 = el2; }
2079 if (hh2 != hh1) { hh0 = hh1; hh1 = hh2; }
2080 if (lh2 != lh1) { lh0 = lh1; lh1 = lh2; }
2081 }
2082 }
2083 if (wells.back() < materials.size() - 2) wells.push_back(materials.size() - 1);
2084 totalqw = 0.;
2085 for (size_t i = 0; i < thicknesses.size(); ++i)
2086 if (isQW(i)) totalqw += thicknesses[i];*/
2087}
2088
2089
2090/*
2091 template <typename Geometry2DType>
2092 const LazyData < Tensor2<double>> DriftDiffusionModel2DSolver<Geometry2DType>::getConductivity(shared_ptr<const MeshD < 2> > dest_mesh, InterpolationMethod) {
2093 this->initCalculation();
2094 this->writelog(LOG_DEBUG, "Getting conductivities");
2095 loadConductivities();
2096 InterpolationFlags flags(this->geometry);
2097 return LazyData < Tensor2<double>>(new LazyDataDelegateImpl < Tensor2<double>>(dest_mesh->size(),
2098 [this, dest_mesh, flags](size_t i) -> Tensor2 < double> {
2099 auto point = flags.wrap(dest_mesh->at(i));
2100 size_t x = this->mesh->axis[0]->findUpIndex(point[0]),
2101 y = this->mesh->axis[1]->findUpIndex(point[1]);
2102 if (x == 0 || y == 0 || x == this->mesh->axis[0]->size() || y == this->mesh->axis[1]->size())
2103 return Tensor2<double>(NAN);
2104 else
2105 return this->conds[this->mesh->elements(x-1, y-1).getIndex()];
2106 }
2107 ));
2108}
2109
2110
2111 template <>
2112 double DriftDiffusionModel2DSolver<Geometry2DCartesian>::getTotalEnergy() {
2113 double W = 0.;
2114 auto T = inTemperature(this->mesh->getElementMesh());
2115 for (auto e: this->mesh->elements) {
2116 size_t ll = e.getLoLoIndex();
2117 size_t lu = e.getUpLoIndex();
2118 size_t ul = e.getLoUpIndex();
2119 size_t uu = e.getUpUpIndex();
2120 double dvx = 0.5e6 * (- potentials[ll] + potentials[lu] - potentials[ul] + potentials[uu])
2121 / (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
2122 double dvy = 0.5e6 * (- potentials[ll] - potentials[lu] + potentials[ul] + potentials[uu])
2123 / (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
2124 double w = this->geometry->getMaterial(e.getMidpoint())->eps(T[e.getIndex()]) * (dvx*dvx + dvy*dvy);
2125 double width = e.getUpper0() - e.getLower0();
2126 double height = e.getUpper1() - e.getLower1();
2127 W += width * height * w;
2128 }
2129 //TODO add outsides of comptational areas
2130 return geometry->getExtrusion()->getLength() * 0.5e-18 * phys::epsilon0 * W; // 1e-18 µm³ -> m³
2131}
2132
2133 template <>
2134 double DriftDiffusionModel2DSolver<Geometry2DCylindrical>::getTotalEnergy() {
2135 double W = 0.;
2136 auto T = inTemperature(this->mesh->getElementMesh());
2137 for (auto e: this->mesh->elements) {
2138 size_t ll = e.getLoLoIndex();
2139 size_t lu = e.getUpLoIndex();
2140 size_t ul = e.getLoUpIndex();
2141 size_t uu = e.getUpUpIndex();
2142 auto midpoint = e.getMidpoint();
2143 double dvx = 0.5e6 * (- potentials[ll] + potentials[lu] - potentials[ul] + potentials[uu])
2144 / (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
2145 double dvy = 0.5e6 * (- potentials[ll] - potentials[lu] + potentials[ul] + potentials[uu])
2146 / (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
2147 double w = this->geometry->getMaterial(midpoint)->eps(T[e.getIndex()]) * (dvx*dvx + dvy*dvy);
2148 double width = e.getUpper0() - e.getLower0();
2149 double height = e.getUpper1() - e.getLower1();
2150 W += width * height * midpoint.rad_r() * w;
2151 }
2152 //TODO add outsides of computational area
2153 return 2.*PI * 0.5e-18 * phys::epsilon0 * W; // 1e-18 µm³ -> m³
2154}
2155
2156
2157 template <typename Geometry2DType>
2158 double DriftDiffusionModel2DSolver<Geometry2DType>::getCapacitance() {
2159
2160 if (this->voltage_boundary.size() != 2) {
2161 throw BadInput(this->getId(), "cannot estimate applied voltage (exactly 2 voltage boundary conditions required)");
2162 }
2163
2164 double U = voltage_boundary[0].value - voltage_boundary[1].value;
2165
2166 return 2e12 * getTotalEnergy() / (U*U); // 1e12 F -> pF
2167}
2168
2169
2170 template <>
2171 double DriftDiffusionModel2DSolver<Geometry2DCartesian>::getTotalHeat() {
2172 double W = 0.;
2173 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
2174 for (auto e: this->mesh->elements) {
2175 double width = e.getUpper0() - e.getLower0();
2176 double height = e.getUpper1() - e.getLower1();
2177 W += width * height * heats[e.getIndex()];
2178 }
2179 return geometry->getExtrusion()->getLength() * 1e-15 * W; // 1e-15 µm³ -> m³, W -> mW
2180}
2181
2182 template <>
2183 double DriftDiffusionModel2DSolver<Geometry2DCylindrical>::getTotalHeat() {
2184 double W = 0.;
2185 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
2186 for (auto e: this->mesh->elements) {
2187 double width = e.getUpper0() - e.getLower0();
2188 double height = e.getUpper1() - e.getLower1();
2189 double r = e.getMidpoint().rad_r();
2190 W += width * height * r * heats[e.getIndex()];
2191 }
2192 return 2e-15*PI * W; // 1e-15 µm³ -> m³, W -> mW
2193}*/
2194
2195
2196 template < > std::string DriftDiffusionModel2DSolver<Geometry2DCartesian>::getClassName() const { return "ddm2d.DriftDiffusion2D"; }
2197 template < > std::string DriftDiffusionModel2DSolver<Geometry2DCylindrical>::getClassName() const { return "ddm2d.DriftDiffusionCyl"; }
2198
2199 PLASK_NO_CONVERSION_WARNING_BEGIN
2200 template struct PLASK_SOLVER_API DriftDiffusionModel2DSolver<Geometry2DCartesian>;
2201 template struct PLASK_SOLVER_API DriftDiffusionModel2DSolver<Geometry2DCylindrical>;
2202 PLASK_NO_WARNING_END
2203
2204}}} // namespace plask::electrical::thermal