source: issm/trunk/src/c/Solverx/Solverx.cpp@ 3226

/*!\file Solverx
 * \brief solver
 */

#include "./Solverx.h"

#undef __FUNCT__ 
#define __FUNCT__ "Solverx"

#include "../shared/shared.h"

#ifdef HAVE_CONFIG_H
        #include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif


void    Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0, char* solver_string){

        /*output: */
        Vec uf=NULL;

        /*intermediary: */
        int local_m,local_n;
        
        /*Solver*/
        KSP ksp=NULL; 
        PC  pc=NULL;
        int iteration_number;
        PetscTruth flag;
        int solver_type;

        /*If initial guess for solution exists, use it to create uf, otherwise, 
         * duplicate right hand side so that solution vector has same structure*/
        if(uf0){
                VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
        }
        else{
                MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVecFromLocalSize(local_n);
        }

        /*Before preparing the solver, add options to the options database*/
        PetscOptionsInsertMultipleString(solver_string);

        /*Process solver_string to see if we are not using special types of external solvers: */
        PetscOptionsDetermineSolverType(&solver_type,solver_string);

        #if _PETSC_VERSION_ == 2 
        if (solver_type==MUMPSPACKAGE_LU){
                /*Convert Kff to MATTAIJMUMPS: */
                MatConvert(Kff,MATAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
        }
        if (solver_type==MUMPSPACKAGE_CHOL){
                /*Convert Kff to MATTSBAIJMUMPS: */
                MatConvert(Kff,MATSBAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
        }
        if (solver_type==SPOOLESPACKAGE_LU){
                /*Convert Kff to MATTSBAIJMUMPS: */
                MatConvert(Kff,MATAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
        }
        if (solver_type==SPOOLESPACKAGE_CHOL){
                /*Convert Kff to MATTSBAIJMUMPS: */
                MatConvert(Kff,MATSBAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
        }
        if (solver_type==SUPERLUDISTPACKAGE){
                /*Convert Kff to MATTSBAIJMUMPS: */
                MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff);
        }
        #endif

        /*Prepare solver*/
        KSPCreate(MPI_COMM_WORLD,&ksp);
        KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
        KSPSetFromOptions(ksp);

        #if _PETSC_VERSION_ == 3 
        /*specific solver?: */
        KSPGetPC(ksp,&pc);
        if (solver_type==MUMPSPACKAGE_LU){
                PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
        }
        #endif


        /*If initial guess for 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);
                }
        }
                

        #ifdef _ISSM_DEBUG_ 
                KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
                KSPGetInitialGuessNonzero(ksp,&flag); 
                _printf_("User provided initial guess solution: %i\n",flag);
        #endif

        KSPSolve(ksp,pf,uf);
        
        /*Check convergence*/
        KSPGetIterationNumber(ksp,&iteration_number);
        if (iteration_number<0){
                throw ErrorException(__FUNCT__,exprintf("%s%i"," Solver diverged at iteration number: ",-iteration_number));
        }
        else{
                #ifdef _ISSM_DEBUG_ 
                _printf_("%s%i%s\n","Solver converged after ",iteration_number," iterations");
                #endif
        }

        /*Free ressources:*/
        KSPFree(&ksp);

        
        /*Assign output pointers:*/
        *puf=uf;
}