SFEMaNS  version 4.1 (work in progress) Reference documentation for SFEMaNS
solver.f90
Go to the documentation of this file.
1 MODULE solve_petsc
2 use my_util
4  INTEGER:: it_max
5  REAL(KIND=8) :: rel_tol, abs_tol
6  CHARACTER(LEN=20) :: solver, precond
7  LOGICAL :: verbose
8  END TYPE solver_param
9
10 CONTAINS
11  SUBROUTINE init_solver(my_par,my_ksp,matrix,communicator,solver,precond, opt_re_init)
13  IMPLICIT NONE
14  LOGICAl, INTENT(IN), OPTIONAL :: opt_re_init
15  type(solver_param) :: my_par
16  CHARACTER(*), OPTIONAL :: solver, precond
17  LOGICAL :: re_init
18  INTEGER :: deb, fin
19 #include "petsc/finclude/petsc.h"
20  mat :: matrix
21  ksp :: my_ksp
22  pc :: prec
23  petscerrorcode :: ierr
24  mpi_comm :: communicator
25
26  IF (.NOT.present(opt_re_init)) THEN
27  re_init=.false.
28  ELSE
29  re_init=opt_re_init
30  END IF
31
32  IF (my_par%it_max.LE.0) THEN
33  my_par%it_max = 100
34  END IF
35  IF (my_par%rel_tol.LE.0.d0) THEN
36  my_par%rel_tol = 1.d-8
37  END IF
38  IF (my_par%abs_tol.LE.0.d0) THEN
39  my_par%abs_tol = 1.d-14
40  END IF
41
42  IF (.NOT.re_init) CALL kspcreate(communicator,my_ksp,ierr)
43  !CALL KSPCreate(communicator,my_ksp,ierr)
44
45  !CALL KSPSetOperators(my_ksp,matrix,matrix,DIFFERENT_NONZERO_PATTERN,ierr)
46  !CALL KSPSetOperators(my_ksp,matrix,matrix,SAME_NONZERO_PATTERN ,ierr) !Petsc 3.4.2
47  CALL kspsetoperators(my_ksp,matrix,matrix,ierr) !Petsc 3.7.2
48
49  IF (present(solver)) THEN
50  deb = start_of_string(solver)
51  fin = last_of_string(solver)
52  IF (solver(deb:fin)=='BCGS') THEN
53  CALL kspsettype(my_ksp, kspbcgs, ierr)
54  ELSE IF (solver(deb:fin)=='GMRES') THEN
55  CALL kspsettype(my_ksp, kspgmres, ierr)
56  ELSE IF (solver(deb:fin)=='FGMRES') THEN
57  CALL kspsettype(my_ksp, kspfgmres, ierr)
58  ELSE IF (solver(deb:fin)=='PCR') THEN
59  CALL kspsettype(my_ksp, kspcr, ierr)
60  ELSE IF (solver(deb:fin)=='CHEBYCHEV') THEN
61  CALL kspsettype(my_ksp, kspchebyshev, ierr)
62  ELSE IF (solver(deb:fin)=='CG') THEN
63  CALL kspsettype(my_ksp, kspcg, ierr)
64  ELSE
65  CALL kspsettype(my_ksp, kspfgmres, ierr)
66  END IF
67  ELSE
68  CALL kspsettype(my_ksp, kspfgmres, ierr)
69  END IF
70  CALL kspsettolerances(my_ksp, my_par%rel_tol, my_par%abs_tol, &
71  petsc_default_real, my_par%it_max, ierr)
72  CALL kspgetpc(my_ksp, prec, ierr)
73  IF (present(precond)) THEN
74  deb = start_of_string(precond)
75  fin = last_of_string(precond)
76  IF (precond(deb:fin)=='JACOBI') THEN
77  call pcsettype(prec, pcbjacobi, ierr)
78  ELSE IF (precond(deb:fin)=='HYPRE') THEN
79  CALL pcsettype(prec, pchypre, ierr)
80  ELSE IF (precond(deb:fin)=='SSOR') THEN
81  CALL pcsettype(prec, pcsor, ierr)
82  ELSE IF (precond(deb:fin)=='MUMPS') THEN
83  CALL pcsettype(prec, pclu, ierr)
84  CALL kspsettype(my_ksp, ksppreonly, ierr)
85  CALL pcfactorsetmatsolverpackage(prec, matsolvermumps, ierr)
86  ELSE
87  CALL pcsettype(prec, pchypre, ierr)
88  END IF
89  ELSE
90  CALL pcsettype(prec, pchypre, ierr)
91  END IF
92  CALL kspsetfromoptions(my_ksp, ierr)
93  END SUBROUTINE init_solver
94
95  SUBROUTINE solver(my_ksp,b,x,reinit,verbose)
96  IMPLICIT NONE
97  LOGICAL, OPTIONAL :: reinit, verbose
98  INTEGER :: its
99 #include "petsc/finclude/petsc.h"
100  ksp :: my_ksp
101  petscerrorcode :: ierr
102  vec :: x, b
103  kspconvergedreason :: reason
104  IF (.NOT.present(reinit)) reinit=.true.
105  IF (.NOT.present(verbose)) verbose=.false.
106
107  IF (reinit) CALL veczeroentries(x,ierr)
108  CALL kspsolve(my_ksp,b,x,ierr)
109  IF (verbose) THEN
110  CALL kspgetiterationnumber(my_ksp, its, ierr)
111  CALL kspgetconvergedreason(my_ksp, reason, ierr)
112  SELECT case(reason)
113  case(2)
114  WRITE(*,*) "KSP_CONVERGED_RTOL, Nb of iterations", its
115  case(3)
116  WRITE(*,*) "KSP_CONVERGED_ATOL, Nb of iterations", its
117  case(4)
118  WRITE(*,*) "Converged after one single iteration of the preconditioner is applied"
119  case(5,6,7,8)
120  WRITE(*,*) "Converge for strange reason:", reason
121  case(-2)
122  WRITE(*,*) "KSP_DIVERGED_NULL"
123  case(-3)
124  WRITE(*,*) "Not converged after it_max", its
125  case(-4)
126  WRITE(*,*) "Not converged: explosion"
127  case(-5,-6,-7)
128  WRITE(*,*) "Not converged for strange reasons", reason
129  case(-8)
130  WRITE(*,*) "Not converged: Indefinite preconditioner"
131  case(-9)
132  WRITE(*,*) "Not converged: NAN"
133  case(-10)
134  WRITE(*,*) "Not converged: Indefinite matrix"
135  case default
136  WRITE(*,*) "Something strange happened", reason
137  END SELECT
138  END IF
139
140  END SUBROUTINE solver
141
142  SUBROUTINE create_local_petsc_matrix(communicator, LA, matrix, clean)
143  USE def_type_mesh
144  IMPLICIT NONE
145  TYPE(petsc_csr_la) :: la
146  LOGICAL, OPTIONAL :: clean
147  REAL(KIND=8),DIMENSION(:), POINTER :: aa
148  INTEGER :: nnzm1, dom_np
149  LOGICAL :: test_clean
150
151 !!$INTEGER, DIMENSION(:), POINTER :: ia, ja 152 153 #include "petsc/finclude/petsc.h" 154 mpi_comm :: communicator 155 mat :: matrix 156 petscerrorcode :: ierr 157 !------------------------------------------------------------------------------ 158 dom_np = SIZE(la%ia)-1 159 nnzm1=la%ia(dom_np)-la%ia(0)-1 160 ALLOCATE(aa(0:nnzm1)) 161 aa =0.d0 162 163 164 !!$ ALLOCATE(ia(0:dom_np),ja(0:nnzm1))
165 !!$ia = LA%ia 166 !!$ ja = LA%ja
167 !!$CALL MatCreateMPIAIJWithArrays(communicator,dom_np,dom_np,PETSC_DECIDE, & 168 !!$ PETSC_DECIDE, ia, ja, aa, matrix, ierr)
169 !!\$DEALLOCATE(ia,ja)
170
171  CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
172  petsc_decide, la%ia, la%ja, aa, matrix, ierr)
173
174  DEALLOCATE(aa)
175  IF (present(clean)) THEN
176  test_clean=clean
177  ELSE
178  test_clean=.true.
179  END IF
180  IF (test_clean) THEN
181  IF (ASSOCIATED(la%ia)) DEALLOCATE(la%ia)
182  IF (ASSOCIATED(la%ja)) DEALLOCATE(la%ja)
183  END IF
184  END SUBROUTINE create_local_petsc_matrix
185
186  SUBROUTINE create_local_petsc_matrix_a_detruire(communicator, aij, i_loc, matrix)
187  USE def_type_mesh
188  IMPLICIT NONE
189  TYPE(aij_type), INTENT(IN) :: aij
190  INTEGER, DIMENSION(2) :: i_loc
191  INTEGER, DIMENSION(:), POINTER :: ia, ja
192  REAL(KIND=8),DIMENSION(:), POINTER :: aa
193  INTEGER :: nnzm1, dom_np, p, i, n
194 #include "petsc/finclude/petsc.h"
195  mpi_comm :: communicator
196  mat :: matrix
197  petscerrorcode :: ierr
198  !------------------------------------------------------------------------------
199  dom_np = i_loc(2) - i_loc(1) + 1
200  nnzm1 = aij%ia(i_loc(2)+1)-aij%ia(i_loc(1))-1
201  ALLOCATE(ia(0:dom_np),ja(0:nnzm1))
202  ia(0)=0
203  DO i = 1, dom_np
204  n = i_loc(1) + i -1
205  ia(i) = aij%ia(n+1)-aij%ia(i_loc(1))
206  DO p=aij%ia(n), aij%ia(n+1)-1
207  ja(p-aij%ia(i_loc(1)))= aij%ja(p)-1
208  end do
209  END DO
210  !------------------------------------------------------------------------------
211  ALLOCATE(aa(0:nnzm1))
212  aa =0
213  CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
214  petsc_decide, ia, ja, aa, matrix, ierr)
215
216  DEALLOCATE(ia,ja,aa)
218
219  SUBROUTINE create_local_petsc_block_matrix(communicator, n_b, aij, i_loc, matrix)
220  USE def_type_mesh
221  IMPLICIT NONE
222  TYPE(aij_type), INTENT(IN) :: aij
223  INTEGER, DIMENSION(2) :: i_loc
224  INTEGER :: n_b
225  INTEGER, DIMENSION(:), POINTER :: ia, ja
226  REAL(KIND=8),DIMENSION(:), POINTER :: aa
227  INTEGER :: nnzm1, dom_np, p, i, n, ib, k
228 #include "petsc/finclude/petsc.h"
229  mpi_comm :: communicator
230  mat :: matrix
231  petscerrorcode :: ierr
232
233  dom_np = i_loc(2) - i_loc(1) + 1
234  nnzm1 = n_b*(aij%ia(i_loc(2)+1)-aij%ia(i_loc(1))-1)
235  ALLOCATE(ia(0:n_b*dom_np),ja(0:nnzm1))
236  ia(0)=0
237  DO k = 1, n_b
238  DO i = 1, dom_np
239  ib = i + (k-1)*dom_np
240  n = i_loc(1) + i - 1
241  ia(i) = n_b*(aij%ia(n+1)-aij%ia(i_loc(1)))
242  DO p=aij%ia(n), aij%ia(n+1)-1
243  ja(p-aij%ia(i_loc(1)))= aij%ja(p)-1
244  end do
245  END DO
246  END DO
247  !------------------------------------------------------------------------------
248  ALLOCATE(aa(0:nnzm1))
249  aa =0
250  CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
251  petsc_decide, ia, ja, aa, matrix, ierr)
252
253  DEALLOCATE(ia,ja,aa)
254  END SUBROUTINE create_local_petsc_block_matrix
255
256
257 END MODULE solve_petsc
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
integer function, public start_of_string(string)
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
integer function, public last_of_string(string)
subroutine create_local_petsc_block_matrix(communicator, n_b, aij, i_loc, matrix)
Definition: solver.f90:219
subroutine create_local_petsc_matrix_a_detruire(communicator, aij, i_loc, matrix)
Definition: solver.f90:186
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:142