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