SFEMaNS  version 4.1 (work in progress) 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 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
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
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
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
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
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
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 subroutine qs_diff_mass_vect_m(type_op, LA, mesh, visco, mass, stab, mode, matrix) Definition: fem_M_axi.f90:11 subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix) Definition: fem_M_axi.f90:127 subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, stab, mode, matrix) Definition: fem_M_axi.f90:198 subroutine qs_11_m(mesh, alpha, ia, ja, a0) subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix) Definition: fem_M_axi.f90:269 subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix) Definition: fem_M_axi.f90:356 subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, mode, matrix) Definition: fem_M_axi.f90:1250 subroutine qs_00_m(mesh, alpha, ia, ja, a0) subroutine error_petsc(string) Definition: my_util.f90:15 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
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