PLaSK library

Loading...
Searching...
No Matches
therm3d.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 <type_traits>
15
16#include "therm3d.hpp"
17
18 namespace plask { namespace thermal { namespace tstatic {
19
20 ThermalFem3DSolver::ThermalFem3DSolver(const std::string& name) :
21 FemSolverWithMaskedMesh<Geometry3D, RectangularMesh<3>>(name),
22 loopno(0),
23 inittemp(300.),
24 maxerr(0.05),
25 outTemperature(this, &ThermalFem3DSolver::getTemperatures),
26 outHeatFlux(this, &ThermalFem3DSolver::getHeatFluxes),
27 outThermalConductivity(this, &ThermalFem3DSolver::getThermalConductivity)
28{
29 temperatures.reset();
30 fluxes.reset();
31 inHeat = 0.;
32 algorithm = ALGORITHM_ITERATIVE;
33}
34
35
36 ThermalFem3DSolver::~ThermalFem3DSolver() {
37}
38
39
40 void ThermalFem3DSolver::loadConfiguration(XMLReader &source, Manager &manager)
41{
42 while (source.requireTagOrEnd())
43 {
44 std::string param = source.getNodeName();
45
46 if (param == "temperature")
47 this->readBoundaryConditions(manager, source, temperature_boundary);
48
49 else if (param == "heatflux")
50 this->readBoundaryConditions(manager, source, heatflux_boundary);
51
52 else if (param == "convection")
53 this->readBoundaryConditions(manager, source, convection_boundary);
54
55 else if (param == "radiation")
56 this->readBoundaryConditions(manager, source, radiation_boundary);
57
58 else if (param == "loop") {
59 inittemp = source.getAttribute<double>("inittemp", inittemp);
60 maxerr = source.getAttribute<double>("maxerr", maxerr);
61 source.requireTagEnd();
62 }
63
64 else if (!parseFemConfiguration(source, manager)) {
65 this->parseStandardConfiguration(source, manager);
66 }
67 }
68}
69
70
71 void ThermalFem3DSolver::onInitialize() {
72 if (!this->geometry) throw NoGeometryException(this->getId());
73 if (!this->mesh) throw NoMeshException(this->getId());
74
75 FemSolverWithMaskedMesh<Geometry3D, RectangularMesh<3>>::onInitialize();
76
77 loopno = 0;
78
79 temperatures.reset(this->maskedMesh->size(), inittemp);
80
81 thickness.reset(this->maskedMesh->getElementsCount(), NAN);
82 // Set stiffness matrix and load vector
83 for (auto elem: this->maskedMesh->elements())
84 {
85 if (!isnan(thickness[elem.getIndex()])) continue;
86 auto material = this->geometry->getMaterial(elem.getMidpoint());
87 double top = elem.getUpper2(), bottom = elem.getLower2();
88 size_t row = elem.getIndex2();
89 size_t itop = row+1, ibottom = row;
90 for (size_t r = row; r > 0; r--) {
91 auto e = this->mesh->element(elem.getIndex0(), elem.getIndex1(), r-1);
92 auto m = this->geometry->getMaterial(e.getMidpoint());
93 if (m == material) { //TODO ignore doping
94 bottom = e.getLower2();
95 ibottom = r-1;
96 }
97 else break;
98 }
99 for (size_t r = elem.getIndex2()+1; r < this->mesh->axis[2]->size()-1; r++) {
100 auto e = this->mesh->element(elem.getIndex0(), elem.getIndex1(), r);
101 auto m = this->geometry->getMaterial(e.getMidpoint());
102 if (m == material) { //TODO ignore doping
103 top = e.getUpper2();
104 itop = r+1;
105 }
106 else break;
107 }
108 double h = top - bottom;
109 for (size_t r = ibottom; r < itop; ++r) {
110 size_t idx = this->maskedMesh->element(elem.getIndex0(), elem.getIndex1(), r).getIndex();
111 if (idx != RectangularMaskedMesh3D::Element::UNKNOWN_ELEMENT_INDEX)
112 thickness[idx] = h;
113 }
114 }
115}
116
117
118 void ThermalFem3DSolver::onInvalidate() {
119 temperatures.reset();
120 fluxes.reset();
121 thickness.reset();
122}
123
124 void ThermalFem3DSolver::setAlgorithm(FemMatrixAlgorithm alg) {
125 //TODO
126 algorithm = alg;
127}
128
140 template <typename ConditionT>
141 static void setBoundaries(const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary,ConditionT>& boundary_conditions,
142 const size_t (&idx)[8], double dx, double dy, double dz, double (&F)[8], double (&K)[8][8],
143 const std::function<double(double,ConditionT,size_t)>& F_function,
144 const std::function<double(double,ConditionT,ConditionT,size_t,size_t,bool)>& K_function
145 )
146{
147 plask::optional<ConditionT> values[8];
148 for (int i = 0; i < 8; ++i) values[i] = boundary_conditions.getValue(idx[i]);
149
150 constexpr int walls[6][4] = { {0,1,2,3}, {4,5,6,7}, {0,2,4,6}, {1,3,5,7}, {0,1,4,5}, {2,3,6,7} };
151 const double areas[3] = { dx*dy, dy*dz, dz*dx };
152
153 for (int side = 0; side < 6; ++side) {
154 const auto& wall = walls[side];
155 if (values[wall[0]] && values[wall[1]] && values[wall[2]] && values[wall[3]]) {
156 double area = areas[side/2];
157 for (int i = 0; i < 4; ++i) {
158 F[i] += F_function(area, *values[wall[i]], wall[i]);
159 for (int j = 0; j <= i; ++j) {
160 int ij = i ^ j; // numbers on the single edge differ by one bit only, so detect different bits
161 bool edge = (ij == 1 || ij == 2 || ij == 4);
162 K[i][j] += K_function(area, *values[wall[i]], *values[wall[j]], wall[i], wall[j], edge);
163 }
164 }
165 }
166
167 }
168}
169
170 void ThermalFem3DSolver::setMatrix(FemMatrix& A, DataVector<double>& B,
171 const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary,double>& btemperature,
172 const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary,double>& bheatflux,
173 const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary,Convection>& bconvection,
174 const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary,Radiation>& bradiation
175 )
176{
177 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
178
179 auto heats = inHeat(maskedMesh->getElementMesh()/*, INTERPOLATION_NEAREST*/);
180
181 // zero the matrix and the load vector
182 A.clear();
183 B.fill(0.);
184
185 // Set stiffness matrix and load vector
186 for (auto elem: maskedMesh->elements())
187 {
188 // nodes numbers for the current element
189 size_t idx[8];
190 idx[0] = elem.getLoLoLoIndex(); // z y 6-----7
191 idx[1] = elem.getUpLoLoIndex(); // |/__x /| /|
192 idx[2] = elem.getLoUpLoIndex(); // 4-----5 |
193 idx[3] = elem.getUpUpLoIndex(); // | 2---|-3
194 idx[4] = elem.getLoLoUpIndex(); // |/ |/
195 idx[5] = elem.getUpLoUpIndex(); // 0-----1
196 idx[6] = elem.getLoUpUpIndex(); //
197 idx[7] = elem.getUpUpUpIndex(); //
198
199 // element size
200 double dx = elem.getUpper0() - elem.getLower0();
201 double dy = elem.getUpper1() - elem.getLower1();
202 double dz = elem.getUpper2() - elem.getLower2();
203
204 // point and material in the middle of the element
205 Vec<3> middle = elem.getMidpoint();
206 auto material = geometry->getMaterial(middle);
207
208 // average temperature on the element
209 double temp = 0.; for (int i = 0; i < 8; ++i) temp += temperatures[idx[i]]; temp *= 0.125;
210
211 // thermal conductivity
212 double kx, ky, kz;
213 std::tie(ky,kz) = std::tuple<double,double>(material->thermk(temp, thickness[elem.getIndex()]));
214
215 ky *= 1e-6; kz *= 1e-6; // W/m -> W/µm
216 kx = ky;
217
218 kx /= dx; kx *= dy; kx *= dz;
219 ky *= dx; ky /= dy; ky *= dz;
220 kz *= dx; kz *= dy; kz /= dz;
221
222 // load vector: heat densities
223 double f = 0.125e-18 * dx * dy * dz * heats[elem.getIndex()]; // 1e-18 -> to transform µm³ into m³
224
225 // set symmetric matrix components
226 double K[8][8];
227 K[0][0] = K[1][1] = K[2][2] = K[3][3] = K[4][4] = K[5][5] = K[6][6] = K[7][7] = (kx + ky + kz) / 9.;
228
229 K[1][0] = K[3][2] = K[5][4] = K[7][6] = (-2.*kx + ky + kz) / 18.;
230 K[2][0] = K[3][1] = K[6][4] = K[7][5] = ( kx - 2.*ky + kz) / 18.;
231 K[4][0] = K[5][1] = K[6][2] = K[7][3] = ( kx + ky - 2.*kz) / 18.;
232
233 K[4][2] = K[5][3] = K[6][0] = K[7][1] = ( kx - 2.*ky - 2.*kz) / 36.;
234 K[4][1] = K[5][0] = K[6][3] = K[7][2] = (-2.*kx + ky - 2.*kz) / 36.;
235 K[2][1] = K[3][0] = K[6][5] = K[7][4] = (-2.*kx - 2.*ky + kz) / 36.;
236
237 K[4][3] = K[5][2] = K[6][1] = K[7][0] = -(kx + ky + kz) / 36.;
238
239 double F[8];
240 std::fill_n(F, 8, f);
241
242 // boundary conditions: heat flux
243 setBoundaries<double>(bheatflux, idx, dx, dy, dz, F, K,
244 [](double area, double value, size_t) { // F
245 return - 0.25e-12 * area * value;
246 },
247 [](double, double, double, size_t, size_t, bool) { return 0.; } // K
248 );
249
250 // boundary conditions: convection
251 setBoundaries<Convection>(bconvection, idx, dx, dy, dz, F, K,
252 [](double area, Convection value, size_t) { // F
253 return 0.25e-12 * area * value.coeff * value.ambient;
254 },
255 [](double area, Convection value1, Convection value2, size_t i1, size_t i2, bool edge) -> double { // K
256 double v = 0.125e-12 * area * (value1.coeff + value2.coeff);
257 return v / (i2==i1? 9. : edge? 18. : 36.);
258 }
259 );
260
261 // boundary conditions: radiation
262 setBoundaries<Radiation>(bradiation, idx, dx, dy, dz, F, K,
263 [this](double area, Radiation value, size_t i) -> double { // F
264 double a = value.ambient; a = a*a;
265 double T = this->temperatures[i]; T = T*T;
266 return - 0.25e-12 * area * value.emissivity * phys::SB * (T*T - a*a);},
267 [](double, Radiation, Radiation, size_t, size_t, bool) {return 0.;} // K
268 );
269
270 for (int i = 0; i < 8; ++i) {
271 for (int j = 0; j <= i; ++j) {
272 A(idx[i],idx[j]) += K[i][j];
273 }
274 B[idx[i]] += F[i];
275 }
276 }
277
278 A.applyBC(btemperature, B);
279}
280
281 double ThermalFem3DSolver::compute(int loops)
282{
283 this->initCalculation();
284
285 fluxes.reset();
286
287 // store boundary conditions for current mesh
288 auto btemperature = temperature_boundary(this->maskedMesh, this->geometry);
289 auto bheatflux = heatflux_boundary(this->maskedMesh, this->geometry);
290 auto bconvection = convection_boundary(this->maskedMesh, this->geometry);
291 auto bradiation = radiation_boundary(this->maskedMesh, this->geometry);
292
293 this->writelog(LOG_INFO, "Running thermal calculations");
294
295 int loop = 0;
296 size_t size = maskedMesh->size();
297
298 std::unique_ptr<FemMatrix> pA(getMatrix());
299 FemMatrix& A = *pA.get();
300
301 double err = 0.;
302 toterr = 0.;
303
304# ifndef NDEBUG
305 if (!temperatures.unique()) this->writelog(LOG_DEBUG, "Temperature data held by something else...");
306# endif
307 temperatures = temperatures.claim();
308
309 DataVector<double> BT(size), temp0(size);
310
311 do {
312 std::copy(temperatures.begin(), temperatures.end(), temp0.begin());
313
314 setMatrix(A, BT, btemperature, bheatflux, bconvection, bradiation);
315 A.solve(BT, temperatures);
316
317 // Update error
318 err = 0.;
319 maxT = 0.;
320 for (auto temp = temperatures.begin(), t = temp0.begin(); t != temp0.end(); ++temp, ++t)
321 {
322 double corr = std::abs(*t - *temp); // for boundary with constant temperature this will be zero anyway
323 if (corr > err) err = corr;
324 if (*temp > maxT) maxT = *temp;
325 }
326 if (err > toterr) toterr = err;
327
328 ++loopno;
329 ++loop;
330
331 // show max correction
332 this->writelog(LOG_RESULT, "Loop {:d}({:d}): max(T) = {:.3f} K, error = {:g} K", loop, loopno, maxT, err);
333
334 } while ((!iter_params.converged || err > maxerr) && (loops == 0 || loop < loops));
335
336 outTemperature.fireChanged();
337 outHeatFlux.fireChanged();
338
339 return toterr;
340}
341
342 void ThermalFem3DSolver::saveHeatFluxes()
343{
344 this->writelog(LOG_DETAIL, "Computing heat fluxes");
345
346 fluxes.reset(this->maskedMesh->getElementsCount());
347
348 for (auto el: this->maskedMesh->elements())
349 {
350 Vec<3,double> midpoint = el.getMidpoint();
351 auto material = this->geometry->getMaterial(midpoint);
352
353 size_t lll = el.getLoLoLoIndex();
354 size_t llu = el.getLoLoUpIndex();
355 size_t lul = el.getLoUpLoIndex();
356 size_t luu = el.getLoUpUpIndex();
357 size_t ull = el.getUpLoLoIndex();
358 size_t ulu = el.getUpLoUpIndex();
359 size_t uul = el.getUpUpLoIndex();
360 size_t uuu = el.getUpUpUpIndex();
361
362 double temp = 0.125 * (temperatures[lll] + temperatures[llu] + temperatures[lul] + temperatures[luu] +
363 temperatures[ull] + temperatures[ulu] + temperatures[uul] + temperatures[uuu]);
364
365 double kxy, kz;
366 auto leaf = dynamic_pointer_cast<const GeometryObjectD<3>>(geometry->getMatchingAt(midpoint, &GeometryObject::PredicateIsLeaf));
367 if (leaf)
368 std::tie(kxy,kz) = std::tuple<double,double>(material->thermk(temp, leaf->getBoundingBox().height()));
369 else
370 std::tie(kxy,kz) = std::tuple<double,double>(material->thermk(temp));
371
372 fluxes[el.getIndex()] = vec(
373 - 0.25e6 * kxy * (- temperatures[lll] - temperatures[llu] - temperatures[lul] - temperatures[luu]
374 + temperatures[ull] + temperatures[ulu] + temperatures[uul] + temperatures[uuu])
375 / (el.getUpper0() - el.getLower0()), // 1e6 - from µm to m
376 - 0.25e6 * kxy * (- temperatures[lll] - temperatures[llu] + temperatures[lul] + temperatures[luu]
377 - temperatures[ull] - temperatures[ulu] + temperatures[uul] + temperatures[uuu])
378 / (el.getUpper1() - el.getLower1()), // 1e6 - from µm to m
379 - 0.25e6 * kz * (- temperatures[lll] + temperatures[llu] - temperatures[lul] + temperatures[luu]
380 - temperatures[ull] + temperatures[ulu] - temperatures[uul] + temperatures[uuu])
381 / (el.getUpper2() - el.getLower2()) // 1e6 - from µm to m
382 );
383 }
384}
385
386
387 const LazyData<double> ThermalFem3DSolver::getTemperatures(const shared_ptr<const MeshD<3>>& dst_mesh, InterpolationMethod method) const {
388 this->writelog(LOG_DEBUG, "Getting temperatures");
389 if (!temperatures) return LazyData<double>(dst_mesh->size(), inittemp); // in case the receiver is connected and no temperature calculated yet
390 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
391 if (this->maskedMesh->full())
392 return SafeData<double>(interpolate(this->mesh, temperatures, dst_mesh, method, this->geometry), 300.);
393 else
394 return SafeData<double>(interpolate(this->maskedMesh, temperatures, dst_mesh, method, this->geometry), 300.);
395}
396
397
398 const LazyData<Vec<3>> ThermalFem3DSolver::getHeatFluxes(const shared_ptr<const MeshD<3>>& dst_mesh, InterpolationMethod method) {
399 this->writelog(LOG_DEBUG, "Getting heat fluxes");
400 if (!temperatures) return LazyData<Vec<3>>(dst_mesh->size(), Vec<3>(0.,0.,0.)); // in case the receiver is connected and no fluxes calculated yet
401 if (!fluxes) saveHeatFluxes(); // we will compute fluxes only if they are needed
402 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
403 if (this->maskedMesh->full())
404 return SafeData<Vec<3>>(interpolate(this->mesh->getElementMesh(), fluxes, dst_mesh, method,
405 InterpolationFlags(this->geometry, InterpolationFlags::Symmetry::NPP, InterpolationFlags::Symmetry::PNP, InterpolationFlags::Symmetry::PPN)),
406 Zero<Vec<3>>());
407 else
408 return SafeData<Vec<3>>(interpolate(this->maskedMesh->getElementMesh(), fluxes, dst_mesh, method,
409 InterpolationFlags(this->geometry,InterpolationFlags::Symmetry::NPP, InterpolationFlags::Symmetry::PNP, InterpolationFlags::Symmetry::PPN)),
410 Zero<Vec<3>>());
411}
412
413
414 ThermalFem3DSolver::
415 ThermalConductivityData::ThermalConductivityData(const ThermalFem3DSolver* solver, const shared_ptr<const MeshD<3>>& dst_mesh):
416 solver(solver), dest_mesh(dst_mesh), flags(solver->geometry)
417{
418 if (solver->temperatures) temps = interpolate(solver->maskedMesh, solver->temperatures, solver->maskedMesh->getElementMesh(), INTERPOLATION_LINEAR);
419 else temps = LazyData<double>(solver->maskedMesh->getElementsCount(), solver->inittemp);
420}
421 Tensor2<double> ThermalFem3DSolver::ThermalConductivityData::at(std::size_t i) const {
422 auto point = flags.wrap(dest_mesh->at(i));
423 std::size_t x = solver->mesh->axis[0]->findUpIndex(point[0]),
424 y = solver->mesh->axis[1]->findUpIndex(point[1]),
425 z = solver->mesh->axis[2]->findUpIndex(point[2]);
426 if (x == 0 || y == 0 || z == 0 || x == solver->mesh->axis[0]->size() || y == solver->mesh->axis[1]->size() || z == solver->mesh->axis[2]->size())
427 return Tensor2<double>(NAN);
428 else {
429 auto elem = solver->maskedMesh->element(x-1, y-1, z-1);
430 size_t idx = elem.getIndex();
431 if (idx == RectangularMaskedMesh3D::Element::UNKNOWN_ELEMENT_INDEX) return Tensor2<double>(NAN);
432 auto material = solver->geometry->getMaterial(elem.getMidpoint());
433 return material->thermk(temps[idx], solver->thickness[idx]);
434 }
435}
436std::size_t ThermalFem3DSolver::ThermalConductivityData::size() const { return dest_mesh->size(); }
437
438 const LazyData<Tensor2<double>> ThermalFem3DSolver::getThermalConductivity(const shared_ptr<const MeshD<3>>& dst_mesh, InterpolationMethod /*method*/) {
439 this->initCalculation();
440 this->writelog(LOG_DEBUG, "Getting thermal conductivities");
441 return LazyData<Tensor2<double>>(new ThermalFem3DSolver::ThermalConductivityData(this, dst_mesh));
442}
443
444
445}}} // namespace plask::thermal::tstatic