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