source: issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp@ 20553

Last change on this file since 20553 was 20553, checked in by Mathieu Morlighem, 9 years ago

CHG: added support for PETSc 3.7
File size: 6.4 KB
Line 
1/*!\file SolverxPetsc
2 * \brief Petsc implementation of solver
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include "./PetscSolver.h"
12#include "../../../shared/Numerics/Verbosity.h"
13#include "../../../shared/MemOps/MemOps.h"
14#include "../../../shared/Exceptions/exceptions.h"
15#include "../../../shared/io/Comm/IssmComm.h"
16#include "../../../shared/Enum/Enum.h"
17
18void PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ /*{{{*/
19
20 PetscVec* uf=new PetscVec();
21
22 Vec uf0_vector = NULL;
23 Vec df_vector = NULL;
24
25 if(uf0) uf0_vector = uf0->vector;
26 if(df) df_vector = df->vector;
27
28 SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
29
30 *puf=uf;
31
32}
33/*}}}*/
34void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ /*{{{*/
35
36 /*Output: */
37 Vec uf = NULL;
38
39 /*Intermediary: */
40 int local_m,local_n,global_m,global_n;
41
42 /*Solver */
43 KSP ksp = NULL;
44 PC pc = NULL;
45 int iteration_number;
46 int solver_type;
47 bool fromlocalsize = true;
48 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
49 PetscTruth flag,flg;
50 #else
51 PetscBool flag,flg;
52 #endif
53
54 /*FS: */
55 IS isv=NULL;
56 IS isp=NULL;
57
58 #if _PETSC_MAJOR_ >= 3
59 char ksp_type[50];
60 #endif
61
62 /*Display message*/
63 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
64 if(VerboseSolver())PetscOptionsPrint(stdout);
65 #else
66 #if _PETSC_MINOR_<7
67 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
68 #else
69 if(VerboseSolver())PetscOptionsView(NULL,PETSC_VIEWER_STDOUT_WORLD);
70 #endif
71 #endif
72
73 /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
74 _assert_(Kff);
75 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);
76 if(!global_n){
77 *puf=NewVec(0,IssmComm::GetComm()); return;
78 }
79
80 /*Initial guess */
81 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
82 #if _PETSC_MAJOR_ >= 3
83 #if _PETSC_MINOR_<7
84 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
85 #else
86 PetscOptionsGetString(NULL,PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
87 #endif
88 if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
89 #endif
90
91 /*If initial guess for the solution exists, use it to create uf, otherwise,
92 * duplicate the right hand side so that the solution vector has the same structure*/
93 if(uf0){
94 VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
95 }
96 else{
97 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
98 }
99
100 /*Process petsc options to see if we are using special types of external solvers*/
101 PetscOptionsDetermineSolverType(&solver_type);
102
103 /*Check the solver is available*/
104 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
105 #if _PETSC_MAJOR_ >=3
106 #ifndef _HAVE_MUMPS_
107 _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
108 #endif
109 #endif
110 }
111
112 /*Prepare solver*/
113 KSPCreate(IssmComm::GetComm(),&ksp);
114 #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5)
115 KSPSetOperators(ksp,Kff,Kff);
116 #else
117 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
118 #endif
119 KSPSetFromOptions(ksp);
120
121 #if _PETSC_MAJOR_==3
122 /*Specific solver?: */
123 KSPGetPC(ksp,&pc);
124 if (solver_type==MUMPSPACKAGE_LU){
125 #if _PETSC_MINOR_==1
126 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
127 #else
128 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
129 #endif
130 }
131
132 /*FS: */
133 if (solver_type==FSSolverEnum){
134 /*Make indices out of doftypes: */
135 if(!df)_error_("need doftypes for FS solver!\n");
136 DofTypesToIndexSet(&isv,&isp,df,FSSolverEnum);
137
138 /*Set field splits: */
139 KSPGetPC(ksp,&pc);
140 #if _PETSC_MINOR_==1
141 PCFieldSplitSetIS(pc,isv);
142 PCFieldSplitSetIS(pc,isp);
143 #else
144 PCFieldSplitSetIS(pc,PETSC_NULL,isv);
145 PCFieldSplitSetIS(pc,PETSC_NULL,isp);
146 #endif
147
148 }
149 #endif
150
151 /*If there is an initial guess for the solution, use it
152 * except if we are using the MUMPS direct solver
153 * where any initial guess will crash Petsc*/
154 if (uf0){
155 if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
156 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
157 }
158 }
159
160 /*Solve: */
161 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
162 KSPSolve(ksp,pf,uf);
163
164 /*Check convergence*/
165 KSPGetIterationNumber(ksp,&iteration_number);
166 if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
167
168 /*Free resources:*/
169 KSPFree(&ksp);
170
171 /*Assign output pointers:*/
172 *puf=uf;
173}
174/*}}}*/
175void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){ /*{{{*/
176
177 /*output: */
178 IS isv=NULL;
179 IS isp=NULL;
180
181 int start,end;
182 IssmDouble* df_local=NULL;
183 int df_local_size;
184
185 int* pressure_indices=NULL;
186 int* velocity_indices=NULL;
187 int pressure_num=0;
188 int velocity_num=0;
189 int pressure_count=0;
190 int velocity_count=0;
191
192 if(typeenum==FSSolverEnum){
193
194 /*Ok, recover doftypes vector values and indices: */
195 VecGetOwnershipRange(df,&start,&end);
196 VecGetLocalSize(df,&df_local_size);
197 VecGetArray(df,&df_local);
198
199 pressure_num=0;
200 velocity_num=0;
201 for(int i=0;i<df_local_size;i++){
202 if (df_local[i]==PressureEnum)pressure_num++;
203 else velocity_num++;
204 }
205
206 /*Allocate indices: */
207 if(pressure_num)pressure_indices=xNew<int>(pressure_num);
208 if(velocity_num)velocity_indices=xNew<int>(velocity_num);
209
210 pressure_count=0;
211 velocity_count=0;
212 for(int i=0;i<df_local_size;i++){
213 if (df_local[i]==PressureEnum){
214 pressure_indices[pressure_count]=start+i;
215 pressure_count++;
216 }
217 if (df_local[i]==VelocityEnum){
218 velocity_indices[velocity_count]=start+i;
219 velocity_count++;
220 }
221 }
222 VecRestoreArray(df,&df_local);
223
224 /*Create indices sets: */
225 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
226 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
227 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
228 #else
229 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
230 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
231 #endif
232 }
233
234 /*Free ressources:*/
235 xDelete<int>(pressure_indices);
236 xDelete<int>(velocity_indices);
237
238 /*Assign output pointers:*/
239 *pisv=isv;
240 *pisp=isp;
241}
242/*}}}*/
Note: See TracBrowser for help on using the repository browser.