SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
fem_M_axi.f90
Go to the documentation of this file.
1 !
2 !
3 !Authors: Jean-Luc Guermond, Raphael Laguerre, Luigi Quartapelle, Copyrights 1996, 2000, 2004
4 !Revised, Jean-Luc Guermond, Franky Luddens, January, 2011
5 !
6 MODULE fem_m_axi
7  USE def_type_mesh
8 CONTAINS
9 
10 
11  SUBROUTINE qs_diff_mass_vect_m (type_op, LA, mesh, visco, mass, stab, mode, matrix)
12  !=================================================
13  USE my_util
14  IMPLICIT NONE
15  TYPE(mesh_type), INTENT(IN) :: mesh
16  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
17  INTEGER, INTENT(IN) :: type_op, mode
18  TYPE(petsc_csr_la) :: la
19  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
20  INTEGER, DIMENSION(:), ALLOCATABLE :: idxm, idxn
21  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
22  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc, b_loc
23  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat_loc
24  REAL(KIND=8) :: ray, eps1, eps2
25  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, k_max, n_w, ix, jx, ki, kj
26  REAL(KIND=8) :: viscolm, xij
27 #include "petsc/finclude/petsc.h"
28  mat :: matrix
29  petscerrorcode :: ierr
30  CALL matzeroentries(matrix, ierr)
31  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
32 
33 
34  IF (type_op == 1) THEN
35  eps1 = 1.d0
36  eps2 = 1.d0
37  k_max = 2 ! 2x2 Structure
38  ELSEIF (type_op == 2) THEN
39  eps1 = 1.d0
40  eps2 = -1.d0
41  k_max = 2 !2x2 Structure
42  ELSEIF (type_op == 3) THEN
43  eps1 = 0.d0
44  eps2 = 0.d0
45  k_max = 1 ! Scalar Structure
46  ELSE
47  eps1 = 0.d0
48  eps2 = 0.d0
49  k_max = 1
50  CALL error_petsc('BUG in qs_diff_mass_vect_M, type_op not correct')
51  ENDIF
52 
53  IF (k_max/=SIZE(la%loc_to_glob,1)) THEN
54  CALL error_petsc('BUG in qs_diff_mass_vect_petsc_M, k_max/=SIZE(LA%loc_to_glob,1)')
55  END IF
56 
57  n_w = mesh%gauss%n_w
58  DO l = 1, mesh%gauss%l_G
59  DO ni = 1, n_w
60  DO nj = 1, n_w
61  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
62  END DO
63  END DO
64  END DO
65 
66  ALLOCATE(mat_loc(k_max*n_w,k_max*n_w),idxm(k_max*n_w),idxn(k_max*n_w))
67  DO m = 1, mesh%dom_me
68  jj_loc = mesh%jj(:,m)
69 
70  a_loc = 0.d0
71  b_loc = 0.d0
72  DO l = 1, mesh%gauss%l_G
73  !Compute radius of Gauss point
74  ray = 0
75  DO ni = 1, n_w; i = jj_loc(ni)
76  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
77  END DO
78  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
79  DO nj = 1, n_w
80  DO ni = 1, n_w
81 
82  !grad(u).grad(v) w.r.t. r and z
83  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
84 
85  !start diagonal block
86  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
87  + viscolm*eps1*wwprod(ni,nj,l)/ray &
88  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
89  + viscolm*mode**2*wwprod(ni,nj,l)/ray
90  !end diagonal block
91  b_loc(ni,nj) = b_loc(ni,nj) + eps2*viscolm*2*mode*wwprod(ni,nj,l)/ray
92  END DO
93  END DO
94  END DO
95 
96  mat_loc = 0.d0
97  DO ki= 1, k_max
98  DO ni = 1, n_w
99  i = jj_loc(ni)
100  iglob = la%loc_to_glob(ki,i)
101  ix = (ki-1)*n_w+ni
102  idxm(ix) = iglob-1
103  DO kj = 1, k_max
104  DO nj = 1, n_w
105  j = jj_loc(nj)
106  jglob = la%loc_to_glob(kj,j)
107  jx = (kj-1)*n_w+nj
108  idxn(jx) = jglob-1
109  IF (ki==kj) THEN
110  mat_loc(ix,jx) = mat_loc(ix,jx) + a_loc(ni,nj)
111  ELSE
112  mat_loc(ix,jx) = mat_loc(ix,jx) + b_loc(ni,nj)
113  END IF
114  END DO
115  END DO
116  END DO
117  END DO
118  CALL matsetvalues(matrix, k_max*n_w, idxm, k_max*n_w, idxn, mat_loc, add_values, ierr)
119  ENDDO
120  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
121  CALL matassemblyend(matrix,mat_final_assembly,ierr)
122 
123  DEALLOCATE(mat_loc,idxm,idxn)
124 
125  END SUBROUTINE qs_diff_mass_vect_m
126 
127  SUBROUTINE qs_diff_mass_scal_m (mesh, LA, visco, mass, stab, mode, matrix)
128  !=================================================
129  IMPLICIT NONE
130  TYPE(mesh_type), INTENT(IN) :: mesh
131  TYPE(petsc_csr_la) :: la
132  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
133  INTEGER, INTENT(IN) :: mode
134  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
135  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
136  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
137  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
138  REAL(KIND=8) :: ray
139  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
140  REAL(KIND=8) :: viscolm, xij
141 #include "petsc/finclude/petsc.h"
142  mat :: matrix
143  petscerrorcode :: ierr
144  CALL matzeroentries(matrix, ierr)
145  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
146 
147  DO l = 1, mesh%gauss%l_G
148  DO ni = 1, mesh%gauss%n_w
149  DO nj = 1, mesh%gauss%n_w
150  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
151  END DO
152  END DO
153  END DO
154 
155  n_w = mesh%gauss%n_w
156  DO m = 1, mesh%dom_me
157  jj_loc = mesh%jj(:,m)
158  a_loc = 0.d0
159 
160  DO nj = 1, mesh%gauss%n_w;
161  j = jj_loc(nj)
162  jglob = la%loc_to_glob(1,j)
163  idxn(nj) = jglob-1
164  DO ni = 1, mesh%gauss%n_w;
165  i = jj_loc(ni)
166  iglob = la%loc_to_glob(1,i)
167  idxm(ni) = iglob-1
168  END DO
169  END DO
170 
171  DO l = 1, mesh%gauss%l_G
172  !Compute radius of Gauss point
173  ray = 0
174  DO ni = 1, n_w; i = jj_loc(ni)
175  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
176  END DO
177  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
178 
179  DO nj = 1, n_w
180  DO ni = 1, n_w
181  !grad(u).grad(v) w.r.t. r and z
182  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
183  !start diagonal block
184  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
185  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
186  + viscolm*mode**2*wwprod(ni,nj,l)/ray
187  !end diagonal block
188  END DO
189  END DO
190  END DO
191  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
192  ENDDO
193  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
194  CALL matassemblyend(matrix,mat_final_assembly,ierr)
195 
196  END SUBROUTINE qs_diff_mass_scal_m
197 
198 SUBROUTINE qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, stab, mode, matrix)
199  !=================================================
200  IMPLICIT NONE
201  TYPE(mesh_type), INTENT(IN) :: mesh
202  TYPE(petsc_csr_la) :: la
203  REAL(KIND=8), INTENT(IN) :: mass, stab
204  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: heat_capa, visco
205  INTEGER, INTENT(IN) :: mode
206  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
207  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
208  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
209  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
210  REAL(KIND=8) :: ray
211  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
212  REAL(KIND=8) :: viscolm, xij
213 #include "petsc/finclude/petsc.h"
214  mat :: matrix
215  petscerrorcode :: ierr
216  CALL matzeroentries(matrix, ierr)
217  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
218 
219  DO l = 1, mesh%gauss%l_G
220  DO ni = 1, mesh%gauss%n_w
221  DO nj = 1, mesh%gauss%n_w
222  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
223  END DO
224  END DO
225  END DO
226 
227  n_w = mesh%gauss%n_w
228  DO m = 1, mesh%dom_me
229  jj_loc = mesh%jj(:,m)
230  a_loc = 0.d0
231 
232  DO nj = 1, mesh%gauss%n_w;
233  j = jj_loc(nj)
234  jglob = la%loc_to_glob(1,j)
235  idxn(nj) = jglob-1
236  DO ni = 1, mesh%gauss%n_w;
237  i = jj_loc(ni)
238  iglob = la%loc_to_glob(1,i)
239  idxm(ni) = iglob-1
240  END DO
241  END DO
242 
243  DO l = 1, mesh%gauss%l_G
244  !Compute radius of Gauss point
245  ray = 0
246  DO ni = 1, n_w; i = jj_loc(ni)
247  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
248  END DO
249  viscolm = (visco(m) + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
250  DO nj = 1, n_w
251  DO ni = 1, n_w
252  !grad(u).grad(v) w.r.t. r and z
253  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
254  !start diagonal block
255  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
256  + heat_capa(m) * mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
257  + viscolm*mode**2*wwprod(ni,nj,l)/ray
258  !end diagonal block
259  END DO
260  END DO
261  END DO
262  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
263  ENDDO
264  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
265  CALL matassemblyend(matrix,mat_final_assembly,ierr)
266 
267  END SUBROUTINE qs_diff_mass_scal_m_variant
268 
269  SUBROUTINE qs_mass_vect_3x3(LA, mesh, mass, matrix)
270  !========================================================================
271  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
272  ! (V1,V4,V5), (V2,V3,V6)
273  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
274  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
275  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
276  ! 2 pour les composantes 2, 3 et 6
277  !------------------------------------------------------------------------
278  USE my_util
279  USE input_data
280  IMPLICIT NONE
281 
282  REAL(KIND=8), INTENT(IN) :: mass
283  TYPE(petsc_csr_la) :: la
284  TYPE(mesh_type) :: mesh
285  INTEGER ::l, m, ni, nj, i, j, ki, kj, n_w
286  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
287  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
288  REAL(KIND=8) :: ray
289  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
290  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
291  INTEGER :: ix, jx, iglob, jglob
292  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
293 #include "petsc/finclude/petsc.h"
294  mat :: matrix
295  petscerrorcode :: ierr
296 
297  CALL matzeroentries(matrix, ierr)
298  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
299 
300  n_w = mesh%gauss%n_w
301 
302 
303  DO l = 1, mesh%gauss%l_G
304  DO ni = 1, mesh%gauss%n_w
305  DO nj = 1, mesh%gauss%n_w
306  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
307  END DO
308  END DO
309  END DO
310 
311  DO m = 1, mesh%me
312  jj_loc = mesh%jj(:,m)
313 
314  aij = 0.d0
315  DO l = 1, mesh%gauss%l_G
316  ray = 0.d0
317  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
318  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
319  END DO
320  DO nj = 1, mesh%gauss%n_w; j = jj_loc(nj)
321  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
322  aij(ni,nj) = aij(ni,nj) + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m)
323  ENDDO
324  ENDDO
325  ENDDO
326 
327  mat_loc = 0.d0
328  DO ki= 1, 3
329  DO ni = 1, n_w
330  i = jj_loc(ni)
331  iglob = la%loc_to_glob(ki,i)
332  ix = (ki-1)*n_w+ni
333  idxn(ix) = iglob-1
334  DO kj = 1, 3
335  DO nj = 1, n_w
336  j = jj_loc(nj)
337  jglob = la%loc_to_glob(kj,j)
338  jx = (kj-1)*n_w+nj
339  jdxn(jx) = jglob-1
340  IF (ki==kj) THEN
341  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
342  END IF
343  END DO
344  END DO
345  END DO
346  END DO
347  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
348  END DO
349 
350  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
351  CALL matassemblyend(matrix,mat_final_assembly,ierr)
352 
353  END SUBROUTINE qs_mass_vect_3x3
354 
355 
356  SUBROUTINE qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
357  !========================================================================
358  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
359  ! (V1,V4,V5), (V2,V3,V6)
360  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
361  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
362  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
363  ! 2 pour les composantes 2, 3 et 6
364  !------------------------------------------------------------------------
365  USE my_util
366  USE input_data
367  IMPLICIT NONE
368  INTEGER , INTENT(IN) :: type_op, mode
369  REAL(KIND=8), INTENT(IN) :: visco, mass, stab, stab_bdy_ns
370  TYPE(petsc_csr_la) :: la
371  TYPE(mesh_type), TARGET :: mesh
372  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
373  REAL(KIND=8) :: xij, viscolm
374  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
375  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
376  REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
377  REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
378  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
379  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
380  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
381  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
382  INTEGER :: ix, jx, iglob, jglob
383  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
384  REAL(KIND=8) :: hm, viscomode, hh
385  !=====DEUX INSERTIONS VALEURS A FAIRE....
386 #include "petsc/finclude/petsc.h"
387  mat :: matrix
388  petscerrorcode :: ierr
389 
390  CALL matzeroentries(matrix, ierr)
391  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
392 
393  np = SIZE(mesh%rr,2)
394  n_w = mesh%gauss%n_w
395  n_ws = mesh%gauss%n_ws
396 
397  IF (type_op == 1) THEN
398  eps1 = 1.d0
399  eps2 = 1.d0
400  k_max = 2 !Structure vectorielle
401  ELSEIF (type_op == 2) THEN
402  eps1 = 1.d0
403  eps2 = -1.d0
404  k_max = 2 !Structure vectorielle
405  ELSEIF (type_op == 3) THEN
406  !cas du laplacien scalaire
407  eps1 = 0.d0
408  eps2 = 0.d0
409  k_max = 1 !Structure scalaire
410  ELSE
411  CALL error_petsc('probleme de type d''operateur')
412  ENDIF
413 
414 
415  DO l = 1, mesh%gauss%l_G
416  DO ni = 1, mesh%gauss%n_w
417  DO nj = 1, mesh%gauss%n_w
418  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
419  END DO
420  END DO
421  END DO
422 
423  DO m = 1, mesh%me
424  jj_loc = mesh%jj(:,m)
425 
426  aij = 0.d0
427  bij = 0.d0
428  cij = 0.d0
429  DO l = 1, mesh%gauss%l_G
430 
431  !--------On calcule le rayon du point gauss
432  ray = 0
433  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
434  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
435  END DO
436  hh=mesh%hloc(m)
437 !!$ hm=MIN(0.5d0/inputs%m_max,hh)!WRONG choice
438  hm=0.5d0/inputs%m_max
439  viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m) ! We add the constant stabilization here
440  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
441 
442  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
443 
444  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
445 
446  !grad(u).grad(v) en r et z
447  xij = 0.d0
448  DO k = 1, mesh%gauss%k_d
449  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
450  END DO
451 
452  !blocs diagonaux
453  z = ray * viscolm* xij &
454  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
455  + viscomode*mode**2*wwprod(ni,nj,l)/ray
456  cij(ni,nj) = cij(ni,nj) + z
457  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
458  !blocs couplant
459  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
460  ENDDO
461  ENDDO
462 
463  ENDDO
464 
465  mat_loc = 0.d0
466  DO ki= 1, 3
467  DO ni = 1, n_w
468  i = jj_loc(ni)
469  iglob = la%loc_to_glob(ki,i)
470  ix = (ki-1)*n_w+ni
471  idxn(ix) = iglob-1
472  DO kj = 1, 3
473  DO nj = 1, n_w
474  j = jj_loc(nj)
475  jglob = la%loc_to_glob(kj,j)
476  jx = (kj-1)*n_w+nj
477  jdxn(jx) = jglob-1
478  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
479  IF (ki==kj) THEN
480  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
481  ELSE
482  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
483  END IF
484  ELSE ! ki=3 OR kj=3
485  IF (ki==kj) THEN
486  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
487  END IF
488  END IF
489  END DO
490  END DO
491  END DO
492  END DO
493  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
494 
495  !==Calcul de visco (grad u)T . (grad v)
496  aij = 0.d0
497  bij = 0.d0
498  cij = 0.d0
499  dij = 0.d0
500  eij = 0.d0
501  fij = 0.d0
502  DO l = 1, mesh%gauss%l_G
503 
504  !--------On calcule le rayon du point gauss
505  ray = 0
506  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
507  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
508  END DO
509 
510  viscolm = visco*mesh%gauss%rj(l,m)*ray
511 
512  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
513  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
514  aij(ni,nj) = aij(ni,nj) &
515  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2)
516  bij(ni,nj) = bij(ni,nj) &
517  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2)
518  cij(ni,nj) = cij(ni,nj) &
519  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m))
520  dij(ni,nj) = dij(ni,nj) &
521  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
522  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray)
523  eij(ni,nj) = eij(ni,nj) &
524  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray)
525  fij(ni,nj) = fij(ni,nj) &
526  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
527  END DO
528  END DO
529  END DO
530 
531  mat_loc=0.d0
532  idxn=0
533  jdxn=0
534  DO ni = 1, n_w
535  DO ki = 1, 3
536  i = jj_loc(ni)
537  iglob = la%loc_to_glob(ki,i)
538  ix = (ki-1)*n_w + ni
539  idxn(ix) = iglob-1
540  END DO
541  END DO
542  jdxn=idxn
543 
544  DO ni = 1, n_w
545  DO nj = 1, n_w
546  !=== Line i 1 (Vr)
547  ix = ni
548  !=== Column j 1 (Vr)
549  jx = nj
550  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
551  !=== Column j 2 (Vt)
552  jx = nj+n_w
553  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
554  !=== Column j 3 (Vz)
555  jx = nj+2*n_w
556  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
557 
558  !=== Line i 2 (Vt)
559  ix = ni+n_w
560  !=== Column j 1 (Vr)
561  jx = nj
562  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
563  !=== Column j 2 (Vt)
564  jx = nj+n_w
565  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
566  !=== Column j 3 (Vz)
567  jx = nj+2*n_w
568  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
569 
570  !=== Line i 3 (Vz)
571  ix = ni+2*n_w
572  !=== Column j 1 (Vr)
573  jx = nj
574  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
575  !=== Column j 2 (Vt)
576  jx = nj+n_w
577  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
578  !=== Column j 3 (Vz)
579  jx = nj+2*n_w
580  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
581  END DO
582  END DO
583  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
584  ENDDO
585  !== Fin du Calcul de visco (grad u)T . (grad v)
586 
587  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
588  !===Surface terms
589  stab_normal = coeff_stab_normal*(1.d0+visco)
590  !===Normalization is not done properly: (1.d0+visco)
591  !===To be revisited some day.
592  DO ms = 1, mesh%mes
593  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
594  aij = 0.d0
595  bij = 0.d0
596  cij = 0.d0
597  dij = 0.d0
598  DO ls = 1, mesh%gauss%l_Gs
599  !--------On calcule le rayon du point gauss
600  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
601  IF (ray.LT.1.d-10) cycle
602  hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
603  x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
604  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
605 
606  DO ni = 1, mesh%gauss%n_ws
607  DO nj = 1, mesh%gauss%n_ws
608  y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
609  aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
610  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
611  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
612  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
613  bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
614  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
615  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
616  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
617  cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
618  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
619  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
620  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
621  dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
622  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
623  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
624  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
625  END DO
626  END DO
627  END DO
628 
629  !++++++++++++++++
630  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
631  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
632  mat_loc_s = 0.d0
633  idxn_s = 0
634  jdxn_s = 0
635  DO ki = 1, 2
636  DO ni = 1, n_ws
637  i = mesh%jjs(ni,ms)
638  iglob = la%loc_to_glob(2*ki-1,i)
639  ix = (ki-1)*n_ws+ni
640  idxn_s(ix) = iglob-1
641  DO kj = 1, 2
642  DO nj = 1, n_ws
643  j = mesh%jjs(nj,ms)
644  jglob = la%loc_to_glob(2*kj-1,j)
645  jx = (kj-1)*n_ws+nj
646  jdxn_s(jx) = jglob-1
647  IF ((ki == 1) .AND. (kj == 1)) THEN
648  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
649  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
650  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
651  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
652  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
653  ELSE
654  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
655  END IF
656  END DO
657  END DO
658  END DO
659  END DO
660  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
661  END DO
662  END IF
663 
664  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
665  CALL matassemblyend(matrix,mat_final_assembly,ierr)
666 
667  END SUBROUTINE qs_diff_mass_vect_3x3
668 
669  SUBROUTINE qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
670  !========================================================================
671  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
672  ! (V1,V4,V5), (V2,V3,V6)
673  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
674  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
675  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
676  ! 2 pour les composantes 2, 3 et 6
677  !------------------------------------------------------------------------
678  USE my_util
679  USE input_data
680  IMPLICIT NONE
681 
682 
683  INTEGER , INTENT(IN) :: type_op, mode
684  REAL(KIND=8), INTENT(IN) :: visco, mass, stab, stab_bdy_ns
685  TYPE(petsc_csr_la) :: la
686  TYPE(mesh_type), TARGET :: mesh
687 
688 !!$ INTEGER, DIMENSION(:,:), POINTER :: jj
689 !!$ INTEGER, POINTER :: me
690 !!$ REAL(KIND=8), DIMENSION(:,:) , POINTER :: rr
691 
692  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
693  REAL(KIND=8) :: xij, viscolm, div_penal
694  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
695  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
696  REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
697  REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
698  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
699  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
700  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
701  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
702  INTEGER :: ix, jx, iglob, jglob
703  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
704  REAL(KIND=8) :: viscomode, hm
705  !=====DEUX INSERTIONS VALEURS A FAIRE....
706 #include "petsc/finclude/petsc.h"
707  mat :: matrix
708  petscerrorcode :: ierr
709 
710  CALL matzeroentries(matrix, ierr)
711  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
712 
713  np = SIZE(mesh%rr,2)
714  n_w = mesh%gauss%n_w
715  n_ws = mesh%gauss%n_ws
716 
717  IF (type_op == 1) THEN
718  eps1 = 1.d0
719  eps2 = 1.d0
720  k_max = 2 !Structure vectorielle
721  ELSEIF (type_op == 2) THEN
722  eps1 = 1.d0
723  eps2 = -1.d0
724  k_max = 2 !Structure vectorielle
725  ELSEIF (type_op == 3) THEN
726  !cas du laplacien scalaire
727  eps1 = 0.d0
728  eps2 = 0.d0
729  k_max = 1 !Structure scalaire
730  ELSE
731  CALL error_petsc('probleme de type d''operateur')
732  ENDIF
733 
734 
735  DO l = 1, mesh%gauss%l_G
736  DO ni = 1, mesh%gauss%n_w
737  DO nj = 1, mesh%gauss%n_w
738  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
739  END DO
740  END DO
741  END DO
742 
743  DO m = 1, mesh%me
744  jj_loc = mesh%jj(:,m)
745 
746  aij = 0.d0
747  bij = 0.d0
748  cij = 0.d0
749  DO l = 1, mesh%gauss%l_G
750 
751  !--------On calcule le rayon du point gauss
752  ray = 0
753  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
754  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
755  END DO
756 
757  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m) ! We add the constant stabilization here
758  hm=0.5d0/inputs%m_max
759  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
760 
761  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
762 
763  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
764 
765  !grad(u).grad(v) en r et z
766  xij = 0.d0
767  DO k = 1, mesh%gauss%k_d
768  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
769  END DO
770 
771  !blocs diagonaux
772  z = ray * viscolm* xij &
773  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
774  + viscomode*mode**2*wwprod(ni,nj,l)/ray
775  cij(ni,nj) = cij(ni,nj) + z
776  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
777  !blocs couplant
778  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
779  ENDDO
780  ENDDO
781 
782  ENDDO
783 
784 
785 
786  !++++++++++++++++
787 !!$ DO ki= 1, k_max
788 !!$ DO ni = 1, mesh%gauss%n_w
789 !!$ i = mesh%jj(ni, m)
790 !!$ ib = i + (ki-1)*np
791 !!$ DO kj = 1, k_max
792 !!$ DO nj = 1, mesh%gauss%n_w
793 !!$ j = mesh%jj(nj, m)
794 !!$ jb = j + (kj-1)*np
795 !!$ DO p = ia(ib), ia(ib+1) - 1
796 !!$ IF (ja(p) == jb) THEN
797 !!$ IF (ki==kj) THEN
798 !!$ a0(p) = a0(p) + aij(ni,nj)
799 !!$ ELSE
800 !!$ a0(p) = a0(p) + bij(ni,nj)
801 !!$ END IF
802 !!$ EXIT
803 !!$ ENDIF
804 !!$ END DO
805 !!$ END DO
806 !!$ END DO
807 !!$ END DO
808 !!$ END DO
809 !!$
810 !!$ IF (type_op /= 3) THEN
811 !!$ ki= 3
812 !!$ DO ni = 1, mesh%gauss%n_w
813 !!$ i = mesh%jj(ni, m)
814 !!$ ib = i + (ki-1)*np
815 !!$ kj = 3
816 !!$ DO nj = 1, mesh%gauss%n_w
817 !!$ j = mesh%jj(nj, m)
818 !!$ jb = j + (kj-1)*np
819 !!$ DO p = ia(ib), ia(ib+1) - 1
820 !!$ IF (ja(p) == jb) THEN
821 !!$ a0(p) = a0(p) + cij(ni,nj)
822 !!$ EXIT
823 !!$ ENDIF
824 !!$ END DO
825 !!$ END DO
826 !!$ END DO
827 !!$ END IF
828 
829  mat_loc = 0.d0
830  DO ki= 1, 3
831  DO ni = 1, n_w
832  i = jj_loc(ni)
833  iglob = la%loc_to_glob(ki,i)
834  ix = (ki-1)*n_w+ni
835  idxn(ix) = iglob-1
836  DO kj = 1, 3
837  DO nj = 1, n_w
838  j = jj_loc(nj)
839  jglob = la%loc_to_glob(kj,j)
840  jx = (kj-1)*n_w+nj
841  jdxn(jx) = jglob-1
842  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
843  IF (ki==kj) THEN
844  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
845  ELSE
846  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
847  END IF
848  ELSE ! ki=3 OR kj=3
849  IF (ki==kj) THEN
850  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
851  END IF
852  END IF
853  END DO
854  END DO
855  END DO
856  END DO
857  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
858 
859  !+++---------------------
860 !!$ IF (type_op==3) CYCLE
861  !==Calcul de visco (grad u)T . (grad v)
862  aij = 0.d0
863  bij = 0.d0
864  cij = 0.d0
865  dij = 0.d0
866  eij = 0.d0
867  fij = 0.d0
868  DO l = 1, mesh%gauss%l_G
869 
870  !--------On calcule le rayon du point gauss
871  ray = 0
872  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
873  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
874  END DO
875 
876  viscolm = visco*mesh%gauss%rj(l,m)*ray
877  div_penal = inputs%div_stab_in_ns*mesh%gauss%rj(l,m)*ray
878 
879  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
880  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
881  aij(ni,nj) = aij(ni,nj) &
882  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
883  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
884  + mesh%gauss%ww(nj,l)/ray)
885  bij(ni,nj) = bij(ni,nj) &
886  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
887  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
888  cij(ni,nj) = cij(ni,nj) &
889  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
890  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
891  dij(ni,nj) = dij(ni,nj) &
892  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
893  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
894  + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
895  eij(ni,nj) = eij(ni,nj) &
896  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
897  + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
898  fij(ni,nj) = fij(ni,nj) &
899  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
900  + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
901  END DO
902  END DO
903  END DO
904  !aij = 0.d0
905  !bij = 0.d0
906  !cij = 0.d0
907  !dij = 0.d0
908  !eij = 0.d0
909  !fij = 0.d0
910  !++++++++++++++++
911  mat_loc=0.d0
912  idxn=0
913  jdxn=0
914  DO ni = 1, n_w
915  DO ki = 1, 3
916  i = jj_loc(ni)
917  iglob = la%loc_to_glob(ki,i)
918  ix = (ki-1)*n_w + ni
919  idxn(ix) = iglob-1
920  END DO
921  END DO
922  jdxn=idxn
923 
924  DO ni = 1, n_w
925  DO nj = 1, n_w
926  !=== Line i 1 (Vr)
927  ix = ni
928  !=== Column j 1 (Vr)
929  jx = nj
930  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
931  !=== Column j 2 (Vt)
932  jx = nj+n_w
933  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
934  !=== Column j 3 (Vz)
935  jx = nj+2*n_w
936  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
937 
938  !=== Line i 2 (Vt)
939  ix = ni+n_w
940  !=== Column j 1 (Vr)
941  jx = nj
942  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
943  !=== Column j 2 (Vt)
944  jx = nj+n_w
945  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
946  !=== Column j 3 (Vz)
947  jx = nj+2*n_w
948  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
949 
950  !=== Line i 3 (Vz)
951  ix = ni+2*n_w
952  !=== Column j 1 (Vr)
953  jx = nj
954  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
955  !=== Column j 2 (Vt)
956  jx = nj+n_w
957  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
958  !=== Column j 3 (Vz)
959  jx = nj+2*n_w
960  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
961  END DO
962  END DO
963  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
964 
965 !!$ DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
966 !!$ DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
967 !!$ ib = i
968 !!$ jb = j
969 !!$ DO p = ia(ib), ia(ib+1) - 1
970 !!$ IF (ja(p) == jb) THEN
971 !!$ a0(p) = a0(p) + aij(ni,nj)
972 !!$ EXIT
973 !!$ ENDIF
974 !!$ END DO
975 !!$ jb= j + np
976 !!$ DO p = ia(ib), ia(ib+1) - 1
977 !!$ IF (ja(p) == jb) THEN
978 !!$ a0(p) = a0(p) + bij(ni,nj)
979 !!$ EXIT
980 !!$ ENDIF
981 !!$ END DO
982 !!$ jb = j + 2*np
983 !!$ DO p = ia(ib), ia(ib+1) - 1
984 !!$ IF (ja(p) == jb) THEN
985 !!$ a0(p) = a0(p) + cij(ni,nj)
986 !!$ EXIT
987 !!$ ENDIF
988 !!$ END DO
989 !!$ ib = i + np
990 !!$ jb = j
991 !!$ DO p = ia(ib), ia(ib+1) - 1
992 !!$ IF (ja(p) == jb) THEN
993 !!$ a0(p) = a0(p) + bij(nj,ni)
994 !!$ EXIT
995 !!$ ENDIF
996 !!$ END DO
997 !!$ jb = j + np
998 !!$ DO p = ia(ib), ia(ib+1) - 1
999 !!$ IF (ja(p) == jb) THEN
1000 !!$ a0(p) = a0(p) + dij(ni,nj)
1001 !!$ EXIT
1002 !!$ ENDIF
1003 !!$ END DO
1004 !!$ jb = j + 2*np
1005 !!$ DO p = ia(ib), ia(ib+1) - 1
1006 !!$ IF (ja(p) == jb) THEN
1007 !!$ a0(p) = a0(p) + eij(ni,nj)
1008 !!$ EXIT
1009 !!$ ENDIF
1010 !!$ END DO
1011 !!$ ib = i + 2*np
1012 !!$ jb = j
1013 !!$ DO p = ia(ib), ia(ib+1) - 1
1014 !!$ IF (ja(p) == jb) THEN
1015 !!$ a0(p) = a0(p) + cij(nj,ni)
1016 !!$ EXIT
1017 !!$ ENDIF
1018 !!$ END DO
1019 !!$ jb = j + np
1020 !!$ DO p = ia(ib), ia(ib+1) - 1
1021 !!$ IF (ja(p) == jb) THEN
1022 !!$ a0(p) = a0(p) + eij(nj,ni)
1023 !!$ EXIT
1024 !!$ ENDIF
1025 !!$ END DO
1026 !!$ jb = j + 2*np
1027 !!$ DO p = ia(ib), ia(ib+1) - 1
1028 !!$ IF (ja(p) == jb) THEN
1029 !!$ a0(p) = a0(p) + fij(ni,nj)
1030 !!$ EXIT
1031 !!$ ENDIF
1032 !!$ END DO
1033 !!$ END DO
1034 !!$ END DO
1035  ENDDO
1036  !== Fin du Calcul de visco (grad u)T . (grad v)
1037  !++++++++++++++++------
1038 
1039  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
1040  !===Surface terms
1041  stab_normal = coeff_stab_normal*(1.d0+visco)
1042  DO ms = 1, mesh%mes
1043  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1044  aij = 0.d0
1045  bij = 0.d0
1046  cij = 0.d0
1047  dij = 0.d0
1048  DO ls = 1, mesh%gauss%l_Gs
1049  !--------On calcule le rayon du point gauss
1050  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1051  IF (ray.LT.1.d-10) cycle
1052  hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1053  x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1054  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1055 
1056  DO ni = 1, mesh%gauss%n_ws
1057  DO nj = 1, mesh%gauss%n_ws
1058  y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1059  aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1060  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1061  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1062  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1063  bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1064  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1065  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1066  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1067  cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1068  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1069  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1070  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1071  dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1072  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1073  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1074  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1075  END DO
1076  END DO
1077  END DO
1078 
1079  !++++++++++++++++
1080  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
1081  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
1082  mat_loc_s = 0.d0
1083  idxn_s = 0
1084  jdxn_s = 0
1085  DO ki = 1, 2
1086  DO ni = 1, n_ws
1087  i = mesh%jjs(ni,ms)
1088  iglob = la%loc_to_glob(2*ki-1,i)
1089  ix = (ki-1)*n_ws+ni
1090  idxn_s(ix) = iglob-1
1091  DO kj = 1, 2
1092  DO nj = 1, n_ws
1093  j = mesh%jjs(nj,ms)
1094  jglob = la%loc_to_glob(2*kj-1,j)
1095  jx = (kj-1)*n_ws+nj
1096  jdxn_s(jx) = jglob-1
1097  IF ((ki == 1) .AND. (kj == 1)) THEN
1098  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1099  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
1100  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1101  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
1102  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1103  ELSE
1104  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1105  END IF
1106  END DO
1107  END DO
1108  END DO
1109  END DO
1110  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1111 
1112 !!$ DO ki= 1, 3, 2
1113 !!$ DO ni = 1, mesh%gauss%n_ws
1114 !!$ i = mesh%jjs(ni, ms)
1115 !!$ ib = i + (ki-1)*np
1116 !!$ DO kj = 1, 3, 2
1117 !!$ DO nj = 1, mesh%gauss%n_ws
1118 !!$ j = mesh%jjs(nj, ms)
1119 !!$ jb = j + (kj-1)*np
1120 !!$ DO p = ia(ib), ia(ib+1) - 1
1121 !!$ IF (ja(p) == jb) THEN
1122 !!$ IF (ki == 1 .AND. kj == 1) THEN
1123 !!$ a0(p) = a0(p) + aij(ni,nj) + aij(nj,ni)
1124 !!$ ELSE IF (ki == 1 .AND. kj == 3) THEN
1125 !!$ a0(p) = a0(p) + bij(ni,nj) + cij(nj,ni)
1126 !!$ ELSE IF (ki == 3 .AND. kj == 1) THEN
1127 !!$ a0(p) = a0(p) + cij(ni,nj) + bij(nj,ni)
1128 !!$ ELSE
1129 !!$ a0(p) = a0(p) + dij(ni,nj) + dij(nj,ni)
1130 !!$ END IF
1131 !!$ EXIT
1132 !!$ ENDIF
1133 !!$ END DO
1134 !!$ END DO
1135 !!$ END DO
1136 !!$ END DO
1137 !!$ END DO
1138  !++++++++++++++++----
1139 
1140  END DO
1141  END IF
1142 
1143  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1144  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1145 
1146 
1147  END SUBROUTINE qs_diff_mass_vect_3x3_divpenal
1148 
1149 
1150 
1151 
1152  SUBROUTINE qs_00_m (mesh, alpha, LA, matrix)
1153  !=================================================
1154  IMPLICIT NONE
1155  TYPE(mesh_type), INTENT(IN) :: mesh
1156  REAL(KIND=8), INTENT(IN) :: alpha
1157  TYPE(petsc_csr_la) :: la
1158  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1159  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1160  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1161  INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1162  REAL(KIND=8) :: al, x, ray
1163 #include "petsc/finclude/petsc.h"
1164  mat :: matrix
1165  petscerrorcode :: ierr
1166  CALL matzeroentries(matrix, ierr)
1167 
1168  DO m = 1, mesh%dom_me
1169  jj_loc = mesh%jj(:,m)
1170  mat_loc = 0.d0
1171  DO l = 1, mesh%gauss%l_G
1172  !Compute radius of Gauss point
1173  ray = 0
1174  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1175  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1176  END DO
1177 
1178  al = alpha *mesh%gauss%rj(l,m)*ray
1179  DO nj = 1, mesh%gauss%n_w;
1180  j = jj_loc(nj)
1181  jglob = la%loc_to_glob(1,j)
1182  idxn(nj) = jglob-1
1183  DO ni = 1, mesh%gauss%n_w;
1184  i = jj_loc(ni)
1185  iglob = la%loc_to_glob(1,i)
1186  idxm(ni) = iglob-1
1187  x = al*mesh%gauss%ww(nj,l)*mesh%gauss%ww(ni,l)
1188  mat_loc(ni,nj) = mat_loc(ni,nj) + x
1189  ENDDO
1190  ENDDO
1191  ENDDO
1192  CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1193  ENDDO
1194  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1195  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1196  END SUBROUTINE qs_00_m
1197 
1198 SUBROUTINE qs_11_m (mesh, visco, LA, matrix)
1199  !=================================================
1200  USE def_type_mesh
1201  IMPLICIT NONE
1202  TYPE(mesh_type), TARGET :: mesh
1203  REAL(KIND=8), INTENT(IN) :: visco
1204  TYPE(petsc_csr_la) :: la
1205  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1206  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1207  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1208  INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1209  REAL(KIND=8) :: al, x, ray
1210 #include "petsc/finclude/petsc.h"
1211  mat :: matrix
1212  petscerrorcode :: ierr
1213  CALL matzeroentries(matrix, ierr)
1214 
1215  DO m = 1, mesh%dom_me
1216  jj_loc = mesh%jj(:,m)
1217  mat_loc = 0.d0
1218 
1219  DO nj = 1, mesh%gauss%n_w;
1220  j = jj_loc(nj)
1221  jglob = la%loc_to_glob(1,j)
1222  idxn(nj) = jglob-1
1223  DO ni = 1, mesh%gauss%n_w;
1224  i = jj_loc(ni)
1225  iglob = la%loc_to_glob(1,i)
1226  idxm(ni) = iglob-1
1227  END DO
1228  END DO
1229 
1230  DO l = 1, mesh%gauss%l_G
1231  !Compute radius of Gauss point
1232  ray = 0
1233  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1234  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1235  END DO
1236  al = visco * mesh%gauss%rj(l,m)
1237  DO nj = 1, mesh%gauss%n_w;
1238  DO ni = 1, mesh%gauss%n_w;
1239  x =al*sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1240  mat_loc(ni,nj) = mat_loc(ni,nj) + x
1241  ENDDO
1242  ENDDO
1243  ENDDO
1244  CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1245  ENDDO
1246  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1247  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1248  END SUBROUTINE qs_11_m
1249 
1250  SUBROUTINE qs_diff_mass_scal_m_level (mesh, LA, visco, mass, stab, mode, matrix)
1251  !=================================================
1252  USE input_data
1253  IMPLICIT NONE
1254  TYPE(mesh_type), INTENT(IN) :: mesh
1255  TYPE(petsc_csr_la) :: la
1256  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
1257  INTEGER, INTENT(IN) :: mode
1258  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1259  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1260  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1261  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
1262  REAL(KIND=8) :: ray
1263  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
1264  REAL(KIND=8) :: viscolm, xij, hm, viscomode, hh
1265 #include "petsc/finclude/petsc.h"
1266  mat :: matrix
1267  petscerrorcode :: ierr
1268  CALL matzeroentries(matrix, ierr)
1269  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
1270 
1271  DO l = 1, mesh%gauss%l_G
1272  DO ni = 1, mesh%gauss%n_w
1273  DO nj = 1, mesh%gauss%n_w
1274  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1275  END DO
1276  END DO
1277  END DO
1278 
1279  n_w = mesh%gauss%n_w
1280  DO m = 1, mesh%dom_me
1281  jj_loc = mesh%jj(:,m)
1282  a_loc = 0.d0
1283 
1284  DO nj = 1, mesh%gauss%n_w;
1285  j = jj_loc(nj)
1286  jglob = la%loc_to_glob(1,j)
1287  idxn(nj) = jglob-1
1288  DO ni = 1, mesh%gauss%n_w;
1289  i = jj_loc(ni)
1290  iglob = la%loc_to_glob(1,i)
1291  idxm(ni) = iglob-1
1292  END DO
1293  END DO
1294 
1295  DO l = 1, mesh%gauss%l_G
1296  !Compute radius of Gauss point
1297  ray = 0
1298  DO ni = 1, n_w; i = jj_loc(ni)
1299  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1300  END DO
1301  hh=mesh%hloc(m)
1302 !!$ hm=MIN(0.5d0/inputs%m_max,hh)!WRONG choice
1303  hm=0.5d0/inputs%m_max
1304  viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
1305  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1306  DO nj = 1, n_w
1307  DO ni = 1, n_w
1308  !grad(u).grad(v) w.r.t. r and z
1309  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1310  !start diagonal block
1311  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
1312  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1313  + viscomode*mode**2*wwprod(ni,nj,l)/ray
1314  !end diagonal block
1315  END DO
1316  END DO
1317  END DO
1318  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
1319  ENDDO
1320  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1321  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1322 
1323  END SUBROUTINE qs_diff_mass_scal_m_level
1324 
1325 END MODULE fem_m_axi
1326 
1327 
1328 
1329 
1330 
1331 
1332 
1333 
1334 
1335 
1336 
1337 
1338 
1339 
fem_m_axi
Definition: fem_M_axi.f90:6
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
fem_m_axi::qs_diff_mass_vect_m
subroutine qs_diff_mass_vect_m(type_op, LA, mesh, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:11
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
fem_m_axi::qs_diff_mass_scal_m
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:127
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
fem_m_axi::qs_diff_mass_scal_m_variant
subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:198
input_data
Definition: read_my_data.f90:200
fem_s_m::qs_11_m
subroutine qs_11_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:171
fem_m_axi::qs_mass_vect_3x3
subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix)
Definition: fem_M_axi.f90:269
fem_m_axi::qs_diff_mass_vect_3x3
subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
Definition: fem_M_axi.f90:356
fem_m_axi::qs_diff_mass_scal_m_level
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:1250
fem_s_m::qs_00_m
subroutine qs_00_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:6
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
z
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic z
Definition: doc_intro.h:193
fem_m_axi::qs_diff_mass_vect_3x3_divpenal
subroutine qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
Definition: fem_M_axi.f90:669