SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
sub_ns_with_velocity.f90
Go to the documentation of this file.
1 !
2 !Authors Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyrights 2005
3 !Revised for PETSC, Jean-Luc Guermond, Franky Luddens, January 2011
4 !
5 MODULE subroutine_ns_with_u
6 USE my_util
7 PUBLIC :: moy, bdf2_ns_stress_bc_with_u
8 PRIVATE
9 CONTAINS
10 
11  SUBROUTINE bdf2_ns_stress_bc_with_u(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, list_mode, pp_mesh, &
12  vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, &
13  chmp_mag, bn_p2, opt_tempn)
14  !==============================
15  USE def_type_mesh
16  USE fem_m_axi
17  USE fem_rhs_axi
18  USE fem_tn_axi
19  USE dir_nodes_petsc
20  USE periodic
21  USE st_matrix
22  USE solve_petsc
23  USE dyn_line
24  USE boundary
25  USE chaine_caractere
26  USE sub_plot
27  USE st_matrix
28  USE input_data
29  USE rhs_gauss_computing
30  USE rhs_para_assembling
31  USE tn_axi
32  USE verbose
33  USE entropy_viscosity
34  USE my_util
35  IMPLICIT NONE
36  REAL(KIND=8) :: time, dt, re
37  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
38  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
39  TYPE(petsc_csr_la) :: vv_3_la, pp_1_la
40  TYPE(periodic_type), INTENT(IN) :: vvz_per, pp_per
41  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: incpn_m1, incpn
42  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: pn_m1, pn
43  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: un_m1, un
44  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), OPTIONAL:: opt_tempn
45  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: chmp_mag, bn_p2
46  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)):: chmp_mag_aux
47  !===Saved variables
48  INTEGER, SAVE :: m_max_c
49  !INTEGER, DIMENSION(:), POINTER, SAVE :: pp_js_D
50  !INTEGER, DIMENSION(:), POINTER, SAVE :: pp_js_axis_D
51  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: pp_global_d
52  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: pp_mode_global_js_d
53  TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_d
54  !INTEGER, DIMENSION(:), POINTER, SAVE :: vv_js_axis_D
55  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: vv_mode_global_js_d
56  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: vel_global_d
57  LOGICAL, SAVE :: once = .true.
58  INTEGER, SAVE :: my_petscworld_rank
59  REAL(KIND=8), DIMENSION(:,:,:,:),ALLOCATABLE, SAVE :: visco_entro_sym_grad_u
60  !===End saved variables
61 
62  !----------Declaration sans save-------------------------------------------------------
63  INTEGER, POINTER, DIMENSION(:) :: pp_1_ifrom, vv_3_ifrom
64 
65  INTEGER :: i, k, m, n, n1, n2, n3, n123
66  INTEGER :: nb_procs, code, nu_mat, mode
67  REAL(KIND=8) :: moyenne
68  !allocations des variables locales
69  REAL(KIND=8), DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: div
70  REAL(KIND=8), DIMENSION(pp_mesh%np, 2) :: pn_p1, phi
71  REAL(KIND=8), DIMENSION(vv_mesh%np, 6) :: un_p1
72  REAL(KIND=8), DIMENSION(vv_mesh%np, 6, SIZE(list_mode)) :: un_m2
73  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rotv_v, rotb_b, rotb_b_aux
74  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: kelvin_force
75  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
76  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: uext
77  REAL(KIND=8), DIMENSION(vv_mesh%np) :: vel_loc, vel_tot
78  REAL(KIND=8) :: tps, tps_tot, tps_cumul, coeff, vloc, cfl, cfl_max, norm
79  !April 17th 2008, JLG
80  REAL(KIND=8) :: one, zero, three
81  DATA zero, one, three/0.d0,1.d0,3.d0/
82 
83  !Communicators for Petsc, in space and Fourier------------------------------
84 #include "petsc/finclude/petsc.h"
85  petscerrorcode :: ierr
86  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
87  mat, DIMENSION(:), POINTER, SAVE :: vel_mat
88  mat, DIMENSION(:), POINTER, SAVE :: press_mat
89  mat, SAVE :: mass_mat, mass_mat0
90  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
91  vec, SAVE :: pb_1, pb_2, px_1, px_1_ghost
92  ksp, DIMENSION(:), POINTER, SAVE :: vel_ksp, press_ksp
93  ksp, SAVE :: mass_ksp, mass_ksp0
94  !------------------------------END OF DECLARATION--------------------------------------
95 
96  IF (once) THEN
97  once = .false.
98  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
99 
100  !===CREATE PETSC VECTORS AND GHOSTS
101  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
102  n = 3*vv_mesh%dom_np
103  CALL veccreateghost(comm_one_d(1), n, &
104  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
105  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
106  CALL vecduplicate(vx_3, vb_3_145, ierr)
107  CALL vecduplicate(vx_3, vb_3_236, ierr)
108 
109  CALL create_my_ghost(pp_mesh,pp_1_la,pp_1_ifrom)
110  n = pp_mesh%dom_np
111  CALL veccreateghost(comm_one_d(1), n, &
112  petsc_determine, SIZE(pp_1_ifrom), pp_1_ifrom, px_1, ierr)
113  CALL vecghostgetlocalform(px_1, px_1_ghost, ierr)
114  CALL vecduplicate(px_1, pb_1, ierr)
115  CALL vecduplicate(px_1, pb_2, ierr)
116  !===End CREATE PETSC VECTORS AND GHOSTS
117 
118  !===Number of modes on proc
119  m_max_c = SIZE(list_mode)
120  !===End umber of modes on proc
121 
122  !===PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
123  !===ATTENTION pressure BCs are no longer implemented
124  CALL scalar_glob_js_d(pp_mesh, list_mode, pp_1_la, pp_mode_global_js_d)
125  ALLOCATE(pp_global_d(m_max_c))
126  DO i = 1, m_max_c
127  ALLOCATE(pp_global_d(i)%DRL(SIZE(pp_mode_global_js_d(i)%DIL)))
128  END DO
129  !===End PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
130 
131  !===PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
132  CALL vector_glob_js_d(vv_mesh, list_mode, vv_3_la, inputs%vv_list_dirichlet_sides, &
133  vv_js_d, vv_mode_global_js_d)
134 
135  ALLOCATE(vel_global_d(m_max_c))
136  DO i = 1, m_max_c
137  ALLOCATE(vel_global_d(i)%DRL(SIZE(vv_mode_global_js_d(i)%DIL)))
138  END DO
139  !===END PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
140 
141  !===ASSEMBLE MASS MATRIX
142  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat, clean=.false.)
143  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat)
144  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
145  CALL periodic_matrix_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, mass_mat, pp_1_la)
146  END IF
147  DO i = 1, m_max_c
148  IF (list_mode(i)==0) cycle
149  CALL dirichlet_m_parallel(mass_mat,pp_mode_global_js_d(i)%DIL)
150  END DO
151  CALL init_solver(inputs%my_par_mass,mass_ksp,mass_mat,comm_one_d(1),&
152  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
153 
154  IF (minval(list_mode)==0) THEN
155  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat0, clean=.false.)
156  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat0)
157  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
158  CALL periodic_matrix_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, mass_mat0, pp_1_la)
159  END IF
160  DO i = 1, m_max_c
161  IF (list_mode(i).NE.0) cycle
162  CALL dirichlet_m_parallel(mass_mat0,pp_mode_global_js_d(i)%DIL)
163  CALL init_solver(inputs%my_par_mass,mass_ksp0,mass_mat0,comm_one_d(1),&
164  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
165  END DO
166  END IF
167  !===END ASSEMBLE MASS MATRIX
168 
169  !===ASSEMBLING VELOCITY MATRICES
170  ALLOCATE(visco_entro_sym_grad_u(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)))
171  visco_entro_sym_grad_u = 0.d0
172  ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
173  ALLOCATE(press_mat(m_max_c),press_ksp(m_max_c))
174  DO i = 1, m_max_c
175  mode = list_mode(i)
176  !===VELOCITY
177  nu_mat = 2*i-1
178  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
179  CALL qs_diff_mass_vect_3x3_divpenal(1, vv_3_la, vv_mesh, one/re, three/(2*dt), &
180  inputs%LES_coeff1, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
181  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
182  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
183  vel_mat(nu_mat), vv_3_la)
184  END IF
185  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
186  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
187  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
188  nu_mat = nu_mat+1
189  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
190  CALL qs_diff_mass_vect_3x3_divpenal(2, vv_3_la, vv_mesh, one/re, three/(2*dt), &
191  inputs%LES_coeff1, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
192  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
193  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
194  vel_mat(nu_mat), vv_3_la)
195  END IF
196  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
197  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
198  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
199  !===End VELOCITY
200 
201  !===PRESSURE
202  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, press_mat(i), clean=.false.)
203  !JLG Jan 2014 (regularize pressure matrix)
204  !CALL qs_diff_mass_scal_M(pp_mesh, pp_1_LA, one, zero, zero, mode, press_mat(i))
205  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, one, 1.d-10, zero, mode, press_mat(i))
206  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
207  CALL periodic_matrix_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, press_mat(i), pp_1_la)
208  END IF
209  CALL dirichlet_m_parallel(press_mat(i),pp_mode_global_js_d(i)%DIL)
210  CALL init_solver(inputs%my_par_pp,press_ksp(i),press_mat(i),comm_one_d(1),&
211  solver=inputs%my_par_pp%solver,precond=inputs%my_par_pp%precond)
212  !===End PRESSURE
213  ENDDO
214  !===End ASSEMBLING VELOCITY MATRICES
215  ENDIF !===End of once
216 
217  !===Compute NL by FFT at Gauss points
218  tps_tot = user_time()
219  tps_cumul = 0
220  tps = user_time()
221  uext = 2*un-un_m1
222  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v,inputs%precession)
223 
224  IF (inputs%type_pb=='mhd') THEN
225  !===Compute Lorentz force if mhd in quasi-static limit
226  IF (inputs%if_quasi_static_approx) THEN
227  DO i = 1, m_max_c
228  mode = list_mode(i)
229  chmp_mag_aux(:,:,i) = h_b_quasi_static('H', vv_mesh%rr, mode)
230  END DO
231  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,chmp_mag_aux,bn_p2,rotb_b_aux)
232  DO i = 1, m_max_c
233  mode = list_mode(i)
234  chmp_mag_aux(:,:,i) = h_b_quasi_static('B', vv_mesh%rr, mode)
235  END DO
236  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,chmp_mag,chmp_mag_aux,rotb_b)
237  rotb_b = rotb_b + rotb_b_aux
238  ELSE !===Compute Lorentz force if mhd
239  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,chmp_mag,bn_p2,rotb_b)
240  END IF
241  rotv_v = rotv_v - rotb_b
242  END IF
243 
244  IF (inputs%type_pb=='fhd') THEN
245  !===Computate Kelvin force if fhd
246  CALL smb_kelvin_force_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,opt_tempn,chmp_mag,kelvin_force)
247  rotv_v = rotv_v - kelvin_force
248  END IF
249 
250  tps = user_time() - tps; tps_cumul=tps_cumul+tps
251  !===End Compute NL by FFT at Gauss points
252 
253  !===Computation of CFL
254  IF (inputs%verbose_CFL) THEN
255  vel_loc = 0.d0
256  DO i = 1, m_max_c
257  IF (list_mode(i)==0) THEN
258  coeff = 1.d0
259  ELSE
260  coeff = .5d0
261  END IF
262  vel_loc = vel_loc + coeff*(un(:,1,i)**2+un(:,2,i)**2+un(:,3,i)**2 &
263  +un(:,4,i)**2+un(:,5,i)**2+un(:,6,i)**2)
264  END DO
265  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
266  CALL mpi_allreduce(vel_loc,vel_tot,vv_mesh%np,mpi_double_precision, mpi_sum, comm_one_d(2), code)
267  vel_tot = sqrt(vel_tot)
268  cfl = 0.d0
269  DO m = 1, vv_mesh%dom_me
270  vloc = maxval(vel_tot(vv_mesh%jj(:,m)))
271  cfl = max(vloc*dt/min(vv_mesh%hloc(m),0.5d0/inputs%m_max),cfl)
272  END DO
273  CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
274  talk_to_me%CFL=cfl_max
275  talk_to_me%time=time
276  END IF
277  !===End Computation of CFL
278 
279  !===Computation of rhs at Gauss points for every mode
280  IF (present(opt_tempn)) THEN
281  CALL rhs_ns_gauss_3x3(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
282  (4*un-un_m1)/(2*inputs%dt), pn, (4.d0*incpn-incpn_m1)/3.d0, &
283  rotv_v, rhs_gauss, opt_tempn=opt_tempn)
284  ELSE
285  CALL rhs_ns_gauss_3x3(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
286  (4*un-un_m1)/(2*inputs%dt), pn, (4.d0*incpn-incpn_m1)/3.d0, &
287  rotv_v, rhs_gauss)
288  END IF
289  !===End Computation of rhs
290 
291  DO i = 1, m_max_c
292  tps = user_time()
293  mode = list_mode(i)
294 
295  !===Compute phi
296  pn_p1(:,:) = pn(:,:,i)
297 
298  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
299  IF (inputs%LES) THEN
300  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236, &
301  opt_tensor=-visco_entro_sym_grad_u(:,:,:,i))
302  ELSE
303  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
304  END IF
305  !===Periodic boundary conditions
306  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
307  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
308  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
309  END IF
310  !===Axis boundary conditions
311  n1 = SIZE(vv_js_d(1)%DIL)
312  n2 = SIZE(vv_js_d(2)%DIL)
313  n3 = SIZE(vv_js_d(3)%DIL)
314  n123 = n1+n2+n3
315  vel_global_d(i)%DRL(1:n1) = vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
316  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
317  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
318  vel_global_d(i)%DRL(n123+1:) = 0.d0
319  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
320  vel_global_d(i)%DRL(1:n1) = vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
321  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
322  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
323  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
324  tps = user_time() - tps; tps_cumul=tps_cumul+tps
325  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
326 
327  !===Solve linear system for momentum equation
328  tps = user_time()
329  !===Solve system 1, ur_c, ut_s, uz_c
330  nu_mat =2*i-1
331  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
332  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
333  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
334  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
335  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
336  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
337 
338  !===Solve system 2, ur_s, ut_c, uz_s
339  nu_mat = nu_mat + 1
340  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
341  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
342  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
343  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
344  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
345  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
346 
347  tps = user_time() - tps; tps_cumul=tps_cumul+tps
348  !===End Solve linear system for momentum equation
349 
350  !===Assemble divergence of velocity in arrays pb_1, pb_2
351  tps = user_time()
352  CALL qs_01_div_hybrid_2006(vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
353  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
354 
355  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
356  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
357  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
358  END IF
359 
360  !===ATENTION BCs are no longer implemented for pressure
361  !===Boundary condition on axis for pressure
362  pp_global_d(i)%DRL = 0.d0
363  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
364  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
365  !===End boundary condition on axis for pressure
366 
367  !===Solve -LAP(PH3) = -3/(2*dt)*DIV(un_p1)
368  CALL solver(press_ksp(i),pb_1,px_1,reinit=.false.,verbose=inputs%my_par_pp%verbose)
369  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
370  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
371  CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,1))
372 
373  CALL solver(press_ksp(i),pb_2,px_1,reinit=.false.,verbose=inputs%my_par_pp%verbose)
374  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
375  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
376  CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,2))
377 
378  phi = -phi*(1.5d0/dt)
379  tps = user_time() - tps; tps_cumul=tps_cumul+tps
380  !===End Solve -LAP(PH3) = -3/(2*dt)*DIV(un_p1)
381 
382  !===Solve mass matrix for pressure correction
383  tps = user_time()
384  IF (mode==0) THEN
385  CALL solver(mass_ksp0,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
386  ELSE
387  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
388  END IF
389  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
390  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
391  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
392  IF (mode==0) THEN
393  CALL solver(mass_ksp0,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
394  ELSE
395  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
396  END IF
397  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
398  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
399  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
400  !===End solve mass matrix for pressure correction
401 
402  !===Pressure computation
403  DO k=1, 2
404  !pn_p1(:,k) = pn_p1(:,k) + phi(:,k) - div(:,k,i)/Re
405  pn_p1(:,k) = pn_p1(:,k) + phi(:,k) - 2*div(:,k,i)/re - inputs%div_stab_in_ns*div(:,k,i)
406  END DO
407  tps = user_time() - tps; tps_cumul=tps_cumul+tps
408  !===End Pressure computation
409 
410  !===UPDATES
411  tps = user_time()
412  !===Handling of mean pressure
413  IF (mode == 0) THEN
414  CALL moy(comm_one_d(1),pp_mesh, pn_p1(:,1),moyenne)
415  pn_p1(:,1) = pn_p1(:,1)-moyenne
416  ENDIF
417  !===End of handling of mean pressure
418 
419  !===Correction of zero mode
420  IF (mode==0) THEN
421  un_p1(:,2) = 0.d0
422  un_p1(:,4) = 0.d0
423  un_p1(:,6) = 0.d0
424  pn_p1(:,2) = 0.d0
425  END IF
426  !===Correction of zero mode
427 
428  IF (inputs%LES) THEN
429  un_m2(:,:,i) = un_m1(:,:,i)
430  END IF
431  !===UPDATES
432  un_m1(:,:,i) = un(:,:,i)
433  un(:,:,i) = un_p1
434  pn_m1(:,:,i) = pn(:,:,i)
435  pn(:,:,i) = pn_p1
436  incpn_m1(:,:,i) = incpn(:,:,i)
437  incpn(:,:,i) = phi
438  tps = user_time() - tps; tps_cumul=tps_cumul+tps
439  !===End UPDATES
440 
441  ENDDO
442 
443  IF (inputs%verbose_divergence) THEN
444  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
445  talk_to_me%div_L2 = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)/norm
446  talk_to_me%weak_div_L2 = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, div)/norm
447  END IF
448 
449  !===Compute entropy viscosity
450  IF (inputs%LES) THEN
451  IF (present(opt_tempn)) THEN
452  CALL compute_entropy_viscosity(comm_one_d, vv_3_la, vv_mesh, pp_mesh, time, list_mode, vvz_per, &
453  un, un_m1, un_m2, pn_m1, rotv_v, visco_entro_sym_grad_u, opt_tempn)
454  ELSE
455  CALL compute_entropy_viscosity(comm_one_d, vv_3_la, vv_mesh, pp_mesh, time, list_mode, vvz_per, &
456  un, un_m1, un_m2, pn_m1, rotv_v, visco_entro_sym_grad_u)
457  END IF
458  END IF
459  !===End Compute entropy viscosity
460 
461  tps_tot = user_time() - tps_tot
462 
463  END SUBROUTINE bdf2_ns_stress_bc_with_u
464  !============================================
465 
466 !!$ SUBROUTINE inject(jj_c, jj_f, pp_c, pp_f)
467 !!$ IMPLICIT NONE
468 !!$ INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_c, jj_f
469 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_c
470 !!$ REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_f
471 !!$ REAL(KIND=8) :: half = 0.5
472 !!$ INTEGER:: m
473 !!$ IF (SIZE(jj_c,1)==3) THEN
474 !!$ DO m = 1, SIZE(jj_f,2)
475 !!$ pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
476 !!$ pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
477 !!$ pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
478 !!$ pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
479 !!$ END DO
480 !!$ ELSE
481 !!$ DO m = 1, SIZE(jj_f,2)
482 !!$ pp_f(jj_f(1:4,m)) = pp_c(jj_c(:,m))
483 !!$ END DO
484 !!$ pp_f(jj_f(5,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(4,:)))*half
485 !!$ pp_f(jj_f(6,:)) = (pp_c(jj_c(4,:)) + pp_c(jj_c(2,:)))*half
486 !!$ pp_f(jj_f(7,:)) = (pp_c(jj_c(2,:)) + pp_c(jj_c(3,:)))*half
487 !!$ pp_f(jj_f(8,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(4,:)))*half
488 !!$ pp_f(jj_f(9,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(1,:)))*half
489 !!$ pp_f(jj_f(10,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(2,:)))*half
490 !!$ END IF
491 !!$ END SUBROUTINE inject
492 
493 
494 ! cas PRECESSION 28/07/09
495  SUBROUTINE smb_cross_prod_gauss_sft_par(communicator,mesh,list_mode,V_in,V_out,precession_in)
496  !=================================
497  USE sft_parallele
498  USE chaine_caractere
499  USE input_data
500  USE def_type_mesh
501  IMPLICIT NONE
502 
503  TYPE(mesh_type), INTENT(IN) :: mesh
504  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
505  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v_in
506  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: v_out
507  LOGICAL, INTENT(IN) :: precession_in
508  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: rotv, w, om
509  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
510  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
511  INTEGER :: m, l , i, mode, index, k
512  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vs, omega_s
513  REAL(KIND=8) :: ray
514  REAL(KIND=8) :: tps, pi
515  REAL(KIND=8), DIMENSION(3) :: temps
516  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: omega
517  LOGICAL, SAVE :: once=.true.
518 !===FOR FFT_PAR_CROSS_PROD_DCL
519  INTEGER :: nb_procs, bloc_size, m_max_pad, code
520 #include "petsc/finclude/petsc.h"
521  mpi_comm :: communicator
522 
523  IF (once) THEN
524  once = .false.
525  ALLOCATE(omega(mesh%np,6,SIZE(list_mode)))
526  omega = 0.d0
527  pi = acos(-1.d0)
528  DO i=1, SIZE(list_mode)
529  IF (list_mode(i) == 1) THEN
530  !precession selon un axe penche d un angle angle_s_pi*PI
531  omega(:,1,i) = inputs%taux_precession * sin(inputs%angle_s_pi*pi)
532  omega(:,4,i) = -inputs%taux_precession * sin(inputs%angle_s_pi*pi)
533  ELSE IF (list_mode(i) == 0) THEN
534  !precession selon un axe penche d un angle angle_s_pi*PI
535  omega(:,5,i) = inputs%taux_precession * cos(inputs%angle_s_pi*pi)
536  ENDIF
537  ENDDO
538  ENDIF ! fin du once
539 
540  tps = user_time()
541  DO i = 1, SIZE(list_mode)
542  mode = list_mode(i)
543  index = 0
544  DO m = 1, mesh%me
545  j_loc = mesh%jj(:,m)
546  DO k = 1, 6
547  vs(:,k) = v_in(j_loc,k,i)
548  omega_s(:,k) = omega(mesh%jj(:,m),k,i)
549  END DO
550 
551  DO l = 1, mesh%gauss%l_G
552  index = index + 1
553  dw_loc = mesh%gauss%dw(:,:,l,m)
554 
555  !===Compute radius of Gauss point
556  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
557 
558  !-----------------vitesse sur les points de Gauss---------------------------
559  w(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
560  w(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
561  w(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
562 
563  w(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
564  w(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
565  w(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
566  !-----------------vecteur rotation sur les points de Gauss---------------------------
567  om(index,1,i) = sum(omega_s(:,1)*mesh%gauss%ww(:,l))
568  om(index,3,i) = sum(omega_s(:,3)*mesh%gauss%ww(:,l))
569  om(index,5,i) = sum(omega_s(:,5)*mesh%gauss%ww(:,l))
570 
571  om(index,2,i) = sum(omega_s(:,2)*mesh%gauss%ww(:,l))
572  om(index,4,i) = sum(omega_s(:,4)*mesh%gauss%ww(:,l))
573  om(index,6,i) = sum(omega_s(:,6)*mesh%gauss%ww(:,l))
574  !-----------------rotational sur les points de Gauss---------------------------
575  !coeff sur les cosinus
576  rotv(index,1,i) = mode/ray*w(index,6,i) &
577  -sum(vs(:,3)*dw_loc(2,:))
578  rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
579  -sum(vs(:,6)*dw_loc(1,:))
580  rotv(index,5,i) = 1/ray*w(index,3,i) &
581  +sum(vs(:,3)*dw_loc(1,:)) &
582  -mode/ray*w(index,2,i)
583  !coeff sur les sinus
584  rotv(index,2,i) =-mode/ray*w(index,5,i) &
585  -sum(vs(:,4)*dw_loc(2,:))
586  rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
587  -sum(vs(:,5)*dw_loc(1,:))
588  rotv(index,6,i) = 1/ray*w(index,4,i) &
589  +sum(vs(:,4)*dw_loc(1,:))&
590  +mode/ray*w(index,1,i)
591  ENDDO
592  ENDDO
593  END DO
594 
595  !JLG, FL, July 23 2010, There was a Bug here
596  IF (.NOT.precession_in) THEN
597  om=0.d0
598  END IF
599  !
600  ! cas PRECESSION 28/07/09
601  !take care to the -F= 2 Om x U term in LHS
602  rotv = rotv + 2.d0 * om
603  temps = 0
604 
605 
606  CALL mpi_comm_size(communicator, nb_procs, code)
607  bloc_size = SIZE(rotv,1)/nb_procs+1
608  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
609  CALL fft_par_cross_prod_dcl(communicator, rotv, w, v_out, nb_procs, bloc_size, m_max_pad, temps)
610  tps = user_time() - tps
611 
612  END SUBROUTINE smb_cross_prod_gauss_sft_par
613 
614  SUBROUTINE smb_curlh_cross_b_gauss_sft_par(communicator,mesh,list_mode,V_in,W_in,V_out)
615  !=================================
616  USE sft_parallele
617  USE chaine_caractere
618  USE input_data
619  USE def_type_mesh
620  IMPLICIT NONE
621  TYPE(mesh_type), INTENT(IN) :: mesh
622  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
623  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v_in, w_in
624  REAL(KIND=8), DIMENSION(:,:,:) :: v_out
625  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: rotv, v, w
626  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
627  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
628  INTEGER :: m, l , i, mode, index, k
629  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vs, ws
630  REAL(KIND=8) :: ray
631  REAL(KIND=8) :: tps
632  REAL(KIND=8), DIMENSION(3) :: temps
633 !===FOR FFT_PAR_CROSS_PROD_DCL
634  INTEGER :: nb_procs, bloc_size, m_max_pad, code
635 #include "petsc/finclude/petsc.h"
636  mpi_comm :: communicator
637 
638  tps = user_time()
639  DO i = 1, SIZE(list_mode)
640  mode = list_mode(i)
641  index = 0
642  DO m = 1, mesh%me
643  j_loc = mesh%jj(:,m)
644  DO k = 1, 6
645  vs(:,k) = v_in(j_loc,k,i)
646  ws(:,k) = w_in(j_loc,k,i)
647  END DO
648 
649  DO l = 1, mesh%gauss%l_G
650  index = index + 1
651  dw_loc = mesh%gauss%dw(:,:,l,m)
652 
653  !===Compute radius of Gauss point
654  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
655 
656  !-----------------vitesse sur les points de Gauss---------------------------
657  w(index,1,i) = sum(ws(:,1)*mesh%gauss%ww(:,l))
658  w(index,3,i) = sum(ws(:,3)*mesh%gauss%ww(:,l))
659  w(index,5,i) = sum(ws(:,5)*mesh%gauss%ww(:,l))
660 
661  w(index,2,i) = sum(ws(:,2)*mesh%gauss%ww(:,l))
662  w(index,4,i) = sum(ws(:,4)*mesh%gauss%ww(:,l))
663  w(index,6,i) = sum(ws(:,6)*mesh%gauss%ww(:,l))
664  v(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
665  v(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
666  v(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
667 
668  v(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
669  v(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
670  v(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
671  !-----------------rotational sur les points de Gauss---------------------------
672  !coeff sur les cosinus
673  rotv(index,1,i) = mode/ray*v(index,6,i) &
674  -sum(vs(:,3)*dw_loc(2,:))
675  rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
676  -sum(vs(:,6)*dw_loc(1,:))
677  rotv(index,5,i) = 1/ray*v(index,3,i) &
678  +sum(vs(:,3)*dw_loc(1,:)) &
679  -mode/ray*v(index,2,i)
680  !coeff sur les sinus
681  rotv(index,2,i) =-mode/ray*v(index,5,i) &
682  -sum(vs(:,4)*dw_loc(2,:))
683  rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
684  -sum(vs(:,5)*dw_loc(1,:))
685  rotv(index,6,i) = 1/ray*v(index,4,i) &
686  +sum(vs(:,4)*dw_loc(1,:))&
687  +mode/ray*v(index,1,i)
688  ENDDO
689  ENDDO
690  END DO
691 
692  !tps = user_time() - tps
693  !WRITE(*,*) ' Tps dans la grande boucle', tps
694  !tps = user_time()
695  temps = 0
696 
697 
698  CALL mpi_comm_size(communicator, nb_procs, code)
699  bloc_size = SIZE(rotv,1)/nb_procs+1
700  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
701 
702  CALL fft_par_cross_prod_dcl(communicator, rotv, w, v_out, nb_procs, bloc_size, m_max_pad, temps)
703  tps = user_time() - tps
704  !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
705  !write(*,*) ' Temps de Comm ', temps(1)
706  !write(*,*) ' Temps de Calc ', temps(2)
707  !write(*,*) ' Temps de Chan ', temps(3)
708  !DEALLOCATE(Om, W, RotV)
709 
710  END SUBROUTINE smb_curlh_cross_b_gauss_sft_par
711 
712  SUBROUTINE smb_kelvin_force_gauss_sft_par(communicator,mesh,list_mode,scal_in,vect_in,vect_out)
713  USE gauss_points
714  USE sft_parallele
715  USE def_type_mesh
716  IMPLICIT NONE
717  TYPE(mesh_type), INTENT(IN) :: mesh
718  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
719  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: scal_in, vect_in
720  REAL(KIND=8), DIMENSION(:,:,:) :: vect_out
721  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: scal_in_gauss
722  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: vect_in_square
723  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: grad_vect_in_square
724  INTEGER :: i, mode, index, m, k, l
725  INTEGER, DIMENSION(:,:), POINTER :: jj
726  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
727  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: scal_in_loc, vect_in_square_loc
728  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
729  REAL(KIND=8) :: rad
730  !===FOR FFT
731  INTEGER :: nb_procs, bloc_size, m_max_pad, code
732 #include "petsc/finclude/petsc.h"
733  mpi_comm :: communicator
734 
735  ! We compute the Kelvin force: Kelvin force = coeff(T) * grad(H**2/2)
736 
737  !===nb_procs and m_max_pad calculus for FFT
738  CALL mpi_comm_size(communicator, nb_procs, code)
739  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
740 
741  !===T on Gauss nodes computation
742  CALL gauss(mesh)
743  jj => mesh%jj
744  DO i = 1, SIZE(list_mode)
745  index = 0 ! global index of Gauss node
746  DO m = 1, mesh%dom_me
747  j_loc = jj(:,m)
748  DO k = 1, 2
749  scal_in_loc(:,k) = scal_in(j_loc,k,i)
750  END DO
751  DO l=1, l_g
752  index = index + 1
753  dw_loc = dw(:,:,l,m)
754  scal_in_gauss(index,1,i) = sum(scal_in_loc(:,1)*ww(:,l))
755  scal_in_gauss(index,2,i) = sum(scal_in_loc(:,2)*ww(:,l))
756  END DO
757  END DO
758  END DO
759 
760  !===H**2 on nodes computation
761  bloc_size = SIZE(vect_in,1)/nb_procs+1
762  CALL fft_par_dot_prod_dcl(communicator, vect_in, vect_in, vect_in_square, nb_procs, bloc_size, m_max_pad)
763 
764  !===grad(H**2) on Gauss nodes computation
765  DO i = 1, SIZE(list_mode)
766  mode = list_mode(i)
767  index = 0
768  DO m = 1, mesh%dom_me
769  j_loc = jj(:,m)
770  DO k = 1, 2
771  vect_in_square_loc(:,k) = vect_in_square(j_loc,k,i)
772  END DO
773  DO l=1, l_g
774  index = index + 1
775  dw_loc = dw(:,:,l,m)
776  rad = sum(mesh%rr(1,jj(:,m))*ww(:,l)) ! radius of Gauss node
777  grad_vect_in_square(index,1,i) = sum(vect_in_square_loc(:,1)*dw_loc(1,:))
778  grad_vect_in_square(index,2,i) = sum(vect_in_square_loc(:,2)*dw_loc(1,:))
779  grad_vect_in_square(index,3,i) = mode/rad * sum(vect_in_square_loc(:,2)*ww(:,l))
780  grad_vect_in_square(index,4,i) = - mode/rad * sum(vect_in_square_loc(:,1)*ww(:,l))
781  grad_vect_in_square(index,5,i) = sum(vect_in_square_loc(:,1)*dw_loc(2,:))
782  grad_vect_in_square(index,6,i) = sum(vect_in_square_loc(:,2)*dw_loc(2,:))
783  END DO
784  END DO
785  END DO
786 
787  !===Kelvin force = coeff(T) * grad(H**2/2) on Gauss nodes computation
788  bloc_size = SIZE(scal_in_gauss,1)/nb_procs+1
789  CALL fft_kelvin_force(communicator, 0.5*grad_vect_in_square, scal_in_gauss, vect_out, nb_procs, bloc_size, m_max_pad)
790 
791  END SUBROUTINE smb_kelvin_force_gauss_sft_par
792 
793  SUBROUTINE moy(communicator,mesh,p,RESLT)
794  !===========================
795  !moyenne
796  USE def_type_mesh
797  IMPLICIT NONE
798  TYPE(mesh_type) :: mesh
799  REAL(KIND=8), DIMENSION(:) , INTENT(IN) :: p
800  REAL(KIND=8) , INTENT(OUT) :: reslt
801  REAL(KIND=8) :: vol_loc, vol_out, r_loc, r_out
802  INTEGER :: m, l , i , ni, code
803  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
804  REAL(KIND=8) :: ray
805 #include "petsc/finclude/petsc.h"
806  mpi_comm :: communicator
807  r_loc = 0.d0
808  vol_loc = 0.d0
809 
810  DO m = 1, mesh%dom_me
811  j_loc = mesh%jj(:,m)
812  DO l = 1, mesh%gauss%l_G
813  ray = 0
814  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
815  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
816  END DO
817 
818  r_loc = r_loc + sum(p(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
819  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
820 
821  ENDDO
822  ENDDO
823 
824  CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
825  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
826  reslt = r_out / vol_out
827 
828  END SUBROUTINE moy
829 
830 END MODULE subroutine_ns_with_u
tn_axi
Definition: tn_axi.f90:5
sft_parallele::fft_par_dot_prod_dcl
subroutine, public fft_par_dot_prod_dcl(communicator, V1_in, V2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:327
sft_parallele::fft_kelvin_force
subroutine, public fft_kelvin_force(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:3916
fem_m_axi
Definition: fem_M_axi.f90:6
fem_rhs_axi
Definition: fem_rhs_axi.f90:1
subroutine_ns_with_m::smb_curlh_cross_b_gauss_sft_par
subroutine smb_curlh_cross_b_gauss_sft_par(communicator, mesh, list_mode, V_in, W_in, V_out)
Definition: sub_ns_with_momentum.f90:711
solve_petsc
Definition: solver.f90:1
subroutine_ns_with_u::smb_kelvin_force_gauss_sft_par
subroutine smb_kelvin_force_gauss_sft_par(communicator, mesh, list_mode, scal_in, vect_in, vect_out)
Definition: sub_ns_with_velocity.f90:712
sft_parallele
Definition: fft_parallel.f90:4
rhs_gauss_computing
Definition: computing_rhs_gauss.f90:1
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
st_matrix::create_my_ghost
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
gauss_points
Definition: associate_gauss.f90:4
sub_plot
Definition: sub_plot.f90:4
dir_nodes_petsc::dirichlet_rhs
subroutine dirichlet_rhs(js_D, bs_D, b)
Definition: dir_nodes_petsc.f90:170
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
dir_nodes_petsc::vector_glob_js_d
subroutine vector_glob_js_d(vv_mesh, list_mode, vv_3_LA, vv_list_dirichlet_sides, vv_js_D, vv_mode_global_js_D)
Definition: dir_nodes_petsc.f90:188
st_matrix
Definition: st_csr.f90:4
dir_nodes_petsc::scalar_glob_js_d
subroutine scalar_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_mode_global_js_D)
Definition: dir_nodes_petsc.f90:312
solve_petsc::init_solver
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
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:5
boundary_generic::h_b_quasi_static
real(kind=8) function, dimension(size(rr, 2), 6), public h_b_quasi_static(char_h_b, rr, m)
Definition: condlim.f90:276
dir_nodes_petsc
Definition: dir_nodes_petsc.f90:4
dyn_line
Definition: dyn_line.f90:4
fem_m_axi::qs_diff_mass_scal_m
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:127
rhs_gauss_computing::rhs_ns_gauss_3x3
subroutine, public rhs_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, pn_inc, rotv_v, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:6
solve_petsc::solver
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
dir_nodes_petsc::dirichlet_m_parallel
subroutine dirichlet_m_parallel(matrix, glob_js_D)
Definition: dir_nodes_petsc.f90:149
fem_s_axi::moy
subroutine moy(mesh, p, RESULT)
Definition: fem_s_direct_axi.f90:34
dyn_line::dyn_int_line
Definition: dyn_line.f90:5
subroutine_ns_with_u
Definition: sub_ns_with_velocity.f90:5
input_data
Definition: read_my_data.f90:200
def_type_mesh::periodic_type
Definition: def_type_mesh.f90:54
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:23
boundary_test_1::vv_exact
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
Definition: condlim_test_1.f90:301
chaine_caractere
Definition: chaine_caractere_par.f90:4
fem_tn_axi
Definition: fem_tn_axi.f90:4
tn_axi::norm_sf
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:38
entropy_viscosity
Definition: entropy_viscosity.f90:1
fem_s_axi::qs_01_div_hybrid_2006
subroutine qs_01_div_hybrid_2006(uu_mesh, pp_mesh, mode, gg, u0_c)
Definition: fem_s_direct_axi.f90:1100
subroutine_ns_with_u::bdf2_ns_stress_bc_with_u
subroutine, public bdf2_ns_stress_bc_with_u(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, chmp_mag, Bn_p2, opt_tempn)
Definition: sub_ns_with_velocity.f90:11
my_util::user_time
real(kind=8) function user_time()
Definition: my_util.f90:7
periodic::periodic_rhs_petsc
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
Definition: prep_mesh_periodic.f90:455
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:6
subroutine_ns_with_u::smb_cross_prod_gauss_sft_par
subroutine smb_cross_prod_gauss_sft_par(communicator, mesh, list_mode, V_in, V_out, precession_in)
Definition: sub_ns_with_velocity.f90:495
rhs_para_assembling
Definition: assembling_rhs.f90:1
dyn_line::dyn_real_line
Definition: dyn_line.f90:9
sft_parallele::fft_par_cross_prod_dcl
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:144
periodic::periodic_matrix_petsc
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
Definition: prep_mesh_periodic.f90:377
st_matrix::extract
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:33
solve_petsc::create_local_petsc_matrix
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:142
fem_m_axi::qs_diff_mass_vect_3x3_divpenal
subroutine qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
Definition: fem_M_axi.f90:669