PLaSK library

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