Line data Source code
1 : // Class FEMPoissonSolver
2 : // Solves the poisson equation using finite element methods and Conjugate
3 : // Gradient + a Preconditioner.
4 :
5 : #ifndef IPPL_PRECONFEMPOISSONSOLVER_H
6 : #define IPPL_PRECONFEMPOISSONSOLVER_H
7 :
8 : // #include "FEM/FiniteElementSpace.h"
9 : #include "LaplaceHelpers.h"
10 : #include "LinearSolvers/PCG.h"
11 : #include "Poisson.h"
12 :
13 : namespace ippl {
14 :
15 : template <typename Tlhs, unsigned Dim, unsigned numElemDOFs>
16 : struct EvalFunctor {
17 : const Vector<Tlhs, Dim> DPhiInvT;
18 : const Tlhs absDetDPhi;
19 :
20 0 : EvalFunctor(Vector<Tlhs, Dim> DPhiInvT, Tlhs absDetDPhi)
21 0 : : DPhiInvT(DPhiInvT)
22 0 : , absDetDPhi(absDetDPhi) {}
23 :
24 0 : KOKKOS_FUNCTION auto operator()(
25 : const size_t& i, const size_t& j,
26 : const Vector<Vector<Tlhs, Dim>, numElemDOFs>& grad_b_q_k) const {
27 0 : return dot((DPhiInvT * grad_b_q_k[j]), (DPhiInvT * grad_b_q_k[i])).apply() * absDetDPhi;
28 : }
29 : };
30 :
31 : /**
32 : * @brief A solver for the poisson equation using finite element methods and
33 : * Conjugate Gradient (CG)
34 : *
35 : * @tparam FieldLHS field type for the left hand side
36 : * @tparam FieldRHS field type for the right hand side
37 : */
38 : template <typename FieldLHS, typename FieldRHS = FieldLHS>
39 : class PreconditionedFEMPoissonSolver : public Poisson<FieldLHS, FieldRHS> {
40 : constexpr static unsigned Dim = FieldLHS::dim;
41 : using Tlhs = typename FieldLHS::value_type;
42 :
43 : public:
44 : using Base = Poisson<FieldLHS, FieldRHS>;
45 : using typename Base::lhs_type, typename Base::rhs_type;
46 : using MeshType = typename FieldRHS::Mesh_t;
47 :
48 : // PCG (Preconditioned Conjugate Gradient) is the solver algorithm used
49 : using PCGSolverAlgorithm_t =
50 : PCG<lhs_type, lhs_type, lhs_type, lhs_type, lhs_type, FieldLHS, FieldRHS>;
51 :
52 : // FEM Space types
53 : using ElementType =
54 : std::conditional_t<Dim == 1, ippl::EdgeElement<Tlhs>,
55 : std::conditional_t<Dim == 2, ippl::QuadrilateralElement<Tlhs>,
56 : ippl::HexahedralElement<Tlhs>>>;
57 :
58 : using QuadratureType = GaussJacobiQuadrature<Tlhs, 5, ElementType>;
59 :
60 : using LagrangeType = LagrangeSpace<Tlhs, Dim, 1, ElementType, QuadratureType, FieldLHS, FieldRHS>;
61 :
62 : // default constructor (compatibility with Alpine)
63 : PreconditionedFEMPoissonSolver()
64 : : Base()
65 : , refElement_m()
66 : , quadrature_m(refElement_m, 0.0, 0.0)
67 : , lagrangeSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)), Vector<Tlhs, Dim>(0),
68 : Vector<Tlhs, Dim>(0))), refElement_m, quadrature_m)
69 : {}
70 :
71 0 : PreconditionedFEMPoissonSolver(lhs_type& lhs, rhs_type& rhs)
72 : : Base(lhs, rhs)
73 : , refElement_m()
74 0 : , quadrature_m(refElement_m, 0.0, 0.0)
75 0 : , lagrangeSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
76 : static_assert(std::is_floating_point<Tlhs>::value, "Not a floating point type");
77 0 : setDefaultParameters();
78 :
79 : // start a timer
80 0 : static IpplTimings::TimerRef init = IpplTimings::getTimer("initFEM");
81 0 : IpplTimings::startTimer(init);
82 :
83 0 : rhs.fillHalo();
84 :
85 0 : lagrangeSpace_m.evaluateLoadVector(rhs);
86 :
87 0 : rhs.fillHalo();
88 :
89 0 : IpplTimings::stopTimer(init);
90 0 : }
91 :
92 0 : void setRhs(rhs_type& rhs) override {
93 0 : Base::setRhs(rhs);
94 :
95 0 : lagrangeSpace_m.initialize(rhs.get_mesh(), rhs.getLayout());
96 :
97 0 : rhs.fillHalo();
98 :
99 0 : lagrangeSpace_m.evaluateLoadVector(rhs);
100 :
101 0 : rhs.fillHalo();
102 0 : }
103 :
104 : /**
105 : * @brief Solve the poisson equation using finite element methods.
106 : * The problem is described by -laplace(lhs) = rhs
107 : */
108 0 : void solve() override {
109 : // start a timer
110 0 : static IpplTimings::TimerRef solve = IpplTimings::getTimer("solve");
111 0 : IpplTimings::startTimer(solve);
112 :
113 0 : const Vector<size_t, Dim> zeroNdIndex = Vector<size_t, Dim>(0);
114 :
115 : // We can pass the zeroNdIndex here, since the transformation jacobian does not depend
116 : // on translation
117 0 : const auto firstElementVertexPoints =
118 : lagrangeSpace_m.getElementMeshVertexPoints(zeroNdIndex);
119 :
120 : // Compute Inverse Transpose Transformation Jacobian ()
121 0 : const Vector<Tlhs, Dim> DPhiInvT =
122 : refElement_m.getInverseTransposeTransformationJacobian(firstElementVertexPoints);
123 :
124 : // Compute absolute value of the determinant of the transformation jacobian (|det D
125 : // Phi_K|)
126 0 : const Tlhs absDetDPhi = Kokkos::abs(
127 : refElement_m.getDeterminantOfTransformationJacobian(firstElementVertexPoints));
128 :
129 0 : EvalFunctor<Tlhs, Dim, LagrangeType::numElementDOFs> poissonEquationEval(
130 : DPhiInvT, absDetDPhi);
131 :
132 : // get BC type of our RHS
133 0 : BConds<FieldRHS, Dim>& bcField = (this->rhs_mp)->getFieldBC();
134 0 : FieldBC bcType = bcField[0]->getBCType();
135 :
136 0 : const auto algoOperator = [poissonEquationEval, &bcField, this](rhs_type field) -> lhs_type {
137 : // set appropriate BCs for the field as the info gets lost in the CG iteration
138 0 : field.setFieldBC(bcField);
139 :
140 0 : field.fillHalo();
141 :
142 0 : auto return_field = lagrangeSpace_m.evaluateAx(field, poissonEquationEval);
143 :
144 0 : return return_field;
145 : };
146 :
147 0 : const auto algoOperatorL = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
148 : // set appropriate BCs for the field as the info gets lost in the CG iteration
149 0 : field.setFieldBC(bcField);
150 :
151 0 : field.fillHalo();
152 :
153 0 : auto return_field = lagrangeSpace_m.evaluateAx_lower(field, poissonEquationEval);
154 :
155 0 : return return_field;
156 : };
157 :
158 0 : const auto algoOperatorU = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
159 : // set appropriate BCs for the field as the info gets lost in the CG iteration
160 0 : field.setFieldBC(bcField);
161 :
162 0 : field.fillHalo();
163 :
164 0 : auto return_field = lagrangeSpace_m.evaluateAx_upper(field, poissonEquationEval);
165 :
166 0 : return return_field;
167 : };
168 :
169 0 : const auto algoOperatorUL = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
170 : // set appropriate BCs for the field as the info gets lost in the CG iteration
171 0 : field.setFieldBC(bcField);
172 :
173 0 : field.fillHalo();
174 :
175 0 : auto return_field = lagrangeSpace_m.evaluateAx_upperlower(field, poissonEquationEval);
176 :
177 0 : return return_field;
178 : };
179 :
180 0 : const auto algoOperatorInvD = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
181 : // set appropriate BCs for the field as the info gets lost in the CG iteration
182 0 : field.setFieldBC(bcField);
183 :
184 0 : field.fillHalo();
185 :
186 0 : auto return_field = lagrangeSpace_m.evaluateAx_inversediag(field, poissonEquationEval);
187 :
188 0 : return return_field;
189 : };
190 :
191 0 : const auto algoOperatorD = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
192 : // set appropriate BCs for the field as the info gets lost in the CG iteration
193 0 : field.setFieldBC(bcField);
194 :
195 0 : field.fillHalo();
196 :
197 0 : auto return_field = lagrangeSpace_m.evaluateAx_diag(field, poissonEquationEval);
198 :
199 0 : return return_field;
200 : };
201 :
202 : // set preconditioner for PCG
203 0 : std::string preconditioner_type =
204 0 : this->params_m.template get<std::string>("preconditioner_type");
205 0 : int level = this->params_m.template get<int>("newton_level");
206 0 : int degree = this->params_m.template get<int>("chebyshev_degree");
207 0 : int inner = this->params_m.template get<int>("gauss_seidel_inner_iterations");
208 0 : int outer = this->params_m.template get<int>("gauss_seidel_outer_iterations");
209 0 : double omega = this->params_m.template get<double>("ssor_omega");
210 : int richardson_iterations =
211 0 : this->params_m.template get<int>("richardson_iterations");
212 :
213 0 : pcg_algo_m.setPreconditioner(algoOperator, algoOperatorL, algoOperatorU, algoOperatorUL,
214 : algoOperatorInvD, algoOperatorD, 0, 0, preconditioner_type,
215 : level, degree, richardson_iterations, inner, outer, omega);
216 :
217 0 : pcg_algo_m.setOperator(algoOperator);
218 :
219 : // send boundary values to RHS (load vector) i.e. lifting (Dirichlet BCs)
220 0 : if (bcType == CONSTANT_FACE) {
221 0 : *(this->rhs_mp) = *(this->rhs_mp) -
222 0 : lagrangeSpace_m.evaluateAx_lift(*(this->rhs_mp), poissonEquationEval);
223 : }
224 :
225 : // start a timer
226 0 : static IpplTimings::TimerRef pcgTimer = IpplTimings::getTimer("pcg");
227 0 : IpplTimings::startTimer(pcgTimer);
228 :
229 : // run PCG -> lhs contains solution
230 0 : pcg_algo_m(*(this->lhs_mp), *(this->rhs_mp), this->params_m);
231 :
232 : // added for BCs to be imposed properly
233 : // (they are not propagated through the preconditioner)
234 0 : if (bcType == CONSTANT_FACE) {
235 0 : bcField.assignGhostToPhysical(*(this->lhs_mp));
236 : }
237 0 : (this->lhs_mp)->fillHalo();
238 :
239 0 : IpplTimings::stopTimer(pcgTimer);
240 :
241 0 : int output = this->params_m.template get<int>("output_type");
242 0 : if (output & Base::GRAD) {
243 0 : *(this->grad_mp) = -grad(*(this->lhs_mp));
244 : }
245 :
246 0 : IpplTimings::stopTimer(solve);
247 0 : }
248 :
249 : /**
250 : * Query how many iterations were required to obtain the solution
251 : * the last time this solver was used
252 : * @return Iteration count of last solve
253 : */
254 0 : int getIterationCount() { return pcg_algo_m.getIterationCount(); }
255 :
256 : /**
257 : * Query the residue
258 : * @return Residue norm from last solve
259 : */
260 0 : Tlhs getResidue() const { return pcg_algo_m.getResidue(); }
261 :
262 : /**
263 : * Query the L2-norm error compared to a given (analytical) sol
264 : * @return L2 error after last solve
265 : */
266 : template <typename F>
267 0 : Tlhs getL2Error(const F& analytic) {
268 0 : Tlhs error_norm = this->lagrangeSpace_m.computeErrorL2(*(this->lhs_mp), analytic);
269 0 : return error_norm;
270 : }
271 :
272 : /**
273 : * Query the average of the solution
274 : * @param vol Boolean indicating whether we divide by volume or not
275 : * @return avg (offset for null space test cases if divided by volume)
276 : */
277 0 : Tlhs getAvg(bool Vol = false) {
278 0 : Tlhs avg = this->lagrangeSpace_m.computeAvg(*(this->lhs_mp));
279 0 : if (Vol) {
280 0 : lhs_type unit((this->lhs_mp)->get_mesh(), (this->lhs_mp)->getLayout());
281 0 : unit = 1.0;
282 0 : Tlhs vol = this->lagrangeSpace_m.computeAvg(unit);
283 0 : return avg/vol;
284 0 : } else {
285 0 : return avg;
286 : }
287 : }
288 :
289 : protected:
290 : PCGSolverAlgorithm_t pcg_algo_m;
291 :
292 0 : virtual void setDefaultParameters() override {
293 0 : this->params_m.add("max_iterations", 1000);
294 0 : this->params_m.add("tolerance", (Tlhs)1e-13);
295 0 : }
296 :
297 : ElementType refElement_m;
298 : QuadratureType quadrature_m;
299 : LagrangeType lagrangeSpace_m;
300 : };
301 :
302 : } // namespace ippl
303 :
304 : #endif
|
