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

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

merged trunk-jpl and trunk for revision 12703

File size: 3.9 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 /*Display message*/
46 if(VerboseModule()) _pprintLine_(" Solving");
47 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
48 if(VerboseSolver())PetscOptionsPrint(stdout);
49 #else
50 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
51 #endif
52
53 /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
54 _assert_(Kff);
55 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
56 if(!global_n){
57 *puf=NULL; return;
58 }
59
60 /*Initial guess */
61 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
62 #if _PETSC_MAJOR_ >= 3
63 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
64 if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
65 #endif
66
67 /*If initial guess for the solution exists, use it to create uf, otherwise,
68 * duplicate the right hand side so that the solution vector has the same structure*/
69 if(uf0){
70 VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
71 }
72 else{
73 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
74 }
75
76 /*Process petsc options to see if we are using special types of external solvers*/
77 PetscOptionsDetermineSolverType(&solver_type);
78
79 /*Check the solver is available*/
80 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
81 #if _PETSC_MAJOR_ >=3
82 #ifndef _HAVE_MUMPS_
83 _error2_("requested MUMPS solver, which was not compiled into ISSM!\n");
84 #endif
85 #endif
86 }
87
88 /*Prepare solver*/
89 KSPCreate(MPI_COMM_WORLD,&ksp);
90 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
91 KSPSetFromOptions(ksp);
92
93 #if _PETSC_MAJOR_==3
94 /*Specific solver?: */
95 KSPGetPC(ksp,&pc);
96 if (solver_type==MUMPSPACKAGE_LU){
97 #if _PETSC_MINOR_==1
98 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
99 #else
100 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
101 #endif
102 }
103
104 /*Stokes: */
105 if (solver_type==StokesSolverEnum){
106 /*Make indices out of doftypes: */
107 if(!df)_error2_("need doftypes for Stokes solver!\n");
108 DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
109
110 /*Set field splits: */
111 KSPGetPC(ksp,&pc);
112 #if _PETSC_MINOR_==1
113 PCFieldSplitSetIS(pc,isv);
114 PCFieldSplitSetIS(pc,isp);
115 #else
116 PCFieldSplitSetIS(pc,PETSC_NULL,isv);
117 PCFieldSplitSetIS(pc,PETSC_NULL,isp);
118 #endif
119
120 }
121 #endif
122
123 /*If there is an initial guess for the solution, use it
124 * except if we are using the MUMPS direct solver
125 * where any initial guess will crash Petsc*/
126 if (uf0){
127 if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
128 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
129 }
130 }
131
132 /*Solve: */
133 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
134 KSPSolve(ksp,pf,uf);
135
136 /*Check convergence*/
137 KSPGetIterationNumber(ksp,&iteration_number);
138 if (iteration_number<0) _error2_("Solver diverged at iteration number: " << -iteration_number);
139
140 /*Free resources:*/
141 KSPFree(&ksp);
142
143 /*Assign output pointers:*/
144 *puf=uf;
145}
Note: See TracBrowser for help on using the repository browser.