7 #include "petsc/finclude/petsc.h" 43 mpi_comm :: fourier_communicator
44 mpi_comm,
DIMENSION(2) :: vizu_grad_curl_comm
48 vec :: vx_phi, vb_phi, vx_phi_ghost
55 list_mode_in, opt_nb_plane)
60 TYPE(
mesh_type),
TARGET,
INTENT(IN) :: vv_mesh_in
65 INTEGER,
DIMENSION(:),
TARGET :: list_mode_in
66 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_nb_plane
68 petscerrorcode :: ierr
69 petscmpiint :: rank, nb_procs, petsc_rank
70 mpi_comm,
DIMENSION(2) :: comm_one_d
71 CALL mpi_comm_rank(petsc_comm_world,petsc_rank,ierr)
72 CALL mpi_comm_rank(comm_one_d(2), rank, ierr)
73 CALL mpi_comm_size(comm_one_d(2), nb_procs, ierr)
74 fourier_communicator = comm_one_d(2)
80 IF ((.NOT.
inputs%if_just_processing) .AND. (
inputs%freq_plot .GT.
inputs%nb_iteration))
THEN 84 IF (
inputs%type_fe_velocity .GE. 3)
THEN 104 CALL mpi_comm_dup(comm_one_d(1),vizu_grad_curl_comm(1) , code)
105 CALL mpi_comm_dup(comm_one_d(2),vizu_grad_curl_comm(2) , code)
116 TYPE(
mesh_type),
INTENT(IN) :: local_mesh
118 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_nb_plane
119 INTEGER,
INTENT(IN) :: rank
120 LOGICAL,
DIMENSION(:),
POINTER :: virgin
121 INTEGER,
DIMENSION(:),
POINTER :: renumber
122 REAL(KIND=8) :: pi, dtheta, theta
123 INTEGER :: i, k, m, n, nb_angle, count, dim, np1, &
124 gap, m_max, m_max_pad
130 mesh_3d%gauss%n_w = 0
131 ALLOCATE(mesh_3d%jj(0,0))
132 ALLOCATE(mesh_3d%rr(0,0))
137 IF (
inputs%select_mode)
THEN 138 IF (
SIZE(
inputs%list_mode_lect)==1)
THEN 139 gap =
inputs%list_mode_lect(1) + 1
141 gap =
inputs%list_mode_lect(2) -
inputs%list_mode_lect(1)
150 IF (
PRESENT(opt_nb_plane))
THEN 151 IF (opt_nb_plane> 2*m_max-1)
THEN 152 m_max_pad = (opt_nb_plane+1)/2
159 nb_angle = 2*m_max_pad-1
165 dtheta = 2*pi/nb_angle
166 IF (mesh%gauss%n_w==3)
THEN 167 mesh_3d%np = local_mesh%np*nb_angle
168 mesh_3d%me = local_mesh%me*nb_angle
169 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
173 DO i = 1, local_mesh%np
175 mesh_3d%rr(1,count) = local_mesh%rr(1,i)*cos(theta)
176 mesh_3d%rr(2,count) = local_mesh%rr(1,i)*sin(theta)
177 mesh_3d%rr(3,count) = local_mesh%rr(2,i)
180 ALLOCATE(mesh_3d%jj(6,mesh_3d%me))
181 mesh_3d%gauss%n_w = 6
184 DO m = 1, local_mesh%me
186 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
187 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)+k*local_mesh%np
190 DO m = 1, local_mesh%me
192 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
193 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
195 ELSE IF (mesh%gauss%n_w==6)
THEN 196 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
200 DO m = 1, local_mesh%me
202 i = local_mesh%jj(n,m)
203 IF (.NOT.virgin(i)) cycle
210 mesh_3d%np = (local_mesh%np+np1)*nb_angle
211 mesh_3d%me = local_mesh%me*nb_angle
212 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
213 ALLOCATE(mesh_3d%jj(15,mesh_3d%me))
214 mesh_3d%gauss%n_w = 15
217 DO m = 1, local_mesh%me
219 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
220 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m) + k*(local_mesh%np+np1)
221 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
222 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
223 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
224 mesh_3d%jj(10,count) = local_mesh%jj(6,m) + k*(local_mesh%np+np1)
225 mesh_3d%jj(11,count) = local_mesh%jj(4,m) + k*(local_mesh%np+np1)
226 mesh_3d%jj(12,count) = local_mesh%jj(5,m) + k*(local_mesh%np+np1)
227 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
228 + (k-1)*(local_mesh%np+np1) + local_mesh%np
232 DO m = 1, local_mesh%me
234 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
235 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
236 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
237 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
238 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
239 mesh_3d%jj(10,count) = local_mesh%jj(6,m)
240 mesh_3d%jj(11,count) = local_mesh%jj(4,m)
241 mesh_3d%jj(12,count) = local_mesh%jj(5,m)
242 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
243 + (k-1)*(local_mesh%np+np1) + local_mesh%np
248 DO i = 1, local_mesh%np
249 n = (k-1)*(local_mesh%np+np1) + i
250 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(theta)
251 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(theta)
252 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
254 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
255 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(theta+dtheta/2)
256 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(theta+dtheta/2)
257 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
261 DEALLOCATE(virgin,renumber)
262 ELSE IF (mesh%gauss%n_w==10)
THEN 264 WRITE(*,*)
'VTK files not done for P3 meshes' 267 CALL error_petsc(
'Bug in prepare_3d_grids: mesh%gauss%n_w= not correct')
271 SUBROUTINE vtu_3d(field, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl, opt_2D, opt_mesh_in)
280 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN),
TARGET :: field
281 REAL(KIND=8),
DIMENSION(:,:,:),
POINTER :: field_in
282 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE,
TARGET :: field_tmp
283 CHARACTER(*),
INTENT(IN) :: name_of_mesh, header, name_of_field, what
284 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
285 CHARACTER(*),
OPTIONAL,
INTENT(IN) :: opt_grad_curl
286 LOGICAL,
OPTIONAL,
INTENT(IN) :: opt_2D
287 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_viz
289 CHARACTER(LEN=3) :: st_mode
290 CHARACTER(LEN=200) :: header_2D
291 CHARACTER(LEN=3) :: name_of_field_2D
295 IF (
PRESENT(opt_grad_curl))
THEN 296 IF ((opt_grad_curl ==
'curl_u').AND.(
PRESENT(opt_mesh_in)))
THEN 297 IF (opt_mesh_in%gauss%n_w==10 .AND.
inputs%type_fe_velocity==3)
THEN 298 ALLOCATE(field_viz(
SIZE(field,1),
SIZE(field,2),
SIZE(field,3)))
306 IF (
PRESENT(opt_mesh_in))
THEN 307 IF (opt_mesh_in%gauss%n_w==10 .AND.
inputs%type_fe_velocity==3)
THEN 309 DO i = 1,
SIZE(field,3)
310 IF (.NOT.(
PRESENT(opt_grad_curl)))
THEN 312 ELSE IF (opt_grad_curl ==
'curl_u')
THEN 315 CALL error_petsc(
'Bug in vtu_3d: P3 and wrong opt_grad_curl')
318 IF (
ALLOCATED(field_viz))
DEALLOCATE(field_viz)
319 field_in => field_tmp
324 IF (name_of_mesh ==
'vv_mesh' .AND.
inputs%type_fe_velocity==3)
THEN 325 CALL error_petsc(
'Bug in vtu_3d: P3 for velocity BUT opt_mesh_in not present')
332 IF (
PRESENT(opt_grad_curl))
THEN 333 IF (opt_grad_curl ==
'grad')
THEN 334 ALLOCATE(field_viz(
SIZE(field_in,1), 3*
SIZE(field_in,2),
SIZE(field_in,3)))
336 ELSE IF (opt_grad_curl ==
'curl_h')
THEN 337 ALLOCATE(field_viz(
SIZE(field_in,1),
SIZE(field_in,2),
SIZE(field_in,3)))
339 ELSE IF ((opt_grad_curl ==
'curl_u').AND.(
PRESENT(opt_mesh_in)))
THEN 340 IF (opt_mesh_in%gauss%n_w==10 .AND.
inputs%type_fe_velocity==3)
THEN 341 ALLOCATE(field_viz(
SIZE(field_in,1),
SIZE(field_in,2),
SIZE(field_in,3)))
344 ALLOCATE(field_viz(
SIZE(field_in,1),
SIZE(field_in,2),
SIZE(field_in,3)))
348 CALL error_petsc(
'Bug in vtu_3d: name_of_opt_grad_curl is not correct')
351 ALLOCATE(field_viz(
SIZE(field_in, 1),
SIZE(field_in,2),
SIZE(field_in,3)))
355 IF (
PRESENT(opt_it))
THEN 356 CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what, opt_it)
358 CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what)
361 IF (
PRESENT(opt_2d))
THEN 365 header_2d =
'J_'//
'mode_'//trim(adjustl(st_mode))
366 name_of_field_2d =
'J' 367 IF (
PRESENT(opt_it))
THEN 369 header_2d, field_viz(:,:,i), name_of_field_2d, what, opt_it=opt_it)
372 header_2d, field_viz(:,:,i), name_of_field_2d, what)
378 DEALLOCATE(field_viz)
379 IF (
ALLOCATED(field_tmp))
DEALLOCATE(field_tmp)
382 SUBROUTINE vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
390 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
391 CHARACTER(*),
INTENT(IN) :: name_of_mesh, header, name_of_field, what
392 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
395 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_FFT
396 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: field_out
397 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_xyz_FFT
398 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: dist_field
399 LOGICAL,
DIMENSION(:),
POINTER :: virgin
400 INTEGER,
DIMENSION(:),
POINTER :: renumber
401 INTEGER,
DIMENSION(:),
POINTER :: list_mode
402 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER :: loc_to_glob
403 INTEGER :: width, mode_min
404 INTEGER,
DIMENSION(1) :: loc
405 REAL(KIND=8) :: pi, dtheta, theta
406 INTEGER :: i, k, m, n, np_b, np_e, np_loc, np_loc_max, &
407 type, nb_angle, count, dim, np1, nb_procs, rank, it
415 ALLOCATE(loc_to_glob(nb_procs))
416 IF (name_of_mesh==
'vv_mesh')
THEN 422 ELSE IF (name_of_mesh==
'pp_mesh')
THEN 428 ELSE IF (name_of_mesh==
'H_mesh')
THEN 434 ELSE IF (name_of_mesh==
'phi_mesh')
THEN 440 ELSE IF (name_of_mesh==
'temp_mesh')
THEN 447 CALL error_petsc(
'Bug in vtu_3d: name_of_mesh is not correct')
453 np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
455 ALLOCATE(dist_field(nb_procs*np_loc_max,
SIZE(field_in,2),width))
460 mode_min = rank*width
462 IF (minval(abs(list_mode - (mode_min+i-1)))/=0) cycle
463 loc = minloc(abs(list_mode - (mode_min+i-1)))
464 DO type = 1, size(field_in,2)
466 np_b = (n-1)*np_loc_max + 1
467 np_loc =
SIZE(loc_to_glob(n)%DIL)
468 np_e = np_b + np_loc - 1
469 dist_field(np_b:np_e,
type,i) = field_in(loc_to_glob(n)%DIL,
type,loc(1))
474 CALL fft_par_real(fourier_communicator, dist_field, field_out_fft, opt_nb_plane=nb_angle)
477 ALLOCATE(field_out_xyz_fft(
SIZE(field_out_fft,1),
SIZE(field_out_fft,2),
SIZE(field_out_fft,3)))
478 IF (
SIZE(field_in,2)==2)
THEN 479 field_out_xyz_fft=field_out_fft
482 dtheta = 2*pi/nb_angle
485 field_out_xyz_fft(k,1,:) = field_out_fft(k,1,:)*cos(theta)-field_out_fft(k,2,:)*sin(theta)
486 field_out_xyz_fft(k,2,:) = field_out_fft(k,1,:)*sin(theta)+field_out_fft(k,2,:)*cos(theta)
487 field_out_xyz_fft(k,3,:) = field_out_fft(k,3,:)
492 IF (
SIZE(field_in,2)==2)
THEN 494 ELSE IF (
SIZE(field_in,2)==6)
THEN 497 CALL error_petsc(
'Bug in vtu_3d: SIZE(field_in,2) not correct')
499 IF (nb_angle /=
SIZE(field_out_xyz_fft,1))
THEN 500 CALL error_petsc(
'Bug in vtu_3d: nb_angle /= SIZE(field_out_xyz_FFT,1')
505 dtheta = 2*pi/nb_angle
506 IF (local_mesh%gauss%n_w==3)
THEN 507 ALLOCATE(field_out(dim,mesh_3d%np))
510 DO i = 1, local_mesh%np
512 field_out(:,count) = field_out_xyz_fft(k,:,i)
516 ELSE IF (local_mesh%gauss%n_w==6)
THEN 517 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
521 DO m = 1, local_mesh%me
523 i = local_mesh%jj(n,m)
524 IF (.NOT.virgin(i)) cycle
532 ALLOCATE(field_out(dim,mesh_3d%np))
534 DO i = 1, local_mesh%np
535 n = (k-1)*(local_mesh%np+np1) + i
536 field_out(:,n) = field_out_xyz_fft(k,:,i)
540 DO i = 1, local_mesh%np
542 DO k = 1, nb_angle - 2
543 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
544 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
545 + (3.d0/4)*field_out_xyz_fft(k+1,:,i) - (1.d0/8)*field_out_xyz_fft(k+2,:,i)
549 DO i = 1, local_mesh%np
551 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
552 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
553 + (3.d0/4)*field_out_xyz_fft(k+1,:,i) - (1.d0/8)*field_out_xyz_fft(1,:,i)
556 DO i = 1, local_mesh%np
558 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
559 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
560 + (3.d0/4)*field_out_xyz_fft(1,:,i) - (1.d0/8)*field_out_xyz_fft(2,:,i)
562 DEALLOCATE(virgin, renumber)
564 CALL error_petsc(
'Bug in vtu_3d: mesh%gauss%n_w is not correct')
567 IF (
PRESENT(opt_it))
THEN 570 field_out, name_of_field, what, opt_it = it)
573 field_out, name_of_field, what)
577 DEALLOCATE(dist_field, field_out, loc_to_glob,field_out_xyz_fft)
578 DEALLOCATE(field_out_fft)
582 header, name_of_field, what, list_mode, opt_it)
591 CHARACTER(*),
INTENT(IN) :: header, name_of_field, what
592 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
593 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
594 REAL(KIND=8),
DIMENSION(:,:,:) :: field_in
596 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_FFT
597 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: field_out
598 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: dist_field
599 LOGICAL,
DIMENSION(:),
POINTER :: virgin
600 INTEGER,
DIMENSION(:),
POINTER :: renumber
602 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER :: loc_to_glob
603 INTEGER :: gap, width, mode_min
604 INTEGER,
DIMENSION(1) :: loc
605 REAL(KIND=8) :: pi, dtheta, theta
606 INTEGER :: i, k, m, n, np_b, np_e, np_loc, np_loc_max, &
607 type, nb_angle, count, dim, np1
608 petscerrorcode :: ierr
609 petscmpiint :: rank, nb_procs
610 mpi_comm :: communicator
611 CALL mpi_comm_rank(communicator, rank, ierr)
612 CALL mpi_comm_size(communicator, nb_procs, ierr)
614 IF (
inputs%select_mode)
THEN 615 IF (
SIZE(
inputs%list_mode_lect)==1)
THEN 616 gap =
inputs%list_mode_lect(1) + 1
618 gap =
inputs%list_mode_lect(2) -
inputs%list_mode_lect(1)
623 width =
SIZE(list_mode)*gap
625 CALL divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
628 np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
630 ALLOCATE(dist_field(nb_procs*np_loc_max,
SIZE(field_in,2),width))
632 mode_min = rank*width
634 IF (minval(abs(list_mode - (mode_min+i-1)))/=0) cycle
635 loc = minloc(abs(list_mode - (mode_min+i-1)))
636 DO type = 1, size(field_in,2)
638 np_b = (n-1)*np_loc_max + 1
639 np_loc =
SIZE(loc_to_glob(n)%DIL)
640 np_e = np_b + np_loc - 1
641 dist_field(np_b:np_e,
type,i) = field_in(loc_to_glob(n)%DIL,
type,loc(1))
646 CALL fft_par_real(communicator, dist_field, field_out_fft, opt_nb_plane=50)
648 IF (
SIZE(field_in,2)==2)
THEN 650 ELSE IF (
SIZE(field_in,2)==6)
THEN 653 CALL error_petsc(
'Bug in transfer_fourier_to_real: SIZE(field_in,2) not correct')
655 nb_angle =
SIZE(field_out_fft,1)
657 dtheta = 2*pi/nb_angle
658 IF (mesh%gauss%n_w==3)
THEN 659 mesh_3d%np = local_mesh%np*nb_angle
660 mesh_3d%me = local_mesh%me*nb_angle
661 ALLOCATE(field_out(dim,mesh_3d%np))
664 DO i = 1, local_mesh%np
666 field_out(:,count) = field_out_fft(k,:,i)
670 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
674 DO i = 1, local_mesh%np
676 mesh_3d%rr(1,count) = local_mesh%rr(1,i)*cos(theta)
677 mesh_3d%rr(2,count) = local_mesh%rr(1,i)*sin(theta)
678 mesh_3d%rr(3,count) = local_mesh%rr(2,i)
681 ALLOCATE(mesh_3d%jj(6,mesh_3d%me))
682 mesh_3d%gauss%n_w = 6
685 DO m = 1, local_mesh%me
687 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
688 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)+k*local_mesh%np
691 DO m = 1, local_mesh%me
693 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
694 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
696 ELSE IF (mesh%gauss%n_w==6)
THEN 697 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
701 DO m = 1, local_mesh%me
703 i = local_mesh%jj(n,m)
704 IF (.NOT.virgin(i)) cycle
711 mesh_3d%np = (local_mesh%np+np1)*nb_angle
712 mesh_3d%me = local_mesh%me*nb_angle
713 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
714 ALLOCATE(mesh_3d%jj(15,mesh_3d%me))
715 mesh_3d%gauss%n_w = 15
718 DO m = 1, local_mesh%me
720 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
721 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m) + k*(local_mesh%np+np1)
722 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
723 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
724 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
725 mesh_3d%jj(10,count) = local_mesh%jj(6,m) + k*(local_mesh%np+np1)
726 mesh_3d%jj(11,count) = local_mesh%jj(4,m) + k*(local_mesh%np+np1)
727 mesh_3d%jj(12,count) = local_mesh%jj(5,m) + k*(local_mesh%np+np1)
728 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
729 + (k-1)*(local_mesh%np+np1) + local_mesh%np
733 DO m = 1, local_mesh%me
735 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
736 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
737 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
738 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
739 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
740 mesh_3d%jj(10,count) = local_mesh%jj(6,m)
741 mesh_3d%jj(11,count) = local_mesh%jj(4,m)
742 mesh_3d%jj(12,count) = local_mesh%jj(5,m)
743 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
744 + (k-1)*(local_mesh%np+np1) + local_mesh%np
749 DO i = 1, local_mesh%np
750 n = (k-1)*(local_mesh%np+np1) + i
751 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(theta)
752 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(theta)
753 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
755 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
756 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(theta+dtheta/2)
757 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(theta+dtheta/2)
758 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
762 ALLOCATE(field_out(dim,mesh_3d%np))
764 DO i = 1, local_mesh%np
765 n = (k-1)*(local_mesh%np+np1) + i
766 field_out(:,n) = field_out_fft(k,:,i)
770 DO i = 1, local_mesh%np
772 DO k = 1, nb_angle - 2
773 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
774 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
775 + (3.d0/4)*field_out_fft(k+1,:,i) - (1.d0/8)*field_out_fft(k+2,:,i)
779 DO i = 1, local_mesh%np
781 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
782 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
783 + (3.d0/4)*field_out_fft(k+1,:,i) - (1.d0/8)*field_out_fft(1,:,i)
786 DO i = 1, local_mesh%np
788 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
789 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
790 + (3.d0/4)*field_out_fft(1,:,i) - (1.d0/8)*field_out_fft(2,:,i)
792 DEALLOCATE(virgin, renumber)
795 CALL error_petsc(
'Bug in transfer_fourier_to_real: mesh%gauss%n_w is not correct')
798 IF (
PRESENT(opt_it))
THEN 800 field_out, name_of_field, what, opt_it=opt_it)
803 field_out, name_of_field, what)
808 IF (
ASSOCIATED(loc_to_glob(n)%DIL))
DEALLOCATE(loc_to_glob(n)%DIL)
810 DEALLOCATE(loc_to_glob)
811 DEALLOCATE(mesh_3d%rr, mesh_3d%jj)
812 DEALLOCATE(local_mesh%rr, local_mesh%jj)
813 DEALLOCATE(dist_field, field_out)
814 DEALLOCATE(field_out_fft)
817 SUBROUTINE divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
822 TYPE(
mesh_type),
INTENT(OUT) :: local_mesh
823 INTEGER :: me_b, me_e, me_loc, np_loc
824 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER :: loc_to_glob
825 INTEGER,
POINTER,
DIMENSION(:) :: glob_to_loc
826 LOGICAL,
DIMENSION(mesh%np) :: virgin
827 INTEGER :: count, m, m_glob, n, i, q, r, loop_rank
828 petscerrorcode :: ierr
829 petscmpiint :: rank, nb_procs
830 mpi_comm :: communicator
836 local_mesh%gauss%n_w = 0
837 ALLOCATE(local_mesh%jj(0,0))
838 ALLOCATE(local_mesh%rr(0,0))
839 ALLOCATE(loc_to_glob(0))
844 CALL mpi_comm_rank(communicator,rank,ierr)
845 CALL mpi_comm_size(communicator,nb_procs,ierr)
846 ALLOCATE(loc_to_glob(nb_procs))
847 ALLOCATE(glob_to_loc(mesh%np))
851 r = mesh%me - q*nb_procs
853 DO loop_rank = 0, nb_procs-1
855 IF (loop_rank+1.LE.r)
THEN 856 me_b = loop_rank*(q+1) + 1
859 me_b = loop_rank*q + r + 1
862 me_loc = me_e - me_b + 1
867 DO n = 1, mesh%gauss%n_w
869 IF (.NOT.virgin(i)) cycle
875 ALLOCATE(loc_to_glob(loop_rank+1)%DIL(np_loc))
879 DO n = 1, mesh%gauss%n_w
881 IF (.NOT.virgin(i)) cycle
884 loc_to_glob(loop_rank+1)%DIL(count) = i
885 IF (loop_rank == rank) glob_to_loc(i) = count
889 IF (loop_rank == rank)
THEN 891 local_mesh%me = me_loc
892 local_mesh%np = np_loc
893 local_mesh%gauss%n_w = mesh%gauss%n_w
894 ALLOCATE(local_mesh%jj(mesh%gauss%n_w,me_loc))
897 local_mesh%jj(:,m) = glob_to_loc(mesh%jj(:,m_glob))
899 ALLOCATE(local_mesh%rr(2,local_mesh%np))
900 local_mesh%rr = mesh%rr(:,loc_to_glob(rank+1)%DIL)
904 DEALLOCATE(glob_to_loc)
914 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
915 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: field_out
917 INTEGER,
POINTER,
DIMENSION(:) :: ifrom
919 REAL(KIND=8),
DIMENSION(SIZE(field_out,1), SIZE(field_out,2)) :: smb_grad_phi
920 REAL(KIND=8),
DIMENSION(6) :: grad_phi_loc
921 INTEGER,
DIMENSION(mesh_grad_phi%gauss%n_w) :: i_loc
922 INTEGER :: i, mode, m, l, j_loc, k, nj
923 REAL(KIND=8) :: ray, drdthetadz
924 petscerrorcode :: ierr
927 param_grad_phi%it_max=1
928 param_grad_phi%rel_tol=1.d-6
929 param_grad_phi%abs_tol=1.d-10
930 param_grad_phi%solver=
'GMRES' 931 param_grad_phi%precond=
'MUMPS' 932 param_grad_phi%verbose=.false.
938 CALL veccreateghost(vizu_grad_curl_comm(1),
mesh_grad_phi%dom_np, &
939 petsc_determine,
SIZE(ifrom), ifrom, vx_phi, ierr)
940 CALL vecghostgetlocalform(vx_phi, vx_phi_ghost, ierr)
941 CALL vecduplicate(vx_phi, vb_phi, ierr)
945 CALL init_solver(param_grad_phi,ksp_grad_phi,mat_grad_phi,vizu_grad_curl_comm(1),&
946 solver=
'GMRES',precond=
'MUMPS')
947 CALL matdestroy(mat_grad_phi,ierr)
963 grad_phi_loc(1) = grad_phi_loc(1) + field_in(j_loc,1,i)*
mesh_grad_phi%gauss%dw(1,nj,l,m)*ray
964 grad_phi_loc(2) = grad_phi_loc(2) + field_in(j_loc,2,i)*
mesh_grad_phi%gauss%dw(1,nj,l,m)*ray
965 grad_phi_loc(3) = grad_phi_loc(3) + field_in(j_loc,2,i)*
mesh_grad_phi%gauss%ww(nj,l)*mode
966 grad_phi_loc(4) = grad_phi_loc(4) - field_in(j_loc,1,i)*
mesh_grad_phi%gauss%ww(nj,l)*mode
967 grad_phi_loc(5) = grad_phi_loc(5) + field_in(j_loc,1,i)*
mesh_grad_phi%gauss%dw(2,nj,l,m)*ray
968 grad_phi_loc(6) = grad_phi_loc(6) + field_in(j_loc,2,i)*
mesh_grad_phi%gauss%dw(2,nj,l,m)*ray
973 smb_grad_phi(i_loc,k) = smb_grad_phi(i_loc,k) + grad_phi_loc(k)* &
980 CALL veczeroentries(vb_phi, ierr)
982 CALL vecassemblybegin(vb_phi,ierr)
983 CALL vecassemblyend(vb_phi,ierr)
985 CALL solver(ksp_grad_phi,vb_phi,vx_phi,reinit=.false.,
verbose=.false.)
987 CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
988 CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
995 field_out(:,2,i) = 0.d0
996 field_out(:,4,i) = 0.d0
997 field_out(:,6,i) = 0.d0
1011 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
1012 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: field_out
1014 INTEGER,
POINTER,
DIMENSION(:) :: ifrom
1016 REAL(KIND=8),
DIMENSION(SIZE(field_out,1), SIZE(field_out,2)) :: smb_rot_h
1017 REAL(KIND=8),
DIMENSION(6) :: rot_h_loc
1018 INTEGER,
DIMENSION(mesh_rot_h%gauss%n_w) :: i_loc
1019 INTEGER :: i, mode, m, l, j_loc, k, nj
1020 REAL(KIND=8) :: ray, drdthetadz
1021 petscerrorcode :: ierr
1023 param_rot_h%it_max=1
1024 param_rot_h%rel_tol=1.d-6
1025 param_rot_h%abs_tol=1.d-10
1026 param_rot_h%solver=
'GMRES' 1027 param_rot_h%precond=
'MUMPS' 1028 param_rot_h%verbose=.false.
1034 CALL veccreateghost(vizu_grad_curl_comm(1),
mesh_rot_h%dom_np, &
1035 petsc_determine,
SIZE(ifrom), ifrom, vx_phi, ierr)
1036 CALL vecghostgetlocalform(vx_phi, vx_phi_ghost, ierr)
1037 CALL vecduplicate(vx_phi, vb_phi, ierr)
1041 CALL init_solver(param_rot_h,ksp_rot_h,mat_rot_h,vizu_grad_curl_comm(1),&
1042 solver=
'GMRES',precond=
'MUMPS')
1043 CALL matdestroy(mat_rot_h,ierr)
1059 rot_h_loc(1) = rot_h_loc(1) + mode*field_in(j_loc,6,i)*
mesh_rot_h%gauss%ww(nj,l) &
1060 -field_in(j_loc,3,i)*
mesh_rot_h%gauss%dw(2,nj,l,m)*ray
1061 rot_h_loc(2) = rot_h_loc(2) - mode*field_in(j_loc,5,i)*
mesh_rot_h%gauss%ww(nj,l) &
1062 -field_in(j_loc,4,i)*
mesh_rot_h%gauss%dw(2,nj,l,m)*ray
1063 rot_h_loc(3) = rot_h_loc(3) + field_in(j_loc,1,i)*
mesh_rot_h%gauss%dw(2,nj,l,m)*ray &
1064 -field_in(j_loc,5,i)*
mesh_rot_h%gauss%dw(1,nj,l,m)*ray
1065 rot_h_loc(4) = rot_h_loc(4) + field_in(j_loc,2,i)*
mesh_rot_h%gauss%dw(2,nj,l,m)*ray &
1066 -field_in(j_loc,6,i)*
mesh_rot_h%gauss%dw(1,nj,l,m)*ray
1067 rot_h_loc(5) = rot_h_loc(5) + field_in(j_loc,3,i)*
mesh_rot_h%gauss%ww(nj,l) &
1068 +field_in(j_loc,3,i)*
mesh_rot_h%gauss%dw(1,nj,l,m)*ray &
1069 - mode*field_in(j_loc,2,i)*
mesh_rot_h%gauss%ww(nj,l)
1071 rot_h_loc(6) = rot_h_loc(6) + field_in(j_loc,4,i)*
mesh_rot_h%gauss%ww(nj,l) &
1072 +field_in(j_loc,4,i)*
mesh_rot_h%gauss%dw(1,nj,l,m)*ray &
1073 + mode*field_in(j_loc,1,i)*
mesh_rot_h%gauss%ww(nj,l)
1078 smb_rot_h(i_loc,k) = smb_rot_h(i_loc,k) + rot_h_loc(k)* &
1085 CALL veczeroentries(vb_phi, ierr)
1087 CALL vecassemblybegin(vb_phi,ierr)
1088 CALL vecassemblyend(vb_phi,ierr)
1090 CALL solver(ksp_rot_h,vb_phi,vx_phi,reinit=.false.,
verbose=.false.)
1092 CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
1093 CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
1100 field_out(:,2,i) = 0.d0
1101 field_out(:,4,i) = 0.d0
1102 field_out(:,6,i) = 0.d0
1116 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
1117 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: field_out
1120 INTEGER,
POINTER,
DIMENSION(:) :: ifrom
1122 REAL(KIND=8),
DIMENSION(SIZE(field_out,1), SIZE(field_out,2)) :: smb_rot_u
1123 REAL(KIND=8),
DIMENSION(6) :: rot_u_loc
1124 INTEGER,
DIMENSION(mesh_in%gauss%n_w) :: i_loc
1125 INTEGER :: i, mode, m, l, j_loc, k, nj
1126 REAL(KIND=8) :: ray, drdthetadz
1127 petscerrorcode :: ierr
1129 param_rot_u%it_max=1
1130 param_rot_u%rel_tol=1.d-6
1131 param_rot_u%abs_tol=1.d-10
1132 param_rot_u%solver=
'GMRES' 1133 param_rot_u%precond=
'MUMPS' 1134 param_rot_u%verbose=.false.
1140 CALL veccreateghost(vizu_grad_curl_comm(1), mesh_in%dom_np, &
1141 petsc_determine,
SIZE(ifrom), ifrom, vx_phi, ierr)
1142 CALL vecghostgetlocalform(vx_phi, vx_phi_ghost, ierr)
1143 CALL vecduplicate(vx_phi, vb_phi, ierr)
1147 CALL init_solver(param_rot_u,ksp_rot_u,mat_rot_u,vizu_grad_curl_comm(1),&
1148 solver=
'GMRES',precond=
'MUMPS')
1149 CALL matdestroy(mat_rot_u,ierr)
1157 DO m = 1, mesh_in%me
1158 DO l = 1, mesh_in%gauss%l_G
1159 ray = sum(mesh_in%rr(1,mesh_in%jj(:,m))*mesh_in%gauss%ww(:,l))
1160 drdthetadz = mesh_in%gauss%rj(l,m)
1162 DO nj = 1,mesh_in%gauss%n_w
1163 j_loc = mesh_in%jj(nj,m)
1165 rot_u_loc(1) = rot_u_loc(1) + mode*field_in(j_loc,6,i)*mesh_in%gauss%ww(nj,l) &
1166 -field_in(j_loc,3,i)*mesh_in%gauss%dw(2,nj,l,m)*ray
1167 rot_u_loc(2) = rot_u_loc(2) - mode*field_in(j_loc,5,i)*mesh_in%gauss%ww(nj,l) &
1168 -field_in(j_loc,4,i)*mesh_in%gauss%dw(2,nj,l,m)*ray
1170 rot_u_loc(3) = rot_u_loc(3) + field_in(j_loc,1,i)*mesh_in%gauss%dw(2,nj,l,m)*ray &
1171 -field_in(j_loc,5,i)*mesh_in%gauss%dw(1,nj,l,m)*ray
1172 rot_u_loc(4) = rot_u_loc(4) + field_in(j_loc,2,i)*mesh_in%gauss%dw(2,nj,l,m)*ray &
1173 -field_in(j_loc,6,i)*mesh_in%gauss%dw(1,nj,l,m)*ray
1175 rot_u_loc(5) = rot_u_loc(5) + field_in(j_loc,3,i)*mesh_in%gauss%ww(nj,l) &
1176 +field_in(j_loc,3,i)*mesh_in%gauss%dw(1,nj,l,m)*ray &
1177 - mode*field_in(j_loc,2,i)*mesh_in%gauss%ww(nj,l)
1178 rot_u_loc(6) = rot_u_loc(6) + field_in(j_loc,4,i)*mesh_in%gauss%ww(nj,l) &
1179 +field_in(j_loc,4,i)*mesh_in%gauss%dw(1,nj,l,m)*ray &
1180 + mode*field_in(j_loc,1,i)*mesh_in%gauss%ww(nj,l)
1182 i_loc = mesh_in%jj(:,m)
1185 smb_rot_u(i_loc,k) = smb_rot_u(i_loc,k) + rot_u_loc(k)* &
1186 mesh_in%gauss%ww(:,l)*drdthetadz
1192 CALL veczeroentries(vb_phi, ierr)
1193 CALL vecsetvalues(vb_phi, mesh_in%np,
vizu_rot_u_la%loc_to_glob(1,:)-1, smb_rot_u(:,k), add_values, ierr)
1194 CALL vecassemblybegin(vb_phi,ierr)
1195 CALL vecassemblyend(vb_phi,ierr)
1197 CALL solver(ksp_rot_u,vb_phi,vx_phi,reinit=.false.,
verbose=.false.)
1199 CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
1200 CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
1201 IF (mesh_in%me/=0)
THEN 1207 field_out(:,2,i) = 0.d0
1208 field_out(:,4,i) = 0.d0
1209 field_out(:,6,i) = 0.d0
1268 TYPE(
mesh_type),
INTENT(OUT) :: mesh_out
1269 LOGICAL,
DIMENSION(mesh_in%np) :: v_virgin
1270 LOGICAL,
DIMENSION(mesh_in%me) :: c_virgin
1271 INTEGER,
DIMENSION(mesh_in%np) :: v_in_to_out
1272 INTEGER :: edge, nb_vertices, nb_edges, n_dof, &
1273 m, m_out, j, m_op, m_out_op, n, n1, n2, n_op, n1_op, n2_op, nn, &
1275 mesh_out%gauss%k_d = 2
1276 mesh_out%gauss%n_w = 6
1277 mesh_out%me = 4*mesh_in%me
1281 DO m = 1, mesh_in%me
1283 IF (mesh_in%neigh(n,m).LE.0) cycle
1288 DO m = 1, mesh_in%me
1290 IF (mesh_in%neigh(n,m).GE.0) cycle
1294 nb_edges = edge + mesh_in%mes
1296 IF (mesh_in%gauss%n_w==10)
THEN 1297 nb_vertices = mesh_in%np - 2*nb_edges - mesh_in%me
1299 CALL error_petsc(
'BUG in divide_mesh_into_four_subcells: mesh to be divided is not P3')
1301 mesh_out%np = nb_vertices + 3*nb_edges + 3*mesh_in%me
1304 ALLOCATE(mesh_out%jj(6,mesh_out%me), mesh_out%rr(2,mesh_out%np))
1307 DO m = 1, mesh_in%me
1310 IF (v_virgin(j))
THEN 1311 v_virgin(j) = .false.
1313 v_in_to_out(j) = n_dof
1318 IF (n_dof.NE.nb_vertices)
THEN 1319 write(*,*) mesh_in%np, nb_edges, mesh_in%me, nb_vertices, n_dof
1320 CALL error_petsc(.NE.
'BUG divide_mesh_into_four_subcells: n_dofnb_vertices')
1323 DO m = 1, mesh_in%me
1327 mesh_out%jj(n,m_out) = v_in_to_out(j)
1328 mesh_out%rr(:,v_in_to_out(j)) = mesh_in%rr(:,j)
1334 DO m = 1, mesh_in%me
1337 m_op = mesh_in%neigh(n,m)
1339 IF (.NOT.c_virgin(m_op)) cycle
1346 IF (mesh_in%neigh(n_op,m_op)==m)
EXIT 1349 n1_op = modulo(n_op,3) + 1
1350 n2_op = modulo(n_op+1,3) + 1
1351 n1 = modulo(n,3) + 1
1352 n2 = modulo(n+1,3) + 1
1353 IF (mesh_in%jj(n1_op,m_op) .NE. mesh_in%jj(n1,m))
THEN 1358 IF (mesh_in%jj(n1,m).NE.mesh_in%jj(n1_op,m_op))
THEN 1359 CALL error_petsc(
'divide_mesh_into_four_subcells: BUG')
1365 mesh_out%rr(:,n_dof) = (mesh_in%rr(:,j1)+mesh_in%rr(:,j2))/2
1368 m_out_op = 4*(m_op-1)+n1_op
1369 mesh_out%jj(n2,m_out) = n_dof
1370 mesh_out%jj(n2_op,m_out_op) = n_dof
1373 m_out_op = 4*(m_op-1)+n2_op
1374 mesh_out%jj(n1,m_out) = n_dof
1375 mesh_out%jj(n1_op,m_out_op) = n_dof
1379 mesh_out%jj(n,m_out) = n_dof
1380 mesh_out%jj(n_op,m_out_op) = n_dof
1383 mesh_out%rr(:,n_dof) = mesh_in%rr(:,j1)*3/4 + mesh_in%rr(:,j2)/4
1386 m_out_op = 4*(m_op-1)+n1_op
1387 mesh_out%jj(n+3,m_out) = n_dof
1388 mesh_out%jj(n_op+3,m_out_op) = n_dof
1391 mesh_out%rr(:,n_dof) = mesh_in%rr(:,j1)/4 + mesh_in%rr(:,j2)*3/4
1394 m_out_op = 4*(m_op-1)+n2_op
1395 mesh_out%jj(n+3,m_out) = n_dof
1396 mesh_out%jj(n_op+3,m_out_op) = n_dof
1402 n1 = modulo(n,3) + 1
1403 n2 = modulo(n+1,3) + 1
1408 mesh_out%rr(:,n_dof) = mesh_in%rr(:,j)/2 + mesh_in%rr(:,j1)/4 + mesh_in%rr(:,j2)/4
1411 mesh_out%jj(n+3,m_out) = n_dof
1413 mesh_out%jj(n+3,m_out) = n_dof
1419 SUBROUTINE interpolate(field_in,mesh_in,field_out,mesh_out)
1421 TYPE(
mesh_type),
INTENT(IN) :: mesh_in, mesh_out
1422 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: field_in
1423 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: field_out
1424 INTEGER :: m_out, m_in, n1in, n2in, nin, nout
1425 REAL(KIND=8) :: mesK, mesin, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, x, y
1426 REAL(KIND=8) :: lambda(3), rout(2), rr(2,3), a(2), b(2)
1427 INTEGER,
DIMENSION(10) :: j10
1428 f1(x,y) = -0.9d1/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y - 0.27d2/0.2d1*x*y**2 &
1429 - 0.9d1/0.2d1*y**3 + 0.9d1*x**2 + 0.18d2*x*y + 0.9d1*y**2 &
1430 - 0.11d2/0.2d1*x - 0.11d2/0.2d1*y + 0.1d1
1431 f2(x,y) = 0.9d1/0.2d1*x**3 - 0.9d1/0.2d1*x**2 + x
1432 f3(x,y) = 0.9d1/0.2d1*y**3 - 0.9d1/0.2d1*y**2 + y
1433 f4(x,y) = 0.27d2/0.2d1*x**2*y - 0.9d1/0.2d1*x*y
1434 f5(x,y) = 0.27d2/0.2d1*x*y**2 - 0.9d1/0.2d1*x*y
1435 f6(x,y) = 0.27d2/0.2d1*x**2*y + 0.27d2*x*y**2 + 0.27d2 / 0.2d1*y**3 &
1436 - 0.45d2/0.2d1*x*y - 0.45d2/0.2d1*y**2 + 0.9d1*y
1437 f7(x,y) = -0.27d2/0.2d1*x*y**2 - 0.27d2/0.2d1*y**3 + 0.9d1/0.2d1*x*y &
1438 + 0.18d2*y**2 - 0.9d1/0.2d1*y
1439 f8(x,y) = 0.27d2/0.2d1*x**3 + 0.27d2*x**2*y + 0.27d2/0.2d1*x*y**2 &
1440 - 0.45d2/0.2d1*x**2 - 0.45d2/0.2d1*x*y + 0.9d1*x
1441 f9(x,y) = -0.27d2/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y + 0.18d2*x**2 &
1442 + 0.9d1/0.2d1*x*y - 0.9d1/0.2d1*x
1443 f10(x,y) = -27*x**2*y - 27*x*y**2 + 27*x*y
1445 DO m_out = 1, mesh_out%me
1446 m_in = (m_out-1)/4 + 1
1447 rr(:,1:3) = mesh_in%rr(:,mesh_in%jj(1:3,m_in))
1452 rout = mesh_out%rr(:,mesh_out%jj(nout,m_out))
1454 n1in = modulo(nin,3) + 1
1455 n2in = modulo(nin+1,3) + 1
1459 lambda(nin) = mesin/mesk
1461 IF (abs(sum(lambda)-1.d0).GT.1.d-13)
THEN 1464 j10 = mesh_in%jj(:,m_in)
1465 field_out(mesh_out%jj(nout,m_out),:) = field_in(j10(1),:)*f1(lambda(2),lambda(3)) &
1466 + field_in(j10(2),:)*f2(lambda(2),lambda(3)) &
1467 + field_in(j10(3),:)*f3(lambda(2),lambda(3)) &
1468 + field_in(j10(4),:)*f4(lambda(2),lambda(3)) &
1469 + field_in(j10(5),:)*f5(lambda(2),lambda(3)) &
1470 + field_in(j10(6),:)*f6(lambda(2),lambda(3)) &
1471 + field_in(j10(7),:)*f7(lambda(2),lambda(3)) &
1472 + field_in(j10(8),:)*f8(lambda(2),lambda(3)) &
1473 + field_in(j10(9),:)*f9(lambda(2),lambda(3)) &
1474 + field_in(j10(10),:)*f10(lambda(2),lambda(3))
1478 FUNCTION mesurek(a,b) RESULT(mes)
1479 REAL(KIND=8),
DIMENSION(2),
INTENT(IN) :: a, b
1481 mes = abs(a(1)*b(2)-a(2)*b(1))
1486 SUBROUTINE make_vtu_file_2d(communicator, mesh_in, header, field_in, field_name, what, opt_it)
1491 TYPE(
mesh_type),
INTENT(IN),
TARGET :: mesh_in
1493 CHARACTER(*),
INTENT(IN) :: header
1494 CHARACTER(*),
INTENT(IN) :: field_name, what
1495 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
1496 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN),
TARGET :: field_in
1497 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE,
TARGET :: field_tmp
1498 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: field
1500 CHARACTER(LEN=200),
DIMENSION(:),
POINTER :: file_list
1501 CHARACTER(LEN=3) :: st_rank, st_it
1502 petscerrorcode :: ierr
1503 petscmpiint :: rank, nb_procs
1504 mpi_comm :: communicator
1505 CALL mpi_comm_rank(communicator, rank, ierr)
1506 CALL mpi_comm_size(communicator, nb_procs, ierr)
1507 ALLOCATE(file_list(nb_procs))
1509 IF (mesh_in%gauss%n_w==10 .AND.
inputs%type_fe_velocity==3)
THEN 1519 IF (
PRESENT(opt_it))
THEN 1521 WRITE(st_it,
'(I3)') it
1523 WRITE(st_rank,
'(I3)') j
1524 file_list(j) = trim(header)//
'_proc_'//trim(adjustl(st_rank))//&
1525 '_it_'//trim(adjustl(st_it))
1529 WRITE(st_rank,
'(I3)') j
1530 file_list(j) = trim(header)//
'_proc_'//trim(adjustl(st_rank))
1534 CALL check_list(communicator, file_list, mesh%np)
1536 IF (
PRESENT(opt_it))
THEN 1544 IF (
SIZE(field,2) == 6)
THEN 1546 ELSE IF (
SIZE(field,2) == 1)
THEN 1548 ELSE IF (
SIZE(field,2) .GT. 0)
THEN 1551 CALL error_petsc(
'Bug in make_vtu_file_2D: field needs at least one component')
1554 IF (
ALLOCATED(field_tmp))
DEALLOCATE(field_tmp)
type(mesh_type), target pp_local_mesh
type(petsc_csr_la), public vizu_grad_phi_la
subroutine compute_rot_h(field_in, field_out)
type(mesh_type), target phi_local_mesh
subroutine solver(my_ksp, b, x, reinit, verbose)
logical, save if_first_rot_u
type(dyn_int_line), dimension(:), pointer vv_loc_to_glob
type(dyn_int_line), dimension(:), pointer temp_loc_to_glob
logical, save if_first_rot_h
type(mesh_type), target vv_mesh_p2_iso_p4
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine, public create_my_ghost(mesh, LA, ifrom)
integer, dimension(:), pointer fourier_list_mode
type(mesh_type), pointer mesh_rot_h
subroutine, public transfer_fourier_to_real(communicator, mesh, field_in, header, name_of_field, what, list_mode, opt_it)
subroutine, public check_list(communicator, file_list, check)
subroutine, public make_vtu_file_2d(communicator, mesh_in, header, field_in, field_name, what, opt_it)
type(dyn_int_line), dimension(:), pointer phi_loc_to_glob
subroutine prepare_3d_grids(mesh, local_mesh, mesh_3d, rank, opt_nb_plane)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
subroutine error_petsc(string)
logical, save if_first_grad_phi
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
type(mesh_type), target h_mesh_3d
type(mesh_type), pointer mesh_grad_phi
subroutine interpolate(field_in, mesh_in, field_out, mesh_out)
type(mesh_type), target pp_mesh_3d
subroutine, public vtu_3d(field, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl, opt_2D, opt_mesh_in)
type(mesh_type), target h_local_mesh
subroutine, public create_xml_vtu_scal_file(field, mesh, file_name, field_name)
real(kind=8) function mesurek(a, b)
type(petsc_csr_la), public vizu_rot_u_la
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine compute_rot_u(field_in, field_out, mesh_in)
subroutine divide_mesh_into_four_subcells(mesh_in, mesh_out)
subroutine, public create_pvd_file(file_list, file_header, time_step, what)
subroutine, public sfemans_initialize_postprocessing(comm_one_d, vv_mesh_in, pp_mesh, H_mesh, phi_mesh, temp_mesh, list_mode_in, opt_nb_plane)
integer, dimension(:), pointer vizu_list_mode
subroutine divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
type(dyn_int_line), dimension(:), pointer h_loc_to_glob
type(mesh_type), target phi_mesh_3d
type(petsc_csr_la), public vizu_rot_h_la
subroutine vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
type(mesh_type), target vv_mesh_3d
type(mesh_type), target vv_local_mesh
subroutine, public make_vtu_file_3d(communicator, mesh, header, field, field_name, what, opt_it)
type(mesh_type), target temp_local_mesh
subroutine compute_grad_phi(field_in, field_out)
subroutine, public create_xml_vtu_vect_file(field, mesh, file_name, opt_st)
subroutine, public fft_par_real(communicator, V1_in, V_out, opt_nb_plane)
type(mesh_type), target temp_mesh_3d
type(dyn_int_line), dimension(:), pointer pp_loc_to_glob