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

Last change on this file since 27157 was 27157, checked in by Mathieu Morlighem, 3 years ago

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