25 INTEGER,
DIMENSION(:),
POINTER :: fourier_list_mode
26 INTEGER :: fourier_nb_angle
27 INTEGER :: fourier_width
28 INTEGER :: fourier_rank
29 INTEGER :: fourier_nb_procs
30 LOGICAL :: if_first_grad_phi=.TRUE.
31 INTEGER,
DIMENSION(:),
POINTER :: vizu_list_mode
32 #include "petsc/finclude/petsc.h"
33 mpi_comm :: fourier_communicator
34 mpi_comm,
DIMENSION(2) :: vizu_grad_curl_comm
36 vec :: vx_phi, vb_phi, vx_phi_ghost
41 list_mode_in, opt_nb_plane)
51 INTEGER,
DIMENSION(:),
TARGET :: list_mode_in
52 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_nb_plane
54 #include "petsc/finclude/petsc.h"
55 petscerrorcode :: ierr
56 petscmpiint :: rank, nb_procs
57 mpi_comm,
DIMENSION(2) :: comm_one_d
58 CALL mpi_comm_rank(comm_one_d(2), rank, ierr)
59 CALL mpi_comm_size(comm_one_d(2), nb_procs, ierr)
60 fourier_communicator = comm_one_d(2)
62 fourier_nb_procs = nb_procs
63 ALLOCATE(fourier_list_mode(
SIZE(list_mode_in)))
64 fourier_list_mode = list_mode_in
65 CALL
divide_mesh(comm_one_d(2), vv_mesh, vv_local_mesh, vv_loc_to_glob)
66 CALL
divide_mesh(comm_one_d(2), pp_mesh, pp_local_mesh, pp_loc_to_glob)
67 CALL
divide_mesh(comm_one_d(2), h_mesh, h_local_mesh, h_loc_to_glob)
68 CALL
divide_mesh(comm_one_d(2), phi_mesh, phi_local_mesh, phi_loc_to_glob)
69 CALL
divide_mesh(comm_one_d(2), temp_mesh, temp_local_mesh, temp_loc_to_glob)
74 CALL
prepare_3d_grids(temp_mesh, temp_local_mesh, temp_mesh_3d, opt_nb_plane)
76 mesh_grad_phi=>phi_mesh
77 CALL mpi_comm_dup(comm_one_d(1),vizu_grad_curl_comm(1) , code)
78 CALL mpi_comm_dup(comm_one_d(2),vizu_grad_curl_comm(2) , code)
79 vizu_list_mode=>list_mode_in
91 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_nb_plane
92 LOGICAL,
DIMENSION(:),
POINTER :: virgin
93 INTEGER,
DIMENSION(:),
POINTER :: renumber
94 REAL(KIND=8) :: pi, dtheta,
theta
95 INTEGER :: i, k, m, n, nb_angle, count, dim, np1, &
102 mesh_3d%gauss%n_w = 0
103 ALLOCATE(mesh_3d%jj(0,0))
104 ALLOCATE(mesh_3d%rr(0,0))
109 IF (inputs%select_mode)
THEN
110 IF (
SIZE(inputs%list_mode_lect)==1)
THEN
111 gap = inputs%list_mode_lect(1) + 1
113 gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
118 fourier_width = (inputs%m_max/fourier_nb_procs)*gap
121 m_max = gap*inputs%m_max
122 IF (present(opt_nb_plane))
THEN
123 IF (opt_nb_plane> 2*m_max-1)
THEN
124 m_max_pad = (opt_nb_plane+1)/2
131 nb_angle = 2*m_max_pad-1
132 fourier_nb_angle = nb_angle
137 dtheta = 2*pi/nb_angle
138 IF (mesh%gauss%n_w==3)
THEN
139 mesh_3d%np = local_mesh%np*nb_angle
140 mesh_3d%me = local_mesh%me*nb_angle
141 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
145 DO i = 1, local_mesh%np
147 mesh_3d%rr(1,count) = local_mesh%rr(1,i)*cos(
theta)
148 mesh_3d%rr(2,count) = local_mesh%rr(1,i)*sin(
theta)
149 mesh_3d%rr(3,count) = local_mesh%rr(2,i)
152 ALLOCATE(mesh_3d%jj(6,mesh_3d%me))
155 DO m = 1, local_mesh%me
157 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
158 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)+k*local_mesh%np
161 DO m = 1, local_mesh%me
163 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
164 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
166 ELSE IF (mesh%gauss%n_w==6)
THEN
167 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
171 DO m = 1, local_mesh%me
173 i = local_mesh%jj(n,m)
174 IF (.NOT.virgin(i)) cycle
181 mesh_3d%np = (local_mesh%np+np1)*nb_angle
182 mesh_3d%me = local_mesh%me*nb_angle
183 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
184 ALLOCATE(mesh_3d%jj(15,mesh_3d%me))
187 DO m = 1, local_mesh%me
189 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
190 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m) + k*(local_mesh%np+np1)
191 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
192 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
193 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
194 mesh_3d%jj(10,count) = local_mesh%jj(6,m) + k*(local_mesh%np+np1)
195 mesh_3d%jj(11,count) = local_mesh%jj(4,m) + k*(local_mesh%np+np1)
196 mesh_3d%jj(12,count) = local_mesh%jj(5,m) + k*(local_mesh%np+np1)
197 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
198 + (k-1)*(local_mesh%np+np1) + local_mesh%np
202 DO m = 1, local_mesh%me
204 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
205 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
206 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
207 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
208 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
209 mesh_3d%jj(10,count) = local_mesh%jj(6,m)
210 mesh_3d%jj(11,count) = local_mesh%jj(4,m)
211 mesh_3d%jj(12,count) = local_mesh%jj(5,m)
212 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
213 + (k-1)*(local_mesh%np+np1) + local_mesh%np
218 DO i = 1, local_mesh%np
219 n = (k-1)*(local_mesh%np+np1) + i
220 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(
theta)
221 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(
theta)
222 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
224 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
225 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(
theta+dtheta/2)
226 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(
theta+dtheta/2)
227 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
231 DEALLOCATE(virgin,renumber)
233 CALL
error_petsc(
'Bug in prepare_3d_grids: mesh%gauss%n_w= not correct')
237 SUBROUTINE vtu_3d(field_in, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl)
245 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
246 CHARACTER(*),
INTENT(IN) :: name_of_mesh, header, name_of_field, what
247 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
248 CHARACTER(*),
OPTIONAL,
INTENT(IN) :: opt_grad_curl
249 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_viz
251 IF (present(opt_grad_curl))
THEN
252 IF (opt_grad_curl ==
'grad')
THEN
253 ALLOCATE(field_viz(
SIZE(field_in,1), 3*
SIZE(field_in,2),
SIZE(field_in,3)))
256 ALLOCATE(field_viz(
SIZE(field_in, 1),
SIZE(field_in,2),
SIZE(field_in,3)))
260 ALLOCATE(field_viz(
SIZE(field_in, 1),
SIZE(field_in,2),
SIZE(field_in,3)))
264 IF (present(opt_it))
THEN
265 CALL
vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what, opt_it)
267 CALL
vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what)
270 DEALLOCATE(field_viz)
274 SUBROUTINE vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
282 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
283 CHARACTER(*),
INTENT(IN) :: name_of_mesh, header, name_of_field, what
284 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
287 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_fft
288 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: field_out
289 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_xyz_fft
290 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: dist_field
295 LOGICAL,
DIMENSION(:),
POINTER :: virgin
296 INTEGER,
DIMENSION(:),
POINTER :: renumber
297 INTEGER,
DIMENSION(:),
POINTER :: list_mode
299 INTEGER :: width, mode_min
300 INTEGER,
DIMENSION(1) :: loc
301 REAL(KIND=8) :: pi, dtheta,
theta
302 INTEGER :: i, k, m, n, np_b, np_e, np_loc, np_loc_max, &
303 type, nb_angle, count, dim, np1, nb_procs, rank, it
306 width = fourier_width
307 nb_procs = fourier_nb_procs
309 nb_angle = fourier_nb_angle
310 list_mode => fourier_list_mode
312 IF (name_of_mesh==
'vv_mesh')
THEN
316 mesh_3d => vv_mesh_3d
317 local_mesh => vv_local_mesh
318 ELSE IF (name_of_mesh==
'pp_mesh')
THEN
322 mesh_3d => pp_mesh_3d
323 local_mesh => pp_local_mesh
324 ELSE IF (name_of_mesh==
'H_mesh')
THEN
329 local_mesh => h_local_mesh
330 ELSE IF (name_of_mesh==
'phi_mesh')
THEN
334 mesh_3d => phi_mesh_3d
335 local_mesh => phi_local_mesh
336 ELSE IF (name_of_mesh==
'temp_mesh')
THEN
340 mesh_3d => temp_mesh_3d
341 local_mesh => temp_local_mesh
343 CALL
error_petsc(
'Bug in vtu_3d: name_of_mesh is not correct')
349 np_loc_max = max(np_loc_max,maxval(
loc_to_glob(n)%DIL))
351 ALLOCATE(dist_field(nb_procs*np_loc_max,
SIZE(field_in,2),width))
354 mode_min = rank*width
356 IF (minval(abs(list_mode - (mode_min+i-1)))/=0) cycle
357 loc = minloc(abs(list_mode - (mode_min+i-1)))
358 DO type = 1, size(field_in,2)
360 np_b = (n-1)*np_loc_max + 1
362 np_e = np_b + np_loc - 1
363 dist_field(np_b:np_e,type,i) = field_in(
loc_to_glob(n)%DIL,type,loc(1))
369 CALL
fft_par_real(fourier_communicator, dist_field, field_out_fft, opt_nb_plane=nb_angle)
372 ALLOCATE(field_out_xyz_fft(
SIZE(field_out_fft,1),
SIZE(field_out_fft,2),
SIZE(field_out_fft,3)))
373 IF (
SIZE(field_in,2)==2)
THEN
374 field_out_xyz_fft=field_out_fft
377 dtheta = 2*pi/nb_angle
380 field_out_xyz_fft(k,1,:) = field_out_fft(k,1,:)*cos(
theta)-field_out_fft(k,2,:)*sin(
theta)
381 field_out_xyz_fft(k,2,:) = field_out_fft(k,1,:)*sin(
theta)+field_out_fft(k,2,:)*cos(
theta)
382 field_out_xyz_fft(k,3,:) = field_out_fft(k,3,:)
387 IF (
SIZE(field_in,2)==2)
THEN
389 ELSE IF (
SIZE(field_in,2)==6)
THEN
392 CALL
error_petsc(
'Bug in vtu_3d: SIZE(field_in,2) not correct')
394 IF (nb_angle /=
SIZE(field_out_xyz_fft,1))
THEN
395 CALL
error_petsc(
'Bug in vtu_3d: nb_angle /= SIZE(field_out_xyz_FFT,1')
400 dtheta = 2*pi/nb_angle
401 IF (local_mesh%gauss%n_w==3)
THEN
402 ALLOCATE(field_out(dim,mesh_3d%np))
405 DO i = 1, local_mesh%np
407 field_out(:,count) = field_out_xyz_fft(k,:,i)
411 ELSE IF (local_mesh%gauss%n_w==6)
THEN
412 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
416 DO m = 1, local_mesh%me
418 i = local_mesh%jj(n,m)
419 IF (.NOT.virgin(i)) cycle
427 ALLOCATE(field_out(dim,mesh_3d%np))
429 DO i = 1, local_mesh%np
430 n = (k-1)*(local_mesh%np+np1) + i
431 field_out(:,n) = field_out_xyz_fft(k,:,i)
435 DO i = 1, local_mesh%np
437 DO k = 1, nb_angle - 2
438 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
439 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
440 + (3.d0/4)*field_out_xyz_fft(k+1,:,i) - (1.d0/8)*field_out_xyz_fft(k+2,:,i)
444 DO i = 1, local_mesh%np
446 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
447 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
448 + (3.d0/4)*field_out_xyz_fft(k+1,:,i) - (1.d0/8)*field_out_xyz_fft(1,:,i)
451 DO i = 1, local_mesh%np
453 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
454 field_out(:,n) = (3.d0/8)*field_out_xyz_fft(k,:,i) &
455 + (3.d0/4)*field_out_xyz_fft(1,:,i) - (1.d0/8)*field_out_xyz_fft(2,:,i)
457 DEALLOCATE(virgin, renumber)
459 CALL
error_petsc(
'Bug in vtu_3d: mesh%gauss%n_w is not correct')
462 IF (present(opt_it))
THEN
465 field_out, name_of_field, what, opt_it = it)
468 field_out, name_of_field, what)
472 DEALLOCATE(dist_field, field_out,
loc_to_glob,field_out_xyz_fft)
475 DEALLOCATE(field_out_fft)
479 header, name_of_field, what, list_mode, opt_it)
488 CHARACTER(*),
INTENT(IN) :: header, name_of_field, what
489 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
490 INTEGER,
OPTIONAL,
INTENT(IN) :: opt_it
491 REAL(KIND=8),
DIMENSION(:,:,:) :: field_in
493 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: field_out_fft
494 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: field_out
495 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE :: dist_field
500 LOGICAL,
DIMENSION(:),
POINTER :: virgin
501 INTEGER,
DIMENSION(:),
POINTER :: renumber
504 INTEGER :: gap, width, mode_min
505 INTEGER,
DIMENSION(1) :: loc
506 REAL(KIND=8) :: pi, dtheta,
theta
507 INTEGER :: i, k, m, n, np_b, np_e, np_loc, np_loc_max, &
508 type, nb_angle, count, dim, np1
509 #include "petsc/finclude/petsc.h"
510 petscerrorcode :: ierr
511 petscmpiint :: rank, nb_procs
512 mpi_comm :: communicator
513 CALL mpi_comm_rank(communicator, rank, ierr)
514 CALL mpi_comm_size(communicator, nb_procs, ierr)
526 IF (inputs%select_mode)
THEN
527 IF (
SIZE(inputs%list_mode_lect)==1)
THEN
528 gap = inputs%list_mode_lect(1) + 1
530 gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
535 width =
SIZE(list_mode)*gap
540 np_loc_max = max(np_loc_max,maxval(
loc_to_glob(n)%DIL))
542 ALLOCATE(dist_field(nb_procs*np_loc_max,
SIZE(field_in,2),width))
547 mode_min = rank*width
549 IF (minval(abs(list_mode - (mode_min+i-1)))/=0) cycle
550 loc = minloc(abs(list_mode - (mode_min+i-1)))
551 DO type = 1, size(field_in,2)
553 np_b = (n-1)*np_loc_max + 1
555 np_e = np_b + np_loc - 1
556 dist_field(np_b:np_e,type,i) = field_in(
loc_to_glob(n)%DIL,type,loc(1))
561 CALL
fft_par_real(communicator, dist_field, field_out_fft, opt_nb_plane=50)
575 IF (
SIZE(field_in,2)==2)
THEN
577 ELSE IF (
SIZE(field_in,2)==6)
THEN
580 CALL
error_petsc(
'Bug in transfer_fourier_to_real: SIZE(field_in,2) not correct')
582 nb_angle =
SIZE(field_out_fft,1)
584 dtheta = 2*pi/nb_angle
585 IF (mesh%gauss%n_w==3)
THEN
587 mesh_3d%np = local_mesh%np*nb_angle
588 mesh_3d%me = local_mesh%me*nb_angle
589 ALLOCATE(field_out(dim,mesh_3d%np))
592 DO i = 1, local_mesh%np
594 field_out(:,count) = field_out_fft(k,:,i)
598 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
602 DO i = 1, local_mesh%np
604 mesh_3d%rr(1,count) = local_mesh%rr(1,i)*cos(
theta)
605 mesh_3d%rr(2,count) = local_mesh%rr(1,i)*sin(
theta)
606 mesh_3d%rr(3,count) = local_mesh%rr(2,i)
609 ALLOCATE(mesh_3d%jj(6,mesh_3d%me))
612 DO m = 1, local_mesh%me
614 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
615 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)+k*local_mesh%np
618 DO m = 1, local_mesh%me
620 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m)+(k-1)*local_mesh%np
621 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
623 ELSE IF (mesh%gauss%n_w==6)
THEN
624 ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
628 DO m = 1, local_mesh%me
630 i = local_mesh%jj(n,m)
631 IF (.NOT.virgin(i)) cycle
638 mesh_3d%np = (local_mesh%np+np1)*nb_angle
639 mesh_3d%me = local_mesh%me*nb_angle
640 ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
641 ALLOCATE(mesh_3d%jj(15,mesh_3d%me))
644 DO m = 1, local_mesh%me
646 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
647 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m) + k*(local_mesh%np+np1)
648 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
649 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
650 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
651 mesh_3d%jj(10,count) = local_mesh%jj(6,m) + k*(local_mesh%np+np1)
652 mesh_3d%jj(11,count) = local_mesh%jj(4,m) + k*(local_mesh%np+np1)
653 mesh_3d%jj(12,count) = local_mesh%jj(5,m) + k*(local_mesh%np+np1)
654 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
655 + (k-1)*(local_mesh%np+np1) + local_mesh%np
659 DO m = 1, local_mesh%me
661 mesh_3d%jj(1:3,count) = local_mesh%jj(1:3,m) + (k-1)*(local_mesh%np+np1)
662 mesh_3d%jj(4:6,count) = local_mesh%jj(1:3,m)
663 mesh_3d%jj(7,count) = local_mesh%jj(6,m) + (k-1)*(local_mesh%np+np1)
664 mesh_3d%jj(8,count) = local_mesh%jj(4,m) + (k-1)*(local_mesh%np+np1)
665 mesh_3d%jj(9,count) = local_mesh%jj(5,m) + (k-1)*(local_mesh%np+np1)
666 mesh_3d%jj(10,count) = local_mesh%jj(6,m)
667 mesh_3d%jj(11,count) = local_mesh%jj(4,m)
668 mesh_3d%jj(12,count) = local_mesh%jj(5,m)
669 mesh_3d%jj(13:15,count) = renumber(local_mesh%jj(1:3,m)) &
670 + (k-1)*(local_mesh%np+np1) + local_mesh%np
675 DO i = 1, local_mesh%np
676 n = (k-1)*(local_mesh%np+np1) + i
677 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(
theta)
678 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(
theta)
679 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
681 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
682 mesh_3d%rr(1,n) = local_mesh%rr(1,i)*cos(
theta+dtheta/2)
683 mesh_3d%rr(2,n) = local_mesh%rr(1,i)*sin(
theta+dtheta/2)
684 mesh_3d%rr(3,n) = local_mesh%rr(2,i)
688 ALLOCATE(field_out(dim,mesh_3d%np))
693 DO i = 1, local_mesh%np
694 n = (k-1)*(local_mesh%np+np1) + i
696 field_out(:,n) = field_out_fft(k,:,i)
708 DO i = 1, local_mesh%np
710 DO k = 1, nb_angle - 2
711 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
712 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
713 + (3.d0/4)*field_out_fft(k+1,:,i) - (1.d0/8)*field_out_fft(k+2,:,i)
717 DO i = 1, local_mesh%np
719 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
720 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
721 + (3.d0/4)*field_out_fft(k+1,:,i) - (1.d0/8)*field_out_fft(1,:,i)
724 DO i = 1, local_mesh%np
726 n = (k-1)*(local_mesh%np+np1) + local_mesh%np + renumber(i)
727 field_out(:,n) = (3.d0/8)*field_out_fft(k,:,i) &
728 + (3.d0/4)*field_out_fft(1,:,i) - (1.d0/8)*field_out_fft(2,:,i)
733 DEALLOCATE(virgin, renumber)
736 CALL
error_petsc(
'Bug in transfer_fourier_to_real: mesh%gauss%n_w is not correct')
739 IF (present(opt_it))
THEN
741 field_out, name_of_field, what, opt_it=opt_it)
744 field_out, name_of_field, what)
752 DEALLOCATE(mesh_3d%rr, mesh_3d%jj)
753 DEALLOCATE(local_mesh%rr, local_mesh%jj)
754 DEALLOCATE(dist_field, field_out)
758 DEALLOCATE(field_out_fft)
761 SUBROUTINE divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
766 TYPE(mesh_type),
INTENT(OUT) :: local_mesh
767 INTEGER :: me_b, me_e, me_loc, np_loc
769 INTEGER,
POINTER,
DIMENSION(:) :: glob_to_loc
770 LOGICAL,
DIMENSION(mesh%np) :: virgin
771 INTEGER :: count, m, m_glob, n, i, q, r, loop_rank
772 #include "petsc/finclude/petsc.h"
773 petscerrorcode :: ierr
774 petscmpiint :: rank, nb_procs
775 mpi_comm :: communicator
781 local_mesh%gauss%n_w = 0
782 ALLOCATE(local_mesh%jj(0,0))
783 ALLOCATE(local_mesh%rr(0,0))
789 CALL mpi_comm_rank(communicator,rank,ierr)
790 CALL mpi_comm_size(communicator,nb_procs,ierr)
792 ALLOCATE(glob_to_loc(mesh%np))
796 r = mesh%me - q*nb_procs
798 DO loop_rank = 0, nb_procs-1
800 IF (loop_rank+1.LE.r)
THEN
801 me_b = loop_rank*(q+1) + 1
804 me_b = loop_rank*q + r + 1
807 me_loc = me_e - me_b + 1
812 DO n = 1, mesh%gauss%n_w
814 IF (.NOT.virgin(i)) cycle
824 DO n = 1, mesh%gauss%n_w
826 IF (.NOT.virgin(i)) cycle
830 IF (loop_rank == rank) glob_to_loc(i) = count
834 IF (loop_rank == rank)
THEN
836 local_mesh%me = me_loc
837 local_mesh%np = np_loc
838 local_mesh%gauss%n_w = mesh%gauss%n_w
839 ALLOCATE(local_mesh%jj(mesh%gauss%n_w,me_loc))
842 local_mesh%jj(:,m) = glob_to_loc(mesh%jj(:,m_glob))
844 ALLOCATE(local_mesh%rr(mesh%gauss%k_d,local_mesh%np))
849 DEALLOCATE(glob_to_loc)
859 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field_in
860 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: field_out
862 INTEGER,
POINTER,
DIMENSION(:) :: ifrom
864 REAL(KIND=8),
DIMENSION(SIZE(field_out,1), SIZE(field_out,2)) :: smb_grad_phi
865 REAL(KIND=8),
DIMENSION(6) :: grad_phi_loc
866 INTEGER,
DIMENSION(mesh_grad_phi%gauss%n_w) :: i_loc
867 INTEGER :: i, mode, m, l, j_loc, k, nj
868 REAL(KIND=8) :: ray, drdthetadz
869 #include "petsc/finclude/petsc.h"
870 petscerrorcode :: ierr
873 param_grad_phi%it_max=1
874 param_grad_phi%rel_tol=1.d-6
875 param_grad_phi%abs_tol=1.d-10
876 param_grad_phi%solver=
'GMRES'
877 param_grad_phi%precond=
'MUMPS'
878 param_grad_phi%verbose=.false.
880 IF (if_first_grad_phi)
THEN
881 if_first_grad_phi=.false.
884 CALL veccreateghost(vizu_grad_curl_comm(1), mesh_grad_phi%dom_np, &
885 petsc_determine,
SIZE(ifrom), ifrom, vx_phi, ierr)
886 CALL vecghostgetlocalform(vx_phi, vx_phi_ghost, ierr)
887 CALL vecduplicate(vx_phi, vb_phi, ierr)
890 CALL
qs_diff_mass_scal_m(mesh_grad_phi, vizu_grad_phi_la, 0.d0, 1.d0, 0.d0, 0, mat_grad_phi)
891 CALL
init_solver(param_grad_phi,ksp_grad_phi,mat_grad_phi,vizu_grad_curl_comm(1),&
892 solver=
'GMRES',precond=
'MUMPS')
893 CALL matdestroy(mat_grad_phi,ierr)
899 DO i = 1,
SIZE(vizu_list_mode)
900 mode = vizu_list_mode(i)
902 DO m = 1, mesh_grad_phi%me
903 DO l = 1, mesh_grad_phi%gauss%l_G
904 ray = sum(mesh_grad_phi%rr(1,mesh_grad_phi%jj(:,m))*mesh_grad_phi%gauss%ww(:,l))
905 drdthetadz = mesh_grad_phi%gauss%rj(l,m)
907 DO nj = 1,mesh_grad_phi%gauss%n_w
908 j_loc = mesh_grad_phi%jj(nj,m)
909 grad_phi_loc(1) = grad_phi_loc(1) + field_in(j_loc,1,i)*mesh_grad_phi%gauss%dw(1,nj,l,m)*ray
910 grad_phi_loc(2) = grad_phi_loc(2) + field_in(j_loc,2,i)*mesh_grad_phi%gauss%dw(1,nj,l,m)*ray
911 grad_phi_loc(3) = grad_phi_loc(3) + field_in(j_loc,2,i)*mesh_grad_phi%gauss%ww(nj,l)*mode
912 grad_phi_loc(4) = grad_phi_loc(4) - field_in(j_loc,1,i)*mesh_grad_phi%gauss%ww(nj,l)*mode
913 grad_phi_loc(5) = grad_phi_loc(5) + field_in(j_loc,1,i)*mesh_grad_phi%gauss%dw(2,nj,l,m)*ray
914 grad_phi_loc(6) = grad_phi_loc(6) + field_in(j_loc,2,i)*mesh_grad_phi%gauss%dw(2,nj,l,m)*ray
916 i_loc = mesh_grad_phi%jj(:,m)
919 smb_grad_phi(i_loc,k) = smb_grad_phi(i_loc,k) + grad_phi_loc(k)* &
920 mesh_grad_phi%gauss%ww(:,l)*drdthetadz
926 CALL veczeroentries(vb_phi, ierr)
927 CALL vecsetvalues(vb_phi, mesh_grad_phi%np, vizu_grad_phi_la%loc_to_glob(1,:)-1, smb_grad_phi(:,k), add_values, ierr)
928 CALL vecassemblybegin(vb_phi,ierr)
929 CALL vecassemblyend(vb_phi,ierr)
931 CALL
solver(ksp_grad_phi,vb_phi,vx_phi,reinit=.false.,
verbose=.false.)
933 CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
934 CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
935 IF (mesh_grad_phi%me/=0)
THEN
936 CALL
extract(vx_phi_ghost,1,1,vizu_grad_phi_la,field_out(:,k,i))
941 field_out(:,2,i) = 0.d0
942 field_out(:,4,i) = 0.d0
943 field_out(:,6,i) = 0.d0
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine, public vtu_3d(field_in, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl)
subroutine, public fft_par_real(communicator, V1_in, V_out, opt_nb_plane)
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic theta
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
subroutine compute_grad_phi(field_in, field_out)
subroutine, public transfer_fourier_to_real(communicator, mesh, field_in, header, name_of_field, what, list_mode, opt_it)
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine prepare_3d_grids(mesh, local_mesh, mesh_3d, opt_nb_plane)
subroutine vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
subroutine, public make_vtu_file_3d(communicator, mesh, header, field, field_name, what, opt_it)
subroutine, public sfemans_initialize_postprocessing(comm_one_d, vv_mesh, pp_mesh, H_mesh, phi_mesh, temp_mesh, list_mode_in, opt_nb_plane)
subroutine error_petsc(string)
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)