SFEMaNS  version 5.3
Reference documentation for SFEMaNS
sub_taylor_navier_stokes.f90
Go to the documentation of this file.
1 !Authors Jean-Luc Guermond, Hugues Faller, 2019
2 !cf https://doi.org/10.1137/18M1209301
4  USE my_util
5  USE input_data
6 #include "petsc/finclude/petsc.h"
7  USE petsc
9  PRIVATE
10 
11 CONTAINS
12 
13  SUBROUTINE navier_stokes_taylor(comm_one_d_ns,time, vv_3_LA, pp_1_LA, &
14  list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un, vvz_per, pp_per)
16  USE periodic
17 
18  USE my_util
19  IMPLICIT NONE
20  TYPE(mesh_type), INTENT(IN) :: vv_mesh, pp_mesh
21  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
22  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: un
23  REAL(KIND=8), DIMENSION(pp_mesh%np,2,SIZE(list_mode)), INTENT(INOUT) :: pn
24  TYPE(dyn_real_array_three), DIMENSION(:) :: der_un
25  TYPE(dyn_real_array_three), DIMENSION(:) :: der_pn
26  REAL(KIND=8), INTENT(IN) :: time
27  TYPE(periodic_type), INTENT(IN) :: vvz_per, pp_per
28  TYPE(petsc_csr_la) :: vv_3_LA, pp_1_LA
29  LOGICAL :: if_navier_stokes_with_taylor = .true. !HF We will nead a 'read_data'
30 #include "petsc/finclude/petsc.h"
31  mpi_comm, DIMENSION(:), POINTER :: comm_one_d_ns
32 
33  IF (if_navier_stokes_with_taylor) THEN
34  IF (inputs%taylor_order==3) THEN
35  CALL update_ns_with_taylor(comm_one_d_ns,time,vv_3_la, pp_1_la, vvz_per, pp_per, &
36  inputs%dt, inputs%Re, inputs%taylor_lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
37  ELSE IF (inputs%taylor_order==4) THEN
38  CALL update_ns_with_taylor_fourth(comm_one_d_ns,time,vv_3_la, pp_1_la, vvz_per, pp_per, &
39  inputs%dt, inputs%Re, inputs%taylor_lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
40  ELSE
41  CALL error_petsc('BUG in navier_stokes_taylor: inputs%taylor_order not well defined')
42  END IF
43  END IF
44  END SUBROUTINE navier_stokes_taylor
45 
46  SUBROUTINE init_velocity_pressure_taylor(vv_mesh, pp_mesh, list_mode, time, pn, der_pn, un, der_un)
49  IMPLICIT NONE
50  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
51  TYPE(mesh_type), INTENT(IN) :: vv_mesh, pp_mesh
52  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: un
53  REAL(KIND=8), DIMENSION(pp_mesh%np,2,SIZE(list_mode)), INTENT(INOUT) :: pn
54  TYPE(dyn_real_array_three), DIMENSION(inputs%taylor_order-1) :: der_un
55  TYPE(dyn_real_array_three), DIMENSION(inputs%taylor_order-1) :: der_pn
56  REAL(KIND=8), INTENT(IN) :: time
57  REAL(KIND=8) :: dt
58  INTEGER :: i, mode, k, kp, m_max_c
59 
60  dt = inputs%dt
61  m_max_c = SIZE(list_mode)
62  DO i= 1, m_max_c
63  mode = list_mode(i)
64  DO k = 1, 6
65  !===initialize velocity at time=t0
66  un(:,k,i) = vv_exact(k,vv_mesh%rr, mode,time)
67  der_un(1)%DRT(:,k,i)= &
68  (-2*vv_exact(k,vv_mesh%rr, mode,time-3*dt)&
69  +9*vv_exact(k,vv_mesh%rr, mode,time-2*dt)&
70  -18*vv_exact(k,vv_mesh%rr, mode,time-dt)&
71  +11*vv_exact(k,vv_mesh%rr, mode,time))/(6*dt)
72  der_un(2)%DRT(:,k,i)= &
73  (-1*vv_exact(k,vv_mesh%rr, mode,time-3*dt)&
74  +4*vv_exact(k,vv_mesh%rr, mode,time-2*dt)&
75  -5*vv_exact(k,vv_mesh%rr, mode,time-dt)&
76  +2*vv_exact(k,vv_mesh%rr, mode,time))/(dt**2)
77  END DO
78  DO k = 1, 2
79  !===pressure at time=t0
80  pn(:,k,i) = pp_exact(k,pp_mesh%rr, mode,time)
81  der_pn(1)%DRT(:,k,i)= &
82  (-2*pp_exact(k,pp_mesh%rr, mode,time-3*dt)&
83  +9*pp_exact(k,pp_mesh%rr, mode,time-2*dt)&
84  -18*pp_exact(k,pp_mesh%rr, mode,time-dt)&
85  +11*pp_exact(k,pp_mesh%rr, mode,time))/(6*dt)
86  der_pn(2)%DRT(:,k,i)= &
87  (-1*pp_exact(k,pp_mesh%rr, mode,time-3*dt)&
88  +4*pp_exact(k,pp_mesh%rr, mode,time-2*dt)&
89  -5*pp_exact(k,pp_mesh%rr, mode,time-dt)&
90  +2*pp_exact(k,pp_mesh%rr, mode,time))/(dt**2)
91  ENDDO
92  IF (inputs%taylor_order==4) THEN
93  DO k = 1, 6
94  der_un(3)%DRT(:,k,i)= &
95  (-1*vv_exact(k,vv_mesh%rr, mode,time-3*dt)&
96  +3*vv_exact(k,vv_mesh%rr, mode,time-2*dt)&
97  -3*vv_exact(k,vv_mesh%rr, mode,time-dt)&
98  +1*vv_exact(k,vv_mesh%rr, mode,time))/(dt**3)
99  END DO
100  DO k = 1, 2
101  der_pn(3)%DRT(:,k,i)= &
102  (-1*pp_exact(k,pp_mesh%rr, mode,time-3*dt)&
103  +3*pp_exact(k,pp_mesh%rr, mode,time-2*dt)&
104  -3*pp_exact(k,pp_mesh%rr, mode,time-dt)&
105  +1*pp_exact(k,pp_mesh%rr, mode,time))/(dt**3)
106  END DO
107  END IF
108  ENDDO
109 
110  !===Set to zero
111  IF (vv_mesh%me/=0) THEN
112  DO kp = 1, inputs%taylor_order-1
113  DO i = 1, m_max_c
114  IF (list_mode(i) == 0) THEN
115  DO k= 1, 3!===velocity
116  der_un(kp)%DRT(:,2*k,i) = 0.d0
117  END DO
118  der_pn(kp)%DRT(:,2,i) = 0.d0!===pressure
119  END IF
120  END DO
121  END DO
122  END IF
123  END SUBROUTINE init_velocity_pressure_taylor
124 
125  SUBROUTINE update_ns_with_taylor(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, &
126  dt, re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
127  !==============================
128  USE def_type_mesh
129  USE fem_m_axi
130  USE fem_rhs_axi
131  USE dir_nodes_petsc
132  USE periodic
133  USE st_matrix
134  USE solve_petsc
135  USE boundary
136  USE st_matrix
137  USE input_data
139  USE tn_axi
140  USE verbose
141  USE my_util
142  IMPLICIT NONE
143  REAL(KIND=8) :: time, dt, Re, lambda
144  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
145  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
146  TYPE(petsc_csr_la) :: vv_3_LA, pp_1_LA
147  TYPE(periodic_type), INTENT(IN) :: vvz_per, pp_per
148  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: un
149  REAL(KIND=8), DIMENSION(pp_mesh%np,2,SIZE(list_mode)), INTENT(INOUT) :: pn
150  TYPE(dyn_real_array_three), DIMENSION(:) :: der_un
151  TYPE(dyn_real_array_three), DIMENSION(:) :: der_pn
152  !===Saved variables
153  INTEGER, SAVE :: m_max_c
154  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: pp_global_D
155  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: pp_mode_global_js_D
156  TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_D
157  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: vv_mode_global_js_D
158  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: vel_global_D
159  LOGICAL, SAVE :: once = .true.
160  INTEGER, SAVE :: my_petscworld_rank
161  !===End saved variables
162  INTEGER, POINTER, DIMENSION(:) :: pp_1_ifrom, vv_3_ifrom
163  INTEGER :: i, m, n, n1, n2, n3, n123
164  INTEGER :: nb_procs, code, nu_mat, mode
165  REAL(KIND=8) :: moyenne
166  REAL(KIND=8), DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: div
167  REAL(KIND=8), DIMENSION(vv_mesh%np, 6) :: un_p1
168  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rotv_v1, rotv_v2, rotv_v3
169  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
170  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: uext, dtt_un_p1, dt_un_p1
171  REAL(KIND=8), DIMENSION(vv_mesh%np) :: vel_loc, vel_tot
172  REAL(KIND=8) :: tps, tps_tot, tps_cumul, coeff, vloc, cfl, cfl_max, norm
173  REAL(KIND=8) :: one, zero, three
174  DATA zero, one, three/0.d0,1.d0,3.d0/
175  !===Communicators for Petsc, in space and Fourier
176  petscerrorcode :: ierr
177  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
178  mat, DIMENSION(:), POINTER, SAVE :: vel_mat
179  mat, SAVE :: mass_mat
180  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
181  vec, SAVE :: pb_1, pb_2, px_1, px_1_ghost
182  ksp, DIMENSION(:), POINTER, SAVE :: vel_ksp
183  ksp, SAVE :: mass_ksp
184  !===END OF DECLARATIONS
185 
186  IF (once) THEN
187  once = .false.
188 
189  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
190  !===CREATE PETSC VECTORS AND GHOSTS
191  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
192  n = 3*vv_mesh%dom_np
193  CALL veccreateghost(comm_one_d(1), n, &
194  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
195  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
196  CALL vecduplicate(vx_3, vb_3_145, ierr)
197  CALL vecduplicate(vx_3, vb_3_236, ierr)
198 
199  CALL create_my_ghost(pp_mesh,pp_1_la,pp_1_ifrom)
200  n = pp_mesh%dom_np
201  CALL veccreateghost(comm_one_d(1), n, &
202  petsc_determine, SIZE(pp_1_ifrom), pp_1_ifrom, px_1, ierr)
203  CALL vecghostgetlocalform(px_1, px_1_ghost, ierr)
204  CALL vecduplicate(px_1, pb_1, ierr)
205  CALL vecduplicate(px_1, pb_2, ierr)
206  !===End CREATE PETSC VECTORS AND GHOSTS
207 
208  !===Number of modes on proc
209  m_max_c = SIZE(list_mode)
210  !===End umber of modes on proc
211 
212 
213  !===PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
214  CALL scalar_without_glob_js_d(pp_mesh, list_mode, pp_1_la, pp_mode_global_js_d)
215  ALLOCATE(pp_global_d(m_max_c))
216  DO i = 1, m_max_c
217  ALLOCATE(pp_global_d(i)%DRL(SIZE(pp_mode_global_js_d(i)%DIL)))
218  END DO
219  !===End PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
220 
221  !===PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
222  CALL vector_glob_js_d(vv_mesh, list_mode, vv_3_la, inputs%vv_list_dirichlet_sides, &
223  vv_js_d, vv_mode_global_js_d)
224 
225  ALLOCATE(vel_global_d(m_max_c))
226  DO i = 1, m_max_c
227  ALLOCATE(vel_global_d(i)%DRL(SIZE(vv_mode_global_js_d(i)%DIL)))
228  END DO
229  !===END PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
230 
231  !===ASSEMBLE MASS MATRIX
232  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat, clean=.false.)
233  CALL qs_diff_mass_scal_m (pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat)
234  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
235  CALL periodic_matrix_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, mass_mat, pp_1_la)
236  END IF
237  DO i = 1, m_max_c
238  CALL dirichlet_m_parallel(mass_mat,pp_mode_global_js_d(i)%DIL)
239  END DO
240  CALL init_solver(inputs%my_par_mass,mass_ksp,mass_mat,comm_one_d(1),&
241  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
242  !===END ASSEMBLE MASS MATRIX
243 
244  !===ASSEMBLING VELOCITY MATRICES
245  ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
246  DO i = 1, m_max_c
247  mode = list_mode(i)
248  !===VELOCITY
249  nu_mat = 2*i-1
250  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
251  CALL qs_diff_mass_vect_3x3_taylor (1, vv_3_la, vv_mesh, one/re, one/dt, & !HF new coeffs
252  lambda, i, mode, vel_mat(nu_mat))
253 
254  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
255  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
256  vel_mat(nu_mat), vv_3_la)
257  END IF
258  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
259  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
260  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
261  nu_mat = nu_mat+1
262  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
263  CALL qs_diff_mass_vect_3x3_taylor (2, vv_3_la, vv_mesh, one/re, one/dt, & !HF new coeffs
264  lambda, i, mode, vel_mat(nu_mat))
265 
266  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
267  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
268  vel_mat(nu_mat), vv_3_la)
269  END IF
270  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
271  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
272  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
273  !===End VELOCITY
274  ENDDO
275  !===End ASSEMBLING VELOCITY MATRICES
276  ENDIF !===End of once
277 
278  !===Compute NL for dtt_un by FFT at Gauss points
279  tps_tot = user_time()
280  tps_cumul = 0
281  tps = user_time()
282 
283  uext = un + 2.d0*dt*der_un(1)%DRT + 2.d0*dt**2 *der_un(2)%DRT !HF 1st NL term
284  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v1,inputs%precession)
285 
286  uext = un + dt*der_un(1)%DRT + 0.5d0*dt**2*der_un(2)%DRT !HF 2nd NL term
287  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v2,inputs%precession)
288 
289  uext = un
290  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v3,inputs%precession)
291  tps = user_time() - tps; tps_cumul=tps_cumul+tps
292  !===End Compute NL by FFT at Gauss points
293 
294  !===Computation of CFL
295  IF (inputs%verbose_CFL) THEN
296  vel_loc = 0.d0
297  DO i = 1, m_max_c
298  IF (list_mode(i)==0) THEN
299  coeff = 1.d0
300  ELSE
301  coeff = .5d0
302  END IF
303  vel_loc = vel_loc + coeff*(un(:,1,i)**2+un(:,2,i)**2+un(:,3,i)**2 &
304  +un(:,4,i)**2+un(:,5,i)**2+un(:,6,i)**2)
305  END DO
306  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
307  CALL mpi_allreduce(vel_loc,vel_tot,vv_mesh%np,mpi_double_precision, mpi_sum, comm_one_d(2), code)
308  vel_tot = sqrt(vel_tot)
309  cfl = 0.d0
310  DO m = 1, vv_mesh%dom_me
311  vloc = maxval(vel_tot(vv_mesh%jj(:,m)))
312  cfl = max(vloc*dt/min(vv_mesh%hloc(m),maxval(vv_mesh%hm)),cfl)
313  END DO
314  CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
315  talk_to_me%CFL=cfl_max
316  talk_to_me%time=time
317  END IF
318  !===End Computation of CFL
319 
320  !===Computation of second order derivatives
321  !===Computation of rhs for dtt_un at Gauss points for every mode
322  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
323  der_un(2)%DRT/dt, der_pn(2)%DRT, (rotv_v1-2.d0*rotv_v2+rotv_v3)/(dt**2), rhs_gauss,2)
324  !===End Computation of rhs for dtt_un
325 
326  DO i = 1, m_max_c
327  tps = user_time()
328  mode = list_mode(i)
329 
330  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
331  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
332  !===Periodic boundary conditions
333  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
334  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
335  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
336  END IF
337  !===Axis boundary conditions
338  n1 = SIZE(vv_js_d(1)%DIL)
339  n2 = SIZE(vv_js_d(2)%DIL)
340  n3 = SIZE(vv_js_d(3)%DIL)
341  n123 = n1+n2+n3
342  vel_global_d(i)%DRL(1:n1) =(vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time+dt)&
343  -2.d0*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
344  +vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt))/(dt**2)
345  vel_global_d(i)%DRL(n1+1:n1+n2) =(vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time+dt)&
346  -2.d0*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
347  +vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt))/(dt**2)
348  vel_global_d(i)%DRL(n1+n2+1:n123)=(vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time+dt)&
349  -2.d0*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
350  +vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt))/(dt**2)
351  vel_global_d(i)%DRL(n123+1:) = 0.d0
352  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
353 
354  vel_global_d(i)%DRL(1:n1) =(vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time+dt)&
355  -2.d0*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
356  +vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt))/(dt**2)
357  vel_global_d(i)%DRL(n1+1:n1+n2) =(vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time+dt)&
358  -2.d0*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
359  +vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt))/(dt**2)
360  vel_global_d(i)%DRL(n1+n2+1:n123)=(vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time+dt)&
361  -2.d0*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
362  +vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt))/(dt**2)
363  vel_global_d(i)%DRL(n123+1:) = 0.d0
364  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
365  tps = user_time() - tps; tps_cumul=tps_cumul+tps
366  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
367 
368  !===Solve linear system for momentum equation
369  tps = user_time()
370  !===Solve system 1, dtt_ur_c, dtt_ut_s, dtt_uz_c
371  nu_mat =2*i-1
372  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
373  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
374  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
375  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
376  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
377  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
378 
379  !===Solve system 2, dtt_ur_s, dtt_ut_c, dtt_uz_s
380  nu_mat = nu_mat + 1
381  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
382  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
383  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
384  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
385  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
386  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
387  tps = user_time() - tps; tps_cumul=tps_cumul+tps
388  !===End Solve linear system for momentum equation
389 
390  !===Assemble divergence of dtt_velocity in arrays pb_1, pb_2
391  tps = user_time()
392  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
393  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
394 
395  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
396  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
397  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
398  END IF
399 
400  !===Boundary condition on axis for dtt_pressure
401  pp_global_d(i)%DRL = 0.d0
402  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
403  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
404  !===End boundary condition on axis for dtt_pressure
405 
406  !===Solve mass matrix for dtt_pressure correction
407  tps = user_time()
408  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
409  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
410  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
411  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
412  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
413  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
414  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
415  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
416  !===End solve mass matrix for dtt_pressure correction
417 
418  !===Compute 2nd order derrivative pressure for this mode
419  der_pn(2)%DRT(:,:,i) = der_pn(2)%DRT(:,:,i) - lambda*div(:,:,i)
420  tps = user_time() - tps; tps_cumul=tps_cumul+tps
421  !===End 2nd order derrivative Pressure computation we will then take mean out
422 
423  !===UPDATES
424  tps = user_time()
425  !===Handling 2nd order derrivative of mean pressure
426  IF (mode == 0) THEN
427  CALL moy(comm_one_d(1),pp_mesh, der_pn(2)%DRT(:,1,i),moyenne)
428  der_pn(2)%DRT(:,1,i) = der_pn(2)%DRT(:,1,i)-moyenne
429  ENDIF
430  !===End of handling of 2nd order derrivative mean pressure
431 
432  !===Correction of zero mode
433  IF (mode==0) THEN
434  un_p1(:,2) = 0.d0
435  un_p1(:,4) = 0.d0
436  un_p1(:,6) = 0.d0
437  der_pn(2)%DRT (:,2,i) = 0.d0
438  END IF
439  !===Correction of zero mode
440 
441  !===UPDATES 2nd order derrivative
442  dtt_un_p1(:,:,i) = un_p1
443  tps = user_time() - tps; tps_cumul=tps_cumul+tps
444  !===End UPDATES 2nd order derrivative
445 
446  ENDDO
447  !===END Computation of second order derivatives
448 
449  !===Computation of first order derivatives
450  !===Computation of rhs for dt_un at Gauss points for every mode
451  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
452  (der_un(1)%DRT-dt/2.d0*(dtt_un_p1-der_un(2)%DRT))/dt, der_pn(1)%DRT+ dt*der_pn(2)%DRT, &
453  (rotv_v1-rotv_v3)/(2.d0*dt), rhs_gauss,1)
454  !===End Computation of rhs for dt_un
455 
456  DO i = 1, m_max_c
457  tps = user_time()
458  mode = list_mode(i)
459 
460  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
461  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
462  !===Periodic boundary conditions
463  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
464  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
465  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
466  END IF
467  !===Axis boundary conditions
468  n1 = SIZE(vv_js_d(1)%DIL)
469  n2 = SIZE(vv_js_d(2)%DIL)
470  n3 = SIZE(vv_js_d(3)%DIL)
471  n123 = n1+n2+n3
472  vel_global_d(i)%DRL(1:n1) =(vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time+dt)&
473  -vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt))/(2.d0*dt)
474  vel_global_d(i)%DRL(n1+1:n1+n2) =(vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time+dt)&
475  -vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt))/(2.d0*dt)
476  vel_global_d(i)%DRL(n1+n2+1:n123)=(vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time+dt)&
477  -vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt))/(2.d0*dt)
478  vel_global_d(i)%DRL(n123+1:) = 0.d0
479  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
480  vel_global_d(i)%DRL(1:n1) =(vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time+dt)&
481  -vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt))/(2.d0*dt)
482  vel_global_d(i)%DRL(n1+1:n1+n2) =(vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time+dt)&
483  -vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt))/(2.d0*dt)
484  vel_global_d(i)%DRL(n1+n2+1:n123)=(vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time+dt)&
485  -vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt))/(2.d0*dt)
486  vel_global_d(i)%DRL(n123+1:) = 0.d0
487  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
488  tps = user_time() - tps; tps_cumul=tps_cumul+tps
489  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
490 
491  !===Solve linear system for momentum equation
492  tps = user_time()
493  !===Solve system 1, dt_ur_c, dt_ut_s, dt_uz_c
494  nu_mat =2*i-1
495  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
496  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
497  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
498  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
499  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
500  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
501 
502  !===Solve system 2, dt_ur_s, dt_ut_c, dt_uz_s
503  nu_mat = nu_mat + 1
504  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
505  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
506  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
507  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
508  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
509  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
510  tps = user_time() - tps; tps_cumul=tps_cumul+tps
511  !===End Solve linear system for momentum equation
512 
513  !===Assemble divergence of dt_velocity in arrays pb_1, pb_2
514  tps = user_time()
515  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
516  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
517 
518  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
519  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
520  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
521  END IF
522 
523  !===Boundary condition on axis for dt_pressure
524  pp_global_d(i)%DRL = 0.d0
525  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
526  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
527  !===End boundary condition on axis for dt_pressure
528 
529  !===Solve mass matrix for dt_pressure correction
530  tps = user_time()
531  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
532  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
533  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
534  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
535  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
536  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
537  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
538  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
539  !===End solve mass matrix for dt_pressure correction
540 
541  !===Compute first order derrivative pressure for this mode
542  !===Pressure first order derrivative computation
543  der_pn(1)%DRT(:,:,i) = der_pn(1)%DRT(:,:,i) +dt*der_pn(2)%DRT(:,:,i)- lambda*div(:,:,i)
544  tps = user_time() - tps; tps_cumul=tps_cumul+tps
545  !===End dt_Pressure computation
546 
547  !===UPDATES
548  tps = user_time()
549  !===Handling of mean first order derrivative pressure
550  IF (mode == 0) THEN
551  CALL moy(comm_one_d(1),pp_mesh, der_pn(1)%DRT(:,1,i),moyenne)
552  der_pn(1)%DRT(:,1,i) = der_pn(1)%DRT(:,1,i)-moyenne
553  ENDIF
554  !===End of handling of mean first order derrivative pressure
555 
556  !===Correction of zero mode
557  IF (mode==0) THEN
558  un_p1(:,2) = 0.d0
559  un_p1(:,4) = 0.d0
560  un_p1(:,6) = 0.d0
561  der_pn(1)%DRT (:,2,i) = 0.d0
562  END IF
563  !===Correction of zero mode
564 
565  !===UPDATE Velocity first order derrivative
566  dt_un_p1(:,:,i) = un_p1
567  tps = user_time() - tps; tps_cumul=tps_cumul+tps
568  !===End UPDATE Velocity first order derrivative
569 
570  ENDDO ! i = 1, m_max_c
571  !===END Computation of first order derivatives
572 
573  !===Computation of 0th order derivatives
574  !===Computation of rhs for un at Gauss points for every mode
575  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
576  (un-dt/2.d0*(dt_un_p1-der_un(1)%DRT)-dt**2/12.d0*(dtt_un_p1-der_un(2)%DRT))/dt, &
577  pn+dt*der_pn(1)%DRT- (dt**2)/2.d0*der_pn(2)%DRT, rotv_v2, rhs_gauss,0)
578  !===End Computation of rhs for un
579 
580  DO i = 1, m_max_c
581  tps = user_time()
582  mode = list_mode(i)
583 
584  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
585  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
586  !===Periodic boundary conditions
587  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
588  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
589  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
590  END IF
591  !===Axis boundary conditions
592  n1 = SIZE(vv_js_d(1)%DIL)
593  n2 = SIZE(vv_js_d(2)%DIL)
594  n3 = SIZE(vv_js_d(3)%DIL)
595  n123 = n1+n2+n3
596  vel_global_d(i)%DRL(1:n1) = vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
597  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
598  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
599  vel_global_d(i)%DRL(n123+1:) = 0.d0
600  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
601  vel_global_d(i)%DRL(1:n1) = vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
602  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
603  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
604  vel_global_d(i)%DRL(n123+1:) = 0.d0
605  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
606  tps = user_time() - tps; tps_cumul=tps_cumul+tps
607  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
608 
609  !===Solve linear system for momentum equation
610  tps = user_time()
611  !===Solve system 1, ur_c, ut_s, uz_c
612  nu_mat =2*i-1
613  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
614  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
615  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
616  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
617  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
618  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
619 
620  !===Solve system 2, ur_s, ut_c, uz_s
621  nu_mat = nu_mat + 1
622  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
623  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
624  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
625  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
626  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
627  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
628 
629  tps = user_time() - tps; tps_cumul=tps_cumul+tps
630  !===End Solve linear system for momentum equation
631 
632  !===Assemble divergence of velocity in arrays pb_1, pb_2
633  tps = user_time()
634  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
635  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
636 
637  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
638  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
639  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
640  END IF
641 
642  !===Boundary condition on axis for pressure
643  pp_global_d(i)%DRL = 0.d0
644  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
645  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
646  !===End boundary condition on axis for pressure
647 
648  !===Solve mass matrix for pressure correction
649  tps = user_time()
650  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
651  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
652  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
653  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
654  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
655  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
656  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
657  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
658  !===End solve mass matrix for pressure correction
659 
660  !===Pressure computation
661  pn(:,:,i) = pn(:,:,i) +dt*der_pn(1)%DRT(:,:,i)-(dt**2)/2.d0*der_pn(2)%DRT(:,:,i)- lambda*div(:,:,i)
662  tps = user_time() - tps; tps_cumul=tps_cumul+tps
663  !===End Pressure computation
664 
665  !===UPDATES
666  tps = user_time()
667  !===Handling of mean pressure
668  IF (mode == 0) THEN
669  CALL moy(comm_one_d(1),pp_mesh, pn(:,1,i),moyenne)
670  pn(:,1,i) = pn(:,1,i)-moyenne
671  ENDIF
672  !===End of handling of mean pressure
673 
674  !===Correction of zero mode
675  IF (mode==0) THEN
676  un_p1(:,2) = 0.d0
677  un_p1(:,4) = 0.d0
678  un_p1(:,6) = 0.d0
679  pn(:,2,i) = 0.d0
680  END IF
681  !===Correction of zero mode
682 
683  !===UPDATES Velocity 0th order derrivative
684  un(:,:,i) = un_p1
685  tps = user_time() - tps; tps_cumul=tps_cumul+tps
686  !===End UPDATES Velocity 0th order derrivative
687 
688  ENDDO ! i = 1, m_max_c
689  !===END Computation of 0th order derivatives
690 
691  !===Update dtt_un un and dt_un once everyone has been compuded
692  DO i = 1, m_max_c
693  der_un(2)%DRT(:,:,i)=dtt_un_p1(:,:,i)
694  der_un(1)%DRT(:,:,i)=dt_un_p1(:,:,i)
695  END DO
696 
697  IF (inputs%verbose_divergence) THEN
698  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
699  talk_to_me%div_L2 = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)/norm
700  talk_to_me%weak_div_L2 = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, div)/norm
701  END IF
702 
703  tps_tot = user_time() - tps_tot
704  END SUBROUTINE update_ns_with_taylor
705  !============================================
706 
707  SUBROUTINE update_ns_with_taylor_fourth(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, &
708  dt, re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
709  !==============================
710  USE def_type_mesh
711  USE fem_m_axi
712  USE fem_rhs_axi
713  USE dir_nodes_petsc
714  USE periodic
715  USE st_matrix
716  USE solve_petsc
717  USE boundary
718  USE st_matrix
719  USE input_data
721  USE tn_axi
722  USE verbose
723  USE my_util
724  IMPLICIT NONE
725  REAL(KIND=8) :: time, dt, Re, lambda
726  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
727  TYPE(mesh_type), INTENT(IN) :: pp_mesh, vv_mesh
728  TYPE(petsc_csr_la) :: vv_3_LA, pp_1_LA
729  TYPE(periodic_type), INTENT(IN) :: vvz_per, pp_per
730  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: un
731  REAL(KIND=8), DIMENSION(pp_mesh%np,2,SIZE(list_mode)), INTENT(INOUT) :: pn
732  TYPE(dyn_real_array_three), DIMENSION(inputs%taylor_order-1) :: der_un
733  TYPE(dyn_real_array_three), DIMENSION(inputs%taylor_order-1) :: der_pn
734  !===Saved variables
735  INTEGER, SAVE :: m_max_c
736  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: pp_global_D
737  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: pp_mode_global_js_D
738  TYPE(dyn_int_line), DIMENSION(3), SAVE :: vv_js_D
739  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: vv_mode_global_js_D
740  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: vel_global_D
741  LOGICAL, SAVE :: once = .true.
742  INTEGER, SAVE :: my_petscworld_rank
743  !===End saved variables
744  INTEGER, POINTER, DIMENSION(:) :: pp_1_ifrom, vv_3_ifrom
745  INTEGER :: i, m, n, n1, n2, n3, n123
746  INTEGER :: nb_procs, code, nu_mat, mode
747  REAL(KIND=8) :: moyenne
748  REAL(KIND=8), DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: div
749  REAL(KIND=8), DIMENSION(vv_mesh%np, 6) :: un_p1
750  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rotv_v1, rotv_v2, rotv_v3, rotv_v4
751  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
752  REAL(KIND=8), DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: uext, dttt_un_p1, dtt_un_p1, dt_un_p1
753  REAL(KIND=8), DIMENSION(vv_mesh%np) :: vel_loc, vel_tot
754  REAL(KIND=8) :: tps, tps_tot, tps_cumul, coeff, vloc, cfl, cfl_max, norm
755  REAL(KIND=8) :: one, zero, three
756  DATA zero, one, three/0.d0,1.d0,3.d0/
757  !===Communicators for Petsc, in space and Fourier
758  petscerrorcode :: ierr
759  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
760  mat, DIMENSION(:), POINTER, SAVE :: vel_mat
761  mat, SAVE :: mass_mat
762  vec, SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
763  vec, SAVE :: pb_1, pb_2, px_1, px_1_ghost
764  ksp, DIMENSION(:), POINTER, SAVE :: vel_ksp
765  ksp, SAVE :: mass_ksp
766  !===END OF DECLARATION
767 
768  IF (once) THEN
769  once = .false.
770 
771  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
772  !===CREATE PETSC VECTORS AND GHOSTS
773  CALL create_my_ghost(vv_mesh,vv_3_la,vv_3_ifrom)
774  n = 3*vv_mesh%dom_np
775  CALL veccreateghost(comm_one_d(1), n, &
776  petsc_determine, SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
777  CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
778  CALL vecduplicate(vx_3, vb_3_145, ierr)
779  CALL vecduplicate(vx_3, vb_3_236, ierr)
780 
781  CALL create_my_ghost(pp_mesh,pp_1_la,pp_1_ifrom)
782  n = pp_mesh%dom_np
783  CALL veccreateghost(comm_one_d(1), n, &
784  petsc_determine, SIZE(pp_1_ifrom), pp_1_ifrom, px_1, ierr)
785  CALL vecghostgetlocalform(px_1, px_1_ghost, ierr)
786  CALL vecduplicate(px_1, pb_1, ierr)
787  CALL vecduplicate(px_1, pb_2, ierr)
788  !===End CREATE PETSC VECTORS AND GHOSTS
789 
790  !===Number of modes on proc
791  m_max_c = SIZE(list_mode)
792  !===End umber of modes on proc
793 
794  !===PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
795  CALL scalar_without_glob_js_d(pp_mesh, list_mode, pp_1_la, pp_mode_global_js_d)
796  ALLOCATE(pp_global_d(m_max_c))
797  DO i = 1, m_max_c
798  ALLOCATE(pp_global_d(i)%DRL(SIZE(pp_mode_global_js_d(i)%DIL)))
799  END DO
800  !===End PREPARE pp_mode_global_js_D ARRAY FOR PRESSURE
801 
802  !===PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
803  CALL vector_glob_js_d(vv_mesh, list_mode, vv_3_la, inputs%vv_list_dirichlet_sides, &
804  vv_js_d, vv_mode_global_js_d)
805 
806  ALLOCATE(vel_global_d(m_max_c))
807  DO i = 1, m_max_c
808  ALLOCATE(vel_global_d(i)%DRL(SIZE(vv_mode_global_js_d(i)%DIL)))
809  END DO
810  !===END PREPARE TYPE OF BOUNDARY CONDITIONS AND js_D ARRAY FOR VELOCITY
811 
812  !===ASSEMBLE MASS MATRIX
813  CALL create_local_petsc_matrix(comm_one_d(1), pp_1_la, mass_mat, clean=.false.)
814  CALL qs_diff_mass_scal_m (pp_mesh, pp_1_la, 0.d0, 1.d0, 0.d0, 0, mass_mat)
815  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
816  CALL periodic_matrix_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, mass_mat, pp_1_la)
817  END IF
818  DO i = 1, m_max_c
819  CALL dirichlet_m_parallel(mass_mat,pp_mode_global_js_d(i)%DIL)
820  END DO
821  CALL init_solver(inputs%my_par_mass,mass_ksp,mass_mat,comm_one_d(1),&
822  solver=inputs%my_par_mass%solver,precond=inputs%my_par_mass%precond)
823  !===END ASSEMBLE MASS MATRIX
824 
825  !===ASSEMBLING VELOCITY MATRICES
826  ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
827  DO i = 1, m_max_c
828  mode = list_mode(i)
829  !===VELOCITY
830  nu_mat = 2*i-1
831  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
832  CALL qs_diff_mass_vect_3x3_taylor (1, vv_3_la, vv_mesh, one/re, one/dt, & !HF new coeffs
833  lambda, i, mode, vel_mat(nu_mat))
834  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
835  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
836  vel_mat(nu_mat), vv_3_la)
837  END IF
838  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
839  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
840  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
841  nu_mat = nu_mat+1
842  CALL create_local_petsc_matrix(comm_one_d(1), vv_3_la, vel_mat(nu_mat), clean=.false.)
843  CALL qs_diff_mass_vect_3x3_taylor (2, vv_3_la, vv_mesh, one/re, one/dt, & !HF new coeffs
844  lambda, i, mode, vel_mat(nu_mat))
845  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
846  CALL periodic_matrix_petsc(vvz_per%n_bord, vvz_per%list,vvz_per%perlist, &
847  vel_mat(nu_mat), vv_3_la)
848  END IF
849  CALL dirichlet_m_parallel(vel_mat(nu_mat),vv_mode_global_js_d(i)%DIL)
850  CALL init_solver(inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
851  solver=inputs%my_par_vv%solver,precond=inputs%my_par_vv%precond)
852  !===End VELOCITY
853  ENDDO
854  !===End ASSEMBLING VELOCITY MATRICES
855 
856  ENDIF !===End of once
857 
858  !===Compute NL for dtt_un by FFT at Gauss points
859  tps_tot = user_time()
860  tps_cumul = 0
861  tps = user_time()
862 
863  uext = un +(4*dt/3)*der_un(1)%DRT + (4*dt/3)**2/2*der_un(2)%DRT + (4*dt/3)**3/6*der_un(3)%DRT !===1st NL term
864  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v1,inputs%precession)
865  !===rotv_v1 is now a non linear term
866 
867  uext = un + dt*der_un(1)%DRT + dt**2/2*der_un(2)%DRT + dt**3/6*der_un(3)%DRT !===2nd NL term
868  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v2,inputs%precession)
869  !===rotv_v2 is now a non linear term
870 
871  uext = un +(dt/2)*der_un(1)%DRT + (dt/2)**2/2*der_un(2)%DRT + (dt/2)**3/6*der_un(3)%DRT !===3rd NL term
872  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v3,inputs%precession)
873  !===rotv_v3 is now a non linear term
874 
875  uext = un !===4th NL term
876  CALL smb_cross_prod_gauss_sft_par(comm_one_d(2),vv_mesh,list_mode,uext,rotv_v4,inputs%precession)
877  !===rotv_v4 is now a non linear term
878  tps = user_time() - tps; tps_cumul=tps_cumul+tps
879  !===End Compute NL by FFT at Gauss points
880 
881  !===Computation of CFL
882  IF (inputs%verbose_CFL) THEN
883  vel_loc = 0.d0
884  DO i = 1, m_max_c
885  IF (list_mode(i)==0) THEN
886  coeff = 1.d0
887  ELSE
888  coeff = 0.5d0
889  END IF
890  vel_loc = vel_loc + coeff*(un(:,1,i)**2+un(:,2,i)**2+un(:,3,i)**2 &
891  +un(:,4,i)**2+un(:,5,i)**2+un(:,6,i)**2)
892  END DO
893  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
894  CALL mpi_allreduce(vel_loc,vel_tot,vv_mesh%np,mpi_double_precision, mpi_sum, comm_one_d(2), code)
895  vel_tot = sqrt(vel_tot)
896  cfl = 0.d0
897  DO m = 1, vv_mesh%dom_me
898  vloc = maxval(vel_tot(vv_mesh%jj(:,m)))
899  cfl = max(vloc*dt/min(vv_mesh%hloc(m),maxval(vv_mesh%hm)),cfl)
900  END DO
901  CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
902  talk_to_me%CFL=cfl_max
903  talk_to_me%time=time
904  END IF
905  !===End Computation of CFL
906 
907  !===Computation of third order derivatives
908  !===Computation of rhs for dttt_un at Gauss points for every mode
909  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
910  der_un(3)%DRT/dt, der_pn(3)%DRT, 9*(9*rotv_v1-20*rotv_v2+16*rotv_v3-5*rotv_v4)/(5*dt**3), rhs_gauss,3)
911  !===End Computation of rhs for dttt_un
912 
913  DO i = 1, m_max_c
914  tps = user_time()
915  mode = list_mode(i)
916 
917  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
918  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
919  !===Periodic boundary conditions
920  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
921  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
922  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
923  END IF
924  !===Axis boundary conditions
925  n1 = SIZE(vv_js_d(1)%DIL)
926  n2 = SIZE(vv_js_d(2)%DIL)
927  n3 = SIZE(vv_js_d(3)%DIL)
928  n123 = n1+n2+n3
929  vel_global_d(i)%DRL(1:n1) = &
930  ( vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
931  -3*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
932  +3*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt) &
933  - vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(dt**3)
934  vel_global_d(i)%DRL(n1+1:n1+n2) = &
935  ( vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
936  -3*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
937  +3*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt) &
938  - vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(dt**3)
939  vel_global_d(i)%DRL(n1+n2+1:n123)= &
940  ( vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
941  -3*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
942  +3*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt) &
943  - vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(dt**3)
944  vel_global_d(i)%DRL(n123+1:) = 0.d0
945  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
946  vel_global_d(i)%DRL(1:n1) = &
947  ( vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
948  -3*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
949  +3*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt) &
950  - vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(dt**3)
951  vel_global_d(i)%DRL(n1+1:n1+n2) = &
952  ( vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
953  -3*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
954  +3*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt) &
955  - vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(dt**3)
956  vel_global_d(i)%DRL(n1+n2+1:n123)= &
957  ( vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
958  -3*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
959  +3*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt) &
960  - vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(dt**3)
961  vel_global_d(i)%DRL(n123+1:) = 0.d0
962  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
963  tps = user_time() - tps; tps_cumul=tps_cumul+tps
964  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
965 
966  !===Solve linear system for momentum equation
967  tps = user_time()
968  !===Solve system 1, dtt_ur_c, dtt_ut_s, dtt_uz_c
969  nu_mat =2*i-1
970  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
971  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
972  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
973  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
974  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
975  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
976 
977  !===Solve system 2, dtt_ur_s, dtt_ut_c, dtt_uz_s
978  nu_mat = nu_mat + 1
979  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
980  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
981  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
982  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
983  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
984  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
985  !===dttt_un_p1 is temporarily in un_p1
986 
987  tps = user_time() - tps; tps_cumul=tps_cumul+tps
988  !===End Solve linear system for momentum equation
989 
990  !===Assemble divergence of dttt_velocity in arrays pb_1, pb_2
991  tps = user_time()
992  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
993  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
994 
995  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
996  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
997  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
998  END IF
999 
1000  !===Boundary condition on axis for dttt_pressure
1001  pp_global_d(i)%DRL = 0.d0
1002  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
1003  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
1004  !===End boundary condition on axis for dttt_pressure
1005 
1006  !===Solve mass matrix for dttt_pressure correction
1007  tps = user_time()
1008  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1009  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1010  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1011  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
1012  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1013  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1014  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1015  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
1016  !===End solve mass matrix for dttt_pressure correction
1017 
1018  !===Compute 3rd order derrivative pressure for this mode
1019  der_pn(3)%DRT(:,:,i) = der_pn(3)%DRT(:,:,i) - lambda*div(:,:,i)
1020  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1021  !===End 3rd order derrivative pressure
1022 
1023  !===UPDATES
1024  tps = user_time()
1025  !===Handling 2nd order derrivative of mean pressure
1026  IF (mode == 0) THEN
1027  CALL moy(comm_one_d(1),pp_mesh, der_pn(3)%DRT(:,1,i),moyenne)
1028  der_pn(3)%DRT(:,1,i) = der_pn(3)%DRT(:,1,i)-moyenne
1029  ENDIF
1030  !===End of handling of 3rd order derrivative mean pressure
1031 
1032  !===Correction of zero mode
1033  IF (mode==0) THEN
1034  un_p1(:,2) = 0.d0
1035  un_p1(:,4) = 0.d0
1036  un_p1(:,6) = 0.d0
1037  der_pn(3)%DRT (:,2,i) = 0.d0
1038  END IF
1039  !===Correction of zero mode
1040 
1041  !===UPDATES 3rd order derrivative
1042  dttt_un_p1(:,:,i) = un_p1
1043  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1044  !===End UPDATES 3rd order derrivative
1045 
1046  ENDDO
1047  !===END Computation of third order derivatives
1048 
1049  !===Computation of second order derivatives
1050  !===Computation of rhs for dtt_un at Gauss points for every mode
1051  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
1052  der_un(2)%DRT/dt-dt/2*(dttt_un_p1-der_un(3)%DRT)/dt, der_pn(2)%DRT+dt*der_pn(3)%DRT, &
1053  (81*rotv_v1-140*rotv_v2+64*rotv_v3-5*rotv_v4)/(10*dt**2), rhs_gauss,2)
1054  !===End Computation of rhs for dtt_un
1055 
1056  DO i = 1, m_max_c
1057  tps = user_time()
1058  mode = list_mode(i)
1059 
1060  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
1061  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
1062  !===Periodic boundary conditions
1063  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1064  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
1065  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
1066  END IF
1067  !===Axis boundary conditions
1068  n1 = SIZE(vv_js_d(1)%DIL)
1069  n2 = SIZE(vv_js_d(2)%DIL)
1070  n3 = SIZE(vv_js_d(3)%DIL)
1071  n123 = n1+n2+n3
1072  vel_global_d(i)%DRL(1:n1) = &
1073  (2*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
1074  -5*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
1075  +4*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt)&
1076  - vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(dt**2)
1077  vel_global_d(i)%DRL(n1+1:n1+n2) = &
1078  (2*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
1079  -5*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
1080  +4*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt)&
1081  - vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(dt**2)
1082  vel_global_d(i)%DRL(n1+n2+1:n123)= &
1083  (2*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
1084  -5*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
1085  +4*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt)&
1086  - vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(dt**2)
1087  vel_global_d(i)%DRL(n123+1:) = 0.d0
1088  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
1089 
1090  vel_global_d(i)%DRL(1:n1) = &
1091  (2*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
1092  -5*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
1093  +4*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt)&
1094  - vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(dt**2)
1095  vel_global_d(i)%DRL(n1+1:n1+n2) = &
1096  (2*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
1097  -5*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
1098  +4*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt)&
1099  - vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(dt**2)
1100  vel_global_d(i)%DRL(n1+n2+1:n123)= &
1101  (2*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
1102  -5*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
1103  +4*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt)&
1104  - vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(dt**2)
1105  vel_global_d(i)%DRL(n123+1:) = 0.d0 !===BUG JLG DEc 11 2019
1106  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
1107  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1108  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
1109 
1110  !===Solve linear system for momentum equation
1111  tps = user_time()
1112  !===Solve system 1, dtt_ur_c, dtt_ut_s, dtt_uz_c
1113  nu_mat =2*i-1
1114  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1115  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1116  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1117  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
1118  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
1119  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
1120 
1121  !===Solve system 2, dtt_ur_s, dtt_ut_c, dtt_uz_s
1122  nu_mat = nu_mat + 1
1123  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1124  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1125  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1126  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
1127  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
1128  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
1129  !rk: dtt_un_p1 is temporarly in un_p1
1130 
1131  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1132  !===End Solve linear system for momentum equation
1133 
1134  !===Assemble divergence of dtt_velocity in arrays pb_1, pb_2
1135  tps = user_time()
1136  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
1137  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
1138 
1139  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1140  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
1141  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
1142  END IF
1143 
1144  !===Boundary condition on axis for dtt_pressure
1145  pp_global_d(i)%DRL = 0.d0
1146  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
1147  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
1148  !===End boundary condition on axis for dtt_pressure
1149 
1150  !===Solve mass matrix for dtt_pressure correction
1151  tps = user_time()
1152  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1153  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1154  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1155  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
1156  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1157  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1158  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1159  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
1160  !===End solve mass matrix for dtt_pressure correction
1161 
1162  !===Compute 2nd order derrivative pressure for this mode
1163  der_pn(2)%DRT(:,:,i) = der_pn(2)%DRT(:,:,i) + dt*der_pn(3)%DRT(:,:,i)- lambda*div(:,:,i)
1164  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1165  !===End 2nd order derrivative Pressure computation we will then take mean out
1166 
1167  !===UPDATES
1168  tps = user_time()
1169  !===Handling 2nd order derrivative of mean pressure
1170  IF (mode == 0) THEN
1171  CALL moy(comm_one_d(1),pp_mesh, der_pn(2)%DRT(:,1,i),moyenne)
1172  der_pn(2)%DRT(:,1,i) = der_pn(2)%DRT(:,1,i)-moyenne
1173  ENDIF
1174  !===End of handling of 2nd order derrivative mean pressure
1175 
1176  !===Correction of zero mode
1177  IF (mode==0) THEN
1178  un_p1(:,2) = 0.d0
1179  un_p1(:,4) = 0.d0
1180  un_p1(:,6) = 0.d0
1181  der_pn(2)%DRT (:,2,i) = 0.d0
1182  END IF
1183  !===Correction of zero mode
1184 
1185  !===UPDATES 2nd order derrivative
1186  dtt_un_p1(:,:,i) = un_p1
1187  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1188  !===End UPDATES 2nd order derrivative
1189 
1190  ENDDO
1191  !===END Computation of second order derivatives
1192 
1193  !===Computation of first order derivatives
1194  !===Computation of rhs for dt_un at Gauss points for every mode
1195  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
1196  der_un(1)%DRT/dt-dt/2*(dtt_un_p1-der_un(2)%DRT)/dt-dt**2/12*(dttt_un_p1-der_un(3)%DRT)/dt, &
1197  der_pn(1)%DRT+dt*der_pn(2)%DRT-dt**2/2*der_pn(3)%DRT, (27*rotv_v1-32*rotv_v3+5*rotv_v4)/(20*dt), rhs_gauss,1)
1198  !===End Computation of rhs for dt_un
1199 
1200  DO i = 1, m_max_c
1201  tps = user_time()
1202  mode = list_mode(i)
1203 
1204  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
1205  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
1206  !===Periodic boundary conditions
1207  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1208  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
1209  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
1210  END IF
1211  !===Axis boundary conditions
1212  n1 = SIZE(vv_js_d(1)%DIL)
1213  n2 = SIZE(vv_js_d(2)%DIL)
1214  n3 = SIZE(vv_js_d(3)%DIL)
1215  n123 = n1+n2+n3
1216  vel_global_d(i)%DRL(1:n1) = &
1217  (11*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
1218  -18*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
1219  + 9*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt)&
1220  - 2*vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(6*dt)
1221  vel_global_d(i)%DRL(n1+1:n1+n2) = &
1222  (11*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
1223  -18*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
1224  + 9*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt)&
1225  - 2*vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(6*dt)
1226  vel_global_d(i)%DRL(n1+n2+1:n123)= &
1227  (11*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
1228  -18*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
1229  + 9*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt)&
1230  - 2*vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(6*dt)
1231  vel_global_d(i)%DRL(n123+1:) = 0.d0
1232  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
1233  vel_global_d(i)%DRL(1:n1) = &
1234  (11*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time) &
1235  -18*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-dt) &
1236  + 9*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-2*dt)&
1237  - 2*vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode,time-3*dt))/(6*dt)
1238  vel_global_d(i)%DRL(n1+1:n1+n2) = &
1239  (11*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time) &
1240  -18*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-dt) &
1241  + 9*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-2*dt)&
1242  - 2*vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode,time-3*dt))/(6*dt)
1243  vel_global_d(i)%DRL(n1+n2+1:n123)= &
1244  (11*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time) &
1245  -18*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-dt) &
1246  + 9*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-2*dt)&
1247  - 2*vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode,time-3*dt))/(6*dt)
1248  vel_global_d(i)%DRL(n123+1:) = 0.d0
1249  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
1250  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1251  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
1252 
1253  !===Solve linear system for momentum equation
1254  tps = user_time()
1255  !===Solve system 1, dt_ur_c, dt_ut_s, dt_uz_c
1256  nu_mat =2*i-1
1257  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1258  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1259  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1260  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
1261  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
1262  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
1263  !===dt_un_p1 is temporarily in un_p1
1264 
1265  !===Solve system 2, dt_ur_s, dt_ut_c, dt_uz_s
1266  nu_mat = nu_mat + 1
1267  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1268  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1269  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1270  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
1271  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
1272  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
1273 
1274  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1275  !===End Solve linear system for momentum equation
1276 
1277  !===Assemble divergence of dt_velocity in arrays pb_1, pb_2
1278  tps = user_time()
1279  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
1280  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
1281 
1282  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1283  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
1284  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
1285  END IF
1286 
1287  !===Boundary condition on axis for dt_pressure
1288  pp_global_d(i)%DRL = 0.d0
1289  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
1290  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
1291  !===End boundary condition on axis for dt_pressure
1292 
1293  !===Solve mass matrix for dt_pressure correction
1294  tps = user_time()
1295  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1296  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1297  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1298  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
1299  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1300  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1301  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1302  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
1303  !===End solve mass matrix for dt_pressure correction
1304 
1305  !===Compute first order derrivative pressure for this mode
1306  !===Pressure first order derrivative computation
1307  der_pn(1)%DRT(:,:,i) = der_pn(1)%DRT(:,:,i) + dt*der_pn(2)%DRT(:,:,i) - dt**2/2*der_pn(3)%DRT(:,:,i) - lambda*div(:,:,i)
1308  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1309  !===End dt_Pressure computation
1310 
1311  !===UPDATES
1312  tps = user_time()
1313  !===Handling of mean first order derrivative pressure
1314  IF (mode == 0) THEN
1315  CALL moy(comm_one_d(1),pp_mesh, der_pn(1)%DRT(:,1,i),moyenne)
1316  der_pn(1)%DRT(:,1,i) = der_pn(1)%DRT(:,1,i)-moyenne
1317  ENDIF
1318  !===End of handling of mean first order derrivative pressure
1319 
1320  !===Correction of zero mode
1321  IF (mode==0) THEN
1322  un_p1(:,2) = 0.d0
1323  un_p1(:,4) = 0.d0
1324  un_p1(:,6) = 0.d0
1325  der_pn(1)%DRT (:,2,i) = 0.d0
1326  END IF
1327  !===Correction of zero mode
1328 
1329  !===UPDATE Velocity first order derrivative
1330  dt_un_p1(:,:,i) = un_p1
1331  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1332  !===End UPDATE Velocity first order derrivative
1333 
1334  ENDDO !=== i = 1, m_max_c
1335  !===END Computation of first order derivatives
1336 
1337  !===Computation of 0th order derivatives
1338  !===Computation of rhs for un at Gauss points for every mode
1339  CALL rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, comm_one_d(2), list_mode, time, &
1340  (un-dt/2*(dt_un_p1-der_un(1)%DRT)-dt**2/12*(dtt_un_p1-der_un(2)%DRT))/dt, &
1341  pn+dt*der_pn(1)%DRT-dt**2/2*der_pn(2)%DRT+dt**3/6*der_pn(3)%DRT,rotv_v2,rhs_gauss,0)
1342  !===End Computation of rhs for un
1343 
1344  DO i = 1, m_max_c
1345  tps = user_time()
1346  mode = list_mode(i)
1347 
1348  !===Assemble vb_3_145, vb_3_236 using rhs_gauss
1349  CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
1350  !===Periodic boundary conditions
1351  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1352  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
1353  CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
1354  END IF
1355  !===Axis boundary conditions
1356  n1 = SIZE(vv_js_d(1)%DIL)
1357  n2 = SIZE(vv_js_d(2)%DIL)
1358  n3 = SIZE(vv_js_d(3)%DIL)
1359  n123 = n1+n2+n3
1360  vel_global_d(i)%DRL(1:n1) = vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
1361  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
1362  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
1363  vel_global_d(i)%DRL(n123+1:) = 0.d0
1364  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
1365  vel_global_d(i)%DRL(1:n1) = vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
1366  vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
1367  vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
1368  vel_global_d(i)%DRL(n123+1:) = 0.d0
1369  CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
1370  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1371  !===End Assemble vb_3_145, vb_3_236 using rhs_gauss
1372 
1373  !===Solve linear system for momentum equation
1374  tps = user_time()
1375  !===Solve system 1, ur_c, ut_s, uz_c
1376  nu_mat =2*i-1
1377  CALL solver(vel_ksp(nu_mat),vb_3_145,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1378  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1379  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1380  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
1381  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
1382  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
1383 
1384  !===Solve system 2, ur_s, ut_c, uz_s
1385  nu_mat = nu_mat + 1
1386  CALL solver(vel_ksp(nu_mat),vb_3_236,vx_3,reinit=.false.,verbose=inputs%my_par_vv%verbose)
1387  CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
1388  CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
1389  CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
1390  CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
1391  CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
1392 
1393  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1394  !===End Solve linear system for momentum equation
1395 
1396  !===Assemble divergence of velocity in arrays pb_1, pb_2
1397  tps = user_time()
1398  CALL qs_01_div_hybrid_generic(inputs%type_fe_velocity,vv_mesh, pp_mesh, pp_1_la, mode, un_p1, pb_1, pb_2)
1399  !===End assembling; pb_1, and pb_2 are petsc vectors for the rhs divergence
1400 
1401  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
1402  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_1, pp_1_la)
1403  CALL periodic_rhs_petsc(pp_per%n_bord, pp_per%list, pp_per%perlist, pb_2, pp_1_la)
1404  END IF
1405 
1406  !===Boundary condition on axis for pressure
1407  pp_global_d(i)%DRL = 0.d0
1408  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
1409  CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
1410  !===End boundary condition on axis for pressure
1411 
1412  !===Solve mass matrix for pressure correction
1413  tps = user_time()
1414  CALL solver(mass_ksp,pb_1,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1415  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1416  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1417  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
1418  CALL solver(mass_ksp,pb_2,px_1,reinit=.false.,verbose=inputs%my_par_mass%verbose)
1419  CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
1420  CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
1421  CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
1422  !===End solve mass matrix for pressure correction
1423 
1424  !===Pressure computation
1425  pn(:,:,i) = pn(:,:,i)+dt*der_pn(1)%DRT(:,:,i)-dt**2/2*der_pn(2)%DRT(:,:,i)+dt**3/6*der_pn(3)%DRT(:,:,i)-lambda*div(:,:,i)
1426  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1427  !===End Pressure computation
1428 
1429  !===UPDATES
1430  tps = user_time()
1431  !===Handling of mean pressure
1432  IF (mode == 0) THEN
1433  CALL moy(comm_one_d(1),pp_mesh, pn(:,1,i),moyenne)
1434  pn(:,1,i) = pn(:,1,i)-moyenne
1435  ENDIF
1436  !===End of handling of mean pressure
1437 
1438  !===Correction of zero mode
1439  IF (mode==0) THEN
1440  un_p1(:,2) = 0.d0
1441  un_p1(:,4) = 0.d0
1442  un_p1(:,6) = 0.d0
1443  pn(:,2,i) = 0.d0
1444  END IF
1445  !===Correction of zero mode
1446 
1447  !===UPDATES Velocity 0th order derrivative
1448  un(:,:,i) = un_p1
1449  tps = user_time() - tps; tps_cumul=tps_cumul+tps
1450  !===End UPDATES Velocity 0th order derrivative
1451 
1452  ENDDO !=== i = 1, m_max_c
1453  !===END Computation of 0th order derivatives
1454 
1455  !===update dttt_un, dtt_un, and dt_un once everyone has been compuded
1456  DO i = 1, m_max_c
1457  der_un(3)%DRT(:,:,i)=dttt_un_p1(:,:,i)
1458  der_un(2)%DRT(:,:,i)= dtt_un_p1(:,:,i)
1459  der_un(1)%DRT(:,:,i)= dt_un_p1(:,:,i)
1460  END DO
1461 
1462  IF (inputs%verbose_divergence) THEN
1463  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
1464  talk_to_me%div_L2 = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)/norm
1465  talk_to_me%weak_div_L2 = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, div)/norm
1466  END IF
1467 
1468  tps_tot = user_time() - tps_tot
1469  END SUBROUTINE update_ns_with_taylor_fourth
1470  !============================================
1471 
1472 
1473  !===PRECESSION 28/07/09
1474  SUBROUTINE smb_cross_prod_gauss_sft_par(communicator,mesh,list_mode,V_in,V_out,precession_in)
1475  !=================================
1476  USE sft_parallele
1477  USE chaine_caractere
1478  USE input_data
1479  USE def_type_mesh
1480  IMPLICIT NONE
1481 
1482  TYPE(mesh_type), INTENT(IN) :: mesh
1483  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1484  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: V_in
1485  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: V_out
1486  LOGICAL, INTENT(IN) :: precession_in
1487  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: RotV, W, Om
1488  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1489  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1490  INTEGER :: m, l , i, mode, index, k
1491  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: Vs, Omega_s
1492  REAL(KIND=8) :: ray
1493  REAL(KIND=8) :: tps, PI
1494  REAL(KIND=8), DIMENSION(3) :: temps
1495  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: Omega
1496  LOGICAL, SAVE :: once=.true.
1497  !===FOR FFT_PAR_CROSS_PROD_DCL
1498  INTEGER :: nb_procs, bloc_size, m_max_pad, code
1499  mpi_comm :: communicator
1500 
1501  IF (once) THEN
1502  once = .false.
1503  ALLOCATE(omega(mesh%np,6,SIZE(list_mode)))
1504  omega = 0.d0
1505  pi = acos(-1.d0)
1506  DO i=1, SIZE(list_mode)
1507  IF (list_mode(i) == 1) THEN
1508  !precession selon un axe penche d un angle angle_s_pi*PI
1509  omega(:,1,i) = inputs%taux_precession * sin(inputs%angle_s_pi*pi)
1510  omega(:,4,i) = -inputs%taux_precession * sin(inputs%angle_s_pi*pi)
1511  ELSE IF (list_mode(i) == 0) THEN
1512  !precession selon un axe penche d un angle angle_s_pi*PI
1513  omega(:,5,i) = inputs%taux_precession * cos(inputs%angle_s_pi*pi)
1514  ENDIF
1515  ENDDO
1516  ENDIF ! fin du once
1517 
1518  tps = user_time()
1519  DO i = 1, SIZE(list_mode)
1520  mode = list_mode(i)
1521  index = 0
1522  DO m = 1, mesh%me
1523  j_loc = mesh%jj(:,m)
1524  DO k = 1, 6
1525  vs(:,k) = v_in(j_loc,k,i)
1526  omega_s(:,k) = omega(mesh%jj(:,m),k,i)
1527  END DO
1528 
1529  DO l = 1, mesh%gauss%l_G
1530  index = index + 1
1531  dw_loc = mesh%gauss%dw(:,:,l,m)
1532 
1533  !===Compute radius of Gauss point
1534  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1535 
1536  !-----------------vitesse sur les points de Gauss---------------------------
1537  w(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
1538  w(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
1539  w(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
1540 
1541  w(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
1542  w(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
1543  w(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
1544  !-----------------vecteur rotation sur les points de Gauss---------------------------
1545  om(index,1,i) = sum(omega_s(:,1)*mesh%gauss%ww(:,l))
1546  om(index,3,i) = sum(omega_s(:,3)*mesh%gauss%ww(:,l))
1547  om(index,5,i) = sum(omega_s(:,5)*mesh%gauss%ww(:,l))
1548 
1549  om(index,2,i) = sum(omega_s(:,2)*mesh%gauss%ww(:,l))
1550  om(index,4,i) = sum(omega_s(:,4)*mesh%gauss%ww(:,l))
1551  om(index,6,i) = sum(omega_s(:,6)*mesh%gauss%ww(:,l))
1552  !-----------------rotational sur les points de Gauss---------------------------
1553  !coeff sur les cosinus
1554  rotv(index,1,i) = mode/ray*w(index,6,i) &
1555  -sum(vs(:,3)*dw_loc(2,:))
1556  rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
1557  -sum(vs(:,6)*dw_loc(1,:))
1558  rotv(index,5,i) = 1/ray*w(index,3,i) &
1559  +sum(vs(:,3)*dw_loc(1,:)) &
1560  -mode/ray*w(index,2,i)
1561  !coeff sur les sinus
1562  rotv(index,2,i) =-mode/ray*w(index,5,i) &
1563  -sum(vs(:,4)*dw_loc(2,:))
1564  rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
1565  -sum(vs(:,5)*dw_loc(1,:))
1566  rotv(index,6,i) = 1/ray*w(index,4,i) &
1567  +sum(vs(:,4)*dw_loc(1,:))&
1568  +mode/ray*w(index,1,i)
1569  ENDDO
1570  ENDDO
1571  END DO
1572 
1573  !JLG, FL, July 23 2010, There was a Bug here
1574  IF (.NOT.precession_in) THEN
1575  om=0.d0
1576  END IF
1577  !
1578  ! cas PRECESSION 28/07/09
1579  !take care to the -F= 2 Om x U term in LHS
1580  rotv = rotv + 2.d0 * om
1581  temps = 0
1582 
1583 
1584  CALL mpi_comm_size(communicator, nb_procs, code)
1585  bloc_size = SIZE(rotv,1)/nb_procs+1
1586  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1587  CALL fft_par_cross_prod_dcl(communicator, rotv, w, v_out, nb_procs, bloc_size, m_max_pad, temps)
1588  tps = user_time() - tps
1589 
1590  END SUBROUTINE smb_cross_prod_gauss_sft_par
1591 
1592  SUBROUTINE moy(communicator,mesh,p,RESLT)
1593  !===========================
1594  USE def_type_mesh
1595  IMPLICIT NONE
1596  TYPE(mesh_type) :: mesh
1597  REAL(KIND=8), DIMENSION(:) , INTENT(IN) :: p
1598  REAL(KIND=8) , INTENT(OUT) :: RESLT
1599  REAL(KIND=8) :: vol_loc, vol_out, r_loc, r_out
1600  INTEGER :: m, l , i , ni, code
1601  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1602  REAL(KIND=8) :: ray
1603  mpi_comm :: communicator
1604  r_loc = 0.d0
1605  vol_loc = 0.d0
1606 
1607  DO m = 1, mesh%dom_me
1608  j_loc = mesh%jj(:,m)
1609  DO l = 1, mesh%gauss%l_G
1610  ray = 0
1611  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
1612  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1613  END DO
1614 
1615  r_loc = r_loc + sum(p(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
1616  vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
1617 
1618  ENDDO
1619  ENDDO
1620 
1621  CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
1622  CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
1623  reslt = r_out / vol_out
1624 
1625  END SUBROUTINE moy
1626 
1627  SUBROUTINE rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, rotv_v, &
1628  rhs_gauss, der_ord)
1629  !=================================
1630  !RHS for Navier-Stokes_taylor : RHS for (d/dt)**der_ord un
1631  USE def_type_mesh
1632  USE my_util
1633  USE input_data
1634  USE fem_rhs_axi
1635  USE sft_parallele
1636  USE boundary
1637  IMPLICIT NONE
1638  TYPE(mesh_type) :: vv_mesh, pp_mesh
1639  REAL(KIND=8), INTENT(IN) :: time
1640  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotv_v
1641  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: V1m
1642  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn
1643  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
1644  INTEGER, INTENT(IN) :: der_ord
1645  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: rhs_gauss
1646  REAL(KIND=8), DIMENSION(6) :: fs, ft
1647  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
1648  REAL(KIND=8), DIMENSION(vv_mesh%np,6) :: ff, imposed_vel
1649  REAL(KIND=8), DIMENSION(vv_mesh%np,2) :: P
1650  REAL(KIND=8), DIMENSION(vv_mesh%gauss%k_d,vv_mesh%gauss%n_w) :: dw_loc
1651  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6) :: imposed_vel_gauss
1652  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6,SIZE(list_mode)) :: fp, rhs_gauss_penal
1653  REAL(KIND=8), DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
1654  REAL(KIND=8) :: ray
1655  INTEGER :: m, l , i, k, index, TYPE
1656  INTEGER :: nb_procs, m_max_pad, bloc_size
1657 #include "petsc/finclude/petsc.h"
1658  petscerrorcode :: ierr
1659  mpi_comm :: communicator
1660 
1661  IF ((der_ord-0)*(der_ord-1)*(der_ord-2)*(der_ord-3)/=0)THEN
1662  CALL error_petsc('BUG in rhs_ns_gauss_3x3_taylor, der_ord not correct')
1663  END IF
1664  DO i = 1, SIZE(list_mode)
1665  DO k= 1, 6 !===external forces
1666  IF (der_ord==0) THEN
1667  ff(:,k) =source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time, inputs%Re,'ns')
1668  ELSEIF (der_ord==1) THEN
1669  ff(:,k) = &
1670  (11*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time, inputs%Re,'ns') &
1671  -18*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-inputs%dt, inputs%Re,'ns') &
1672  + 9*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-2*inputs%dt, inputs%Re,'ns') &
1673  - 2*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-3*inputs%dt, inputs%Re,'ns'))/(6*inputs%dt)
1674  ELSEIF (der_ord==2) THEN
1675  ff(:,k) = &
1676  (2*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time, inputs%Re,'ns') &
1677  -5*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-inputs%dt, inputs%Re,'ns') &
1678  +4*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-2*inputs%dt, inputs%Re,'ns') &
1679  - source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-3*inputs%dt, inputs%Re,'ns'))/(inputs%dt**2)
1680  ELSEIF (der_ord==3) THEN
1681  ff(:,k) = &
1682  ( source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time, inputs%Re,'ns') &
1683  -3*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-inputs%dt, inputs%Re,'ns') &
1684  +3*source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-2*inputs%dt, inputs%Re,'ns') &
1685  - source_in_ns_momentum(k, vv_mesh%rr, list_mode(i),i, time-3*inputs%dt, inputs%Re,'ns'))/(inputs%dt**3)
1686  END IF
1687  END DO
1688  DO k = 1, 2 !===pressure
1689  !===JLG+HF July 24th, 2019
1690  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn(:,k,i), p(:,k))
1691  ENDDO
1692 
1693  index = 0
1694  DO m = 1, vv_mesh%dom_me
1695  j_loc = vv_mesh%jj(:,m)
1696 
1697  DO l = 1, vv_mesh%gauss%l_G !===compute values on gauss pts
1698  index = index +1
1699  dw_loc = vv_mesh%gauss%dw(:,:,l,m)
1700 
1701  !===Compute radius of Gauss point
1702  ray = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
1703  rr_gauss(1,index) = ray
1704  rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
1705 
1706  !=== u0(1,:) <--> f(r,m,c)
1707  fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l)) !external forcing
1708  ft(1) = sum(v1m(j_loc,1,i) * vv_mesh%gauss%ww(:,l)) !inertia
1709  fp(index,1,i) = -sum(p(j_loc,1)*dw_loc(1,:)) !pressure term
1710  !=== u0(2,:) <--> f(r,m,s)
1711  fs(2) = sum(ff(j_loc,2) * vv_mesh%gauss%ww(:,l))
1712  ft(2) = sum(v1m(j_loc,2,i) * vv_mesh%gauss%ww(:,l))
1713  fp(index,2,i) = -sum(p(j_loc,2)*dw_loc(1,:))
1714  !=== u0(3,:) <--> f(th,m,c)
1715  fs(3) = sum(ff(j_loc,3) * vv_mesh%gauss%ww(:,l))
1716  ft(3) = sum(v1m(j_loc,3,i) * vv_mesh%gauss%ww(:,l))
1717  fp(index,3,i) = -sum(p(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
1718  !=== u0(4,:) <--> f(th,m,s)
1719  fs(4) = sum(ff(j_loc,4) * vv_mesh%gauss%ww(:,l))
1720  ft(4) = sum(v1m(j_loc,4,i) * vv_mesh%gauss%ww(:,l))
1721  fp(index,4,i) = sum(p(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
1722  !=== u0(5,:) <--> f(z,m,c)
1723  fs(5) = sum(ff(j_loc,5) * vv_mesh%gauss%ww(:,l))
1724  ft(5) = sum(v1m(j_loc,5,i) * vv_mesh%gauss%ww(:,l))
1725  fp(index,5,i) = -sum(p(j_loc,1)*dw_loc(2,:))
1726  !=== u0(6,:) <--> f(z,m,s)
1727  fs(6) = sum(ff(j_loc,6) * vv_mesh%gauss%ww(:,l))
1728  ft(6) = sum(v1m(j_loc,6,i) * vv_mesh%gauss%ww(:,l))
1729  fp(index,6,i) = -sum(p(j_loc,2)*dw_loc(2,:))
1730 
1731  rhs_gauss(index,:,i) = (ft+fs-rotv_v(index,:,i))
1732 
1733  ENDDO
1734  ENDDO
1735  IF (inputs%if_ns_penalty) THEN
1736  IF(inputs%if_impose_vel_in_solids) THEN
1737  IF (list_mode(i)==0) THEN
1738  IF (der_ord==0) THEN
1739  imposed_vel(:,:) = imposed_velocity_by_penalty(vv_mesh%rr,time)/(inputs%dt)
1740  index = 0
1741  DO m = 1, vv_mesh%dom_me
1742  j_loc = vv_mesh%jj(:,m)
1743  DO l = 1, vv_mesh%gauss%l_G
1744  index = index +1
1745  DO TYPE = 1, 6
1746  imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1747  * vv_mesh%gauss%ww(:,l))
1748  END DO
1749  END DO
1750  END DO
1751  rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1752  ELSE IF (der_ord==1) THEN
1753  !===CN HF adapt to 4th order
1754  ! imposed_vel(:,:) = (imposed_velocity_by_penalty(vv_mesh%rr,time+inputs%dt) &
1755  ! - imposed_velocity_by_penalty(vv_mesh%rr,time-inputs%dt))/(2.d0*inputs%dt**2)
1756  imposed_vel(:,:) = &
1757  (-2*imposed_velocity_by_penalty(vv_mesh%rr,time-3*inputs%dt)&
1758  +9*imposed_velocity_by_penalty(vv_mesh%rr,time-2*inputs%dt)&
1759  -18*imposed_velocity_by_penalty(vv_mesh%rr,time-1*inputs%dt)&
1760  +11*imposed_velocity_by_penalty(vv_mesh%rr,time))/(6*inputs%dt**2)
1761  !===CN HF adapt to 4th order
1762  index = 0
1763  DO m = 1, vv_mesh%dom_me
1764  j_loc = vv_mesh%jj(:,m)
1765  DO l = 1, vv_mesh%gauss%l_G
1766  index = index +1
1767  DO TYPE = 1, 6
1768  imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1769  * vv_mesh%gauss%ww(:,l))
1770  END DO
1771  END DO
1772  END DO
1773  rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1774  ELSE IF (der_ord==2) THEN
1775  !===CN HF adapt to 4th order: nothing to do
1776  imposed_vel(:,:) = &
1777  (-1*imposed_velocity_by_penalty(vv_mesh%rr,time-3*inputs%dt)&
1778  +4*imposed_velocity_by_penalty(vv_mesh%rr,time-2*inputs%dt)&
1779  -5*imposed_velocity_by_penalty(vv_mesh%rr,time-1*inputs%dt)&
1780  +2*imposed_velocity_by_penalty(vv_mesh%rr,time))/(inputs%dt**3)
1781  !===CN HF adapt to 4th order: nothing to do
1782  index = 0
1783  DO m = 1, vv_mesh%dom_me
1784  j_loc = vv_mesh%jj(:,m)
1785  DO l = 1, vv_mesh%gauss%l_G
1786  index = index +1
1787  DO TYPE = 1, 6
1788  imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1789  * vv_mesh%gauss%ww(:,l))
1790  END DO
1791  END DO
1792  END DO
1793  rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1794  !===CN HF adapt to 4th order
1795  ELSE IF (der_ord==3) THEN
1796  imposed_vel(:,:) = &
1797  (-1*imposed_velocity_by_penalty(vv_mesh%rr,time-3*inputs%dt)&
1798  +3*imposed_velocity_by_penalty(vv_mesh%rr,time-2*inputs%dt)&
1799  -3*imposed_velocity_by_penalty(vv_mesh%rr,time-1*inputs%dt)&
1800  +1*imposed_velocity_by_penalty(vv_mesh%rr,time))/(inputs%dt**4)
1801  index = 0
1802  DO m = 1, vv_mesh%dom_me
1803  j_loc = vv_mesh%jj(:,m)
1804  DO l = 1, vv_mesh%gauss%l_G
1805  index = index +1
1806  DO TYPE = 1, 6
1807  imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1808  * vv_mesh%gauss%ww(:,l))
1809  END DO
1810  END DO
1811  END DO
1812  rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1813  !===CN HF adapt to 4th order
1814  ENDIF
1815  END IF
1816  END IF
1817  END IF
1818  END DO
1819 
1820  IF (inputs%if_ns_penalty) THEN
1821  IF(inputs%if_impose_vel_in_solids) THEN
1822  CALL mpi_comm_size(communicator, nb_procs, ierr)
1823  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1824  bloc_size = SIZE(rhs_gauss,1)/nb_procs+1
1825 
1826  CALL fft_par_var_eta_prod_gauss_dcl(communicator, penal_in_real_space, vv_mesh, &
1827  rhs_gauss, rhs_gauss_penal, nb_procs, bloc_size, m_max_pad, rr_gauss, time)
1828  DO i = 1, SIZE(list_mode)
1829  IF (list_mode(i)==0) THEN
1830  rhs_gauss(:,:,i) = rhs_gauss_penal(:,:,i) + imposed_vel_gauss(:,:)
1831  ELSE
1832  rhs_gauss(:,:,i) = rhs_gauss_penal(:,:,i)
1833  END IF
1834  END DO
1835  END IF
1836  END IF
1837 
1838  rhs_gauss = rhs_gauss + fp
1839 
1840  END SUBROUTINE rhs_ns_gauss_3x3_taylor
1841 
1842  SUBROUTINE qs_diff_mass_vect_3x3_taylor(type_op, LA, mesh, visco, mass, lambda, i_mode, mode, matrix)
1844 #include "petsc/finclude/petsc.h"
1845  USE petsc
1846  USE my_util
1847  USE input_data
1848  IMPLICIT NONE
1849 
1850  INTEGER , INTENT(IN) :: type_op, mode, i_mode
1851  REAL(KIND=8), INTENT(IN) :: visco, mass, lambda
1852  TYPE(petsc_csr_la) :: LA
1853  TYPE(mesh_type), TARGET :: mesh
1854 
1855  INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
1856  REAL(KIND=8) :: xij, viscolm, div_penal
1857  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1858  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
1859  REAL(KIND=8) :: ray, eps1, eps2, z
1860  REAL(KIND=8) :: two = 2.d0
1861  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
1862  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
1863  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
1864  INTEGER, DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
1865  INTEGER :: ix, jx, iglob, jglob
1866  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
1867  REAL(KIND=8) :: viscomode, hm
1868  !#include "petsc/finclude/petsc.h"
1869  mat :: matrix
1870  petscerrorcode :: ierr
1871 
1872  CALL matzeroentries (matrix, ierr)
1873  CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1874 
1875  np = SIZE(mesh%rr,2)
1876  n_w = mesh%gauss%n_w
1877  n_ws = mesh%gauss%n_ws
1878 
1879  IF (type_op == 1) THEN
1880  eps1 = 1.d0
1881  eps2 = 1.d0
1882  k_max = 2 !Structure vectorielle
1883  ELSEIF (type_op == 2) THEN
1884  eps1 = 1.d0
1885  eps2 = -1.d0
1886  k_max = 2 !Structure vectorielle
1887  ELSEIF (type_op == 3) THEN
1888  !cas du laplacien scalaire
1889  eps1 = 0.d0
1890  eps2 = 0.d0
1891  k_max = 1 !Structure scalaire
1892  ELSE
1893  eps1 = 0.d0
1894  eps2 = 0.d0
1895  k_max = 1
1896  CALL error_petsc('probleme de type d''operateur')
1897  ENDIF
1898 
1899 
1900  DO l = 1, mesh%gauss%l_G
1901  DO ni = 1, mesh%gauss%n_w
1902  DO nj = 1, mesh%gauss%n_w
1903  wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1904  END DO
1905  END DO
1906  END DO
1907 
1908  DO m = 1, mesh%me
1909  jj_loc = mesh%jj(:,m)
1910 
1911  aij = 0.d0
1912  bij = 0.d0
1913  cij = 0.d0
1914  DO l = 1, mesh%gauss%l_G
1915 
1916  !--------On calcule le rayon du point gauss
1917  ray = 0
1918  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1919  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1920  END DO
1921 
1922  viscolm = visco*mesh%gauss%rj(l,m)
1923  !hm=0.5d0/inputs%m_max
1924  !hm=mesh%hm(i_mode) !(JLG April 7 2017)
1925  hm=min(mesh%hm(i_mode),mesh%hloc(m))!WRONG choice
1926  viscomode = visco*mesh%gauss%rj(l,m)
1927 
1928  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1929 
1930  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1931 
1932  !grad(u).grad(v) en r et z
1933  xij = 0.d0
1934  DO k = 1, mesh%gauss%k_d
1935  xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
1936  END DO
1937 
1938  !blocs diagonaux
1939  z = ray * viscolm* xij &
1940  + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1941  + viscomode*mode**2*wwprod(ni,nj,l)/ray
1942  cij(ni,nj) = cij(ni,nj) + z
1943  aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
1944  !blocs couplant
1945  bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
1946  ENDDO
1947  ENDDO
1948 
1949  ENDDO
1950 
1951  mat_loc = 0.d0
1952  DO ki= 1, 3
1953  DO ni = 1, n_w
1954  i = jj_loc(ni)
1955  iglob = la%loc_to_glob(ki,i)
1956  ix = (ki-1)*n_w+ni
1957  idxn(ix) = iglob-1
1958  DO kj = 1, 3
1959  DO nj = 1, n_w
1960  j = jj_loc(nj)
1961  jglob = la%loc_to_glob(kj,j)
1962  jx = (kj-1)*n_w+nj
1963  jdxn(jx) = jglob-1
1964  IF ((ki .LT. 3) .AND. (kj .LT. 3)) THEN
1965  IF (ki==kj) THEN
1966  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1967  ELSE
1968  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1969  END IF
1970  ELSE ! ki=3 OR kj=3
1971  IF (ki==kj) THEN
1972  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1973  END IF
1974  END IF
1975  END DO
1976  END DO
1977  END DO
1978  END DO
1979  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1980 
1981  !+++---------------------
1982 !!$ IF (type_op==3) CYCLE
1983  !==Calcul de visco (grad u)T . (grad v)
1984  aij = 0.d0
1985  bij = 0.d0
1986  cij = 0.d0
1987  dij = 0.d0
1988  eij = 0.d0
1989  fij = 0.d0
1990  DO l = 1, mesh%gauss%l_G
1991 
1992  !--------On calcule le rayon du point gauss
1993  ray = 0
1994  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1995  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1996  END DO
1997 
1998  viscolm = visco*mesh%gauss%rj(l,m)*ray
1999  !HF 2019/04/11
2000  div_penal = lambda*mesh%gauss%rj(l,m)*ray
2001  !HF 2019/04/11
2002 
2003  DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
2004  DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
2005  aij(ni,nj) = aij(ni,nj) &
2006  + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
2007  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
2008  + mesh%gauss%ww(nj,l)/ray)
2009  bij(ni,nj) = bij(ni,nj) &
2010  + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
2011  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
2012  cij(ni,nj) = cij(ni,nj) &
2013  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
2014  + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
2015  dij(ni,nj) = dij(ni,nj) &
2016  + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
2017  -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
2018  + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
2019  eij(ni,nj) = eij(ni,nj) &
2020  + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
2021  + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
2022  fij(ni,nj) = fij(ni,nj) &
2023  + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
2024  + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
2025  END DO
2026  END DO
2027  END DO
2028  !++++++++++++++++
2029  mat_loc=0.d0
2030  idxn=0
2031  jdxn=0
2032  DO ni = 1, n_w
2033  DO ki = 1, 3
2034  i = jj_loc(ni)
2035  iglob = la%loc_to_glob(ki,i)
2036  ix = (ki-1)*n_w + ni
2037  idxn(ix) = iglob-1
2038  END DO
2039  END DO
2040  jdxn=idxn
2041 
2042  DO ni = 1, n_w
2043  DO nj = 1, n_w
2044  !=== Line i 1 (Vr)
2045  ix = ni
2046  !=== Column j 1 (Vr)
2047  jx = nj
2048  mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
2049  !=== Column j 2 (Vt)
2050  jx = nj+n_w
2051  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
2052  !=== Column j 3 (Vz)
2053  jx = nj+2*n_w
2054  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
2055 
2056  !=== Line i 2 (Vt)
2057  ix = ni+n_w
2058  !=== Column j 1 (Vr)
2059  jx = nj
2060  mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
2061  !=== Column j 2 (Vt)
2062  jx = nj+n_w
2063  mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
2064  !=== Column j 3 (Vz)
2065  jx = nj+2*n_w
2066  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
2067 
2068  !=== Line i 3 (Vz)
2069  ix = ni+2*n_w
2070  !=== Column j 1 (Vr)
2071  jx = nj
2072  mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
2073  !=== Column j 2 (Vt)
2074  jx = nj+n_w
2075  mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
2076  !=== Column j 3 (Vz)
2077  jx = nj+2*n_w
2078  mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
2079  END DO
2080  END DO
2081  CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
2082 
2083  ENDDO
2084  !== Fin du Calcul de visco (grad u)T . (grad v)
2085  !++++++++++++++++------
2086 
2087  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
2088  !===Surface terms
2089  DO ms = 1, mesh%mes
2090  IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
2091  aij = 0.d0
2092  bij = 0.d0
2093  cij = 0.d0
2094  dij = 0.d0
2095  DO ls = 1, mesh%gauss%l_Gs
2096  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
2097  IF (ray.LT.1.d-10) cycle
2098  z = two*mesh%gauss%rjs(ls,ms)*ray*visco
2099 
2100  DO ni = 1, mesh%gauss%n_ws
2101  DO nj = 1, mesh%gauss%n_ws
2102  aij(ni,nj) = aij(ni,nj) - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
2103  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
2104  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
2105  bij(ni,nj) = bij(ni,nj) - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
2106  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
2107  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
2108  cij(ni,nj) = cij(ni,nj) - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
2109  *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
2110  + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
2111  dij(ni,nj) = dij(ni,nj) - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
2112  *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
2113  + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
2114  END DO
2115  END DO
2116  END DO
2117 
2118  !++++++++++++++++
2119  !=== In the following loops, ki=1 for Vr and ki=2 for Vz
2120  !=== ==> 2ki-1 = 1 for Vr, 2ki-1 = 3 for Vz (LA%loc_to_glob(2*ki-1,i))
2121  mat_loc_s = 0.d0
2122  idxn_s = 0
2123  jdxn_s = 0
2124  DO ki = 1, 2
2125  DO ni = 1, n_ws
2126  i = mesh%jjs(ni,ms)
2127  iglob = la%loc_to_glob(2*ki-1,i)
2128  ix = (ki-1)*n_ws+ni
2129  idxn_s(ix) = iglob-1
2130  DO kj = 1, 2
2131  DO nj = 1, n_ws
2132  j = mesh%jjs(nj,ms)
2133  jglob = la%loc_to_glob(2*kj-1,j)
2134  jx = (kj-1)*n_ws+nj
2135  jdxn_s(jx) = jglob-1
2136  IF ((ki == 1) .AND. (kj == 1)) THEN
2137  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
2138  ELSE IF ((ki == 1) .AND. (kj==2)) THEN
2139  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
2140  ELSE IF ((ki == 2) .AND. (kj==1)) THEN
2141  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
2142  ELSE
2143  mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
2144  END IF
2145  END DO
2146  END DO
2147  END DO
2148  END DO
2149  CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
2150  END DO
2151  END IF
2152 
2153  CALL matassemblybegin(matrix,mat_final_assembly,ierr)
2154  CALL matassemblyend(matrix,mat_final_assembly,ierr)
2155 
2156 
2157  END SUBROUTINE qs_diff_mass_vect_3x3_taylor
2158 END MODULE update_taylor_navier_stokes
subroutine update_ns_with_taylor_fourth(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:99
subroutine smb_cross_prod_gauss_sft_par(communicator, mesh, list_mode, V_in, V_out, precession_in)
subroutine update_ns_with_taylor(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:34
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:15
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
Definition: condlim.f90:191
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:147
subroutine, public fft_par_var_eta_prod_gauss_dcl(communicator, eta, H_mesh, c_in, c_out, nb_procs, bloc_size, m_max_pad, rr_gauss, time, temps)
subroutine error_petsc(string)
Definition: my_util.f90:16
type(my_data), public inputs
real(kind=8) function user_time()
Definition: my_util.f90:8
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:130
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
Definition: condlim.f90:162
subroutine qs_01_div_hybrid_generic(type_fe_velocity, vv_mesh, pp_mesh, pp_LA, mode, gg, pb_1, pb_2)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:12
subroutine dirichlet_rhs(js_D, bs_D, b)
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, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine qs_diff_mass_vect_3x3_taylor(type_op, LA, mesh, visco, mass, lambda, i_mode, mode, matrix)
Definition: tn_axi.f90:5
subroutine dirichlet_m_parallel(matrix, glob_js_D)
type(my_verbose), public talk_to_me
Definition: verbose.f90:27
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:40
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
subroutine rhs_ns_gauss_3x3_taylor(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, rotv_v, rhs_gauss, der_ord)
subroutine moy(communicator, mesh, p, RESLT)
subroutine, public navier_stokes_taylor(comm_one_d_ns, time, vv_3_LA, pp_1_LA, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un, vvz_per, pp_per)
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
subroutine, public init_velocity_pressure_taylor(vv_mesh, pp_mesh, list_mode, time, pn, der_pn, un, der_un)
subroutine scalar_without_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_mode_global_js_D)
subroutine inject_generic(type_fe_velocity, jj_c, jj_f, pp_c, pp_f)