1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Eric1; Nguyen, Manh-Ha; Chambeyron, Colin
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, either version 3 of the License, or
8 (at your option) any later version.
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 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18
19 \file dev_Eric1.cpp
20 \author E. Lunéville
21 \since 1 March 2014
22 \date 1 March 2014
23 */
24
25 #include "xlife++-libs.h"
26 #include "testUtils.hpp"
27
28 using namespace xlifepp;
29
30 namespace dev_Eric1
31 {
un(const Point & P,Parameters & pa=defaultParameters)32 Real un(const Point& P, Parameters& pa = defaultParameters)
33 {
34 return 1.;
35 }
36 /********************************************************************************************************/
37
nonZeros(const TermMatrix & A)38 Number nonZeros(const TermMatrix& A)
39 {
40 Number NZ=0;
41 SuTermMatrix* suterm = A.begin()->second;
42 LargeMatrix<Complex>* cA = suterm->entries()->cEntries_p;
43 if (0 != cA) NZ = cA->nbNonZero();
44
45 LargeMatrix<Real>* rA = suterm->entries()->rEntries_p;
46 if (0 != rA) NZ = rA->nbNonZero();
47
48 LargeMatrix<Matrix<Real> >* rmA = suterm->entries()->rmEntries_p;
49 if (0 != rmA) NZ = rmA->nbNonZero();
50
51 return NZ;
52 }
53
printResultsOfCV(const TermMatrix & MAT,const TermVector & B,const TermVector & SOL,const TermVector & Uex,Real epsilon,std::ostream & out)54 void printResultsOfCV( const TermMatrix& MAT, const TermVector& B, const TermVector& SOL, const TermVector& Uex, Real epsilon, std::ostream& out)
55 {
56 TermVector E, D;
57 E=B;
58 E -= MAT*SOL;
59 D = SOL; D -= Uex;
60 if ( (norminfty(E) < epsilon) && (norminfty(D) < epsilon))
61 {
62 std::cout << "CV with epsilon = " << epsilon << eol;
63 std::cout << "||residu|| = " << norminfty(E) << eol << "||Uex-Unum|| = " << norminfty(D) << _nbIterations<<eol << eol;
64 }
65 else
66 {
67 std::cout << "NO CV with epsilon = " << epsilon << eol;
68 std::cout << "||residu|| = " << norminfty(E) << eol << "||Uex-Unum|| = " << norminfty(D) << _nbIterations<< eol << eol;
69 }
70 }
71
72
ISEpsilon_Pk(const GeomDomain & omg,int k,Real h,std::ostream & out)73 void ISEpsilon_Pk(const GeomDomain& omg, int k, Real h, std::ostream& out)
74 {
75 std::stringstream s;
76 trace_p->push("Epsilon_Pk");
77 verboseLevel(20);
78
79 // Real epsil=0.0001;
80 //create Lagrange Pind space and unknown
81 Interpolation& Pk=interpolation(Lagrange, standard, k, H1);
82 Space Vk(omg, Pk, "Vk", true);
83 // Create unknows and test function
84 Unknown u(Vk, "u", 2);
85 TestFunction v(u, "v");
86
87 // Create bilinear form
88 BilinearForm auv=intg( omg, epsilon(u)%epsilon(v) ) + intg(omg, div(u) * div(v) ) + intg(omg, u | v);
89 // BilinearForm auv=intg( omg, grad(u)|grad(v))+ intg(omg, u * v);
90 // BilinearForm muv= intg(omg, u | v);
91 // BilinearForm euv= intg(omg, epsilon(u)%epsilon(v));
92 // BilinearForm duv= intg(omg, div(u)*div(v));
93
94 // Create and compute term matrix
95 TermMatrix A(auv, "a(u,v)");
96 TermMatrix Atmp(A);
97 TermMatrix ACS(auv, _csDualStorage, "a(u,v)");
98 TermMatrix ACStmp(ACS);
99 TermMatrix ACSS(auv, _csSymStorage, "a(u,v) sym");
100 TermMatrix ACSStmp(ACSS);
101 //
102 Number NbP = A.numberOfRows();//begin()->second->entries()->rmEntries_p->nbRows;
103 Number NbPs = A.numberOfRows();//begin()->second->entries()->rmEntries_p->nbRowsSub;
104 Number NbCoef=NbP*NbPs;
105 std::cout<<"Matrice de "<<NbP<<"x"<<NbP<<" blocs:"<<NbPs<<"x"<<NbPs <<eol;
106 std::cout<<"nb non-zero coeficients of the matrix: "<<nonZeros(A)<<eol<<eol;
107
108 //================================================================================================
109
110 TermVector rvX0(u,omg,0.,"vecteur initial");
111 TermVector rvX1(u,omg,1.5,"vecteur initial");
112 TermVector S, UN(u,omg,2.,"un"),SS;
113 TermMatrix III(UN, "identite"); //compute(III);
114 S = A*UN;
115 SS = ACS*UN;
116 // Parameters
117 const Number krylovDim = NbP;
118 const Real epsilon = 1.e-6;
119 const Real epsilonP = 1.e-6;
120 const Number nbIteration = 1000;
121 //const Number vb = 2;
122
123
124 /////////////////////////////////////////////////
125 // BicG
126 /////////////////////////////////////////////////
127 std::cout <<cpuTime();
128 std::cout << "* BiCG without Preconditionner" << eol;
129 TermVector rvBicg = bicgSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=nbIteration);
130 std::cout << cpuTime("BiGG without preconditioner")<<std::endl<<std::endl;
131 printResultsOfCV(A, S, rvBicg, UN, epsilon, out);
132 clear(rvBicg);
133
134 A=Atmp;
135 std::cout << "* BiCG with Preconditionner (diag)" << eol;
136 cpuTime("");
137 PreconditionerTerm preMatBi(A, _diagPrec);
138 TermVector rvBicgP = bicgSolve(A, S, rvX0, preMatBi, _tolerance=epsilonP, _maxIt=nbIteration);
139 std::cout<< cpuTime("BiGG with preconditioner diag")<<std::endl;
140 printResultsOfCV(A, S, rvBicgP, UN, epsilonP, out);
141 clear(rvBicgP);
142
143 A=Atmp;
144 std::cout << "* BiCG with Preconditionner (LDLT)" << eol;
145 std::cout <<cpuTime();
146 PreconditionerTerm preMatBi2(A, _ldltPrec);
147 TermVector rvBicgP2 = bicgSolve(A, S, rvX0, preMatBi2, _tolerance=2*epsilonP, _maxIt=nbIteration);
148 std::cout<< cpuTime("BiGG with preconditioner LDLt")<<std::endl;;
149 printResultsOfCV(A, S, rvBicgP2, UN, 2*epsilonP, out);
150 clear(rvBicgP2);
151
152 A=Atmp;
153 std::cout << "* BiCG with Preconditionner (iLDLT)" << eol;
154 std::cout <<cpuTime();
155 PreconditionerTerm preMatBi4(A, _ildltPrec);
156 TermVector rvBicgP4 = bicgSolve(A, S, rvX0, preMatBi4, _tolerance=2*epsilonP, _maxIt=nbIteration);
157 std::cout << cpuTime("BiGG with preconditioner iLDLt")<<std::endl;
158 printResultsOfCV(A, S, rvBicgP4, UN, 2*epsilonP,out);
159 clear(rvBicgP4);
160
161 ACSS=ACSStmp;
162 std::cout << "* BiCG with incomplete Preconditionner (iLDLT) SymCs" << eol;
163 std::cout <<cpuTime();
164 PreconditionerTerm preMatBi3(ACSS, _ildltPrec);
165 TermVector rvBicgP3 = bicgSolve(ACSS, S, rvX0, preMatBi3, _tolerance=2*epsilonP, _maxIt=nbIteration);
166 std::cout << cpuTime("BiGG with preconditioner iLDLt")<<std::endl;
167 printResultsOfCV(ACSS, S, rvBicgP3, UN, 2*epsilonP,out);
168 clear(rvBicgP3);
169
170 ACSS=ACSStmp;
171 std::cout << "* BiCG with incomplete Preconditionner (iLLt)" << eol;
172 std::cout <<cpuTime();
173 PreconditionerTerm preMatBi5(ACSS, _illtPrec);
174 TermVector rvBicgP5 = bicgSolve(ACSS, S, rvX0, preMatBi5, _tolerance=2*epsilonP, _maxIt=nbIteration);
175 std::cout << cpuTime("BiGG with preconditioner iLLt")<<std::endl;
176 printResultsOfCV(ACSS, S, rvBicgP5, UN, 2*epsilonP, out);
177 clear(rvBicgP5);
178
179 // no transpose solver
180 // ACSS=ACSStmp;
181 // out << "* BiCG with incomplete Preconditionner (iLU)" << eol;
182 // std::cout << "* BiCG with incomplete Preconditionner (iLU)" << eol;
183 // PreconditionerTerm preMatBi4(ACSS, _iluPrec);
184 // TermVector rvBicgP4 = bicgSolve(ACSS, S, rvX0, preMatBi4, _tolerance=epsilonP, _maxIt=nbIteration);
185 // printResultsOfCV(ACSS, S, rvBicgP4, UN, epsilonP, out);
186 // clear(rvBicgP4);
187
188 //
189 /////////////////////////////////////////////////
190 // CG
191 /////////////////////////////////////////////////
192
193 A=Atmp;
194 std::cout << "* CG without Preconditionner" << eol;
195 TermVector rvCg = cgSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=nbIteration);
196 std::cout << cpuTime("CG without preconditioner") << std::endl;
197 printResultsOfCV(A, S, rvCg, UN, epsilon, out);
198 clear(rvCg);
199
200 A=Atmp;
201 std::cout << "* CG with Preconditionner (diag)" << eol;
202 PreconditionerTerm preMat(A, _diagPrec);
203 TermVector rvCgP = cgSolve(A, S, rvX0, preMat, _tolerance=epsilonP, _maxIt=nbIteration);
204 std::cout << cpuTime("CG with preconditioner diag") << std::endl;
205 printResultsOfCV(A, S, rvCgP, UN, epsilonP, out);
206 clear(rvCgP);
207
208 A=Atmp;
209 std::cout << "* CG with Preconditionner (LDLT)" << eol;
210 PreconditionerTerm preMat2(A, _ldltPrec);
211 TermVector rvCgP2 = cgSolve(A, S, rvX0, preMat2, _tolerance=epsilonP, _maxIt=nbIteration);
212 std::cout<<cpuTime("CG with preconditioner LDLT")<< std::endl;
213 printResultsOfCV(A, S, rvCgP2, UN, epsilonP, out);
214 clear(rvCgP2);
215
216 ACSS=ACSStmp;
217 std::cout << "* CG with incomplete Preconditionner (iLDLT)" << eol;
218 PreconditionerTerm preMat3(ACSS, _ildltPrec);
219 TermVector rvCgP3 = cgSolve(ACSS, S, rvX0, preMat3, _tolerance=epsilonP, _maxIt=nbIteration);
220 std::cout << cpuTime("CG with preconditioner iLDLT mat sym") << std::endl;
221 printResultsOfCV(ACSS, S, rvCgP3, UN, epsilonP, out);
222 clear(rvCgP3);
223 /*pourquoi? -> matrice non symq*/
224 // ACS=ACStmp;
225 // out << "* CG with incomplete Preconditionner (iLDLT)" << eol;
226 // std::cout << "* CG with incomplete Preconditionner (iLDLT)" << ACS.symmetry() << eol;
227 // PreconditionerTerm preMat3bis(ACS, _ildltPrec);
228 // TermVector rvCgP3bis = cgSolve(ACS, S, rvX0, preMat3bis, _tolerance=epsilonP, _maxIt=nbIteration);
229 // printResultsOfCV(ACS, S, rvCgP3bis, UN, epsilonP, out);
230 // clear(rvCgP3);
231
232 ACS=ACStmp;
233 std::cout << "* CG with incomplete Preconditionner (iLU)" << eol;
234 PreconditionerTerm preMat4(ACS, _iluPrec);
235 TermVector rvCgP4 = cgSolve(ACS, S, rvX0, preMat4, _tolerance=epsilonP, _maxIt=nbIteration);
236 std::cout << cpuTime("CG with preconditioner iLU") << std::endl;
237 printResultsOfCV(ACS, S, rvCgP4, UN, epsilonP, out);
238 clear(rvCgP4);
239
240 ACSS=ACSStmp;
241 std::cout << "* CG with incomplete Preconditionner (iLLt)" << eol;
242 PreconditionerTerm preMat5(ACSS, _illtPrec);
243 TermVector rvCgP5 = cgSolve(ACSS, S, rvX0, preMat5, _tolerance=epsilonP, _maxIt=nbIteration);
244 std::cout << cpuTime("CG with preconditioner iLLT") << std::endl;
245 printResultsOfCV(ACSS, S, rvCgP5, UN, epsilonP, out);
246 clear(rvCgP5);
247
248 /////////////////////////////////////////////////
249 // CGS
250 /////////////////////////////////////////////////
251
252 // ##### CGS DOES NOT CONVERGE, perhaps a case of erratic behaviour
253
254 A=Atmp;
255 std::cout << "* CGS without Preconditionner" << eol;
256 TermVector rvCgS = cgsSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=nbIteration);
257 std::cout << cpuTime("CGS without preconditioner ") << std::endl;
258 printResultsOfCV(A, S, rvCgS, UN, epsilon, out);
259 clear( rvCgS);
260
261 A=Atmp;
262 std::cout << "* CGS with Preconditionner (diag)" << eol;
263 PreconditionerTerm preMatS(A, _diagPrec);
264 TermVector rvCgSP = cgsSolve(A, S, rvX0, preMatS, _tolerance=epsilonP, _maxIt=nbIteration);
265 std::cout << cpuTime("CGS with preconditioner diag ") << std::endl;
266 printResultsOfCV(A, S, rvCgSP, UN, epsilonP, out);
267 clear(rvCgSP);
268
269 A=Atmp;
270 std::cout << "* CGS with Preconditionner (LDLT)" << eol;
271 PreconditionerTerm preMatS2(A, _ldltPrec);
272 TermVector rvCgSP2 = cgsSolve(A, S, rvX1, preMatS2, _tolerance=epsilonP, _maxIt=nbIteration);
273 std::cout << cpuTime("CGS with preconditioner LDLt ") << std::endl;
274 printResultsOfCV(A, S, rvCgSP2, UN, epsilonP, out);
275 clear(rvCgSP2);
276
277
278 ACSS=ACSStmp;
279 std::cout << "* CGS with incomplete Preconditionner (iLDLT)" << eol;
280 PreconditionerTerm preMatS3(ACSS, _ildltPrec);
281 TermVector rvCgSP3 = cgsSolve(ACSS, S, rvX1, preMatS3, _tolerance=epsilonP, _maxIt=nbIteration);
282 std::cout << cpuTime("CGS with preconditioner iLDLt ") << std::endl;
283 printResultsOfCV(ACSS, S, rvCgSP3, UN, epsilonP, out);
284 clear(rvCgSP3);
285
286 ACS=ACStmp;
287 std::cout << "* CGS with incomplete Preconditionner (iLU)" << eol;
288 PreconditionerTerm preMatS4(ACS, _iluPrec);
289 TermVector rvCgSP4 = cgsSolve(ACS, S, rvX1, preMatS4, _tolerance=epsilonP, _maxIt=nbIteration);
290 std::cout << cpuTime("CGS with preconditioner iLU ") << std::endl;
291 printResultsOfCV(ACS, S, rvCgSP4, UN, epsilonP, out);
292 clear(rvCgSP4);
293
294 ACSS=ACSStmp;
295 std::cout << "* CGS with incomplete Preconditionner (iLLt)" << eol;
296 PreconditionerTerm preMatS5(ACSS, _illtPrec);
297 TermVector rvCgSP5 = cgsSolve(ACSS, S, rvX1, preMatS5, _tolerance=epsilonP, _maxIt=nbIteration);
298 std::cout << cpuTime("CGS with preconditioner iLLt ") << std::endl;
299 printResultsOfCV(ACSS, S, rvCgSP5, UN, epsilonP, out);
300 clear(rvCgSP5);
301
302 /////////////////////////////////////////////////
303 // GMRES
304 /////////////////////////////////////////////////
305
306 A=Atmp;
307 std::cout << "* Gmres without Preconditionner" << eol;
308 TermVector rvG = gmresSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=nbIteration, _krylovDim=krylovDim);
309 std::cout << cpuTime("GMRES without preconditioner ") << std::endl;
310 printResultsOfCV(A, S, rvG, UN, epsilon, out);
311 clear(rvG);
312
313 A=Atmp;
314 std::cout << "* Gmres with Preconditionner (diag)" << eol;
315 PreconditionerTerm preMatGmres(A, _diagPrec);
316 TermVector rvGP = gmresSolve(A, S, rvX1, preMatGmres, _tolerance=epsilonP, _maxIt=nbIteration, _krylovDim=krylovDim);
317 std::cout << cpuTime("GMRES with preconditioner diag") << std::endl;
318 printResultsOfCV(A, S, rvGP, UN, epsilonP, out);
319 clear(rvGP);
320
321 A=Atmp;
322 std::cout << "* Gmres with Preconditionner (LDLT)" << eol;
323 PreconditionerTerm preMatGmres2(A, _ldltPrec);
324 TermVector rvGP2 = gmresSolve(A, S, rvX0, preMatGmres2, _tolerance=epsilonP, _maxIt=nbIteration, _krylovDim=krylovDim);
325 std::cout << cpuTime("GMRES with preconditioner LDLt") << std::endl;
326 printResultsOfCV(A, S, rvGP2, UN, epsilonP, out);
327 clear(rvGP2);
328
329 ACSS=ACSStmp;
330 std::cout << "* GMRES with incomplete Preconditionner (iLDLT)" << eol;
331 PreconditionerTerm preMatGmres3(ACSS, _ildltPrec);
332 TermVector rvGP3 = gmresSolve(ACSS, S, rvX1, preMatGmres3, _tolerance=epsilonP, _maxIt=nbIteration);
333 std::cout << cpuTime("GMRES with preconditioner iLDLt") << std::endl;
334 printResultsOfCV(ACSS, S, rvGP3, UN, epsilonP, out);
335 clear(rvGP3);
336
337 ACS=ACStmp;
338 std::cout << "* GMRES with incomplete Preconditionner (iLU)" << eol;
339 PreconditionerTerm preMatGmres4(ACS, _iluPrec);
340 TermVector rvGP4 = gmresSolve(ACS, S, rvX1, preMatGmres4, _tolerance=epsilonP, _maxIt=nbIteration);
341 std::cout << cpuTime("GMRES with preconditioner iLU") <<std::endl;
342 printResultsOfCV(ACS, S, rvGP4, UN, epsilonP, out);
343 clear(rvGP4);
344
345 ACSS=ACSStmp;
346 std::cout << "* GMRES with incomplete Preconditionner (iLLt)" << eol;
347 PreconditionerTerm preMatGmres5(ACSS, _illtPrec);
348 TermVector rvGP5 = gmresSolve(ACSS, S, rvX1, preMatGmres5, _tolerance=epsilonP, _maxIt=nbIteration);
349 std::cout << cpuTime("GMRES with preconditioner iLLt")<<std::endl;
350 printResultsOfCV(ACSS, S, rvGP5, UN, epsilonP, out);
351 clear(rvGP5);
352
353
354
355 /////////////////////////////////////////////////
356 // BiCGStab
357 /////////////////////////////////////////////////
358
359 A=Atmp;
360 std::cout << "* BiCGstab without Preconditionner" << eol;
361 TermVector rvBicgStab = bicgStabSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=10*nbIteration);
362 std::cout << cpuTime("BiCgSrab without preconditioner ") << std::endl;
363 printResultsOfCV(A, S, rvBicgStab, UN, epsilon, out);
364 clear(rvBicgStab);
365
366 A=Atmp;
367 std::cout << "* BiCGstab with Preconditionner (diag)" << eol;
368 PreconditionerTerm preMatBiS(A, _diagPrec);
369 TermVector rvBicgStabP = bicgStabSolve(A, S, rvX0, preMatBiS, _tolerance=epsilonP, _maxIt=nbIteration);
370 std::cout << cpuTime("BiCgSrab with preconditioner diag ")<< std::endl;
371 printResultsOfCV(A, S, rvBicgStabP, UN, epsilonP, out);
372 clear(rvBicgStabP);
373
374 A=Atmp;
375 std::cout << "* BiCGstab with Preconditionner (LDLT)" << eol;
376 PreconditionerTerm preMatBiS2(A, _ldltPrec);
377 TermVector rvBicgStabP2 = bicgStabSolve(A, S, rvX0, preMatBiS2, _tolerance=epsilonP, _maxIt=nbIteration);
378 std::cout << cpuTime("BiCgSrab with preconditioner LDLt ") << std::endl;
379 printResultsOfCV(A, S, rvBicgStabP2, UN, epsilonP, out);
380 clear(rvBicgStabP2);
381
382
383 ACSS=ACSStmp;
384 std::cout << "* BiCGstab with incomplete Preconditionner (iLDLT)" << eol;
385 PreconditionerTerm preMatBiS3(ACSS, _ildltPrec);
386 TermVector rvBicgStabP3 = bicgStabSolve(ACSS, S, rvX1, preMatBiS3, _tolerance=2*epsilonP, _maxIt=2*nbIteration);
387 std::cout << cpuTime("BiCgSrab with preconditioner iLDLt ") << std::endl;
388 printResultsOfCV(ACSS, S, rvBicgStabP3, UN, 2*epsilonP, out);
389 clear(rvBicgStabP3);
390
391
392 ACS=ACStmp;
393 std::cout << "* BiCGstab with incomplete Preconditionner (iLU)" << eol;
394 PreconditionerTerm preMatBiS4(ACS, _iluPrec);
395 TermVector rvBicgStabP4 = bicgStabSolve(ACS, S, rvX1, preMatBiS4, _tolerance=2*epsilonP, _maxIt=nbIteration);
396 std::cout << cpuTime("BiCgSrab with preconditioner iLU ") << std::endl;
397 printResultsOfCV(ACS, S, rvBicgStabP4, UN, 2*epsilonP, out);
398 clear(rvBicgStabP4);
399
400
401 ACSS=ACSStmp;
402 std::cout << "* BiCGstab with incomplete Preconditionner (iLLt)" << eol;
403 PreconditionerTerm preMatBiS5(ACSS, _illtPrec);
404 TermVector rvBicgStabP5 = bicgStabSolve(ACSS, S, rvX1, preMatBiS5, _tolerance=2*epsilonP, _maxIt=nbIteration);
405 std::cout << cpuTime("BiCgSrab with preconditioner iLLt ") << std::endl;
406 printResultsOfCV(ACS, S, rvBicgStabP5, UN, 2*epsilonP, out);
407 clear(rvBicgStabP5);
408
409 /////////////////////////////////////////////////
410 // QMR //
411 /////////////////////////////////////////////////
412
413 A=Atmp;
414 std::cout << "* QMR without Preconditionner" << eol;
415 TermVector rvQmr = qmrSolve(A, S, rvX0, _tolerance=epsilon, _maxIt=nbIteration);
416 std::cout << cpuTime("QMR without preconditioner ") << std::endl;
417 printResultsOfCV(A, S, rvQmr, UN, epsilon, out);
418 clear(rvQmr);
419
420 A=Atmp;
421 std::cout << "* QMR with Preconditionner (diag)" << eol;
422 PreconditionerTerm preQ(A, _diagPrec);
423 TermVector rvQmrP = qmrSolve(A, S, rvX0, preQ, _tolerance=epsilonP, _maxIt=nbIteration);
424 std::cout << cpuTime("QMR with preconditioner diag ") << std::endl;
425 printResultsOfCV(A, S, rvQmrP, UN, epsilonP, out);
426 clear(rvQmrP);
427
428 A=Atmp;
429 std::cout << "* QMR with Preconditionner (LDLT)" << eol;
430 PreconditionerTerm preQ2(A, _ldltPrec);
431 TermVector rvQmrP2 = qmrSolve(A, S, rvX0, preQ2, _tolerance=epsilonP, _maxIt=nbIteration);
432 std::cout << cpuTime("QMR with preconditioner LDLt ") << std::endl;
433 printResultsOfCV(A, S, rvQmrP2, UN, epsilonP, out);
434 clear(rvQmrP2);
435
436 ACSS=ACSStmp;
437 std::cout << "* QMR with incomplete Preconditionner (iLDLT)" << eol;
438 PreconditionerTerm preQ3(ACSS, _ildltPrec);
439 TermVector rvQmrP3 = qmrSolve(ACSS, S, rvX0, preQ3, _tolerance=epsilonP, _maxIt=nbIteration);
440 std::cout << cpuTime("QMR SYM with preconditioner iLDLt ") << std::endl;
441 printResultsOfCV(ACSS, S, rvQmrP3, UN, epsilonP, out);
442 clear(rvQmrP3);
443
444 A=Atmp;
445 std::cout << "* QMR with Preconditionner (iLDLT)" << eol;
446 PreconditionerTerm preQ4(A, _ildltPrec);
447 TermVector rvQmrP4 = qmrSolve(A, S, rvX0, preQ4, _tolerance=epsilonP, _maxIt=nbIteration);
448 std::cout << cpuTime("QMR with preconditioner iLDLt ") << std::endl;
449 printResultsOfCV(A, S, rvQmrP4, UN, epsilonP, out);
450 clear(rvQmrP4);
451
452 ACSS=ACSStmp;
453 std::cout << "* QMR with incomplete Preconditionner (iLLt)" << eol;
454 PreconditionerTerm preQ5(ACSS, _illtPrec);
455 TermVector rvQmrP5 = qmrSolve(ACSS, S, rvX0, preQ5, _tolerance=epsilonP, _maxIt=nbIteration);
456 std::cout << cpuTime("QMR with preconditioner iLLt ") << std::endl;
457 printResultsOfCV(ACSS, S, rvQmrP5, UN, epsilonP, out);
458 clear(rvQmrP5);
459 //
460 // no transpose solver
461 // ACS=ACStmp;
462 // out << "* QMR with incomplete Preconditionner (iLU)" << eol;
463 // std::cout << "* QMR with incomplete Preconditionner (iLU)" << eol;
464 // PreconditionerTerm preQ4(ACS, _iluPrec);
465 // TermVector rvQmrP4 = qmrSolve(ACS, S, rvX1, preQ4, _tolerance=epsilonP, _maxIt=nbIteration);
466 // printResultsOfCV(ACS, S, rvQmrP4, UN, epsilonP, out);
467 // clear(rvQmrP4);
468
469
470 }
471
472
dev_Eric1(bool check)473 String dev_Eric1(bool check)
474 {
475 String rootname = "dev_Eric1";
476 trace_p->push(rootname);
477 std::stringstream out;
478 out.precision(testPrec);
479 verboseLevel(30);
480 numberOfThreads(1);
481 bool withTriangles = true;
482 bool withQuadrangles = true;
483 bool withTetrahedra = true;
484 bool withHexahedra = true;
485 verboseLevel(9);;
486
487 // std::stringstream out;
488 // out.precision(testPrec);
489 std::string meshInfo;
490 verboseLevel(2);
491 Number nbDivMin , nbDivMax, dd=1;;
492 Number ordMin, ordMax;
493 std::vector<String> sidenames(4, "");
494 Real h;
495 std::string nmFich;
496 Number NbLignes;
497 Number NbSolvers = 2;
498
499
500 std::cout << "==========================================" << eol;
501 std::cout << "= TEST 2D =" << eol;
502 std::cout << "==========================================" << eol<<eol;;
503 if (withTriangles) {
504 //create a mesh and Domains TRIANGLES
505 std::cout << "# -> Triangles P" <<1<< eol;
506 std::stringstream NmEps1;
507 meshInfo= " Triangle mesh of [0,1]x[0,1] with nbStep = ";
508 nbDivMin =5, nbDivMax=5;
509 ordMin = 1, ordMax=2;
510 NbLignes= (nbDivMax - nbDivMin + 1) * ( ordMax - ordMin +1) * (NbSolvers+2)+1;
511 std::cout<<"* "<<"triangles "<<NbLignes<<" "<<NbSolvers<<eol;
512 for(Number div = nbDivMin; div <= nbDivMax; div+=1)
513 {
514 meshInfo += tostring(div);
515 Real h=1./(dd*div);
516 Number nn=std::pow(2,div)+1;
517 Mesh mesh2d(Square(_origin=Point(0.,0.), _length=1, _nnodes=nn), _triangle, 1, _structured, "P1 mesh of [0,1]x[1,2]");
518 // string fn("xlifepp2D.geo");
519 // Mesh mesh2d=Mesh(inputsPathTo(fn),"Rectangle", _geo);
520 Domain omega = mesh2d.domain("Omega");
521 std::cout <<eol<<eol<< "====================================" << eol;
522 std::cout << "# Triangle mesh: " << div << " sudivisions, h = " <<h <<eol;
523 std::cout << "===========================================" << eol;
524 for (Number ord=ordMin; ord<=ordMax; ord++)
525 {
526 std::cout<<"--------------------------" <<eol;
527 std::cout << " -> interpolation P" <<ord<< eol;
528 std::cout<<"--------------------------" <<eol;
529 std::cout <<"# -> Ordre P" <<std::endl<<ord<<std::endl;
530 nmFich = NmEps1.str();
531 std::cout<<mesh2d.nbOfNodes()<<eol;
532 std::cout<<"* triangle"<<std::endl;
533 ISEpsilon_Pk(omega, ord, h, out);
534 }
535 }
536 }
537
538 if (withQuadrangles) {
539 //create a mesh and Domains Quadrangles
540 meshInfo= " quadrangle mesh of [0,1]x[0,1] with nbStep = ";
541 nbDivMin =5, nbDivMax=5;
542 ordMin = 1, ordMax=2;
543 std::stringstream NmEps2;
544 NmEps2 << "Epsi";
545 NbLignes= (nbDivMax - nbDivMin + 1) * ( ordMax - ordMin +1) * (NbSolvers+2)+1;
546 for(Number div = nbDivMin; div <= nbDivMax; div++)
547 {
548 meshInfo += tostring(div);
549 Real h=1./(dd*div);
550 Number nn=std::pow(2,div)+1;
551 Mesh mesh2d(Square(_origin=Point(0.,0.), _length=1, _nnodes=nn), _quadrangle, 1, _structured, "Q1 mesh of [0,1]x[0,1]");
552 Domain omega = mesh2d.domain("Omega");
553
554 std::cout <<std::endl;out <<eol<<eol<< "================" << eol;
555 std::cout <<std::endl;out << "# Quadrangle mesh: " << div << " sudivisions, h = " <<h <<eol;
556 std::cout <<std::endl;out << "===================================" << eol;
557 for(Number ord=ordMin; ord<=ordMax; ord++)
558 {
559 std::cout<<"--------------------------" <<eol;
560 std::cout << " -> interpolation P" <<ord<< eol;
561 std::cout<<"--------------------------" <<eol;
562 std::cout<<"* quadrangle"<<std::endl;
563 // ISEpsilon_Pk(omega, ord, h,NmEps2, out);
564 }
565 }
566 }
567
568 // create a mesh and Domains TETRAHEDRON
569 if (withTetrahedra)
570 {
571 std::cout << eol<<"========================================== " << eol;
572 std::cout << "= TEST 3D = " << eol;
573 std::cout << "========================================== " << eol;
574 meshInfo=" TETRAHEDRONS mesh of cube ";
575 nbDivMin =2, nbDivMax=2;
576 ordMin = 1, ordMax=2;
577
578 std::stringstream NmEps3;
579 NmEps3 << "Epsi";
580 NbLignes= (nbDivMax - nbDivMin + 1) * ( ordMax - ordMin +1) * (NbSolvers+2)+1;
581 for(Number div = nbDivMin; div <= nbDivMax; div++)
582 {
583 Number nn=std::pow(2,div)+1;
584 Mesh mesh3d(Cube(_origin=Point(0.,0.,0.), _length=1., _nnodes=nn), _tetrahedron, 1, _subdiv, meshInfo);
585 Domain omega = mesh3d.domain("Omega");
586 h = 2. / std::pow(2,div);
587 meshInfo += tostring(div);
588 std::cout <<eol<<eol<< "=========================================" << eol;
589 std::cout << "# TETRAHEDRONS mesh : " << div << " sudivisions, h = " <<h <<eol;
590 std::cout << "=========================================" << eol;
591 for(Number ord = ordMin; ord <= ordMax; ord++)
592 {
593 std::cout<<"--------------------------" <<eol;
594 std::cout<< " -> interpolation P" <<ord<< eol;
595 std::cout<<"--------------------------" <<eol;
596 std::cout<<"* tetrahedre"<<std::endl;
597 // ISEpsilon_Pk(omega, ord, h, NmEps3, out);
598 }
599 }
600 }
601
602 // create a mesh and Domains HEXAHEDRON
603 if (withHexahedra)
604 {
605 meshInfo=" HEXAHEDRONS mesh of cube ";
606 nbDivMin =2, nbDivMax=2;
607 ordMin = 1, ordMax=2;
608
609 std::stringstream NmEps4;
610 NmEps4 << "Epsi";
611 NbLignes= (nbDivMax - nbDivMin + 1) * ( ordMax - ordMin +1) * (NbSolvers+2)+1;
612 for(Number div = nbDivMin; div <= nbDivMax; div++)
613 {
614 Number nn=std::pow(2,div)+1;
615 Mesh mesh3d(Cube(_origin=Point(0.,0.,0.), _length=1., _nnodes=nn), _hexahedron, 1, _structured, meshInfo);
616 Domain omega = mesh3d.domain("Omega");
617 h = 2. / std::pow(2,div);
618 meshInfo += tostring(div);
619 std::cout <<eol<<eol<< "=========================================" << eol;
620 std::cout << "# HEXAHEDRONS mesh : " << div << " sudivisions, h = " <<h <<eol;
621 std::cout << "=========================================" << eol;
622 for(Number ord = ordMin; ord <= ordMax; ord++)
623 {
624 std::cout<<"--------------------------" <<eol;
625 std::cout << " -> interpolation P" <<ord<< eol;
626 std::cout<<"--------------------------" <<eol;
627 std::cout<<"* hexahedre"<<std::endl;
628 // ISEpsilon_Pk(omega, ord, h, NmEps4, out);
629 }
630 }
631 }
632
633
634
635 //------------------------------------------------------------------------------------
636 // save results in a file or compare results with some references value in a file
637 //------------------------------------------------------------------------------------
638 trace_p->pop();
639 if(check) { return diffResFile(out, rootname); }
640 else { return saveResToFile(out, rootname); }
641 }
642
643 }
644