PLaSK library
Loading...
Searching...
No Matches
diffusion1d.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 "diffusion1d.hpp"
15
16#define dpbtrf_ F77_GLOBAL(dpbtrf, DPBTRF)
17#define dpbtrs_ F77_GLOBAL(dpbtrs, DPBTRS)
18
19F77SUB dpbtrf_(const char& uplo, const int& n, const int& kd, double* ab, const int& ldab, int& info);
20F77SUB dpbtrs_(const char& uplo,
21 const int& n,
22 const int& kd,
23 const int& nrhs,
24 double* ab,
25 const int& ldab,
26 double* b,
27 const int& ldb,
28 int& info);
29
30namespace plask { namespace electrical { namespace diffusion1d {
31
32constexpr double inv_hc = 1.0e-9 / (phys::c * phys::h_J);
33using phys::Z0;
34
35template <typename Geometry2DType>
36void DiffusionFem2DSolver<Geometry2DType>::loadConfiguration(XMLReader& reader, Manager& manager) {
37 while (reader.requireTagOrEnd()) {
38 const std::string& param = reader.getNodeName();
39
40 if (param == "config") {
41 fem_method = reader.enumAttribute<FemMethod>("fem-method")
42 .value("linear", FEM_LINEAR)
43 .value("parabolic", FEM_PARABOLIC)
44 .get(fem_method);
45 relative_accuracy = reader.getAttribute<double>("accuracy", relative_accuracy);
46 minor_concentration = reader.getAttribute<double>("abs-accuracy", minor_concentration);
47 interpolation_method = reader.getAttribute<InterpolationMethod>("interpolation", interpolation_method);
48 max_mesh_changes = reader.getAttribute<int>("maxrefines", max_mesh_changes);
49 max_iterations = reader.getAttribute<int>("maxiters", max_iterations);
50 do_initial = reader.getAttribute<bool>("initial", do_initial);
51
52 reader.requireTagEnd();
53 } else if (param == "mesh" && (reader.hasAttribute("start") || reader.hasAttribute("stop") || reader.hasAttribute("num"))) {
54 this->writelog(LOG_WARNING,
55 "Setting mesh directly in the solver is obsolete. Please use 'ref' attribute to refer to the mesh "
56 "defined in the <grids> section");
57 double r_min = reader.requireAttribute<double>("start");
58 double r_max = reader.requireAttribute<double>("stop");
59 size_t no_points = reader.requireAttribute<size_t>("num");
60 this->mesh = plask::make_shared<RegularMesh1D>(r_min, r_max, no_points);
61 reader.requireTagEnd();
62 } else
63 this->parseStandardConfiguration(reader, manager);
64 }
65}
66
67template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::compute_initial() { compute(COMPUTATION_INITIAL); }
68
69template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::compute_threshold() {
70 compute(COMPUTATION_THRESHOLD);
71}
72
73template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::compute_overthreshold() {
74 compute(COMPUTATION_OVERTHRESHOLD);
75}
76
77template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::onInitialize() {
78 iterations = 0;
79 detected_QW = detectQuantumWells();
80 determineQwWidth();
81
82 // reset mesh to original value
83 z = getZQWCoordinate();
84 if (!this->mesh) {
85 double left = INFINITY, right = -INFINITY;
86 for (const auto& box : detected_QW) {
87 if (box.left() < left) left = box.left();
88 if (box.right() > left) right = box.right();
89 }
90 size_t num = static_cast<size_t>(std::round((right - left) * 100)) + 1;
91 this->writelog(LOG_DETAIL, "Making default mesh with {} points", num);
92 this->setMesh(make_shared<RegularMesh1D>(left, right, num));
93 }
94 mesh2->setAxis0(this->mesh);
95 mesh2->setAxis1(plask::make_shared<plask::RegularAxis>(z, z, 1));
96 if (current_mesh().size() % 2 == 0)
97 current_mesh().reset(current_mesh().first(), current_mesh().last(), current_mesh().size() + 1);
98
99 n_present.reset(current_mesh().size(), 0.0);
100}
101
102template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::onInvalidate() {
103 j_on_the_mesh.reset();
104 T_on_the_mesh.reset();
105
106 PM.reset();
107 overthreshold_dgdn.reset();
108 // overthreshold_g.reset();
109
110 n_present.reset();
111 n_previous.reset();
112}
113
114template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::compute(ComputationType type) {
115 initial_computation = type == COMPUTATION_INITIAL || (this->initCalculation() && do_initial);
116 threshold_computation = type == COMPUTATION_THRESHOLD;
117 overthreshold_computation = type == COMPUTATION_OVERTHRESHOLD;
118
119 this->writelog(LOG_INFO, "Computing lateral carriers diffusion using {0} FEM method",
120 fem_method == FEM_LINEAR ? "linear" : "parabolic");
121
122 T_on_the_mesh = inTemperature(mesh2, interpolation_method); // data temperature vector provided by inTemperature reciever
123 j_on_the_mesh =
124 inCurrentDensity(mesh2, interpolation_method); // data current density vector provided by (in|out)Current reciever
125
126 int mesh_changes = 0;
127 bool convergence = true;
128
129 auto& current_mesh = this->current_mesh();
130 do {
131 if (!convergence) {
132 mesh_changes += 1;
133 if (mesh_changes > max_mesh_changes)
134 throw ComputationError(this->getId(), "maximum number of mesh refinements ({0}) reached", max_mesh_changes);
135 size_t new_size = 2 * current_mesh.size() - 1;
136 writelog(LOG_DETAIL, "Refining mesh (new size: {0})", new_size);
137
138 plask::DataVector<double> old_n = n_present;
139 size_t nm = old_n.size() - 1;
140
141 current_mesh.reset(current_mesh.first(), current_mesh.last(), new_size); // this calls invalidate
142 T_on_the_mesh =
143 inTemperature(mesh2, interpolation_method); // data temperature vector provided by inTemperature reciever
144 j_on_the_mesh =
145 inCurrentDensity(mesh2, interpolation_method); // data current density vector provided by inCurrentDensity reciever
146
147 n_present.reset(current_mesh.size(), 0.0);
148 // n in new points is average of the surrounding ones
149
150 for (size_t i = 0; i != nm; ++i) {
151 n_present[2*i] = old_n[i];
152 n_present[2*i+1] = 0.5 * (old_n[i] + old_n[i+1]);
153 }
154 n_present[2*nm] = old_n[nm];
155 }
156 if (initial_computation) {
157 this->writelog(LOG_DETAIL, "Conducting initial computations");
158 convergence = MatrixFEM();
159 if (convergence) initial_computation = false;
160 }
161 if (threshold_computation) {
162 this->writelog(LOG_DETAIL, "Conducting threshold computations");
163 convergence = MatrixFEM();
164 if (convergence) threshold_computation = false;
165 }
166 if (overthreshold_computation) {
167 // write_debug("Git0a!");
168
169 this->writelog(LOG_DETAIL, "Conducting overthreshold computations");
170 convergence = MatrixFEM();
171 if (convergence) overthreshold_computation = false;
172 }
173 } while (initial_computation || threshold_computation || overthreshold_computation);
174
175 this->writelog(LOG_DETAIL, "Converged after {0} mesh refinements and {1} computational loops", mesh_changes, iterations);
176}
177
178template <typename Geometry2DType> bool DiffusionFem2DSolver<Geometry2DType>::MatrixFEM() {
179 // Computation of K*n" - E*n = -F
180 bool _convergence;
181 iterations = 0;
182
183 // LAPACK factorization (dpbtrf) and equation solver (dpbtrs) info variables:
184 int info_f = 0;
185 int info_s = 0;
186
187 double A = 0.0;
188 double B = 0.0;
189 double C = 0.0;
190
191 double L = 0.0;
192 double R = 0.0;
193
194 // equation AX = RHS
195
196 DataVector<double> A_matrix;
197 DataVector<double> B_vector;
198
199 do {
200 if (fem_method == FEM_LINEAR)
201 A_matrix.reset(2 * current_mesh().size(), 0.0);
202 else if (fem_method == FEM_PARABOLIC)
203 A_matrix.reset(3 * current_mesh().size(), 0.0);
204
205 B_vector.reset(current_mesh().size(), 0.0);
206
207 n_previous = n_present.copy();
208
209 if (initial_computation) {
210 double T = 0.0;
211
212 for (std::size_t i = 0; i < current_mesh().size(); i++) {
213 T = T_on_the_mesh[i];
214
215 A = this->QW_material->A(T);
216 B = this->QW_material->B(T);
217 C = this->QW_material->C(T);
218
219 double RS = rightSide(i); // instead of rightSide(i, T, 0) ?
220 B_vector[i] =
221 pow((sqrt(27 * C * C * RS * RS + (4 * B * B * B - 18 * A * B * C) * RS + 4 * A * A * A * C - A * A * B * B) /
222 (2 * pow(3.0, 3. / 2.) * C * C) +
223 (-27 * C * C * RS + 9 * A * B * C - 2 * B * B * B) / (54 * C * C * C)),
224 1. / 3.) -
225 (3 * A * C - B * B) / pow(9 * C * C *
226 (sqrt(27 * C * C * RS * RS + (4 * B * B * B - 18 * A * B * C) * RS +
227 4 * A * A * A * C - A * A * B * B) /
228 (2 * pow(3.0, 3. / 2.) * C * C) +
229 (-27 * C * C * RS + 9 * A * B * C - 2 * B * B * B) / (54 * C * C * C)),
230 1. / 3.) -
231 B / (3 * C);
232 // X_vector[i] = -B/(3.0*C) - (pow(2.0,1.0/3.0)*(- B*B + 3.0*A*C))/(3.0*C*pow(-2.0*B*B*B +
233 // 9.0*A*B*C + 27.0*C*C*RS + sqrt(4.0*pow(- B*B + 3.0*A*C,3.0) +
234 // pow(-2.0*B*B*B + 9.0*A*B*C + 27.0*C*C*RS,2.0)),1.0/3.0)) +
235 // pow(-2.0*B*B*B + 9.0*A*B*C + 27.0*C*C*RS + sqrt(4.0*pow(-B*B +
236 // 3.0*A*C,3.0) + pow(-2.0*B*B*B + 9.0*A*B*C +
237 // 27.0*C*C*RS,2.0)),1.0/3.0)/(3.0*pow(2.0,1.0/3.0)*C);
238 // double X_part = 9*C*sqrt(27*C*C*rightSide(i)*rightSide(i) + (4*B*B*B-18*A*B*C)*rightSide(i) +
239 // 4*A*A*A*C - A*A*B*B); X_vector[i] = (pow(X_part, 2./3.) -
240 // pow(2,1./3.)*pow(3,1./6.)*B*pow(X_part, 1./3.) - pow(2,2./3.)*pow(3,4./3.)*A*C +
241 // pow(2,2./3.)*pow(3,1./3.)*B*B)/(pow(2,1./3.)*pow(3,7./6.)*C*pow(X_part, 1./3.)); double X_part =
242 // pow(sqrt(27*C*C*rightSide(i)*rightSide(i) + (4*B*B*B-18*A*B*C)*rightSide(i) + 4*A*A*A*C -
243 // A*A*B*B)/(2*pow(3,3/2)*C*C) + (-27*C*C*rightSide(i) + 9*A*B*C - 2*B*B*B)/(54*C*C*C),1/3);
244 // X_vector[i] = X_part - (3*A*C - B*B)/X_part - B/(3*C);
245 }
246 _convergence = true;
247 n_present = B_vector.copy();
248 } else {
249 const char uplo = 'U';
250 int lapack_n = 0;
251 int lapack_kd = 0;
252 int lapack_ldab = 0;
253 int lapack_nrhs = 0;
254 int lapack_ldb = 0;
255
256 if (overthreshold_computation) {
257 // Compute E and F components for overthreshold computations
258 if (inWavelength.size() != inLightE.size())
259 throw BadInput(this->getId(), "number of modes in inWavelength ({}) and inLightE ({}) differ",
260 inWavelength.size(), inLightE.size());
261
262 // Sum all modes
263 PM = DataVector<double>(mesh2->size(), 0.);
264 overthreshold_dgdn = DataVector<double>(mesh2->size(), 0.);
265 // overthreshold_g = DataVector<double>(mesh2->size(), 0.);
266 modesP.assign(inWavelength.size(), 0.);
267
268 double D;
269 if (current_mesh().size() > 1) D = current_mesh()[1] - current_mesh()[0]; // only ok for regular mesh
270
271 for (size_t n = 0; n != inWavelength.size(); ++n) {
272 double wavelength = real(inWavelength(n));
273 write_debug("wavelength: {0} nm", wavelength);
274
275 auto mesh_Li = plask::make_shared<plask::RectangularMesh<2>>();
276
277 mesh_Li->setAxis0(current_mesh_ptr());
278 mesh_Li->setAxis1(plask::make_shared<plask::OrderedAxis>(getZQWCoordinates()));
279
280 // auto Li = inLightMagnitude(n, mesh2, interpolation_method);
281 auto Li = averageLi(inLightE(n, mesh_Li, interpolation_method), *mesh_Li);
282
283 write_debug("Li[0]: {0} W/cm2", Li[0] * 1e-4);
284 int ile = 0;
285 for (auto n : n_present) {
286 if (n <= 0.) ile++;
287 }
288 // write_debug("n < 0: {0} times", ile);
289 auto g = inGain(mesh2, wavelength, interpolation_method);
290 // write_debug("g[0]: {0} cm(-1)", g[0]);
291 auto dgdn = inGain(Gain::DGDN, mesh2, wavelength, interpolation_method);
292 // write_debug("dgdn[0]: {0} cm(-4)", dgdn[0]);
293 auto factor = inv_hc * wavelength * 1e-4; // inverse one photon energy
294 for (size_t i = 0; i != mesh2->size(); ++i) {
295 auto common = factor * this->QW_material->nr(wavelength, T_on_the_mesh[i]) * Li[i];
296 double pm = common.c00 * g[i].c00 + common.c11 * g[i].c11;
297 PM[i] += pm;
298 overthreshold_dgdn[i] += common.c00 * dgdn[i].c00 + common.c11 * dgdn[i].c11;
299 ;
300 // overthreshold_g[i] += common * g[i];
301 modesP[n] += pm * D * jacobian(current_mesh()[i]);
302 }
303 modesP[n] *= 1.0e-5 * global_QW_width * (1e9 * phys::h_J * phys::c / wavelength);
304 // 1.0e-5 from µm to cm conversion and conversion to mW (r*dr), (...) - photon energy
305 }
306 }
307
308 DiffusionFem2DSolver<Geometry2DType>::createMatrices(A_matrix, B_vector);
309
310 if (fem_method == FEM_LINEAR) {
311 lapack_n = lapack_ldb = (int)current_mesh().size();
312 lapack_kd = 1;
313 lapack_ldab = 2;
314 lapack_nrhs = 1;
315 } else if (fem_method == FEM_PARABOLIC) {
316 DiffusionFem2DSolver<Geometry2DType>::createMatrices(A_matrix, B_vector);
317
318 lapack_n = lapack_ldb = (int)current_mesh().size();
319 lapack_kd = 2;
320 lapack_ldab = 3;
321 lapack_nrhs = 1;
322 }
323
324 dpbtrf_(uplo, lapack_n, lapack_kd, A_matrix.begin(), lapack_ldab, info_f); // factorize A
325 dpbtrs_(uplo, lapack_n, lapack_kd, lapack_nrhs, A_matrix.begin(), lapack_ldab, B_vector.begin(), lapack_ldb,
326 info_s); // solve Ax = B
327
328 n_present = B_vector.copy();
329
330 double absolute_error = 0.0;
331 double relative_error = 0.0;
332 double absolute_concentration_error =
333 abs(this->QW_material->A(300.0) * minor_concentration +
334 this->QW_material->B(300.0) * minor_concentration * minor_concentration +
335 this->QW_material->C(300.0) * minor_concentration * minor_concentration * minor_concentration);
336
337 /****************** RPSFEM: ******************/
338 // double tolerance = 5e+15;
339 // double n_tolerance = this->QW_material->A(300)*tolerance + this->QW_material->B(300)*tolerance*tolerance
340 // + this->QW_material->C(300)*tolerance*tolerance*tolerance;
341 /****************** end RPSFEM: ******************/
342
343 if (fem_method == FEM_LINEAR) {
344 _convergence = true;
345 for (std::size_t i = 0; i < current_mesh().size(); i++) {
346 L = leftSide(i);
347 R = rightSide(i);
348
349 absolute_error = L - R;
350 relative_error = absolute_error / R;
351 if (relative_accuracy < relative_error) _convergence = false;
352 }
353 } else if (fem_method == FEM_PARABOLIC) {
354 // double max_error_absolute = 0.0;
355 double max_error_relative = 0.0;
356 // int max_error_point = 0.0;
357 // double max_error_R = 0.0;
358
359 _convergence = true;
360 for (std::size_t i = 0; std::ptrdiff_t(i) < (std::ptrdiff_t(current_mesh().size()) - 1) / 2; i++) {
361 L = leftSide(2 * i + 1);
362 R = rightSide(2 * i + 1);
363
364 absolute_error = abs(L - R);
365 relative_error = abs(absolute_error / R);
366 // write_debug("absolute error: {0}", absolute_error);
367 // write_debug("relative error: {0}", relative_error);
368 // write_debug("n_present[0]: {0}", n_present[0]);
369 // write_debug("F(0): {0}", F(0));
370
371 if (max_error_relative < relative_error) {
372 max_error_relative = relative_error;
373 // max_error_absolute = absolute_error;
374 // max_error_point = current_mesh()[2*i+1];
375 // max_error_R = R;
376 }
377
378 if ((relative_accuracy < relative_error) && (absolute_concentration_error < absolute_error)) {
379 _convergence = false;
380 break;
381 }
382 }
383 }
384 iterations += 1;
385#ifndef NDEBUG
386 writelog(LOG_DEBUG, "iteration: {0}", iterations);
387 if (overthreshold_computation)
388 writelog(LOG_DEBUG, "Integral of overthreshold loses: {0} mW, qw_width: {1} cm", burning_integral(),
389 global_QW_width);
390#endif
391 }
392
393 } while (!_convergence && (iterations < max_iterations));
394
395 outCarriersConcentration.fireChanged();
396
397 return _convergence;
398}
399
400template <>
401void DiffusionFem2DSolver<Geometry2DCartesian>::createMatrices(DataVector<double> A_matrix, DataVector<double> B_vector) {
402 // linear FEM elements variables:
403
404 double r1 = 0.0, r2 = 0.0;
405 double k11e = 0, k12e = 0, k22e = 0;
406 double p1e = 0, p2e = 0;
407
408 // parabolic FEM elements variables:
409
410 double K = 0.0;
411 double E = 0.0;
412 double F = 0.0;
413
414 if (fem_method == FEM_LINEAR) // 02.10.2012 Marcin Gebski
415 {
416 double j1 = 0.0;
417 double j2 = 0.0;
418
419 auto& current_mesh = this->current_mesh();
420 for (int i = 0; i < int(current_mesh.size()) - 1; i++) // loop over all elements
421 {
422 r1 = current_mesh[i] * 1e-4;
423 r2 = current_mesh[i+1] * 1e-4;
424
425 j1 = abs(j_on_the_mesh[i][1] * 1e+3);
426 j2 = abs(j_on_the_mesh[i+1][1] * 1e+3);
427
428 K = this->K(i);
429 F = this->F(i);
430 E = this->E(i);
431
432 k11e = K / (r2 - r1) + E * (r2 - r1) / 3;
433 k12e = -K / (r2 - r1) + E * (r2 - r1) / 6;
434 k22e = K / (r2 - r1) + E * (r2 - r1) / 3;
435
436 p1e = ((r2 - r1) / 2) * (F + (2 * j1 + j2) / (3 * plask::phys::qe * global_QW_width));
437 p2e = ((r2 - r1) / 2) * (F + (2 * j2 + j1) / (3 * plask::phys::qe * global_QW_width));
438
439 A_matrix[2*i+1] += k11e;
440 A_matrix[2*i+2] += k12e;
441 A_matrix[2*i+3] += k22e;
442
443 B_vector[i] += p1e;
444 B_vector[i+1] += p2e;
445
446 } // end loop over all elements
447
448 } else if (fem_method == FEM_PARABOLIC) {
449 double r3 = 0.0;
450 double k13e = 0.0, k23e = 0.0, k33e = 0.0; // local stiffness matrix elements
451 double p3e = 0.0;
452
453 auto& current_mesh = this->current_mesh();
454 for (int i = 0; i < (int(current_mesh.size()) - 1) / 2; i++) // loop over all elements
455 {
456 r1 = current_mesh[2*i] * 1e-4;
457 r3 = current_mesh[2*i+2] * 1e-4;
458
459 K = this->K(2 * i + 1);
460 F = this->F(2 * i + 1);
461 E = this->E(2 * i + 1);
462
463 K = K / ((r3 - r1) * (r3 - r1));
464 double Cnst = (r3 - r1) / 30;
465
466 k11e = Cnst * (70 * K + 4 * E);
467 k12e = Cnst * (-80 * K + 2 * E); // = k21e
468 k13e = Cnst * (10 * K - E); // = k31e
469 k22e = Cnst * (160 * K + 16 * E);
470 k23e = k12e; // = k32e
471 k33e = Cnst * (70 * K + 4 * E);
472
473 p1e = F * (r3 - r1) / 6;
474 p2e = 2 * F * (r3 - r1) / 3;
475 p3e = p1e;
476
477 // Fill matrix A_matrix columnwise: //29.06.2012 r. Marcin Gebski
478
479 A_matrix[6*i+2] += k11e;
480 A_matrix[6*i+4] += k12e;
481 A_matrix[6*i+6] += k13e;
482 A_matrix[6*i+5] += k22e;
483 A_matrix[6*i+7] += k23e;
484 A_matrix[6*i+8] += k33e;
485 A_matrix[6*i+3] += 0; // k24 = 0 - fill top band
486
487 B_vector[2*i] += p1e;
488 B_vector[2*i+1] += p2e;
489 B_vector[2*i+2] += p3e;
490
491 } // end loop over all elements
492 }
493}
494
495template <>
496void DiffusionFem2DSolver<Geometry2DCylindrical>::createMatrices(DataVector<double> A_matrix, DataVector<double> B_vector) {
497 // linear FEM elements variables:
498
499 double r1 = 0.0, r2 = 0.0;
500 double k11e = 0, k12e = 0, k22e = 0;
501 double p1e = 0, p2e = 0;
502
503 // parabolic FEM elements variables:
504
505 double K = 0.0;
506 double E = 0.0;
507 double F = 0.0;
508
509 if (fem_method == FEM_LINEAR) // 02.10.2012 Marcin Gebski
510 {
511 double j1 = 0.0;
512 double j2 = 0.0;
513
514 auto& current_mesh = this->current_mesh();
515 for (int i = 0; i < int(current_mesh.size()) - 1; i++) // loop over all elements
516 {
517 r1 = current_mesh[i] * 1e-4;
518 r2 = current_mesh[i+1] * 1e-4;
519
520 j1 = abs(j_on_the_mesh[i][1] * 1e+3);
521 j2 = abs(j_on_the_mesh[i+1][1] * 1e+3);
522
523 K = this->K(i);
524 F = this->F(i);
525 E = this->E(i);
526
527 K = 4 * K / ((r2 - r1) * (r2 - r1));
528
529 k11e = PI * (r2 - r1) / 4 * ((K + E) * (r1 + r2) + E * (3 * r1 - r2) / 3);
530 k12e = PI * (r2 - r1) / 4 * ((-K + E) * (r1 + r2) - E * (r1 + r2) / 3);
531 k22e = PI * (r2 - r1) / 4 * ((K + E) * (r1 + r2) + E * (3 * r2 - r1) / 3);
532
533 p1e = PI * (r2 - r1) *
534 (F / 3 * (2 * r1 + r2) +
535 (1 / (6 * plask::phys::qe * global_QW_width)) * (3 * j1 * r1 + j1 * r2 + j2 * r1 + r2 * j2));
536 p2e = PI * (r2 - r1) *
537 (F / 3 * (r1 + 2 * r2) +
538 (1 / (6 * plask::phys::qe * global_QW_width)) * (3 * j2 * r2 + j1 * r2 + j2 * r1 + r1 * j1));
539
540 A_matrix[2*i+1] += k11e;
541 A_matrix[2*i+2] += k12e;
542 A_matrix[2*i+3] += k22e;
543
544 B_vector[i] += p1e;
545 B_vector[i+1] += p2e;
546
547 } // end loop over all elements
548 } else if (fem_method == FEM_PARABOLIC) // 02.10.2012 Marcin Gebski
549 {
550 double r3 = 0.0;
551 double k13e = 0.0, k23e = 0.0, k33e = 0.0; // local stiffness matrix elements
552 double p3e = 0.0;
553
554 auto& current_mesh = this->current_mesh();
555 for (int i = 0; i < (int(current_mesh.size()) - 1) / 2; i++) // loop over all elements
556 {
557 r1 = current_mesh[2*i] * 1e-4;
558 r3 = current_mesh[2*i+2] * 1e-4;
559
560 K = this->K(2 * i + 1);
561 F = this->F(2 * i + 1);
562 E = this->E(2 * i + 1);
563
564 double Cnst = PI * (r3 - r1) / 30;
565
566 k11e = Cnst * (10 * K * (11 * r1 + 3 * r3) / ((r3 - r1) * (r3 - r1)) + E * (7 * r1 + r3));
567 k12e = Cnst * (-40 * K * (3 * r1 + r3) / ((r3 - r1) * (r3 - r1)) + 4 * E * r1); // = k21e
568 k13e = Cnst * (10 * K * (r1 + r3) / ((r3 - r1) * (r3 - r1)) - E * (r1 + r3)); // = k31e
569 k22e = Cnst * (160 * K * (r1 + r3) / ((r3 - r1) * (r3 - r1)) + 16 * E * (r1 + r3));
570 k23e = Cnst * (-40 * K * (r1 + 3 * r3) / ((r3 - r1) * (r3 - r1)) + 4 * E * r3); // = k32e
571 k33e = Cnst * (10 * K * (3 * r1 + 11 * r3) / ((r3 - r1) * (r3 - r1)) + E * (r1 + 7 * r3));
572
573 p1e = Cnst * 10 * F * r1;
574 p2e = Cnst * 20 * F * (r1 + r3);
575 p3e = Cnst * 10 * F * r3;
576
577 // Fill matrix A_matrix columnwise: //29.06.2012 r. Marcin Gebski
578
579 A_matrix[6*i+2] += k11e;
580 A_matrix[6*i+4] += k12e;
581 A_matrix[6*i+6] += k13e;
582 A_matrix[6*i+5] += k22e;
583 A_matrix[6*i+7] += k23e;
584 A_matrix[6*i+8] += k33e;
585 A_matrix[6*i+3] += 0; // k24 = 0 - fill top band
586
587 B_vector[2*i] += p1e;
588 B_vector[2*i+1] += p2e;
589 B_vector[2*i+2] += p3e;
590
591 } // end loop over all elements
592 }
593}
594
595template <typename Geometry2DType>
596DiffusionFem2DSolver<Geometry2DType>::ConcentrationDataImpl::ConcentrationDataImpl(const DiffusionFem2DSolver* solver,
597 shared_ptr<const plask::MeshD<2>> dest_mesh,
598 InterpolationMethod interp)
599 : solver(solver),
600 destination_mesh(dest_mesh),
601 interpolationFlags(InterpolationFlags(solver->geometry)),
602 concentration(interpolate(solver->mesh2,
603 solver->n_present,
604 dest_mesh,
605 getInterpolationMethod<INTERPOLATION_SPLINE>(interp),
606 interpolationFlags)) {}
607
608template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::ConcentrationDataImpl::at(size_t i) const {
609 // Make sure we have concentration only in the quantum wells
610 // TODO maybe more optimal approach would be reasonable?
611 auto point = interpolationFlags.wrap(destination_mesh->at(i));
612 bool inqw = false;
613 for (auto QW : solver->detected_QW)
614 if (QW.contains(point)) {
615 inqw = true;
616 break;
617 }
618 if (!inqw) return 0.;
619 return concentration[i];
620}
621
622template <typename Geometry2DType>
623const LazyData<double> DiffusionFem2DSolver<Geometry2DType>::getConcentration(CarriersConcentration::EnumType what,
624 shared_ptr<const plask::MeshD<2>> dest_mesh,
625 InterpolationMethod interpolation) const {
626 if (what != CarriersConcentration::MAJORITY && what != CarriersConcentration::PAIRS) {
627 return LazyData<double>(dest_mesh->size(), NAN);
628 }
629 if (!n_present.data()) throw NoValue("carriers concentration");
630 return LazyData<double>(new DiffusionFem2DSolver<Geometry2DType>::ConcentrationDataImpl(this, dest_mesh, interpolation));
631}
632
633template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::K(int i) {
634 double T = T_on_the_mesh[i];
635 if (threshold_computation || overthreshold_computation)
636 return this->QW_material->D(T);
637 else
638 return 0.; // for initial distribution there is no diffusion
639}
640
641template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::E(int i) {
642 double T = T_on_the_mesh[i];
643 double n0 = n_previous[i];
644 double result = 0.0;
645
646 result = (this->QW_material->A(T) + 2 * this->QW_material->B(T) * n0 + 3 * this->QW_material->C(T) * n0 * n0);
647
648 if (overthreshold_computation) {
649 result += overthreshold_dgdn[i];
650 }
651
652 return result;
653}
654
655template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::F(int i) {
656 double T = T_on_the_mesh[i];
657 double n0 = n_previous[i];
658 double product = 0.0;
659
660 product = (this->QW_material->B(T) * n0 * n0 + 2 * this->QW_material->C(T) * n0 * n0 * n0);
661
662 if (fem_method == FEM_PARABOLIC) // 02.10.2012 Marcin Gebski
663 product += abs(j_on_the_mesh[i][1] * 1e+3) / (plask::phys::qe * global_QW_width);
664
665 if (overthreshold_computation) {
666 product += overthreshold_dgdn[i] * n0 - PM[i];
667 }
668
669 return product;
670}
671
672template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::nSecondDeriv(std::size_t i) {
673 double n_second_deriv = 0.0; // second derivative with respect to r
674 double dr = 0.0;
675
676 auto& current_mesh = this->current_mesh();
677 if (fem_method == FEM_LINEAR) // 02.10.2012 Marcin Gebski
678 {
679 dr = (current_mesh.last() - current_mesh.first()) * 1e-4 / (double)current_mesh.size();
680 double n_right = 0, n_left = 0, n_central = 0; // n values for derivative: right-side, left-side, central
681
682 if ((i > 0) && (i + 1 < current_mesh.size())) // middle of the range
683 {
684 n_right = n_present[i+1];
685 n_left = n_present[i-1];
686 n_central = n_present[i];
687
688 if (std::is_same<Geometry2DType, Geometry2DCylindrical>::value)
689 n_second_deriv =
690 (n_right - 2 * n_central + n_left) / (dr * dr) + 1.0 / (current_mesh[i] * 1e-4) * (n_right - n_left) / (2 * dr);
691 else
692 n_second_deriv = (n_right - 2 * n_central + n_left) / (dr * dr); // + 1.0/(current_mesh()[i]*1e-4);
693 } else if (i == 0) // punkt r = 0
694 {
695 n_right = n_present[i+1];
696 n_left = n_present[i+1]; // w r = 0 jest maksimum(warunki brzegowe), zalozylem nP = nL w blisko maksimum
697 n_central = n_present[i];
698
699 n_second_deriv = 2 * (n_right - 2 * n_central + n_left) / (dr * dr); // wymyslone przez Wasiak 2011.02.07
700 } else // prawy brzeg
701 {
702 n_right = n_present[i-1];
703 n_left = n_present[i-1]; // podobnie jak we wczesniejszym warunku
704 n_central = n_present[i];
705
706 if (std::is_same<Geometry2DType, Geometry2DCylindrical>::value)
707 n_second_deriv =
708 (n_right - 2 * n_central + n_left) / (dr * dr) + 1.0 / (current_mesh[i] * 1e-4) * (n_right - n_left) / (2 * dr);
709 else
710 n_second_deriv = (n_right - 2 * n_central + n_left) / (dr * dr); // + 1.0/(current_mesh()[i]*1e-4);
711 }
712 } else if (fem_method == FEM_PARABOLIC) // 02.10.2012 Marcin Gebski
713 {
714 dr = (current_mesh[i+1] - current_mesh[i-1]) * 1e-4;
715 n_second_deriv = (n_present[i-1] + n_present[i+1] - 2.0 * n_present[i]) * (4.0 / (dr * dr));
716 if (std::is_same<Geometry2DType, Geometry2DCylindrical>::value)
717 n_second_deriv += (1.0 / (current_mesh[i] * 1e-4)) * (1.0 / dr) *
718 (n_present[i+1] - n_present[i-1]); // adding cylindrical component of laplace operator
719 }
720
721 return n_second_deriv;
722}
723
724template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::leftSide(std::size_t i) {
725 double T = T_on_the_mesh[i];
726 double n = n_present[i];
727
728 double product = -(this->QW_material->A(T) * n + this->QW_material->B(T) * n * n + this->QW_material->C(T) * n * n * n);
729
730 if (threshold_computation || overthreshold_computation) product += this->QW_material->D(T) * nSecondDeriv(i);
731
732 if (overthreshold_computation) {
733 // write_debug("product before: {0}", product);
734 product -= PM[i];
735 // write_debug("product after: {0}", product);
736 // write_debug("PM[i]: {0}", PM[i]);
737 }
738
739 return product;
740}
741
742template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::burning_integral() {
743 if (modesP.size() == 0)
744 throw Exception("{0}: You must run over-threshold computations first before getting burring integral.", this->getId());
745 double int_val = 0.0;
746 for (double p : modesP) int_val += p;
747 return int_val;
748}
749
750template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::rightSide(std::size_t i) {
751 return -abs(j_on_the_mesh[i][1]) * 1e+3 / (plask::phys::qe * global_QW_width);
752}
753
754template <typename Geometry2DType> std::vector<Box2D> DiffusionFem2DSolver<Geometry2DType>::detectQuantumWells() {
755 if (!this->geometry) throw NoGeometryException(this->getId());
756
757 shared_ptr<RectangularMesh<2>> grid =
758 RectangularMesh2DSimpleGenerator().generate_t<RectangularMesh<2>>(this->geometry->getChild());
759 shared_ptr<RectangularMesh<2>> points = grid->getElementMesh();
760
761 std::vector<Box2D> results;
762
763 // Compact each row (it can contain only one QW and each must start and end in the same point)
764 double left = 0., right = 0.;
765 bool foundQW = false;
766 bool had_active = false, after_active = false;
767 for (std::size_t r = 0; r < points->axis[1]->size(); ++r) {
768 bool inQW = false;
769 bool active_row = false;
770 for (std::size_t c = 0; c < points->axis[0]->size(); ++c) {
771 auto point = points->at(c, r);
772 auto tags = this->geometry->getRolesAt(point);
773 bool QW = tags.find("QW") != tags.end() || tags.find("QD") != tags.end();
774 bool active = false;
775 for (const std::string& tag : tags)
776 if (tag.substr(0, 6) == "active") {
777 active = true;
778 break;
779 }
780 if (QW && !active) throw Exception("{0}: All marked quantum wells must belong to marked active region.", this->getId());
781 if (QW && !inQW) // QW start
782 {
783 if (foundQW) {
784 if (left != grid->axis[0]->at(c))
785 throw Exception("{0}: Left edge of quantum wells not vertically aligned.", this->getId());
786 if (*this->geometry->getMaterial(point) != *QW_material)
787 throw Exception("{0}: Quantum wells of multiple materials not supported.", this->getId());
788 } else {
789 QW_material = this->geometry->getMaterial(point);
790 }
791 left = grid->axis[0]->at(c);
792 inQW = true;
793 }
794 if (!QW && inQW) // QW end
795 {
796 if (foundQW && right != grid->axis[0]->at(c))
797 throw Exception("{0}: Right edge of quantum wells not vertically aligned.", this->getId());
798 right = grid->axis[0]->at(c);
799 results.push_back(Box2D(left, grid->axis[1]->at(r), right, grid->axis[1]->at(r + 1)));
800 foundQW = true;
801 inQW = false;
802 }
803 if (active) {
804 active_row = had_active = true;
805 if (after_active) throw Exception("{0}: Multiple active regions not supported.", this->getId());
806 }
807 }
808 if (inQW) { // handle situation when QW spans to the end of the structure
809 if (foundQW && right != grid->axis[0]->at(points->axis[0]->size()))
810 throw Exception("{0}: Right edge of quantum wells not vertically aligned.", this->getId());
811 right = grid->axis[0]->at(points->axis[0]->size());
812 results.push_back(Box2D(left, grid->axis[1]->at(r), right, grid->axis[1]->at(r + 1)));
813 foundQW = true;
814 inQW = false;
815 }
816 if (!active_row && had_active) after_active = true;
817 }
818
819 // Compact results in vertical direction
820 // TODO
821
822 return results;
823}
824
825template <typename Geometry2DType> double DiffusionFem2DSolver<Geometry2DType>::getZQWCoordinate() {
826 double coordinate = 0.0;
827 std::size_t no_QW = detected_QW.size();
828 if (no_QW == 0) throw Exception("no quantum wells defined");
829 // here no_QW > 0:
830 std::size_t no_Box;
831 if (no_QW % 2 == 0) {
832 no_Box = no_QW / 2 - 1;
833 coordinate = (detected_QW[no_Box].lower[1] + detected_QW[no_Box].upper[1]) / 2.0;
834 } else {
835 no_Box = (no_QW - 1) / 2;
836 coordinate = (detected_QW[no_Box].lower[1] + detected_QW[no_Box].upper[1]) / 2.0;
837 }
838 return coordinate;
839}
840
841template <typename Geometry2DType> std::vector<double> DiffusionFem2DSolver<Geometry2DType>::getZQWCoordinates() {
842 std::size_t no_QW = detected_QW.size();
843 if (no_QW == 0) throw Exception("no quantum wells defined");
844
845 std::vector<double> coordinates(no_QW);
846
847 for (std::size_t i = 0; i < no_QW; i++) coordinates[i] = (detected_QW[i].lower[1] + detected_QW[i].upper[1]) / 2.0;
848
849 return coordinates;
850}
851
852template <typename Geometry2DType> void DiffusionFem2DSolver<Geometry2DType>::determineQwWidth() {
853 global_QW_width = 0.;
854 for (std::size_t i = 0; i < detected_QW.size(); i++) {
855 global_QW_width += (this->detected_QW[i].upper[1] - this->detected_QW[i].lower[1]);
856 }
857 global_QW_width *= 1e-4;
858}
859
860template <typename Geometry2DType>
861plask::DataVector<const Tensor2<double>> DiffusionFem2DSolver<Geometry2DType>::averageLi(plask::LazyData<Vec<3, dcomplex>> Le,
862 const plask::RectangularMesh<2>& mesh_Li) {
863 plask::DataVector<Tensor2<double>> Li(current_mesh().size());
864
865 for (std::size_t i = 0; i < current_mesh().size(); ++i) {
866 Tensor2<double> current_Li(0., 0.);
867
868 for (std::size_t j = 0; j < detected_QW.size(); ++j) {
869 std::size_t k = mesh_Li.index(i, j);
870 current_Li +=
871 Tensor2<double>(real(Le[k].c0 * conj(Le[k].c0) + Le[k].c1 * conj(Le[k].c1)), real(Le[k].c2 * conj(Le[k].c2)));
872 }
873 Li[i] = 0.5 / Z0 * current_Li / double(detected_QW.size());
874 }
875 return Li;
876}
877
878template <> std::string DiffusionFem2DSolver<Geometry2DCartesian>::getClassName() const { return "electrical.OldDiffusion2D"; }
879template <> std::string DiffusionFem2DSolver<Geometry2DCylindrical>::getClassName() const { return "electrical.OldDiffusionCyl"; }
880
881template class PLASK_SOLVER_API DiffusionFem2DSolver<Geometry2DCartesian>;
882template class PLASK_SOLVER_API DiffusionFem2DSolver<Geometry2DCylindrical>;
883
884}}} // namespace plask::electrical::diffusion1d