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();
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);
117 Vec m_petsc_vector =
nullptr;
133 void AlephMatrixCreate(
void)
override;
134 void AlephMatrixSetFilled(
bool)
override;
135 int AlephMatrixAssemble(
void)
override;
142 void writeToFile(
const String)
override;
146 Mat m_petsc_matrix =
nullptr;
147 KSP m_ksp_solver =
nullptr;
171 for (
auto* v : m_IAlephVectors)
173 for (
auto* v : m_IAlephMatrixs)
175 for (
auto* v : m_IAlephTopologys)
189 void initialize()
override {}
193 Integer nb_row_size)
override
201 Integer index)
override
210 Integer index)
override
235 if (!m_participating_in_solver) {
236 debug() <<
"\33[1;32m\t[AlephTopologyPETSc] Not concerned with this solver, returning\33[0m";
239 debug() <<
"\33[1;32m\t\t[AlephTopologyPETSc] @" <<
this <<
"\33[0m";
240 if (!m_kernel->isParallel()) {
241 PETSC_COMM_WORLD = PETSC_COMM_SELF;
253void AlephTopologyPETSc::
257 PetscBool is_petsc_initialized = PETSC_FALSE;
258 PetscInitialized(&is_petsc_initialized);
259 if (is_petsc_initialized==PETSC_TRUE)
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;
267 const CommandLineArguments& args = init_args.commandLineArguments();
268 PetscInitialize(args.commandLineArgc(),args.commandLineArgv(),
nullptr,
nullptr);
271 PetscInitializeNoArguments();
273 global_need_petsc_finalize =
true;
283AlephVectorPETSc(ITraceMng *tm,AlephKernel *kernel,Integer index)
284: IAlephVector(tm,kernel,index)
290 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorPETSc] new SolverVectorPETSc";
299 VecDestroy(&m_petsc_vector);
308void AlephVectorPETSc::
309AlephVectorCreate(
void)
311 for(
int iCpu=0;iCpu<m_kernel->size();++iCpu){
312 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
315 jLower=m_kernel->topology()->gathered_nb_row(iCpu);
316 jUpper=m_kernel->topology()->gathered_nb_row(iCpu+1);
320 VecCreateMPI(MPI_COMM_SUB,
322 m_kernel->topology()->nb_row_size(),
324 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorCreate] PETSC VectorCreate"
325 <<
" of local size=" <<jSize<<
"/"<<m_kernel->topology()->nb_row_size()
332void AlephVectorPETSc::
333AlephVectorSet(
const double *bfr_val,
const AlephInt *bfr_idx, Integer size)
335 VecSetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val,INSERT_VALUES);
336 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorSet] "<<size<<
" values inserted!";
338 AlephVectorAssemble();
345int AlephVectorPETSc::
346AlephVectorAssemble(
void)
348 VecAssemblyBegin(m_petsc_vector);
349 VecAssemblyEnd(m_petsc_vector);
350 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorAssemble]";
358void AlephVectorPETSc::
359AlephVectorGet(
double* bfr_val,
const AlephInt* bfr_idx, Integer size)
361 VecGetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val);
362 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorGet] fetched "<< size <<
" values!";
369void AlephVectorPETSc::
370writeToFile(
const String filename)
372 ARCANE_UNUSED(filename);
373 debug()<<
"\t\t[AlephVectorPETSc::writeToFile]";
379Real AlephVectorPETSc::
383 VecNorm(m_petsc_vector,NORM_INFINITY,&val);
395AlephMatrixPETSc(ITraceMng *tm, AlephKernel *kernel, Integer index)
396: IAlephMatrix(tm,kernel,index)
398 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixPETSc] new AlephMatrixPETSc";
408 MatDestroy(&m_petsc_matrix);
410 KSPDestroy(&m_ksp_solver);
419void AlephMatrixPETSc::
420AlephMatrixCreate(
void)
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);
434#if PETSC_VERSION_GE(3,3,0)
435 #define PETSC_VERSION_MatCreate MatCreateAIJ
436#elif PETSC_VERSION_(3,0,0)
437 #define PETSC_VERSION_MatCreate MatCreateMPIAIJ
439 PetscInt* row_data = _toPetscInt(m_kernel->topology()->gathered_nb_row_elements().subView(ilower,size).data());
440 PETSC_VERSION_MatCreate(MPI_COMM_SUB,
443 m_kernel->topology()->nb_row_size(),
444 m_kernel->topology()->nb_row_size(),
452 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
453 MatSetUp(m_petsc_matrix);
454 debug()<<
"\t\t[AlephMatrixPetsc::AlephMatrixCreate] PETSC MatrixCreate idx:"<<m_index
455 <<
", ("<<ilower<<
"->"<<(iupper-1)<<
")";
462void AlephMatrixPETSc::
463AlephMatrixSetFilled(
bool)
465 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixSetFilled] done";
466 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
473int AlephMatrixPETSc::
474AlephMatrixAssemble(
void)
476 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixAssemble] AlephMatrixAssemble";
477 MatAssemblyBegin(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
478 MatAssemblyEnd(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
486void AlephMatrixPETSc::
489 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] size="<<size;
490 for(
int i=0;i<size;i++){
492 MatSetValue(m_petsc_matrix, rows[i], cols[i], values[i], INSERT_VALUES);
494 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] done";
496 AlephMatrixAssemble();
504Real AlephMatrixPETSc::
505LinftyNormVectorProductAndSub(AlephVector* x,AlephVector* b)
509 throw FatalErrorException(
"LinftyNormVectorProductAndSub",
"error");
516bool AlephMatrixPETSc::
517isAlreadySolved(AlephVectorPETSc* x,
519 AlephVectorPETSc* tmp,
526 ARCANE_UNUSED(residual_norm);
527 ARCANE_UNUSED(params);
535int AlephMatrixPETSc::
536AlephMatrixSolve(AlephVector* x,AlephVector* b,AlephVector* t,
537 Integer& nb_iteration,
Real* residual_norm,
538 AlephParams* solver_param)
545 KSPConvergedReason reason;
547 AlephVectorPETSc* x_petsc =
dynamic_cast<AlephVectorPETSc*
> (x->implementation());
548 AlephVectorPETSc* b_petsc =
dynamic_cast<AlephVectorPETSc*
> (b->implementation());
552 const bool is_parallel = m_kernel->subParallelMng(m_index)->
isParallel();
554 Vec solution = x_petsc->m_petsc_vector;
555 Vec RHS = b_petsc->m_petsc_vector;
557 debug()<<
"[AlephMatrixSolve]";
559 KSPDestroy(&m_ksp_solver);
560 KSPCreate(MPI_COMM_SUB,&m_ksp_solver);
561#if PETSC_VERSION_GE(3,6,1)
562 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix);
564 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix,SAME_NONZERO_PATTERN);
567 switch(solver_param->method()){
569 case TypesSolver::PCG : KSPSetType(m_ksp_solver,KSPCG);
break;
571 case TypesSolver::BiCGStab : KSPSetType(m_ksp_solver,KSPBCGS);
break;
574 case TypesSolver::BiCGStab2: KSPSetType(m_ksp_solver,KSPIBCGS);
break;
576 case TypesSolver::GMRES : KSPSetType(m_ksp_solver,KSPGMRES);
break;
577 default :
throw ArgumentException(
"AlephMatrixPETSc::AlephMatrixSolve",
"Unknown solver method");
580 KSPGetPC(m_ksp_solver,&prec);
582 switch(solver_param->precond()){
583 case TypesSolver::NONE : PCSetType(prec,PCNONE);
break;
585 case TypesSolver::DIAGONAL : PCSetType(prec,PCJACOBI);
break;
587 case TypesSolver::ILU:
591 PCSetType(prec,PCBJACOBI);
593 PCSetType(prec,PCILU);
596 case TypesSolver::IC:
599 PCSetType(prec,PCBJACOBI);
600 PetscInt nlocal, first;
603 KSPSetUp(m_ksp_solver);
604 PCBJacobiGetSubKSP(prec, &nlocal, &first, &subksp);
605 for (
int i = 0; i < nlocal; i++) {
606 KSPGetPC(subksp[i], &subpc);
607 PCSetType(subpc, PCICC);
611 PCSetType(prec,PCICC);
614 case TypesSolver::SPAIstat : PCSetType(prec,PCSPAI);
break;
618 case TypesSolver::AMG:
622 KSPSetType(m_ksp_solver,KSPFGMRES);
623 PCSetType(prec,PCMG);
625 PCMGSetLevels(prec,1,
628 PCMGSetType(prec,PC_MG_MULTIPLICATIVE);
630 PCMGSetCycleType(prec,PC_MG_CYCLE_V);
632#if PETSC_VERSION_GE(3,10,2)
633 PCMGSetNumberSmooth(prec, 1);
635 PCMGSetNumberSmoothDown(prec, 1);
636 PCMGSetNumberSmoothUp(prec, 1);
647 case TypesSolver::AINV :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement AINV indisponible");
648 case TypesSolver::SPAIdyn :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement SPAIdyn indisponible");
649 case TypesSolver::ILUp :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement ILUp indisponible");
650 case TypesSolver::POLY :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement POLY indisponible");
651 default :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement inconnu");
654 if (solver_param->xoUser())
655 KSPSetInitialGuessNonzero(m_ksp_solver,PETSC_TRUE);
657 KSPSetTolerances(m_ksp_solver,
658 solver_param->epsilon(),
661 solver_param->maxIter());
664 KSPSetFromOptions(m_ksp_solver);
667 PCSetFromOptions(prec) ;
670 KSPSetUp(m_ksp_solver);
675 KSPGetType(m_ksp_solver,&kt);
678 info(4) <<
"[AlephMatrixSolve] Will use KSP type :" << kt <<
" PC type : "<< pt
679 <<
" is_parallel=" << is_parallel;
682 debug()<<
"[AlephMatrixSolve] All set up, now solving";
683 petscCheck(KSPSolve(m_ksp_solver,RHS,solution));
684 debug()<<
"[AlephMatrixSolve] solved";
685 petscCheck(KSPGetConvergedReason(m_ksp_solver,&reason));
688#if !PETSC_VERSION_(3,0,0)
689 case(KSP_CONVERGED_RTOL_NORMAL):{
break;}
690 case(KSP_CONVERGED_ATOL_NORMAL):{
break;}
692 case(KSP_CONVERGED_RTOL):{
break;}
693 case(KSP_CONVERGED_ATOL):{
break;}
694 case(KSP_CONVERGED_ITS):{
break;}
695#if !PETSC_VERSION_GE(3,19,0)
696 case(KSP_CONVERGED_CG_NEG_CURVE):{
break;}
697 case(KSP_CONVERGED_CG_CONSTRAINED):{
break;}
699 case(KSP_CONVERGED_STEP_LENGTH):{
break;}
700 case(KSP_CONVERGED_HAPPY_BREAKDOWN):{
break;}
702 case(KSP_DIVERGED_NULL):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NULL");}
703 case(KSP_DIVERGED_ITS):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_ITS");}
704 case(KSP_DIVERGED_DTOL):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_DTOL");}
705 case(KSP_DIVERGED_BREAKDOWN):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_BREAKDOWN");}
706 case(KSP_DIVERGED_BREAKDOWN_BICG):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_BREAKDOWN_BICG");}
707 case(KSP_DIVERGED_NONSYMMETRIC):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NONSYMMETRIC");}
708 case(KSP_DIVERGED_INDEFINITE_PC):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_INDEFINITE_PC");}
709#if PETSC_VERSION_GE(3,6,1)
710 case(KSP_DIVERGED_NANORINF):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NANORINF");}
712 case(KSP_DIVERGED_NAN):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NAN");}
714 case(KSP_DIVERGED_INDEFINITE_MAT):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_INDEFINITE_MAT");}
716 case(KSP_CONVERGED_ITERATING):{
break;}
717 default:
throw Exception(
"AlephMatrixPETSc::Solve",
"");
720 KSPGetIterationNumber(m_ksp_solver,&its);
723 KSPGetResidualNorm(m_ksp_solver,&norm);
734void AlephMatrixPETSc::
735writeToFile(
const String filename)
737 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.