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

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 13393

File size: 4.2 KB
RevLine 
[11679]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
[13395]16void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){
17
18 PetscVec* uf=new PetscVec();
19
20 Vec uf0_vector = NULL;
21 Vec df_vector = NULL;
22
23 if(uf0) uf0_vector = uf0->vector;
24 if(df) df_vector = df->vector;
25
26 SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
27
28 *puf=uf;
29
30}
[11679]31void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
32
33 /*Output: */
34 Vec uf = NULL;
35
36 /*Intermediary: */
37 int local_m,local_n,global_m,global_n;
38 int analysis_type;
39
40 /*Solver */
41 KSP ksp = NULL;
42 PC pc = NULL;
43 int iteration_number;
44 int solver_type;
45 bool fromlocalsize = true;
46 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
47 PetscTruth flag,flg;
48 #else
49 PetscBool flag,flg;
50 #endif
51
52 /*Stokes: */
53 IS isv=NULL;
54 IS isp=NULL;
55
56 #if _PETSC_MAJOR_ >= 3
57 char ksp_type[50];
58 #endif
59
60 /*Display message*/
[12706]61 if(VerboseModule()) _pprintLine_(" Solving");
[11679]62 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
63 if(VerboseSolver())PetscOptionsPrint(stdout);
64 #else
65 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
66 #endif
67
[12330]68 /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
[11679]69 _assert_(Kff);
70 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
71 if(!global_n){
[13395]72 *puf=NewVec(0); return;
[11679]73 }
[12330]74
75 /*Initial guess */
[11679]76 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
77 #if _PETSC_MAJOR_ >= 3
78 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
79 if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
80 #endif
81
82 /*If initial guess for the solution exists, use it to create uf, otherwise,
83 * duplicate the right hand side so that the solution vector has the same structure*/
84 if(uf0){
85 VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
86 }
87 else{
88 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
89 }
[12330]90
91 /*Process petsc options to see if we are using special types of external solvers*/
[11679]92 PetscOptionsDetermineSolverType(&solver_type);
93
[12330]94 /*Check the solver is available*/
[11679]95 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
[12330]96 #if _PETSC_MAJOR_ >=3
97 #ifndef _HAVE_MUMPS_
[13395]98 _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
[12330]99 #endif
[11679]100 #endif
[12330]101 }
[11679]102
[12330]103 /*Prepare solver*/
[11679]104 KSPCreate(MPI_COMM_WORLD,&ksp);
105 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
106 KSPSetFromOptions(ksp);
107
[12330]108 #if _PETSC_MAJOR_==3
[11679]109 /*Specific solver?: */
110 KSPGetPC(ksp,&pc);
111 if (solver_type==MUMPSPACKAGE_LU){
112 #if _PETSC_MINOR_==1
113 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
114 #else
115 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
116 #endif
117 }
118
119 /*Stokes: */
120 if (solver_type==StokesSolverEnum){
121 /*Make indices out of doftypes: */
[13395]122 if(!df)_error_("need doftypes for Stokes solver!\n");
[11679]123 DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
124
125 /*Set field splits: */
126 KSPGetPC(ksp,&pc);
127 #if _PETSC_MINOR_==1
128 PCFieldSplitSetIS(pc,isv);
129 PCFieldSplitSetIS(pc,isp);
130 #else
131 PCFieldSplitSetIS(pc,PETSC_NULL,isv);
132 PCFieldSplitSetIS(pc,PETSC_NULL,isp);
133 #endif
134
135 }
136 #endif
137
[12330]138 /*If there is an initial guess for the solution, use it
139 * except if we are using the MUMPS direct solver
140 * where any initial guess will crash Petsc*/
[11679]141 if (uf0){
[12330]142 if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
[11679]143 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
144 }
145 }
146
147 /*Solve: */
[12330]148 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
[11679]149 KSPSolve(ksp,pf,uf);
150
151 /*Check convergence*/
152 KSPGetIterationNumber(ksp,&iteration_number);
[13395]153 if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
[11679]154
155 /*Free resources:*/
156 KSPFree(&ksp);
157
158 /*Assign output pointers:*/
159 *puf=uf;
160}
Note: See TracBrowser for help on using the repository browser.