PLaSK library
Loading...
Searching...
No Matches
ddm2d.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK__MODULE_ELECTRICAL_DDM2D_H
15#define PLASK__MODULE_ELECTRICAL_DDM2D_H
16
17#include <fstream>
18#include <limits>
19
20#include <plask/common/fem.hpp>
21#include <plask/plask.hpp>
22
23PLASK_NO_CONVERSION_WARNING_BEGIN
24#include <Eigen/Eigen>
25PLASK_NO_WARNING_END
26#include "fd.hpp"
27
28namespace plask { namespace electrical { namespace drift_diffusion {
29
35
43
49
53template <typename Geometry2DType>
54struct PLASK_SOLVER_API DriftDiffusionModel2DSolver : public FemSolverWithMesh<Geometry2DType, RectangularMesh<2>> {
55 protected:
56 std::size_t size;
57
58 // scalling parameters
59 double mTx;
60 double mEx;
61 double mNx;
62 double mEpsRx;
63 double mXx;
64 // double mKx; ///< thermal conductivity (W/(m*K))
65 double mMix;
66 double mRx;
67 double mJx;
68 // double mtx; ///< SRH recombination lifetime (s)
69 double mAx;
70 double mBx;
71 double mCx;
72 // double mHx; ///< heat source (W/(m^3))
73 double mPx;
74
75 double dU;
76 double maxDelPsi0;
77 double maxDelPsi;
78 double maxDelFn;
79 double maxDelFp;
80
81 Stat stat;
82 ContType conttype;
83
84 // int loopno; ///< Number of completed loops
85 // double toterr; ///< Maximum estimated error during all iterations (useful for single calculations managed by
86 // external python script) Vec<2,double> maxcur; ///< Maximum current in the structure
87
88 DataVector<double> dveN;
89 DataVector<double> dveP;
90 DataVector<double> dvePsi;
91 DataVector<double> dveFnEta;
92 DataVector<double> dveFpKsi;
93
94 DataVector<double> dvnPsi0;
95 DataVector<double> dvnPsi;
96 DataVector<double> dvnFnEta;
97 DataVector<double> dvnFpKsi;
98 DataVector<Vec<2, double>> currentsN;
99 DataVector<Vec<2, double>> currentsP;
100 DataVector<double> heats;
101
102 bool needPsi0;
103
104 double T0;
105
106 bool strained;
107
109 void onInitialize() override;
110
112 void onInvalidate() override;
113
115 size_t getActiveRegionMeshIndex(size_t actnum) const;
116
117 void onMeshChange(const typename RectangularMesh<2>::Event& evt) override {
118 SolverWithMesh<Geometry2DType, RectangularMesh<2>>::onMeshChange(evt);
119 }
120
121 void onGeometryChange(const Geometry::Event& evt) override {
122 SolverWithMesh<Geometry2DType, RectangularMesh<2>>::onGeometryChange(evt);
123 }
124
128 void computePsiI();
129
134 size_t isActive(const Vec<2>& point) const {
135 size_t no(0);
136 auto roles = this->geometry->getRolesAt(point);
137 for (auto role : roles) {
138 size_t l = 0;
139 if (role.substr(0, 6) == "active")
140 l = 6;
141 else if (role.substr(0, 8) == "junction")
142 l = 8;
143 else
144 continue;
145 if (no != 0) throw BadInput(this->getId(), "multiple 'active'/'junction' roles specified");
146 if (role.size() == l)
147 no = 1;
148 else {
149 try {
150 no = boost::lexical_cast<size_t>(role.substr(l)) + 1;
151 } catch (boost::bad_lexical_cast&) {
152 throw BadInput(this->getId(), "bad junction number in role '{0}'", role);
153 }
154 }
155 }
156 return no;
157 }
158
160 size_t isActive(const RectangularMesh<2>::Element& element) const { return isActive(element.getMidpoint()); }
161
162 private:
164 void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) { needPsi0 = true; }
165
167 double findPsiI(double iEc0,
168 double iEv0,
169 double iNc,
170 double iNv,
171 double iNd,
172 double iNa,
173 double iEd,
174 double iEa,
175 double iFnEta,
176 double iFpKsi,
177 double iT,
178 std::size_t& loop) const;
179
189 double calcN(double iNc, double iFnEta, double iPsi, double iEc0, double iT) const {
190 switch (stat) {
191 // case STAT_MB: return ( iNc * iFnEta * exp(iPsi-iEc0) );
192 case STAT_MB:
193 // this->writelog(LOG_INFO, "Maxwell-Boltzmann statistics");
194 return (iNc * pow(iFnEta, 1. / iT) * exp((iPsi - iEc0) / iT));
195 // case STAT_FD: return ( iNc * fermiDiracHalf(log(iFnEta) + iPsi - iEc0) );
196 case STAT_FD:
197 // this->writelog(LOG_INFO, "Fermi-Dirac statistics");
198 return (iNc * fermiDiracHalf((log(iFnEta) + iPsi - iEc0) / iT));
199 }
200 return NAN;
201 }
202
212 double calcP(double iNv, double iFpKsi, double iPsi, double iEv0, double iT) const {
213 switch (stat) {
214 // case STAT_MB: return ( iNv * iFpKsi * exp(iEv0-iPsi) );
215 case STAT_MB: return (iNv * pow(iFpKsi, 1. / iT) * exp((iEv0 - iPsi) / iT));
216 // case STAT_FD: return ( iNv * fermiDiracHalf(log(iFpKsi) - iPsi + iEv0) );
217 case STAT_FD: return (iNv * fermiDiracHalf((log(iFpKsi) - iPsi + iEv0) / iT));
218 }
219 return NAN;
220 }
221
222 void divideByElements(DataVector<double>& values) {
223 size_t majs = this->mesh->majorAxis()->size(), mins = this->mesh->minorAxis()->size();
224 if (mins == 0 || majs == 0) return;
225 for (size_t j = 1, jend = mins - 1; j < jend; ++j) values[j] *= 0.5;
226 for (size_t i = 1, iend = majs - 1; i < iend; ++i) {
227 values[mins * i] *= 0.5;
228 for (size_t j = 1, jend = mins - 1; j < jend; ++j) values[mins * i + j] *= 0.25;
229 values[mins * (i + 1) - 1] *= 0.5;
230 }
231 for (size_t j = mins * (majs - 1) + 1, jend = this->mesh->size() - 1; j < jend; ++j) values[j] *= 0.5;
232 }
233
234 void savePsi0();
235 void savePsi();
236 void saveFnEta();
237 void saveFpKsi();
238 void saveN();
239 void saveP();
240
242 template <CalcType calctype>
243 double addCorr(DataVector<double>& corr, const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary, double>& vconst);
244
246 void saveHeatDensities();
247
249 inline void addCurvature(double& k44,
250 double& k33,
251 double& k22,
252 double& k11,
253 double& k43,
254 double& k21,
255 double& k42,
256 double& k31,
257 double& k32,
258 double& k41,
259 double ky,
260 double width,
261 const Vec<2, double>& midpoint);
262
264 template <CalcType calctype>
265 void setMatrix(FemMatrix& A,
266 DataVector<double>& B,
267 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary, double>& bvoltage);
268
271 struct ActiveRegionInfo {
272 shared_ptr<StackContainer<2>> layers;
273 Vec<2> origin;
274
275 ActiveRegionInfo(Vec<2> origin) : layers(plask::make_shared<StackContainer<2>>()), origin(origin) {}
276
278 size_t size() const { return layers->getChildrenCount(); }
279
281 shared_ptr<Material> getLayerMaterial(size_t n) const {
282 auto block = static_cast<Block<2>*>(static_cast<Translation<2>*>(layers->getChildNo(n).get())->getChild().get());
283 if (auto m = block->singleMaterial()) return m;
284 throw plask::Exception("freeCarrierGainSolver requires solid layers.");
285 }
286
288 Box2D getLayerBox(size_t n) const {
289 return static_cast<GeometryObjectD<2>*>(layers->getChildNo(n).get())->getBoundingBox() + origin;
290 }
291
293 bool isQW(size_t n) const { return static_cast<Translation<2>*>(layers->getChildNo(n).get())->getChild()->hasRole("QW"); }
294
296 Box2D getBoundingBox() const { return layers->getBoundingBox() + origin; }
297
299 bool contains(const Vec<2>& point) const { return getBoundingBox().contains(point); }
300
302 bool inQW(const Vec<2>& point) const {
303 if (!contains(point)) return false;
304 assert(layers->getChildForHeight(point.c1 - origin.c1));
305 return layers->getChildForHeight(point.c1 - origin.c1)->getChild()->hasRole("QW");
306 }
307
308 double averageNr(double lam, double T, double conc = 0.) const {
309 double nr = 0.;
310 for (size_t i = 0; i != materials.size(); ++i)
311 if (isQW(i)) nr += thicknesses[i] * materials[i]->Nr(lam, T, conc).real();
312 return nr / totalqw;
313 }
314
315 std::vector<shared_ptr<Material>> materials;
316 std::vector<double> thicknesses;
317 std::vector<size_t> wells;
318
319 double total;
320 double totalqw;
321 double bottom;
322 double top;
323
324 enum ConsideredHoles : unsigned {
325 NO_HOLES = 0,
326 HEAVY_HOLES = 1,
327 LIGHT_HOLES = 2,
328 BOTH_HOLES = 3
329 } holes;
330
335 void summarize(const DriftDiffusionModel2DSolver<Geometry2DType>* solver);
336 };
337
339 shared_ptr<Material> substrateMaterial;
340
342 std::vector<ActiveRegionInfo> regions;
343
345 void detectActiveRegions();
346
348 shared_ptr<RectangularMesh<2>> meshAct;
349
351 shared_ptr<RectangularMesh<2>> meshActMid;
352
354 double zSub;
355 std::vector<double> z;
356 double dz;
357 int nn;
358 int ne;
359 double hh2m;
360
361 // std::vector<double> CBel; /// vector with energy band diagram for electrons from CB (nm)
362 std::vector<double> CBelLev;
363
364 std::vector<double> lev_el;
365 int n_lev_el;
366 double CBelMin, CBelMax;
367 double T;
368 double Eshift;
369
370 public:
371 double maxerr;
372
374 BoundaryConditions<RectangularMesh<2>::Boundary, double> voltage_boundary;
375
376 typename ProviderFor<Potential, Geometry2DType>::Delegate outPotential;
377
378 typename ProviderFor<FermiLevels, Geometry2DType>::Delegate outFermiLevels;
379
380 typename ProviderFor<BandEdges, Geometry2DType>::Delegate outBandEdges;
381
382 typename ProviderFor<CurrentDensity, Geometry2DType>::Delegate outCurrentDensityForElectrons;
383
384 typename ProviderFor<CurrentDensity, Geometry2DType>::Delegate outCurrentDensityForHoles;
385
386 typename ProviderFor<CarriersConcentration, Geometry2DType>::Delegate outCarriersConcentration;
387
388 typename ProviderFor<Heat, Geometry2DType>::Delegate outHeat;
389 /*
390 typename ProviderFor<Conductivity, Geometry2DType>::Delegate outConductivity;
391 */
392 ReceiverFor<Temperature, Geometry2DType> inTemperature;
393
394 bool mRsrh;
395 bool mRrad;
396 bool mRaug;
397 bool mPol;
398 bool mFullIon;
399
400 double mSchottkyP;
401 double mSchottkyN;
402
403 double maxerrPsiI;
405 double maxerrPsi0;
407 double maxerrPsi;
409 double maxerrFn;
411 double maxerrFp;
413 size_t loopsPsiI;
414 size_t loopsPsi0;
415 size_t loopsPsi;
416 size_t loopsFn;
417 size_t loopsFp;
418
423 double compute(unsigned loops = 1);
424
429 double findEnergyLevels();
430
435 bool checkWell(std::string _carrier);
436
443 double integrateCurrent(size_t vindex, bool onlyactive = false);
444
450 double getTotalCurrent(size_t nact = 0);
451
456 // double getTotalEnergy();// LP_09.2015
457
462 // double getCapacitance();// LP_09.2015
463
468 // double getTotalHeat();// LP_09.2015
469
471 // double getErr() const { return toterr; }
472
473 void loadConfiguration(XMLReader& source,
474 Manager& manager) override; // for solver configuration (see: *.xpl file with structures)
475
476 DriftDiffusionModel2DSolver(const std::string& name = "");
477
478 std::string getClassName() const override;
479
480 ~DriftDiffusionModel2DSolver();
481
482 protected:
483 const LazyData<double> getPotentials(shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method) const;
484
485 const LazyData<double> getFermiLevels(FermiLevels::EnumType what,
486 shared_ptr<const MeshD<2>> dest_mesh,
487 InterpolationMethod method) const;
488
489 const LazyData<double> getBandEdges(BandEdges::EnumType what, shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method);
490
491 const LazyData<double> getHeatDensities(shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method);
492
493 const LazyData<Vec<2>> getCurrentDensitiesForElectrons(shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method);
494
495 const LazyData<Vec<2>> getCurrentDensitiesForHoles(shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method);
496
497 const LazyData<double> getCarriersConcentration(CarriersConcentration::EnumType what,
498 shared_ptr<const MeshD<2>> dest_mesh,
499 InterpolationMethod method);
500
501 /*
502 const LazyData<Tensor2<double>> getConductivity(shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod method);
503 */
504};
505
506}}} // namespace plask::electrical::drift_diffusion
507
508#endif