[11679] | 1 | /*!\file SolverxPetsc
|
---|
| 2 | * \brief Petsc implementation of solver
|
---|
| 3 | */
|
---|
| 4 |
|
---|
| 5 | #include "./Solverx.h"
|
---|
| 6 | #include "../../shared/shared.h"
|
---|
| 7 | #include "../../include/include.h"
|
---|
| 8 | #include "../../io/io.h"
|
---|
| 9 |
|
---|
| 10 | #ifdef HAVE_CONFIG_H
|
---|
| 11 | #include <config.h>
|
---|
| 12 | #else
|
---|
| 13 | #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
|
---|
| 14 | #endif
|
---|
| 15 |
|
---|
[13395] | 16 | void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){
|
---|
| 17 |
|
---|
| 18 | PetscVec* uf=new PetscVec();
|
---|
| 19 |
|
---|
| 20 | Vec uf0_vector = NULL;
|
---|
| 21 | Vec df_vector = NULL;
|
---|
| 22 |
|
---|
| 23 | if(uf0) uf0_vector = uf0->vector;
|
---|
| 24 | if(df) df_vector = df->vector;
|
---|
| 25 |
|
---|
| 26 | SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
|
---|
| 27 |
|
---|
| 28 | *puf=uf;
|
---|
| 29 |
|
---|
| 30 | }
|
---|
[11679] | 31 | void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
|
---|
| 32 |
|
---|
| 33 | /*Output: */
|
---|
| 34 | Vec uf = NULL;
|
---|
| 35 |
|
---|
| 36 | /*Intermediary: */
|
---|
| 37 | int local_m,local_n,global_m,global_n;
|
---|
| 38 | int analysis_type;
|
---|
| 39 |
|
---|
| 40 | /*Solver */
|
---|
| 41 | KSP ksp = NULL;
|
---|
| 42 | PC pc = NULL;
|
---|
| 43 | int iteration_number;
|
---|
| 44 | int solver_type;
|
---|
| 45 | bool fromlocalsize = true;
|
---|
| 46 | #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
|
---|
| 47 | PetscTruth flag,flg;
|
---|
| 48 | #else
|
---|
| 49 | PetscBool flag,flg;
|
---|
| 50 | #endif
|
---|
| 51 |
|
---|
| 52 | /*Stokes: */
|
---|
| 53 | IS isv=NULL;
|
---|
| 54 | IS isp=NULL;
|
---|
| 55 |
|
---|
| 56 | #if _PETSC_MAJOR_ >= 3
|
---|
| 57 | char ksp_type[50];
|
---|
| 58 | #endif
|
---|
| 59 |
|
---|
| 60 | /*Display message*/
|
---|
[12706] | 61 | if(VerboseModule()) _pprintLine_(" Solving");
|
---|
[11679] | 62 | #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
|
---|
| 63 | if(VerboseSolver())PetscOptionsPrint(stdout);
|
---|
| 64 | #else
|
---|
| 65 | if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
|
---|
| 66 | #endif
|
---|
| 67 |
|
---|
[12330] | 68 | /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
|
---|
[11679] | 69 | _assert_(Kff);
|
---|
| 70 | MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
|
---|
| 71 | if(!global_n){
|
---|
[13395] | 72 | *puf=NewVec(0); return;
|
---|
[11679] | 73 | }
|
---|
[12330] | 74 |
|
---|
| 75 | /*Initial guess */
|
---|
[11679] | 76 | /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
|
---|
| 77 | #if _PETSC_MAJOR_ >= 3
|
---|
| 78 | PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
|
---|
| 79 | if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
|
---|
| 80 | #endif
|
---|
| 81 |
|
---|
| 82 | /*If initial guess for the solution exists, use it to create uf, otherwise,
|
---|
| 83 | * duplicate the right hand side so that the solution vector has the same structure*/
|
---|
| 84 | if(uf0){
|
---|
| 85 | VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
|
---|
| 86 | }
|
---|
| 87 | else{
|
---|
| 88 | MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
|
---|
| 89 | }
|
---|
[12330] | 90 |
|
---|
| 91 | /*Process petsc options to see if we are using special types of external solvers*/
|
---|
[11679] | 92 | PetscOptionsDetermineSolverType(&solver_type);
|
---|
| 93 |
|
---|
[12330] | 94 | /*Check the solver is available*/
|
---|
[11679] | 95 | if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
|
---|
[12330] | 96 | #if _PETSC_MAJOR_ >=3
|
---|
| 97 | #ifndef _HAVE_MUMPS_
|
---|
[13395] | 98 | _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
|
---|
[12330] | 99 | #endif
|
---|
[11679] | 100 | #endif
|
---|
[12330] | 101 | }
|
---|
[11679] | 102 |
|
---|
[12330] | 103 | /*Prepare solver*/
|
---|
[11679] | 104 | KSPCreate(MPI_COMM_WORLD,&ksp);
|
---|
| 105 | KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
|
---|
| 106 | KSPSetFromOptions(ksp);
|
---|
| 107 |
|
---|
[12330] | 108 | #if _PETSC_MAJOR_==3
|
---|
[11679] | 109 | /*Specific solver?: */
|
---|
| 110 | KSPGetPC(ksp,&pc);
|
---|
| 111 | if (solver_type==MUMPSPACKAGE_LU){
|
---|
| 112 | #if _PETSC_MINOR_==1
|
---|
| 113 | PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
|
---|
| 114 | #else
|
---|
| 115 | PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
|
---|
| 116 | #endif
|
---|
| 117 | }
|
---|
| 118 |
|
---|
| 119 | /*Stokes: */
|
---|
| 120 | if (solver_type==StokesSolverEnum){
|
---|
| 121 | /*Make indices out of doftypes: */
|
---|
[13395] | 122 | if(!df)_error_("need doftypes for Stokes solver!\n");
|
---|
[11679] | 123 | DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
|
---|
| 124 |
|
---|
| 125 | /*Set field splits: */
|
---|
| 126 | KSPGetPC(ksp,&pc);
|
---|
| 127 | #if _PETSC_MINOR_==1
|
---|
| 128 | PCFieldSplitSetIS(pc,isv);
|
---|
| 129 | PCFieldSplitSetIS(pc,isp);
|
---|
| 130 | #else
|
---|
| 131 | PCFieldSplitSetIS(pc,PETSC_NULL,isv);
|
---|
| 132 | PCFieldSplitSetIS(pc,PETSC_NULL,isp);
|
---|
| 133 | #endif
|
---|
| 134 |
|
---|
| 135 | }
|
---|
| 136 | #endif
|
---|
| 137 |
|
---|
[12330] | 138 | /*If there is an initial guess for the solution, use it
|
---|
| 139 | * except if we are using the MUMPS direct solver
|
---|
| 140 | * where any initial guess will crash Petsc*/
|
---|
[11679] | 141 | if (uf0){
|
---|
[12330] | 142 | if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
|
---|
[11679] | 143 | KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
|
---|
| 144 | }
|
---|
| 145 | }
|
---|
| 146 |
|
---|
| 147 | /*Solve: */
|
---|
[12330] | 148 | if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
|
---|
[11679] | 149 | KSPSolve(ksp,pf,uf);
|
---|
| 150 |
|
---|
| 151 | /*Check convergence*/
|
---|
| 152 | KSPGetIterationNumber(ksp,&iteration_number);
|
---|
[13395] | 153 | if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
|
---|
[11679] | 154 |
|
---|
| 155 | /*Free resources:*/
|
---|
| 156 | KSPFree(&ksp);
|
---|
| 157 |
|
---|
| 158 | /*Assign output pointers:*/
|
---|
| 159 | *puf=uf;
|
---|
| 160 | }
|
---|