Ice Sheet System Model  4.18
Code documentation
Functions
PetscSolver.h File Reference
#include "../petscincludes.h"

Go to the source code of this file.

Functions

void PetscSolve (PetscVec **puf, PetscMat *Kff, PetscVec *pf, PetscVec *uf0, PetscVec *df, Parameters *parameters)
 
void SolverxPetsc (Vec *puf, Mat Kff, Vec pf, Vec uf0, Vec df, Parameters *parameters)
 
void DofTypesToIndexSet (IS *pisv, IS *pisp, Vec df, int typeenum)
 

Function Documentation

◆ PetscSolve()

void PetscSolve ( PetscVec **  puf,
PetscMat Kff,
PetscVec pf,
PetscVec uf0,
PetscVec df,
Parameters parameters 
)

Definition at line 19 of file PetscSolver.cpp.

19  { /*{{{*/
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 }

◆ SolverxPetsc()

void SolverxPetsc ( Vec *  puf,
Mat  Kff,
Vec  pf,
Vec  uf0,
Vec  df,
Parameters parameters 
)

Definition at line 35 of file PetscSolver.cpp.

35  { /*{{{*/
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  /*FS: */
56  IS isv=NULL;
57  IS isp=NULL;
58  char ksp_type[50];
59 
60  /*Display message*/
61  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
62  if(VerboseSolver())PetscOptionsPrint(stdout);
63  #else
64  #if _PETSC_MINOR_<7
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_MINOR_<7
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 (strcmp(ksp_type,"preonly")==0)uf0=NULL;
86 
87  /*If initial guess for the solution exists, use it to create uf, otherwise,
88  * duplicate the right hand side so that the solution vector has the same structure*/
89  if(uf0){
90  VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
91  }
92  else{
93  MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
94  }
95 
96  /*Process petsc options to see if we are using special types of external solvers*/
97  PetscOptionsDetermineSolverType(&solver_type);
98 
99  /*Check the solver is available*/
100  if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
101  #ifndef _HAVE_MUMPS_
102  _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
103  #endif
104  }
105 
106  /*Prepare solver*/
107  KSPCreate(IssmComm::GetComm(),&ksp);
108  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=5)
109  KSPSetOperators(ksp,Kff,Kff);
110  #else
111  KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
112  #endif
113  KSPSetFromOptions(ksp);
114 
115  /*Specific solver?: */
116  KSPGetPC(ksp,&pc);
117  if (solver_type==MUMPSPACKAGE_LU){
118  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=9)
119  PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
120  #else
121  PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
122  #endif
123  }
124 
125  /*FS: */
126  if (solver_type==FSSolverEnum){
127  /*Make indices out of doftypes: */
128  if(!df)_error_("need doftypes for FS solver!\n");
129  DofTypesToIndexSet(&isv,&isp,df,FSSolverEnum);
130 
131  /*Set field splits: */
132  KSPGetPC(ksp,&pc);
133  #if _PETSC_MINOR_==1
134  PCFieldSplitSetIS(pc,isv);
135  PCFieldSplitSetIS(pc,isp);
136  #else
137  PCFieldSplitSetIS(pc,PETSC_NULL,isv);
138  PCFieldSplitSetIS(pc,PETSC_NULL,isp);
139  #endif
140 
141  }
142 
143  /*If there is an initial guess for the solution, use it
144  * except if we are using the MUMPS direct solver
145  * where any initial guess will crash Petsc*/
146  if (uf0){
147  if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
148  KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
149  }
150  }
151 
152  /*Solve: */
153  if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
154  KSPSolve(ksp,pf,uf);
155 
156  /*Check convergence*/
157  KSPGetIterationNumber(ksp,&iteration_number);
158  if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
159  if (VerboseSolver()) _printf0_("Petsc: "<< iteration_number << " KSP iterations\n");
160 
161  /*Free resources:*/
162  KSPFree(&ksp);
163 
164  /*Assign output pointers:*/
165  *puf=uf;
166 }

◆ DofTypesToIndexSet()

