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

Last change on this file since 14792 was 14792, checked in by Eric.Larour, 12 years ago

CHG: finished implementation of IssmMpiVec and IssmMpiDenseMat objects. Still need to V&V.
CHG: moved the src/c/classes/matrix/Matrix and Vector objects to src/c/classes/toolkits.
CHG and NEW: offloaded the solver code to the toolkits. To do so, created a Solver class
in the src/c/classes/toolkits, which connects the toolkits to the abstract classes Matrix
and Vector.

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