PLaSK library

Loading...
Searching...
No Matches
diffusion3d.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_DIFFUSION3D_H
15#define PLASK__MODULE_ELECTRICAL_DIFFUSION3D_H
16
17#include <plask/common/fem.hpp>
18#include <plask/plask.hpp>
19
20 namespace plask { namespace electrical { namespace diffusion {
21
22 template <typename MeshT> struct PLASK_SOLVER_API QwsLateralMesh3D : MultiLateralMesh3D<MeshT> {
23
24 QwsLateralMesh3D(const shared_ptr<MeshT>& lateral, std::vector<double> vert) :
25 MultiLateralMesh3D<MeshT>(lateral, shared_ptr<OrderedAxis>(new OrderedAxis(vert, 0.))) {}
26
27 QwsLateralMesh3D(const shared_ptr<MeshT>& lateral, const shared_ptr<MeshAxis>& vert) :
28 MultiLateralMesh3D<MeshT>(lateral, vert) {}
29
30 shared_ptr<QwsLateralMesh3D<typename MeshT::ElementMesh>> getElementMesh() const {
31 return shared_ptr<QwsLateralMesh3D<typename MeshT::ElementMesh>>(
32 new QwsLateralMesh3D<typename MeshT::ElementMesh>(this->lateralMesh->getElementMesh(), this->vertAxis));
33 }
34};
35
36 struct ActiveRegion3D {
37 struct Region {
38 size_t bottom, top, lon, tra;
39 bool warn;
40 std::vector<bool> isQW;
41 Region() {}
42 Region(size_t b, size_t t, size_t x, size_t y, const std::vector<bool>& isQW)
43 : bottom(b), top(t), lon(x), tra(y), warn(true), isQW(isQW) {}
44 };
45
46 size_t bottom, top;
47
48 double QWheight;
49
50 shared_ptr<LateralMesh3D<RectangularMaskedMesh2D>> mesh2;
51 shared_ptr<LateralMesh3D<RectangularMaskedMesh2D::ElementMesh>> emesh2;
52 shared_ptr<QwsLateralMesh3D<RectangularMaskedMesh2D>> mesh3;
53 shared_ptr<QwsLateralMesh3D<RectangularMaskedMesh2D::ElementMesh>> emesh3;
54
55 std::vector<std::pair<double, double>> QWs;
56
57 DataVector<double> U;
58
59 std::vector<double> modesP;
60
61 template <typename SolverT>
62 ActiveRegion3D(const SolverT* solver,
63 size_t bottom,
64 size_t top,
65 double h,
66 std::vector<double> QWz,
67 std::vector<std::pair<size_t, size_t>> QWbt)
68 : bottom(bottom), top(top), QWheight(h) {
69 QWs.reserve(QWbt.size());
70 for (auto& bt : QWbt) QWs.emplace_back(solver->getMesh()->vert()->at(bt.first), solver->getMesh()->vert()->at(bt.second));
71
72 double z = QWz[(QWz.size() + 1) / 2 - 1];
73 auto lateral_mesh_unmasked = make_shared<RectangularMesh2D>(solver->getMesh()->lon(), solver->getMesh()->tran());
74 auto lateral_mesh = make_shared<RectangularMaskedMesh2D>(
75 *lateral_mesh_unmasked, [solver, z](const RectangularMesh2D::Element& element) -> bool {
76 auto point = element.getMidpoint();
77 auto roles = solver->getGeometry()->getRolesAt(Vec<3>(point.c0, point.c1, z));
78 return roles.find("QW") != roles.end() || roles.find("QD") != roles.end() || roles.find("carriers") != roles.end();
79 });
80
81 mesh2 = make_shared<LateralMesh3D<RectangularMaskedMesh2D>>(lateral_mesh, z);
82 emesh2 = mesh2->getElementMesh();
83 // mesh3 = make_shared<QwsLateralMesh3D<RectangularMaskedMesh2D>>(lateral_mesh, std::move(QWz));
84 mesh3.reset(new QwsLateralMesh3D<RectangularMaskedMesh2D>(lateral_mesh, std::move(QWz)));
85 emesh3 = mesh3->getElementMesh();
86 }
87
88 double vert() const { return mesh2->vert; }
89
90 template <typename ReceiverType, typename MeshType>
91 LazyData<typename ReceiverType::ValueType> verticallyAverage(
92 const ReceiverType& receiver,
93 const shared_ptr<QwsLateralMesh3D<MeshType>>& mesh,
94 InterpolationMethod interp = InterpolationMethod::INTERPOLATION_DEFAULT) const {
95 auto data = receiver(mesh, interp);
96 const size_t n = mesh->vertAxis->size();
97 return LazyData<typename ReceiverType::ValueType>(
98 mesh->lateralMesh->size(), [this, data, n](size_t i) -> typename ReceiverType::ValueType {
99 typename ReceiverType::ValueType val(Zero<typename ReceiverType::ValueType>());
100 for (size_t j = n * i, end = n * (i + 1); j < end; ++j) val += data[j];
101 return val / double(n);
102 });
103 }
104};
105
106 struct ElementParams3D {
107 const size_t n00, n01, n10, n11, i00, i01, i10, i02, i03, i12, i20, i21, i30, i22, i23, i32;
108 const double X, Y;
109 ElementParams3D(double X, double Y, size_t n00, size_t n01, size_t n10, size_t n11)
110 : n00(n00),
111 n01(n01),
112 n10(n10),
113 n11(n11),
114 i00(3 * n00),
115 i01(i00 + 1),
116 i10(i00 + 2),
117 i02(3 * n01),
118 i03(i02 + 1),
119 i12(i02 + 2),
120 i20(3 * n10),
121 i21(i20 + 1),
122 i30(i20 + 2),
123 i22(3 * n11),
124 i23(i22 + 1),
125 i32(i22 + 2),
126 X(X),
127 Y(Y) {}
128 ElementParams3D(const RectangularMaskedMesh2D::Element& element)
129 : ElementParams3D(element.getSize0(),
130 element.getSize1(),
131 element.getLoLoIndex(),
132 element.getLoUpIndex(),
133 element.getUpLoIndex(),
134 element.getUpUpIndex()) {}
135 ElementParams3D(const ActiveRegion3D& active, size_t ie) : ElementParams3D(active.mesh2->lateralMesh->element(ie)) {}
136};
137
141 struct PLASK_SOLVER_API Diffusion3DSolver : public FemSolverWithMesh<Geometry3D, RectangularMesh<3>> {
142 protected:
144 int loopno;
145 double toterr;
146
147 struct ConcentrationDataImpl : public LazyDataImpl<double> {
148 const Diffusion3DSolver* solver;
149 shared_ptr<const MeshD<3>> destination_mesh;
150 InterpolationFlags interpolationFlags;
151 std::vector<LazyData<double>> concentrations;
152 ConcentrationDataImpl(const Diffusion3DSolver* solver, shared_ptr<const MeshD<3>> dest_mesh, InterpolationMethod interp);
153 double at(size_t i) const override;
154 size_t size() const override { return destination_mesh->size(); }
155 };
156
157 std::map<size_t, ActiveRegion3D> active;
158
160 inline void setLocalMatrix(FemMatrix& K,
161 DataVector<double>& F,
162 const ElementParams3D e,
163 const double A,
164 const double B,
165 const double C,
166 const double D,
167 const double* U,
168 const double* J);
169
171 inline void addLocalBurningMatrix(FemMatrix& K,
172 DataVector<double>& F,
173 const ElementParams3D e,
174 const Tensor2<double> G,
175 const Tensor2<double> dG,
176 const double Ug,
177 const Tensor2<double>* P);
178
180 template <typename T> inline T integrateBilinear(const double Lx, const double Ly, const T* P) {
181 return 0.25 * (P[0] + P[1] + P[2] + P[3]) * Lx * Lx;
182 }
183
185 void onInitialize() override;
186
188 void onInvalidate() override;
189
191 void setupActiveRegions();
192
193 // void computeInitial(ActiveRegion3D& active);
194
195 private:
196 void summarizeActiveRegion(std::map<size_t, ActiveRegion3D::Region>& regions,
197 size_t num,
198 size_t start,
199 size_t ver,
200 size_t lon,
201 size_t tra,
202 const std::vector<bool>& isQW,
203 const shared_ptr<RectangularMesh3D::ElementMesh>& points);
204
205 public:
206 double maxerr;
207
209 BoundaryConditions<RectangularMesh<2>::Boundary, double> voltage_boundary;
210
211 ReceiverFor<CurrentDensity, Geometry3D> inCurrentDensity;
212
213 ReceiverFor<Temperature, Geometry3D> inTemperature;
214
215 ReceiverFor<Gain, Geometry3D> inGain;
216
217 ReceiverFor<ModeWavelength> inWavelength;
218
219 ReceiverFor<ModeLightE, Geometry3D> inLightE;
220
221 typename ProviderFor<CarriersConcentration, Geometry3D>::Delegate outCarriersConcentration;
222
223 std::string getClassName() const override { return "electrical.Diffusion3D"; }
224
232 double compute(unsigned loops, bool shb, size_t act);
233
240 double compute(unsigned loops = 0, bool shb = false) {
241 this->initCalculation();
242 double maxerr = 0.;
243 for (const auto& act : active) maxerr = max(maxerr, compute(loops, shb, act.first));
244 return maxerr;
245 }
246
247 void loadConfiguration(XMLReader& source, Manager& manager) override;
248
249 void parseConfiguration(XMLReader& source, Manager& manager);
250
251 Diffusion3DSolver(const std::string& name = "");
252
253 ~Diffusion3DSolver();
254
255 size_t activeRegionsCount() const { return active.size(); }
256
257 double get_burning_integral_for_mode(size_t mode) const;
258
259 double get_burning_integral() const;
260
261 protected:
262 const LazyData<double> getConcentration(CarriersConcentration::EnumType what,
263 shared_ptr<const MeshD<3>> dest_mesh,
264 InterpolationMethod interpolation = INTERPOLATION_DEFAULT) const;
265};
266
267}}} // namespace plask::electrical::diffusion
268
269#endif