SFEMaNS  version 5.3
Reference documentation for SFEMaNS
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 #include "petsc/finclude/petsc.h"
9  USE petsc
10 CONTAINS
11 
12 
13  SUBROUTINE qs_diff_mass_vect_m (type_op, LA, mesh, visco, mass, stab, mode, matrix)
14  !=================================================
15  USE my_util
16  IMPLICIT NONE
17  TYPE(mesh_type), INTENT(IN) :: mesh
18  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
19  INTEGER, INTENT(IN) :: type_op, mode
20  TYPE(petsc_csr_la) :: LA
21  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
22  INTEGER, DIMENSION(:), ALLOCATABLE :: idxm, idxn
23  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
24  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc, b_loc
25  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat_loc
26  REAL(KIND=8) :: ray, eps1, eps2
27  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, k_max, n_w, ix, jx, ki, kj
28  REAL(KIND=8) :: viscolm, xij
29 !#include "petsc/finclude/petsc.h"
30  mat :: matrix
31  petscerrorcode :: ierr
32  CALL matzeroentries (matrix, ierr)
33  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
34 
35 
36  IF (type_op == 1) THEN
37  eps1 = 1.d0
38  eps2 = 1.d0
39  k_max = 2 ! 2x2 Structure
40  ELSEIF (type_op == 2) THEN
41  eps1 = 1.d0
42  eps2 = -1.d0
43  k_max = 2 !2x2 Structure
44  ELSEIF (type_op == 3) THEN
45  eps1 = 0.d0
46  eps2 = 0.d0
47  k_max = 1 ! Scalar Structure
48  ELSE
49  eps1 = 0.d0
50  eps2 = 0.d0
51  k_max = 1
52  CALL error_petsc('BUG in qs_diff_mass_vect_M, type_op not correct')
53  ENDIF
54 
55  IF (k_max/=SIZE(la%loc_to_glob,1)) THEN
56  CALL error_petsc('BUG in qs_diff_mass_vect_petsc_M, k_max/=SIZE(LA%loc_to_glob,1)')
57  END IF
58 
59  n_w = mesh%gauss%n_w
60  DO l = 1, mesh%gauss%l_G
61  DO ni = 1, n_w
62  DO nj = 1, n_w
63  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
64  END DO
65  END DO
66  END DO
67 
68  ALLOCATE(mat_loc(k_max*n_w,k_max*n_w),idxm(k_max*n_w),idxn(k_max*n_w))
69  DO m = 1, mesh%dom_me
70  jj_loc = mesh%jj(:,m)
71 
72  a_loc = 0.d0
73  b_loc = 0.d0
74  DO l = 1, mesh%gauss%l_G
75  !Compute radius of Gauss point
76  ray = 0
77  DO ni = 1, n_w; i = jj_loc(ni)
78  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
79  END DO
80  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
81  DO nj = 1, n_w
82  DO ni = 1, n_w
83 
84  !grad(u).grad(v) w.r.t. r and z
85  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
86 
87  !start diagonal block
88  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
89  + viscolm*eps1*wwprod(ni,nj,l)/ray &
90  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
91  + viscolm*mode**2*wwprod(ni,nj,l)/ray
92  !end diagonal block
93  b_loc(ni,nj) = b_loc(ni,nj) + eps2*viscolm*2*mode*wwprod(ni,nj,l)/ray
94  END DO
95  END DO
96  END DO
97 
98  mat_loc = 0.d0
99  DO ki= 1, k_max
100  DO ni = 1, n_w
101  i = jj_loc(ni)
102  iglob = la%loc_to_glob(ki,i)
103  ix = (ki-1)*n_w+ni
104  idxm(ix) = iglob-1
105  DO kj = 1, k_max
106  DO nj = 1, n_w
107  j = jj_loc(nj)
108  jglob = la%loc_to_glob(kj,j)
109  jx = (kj-1)*n_w+nj
110  idxn(jx) = jglob-1
111  IF (ki==kj) THEN
112  mat_loc(ix,jx) = mat_loc(ix,jx) + a_loc(ni,nj)
113  ELSE
114  mat_loc(ix,jx) = mat_loc(ix,jx) + b_loc(ni,nj)
115  END IF
116  END DO
117  END DO
118  END DO
119  END DO
120  CALL matsetvalues(matrix, k_max*n_w, idxm, k_max*n_w, idxn, mat_loc, add_values, ierr)
121  ENDDO
122  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
123  CALL matassemblyend(matrix,mat_final_assembly,ierr)
124 
125  DEALLOCATE(mat_loc,idxm,idxn)
126 
127  END SUBROUTINE qs_diff_mass_vect_m
128 
129  SUBROUTINE qs_diff_mass_scal_m (mesh, LA, visco, mass, stab, mode, matrix)
130  !=================================================
131  IMPLICIT NONE
132  TYPE(mesh_type), INTENT(IN) :: mesh
133  TYPE(petsc_csr_la) :: LA
134  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
135  INTEGER, INTENT(IN) :: mode
136  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
137  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
138  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
139  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
140  REAL(KIND=8) :: ray
141  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
142  REAL(KIND=8) :: viscolm, xij
143 !#include "petsc/finclude/petsc.h"
144  mat :: matrix
145  petscerrorcode :: ierr
146  CALL matzeroentries (matrix, ierr)
147  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
148 
149  DO l = 1, mesh%gauss%l_G
150  DO ni = 1, mesh%gauss%n_w
151  DO nj = 1, mesh%gauss%n_w
152  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
153  END DO
154  END DO
155  END DO
156 
157  n_w = mesh%gauss%n_w
158  DO m = 1, mesh%dom_me
159  jj_loc = mesh%jj(:,m)
160  a_loc = 0.d0
161 
162  DO nj = 1, mesh%gauss%n_w;
163  j = jj_loc(nj)
164  jglob = la%loc_to_glob(1,j)
165  idxn(nj) = jglob-1
166  DO ni = 1, mesh%gauss%n_w;
167  i = jj_loc(ni)
168  iglob = la%loc_to_glob(1,i)
169  idxm(ni) = iglob-1
170  END DO
171  END DO
172 
173  DO l = 1, mesh%gauss%l_G
174  !Compute radius of Gauss point
175  ray = 0
176  DO ni = 1, n_w; i = jj_loc(ni)
177  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
178  END DO
179  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
180 
181  DO nj = 1, n_w
182  DO ni = 1, n_w
183  !grad(u).grad(v) w.r.t. r and z
184  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
185  !start diagonal block
186  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
187  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
188  + viscolm*mode**2*wwprod(ni,nj,l)/ray
189  !end diagonal block
190  END DO
191  END DO
192  END DO
193  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
194  ENDDO
195  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
196  CALL matassemblyend(matrix,mat_final_assembly,ierr)
197 
198  END SUBROUTINE qs_diff_mass_scal_m
199 
200  SUBROUTINE qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, temp_list_robin_sides, &
201  convection_coeff, stab, mode, matrix)
203  !=================================================
204  IMPLICIT NONE
205  TYPE(mesh_type), INTENT(IN) :: mesh
206  TYPE(petsc_csr_la) :: LA
207  REAL(KIND=8), INTENT(IN) :: mass, stab
208  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: heat_capa, visco, convection_coeff ! MODIFICATION: Robin coeff = convection coeff h for temperature case
209  INTEGER, DIMENSION(:), INTENT(IN) :: temp_list_robin_sides ! MODIFICATION: robin sides added to build the matrix
210  INTEGER, INTENT(IN) :: mode
211  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
212  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
213  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
214  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
215  REAL(KIND=8) :: ray
216  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w, n_ws, ms, ls, ib ! MODIFICATION: for Robin
217  REAL(KIND=8) :: viscolm, xij, x ! MODIFICATION: for Robin
218  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,mesh%gauss%n_ws) :: h_phii_phij ! MODIFICATION: terms for Robin
219  INTEGER, DIMENSION(1:1) :: coeff_index ! MODIFICATION: index of robin coefficient in the list
220 
221 !#include "petsc/finclude/petsc.h"
222  mat :: matrix
223  petscerrorcode :: ierr
224  CALL matzeroentries (matrix, ierr)
225  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
226 
227  DO l = 1, mesh%gauss%l_G
228  DO ni = 1, mesh%gauss%n_w
229  DO nj = 1, mesh%gauss%n_w
230  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
231  END DO
232  END DO
233  END DO
234 
235  n_w = mesh%gauss%n_w
236  DO m = 1, mesh%dom_me
237  jj_loc = mesh%jj(:,m)
238  a_loc = 0.d0
239 
240  DO nj = 1, mesh%gauss%n_w;
241  j = jj_loc(nj)
242  jglob = la%loc_to_glob(1,j)
243  idxn(nj) = jglob-1
244  DO ni = 1, mesh%gauss%n_w;
245  i = jj_loc(ni)
246  iglob = la%loc_to_glob(1,i)
247  idxm(ni) = iglob-1
248  END DO
249  END DO
250 
251  DO l = 1, mesh%gauss%l_G
252  !Compute radius of Gauss point
253  ray = 0
254  DO ni = 1, n_w; i = jj_loc(ni)
255  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
256  END DO
257  viscolm = (visco(m) + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
258  DO nj = 1, n_w
259  DO ni = 1, n_w
260  !grad(u).grad(v) w.r.t. r and z
261  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
262  !start diagonal block
263  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
264  + heat_capa(m) * mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
265  + viscolm*mode**2*wwprod(ni,nj,l)/ray
266  !end diagonal block
267  END DO
268  END DO
269  END DO
270  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
271  ENDDO
272 
273  !===Robin conditions ! MODIFICATION: Robin: addition of the term int_(partial Omega) h*u*v, with h the convection coefficient
274 
275  IF (SIZE(temp_list_robin_sides) > 0) THEN
276 
277  n_ws = mesh%gauss%n_ws
278  DO ms = 1, mesh%mes
279  IF (minval(abs(temp_list_robin_sides - mesh%sides(ms))) > 0) cycle
280  h_phii_phij = 0d0
281  coeff_index = minloc(abs(temp_list_robin_sides - mesh%sides(ms)))
282  DO ls = 1, mesh%gauss%l_Gs
283  !===Compute radius of Gauss point
284  ray = sum(mesh%rr(1,mesh%jjs(:,ms))* mesh%gauss%wws(:,ls))
285  x = convection_coeff(coeff_index(1)) * ray * mesh%gauss%rjs(ls,ms)
286  DO ni=1, n_ws
287  DO nj=1, n_ws
288  h_phii_phij(ni,nj) = h_phii_phij(ni,nj) + &
289  x * mesh%gauss%wws(ni,ls) * mesh%gauss%wws(nj,ls)
290  ENDDO
291  ENDDO
292  ENDDO
293  DO ni = 1, n_ws
294  i = mesh%jjs(ni,ms)
295  ib = la%loc_to_glob(1,i)
296  idxn(ni) = ib - 1
297  END DO
298  CALL matsetvalues(matrix, n_ws, idxn(1:n_ws), n_ws, idxn(1:n_ws), &
299  h_phii_phij, add_values, ierr)
300  END DO
301 
302  END IF
303 
304  !===End Robin conditions
305 
306  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
307  CALL matassemblyend(matrix,mat_final_assembly,ierr)
308 
309  END SUBROUTINE qs_diff_mass_scal_m_variant
310 
311  SUBROUTINE qs_mass_vect_3x3(LA, mesh, mass, matrix)
312  !========================================================================
313  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
314  ! (V1,V4,V5), (V2,V3,V6)
315  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
316  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
317  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
318  ! 2 pour les composantes 2, 3 et 6
319  !------------------------------------------------------------------------
320  USE my_util
321  USE input_data
322  IMPLICIT NONE
323 
324  REAL(KIND=8), INTENT(IN) :: mass
325  TYPE(petsc_csr_la) :: LA
326  TYPE(mesh_type) :: mesh
327  INTEGER ::l, m, ni, nj, i, j, ki, kj, n_w
328  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
329  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
330  REAL(KIND=8) :: ray
331  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
332  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
333  INTEGER :: ix, jx, iglob, jglob
334  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
335 !#include "petsc/finclude/petsc.h"
336  mat :: matrix
337  petscerrorcode :: ierr
338 
339  CALL matzeroentries (matrix, ierr)
340  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
341 
342  n_w = mesh%gauss%n_w
343 
344 
345  DO l = 1, mesh%gauss%l_G
346  DO ni = 1, mesh%gauss%n_w
347  DO nj = 1, mesh%gauss%n_w
348  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
349  END DO
350  END DO
351  END DO
352 
353  DO m = 1, mesh%me
354  jj_loc = mesh%jj(:,m)
355 
356  aij = 0.d0
357  DO l = 1, mesh%gauss%l_G
358  ray = 0.d0
359  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
360  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
361  END DO
362  DO nj = 1, mesh%gauss%n_w; j = jj_loc(nj)
363  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
364  aij(ni,nj) = aij(ni,nj) + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m)
365  ENDDO
366  ENDDO
367  ENDDO
368 
369  mat_loc = 0.d0
370  DO ki= 1, 3
371  DO ni = 1, n_w
372  i = jj_loc(ni)
373  iglob = la%loc_to_glob(ki,i)
374  ix = (ki-1)*n_w+ni
375  idxn(ix) = iglob-1
376  DO kj = 1, 3
377  DO nj = 1, n_w
378  j = jj_loc(nj)
379  jglob = la%loc_to_glob(kj,j)
380  jx = (kj-1)*n_w+nj
381  jdxn(jx) = jglob-1
382  IF (ki==kj) THEN
383  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
384  END IF
385  END DO
386  END DO
387  END DO
388  END DO
389  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
390  END DO
391 
392  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
393  CALL matassemblyend(matrix,mat_final_assembly,ierr)
394 
395  END SUBROUTINE qs_mass_vect_3x3
396 
397 
398  SUBROUTINE qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
399  !========================================================================
400  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
401  ! (V1,V4,V5), (V2,V3,V6)
402  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
403  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
404  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
405  ! 2 pour les composantes 2, 3 et 6
406  !------------------------------------------------------------------------
407  USE my_util
408  USE input_data
409  IMPLICIT NONE
410  INTEGER , INTENT(IN) :: type_op, mode, i_mode
411  REAL(KIND=8), INTENT(IN) :: visco, mass, stab, stab_bdy_ns
412  TYPE(petsc_csr_la) :: LA
413  TYPE(mesh_type), TARGET :: mesh
414  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
415  REAL(KIND=8) :: xij, viscolm
416  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
417  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
418  REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
419  REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
420  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
421  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
422  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
423  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
424  INTEGER :: ix, jx, iglob, jglob
425  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
426  REAL(KIND=8) :: hm, viscomode, hh
427  !=====DEUX INSERTIONS VALEURS A FAIRE....
428 !#include "petsc/finclude/petsc.h"
429  mat :: matrix
430  petscerrorcode :: ierr
431 
432  CALL matzeroentries (matrix, ierr)
433  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
434 
435  np = SIZE(mesh%rr,2)
436  n_w = mesh%gauss%n_w
437  n_ws = mesh%gauss%n_ws
438 
439  IF (type_op == 1) THEN
440  eps1 = 1.d0
441  eps2 = 1.d0
442  k_max = 2 !Structure vectorielle
443  ELSEIF (type_op == 2) THEN
444  eps1 = 1.d0
445  eps2 = -1.d0
446  k_max = 2 !Structure vectorielle
447  ELSEIF (type_op == 3) THEN
448  !cas du laplacien scalaire
449  eps1 = 0.d0
450  eps2 = 0.d0
451  k_max = 1 !Structure scalaire
452  ELSE
453  eps1 = 0.d0
454  eps2 = 0.d0
455  k_max = 1
456  CALL error_petsc('probleme de type d''operateur')
457  ENDIF
458 
459 
460  DO l = 1, mesh%gauss%l_G
461  DO ni = 1, mesh%gauss%n_w
462  DO nj = 1, mesh%gauss%n_w
463  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
464  END DO
465  END DO
466  END DO
467 
468  DO m = 1, mesh%me
469  jj_loc = mesh%jj(:,m)
470 
471  aij = 0.d0
472  bij = 0.d0
473  cij = 0.d0
474  DO l = 1, mesh%gauss%l_G
475 
476  !--------On calcule le rayon du point gauss
477  ray = 0
478  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
479  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
480  END DO
481  hh=mesh%hloc(m)
482  hm=min(mesh%hm(i_mode),hh)!WRONG choice
483  !hm=0.5d0/inputs%m_max
484  !hm=mesh%hm(i_mode) !(JLG April 7 2017)
485  viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m) ! We add the constant stabilization here
486  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
487 
488  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
489 
490  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
491 
492  !grad(u).grad(v) en r et z
493  xij = 0.d0
494  DO k = 1, mesh%gauss%k_d
495  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
496  END DO
497 
498  !blocs diagonaux
499  z = ray * viscolm* xij &
500  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
501  + viscomode*mode**2*wwprod(ni,nj,l)/ray
502  cij(ni,nj) = cij(ni,nj) + z
503  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
504  !blocs couplant
505  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
506  ENDDO
507  ENDDO
508 
509  ENDDO
510 
511  mat_loc = 0.d0
512  DO ki= 1, 3
513  DO ni = 1, n_w
514  i = jj_loc(ni)
515  iglob = la%loc_to_glob(ki,i)
516  ix = (ki-1)*n_w+ni
517  idxn(ix) = iglob-1
518  DO kj = 1, 3
519  DO nj = 1, n_w
520  j = jj_loc(nj)
521  jglob = la%loc_to_glob(kj,j)
522  jx = (kj-1)*n_w+nj
523  jdxn(jx) = jglob-1
524  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
525  IF (ki==kj) THEN
526  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
527  ELSE
528  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
529  END IF
530  ELSE ! ki=3 OR kj=3
531  IF (ki==kj) THEN
532  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
533  END IF
534  END IF
535  END DO
536  END DO
537  END DO
538  END DO
539  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
540 
541  !==Calcul de visco (grad u)T . (grad v)
542  aij = 0.d0
543  bij = 0.d0
544  cij = 0.d0
545  dij = 0.d0
546  eij = 0.d0
547  fij = 0.d0
548  DO l = 1, mesh%gauss%l_G
549 
550  !--------On calcule le rayon du point gauss
551  ray = 0
552  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
553  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
554  END DO
555 
556  viscolm = visco*mesh%gauss%rj(l,m)*ray
557 
558  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
559  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
560  aij(ni,nj) = aij(ni,nj) &
561  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2)
562  bij(ni,nj) = bij(ni,nj) &
563  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2)
564  cij(ni,nj) = cij(ni,nj) &
565  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m))
566  dij(ni,nj) = dij(ni,nj) &
567  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
568  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray)
569  eij(ni,nj) = eij(ni,nj) &
570  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray)
571  fij(ni,nj) = fij(ni,nj) &
572  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
573  END DO
574  END DO
575  END DO
576 
577  mat_loc=0.d0
578  idxn=0
579  jdxn=0
580  DO ni = 1, n_w
581  DO ki = 1, 3
582  i = jj_loc(ni)
583  iglob = la%loc_to_glob(ki,i)
584  ix = (ki-1)*n_w + ni
585  idxn(ix) = iglob-1
586  END DO
587  END DO
588  jdxn=idxn
589 
590  DO ni = 1, n_w
591  DO nj = 1, n_w
592  !=== Line i 1 (Vr)
593  ix = ni
594  !=== Column j 1 (Vr)
595  jx = nj
596  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
597  !=== Column j 2 (Vt)
598  jx = nj+n_w
599  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
600  !=== Column j 3 (Vz)
601  jx = nj+2*n_w
602  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
603 
604  !=== Line i 2 (Vt)
605  ix = ni+n_w
606  !=== Column j 1 (Vr)
607  jx = nj
608  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
609  !=== Column j 2 (Vt)
610  jx = nj+n_w
611  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
612  !=== Column j 3 (Vz)
613  jx = nj+2*n_w
614  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
615 
616  !=== Line i 3 (Vz)
617  ix = ni+2*n_w
618  !=== Column j 1 (Vr)
619  jx = nj
620  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
621  !=== Column j 2 (Vt)
622  jx = nj+n_w
623  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
624  !=== Column j 3 (Vz)
625  jx = nj+2*n_w
626  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
627  END DO
628  END DO
629  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
630  ENDDO
631  !== Fin du Calcul de visco (grad u)T . (grad v)
632 
633  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
634  !===Surface terms
635  stab_normal = coeff_stab_normal*(1.d0+visco)
636  !===Normalization is not done properly: (1.d0+visco)
637  !===To be revisited some day.
638  DO ms = 1, mesh%mes
639  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
640  aij = 0.d0
641  bij = 0.d0
642  cij = 0.d0
643  dij = 0.d0
644  DO ls = 1, mesh%gauss%l_Gs
645  !--------On calcule le rayon du point gauss
646  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
647  IF (ray.LT.1.d-10) cycle
648  hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
649  x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
650  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
651 
652  DO ni = 1, mesh%gauss%n_ws
653  DO nj = 1, mesh%gauss%n_ws
654  y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
655  aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
656  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
657  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
658  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
659  bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
660  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
661  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
662  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
663  cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
664  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
665  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
666  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
667  dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
668  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
669  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
670  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
671  END DO
672  END DO
673  END DO
674 
675  !++++++++++++++++
676  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
677  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
678  mat_loc_s = 0.d0
679  idxn_s = 0
680  jdxn_s = 0
681  DO ki = 1, 2
682  DO ni = 1, n_ws
683  i = mesh%jjs(ni,ms)
684  iglob = la%loc_to_glob(2*ki-1,i)
685  ix = (ki-1)*n_ws+ni
686  idxn_s(ix) = iglob-1
687  DO kj = 1, 2
688  DO nj = 1, n_ws
689  j = mesh%jjs(nj,ms)
690  jglob = la%loc_to_glob(2*kj-1,j)
691  jx = (kj-1)*n_ws+nj
692  jdxn_s(jx) = jglob-1
693  IF ((ki == 1) .AND. (kj == 1)) THEN
694  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
695  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
696  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
697  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
698  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
699  ELSE
700  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
701  END IF
702  END DO
703  END DO
704  END DO
705  END DO
706  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
707  END DO
708  END IF
709 
710  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
711  CALL matassemblyend(matrix,mat_final_assembly,ierr)
712 
713  END SUBROUTINE qs_diff_mass_vect_3x3
714 
715  SUBROUTINE qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
716  !========================================================================
717  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
718  ! (V1,V4,V5), (V2,V3,V6)
719  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
720  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
721  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
722  ! 2 pour les composantes 2, 3 et 6
723  !------------------------------------------------------------------------
724  USE my_util
725  USE input_data
726  IMPLICIT NONE
727 
728 
729  INTEGER , INTENT(IN) :: type_op, mode, i_mode
730  REAL(KIND=8), INTENT(IN) :: visco, mass, stab, stab_bdy_ns
731  TYPE(petsc_csr_la) :: LA
732  TYPE(mesh_type), TARGET :: mesh
733 
734 !!$ INTEGER, DIMENSION(:,:), POINTER :: jj
735 !!$ INTEGER, POINTER :: me
736 !!$ REAL(KIND=8), DIMENSION(:,:) , POINTER :: rr
737 
738  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
739  REAL(KIND=8) :: xij, viscolm, div_penal
740  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
741  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
742  REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
743  REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
744  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
745  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
746  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
747  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
748  INTEGER :: ix, jx, iglob, jglob
749  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
750  REAL(KIND=8) :: viscomode, hm
751  !=====DEUX INSERTIONS VALEURS A FAIRE....
752 !#include "petsc/finclude/petsc.h"
753  mat :: matrix
754  petscerrorcode :: ierr
755 
756  CALL matzeroentries (matrix, ierr)
757  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
758 
759  np = SIZE(mesh%rr,2)
760  n_w = mesh%gauss%n_w
761  n_ws = mesh%gauss%n_ws
762 
763  IF (type_op == 1) THEN
764  eps1 = 1.d0
765  eps2 = 1.d0
766  k_max = 2 !Structure vectorielle
767  ELSEIF (type_op == 2) THEN
768  eps1 = 1.d0
769  eps2 = -1.d0
770  k_max = 2 !Structure vectorielle
771  ELSEIF (type_op == 3) THEN
772  !cas du laplacien scalaire
773  eps1 = 0.d0
774  eps2 = 0.d0
775  k_max = 1 !Structure scalaire
776  ELSE
777  eps1 = 0.d0
778  eps2 = 0.d0
779  k_max = 1
780  CALL error_petsc('probleme de type d''operateur')
781  ENDIF
782 
783 
784  DO l = 1, mesh%gauss%l_G
785  DO ni = 1, mesh%gauss%n_w
786  DO nj = 1, mesh%gauss%n_w
787  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
788  END DO
789  END DO
790  END DO
791 
792  DO m = 1, mesh%me
793  jj_loc = mesh%jj(:,m)
794 
795  aij = 0.d0
796  bij = 0.d0
797  cij = 0.d0
798  DO l = 1, mesh%gauss%l_G
799 
800  !--------On calcule le rayon du point gauss
801  ray = 0
802  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
803  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
804  END DO
805 
806  viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m) ! We add the constant stabilization here
807  !hm=0.5d0/inputs%m_max
808  !hm=mesh%hm(i_mode) !(JLG April 7 2017)
809  hm=min(mesh%hm(i_mode),mesh%hloc(m))!WRONG choice
810  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
811 
812  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
813 
814  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
815 
816  !grad(u).grad(v) en r et z
817  xij = 0.d0
818  DO k = 1, mesh%gauss%k_d
819  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
820  END DO
821 
822  !blocs diagonaux
823  z = ray * viscolm* xij &
824  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
825  + viscomode*mode**2*wwprod(ni,nj,l)/ray
826  cij(ni,nj) = cij(ni,nj) + z
827  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
828  !blocs couplant
829  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
830  ENDDO
831  ENDDO
832 
833  ENDDO
834 
835 
836 
837  !++++++++++++++++
838 !!$ DO ki= 1, k_max
839 !!$ DO ni = 1, mesh%gauss%n_w
840 !!$ i = mesh%jj(ni, m)
841 !!$ ib = i + (ki-1)*np
842 !!$ DO kj = 1, k_max
843 !!$ DO nj = 1, mesh%gauss%n_w
844 !!$ j = mesh%jj(nj, m)
845 !!$ jb = j + (kj-1)*np
846 !!$ DO p = ia(ib), ia(ib+1) - 1
847 !!$ IF (ja(p) == jb) THEN
848 !!$ IF (ki==kj) THEN
849 !!$ a0(p) = a0(p) + aij(ni,nj)
850 !!$ ELSE
851 !!$ a0(p) = a0(p) + bij(ni,nj)
852 !!$ END IF
853 !!$ EXIT
854 !!$ ENDIF
855 !!$ END DO
856 !!$ END DO
857 !!$ END DO
858 !!$ END DO
859 !!$ END DO
860 !!$
861 !!$ IF (type_op /= 3) THEN
862 !!$ ki= 3
863 !!$ DO ni = 1, mesh%gauss%n_w
864 !!$ i = mesh%jj(ni, m)
865 !!$ ib = i + (ki-1)*np
866 !!$ kj = 3
867 !!$ DO nj = 1, mesh%gauss%n_w
868 !!$ j = mesh%jj(nj, m)
869 !!$ jb = j + (kj-1)*np
870 !!$ DO p = ia(ib), ia(ib+1) - 1
871 !!$ IF (ja(p) == jb) THEN
872 !!$ a0(p) = a0(p) + cij(ni,nj)
873 !!$ EXIT
874 !!$ ENDIF
875 !!$ END DO
876 !!$ END DO
877 !!$ END DO
878 !!$ END IF
879 
880  mat_loc = 0.d0
881  DO ki= 1, 3
882  DO ni = 1, n_w
883  i = jj_loc(ni)
884  iglob = la%loc_to_glob(ki,i)
885  ix = (ki-1)*n_w+ni
886  idxn(ix) = iglob-1
887  DO kj = 1, 3
888  DO nj = 1, n_w
889  j = jj_loc(nj)
890  jglob = la%loc_to_glob(kj,j)
891  jx = (kj-1)*n_w+nj
892  jdxn(jx) = jglob-1
893  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
894  IF (ki==kj) THEN
895  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
896  ELSE
897  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
898  END IF
899  ELSE ! ki=3 OR kj=3
900  IF (ki==kj) THEN
901  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
902  END IF
903  END IF
904  END DO
905  END DO
906  END DO
907  END DO
908  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
909 
910  !+++---------------------
911 !!$ IF (type_op==3) CYCLE
912  !==Calcul de visco (grad u)T . (grad v)
913  aij = 0.d0
914  bij = 0.d0
915  cij = 0.d0
916  dij = 0.d0
917  eij = 0.d0
918  fij = 0.d0
919  DO l = 1, mesh%gauss%l_G
920 
921  !--------On calcule le rayon du point gauss
922  ray = 0
923  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
924  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
925  END DO
926 
927  viscolm = visco*mesh%gauss%rj(l,m)*ray
928 !LC 2017/01/27
929  !div_penal = inputs%div_stab_in_ns*mesh%gauss%rj(l,m)*ray
930  div_penal = inputs%div_stab_in_ns/inputs%Re*mesh%gauss%rj(l,m)*ray
931 
932  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
933  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
934  aij(ni,nj) = aij(ni,nj) &
935  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
936  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
937  + mesh%gauss%ww(nj,l)/ray)
938  bij(ni,nj) = bij(ni,nj) &
939  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
940  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
941  cij(ni,nj) = cij(ni,nj) &
942  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
943  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
944  dij(ni,nj) = dij(ni,nj) &
945  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
946  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
947  + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
948  eij(ni,nj) = eij(ni,nj) &
949  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
950  + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
951  fij(ni,nj) = fij(ni,nj) &
952  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
953  + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
954  END DO
955  END DO
956  END DO
957  !aij = 0.d0
958  !bij = 0.d0
959  !cij = 0.d0
960  !dij = 0.d0
961  !eij = 0.d0
962  !fij = 0.d0
963  !++++++++++++++++
964  mat_loc=0.d0
965  idxn=0
966  jdxn=0
967  DO ni = 1, n_w
968  DO ki = 1, 3
969  i = jj_loc(ni)
970  iglob = la%loc_to_glob(ki,i)
971  ix = (ki-1)*n_w + ni
972  idxn(ix) = iglob-1
973  END DO
974  END DO
975  jdxn=idxn
976 
977  DO ni = 1, n_w
978  DO nj = 1, n_w
979  !=== Line i 1 (Vr)
980  ix = ni
981  !=== Column j 1 (Vr)
982  jx = nj
983  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
984  !=== Column j 2 (Vt)
985  jx = nj+n_w
986  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
987  !=== Column j 3 (Vz)
988  jx = nj+2*n_w
989  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
990 
991  !=== Line i 2 (Vt)
992  ix = ni+n_w
993  !=== Column j 1 (Vr)
994  jx = nj
995  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
996  !=== Column j 2 (Vt)
997  jx = nj+n_w
998  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
999  !=== Column j 3 (Vz)
1000  jx = nj+2*n_w
1001  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
1002 
1003  !=== Line i 3 (Vz)
1004  ix = ni+2*n_w
1005  !=== Column j 1 (Vr)
1006  jx = nj
1007  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
1008  !=== Column j 2 (Vt)
1009  jx = nj+n_w
1010  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
1011  !=== Column j 3 (Vz)
1012  jx = nj+2*n_w
1013  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
1014  END DO
1015  END DO
1016  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1017 
1018 !!$ DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1019 !!$ DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1020 !!$ ib = i
1021 !!$ jb = j
1022 !!$ DO p = ia(ib), ia(ib+1) - 1
1023 !!$ IF (ja(p) == jb) THEN
1024 !!$ a0(p) = a0(p) + aij(ni,nj)
1025 !!$ EXIT
1026 !!$ ENDIF
1027 !!$ END DO
1028 !!$ jb= j + np
1029 !!$ DO p = ia(ib), ia(ib+1) - 1
1030 !!$ IF (ja(p) == jb) THEN
1031 !!$ a0(p) = a0(p) + bij(ni,nj)
1032 !!$ EXIT
1033 !!$ ENDIF
1034 !!$ END DO
1035 !!$ jb = j + 2*np
1036 !!$ DO p = ia(ib), ia(ib+1) - 1
1037 !!$ IF (ja(p) == jb) THEN
1038 !!$ a0(p) = a0(p) + cij(ni,nj)
1039 !!$ EXIT
1040 !!$ ENDIF
1041 !!$ END DO
1042 !!$ ib = i + np
1043 !!$ jb = j
1044 !!$ DO p = ia(ib), ia(ib+1) - 1
1045 !!$ IF (ja(p) == jb) THEN
1046 !!$ a0(p) = a0(p) + bij(nj,ni)
1047 !!$ EXIT
1048 !!$ ENDIF
1049 !!$ END DO
1050 !!$ jb = j + np
1051 !!$ DO p = ia(ib), ia(ib+1) - 1
1052 !!$ IF (ja(p) == jb) THEN
1053 !!$ a0(p) = a0(p) + dij(ni,nj)
1054 !!$ EXIT
1055 !!$ ENDIF
1056 !!$ END DO
1057 !!$ jb = j + 2*np
1058 !!$ DO p = ia(ib), ia(ib+1) - 1
1059 !!$ IF (ja(p) == jb) THEN
1060 !!$ a0(p) = a0(p) + eij(ni,nj)
1061 !!$ EXIT
1062 !!$ ENDIF
1063 !!$ END DO
1064 !!$ ib = i + 2*np
1065 !!$ jb = j
1066 !!$ DO p = ia(ib), ia(ib+1) - 1
1067 !!$ IF (ja(p) == jb) THEN
1068 !!$ a0(p) = a0(p) + cij(nj,ni)
1069 !!$ EXIT
1070 !!$ ENDIF
1071 !!$ END DO
1072 !!$ jb = j + np
1073 !!$ DO p = ia(ib), ia(ib+1) - 1
1074 !!$ IF (ja(p) == jb) THEN
1075 !!$ a0(p) = a0(p) + eij(nj,ni)
1076 !!$ EXIT
1077 !!$ ENDIF
1078 !!$ END DO
1079 !!$ jb = j + 2*np
1080 !!$ DO p = ia(ib), ia(ib+1) - 1
1081 !!$ IF (ja(p) == jb) THEN
1082 !!$ a0(p) = a0(p) + fij(ni,nj)
1083 !!$ EXIT
1084 !!$ ENDIF
1085 !!$ END DO
1086 !!$ END DO
1087 !!$ END DO
1088  ENDDO
1089  !== Fin du Calcul de visco (grad u)T . (grad v)
1090  !++++++++++++++++------
1091 
1092  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
1093  !===Surface terms
1094  stab_normal = coeff_stab_normal*(1.d0+visco)
1095  DO ms = 1, mesh%mes
1096  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1097  aij = 0.d0
1098  bij = 0.d0
1099  cij = 0.d0
1100  dij = 0.d0
1101  DO ls = 1, mesh%gauss%l_Gs
1102  !--------On calcule le rayon du point gauss
1103  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1104  IF (ray.LT.1.d-10) cycle
1105  hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1106  x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1107  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1108 
1109  DO ni = 1, mesh%gauss%n_ws
1110  DO nj = 1, mesh%gauss%n_ws
1111  y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1112  aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1113  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1114  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1115  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1116  bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1117  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1118  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1119  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1120  cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1121  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1122  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1123  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1124  dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1125  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1126  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1127  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1128  END DO
1129  END DO
1130  END DO
1131 
1132  !++++++++++++++++
1133  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
1134  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
1135  mat_loc_s = 0.d0
1136  idxn_s = 0
1137  jdxn_s = 0
1138  DO ki = 1, 2
1139  DO ni = 1, n_ws
1140  i = mesh%jjs(ni,ms)
1141  iglob = la%loc_to_glob(2*ki-1,i)
1142  ix = (ki-1)*n_ws+ni
1143  idxn_s(ix) = iglob-1
1144  DO kj = 1, 2
1145  DO nj = 1, n_ws
1146  j = mesh%jjs(nj,ms)
1147  jglob = la%loc_to_glob(2*kj-1,j)
1148  jx = (kj-1)*n_ws+nj
1149  jdxn_s(jx) = jglob-1
1150  IF ((ki == 1) .AND. (kj == 1)) THEN
1151  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1152  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
1153  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1154  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
1155  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1156  ELSE
1157  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1158  END IF
1159  END DO
1160  END DO
1161  END DO
1162  END DO
1163  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1164 
1165 !!$ DO ki= 1, 3, 2
1166 !!$ DO ni = 1, mesh%gauss%n_ws
1167 !!$ i = mesh%jjs(ni, ms)
1168 !!$ ib = i + (ki-1)*np
1169 !!$ DO kj = 1, 3, 2
1170 !!$ DO nj = 1, mesh%gauss%n_ws
1171 !!$ j = mesh%jjs(nj, ms)
1172 !!$ jb = j + (kj-1)*np
1173 !!$ DO p = ia(ib), ia(ib+1) - 1
1174 !!$ IF (ja(p) == jb) THEN
1175 !!$ IF (ki == 1 .AND. kj == 1) THEN
1176 !!$ a0(p) = a0(p) + aij(ni,nj) + aij(nj,ni)
1177 !!$ ELSE IF (ki == 1 .AND. kj == 3) THEN
1178 !!$ a0(p) = a0(p) + bij(ni,nj) + cij(nj,ni)
1179 !!$ ELSE IF (ki == 3 .AND. kj == 1) THEN
1180 !!$ a0(p) = a0(p) + cij(ni,nj) + bij(nj,ni)
1181 !!$ ELSE
1182 !!$ a0(p) = a0(p) + dij(ni,nj) + dij(nj,ni)
1183 !!$ END IF
1184 !!$ EXIT
1185 !!$ ENDIF
1186 !!$ END DO
1187 !!$ END DO
1188 !!$ END DO
1189 !!$ END DO
1190 !!$ END DO
1191  !++++++++++++++++----
1192 
1193  END DO
1194  END IF
1195 
1196  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1197  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1198 
1199 
1200  END SUBROUTINE qs_diff_mass_vect_3x3_divpenal
1201 
1202 
1203  SUBROUTINE qs_diff_mass_vect_3x3_divpenal_art_comp(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, &
1204  stab_art_comp, i_mode, mode, matrix)
1205  !========================================================================
1206  ! laplacien vectoriel qui renvoie une matrice 3np * 3np, pour trois composantes
1207  ! (V1,V4,V5), (V2,V3,V6)
1208  ! calcul de l'operateur mass-visco(lap-eps1*1/r^2) et eps2*2m*visco/r^2
1209  ! eps1 et eps2 sont des parametres qui definissent le type d'operateur
1210  ! type de l'operateur : 1 pour les composantes 1, 4 et 5
1211  ! 2 pour les composantes 2, 3 et 6
1212  !------------------------------------------------------------------------
1213  USE my_util
1214  USE input_data
1215  IMPLICIT NONE
1216 
1217 
1218  INTEGER , INTENT(IN) :: type_op, mode, i_mode
1219  REAL(KIND=8), INTENT(IN) :: visco, mass, stab, stab_bdy_ns, stab_art_comp
1220  TYPE(petsc_csr_la) :: LA
1221  TYPE(mesh_type), TARGET :: mesh
1222  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
1223  REAL(KIND=8) :: xij, viscolm, div_penal
1224  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1225  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
1226  REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
1227  REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
1228  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
1229  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
1230  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
1231  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
1232  INTEGER :: ix, jx, iglob, jglob
1233  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1234  REAL(KIND=8) :: viscomode, hh, hm
1235  !=====DEUX INSERTIONS VALEURS A FAIRE....
1236 !#include "petsc/finclude/petsc.h"
1237  mat :: matrix
1238  petscerrorcode :: ierr
1239 
1240  CALL matzeroentries (matrix, ierr)
1241  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1242 
1243  np = SIZE(mesh%rr,2)
1244  n_w = mesh%gauss%n_w
1245  n_ws = mesh%gauss%n_ws
1246 
1247  IF (type_op == 1) THEN
1248  eps1 = 1.d0
1249  eps2 = 1.d0
1250  k_max = 2 !Structure vectorielle
1251  ELSEIF (type_op == 2) THEN
1252  eps1 = 1.d0
1253  eps2 = -1.d0
1254  k_max = 2 !Structure vectorielle
1255  ELSEIF (type_op == 3) THEN
1256  !cas du laplacien scalaire
1257  eps1 = 0.d0
1258  eps2 = 0.d0
1259  k_max = 1 !Structure scalaire
1260  ELSE
1261  eps1 = 0.d0
1262  eps2 = 0.d0
1263  k_max = 1
1264  CALL error_petsc('probleme de type d''operateur')
1265  ENDIF
1266 
1267 
1268  DO l = 1, mesh%gauss%l_G
1269  DO ni = 1, mesh%gauss%n_w
1270  DO nj = 1, mesh%gauss%n_w
1271  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1272  END DO
1273  END DO
1274  END DO
1275 
1276  DO m = 1, mesh%me
1277  jj_loc = mesh%jj(:,m)
1278 
1279  aij = 0.d0
1280  bij = 0.d0
1281  cij = 0.d0
1282  DO l = 1, mesh%gauss%l_G
1283 
1284  !--------On calcule le rayon du point gauss
1285  ray = 0
1286  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1287  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1288  END DO
1289 
1290  hh=mesh%hloc(m)
1291  hm=min(mesh%hm(i_mode),hh)!WRONG choice
1292  !hm=0.5d0/inputs%m_max
1293  !hm=mesh%hm(i_mode) !(JLG April 7 2017)
1294  viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m) ! We add the constant stabilization here
1295  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1296 
1297  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1298 
1299  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1300 
1301  !grad(u).grad(v) en r et z
1302  xij = 0.d0
1303  DO k = 1, mesh%gauss%k_d
1304  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
1305  END DO
1306 
1307  !blocs diagonaux
1308  z = ray * viscolm* xij &
1309  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1310  + viscomode*mode**2*wwprod(ni,nj,l)/ray
1311  cij(ni,nj) = cij(ni,nj) + z
1312  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
1313  !blocs couplant
1314  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
1315  ENDDO
1316  ENDDO
1317 
1318  ENDDO
1319 
1320  mat_loc = 0.d0
1321  DO ki= 1, 3
1322  DO ni = 1, n_w
1323  i = jj_loc(ni)
1324  iglob = la%loc_to_glob(ki,i)
1325  ix = (ki-1)*n_w+ni
1326  idxn(ix) = iglob-1
1327  DO kj = 1, 3
1328  DO nj = 1, n_w
1329  j = jj_loc(nj)
1330  jglob = la%loc_to_glob(kj,j)
1331  jx = (kj-1)*n_w+nj
1332  jdxn(jx) = jglob-1
1333  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
1334  IF (ki==kj) THEN
1335  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1336  ELSE
1337  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1338  END IF
1339  ELSE ! ki=3 OR kj=3
1340  IF (ki==kj) THEN
1341  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1342  END IF
1343  END IF
1344  END DO
1345  END DO
1346  END DO
1347  END DO
1348  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1349 
1350  !+++---------------------
1351  !==Calcul de visco (grad u)T . (grad v)
1352  aij = 0.d0
1353  bij = 0.d0
1354  cij = 0.d0
1355  dij = 0.d0
1356  eij = 0.d0
1357  fij = 0.d0
1358  DO l = 1, mesh%gauss%l_G
1359 
1360  !--------On calcule le rayon du point gauss
1361  ray = 0
1362  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1363  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1364  END DO
1365 
1366  viscolm = visco*mesh%gauss%rj(l,m)*ray
1367  div_penal = stab_art_comp*mesh%gauss%rj(l,m)*ray
1368  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1369  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1370  aij(ni,nj) = aij(ni,nj) &
1371  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
1372  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
1373  + mesh%gauss%ww(nj,l)/ray)
1374  bij(ni,nj) = bij(ni,nj) &
1375  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
1376  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
1377  cij(ni,nj) = cij(ni,nj) &
1378  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
1379  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
1380  dij(ni,nj) = dij(ni,nj) &
1381  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
1382  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
1383  + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
1384  eij(ni,nj) = eij(ni,nj) &
1385  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
1386  + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
1387  fij(ni,nj) = fij(ni,nj) &
1388  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
1389  + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
1390  END DO
1391  END DO
1392  END DO
1393 
1394  mat_loc=0.d0
1395  idxn=0
1396  jdxn=0
1397  DO ni = 1, n_w
1398  DO ki = 1, 3
1399  i = jj_loc(ni)
1400  iglob = la%loc_to_glob(ki,i)
1401  ix = (ki-1)*n_w + ni
1402  idxn(ix) = iglob-1
1403  END DO
1404  END DO
1405  jdxn=idxn
1406 
1407  DO ni = 1, n_w
1408  DO nj = 1, n_w
1409  !=== Line i 1 (Vr)
1410  ix = ni
1411  !=== Column j 1 (Vr)
1412  jx = nj
1413  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1414  !=== Column j 2 (Vt)
1415  jx = nj+n_w
1416  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1417  !=== Column j 3 (Vz)
1418  jx = nj+2*n_w
1419  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1420 
1421  !=== Line i 2 (Vt)
1422  ix = ni+n_w
1423  !=== Column j 1 (Vr)
1424  jx = nj
1425  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
1426  !=== Column j 2 (Vt)
1427  jx = nj+n_w
1428  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
1429  !=== Column j 3 (Vz)
1430  jx = nj+2*n_w
1431  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
1432 
1433  !=== Line i 3 (Vz)
1434  ix = ni+2*n_w
1435  !=== Column j 1 (Vr)
1436  jx = nj
1437  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
1438  !=== Column j 2 (Vt)
1439  jx = nj+n_w
1440  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
1441  !=== Column j 3 (Vz)
1442  jx = nj+2*n_w
1443  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
1444  END DO
1445  END DO
1446  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1447  ENDDO
1448  !== Fin du Calcul de visco (grad u)T . (grad v)
1449  !++++++++++++++++------
1450 
1451  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
1452  !===Surface terms
1453  stab_normal = coeff_stab_normal*(1.d0+visco)
1454  DO ms = 1, mesh%mes
1455  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1456  aij = 0.d0
1457  bij = 0.d0
1458  cij = 0.d0
1459  dij = 0.d0
1460  DO ls = 1, mesh%gauss%l_Gs
1461  !--------On calcule le rayon du point gauss
1462  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1463  IF (ray.LT.1.d-10) cycle
1464  hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1465  x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1466  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1467 
1468  DO ni = 1, mesh%gauss%n_ws
1469  DO nj = 1, mesh%gauss%n_ws
1470  y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1471  aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1472  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1473  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1474  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1475  bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1476  - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1477  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1478  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1479  cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1480  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1481  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1482  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1483  dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1484  - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1485  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1486  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1487  END DO
1488  END DO
1489  END DO
1490 
1491  !++++++++++++++++
1492  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
1493  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
1494  mat_loc_s = 0.d0
1495  idxn_s = 0
1496  jdxn_s = 0
1497  DO ki = 1, 2
1498  DO ni = 1, n_ws
1499  i = mesh%jjs(ni,ms)
1500  iglob = la%loc_to_glob(2*ki-1,i)
1501  ix = (ki-1)*n_ws+ni
1502  idxn_s(ix) = iglob-1
1503  DO kj = 1, 2
1504  DO nj = 1, n_ws
1505  j = mesh%jjs(nj,ms)
1506  jglob = la%loc_to_glob(2*kj-1,j)
1507  jx = (kj-1)*n_ws+nj
1508  jdxn_s(jx) = jglob-1
1509  IF ((ki == 1) .AND. (kj == 1)) THEN
1510  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1511  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
1512  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1513  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
1514  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1515  ELSE
1516  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1517  END IF
1518  END DO
1519  END DO
1520  END DO
1521  END DO
1522  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1523  END DO
1524  END IF
1525 
1526  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1527  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1528 
1529 
1531 
1532  SUBROUTINE qs_00_m (mesh, alpha, LA, matrix)
1533  !=================================================
1534  IMPLICIT NONE
1535  TYPE(mesh_type), INTENT(IN) :: mesh
1536  REAL(KIND=8), INTENT(IN) :: alpha
1537  TYPE(petsc_csr_la) :: LA
1538  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1539  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1540  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1541  INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1542  REAL(KIND=8) :: al, x, ray
1543 !#include "petsc/finclude/petsc.h"
1544  mat :: matrix
1545  petscerrorcode :: ierr
1546  CALL matzeroentries (matrix, ierr)
1547 
1548  DO m = 1, mesh%dom_me
1549  jj_loc = mesh%jj(:,m)
1550  mat_loc = 0.d0
1551  DO l = 1, mesh%gauss%l_G
1552  !Compute radius of Gauss point
1553  ray = 0
1554  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1555  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1556  END DO
1557 
1558  al = alpha *mesh%gauss%rj(l,m)*ray
1559  DO nj = 1, mesh%gauss%n_w;
1560  j = jj_loc(nj)
1561  jglob = la%loc_to_glob(1,j)
1562  idxn(nj) = jglob-1
1563  DO ni = 1, mesh%gauss%n_w;
1564  i = jj_loc(ni)
1565  iglob = la%loc_to_glob(1,i)
1566  idxm(ni) = iglob-1
1567  x = al*mesh%gauss%ww(nj,l)*mesh%gauss%ww(ni,l)
1568  mat_loc(ni,nj) = mat_loc(ni,nj) + x
1569  ENDDO
1570  ENDDO
1571  ENDDO
1572  CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1573  ENDDO
1574  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1575  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1576  END SUBROUTINE qs_00_m
1577 
1578 SUBROUTINE qs_11_m (mesh, visco, LA, matrix)
1579  !=================================================
1580  USE def_type_mesh
1581  IMPLICIT NONE
1582  TYPE(mesh_type), TARGET :: mesh
1583  REAL(KIND=8), INTENT(IN) :: visco
1584  TYPE(petsc_csr_la) :: LA
1585  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1586  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1587  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1588  INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1589  REAL(KIND=8) :: al, x, ray
1590 !#include "petsc/finclude/petsc.h"
1591  mat :: matrix
1592  petscerrorcode :: ierr
1593  CALL matzeroentries (matrix, ierr)
1594 
1595  DO m = 1, mesh%dom_me
1596  jj_loc = mesh%jj(:,m)
1597  mat_loc = 0.d0
1598 
1599  DO nj = 1, mesh%gauss%n_w;
1600  j = jj_loc(nj)
1601  jglob = la%loc_to_glob(1,j)
1602  idxn(nj) = jglob-1
1603  DO ni = 1, mesh%gauss%n_w;
1604  i = jj_loc(ni)
1605  iglob = la%loc_to_glob(1,i)
1606  idxm(ni) = iglob-1
1607  END DO
1608  END DO
1609 
1610  DO l = 1, mesh%gauss%l_G
1611  !Compute radius of Gauss point
1612  ray = 0
1613  DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1614  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1615  END DO
1616  al = visco * mesh%gauss%rj(l,m)
1617  DO nj = 1, mesh%gauss%n_w;
1618  DO ni = 1, mesh%gauss%n_w;
1619  x =al*sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1620  mat_loc(ni,nj) = mat_loc(ni,nj) + x
1621  ENDDO
1622  ENDDO
1623  ENDDO
1624  CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1625  ENDDO
1626  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1627  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1628  END SUBROUTINE qs_11_m
1629 
1630  SUBROUTINE qs_diff_mass_scal_m_level (mesh, LA, visco, mass, stab, i_mode, mode, matrix)
1631  !=================================================
1632  USE input_data
1633  IMPLICIT NONE
1634  TYPE(mesh_type), INTENT(IN) :: mesh
1635  TYPE(petsc_csr_la) :: LA
1636  REAL(KIND=8), INTENT(IN) :: visco, mass, stab
1637  INTEGER, INTENT(IN) :: mode, i_mode
1638  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1639  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1640  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1641  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
1642  REAL(KIND=8) :: ray
1643  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
1644  REAL(KIND=8) :: viscolm, xij, hm, viscomode, hh
1645 !#include "petsc/finclude/petsc.h"
1646  mat :: matrix
1647  petscerrorcode :: ierr
1648  CALL matzeroentries (matrix, ierr)
1649  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1650 
1651  DO l = 1, mesh%gauss%l_G
1652  DO ni = 1, mesh%gauss%n_w
1653  DO nj = 1, mesh%gauss%n_w
1654  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1655  END DO
1656  END DO
1657  END DO
1658 
1659  n_w = mesh%gauss%n_w
1660  DO m = 1, mesh%dom_me
1661  jj_loc = mesh%jj(:,m)
1662  a_loc = 0.d0
1663 
1664  DO nj = 1, mesh%gauss%n_w;
1665  j = jj_loc(nj)
1666  jglob = la%loc_to_glob(1,j)
1667  idxn(nj) = jglob-1
1668  DO ni = 1, mesh%gauss%n_w;
1669  i = jj_loc(ni)
1670  iglob = la%loc_to_glob(1,i)
1671  idxm(ni) = iglob-1
1672  END DO
1673  END DO
1674 
1675  DO l = 1, mesh%gauss%l_G
1676  !Compute radius of Gauss point
1677  ray = 0
1678  DO ni = 1, n_w; i = jj_loc(ni)
1679  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1680  END DO
1681  hh=mesh%hloc(m)
1682  hm=min(mesh%hm(i_mode),hh)!WRONG choice
1683  !hm=0.5d0/inputs%m_max
1684  !hm=mesh%hm(i_mode) !(JLG April 7 2017)
1685  viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
1686  viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1687  DO nj = 1, n_w
1688  DO ni = 1, n_w
1689  !grad(u).grad(v) w.r.t. r and z
1690  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1691  !start diagonal block
1692  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
1693  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1694  + viscomode*mode**2*wwprod(ni,nj,l)/ray
1695  !end diagonal block
1696  END DO
1697  END DO
1698  END DO
1699  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
1700  ENDDO
1701  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1702  CALL matassemblyend(matrix,mat_final_assembly,ierr)
1703 
1704  END SUBROUTINE qs_diff_mass_scal_m_level
1705 
1706 END MODULE fem_m_axi
1707 
1708 
1709 
1710 
1711 
1712 
1713 
1714 
1715 
1716 
1717 
1718 
1719 
1720 
subroutine qs_11_m(mesh, visco, LA, matrix)
Definition: fem_M_axi.f90:1579
subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, temp_list_robin_sides, convection_coeff, stab, mode, matrix)
Definition: fem_M_axi.f90:202
subroutine error_petsc(string)
Definition: my_util.f90:16
type(my_data), public inputs
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:130
subroutine qs_diff_mass_vect_m(type_op, LA, mesh, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:14
subroutine qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
Definition: fem_M_axi.f90:716
subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix)
Definition: fem_M_axi.f90:312
subroutine qs_00_m(mesh, alpha, LA, matrix)
Definition: fem_M_axi.f90:1533
subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
Definition: fem_M_axi.f90:399
subroutine qs_diff_mass_vect_3x3_divpenal_art_comp(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, stab_art_comp, i_mode, mode, matrix)
Definition: fem_M_axi.f90:1205
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, i_mode, mode, matrix)
Definition: fem_M_axi.f90:1631