SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
sub_level_set.f90
Go to the documentation of this file.
1 !
2 !Authors Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyrights 2005
3 !Revised for PETSC, Jean-Luc Guermond, Franky Luddens, January 2011
4 !
5 MODULE subroutine_level_set
6  USE my_util
7  PUBLIC :: three_level_level_set
8  PRIVATE
9  REAL(KIND=8) :: max_velocity_at_tn
10 !===JLG Sept 27, 2016
11 ! LOGICAL :: compression_mthd_JLG=.TRUE.
12 !===JLG Sept 27, 2016
13 
14 CONTAINS
15 
16  SUBROUTINE three_level_level_set(comm_one_d,time, cc_1_LA, dt, list_mode, cc_mesh, cn_m1, cn, &
17  chmp_vit, max_vel, my_par_cc, cc_list_dirichlet_sides, cc_per, nb_inter, &
18  visc_entro_level)
19  !==============================
20  USE def_type_mesh
21  USE fem_m_axi
22  USE fem_rhs_axi
23  USE fem_tn_axi
24  USE dir_nodes_petsc
25  USE periodic
26  USE st_matrix
27  USE solve_petsc
28  USE dyn_line
29  USE boundary
30  USE chaine_caractere
31  USE sub_plot
32  USE st_matrix
33  USE sft_parallele
34  USE input_data
35  IMPLICIT NONE
36  REAL(KIND=8) :: time, dt
37  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
38  INTEGER, INTENT(IN) :: nb_inter
39  TYPE(mesh_type), INTENT(IN) :: cc_mesh
40  type(petsc_csr_la) :: cc_1_la
41  TYPE(periodic_type), INTENT(IN) :: cc_per
42  TYPE(solver_param), INTENT(IN) :: my_par_cc
43  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: cn_m1, cn
44  INTEGER, DIMENSION(:), INTENT(IN) :: cc_list_dirichlet_sides
45  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: chmp_vit
46  REAL(KIND=8), INTENT(INOUT) :: max_vel
47  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro_level
48  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: cc_global_d
49  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: cc_mode_global_js_d
50  LOGICAL, SAVE :: once = .true., once_vel=.true.
51  LOGICAL, SAVE :: re_init=.false.
52  INTEGER, SAVE :: m_max_c
53  INTEGER, DIMENSION(:), POINTER, SAVE :: cc_js_d ! Dirichlet nodes
54  INTEGER, SAVE :: my_petscworld_rank, nb_procs
55  REAL(KIND=8), SAVE :: les_coeff1_in_level
56  !----------FIN SAVE--------------------------------------------------------------------
57 
58  !----------Declaration sans save-------------------------------------------------------
59  INTEGER, POINTER, DIMENSION(:) :: cc_1_ifrom
60  INTEGER :: i, n
61  INTEGER :: code, mode
62  INTEGER :: bloc_size, m_max_pad
63  !allocations des variables locales
64  REAL(KIND=8), DIMENSION(cc_mesh%np) :: ff
65  REAL(KIND=8), DIMENSION(cc_mesh%np, 2) :: cn_p1
66  REAL(KIND=8), DIMENSION(cc_mesh%np,2,SIZE(list_mode)) :: cext
67  REAL(KIND=8), DIMENSION(cc_mesh%np,2,SIZE(list_mode)) :: cext_reg
68  REAL(KIND=8), DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_conv
69  REAL(KIND=8), DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_comp
70  REAL(KIND=8), DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 6, SIZE(list_mode)) :: ff_entro
71  REAL(KIND=8), DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_phi_1mphi
72  REAL(KIND=8) :: int_mass_correct
73  REAL(KIND=8) :: tps, tps_tot, tps_cumul
74  REAL(KIND=8) :: one, zero, three
75  DATA zero, one, three/0.d0,1.d0,3.d0/
76  !Communicators for Petsc, in space and Fourier------------------------------
77 #include "petsc/finclude/petsc.h"
78  petscerrorcode :: ierr
79  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
80  mat, DIMENSION(:), POINTER, SAVE :: cc_mat
81  vec, SAVE :: cb_1, cb_2, cx_1, cx_1_ghost
82  ksp, DIMENSION(:), POINTER, SAVE :: cc_ksp
83  vec, SAVE :: cb_reg_1, cb_reg_2 !vectors for level set regularization
84  mat, DIMENSION(:), POINTER, SAVE :: cc_reg_mat
85  ksp, DIMENSION(:), POINTER, SAVE :: cc_reg_ksp
86  !------------------------------END OF DECLARATION--------------------------------------
87 
88  IF (once) THEN
89 
90  once = .false.
91 
92  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
93  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
94 
95  !-----CREATE PETSC VECTORS AND GHOSTS-----------------------------------------
96  CALL create_my_ghost(cc_mesh,cc_1_la,cc_1_ifrom)
97  n = cc_mesh%dom_np
98  CALL veccreateghost(comm_one_d(1), n, &
99  petsc_determine, SIZE(cc_1_ifrom), cc_1_ifrom, cx_1, ierr)
100  CALL vecghostgetlocalform(cx_1, cx_1_ghost, ierr)
101  CALL vecduplicate(cx_1, cb_1, ierr)
102  CALL vecduplicate(cx_1, cb_2, ierr)
103  CALL vecduplicate(cx_1, cb_reg_1, ierr)
104  CALL vecduplicate(cx_1, cb_reg_2, ierr)
105  !------------------------------------------------------------------------------
106 
107  !-------------DIMENSIONS-------------------------------------------------------
108  m_max_c = SIZE(list_mode)
109  !------------------------------------------------------------------------------
110 
111  !---------PREPARE cc_js_D ARRAY FOR PHASE--------------------------------------
112  CALL dirichlet_nodes_parallel(cc_mesh, cc_list_dirichlet_sides, cc_js_d)
113  CALL scalar_glob_js_d(cc_mesh, list_mode, cc_1_la, cc_mode_global_js_d)
114  ALLOCATE(cc_global_d(m_max_c))
115  DO i = 1, m_max_c
116  ALLOCATE(cc_global_d(i)%DRL(SIZE(cc_mode_global_js_d(i)%DIL)))
117  END DO
118  !------------------------------------------------------------------------------
119 
120  !-------------ASSEMBLE PHASE MATRICES------------------------------------------
121  ALLOCATE(cc_mat(m_max_c),cc_ksp(m_max_c))
122 
123  IF (inputs%if_compression_mthd_JLG) THEN
124  ALLOCATE(cc_reg_mat(m_max_c),cc_reg_ksp(m_max_c))
125  DO i = 1, m_max_c
126  mode = list_mode(i)
127 
128  !---PHASE MATRIX for level set regularization
129  CALL create_local_petsc_matrix(comm_one_d(1), cc_1_la, cc_reg_mat(i), clean=.false.)
130  CALL qs_regul_m(cc_mesh, cc_1_la, 3.d0, mode, cc_reg_mat(i))
131  IF (cc_per%n_bord/=0) THEN
132  CALL periodic_matrix_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cc_reg_mat(i), cc_1_la)
133  END IF
134  CALL dirichlet_m_parallel(cc_reg_mat(i),cc_mode_global_js_d(i)%DIL)
135  CALL init_solver(my_par_cc,cc_reg_ksp(i),cc_reg_mat(i),comm_one_d(1),&
136  solver=my_par_cc%solver,precond=my_par_cc%precond)
137  ENDDO
138  END IF
139  max_velocity_at_tn = 0.d0
140 
141  IF (inputs%if_level_set_P2) THEN
142  les_coeff1_in_level=inputs%LES_coeff1
143  ELSE
144  les_coeff1_in_level=4.d0*inputs%LES_coeff1
145  END IF
146  ENDIF
147 
148  !===Assembling left hand side
149  IF (once_vel) THEN
150  re_init=.false.
151  once_vel = .false.
152 
153  ! Take care of once_vel
154  max_vel = max(1.1d0*max_velocity_at_tn,max_vel)
155  IF (my_petscworld_rank==0) WRITE(*,*) ' Recomputing matrix for level set function'
156  IF (my_petscworld_rank==0) WRITE(*,*) ' NEW MAX VEL test 1', time, max_vel
157 
158  DO i = 1, m_max_c
159  mode = list_mode(i)
160 
161  !---PHASE MATRIX
162  CALL create_local_petsc_matrix(comm_one_d(1), cc_1_la, cc_mat(i), clean=.false.)
163  IF (inputs%if_moment_bdf2) THEN
164  CALL qs_diff_mass_scal_m_level(cc_mesh, cc_1_la, 0.d0, 1.5d0/dt, &
165  les_coeff1_in_level, mode, cc_mat(i)) !===Coefficient 4 because P1 approximation
166  ELSE
167  CALL qs_diff_mass_scal_m_level(cc_mesh, cc_1_la, 0.d0, 1.d0/dt, &
168  les_coeff1_in_level, mode, cc_mat(i)) !===Coefficient 4 because P1 approximation
169  END IF
170  IF (cc_per%n_bord/=0) THEN
171  CALL periodic_matrix_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cc_mat(i), cc_1_la)
172  END IF
173  CALL dirichlet_m_parallel(cc_mat(i),cc_mode_global_js_d(i)%DIL)
174  CALL init_solver(my_par_cc,cc_ksp(i),cc_mat(i),comm_one_d(1),&
175  solver=my_par_cc%solver,precond=my_par_cc%precond, opt_re_init=re_init)
176  END DO
177  END IF
178  !===End assembling left hand side
179 
180  tps_tot = user_time()
181  tps_cumul = 0
182 
183  IF (inputs%if_moment_bdf2) THEN
184  cext = 2*cn-cn_m1
185  ELSE
186  cext = cn
187  END IF
188 
189  IF (inputs%if_compression_mthd_JLG) THEN
190  !---------------REGULARIZATION OF LEVEL SET FOR COMPRESSION---------------------------
191  DO i = 1, m_max_c
192  !===Compute rhs for level set regularization
193  CALL qs_00(cc_mesh,cc_1_la, cext(:,1,i), cb_reg_1)
194  CALL qs_00(cc_mesh,cc_1_la, cext(:,2,i), cb_reg_2)
195 
196  !===RHS periodicity
197  IF (cc_per%n_bord/=0) THEN
198  CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_reg_1, cc_1_la)
199  CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_reg_2, cc_1_la)
200  END IF
201 
202  !===RHS Dirichlet
203  n = SIZE(cc_js_d)
204  cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,1,cc_mesh%rr(:,cc_js_d), mode, time)
205  cc_global_d(i)%DRL(n+1:) = 0.d0
206  CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_reg_1)
207  cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,2,cc_mesh%rr(:,cc_js_d), mode, time)
208  cc_global_d(i)%DRL(n+1:) = 0.d0
209  CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_reg_2)
210 
211  !===Solve level set regularization equation
212  tps = user_time()
213  !Solve system cc_c
214  CALL solver(cc_reg_ksp(i),cb_reg_1,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
215  CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
216  CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
217  CALL extract(cx_1_ghost,1,1,cc_1_la,cext_reg(:,1,i))
218 
219  !Solve system cc_s
220  CALL solver(cc_reg_ksp(i),cb_reg_2,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
221  CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
222  CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
223  CALL extract(cx_1_ghost,1,1,cc_1_la,cext_reg(:,2,i))
224  tps = user_time() - tps; tps_cumul=tps_cumul+tps
225  !WRITE(*,*) ' Tps solution des pb de vitesse', tps, 'for mode ', mode
226 
227  END DO
228  !--------------- END REGULARIZATION OF LEVEL SET FOR COMPRESSION----------------------
229  END IF
230 
231  !===Compute rhs by FFT at Gauss points
232  CALL smb_ugradc_gauss_fft_par(comm_one_d(2), cc_mesh, list_mode, chmp_vit, cext, nb_procs, ff_conv)
233 
234  IF (inputs%if_compression_mthd_JLG) THEN
235  CALL smb_compr_visc_entro_gauss_fft_par(comm_one_d, cc_mesh, &
236  list_mode, cext, cext_reg, visc_entro_level,&
237  les_coeff1_in_level, nb_procs, ff_entro, ff_phi_1mphi)
238  ELSE
239  CALL smb_compression_gauss_fft_par(comm_one_d, cc_mesh, list_mode, &
240  chmp_vit, cext, nb_procs, ff_comp)
241  ff_conv=ff_conv-ff_comp
242 
243  CALL smb_visc_entro_gauss_fft_par(comm_one_d, cc_mesh, list_mode, cext, visc_entro_level,&
244  les_coeff1_in_level, nb_procs, ff_entro, ff_phi_1mphi)
245  END IF
246 
247  IF (inputs%if_mass_correction) THEN
248  CALL compute_int_mass_correct(comm_one_d, cc_mesh, list_mode, ff_conv, ff_phi_1mphi, int_mass_correct)
249  ELSE
250  int_mass_correct = 0.d0
251  END IF
252 
253  !------------DEBUT BOUCLE SUR LES MODES----------------
254  DO i = 1, m_max_c
255  mode = list_mode(i)
256 
257  !===RHS phase
258  IF (inputs%if_moment_bdf2) THEN
259  ff = (2/dt)*cn(:,1,i) - (1/(2*dt))*cn_m1(:,1,i) &
260  +source_in_level_set(nb_inter,1, cc_mesh%rr, mode, time)
261  CALL qs_00_level_set_gauss(cc_mesh, cc_1_la, ff, -ff_conv(:,1,i), &
262  mode, 1, cb_1, cext(:,1,i), &
263  ff_entro(:,:,i), -ff_phi_1mphi(:,1,i), int_mass_correct)
264 
265  ff = (2/dt)*cn(:,2,i) - (1/(2*dt))*cn_m1(:,2,i) &
266  +source_in_level_set(nb_inter,2, cc_mesh%rr, mode, time)
267  CALL qs_00_level_set_gauss(cc_mesh, cc_1_la, ff, -ff_conv(:,2,i),&
268  mode, 2, cb_2, cext(:,2,i), &
269  ff_entro(:,:,i), -ff_phi_1mphi(:,2,i), int_mass_correct)
270  ELSE !BDF1
271  ff = (1/dt)*cn(:,1,i) + source_in_level_set(nb_inter,1, cc_mesh%rr, mode, time)
272  CALL qs_00_level_set_gauss(cc_mesh, cc_1_la, ff, -ff_conv(:,1,i), &
273  mode, 1, cb_1, cext(:,1,i), &
274  ff_entro(:,:,i), -ff_phi_1mphi(:,1,i), int_mass_correct)
275 
276  ff = (1/dt)*cn(:,2,i) + source_in_level_set(nb_inter,2, cc_mesh%rr, mode, time)
277  CALL qs_00_level_set_gauss(cc_mesh, cc_1_la, ff, -ff_conv(:,2,i),&
278  mode, 2, cb_2, cext(:,2,i), &
279  ff_entro(:,:,i), -ff_phi_1mphi(:,2,i), int_mass_correct)
280  END IF
281 
282  !===RHS periodicity
283  IF (cc_per%n_bord/=0) THEN
284  CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_1, cc_1_la)
285  CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_2, cc_1_la)
286  END IF
287  !----------------------------------------------------------
288 
289  n = SIZE(cc_js_d)
290  cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,1,cc_mesh%rr(:,cc_js_d), mode, time)
291  cc_global_d(i)%DRL(n+1:) = 0.d0
292  CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_1)
293  cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,2,cc_mesh%rr(:,cc_js_d), mode, time)
294  cc_global_d(i)%DRL(n+1:) = 0.d0
295  CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_2)
296 
297  tps = user_time() - tps; tps_cumul=tps_cumul+tps
298  !WRITE(*,*) ' Tps second membre vitesse', tps
299  !-------------------------------------------------------------------------------------
300 
301  !--------------------INVERSION DES OPERATEURS--------------
302  tps = user_time()
303  !Solve system cc_c
304  CALL solver(cc_ksp(i),cb_1,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
305  CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
306  CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
307  CALL extract(cx_1_ghost,1,1,cc_1_la,cn_p1(:,1))
308 
309  !Solve system cc_s
310  CALL solver(cc_ksp(i),cb_2,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
311  CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
312  CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
313  CALL extract(cx_1_ghost,1,1,cc_1_la,cn_p1(:,2))
314  tps = user_time() - tps; tps_cumul=tps_cumul+tps
315  !WRITE(*,*) ' Tps solution des pb de vitesse', tps, 'for mode ', mode
316  !-------------------------------------------------------------------------------------
317 
318  !---------------UPDATES-----------------------
319  tps = user_time()
320 
321  IF (mode==0) THEN
322  cn_p1(:,2) = 0.d0
323  END IF
324 
325  cn_m1(:,:,i) = cn(:,:,i)
326  cn(:,:,i) = cn_p1
327 
328  tps = user_time() - tps; tps_cumul=tps_cumul+tps
329  !WRITE(*,*) ' Tps des updates', tps
330  !-------------------------------------------------------------------------------------
331  ENDDO
332 
333 !!$ IF (inputs%if_kill_overshoot) THEN
334 !!$ IF (nb_procs==1.AND.SIZE(list_mode)==1.AND.list_mode(1)==0) THEN !case axisym
335 !!$ cn = MIN(1.d0, cn)
336 !!$ cn = MAX(0.d0, cn)
337 !!$ ELSE !level set depends of theta
338 !!$ bloc_size = SIZE(cn,1)/nb_procs+1
339 !!$ m_max_pad = 3*SIZE(list_mode)*nb_procs/2
340 !!$ CALL FFT_NO_OVERSHOOT_LEVEL_SET(comm_one_d(2), cn, nb_procs, bloc_size, m_max_pad)
341 !!$ END IF
342 !!$ END IF
343 
344  tps_tot = user_time() - tps_tot
345  !WRITE(*,'(A,2(f13.3,2x))') ' Tps boucle en temps Navier_stokes', tps_tot, tps_cumul
346  !WRITE(*,*) ' TIME = ', time, '========================================'
347 
348  END SUBROUTINE three_level_level_set
349  !============================================
350 
351  SUBROUTINE smb_ugradc_gauss_fft_par(communicator,mesh,list_mode,V_in,c_in, nb_procs, c_out)
352  !=================================
353  USE gauss_points
354  USE sft_parallele
355  USE chaine_caractere
356  USE boundary
357  IMPLICIT NONE
358  TYPE(mesh_type), TARGET :: mesh
359  INTEGER, INTENT(IN) :: nb_procs
360  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
361  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v_in, c_in
362  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: c_out
363  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: gradc, w
364  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: div, cgauss
365 !!$ REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: cint
366  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
367  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
368  INTEGER :: m, l , i, mode, index, k
369  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vs
370  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: cs
371  INTEGER, DIMENSION(:,:), POINTER :: jj
372  INTEGER, POINTER :: me
373  REAL(KIND=8) :: ray, tps
374  REAL(KIND=8), DIMENSION(3) :: temps
375  INTEGER :: m_max_pad, bloc_size
376 #include "petsc/finclude/petsc.h"
377  mpi_comm :: communicator
378 
379  CALL gauss(mesh)
380  jj => mesh%jj
381  me => mesh%me
382 
383  tps = user_time()
384  DO i = 1, SIZE(list_mode)
385  mode = list_mode(i)
386  index = 0
387  DO m = 1, mesh%dom_me
388  j_loc = jj(:,m)
389  DO k = 1, 6
390  vs(:,k) = v_in(j_loc,k,i)
391  END DO
392  DO k = 1, 2
393  cs(:,k) = c_in(j_loc,k,i)
394  END DO
395  DO l = 1, l_g
396  index = index + 1
397  dw_loc = dw(:,:,l,m)
398 
399  !===Compute radius of Gauss point
400  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
401 
402  !-----------------vitesse sur les points de Gauss---------------------------
403  w(index,1,i) = sum(vs(:,1)*ww(:,l))
404  w(index,3,i) = sum(vs(:,3)*ww(:,l))
405  w(index,5,i) = sum(vs(:,5)*ww(:,l))
406 
407  w(index,2,i) = sum(vs(:,2)*ww(:,l))
408  w(index,4,i) = sum(vs(:,4)*ww(:,l))
409  w(index,6,i) = sum(vs(:,6)*ww(:,l))
410 
411  div(index,1,i) = sum(vs(:,1)*dw_loc(1,:)) + sum(vs(:,1)*ww(:,l))/ray &
412  + (mode/ray)*sum(vs(:,4)*ww(:,l)) + sum(vs(:,5)*dw_loc(2,:))
413  div(index,2,i) = sum(vs(:,2)*dw_loc(1,:)) + sum(vs(:,2)*ww(:,l))/ray &
414  - (mode/ray)*sum(vs(:,3)*ww(:,l)) + sum(vs(:,6)*dw_loc(2,:))
415 
416  !-----------------gradient de c sur les points de Gauss---------------------------
417  !coeff sur les cosinus et sinus
418  gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
419  gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
420  gradc(index,3,i) = mode/ray*sum(cs(:,2)*ww(:,l))
421  gradc(index,4,i) = -mode/ray*sum(cs(:,1)*ww(:,l))
422  gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
423  gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))
424 
425  cgauss(index,1,i) = sum(cs(:,1)*ww(:,l))
426  cgauss(index,2,i) = sum(cs(:,2)*ww(:,l))
427 
428  ENDDO
429  ENDDO
430  END DO
431 
432  !tps = user_time() - tps
433  !WRITE(*,*) ' Tps dans la grande boucle', tps
434  !tps = user_time()
435  temps = 0
436 
437  bloc_size = SIZE(gradc,1)/nb_procs+1
438  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
439  CALL fft_par_dot_prod_dcl(communicator, gradc, w, c_out, nb_procs, bloc_size, m_max_pad, temps)
440 !!$ !===Conservative term (phi div(u))
441 !!$ bloc_size = SIZE(Div,1)/nb_procs+1
442 !!$ CALL FFT_PAR_PROD_DCL(communicator, Div, Cgauss, cint, nb_procs, bloc_size, m_max_pad, temps)
443 !!$ c_out = c_out + cint
444 !!$ !===End Conservative term (phi div(u))
445  tps = user_time() - tps
446  !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
447  !write(*,*) ' Temps de Comm ', temps(1)
448  !write(*,*) ' Temps de Calc ', temps(2)
449  !write(*,*) ' Temps de Chan ', temps(3)
450 
451  END SUBROUTINE smb_ugradc_gauss_fft_par
452 
453  SUBROUTINE smb_compression_gauss_fft_par(communicator,mesh,list_mode,V_in,c_in, nb_procs, c_out)
454  !=================================
455  USE gauss_points
456  USE sft_parallele
457  USE chaine_caractere
458  USE boundary
459  USE tn_axi
460  USE input_data
461  IMPLICIT NONE
462  TYPE(mesh_type), TARGET :: mesh
463  INTEGER, INTENT(IN) :: nb_procs
464  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
465  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v_in, c_in
466  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: c_out
467  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: gradc, w
468  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: cgauss
469  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
470  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
471  INTEGER :: m, l , i, mode, index, k
472  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vs
473  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: cs
474  INTEGER, DIMENSION(:,:), POINTER :: jj
475  INTEGER, POINTER :: me
476  REAL(KIND=8) :: ray, tps, norm_vel_l2, volume_3d
477  INTEGER :: m_max_pad, bloc_size
478 #include "petsc/finclude/petsc.h"
479  mpi_comm, DIMENSION(2) :: communicator
480 
481  CALL gauss(mesh)
482  jj => mesh%jj
483  me => mesh%me
484 
485  tps = user_time()
486  DO i = 1, SIZE(list_mode)
487  mode = list_mode(i)
488  index = 0
489  DO m = 1, mesh%dom_me
490  j_loc = jj(:,m)
491  DO k = 1, 6
492  vs(:,k) = v_in(j_loc,k,i)
493  END DO
494  DO k = 1, 2
495  cs(:,k) = c_in(j_loc,k,i)
496  END DO
497  DO l = 1, l_g
498  index = index + 1
499  dw_loc = dw(:,:,l,m)
500 
501  !===Compute radius of Gauss point
502  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
503 
504  !-----------------vitesse sur les points de Gauss---------------------------
505  w(index,1,i) = sum(vs(:,1)*ww(:,l))
506  w(index,3,i) = sum(vs(:,3)*ww(:,l))
507  w(index,5,i) = sum(vs(:,5)*ww(:,l))
508 
509  w(index,2,i) = sum(vs(:,2)*ww(:,l))
510  w(index,4,i) = sum(vs(:,4)*ww(:,l))
511  w(index,6,i) = sum(vs(:,6)*ww(:,l))
512 
513  !-----------------gradient de c sur les points de Gauss---------------------------
514  !coeff sur les cosinus et sinus
515  gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
516  gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
517  gradc(index,3,i) = mode/ray*sum(cs(:,2)*ww(:,l))
518  gradc(index,4,i) = -mode/ray*sum(cs(:,1)*ww(:,l))
519  gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
520  gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))
521 
522  cgauss(index,1,i) = sum(cs(:,1)*ww(:,l))
523  cgauss(index,2,i) = sum(cs(:,2)*ww(:,l))
524  ENDDO
525  ENDDO
526  END DO
527 
528  bloc_size = mesh%gauss%l_G*mesh%me/nb_procs+1
529  bloc_size = mesh%gauss%l_G*(bloc_size/mesh%gauss%l_G)+mesh%gauss%l_G
530  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
531 
532 !TEST JLG LC vel L2
533  CALL twod_volume(communicator(1),mesh,volume_3d)
534  volume_3d = volume_3d*2*acos(-1.d0)
535  norm_vel_l2 = norm_sf(communicator, 'L2', mesh, list_mode, v_in)/volume_3d
536 !TEST JLG LC vel L2
537 
538  CALL fft_compression_level_set_dcl(communicator(2), communicator(1),gradc, w, cgauss, c_out, &
539  mesh%hloc_gauss, mesh%gauss%l_G, nb_procs, bloc_size, m_max_pad)
540 
541  tps = user_time() - tps
542  !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
543  !write(*,*) ' Temps de Comm ', temps(1)
544  !write(*,*) ' Temps de Calc ', temps(2)
545  !write(*,*) ' Temps de Chan ', temps(3)
546 
547  END SUBROUTINE smb_compression_gauss_fft_par
548 
549  SUBROUTINE twod_volume(communicator,mesh,RESLT)
550  !===========================
551  USE def_type_mesh
552  IMPLICIT NONE
553  TYPE(mesh_type) :: mesh
554  REAL(KIND=8), INTENT(OUT) :: reslt
555  REAL(KIND=8) :: vol_loc, vol_out
556  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
557  INTEGER :: m, l , i , ni, code
558  REAL(KIND=8) :: ray
559 #include "petsc/finclude/petsc.h"
560  mpi_comm :: communicator
561 
562  vol_loc = 0.d0
563  DO m = 1, mesh%dom_me
564  j_loc = mesh%jj(:,m)
565  DO l = 1, mesh%gauss%l_G
566  ray = 0
567  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
568  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
569  END DO
570  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
571  ENDDO
572  ENDDO
573  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
574  reslt = vol_out
575  END SUBROUTINE twod_volume
576 
577  SUBROUTINE qs_regul_m (mesh, LA, stab, mode, matrix)
578  !=================================================
579  USE def_type_mesh
580  USE input_data
581  IMPLICIT NONE
582  TYPE(mesh_type), INTENT(IN) :: mesh
583  type(petsc_csr_la) :: la
584  REAL(KIND=8), INTENT(IN) :: stab
585  INTEGER, INTENT(IN) :: mode
586  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
587  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm, idxn
588  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
589  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
590  REAL(KIND=8) :: ray
591  INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
592  REAL(KIND=8) :: viscolm, xij, viscomode, hm, hh
593 #include "petsc/finclude/petsc.h"
594  mat :: matrix
595  petscerrorcode :: ierr
596 
597  CALL matzeroentries(matrix, ierr)
598  CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
599 
600  DO l = 1, mesh%gauss%l_G
601  DO ni = 1, mesh%gauss%n_w
602  DO nj = 1, mesh%gauss%n_w
603  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
604  END DO
605  END DO
606  END DO
607 
608  n_w = mesh%gauss%n_w
609  DO m = 1, mesh%dom_me
610  jj_loc = mesh%jj(:,m)
611  a_loc = 0.d0
612 
613  DO nj = 1, mesh%gauss%n_w;
614  j = jj_loc(nj)
615  jglob = la%loc_to_glob(1,j)
616  idxn(nj) = jglob-1
617  DO ni = 1, mesh%gauss%n_w;
618  i = jj_loc(ni)
619  iglob = la%loc_to_glob(1,i)
620  idxm(ni) = iglob-1
621  END DO
622  END DO
623 
624  DO l = 1, mesh%gauss%l_G
625  !Compute radius of Gauss point
626  ray = 0
627  DO ni = 1, n_w; i = jj_loc(ni)
628  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
629  END DO
630  hh=mesh%hloc(m)
631 !!$ hm=MIN(0.5d0/inputs%m_max,hh)!WRONG choice
632  hm=0.5d0/inputs%m_max
633  viscolm = (stab*hh)**2*mesh%gauss%rj(l,m)
634  viscomode = (stab*hm)**2*mesh%gauss%rj(l,m)
635  DO nj = 1, n_w
636  DO ni = 1, n_w
637  !grad(u).grad(v) w.r.t. r and z
638  xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
639  !start diagonal block
640  a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
641  + ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
642  + viscomode*mode**2*wwprod(ni,nj,l)/ray
643  !end diagonal block
644  END DO
645  END DO
646  END DO
647  CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
648  ENDDO
649  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
650  CALL matassemblyend(matrix,mat_final_assembly,ierr)
651 
652  END SUBROUTINE qs_regul_m
653 
654  SUBROUTINE smb_compr_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, c_reg_in, visc_entro_real,&
655  coeff1_in_level, nb_procs, v_out, c_out)
656  !=================================
657  USE gauss_points
658  USE sft_parallele
659  USE chaine_caractere
660  USE boundary
661  USE input_data
662  IMPLICIT NONE
663  TYPE(mesh_type), TARGET :: mesh
664  INTEGER, INTENT(IN) :: nb_procs
665  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
666  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: c_in
667  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: c_reg_in
668  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro_real
669  REAL(KIND=8), INTENT(IN) :: coeff1_in_level
670  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: v_out
671  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)), INTENT(OUT) :: c_out
672  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: gradc_in, gradc_reg_in
673  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)) :: c_in_gauss
674  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
675  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
676  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: c_in_loc
677  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: c_reg_in_loc
678  INTEGER, DIMENSION(:,:), POINTER :: jj
679  INTEGER, POINTER :: me
680  REAL(KIND=8) :: ray, hh, hm
681  INTEGER :: m, l , i, mode, index, k
682  INTEGER :: m_max_pad, bloc_size
683 #include "petsc/finclude/petsc.h"
684  mpi_comm, DIMENSION(2) :: communicator
685 
686  CALL gauss(mesh)
687  jj => mesh%jj
688  me => mesh%me
689 
690  DO i = 1, SIZE(list_mode)
691  mode = list_mode(i)
692  index = 0
693  DO m = 1, mesh%dom_me
694  j_loc = jj(:,m)
695  DO k = 1, 2
696  c_in_loc(:,k) = c_in(j_loc,k,i)
697  c_reg_in_loc(:,k) = c_reg_in(j_loc,k,i)
698  END DO
699 
700  DO l = 1, l_g
701  index = index + 1
702  dw_loc = dw(:,:,l,m)
703 
704  !===Compute local mesh sizes
705  hh=mesh%hloc_gauss(index)
706 !!$ !hm=MIN(0.5d0/inputs%m_max,hh)!WRONG choice
707  hm=0.5d0/inputs%m_max
708 
709  !===Compute radius of Gauss point
710  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
711 
712  !-----------------c_in on Gauss points------------------
713  c_in_gauss(index,1,i) = sum(c_in_loc(:,1)*ww(:,l))
714  c_in_gauss(index,2,i) = sum(c_in_loc(:,2)*ww(:,l))
715 
716  !----------Gradient of c_in on Gauss points---------
717  gradc_in(index,1,i) = sum(c_in_loc(:,1)*dw_loc(1,:))*hh
718  gradc_in(index,2,i) = sum(c_in_loc(:,2)*dw_loc(1,:))*hh
719  gradc_in(index,3,i) = mode/ray*sum(c_in_loc(:,2)*ww(:,l))*hm
720  gradc_in(index,4,i) = -mode/ray*sum(c_in_loc(:,1)*ww(:,l))*hm
721  gradc_in(index,5,i) = sum(c_in_loc(:,1)*dw_loc(2,:))*hh
722  gradc_in(index,6,i) = sum(c_in_loc(:,2)*dw_loc(2,:))*hh
723 
724  !----------Gradient of c_reg_in on Gauss points---------
725  gradc_reg_in(index,1,i) = sum(c_reg_in_loc(:,1)*dw_loc(1,:))*hh
726  gradc_reg_in(index,2,i) = sum(c_reg_in_loc(:,2)*dw_loc(1,:))*hh
727  gradc_reg_in(index,3,i) = mode/ray*sum(c_reg_in_loc(:,2)*ww(:,l))*hm
728  gradc_reg_in(index,4,i) = -mode/ray*sum(c_reg_in_loc(:,1)*ww(:,l))*hm
729  gradc_reg_in(index,5,i) = sum(c_reg_in_loc(:,1)*dw_loc(2,:))*hh
730  gradc_reg_in(index,6,i) = sum(c_reg_in_loc(:,2)*dw_loc(2,:))*hh
731  END DO
732  END DO
733  END DO
734  bloc_size = SIZE(c_in_gauss,1)/nb_procs+1
735  bloc_size = l_g*(bloc_size/l_g)+l_g
736  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
737 
738  !===Compute c_out = max(0,phi*(1-phi))
739  !===Compute V_out = (coeff1_in_level-visc_entro)*Grad(phi) + comp*visc_entro/(h_loc*|Grad(phi_reg)(x,t)|)*c_out*Grad(phi_reg)
740  CALL fft_par_compr_entro_visc_dcl(communicator(2), gradc_in, gradc_reg_in, c_in_gauss, visc_entro_real, &
741  mesh%hloc_gauss, coeff1_in_level, v_out, c_out, nb_procs, bloc_size, m_max_pad)
742 
743  END SUBROUTINE smb_compr_visc_entro_gauss_fft_par
744 
745  SUBROUTINE smb_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, visc_entro_real,&
746  coeff1_in_level, nb_procs, v_out, c_out)
747  !=================================
748  USE gauss_points
749  USE sft_parallele
750  USE chaine_caractere
751  USE boundary
752  USE input_data
753  IMPLICIT NONE
754  TYPE(mesh_type), TARGET :: mesh
755  INTEGER, INTENT(IN) :: nb_procs
756  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
757  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: c_in
758  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro_real
759  REAL(KIND=8), INTENT(IN) :: coeff1_in_level
760  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: v_out
761  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)), INTENT(OUT) :: c_out
762  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: gradc_in
763  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)) :: c_in_gauss
764  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
765  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
766  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: c_in_loc
767  INTEGER, DIMENSION(:,:), POINTER :: jj
768  INTEGER, POINTER :: me
769  REAL(KIND=8) :: ray, hh, hm
770  INTEGER :: m, l , i, mode, index, k
771  INTEGER :: m_max_pad, bloc_size
772 #include "petsc/finclude/petsc.h"
773  mpi_comm, DIMENSION(2) :: communicator
774 
775  CALL gauss(mesh)
776  jj => mesh%jj
777  me => mesh%me
778 
779  DO i = 1, SIZE(list_mode)
780  mode = list_mode(i)
781  index = 0
782  DO m = 1, mesh%dom_me
783  j_loc = jj(:,m)
784  DO k = 1, 2
785  c_in_loc(:,k) = c_in(j_loc,k,i)
786  END DO
787 
788  DO l = 1, l_g
789  index = index + 1
790  dw_loc = dw(:,:,l,m)
791 
792  !===Compute local mesh sizes
793  hh=mesh%hloc_gauss(index)
794 !!$ !hm=MIN(0.5d0/inputs%m_max,hh)!WRONG choice
795  hm=0.5d0/inputs%m_max
796 
797  !===Compute radius of Gauss point
798  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
799 
800  !-----------------c_in on Gauss points------------------
801  c_in_gauss(index,1,i) = sum(c_in_loc(:,1)*ww(:,l))
802  c_in_gauss(index,2,i) = sum(c_in_loc(:,2)*ww(:,l))
803 
804  !----------Gradient of c_in on Gauss points---------
805  gradc_in(index,1,i) = sum(c_in_loc(:,1)*dw_loc(1,:))*hh
806  gradc_in(index,2,i) = sum(c_in_loc(:,2)*dw_loc(1,:))*hh
807  gradc_in(index,3,i) = mode/ray*sum(c_in_loc(:,2)*ww(:,l))*hm
808  gradc_in(index,4,i) = -mode/ray*sum(c_in_loc(:,1)*ww(:,l))*hm
809  gradc_in(index,5,i) = sum(c_in_loc(:,1)*dw_loc(2,:))*hh
810  gradc_in(index,6,i) = sum(c_in_loc(:,2)*dw_loc(2,:))*hh
811  END DO
812  END DO
813  END DO
814  !CALL MPI_COMM_SIZE(communicator, nb_procs, code)
815  bloc_size = SIZE(c_in_gauss,1)/nb_procs+1
816  bloc_size = l_g*(bloc_size/l_g)+l_g
817  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
818 
819  !===Compute c_out = max(0,phi*(1-phi))
820  !===Compute V_out = (coeff1_in_level-visc_entro)*Grad(phi)
821  CALL fft_par_entro_visc_dcl(communicator(2), gradc_in, c_in_gauss, visc_entro_real, &
822  mesh%hloc_gauss, coeff1_in_level, v_out, c_out, nb_procs, bloc_size, m_max_pad)
823 
824  END SUBROUTINE smb_visc_entro_gauss_fft_par
825 
826  SUBROUTINE compute_int_mass_correct(communicator, mesh, list_mode, c1_in, c2_in, int_out)
827  !=================================
828  USE gauss_points
829  IMPLICIT NONE
830  TYPE(mesh_type), TARGET :: mesh
831  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
832  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: c1_in !=ff_conv
833  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: c2_in !=ff_phi_1mphi
834  REAL(KIND=8), INTENT(OUT) :: int_out
835  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
836  REAL(KIND=8) :: ray
837  REAL(KIND=8) :: int_c1_in, int_c1_in_loc, int_c1_in_f
838  REAL(KIND=8) :: int_c2_in, int_c2_in_loc, int_c2_in_f
839  INTEGER :: m, l , i, index, code
840 #include "petsc/finclude/petsc.h"
841  mpi_comm, DIMENSION(2) :: communicator
842 
843  int_c1_in_loc = 0.d0
844  int_c2_in_loc = 0.d0
845 
846  DO i = 1, SIZE(list_mode)
847  index = 0
848  IF (list_mode(i)==0) THEN
849  DO m = 1, mesh%dom_me
850  j_loc = mesh%jj(:,m)
851  DO l = 1, mesh%gauss%l_G
852  index = index + 1
853  !===Compute radius of Gauss point
854  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
855 
856  !===Compute integral of c1_in and c2_in
857  int_c1_in_loc = int_c1_in_loc + c1_in(index,1,i)*ray*mesh%gauss%rj(l,m)
858  int_c2_in_loc = int_c2_in_loc + c2_in(index,1,i)*ray*mesh%gauss%rj(l,m)
859  END DO
860  END DO
861  END IF
862  END DO
863 
864  int_c1_in_loc = int_c1_in_loc*2*acos(-1.d0)
865  CALL mpi_allreduce(int_c1_in_loc, int_c1_in_f, 1, mpi_double_precision, mpi_sum, &
866  communicator(2), code)
867  CALL mpi_allreduce(int_c1_in_f, int_c1_in, 1, mpi_double_precision, mpi_sum, &
868  communicator(1), code)
869 
870  int_c2_in_loc = int_c2_in_loc*2*acos(-1.d0)
871  CALL mpi_allreduce(int_c2_in_loc, int_c2_in_f, 1, mpi_double_precision, mpi_sum, &
872  communicator(2), code)
873  CALL mpi_allreduce(int_c2_in_f, int_c2_in, 1, mpi_double_precision, mpi_sum, &
874  communicator(1), code)
875 
876  IF (int_c2_in.LE.1.d-14) THEN
877  int_out=0.d0
878  ELSE
879  int_out=-int_c1_in/int_c2_in
880  END IF
881 
882  END SUBROUTINE compute_int_mass_correct
883 
884  SUBROUTINE qs_00_level_set_gauss (mesh, LA, ff, ff_gauss, mode, type, vect, level_set_ext, &
885  fcompr, ff_phi_1mphi, stab_mass)
886  !=================================
887  USE def_type_mesh
888  USE input_data
889  IMPLICIT NONE
890  TYPE(mesh_type), TARGET :: mesh
891  TYPE(petsc_csr_la) :: la
892  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff, ff_gauss
893  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: level_set_ext
894  INTEGER , INTENT(IN) :: mode
895  INTEGER , INTENT(IN) :: type ! 1 = cosine, 2 = sine
896  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: fcompr
897  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff_phi_1mphi
898  REAL(KIND=8), INTENT(IN) :: stab_mass
899  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: ff_loc, level_set_loc
900  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
901  REAL(KIND=8), DIMENSION(3) :: fcomprl
902  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v_loc
903  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm
904  INTEGER :: i, m, l, ni, iglob, index
905  REAL(KIND=8) :: fl, ray
906 #include "petsc/finclude/petsc.h"
907  vec :: vect
908  petscerrorcode :: ierr
909 
910  CALL vecset(vect, 0.d0, ierr)
911 
912  index = 0
913  DO m = 1, mesh%dom_me
914  jj_loc = mesh%jj(:,m)
915  ff_loc = ff(jj_loc)
916  level_set_loc = level_set_ext(jj_loc)
917 
918  DO ni = 1, mesh%gauss%n_w
919  i = mesh%jj(ni,m)
920  iglob = la%loc_to_glob(1,i)
921  idxm(ni) = iglob-1
922  END DO
923 
924  v_loc = 0.d0
925  DO l = 1, mesh%gauss%l_G
926  index = index + 1
927  ray = 0
928  DO ni = 1, mesh%gauss%n_w
929  i = jj_loc(ni)
930  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
931  END DO
932 
933  ! Compute ff on gauss points + ff_gauss
934  fl = (sum(ff_loc(:)*mesh%gauss%ww(:,l)) + ff_gauss(index)+stab_mass*ff_phi_1mphi(index))*mesh%gauss%rj(l,m)*ray
935 
936  ! Compute compressive term on gauss points
937  IF (type==1) THEN
938  fcomprl(1) = fcompr(index,1)*mesh%gauss%rj(l,m)*ray
939  fcomprl(2) = -mode*fcompr(index,4)*mesh%gauss%rj(l,m)
940  fcomprl(3) = fcompr(index,5)*mesh%gauss%rj(l,m)*ray
941 
942  ELSE IF (type==2) THEN
943  fcomprl(1) = fcompr(index,2)*mesh%gauss%rj(l,m)*ray
944  fcomprl(2) = mode*fcompr(index,3)*mesh%gauss%rj(l,m)
945  fcomprl(3) = fcompr(index,6)*mesh%gauss%rj(l,m)*ray
946  ELSE
947  CALL error_petsc('error in type while calling qs_00_level_set_gauss')
948  END IF
949 
950  DO ni = 1, mesh%gauss%n_w
951  ! Add time derivative, advection and forcing term
952  v_loc(ni) = v_loc(ni) + mesh%gauss%ww(ni,l) *fl
953  ! Add fcompr=(c1_LES-nu_entro)*Grad(phi) + compression
954  v_loc(ni) = v_loc(ni) + (fcomprl(1)*mesh%gauss%dw(1,ni,l,m) + fcomprl(2)*mesh%gauss%ww(ni,l) &
955  + fcomprl(3)*mesh%gauss%dw(2,ni,l,m))
956  END DO
957  END DO
958 
959  CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
960  ENDDO
961  CALL vecassemblybegin(vect,ierr)
962  CALL vecassemblyend(vect,ierr)
963  END SUBROUTINE qs_00_level_set_gauss
964 
965 END MODULE subroutine_level_set
tn_axi
Definition: tn_axi.f90:5
sft_parallele::fft_par_dot_prod_dcl
subroutine, public fft_par_dot_prod_dcl(communicator, V1_in, V2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:327
fem_m_axi
Definition: fem_M_axi.f90:6
fem_rhs_axi
Definition: fem_rhs_axi.f90:1
solve_petsc
Definition: solver.f90:1
sft_parallele
Definition: fft_parallel.f90:4
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
st_matrix::create_my_ghost
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
gauss_points
Definition: associate_gauss.f90:4
sub_plot
Definition: sub_plot.f90:4
dir_nodes_petsc::dirichlet_rhs
subroutine dirichlet_rhs(js_D, bs_D, b)
Definition: dir_nodes_petsc.f90:170
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
subroutine_level_set::three_level_level_set
subroutine, public three_level_level_set(comm_one_d, time, cc_1_LA, dt, list_mode, cc_mesh, cn_m1, cn, chmp_vit, max_vel, my_par_cc, cc_list_dirichlet_sides, cc_per, nb_inter, visc_entro_level)
Definition: sub_level_set.f90:16
st_matrix
Definition: st_csr.f90:4
boundary_test_19::source_in_level_set
real(kind=8) function, dimension(size(rr, 2)), public source_in_level_set(interface_nb, TYPE, rr, m, t)
Definition: condlim_test_19.f90:170
dir_nodes_petsc::scalar_glob_js_d
subroutine scalar_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_mode_global_js_D)
Definition: dir_nodes_petsc.f90:312
solve_petsc::init_solver
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
sft_parallele::fft_par_entro_visc_dcl
subroutine, public fft_par_entro_visc_dcl(communicator, V1_in, c1_in, c_in_real, hloc_gauss, coeff1_in_level, V_out, c_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:3677
subroutine_level_set::smb_compr_visc_entro_gauss_fft_par
subroutine smb_compr_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, c_reg_in, visc_entro_real, coeff1_in_level, nb_procs, V_out, c_out)
Definition: sub_level_set.f90:654
dir_nodes_petsc::dirichlet_nodes_parallel
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
Definition: dir_nodes_petsc.f90:52
subroutine_level_set::compute_int_mass_correct
subroutine compute_int_mass_correct(communicator, mesh, list_mode, c1_in, c2_in, int_out)
Definition: sub_level_set.f90:826
subroutine_level_set
Definition: sub_level_set.f90:5
dir_nodes_petsc
Definition: dir_nodes_petsc.f90:4
dyn_line
Definition: dyn_line.f90:4
solve_petsc::solver
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
dir_nodes_petsc::dirichlet_m_parallel
subroutine dirichlet_m_parallel(matrix, glob_js_D)
Definition: dir_nodes_petsc.f90:149
dyn_line::dyn_int_line
Definition: dyn_line.f90:5
sft_parallele::fft_par_compr_entro_visc_dcl
subroutine, public fft_par_compr_entro_visc_dcl(communicator, V1_in, V2_in, c1_in, c_in_real, hloc_gauss, coeff1_in_level, V_out, c_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:2775
input_data
Definition: read_my_data.f90:200
def_type_mesh::periodic_type
Definition: def_type_mesh.f90:54
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:23
sft_parallele::fft_compression_level_set_dcl
subroutine, public fft_compression_level_set_dcl(communicator_F, communicator_S, V1_in, V2_in, c_in, c_out, hloc_gauss, l_G, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:3436
chaine_caractere
Definition: chaine_caractere_par.f90:4
fem_tn_axi
Definition: fem_tn_axi.f90:4
fem_s_axi::qs_00
subroutine qs_00(mesh, ff, u0)
Definition: fem_s_direct_axi.f90:127
tn_axi::norm_sf
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:38
subroutine_level_set::qs_regul_m
subroutine qs_regul_m(mesh, LA, stab, mode, matrix)
Definition: sub_level_set.f90:577
fem_m_axi::qs_diff_mass_scal_m_level
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:1250
my_util::user_time
real(kind=8) function user_time()
Definition: my_util.f90:7
periodic::periodic_rhs_petsc
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
Definition: prep_mesh_periodic.f90:455
solve_petsc::solver_param
Definition: solver.f90:3
subroutine_level_set::qs_00_level_set_gauss
subroutine qs_00_level_set_gauss(mesh, LA, ff, ff_gauss, mode, type, vect, level_set_ext, fcompr, ff_phi_1mphi, stab_mass)
Definition: sub_level_set.f90:884
subroutine_level_set::smb_visc_entro_gauss_fft_par
subroutine smb_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, visc_entro_real, coeff1_in_level, nb_procs, V_out, c_out)
Definition: sub_level_set.f90:745
dyn_line::dyn_real_line
Definition: dyn_line.f90:9
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
periodic::periodic_matrix_petsc
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
Definition: prep_mesh_periodic.f90:377
subroutine_level_set::smb_compression_gauss_fft_par
subroutine smb_compression_gauss_fft_par(communicator, mesh, list_mode, V_in, c_in, nb_procs, c_out)
Definition: sub_level_set.f90:453
st_matrix::extract
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:33
boundary_test_19::level_set_exact
real(kind=8) function, dimension(size(rr, 2)), public level_set_exact(interface_nb, TYPE, rr, m, t)
Definition: condlim_test_19.f90:260
solve_petsc::create_local_petsc_matrix
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:142
entropy_viscosity::twod_volume
subroutine twod_volume(communicator, mesh, RESLT)
Definition: entropy_viscosity.f90:1014
subroutine_level_set::smb_ugradc_gauss_fft_par
subroutine smb_ugradc_gauss_fft_par(communicator, mesh, list_mode, V_in, c_in, nb_procs, c_out)
Definition: sub_level_set.f90:351