PLaSK library

Loading...
Searching...
No Matches
diffusion2d.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_DIFFUSION2D_H
15#define PLASK__MODULE_ELECTRICAL_DIFFUSION2D_H
16
17#include <plask/common/fem.hpp>
18#include <plask/plask.hpp>
19
20 namespace plask { namespace electrical { namespace diffusion {
21
22 struct ActiveRegion2D {
23 struct Region {
24 size_t left, right, bottom, top;
25 size_t rowl, rowr;
26 bool warn;
27 Region()
28 : left(0),
29 right(0),
30 bottom(std::numeric_limits<size_t>::max()),
31 top(std::numeric_limits<size_t>::max()),
32 rowl(std::numeric_limits<size_t>::max()),
33 rowr(0),
34 warn(true) {}
35 };
36
37 size_t left, right, bottom, top;
38 double QWheight;
39
40 shared_ptr<RectangularMesh<2>> mesh2, emesh2, mesh1, emesh1;
41 std::vector<std::pair<double, double>> QWs;
42
43 DataVector<double> U;
44
45 std::vector<double> modesP;
46
47 template <typename SolverT>
48 ActiveRegion2D(const SolverT* solver,
49 size_t l,
50 size_t r,
51 size_t b,
52 size_t t,
53 double h,
54 std::vector<double> QWz,
55 std::vector<std::pair<size_t, size_t>> QWbt)
56 : left(l),
57 right(r),
58 bottom(b),
59 top(t),
60 QWheight(h),
61 mesh2(new RectangularMesh<2>(shared_ptr<OrderedAxis>(new OrderedAxis()),
62 shared_ptr<OrderedAxis>(new OrderedAxis(QWz)),
63 RectangularMesh<2>::ORDER_01)) {
64 auto begin = solver->getMesh()->tran()->begin();
65 dynamic_pointer_cast<OrderedAxis>(mesh2->tran())->addOrderedPoints(begin + left, begin + right + 1);
66
67 QWs.reserve(QWbt.size());
68 for (auto& bt : QWbt) QWs.emplace_back(solver->getMesh()->vert()->at(bt.first), solver->getMesh()->vert()->at(bt.second));
69
70 shared_ptr<OrderedAxis> mesh = this->mesh();
71 emesh2.reset(new RectangularMesh<2>(mesh->getMidpointAxis(), mesh2->vert(), RectangularMesh<2>::ORDER_01));
72 mesh1.reset(new RectangularMesh<2>(mesh, make_shared<OnePointAxis>(this->vert())));
73 emesh1.reset(new RectangularMesh<2>(emesh2->tran(), mesh1->vert()));
74
75 // OrderedAxis faxis0;
76 // faxis0.addOrderedPoints(mesh->begin(), mesh->end(), mesh->size(), 0.);
77 // faxis0.addOrderedPoints(emesh1->tran()->begin(), emesh1->tran()->end(), emesh1->tran()->size(), 0.);
78 // fmesh1.reset(new RectangularMesh<2>(faxis0.getMidpointAxis(), mesh1->vert(), RectangularMesh<2>::ORDER_01));
79 // fmesh1.reset(new RectangularMesh<2>(fmesh1->tran(), emesh1->vert()));
80 // assert(fmesh1->size() == 2 * emesh1->size());
81 }
82
83 shared_ptr<OrderedAxis> mesh() const {
84 assert(dynamic_pointer_cast<OrderedAxis>(mesh2->tran()));
85 return static_pointer_cast<OrderedAxis>(mesh2->tran());
86 }
87
88 double vert() const { return mesh2->vert()->at((mesh2->vert()->size() + 1) / 2 - 1); }
89
90 template <typename ReceiverType>
91 LazyData<typename ReceiverType::ValueType> verticallyAverage(
92 const ReceiverType& receiver,
93 const shared_ptr<const RectangularMesh<2>>& mesh,
94 InterpolationMethod interp = InterpolationMethod::INTERPOLATION_DEFAULT) const {
95 assert(mesh->getIterationOrder() == RectangularMesh<2>::ORDER_01);
96 auto data = receiver(mesh, interp);
97 const size_t n = mesh->vert()->size();
98 return LazyData<typename ReceiverType::ValueType>(
99 mesh->tran()->size(), [this, data, n](size_t i) -> typename ReceiverType::ValueType {
100 typename ReceiverType::ValueType val(Zero<typename ReceiverType::ValueType>());
101 for (size_t j = n * i, end = n * (i + 1); j < end; ++j) val += data[j];
102 return val / double(n);
103 });
104 }
105};
106
110 template <typename Geometry2DType>
111 struct PLASK_SOLVER_API Diffusion2DSolver : public FemSolverWithMesh<Geometry2DType, RectangularMesh<2>> {
112 protected:
114 int loopno;
115 double toterr;
116
117 struct ConcentrationDataImpl : public LazyDataImpl<double> {
118 const Diffusion2DSolver* solver;
119 shared_ptr<const MeshD<2>> destination_mesh;
120 InterpolationFlags interpolationFlags;
121 std::vector<LazyData<double>> concentrations;
122 ConcentrationDataImpl(const Diffusion2DSolver* solver, shared_ptr<const MeshD<2>> dest_mesh, InterpolationMethod interp);
123 double at(size_t i) const override;
124 size_t size() const override { return destination_mesh->size(); }
125 };
126
127 std::map<size_t, ActiveRegion2D> active;
128
130 inline void setLocalMatrix(const double R,
131 const double L,
132 const double L2,
133 const double L3,
134 const double L4,
135 const double L5,
136 const double L6,
137 const double A,
138 const double B,
139 const double C,
140 const double D,
141 const double* U,
142 const double* J,
143 double& K00,
144 double& K01,
145 double& K02,
146 double& K03,
147 double& K11,
148 double& K12,
149 double& K13,
150 double& K22,
151 double& K23,
152 double& K33,
153 double& F0,
154 double& F1,
155 double& F2,
156 double& F3);
157
159 inline void addLocalBurningMatrix(const double R,
160 const double L,
161 const double L2,
162 const double L3,
163 const double* P,
164 const double* g,
165 const double* dg,
166 const double ug,
167 double& K00,
168 double& K01,
169 double& K02,
170 double& K03,
171 double& K11,
172 double& K12,
173 double& K13,
174 double& K22,
175 double& K23,
176 double& K33,
177 double& F0,
178 double& F1,
179 double& F2,
180 double& F3);
181
183 template <typename T> inline T integrateLinear(const double R, const double L, const T* P);
184
186 void onInitialize() override;
187
189 void onInvalidate() override;
190
192 void setupActiveRegion2Ds();
193
194 // void computeInitial(ActiveRegion2D& active);
195
200 size_t isActive(const Vec<2>& point) const {
201 size_t no(0);
202 auto roles = this->geometry->getRolesAt(point);
203 for (auto role : roles) {
204 size_t l = 0;
205 if (role.substr(0, 6) == "active")
206 l = 6;
207 else if (role.substr(0, 8) == "junction")
208 l = 8;
209 else
210 continue;
211 if (no != 0) throw BadInput(this->getId(), "multiple 'active'/'junction' roles specified");
212 if (role.size() == l)
213 no = 1;
214 else {
215 try {
216 no = boost::lexical_cast<size_t>(role.substr(l)) + 1;
217 } catch (boost::bad_lexical_cast&) {
218 throw BadInput(this->getId(), "bad junction number in role '{0}'", role);
219 }
220 }
221 }
222 return no;
223 }
224
226 size_t isActive(const RectangularMesh2D::Element& element) const { return isActive(element.getMidpoint()); }
227
228 struct PLASK_SOLVER_API From1DGenerator : public MeshGeneratorD<2> {
229 shared_ptr<MeshGeneratorD<1>> generator;
230
231 From1DGenerator(const shared_ptr<MeshGeneratorD<1>>& generator) : generator(generator) {}
232
233 shared_ptr<MeshD<2>> generate(const shared_ptr<GeometryObjectD<2>>& geometry) override {
234 auto simple_mesh = makeGeometryGrid(geometry);
235 auto mesh1d = (*generator)(geometry);
236 if (shared_ptr<MeshAxis> axis = dynamic_pointer_cast<MeshAxis>(mesh1d))
237 return make_shared<RectangularMesh<2>>(axis, simple_mesh->vert());
238 throw BadInput("generator1D", "1D mesh must be MeshAxis");
239 }
240 };
241
242 public:
243 double maxerr;
244
246 BoundaryConditions<RectangularMesh<2>::Boundary, double> voltage_boundary;
247
248 ReceiverFor<CurrentDensity, Geometry2DType> inCurrentDensity;
249
250 ReceiverFor<Temperature, Geometry2DType> inTemperature;
251
252 ReceiverFor<Gain, Geometry2DType> inGain;
253
254 ReceiverFor<ModeWavelength> inWavelength;
255
256 ReceiverFor<ModeLightE, Geometry2DType> inLightE;
257
258 typename ProviderFor<CarriersConcentration, Geometry2DType>::Delegate outCarriersConcentration;
259
260 std::string getClassName() const override;
261
269 double compute(unsigned loops, bool shb, size_t act);
270
277 double compute(unsigned loops = 0, bool shb = false) {
278 this->initCalculation();
279 double maxerr = 0.;
280 for (const auto& act: active) maxerr = max(maxerr, compute(loops, shb, act.first));
281 return maxerr;
282 }
283
284 void loadConfiguration(XMLReader& source, Manager& manager) override;
285
286 void parseConfiguration(XMLReader& source, Manager& manager);
287
288 Diffusion2DSolver(const std::string& name = "");
289
290 ~Diffusion2DSolver();
291
292 using SolverWithMesh<Geometry2DType, RectangularMesh<2>>::setMesh;
293
294 void setMesh(shared_ptr<MeshD<1>> mesh) {
295 auto simple_mesh = makeGeometryGrid(this->geometry);
296 if (shared_ptr<MeshAxis> axis = dynamic_pointer_cast<MeshAxis>(mesh)) {
297 shared_ptr<RectangularMesh<2>> mesh2d(new RectangularMesh<2>(axis, simple_mesh->vert()));
298 this->setMesh(mesh2d);
299 } else
300 throw BadInput(this->getId(), "1D mesh must be MeshAxis");
301 }
302
303 void setMesh(shared_ptr<MeshGeneratorD<1>> generator) { this->setMesh(make_shared<From1DGenerator>(generator)); }
304
305 size_t activeRegionsCount() const { return active.size(); }
306
307 double get_burning_integral_for_mode(size_t mode) const;
308
309 double get_burning_integral() const;
310
311 protected:
312 const LazyData<double> getConcentration(CarriersConcentration::EnumType what,
313 shared_ptr<const MeshD<2>> dest_mesh,
314 InterpolationMethod interpolation = INTERPOLATION_DEFAULT) const;
315};
316
317}}} // namespace plask::electrical::diffusion
318
319#endif