/*!\file SolverxPetsc * \brief Petsc implementation of solver */ #include "./Solverx.h" #include "../../shared/shared.h" #include "../../include/include.h" #include "../../io/io.h" #ifdef HAVE_CONFIG_H #include #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ PetscVec* uf=new PetscVec(); Vec uf0_vector = NULL; Vec df_vector = NULL; if(uf0) uf0_vector = uf0->vector; if(df) df_vector = df->vector; SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters); *puf=uf; } void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ /*Output: */ Vec uf = NULL; /*Intermediary: */ int local_m,local_n,global_m,global_n; int analysis_type; /*Solver */ KSP ksp = NULL; PC pc = NULL; int iteration_number; int solver_type; bool fromlocalsize = true; #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) PetscTruth flag,flg; #else PetscBool flag,flg; #endif /*Stokes: */ IS isv=NULL; IS isp=NULL; #if _PETSC_MAJOR_ >= 3 char ksp_type[50]; #endif /*Display message*/ if(VerboseModule()) _pprintLine_(" Solving"); #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) if(VerboseSolver())PetscOptionsPrint(stdout); #else if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD); #endif /*First, check that f-set is not NULL, i.e. model is fully constrained:*/ _assert_(Kff); MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m); if(!global_n){ *puf=NewVec(0); return; } /*Initial guess */ /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */ #if _PETSC_MAJOR_ >= 3 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg); if (strcmp(ksp_type,"preonly")==0)uf0=NULL; #endif /*If initial guess for the solution exists, use it to create uf, otherwise, * duplicate the right hand side so that the solution vector has the same structure*/ if(uf0){ VecDuplicate(uf0,&uf); VecCopy(uf0,uf); } else{ MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize); } /*Process petsc options to see if we are using special types of external solvers*/ PetscOptionsDetermineSolverType(&solver_type); /*Check the solver is available*/ if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){ #if _PETSC_MAJOR_ >=3 #ifndef _HAVE_MUMPS_ _error_("requested MUMPS solver, which was not compiled into ISSM!\n"); #endif #endif } /*Prepare solver*/ KSPCreate(MPI_COMM_WORLD,&ksp); KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN); KSPSetFromOptions(ksp); #if _PETSC_MAJOR_==3 /*Specific solver?: */ KSPGetPC(ksp,&pc); if (solver_type==MUMPSPACKAGE_LU){ #if _PETSC_MINOR_==1 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS); #else PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS); #endif } /*Stokes: */ if (solver_type==StokesSolverEnum){ /*Make indices out of doftypes: */ if(!df)_error_("need doftypes for Stokes solver!\n"); DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum); /*Set field splits: */ KSPGetPC(ksp,&pc); #if _PETSC_MINOR_==1 PCFieldSplitSetIS(pc,isv); PCFieldSplitSetIS(pc,isp); #else PCFieldSplitSetIS(pc,PETSC_NULL,isv); PCFieldSplitSetIS(pc,PETSC_NULL,isp); #endif } #endif /*If there is an initial guess for the solution, use it * except if we are using the MUMPS direct solver * where any initial guess will crash Petsc*/ if (uf0){ if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){ KSPSetInitialGuessNonzero(ksp,PETSC_TRUE); } } /*Solve: */ if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); KSPSolve(ksp,pf,uf); /*Check convergence*/ KSPGetIterationNumber(ksp,&iteration_number); if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number); /*Free resources:*/ KSPFree(&ksp); /*Assign output pointers:*/ *puf=uf; }