Arcane  v3.16.2.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-2025 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-2025 */
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
77class AlephTopologyPETSc
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
101class AlephVectorPETSc
102: public IAlephVector
103{
104 public:
105
106 AlephVectorPETSc(ITraceMng*, AlephKernel*, Integer);
107 ~AlephVectorPETSc();
108 void AlephVectorCreate(void) override;
109 void AlephVectorSet(const double*, const AlephInt*, Integer) override;
110 int AlephVectorAssemble(void) override;
111 void AlephVectorGet(double*, const AlephInt*, Integer) override;
112 void writeToFile(const String) override;
113 Real LinftyNorm(void);
114
115 public:
116
117 Vec m_petsc_vector = nullptr;
118 PetscInt jSize = 0;
119 PetscInt jUpper = 0;
120 PetscInt jLower = -1;
121};
122
123/*---------------------------------------------------------------------------*/
124/*---------------------------------------------------------------------------*/
125
126class AlephMatrixPETSc
127: public IAlephMatrix
128{
129 public:
130
131 AlephMatrixPETSc(ITraceMng*, AlephKernel*, Integer);
132 ~AlephMatrixPETSc();
133 void AlephMatrixCreate(void) override;
134 void AlephMatrixSetFilled(bool) override;
135 int AlephMatrixAssemble(void) override;
136 void AlephMatrixFill(int, AlephInt*, AlephInt*, double*) override;
137 Real LinftyNormVectorProductAndSub(AlephVector*, AlephVector*);
138 bool isAlreadySolved(AlephVectorPETSc*, AlephVectorPETSc*,
140 int AlephMatrixSolve(AlephVector*, AlephVector*,
141 AlephVector*, Integer&, Real*, AlephParams*) override;
142 void writeToFile(const String) override;
143
144 private:
145
146 Mat m_petsc_matrix = nullptr;
147 KSP m_ksp_solver = nullptr;
148};
149
150/*---------------------------------------------------------------------------*/
151/*---------------------------------------------------------------------------*/
152
153namespace
154{
155bool global_need_petsc_finalize = false;
156}
157
158/*---------------------------------------------------------------------------*/
159/*---------------------------------------------------------------------------*/
160
161class PETScAlephFactoryImpl
162: public AbstractService
163, public IAlephFactoryImpl
164{
165 public:
166 explicit PETScAlephFactoryImpl(const ServiceBuildInfo& sbi)
167 : AbstractService(sbi)
168 {}
169 ~PETScAlephFactoryImpl() override
170 {
171 for (auto* v : m_IAlephVectors)
172 delete v;
173 for (auto* v : m_IAlephMatrixs)
174 delete v;
175 for (auto* v : m_IAlephTopologys)
176 delete v;
177 if (global_need_petsc_finalize){
178 // TODO: Il faudrait appeler PETscFinalize() mais il
179 // faut le faire faire avant le MPI_Finalize() ce qui
180 // n'est pas toujours le cas suivant quand cette fabrique
181 // est détruite. En attendant, on n'appelle pas la fonction
182 //PetscFinalize();
183 global_need_petsc_finalize = false;
184 }
185 }
186
187 public:
188
189 void initialize() override {}
190 IAlephTopology* createTopology(ITraceMng* tm,
191 AlephKernel* kernel,
192 Integer index,
193 Integer nb_row_size) override
194 {
195 IAlephTopology* new_topology = new AlephTopologyPETSc(tm, kernel, index, nb_row_size);
196 m_IAlephTopologys.add(new_topology);
197 return new_topology;
198 }
199 IAlephVector* createVector(ITraceMng* tm,
200 AlephKernel* kernel,
201 Integer index) override
202 {
203 IAlephVector* new_vector = new AlephVectorPETSc(tm, kernel, index);
204 m_IAlephVectors.add(new_vector);
205 return new_vector;
206 }
207
208 IAlephMatrix* createMatrix(ITraceMng* tm,
209 AlephKernel* kernel,
210 Integer index) override
211 {
212 IAlephMatrix* new_matrix = new AlephMatrixPETSc(tm, kernel, index);
213 m_IAlephMatrixs.add(new_matrix);
214 return new_matrix;
215 }
216
217 private:
218
219 UniqueArray<IAlephVector*> m_IAlephVectors;
220 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
221 UniqueArray<IAlephTopology*> m_IAlephTopologys;
222};
223
224/*---------------------------------------------------------------------------*/
225/*---------------------------------------------------------------------------*/
226
227/*---------------------------------------------------------------------------*/
228/*---------------------------------------------------------------------------*/
229
230AlephTopologyPETSc::
231AlephTopologyPETSc(ITraceMng* tm, AlephKernel* kernel, Integer index, Integer nb_row_size)
232: IAlephTopology(tm, kernel, index, nb_row_size)
233{
234 ARCANE_CHECK_POINTER(kernel);
235 if (!m_participating_in_solver) {
236 debug() << "\33[1;32m\t[AlephTopologyPETSc] Not concerned with this solver, returning\33[0m";
237 return;
238 }
239 debug() << "\33[1;32m\t\t[AlephTopologyPETSc] @" << this << "\33[0m";
240 if (!m_kernel->isParallel()) {
241 PETSC_COMM_WORLD = PETSC_COMM_SELF;
242 }
243 else {
244 PETSC_COMM_WORLD = *(MPI_Comm*)(kernel->subParallelMng(index)->getMPICommunicator());
245 }
246
247 _checkInitPETSc();
248}
249
250/*---------------------------------------------------------------------------*/
251/*---------------------------------------------------------------------------*/
252
253void AlephTopologyPETSc::
254_checkInitPETSc()
255{
256 // Ne fait rien si PETSc a déjà été initialisé.
257 PetscBool is_petsc_initialized = PETSC_FALSE;
258 PetscInitialized(&is_petsc_initialized);
259 if (is_petsc_initialized==PETSC_TRUE)
260 return;
261
262 AlephKernelSolverInitializeArguments& init_args = m_kernel->solverInitializeArgs();
263 bool has_args = init_args.hasValues();
264 debug() << Trace::Color::cyan() << "[AlephTopologyPETSc] Initializing PETSc. UseArg=" << has_args;
265
266 if (has_args){
267 const CommandLineArguments& args = init_args.commandLineArguments();
268 PetscInitialize(args.commandLineArgc(),args.commandLineArgv(),nullptr,nullptr);
269 }
270 else{
271 PetscInitializeNoArguments();
272 }
273 global_need_petsc_finalize = true;
274}
275
276/*---------------------------------------------------------------------------*/
277/*---------------------------------------------------------------------------*/
278
279// ****************************************************************************
280// * AlephVectorPETSc
281// ****************************************************************************
282AlephVectorPETSc::
283AlephVectorPETSc(ITraceMng *tm,AlephKernel *kernel,Integer index)
284: IAlephVector(tm,kernel,index)
285, m_petsc_vector()
286, jSize(0)
287, jUpper(0)
288, jLower(-1)
289{
290 debug()<<"\t\t[AlephVectorPETSc::AlephVectorPETSc] new SolverVectorPETSc";
291}
292
293/*---------------------------------------------------------------------------*/
294/*---------------------------------------------------------------------------*/
295
296AlephVectorPETSc::
297~AlephVectorPETSc()
298{
299 VecDestroy(&m_petsc_vector);
300}
301
302/*---------------------------------------------------------------------------*/
303/*---------------------------------------------------------------------------*/
304
305// ****************************************************************************
306// * AlephVectorCreate
307// ****************************************************************************
308void AlephVectorPETSc::
309AlephVectorCreate(void)
310{
311 for( int iCpu=0;iCpu<m_kernel->size();++iCpu){
312 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
313 continue;
314 if (jLower==-1)
315 jLower=m_kernel->topology()->gathered_nb_row(iCpu);
316 jUpper=m_kernel->topology()->gathered_nb_row(iCpu+1);
317 }
318 // Mise à jour de la taille locale du buffer pour le calcul plus tard de la norme max, par exemple
319 jSize=jUpper-jLower;
320 VecCreateMPI(MPI_COMM_SUB,
321 jSize, // n
322 m_kernel->topology()->nb_row_size(), // N
323 &m_petsc_vector);
324 debug()<<"\t\t[AlephVectorPETSc::AlephVectorCreate] PETSC VectorCreate"
325 <<" of local size=" <<jSize<<"/"<<m_kernel->topology()->nb_row_size()
326 <<", done";
327}
328
329// ****************************************************************************
330// * AlephVectorSet
331// ****************************************************************************
332void AlephVectorPETSc::
333AlephVectorSet(const double *bfr_val, const AlephInt *bfr_idx, Integer size)
334{
335 VecSetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val,INSERT_VALUES);
336 debug()<<"\t\t[AlephVectorPETSc::AlephVectorSet] "<<size<<" values inserted!";
337 // En séquentiel, il faut le fair aussi avec PETSc
338 AlephVectorAssemble();
339}
340
341
342// ****************************************************************************
343// * AlephVectorAssemble
344// ****************************************************************************
345int AlephVectorPETSc::
346AlephVectorAssemble(void)
347{
348 VecAssemblyBegin(m_petsc_vector);
349 VecAssemblyEnd(m_petsc_vector);
350 debug()<<"\t\t[AlephVectorPETSc::AlephVectorAssemble]";
351 return 0;
352}
353
354
355// ****************************************************************************
356// * AlephVectorGet
357// ****************************************************************************
358void AlephVectorPETSc::
359AlephVectorGet(double* bfr_val, const AlephInt* bfr_idx, Integer size)
360{
361 VecGetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val);
362 debug()<<"\t\t[AlephVectorPETSc::AlephVectorGet] fetched "<< size <<" values!";
363}
364
365
366// ****************************************************************************
367// * writeToFile
368// ****************************************************************************
369void AlephVectorPETSc::
370writeToFile(const String filename)
371{
372 ARCANE_UNUSED(filename);
373 debug()<<"\t\t[AlephVectorPETSc::writeToFile]";
374}
375
376// ****************************************************************************
377// * LinftyNorm
378// ****************************************************************************
379Real AlephVectorPETSc::
380LinftyNorm(void)
381{
382 PetscReal val;
383 VecNorm(m_petsc_vector,NORM_INFINITY,&val);
384 return val;
385}
386
387// ****************************************************************************
388// * AlephMatrixPETSc
389// ****************************************************************************
390
391/*---------------------------------------------------------------------------*/
392/*---------------------------------------------------------------------------*/
393
394AlephMatrixPETSc::
395AlephMatrixPETSc(ITraceMng *tm, AlephKernel *kernel, Integer index)
396: IAlephMatrix(tm,kernel,index)
397{
398 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixPETSc] new AlephMatrixPETSc";
399}
400
401/*---------------------------------------------------------------------------*/
402/*---------------------------------------------------------------------------*/
403
404AlephMatrixPETSc::
405~AlephMatrixPETSc()
406{
407 if (m_petsc_matrix)
408 MatDestroy(&m_petsc_matrix);
409 if (m_ksp_solver)
410 KSPDestroy(&m_ksp_solver);
411}
412
413/*---------------------------------------------------------------------------*/
414/*---------------------------------------------------------------------------*/
415
416// ****************************************************************************
417// * AlephMatrixCreate
418// ****************************************************************************
419void AlephMatrixPETSc::
420AlephMatrixCreate(void)
421{
422 AlephInt ilower=-1;
423 AlephInt iupper=0;
424
425 for( int iCpu=0;iCpu<m_kernel->size();++iCpu ){
426 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu]) continue;
427 if (ilower==-1) ilower=m_kernel->topology()->gathered_nb_row(iCpu);
428 iupper=m_kernel->topology()->gathered_nb_row(iCpu+1);
429 }
430
431 AlephInt size = iupper-ilower;
432
433 PetscInt* row_data = _toPetscInt(m_kernel->topology()->gathered_nb_row_elements().subView(ilower,size).data());
434
435 PetscInt m = size; // Local number of rows
436 PetscInt n = size; // Local number of columns (assuming square)
437 PetscInt M = m_kernel->topology()->nb_row_size(); // Global number of rows
438 PetscInt N = M; // Global number of columns (assuming square)
439
440 MatCreate(PETSC_COMM_WORLD, &m_petsc_matrix);
441 MatSetSizes(m_petsc_matrix, m, n, M, N);
442
443 // Preallocate (AIJ is sparse, so preallocation improves performance)
444 MatSetType(m_petsc_matrix, MATAIJ);
445 MatMPIAIJSetPreallocation(m_petsc_matrix, 0, row_data, 0, row_data);
446
447 // Allow matrix type selection at runtime (Flexible, Runtime Configurable with -mat_type)
448 MatSetFromOptions(m_petsc_matrix);
449#if PETSC_VERSION_LE(3,0,0)
450 AlephInt jlower=ilower;
451 AlephInt jupper=iupper;
452
453 PETSC_VERSION_MatCreate(MPI_COMM_SUB,
454 iupper-ilower, // m = number of rows
455 jupper-jlower, // n = number of columns
456 m_kernel->topology()->nb_row_size(), // M = number of global rows
457 m_kernel->topology()->nb_row_size(), // N = number of global columns
458 0, // ignored (number of nonzeros per row in DIAGONAL portion of local submatrix)
459 // array containing the number of nonzeros in the various rows of the DIAGONAL portion of local submatrix
460 row_data,
461 0, // ignored (number of nonzeros per row in the OFF-DIAGONAL portion of local submatrix)
462 // array containing the number of nonzeros in the various rows of the OFF-DIAGONAL portion of local submatrix
463 row_data,
464 &m_petsc_matrix);
465#endif
466
467 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
468 MatSetUp(m_petsc_matrix);
469 debug()<<"\t\t[AlephMatrixPetsc::AlephMatrixCreate] PETSC MatrixCreate idx:"<<m_index
470 <<", ("<<ilower<<"->"<<(iupper-1)<<")";
471}
472
473
474// ****************************************************************************
475// * AlephMatrixSetFilled
476// ****************************************************************************
477void AlephMatrixPETSc::
478AlephMatrixSetFilled(bool)
479{
480 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixSetFilled] done";
481 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
482}
483
484
485// ****************************************************************************
486// * AlephMatrixAssemble
487// ****************************************************************************
488int AlephMatrixPETSc::
489AlephMatrixAssemble(void)
490{
491 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixAssemble] AlephMatrixAssemble";
492 MatAssemblyBegin(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
493 MatAssemblyEnd(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
494 return 0;
495}
496
497
498// ****************************************************************************
499// * AlephMatrixFill
500// ****************************************************************************
501void AlephMatrixPETSc::
502AlephMatrixFill(int size, AlephInt* rows, AlephInt* cols, double *values)
503{
504 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] size="<<size;
505 for(int i=0;i<size;i++){
506 //debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] i="<<i;
507 MatSetValue(m_petsc_matrix, rows[i], cols[i], values[i], INSERT_VALUES);
508 }
509 debug()<<"\t\t[AlephMatrixPETSc::AlephMatrixFill] done";
510 // PETSc réclame systématiquement l'assemblage
511 AlephMatrixAssemble();
512}
513
514
515
516// ****************************************************************************
517// * LinftyNormVectorProductAndSub
518// ****************************************************************************
519Real AlephMatrixPETSc::
520LinftyNormVectorProductAndSub(AlephVector* x,AlephVector* b)
521{
522 ARCANE_UNUSED(x);
523 ARCANE_UNUSED(b);
524 throw FatalErrorException("LinftyNormVectorProductAndSub", "error");
525}
526
527
528// ****************************************************************************
529// * isAlreadySolved
530// ****************************************************************************
531bool AlephMatrixPETSc::
532isAlreadySolved(AlephVectorPETSc* x,
534 AlephVectorPETSc* tmp,
535 Real* residual_norm,
536 AlephParams* params)
537{
538 ARCANE_UNUSED(x);
539 ARCANE_UNUSED(b);
540 ARCANE_UNUSED(tmp);
541 ARCANE_UNUSED(residual_norm);
542 ARCANE_UNUSED(params);
543 return false;
544}
545
546
547// ****************************************************************************
548// * AlephMatrixSolve
549// ****************************************************************************
550int AlephMatrixPETSc::
551AlephMatrixSolve(AlephVector* x,AlephVector* b,AlephVector* t,
552 Integer& nb_iteration,Real* residual_norm,
553 AlephParams* solver_param)
554{
555 ARCANE_UNUSED(t);
556
557 PC prec;
558 PetscInt its;
559 PetscReal norm;
560 KSPConvergedReason reason;
561
562 AlephVectorPETSc* x_petsc = dynamic_cast<AlephVectorPETSc*> (x->implementation());
563 AlephVectorPETSc* b_petsc = dynamic_cast<AlephVectorPETSc*> (b->implementation());
564
565 ARCANE_CHECK_POINTER(x_petsc);
566 ARCANE_CHECK_POINTER(b_petsc);
567 const bool is_parallel = m_kernel->subParallelMng(m_index)->isParallel();
568
569 Vec solution = x_petsc->m_petsc_vector;
570 Vec RHS = b_petsc->m_petsc_vector;
571
572 debug()<<"[AlephMatrixSolve]";
573 if (m_ksp_solver)
574 KSPDestroy(&m_ksp_solver);
575 KSPCreate(MPI_COMM_SUB,&m_ksp_solver);
576#if PETSC_VERSION_GE(3,6,1)
577 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix);//SAME_PRECONDITIONER);
578#else
579 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix,SAME_NONZERO_PATTERN);//SAME_PRECONDITIONER);
580#endif
581 // Builds KSP for a particular solver
582 switch(solver_param->method()){
583 // Preconditioned conjugate gradient (PCG) iterative method
584 case TypesSolver::PCG : KSPSetType(m_ksp_solver,KSPCG); break;
585 // BiCGStab (Stabilized version of BiConjugate Gradient Squared) method
586 case TypesSolver::BiCGStab : KSPSetType(m_ksp_solver,KSPBCGS); break;
587 // IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared) method
588 // in an alternative form to have only a single global reduction operation instead of the usual 3 (or 4)
589 case TypesSolver::BiCGStab2: KSPSetType(m_ksp_solver,KSPIBCGS); break;
590 // Generalized Minimal Residual method. (Saad and Schultz, 1986) with restart
591 case TypesSolver::GMRES : KSPSetType(m_ksp_solver,KSPGMRES); break;
592 default : throw ArgumentException("AlephMatrixPETSc::AlephMatrixSolve", "Unknown solver method");
593 }
594
595 KSPGetPC(m_ksp_solver,&prec);
596
597 switch(solver_param->precond()){
598 case TypesSolver::NONE : PCSetType(prec,PCNONE); break;
599 // Jacobi (i.e. diagonal scaling preconditioning)
600 case TypesSolver::DIAGONAL : PCSetType(prec,PCJACOBI); break;
601 // Incomplete factorization preconditioners.
602 case TypesSolver::ILU:
603 // ILU ne fonctionne pas en parallèle avec PETSc (sauf si compilé avec Hypre)
604 // ILU can be set up in parallel via a sub-preconditioner to Block-Jacobi preconditioner
605 if (is_parallel)
606 PCSetType(prec,PCBJACOBI);
607 else
608 PCSetType(prec,PCILU);
609 break;
610 // Incomplete Cholesky factorization preconditioners.
611 case TypesSolver::IC:
612 // IC can be set up in parallel via a sub-preconditioner to Block-Jacobi preconditioner
613 if (is_parallel) {
614 PCSetType(prec,PCBJACOBI); // set PC to Block Jacobi
615 PetscInt nlocal, first;
616 KSP *subksp; // KSP context for subdomain
617 PC subpc; // PC context for subdomain
618 KSPSetUp(m_ksp_solver);
619 PCBJacobiGetSubKSP(prec, &nlocal, &first, &subksp); // Extract the KSP context array for the local blocks
620 for (int i = 0; i < nlocal; i++) {
621 KSPGetPC(subksp[i], &subpc);
622 PCSetType(subpc, PCICC); // set subdomain PC to IC
623 }
624 }
625 else
626 PCSetType(prec,PCICC);
627 break;
628 // Sparse Approximate Inverse method of Grote and Barnard as a preconditioner (SIAM J. Sci. Comput.; vol 18, nr 3)
629 case TypesSolver::SPAIstat : PCSetType(prec,PCSPAI); break;
630 // Use multigrid preconditioning.
631 // This preconditioner requires you provide additional information
632 // about the coarser grid matrices and restriction/interpolation operators.
633 case TypesSolver::AMG:
634 // By default PCMG uses GMRES on the fine grid smoother so this should be used with KSPFGMRES
635 // or the smoother changed to not use GMRES
636 // Implements the Flexible Generalized Minimal Residual method. developed by Saad with restart
637 KSPSetType(m_ksp_solver,KSPFGMRES);
638 PCSetType(prec,PCMG);
639 //Sets the number of levels to use with MG. Must be called before any other MG routine.
640 PCMGSetLevels(prec,1,
641 (MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator()));
642 // Determines the form of multigrid to use: multiplicative, additive, full, or the Kaskade algorithm.
643 PCMGSetType(prec,PC_MG_MULTIPLICATIVE); // PC_MG_MULTIPLICATIVE,PC_MG_ADDITIVE,PC_MG_FULL,PC_MG_KASKADE
644 // Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more complicated cycling.
645 PCMGSetCycleType(prec,PC_MG_CYCLE_V); // PC_MG_CYCLE_V or PC_MG_CYCLE_W
646 //Sets the number of pre-smoothing steps to use on all levels
647#if PETSC_VERSION_GE(3,10,2)
648 PCMGSetNumberSmooth(prec, 1);
649#else
650 PCMGSetNumberSmoothDown(prec, 1);
651 PCMGSetNumberSmoothUp(prec, 1);
652#endif
653 // PCMGSetResidual: Sets the function to be used to calculate the residual on the lth level
654 // PCMGSetInterpolation: Sets the function to be used to calculate the interpolation from l-1 to the lth level
655 // PCMGSetRestriction: Sets the function to be used to restrict vector from level l to l-1
656 // PCMGSetRhs: Sets the vector space to be used to store the right-hand side on a particular level
657 // PCMGSetX: Sets the vector space to be used to store the solution on a particular level
658 // PCMGSetR: Sets the vector space to be used to store the residual on a particular level
659 //PCMGSetRestriction, PCMGSetRhs, PCMGSetX, PCMGSetR, etc.
660 //PCMGSetLevels
661 break;
662 case TypesSolver::AINV : throw ArgumentException("AlephMatrixPETSc", "preconditionnement AINV indisponible");
663 case TypesSolver::SPAIdyn : throw ArgumentException("AlephMatrixPETSc", "preconditionnement SPAIdyn indisponible");
664 case TypesSolver::ILUp : throw ArgumentException("AlephMatrixPETSc", "preconditionnement ILUp indisponible");
665 case TypesSolver::POLY : throw ArgumentException("AlephMatrixPETSc", "preconditionnement POLY indisponible");
666 default : throw ArgumentException("AlephMatrixPETSc", "preconditionnement inconnu");
667 }
668
669 if (solver_param->xoUser())
670 KSPSetInitialGuessNonzero(m_ksp_solver,PETSC_TRUE);
671
672 KSPSetTolerances(m_ksp_solver,
673 solver_param->epsilon(),
674 PETSC_DEFAULT,
675 PETSC_DEFAULT,
676 solver_param->maxIter());
677
678 // override KSP options from commandline if present//
679 KSPSetFromOptions(m_ksp_solver);
680
681 // override PC options from commandline if present//
682 PCSetFromOptions(prec) ;
683
684 // setup KSP solver and precondtioner //
685 KSPSetUp(m_ksp_solver);
686 PCSetUp(prec);
687
688 // print info which KSP solver and PC used //
689 KSPType kt;
690 KSPGetType(m_ksp_solver,&kt);
691 PCType pt;
692 PCGetType(prec,&pt);
693 info(4) << "[AlephMatrixSolve] Will use KSP type :" << kt << " PC type : "<< pt
694 << " is_parallel=" << is_parallel;
695
696 // solve the linear system //
697 debug()<<"[AlephMatrixSolve] All set up, now solving";
698 petscCheck(KSPSolve(m_ksp_solver,RHS,solution));
699 debug()<<"[AlephMatrixSolve] solved";
700 petscCheck(KSPGetConvergedReason(m_ksp_solver,&reason));
701
702 switch(reason){
703#if !PETSC_VERSION_(3,0,0)
704 case(KSP_CONVERGED_RTOL_NORMAL):{break;}
705 case(KSP_CONVERGED_ATOL_NORMAL):{break;}
706#endif
707 case(KSP_CONVERGED_RTOL):{break;}
708 case(KSP_CONVERGED_ATOL):{break;}
709 case(KSP_CONVERGED_ITS):{break;}
710#if !PETSC_VERSION_GE(3,19,0)
711 case(KSP_CONVERGED_CG_NEG_CURVE):{break;}
712 case(KSP_CONVERGED_CG_CONSTRAINED):{break;}
713#endif
714 case(KSP_CONVERGED_STEP_LENGTH):{break;}
715 case(KSP_CONVERGED_HAPPY_BREAKDOWN):{break;}
716 /* diverged */
717 case(KSP_DIVERGED_NULL):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NULL");}
718 case(KSP_DIVERGED_ITS):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_ITS");}
719 case(KSP_DIVERGED_DTOL):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_DTOL");}
720 case(KSP_DIVERGED_BREAKDOWN):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_BREAKDOWN");}
721 case(KSP_DIVERGED_BREAKDOWN_BICG):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_BREAKDOWN_BICG");}
722 case(KSP_DIVERGED_NONSYMMETRIC):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NONSYMMETRIC");}
723 case(KSP_DIVERGED_INDEFINITE_PC):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_INDEFINITE_PC");}
724#if PETSC_VERSION_GE(3,6,1)
725 case(KSP_DIVERGED_NANORINF):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NANORINF");}
726#else
727 case(KSP_DIVERGED_NAN):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_NAN");}
728#endif
729 case(KSP_DIVERGED_INDEFINITE_MAT):{ throw Exception("AlephMatrixPETSc::Solve", "KSP_DIVERGED_INDEFINITE_MAT");}
730
731 case(KSP_CONVERGED_ITERATING):{break;}
732 default: throw Exception("AlephMatrixPETSc::Solve", "");
733 }
734
735 KSPGetIterationNumber(m_ksp_solver,&its);
736 nb_iteration=its;
737
738 KSPGetResidualNorm(m_ksp_solver,&norm);
739 *residual_norm=norm;
740
741 return 0;
742}
743
744
745// ****************************************************************************
746// * writeToFile
747// ****************************************************************************
748
749void AlephMatrixPETSc::
750writeToFile(const String filename)
751{
752 ARCANE_UNUSED(filename);
753}
754
755/*---------------------------------------------------------------------------*/
756/*---------------------------------------------------------------------------*/
757
759
760/*---------------------------------------------------------------------------*/
761/*---------------------------------------------------------------------------*/
762
763}
764
765/*---------------------------------------------------------------------------*/
766/*---------------------------------------------------------------------------*/
#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.
AbstractService(const ServiceBuildInfo &)
Constructeur à partir d'un ServiceBuildInfo.
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 void * getMPICommunicator()=0
Adresse du communicateur MPI associé à ce gestionnaire.
Interface du gestionnaire de traces.
Structure contenant les informations pour créer un service.
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.
Vecteur 1D de données avec sémantique par valeur (style STL).
-*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
Int32 Integer
Type représentant un entier.
double Real
Type représentant un réel.
int AlephInt
Type par défaut pour indexer les lignes et les colonnes des matrices et vecteurs.
Definition AlephGlobal.h:50