void DofTypesToIndexSet ( IS *  pisv,
IS *  pisp,
Vec  df,
int  typeenum 
)

Definition at line 168 of file PetscSolver.cpp.

168  { /*{{{*/
169 
170  /*output: */
171  IS isv=NULL;
172  IS isp=NULL;
173 
174  int start,end;
175  IssmDouble* df_local=NULL;
176  int df_local_size;
177 
178  int* pressure_indices=NULL;
179  int* velocity_indices=NULL;
180  int pressure_num=0;
181  int velocity_num=0;
182  int pressure_count=0;
183  int velocity_count=0;
184 
185  if(typeenum==FSSolverEnum){
186 
187  /*Ok, recover doftypes vector values and indices: */
188  VecGetOwnershipRange(df,&start,&end);
189  VecGetLocalSize(df,&df_local_size);
190  VecGetArray(df,&df_local);
191 
192  pressure_num=0;
193  velocity_num=0;
194  for(int i=0;i<df_local_size;i++){
195  if (df_local[i]==PressureEnum)pressure_num++;
196  else velocity_num++;
197  }
198 
199  /*Allocate indices: */
200  if(pressure_num)pressure_indices=xNew<int>(pressure_num);
201  if(velocity_num)velocity_indices=xNew<int>(velocity_num);
202 
203  pressure_count=0;
204  velocity_count=0;
205  for(int i=0;i<df_local_size;i++){
206  if (df_local[i]==PressureEnum){
207  pressure_indices[pressure_count]=start+i;
208  pressure_count++;
209  }
210  if (df_local[i]==VelocityEnum){
211  velocity_indices[velocity_count]=start+i;
212  velocity_count++;
213  }
214  }
215  VecRestoreArray(df,&df_local);
216 
217  /*Create indices sets: */
218  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
219  ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
220  ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
221  #else
222  ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
223  ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
224  #endif
225  }
226 
227  /*Free ressources:*/
228  xDelete<int>(pressure_indices);
229  xDelete<int>(velocity_indices);
230 
231  /*Assign output pointers:*/
232  *pisv=isv;
233  *pisp=isp;
234 }
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
SPOOLESPACKAGE_CHOL
@ SPOOLESPACKAGE_CHOL
Definition: SolverEnum.h:13
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
PetscOptionsDetermineSolverType
void PetscOptionsDetermineSolverType(int *psolver_type)
Definition: PetscOptionsDetermineSolverType.cpp:20
SUPERLUDISTPACKAGE
@ SUPERLUDISTPACKAGE
Definition: SolverEnum.h:14
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
SPOOLESPACKAGE_LU
@ SPOOLESPACKAGE_LU
Definition: SolverEnum.h:12
KSPFree
void KSPFree(KSP *pksp)
Definition: KSPFree.cpp:16
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
FSSolverEnum
@ FSSolverEnum
Definition: EnumDefinitions.h:1061
VelocityEnum
@ VelocityEnum
Definition: EnumDefinitions.h:458
PetscVec::vector
Vec vector
Definition: PetscVec.h:25
SolverxPetsc
void SolverxPetsc(Vec *puf, Mat Kff, Vec pf, Vec uf0, Vec df, Parameters *parameters)
Definition: PetscSolver.cpp:35
NewVec
Vec NewVec(int size, ISSM_MPI_Comm comm, bool fromlocalsize)
Definition: NewVec.cpp:19
DofTypesToIndexSet
void DofTypesToIndexSet(IS *pisv, IS *pisp, Vec df, int typeenum)
Definition: PetscSolver.cpp:168
VerboseSolver
bool VerboseSolver(void)
Definition: Verbosity.cpp:25
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
PetscVec
Definition: PetscVec.h:22
PetscMat::matrix
Mat matrix
Definition: PetscMat.h:27
MUMPSPACKAGE_CHOL
@ MUMPSPACKAGE_CHOL
Definition: SolverEnum.h:11
MUMPSPACKAGE_LU
@ MUMPSPACKAGE_LU
Definition: SolverEnum.h:10