Arcane  v3.14.10.0
Documentation développeur
Chargement...
Recherche...
Aucune correspondance
AlephPETSc.cc
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2024 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7/*---------------------------------------------------------------------------*/
8/* AlephPETSc.cc (C) 2000-2024 */
9/* */
10/* Implémentation PETSc de Aleph. */
11/*---------------------------------------------------------------------------*/
12/*---------------------------------------------------------------------------*/
13
14#include "arcane/aleph/AlephArcane.h"
15
16#include "arcane/utils/StringList.h"
17
18#define MPI_COMM_SUB *(MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator())
19
20#include "petscksp.h"
21#include "petscsys.h"
22
23// TODO: mieux gérer les sous-versions de PETSC
24#if PETSC_VERSION_GE(3,6,1)
25#include "petscmat.h"
26#elif PETSC_VERSION_(3,3,0)
27#include "petscpcmg.h"
28#elif PETSC_VERSION_(3,0,0)
29#include "petscmg.h"
30#else
31#error "Bad version of PETSc: should be 3.0.0, 3.3.0 or greater than 3.6.1."
32#endif
33
34/*---------------------------------------------------------------------------*/
35/*---------------------------------------------------------------------------*/
36
37namespace Arcane
38{
39
40namespace
41{
42// To convert from AlephInt data to PetscInt data
43PetscInt* _toPetscInt(AlephInt* v)
44{
45 static_assert(sizeof(PetscInt)==sizeof(AlephInt),"AlephInt and PetscInt should have the same size");
46 return reinterpret_cast<PetscInt*>(v);
47}
48// To convert from AlephInt data to PetscInt data
49const PetscInt* _toPetscInt(const AlephInt* v)
50{
51 static_assert(sizeof(PetscInt)==sizeof(AlephInt),"AlephInt and PetscInt should have the same size");
52 return reinterpret_cast<const PetscInt*>(v);
53}
54
55inline void
56petscCheck(PetscErrorCode error_code)
57{
58 // Note GG: Je ne sais pas vraiment à partir de quelle version de PETSc
59 // la valeur 'PETSC_SUCCESS' la fonction 'PetscCallVoid' sont
60 // disponibles mais elles le sont dans la 3.18.0.
61#if PETSC_VERSION_GE(3,19,0)
62 if (error_code==PETSC_SUCCESS)
63 return;
64 PetscCallVoid(error_code);
65 ARCANE_FATAL("Error in PETSc backend errcode={0}",error_code);
66#else
67 if (error_code!=0)
68 ARCANE_FATAL("Error in PETSc backend errcode={0}",error_code);
69#endif
70}
71
72}
73
74/*---------------------------------------------------------------------------*/
75/*---------------------------------------------------------------------------*/
76
78: public IAlephTopology
79{
80 public:
81
82 AlephTopologyPETSc(ITraceMng* tm, AlephKernel* kernel, Integer index, Integer nb_row_size);
83 ~AlephTopologyPETSc() override
84 {
85 debug() << "\33[1;5;32m\t\t\t[~AlephTopologyPETSc]\33[0m";
86 }
87
88 public:
89
90 void backupAndInitialize() override {}
91 void restore() override {}
92
93 private:
94
95 void _checkInitPETSc();
96};
97
98/*---------------------------------------------------------------------------*/
99/*---------------------------------------------------------------------------*/
100
102: public IAlephVector
103{
104 public:
105
107 void AlephVectorCreate(void) override;
108 void AlephVectorSet(const double*, const AlephInt*, Integer) override;
109 int AlephVectorAssemble(void) override;
110 void AlephVectorGet(double*, const AlephInt*, Integer) override;
111 void writeToFile(const String) override;
112 Real LinftyNorm(void);
113
114 public:
115
116 Vec m_petsc_vector;
117 PetscInt jSize, jUpper, jLower;
118};
119
120/*---------------------------------------------------------------------------*/
121/*---------------------------------------------------------------------------*/
122
124: public IAlephMatrix
125{
126 public:
127
129 void AlephMatrixCreate(void) override;
130 void AlephMatrixSetFilled(bool) override;
131 int AlephMatrixAssemble(void) override;
132 void AlephMatrixFill(int, AlephInt*, AlephInt*, double*) override;
133 Real LinftyNormVectorProductAndSub(AlephVector*, AlephVector*);
134 bool isAlreadySolved(AlephVectorPETSc*, AlephVectorPETSc*,
135 AlephVectorPETSc*, Real*, AlephParams*);
136 int AlephMatrixSolve(AlephVector*, AlephVector*,
137 AlephVector*, Integer&, Real*, AlephParams*) override;
138 void writeToFile(const String) override;
139
140 private:
141
142 Mat m_petsc_matrix;
143 KSP m_ksp_solver;
144};
145
146/*---------------------------------------------------------------------------*/
147/*---------------------------------------------------------------------------*/
148
150: public AbstractService
151, public IAlephFactoryImpl
152{
153 public:
154
157 {}
158 ~PETScAlephFactoryImpl() override
159 {
160 for (auto* v : m_IAlephVectors)
161 delete v;
162 for (auto* v : m_IAlephMatrixs)
163 delete v;
164 for (auto* v : m_IAlephTopologys)
165 delete v;
166 }
167
168 public:
169
170 void initialize() override {}
171 IAlephTopology* createTopology(ITraceMng* tm,
172 AlephKernel* kernel,
173 Integer index,
174 Integer nb_row_size) override
175 {
176 IAlephTopology* new_topology = new AlephTopologyPETSc(tm, kernel, index, nb_row_size);
177 m_IAlephTopologys.add(new_topology);
178 return new_topology;
179 }
180 IAlephVector* createVector(ITraceMng* tm,
181 AlephKernel* kernel,
182 Integer index) override
183 {
184 IAlephVector* new_vector = new AlephVectorPETSc(tm, kernel, index);
185 m_IAlephVectors.add(new_vector);
186 return new_vector;
187 }
188
189 IAlephMatrix* createMatrix(ITraceMng* tm,
190 AlephKernel* kernel,
191 Integer index) override
192 {
193 IAlephMatrix* new_matrix = new AlephMatrixPETSc(tm, kernel, index);
194 m_IAlephMatrixs.add(new_matrix);
195 return new_matrix;
196 }
197
198 private:
199
200 UniqueArray<IAlephVector*> m_IAlephVectors;
201 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
202 UniqueArray<IAlephTopology*> m_IAlephTopologys;
203};
204
205/*---------------------------------------------------------------------------*/
206/*---------------------------------------------------------------------------*/
207
208/*---------------------------------------------------------------------------*/
209/*---------------------------------------------------------------------------*/
210
211AlephTopologyPETSc::
212AlephTopologyPETSc(ITraceMng* tm,
213 AlephKernel* kernel,
214 Integer index,
215 Integer nb_row_size)
216: IAlephTopology(tm, kernel, index, nb_row_size)
217{
218 ARCANE_CHECK_POINTER(kernel);
219 if (!m_participating_in_solver) {
220 debug() << "\33[1;32m\t[AlephTopologyPETSc] Not concerned with this solver, returning\33[0m";
221 return;
222 }
223 debug() << "\33[1;32m\t\t[AlephTopologyPETSc] @" << this << "\33[0m";
224 if (!m_kernel->isParallel()) {
225 PETSC_COMM_WORLD = PETSC_COMM_SELF;
226 }
227 else {
228 PETSC_COMM_WORLD = *(MPI_Comm*)(kernel->subParallelMng(index)->getMPICommunicator());
229 }
230
231 _checkInitPETSc();
232}
233
234/*---------------------------------------------------------------------------*/
235/*---------------------------------------------------------------------------*/
236
237void AlephTopologyPETSc::
238_checkInitPETSc()
239{
240 // Ne fait rien si PETSc a déjà été initialisé.
241 PetscBool is_petsc_initialized = PETSC_FALSE;
242 PetscInitialized(&is_petsc_initialized);
243 if (is_petsc_initialized==PETSC_TRUE)
244 return;
245
246 AlephKernelSolverInitializeArguments& init_args = m_kernel->solverInitializeArgs();
247 bool has_args = init_args.hasValues();
248 debug() << Trace::Color::cyan() << "[AlephTopologyPETSc] Initializing PETSc. UseArg=" << has_args;
249 if (has_args){
250 const CommandLineArguments& args = init_args.commandLineArguments();
251 PetscInitialize(args.commandLineArgc(),args.commandLineArgv(),nullptr,nullptr);
252 }
253 else{
254 PetscInitializeNoArguments();
255 }
256}
257
258/*---------------------------------------------------------------------------*/
259/*---------------------------------------------------------------------------*/
260
261// ****************************************************************************
262// * AlephVectorPETSc
263// ****************************************************************************
264AlephVectorPETSc::
265AlephVectorPETSc(ITraceMng *tm,AlephKernel *kernel,Integer index)
266: IAlephVector(tm,kernel,index)
267, m_petsc_vector()
268, jSize(0)
269, jUpper(0)
270, jLower(-1)
271{
272 debug()<<"\t\t[AlephVectorPETSc::AlephVectorPETSc] new SolverVectorPETSc";
273}
274
275// ****************************************************************************
276// * AlephVectorCreate
277// ****************************************************************************
278void AlephVectorPETSc::
279AlephVectorCreate(void)
280{
281 for( int iCpu=0;iCpu<m_kernel->size();++iCpu){
282 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
283 continue;
284 if (jLower==-1)
285 jLower=m_kernel->topology()->gathered_nb_row(iCpu);
286 jUpper=m_kernel->topology()->gathered_nb_row(iCpu+1);
287 }
288 // Mise à jour de la taille locale du buffer pour le calcul plus tard de la norme max, par exemple
289 jSize=jUpper-jLower;
290 VecCreateMPI(MPI_COMM_SUB,
291 jSize, // n
292 m_kernel->topology()->nb_row_size(), // N
293 &m_petsc_vector);
294 debug()<<"\t\t[AlephVectorPETSc::AlephVectorCreate] PETSC VectorCreate"
295 <<" of local size=" <<jSize<<"/"<<m_kernel->topology()->nb_row_size()
296 <<", done";
297}
298
299// ****************************************************************************
300// * AlephVectorSet
301// ****************************************************************************
302void AlephVectorPETSc::
303AlephVectorSet(const double *bfr_val, const AlephInt *bfr_idx, Integer size)
304{
305 VecSetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val,INSERT_VALUES);
306 debug()<<"\t\t[AlephVectorPETSc::AlephVectorSet] "<<size<<" values inserted!";
307 // En séquentiel, il faut le fair aussi avec PETSc
308 AlephVectorAssemble();
309}
310
311
312// ****************************************************************************
313// * AlephVectorAssemble
314// ****************************************************************************
315int AlephVectorPETSc::
316AlephVectorAssemble(void)
317{
318 VecAssemblyBegin(m_petsc_vector);
319 VecAssemblyEnd(m_petsc_vector);
320 debug()<<"\t\t[AlephVectorPETSc::AlephVectorAssemble]";
321 return 0;
322}
323
324
325// ****************************************************************************
326// * AlephVectorGet
327// ****************************************************************************
328void AlephVectorPETSc::
329AlephVectorGet(double* bfr_val, const AlephInt* bfr_idx, Integer size)
330{
331 VecGetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val);
332 debug()<<"\t\t[AlephVectorPETSc::AlephVectorGet] fetched "<< size <<" values!";
333}
334
335
336// ****************************************************************************
337// * writeToFile
338// ****************************************************************************
339void AlephVectorPETSc::
340writeToFile(const String filename)
341{
342 ARCANE_UNUSED(filename);
343 debug()<<"\t\t[AlephVectorPETSc::writeToFile]";
344}
345
346// ****************************************************************************
347// * LinftyNorm
348// ****************************************************************************
349Real AlephVectorPETSc::
350LinftyNorm(void)
351{
352 PetscReal val;
353 VecNorm(m_petsc_vector,NORM_INFINITY,&val);
354 return val;
355}
356
357// ****************************************************************************
358// * AlephMatrixPETSc
359// ****************************************************************************
360AlephMatrixPETSc::
361AlephMatrixPETSc(ITraceMng *tm,
362 AlephKernel *kernel,
363 Integer index)
364: IAlephMatrix(tm,kernel,index)
365, m_petsc_matrix()
366, m_ksp_solver()
367{
368 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixPETSc] new AlephMatrixPETSc";
369}
370
371
372// ****************************************************************************
373// * AlephMatrixCreate
374// ****************************************************************************
375void AlephMatrixPETSc::
376AlephMatrixCreate(void)
377{
378 AlephInt ilower=-1;
379 AlephInt iupper=0;
380
381 for( int iCpu=0;iCpu<m_kernel->size();++iCpu ){
382 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu]) continue;
383 if (ilower==-1) ilower=m_kernel->topology()->gathered_nb_row(iCpu);
384 iupper=m_kernel->topology()->gathered_nb_row(iCpu+1);
385 }
386 AlephInt size = iupper-ilower;
387 AlephInt jlower=ilower;
388 AlephInt jupper=iupper;
389
390#if PETSC_VERSION_GE(3,3,0)
391 #define PETSC_VERSION_MatCreate MatCreateAIJ
392#elif PETSC_VERSION_(3,0,0)
393 #define PETSC_VERSION_MatCreate MatCreateMPIAIJ
394#endif
395 PetscInt* row_data = _toPetscInt(m_kernel->topology()->gathered_nb_row_elements().subView(ilower,size).data());
396 PETSC_VERSION_MatCreate(MPI_COMM_SUB,
397 iupper-ilower, // m = number of rows
398 jupper-jlower, // n = number of columns
399 m_kernel->topology()->nb_row_size(), // M = number of global rows
400 m_kernel->topology()->nb_row_size(), // N = number of global columns
401 0, // ignored (number of nonzeros per row in DIAGONAL portion of local submatrix)
402 // array containing the number of nonzeros in the various rows of the DIAGONAL portion of local submatrix
403 row_data,
404 0, // ignored (number of nonzeros per row in the OFF-DIAGONAL portion of local submatrix)
405 // array containing the number of nonzeros in the various rows of the OFF-DIAGONAL portion of local submatrix
406 row_data,
407 &m_petsc_matrix);
408 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
409 MatSetUp(m_petsc_matrix);
410 debug()<<"\t\t[AlephMatrixPetsc::AlephMatrixCreate] PETSC MatrixCreate idx:"<<m_index
411 <<", ("<<ilower<<"->"<<(iupper-1)<<")";
412}
413
414
415// ****************************************************************************
416// * AlephMatrixSetFilled
417// ****************************************************************************
418void AlephMatrixPETSc::
419AlephMatrixSetFilled(bool)
420{
421 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixSetFilled] done";
422 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
423}
424
425
426// ****************************************************************************
427// * AlephMatrixAssemble
428// ****************************************************************************
429int AlephMatrixPETSc::
430AlephMatrixAssemble(void)
431{
432 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixAssemble] AlephMatrixAssemble";
433 MatAssemblyBegin(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
434 MatAssemblyEnd(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
435 return 0;
436}
437
438
439// ****************************************************************************
440// * AlephMatrixFill
441// ****************************************************************************
442void AlephMatrixPETSc::
443AlephMatrixFill(int size, AlephInt* rows, AlephInt* cols, double *values)
444{
445 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] size="<<size;
446 for(int i=0;i<size;i++){
447 //debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] i="<<i;
448 MatSetValue(m_petsc_matrix, rows[i], cols[i], values[i], INSERT_VALUES);
449 }
450 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] done";
451 // PETSc réclame systématiquement l'assemblage
452 AlephMatrixAssemble();
453}
454
455
456
457// ****************************************************************************
458// * LinftyNormVectorProductAndSub
459// ****************************************************************************
460Real AlephMatrixPETSc::
461LinftyNormVectorProductAndSub(AlephVector* x,AlephVector* b)
462{
463 ARCANE_UNUSED(x);
464 ARCANE_UNUSED(b);
465 throw FatalErrorException("LinftyNormVectorProductAndSub", "error");
466}
467
468
469// ****************************************************************************
470// * isAlreadySolved
471// ****************************************************************************
472bool AlephMatrixPETSc::
473isAlreadySolved(AlephVectorPETSc* x,
474 AlephVectorPETSc* b,
475 AlephVectorPETSc* tmp,
476 Real* residual_norm,
477 AlephParams* params)
478{
479 ARCANE_UNUSED(x);
480 ARCANE_UNUSED(b);
481 ARCANE_UNUSED(tmp);
482 ARCANE_UNUSED(residual_norm);
483 ARCANE_UNUSED(params);
484 return false;
485}
486
487
488// ****************************************************************************
489// * AlephMatrixSolve
490// ****************************************************************************
491int AlephMatrixPETSc::
492AlephMatrixSolve(AlephVector* x,AlephVector* b,AlephVector* t,
493 Integer& nb_iteration,Real* residual_norm,
494 AlephParams* solver_param)
495{
496 ARCANE_UNUSED(t);
497
498 PC prec;
499 PetscInt its;
500 PetscReal norm;
501 KSPConvergedReason reason;
502
503 AlephVectorPETSc* x_petsc = dynamic_cast<AlephVectorPETSc*> (x->implementation());
504 AlephVectorPETSc* b_petsc = dynamic_cast<AlephVectorPETSc*> (b->implementation());
505
506 ARCANE_CHECK_POINTER(x_petsc);
507 ARCANE_CHECK_POINTER(b_petsc);
508 const bool is_parallel = m_kernel->subParallelMng(m_index)->isParallel();
509
510 Vec solution = x_petsc->m_petsc_vector;
511 Vec RHS = b_petsc->m_petsc_vector;
512
513 debug()<<"[AlephMatrixSolve]";
514 KSPCreate(MPI_COMM_SUB,&m_ksp_solver);
515#if PETSC_VERSION_GE(3,6,1)
516 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix);//SAME_PRECONDITIONER);
517#else
518 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix,SAME_NONZERO_PATTERN);//SAME_PRECONDITIONER);
519#endif
520 // Builds KSP for a particular solver
521 switch(solver_param->method()){
522 // Preconditioned conjugate gradient (PCG) iterative method
523 case TypesSolver::PCG : KSPSetType(m_ksp_solver,KSPCG); break;
524 // BiCGStab (Stabilized version of BiConjugate Gradient Squared) method
525 case TypesSolver::BiCGStab : KSPSetType(m_ksp_solver,KSPBCGS); break;
526 // IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared) method
527 // in an alternative form to have only a single global reduction operation instead of the usual 3 (or 4)
528 case TypesSolver::BiCGStab2: KSPSetType(m_ksp_solver,KSPIBCGS); break;
529 // Generalized Minimal Residual method. (Saad and Schultz, 1986) with restart
530 case TypesSolver::GMRES : KSPSetType(m_ksp_solver,KSPGMRES); break;
531 default : throw ArgumentException("AlephMatrixPETSc::AlephMatrixSolve", "Unknown solver method");
532 }
533
534 KSPGetPC(m_ksp_solver,&prec);
535
536 switch(solver_param->precond()){
537 case TypesSolver::NONE : PCSetType(prec,PCNONE); break;
538 // Jacobi (i.e. diagonal scaling preconditioning)
539 case TypesSolver::DIAGONAL : PCSetType(prec,PCJACOBI); break;
540 // Incomplete factorization preconditioners.
541 case TypesSolver::ILU:
542 // ILU ne fonctionne pas en parallèle avec PETSc (sauf si compilé avec Hypre)
543 // ILU can be set up in parallel via a sub-preconditioner to Block-Jacobi preconditioner
544 if (is_parallel)
545 PCSetType(prec,PCBJACOBI);
546 else
547 PCSetType(prec,PCILU);
548 break;
549 // Incomplete Cholesky factorization preconditioners.
550 case TypesSolver::IC:
551 // IC can be set up in parallel via a sub-preconditioner to Block-Jacobi preconditioner
552 if (is_parallel) {
553 PCSetType(prec,PCBJACOBI); // set PC to Block Jacobi
554 PetscInt nlocal, first;
555 KSP *subksp; // KSP context for subdomain
556 PC subpc; // PC context for subdomain
557 KSPSetUp(m_ksp_solver);
558 PCBJacobiGetSubKSP(prec, &nlocal, &first, &subksp); // Extract the KSP context array for the local blocks
559 for (int i = 0; i < nlocal; i++) {
560 KSPGetPC(subksp[i], &subpc);
561 PCSetType(subpc, PCICC); // set subdomain PC to IC
562 }
563 }
564 else
565 PCSetType(prec,PCICC);
566 break;
567 // Sparse Approximate Inverse method of Grote and Barnard as a preconditioner (SIAM J. Sci. Comput.; vol 18, nr 3)
568 case TypesSolver::SPAIstat : PCSetType(prec,PCSPAI); break;
569 // Use multigrid preconditioning.
570 // This preconditioner requires you provide additional information
571 // about the coarser grid matrices and restriction/interpolation operators.
572 case TypesSolver::AMG:
573 // By default PCMG uses GMRES on the fine grid smoother so this should be used with KSPFGMRES
574 // or the smoother changed to not use GMRES
575 // Implements the Flexible Generalized Minimal Residual method. developed by Saad with restart
576 KSPSetType(m_ksp_solver,KSPFGMRES);
577 PCSetType(prec,PCMG);
578 //Sets the number of levels to use with MG. Must be called before any other MG routine.
579 PCMGSetLevels(prec,1,
580 (MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator()));
581 // Determines the form of multigrid to use: multiplicative, additive, full, or the Kaskade algorithm.
582 PCMGSetType(prec,PC_MG_MULTIPLICATIVE); // PC_MG_MULTIPLICATIVE,PC_MG_ADDITIVE,PC_MG_FULL,PC_MG_KASKADE
583 // Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more complicated cycling.
584 PCMGSetCycleType(prec,PC_MG_CYCLE_V); // PC_MG_CYCLE_V or PC_MG_CYCLE_W
585 //Sets the number of pre-smoothing steps to use on all levels
586#if PETSC_VERSION_GE(3,10,2)
587 PCMGSetNumberSmooth(prec, 1);
588#else
589 PCMGSetNumberSmoothDown(prec, 1);
590 PCMGSetNumberSmoothUp(prec, 1);
591#endif
592 // PCMGSetResidual: Sets the function to be used to calculate the residual on the lth level
593 // PCMGSetInterpolation: Sets the function to be used to calculate the interpolation from l-1 to the lth level
594 // PCMGSetRestriction: Sets the function to be used to restrict vector from level l to l-1
595 // PCMGSetRhs: Sets the vector space to be used to store the right-hand side on a particular level
596 // PCMGSetX: Sets the vector space to be used to store the solution on a particular level
597 // PCMGSetR: Sets the vector space to be used to store the residual on a particular level
598 //PCMGSetRestriction, PCMGSetRhs, PCMGSetX, PCMGSetR, etc.
599 //PCMGSetLevels
600 break;
601 case TypesSolver::AINV : throw ArgumentException("AlephMatrixPETSc", "preconditionnement AINV indisponible");
602 case TypesSolver::SPAIdyn : throw ArgumentException("AlephMatrixPETSc", "preconditionnement SPAIdyn indisponible");
603 case TypesSolver::ILUp : throw ArgumentException("AlephMatrixPETSc", "preconditionnement ILUp indisponible");
604 case TypesSolver::POLY : throw ArgumentException("AlephMatrixPETSc", "preconditionnement POLY indisponible");
605 default : throw ArgumentException("AlephMatrixPETSc", "preconditionnement inconnu");
606 }
607
608 if (solver_param->xoUser())
609 KSPSetInitialGuessNonzero(m_ksp_solver,PETSC_TRUE);
610
611 KSPSetTolerances(m_ksp_solver,
612 solver_param->epsilon(),
613 PETSC_DEFAULT,
614 PETSC_DEFAULT,
615 solver_param->maxIter());
616
617 // override KSP options from commandline if present//
618 KSPSetFromOptions(m_ksp_solver);
619
620 // override PC options from commandline if present//
621 PCSetFromOptions(prec) ;
622
623 // setup KSP solver and precondtioner //
624 KSPSetUp(m_ksp_solver);
625 PCSetUp(prec);
626
627 // print info which KSP solver and PC used //
628 KSPType kt;
629 KSPGetType(m_ksp_solver,&kt);
630 PCType pt;
631 PCGetType(prec,&pt);
632 info(4) << "[AlephMatrixSolve] Will use KSP type :" << kt << " PC type : "<< pt
633 << " is_parallel=" << is_parallel;
634
635 // solve the linear system //
636 debug()<<"[AlephMatrixSolve] All set up, now solving";
637 petscCheck(KSPSolve(m_ksp_solver,RHS,solution));
638 debug()<<"[AlephMatrixSolve] solved";
639 petscCheck(KSPGetConvergedReason(m_ksp_solver,&reason));
640
641 switch(reason){
642#if !PETSC_VERSION_(3,0,0)
643 case(KSP_CONVERGED_RTOL_NORMAL):{break;}
644 case(KSP_CONVERGED_ATOL_NORMAL):{break;}
645#endif
646 case(KSP_CONVERGED_RTOL):{break;}
647 case(KSP_CONVERGED_ATOL):{break;}
648 case(KSP_CONVERGED_ITS):{break;}
649#if !PETSC_VERSION_GE(3,19,0)
650 case(KSP_CONVERGED_CG_NEG_CURVE):{break;}
651 case(KSP_CONVERGED_CG_CONSTRAINED):{break;}
652#endif
653 case(KSP_CONVERGED_STEP_LENGTH):{break;}
654 case(KSP_CONVERGED_HAPPY_BREAKDOWN):{break;}
655 /* diverged */
656 case(KSP_DIVERGED_NULL):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NULL");}
657 case(KSP_DIVERGED_ITS):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_ITS");}
658 case(KSP_DIVERGED_DTOL):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_DTOL");}
659 case(KSP_DIVERGED_BREAKDOWN):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_BREAKDOWN");}
660 case(KSP_DIVERGED_BREAKDOWN_BICG):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_BREAKDOWN_BICG");}
661 case(KSP_DIVERGED_NONSYMMETRIC):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NONSYMMETRIC");}
662 case(KSP_DIVERGED_INDEFINITE_PC):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_INDEFINITE_PC");}
663#if PETSC_VERSION_GE(3,6,1)
664 case(KSP_DIVERGED_NANORINF):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NANORINF");}
665#else
666 case(KSP_DIVERGED_NAN):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NAN");}
667#endif
668 case(KSP_DIVERGED_INDEFINITE_MAT):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_INDEFINITE_MAT");}
669
670 case(KSP_CONVERGED_ITERATING):{break;}
671 default: throw Exception("AlephMatrixPETSc::Solve", "");
672 }
673
674 KSPGetIterationNumber(m_ksp_solver,&its);
675 nb_iteration=its;
676
677 KSPGetResidualNorm(m_ksp_solver,&norm);
678 *residual_norm=norm;
679
680 return 0;
681}
682
683
684// ****************************************************************************
685// * writeToFile
686// ****************************************************************************
687
688void AlephMatrixPETSc::
689writeToFile(const String filename)
690{
691 ARCANE_UNUSED(filename);
692}
693
694/*---------------------------------------------------------------------------*/
695/*---------------------------------------------------------------------------*/
696
697ARCANE_REGISTER_APPLICATION_FACTORY(PETScAlephFactoryImpl,IAlephFactoryImpl,PETScAlephFactory);
698
699/*---------------------------------------------------------------------------*/
700/*---------------------------------------------------------------------------*/
701
702}
703
704/*---------------------------------------------------------------------------*/
705/*---------------------------------------------------------------------------*/
#define ARCANE_CHECK_POINTER(ptr)
Macro retournant le pointeur ptr s'il est non nul ou lancant une exception s'il est nul.
#define ARCANE_FATAL(...)
Macro envoyant une exception FatalErrorException.
#define ARCANE_REGISTER_APPLICATION_FACTORY(aclass, ainterface, aname)
Enregistre un service de fabrique pour la classe aclass.
Classe de base d'un service.
Paramètres d'un système linéraire.
Definition AlephParams.h:34
Vecteur d'un système linéaire.
Definition AlephVector.h:33
Interface d'une fabrique d'implémentation pour Aleph.
virtual bool isParallel() const =0
Retourne true si l'exécution est parallèle.
virtual void * getMPICommunicator()=0
Adresse du communicateur MPI associé à ce gestionnaire.
Lecteur des fichiers de maillage via la bibliothèque LIMA.
Definition Lima.cc:120
Structure contenant les informations pour créer un service.
Interface du gestionnaire de traces.
Chaîne de caractères unicode.
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flot pour un message de debug.
TraceMessage info() const
Flot pour un message d'information.
-*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
int AlephInt
Type par défaut pour indexer les lignes et les colonnes des matrices et vecteurs.
Definition AlephGlobal.h:50