source: issm/trunk/src/c/modules/Solverx/SolverxPetsc.cpp@ 11995

Last change on this file since 11995 was 11995, checked in by Mathieu Morlighem, 13 years ago

merged trunk-jpl and trunk for revision 11994M
File size: 5.0 KB
Line 
1/*!\file SolverxPetsc
2 * \brief Petsc implementation of solver
3 */
4
5#include "./Solverx.h"
6#include "../../shared/shared.h"
7#include "../../include/include.h"
8#include "../../io/io.h"
9
10#ifdef HAVE_CONFIG_H
11 #include <config.h>
12#else
13#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
14#endif
15
16void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
17
18 /*Output: */
19 Vec uf = NULL;
20
21 /*Intermediary: */
22 int local_m,local_n,global_m,global_n;
23 int analysis_type;
24
25 /*Solver */
26 KSP ksp = NULL;
27 PC pc = NULL;
28 int iteration_number;
29 int solver_type;
30 bool fromlocalsize = true;
31 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
32 PetscTruth flag,flg;
33 #else
34 PetscBool flag,flg;
35 #endif
36
37 /*Stokes: */
38 IS isv=NULL;
39 IS isp=NULL;
40
41 #if _PETSC_MAJOR_ >= 3
42 char ksp_type[50];
43 #endif
44
45
46 /*Display message*/
47 _printf_(VerboseModule()," Solving\n");
48 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
49 if(VerboseSolver())PetscOptionsPrint(stdout);
50 #else
51 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
52 #endif
53
54 /*First, check that f-set is not NULL, i.e. model is fully constrained: {{{*/
55 _assert_(Kff);
56 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
57 if(!global_n){
58 *puf=NULL; return;
59 }
60 /*}}}*/
61 /*Initial guess logic here: {{{1*/
62 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
63 #if _PETSC_MAJOR_ >= 3
64 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
65 if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
66 #endif
67
68 /*If initial guess for the solution exists, use it to create uf, otherwise,
69 * duplicate the right hand side so that the solution vector has the same structure*/
70 if(uf0){
71 VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
72 }
73 else{
74 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
75 }
76 /*}}}*/
77 /*Process petsc options to see if we are using special types of external solvers: {{{1*/
78 PetscOptionsDetermineSolverType(&solver_type);
79
80 /*In serial mode, the matrices have been loaded as MPIAIJ or AIJ matrices.
81 We need to convert them if we are going to run the solvers successfully: */
82 #ifdef _SERIAL_
83 #if _PETSC_MAJOR_ == 2
84 if (solver_type==MUMPSPACKAGE_LU){
85 /*Convert Kff to MATTAIJMUMPS: */
86 MatConvert(Kff,MATAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
87 }
88 if (solver_type==MUMPSPACKAGE_CHOL){
89 /*Convert Kff to MATTSBAIJMUMPS: */
90 MatConvert(Kff,MATSBAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
91 }
92 if (solver_type==SPOOLESPACKAGE_LU){
93 /*Convert Kff to MATTSBAIJMUMPS: */
94 MatConvert(Kff,MATAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
95 }
96 if (solver_type==SPOOLESPACKAGE_CHOL){
97 /*Convert Kff to MATTSBAIJMUMPS: */
98 MatConvert(Kff,MATSBAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
99 }
100 if (solver_type==SUPERLUDISTPACKAGE){
101 /*Convert Kff to MATTSBAIJMUMPS: */
102 MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff);
103 }
104 if (solver_type==StokesSolverEnum){
105 _error_("Petsc 2 does not support multi-physics solvers");
106 }
107 #endif
108 #endif
109 /*}}}*/
110 /*Check the solver is available: {{{1*/
111 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
112 #if _PETSC_MAJOR_ >=3
113 #ifndef _HAVE_MUMPS_
114 _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
115 #endif
116
117 #endif
118 }
119 /*}}}*/
120 /*Prepare solver:{{{1*/
121 KSPCreate(MPI_COMM_WORLD,&ksp);
122 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
123 KSPSetFromOptions(ksp);
124
125 #if defined(_SERIAL_) && _PETSC_MAJOR_==3
126 /*Specific solver?: */
127 KSPGetPC(ksp,&pc);
128 if (solver_type==MUMPSPACKAGE_LU){
129 #if _PETSC_MINOR_==1
130 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
131 #else
132 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
133 #endif
134 }
135 #endif
136
137 #if defined(_PARALLEL_) && _PETSC_MAJOR_==3
138 /*Stokes: */
139 if (solver_type==StokesSolverEnum){
140 /*Make indices out of doftypes: */
141 if(!df)_error_("need doftypes for Stokes solver!\n");
142 DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
143
144 /*Set field splits: */
145 KSPGetPC(ksp,&pc);
146 #if _PETSC_MINOR_==1
147 PCFieldSplitSetIS(pc,isv);
148 PCFieldSplitSetIS(pc,isp);
149 #else
150 PCFieldSplitSetIS(pc,PETSC_NULL,isv);
151 PCFieldSplitSetIS(pc,PETSC_NULL,isp);
152 #endif
153
154 }
155 #endif
156
157 /*}}}*/
158 /*If there is an initial guess for the solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/
159 if (uf0){
160 if( (solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
161 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
162 }
163 }
164 /*}}}*/
165
166 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
167
168 /*Solve: */
169 KSPSolve(ksp,pf,uf);
170
171 /*Check convergence*/
172 KSPGetIterationNumber(ksp,&iteration_number);
173 if (iteration_number<0) _error_("%s%i"," Solver diverged at iteration number: ",-iteration_number);
174
175 /*Free resources:*/
176 KSPFree(&ksp);
177
178 /*Assign output pointers:*/
179 *puf=uf;
180}
Note: See TracBrowser for help on using the repository browser.