6 #include "petsc/finclude/petsc.h" 14 list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un, vvz_per, pp_per)
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
26 REAL(KIND=8),
INTENT(IN) :: time
29 LOGICAL :: if_navier_stokes_with_taylor = .true.
30 #include "petsc/finclude/petsc.h" 31 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d_ns
33 IF (if_navier_stokes_with_taylor)
THEN 34 IF (
inputs%taylor_order==3)
THEN 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 39 inputs%dt,
inputs%Re,
inputs%taylor_lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
41 CALL error_petsc(
'BUG in navier_stokes_taylor: inputs%taylor_order not well defined')
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
56 REAL(KIND=8),
INTENT(IN) :: time
58 INTEGER :: i, mode, k, kp, m_max_c
61 m_max_c =
SIZE(list_mode)
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)
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)
92 IF (
inputs%taylor_order==4)
THEN 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)
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)
111 IF (vv_mesh%me/=0)
THEN 112 DO kp = 1,
inputs%taylor_order-1
114 IF (list_mode(i) == 0)
THEN 116 der_un(kp)%DRT(:,2*k,i) = 0.d0
118 der_pn(kp)%DRT(:,2,i) = 0.d0
126 dt, re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
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
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
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
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
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/
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
189 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
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)
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)
209 m_max_c =
SIZE(list_mode)
215 ALLOCATE(pp_global_d(m_max_c))
217 ALLOCATE(pp_global_d(i)%DRL(
SIZE(pp_mode_global_js_d(i)%DIL)))
223 vv_js_d, vv_mode_global_js_d)
225 ALLOCATE(vel_global_d(m_max_c))
227 ALLOCATE(vel_global_d(i)%DRL(
SIZE(vv_mode_global_js_d(i)%DIL)))
234 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 245 ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
252 lambda, i, mode, vel_mat(nu_mat))
254 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 256 vel_mat(nu_mat), vv_3_la)
259 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
264 lambda, i, mode, vel_mat(nu_mat))
266 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 268 vel_mat(nu_mat), vv_3_la)
271 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
283 uext = un + 2.d0*dt*der_un(1)%DRT + 2.d0*dt**2 *der_un(2)%DRT
286 uext = un + dt*der_un(1)%DRT + 0.5d0*dt**2*der_un(2)%DRT
291 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
295 IF (
inputs%verbose_CFL)
THEN 298 IF (list_mode(i)==0)
THEN 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)
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)
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)
314 CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
323 der_un(2)%DRT/dt, der_pn(2)%DRT, (rotv_v1-2.d0*rotv_v2+rotv_v3)/(dt**2), rhs_gauss,2)
331 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
338 n1 =
SIZE(vv_js_d(1)%DIL)
339 n2 =
SIZE(vv_js_d(2)%DIL)
340 n3 =
SIZE(vv_js_d(3)%DIL)
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)
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
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))
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
395 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
419 der_pn(2)%DRT(:,:,i) = der_pn(2)%DRT(:,:,i) - lambda*div(:,:,i)
420 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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
437 der_pn(2)%DRT (:,2,i) = 0.d0
442 dtt_un_p1(:,:,i) = un_p1
443 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
461 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
468 n1 =
SIZE(vv_js_d(1)%DIL)
469 n2 =
SIZE(vv_js_d(2)%DIL)
470 n3 =
SIZE(vv_js_d(3)%DIL)
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
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))
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
518 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
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
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
561 der_pn(1)%DRT (:,2,i) = 0.d0
566 dt_un_p1(:,:,i) = un_p1
567 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
585 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
592 n1 =
SIZE(vv_js_d(1)%DIL)
593 n2 =
SIZE(vv_js_d(2)%DIL)
594 n3 =
SIZE(vv_js_d(3)%DIL)
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
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))
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))
629 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
637 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
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
669 CALL moy(comm_one_d(1),pp_mesh, pn(:,1,i),moyenne)
670 pn(:,1,i) = pn(:,1,i)-moyenne
685 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
693 der_un(2)%DRT(:,:,i)=dtt_un_p1(:,:,i)
694 der_un(1)%DRT(:,:,i)=dt_un_p1(:,:,i)
697 IF (
inputs%verbose_divergence)
THEN 698 norm =
norm_sf(comm_one_d,
'H1', vv_mesh, list_mode, un)
708 dt, re, lambda, list_mode, pp_mesh, vv_mesh, pn, der_pn, un, der_un)
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
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
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
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
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/
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
771 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
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)
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)
791 m_max_c =
SIZE(list_mode)
796 ALLOCATE(pp_global_d(m_max_c))
798 ALLOCATE(pp_global_d(i)%DRL(
SIZE(pp_mode_global_js_d(i)%DIL)))
804 vv_js_d, vv_mode_global_js_d)
806 ALLOCATE(vel_global_d(m_max_c))
808 ALLOCATE(vel_global_d(i)%DRL(
SIZE(vv_mode_global_js_d(i)%DIL)))
815 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 826 ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
833 lambda, i, mode, vel_mat(nu_mat))
834 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 836 vel_mat(nu_mat), vv_3_la)
839 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
844 lambda, i, mode, vel_mat(nu_mat))
845 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 847 vel_mat(nu_mat), vv_3_la)
850 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
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
867 uext = un + dt*der_un(1)%DRT + dt**2/2*der_un(2)%DRT + dt**3/6*der_un(3)%DRT
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
878 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
882 IF (
inputs%verbose_CFL)
THEN 885 IF (list_mode(i)==0)
THEN 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)
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)
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)
901 CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
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)
918 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
925 n1 =
SIZE(vv_js_d(1)%DIL)
926 n2 =
SIZE(vv_js_d(2)%DIL)
927 n3 =
SIZE(vv_js_d(3)%DIL)
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
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))
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))
987 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
995 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
1019 der_pn(3)%DRT(:,:,i) = der_pn(3)%DRT(:,:,i) - lambda*div(:,:,i)
1020 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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
1037 der_pn(3)%DRT (:,2,i) = 0.d0
1042 dttt_un_p1(:,:,i) = un_p1
1043 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
1061 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
1068 n1 =
SIZE(vv_js_d(1)%DIL)
1069 n2 =
SIZE(vv_js_d(2)%DIL)
1070 n3 =
SIZE(vv_js_d(3)%DIL)
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)
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
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
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))
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))
1131 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
1139 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
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
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
1181 der_pn(2)%DRT (:,2,i) = 0.d0
1186 dtt_un_p1(:,:,i) = un_p1
1187 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
1205 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
1212 n1 =
SIZE(vv_js_d(1)%DIL)
1213 n2 =
SIZE(vv_js_d(2)%DIL)
1214 n3 =
SIZE(vv_js_d(3)%DIL)
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
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))
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))
1274 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
1282 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
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
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
1325 der_pn(1)%DRT (:,2,i) = 0.d0
1330 dt_un_p1(:,:,i) = un_p1
1331 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
1349 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
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)
1356 n1 =
SIZE(vv_js_d(1)%DIL)
1357 n2 =
SIZE(vv_js_d(2)%DIL)
1358 n3 =
SIZE(vv_js_d(3)%DIL)
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
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))
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))
1393 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
1401 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 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)
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))
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))
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
1433 CALL moy(comm_one_d(1),pp_mesh, pn(:,1,i),moyenne)
1434 pn(:,1,i) = pn(:,1,i)-moyenne
1449 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
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)
1462 IF (
inputs%verbose_divergence)
THEN 1463 norm =
norm_sf(comm_one_d,
'H1', vv_mesh, list_mode, un)
1465 talk_to_me%weak_div_L2 =
norm_sf(comm_one_d,
'L2', pp_mesh, list_mode, div)/norm
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
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.
1498 INTEGER :: nb_procs, bloc_size, m_max_pad, code
1499 mpi_comm :: communicator
1503 ALLOCATE(omega(mesh%np,6,
SIZE(list_mode)))
1506 DO i=1,
SIZE(list_mode)
1507 IF (list_mode(i) == 1)
THEN 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 1513 omega(:,5,i) =
inputs%taux_precession * cos(
inputs%angle_s_pi*pi)
1519 DO i = 1,
SIZE(list_mode)
1523 j_loc = mesh%jj(:,m)
1525 vs(:,k) = v_in(j_loc,k,i)
1526 omega_s(:,k) = omega(mesh%jj(:,m),k,i)
1529 DO l = 1, mesh%gauss%l_G
1531 dw_loc = mesh%gauss%dw(:,:,l,m)
1534 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
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))
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))
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))
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))
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)
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)
1574 IF (.NOT.precession_in)
THEN 1580 rotv = rotv + 2.d0 * om
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
1592 SUBROUTINE moy(communicator,mesh,p,RESLT)
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
1603 mpi_comm :: communicator
1607 DO m = 1, mesh%dom_me
1608 j_loc = mesh%jj(:,m)
1609 DO l = 1, mesh%gauss%l_G
1611 DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
1612 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
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)
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
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
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
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')
1664 DO i = 1,
SIZE(list_mode)
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 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 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 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)
1694 DO m = 1, vv_mesh%dom_me
1695 j_loc = vv_mesh%jj(:,m)
1697 DO l = 1, vv_mesh%gauss%l_G
1699 dw_loc = vv_mesh%gauss%dw(:,:,l,m)
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))
1707 fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l))
1708 ft(1) = sum(v1m(j_loc,1,i) * vv_mesh%gauss%ww(:,l))
1709 fp(index,1,i) = -sum(p(j_loc,1)*dw_loc(1,:))
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,:))
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)
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)
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,:))
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,:))
1731 rhs_gauss(index,:,i) = (ft+fs-rotv_v(index,:,i))
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)
1741 DO m = 1, vv_mesh%dom_me
1742 j_loc = vv_mesh%jj(:,m)
1743 DO l = 1, vv_mesh%gauss%l_G
1746 imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1747 * vv_mesh%gauss%ww(:,l))
1751 rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1752 ELSE IF (der_ord==1)
THEN 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)
1763 DO m = 1, vv_mesh%dom_me
1764 j_loc = vv_mesh%jj(:,m)
1765 DO l = 1, vv_mesh%gauss%l_G
1768 imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1769 * vv_mesh%gauss%ww(:,l))
1773 rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
1774 ELSE IF (der_ord==2)
THEN 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)
1783 DO m = 1, vv_mesh%dom_me
1784 j_loc = vv_mesh%jj(:,m)
1785 DO l = 1, vv_mesh%gauss%l_G
1788 imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1789 * vv_mesh%gauss%ww(:,l))
1793 rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
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)
1802 DO m = 1, vv_mesh%dom_me
1803 j_loc = vv_mesh%jj(:,m)
1804 DO l = 1, vv_mesh%gauss%l_G
1807 imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
1808 * vv_mesh%gauss%ww(:,l))
1812 rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
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
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(:,:)
1832 rhs_gauss(:,:,i) = rhs_gauss_penal(:,:,i)
1838 rhs_gauss = rhs_gauss + fp
1844 #include "petsc/finclude/petsc.h" 1850 INTEGER ,
INTENT(IN) :: type_op, mode, i_mode
1851 REAL(KIND=8),
INTENT(IN) :: visco, mass, lambda
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
1870 petscerrorcode :: ierr
1872 CALL matzeroentries (matrix, ierr)
1873 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1875 np =
SIZE(mesh%rr,2)
1876 n_w = mesh%gauss%n_w
1877 n_ws = mesh%gauss%n_ws
1879 IF (type_op == 1)
THEN 1883 ELSEIF (type_op == 2)
THEN 1887 ELSEIF (type_op == 3)
THEN 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)
1909 jj_loc = mesh%jj(:,m)
1914 DO l = 1, mesh%gauss%l_G
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)
1922 viscolm = visco*mesh%gauss%rj(l,m)
1925 hm=min(mesh%hm(i_mode),mesh%hloc(m))
1926 viscomode = visco*mesh%gauss%rj(l,m)
1928 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1930 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
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)
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
1945 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
1955 iglob = la%loc_to_glob(ki,i)
1961 jglob = la%loc_to_glob(kj,j)
1964 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN 1966 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1968 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1972 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1979 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1990 DO l = 1, mesh%gauss%l_G
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)
1998 viscolm = visco*mesh%gauss%rj(l,m)*ray
2000 div_penal = lambda*mesh%gauss%rj(l,m)*ray
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))
2035 iglob = la%loc_to_glob(ki,i)
2036 ix = (ki-1)*n_w + ni
2048 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
2051 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
2054 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
2060 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
2063 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
2066 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
2072 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
2075 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
2078 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
2081 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
2087 IF (
inputs%vv_nb_dirichlet_normal_velocity>0)
THEN 2090 IF (minval(abs(mesh%sides(ms)-
inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
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
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))
2127 iglob = la%loc_to_glob(2*ki-1,i)
2129 idxn_s(ix) = iglob-1
2133 jglob = la%loc_to_glob(2*kj-1,j)
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)
2143 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
2149 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
2153 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
2154 CALL matassemblyend(matrix,mat_final_assembly,ierr)
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)
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)
subroutine, public create_my_ghost(mesh, LA, ifrom)
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
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)
real(kind=8) function user_time()
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)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
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)
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)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
type(my_verbose), public talk_to_me
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
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)