14#include "arcane/aleph/AlephArcane.h"
16#include "arcane/utils/StringList.h"
18#define MPI_COMM_SUB *(MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator())
24#if PETSC_VERSION_GE(3,6,1)
26#elif PETSC_VERSION_(3,3,0)
28#elif PETSC_VERSION_(3,0,0)
31#error "Bad version of PETSc: should be 3.0.0, 3.3.0 or greater than 3.6.1."
45 static_assert(
sizeof(PetscInt)==
sizeof(
AlephInt),
"AlephInt and PetscInt should have the same size");
46 return reinterpret_cast<PetscInt*
>(v);
49const PetscInt* _toPetscInt(
const AlephInt* v)
51 static_assert(
sizeof(PetscInt)==
sizeof(
AlephInt),
"AlephInt and PetscInt should have the same size");
52 return reinterpret_cast<const PetscInt*
>(v);
56petscCheck(PetscErrorCode error_code)
61#if PETSC_VERSION_GE(3,19,0)
62 if (error_code==PETSC_SUCCESS)
64 PetscCallVoid(error_code);
65 ARCANE_FATAL(
"Error in PETSc backend errcode={0}",error_code);
68 ARCANE_FATAL(
"Error in PETSc backend errcode={0}",error_code);
85 debug() <<
"\33[1;5;32m\t\t\t[~AlephTopologyPETSc]\33[0m";
90 void backupAndInitialize()
override {}
91 void restore()
override {}
95 void _checkInitPETSc();
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);
129 void AlephMatrixCreate(
void)
override;
130 void AlephMatrixSetFilled(
bool)
override;
131 int AlephMatrixAssemble(
void)
override;
138 void writeToFile(
const String)
override;
160 for (
auto* v : m_IAlephVectors)
162 for (
auto* v : m_IAlephMatrixs)
164 for (
auto* v : m_IAlephTopologys)
170 void initialize()
override {}
174 Integer nb_row_size)
override
182 Integer index)
override
191 Integer index)
override
219 if (!m_participating_in_solver) {
220 debug() <<
"\33[1;32m\t[AlephTopologyPETSc] Not concerned with this solver, returning\33[0m";
223 debug() <<
"\33[1;32m\t\t[AlephTopologyPETSc] @" <<
this <<
"\33[0m";
224 if (!m_kernel->isParallel()) {
225 PETSC_COMM_WORLD = PETSC_COMM_SELF;
237void AlephTopologyPETSc::
241 PetscBool is_petsc_initialized = PETSC_FALSE;
242 PetscInitialized(&is_petsc_initialized);
243 if (is_petsc_initialized==PETSC_TRUE)
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;
250 const CommandLineArguments& args = init_args.commandLineArguments();
251 PetscInitialize(args.commandLineArgc(),args.commandLineArgv(),
nullptr,
nullptr);
254 PetscInitializeNoArguments();
265AlephVectorPETSc(ITraceMng *tm,AlephKernel *kernel,Integer index)
266: IAlephVector(tm,kernel,index)
272 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorPETSc] new SolverVectorPETSc";
278void AlephVectorPETSc::
279AlephVectorCreate(
void)
281 for(
int iCpu=0;iCpu<m_kernel->size();++iCpu){
282 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
285 jLower=m_kernel->topology()->gathered_nb_row(iCpu);
286 jUpper=m_kernel->topology()->gathered_nb_row(iCpu+1);
290 VecCreateMPI(MPI_COMM_SUB,
292 m_kernel->topology()->nb_row_size(),
294 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorCreate] PETSC VectorCreate"
295 <<
" of local size=" <<jSize<<
"/"<<m_kernel->topology()->nb_row_size()
302void AlephVectorPETSc::
303AlephVectorSet(
const double *bfr_val,
const AlephInt *bfr_idx, Integer size)
305 VecSetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val,INSERT_VALUES);
306 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorSet] "<<size<<
" values inserted!";
308 AlephVectorAssemble();
315int AlephVectorPETSc::
316AlephVectorAssemble(
void)
318 VecAssemblyBegin(m_petsc_vector);
319 VecAssemblyEnd(m_petsc_vector);
320 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorAssemble]";
328void AlephVectorPETSc::
329AlephVectorGet(
double* bfr_val,
const AlephInt* bfr_idx, Integer size)
331 VecGetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val);
332 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorGet] fetched "<< size <<
" values!";
339void AlephVectorPETSc::
340writeToFile(
const String filename)
342 ARCANE_UNUSED(filename);
343 debug()<<
"\t\t[AlephVectorPETSc::writeToFile]";
349Real AlephVectorPETSc::
353 VecNorm(m_petsc_vector,NORM_INFINITY,&val);
361AlephMatrixPETSc(ITraceMng *tm,
364: IAlephMatrix(tm,kernel,index)
368 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixPETSc] new AlephMatrixPETSc";
375void AlephMatrixPETSc::
376AlephMatrixCreate(
void)
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);
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
395 PetscInt* row_data = _toPetscInt(m_kernel->topology()->gathered_nb_row_elements().subView(ilower,size).data());
396 PETSC_VERSION_MatCreate(MPI_COMM_SUB,
399 m_kernel->topology()->nb_row_size(),
400 m_kernel->topology()->nb_row_size(),
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)<<
")";
418void AlephMatrixPETSc::
419AlephMatrixSetFilled(
bool)
421 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixSetFilled] done";
422 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
429int AlephMatrixPETSc::
430AlephMatrixAssemble(
void)
432 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixAssemble] AlephMatrixAssemble";
433 MatAssemblyBegin(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
434 MatAssemblyEnd(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
442void AlephMatrixPETSc::
445 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] size="<<size;
446 for(
int i=0;i<size;i++){
448 MatSetValue(m_petsc_matrix, rows[i], cols[i], values[i], INSERT_VALUES);
450 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] done";
452 AlephMatrixAssemble();
460Real AlephMatrixPETSc::
461LinftyNormVectorProductAndSub(AlephVector* x,AlephVector* b)
465 throw FatalErrorException(
"LinftyNormVectorProductAndSub",
"error");
472bool AlephMatrixPETSc::
473isAlreadySolved(AlephVectorPETSc* x,
475 AlephVectorPETSc* tmp,
482 ARCANE_UNUSED(residual_norm);
483 ARCANE_UNUSED(params);
491int AlephMatrixPETSc::
492AlephMatrixSolve(AlephVector* x,AlephVector* b,AlephVector* t,
493 Integer& nb_iteration,
Real* residual_norm,
494 AlephParams* solver_param)
501 KSPConvergedReason reason;
503 AlephVectorPETSc* x_petsc =
dynamic_cast<AlephVectorPETSc*
> (x->implementation());
504 AlephVectorPETSc* b_petsc =
dynamic_cast<AlephVectorPETSc*
> (b->implementation());
508 const bool is_parallel = m_kernel->subParallelMng(m_index)->
isParallel();
510 Vec solution = x_petsc->m_petsc_vector;
511 Vec RHS = b_petsc->m_petsc_vector;
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);
518 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix,SAME_NONZERO_PATTERN);
521 switch(solver_param->method()){
523 case TypesSolver::PCG : KSPSetType(m_ksp_solver,KSPCG);
break;
525 case TypesSolver::BiCGStab : KSPSetType(m_ksp_solver,KSPBCGS);
break;
528 case TypesSolver::BiCGStab2: KSPSetType(m_ksp_solver,KSPIBCGS);
break;
530 case TypesSolver::GMRES : KSPSetType(m_ksp_solver,KSPGMRES);
break;
531 default :
throw ArgumentException(
"AlephMatrixPETSc::AlephMatrixSolve",
"Unknown solver method");
534 KSPGetPC(m_ksp_solver,&prec);
536 switch(solver_param->precond()){
537 case TypesSolver::NONE : PCSetType(prec,PCNONE);
break;
539 case TypesSolver::DIAGONAL : PCSetType(prec,PCJACOBI);
break;
541 case TypesSolver::ILU:
545 PCSetType(prec,PCBJACOBI);
547 PCSetType(prec,PCILU);
550 case TypesSolver::IC:
553 PCSetType(prec,PCBJACOBI);
554 PetscInt nlocal, first;
557 KSPSetUp(m_ksp_solver);
558 PCBJacobiGetSubKSP(prec, &nlocal, &first, &subksp);
559 for (
int i = 0; i < nlocal; i++) {
560 KSPGetPC(subksp[i], &subpc);
561 PCSetType(subpc, PCICC);
565 PCSetType(prec,PCICC);
568 case TypesSolver::SPAIstat : PCSetType(prec,PCSPAI);
break;
572 case TypesSolver::AMG:
576 KSPSetType(m_ksp_solver,KSPFGMRES);
577 PCSetType(prec,PCMG);
579 PCMGSetLevels(prec,1,
582 PCMGSetType(prec,PC_MG_MULTIPLICATIVE);
584 PCMGSetCycleType(prec,PC_MG_CYCLE_V);
586#if PETSC_VERSION_GE(3,10,2)
587 PCMGSetNumberSmooth(prec, 1);
589 PCMGSetNumberSmoothDown(prec, 1);
590 PCMGSetNumberSmoothUp(prec, 1);
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");
608 if (solver_param->xoUser())
609 KSPSetInitialGuessNonzero(m_ksp_solver,PETSC_TRUE);
611 KSPSetTolerances(m_ksp_solver,
612 solver_param->epsilon(),
615 solver_param->maxIter());
618 KSPSetFromOptions(m_ksp_solver);
621 PCSetFromOptions(prec) ;
624 KSPSetUp(m_ksp_solver);
629 KSPGetType(m_ksp_solver,&kt);
632 info(4) <<
"[AlephMatrixSolve] Will use KSP type :" << kt <<
" PC type : "<< pt
633 <<
" is_parallel=" << is_parallel;
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));
642#if !PETSC_VERSION_(3,0,0)
643 case(KSP_CONVERGED_RTOL_NORMAL):{
break;}
644 case(KSP_CONVERGED_ATOL_NORMAL):{
break;}
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;}
653 case(KSP_CONVERGED_STEP_LENGTH):{
break;}
654 case(KSP_CONVERGED_HAPPY_BREAKDOWN):{
break;}
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");}
666 case(KSP_DIVERGED_NAN):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NAN");}
668 case(KSP_DIVERGED_INDEFINITE_MAT):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_INDEFINITE_MAT");}
670 case(KSP_CONVERGED_ITERATING):{
break;}
671 default:
throw Exception(
"AlephMatrixPETSc::Solve",
"");
674 KSPGetIterationNumber(m_ksp_solver,&its);
677 KSPGetResidualNorm(m_ksp_solver,&norm);
688void AlephMatrixPETSc::
689writeToFile(
const String filename)
691 ARCANE_UNUSED(filename);
#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.
Vecteur d'un système linéaire.
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.
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.