SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
sub_ns_with_momentum.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 !
6  USE my_util
8  PRIVATE
9 CONTAINS
10 
11  SUBROUTINE three_level_ns_tensor_sym_with_m(comm_one_d, time, vv_3_LA, pp_1_LA, &
12  dt, re, list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, &
13  hn_p2, bn_p2, density_m1, density, density_p1, visco_dyn, tempn, level_set, level_set_p1, &
14  visc_entro_level)
15 
16  !==============================
17  USE def_type_mesh
18  USE fem_m_axi
19  USE fem_rhs_axi
20  USE fem_tn_axi
21  USE dir_nodes_petsc
22  USE periodic
23  USE st_matrix
24  USE solve_petsc
25  USE dyn_line
26  USE boundary
28  USE sub_plot
29  USE st_matrix
30  USE input_data
31  USE sft_parallele
32  USE tn_axi
33  USE verbose
35  USE subroutine_mass
36  IMPLICIT NONE
37  REAL(KIND=8) :: time, dt, re
38  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
39  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
40  TYPE(petsc_csr_la) :: vv_3_la, pp_1_la
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) :: density_m1, density, density_p1
45  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: level_set, level_set_p1
46  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visco_dyn
47  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tempn
48  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: hn_p2, bn_p2
49  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: hn_p2_aux
50  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: visc_entro_level
51  INTEGER, SAVE :: m_max_c
52  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: pp_global_d
53  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: pp_mode_global_js_d
54  TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_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), SAVE :: mu_bar, nu_bar, rho_bar, sqrt_2d_vol
60  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: momentum, momentum_m1, momentum_m2
61  INTEGER, SAVE :: bloc_size, m_max_pad, nb_procs
62  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: rotb_b_m1
63  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: visc_grad_vel_m1
64  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: tensor_m1
65  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:), SAVE :: visc_entro_real
66  !----------FIN SAVE--------------------------------------------------------------------
67 
68  !----------Declaration sans save-------------------------------------------------------
69  INTEGER, POINTER, DIMENSION(:) :: pp_1_ifrom, vv_3_ifrom
70  INTEGER :: i, k, m, n
71  INTEGER :: code,nu_mat, mode
72  REAL(KIND=8) :: moyenne
73  !allocations des variables locales
74  REAL(KIND=8), DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: div
75  REAL(KIND=8), DIMENSION(pp_mesh%np, 2) :: pn_p1, phi
76  REAL(KIND=8), DIMENSION(vv_mesh%np, 2) :: p_p2
77  REAL(KIND=8), DIMENSION(vv_mesh%np, 6) :: un_p1, src
78  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rotb_b, rotb_b_aux
79  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: visc_grad_vel
80  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: visc_entro_grad_mom
81  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: tensor_surface_gauss
82  REAL(KIND=8), DIMENSION(3,vv_mesh%np,6,SIZE(list_mode)) :: tensor
83  REAL(KIND=8), DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6) :: visc_grad_vel_ext
84  REAL(KIND=8), DIMENSION(3,vv_mesh%np,6) :: tensor_ext
85  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: uext, momentumext, momentum_exact
86  REAL(KIND=8), DIMENSION(inputs%nb_fluid-1, vv_mesh%np, 2, SIZE(list_mode)) :: level_set_fem_p2
87  REAL(KIND=8), DIMENSION(vv_mesh%np) :: vel_loc, vel_tot
88  REAL(KIND=8), DIMENSION(SIZE(list_mode)) :: normalization_mt
89  REAL(KINd=8) :: vel_tot_max_s,vel_tot_max
90  INTEGER :: n1, n2, n3, n123, nb_inter
91  REAL(KIND=8) :: tps, tps_tot, tps_cumul, coeff, cfl, cfl_max, norm, coeff1_in_momemtum
92  INTEGER :: nb_procs_les, bloc_size_les, m_max_pad_les
93  !April 17th 2008, JLG
94  REAL(KIND=8) :: one, zero, three
95  DATA zero, one, three/0.d0,1.d0,3.d0/
96 
97  !Communicators for Petsc, in space and Fourier------------------------------
98 #include "petsc/finclude/petsc.h"
99  petscerrorcode :: ierr
100  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
101  mat, DIMENSION(:), POINTER, SAVE :: vel_mat
102  mat, DIMENSION(:), POINTER, SAVE :: press_mat
103  mat, SAVE :: mass_mat, mass_mat0
104  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
105  vec, SAVE :: pb_1, pb_2, px_1, px_1_ghost
106  ksp, DIMENSION(:), POINTER, SAVE :: vel_ksp, press_ksp
107  ksp, SAVE :: mass_ksp, mass_ksp0
108  !------------------------------END OF DECLARATION--------------------------------------
109 
110 
111  IF (once) THEN
112 
113  once = .false.
114 
115  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
116 
117  !-----CREATE PETSC VECTORS AND GHOSTS-----------------------------------------
118  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
119  n = 3*vv_mesh%dom_np
120  CALL veccreateghost(comm_one_d(1), n, &
121  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
122  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
123  CALL vecduplicate(vx_3, vb_3_145, ierr)
124  CALL vecduplicate(vx_3, vb_3_236, ierr)
125 
126  CALL create_my_ghost(pp_mesh,pp_1_la,pp_1_ifrom)
127  n = pp_mesh%dom_np
128  CALL veccreateghost(comm_one_d(1), n, &
129  petsc_determine, SIZE(pp_1_ifrom), pp_1_ifrom, px_1, ierr)
130  CALL vecghostgetlocalform(px_1, px_1_ghost, ierr)
131  CALL vecduplicate(px_1, pb_1, ierr)
132  CALL vecduplicate(px_1, pb_2, ierr)
133  !------------------------------------------------------------------------------
134 
135  !-------------DIMENSIONS-------------------------------------------------------
136  m_max_c = SIZE(list_mode)
137  !------------------------------------------------------------------------------
138 
139  !-------------MOMENTUM INITIALIZATION------------------------------------------
140  ALLOCATE(momentum(SIZE(un,1),SIZE(un,2), SIZE(un,3)),&
141  momentum_m1(SIZE(un,1),SIZE(un,2), SIZE(un,3)), &
142  momentum_m2(SIZE(un,1),SIZE(un,2), SIZE(un,3)))
143  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
144  bloc_size = SIZE(un,1)/nb_procs+1
145  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
146 
147  IF (inputs%if_level_set) THEN
148  CALL fft_scalar_vect_dcl(comm_one_d(2), un_m1, density_m1, momentum_m1, 1, nb_procs, &
149  bloc_size, m_max_pad)
150  CALL fft_scalar_vect_dcl(comm_one_d(2), un, density, momentum, 1, nb_procs, &
151  bloc_size, m_max_pad)
152  ELSE
153  momentum_m1 = un_m1
154  momentum = un
155  END IF
156 
157  ALLOCATE(rotb_b_m1(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)))
158  rotb_b_m1 = 0.d0
159 
160  ALLOCATE(tensor_m1(3,vv_mesh%np,6,SIZE(list_mode)))
161  bloc_size = SIZE(un,1)/nb_procs+1
162  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
163  CALL fft_tensor_dcl(comm_one_d(2), momentum_m1, un_m1, tensor_m1, nb_procs, bloc_size, m_max_pad)
164 
165  ALLOCATE(visc_grad_vel_m1(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)))
166  CALL smb_explicit_diffu_sym(comm_one_d(2), vv_mesh, list_mode, nb_procs, &
167  visco_dyn/re, un_m1, visc_grad_vel_m1)
168 
169  CALL mpi_comm_size(comm_one_d(2), nb_procs_les, code)
170  bloc_size_les = vv_mesh%gauss%l_G*vv_mesh%dom_me/nb_procs_les+1
171  bloc_size_les = vv_mesh%gauss%l_G*(bloc_size_les/vv_mesh%gauss%l_G)+vv_mesh%gauss%l_G
172  m_max_pad_les = 3*SIZE(list_mode)*nb_procs_les/2
173  ALLOCATE(visc_entro_real(2*m_max_pad_les-1,bloc_size_les))
174  visc_entro_real = 0.d0
175 
176  !---------PREPARE pp_js_D ARRAY FOR PRESSURE-----------------------------------
177  !CALL dirichlet_nodes_parallel(pp_mesh, inputs%pp_list_dirichlet_sides, pp_js_D)
178  !------------------------------------------------------------------------------
179  !===PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
180  !===ATTENTION pressure BCs are no longer implemented
181  CALL scalar_glob_js_d(pp_mesh, list_mode, pp_1_la, pp_mode_global_js_d)
182  ALLOCATE(pp_global_d(m_max_c))
183  DO i = 1, m_max_c
184  ALLOCATE(pp_global_d(i)%DRL(SIZE(pp_mode_global_js_d(i)%DIL)))
185  END DO
186 
187  !===PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
188  CALL vector_glob_js_d(vv_mesh, list_mode, vv_3_la, inputs%vv_list_dirichlet_sides, &
189  vv_js_d, vv_mode_global_js_d)
190 
191  ALLOCATE(vel_global_d(m_max_c))
192  DO i = 1, m_max_c
193  ALLOCATE(vel_global_d(i)%DRL(SIZE(vv_mode_global_js_d(i)%DIL)))
194  END DO
195  !===END PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
196  !------------------------------------------------------------------------------
197 
198  !---------------------------------------------------------------------------------
199  !===JLG LC, WE should introduce inputs%LES_coeff1_in_momemtum and inputs%LES_coeff1_in_level
200  IF (inputs%LES.AND.inputs%if_LES_in_momentum) THEN
201  coeff1_in_momemtum = inputs%LES_coeff1
202  ELSE
203  coeff1_in_momemtum = 0.d0
204  END IF
205 
206  !===ASSEMBLE MASS MATRIX
207  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat, clean=.false.)
208  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat)
209  DO i = 1, m_max_c
210  IF (list_mode(i)==0) cycle
211  CALL dirichlet_m_parallel(mass_mat,pp_mode_global_js_d(i)%DIL)
212  END DO
213  CALL init_solver(inputs%my_par_mass,mass_ksp,mass_mat,comm_one_d(1),&
214  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
215 
216  IF (minval(list_mode)==0) THEN
217  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat0, clean=.false.)
218  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat0)
219  DO i = 1, m_max_c
220  IF (list_mode(i).NE.0) cycle
221  CALL dirichlet_m_parallel(mass_mat0,pp_mode_global_js_d(i)%DIL)
222  CALL init_solver(inputs%my_par_mass,mass_ksp0,mass_mat0,comm_one_d(1),&
223  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
224  END DO
225  END IF
226  !===END ASSEMBLE MASS MATRIX
227 
228  !-------------ASSEMBLE VELOCITY MATRICES----------------------
229  ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
230  ALLOCATE(press_mat(m_max_c),press_ksp(m_max_c))
231 
232  ! Definition of nu_bar
233  IF (inputs%if_level_set) THEN
234  !To be tested thoroughly
235  mu_bar = minval(inputs%dyna_visc_fluid)
236  rho_bar = minval(inputs%density_fluid)
237  nu_bar = 0.d0
238  DO n = 1, inputs%nb_fluid
239  nu_bar = max(nu_bar,inputs%dyna_visc_fluid(n)/inputs%density_fluid(n))
240  END DO
241  nu_bar = 2.d0*nu_bar
242  ELSE
243  mu_bar = 1.d0
244  nu_bar = 1.d0
245  rho_bar = 1.d0
246  END IF
247 
248  DO i = 1, m_max_c
249  mode = list_mode(i)
250  !---VELOCITY
251  nu_mat = 2*i-1
252  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
253  IF (inputs%if_moment_bdf2) THEN
254  CALL qs_diff_mass_vect_3x3(1, vv_3_la, vv_mesh, nu_bar/re, three/(2*dt), &
255  coeff1_in_momemtum, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
256  ELSE
257  CALL qs_diff_mass_vect_3x3(1, vv_3_la, vv_mesh, nu_bar/re, one/dt, &
258  coeff1_in_momemtum, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
259  END IF
260  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
261  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
262  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
263 
264  nu_mat = nu_mat+1
265  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
266  IF (inputs%if_moment_bdf2) THEN
267  CALL qs_diff_mass_vect_3x3(2, vv_3_la, vv_mesh, nu_bar/re, three/(2*dt), &
268  coeff1_in_momemtum, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
269  ELSE
270  CALL qs_diff_mass_vect_3x3(2, vv_3_la, vv_mesh, nu_bar/re, one/dt, &
271  coeff1_in_momemtum, inputs%stab_bdy_ns, mode, vel_mat(nu_mat))
272  END IF
273  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
274  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
275  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
276 
277  !---PRESSURE
278  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, press_mat(i), clean=.false.)
279  CALL qs_diff_mass_scal_m(pp_mesh, pp_1_la, one, 1.d-10, zero, mode, press_mat(i))
280  CALL dirichlet_m_parallel(press_mat(i),pp_mode_global_js_d(i)%DIL)
281  CALL init_solver(inputs%my_par_pp,press_ksp(i),press_mat(i),comm_one_d(1),&
282  solver=inputs%my_par_pp%solver,precond=inputs%my_par_pp%precond)
283 
284  ENDDO
285 
286  CALL twod_volume(comm_one_d(1),vv_mesh,sqrt_2d_vol)
287  sqrt_2d_vol = sqrt(sqrt_2d_vol)
288 
289  ENDIF ! Fin du once
290  tps_tot = user_time()
291  tps_cumul = 0
292 
293 
294  !===Compute rhs by FFT at Gauss points
295  tps = user_time()
296 
297  IF (inputs%if_moment_bdf2) THEN
298  uext = 2*un-un_m1
299  IF (inputs%if_level_set) THEN
300  ! BDF2: momentumext = rho_np1*(2*un-un_m1)
301  CALL fft_scalar_vect_dcl(comm_one_d(2), uext, density_p1, momentumext, 1,nb_procs, &
302  bloc_size, m_max_pad)
303  ELSE
304  momentumext = 2*momentum - momentum_m1
305  END IF
306  ELSE
307  uext = un
308  IF (inputs%if_level_set) THEN
309  ! BDF1: momentumext = rho_np1*un
310  CALL fft_scalar_vect_dcl(comm_one_d(2), un, density_p1, momentumext, 1, nb_procs, &
311  bloc_size, m_max_pad)
312  ELSE
313  momentumext = momentum
314  END IF
315  END IF
316 
317  !===Precession should be in condlim
318  IF (inputs%precession) THEN
319  CALL error_petsc('for momentum ns : precession should be in condlim')
320  END IF
321 
322  !===Compute Lorentz force if mhd
323  IF (inputs%type_pb=='mhd') THEN
324  !===Compute Lorentz force if mhd in quasi-static limit
325  IF (inputs%if_quasi_static_approx) THEN
326  DO i = 1, m_max_c
327  mode = list_mode(i)
328  hn_p2_aux(:,:,i) = h_b_quasi_static('H', vv_mesh%rr, mode)
329  END DO
330  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,hn_p2_aux,bn_p2,rotb_b_aux)
331  DO i = 1, m_max_c
332  mode = list_mode(i)
333  hn_p2_aux(:,:,i) = h_b_quasi_static('B', vv_mesh%rr, mode)
334  END DO
335  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,hn_p2,hn_p2_aux,rotb_b)
336  rotb_b = rotb_b + rotb_b_aux
337  ELSE !===Compute Lorentz force if mhd
338  CALL smb_curlh_cross_b_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,hn_p2,bn_p2,rotb_b)
339  END IF
340  ELSE
341  rotb_b = 0.d0
342  END IF
343  !===End compute Lorentz force if mhd
344 
345  IF (inputs%if_level_set) THEN
346  !===Compute visco_dyn/Re*Grad(un)
347  CALL smb_explicit_diffu_sym(comm_one_d(2), vv_mesh, list_mode, nb_procs, &
348  visco_dyn/re, un, visc_grad_vel)
349  !===End compute visco_dyn/Re*Grad(un)
350 
351  !===Compute (-LES_coeff1+visc_entro)*Grad(momentumext)
352  IF (inputs%LES.AND.inputs%if_LES_in_momentum) THEN
353  CALL smb_explicit_les(comm_one_d(2), vv_mesh, list_mode, nb_procs, &
354  visc_entro_real, momentumext, visc_entro_grad_mom)
355  ELSE
356  visc_entro_grad_mom=0.d0
357  END IF
358  !===End compute (-LES_coeff1+visc_entro)*Grad(momentumext)
359 
360  IF (inputs%if_surface_tension) THEN
361  !===Compute coeff_surface*Grad(level_set):Grad(level_set)
362  IF (inputs%if_level_set_P2) THEN
363  CALL smb_surface_tension(comm_one_d(2), vv_mesh, list_mode, nb_procs, &
364  level_set_p1, tensor_surface_gauss)
365  ELSE
366  DO nb_inter = 1, inputs%nb_fluid-1
367  DO i = 1, SIZE(list_mode)
368  DO k = 1, 2
369  CALL inject_p1_p2(pp_mesh%jj, vv_mesh%jj, level_set_p1(nb_inter,:,k,i), &
370  level_set_fem_p2(nb_inter,:,k,i))
371  END DO
372  END DO
373  END DO
374  CALL smb_surface_tension(comm_one_d(2), vv_mesh, list_mode, nb_procs, &
375  level_set_fem_p2, tensor_surface_gauss)
376  END IF
377  ELSE
378  tensor_surface_gauss = 0.d0
379  END IF
380  ELSE
381  visc_grad_vel = 0.d0
382  tensor_surface_gauss = 0.d0
383  END IF
384 
385  !===Compute tensor product of momentum by velocity
386  CALL fft_tensor_dcl(comm_one_d(2), momentum, un, tensor, nb_procs, bloc_size, m_max_pad)
387 
388  !===PREPARE BOUNDARY CONDITION FOR MOMENTUM
389  CALL momentum_dirichlet(comm_one_d(2), vv_mesh, list_mode, time, nb_procs,density_p1, &
390  momentum_exact, vv_js_d)
391 
392  tps = user_time() - tps; tps_cumul=tps_cumul+tps
393  !WRITE(*,*) ' Tps fft vitesse', tps
394 
395  !------------DEBUT BOUCLE SUR LES MODES----------------
396  DO i = 1, m_max_c
397  mode = list_mode(i)
398  !===Compute phi
399  tps = user_time()
400  !jan 29 2007
401  pn_p1(:,:) = pn(:,:,i)
402  IF (inputs%if_moment_bdf2) THEN
403  phi = pn_p1(:,:) + (4.d0 * incpn(:,:,i) - incpn_m1(:,:,i))/3.d0
404  ELSE
405  phi = pn_p1(:,:) + incpn(:,:,i)
406  END IF
407  !jan 29 2007
408 
409  !===Inject pressure P1 -> P2
410  DO k = 1, 2
411  CALL inject(pp_mesh%jj, vv_mesh%jj, phi(:,k), p_p2(:,k))
412  ENDDO
413 
414  !===Prediction step
415  DO k=1,6
416  IF (inputs%if_temperature) THEN
417  src(:,k) = source_in_ns_momentum(k, vv_mesh%rr, mode, i, time, re, 'ns', &
418  opt_density=density_p1, opt_tempn=tempn)
419  ELSE
420  src(:,k) = source_in_ns_momentum(k, vv_mesh%rr, mode, i, time, re, 'ns', &
421  opt_density=density_p1)
422  END IF
423  END DO
424 
425  IF (inputs%if_moment_bdf2) THEN
426  tensor_ext = 2*tensor(:,:,:,i) - tensor_m1(:,:,:,i)
427  visc_grad_vel_ext = 2*visc_grad_vel(:,:,:,i) - visc_grad_vel_m1(:,:,:,i)
428 
429  !===Update terms at m1
430  rotb_b_m1(:,:,i) = rotb_b(:,:,i)
431  tensor_m1(:,:,:,i) = tensor(:,:,:,i)
432  visc_grad_vel_m1(:,:,:,i) = visc_grad_vel(:,:,:,i)
433  CALL qs_ns_momentum_stab_3x3(vv_mesh, vv_3_la, mode, src, &
434  (4*momentum(:,:,i)-momentum_m1(:,:,i))/(2*dt), p_p2(:,:), &
435  vb_3_145, vb_3_236, rotb_b(:,:,i), tensor_ext,&
436  tensor_surface_gauss(:,:,:,i), nu_bar/re, momentumext(:,:,i), &
437  -visc_grad_vel_ext-visc_entro_grad_mom(:,:,:,i))
438  ELSE
439  CALL qs_ns_momentum_stab_3x3(vv_mesh, vv_3_la, mode, src, &
440  momentum(:,:,i)/dt, p_p2(:,:), &
441  vb_3_145, vb_3_236, rotb_b(:,:,i), tensor(:,:,:,i),&
442  tensor_surface_gauss(:,:,:,i), nu_bar/re, momentumext(:,:,i), &
443  -visc_grad_vel(:,:,:,i)-visc_entro_grad_mom(:,:,:,i))
444  END IF
445 
446  !===Axis boundary conditions
447  n1 = SIZE(vv_js_d(1)%DIL)
448  n2 = SIZE(vv_js_d(2)%DIL)
449  n3 = SIZE(vv_js_d(3)%DIL)
450  n123 = n1+n2+n3
451  vel_global_d(i)%DRL(1:n1) = momentum_exact(vv_js_d(1)%DIL,1,i)
452  vel_global_d(i)%DRL(n1+1:n1+n2) = momentum_exact(vv_js_d(2)%DIL,4,i)
453  vel_global_d(i)%DRL(n1+n2+1:n123)= momentum_exact(vv_js_d(3)%DIL,5,i)
454  vel_global_d(i)%DRL(n123+1:) = 0.d0
455  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
456  vel_global_d(i)%DRL(1:n1) = momentum_exact(vv_js_d(1)%DIL,2,i)
457  vel_global_d(i)%DRL(n1+1:n1+n2) = momentum_exact(vv_js_d(2)%DIL,3,i)
458  vel_global_d(i)%DRL(n1+n2+1:n123)= momentum_exact(vv_js_d(3)%DIL,6,i)
459  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
460  tps = user_time() - tps; tps_cumul=tps_cumul+tps
461  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
462 
463 
464  tps = user_time() - tps; tps_cumul=tps_cumul+tps
465  !WRITE(*,*) ' Tps second membre vitesse', tps
466  !-------------------------------------------------------------------------------------
467 
468  !--------------------INVERSION DE L'OPERATEUR 1 --------------
469  tps = user_time()
470  !Solve system 1, ur_c, ut_s, uz_c
471  nu_mat =2*i-1
472  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
473  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
474  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
475  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
476  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
477  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
478 
479  !Solve system 2, ur_s, ut_c, uz_s
480  nu_mat = nu_mat + 1
481  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
482  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
483  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
484  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
485  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
486  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
487 
488  tps = user_time() - tps; tps_cumul=tps_cumul+tps
489  !WRITE(*,*) ' Tps solution des pb de vitesse', tps, 'for mode ', mode
490  !-------------------------------------------------------------------------------------
491 
492  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
493  IF (mode==0) THEN
494  un_p1(:,2) = 0.d0
495  un_p1(:,4) = 0.d0
496  un_p1(:,6) = 0.d0
497  pn_p1(:,2) = 0.d0
498  END IF
499  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
500 
501  momentum_m2(:,:,i) = momentum_m1(:,:,i)
502  momentum_m1(:,:,i) = momentum(:,:,i)
503  momentum(:,:,i) = un_p1
504 
505  END DO
506 
507  !---------------UPDATE VELOCITY---------------------
508  un_m1 = un
509  IF (inputs%if_level_set) THEN
510  CALL fft_scalar_vect_dcl(comm_one_d(2), momentum, density_p1, un, 2, nb_procs, &
511  bloc_size, m_max_pad)
512  ELSE
513  un = momentum
514  END IF
515 
516 
517  DO i = 1, m_max_c
518  mode = list_mode(i)
519  !===Compute phi
520  tps = user_time()
521  !jan 29 2007
522  pn_p1(:,:) = pn(:,:,i)
523  !jan 29 2007
524 
525  !---------------SOLUTION OF THE POISSON EQUATION--------------
526  ! BDF2 : solve -LAP(PH3) = -3*rho_bar/(2*dt)*DIV(un_p1)
527  ! BDF1 : solve -LAP(PH3) = -rho_bar/dt*DIV(un_p1)
528  tps = user_time()
529  CALL qs_01_div_hybrid_2006(vv_mesh, pp_mesh, pp_1_la, mode, un(:,:,i), pb_1, pb_2)
530  !pb_1, and pb_2 are petsc vectors for the rhs divergence
531 
532  !===ATENTION BCs are no longer implemented for pressure
533  !===Boundary condition on axis for pressure
534  pp_global_d(i)%DRL = 0.d0
535  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
536  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
537  !===End boundary condition on axis for pressure
538 
539  CALL solver(press_ksp(i),pb_1,px_1,reinit=.false.,verbose=inputs%my_par_pp%verbose)
540  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
541  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
542  CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,1))
543 
544  CALL solver(press_ksp(i),pb_2,px_1,reinit=.false.,verbose=inputs%my_par_pp%verbose)
545  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
546  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
547  CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,2))
548 
549  !Don't forget the -(1.5d0/dt)*min(rh0) coefficient
550  IF (inputs%if_moment_bdf2) THEN
551  phi = -phi*(1.5d0/dt)*rho_bar
552  ELSE
553  phi = -phi/dt*rho_bar
554  END IF
555  tps = user_time() - tps; tps_cumul=tps_cumul+tps
556  !WRITE(*,*) ' Tps solution des pb de pression', tps, 'for mode ', mode
557  !-------------------------------------------------------------------------------------
558 
559 
560  !---------------CORRECTION DE LA PRESSION-----------------------
561  tps = user_time()
562  !CALL solver(mass_ksp,pb_1,px_1,reinit=.FALSE.,verbose=inputs%my_par_mass%verbose)
563  IF (mode==0) THEN
564  CALL solver(mass_ksp0,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
565  ELSE
566  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
567  END IF
568  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
569  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
570  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
571 
572  !CALL solver(mass_ksp,pb_2,px_1,reinit=.FALSE.,verbose=inputs%my_par_mass%verbose)
573  IF (mode==0) THEN
574  CALL solver(mass_ksp0,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
575  ELSE
576  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
577  END IF
578  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
579  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
580  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
581 
582  !Pressure computation
583  !jan 29 2007
584  DO k=1, 2
585  pn_p1(:,k) = pn_p1(:,k) + phi(:,k) - div(:,k,i)*(mu_bar/re)
586  END DO
587  !jan 29 2007
588  tps = user_time() - tps; tps_cumul=tps_cumul+tps
589  !WRITE(*,*) ' Tps correction de divergence', tps
590  !-------------------------------------------------------------------------------------
591 
592  !---------------UPDATE PRESSURE---------------------
593  tps = user_time()
594  IF (mode == 0) THEN
595  CALL moy(comm_one_d(1),pp_mesh, pn_p1(:,1),moyenne)
596  pn_p1(:,1) = pn_p1(:,1)-moyenne
597  ENDIF
598 
599  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
600  IF (mode==0) THEN
601  pn_p1(:,2) = 0.d0
602  END IF
603  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
604 
605  pn_m1(:,:,i) = pn(:,:,i)
606  pn(:,:,i) = pn_p1
607 
608  incpn_m1(:,:,i) = incpn(:,:,i)
609  incpn(:,:,i) = phi
610 
611  tps = user_time() - tps; tps_cumul=tps_cumul+tps
612  !WRITE(*,*) ' Tps des updates', tps
613  !-------------------------------------------------------------------------------------
614 
615  ENDDO
616  tps_tot = user_time() - tps_tot
617 
618  !===Verbose divergence of velocity
619  IF (inputs%verbose_divergence) THEN
620  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
621  talk_to_me%div_L2 = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)/norm
622  talk_to_me%weak_div_L2 = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, div)/norm
623  END IF
624  !===End verbose divergence of velocity
625 
626  !===Compute vel max and CFL
627  IF (inputs%verbose_CFL.OR.inputs%LES) THEN
628  vel_loc = 0.d0
629  DO i = 1, m_max_c
630  IF (list_mode(i)==0) THEN
631  coeff = 1.d0
632  ELSE
633  coeff = .5d0
634  END IF
635  vel_loc = vel_loc + coeff*(un(:,1,i)**2+un(:,2,i)**2+un(:,3,i)**2+&
636  un(:,4,i)**2+un(:,5,i)**2+un(:,6,i)**2)
637  END DO
638  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
639  CALL mpi_allreduce(vel_loc,vel_tot,vv_mesh%np,mpi_double_precision, mpi_sum, comm_one_d(2), code)
640 
641  vel_tot = sqrt(vel_tot)
642  vel_tot_max_s = maxval(vel_tot)
643  CALL mpi_allreduce(vel_tot_max_s,vel_tot_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
644 
645  !===New normalization (JLG; April 24, 2015)
646  DO i = 1, m_max_c
647  normalization_mt(i) = norm_s(comm_one_d, 'L2', vv_mesh, list_mode, momentum)/(sqrt(2.d0)*sqrt_2d_vol)
648  END DO
649  !===New normalization (JLG; April 24, 2015)
650 
651  !===Computation of CFL
652  IF (inputs%verbose_CFL) THEN
653  cfl = 0
654  DO m = 1, vv_mesh%dom_me
655  cfl = max(vel_tot_max*dt/vv_mesh%hloc(m),cfl)
656  END DO
657  CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
658  talk_to_me%CFL=cfl_max
659  talk_to_me%time=time
660  !IF (my_petscworld_rank==0) THEN
661  ! WRITE(*,'(3(A,e10.3))') ' Time = ', time, ', CFL = ', cfl_max, 'vel_tot_max', vel_tot_max
662  !END IF
663  END IF
664  END IF
665 
666  !===Compute entropy viscosity
667  IF (inputs%LES) THEN
668  CALL compute_entropy_viscosity_mom(comm_one_d, vv_3_la, vv_mesh, pp_mesh, time, list_mode, &
669  momentum, momentum_m1, momentum_m2, pn_m1, un_m1, tensor, visc_grad_vel, tensor_surface_gauss, &
670  rotb_b, visco_dyn, density_m1, density, density_p1, tempn, visc_entro_real, visc_entro_level)
671  ELSE
672  visc_entro_real = 0.d0
673  visc_entro_level = 0.d0
674  END IF
675  !===End compute entropy viscosity
676 
677  END SUBROUTINE three_level_ns_tensor_sym_with_m
678  !============================================
679 
680  SUBROUTINE inject(jj_c, jj_f, pp_c, pp_f)
681  IMPLICIT NONE
682 
683  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_c, jj_f
684  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_c
685  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_f
686  REAL(KIND=8) :: half = 0.5
687  INTEGER:: m
688  IF (SIZE(jj_c,1)==3) THEN
689  DO m = 1, SIZE(jj_f,2)
690  pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
691  pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
692  pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
693  pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
694  END DO
695 
696  ELSE
697  DO m = 1, SIZE(jj_f,2)
698  pp_f(jj_f(1:4,m)) = pp_c(jj_c(:,m))
699  END DO
700  pp_f(jj_f(5,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(4,:)))*half
701  pp_f(jj_f(6,:)) = (pp_c(jj_c(4,:)) + pp_c(jj_c(2,:)))*half
702  pp_f(jj_f(7,:)) = (pp_c(jj_c(2,:)) + pp_c(jj_c(3,:)))*half
703  pp_f(jj_f(8,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(4,:)))*half
704  pp_f(jj_f(9,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(1,:)))*half
705  pp_f(jj_f(10,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(2,:)))*half
706 
707  END IF
708 
709  END SUBROUTINE inject
710 
711  SUBROUTINE smb_curlh_cross_b_gauss_sft_par(communicator,mesh,list_mode,V_in,W_in,V_out)
712  !=================================
713  USE sft_parallele
714  USE chaine_caractere
715  USE input_data
716  USE def_type_mesh
717  IMPLICIT NONE
718  TYPE(mesh_type), INTENT(IN) :: mesh
719  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
720  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v_in, w_in
721  REAL(KIND=8), DIMENSION(:,:,:) :: v_out
722  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: rotv, v, w
723  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
724  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
725  INTEGER :: m, l , i, mode, index, k
726  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vs, ws
727  REAL(KIND=8) :: ray
728  REAL(KIND=8) :: tps
729  REAL(KIND=8), DIMENSION(3) :: temps
730 !===FOR FFT_PAR_CROSS_PROD_DCL
731  INTEGER :: nb_procs, bloc_size, m_max_pad, code
732 #include "petsc/finclude/petsc.h"
733  mpi_comm :: communicator
734 
735  tps = user_time()
736  DO i = 1, SIZE(list_mode)
737  mode = list_mode(i)
738  index = 0
739  DO m = 1, mesh%me
740  j_loc = mesh%jj(:,m)
741  DO k = 1, 6
742  vs(:,k) = v_in(j_loc,k,i)
743  ws(:,k) = w_in(j_loc,k,i)
744  END DO
745 
746  DO l = 1, mesh%gauss%l_G
747  index = index + 1
748  dw_loc = mesh%gauss%dw(:,:,l,m)
749 
750  !===Compute radius of Gauss point
751  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
752 
753  !-----------------vitesse sur les points de Gauss---------------------------
754  w(index,1,i) = sum(ws(:,1)*mesh%gauss%ww(:,l))
755  w(index,3,i) = sum(ws(:,3)*mesh%gauss%ww(:,l))
756  w(index,5,i) = sum(ws(:,5)*mesh%gauss%ww(:,l))
757 
758  w(index,2,i) = sum(ws(:,2)*mesh%gauss%ww(:,l))
759  w(index,4,i) = sum(ws(:,4)*mesh%gauss%ww(:,l))
760  w(index,6,i) = sum(ws(:,6)*mesh%gauss%ww(:,l))
761  v(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
762  v(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
763  v(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
764 
765  v(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
766  v(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
767  v(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
768  !-----------------rotational sur les points de Gauss---------------------------
769  !coeff sur les cosinus
770  rotv(index,1,i) = mode/ray*v(index,6,i) &
771  -sum(vs(:,3)*dw_loc(2,:))
772  rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
773  -sum(vs(:,6)*dw_loc(1,:))
774  rotv(index,5,i) = 1/ray*v(index,3,i) &
775  +sum(vs(:,3)*dw_loc(1,:)) &
776  -mode/ray*v(index,2,i)
777  !coeff sur les sinus
778  rotv(index,2,i) =-mode/ray*v(index,5,i) &
779  -sum(vs(:,4)*dw_loc(2,:))
780  rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
781  -sum(vs(:,5)*dw_loc(1,:))
782  rotv(index,6,i) = 1/ray*v(index,4,i) &
783  +sum(vs(:,4)*dw_loc(1,:))&
784  +mode/ray*v(index,1,i)
785  ENDDO
786  ENDDO
787  END DO
788 
789  !tps = user_time() - tps
790  !WRITE(*,*) ' Tps dans la grande boucle', tps
791  !tps = user_time()
792  temps = 0
793 
794 
795  CALL mpi_comm_size(communicator, nb_procs, code)
796  bloc_size = SIZE(rotv,1)/nb_procs+1
797  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
798 
799  CALL fft_par_cross_prod_dcl(communicator, rotv, w, v_out, nb_procs, bloc_size, m_max_pad, temps)
800  tps = user_time() - tps
801  !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
802  !write(*,*) ' Temps de Comm ', temps(1)
803  !write(*,*) ' Temps de Calc ', temps(2)
804  !write(*,*) ' Temps de Chan ', temps(3)
805  !DEALLOCATE(Om, W, RotV)
806 
807  END SUBROUTINE smb_curlh_cross_b_gauss_sft_par
808 
809  SUBROUTINE smb_explicit_diffu_sym(communicator, mesh, list_mode, nb_procs, visc_dyn, vel, V_out)
810  USE gauss_points
811  USE sft_parallele
812  USE chaine_caractere
813  USE boundary
814  USE input_data
815  IMPLICIT NONE
816  TYPE(mesh_type), TARGET :: mesh
817  INTEGER, INTENT(IN) :: nb_procs
818  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
819  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visc_dyn
820  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vel
821  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: v_out
822  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: grad1_vel, grad2_vel, grad3_vel
823  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: gradt1_vel, gradt2_vel, gradt3_vel
824  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: part_tensor_1, part_tensor_2
825  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: part_visc_sym_grad_1
826  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: part_visc_sym_grad_2
827  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)) :: visc_dyn_gauss
828  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
829  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
830  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: vel_loc
831  REAL(KIND=8) :: ray
832  INTEGER :: m, l , i, mode, index, k
833  INTEGER :: m_max_pad, bloc_size
834 #include "petsc/finclude/petsc.h"
835  mpi_comm :: communicator
836 
837  CALL gauss(mesh)
838 
839  DO i = 1, SIZE(list_mode)
840  mode = list_mode(i)
841  index = 0
842  DO m = 1, mesh%dom_me
843  j_loc = mesh%jj(:,m)
844  DO k = 1, 6
845  vel_loc(:,k) = vel(j_loc,k,i)
846  END DO
847  DO l = 1, l_g
848  index = index + 1
849  dw_loc = dw(:,:,l,m)
850 
851  !===Compute radius of Gauss point
852  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
853 
854  !-----------------Dynamic viscosity on Gauss points---------------------------
855  visc_dyn_gauss(index,1,i) = sum(visc_dyn(j_loc,1,i)*ww(:,l))
856  visc_dyn_gauss(index,2,i) = sum(visc_dyn(j_loc,2,i)*ww(:,l))
857 
858  !-----------------Grad u_r on Gauss points------------------------------------
859  grad1_vel(index,1,i) = sum(vel_loc(:,1)*dw_loc(1,:))
860  grad1_vel(index,2,i) = sum(vel_loc(:,2)*dw_loc(1,:))
861  grad1_vel(index,3,i) = (mode*sum(vel_loc(:,2)*ww(:,l)) - sum(vel_loc(:,3)*ww(:,l)))/ray
862  grad1_vel(index,4,i) = (-mode*sum(vel_loc(:,1)*ww(:,l)) - sum(vel_loc(:,4)*ww(:,l)))/ray
863  grad1_vel(index,5,i) = sum(vel_loc(:,1)*dw_loc(2,:))
864  grad1_vel(index,6,i) = sum(vel_loc(:,2)*dw_loc(2,:))
865 
866  !-----------------Grad u_th on Gauss points-----------------------------------
867  grad2_vel(index,1,i) = sum(vel_loc(:,3)*dw_loc(1,:))
868  grad2_vel(index,2,i) = sum(vel_loc(:,4)*dw_loc(1,:))
869  grad2_vel(index,3,i) = (mode*sum(vel_loc(:,4)*ww(:,l)) + sum(vel_loc(:,1)*ww(:,l)))/ray
870  grad2_vel(index,4,i) = (-mode*sum(vel_loc(:,3)*ww(:,l)) + sum(vel_loc(:,2)*ww(:,l)))/ray
871  grad2_vel(index,5,i) = sum(vel_loc(:,3)*dw_loc(2,:))
872  grad2_vel(index,6,i) = sum(vel_loc(:,4)*dw_loc(2,:))
873 
874  !-----------------Grad u_z on Gauss points------------------------------------
875  grad3_vel(index,1,i) = sum(vel_loc(:,5)*dw_loc(1,:))
876  grad3_vel(index,2,i) = sum(vel_loc(:,6)*dw_loc(1,:))
877  grad3_vel(index,3,i) = mode*sum(vel_loc(:,6)*ww(:,l))/ray
878  grad3_vel(index,4,i) = -mode*sum(vel_loc(:,5)*ww(:,l))/ray
879  grad3_vel(index,5,i) = sum(vel_loc(:,5)*dw_loc(2,:))
880  grad3_vel(index,6,i) = sum(vel_loc(:,6)*dw_loc(2,:))
881  ENDDO
882  ENDDO
883  END DO
884 
885  IF (inputs%if_tensor_sym) THEN
886  !-----------------GradT u_r on Gauss points------------------------------------
887  gradt1_vel(:,1,:) = grad1_vel(:,1,:)
888  gradt1_vel(:,2,:) = grad1_vel(:,2,:)
889  gradt1_vel(:,3,:) = grad2_vel(:,1,:)
890  gradt1_vel(:,4,:) = grad2_vel(:,2,:)
891  gradt1_vel(:,5,:) = grad3_vel(:,1,:)
892  gradt1_vel(:,6,:) = grad3_vel(:,2,:)
893  !-----------------GradT u_th on Gauss points-----------------------------------
894  gradt2_vel(:,1,:) = grad1_vel(:,3,:)
895  gradt2_vel(:,2,:) = grad1_vel(:,4,:)
896  gradt2_vel(:,3,:) = grad2_vel(:,3,:)
897  gradt2_vel(:,4,:) = grad2_vel(:,4,:)
898  gradt2_vel(:,5,:) = grad3_vel(:,3,:)
899  gradt2_vel(:,6,:) = grad3_vel(:,4,:)
900  !-----------------GradT u_z on Gauss points------------------------------------
901  gradt3_vel(:,1,:) = grad1_vel(:,5,:)
902  gradt3_vel(:,2,:) = grad1_vel(:,6,:)
903  gradt3_vel(:,3,:) = grad2_vel(:,5,:)
904  gradt3_vel(:,4,:) = grad2_vel(:,6,:)
905  gradt3_vel(:,5,:) = grad3_vel(:,5,:)
906  gradt3_vel(:,6,:) = grad3_vel(:,6,:)
907 
908  !===Grad = Grad + Grad^T
909  grad1_vel = grad1_vel + gradt1_vel
910  grad2_vel = grad2_vel + gradt2_vel
911  grad3_vel = grad3_vel + gradt3_vel
912  END IF
913 
914  !===Compute (visc_dyn + nu_entro)*Grad^s(vel)
915  part_tensor_1 = grad1_vel(:,:,:)
916  part_tensor_2(:,1:4,:) = grad2_vel(:,3:6,:)
917  part_tensor_2(:,5:6,:) = grad3_vel(:,5:6,:)
918 
919  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
920  bloc_size = SIZE(grad1_vel,1)/nb_procs+1
921  CALL fft_compute_diffu_mom(communicator,part_tensor_1, part_tensor_2, visc_dyn_gauss, &
922  part_visc_sym_grad_1, part_visc_sym_grad_2, nb_procs, bloc_size, m_max_pad)
923 
924  !===Using strain rate tensor symmetry to get V1_out= visc_dyn*Grad^s(vel)
925  v_out(1,:,:,:) = part_visc_sym_grad_1
926  v_out(2,:,1:2,:) = part_visc_sym_grad_1(:,3:4,:)
927  v_out(2,:,3:6,:) = part_visc_sym_grad_2(:,1:4,:)
928  v_out(3,:,1:2,:) = part_visc_sym_grad_1(:,5:6,:)
929  v_out(3,:,3:4,:) = part_visc_sym_grad_2(:,3:4,:)
930  v_out(3,:,5:6,:) = part_visc_sym_grad_2(:,5:6,:)
931 
932  END SUBROUTINE smb_explicit_diffu_sym
933 
934  SUBROUTINE smb_explicit_les(communicator, mesh, list_mode, nb_procs, visc_entro_real, mom, V_out)
935  USE gauss_points
936  USE sft_parallele
937  USE chaine_caractere
938  USE boundary
939  USE input_data
940  IMPLICIT NONE
941  TYPE(mesh_type), TARGET :: mesh
942  INTEGER, INTENT(IN) :: nb_procs
943  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
944  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro_real
945  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: mom
946  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: v_out
947  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: grad1_mom, grad2_mom, grad3_mom
948  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
949  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
950  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: mom_loc
951  REAL(KIND=8) :: ray, hh, hm
952  INTEGER :: m, l , i, mode, index, k
953  INTEGER :: m_max_pad, bloc_size
954 #include "petsc/finclude/petsc.h"
955  mpi_comm :: communicator
956 
957  CALL gauss(mesh)
958 
959  DO i = 1, SIZE(list_mode)
960  mode = list_mode(i)
961  index = 0
962  DO m = 1, mesh%dom_me
963  j_loc = mesh%jj(:,m)
964  DO k = 1, 6
965  mom_loc(:,k) = mom(j_loc,k,i)
966  END DO
967  DO l = 1, l_g
968  index = index + 1
969  dw_loc = dw(:,:,l,m)
970 
971  !===Compute local mesh sizes
972  hh=mesh%hloc_gauss(index)
973 !!$ hm=MIN(0.5d0/inputs%m_max,hh) !WRONG choice
974  hm=0.5d0/inputs%m_max
975 
976  !===Compute radius of Gauss point
977  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
978 
979  !-----------------Grad mom_r on Gauss points------------------------------------
980  grad1_mom(index,1,i) = sum(mom_loc(:,1)*dw_loc(1,:))*hh
981  grad1_mom(index,2,i) = sum(mom_loc(:,2)*dw_loc(1,:))*hh
982  grad1_mom(index,3,i) = (mode*sum(mom_loc(:,2)*ww(:,l)) - &
983  sum(mom_loc(:,3)*ww(:,l)))/ray*hm
984  grad1_mom(index,4,i) = (-mode*sum(mom_loc(:,1)*ww(:,l)) - &
985  sum(mom_loc(:,4)*ww(:,l)))/ray*hm
986  grad1_mom(index,5,i) = sum(mom_loc(:,1)*dw_loc(2,:))*hh
987  grad1_mom(index,6,i) = sum(mom_loc(:,2)*dw_loc(2,:))*hh
988 
989  !-----------------Grad mom_th on Gauss points-----------------------------------
990  grad2_mom(index,1,i) = sum(mom_loc(:,3)*dw_loc(1,:))*hh
991  grad2_mom(index,2,i) = sum(mom_loc(:,4)*dw_loc(1,:))*hh
992  grad2_mom(index,3,i) = (mode*sum(mom_loc(:,4)*ww(:,l)) + &
993  sum(mom_loc(:,1)*ww(:,l)))/ray*hm
994  grad2_mom(index,4,i) = (-mode*sum(mom_loc(:,3)*ww(:,l)) + &
995  sum(mom_loc(:,2)*ww(:,l)))/ray*hm
996  grad2_mom(index,5,i) = sum(mom_loc(:,3)*dw_loc(2,:))*hh
997  grad2_mom(index,6,i) = sum(mom_loc(:,4)*dw_loc(2,:))*hh
998 
999  !-----------------Grad mom_z on Gauss points------------------------------------
1000  grad3_mom(index,1,i) = sum(mom_loc(:,5)*dw_loc(1,:))*hh
1001  grad3_mom(index,2,i) = sum(mom_loc(:,6)*dw_loc(1,:))*hh
1002  grad3_mom(index,3,i) = mode*sum(mom_loc(:,6)*ww(:,l))/ray*hm
1003  grad3_mom(index,4,i) = -mode*sum(mom_loc(:,5)*ww(:,l))/ray*hm
1004  grad3_mom(index,5,i) = sum(mom_loc(:,5)*dw_loc(2,:))*hh
1005  grad3_mom(index,6,i) = sum(mom_loc(:,6)*dw_loc(2,:))*hh
1006  ENDDO
1007  ENDDO
1008  END DO
1009 
1010  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1011  bloc_size = SIZE(grad1_mom,1)/nb_procs+1
1012  bloc_size = mesh%gauss%l_G*(bloc_size/mesh%gauss%l_G)+mesh%gauss%l_G
1013  CALL fft_compute_entropy_visc_grad_mom(communicator,grad1_mom, grad2_mom, grad3_mom, visc_entro_real, &
1014  v_out, nb_procs, bloc_size, m_max_pad)
1015 
1016  END SUBROUTINE smb_explicit_les
1017 
1018  SUBROUTINE smb_surface_tension(communicator, mesh, list_mode, nb_procs, level_set, tensor_surface_gauss)
1019  USE gauss_points
1020  USE sft_parallele
1021  USE chaine_caractere
1022  USE boundary
1023  USE input_data
1024  IMPLICIT NONE
1025  TYPE(mesh_type), TARGET :: mesh
1026  INTEGER, INTENT(IN) :: nb_procs
1027  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1028  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: level_set
1029  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: tensor_surface_gauss
1030  REAL(KIND=8), DIMENSION(3,mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: tensor
1031  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: grad_level
1032  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1033  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1034  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: level_set_loc
1035  REAL(KIND=8) :: ray
1036  INTEGER :: m, l , i, mode, index, k, n
1037  INTEGER :: m_max_pad, bloc_size
1038 #include "petsc/finclude/petsc.h"
1039  mpi_comm :: communicator
1040 
1041  CALL gauss(mesh)
1042 
1043  tensor_surface_gauss = 0.d0
1044  DO n = 1, inputs%nb_fluid-1
1045  DO i = 1, SIZE(list_mode)
1046  mode = list_mode(i)
1047  index = 0
1048  DO m = 1, mesh%dom_me
1049  j_loc = mesh%jj(:,m)
1050  DO k = 1, 2
1051  level_set_loc(:,k) = level_set(n,j_loc,k,i)
1052  END DO
1053  DO l = 1, l_g
1054  index = index + 1
1055  dw_loc = dw(:,:,l,m)
1056 
1057  !===Compute radius of Gauss point
1058  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
1059 
1060  !-----------------Grad u_r on Gauss points------------------------------------
1061  grad_level(index,1,i) = sum(level_set_loc(:,1)*dw_loc(1,:))
1062  grad_level(index,2,i) = sum(level_set_loc(:,2)*dw_loc(1,:))
1063  grad_level(index,3,i) = mode/ray*sum(level_set_loc(:,2)*ww(:,l))
1064  grad_level(index,4,i) = -mode/ray*sum(level_set_loc(:,1)*ww(:,l))
1065  grad_level(index,5,i) = sum(level_set_loc(:,1)*dw_loc(2,:))
1066  grad_level(index,6,i) = sum(level_set_loc(:,2)*dw_loc(2,:))
1067  END DO
1068  ENDDO
1069  ENDDO
1070  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1071  bloc_size = SIZE(grad_level,1)/nb_procs+1
1072  CALL fft_tensor_dcl(communicator, grad_level, grad_level, tensor, nb_procs, bloc_size, m_max_pad, opt_tension=.true.)
1073  tensor_surface_gauss = tensor_surface_gauss + inputs%coeff_surface(n)*tensor
1074  END DO
1075 
1076  END SUBROUTINE smb_surface_tension
1077 
1078  SUBROUTINE momentum_dirichlet(communicator, mesh, list_mode, t, nb_procs, density, momentum_exact, vv_js_D)
1079  USE gauss_points
1080  USE sft_parallele
1081  USE chaine_caractere
1082  USE boundary
1083  USE input_data
1084  IMPLICIT NONE
1085  TYPE(mesh_type), TARGET :: mesh
1086  INTEGER, INTENT(IN) :: nb_procs
1087  REAL(KIND=8), INTENT(IN) :: t
1088  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1089  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: density
1090  TYPE(dyn_int_line), DIMENSION(3),INTENT(IN) :: vv_js_d
1091  REAL(KIND=8), DIMENSION(mesh%np,6,SIZE(list_mode)), INTENT(OUT) :: momentum_exact
1092  REAL(KIND=8), DIMENSION(SIZE(vv_js_D(1)%DIL),2,SIZE(list_mode)) :: vel_exact_r, mr
1093  REAL(KIND=8), DIMENSION(SIZE(vv_js_D(2)%DIL),2,SIZE(list_mode)) :: vel_exact_t, mt
1094  REAL(KIND=8), DIMENSION(SIZE(vv_js_D(3)%DIL),2,SIZE(list_mode)) :: vel_exact_z, mz
1095  INTEGER :: i, k, kk
1096  INTEGER :: m_max_pad, bloc_size
1097 #include "petsc/finclude/petsc.h"
1098  mpi_comm :: communicator
1099 
1100  IF (inputs%if_level_set) THEN
1101  DO i = 1, SIZE(list_mode)
1102  DO k = 1, 2
1103  vel_exact_r(:,k,i) = vv_exact(k, mesh%rr(:,vv_js_d(1)%DIL),list_mode(i),t)
1104  vel_exact_t(:,k,i) = vv_exact(k+2,mesh%rr(:,vv_js_d(2)%DIL),list_mode(i),t)
1105  vel_exact_z(:,k,i) = vv_exact(k+4,mesh%rr(:,vv_js_d(3)%DIL),list_mode(i),t)
1106  END DO
1107  END DO
1108 
1109  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1110  bloc_size = SIZE(vv_js_d(1)%DIL)/nb_procs+1
1111  CALL fft_par_prod_dcl(communicator, density(vv_js_d(1)%DIL,:,:), vel_exact_r, &
1112  mr, nb_procs, bloc_size, m_max_pad)
1113  momentum_exact(vv_js_d(1)%DIL,1:2,:) = mr
1114 
1115  bloc_size = SIZE(vv_js_d(2)%DIL)/nb_procs+1
1116  CALL fft_par_prod_dcl(communicator, density(vv_js_d(2)%DIL,:,:), vel_exact_t, &
1117  mt, nb_procs, bloc_size, m_max_pad)
1118  momentum_exact(vv_js_d(2)%DIL,3:4,:) = mt
1119 
1120  bloc_size = SIZE(vv_js_d(3)%DIL)/nb_procs+1
1121  CALL fft_par_prod_dcl(communicator, density(vv_js_d(3)%DIL,:,:), vel_exact_z, &
1122  mz, nb_procs, bloc_size, m_max_pad)
1123  momentum_exact(vv_js_d(3)%DIL,5:6,:) = mz
1124  ELSE
1125  DO i = 1, SIZE(list_mode)
1126  DO k = 1, 6
1127  kk = (k-1)/2 + 1
1128  momentum_exact(vv_js_d(kk)%DIL,k,i) = vv_exact(k,mesh%rr(:,vv_js_d(kk)%DIL),list_mode(i),t)
1129  END DO
1130  END DO
1131  END IF
1132 
1133  END SUBROUTINE momentum_dirichlet
1134 
1135  SUBROUTINE twod_volume(communicator,mesh,RESLT)
1136  !===========================
1137  USE def_type_mesh
1138  IMPLICIT NONE
1139  TYPE(mesh_type) :: mesh
1140  REAL(KIND=8), INTENT(OUT) :: reslt
1141  REAL(KIND=8) :: vol_loc, vol_out
1142  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1143  INTEGER :: m, l , i , ni, code
1144  REAL(KIND=8) :: ray
1145 #include "petsc/finclude/petsc.h"
1146  mpi_comm :: communicator
1147  vol_loc = 0.d0
1148  DO m = 1, mesh%dom_me
1149  j_loc = mesh%jj(:,m)
1150  DO l = 1, mesh%gauss%l_G
1151  ray = 0
1152  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
1153  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1154  END DO
1155  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
1156  ENDDO
1157  ENDDO
1158  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
1159  reslt = vol_out
1160  END SUBROUTINE twod_volume
1161 
1162  SUBROUTINE moy(communicator,mesh,p,RESLT)
1163  !===========================
1164  !===average of p
1165  USE def_type_mesh
1166  IMPLICIT NONE
1167  TYPE(mesh_type) :: mesh
1168  REAL(KIND=8), DIMENSION(:) , INTENT(IN) :: p
1169  REAL(KIND=8) , INTENT(OUT) :: reslt
1170  REAL(KIND=8) :: vol_loc, vol_out, r_loc, r_out
1171  INTEGER :: m, l , i , ni, code
1172  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1173  REAL(KIND=8) :: ray
1174 #include "petsc/finclude/petsc.h"
1175  mpi_comm :: communicator
1176  r_loc = 0.d0
1177  vol_loc = 0.d0
1178 
1179  DO m = 1, mesh%dom_me
1180  j_loc = mesh%jj(:,m)
1181  DO l = 1, mesh%gauss%l_G
1182  ray = 0
1183  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
1184  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1185  END DO
1186 
1187  r_loc = r_loc + sum(p(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
1188  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
1189 
1190  ENDDO
1191  ENDDO
1192 
1193  CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
1194  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
1195  reslt = r_out / vol_out
1196 
1197  END SUBROUTINE moy
1198 
1199 END MODULE subroutine_ns_with_m
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic t
Definition: doc_intro.h:199
Definition: tn_axi.f90:5
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
subroutine smb_curlh_cross_b_gauss_sft_par(communicator, mesh, list_mode, V_in, W_in, V_out)
subroutine smb_surface_tension(communicator, mesh, list_mode, nb_procs, level_set, tensor_surface_gauss)
subroutine smb_explicit_les(communicator, mesh, list_mode, nb_procs, visc_entro_real, mom, V_out)
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
subroutine dirichlet_rhs(js_D, bs_D, b)
subroutine, public fft_par_prod_dcl(communicator, c1_in, c2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine vector_glob_js_d(vv_mesh, list_mode, vv_3_LA, vv_list_dirichlet_sides, vv_js_D, vv_mode_global_js_D)
subroutine scalar_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_mode_global_js_D)
subroutine, public inject_p1_p2(jj_c, jj_f, pp_c, pp_f)
Definition: sub_mass.f90:242
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
subroutine qs_ns_momentum_stab_3x3(mesh, vv_3_LA, mode, ff, V1m, P, vb_3_145, vb_3_236, rotb_b, tensor, tensor_surface_gauss, stab, momentum, visc_grad_vel)
real(kind=8) function, dimension(size(rr, 2), 6), public h_b_quasi_static(char_h_b, rr, m)
Definition: condlim.f90:276
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:127
real(kind=8) function norm_s(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:85
subroutine, public fft_compute_entropy_visc_grad_mom(communicator, V1_in, V2_in, V3_in, c_in_real, V_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
subroutine dirichlet_m_parallel(matrix, glob_js_D)
subroutine moy(mesh, p, RESULT)
subroutine inject(jj_c, jj_f, pp_c, pp_f)
subroutine, public three_level_ns_tensor_sym_with_m(comm_one_d, time, vv_3_LA, pp_1_LA, dt, Re, list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, Hn_p2, Bn_p2, density_m1, density, density_p1, visco_dyn, tempn, level_set, level_set_p1, visc_entro_level)
subroutine gauss(mesh)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
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 qs_01_div_hybrid_2006(uu_mesh, pp_mesh, mode, gg, u0_c)
subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
Definition: fem_M_axi.f90:356
real(kind=8) function user_time()
Definition: my_util.f90:7
subroutine momentum_dirichlet(communicator, mesh, list_mode, t, nb_procs, density, momentum_exact, vv_js_D)
subroutine, public fft_compute_diffu_mom(communicator, V1_in, V2_in, V3_in, V1_out, V2_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public fft_tensor_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps, opt_tension)
subroutine smb_explicit_diffu_sym(communicator, mesh, list_mode, nb_procs, visc_dyn, vel, V_out)
subroutine error_petsc(string)
Definition: my_util.f90:15
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
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)