PLaSK library

Loading...
Searching...
No Matches
electr2d.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 "electr2d.hpp"
15
16 namespace plask { namespace electrical { namespace shockley {
17
18 template <typename Geometry2DType>
19 ElectricalFem2DSolver<Geometry2DType>::ElectricalFem2DSolver(const std::string& name)
20 : FemSolverWithMaskedMesh<Geometry2DType, RectangularMesh<2>>(name),
21 pcond(5.),
22 ncond(50.),
23 loopno(0),
24 default_junction_conductivity(Tensor2<double>(0., 5.)),
25 maxerr(0.05),
26 outVoltage(this, &ElectricalFem2DSolver<Geometry2DType>::getVoltage),
27 outCurrentDensity(this, &ElectricalFem2DSolver<Geometry2DType>::getCurrentDensities),
28 outHeat(this, &ElectricalFem2DSolver<Geometry2DType>::getHeatDensities),
29 outConductivity(this, &ElectricalFem2DSolver<Geometry2DType>::getConductivity),
30 convergence(CONVERGENCE_FAST) {
31 onInvalidate();
32 inTemperature = 300.;
33 junction_conductivity.reset(1, default_junction_conductivity);
34}
35
36 template <typename Geometry2DType>
37 void ElectricalFem2DSolver<Geometry2DType>::loadConfiguration(XMLReader& source, Manager& manager) {
38 while (source.requireTagOrEnd()) parseConfiguration(source, manager);
39}
40
41 template <typename Geometry2DType>
42 void ElectricalFem2DSolver<Geometry2DType>::parseConfiguration(XMLReader& source, Manager& manager) {
43 std::string param = source.getNodeName();
44
45 if (param == "potential") source.throwException("<potential> boundary conditions have been permanently renamed to <voltage>");
46
47 if (param == "voltage")
48 this->readBoundaryConditions(manager, source, voltage_boundary);
49
50 else if (param == "loop") {
51 if (source.hasAttribute("start-cond") || source.hasAttribute("start-cond-inplane")) {
52 double c0 = default_junction_conductivity.c00, c1 = default_junction_conductivity.c11;
53 auto vert = source.getAttribute<std::string>("start-cond");
54 if (vert) {
55 if (vert->find(',') == std::string::npos) {
56 c1 = boost::lexical_cast<double>(*vert);
57 } else {
58 if (source.hasAttribute("start-cond-inplane"))
59 throw XMLException(
60 source, "tag attribute 'start-cond' has two values, but attribute 'start-cond-inplane' is also provided");
61 auto values = splitString2(*vert, ',');
62 c0 = boost::lexical_cast<double>(values.first);
63 c1 = boost::lexical_cast<double>(values.second);
64 }
65 }
66 c0 = source.getAttribute<double>("start-cond-inplane", c0);
67 this->setCondJunc(Tensor2<double>(c0, c1));
68 }
69 convergence = source.enumAttribute<Convergence>("convergence")
70 .value("fast", CONVERGENCE_FAST)
71 .value("stable", CONVERGENCE_STABLE)
72 .get(convergence);
73 maxerr = source.getAttribute<double>("maxerr", maxerr);
74 source.requireTagEnd();
75 }
76
77 else if (param == "contacts") {
78 pcond = source.getAttribute<double>("pcond", pcond);
79 ncond = source.getAttribute<double>("ncond", ncond);
80 source.requireTagEnd();
81 }
82
83 else if (!this->parseFemConfiguration(source, manager)) {
84 this->parseStandardConfiguration(source, manager);
85 }
86}
87
88 template <typename Geometry2DType> ElectricalFem2DSolver<Geometry2DType>::~ElectricalFem2DSolver() {}
89
90 template <typename Geometry2DType> void ElectricalFem2DSolver<Geometry2DType>::setupActiveRegions() {
91 this->invalidate();
92
93 if (!this->geometry || !this->mesh) {
94 if (junction_conductivity.size() != 1) {
95 Tensor2<double> condy(0., 0.);
96 for (auto cond : junction_conductivity) condy += cond;
97 junction_conductivity.reset(1, condy / double(junction_conductivity.size()));
98 }
99 return;
100 }
101
102 this->setupMaskedMesh();
103
104 auto points = this->mesh->getElementMesh();
105
106 std::vector<typename Active::Region> regions;
107
108 for (size_t r = 0; r < points->axis[1]->size(); ++r) {
109 size_t prev = 0;
110 shared_ptr<Material> material;
111 for (size_t c = 0; c < points->axis[0]->size(); ++c) { // In the (possible) active region
112 auto point = points->at(c, r);
113 size_t num = isActive(point);
114
115 if (num) { // here we are inside the active region
116 if (regions.size() >= num && regions[num - 1].warn) {
117 if (!material)
118 material = this->geometry->getMaterial(points->at(c, r));
119 else if (*material != *this->geometry->getMaterial(points->at(c, r))) {
120 writelog(LOG_WARNING, "Junction {} is laterally non-uniform", num - 1);
121 regions[num - 1].warn = false;
122 }
123 }
124 regions.resize(max(regions.size(), num));
125 auto& reg = regions[num - 1];
126 if (prev != num) { // this region starts in the current row
127 if (reg.top < r) {
128 throw Exception("{0}: Junction {1} is disjoint", this->getId(), num - 1);
129 }
130 if (reg.bottom >= r)
131 reg.bottom = r; // first row
132 else if (reg.rowr <= c)
133 throw Exception("{0}: Junction {1} is disjoint", this->getId(), num - 1);
134 reg.top = r + 1;
135 reg.rowl = c;
136 if (reg.left > reg.rowl) reg.left = reg.rowl;
137 }
138 }
139 if (prev && prev != num) { // previous region ended
140 auto& reg = regions[prev - 1];
141 if (reg.bottom < r && reg.rowl >= c) throw Exception("{0}: Junction {1} is disjoint", this->getId(), prev - 1);
142 reg.rowr = c;
143 if (reg.right < reg.rowr) reg.right = reg.rowr;
144 }
145 prev = num;
146 }
147 if (prev) // junction reached the edge
148 regions[prev - 1].rowr = regions[prev - 1].right = points->axis[0]->size();
149 }
150
151 size_t condsize = 0;
152 active.clear();
153 active.reserve(regions.size());
154 size_t i = 0;
155 for (auto& reg : regions) {
156 if (reg.bottom == std::numeric_limits<size_t>::max()) reg.bottom = reg.top = 0;
157 active.emplace_back(condsize, reg.left, reg.right, reg.bottom, reg.top,
158 this->mesh->axis[1]->at(reg.top) - this->mesh->axis[1]->at(reg.bottom));
159 condsize += reg.right - reg.left;
160 this->writelog(LOG_DETAIL, "Detected junction {0} thickness = {1}nm", i++, 1e3 * active.back().height);
161 this->writelog(LOG_DEBUG, "Junction {0} span: [{1},{3}]-[{2},{4}]", i - 1, reg.left, reg.right, reg.bottom, reg.top);
162 }
163
164 if (junction_conductivity.size() != condsize) {
165 Tensor2<double> condy(0., 0.);
166 for (auto cond : junction_conductivity) condy += cond;
167 junction_conductivity.reset(max(condsize, size_t(1)), condy / double(junction_conductivity.size()));
168 }
169}
170
171 template <typename Geometry2DType> void ElectricalFem2DSolver<Geometry2DType>::onInitialize() {
172 if (!this->geometry) throw NoGeometryException(this->getId());
173 if (!this->mesh) throw NoMeshException(this->getId());
174 setupActiveRegions();
175 loopno = 0;
176 potentials.reset(this->maskedMesh->size(), 0.);
177 currents.reset(this->maskedMesh->getElementsCount(), vec(0., 0.));
178 conds.reset(this->maskedMesh->getElementsCount());
179}
180
181 template <typename Geometry2DType> void ElectricalFem2DSolver<Geometry2DType>::onInvalidate() {
182 conds.reset();
183 potentials.reset();
184 currents.reset();
185 heats.reset();
186 junction_conductivity.reset(1, default_junction_conductivity);
187}
188
189 template <>
190 inline void ElectricalFem2DSolver<Geometry2DCartesian>::setLocalMatrix(double&,
191 double&,
192 double&,
193 double&,
194 double&,
195 double&,
196 double&,
197 double&,
198 double&,
199 double&,
200 double,
201 double,
202 const Vec<2, double>&) {
203 return;
204}
205
206 template <>
207 inline void ElectricalFem2DSolver<Geometry2DCylindrical>::setLocalMatrix(double& k44,
208 double& k33,
209 double& k22,
210 double& k11,
211 double& k43,
212 double& k21,
213 double& k42,
214 double& k31,
215 double& k32,
216 double& k41,
217 double,
218 double,
219 const Vec<2, double>& midpoint) {
220 double r = midpoint.rad_r();
221 k44 = r * k44;
222 k33 = r * k33;
223 k22 = r * k22;
224 k11 = r * k11;
225 k43 = r * k43;
226 k21 = r * k21;
227 k42 = r * k42;
228 k31 = r * k31;
229 k32 = r * k32;
230 k41 = r * k41;
231}
232
234 template <typename Geometry2DType>
235 void ElectricalFem2DSolver<Geometry2DType>::setMatrix(
236 FemMatrix& A,
237 DataVector<double>& B,
238 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary, double>& bvoltage,
239 const LazyData<double>& temperature) {
240 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
241
242 // Update junction conductivities
243 if (loopno != 0) {
244 for (auto e : this->maskedMesh->elements()) {
245 if (size_t nact = isActive(e)) {
246 size_t i = e.getIndex();
247 size_t left = this->maskedMesh->index0(e.getLoLoIndex());
248 size_t right = this->maskedMesh->index0(e.getUpLoIndex());
249 const Active& act = active[nact - 1];
250 double U =
251 0.5 *
252 (potentials[this->maskedMesh->index(left, act.top)] - potentials[this->maskedMesh->index(left, act.bottom)] +
253 potentials[this->maskedMesh->index(right, act.top)] - potentials[this->maskedMesh->index(right, act.bottom)]);
254 double jy = 0.1 * conds[i].c11 * U / act.height; // [j] = kA/cm²
255 size_t ti = this->maskedMesh->element(e.getIndex0(), (act.top + act.bottom) / 2).getIndex();
256 Tensor2<double> cond = activeCond(nact - 1, U, jy, temperature[ti]);
257 switch (convergence) {
258 case CONVERGENCE_STABLE:
259 cond = 0.5 * (conds[i] + cond);
260 case CONVERGENCE_FAST:
261 conds[i] = cond;
262 }
263 if (isnan(conds[i].c11) || abs(conds[i].c11) < 1e-16) conds[i].c11 = 1e-16;
264 }
265 }
266 }
267
268 A.clear();
269 B.fill(0.);
270
271 // Set stiffness matrix and load vector
272 for (auto e : this->maskedMesh->elements()) {
273 size_t i = e.getIndex();
274
275 // nodes numbers for the current element
276 size_t loleftno = e.getLoLoIndex();
277 size_t lorghtno = e.getUpLoIndex();
278 size_t upleftno = e.getLoUpIndex();
279 size_t uprghtno = e.getUpUpIndex();
280
281 // element size
282 double elemwidth = e.getUpper0() - e.getLower0();
283 double elemheight = e.getUpper1() - e.getLower1();
284
285 Vec<2, double> midpoint = e.getMidpoint();
286
287 double kx = conds[i].c00;
288 double ky = conds[i].c11;
289
290 kx *= elemheight;
291 kx /= elemwidth;
292 ky *= elemwidth;
293 ky /= elemheight;
294
295 // set symmetric matrix components
296 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
297
298 k44 = k33 = k22 = k11 = (kx + ky) / 3.;
299 k43 = k21 = (-2. * kx + ky) / 6.;
300 k42 = k31 = -(kx + ky) / 6.;
301 k32 = k41 = (kx - 2. * ky) / 6.;
302
303 // set stiffness matrix
304 setLocalMatrix(k44, k33, k22, k11, k43, k21, k42, k31, k32, k41, ky, elemwidth, midpoint);
305
306 A(loleftno, loleftno) += k11;
307 A(lorghtno, lorghtno) += k22;
308 A(uprghtno, uprghtno) += k33;
309 A(upleftno, upleftno) += k44;
310
311 A(lorghtno, loleftno) += k21;
312 A(uprghtno, loleftno) += k31;
313 A(upleftno, loleftno) += k41;
314 A(uprghtno, lorghtno) += k32;
315 A(upleftno, lorghtno) += k42;
316 A(upleftno, uprghtno) += k43;
317 }
318
319 A.applyBC(bvoltage, B);
320
321#ifndef NDEBUG
322 double* aend = A.data + A.size;
323 for (double* pa = A.data; pa != aend; ++pa) {
324 if (isnan(*pa) || isinf(*pa))
325 throw ComputationError(this->getId(), "error in stiffness matrix at position {0} ({1})", pa - A.data,
326 isnan(*pa) ? "nan" : "inf");
327 }
328#endif
329}
330
331 template <typename Geometry2DType> LazyData<double> ElectricalFem2DSolver<Geometry2DType>::loadConductivities() {
332 auto midmesh = this->maskedMesh->getElementMesh();
333 auto temperature = inTemperature(midmesh);
334
335 for (auto e : this->maskedMesh->elements()) {
336 size_t i = e.getIndex();
337 Vec<2, double> midpoint = e.getMidpoint();
338
339 auto roles = this->geometry->getRolesAt(midpoint);
340 if (size_t actn = isActive(midpoint)) {
341 const auto& act = active[actn - 1];
342 conds[i] = junction_conductivity[act.offset + e.getIndex0()];
343 if (isnan(conds[i].c11) || abs(conds[i].c11) < 1e-16) conds[i].c11 = 1e-16;
344 } else if (roles.find("p-contact") != roles.end()) {
345 conds[i] = Tensor2<double>(pcond, pcond);
346 } else if (roles.find("n-contact") != roles.end()) {
347 conds[i] = Tensor2<double>(ncond, ncond);
348 } else
349 conds[i] = this->geometry->getMaterial(midpoint)->cond(temperature[i]);
350 }
351
352 return temperature;
353}
354
355 template <typename Geometry2DType> void ElectricalFem2DSolver<Geometry2DType>::saveConductivities() {
356 for (size_t n = 0; n < active.size(); ++n) {
357 const auto& act = active[n];
358 for (size_t i = act.left, r = (act.top + act.bottom) / 2; i != act.right; ++i)
359 junction_conductivity[act.offset + i] = conds[this->maskedMesh->element(i, r).getIndex()];
360 }
361}
362
363 template <typename Geometry2DType> double ElectricalFem2DSolver<Geometry2DType>::compute(unsigned loops) {
364 this->initCalculation();
365
366 heats.reset();
367
368 // Store boundary conditions for current mesh
369 auto vconst = voltage_boundary(this->maskedMesh, this->geometry);
370
371 this->writelog(LOG_INFO, "Running electrical calculations");
372
373 unsigned loop = 0;
374
375 std::unique_ptr<FemMatrix> pA(this->getMatrix());
376 FemMatrix& A = *pA.get();
377
378 double err = 0.;
379 toterr = 0.;
380
381#ifndef NDEBUG
382 if (!potentials.unique()) this->writelog(LOG_DEBUG, "Voltage data held by something else...");
383#endif
384 potentials = potentials.claim();
385
386 DataVector<double> rhs(potentials.size());
387
388 auto temperature = loadConductivities();
389
390 bool noactive = (active.size() == 0);
391 double minj = 100e-7; // assume no significant heating below this current
392
393 do {
394 setMatrix(A, rhs, vconst, temperature);
395 A.solve(rhs, potentials);
396
397 err = 0.;
398 double mcur = 0.;
399 for (auto el : this->maskedMesh->elements()) {
400 size_t i = el.getIndex();
401 size_t loleftno = el.getLoLoIndex();
402 size_t lorghtno = el.getUpLoIndex();
403 size_t upleftno = el.getLoUpIndex();
404 size_t uprghtno = el.getUpUpIndex();
405 double dvx = -0.05 * (-potentials[loleftno] + potentials[lorghtno] - potentials[upleftno] + potentials[uprghtno]) /
406 (el.getUpper0() - el.getLower0()); // [j] = kA/cm²
407 double dvy = -0.05 * (-potentials[loleftno] - potentials[lorghtno] + potentials[upleftno] + potentials[uprghtno]) /
408 (el.getUpper1() - el.getLower1()); // [j] = kA/cm²
409 auto cur = vec(conds[i].c00 * dvx, conds[i].c11 * dvy);
410 if (noactive || isActive(el)) {
411 double acur = abs2(cur);
412 if (acur > mcur) {
413 mcur = acur;
414 maxcur = cur;
415 }
416 }
417 double delta = abs2(currents[i] - cur);
418 if (delta > err) err = delta;
419 currents[i] = cur;
420 }
421 mcur = sqrt(mcur);
422 err = 100. * sqrt(err) / max(mcur, minj);
423 if ((loop != 0 || mcur >= minj) && err > toterr) toterr = err;
424
425 ++loopno;
426 ++loop;
427
428 this->writelog(LOG_RESULT, "Loop {:d}({:d}): max(j{}) = {:g} kA/cm2, error = {:g}%", loop, loopno, noactive ? "" : "@junc",
429 mcur, err);
430
431 } while ((!this->iter_params.converged || err > maxerr) && (loops == 0 || loop < loops));
432
433 saveConductivities();
434
435 outVoltage.fireChanged();
436 outCurrentDensity.fireChanged();
437 outHeat.fireChanged();
438
439 return toterr;
440}
441
442 template <typename Geometry2DType> void ElectricalFem2DSolver<Geometry2DType>::saveHeatDensities() {
443 this->writelog(LOG_DETAIL, "Computing heat densities");
444
445 heats.reset(this->maskedMesh->getElementsCount());
446
447 for (auto e : this->maskedMesh->elements()) {
448 size_t i = e.getIndex();
449 size_t loleftno = e.getLoLoIndex();
450 size_t lorghtno = e.getUpLoIndex();
451 size_t upleftno = e.getLoUpIndex();
452 size_t uprghtno = e.getUpUpIndex();
453 auto midpoint = e.getMidpoint();
454 if (this->geometry->getMaterial(midpoint)->kind() == Material::EMPTY || this->geometry->hasRoleAt("noheat", midpoint))
455 heats[i] = 0.;
456 else {
457 double dvx = 0.5e6 * (-potentials[loleftno] + potentials[lorghtno] - potentials[upleftno] + potentials[uprghtno]) /
458 (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
459 double dvy = 0.5e6 * (-potentials[loleftno] - potentials[lorghtno] + potentials[upleftno] + potentials[uprghtno]) /
460 (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
461 heats[i] = conds[i].c00 * dvx * dvx + conds[i].c11 * dvy * dvy;
462 }
463 }
464}
465
466 template <> double ElectricalFem2DSolver<Geometry2DCartesian>::integrateCurrent(size_t vindex, bool onlyactive) {
467 if (!potentials) throw NoValue("current densities");
468 this->writelog(LOG_DETAIL, "Computing total current");
469 double result = 0.;
470 for (size_t i = 0; i < mesh->axis[0]->size() - 1; ++i) {
471 auto element = maskedMesh->element(i, vindex);
472 if (!onlyactive || isActive(element.getMidpoint())) {
473 size_t index = element.getIndex();
474 if (index != RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX) result += currents[index].c1 * element.getSize0();
475 }
476 }
477 if (this->getGeometry()->isSymmetric(Geometry::DIRECTION_TRAN)) result *= 2.;
478 return result * geometry->getExtrusion()->getLength() * 0.01; // kA/cm² µm² --> mA;
479}
480
481 template <> double ElectricalFem2DSolver<Geometry2DCylindrical>::integrateCurrent(size_t vindex, bool onlyactive) {
482 if (!potentials) throw NoValue("current densities");
483 this->writelog(LOG_DETAIL, "Computing total current");
484 double result = 0.;
485 for (size_t i = 0; i < mesh->axis[0]->size() - 1; ++i) {
486 auto element = maskedMesh->element(i, vindex);
487 if (!onlyactive || isActive(element.getMidpoint())) {
488 size_t index = element.getIndex();
489 if (index != RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX) {
490 double rin = element.getLower0(), rout = element.getUpper0();
491 result += currents[index].c1 * (rout * rout - rin * rin);
492 }
493 }
494 }
495 return result * plask::PI * 0.01; // kA/cm² µm² --> mA
496}
497
498 template <typename Geometry2DType> double ElectricalFem2DSolver<Geometry2DType>::getTotalCurrent(size_t nact) {
499 if (!potentials) throw NoValue("current");
500 if (nact >= active.size()) throw BadInput(this->getId(), "wrong active region number");
501 const auto& act = active[nact];
502 // Find the average of the active region
503 size_t level = (act.bottom + act.top) / 2;
504 return integrateCurrent(level, true);
505}
506
507 template <typename Geometry2DType>
508 const LazyData<double> ElectricalFem2DSolver<Geometry2DType>::getVoltage(shared_ptr<const MeshD<2>> dest_mesh,
509 InterpolationMethod method) const {
510 if (!potentials) throw NoValue("voltage");
511 this->writelog(LOG_DEBUG, "Getting voltage");
512 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
513 if (this->maskedMesh->full())
514 return interpolate(this->mesh, potentials, dest_mesh, method, this->geometry);
515 else
516 return interpolate(this->maskedMesh, potentials, dest_mesh, method, this->geometry);
517}
518
519 template <typename Geometry2DType>
520 const LazyData<Vec<2>> ElectricalFem2DSolver<Geometry2DType>::getCurrentDensities(shared_ptr<const MeshD<2>> dest_mesh,
521 InterpolationMethod method) {
522 if (!potentials) throw NoValue("current density");
523 this->writelog(LOG_DEBUG, "Getting current densities");
524 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
525 InterpolationFlags flags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN);
526 if (this->maskedMesh->full()) {
527 auto result = interpolate(this->mesh->getElementMesh(), currents, dest_mesh, method, flags);
528 return LazyData<Vec<2>>(result.size(), [result, this, flags, dest_mesh](size_t i) {
529 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i))) ? result[i] : Vec<2>(0., 0.);
530 });
531 } else {
532 auto result = interpolate(this->maskedMesh->getElementMesh(), currents, dest_mesh, method, flags);
533 return LazyData<Vec<2>>(result.size(), [result](size_t i) {
534 // Masked mesh always returns NaN outside of itself
535 auto val = result[i];
536 return isnan(val) ? Vec<2>(0., 0.) : val;
537 });
538 }
539}
540
541 template <typename Geometry2DType>
542 const LazyData<double> ElectricalFem2DSolver<Geometry2DType>::getHeatDensities(shared_ptr<const MeshD<2>> dest_mesh,
543 InterpolationMethod method) {
544 if (!potentials) throw NoValue("heat density");
545 this->writelog(LOG_DEBUG, "Getting heat density");
546 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
547 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
548 InterpolationFlags flags(this->geometry);
549 if (this->maskedMesh->full()) {
550 auto result = interpolate(this->mesh->getElementMesh(), heats, dest_mesh, method, flags);
551 return LazyData<double>(result.size(), [result, this, flags, dest_mesh](size_t i) {
552 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i))) ? result[i] : 0.;
553 });
554 } else {
555 auto result = interpolate(this->maskedMesh->getElementMesh(), heats, dest_mesh, method, flags);
556 return LazyData<double>(result.size(), [result](size_t i) {
557 // Masked mesh always returns NaN outside of itself
558 auto val = result[i];
559 return isnan(val) ? 0. : val;
560 });
561 }
562}
563
564 template <typename Geometry2DType>
565 const LazyData<Tensor2<double>> ElectricalFem2DSolver<Geometry2DType>::getConductivity(shared_ptr<const MeshD<2>> dest_mesh,
566 InterpolationMethod) {
567 this->initCalculation();
568 this->writelog(LOG_DEBUG, "Getting conductivities");
569 loadConductivities();
570 InterpolationFlags flags(this->geometry);
571 return interpolate(this->maskedMesh->getElementMesh(), conds, dest_mesh, INTERPOLATION_NEAREST, flags);
572}
573
574 template <> double ElectricalFem2DSolver<Geometry2DCartesian>::getTotalEnergy() {
575 double W = 0.;
576 auto T = inTemperature(this->maskedMesh->getElementMesh());
577 for (auto e : this->maskedMesh->elements()) {
578 size_t ll = e.getLoLoIndex();
579 size_t lu = e.getUpLoIndex();
580 size_t ul = e.getLoUpIndex();
581 size_t uu = e.getUpUpIndex();
582 double dvx = 0.5e6 * (-potentials[ll] + potentials[lu] - potentials[ul] + potentials[uu]) /
583 (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
584 double dvy = 0.5e6 * (-potentials[ll] - potentials[lu] + potentials[ul] + potentials[uu]) /
585 (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
586 double w = this->geometry->getMaterial(e.getMidpoint())->eps(T[e.getIndex()]) * (dvx * dvx + dvy * dvy);
587 double width = e.getUpper0() - e.getLower0();
588 double height = e.getUpper1() - e.getLower1();
589 W += width * height * w;
590 }
591 // TODO add outsides of comptational areas
592 return geometry->getExtrusion()->getLength() * 0.5e-18 * phys::epsilon0 * W; // 1e-18 µm³ -> m³
593}
594
595 template <> double ElectricalFem2DSolver<Geometry2DCylindrical>::getTotalEnergy() {
596 double W = 0.;
597 auto T = inTemperature(this->maskedMesh->getElementMesh());
598 for (auto e : this->maskedMesh->elements()) {
599 size_t ll = e.getLoLoIndex();
600 size_t lu = e.getUpLoIndex();
601 size_t ul = e.getLoUpIndex();
602 size_t uu = e.getUpUpIndex();
603 auto midpoint = e.getMidpoint();
604 double dvx = 0.5e6 * (-potentials[ll] + potentials[lu] - potentials[ul] + potentials[uu]) /
605 (e.getUpper0() - e.getLower0()); // [grad(dV)] = V/m
606 double dvy = 0.5e6 * (-potentials[ll] - potentials[lu] + potentials[ul] + potentials[uu]) /
607 (e.getUpper1() - e.getLower1()); // [grad(dV)] = V/m
608 double w = this->geometry->getMaterial(midpoint)->eps(T[e.getIndex()]) * (dvx * dvx + dvy * dvy);
609 double width = e.getUpper0() - e.getLower0();
610 double height = e.getUpper1() - e.getLower1();
611 W += width * height * midpoint.rad_r() * w;
612 }
613 // TODO add outsides of computational area
614 return 2. * plask::PI * 0.5e-18 * phys::epsilon0 * W; // 1e-18 µm³ -> m³
615}
616
617 template <typename Geometry2DType> double ElectricalFem2DSolver<Geometry2DType>::getCapacitance() {
618 if (this->voltage_boundary.size() != 2) {
619 throw BadInput(this->getId(), "cannot estimate applied voltage (exactly 2 voltage boundary conditions required)");
620 }
621
622 double U = voltage_boundary[0].value - voltage_boundary[1].value;
623
624 return 2e12 * getTotalEnergy() / (U * U); // 1e12 F -> pF
625}
626
627 template <> double ElectricalFem2DSolver<Geometry2DCartesian>::getTotalHeat() {
628 double W = 0.;
629 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
630 for (auto e : this->maskedMesh->elements()) {
631 double width = e.getUpper0() - e.getLower0();
632 double height = e.getUpper1() - e.getLower1();
633 W += width * height * heats[e.getIndex()];
634 }
635 if (this->getGeometry()->isSymmetric(Geometry::DIRECTION_TRAN)) W *= 2.;
636 return geometry->getExtrusion()->getLength() * 1e-15 * W; // 1e-15 µm³ -> m³, W -> mW
637}
638
639 template <> double ElectricalFem2DSolver<Geometry2DCylindrical>::getTotalHeat() {
640 double W = 0.;
641 if (!heats) saveHeatDensities(); // we will compute heats only if they are needed
642 for (auto e : this->maskedMesh->elements()) {
643 double width = e.getUpper0() - e.getLower0();
644 double height = e.getUpper1() - e.getLower1();
645 double r = e.getMidpoint().rad_r();
646 W += width * height * r * heats[e.getIndex()];
647 }
648 return 2e-15 * plask::PI * W; // 1e-15 µm³ -> m³, W -> mW
649}
650
651 template struct PLASK_SOLVER_API ElectricalFem2DSolver<Geometry2DCartesian>;
652 template struct PLASK_SOLVER_API ElectricalFem2DSolver<Geometry2DCylindrical>;
653
654}}} // namespace plask::electrical::shockley