SFEMaNS  version 5.3
Reference documentation for SFEMaNS
entropy_viscosity.f90
Go to the documentation of this file.
1 MODULE entropy_viscosity
2  USE my_util
3  PUBLIC :: compute_entropy_viscosity, compute_entropy_viscosity_mom, &
4  compute_entropy_viscosity_mom_no_level_set
5  PRIVATE
6 CONTAINS
7  SUBROUTINE compute_entropy_viscosity(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, vvz_per, &
8  un, un_m1, un_m2, pn_m1, rotv_v_m1, visco_entro_grad_u, opt_tempn)
9  USE def_type_mesh
10  USE fem_m_axi
11  USE solve_petsc
12  USE periodic
13  USE rhs_gauss_computing
14  USE rhs_para_assembling
15  USE dir_nodes_petsc
16  USE st_matrix
17  USE sft_parallele
18  USE tn_axi
19  USE input_data
20  USE my_util
21  USE boundary
22  USE fourier_to_real_for_vtu
23 #include "petsc/finclude/petsc.h"
24  USE petsc
25  IMPLICIT NONE
26  TYPE(petsc_csr_la) :: vv_3_LA
27  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
28  TYPE(periodic_type), INTENT(IN) :: vvz_per
29  REAL(KIND=8), INTENT(IN) :: time
30  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
31  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: un, un_m1, un_m2
32  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: pn_m1
33  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotv_v_m1
34  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: opt_tempn
35  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(OUT) :: visco_entro_grad_u
36  TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_D
37  LOGICAL, SAVE :: once = .true.
38  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: zero_dir1
39  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: zero_dir2
40  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: zero_dir3
41  REAL(KIND=8), SAVE :: Volume_3D
42  INTEGER, SAVE :: iteration
43 
44  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: strain_rate_tensor
45  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_Gs*vv_mesh%dom_mes,6,SIZE(list_mode)) :: strain_rate_tensor_scal_n_bdy
46  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
47  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: un_m1_gauss, res_ns_gauss
48  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: res_ns
49  INTEGER, POINTER, DIMENSION(:) :: vv_3_ifrom
50  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
51  REAL(KIND=8) :: norm_vel_L2
52  !REAL(KIND=8) :: norm_res_ns_L2
53  !INTEGER :: rank
54  INTEGER :: i, k, nu_mat, mode, m, l, TYPE, index, n
55  INTEGER :: nb_procs, code, bloc_size, m_max_pad
56  petscerrorcode :: ierr
57  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
58  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
59  mat, DIMENSION(:), POINTER, SAVE :: les_mat
60  ksp, DIMENSION(:), POINTER, SAVE :: les_ksp
61 
62  IF (once) THEN
63 
64  once = .false.
65 
66  !===CREATE PETSC VECTORS AND GHOSTS
67  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
68  n = 3*vv_mesh%dom_np
69  CALL veccreateghost(comm_one_d(1), n, &
70  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
71  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
72  CALL vecduplicate(vx_3, vb_3_145, ierr)
73  CALL vecduplicate(vx_3, vb_3_236, ierr)
74 
75  !===Compute Volume
76  CALL twod_volume(comm_one_d(1),vv_mesh,volume_3d)
77  volume_3d = volume_3d*2*acos(-1.d0)
78 
79  !===PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
80  DO k = 1, 3
81  CALL dirichlet_nodes_parallel(vv_mesh, inputs%vv_list_dirichlet_sides(k)%DIL, vv_js_d(k)%DIL)
82  END DO
83  !===End PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
84 
85  !===ASSEMBLING NS RESIDUAL MATRICES
86  ALLOCATE(les_mat(2),les_ksp(2))
87  ALLOCATE(zero_dir1(SIZE(vv_js_d(1)%DIL)))
88  ALLOCATE(zero_dir2(SIZE(vv_js_d(2)%DIL)))
89  ALLOCATE(zero_dir3(SIZE(vv_js_d(3)%DIL)))
90 
91  DO k = 1, 2
92  nu_mat = k
93  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la,les_mat(nu_mat), clean=.false.)
94  CALL qs_mass_vect_3x3(vv_3_la, vv_mesh, 1.d0, les_mat(nu_mat))
95  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
96  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
97  les_mat(nu_mat), vv_3_la)
98  END IF
99  CALL init_solver(inputs%my_par_vv,les_ksp(nu_mat),les_mat(nu_mat),comm_one_d(1),&
100  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
101  END DO
102  iteration = 0
103  END IF !end of once
104 
105  iteration = iteration + 1
106 
107  !===Compute Strain rate tensor
108  CALL smb_explicit_strain_rate_tensor(vv_mesh, list_mode, un_m1, strain_rate_tensor)
109  strain_rate_tensor = 1/inputs%Re*strain_rate_tensor
110  CALL smb_explicit_strain_rate_tensor_bdy(vv_mesh, list_mode, un_m1, strain_rate_tensor_scal_n_bdy)
111  strain_rate_tensor_scal_n_bdy = -1/inputs%Re*strain_rate_tensor_scal_n_bdy
112 
113  !===Computation of rhs at Gauss points for every mode without the strain rate tensor
114  IF (PRESENT(opt_tempn)) THEN
115  CALL rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time-inputs%dt, &
116  (un-un_m2)/(2*inputs%dt), pn_m1, rotv_v_m1, rhs_gauss, opt_tempn=opt_tempn)
117  ELSE
118  CALL rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time-inputs%dt, &
119  (un-un_m2)/(2*inputs%dt), pn_m1, rotv_v_m1, rhs_gauss)
120  END IF
121  !===End Computation of rhs
122 
123  DO i = 1, SIZE(list_mode)
124  mode = list_mode(i)
125 
126  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
127  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236,&
128  opt_tensor=strain_rate_tensor(:,:,:,i), opt_tensor_scaln_bdy=strain_rate_tensor_scal_n_bdy(:,:,i))
129 
130  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
131  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
132  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
133  END IF
134  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
135 
136 
137  !===Solve linear system for momentum equation
138  !Solve system 1, ur_c, ut_s, uz_c
139  nu_mat = 1
140  CALL solver(les_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
141  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
142  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
143  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,1,i))
144  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,4,i))
145  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,5,i))
146 
147  !Solve system 2, ur_s, ut_c, uz_s
148  nu_mat = 2
149  CALL solver(les_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
150  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
151  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
152  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,2,i))
153  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,3,i))
154  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,6,i))
155  !===End Solve linear system for momentum equation
156  END DO
157 
158  !norm_res_ns_L2 = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, res_ns)
159  !CALL vtu_3d(res_ns, 'vv_mesh', 'Resns', 'resns', 'new')
160  !WRITE(*,*) 'norm_res_ns_L2 = ', norm_res_ns_L2
161 !!$ IF (MOD(iteration,inputs%freq_en) == 0) THEN
162 !!$ CALL MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
163 !!$ norm_res_ns_L2 = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, res_ns)
164 !!$ !CALL vtu_3d(res_ns, 'vv_mesh', 'Resns', 'resns', 'new')
165 !!$ IF (rank==0) THEN
166 !!$ WRITE(*,*) 'norm_res_ns_L2 = ', norm_res_ns_L2
167 !!$ END IF
168 !!$ END IF
169 
170  !===Compute un_m1 and res_ns on gauss points
171  DO i = 1, SIZE(list_mode)
172  index = 0
173  DO m = 1, vv_mesh%dom_me
174  j_loc = vv_mesh%jj(:,m)
175  DO l = 1, vv_mesh%gauss%l_G
176  index = index + 1
177  DO TYPE = 1, 6
178  un_m1_gauss(index,TYPE,i) = sum(un_m1(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
179  res_ns_gauss(index,TYPE,i) = sum(res_ns(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
180  END DO
181  END DO
182  END DO
183  END DO
184  !===End compute un_m1 and res_ns on gauss points
185 
186  !===Compute entropy viscosity times gradient of 2*un-un_m1 in real space by FFT
187  CALL smb_explicit_grad_vel_les(vv_mesh, list_mode, 2*un-un_m1, strain_rate_tensor)
188 
189  norm_vel_l2 = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_m1)
190 
191  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
192  bloc_size = SIZE(un_m1_gauss,1)/nb_procs+1
193  bloc_size = vv_mesh%gauss%l_G*(bloc_size/vv_mesh%gauss%l_G)+vv_mesh%gauss%l_G
194  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
195  CALL fft_compute_entropy_visc(comm_one_d(2), comm_one_d(1), un_m1_gauss, res_ns_gauss, &
196  strain_rate_tensor(1,:,:,:), strain_rate_tensor(2,:,:,:), strain_rate_tensor(3,:,:,:), &
197  vv_mesh%hloc_gauss, visco_entro_grad_u(1,:,:,:), visco_entro_grad_u(2,:,:,:), &
198  visco_entro_grad_u(3,:,:,:), nb_procs, bloc_size, m_max_pad, norm_vel_l2**2/volume_3d, vv_mesh%gauss%l_G)
199  !===End compute entropy viscosity times gradient of 2*un-un_m1 in real space by FFT
200 
201  END SUBROUTINE compute_entropy_viscosity
202 
203  SUBROUTINE compute_entropy_viscosity_mom(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, &
204  momentum, momentum_m1, momentum_m2, pn_m1, un_m1, tensor_m1, visc_grad_vel_m1, tensor_surface_gauss, &
205  rotb_b_m1, visco_dyn_m1, density_m2, density_m1, density, tempn, visc_entro_real, visc_entro_level_real)
206  USE def_type_mesh
207  USE fem_m_axi
208  USE solve_petsc
209  USE periodic
210  USE rhs_gauss_computing
211  USE rhs_para_assembling
212  USE dir_nodes_petsc
213  USE st_matrix
214  USE sft_parallele
215  USE tn_axi
216  USE input_data
217  USE my_util
218  USE subroutine_mass
219 #include "petsc/finclude/petsc.h"
220  USE petsc
221  IMPLICIT NONE
222  TYPE(petsc_csr_la) :: vv_3_LA
223  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
224  REAL(KIND=8), INTENT(IN) :: time
225  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
226  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: momentum, momentum_m1, momentum_m2
227  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: un_m1
228  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn_m1
229  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: tensor_m1, visc_grad_vel_m1
230  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: tensor_surface_gauss !t=tn_m1 in bdf1 / tn if bdf2)
231  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotb_b_m1
232  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visco_dyn_m1
233  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: density_m2, density_m1, density
234  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tempn
235  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: visc_entro_real
236  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: visc_entro_level_real
237  !TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_D
238  LOGICAL, SAVE :: once = .true.
239  REAL(KIND=8), SAVE :: Volume_3D
240  INTEGER, SAVE :: iteration
241  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: tensor_m1_gauss
242  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_Gs*vv_mesh%dom_mes,6,SIZE(list_mode)) :: visc_grad_vel_m1_scal_n_bdy
243  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_Gs*vv_mesh%dom_mes,6,SIZE(list_mode)) :: tensor_m1_scal_n_bdy
244  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
245  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: un_m1_gauss
246  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: momentum_m1_gauss
247  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: res_ns_gauss
248  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: res_ns
249  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,2,SIZE(list_mode)) :: res_mass_gauss
250  REAL(KIND=8), DIMENSION(pp_mesh%np,6,SIZE(list_mode)) :: un_m1_P1, momentum_m1_P1, res_ns_P1
251  REAL(KIND=8), DIMENSION(pp_mesh%np,2,SIZE(list_mode)) :: density_P1, density_m2_P1
252  REAL(KIND=8), DIMENSION(pp_mesh%gauss%l_G*pp_mesh%dom_me,6,SIZE(list_mode)) :: un_m1_P1_gauss
253  REAL(KIND=8), DIMENSION(pp_mesh%gauss%l_G*pp_mesh%dom_me,6,SIZE(list_mode)) :: momentum_m1_P1_gauss
254  REAL(KIND=8), DIMENSION(pp_mesh%gauss%l_G*pp_mesh%dom_me,6,SIZE(list_mode)) :: res_ns_P1_gauss
255  REAL(KIND=8), DIMENSION(pp_mesh%gauss%l_G*pp_mesh%dom_me,2,SIZE(list_mode)) :: res_mass_P1_gauss
256  INTEGER, POINTER, DIMENSION(:) :: vv_3_ifrom
257  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
258  INTEGER, DIMENSION(pp_mesh%gauss%n_w) :: j_loc_P1
259  REAL(KIND=8) :: norm_vel_L2, norm_mom_L2
260  !REAL(KIND=8) :: norm_res_ns_L2
261  !INTEGER :: rank
262  INTEGER :: i, k, nu_mat, mode, m, l, TYPE, index, n
263  INTEGER :: nb_procs, code, bloc_size, m_max_pad
264  petscerrorcode :: ierr
265  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
266  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
267  mat, DIMENSION(:), POINTER, SAVE :: les_mat
268  ksp, DIMENSION(:), POINTER, SAVE :: les_ksp
269 
270  IF (once) THEN
271 
272  once = .false.
273 
274  !===CREATE PETSC VECTORS AND GHOSTS
275  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
276  n = 3*vv_mesh%dom_np
277  CALL veccreateghost(comm_one_d(1), n, &
278  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
279  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
280  CALL vecduplicate(vx_3, vb_3_145, ierr)
281  CALL vecduplicate(vx_3, vb_3_236, ierr)
282 
283  !===Compute Volume
284  CALL twod_volume(comm_one_d(1),vv_mesh,volume_3d)
285  volume_3d = volume_3d*2*acos(-1.d0)
286 
287  !===ASSEMBLING NS RESIDUAL MATRICES
288  ALLOCATE(les_mat(2),les_ksp(2))
289  DO k = 1, 2
290  nu_mat = k
291  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la,les_mat(nu_mat), clean=.false.)
292  CALL qs_mass_vect_3x3(vv_3_la, vv_mesh, 1.d0, les_mat(nu_mat))
293 !!$ IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
294 !!$ CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
295 !!$ LES_mat(nu_mat), vv_3_LA)
296 !!$ END IF
297  CALL init_solver(inputs%my_par_vv,les_ksp(nu_mat),les_mat(nu_mat),comm_one_d(1),&
298  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
299  END DO
300 
301  iteration = 0
302  END IF !end of once
303 
304  iteration = iteration + 1
305 
306  !===Compute Strain rate tensor scalar normal on bdy
307  CALL smb_explicit_strain_rate_tensor_bdy_mom(comm_one_d(2), vv_mesh, list_mode, visco_dyn_m1, un_m1, &
308  visc_grad_vel_m1_scal_n_bdy)
309  visc_grad_vel_m1_scal_n_bdy = -1/inputs%Re*visc_grad_vel_m1_scal_n_bdy
310 
311  !===Compute tensor scalar normal on bdy
312  CALL smb_explicit_tensor_bdy(vv_mesh, list_mode, tensor_m1, tensor_m1_scal_n_bdy)
313 
314  !===Computation of rhs at Gauss points for every mode without tensors
315  CALL rhs_residual_ns_gauss_3x3_mom(vv_mesh, pp_mesh, list_mode, time-inputs%dt, &
316  (momentum-momentum_m2)/(2*inputs%dt), pn_m1, density_m1, rotb_b_m1, rhs_gauss, opt_tempn=tempn)
317  !===End Computation of rhs
318 
319  !===Computation of tensor_m1(=m x u) on gauss points
320  DO i = 1, SIZE(list_mode)
321  index = 0
322  DO m = 1, vv_mesh%dom_me
323  j_loc = vv_mesh%jj(:,m)
324  DO l = 1, vv_mesh%gauss%l_G
325  index = index + 1
326  DO TYPE = 1, 6
327  DO k = 1, 3
328  tensor_m1_gauss(k,index,TYPE,i)=sum(tensor_m1(k,j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
329  END DO
330  END DO
331  END DO
332  END DO
333  END DO
334  !===End computation of tensor_m1(=m x u) on gauss points
335 
336  DO i = 1, SIZE(list_mode)
337  mode = list_mode(i)
338 
339  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
340  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236,&
341  opt_tensor=-tensor_m1_gauss(:,:,:,i)+visc_grad_vel_m1(:,:,:,i)+tensor_surface_gauss(:,:,:,i), &
342  opt_tensor_scaln_bdy=visc_grad_vel_m1_scal_n_bdy(:,:,i)+tensor_m1_scal_n_bdy(:,:,i))
343 
344 !!$ IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
345 !!$ CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_LA)
346 !!$ CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_LA)
347 !!$ END IF
348  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
349 
350 
351  !===Solve linear system for momentum equation
352  !Solve system 1, ur_c, ut_s, uz_c
353  nu_mat = 1
354  CALL solver(les_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
355  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
356  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
357  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,1,i))
358  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,4,i))
359  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,5,i))
360 
361  !Solve system 2, ur_s, ut_c, uz_s
362  nu_mat = 2
363  CALL solver(les_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
364  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
365  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
366  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,2,i))
367  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,3,i))
368  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,6,i))
369  !===End Solve linear system for momentum equation
370  END DO
371 
372 !!$ IF (MOD(iteration,inputs%freq_en) == 0) THEN
373 !!$ CALL MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
374 !!$ norm_res_ns_L2 = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, res_ns)
375 !!$ !CALL vtu_3d(res_ns, 'vv_mesh', 'Resns', 'resns', 'new')
376 !!$ IF (rank==0) THEN
377 !!$ WRITE(*,*) 'norm_res_ns_L2 = ', norm_res_ns_L2
378 !!$ END IF
379 !!$ END IF
380 
381  !===Compute un_m1 and res_ns on gauss points
382  DO i = 1, SIZE(list_mode)
383  index = 0
384  DO m = 1, vv_mesh%dom_me
385  j_loc = vv_mesh%jj(:,m)
386  DO l = 1, vv_mesh%gauss%l_G
387  index = index + 1
388  DO TYPE = 1, 6
389  un_m1_gauss(index,TYPE,i) = sum(un_m1(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
390  momentum_m1_gauss(index,TYPE,i) = sum(momentum_m1(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
391  res_ns_gauss(index,TYPE,i) = sum(res_ns(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
392  END DO
393  END DO
394  END DO
395  END DO
396  !===End compute un_m1 and res_ns on gauss points
397 
398  !===Compute res_mass on gauss points
399  CALL compute_res_mass_gauss(vv_mesh, list_mode, density_m2, density, momentum_m1, res_mass_gauss)
400  !===End compute res_mass on gauss points
401 
402  !===Compute entropy viscosity for momentum and level set equations on gauss points
403  norm_vel_l2 = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_m1)
404  norm_mom_l2 = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, momentum_m1)
405 
406  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
407  bloc_size = SIZE(un_m1_gauss,1)/nb_procs+1
408  bloc_size = vv_mesh%gauss%l_G*(bloc_size/vv_mesh%gauss%l_G)+vv_mesh%gauss%l_G
409  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
410 
411  IF (inputs%if_level_set_P2) THEN
412  !===Entropy viscosity for momentum
413  CALL fft_compute_entropy_visc_mom(comm_one_d(2), comm_one_d(1), un_m1_gauss, momentum_m1_gauss, res_ns_gauss, &
414  res_mass_gauss, vv_mesh%hloc_gauss, visc_entro_real, nb_procs, bloc_size, m_max_pad,&
415  vv_mesh%gauss%l_G, opt_c2_real_out=visc_entro_level_real)
416  ELSE
417  !===Entropy viscosity for momentum
418  CALL fft_compute_entropy_visc_mom(comm_one_d(2), comm_one_d(1), un_m1_gauss, momentum_m1_gauss, res_ns_gauss, &
419  res_mass_gauss, vv_mesh%hloc_gauss, visc_entro_real, nb_procs, bloc_size, m_max_pad,&
420  vv_mesh%gauss%l_G)
421 
422  !===Compute un_m1, momentum_m1, res_ns, density and density_m2 on P1 nodes
423  DO i = 1, SIZE(list_mode)
424  DO k = 1, 6
425  CALL project_p2_p1(vv_mesh%jj, pp_mesh%jj, un_m1(:,k,i), un_m1_p1(:,k,i))
426  CALL project_p2_p1(vv_mesh%jj, pp_mesh%jj, momentum_m1(:,k,i), momentum_m1_p1(:,k,i))
427  CALL project_p2_p1(vv_mesh%jj, pp_mesh%jj, res_ns(:,k,i), res_ns_p1(:,k,i))
428  END DO
429  DO k = 1, 2
430  CALL project_p2_p1(vv_mesh%jj, pp_mesh%jj, density(:,k,i), density_p1(:,k,i))
431  CALL project_p2_p1(vv_mesh%jj, pp_mesh%jj, density_m2(:,k,i), density_m2_p1(:,k,i))
432  END DO
433  END DO
434  !===End compute un_m1, momentum_m1, res_ns, density and density_m2 on P1 nodes
435 
436  !===Compute un_m1 and res_ns on P1 gauss points
437  DO i = 1, SIZE(list_mode)
438  index = 0
439  DO m = 1, pp_mesh%dom_me
440  j_loc_p1 = pp_mesh%jj(:,m)
441  DO l = 1, pp_mesh%gauss%l_G
442  index = index + 1
443  DO TYPE = 1, 6
444  un_m1_p1_gauss(index,TYPE,i) = sum(un_m1_p1(j_loc_p1,TYPE,i)*pp_mesh%gauss%ww(:,l))
445  momentum_m1_p1_gauss(index,TYPE,i) = sum(momentum_m1_p1(j_loc_p1,TYPE,i)*pp_mesh%gauss%ww(:,l))
446  res_ns_p1_gauss(index,TYPE,i) = sum(res_ns_p1(j_loc_p1,TYPE,i)*pp_mesh%gauss%ww(:,l))
447  END DO
448  END DO
449  END DO
450  END DO
451  !===End compute un_m1 and res_ns on P1 gauss points
452 
453  !===Compute res_mass on P1 gauss points
454  CALL compute_res_mass_gauss(pp_mesh, list_mode, density_m2_p1, density_p1, momentum_m1_p1, res_mass_p1_gauss)
455  !===End compute res_mass on P1 gauss points
456 
457  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
458  bloc_size = SIZE(un_m1_p1_gauss,1)/nb_procs+1
459  bloc_size = pp_mesh%gauss%l_G*(bloc_size/pp_mesh%gauss%l_G)+pp_mesh%gauss%l_G
460  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
461 
462  !===Entropy viscosity for level set
463  CALL fft_compute_entropy_visc_mom(comm_one_d(2), comm_one_d(1), un_m1_p1_gauss, momentum_m1_p1_gauss,&
464  res_ns_p1_gauss, res_mass_p1_gauss, pp_mesh%hloc_gauss, visc_entro_level_real, nb_procs, bloc_size, m_max_pad,&
465  pp_mesh%gauss%l_G)
466  END IF
467 
468 !!$ IF (MOD(iteration,inputs%freq_en) == 0) THEN
469 !!$ CALL MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
470 !!$ norm_res_ns_L2 = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, visc_entro)
471 !!$ IF (rank==0) THEN
472 !!$ WRITE(*,*) 'norm_visc_entro_L2 = ', norm_res_ns_L2
473 !!$ END IF
474 !!$ norm_res_ns_L2 = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, visc_entro_level)
475 !!$ IF (rank==0) THEN
476 !!$ WRITE(*,*) 'norm_visc_entro_level_L2 = ', norm_res_ns_L2
477 !!$ END IF
478 !!$ END IF
479 
480  END SUBROUTINE compute_entropy_viscosity_mom
481 
482  SUBROUTINE compute_entropy_viscosity_mom_no_level_set(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, &
483  momentum, momentum_m1, momentum_m2, pn_m1, un_m1, tensor_m1, &
484  rotb_b_m1, tempn, visc_entro_real)
485  USE def_type_mesh
486  USE fem_m_axi
487  USE solve_petsc
488  USE periodic
489  USE rhs_gauss_computing
490  USE rhs_para_assembling
491  USE dir_nodes_petsc
492  USE st_matrix
493  USE sft_parallele
494  USE tn_axi
495  USE input_data
496  USE my_util
497  USE subroutine_mass
498 #include "petsc/finclude/petsc.h"
499  USE petsc
500  IMPLICIT NONE
501  TYPE(petsc_csr_la) :: vv_3_LA
502  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
503  REAL(KIND=8), INTENT(IN) :: time
504  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
505  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: momentum, momentum_m1, momentum_m2
506  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: un_m1
507  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn_m1
508  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: tensor_m1
509  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotb_b_m1
510  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tempn
511  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: visc_entro_real
512  !TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_D
513  LOGICAL, SAVE :: once = .true.
514  REAL(KIND=8), SAVE :: Volume_3D
515  INTEGER, SAVE :: iteration
516  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: tensor_m1_gauss
517  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_Gs*vv_mesh%dom_mes,6,SIZE(list_mode)) :: tensor_m1_scal_n_bdy
518  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: visc_grad_vel_m1
519  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_Gs*vv_mesh%dom_mes,6,SIZE(list_mode)) :: visc_grad_vel_m1_scal_n_bdy
520  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
521  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,2,SIZE(list_mode)) :: res_mass_gauss
522  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: un_m1_gauss
523  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: momentum_m1_gauss
524  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: res_ns_gauss
525  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: res_ns
526  INTEGER, POINTER, DIMENSION(:) :: vv_3_ifrom
527  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
528  REAL(KIND=8) :: norm_vel_L2, norm_mom_L2
529  !REAL(KIND=8) :: norm_res_ns_L2
530  !INTEGER :: rank
531  INTEGER :: i, k, nu_mat, mode, m, l, TYPE, index, n
532  INTEGER :: nb_procs, code, bloc_size, m_max_pad
533  petscerrorcode :: ierr
534  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
535  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
536  mat, DIMENSION(:), POINTER, SAVE :: les_mat
537  ksp, DIMENSION(:), POINTER, SAVE :: les_ksp
538 
539  IF (once) THEN
540 
541  once = .false.
542 
543  !===CREATE PETSC VECTORS AND GHOSTS
544  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
545  n = 3*vv_mesh%dom_np
546  CALL veccreateghost(comm_one_d(1), n, &
547  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
548  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
549  CALL vecduplicate(vx_3, vb_3_145, ierr)
550  CALL vecduplicate(vx_3, vb_3_236, ierr)
551 
552  !===Compute Volume
553  CALL twod_volume(comm_one_d(1),vv_mesh,volume_3d)
554  volume_3d = volume_3d*2*acos(-1.d0)
555 
556  !===ASSEMBLING NS RESIDUAL MATRICES
557  ALLOCATE(les_mat(2),les_ksp(2))
558  DO k = 1, 2
559  nu_mat = k
560  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la,les_mat(nu_mat), clean=.false.)
561  CALL qs_mass_vect_3x3(vv_3_la, vv_mesh, 1.d0, les_mat(nu_mat))
562  CALL init_solver(inputs%my_par_vv,les_ksp(nu_mat),les_mat(nu_mat),comm_one_d(1),&
563  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
564  END DO
565 
566  iteration = 0
567  END IF !end of once
568 
569  iteration = iteration + 1
570 
571  !===Compute Strain rate tensor
572  CALL smb_explicit_strain_rate_tensor(vv_mesh, list_mode, un_m1, visc_grad_vel_m1)
573  visc_grad_vel_m1 = 1/inputs%Re* visc_grad_vel_m1
574  CALL smb_explicit_strain_rate_tensor_bdy(vv_mesh, list_mode, un_m1, visc_grad_vel_m1_scal_n_bdy)
575  visc_grad_vel_m1_scal_n_bdy = -1/inputs%Re*visc_grad_vel_m1_scal_n_bdy
576  !===End Compute Strain rate tensor
577 
578  !===Compute tensor scalar normal on bdy
579  CALL smb_explicit_tensor_bdy(vv_mesh, list_mode, tensor_m1, tensor_m1_scal_n_bdy)
580 
581  !===Computation of rhs at Gauss points for every mode without tensors
582  CALL rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time-inputs%dt, &
583  (momentum-momentum_m2)/(2*inputs%dt), pn_m1, -rotb_b_m1, rhs_gauss, opt_tempn=tempn)
584  !===End Computation of rhs
585 
586  !===Computation of tensor_m1(=m x u) on gauss points
587  DO i = 1, SIZE(list_mode)
588  index = 0
589  DO m = 1, vv_mesh%dom_me
590  j_loc = vv_mesh%jj(:,m)
591  DO l = 1, vv_mesh%gauss%l_G
592  index = index + 1
593  DO TYPE = 1, 6
594  DO k = 1, 3
595  tensor_m1_gauss(k,index,TYPE,i)=sum(tensor_m1(k,j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
596  END DO
597  END DO
598  END DO
599  END DO
600  END DO
601  !===End computation of tensor_m1(=m x u) on gauss points
602 
603  DO i = 1, SIZE(list_mode)
604  mode = list_mode(i)
605 
606  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
607  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236,&
608  opt_tensor=-tensor_m1_gauss(:,:,:,i)+visc_grad_vel_m1(:,:,:,i), &
609  opt_tensor_scaln_bdy=visc_grad_vel_m1_scal_n_bdy(:,:,i)+tensor_m1_scal_n_bdy(:,:,i))
610  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
611 
612 
613  !===Solve linear system for momentum equation
614  !Solve system 1, ur_c, ut_s, uz_c
615  nu_mat = 1
616  CALL solver(les_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
617  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
618  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
619  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,1,i))
620  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,4,i))
621  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,5,i))
622 
623  !Solve system 2, ur_s, ut_c, uz_s
624  nu_mat = 2
625  CALL solver(les_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
626  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
627  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
628  CALL extract(vx_3_ghost,1,1,vv_3_la,res_ns(:,2,i))
629  CALL extract(vx_3_ghost,2,2,vv_3_la,res_ns(:,3,i))
630  CALL extract(vx_3_ghost,3,3,vv_3_la,res_ns(:,6,i))
631  !===End Solve linear system for momentum equation
632  END DO
633 
634  !===Compute un_m1 and res_ns on gauss points
635  DO i = 1, SIZE(list_mode)
636  index = 0
637  DO m = 1, vv_mesh%dom_me
638  j_loc = vv_mesh%jj(:,m)
639  DO l = 1, vv_mesh%gauss%l_G
640  index = index + 1
641  DO TYPE = 1, 6
642  un_m1_gauss(index,TYPE,i) = sum(un_m1(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
643  momentum_m1_gauss(index,TYPE,i) = sum(momentum_m1(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
644  res_ns_gauss(index,TYPE,i) = sum(res_ns(j_loc,TYPE,i)*vv_mesh%gauss%ww(:,l))
645  END DO
646  END DO
647  END DO
648  END DO
649  !===End compute un_m1 and res_ns on gauss points
650 
651  !===Compute entropy viscosity for momentum and level set equations on gauss points
652  norm_vel_l2 = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_m1)
653  norm_mom_l2 = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, momentum_m1)
654 
655  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
656  bloc_size = SIZE(un_m1_gauss,1)/nb_procs+1
657  bloc_size = vv_mesh%gauss%l_G*(bloc_size/vv_mesh%gauss%l_G)+vv_mesh%gauss%l_G
658  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
659 
660  !===Entropy viscosity for momentum
661  res_mass_gauss=0.d0
662  CALL fft_compute_entropy_visc_mom(comm_one_d(2), comm_one_d(1), un_m1_gauss, momentum_m1_gauss, res_ns_gauss, &
663  res_mass_gauss, vv_mesh%hloc_gauss, visc_entro_real, nb_procs, bloc_size, m_max_pad,&
664  vv_mesh%gauss%l_G)
665 
666  END SUBROUTINE compute_entropy_viscosity_mom_no_level_set
667 
668  SUBROUTINE smb_explicit_grad_vel_les(mesh, list_mode, vel, V_out)
669  USE gauss_points
670  USE sft_parallele
671  USE chaine_caractere
672  USE boundary
673  USE input_data
674  IMPLICIT NONE
675  TYPE(mesh_type), TARGET :: mesh
676  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
677  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vel
678  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: V_out
679  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: grad1_vel, grad2_vel, grad3_vel
680  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
681  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
682  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vel_loc
683  REAL(KIND=8) :: ray, hh, hm
684  INTEGER :: m, l , i, mode, index, k
685 
686  DO i = 1, SIZE(list_mode)
687  mode = list_mode(i)
688  index = 0
689  DO m = 1, mesh%dom_me
690  j_loc = mesh%jj(:,m)
691  DO k = 1, 6
692  vel_loc(:,k) = vel(j_loc,k,i)
693  END DO
694  DO l = 1, mesh%gauss%l_G
695  index = index + 1
696  dw_loc = mesh%gauss%dw(:,:,l,m)
697 
698  !===Compute local mesh sizes
699  hh=mesh%hloc_gauss(index)
700  hm=min(mesh%hm(i),hh) !(JLG April 7 2017)
701  !hm=0.5d0/inputs%m_max
702 
703  !===Compute radius of Gauss point
704  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
705 
706  !-----------------Grad u_r on Gauss points------------------------------------
707  grad1_vel(index,1,i) = sum(vel_loc(:,1)*dw_loc(1,:))*hh
708  grad1_vel(index,2,i) = sum(vel_loc(:,2)*dw_loc(1,:))*hh
709  grad1_vel(index,3,i) = (mode*sum(vel_loc(:,2)*mesh%gauss%ww(:,l)) - &
710  sum(vel_loc(:,3)*mesh%gauss%ww(:,l)))/ray*hm
711  grad1_vel(index,4,i) = (-mode*sum(vel_loc(:,1)*mesh%gauss%ww(:,l)) - &
712  sum(vel_loc(:,4)*mesh%gauss%ww(:,l)))/ray*hm
713  grad1_vel(index,5,i) = sum(vel_loc(:,1)*dw_loc(2,:))*hh
714  grad1_vel(index,6,i) = sum(vel_loc(:,2)*dw_loc(2,:))*hh
715 
716  !-----------------Grad u_th on Gauss points-----------------------------------
717  grad2_vel(index,1,i) = sum(vel_loc(:,3)*dw_loc(1,:))*hh
718  grad2_vel(index,2,i) = sum(vel_loc(:,4)*dw_loc(1,:))*hh
719  grad2_vel(index,3,i) = (mode*sum(vel_loc(:,4)*mesh%gauss%ww(:,l)) + &
720  sum(vel_loc(:,1)*mesh%gauss%ww(:,l)))/ray*hm
721  grad2_vel(index,4,i) = (-mode*sum(vel_loc(:,3)*mesh%gauss%ww(:,l)) + &
722  sum(vel_loc(:,2)*mesh%gauss%ww(:,l)))/ray*hm
723  grad2_vel(index,5,i) = sum(vel_loc(:,3)*dw_loc(2,:))*hh
724  grad2_vel(index,6,i) = sum(vel_loc(:,4)*dw_loc(2,:))*hh
725 
726  !-----------------Grad u_z on Gauss points------------------------------------
727  grad3_vel(index,1,i) = sum(vel_loc(:,5)*dw_loc(1,:))*hh
728  grad3_vel(index,2,i) = sum(vel_loc(:,6)*dw_loc(1,:))*hh
729  grad3_vel(index,3,i) = mode*sum(vel_loc(:,6)*mesh%gauss%ww(:,l))/ray*hm
730  grad3_vel(index,4,i) = -mode*sum(vel_loc(:,5)*mesh%gauss%ww(:,l))/ray*hm
731  grad3_vel(index,5,i) = sum(vel_loc(:,5)*dw_loc(2,:))*hh
732  grad3_vel(index,6,i) = sum(vel_loc(:,6)*dw_loc(2,:))*hh
733  ENDDO
734  ENDDO
735  END DO
736 
737  v_out(1,:,:,:) = grad1_vel
738  v_out(2,:,:,:) = grad2_vel
739  v_out(3,:,:,:) = grad3_vel
740 
741  END SUBROUTINE smb_explicit_grad_vel_les
742 
743  SUBROUTINE smb_explicit_strain_rate_tensor(mesh, list_mode, vel, V_out)
744  USE gauss_points
745  USE sft_parallele
746  USE chaine_caractere
747  USE boundary
748  USE input_data
749  IMPLICIT NONE
750  TYPE(mesh_type), TARGET :: mesh
751  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
752  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vel
753  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: V_out
754  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: grad1_vel, grad2_vel, grad3_vel
755  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: gradT1_vel, gradT2_vel, gradT3_vel
756  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
757  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
758  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vel_loc
759  REAL(KIND=8) :: ray
760  INTEGER :: m, l , i, mode, index, k
761 
762  DO i = 1, SIZE(list_mode)
763  mode = list_mode(i)
764  index = 0
765  DO m = 1, mesh%dom_me
766  j_loc = mesh%jj(:,m)
767  DO k = 1, 6
768  vel_loc(:,k) = vel(j_loc,k,i)
769  END DO
770  DO l = 1, mesh%gauss%l_G
771  index = index + 1
772  dw_loc = mesh%gauss%dw(:,:,l,m)
773 
774  !===Compute radius of Gauss point
775  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
776 
777  !-----------------Grad u_r on Gauss points------------------------------------
778  grad1_vel(index,1,i) = sum(vel_loc(:,1)*dw_loc(1,:))
779  grad1_vel(index,2,i) = sum(vel_loc(:,2)*dw_loc(1,:))
780  grad1_vel(index,3,i) = (mode*sum(vel_loc(:,2)*mesh%gauss%ww(:,l)) - &
781  sum(vel_loc(:,3)*mesh%gauss%ww(:,l)))/ray
782  grad1_vel(index,4,i) = (-mode*sum(vel_loc(:,1)*mesh%gauss%ww(:,l)) - &
783  sum(vel_loc(:,4)*mesh%gauss%ww(:,l)))/ray
784  grad1_vel(index,5,i) = sum(vel_loc(:,1)*dw_loc(2,:))
785  grad1_vel(index,6,i) = sum(vel_loc(:,2)*dw_loc(2,:))
786 
787  !-----------------Grad u_th on Gauss points-----------------------------------
788  grad2_vel(index,1,i) = sum(vel_loc(:,3)*dw_loc(1,:))
789  grad2_vel(index,2,i) = sum(vel_loc(:,4)*dw_loc(1,:))
790  grad2_vel(index,3,i) = (mode*sum(vel_loc(:,4)*mesh%gauss%ww(:,l)) + &
791  sum(vel_loc(:,1)*mesh%gauss%ww(:,l)))/ray
792  grad2_vel(index,4,i) = (-mode*sum(vel_loc(:,3)*mesh%gauss%ww(:,l)) + &
793  sum(vel_loc(:,2)*mesh%gauss%ww(:,l)))/ray
794  grad2_vel(index,5,i) = sum(vel_loc(:,3)*dw_loc(2,:))
795  grad2_vel(index,6,i) = sum(vel_loc(:,4)*dw_loc(2,:))
796 
797  !-----------------Grad u_z on Gauss points------------------------------------
798  grad3_vel(index,1,i) = sum(vel_loc(:,5)*dw_loc(1,:))
799  grad3_vel(index,2,i) = sum(vel_loc(:,6)*dw_loc(1,:))
800  grad3_vel(index,3,i) = mode*sum(vel_loc(:,6)*mesh%gauss%ww(:,l))/ray
801  grad3_vel(index,4,i) = -mode*sum(vel_loc(:,5)*mesh%gauss%ww(:,l))/ray
802  grad3_vel(index,5,i) = sum(vel_loc(:,5)*dw_loc(2,:))
803  grad3_vel(index,6,i) = sum(vel_loc(:,6)*dw_loc(2,:))
804  ENDDO
805  ENDDO
806  END DO
807 
808  !-----------------GradT u_r on Gauss points------------------------------------
809  gradt1_vel(:,1,:) = grad1_vel(:,1,:)
810  gradt1_vel(:,2,:) = grad1_vel(:,2,:)
811  gradt1_vel(:,3,:) = grad2_vel(:,1,:)
812  gradt1_vel(:,4,:) = grad2_vel(:,2,:)
813  gradt1_vel(:,5,:) = grad3_vel(:,1,:)
814  gradt1_vel(:,6,:) = grad3_vel(:,2,:)
815  !-----------------GradT u_th on Gauss points-----------------------------------
816  gradt2_vel(:,1,:) = grad1_vel(:,3,:)
817  gradt2_vel(:,2,:) = grad1_vel(:,4,:)
818  gradt2_vel(:,3,:) = grad2_vel(:,3,:)
819  gradt2_vel(:,4,:) = grad2_vel(:,4,:)
820  gradt2_vel(:,5,:) = grad3_vel(:,3,:)
821  gradt2_vel(:,6,:) = grad3_vel(:,4,:)
822  !-----------------GradT u_z on Gauss points------------------------------------
823  gradt3_vel(:,1,:) = grad1_vel(:,5,:)
824  gradt3_vel(:,2,:) = grad1_vel(:,6,:)
825  gradt3_vel(:,3,:) = grad2_vel(:,5,:)
826  gradt3_vel(:,4,:) = grad2_vel(:,6,:)
827  gradt3_vel(:,5,:) = grad3_vel(:,5,:)
828  gradt3_vel(:,6,:) = grad3_vel(:,6,:)
829 
830  !===Grad = Grad + Grad^T
831  v_out(1,:,:,:) = grad1_vel + gradt1_vel
832  v_out(2,:,:,:) = grad2_vel + gradt2_vel
833  v_out(3,:,:,:) = grad3_vel + gradt3_vel
834 
835  END SUBROUTINE smb_explicit_strain_rate_tensor
836 
837  SUBROUTINE smb_explicit_strain_rate_tensor_bdy(mesh, list_mode, vel, V_out)
838  USE gauss_points
839  USE sft_parallele
840  USE chaine_caractere
841  USE boundary
842  USE input_data
843  IMPLICIT NONE
844  TYPE(mesh_type), TARGET :: mesh
845  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
846  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vel
847  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)), INTENT(OUT) :: V_out
848  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)) :: grad1_vel, grad2_vel, grad3_vel
849  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)) :: gradT1_vel, gradT2_vel, gradT3_vel
850  INTEGER, DIMENSION(mesh%gauss%n_ws) :: js_loc
851  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
852  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dws_loc
853  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vel_loc
854  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,6) :: vels_loc
855  REAL(KIND=8) :: ray
856  INTEGER :: m, ms, ls , i, mode, indexs, k
857 
858  DO i = 1, SIZE(list_mode)
859  mode = list_mode(i)
860  indexs = 0
861  DO ms = 1, mesh%dom_mes
862  js_loc = mesh%jjs(:,ms)
863  m = mesh%neighs(ms)
864  j_loc = mesh%jj(:,m)
865  DO k = 1, 6
866  vels_loc(:,k) = vel(js_loc,k,i)
867  vel_loc(:,k) = vel(j_loc,k,i)
868  END DO
869  DO ls = 1, mesh%gauss%l_Gs
870  indexs = indexs + 1
871  dws_loc = mesh%gauss%dw_s(:,:,ls,ms)
872 
873  !===Compute radius of Gauss point
874  ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
875 
876  !===Dont compute on r = 0
877  IF (ray.LT.1.d-10) THEN
878  v_out(indexs,:,i) = 0.d0
879  cycle
880  END IF
881 
882  !-----------------Grad u_r on Gauss points------------------------------------
883  grad1_vel(indexs,1,i) = sum(vel_loc(:,1)*dws_loc(1,:))
884  grad1_vel(indexs,2,i) = sum(vel_loc(:,2)*dws_loc(1,:))
885  grad1_vel(indexs,3,i) = (mode*sum(vels_loc(:,2)*mesh%gauss%wws(:,ls)) - &
886  sum(vels_loc(:,3)*mesh%gauss%wws(:,ls)))/ray
887  grad1_vel(indexs,4,i) = (-mode*sum(vels_loc(:,1)*mesh%gauss%wws(:,ls)) - &
888  sum(vels_loc(:,4)*mesh%gauss%wws(:,ls)))/ray
889  grad1_vel(indexs,5,i) = sum(vel_loc(:,1)*dws_loc(2,:))
890  grad1_vel(indexs,6,i) = sum(vel_loc(:,2)*dws_loc(2,:))
891 
892  !-----------------Grad u_th on Gauss points-----------------------------------
893  grad2_vel(indexs,1,i) = sum(vel_loc(:,3)*dws_loc(1,:))
894  grad2_vel(indexs,2,i) = sum(vel_loc(:,4)*dws_loc(1,:))
895  grad2_vel(indexs,3,i) = (mode*sum(vels_loc(:,4)*mesh%gauss%wws(:,ls)) + &
896  sum(vels_loc(:,1)*mesh%gauss%wws(:,ls)))/ray
897  grad2_vel(indexs,4,i) = (-mode*sum(vels_loc(:,3)*mesh%gauss%wws(:,ls)) + &
898  sum(vels_loc(:,2)*mesh%gauss%wws(:,ls)))/ray
899  grad2_vel(indexs,5,i) = sum(vel_loc(:,3)*dws_loc(2,:))
900  grad2_vel(indexs,6,i) = sum(vel_loc(:,4)*dws_loc(2,:))
901 
902  !-----------------Grad u_z on Gauss points------------------------------------
903  grad3_vel(indexs,1,i) = sum(vel_loc(:,5)*dws_loc(1,:))
904  grad3_vel(indexs,2,i) = sum(vel_loc(:,6)*dws_loc(1,:))
905  grad3_vel(indexs,3,i) = mode*sum(vels_loc(:,6)*mesh%gauss%wws(:,ls))/ray
906  grad3_vel(indexs,4,i) = -mode*sum(vels_loc(:,5)*mesh%gauss%wws(:,ls))/ray
907  grad3_vel(indexs,5,i) = sum(vel_loc(:,5)*dws_loc(2,:))
908  grad3_vel(indexs,6,i) = sum(vel_loc(:,6)*dws_loc(2,:))
909 
910  !-----------------GradT u_r on Gauss points------------------------------------
911  gradt1_vel(indexs,1,i) = grad1_vel(indexs,1,i)
912  gradt1_vel(indexs,2,i) = grad1_vel(indexs,2,i)
913  gradt1_vel(indexs,3,i) = grad2_vel(indexs,1,i)
914  gradt1_vel(indexs,4,i) = grad2_vel(indexs,2,i)
915  gradt1_vel(indexs,5,i) = grad3_vel(indexs,1,i)
916  gradt1_vel(indexs,6,i) = grad3_vel(indexs,2,i)
917  !-----------------GradT u_th on Gauss points-----------------------------------
918  gradt2_vel(indexs,1,i) = grad1_vel(indexs,3,i)
919  gradt2_vel(indexs,2,i) = grad1_vel(indexs,4,i)
920  gradt2_vel(indexs,3,i) = grad2_vel(indexs,3,i)
921  gradt2_vel(indexs,4,i) = grad2_vel(indexs,4,i)
922  gradt2_vel(indexs,5,i) = grad3_vel(indexs,3,i)
923  gradt2_vel(indexs,6,i) = grad3_vel(indexs,4,i)
924  !-----------------GradT u_z on Gauss points------------------------------------
925  gradt3_vel(indexs,1,i) = grad1_vel(indexs,5,i)
926  gradt3_vel(indexs,2,i) = grad1_vel(indexs,6,i)
927  gradt3_vel(indexs,3,i) = grad2_vel(indexs,5,i)
928  gradt3_vel(indexs,4,i) = grad2_vel(indexs,6,i)
929  gradt3_vel(indexs,5,i) = grad3_vel(indexs,5,i)
930  gradt3_vel(indexs,6,i) = grad3_vel(indexs,6,i)
931 
932  !-----------------Grad_sym_u scalar normal-------------------------------------
933  v_out(indexs,1,i) = (grad1_vel(indexs,1,i)+gradt1_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
934  + (grad1_vel(indexs,5,i)+gradt1_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
935  v_out(indexs,2,i) = (grad1_vel(indexs,2,i)+gradt1_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
936  + (grad1_vel(indexs,6,i)+gradt1_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
937  v_out(indexs,3,i) = (grad2_vel(indexs,1,i)+gradt2_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
938  + (grad2_vel(indexs,5,i)+gradt2_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
939  v_out(indexs,4,i) = (grad2_vel(indexs,2,i)+gradt2_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
940  + (grad2_vel(indexs,6,i)+gradt2_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
941  v_out(indexs,5,i) = (grad3_vel(indexs,1,i)+gradt3_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
942  + (grad3_vel(indexs,5,i)+gradt3_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
943  v_out(indexs,6,i) = (grad3_vel(indexs,2,i)+gradt3_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
944  + (grad3_vel(indexs,6,i)+gradt3_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
945  ENDDO !l_Gs
946  ENDDO !mes
947  END DO !i
948 
949  END SUBROUTINE smb_explicit_strain_rate_tensor_bdy
950 
951  SUBROUTINE smb_explicit_strain_rate_tensor_bdy_mom(communicator, mesh, list_mode, visc_dyn, vel, V_out)
952  USE gauss_points
953  USE sft_parallele
954  USE chaine_caractere
955  USE boundary
956  USE input_data
957 #include "petsc/finclude/petsc.h"
958  USE petsc
959  IMPLICIT NONE
960  TYPE(mesh_type), TARGET :: mesh
961  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
962  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visc_dyn
963  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vel
964  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)), INTENT(OUT) :: V_out
965  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)) :: grad1_vel, grad2_vel, grad3_vel
966  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)) :: gradT1_vel, gradT2_vel, gradT3_vel
967  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)) :: V_bdy
968  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,2,SIZE(list_mode)) :: visc_dyn_gauss
969  INTEGER, DIMENSION(mesh%gauss%n_ws) :: js_loc
970  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
971  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dws_loc
972  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vel_loc
973  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,6) :: vels_loc
974  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,2) :: visc_dyns_loc
975  REAL(KIND=8) :: ray
976  INTEGER :: m, ms, ls , i, mode, indexs, k
977  INTEGER :: code, bloc_size, m_max_pad, nb_procs
978  mpi_comm :: communicator
979 
980  DO i = 1, SIZE(list_mode)
981  mode = list_mode(i)
982  indexs = 0
983  DO ms = 1, mesh%dom_mes
984  js_loc = mesh%jjs(:,ms)
985  m = mesh%neighs(ms)
986  j_loc = mesh%jj(:,m)
987  DO k = 1, 6
988  vels_loc(:,k) = vel(js_loc,k,i)
989  vel_loc(:,k) = vel(j_loc,k,i)
990  END DO
991  DO k = 1, 2
992  visc_dyns_loc(:,k) = visc_dyn(js_loc,k,i)
993  END DO
994  DO ls = 1, mesh%gauss%l_Gs
995  indexs = indexs + 1
996  dws_loc = mesh%gauss%dw_s(:,:,ls,ms)
997 
998  !===Compute radius of Gauss point
999  ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
1000 
1001  !===Dont compute on r = 0
1002  IF (ray.LT.1.d-10) THEN
1003  visc_dyn_gauss(indexs,:,i)=0.d0
1004  v_bdy(indexs,:,i) =0.d0
1005  cycle
1006  END IF
1007 
1008  !-----------------Visc_dyn on Gauss points------------------------------------
1009  visc_dyn_gauss(indexs,1,i) = sum(visc_dyns_loc(:,1)*mesh%gauss%wws(:,ls))
1010  visc_dyn_gauss(indexs,2,i) = sum(visc_dyns_loc(:,2)*mesh%gauss%wws(:,ls))
1011 
1012  !-----------------Grad u_r on Gauss points------------------------------------
1013  grad1_vel(indexs,1,i) = sum(vel_loc(:,1)*dws_loc(1,:))
1014  grad1_vel(indexs,2,i) = sum(vel_loc(:,2)*dws_loc(1,:))
1015  grad1_vel(indexs,3,i) = (mode*sum(vels_loc(:,2)*mesh%gauss%wws(:,ls)) - &
1016  sum(vels_loc(:,3)*mesh%gauss%wws(:,ls)))/ray
1017  grad1_vel(indexs,4,i) = (-mode*sum(vels_loc(:,1)*mesh%gauss%wws(:,ls)) - &
1018  sum(vels_loc(:,4)*mesh%gauss%wws(:,ls)))/ray
1019  grad1_vel(indexs,5,i) = sum(vel_loc(:,1)*dws_loc(2,:))
1020  grad1_vel(indexs,6,i) = sum(vel_loc(:,2)*dws_loc(2,:))
1021 
1022  !-----------------Grad u_th on Gauss points-----------------------------------
1023  grad2_vel(indexs,1,i) = sum(vel_loc(:,3)*dws_loc(1,:))
1024  grad2_vel(indexs,2,i) = sum(vel_loc(:,4)*dws_loc(1,:))
1025  grad2_vel(indexs,3,i) = (mode*sum(vels_loc(:,4)*mesh%gauss%wws(:,ls)) + &
1026  sum(vels_loc(:,1)*mesh%gauss%wws(:,ls)))/ray
1027  grad2_vel(indexs,4,i) = (-mode*sum(vels_loc(:,3)*mesh%gauss%wws(:,ls)) + &
1028  sum(vels_loc(:,2)*mesh%gauss%wws(:,ls)))/ray
1029  grad2_vel(indexs,5,i) = sum(vel_loc(:,3)*dws_loc(2,:))
1030  grad2_vel(indexs,6,i) = sum(vel_loc(:,4)*dws_loc(2,:))
1031 
1032  !-----------------Grad u_z on Gauss points------------------------------------
1033  grad3_vel(indexs,1,i) = sum(vel_loc(:,5)*dws_loc(1,:))
1034  grad3_vel(indexs,2,i) = sum(vel_loc(:,6)*dws_loc(1,:))
1035  grad3_vel(indexs,3,i) = mode*sum(vels_loc(:,6)*mesh%gauss%wws(:,ls))/ray
1036  grad3_vel(indexs,4,i) = -mode*sum(vels_loc(:,5)*mesh%gauss%wws(:,ls))/ray
1037  grad3_vel(indexs,5,i) = sum(vel_loc(:,5)*dws_loc(2,:))
1038  grad3_vel(indexs,6,i) = sum(vel_loc(:,6)*dws_loc(2,:))
1039 
1040  !-----------------GradT u_r on Gauss points------------------------------------
1041  gradt1_vel(indexs,1,i) = grad1_vel(indexs,1,i)
1042  gradt1_vel(indexs,2,i) = grad1_vel(indexs,2,i)
1043  gradt1_vel(indexs,3,i) = grad2_vel(indexs,1,i)
1044  gradt1_vel(indexs,4,i) = grad2_vel(indexs,2,i)
1045  gradt1_vel(indexs,5,i) = grad3_vel(indexs,1,i)
1046  gradt1_vel(indexs,6,i) = grad3_vel(indexs,2,i)
1047  !-----------------GradT u_th on Gauss points-----------------------------------
1048  gradt2_vel(indexs,1,i) = grad1_vel(indexs,3,i)
1049  gradt2_vel(indexs,2,i) = grad1_vel(indexs,4,i)
1050  gradt2_vel(indexs,3,i) = grad2_vel(indexs,3,i)
1051  gradt2_vel(indexs,4,i) = grad2_vel(indexs,4,i)
1052  gradt2_vel(indexs,5,i) = grad3_vel(indexs,3,i)
1053  gradt2_vel(indexs,6,i) = grad3_vel(indexs,4,i)
1054  !-----------------GradT u_z on Gauss points------------------------------------
1055  gradt3_vel(indexs,1,i) = grad1_vel(indexs,5,i)
1056  gradt3_vel(indexs,2,i) = grad1_vel(indexs,6,i)
1057  gradt3_vel(indexs,3,i) = grad2_vel(indexs,5,i)
1058  gradt3_vel(indexs,4,i) = grad2_vel(indexs,6,i)
1059  gradt3_vel(indexs,5,i) = grad3_vel(indexs,5,i)
1060  gradt3_vel(indexs,6,i) = grad3_vel(indexs,6,i)
1061 
1062  !-----------------Grad_sym_u scalar normal-------------------------------------
1063  v_bdy(indexs,1,i) = (grad1_vel(indexs,1,i)+gradt1_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
1064  + (grad1_vel(indexs,5,i)+gradt1_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
1065  v_bdy(indexs,2,i) = (grad1_vel(indexs,2,i)+gradt1_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
1066  + (grad1_vel(indexs,6,i)+gradt1_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
1067  v_bdy(indexs,3,i) = (grad2_vel(indexs,1,i)+gradt2_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
1068  + (grad2_vel(indexs,5,i)+gradt2_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
1069  v_bdy(indexs,4,i) = (grad2_vel(indexs,2,i)+gradt2_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
1070  + (grad2_vel(indexs,6,i)+gradt2_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
1071  v_bdy(indexs,5,i) = (grad3_vel(indexs,1,i)+gradt3_vel(indexs,1,i))*mesh%gauss%rnorms(1,ls,ms) &
1072  + (grad3_vel(indexs,5,i)+gradt3_vel(indexs,5,i))*mesh%gauss%rnorms(2,ls,ms)
1073  v_bdy(indexs,6,i) = (grad3_vel(indexs,2,i)+gradt3_vel(indexs,2,i))*mesh%gauss%rnorms(1,ls,ms) &
1074  + (grad3_vel(indexs,6,i)+gradt3_vel(indexs,6,i))*mesh%gauss%rnorms(2,ls,ms)
1075  ENDDO !l_Gs
1076  ENDDO !mes
1077  END DO !i
1078 
1079  !===Compute visc_dyn*(Grad_sym_u scalar normal)
1080  CALL mpi_comm_size(communicator, nb_procs, code)
1081  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1082  bloc_size = SIZE(v_bdy,1)/nb_procs+1
1083  CALL fft_scalar_vect_dcl(communicator, v_bdy, visc_dyn_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
1084  !===End Compute visc_dyn*(Grad_sym_u scalar normal)
1085 
1086  END SUBROUTINE smb_explicit_strain_rate_tensor_bdy_mom
1087 
1088  SUBROUTINE smb_explicit_tensor_bdy(mesh, list_mode, tensor, V_out)
1089  USE gauss_points
1090  USE sft_parallele
1091  USE chaine_caractere
1092  USE boundary
1093  USE input_data
1094  IMPLICIT NONE
1095  TYPE(mesh_type), TARGET :: mesh
1096  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1097  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: tensor
1098  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*mesh%dom_mes,6,SIZE(list_mode)), INTENT(OUT) :: V_out
1099  INTEGER, DIMENSION(mesh%gauss%n_ws) :: js_loc
1100  REAL(KIND=8), DIMENSION(3,mesh%gauss%n_ws,6) :: tensors_loc
1101  REAL(KIND=8) :: ray
1102  INTEGER :: m, ms, ls , i, mode, indexs, k, n
1103 
1104  DO i = 1, SIZE(list_mode)
1105  mode = list_mode(i)
1106  indexs = 0
1107  DO ms = 1, mesh%dom_mes
1108  js_loc = mesh%jjs(:,ms)
1109  m = mesh%neighs(ms)
1110  DO k = 1, 6
1111  DO n = 1, 3
1112  tensors_loc(n,:,k) = tensor(n,js_loc,k,i)
1113  END DO
1114  END DO
1115  DO ls = 1, mesh%gauss%l_Gs
1116  indexs = indexs + 1
1117 
1118  !===Compute radius of Gauss point
1119  ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
1120 
1121  !===Dont compute on r = 0
1122  IF (ray.LT.1.d-10) THEN
1123  v_out(indexs,:,i) = 0.d0
1124  cycle
1125  END IF
1126  !-----------------Tensor scalar normal on gauss points------------------------
1127  v_out(indexs,1,i) = sum(tensors_loc(1,:,1)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1128  + sum(tensors_loc(1,:,5)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1129  v_out(indexs,2,i) = sum(tensors_loc(1,:,2)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1130  + sum(tensors_loc(1,:,6)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1131  v_out(indexs,3,i) = sum(tensors_loc(2,:,1)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1132  + sum(tensors_loc(2,:,5)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1133  v_out(indexs,4,i) = sum(tensors_loc(2,:,2)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1134  + sum(tensors_loc(2,:,6)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1135  v_out(indexs,5,i) = sum(tensors_loc(3,:,1)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1136  + sum(tensors_loc(3,:,5)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1137  v_out(indexs,6,i) = sum(tensors_loc(3,:,2)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(1,ls,ms) &
1138  + sum(tensors_loc(3,:,6)*mesh%gauss%wws(:,ls))*mesh%gauss%rnorms(2,ls,ms)
1139  ENDDO !l_Gs
1140  ENDDO !mes
1141  END DO !i
1142 
1143  END SUBROUTINE smb_explicit_tensor_bdy
1144 
1145  SUBROUTINE compute_res_mass_gauss(mesh, list_mode, density_m2, density, momentum_m1, c_out)
1146  USE gauss_points
1147  USE sft_parallele
1148  USE chaine_caractere
1149  USE boundary
1150  USE input_data
1151  IMPLICIT NONE
1152  TYPE(mesh_type), TARGET :: mesh
1153  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1154  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: density_m2, density
1155  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: momentum_m1
1156  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)), INTENT(OUT) :: c_out
1157  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: Div
1158  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: dens_m2_gauss, dens_gauss
1159  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1160  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1161  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: mom_loc
1162  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: dens_m2_loc, dens_loc
1163  REAL(KIND=8) :: ray
1164  INTEGER :: m, l , i, mode, index, k
1165 
1166  DO i = 1, SIZE(list_mode)
1167  mode = list_mode(i)
1168  index = 0
1169  DO m = 1, mesh%dom_me
1170  j_loc = mesh%jj(:,m)
1171  DO k = 1, 6
1172  mom_loc(:,k) = momentum_m1(j_loc,k,i)
1173  END DO
1174  DO k = 1, 2
1175  dens_m2_loc(:,k) = density_m2(j_loc,k,i)
1176  dens_loc(:,k) = density(j_loc,k,i)
1177  END DO
1178  DO l = 1, mesh%gauss%l_G
1179  index = index + 1
1180  dw_loc = mesh%gauss%dw(:,:,l,m)
1181 
1182  !===Compute radius of Gauss point
1183  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1184 
1185  !===Compute density_m2 and density on gauss point
1186  dens_m2_gauss(index,1,i) = sum(dens_m2_loc(:,1)*mesh%gauss%ww(:,l))
1187  dens_m2_gauss(index,2,i) = sum(dens_m2_loc(:,2)*mesh%gauss%ww(:,l))
1188 
1189  dens_gauss(index,1,i) = sum(dens_loc(:,1)*mesh%gauss%ww(:,l))
1190  dens_gauss(index,2,i) = sum(dens_loc(:,2)*mesh%gauss%ww(:,l))
1191 
1192  !===Compute divergence of momentum on gauss point
1193  div(index,1,i) = sum(mom_loc(:,1)*dw_loc(1,:)) + sum(mom_loc(:,1)*mesh%gauss%ww(:,l))/ray &
1194  + (mode/ray)*sum(mom_loc(:,4)*mesh%gauss%ww(:,l)) + sum(mom_loc(:,5)*dw_loc(2,:))
1195  div(index,2,i) = sum(mom_loc(:,2)*dw_loc(1,:)) + sum(mom_loc(:,2)*mesh%gauss%ww(:,l))/ray &
1196  - (mode/ray)*sum(mom_loc(:,3)*mesh%gauss%ww(:,l)) + sum(mom_loc(:,6)*dw_loc(2,:))
1197  END DO
1198  END DO
1199  END DO
1200 
1201  c_out = (dens_gauss-dens_m2_gauss)/(2*inputs%dt) + div
1202 
1203  END SUBROUTINE compute_res_mass_gauss
1204 
1205  SUBROUTINE twod_volume(communicator,mesh,RESLT)
1206  !===========================
1207  USE def_type_mesh
1208 #include "petsc/finclude/petsc.h"
1209  USE petsc
1210  IMPLICIT NONE
1211  TYPE(mesh_type) :: mesh
1212  REAL(KIND=8), INTENT(OUT) :: RESLT
1213  REAL(KIND=8) :: vol_loc, vol_out
1214  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1215  INTEGER :: m, l , i , ni, code
1216  REAL(KIND=8) :: ray
1217  mpi_comm :: communicator
1218  vol_loc = 0.d0
1219  DO m = 1, mesh%dom_me
1220  j_loc = mesh%jj(:,m)
1221  DO l = 1, mesh%gauss%l_G
1222  ray = 0
1223  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
1224  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1225  END DO
1226  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
1227  ENDDO
1228  ENDDO
1229  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
1230  reslt = vol_out
1231  END SUBROUTINE twod_volume
1232 
1233 END MODULE entropy_viscosity
entropy_viscosity::smb_explicit_tensor_bdy
subroutine smb_explicit_tensor_bdy(mesh, list_mode, tensor, V_out)
Definition: entropy_viscosity.f90:1089
entropy_viscosity::smb_explicit_grad_vel_les
subroutine smb_explicit_grad_vel_les(mesh, list_mode, vel, V_out)
Definition: entropy_viscosity.f90:669
solve_petsc::solver
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:99
sft_parallele::fft_compute_entropy_visc
subroutine, public fft_compute_entropy_visc(communicator, communicator_S, V1_in, V2_in, V3_in, V4_in, V5_in, hloc_gauss, V1_out, V2_out, V3_out, nb_procs, bloc_size, m_max_pad, residual_normalization, l_G, temps)
Definition: fft_parallel.f90:2255
entropy_viscosity::smb_explicit_strain_rate_tensor_bdy_mom
subroutine smb_explicit_strain_rate_tensor_bdy_mom(communicator, mesh, list_mode, visc_dyn, vel, V_out)
Definition: entropy_viscosity.f90:952
gauss_points
Definition: associate_gauss.f90:4
st_matrix::extract
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:34
entropy_viscosity::compute_res_mass_gauss
subroutine compute_res_mass_gauss(mesh, list_mode, density_m2, density, momentum_m1, c_out)
Definition: entropy_viscosity.f90:1146
entropy_viscosity::smb_explicit_strain_rate_tensor_bdy
subroutine smb_explicit_strain_rate_tensor_bdy(mesh, list_mode, vel, V_out)
Definition: entropy_viscosity.f90:838
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
entropy_viscosity::twod_volume
subroutine twod_volume(communicator, mesh, RESLT)
Definition: entropy_viscosity.f90:1206
subroutine_mass::project_p2_p1
subroutine, public project_p2_p1(jj_P2, jj_P1, pp_P2, pp_P1)
Definition: sub_mass.f90:291
entropy_viscosity::smb_explicit_strain_rate_tensor
subroutine smb_explicit_strain_rate_tensor(mesh, list_mode, vel, V_out)
Definition: entropy_viscosity.f90:744
entropy_viscosity::compute_entropy_viscosity_mom
subroutine, public compute_entropy_viscosity_mom(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, momentum, momentum_m1, momentum_m2, pn_m1, un_m1, tensor_m1, visc_grad_vel_m1, tensor_surface_gauss, rotb_b_m1, visco_dyn_m1, density_m2, density_m1, density, tempn, visc_entro_real, visc_entro_level_real)
Definition: entropy_viscosity.f90:206
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
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
rhs_para_assembling
Definition: assembling_rhs.f90:1
dir_nodes_petsc
Definition: dir_nodes_petsc.f90:4
fem_m_axi::qs_mass_vect_3x3
subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix)
Definition: fem_M_axi.f90:312
sft_parallele
Definition: fft_parallel.f90:4
entropy_viscosity
Definition: entropy_viscosity.f90:1
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_nodes_parallel
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
Definition: dir_nodes_petsc.f90:53
sft_parallele::fft_compute_entropy_visc_mom
subroutine, public fft_compute_entropy_visc_mom(communicator, communicator_S, V1_in, V2_in, V3_in, c1_in, hloc_gauss, c1_real_out, nb_procs, bloc_size, m_max_pad, l_G, opt_c2_real_out, temps)
Definition: fft_parallel.f90:2537
entropy_viscosity::compute_entropy_viscosity_mom_no_level_set
subroutine, public compute_entropy_viscosity_mom_no_level_set(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, momentum, momentum_m1, momentum_m2, pn_m1, un_m1, tensor_m1, rotb_b_m1, tempn, visc_entro_real)
Definition: entropy_viscosity.f90:485
periodic::periodic_matrix_petsc
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
Definition: prep_mesh_periodic.f90:378
tn_axi
Definition: tn_axi.f90:5
fourier_to_real_for_vtu
Definition: fourier_real.f90:1
rhs_gauss_computing
Definition: computing_rhs_gauss.f90:1
solve_petsc
Definition: solver.f90:1
rhs_gauss_computing::rhs_residual_ns_gauss_3x3_mom
subroutine, public rhs_residual_ns_gauss_3x3_mom(vv_mesh, pp_mesh, list_mode, time, du_dt, pn, density, rotb_b, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:260
tn_axi::norm_sf
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:40
rhs_gauss_computing::rhs_residual_ns_gauss_3x3
subroutine, public rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, du_dt, pn, rotv_v, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:155
sft_parallele::fft_scalar_vect_dcl
subroutine, public fft_scalar_vect_dcl(communicator, V1_in, V2_in, V_out, pb, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:857
rhs_para_assembling::rhs_3x3
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
Definition: assembling_rhs.f90:6
subroutine_mass
Definition: sub_mass.f90:1
entropy_viscosity::compute_entropy_viscosity
subroutine, public compute_entropy_viscosity(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, vvz_per, un, un_m1, un_m2, pn_m1, rotv_v_m1, visco_entro_grad_u, opt_tempn)
Definition: entropy_viscosity.f90:9
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
def_type_mesh
Definition: def_type_mesh.f90:4
gauss_points::ww
real(kind=8), dimension(:,:), pointer ww
Definition: associate_gauss.f90:10
periodic::periodic_rhs_petsc
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
Definition: prep_mesh_periodic.f90:459