SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
fourier_real.f90
Go to the documentation of this file.
2  USE dyn_line
3  USE def_type_mesh
6  PUBLIC :: vtu_3d
7  TYPE(petsc_csr_la), PUBLIC :: vizu_grad_phi_LA
8  TYPE(mesh_type), POINTER :: mesh_grad_phi
9  PRIVATE
10  TYPE(mesh_type), TARGET :: vv_mesh_3d
11  TYPE(mesh_type), TARGET :: pp_mesh_3d
12  TYPE(mesh_type), TARGET :: H_mesh_3d
13  TYPE(mesh_type), TARGET :: phi_mesh_3d
14  TYPE(mesh_type), TARGET :: temp_mesh_3d
15  TYPE(mesh_type), TARGET :: vv_local_mesh
16  TYPE(mesh_type), TARGET :: pp_local_mesh
17  TYPE(mesh_type), TARGET :: H_local_mesh
18  TYPE(mesh_type), TARGET :: phi_local_mesh
19  TYPE(mesh_type), TARGET :: temp_local_mesh
20  TYPE(dyn_int_line), DIMENSION(:), POINTER :: vv_loc_to_glob
21  TYPE(dyn_int_line), DIMENSION(:), POINTER :: pp_loc_to_glob
22  TYPE(dyn_int_line), DIMENSION(:), POINTER :: H_loc_to_glob
23  TYPE(dyn_int_line), DIMENSION(:), POINTER :: phi_loc_to_glob
24  TYPE(dyn_int_line), DIMENSION(:), POINTER :: temp_loc_to_glob
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
35  mat :: mat_grad_phi
36  vec :: vx_phi, vb_phi, vx_phi_ghost
37  ksp :: ksp_grad_phi
38 CONTAINS
39 
40  SUBROUTINE sfemans_initialize_postprocessing(comm_one_d, vv_mesh, pp_mesh, H_mesh, phi_mesh, temp_mesh, &
41  list_mode_in, opt_nb_plane)
42  USE my_util
43  IMPLICIT NONE
44  TYPE(mesh_type), INTENT(IN) :: vv_mesh
45  TYPE(mesh_type), INTENT(IN) :: pp_mesh
46  TYPE(mesh_type), INTENT(IN) :: h_mesh
47 !!$ TYPE(mesh_type), INTENT(IN) :: phi_mesh
48  TYPE(mesh_type), TARGET :: phi_mesh
49  TYPE(mesh_type), INTENT(IN) :: temp_mesh
50 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode_in
51  INTEGER, DIMENSION(:), TARGET :: list_mode_in
52  INTEGER, OPTIONAL, INTENT(IN) :: opt_nb_plane
53  INTEGER :: code
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)
61  fourier_rank = rank
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)
70  CALL prepare_3d_grids(vv_mesh, vv_local_mesh, vv_mesh_3d, opt_nb_plane)
71  CALL prepare_3d_grids(pp_mesh, pp_local_mesh, pp_mesh_3d, opt_nb_plane)
72  CALL prepare_3d_grids(h_mesh, h_local_mesh, h_mesh_3d, opt_nb_plane)
73  CALL prepare_3d_grids(phi_mesh, phi_local_mesh, phi_mesh_3d, opt_nb_plane)
74  CALL prepare_3d_grids(temp_mesh, temp_local_mesh, temp_mesh_3d, opt_nb_plane)
75 
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
81 
82  SUBROUTINE prepare_3d_grids(mesh, local_mesh, mesh_3d, opt_nb_plane)
83  USE dyn_line
84  USE def_type_mesh
85  USE input_data
86  USE my_util
87  IMPLICIT NONE
88  TYPE(mesh_type), INTENT(IN) :: mesh
89  TYPE(mesh_type), INTENT(IN) :: local_mesh
90  TYPE(mesh_type), INTENT(OUT):: mesh_3d
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, &
96  gap, m_max, m_max_pad
97 
98  !===Deal with empty mesh
99  IF (mesh%me==0) THEN
100  mesh_3d%me = 0
101  mesh_3d%np = 0
102  mesh_3d%gauss%n_w = 0
103  ALLOCATE(mesh_3d%jj(0,0))
104  ALLOCATE(mesh_3d%rr(0,0))
105  RETURN
106  END IF
107 
108  !===Compute gap
109  IF (inputs%select_mode) THEN
110  IF (SIZE(inputs%list_mode_lect)==1) THEN
111  gap = inputs%list_mode_lect(1) + 1
112  ELSE
113  gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
114  END IF
115  ELSE
116  gap = 1
117  END IF
118  fourier_width = (inputs%m_max/fourier_nb_procs)*gap
119 
120  !===Compute nb_angle
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
125  ELSE
126  m_max_pad = m_max
127  END IF
128  ELSE
129  m_max_pad = m_max
130  END IF
131  nb_angle = 2*m_max_pad-1
132  fourier_nb_angle = nb_angle
133 
134  !===Create 3d mesh
135  dim = 3
136  pi = acos(-1.d0)
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))
142  count = 0
143  DO k = 1, nb_angle
144  theta = (k-1)*dtheta
145  DO i = 1, local_mesh%np
146  count = count + 1
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)
150  ENDDO
151  END DO
152  ALLOCATE(mesh_3d%jj(6,mesh_3d%me)) !P1 3d VTK wedge 13
153  count = 0
154  DO k = 1, nb_angle-1
155  DO m = 1, local_mesh%me
156  count = count + 1
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
159  END DO
160  END DO
161  DO m = 1, local_mesh%me
162  count = count + 1
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)
165  END DO
166  ELSE IF (mesh%gauss%n_w==6) THEN
167  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
168  virgin = .true.
169  renumber = -1
170  count = 0
171  DO m = 1, local_mesh%me
172  DO n = 1, 3
173  i = local_mesh%jj(n,m)
174  IF (.NOT.virgin(i)) cycle
175  virgin(i) = .false.
176  count = count + 1 !===Count P1 nodes on local_mesh
177  renumber(i) = count
178  END DO
179  END DO
180  np1 = count
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)) !P2 3d VTK wedge 26
185  count = 0
186  DO k = 1, nb_angle-1
187  DO m = 1, local_mesh%me
188  count = count + 1
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
199  END DO
200  END DO
201  k = nb_angle
202  DO m = 1, local_mesh%me
203  count = count + 1
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
214  END DO
215  !===Create 3d mesh
216  DO k = 1, nb_angle
217  theta = (k-1)*dtheta
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)
223  IF (virgin(i)) cycle !===This is a vertex
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)
228  ENDDO
229  END DO
230  !===Cleanup
231  DEALLOCATE(virgin,renumber)
232  ELSE
233  CALL error_petsc('Bug in prepare_3d_grids: mesh%gauss%n_w= not correct')
234  END IF
235  END SUBROUTINE prepare_3d_grids
236 
237  SUBROUTINE vtu_3d(field_in, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl)
238  USE dyn_line
239  USE def_type_mesh
240  USE vtk_viz
241  USE sft_parallele
242  USE input_data
243  USE my_util
244  IMPLICIT NONE
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
250 
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)))
254  CALL compute_grad_phi(field_in, field_viz)
255  ELSE
256  ALLOCATE(field_viz(SIZE(field_in, 1), SIZE(field_in,2), SIZE(field_in,3)))
257  field_viz = field_in
258  END IF
259  ELSE
260  ALLOCATE(field_viz(SIZE(field_in, 1), SIZE(field_in,2), SIZE(field_in,3)))
261  field_viz = field_in
262  END IF
263 
264  IF (present(opt_it)) THEN
265  CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what, opt_it)
266  ELSE
267  CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what)
268  END IF
269 
270  DEALLOCATE(field_viz)
271 
272  END SUBROUTINE vtu_3d
273 
274  SUBROUTINE vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
275  USE dyn_line
276  USE def_type_mesh
277  USE vtk_viz
278  USE sft_parallele
279  USE input_data
280  USE my_util
281  IMPLICIT NONE
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
285  TYPE(mesh_type), POINTER :: mesh_3d
286  TYPE(mesh_type), POINTER :: local_mesh
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
291 ! FL-CN possible faute 19/03/2013
292  !REAL(KIND=8), DIMENSION(:,:,:), POINTER :: field_out_FFT
293  !REAL(KIND=8), DIMENSION(:,:), POINTER :: field_out
294  !REAL(KIND=8), DIMENSION(:,:,:), POINTER :: dist_field
295  LOGICAL, DIMENSION(:), POINTER :: virgin
296  INTEGER, DIMENSION(:), POINTER :: renumber
297  INTEGER, DIMENSION(:), POINTER :: list_mode
298  TYPE(dyn_int_line), DIMENSION(:), POINTER :: loc_to_glob
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
304 
305  !===Initialization
306  width = fourier_width
307  nb_procs = fourier_nb_procs
308  rank = fourier_rank
309  nb_angle = fourier_nb_angle
310  list_mode => fourier_list_mode
311  ALLOCATE(loc_to_glob(nb_procs))
312  IF (name_of_mesh=='vv_mesh') THEN
313  DO n = 1, nb_procs
314  loc_to_glob(n)%DIL => vv_loc_to_glob(n)%DIL
315  END DO
316  mesh_3d => vv_mesh_3d
317  local_mesh => vv_local_mesh
318  ELSE IF (name_of_mesh=='pp_mesh') THEN
319  DO n = 1, nb_procs
320  loc_to_glob(n)%DIL => pp_loc_to_glob(n)%DIL
321  END DO
322  mesh_3d => pp_mesh_3d
323  local_mesh => pp_local_mesh
324  ELSE IF (name_of_mesh=='H_mesh') THEN
325  DO n = 1, nb_procs
326  loc_to_glob(n)%DIL => h_loc_to_glob(n)%DIL
327  END DO
328  mesh_3d => h_mesh_3d
329  local_mesh => h_local_mesh
330  ELSE IF (name_of_mesh=='phi_mesh') THEN
331  DO n = 1, nb_procs
332  loc_to_glob(n)%DIL => phi_loc_to_glob(n)%DIL
333  END DO
334  mesh_3d => phi_mesh_3d
335  local_mesh => phi_local_mesh
336  ELSE IF (name_of_mesh=='temp_mesh') THEN
337  DO n = 1, nb_procs
338  loc_to_glob(n)%DIL => temp_loc_to_glob(n)%DIL
339  END DO
340  mesh_3d => temp_mesh_3d
341  local_mesh => temp_local_mesh
342  ELSE
343  CALL error_petsc('Bug in vtu_3d: name_of_mesh is not correct')
344  END IF
345 
346  !===Compute np_loc_max
347  np_loc_max = -1
348  DO n = 1, nb_procs
349  np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
350  END DO
351  ALLOCATE(dist_field(nb_procs*np_loc_max,SIZE(field_in,2),width))
352 
353  dist_field = 0.d0
354  mode_min = rank*width
355  DO i = 1, 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)
359  DO n = 1, nb_procs
360  np_b = (n-1)*np_loc_max + 1
361  np_loc = SIZE(loc_to_glob(n)%DIL)
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))
364  END DO
365  END DO
366  END DO
367 
368  !===Compute field_out_FFT in real space on local_mesh
369  CALL fft_par_real(fourier_communicator, dist_field, field_out_fft, opt_nb_plane=nb_angle)
370 
371 !++++++++++FL 29/03/2013
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
375  ELSE
376  pi = acos(-1.d0)
377  dtheta = 2*pi/nb_angle
378  DO k = 1, nb_angle
379  theta = (k-1)*dtheta
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,:)
383  END DO
384  END IF
385 !++++++++++FL 29/03/2013
386 
387  IF (SIZE(field_in,2)==2) THEN
388  dim = 1
389  ELSE IF (SIZE(field_in,2)==6) THEN
390  dim = 3
391  ELSE
392  CALL error_petsc('Bug in vtu_3d: SIZE(field_in,2) not correct')
393  END IF
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')
396  END IF
397 
398  !===Reconstruct field_out in real space on mesh_3d
399  pi = acos(-1.d0)
400  dtheta = 2*pi/nb_angle
401  IF (local_mesh%gauss%n_w==3) THEN
402  ALLOCATE(field_out(dim,mesh_3d%np))
403  count = 0
404  DO k = 1, nb_angle
405  DO i = 1, local_mesh%np
406  count = count + 1
407  field_out(:,count) = field_out_xyz_fft(k,:,i)
408  ENDDO
409  END DO
410 
411  ELSE IF (local_mesh%gauss%n_w==6) THEN
412  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
413  virgin = .true.
414  renumber = -1
415  count = 0
416  DO m = 1, local_mesh%me
417  DO n = 1, 3
418  i = local_mesh%jj(n,m)
419  IF (.NOT.virgin(i)) cycle
420  virgin(i) = .false.
421  count = count + 1 !===Count P1 nodes on local_mesh
422  renumber(i) = count
423  END DO
424  END DO
425  np1 = count
426 
427  ALLOCATE(field_out(dim,mesh_3d%np))
428  DO k = 1, nb_angle
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)
432  END DO
433  END DO
434 
435  DO i = 1, local_mesh%np
436  IF (virgin(i)) cycle !===This is a vertex
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)
441  END DO
442  END DO
443  k = nb_angle - 1
444  DO i = 1, local_mesh%np
445  IF (virgin(i)) cycle !===This is a vertex
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)
449  END DO
450  k = nb_angle
451  DO i = 1, local_mesh%np
452  IF (virgin(i)) cycle !===This is a vertex
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)
456  END DO
457  DEALLOCATE(virgin, renumber)
458  ELSE
459  CALL error_petsc('Bug in vtu_3d: mesh%gauss%n_w is not correct')
460  END IF
461 
462  IF (present(opt_it)) THEN
463  it = opt_it
464  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
465  field_out, name_of_field, what, opt_it = it)
466  ELSE
467  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
468  field_out, name_of_field, what)
469  END IF
470 
471  !===Cleanup
472  DEALLOCATE(dist_field, field_out, loc_to_glob,field_out_xyz_fft)
473 ! FL-CN possible faute ?? 19/03/2013
474  !NULLIFY(field_out_FFT)
475  DEALLOCATE(field_out_fft)
476  END SUBROUTINE vtu_3d_base
477 
478  SUBROUTINE transfer_fourier_to_real(communicator, mesh, field_in, &
479  header, name_of_field, what, list_mode, opt_it)
480  USE dyn_line
481  USE def_type_mesh
482  USE vtk_viz
483  USE sft_parallele
484  USE input_data
485  USE my_util
486  IMPLICIT NONE
487  TYPE(mesh_type), INTENT(IN) :: mesh
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
492  TYPE(mesh_type) :: mesh_3d
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
496 ! FL-CN possible faute 19/03/2013
497  !REAL(KIND=8), DIMENSION(:,:,:), POINTER :: field_out_FFT
498  !REAL(KIND=8), DIMENSION(:,:), POINTER :: field_out
499  !REAL(KIND=8), DIMENSION(:,:,:), POINTER :: dist_field
500  LOGICAL, DIMENSION(:), POINTER :: virgin
501  INTEGER, DIMENSION(:), POINTER :: renumber
502  TYPE(mesh_type) :: local_mesh
503  TYPE(dyn_int_line), DIMENSION(:), POINTER :: loc_to_glob
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)
515 
516  !field_in =0.d0
517  !IF (rank==2) THEN
518  ! field_in(:,1,1) = mesh%rr(1,:)
519  ! field_in(:,6,1) = mesh%rr(2,:)
520  !END IF
521  !IF (rank==1) THEN
522  ! field_in(:,3,1) = mesh%rr(1,:)
523  ! field_in(:,4,1) = mesh%rr(2,:)
524  !END IF
525 
526  IF (inputs%select_mode) THEN
527  IF (SIZE(inputs%list_mode_lect)==1) THEN
528  gap = inputs%list_mode_lect(1) + 1
529  ELSE
530  gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
531  END IF
532  ELSE
533  gap = 1
534  END IF
535  width = SIZE(list_mode)*gap
536 
537  CALL divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
538  np_loc_max = -1
539  DO n = 1, nb_procs
540  np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
541  END DO
542  ALLOCATE(dist_field(nb_procs*np_loc_max,SIZE(field_in,2),width))
543  !TEST
544  !write(*,*) ' TEST', local_mesh%np, MAXVAL(local_mesh%jj)
545  !TEST
546  dist_field = 0.d0
547  mode_min = rank*width
548  DO i = 1, 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)
552  DO n = 1, nb_procs
553  np_b = (n-1)*np_loc_max + 1
554  np_loc = SIZE(loc_to_glob(n)%DIL)
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))
557  END DO
558  END DO
559  END DO
560 
561  CALL fft_par_real(communicator, dist_field, field_out_fft, opt_nb_plane=50)
562 
563  !i = SIZE(field_out_FFT,1)
564  !write(*,*) ' i ', i
565  !write(*,*) 'proc', rank, MAXVAL(ABS(field_out_FFT(2,2,1:local_mesh%np) &
566  ! - local_mesh%rr(1,:)*COS(4*(2*ACOS(-1.d0))/i) &
567  ! - local_mesh%rr(2,:)*SIN(4*(2*ACOS(-1.d0))/i))) &
568  ! + MAXVAL(ABS(field_out_FFT(2,1,1:local_mesh%np) &
569  ! - local_mesh%rr(1,:)*COS(8*(2*ACOS(-1.d0))/i))) &
570  ! + MAXVAL(ABS(field_out_FFT(2,3,1:local_mesh%np) &
571  ! - local_mesh%rr(2,:)*SIN(8*(2*ACOS(-1.d0))/i)))
572 
573  !write(*,*) 'rank', rank, index, field_out(:,1,local_mesh%np/2)
574 
575  IF (SIZE(field_in,2)==2) THEN
576  dim = 1
577  ELSE IF (SIZE(field_in,2)==6) THEN
578  dim = 3
579  ELSE
580  CALL error_petsc('Bug in transfer_fourier_to_real: SIZE(field_in,2) not correct')
581  END IF
582  nb_angle = SIZE(field_out_fft,1)
583  pi = acos(-1.d0)
584  dtheta = 2*pi/nb_angle
585  IF (mesh%gauss%n_w==3) THEN
586  !IF (.TRUE.) 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))
590  count = 0
591  DO k = 1, nb_angle
592  DO i = 1, local_mesh%np
593  count = count + 1
594  field_out(:,count) = field_out_fft(k,:,i)
595  ENDDO
596  END DO
597  !===Create 3d mesh
598  ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
599  count = 0
600  DO k = 1, nb_angle
601  theta = (k-1)*dtheta
602  DO i = 1, local_mesh%np
603  count = count + 1
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)
607  ENDDO
608  END DO
609  ALLOCATE(mesh_3d%jj(6,mesh_3d%me)) !P1 3d VTK wedge 13
610  count = 0
611  DO k = 1, nb_angle-1
612  DO m = 1, local_mesh%me
613  count = count + 1
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
616  END DO
617  END DO
618  DO m = 1, local_mesh%me
619  count = count + 1
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)
622  END DO
623  ELSE IF (mesh%gauss%n_w==6) THEN
624  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
625  virgin = .true.
626  renumber = -1
627  count = 0
628  DO m = 1, local_mesh%me
629  DO n = 1, 3
630  i = local_mesh%jj(n,m)
631  IF (.NOT.virgin(i)) cycle
632  virgin(i) = .false.
633  count = count + 1 !===Count P1 nodes on local_mesh
634  renumber(i) = count
635  END DO
636  END DO
637  np1 = count
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)) !P2 3d VTK wedge 26
642  count = 0
643  DO k = 1, nb_angle-1
644  DO m = 1, local_mesh%me
645  count = count + 1
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
656  END DO
657  END DO
658  k = nb_angle
659  DO m = 1, local_mesh%me
660  count = count + 1
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
671  END DO
672  !===Create 3d mesh
673  DO k = 1, nb_angle
674  theta = (k-1)*dtheta
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)
680  IF (virgin(i)) cycle !===This is a vertex
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)
685  ENDDO
686  END DO
687 
688  ALLOCATE(field_out(dim,mesh_3d%np))
689  DO k = 1, nb_angle
690 !++++++++++++++++
691 !!$ theta = (k-1)*dtheta
692 !++++++++++++++++
693  DO i = 1, local_mesh%np
694  n = (k-1)*(local_mesh%np+np1) + i
695 !++++++++++++++++
696  field_out(:,n) = field_out_fft(k,:,i)
697 !!$ IF (SIZE(field_out,1)==1) THEN
698 !!$ field_out(:,n) = field_out_FFT(k,:,i)
699 !!$ ELSE
700 !!$ field_out(1,n) = field_out_FFT(k,1,i)*COS(theta)-field_out_FFT(k,2,i)*SIN(theta)
701 !!$ field_out(2,n) = field_out_FFT(k,1,i)*SIN(theta)+field_out_FFT(k,2,i)*COS(theta)
702 !!$ field_out(3,n) = field_out_FFT(k,3,i)
703 !!$ END IF
704 !++++++++++++++++
705  END DO
706  END DO
707 
708  DO i = 1, local_mesh%np
709  IF (virgin(i)) cycle !===This is a vertex
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)
714  END DO
715  END DO
716  k = nb_angle - 1
717  DO i = 1, local_mesh%np
718  IF (virgin(i)) cycle !===This is a vertex
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)
722  END DO
723  k = nb_angle
724  DO i = 1, local_mesh%np
725  IF (virgin(i)) cycle !===This is a vertex
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)
729  END DO
730  !field_out(1,:) = COS(8*ATAN2(mesh_3d%rr(2,:),mesh_3d%rr(1,:)))
731  !field_out(2,:) = mesh_3d%rr(2,:)
732  !field_out(3,:) = mesh_3d%rr(3,:)
733  DEALLOCATE(virgin, renumber)
734 
735  ELSE
736  CALL error_petsc('Bug in transfer_fourier_to_real: mesh%gauss%n_w is not correct')
737  END IF
738 
739  IF (present(opt_it)) THEN
740  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
741  field_out, name_of_field, what, opt_it=opt_it)
742  ELSE
743  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
744  field_out, name_of_field, what)
745  END IF
746 
747  !===Cleanup
748  DO n = 1, nb_procs
749  IF (ASSOCIATED(loc_to_glob(n)%DIL)) DEALLOCATE(loc_to_glob(n)%DIL)
750  END DO
751  DEALLOCATE(loc_to_glob)
752  DEALLOCATE(mesh_3d%rr, mesh_3d%jj)
753  DEALLOCATE(local_mesh%rr, local_mesh%jj)
754  DEALLOCATE(dist_field, field_out)
755 
756 ! FL-CN possible faute ?? 19/03/2013
757  !NULLIFY(field_out_FFT)
758  DEALLOCATE(field_out_fft)
759  END SUBROUTINE transfer_fourier_to_real
760 
761  SUBROUTINE divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
762  USE def_type_mesh
763  USE vtk_viz
764  IMPLICIT NONE
765  TYPE(mesh_type), INTENT(IN) :: mesh
766  TYPE(mesh_type), INTENT(OUT) :: local_mesh
767  INTEGER :: me_b, me_e, me_loc, np_loc
768  TYPE(dyn_int_line), DIMENSION(:), POINTER :: loc_to_glob
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
776 
777  !===Deal with empty mesh
778  IF (mesh%me==0) THEN
779  local_mesh%me = 0
780  local_mesh%np = 0
781  local_mesh%gauss%n_w = 0
782  ALLOCATE(local_mesh%jj(0,0))
783  ALLOCATE(local_mesh%rr(0,0))
784  ALLOCATE(loc_to_glob(0))
785  RETURN
786  END IF
787 
788  !===Continue with non empty mesh
789  CALL mpi_comm_rank(communicator,rank,ierr)
790  CALL mpi_comm_size(communicator,nb_procs,ierr)
791  ALLOCATE(loc_to_glob(nb_procs))
792  ALLOCATE(glob_to_loc(mesh%np))
793 
794  !===Quotient and remainder
795  q = mesh%me/nb_procs
796  r = mesh%me - q*nb_procs
797 
798  DO loop_rank = 0, nb_procs-1
799  !===Compute me_b and me_b
800  IF (loop_rank+1.LE.r) THEN
801  me_b = loop_rank*(q+1) + 1
802  me_e = me_b + q
803  ELSE
804  me_b = loop_rank*q + r + 1
805  me_e = me_b + q - 1
806  END IF
807  me_loc = me_e - me_b + 1
808  !===Create local mesh
809  count = 0
810  virgin = .true.
811  DO m = me_b, me_e
812  DO n = 1, mesh%gauss%n_w
813  i = mesh%jj(n,m)
814  IF (.NOT.virgin(i)) cycle
815  virgin(i)=.false.
816  count = count + 1
817  END DO
818  END DO
819  np_loc = count
820  ALLOCATE(loc_to_glob(loop_rank+1)%DIL(np_loc))
821  count = 0
822  virgin = .true.
823  DO m = me_b, me_e
824  DO n = 1, mesh%gauss%n_w
825  i = mesh%jj(n,m)
826  IF (.NOT.virgin(i)) cycle
827  virgin(i)=.false.
828  count = count + 1
829  loc_to_glob(loop_rank+1)%DIL(count) = i
830  IF (loop_rank == rank) glob_to_loc(i) = count
831  END DO
832  END DO
833 
834  IF (loop_rank == rank) THEN
835  !===Create local mesh
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))
840  DO m = 1, me_loc
841  m_glob = me_b + m -1
842  local_mesh%jj(:,m) = glob_to_loc(mesh%jj(:,m_glob))
843  END DO
844  ALLOCATE(local_mesh%rr(mesh%gauss%k_d,local_mesh%np))
845  local_mesh%rr = mesh%rr(:,loc_to_glob(rank+1)%DIL)
846  END IF
847  END DO
848 
849  DEALLOCATE(glob_to_loc)
850  END SUBROUTINE divide_mesh
851 
852  SUBROUTINE compute_grad_phi(field_in, field_out)
853  USE def_type_mesh
854  USE my_util
855  USE fem_m_axi
856  USE st_matrix
857  USE solve_petsc
858  IMPLICIT NONE
859  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field_in
860  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: field_out
861 
862  INTEGER, POINTER, DIMENSION(:) :: ifrom
863  TYPE(solver_param) :: param_grad_phi
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
871 
872 
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.
879 
880  IF (if_first_grad_phi) THEN
881  if_first_grad_phi=.false.
882 
883  CALL create_my_ghost(mesh_grad_phi, vizu_grad_phi_la,ifrom)
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)
888 
889  CALL create_local_petsc_matrix(vizu_grad_curl_comm(1),vizu_grad_phi_la , mat_grad_phi, clean=.true.)
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)
894 
895  END IF
896 
897  smb_grad_phi = 0
898  field_out = 0.d0
899  DO i = 1, SIZE(vizu_list_mode)
900  mode = vizu_list_mode(i)
901  smb_grad_phi=0.d0
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)
906  grad_phi_loc = 0.d0
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
915  END DO
916  i_loc = mesh_grad_phi%jj(:,m)
917 
918  DO k = 1, 6
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
921  END DO
922  END DO
923  END DO
924 
925  DO k = 1, 6
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)
930 
931  CALL solver(ksp_grad_phi,vb_phi,vx_phi,reinit=.false.,verbose=.false.)
932 
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))
937  END IF
938  END DO
939 
940  IF (mode == 0) THEN
941  field_out(:,2,i) = 0.d0
942  field_out(:,4,i) = 0.d0
943  field_out(:,6,i) = 0.d0
944  END IF
945 
946  END DO
947 
948 
949  END SUBROUTINE compute_grad_phi
950 
951 END MODULE fourier_to_real_for_vtu
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
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
Definition: doc_intro.h:193
subroutine loc_to_glob(mesh_loc, mesh_glob, l_t_g)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
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)
Definition: fem_M_axi.f90:127
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
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)
Definition: plot_vtk.f90:637
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)
Definition: my_util.f90:15
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:33
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:142