Line data Source code
1 : // Class FEMMaxwellDifussionSolver
2 : // Solves the electric diffusion probelm given by curl(curl(E)) + E = f in
3 : // the domain and n x E = 0 on the boundary.
4 :
5 : #ifndef IPPL_FEM_MAXWELL_DIFUSSION_SOLVER_H
6 : #define IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
7 :
8 : #include "LinearSolvers/PCG.h"
9 : #include "Maxwell.h"
10 : #include <iomanip>
11 : #include <iostream>
12 : #include <fstream>
13 :
14 : namespace ippl {
15 :
16 : /**
17 : * @brief Representation of the lhs of the problem we are trying to solve.
18 : *
19 : * In our case this corresponds to the variational formulation of the
20 : * curl(curl(E)) + E and is curl(b_i)*curl(b_j) + b_i*b_j.
21 : *
22 : * @tparam T The type we are working with.
23 : * @tparam Dim the dimension of the space.
24 : * @tparam numElementDOFs the number of DOFs per element that we have.
25 : */
26 : template <typename T, unsigned Dim, unsigned numElementDOFs>
27 : struct EvalFunctor {
28 : /**
29 : * @brief The inverse transpose Jacobian.
30 : *
31 : * As we have a unirectangular grid it is the same for all the differnt
32 : * Elements and we therefore have to store it only once.
33 : */
34 : const Vector<T, Dim> DPhiInvT;
35 :
36 : /**
37 : * @brief The determinant of the Jacobian.
38 : *
39 : * As we have a unirectangular grid it is the same for all the differnt
40 : * Elements and we therefore have to store it only once.
41 : */
42 : const T absDetDPhi;
43 :
44 : /**
45 : * @brief Constructor.
46 : */
47 0 : EvalFunctor(Vector<T, Dim> DPhiInvT, T absDetDPhi)
48 0 : : DPhiInvT(DPhiInvT)
49 0 : , absDetDPhi(absDetDPhi) {}
50 :
51 : /**
52 : * @brief Returns the evaluation of
53 : * (curl(b_i)*curl(b_j) + b_i*b_j)*absDetDPhi.
54 : *
55 : * This function takes as input the basis function values and their curl
56 : * for the different DOFs and returns the evaluation of the inner part
57 : * of the integral of the variational formuation, which corresponds to
58 : * (curl(b_i)*curl(b_j) + b_i*b_j), but note that we additionally also
59 : * multiply this with absDetDPhi, which is required by the quadrature
60 : * rule. In theroy this could also be done outside of this.
61 : *
62 : * @param i The first DOF index.
63 : * @param j The second DOF index.
64 : * @param curl_b_q_k The curl of the DOFs.
65 : * @param val_b_q_k The values of the DOFs.
66 : *
67 : * @returns (curl(b_i)*curl(b_j) + b_i*b_j)*absDetDPhi
68 : */
69 0 : KOKKOS_FUNCTION const auto operator()(size_t i, size_t j,
70 : const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& curl_b_q_k,
71 : const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& val_b_q_k) const {
72 :
73 0 : T curlTerm = dot(DPhiInvT*curl_b_q_k[j], DPhiInvT*curl_b_q_k[i]).apply();
74 0 : T massTerm = dot(val_b_q_k[j], val_b_q_k[i]).apply();
75 0 : return (curlTerm + massTerm)*absDetDPhi;
76 : }
77 : };
78 :
79 : /**
80 : * @brief A solver for the electric diffusion equation given by
81 : * \f$ \nabla \times \nabla \times E + E = f \text{ in } \Omega\f$ and
82 : * \f$ n \times E = 0 \text{ on } \partial \Omega\f$ using the Nédélec basis
83 : * functions.
84 : *
85 : * @tparam FieldType The type used to represent a field on a mesh.
86 : */
87 : template <typename FieldType>
88 : class FEMMaxwellDiffusionSolver : public Maxwell<FieldType, FieldType> {
89 : constexpr static unsigned Dim = FieldType::dim;
90 :
91 : // we call value_type twice as in theory we expect the field to store
92 : // vector data represented by an ippl::Vector
93 : using T = typename FieldType::value_type::value_type;
94 :
95 : typedef Vector<T,Dim> point_t;
96 :
97 : public:
98 : using Base = Maxwell<FieldType, FieldType>;
99 : using MeshType = typename FieldType::Mesh_t;
100 :
101 : // PCG (Preconditioned Conjugate Gradient) is the solver algorithm used
102 : using PCGSolverAlgorithm_t = CG<FEMVector<T>, FEMVector<T>, FEMVector<T>, FEMVector<T>,
103 : FEMVector<T>, FEMVector<T>, FEMVector<T>>;
104 :
105 : // FEM Space types
106 : using ElementType = std::conditional_t<Dim == 2, ippl::QuadrilateralElement<T>,
107 : ippl::HexahedralElement<T>>;
108 :
109 : using QuadratureType = GaussJacobiQuadrature<T, 5, ElementType>;
110 :
111 : using NedelecType = NedelecSpace<T, Dim, 1, ElementType, QuadratureType, FieldType>;
112 :
113 : // default constructor (compatibility with Alpine)
114 : FEMMaxwellDiffusionSolver()
115 : : Base()
116 : , rhsVector_m(nullptr)
117 : , refElement_m()
118 : , quadrature_m(refElement_m, 0.0, 0.0)
119 : , nedelecSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)), Vector<T, Dim>(0),
120 : Vector<T, Dim>(0))), refElement_m, quadrature_m)
121 : {}
122 :
123 0 : FEMMaxwellDiffusionSolver(FieldType& lhs, FieldType& rhs, const FEMVector<point_t>& rhsVector)
124 : : Base(lhs, lhs, rhs)
125 0 : , rhsVector_m(nullptr)
126 : , refElement_m()
127 0 : , quadrature_m(refElement_m, 0.0, 0.0)
128 0 : , nedelecSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
129 :
130 : static_assert(std::is_floating_point<T>::value, "Not a floating point type");
131 0 : setDefaultParameters();
132 :
133 : // Calcualte the rhs, using the Nedelec space
134 0 : rhsVector_m =
135 0 : std::make_unique<FEMVector<T>>(nedelecSpace_m.evaluateLoadVector(rhsVector));
136 :
137 0 : rhsVector_m->accumulateHalo();
138 0 : rhsVector_m->fillHalo();
139 :
140 0 : }
141 :
142 : void setRhs(FieldType& rhs, const FEMVector<point_t>& rhsVector) {
143 :
144 : Base::setRhs(rhs);
145 :
146 : // Calcualte the rhs, using the Nedelec space
147 : rhsVector_m =
148 : std::make_unique<FEMVector<T>>(nedelecSpace_m.evaluateLoadVector(rhsVector));
149 :
150 : rhsVector_m->accumulateHalo();
151 : rhsVector_m->fillHalo();
152 : }
153 :
154 : /**
155 : * @brief Solve the equation using finite element methods.
156 : */
157 0 : void solve() override {
158 :
159 0 : const Vector<size_t, Dim> zeroNdIndex = Vector<size_t, Dim>(0);
160 :
161 : // We can pass the zeroNdIndex here, since the transformation
162 : // jacobian does not depend on translation
163 0 : const auto firstElementVertexPoints =
164 : nedelecSpace_m.getElementMeshVertexPoints(zeroNdIndex);
165 :
166 : // Compute Inverse Transpose Transformation Jacobian ()
167 0 : const Vector<T, Dim> DPhiInvT =
168 : refElement_m.getInverseTransposeTransformationJacobian(firstElementVertexPoints);
169 :
170 : // Compute absolute value of the determinant of the transformation
171 : // jacobian (|det D Phi_K|)
172 0 : const T absDetDPhi = Kokkos::abs(
173 : refElement_m.getDeterminantOfTransformationJacobian(firstElementVertexPoints));
174 :
175 : // Create the functor object which stores the function we have to
176 : // solve for the lhs
177 0 : EvalFunctor<T, Dim, this->nedelecSpace_m.numElementDOFs> maxwellDiffusionEval(
178 : DPhiInvT, absDetDPhi);
179 :
180 : // The Ax operator
181 0 : const auto algoOperator = [maxwellDiffusionEval, this](FEMVector<T> vector)
182 : -> FEMVector<T> {
183 :
184 0 : vector.fillHalo();
185 :
186 0 : FEMVector<T> return_vector = nedelecSpace_m.evaluateAx(vector,maxwellDiffusionEval);
187 :
188 0 : return_vector.accumulateHalo();
189 :
190 0 : return return_vector;
191 0 : };
192 :
193 : // setup the CG solver
194 0 : pcg_algo_m.setOperator(algoOperator);
195 :
196 : // Create the coefficient vector for the solution
197 0 : FEMVector<T> lhsVector = rhsVector_m->deepCopy();
198 :
199 : // Solve the system using CG
200 : try {
201 0 : pcg_algo_m(lhsVector, *rhsVector_m, this->params_m);
202 0 : } catch (IpplException& e) {
203 0 : std::string msg = e.where() + ": " + e.what() + "\n";
204 0 : Kokkos::abort(msg.c_str());
205 0 : }
206 :
207 : // store solution.
208 0 : lhsVector_m = std::make_unique<FEMVector<T>>(lhsVector);
209 :
210 : // set the boundary values to the correct values.
211 0 : lhsVector.fillHalo();
212 :
213 0 : }
214 :
215 : /**
216 : * Query how many iterations were required to obtain the solution
217 : * the last time this solver was used
218 : * @return Iteration count of last solve
219 : */
220 0 : int getIterationCount() { return pcg_algo_m.getIterationCount(); }
221 :
222 : /**
223 : * Query the residue
224 : * @return Residue norm from last solve
225 : */
226 0 : T getResidue() const { return pcg_algo_m.getResidue(); }
227 :
228 : /**
229 : * @brief Reconstructs function values at arbitrary points in the mesh.
230 : *
231 : * This function can be used to retrieve the values of a solution
232 : * function at arbitrary points inside of the mesh.
233 : *
234 : * @note Currently the function is able to handle both cases, where we
235 : * have that \p positions only contains positions which are inside of
236 : * local domain of this MPI rank (i.e. each rank gets its own unique
237 : * \p position ) and where \p positions contains the positions of all
238 : * ranks (i.e. \p positions is the same for all ranks). If in the future
239 : * it can be guaranteed, that each rank will get its own \p positions
240 : * then certain parts of the function implementation can be removed.
241 : * Instructions for this are given in the implementation itself.
242 : *
243 : * @param positions The points at which the function should be
244 : * evaluated. A \c Kokkos::View which stores in each element a 2D/3D
245 : * point.
246 : *
247 : * @return The function evaluated at the given points, stored inside of
248 : * \c Kokkos::View where each element corresponts to the function value
249 : * at the point described by the same element inside of \p positions.
250 : */
251 0 : Kokkos::View<point_t*> reconstructToPoints(const Kokkos::View<point_t*>& positions) {
252 0 : return this->nedelecSpace_m.reconstructToPoints(positions, *lhsVector_m);
253 : }
254 :
255 :
256 :
257 : /**
258 : * @brief Given an analytical solution computes the L2 norm error.
259 : *
260 : * @param analytical The analytical solution (functor)
261 : *
262 : * @return error - The error ||u - analytical||_L2
263 : */
264 : template <typename F>
265 0 : T getL2Error(const F& analytic) {
266 0 : T error_norm = this->nedelecSpace_m.computeError(*lhsVector_m, analytic);
267 0 : return error_norm;
268 : }
269 :
270 :
271 : protected:
272 :
273 : /**
274 : * @brief The CG Solver we use
275 : */
276 : PCGSolverAlgorithm_t pcg_algo_m;
277 :
278 : /**
279 : * @brief Sets the default values for the CG solver.
280 : * Defaults are: max Iterations = 10, tolerance = 1e-13
281 : */
282 0 : virtual void setDefaultParameters() {
283 0 : this->params_m.add("max_iterations", 10);
284 0 : this->params_m.add("tolerance", (T)1e-13);
285 0 : }
286 :
287 : /**
288 : * @brief FEM represenation of the rhs
289 : * We use this to store the rhs b of the System Ax = b used in the
290 : * Galerkin FEM scheme.
291 : */
292 : std::unique_ptr<FEMVector<T>> rhsVector_m;
293 :
294 : /**
295 : * @brief FEM represenation of the solution vector
296 : * We use this to store the solution x of the System Ax = b used in the
297 : * Galerkin FEM scheme.
298 : */
299 : std::unique_ptr<FEMVector<T>> lhsVector_m;
300 :
301 : /**
302 : * @brief the reference element we have.
303 : */
304 : ElementType refElement_m;
305 :
306 : /**
307 : * @brief The quadrature rule we use.
308 : */
309 : QuadratureType quadrature_m;
310 :
311 : /**
312 : * @brief The Nedelec Space object.
313 : *
314 : * This is the representation of the Nedelec space that we have and
315 : * which we use to interact with all the Nedelec stuff.
316 : */
317 : NedelecType nedelecSpace_m;
318 : };
319 :
320 : } // namespace ippl
321 :
322 :
323 :
324 : #endif // IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
|
