PLaSK library
Loading...
Searching...
No Matches
fem_solver.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) 2023 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_COMMON_FEM_FEM_SOLVER_HPP
15#define PLASK_COMMON_FEM_FEM_SOLVER_HPP
16
17#include <plask/plask.hpp>
18#include "cholesky_matrix.hpp"
19#include "gauss_matrix.hpp"
20#include "iterative_matrix.hpp"
21
22namespace plask {
23
30
31template <typename SpaceT, typename MeshT> struct FemSolverWithMesh : public SolverWithMesh<SpaceT, MeshT> {
32 FemMatrixAlgorithm algorithm = ALGORITHM_CHOLESKY;
33
34 IterativeMatrixParams iter_params;
35
36 FemSolverWithMesh(const std::string& name = "") : SolverWithMesh<SpaceT, MeshT>(name) {}
37
38 bool parseFemConfiguration(XMLReader& reader, Manager& manager) {
39 if (reader.getNodeName() == "matrix") {
40 algorithm = reader.enumAttribute<FemMatrixAlgorithm>("algorithm")
41 .value("cholesky", ALGORITHM_CHOLESKY)
42 .value("gauss", ALGORITHM_GAUSS)
43 .value("iterative", ALGORITHM_ITERATIVE)
44 .get(algorithm);
45
46 if (reader.requireTagOrEnd("iterative")) {
47 iter_params.accelerator = reader.enumAttribute<IterativeMatrixParams::Accelelator>("accelerator")
48 .value("cg", IterativeMatrixParams::ACCEL_CG)
49 .value("si", IterativeMatrixParams::ACCEL_SI)
50 .value("sor", IterativeMatrixParams::ACCEL_SOR)
51 .value("srcg", IterativeMatrixParams::ACCEL_SRCG)
52 .value("srsi", IterativeMatrixParams::ACCEL_SRSI)
53 .value("basic", IterativeMatrixParams::ACCEL_BASIC)
54 .value("me", IterativeMatrixParams::ACCEL_ME)
55 .value("cgnr", IterativeMatrixParams::ACCEL_CGNR)
56 .value("lsqr", IterativeMatrixParams::ACCEL_LSQR)
57 .value("odir", IterativeMatrixParams::ACCEL_ODIR)
58 .value("omin", IterativeMatrixParams::ACCEL_OMIN)
59 .value("ores", IterativeMatrixParams::ACCEL_ORES)
60 .value("iom", IterativeMatrixParams::ACCEL_IOM)
61 .value("gmres", IterativeMatrixParams::ACCEL_GMRES)
62 .value("usymlq", IterativeMatrixParams::ACCEL_USYMLQ)
63 .value("usymqr", IterativeMatrixParams::ACCEL_USYMQR)
64 .value("landir", IterativeMatrixParams::ACCEL_LANDIR)
65 .value("lanmin", IterativeMatrixParams::ACCEL_LANMIN)
66 .value("lanres", IterativeMatrixParams::ACCEL_LANRES)
67 .value("cgcr", IterativeMatrixParams::ACCEL_CGCR)
68 .value("bcgs", IterativeMatrixParams::ACCEL_BCGS)
69 .get(iter_params.accelerator);
70 iter_params.preconditioner = reader.enumAttribute<IterativeMatrixParams::Preconditioner>("preconditioner")
71 .value("rich", IterativeMatrixParams::PRECOND_RICH)
72 .value("jac", IterativeMatrixParams::PRECOND_JAC)
73 .value("ljac", IterativeMatrixParams::PRECOND_LJAC)
74 .value("ljacx", IterativeMatrixParams::PRECOND_LJACX)
75 .value("sor", IterativeMatrixParams::PRECOND_SOR)
76 .value("ssor", IterativeMatrixParams::PRECOND_SSOR)
77 .value("ic", IterativeMatrixParams::PRECOND_IC)
78 .value("mic", IterativeMatrixParams::PRECOND_MIC)
79 .value("lsp", IterativeMatrixParams::PRECOND_LSP)
80 .value("neu", IterativeMatrixParams::PRECOND_NEU)
81 .value("lsor", IterativeMatrixParams::PRECOND_LSOR)
82 .value("lssor", IterativeMatrixParams::PRECOND_LSSOR)
83 .value("llsp", IterativeMatrixParams::PRECOND_LLSP)
84 .value("lneu", IterativeMatrixParams::PRECOND_LNEU)
85 .value("bic", IterativeMatrixParams::PRECOND_BIC)
86 .value("bicx", IterativeMatrixParams::PRECOND_BICX)
87 .value("mbic", IterativeMatrixParams::PRECOND_MBIC)
88 .value("mbicx", IterativeMatrixParams::PRECOND_MBICX)
89 .get(iter_params.preconditioner);
90 iter_params.no_convergence_behavior = reader.enumAttribute<IterativeMatrixParams::NoConvergenceBehavior>("noconv")
91 .value("error", IterativeMatrixParams::NO_CONVERGENCE_ERROR)
92 .value("warning", IterativeMatrixParams::NO_CONVERGENCE_WARNING)
93 .value("continue", IterativeMatrixParams::NO_CONVERGENCE_CONTINUE)
94 .get(iter_params.no_convergence_behavior);
95 iter_params.maxit = reader.getAttribute<int>("maxit", iter_params.maxit);
96 iter_params.maxerr = reader.getAttribute<double>("maxerr", iter_params.maxerr);
97 iter_params.nfact = reader.getAttribute<int>("nfact", iter_params.nfact);
98 iter_params.omega = reader.getAttribute<double>("omega", iter_params.omega);
99 iter_params.ndeg = reader.getAttribute<int>("ndeg", iter_params.ndeg);
100 iter_params.lvfill = reader.getAttribute<int>("lvfill", iter_params.lvfill);
101 iter_params.ltrunc = reader.getAttribute<int>("ltrunc", iter_params.ltrunc);
102 iter_params.ns1 = reader.getAttribute<int>("nsave", iter_params.ns1);
103 iter_params.ns2 = reader.getAttribute<int>("nrestart", iter_params.ns2);
104 reader.requireTagEnd(); // </iterative>
105 reader.requireTagEnd(); // </matrix>
106 }
107
108 return true;
109 }
110 return false;
111 }
112
113 inline FemMatrix* getMatrix();
114};
115
116template <typename SpaceT, typename MeshT> inline FemMatrix* FemSolverWithMesh<SpaceT, MeshT>::getMatrix() {
117 switch (algorithm) {
118 case ALGORITHM_CHOLESKY: return new DpbMatrix(this, this->mesh->size(), this->mesh->minorAxis()->size() + 1);
119 case ALGORITHM_GAUSS: return new DgbMatrix(this, this->mesh->size(), this->mesh->minorAxis()->size() + 1);
120 case ALGORITHM_ITERATIVE: return new SparseBandMatrix(this, this->mesh->size(), this->mesh->minorAxis()->size());
121 }
122 return nullptr;
123}
124
125template <> inline FemMatrix* FemSolverWithMesh<Geometry3D, RectangularMesh<3>>::getMatrix() {
126 size_t band = this->mesh->minorAxis()->size() * (this->mesh->mediumAxis()->size() + 1) + 1;
127 switch (algorithm) {
128 case ALGORITHM_CHOLESKY: return new DpbMatrix(this, this->mesh->size(), band);
129 case ALGORITHM_GAUSS: return new DgbMatrix(this, this->mesh->size(), band);
130 case ALGORITHM_ITERATIVE:
131 return new SparseBandMatrix(this, this->mesh->size(), mesh->mediumAxis()->size() * mesh->minorAxis()->size(),
132 mesh->minorAxis()->size());
133 }
134 return nullptr;
135}
136
138
144
145template <typename SpaceT, typename MeshT> struct FemSolverWithMaskedMesh : public FemSolverWithMesh<SpaceT, MeshT> {
146 static_assert(std::is_base_of<RectangularMesh<MeshT::DIM>, MeshT>::value,
147 "FemSolverWithMaskedMesh only works with RectangularMesh");
148
149 FemSolverWithMaskedMesh(const std::string& name = "") : FemSolverWithMesh<SpaceT, MeshT>(name) {}
150
151 protected:
152 plask::shared_ptr<RectangularMaskedMesh<MeshT::DIM>> maskedMesh = plask::make_shared<RectangularMaskedMesh<MeshT::DIM>>();
153 EmptyElementsHandling empty_elements = EMPTY_ELEMENTS_DEFAULT;
154
155 public:
157 EmptyElementsHandling getEmptyElements() const { return empty_elements; }
159 void setEmptyElements(EmptyElementsHandling val) {
160 empty_elements = val;
161 this->invalidate();
162 }
163
164 bool parseFemConfiguration(XMLReader& reader, Manager& manager) {
165 if (reader.getNodeName() == "mesh") {
166 if (reader.hasAttribute("include-empty")) {
167 this->writelog(LOG_WARNING, this->getId(), "Attribute 'include-empty' is deprecated, use 'empty-elements' instead");
168 empty_elements = reader.requireAttribute<bool>("include-empty") ? EMPTY_ELEMENTS_INCLUDED : EMPTY_ELEMENTS_EXCLUDED;
169 }
170 empty_elements = reader.enumAttribute<EmptyElementsHandling>("empty-elements")
171 .value("default", EMPTY_ELEMENTS_DEFAULT)
172 .value("exclude", EMPTY_ELEMENTS_EXCLUDED)
173 .value("include", EMPTY_ELEMENTS_INCLUDED)
174 .value("excluded", EMPTY_ELEMENTS_EXCLUDED)
175 .value("included", EMPTY_ELEMENTS_INCLUDED)
176 .get(empty_elements);
177 return false;
178 }
179 return FemSolverWithMesh<SpaceT, MeshT>::parseFemConfiguration(reader, manager);
180 }
181
182 void setupMaskedMesh() {
183 if (empty_elements == EMPTY_ELEMENTS_INCLUDED ||
184 (this->algorithm == ALGORITHM_ITERATIVE && empty_elements == EMPTY_ELEMENTS_DEFAULT)) {
185 maskedMesh->selectAll(*this->mesh);
186 } else {
187 maskedMesh->reset(*this->mesh, *this->geometry, ~plask::Material::EMPTY);
188 }
189 }
190
191 void onInitialize() { setupMaskedMesh(); }
192
193 inline FemMatrix* getMatrix();
194};
195
196template <typename SpaceT, typename MeshT> inline FemMatrix* FemSolverWithMaskedMesh<SpaceT, MeshT>::getMatrix() {
197 size_t band;
198 if (empty_elements == EMPTY_ELEMENTS_INCLUDED || this->algorithm == ALGORITHM_ITERATIVE) {
199 band = this->mesh->minorAxis()->size() + 1;
200 } else {
201 band = 0;
202 for (auto element : this->maskedMesh->elements()) {
203 size_t span = element.getUpUpIndex() - element.getLoLoIndex();
204 if (span > band) band = span;
205 }
206 }
207 switch (this->algorithm) {
208 case ALGORITHM_CHOLESKY: return new DpbMatrix(this, this->maskedMesh->size(), band);
209 case ALGORITHM_GAUSS: return new DgbMatrix(this, this->maskedMesh->size(), band);
210 case ALGORITHM_ITERATIVE:
211 if (empty_elements != EMPTY_ELEMENTS_EXCLUDED)
212 return new SparseBandMatrix(this, this->maskedMesh->size(), this->mesh->minorAxis()->size());
213 else
214 return new SparseFreeMatrix(this, this->maskedMesh->size(), this->maskedMesh->elements().size() * 10);
215 }
216 return nullptr;
217}
218
219template <> inline FemMatrix* FemSolverWithMaskedMesh<Geometry3D, RectangularMesh<3>>::getMatrix() {
220 size_t band;
221 if (empty_elements || algorithm == ALGORITHM_ITERATIVE) {
222 band = this->mesh->minorAxis()->size() * (this->mesh->mediumAxis()->size() + 1) + 1;
223 } else {
224 band = 0;
225 for (auto element : this->maskedMesh->elements()) {
226 size_t span = element.getUpUpUpIndex() - element.getLoLoLoIndex();
227 if (span > band) band = span;
228 }
229 }
230 switch (algorithm) {
231 case ALGORITHM_CHOLESKY: return new DpbMatrix(this, this->maskedMesh->size(), band);
232 case ALGORITHM_GAUSS: return new DgbMatrix(this, this->maskedMesh->size(), band);
233 case ALGORITHM_ITERATIVE:
234 if (empty_elements != EMPTY_ELEMENTS_EXCLUDED)
235 return new SparseBandMatrix(this, this->maskedMesh->size(), mesh->mediumAxis()->size() * mesh->minorAxis()->size(),
236 mesh->minorAxis()->size());
237 else
238 return new SparseFreeMatrix(this, this->maskedMesh->size(), this->maskedMesh->elements().size() * 36);
239 }
240 return nullptr;
241}
242
243} // namespace plask
244
245#endif