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

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

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