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