SFEMaNS  version 5.3
Reference documentation for SFEMaNS
fourier_real.f90
Go to the documentation of this file.
2  USE my_util
3  USE dyn_line
4  USE def_type_mesh
5  USE sub_plot
6  USE input_data
7 #include "petsc/finclude/petsc.h"
8  USE petsc
11  PUBLIC :: vtu_3d
12  TYPE(petsc_csr_la), PUBLIC :: vizu_grad_phi_la
13  TYPE(petsc_csr_la), PUBLIC :: vizu_rot_h_la
14  TYPE(petsc_csr_la), PUBLIC :: vizu_rot_u_la !ajout rot
15  TYPE(mesh_type), POINTER :: mesh_grad_phi
16  TYPE(mesh_type), POINTER :: mesh_rot_h
17  PRIVATE
18  TYPE(mesh_type), TARGET :: vv_mesh_3d
19  TYPE(mesh_type), TARGET :: pp_mesh_3d
20  TYPE(mesh_type), TARGET :: h_mesh_3d
21  TYPE(mesh_type), TARGET :: phi_mesh_3d
22  TYPE(mesh_type), TARGET :: temp_mesh_3d
23  TYPE(mesh_type), TARGET :: vv_mesh_p2_iso_p4
24  TYPE(mesh_type), TARGET :: vv_local_mesh
25  TYPE(mesh_type), TARGET :: pp_local_mesh
26  TYPE(mesh_type), TARGET :: h_local_mesh
27  TYPE(mesh_type), TARGET :: phi_local_mesh
28  TYPE(mesh_type), TARGET :: temp_local_mesh
29  TYPE(dyn_int_line), DIMENSION(:), POINTER :: vv_loc_to_glob
30  TYPE(dyn_int_line), DIMENSION(:), POINTER :: pp_loc_to_glob
31  TYPE(dyn_int_line), DIMENSION(:), POINTER :: h_loc_to_glob
32  TYPE(dyn_int_line), DIMENSION(:), POINTER :: phi_loc_to_glob
33  TYPE(dyn_int_line), DIMENSION(:), POINTER :: temp_loc_to_glob
34  INTEGER, DIMENSION(:), POINTER :: fourier_list_mode
35  INTEGER :: fourier_nb_angle
36  INTEGER :: fourier_width
37  INTEGER :: fourier_rank
38  INTEGER :: fourier_nb_procs
39  LOGICAL, SAVE :: if_first_grad_phi=.true.
40  LOGICAL, SAVE :: if_first_rot_h=.true.
41  LOGICAL, SAVE :: if_first_rot_u=.true. !===Add rot
42  INTEGER, DIMENSION(:), POINTER :: vizu_list_mode
43  mpi_comm :: fourier_communicator
44  mpi_comm, DIMENSION(2) :: vizu_grad_curl_comm
45  mat :: mat_grad_phi
46  mat :: mat_rot_h
47  mat :: mat_rot_u !===Add rot
48  vec :: vx_phi, vb_phi, vx_phi_ghost
49  ksp :: ksp_grad_phi
50  ksp :: ksp_rot_h
51  ksp :: ksp_rot_u !===Add rot
52 CONTAINS
53 
54  SUBROUTINE sfemans_initialize_postprocessing(comm_one_d, vv_mesh_in, pp_mesh, H_mesh, phi_mesh, temp_mesh, &
55  list_mode_in, opt_nb_plane)
57  USE sub_plot
58  IMPLICIT NONE
59  TYPE(mesh_type), POINTER :: vv_mesh
60  TYPE(mesh_type), TARGET, INTENT(IN) :: vv_mesh_in
61  TYPE(mesh_type), INTENT(IN) :: pp_mesh
62  TYPE(mesh_type), TARGET :: H_mesh
63  TYPE(mesh_type), TARGET :: phi_mesh
64  TYPE(mesh_type), INTENT(IN) :: temp_mesh
65  INTEGER, DIMENSION(:), TARGET :: list_mode_in
66  INTEGER, OPTIONAL, INTENT(IN) :: opt_nb_plane
67  INTEGER :: code
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)
75  fourier_rank = rank
76  fourier_nb_procs = nb_procs
77  ALLOCATE(fourier_list_mode(SIZE(list_mode_in)))
78  fourier_list_mode = list_mode_in
79  !===JLG HF July 25 2019
80  IF ((.NOT.inputs%if_just_processing) .AND. (inputs%freq_plot .GT. inputs%nb_iteration)) THEN
81  !===Not preparation of vizu meshes
82  RETURN
83  END IF
84  IF (inputs%type_fe_velocity .GE. 3) THEN
85  !===divide each cell of vv_mesh into 4 cells
87  vv_mesh => vv_mesh_p2_iso_p4
88  ELSE
89  vv_mesh => vv_mesh_in
90  END IF
91  !===JLG HF July 25 2019
92  CALL divide_mesh(comm_one_d(2), vv_mesh, vv_local_mesh, vv_loc_to_glob)
93  CALL divide_mesh(comm_one_d(2), pp_mesh, pp_local_mesh, pp_loc_to_glob)
94  CALL divide_mesh(comm_one_d(2), h_mesh, h_local_mesh, h_loc_to_glob)
95  CALL divide_mesh(comm_one_d(2), phi_mesh, phi_local_mesh, phi_loc_to_glob)
96  CALL divide_mesh(comm_one_d(2), temp_mesh, temp_local_mesh, temp_loc_to_glob)
97  CALL prepare_3d_grids(vv_mesh, vv_local_mesh, vv_mesh_3d, petsc_rank, opt_nb_plane)
98  CALL prepare_3d_grids(pp_mesh, pp_local_mesh, pp_mesh_3d, petsc_rank, opt_nb_plane)
99  CALL prepare_3d_grids(h_mesh, h_local_mesh, h_mesh_3d, petsc_rank, opt_nb_plane)
100  CALL prepare_3d_grids(phi_mesh, phi_local_mesh, phi_mesh_3d, petsc_rank, opt_nb_plane)
101  CALL prepare_3d_grids(temp_mesh, temp_local_mesh, temp_mesh_3d, petsc_rank, opt_nb_plane)
102  mesh_rot_h => h_mesh
103  mesh_grad_phi => phi_mesh
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)
106  vizu_list_mode => list_mode_in
107  END SUBROUTINE sfemans_initialize_postprocessing
108 
109  SUBROUTINE prepare_3d_grids(mesh, local_mesh, mesh_3d, rank, opt_nb_plane)
111  USE def_type_mesh
112  USE input_data
113  USE my_util
114  IMPLICIT NONE
115  TYPE(mesh_type), INTENT(IN) :: mesh
116  TYPE(mesh_type), INTENT(IN) :: local_mesh
117  TYPE(mesh_type), INTENT(OUT):: mesh_3d
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
125 
126  !===Deal with empty mesh
127  IF (mesh%me==0) THEN
128  mesh_3d%me = 0
129  mesh_3d%np = 0
130  mesh_3d%gauss%n_w = 0
131  ALLOCATE(mesh_3d%jj(0,0))
132  ALLOCATE(mesh_3d%rr(0,0))
133  RETURN
134  END IF
135 
136  !===Compute gap
137  IF (inputs%select_mode) THEN
138  IF (SIZE(inputs%list_mode_lect)==1) THEN
139  gap = inputs%list_mode_lect(1) + 1
140  ELSE
141  gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
142  END IF
143  ELSE
144  gap = 1
145  END IF
147 
148  !===Compute nb_angle
149  m_max = gap*inputs%m_max
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
153  ELSE
154  m_max_pad = m_max
155  END IF
156  ELSE
157  m_max_pad = m_max
158  END IF
159  nb_angle = 2*m_max_pad-1
160  fourier_nb_angle = nb_angle
161 
162  !===Create 3d mesh
163  dim = 3
164  pi = acos(-1.d0)
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))
170  count = 0
171  DO k = 1, nb_angle
172  theta = (k-1)*dtheta
173  DO i = 1, local_mesh%np
174  count = count + 1
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)
178  ENDDO
179  END DO
180  ALLOCATE(mesh_3d%jj(6,mesh_3d%me)) !===P1 3d VTK wedge 13
181  mesh_3d%gauss%n_w = 6
182  count = 0
183  DO k = 1, nb_angle-1
184  DO m = 1, local_mesh%me
185  count = count + 1
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
188  END DO
189  END DO
190  DO m = 1, local_mesh%me
191  count = count + 1
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)
194  END DO
195  ELSE IF (mesh%gauss%n_w==6) THEN
196  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
197  virgin = .true.
198  renumber = -1
199  count = 0
200  DO m = 1, local_mesh%me
201  DO n = 1, 3
202  i = local_mesh%jj(n,m)
203  IF (.NOT.virgin(i)) cycle
204  virgin(i) = .false.
205  count = count + 1 !===Count P1 nodes on local_mesh
206  renumber(i) = count
207  END DO
208  END DO
209  np1 = count
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)) !===P2 3d VTK wedge 26
214  mesh_3d%gauss%n_w = 15
215  count = 0
216  DO k = 1, nb_angle-1
217  DO m = 1, local_mesh%me
218  count = count + 1
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
229  END DO
230  END DO
231  k = nb_angle
232  DO m = 1, local_mesh%me
233  count = count + 1
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
244  END DO
245  !===Create 3d mesh
246  DO k = 1, nb_angle
247  theta = (k-1)*dtheta
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)
253  IF (virgin(i)) cycle !===This is a vertex
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)
258  ENDDO
259  END DO
260  !===Cleanup
261  DEALLOCATE(virgin,renumber)
262  ELSE IF (mesh%gauss%n_w==10) THEN
263  IF (rank==0) THEN
264  WRITE(*,*) 'VTK files not done for P3 meshes'
265  END IF
266  ELSE
267  CALL error_petsc('Bug in prepare_3d_grids: mesh%gauss%n_w= not correct')
268  END IF
269  END SUBROUTINE prepare_3d_grids
270 
271  SUBROUTINE vtu_3d(field, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl, opt_2D, opt_mesh_in)
273  USE def_type_mesh
274  USE vtk_viz
275  USE sft_parallele
276  USE input_data
277  USE my_util
278  IMPLICIT NONE
279  TYPE(mesh_type), OPTIONAL :: 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
288  !===VTU 2d======================================================================
289  CHARACTER(LEN=3) :: st_mode
290  CHARACTER(LEN=200) :: header_2D
291  CHARACTER(LEN=3) :: name_of_field_2D
292  INTEGER :: i
293 
294  !===HF-CN 20/01/2020
295  IF (PRESENT(opt_grad_curl)) THEN !=== compute rot_u for P3
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)))
299  CALL compute_rot_u(field, field_viz, opt_mesh_in) !=== only opt_mesh_in
300  END IF
301  END IF
302  END IF
303  !===HF-CN 20/01/2020
304 
305  !===JLG HF July 26, 2019
306  IF (PRESENT(opt_mesh_in)) THEN
307  IF (opt_mesh_in%gauss%n_w==10 .AND. inputs%type_fe_velocity==3) THEN
308  ALLOCATE(field_tmp(vv_mesh_p2_iso_p4%np,SIZE(field,2),SIZE(field,3)))
309  DO i = 1, SIZE(field,3)
310  IF (.NOT.(PRESENT(opt_grad_curl))) THEN !=== visualization velocity
311  CALL interpolate(field(:,:,i),opt_mesh_in,field_tmp(:,:,i),vv_mesh_p2_iso_p4) !===field_tmp=u
312  ELSE IF (opt_grad_curl == 'curl_u') THEN !=== visualization curl
313  CALL interpolate(field_viz(:,:,i),opt_mesh_in,field_tmp(:,:,i),vv_mesh_p2_iso_p4) !===field_tmp=curl_u
314  ELSE
315  CALL error_petsc('Bug in vtu_3d: P3 and wrong opt_grad_curl')
316  END IF
317  END DO
318  IF (ALLOCATED(field_viz)) DEALLOCATE(field_viz)
319  field_in => field_tmp !=== if P3: links to u or rot_u on vv_mesh_p2_iso_p
320  ELSE
321  field_in => field !=== not P3 but (present(opt_mesh(in)))
322  END IF
323  ELSE
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')
326  ELSE
327  field_in => field !=== not present (opt_mesh_in)
328  ENDIF
329  END IF
330  !===JLG HF July 26, 2019
331 
332  IF (PRESENT(opt_grad_curl)) THEN
333  IF (opt_grad_curl == 'grad') THEN !===grad phi
334  ALLOCATE(field_viz(SIZE(field_in,1), 3*SIZE(field_in,2), SIZE(field_in,3)))
335  CALL compute_grad_phi(field_in, field_viz)
336  ELSE IF (opt_grad_curl == 'curl_h') THEN !===curl h
337  ALLOCATE(field_viz(SIZE(field_in,1), SIZE(field_in,2), SIZE(field_in,3)))
338  CALL compute_rot_h(field_in, field_viz)
339  ELSE IF ((opt_grad_curl == 'curl_u').AND.(PRESENT(opt_mesh_in))) THEN !===curl u
340  IF (opt_mesh_in%gauss%n_w==10 .AND. inputs%type_fe_velocity==3) THEN !===P3
341  ALLOCATE(field_viz(SIZE(field_in,1), SIZE(field_in,2), SIZE(field_in,3)))
342  field_viz = field_in !===field_in = curl u
343  ELSE
344  ALLOCATE(field_viz(SIZE(field_in,1), SIZE(field_in,2), SIZE(field_in,3)))
345  CALL compute_rot_u(field_in, field_viz, opt_mesh_in) !===field_in = u
346  END IF
347  ELSE
348  CALL error_petsc('Bug in vtu_3d: name_of_opt_grad_curl is not correct')
349  END IF
350  ELSE
351  ALLOCATE(field_viz(SIZE(field_in, 1), SIZE(field_in,2), SIZE(field_in,3)))
352  field_viz = field_in
353  END IF
354 
355  IF (PRESENT(opt_it)) THEN
356  CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what, opt_it)
357  ELSE
358  CALL vtu_3d_base(field_viz, name_of_mesh, header, name_of_field, what)
359  END IF
360 
361  IF (PRESENT(opt_2d)) THEN
362  IF (opt_2d) THEN
363  DO i = 1, SIZE(vizu_list_mode)
364  WRITE(st_mode,'(I3)') vizu_list_mode(i) !=== (CHARACTER(LEN=3) :: st_mode)
365  header_2d = 'J_'//'mode_'//trim(adjustl(st_mode)) !=== (CHARACTER(LEN=200) :: header)
366  name_of_field_2d = 'J' !===(for instance) (CHARACTER(LEN=3) :: name_of_field)
367  IF (PRESENT(opt_it)) THEN
368  CALL make_vtu_file_2d(vizu_grad_curl_comm(1), mesh_rot_h, &
369  header_2d, field_viz(:,:,i), name_of_field_2d, what, opt_it=opt_it)
370  ELSE
371  CALL make_vtu_file_2d(vizu_grad_curl_comm(1), mesh_rot_h, &
372  header_2d, field_viz(:,:,i), name_of_field_2d, what)
373  END IF
374  END DO
375  END IF
376  END IF
377 
378  DEALLOCATE(field_viz)
379  IF (ALLOCATED(field_tmp)) DEALLOCATE(field_tmp)
380  END SUBROUTINE vtu_3d
381 
382  SUBROUTINE vtu_3d_base(field_in, name_of_mesh, header, name_of_field, what, opt_it)
384  USE def_type_mesh
385  USE vtk_viz
386  USE sft_parallele
387  USE input_data
388  USE my_util
389  IMPLICIT NONE
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
393  TYPE(mesh_type), POINTER :: mesh_3d
394  TYPE(mesh_type), POINTER :: local_mesh
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
408 
409  !===Initialization
410  width = fourier_width
411  nb_procs = fourier_nb_procs
412  rank = fourier_rank
413  nb_angle = fourier_nb_angle
414  list_mode => fourier_list_mode
415  ALLOCATE(loc_to_glob(nb_procs))
416  IF (name_of_mesh=='vv_mesh') THEN
417  DO n = 1, nb_procs
418  loc_to_glob(n)%DIL => vv_loc_to_glob(n)%DIL
419  END DO
420  mesh_3d => vv_mesh_3d
421  local_mesh => vv_local_mesh
422  ELSE IF (name_of_mesh=='pp_mesh') THEN
423  DO n = 1, nb_procs
424  loc_to_glob(n)%DIL => pp_loc_to_glob(n)%DIL
425  END DO
426  mesh_3d => pp_mesh_3d
427  local_mesh => pp_local_mesh
428  ELSE IF (name_of_mesh=='H_mesh') THEN
429  DO n = 1, nb_procs
430  loc_to_glob(n)%DIL => h_loc_to_glob(n)%DIL
431  END DO
432  mesh_3d => h_mesh_3d
433  local_mesh => h_local_mesh
434  ELSE IF (name_of_mesh=='phi_mesh') THEN
435  DO n = 1, nb_procs
436  loc_to_glob(n)%DIL => phi_loc_to_glob(n)%DIL
437  END DO
438  mesh_3d => phi_mesh_3d
439  local_mesh => phi_local_mesh
440  ELSE IF (name_of_mesh=='temp_mesh') THEN
441  DO n = 1, nb_procs
442  loc_to_glob(n)%DIL => temp_loc_to_glob(n)%DIL
443  END DO
444  mesh_3d => temp_mesh_3d
445  local_mesh => temp_local_mesh
446  ELSE
447  CALL error_petsc('Bug in vtu_3d: name_of_mesh is not correct')
448  END IF
449 
450  !===Compute np_loc_max
451  np_loc_max = -1
452  DO n = 1, nb_procs
453  np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
454  END DO
455  ALLOCATE(dist_field(nb_procs*np_loc_max,SIZE(field_in,2),width))
456 
457  !===Navier subdomain is divided into fourier_rank parts
458  !===This loop ditributes fourier modes to the sub-meshes (1 to fourier_rank)
459  dist_field = 0.d0
460  mode_min = rank*width
461  DO i = 1, 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)
465  DO n = 1, nb_procs
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))
470  END DO
471  END DO
472  END DO
473  !===Compute field_out_FFT in real space on local_mesh
474  CALL fft_par_real(fourier_communicator, dist_field, field_out_fft, opt_nb_plane=nb_angle)
475 
476  !===FL 29/03/2013
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
480  ELSE
481  pi = acos(-1.d0)
482  dtheta = 2*pi/nb_angle
483  DO k = 1, nb_angle
484  theta = (k-1)*dtheta
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,:)
488  END DO
489  END IF
490  !===FL 29/03/2013
491 
492  IF (SIZE(field_in,2)==2) THEN
493  dim = 1
494  ELSE IF (SIZE(field_in,2)==6) THEN
495  dim = 3
496  ELSE
497  CALL error_petsc('Bug in vtu_3d: SIZE(field_in,2) not correct')
498  END IF
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')
501  END IF
502 
503  !===Reconstruct field_out in real space on mesh_3d
504  pi = acos(-1.d0)
505  dtheta = 2*pi/nb_angle
506  IF (local_mesh%gauss%n_w==3) THEN
507  ALLOCATE(field_out(dim,mesh_3d%np))
508  count = 0
509  DO k = 1, nb_angle
510  DO i = 1, local_mesh%np
511  count = count + 1
512  field_out(:,count) = field_out_xyz_fft(k,:,i)
513  ENDDO
514  END DO
515 
516  ELSE IF (local_mesh%gauss%n_w==6) THEN
517  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
518  virgin = .true.
519  renumber = -1
520  count = 0
521  DO m = 1, local_mesh%me
522  DO n = 1, 3
523  i = local_mesh%jj(n,m)
524  IF (.NOT.virgin(i)) cycle
525  virgin(i) = .false.
526  count = count + 1 !===Count P1 nodes on local_mesh
527  renumber(i) = count
528  END DO
529  END DO
530  np1 = count
531 
532  ALLOCATE(field_out(dim,mesh_3d%np))
533  DO k = 1, nb_angle
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)
537  END DO
538  END DO
539 
540  DO i = 1, local_mesh%np
541  IF (virgin(i)) cycle !===This is a vertex
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)
546  END DO
547  END DO
548  k = nb_angle - 1
549  DO i = 1, local_mesh%np
550  IF (virgin(i)) cycle !===This is a vertex
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)
554  END DO
555  k = nb_angle
556  DO i = 1, local_mesh%np
557  IF (virgin(i)) cycle !===This is a vertex
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)
561  END DO
562  DEALLOCATE(virgin, renumber)
563  ELSE
564  CALL error_petsc('Bug in vtu_3d: mesh%gauss%n_w is not correct')
565  END IF
566 
567  IF (PRESENT(opt_it)) THEN
568  it = opt_it
569  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
570  field_out, name_of_field, what, opt_it = it)
571  ELSE
572  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
573  field_out, name_of_field, what)
574  END IF
575 
576  !===Cleanup
577  DEALLOCATE(dist_field, field_out, loc_to_glob,field_out_xyz_fft)
578  DEALLOCATE(field_out_fft)
579  END SUBROUTINE vtu_3d_base
580 
581  SUBROUTINE transfer_fourier_to_real(communicator, mesh, field_in, &
582  header, name_of_field, what, list_mode, opt_it)
584  USE def_type_mesh
585  USE vtk_viz
586  USE sft_parallele
587  USE input_data
588  USE my_util
589  IMPLICIT NONE
590  TYPE(mesh_type), INTENT(IN) :: mesh
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
595  TYPE(mesh_type) :: mesh_3d
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
601  TYPE(mesh_type) :: local_mesh
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)
613 
614  IF (inputs%select_mode) THEN
615  IF (SIZE(inputs%list_mode_lect)==1) THEN
616  gap = inputs%list_mode_lect(1) + 1
617  ELSE
618  gap = inputs%list_mode_lect(2) - inputs%list_mode_lect(1)
619  END IF
620  ELSE
621  gap = 1
622  END IF
623  width = SIZE(list_mode)*gap
624 
625  CALL divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
626  np_loc_max = -1
627  DO n = 1, nb_procs
628  np_loc_max = max(np_loc_max,maxval(loc_to_glob(n)%DIL))
629  END DO
630  ALLOCATE(dist_field(nb_procs*np_loc_max,SIZE(field_in,2),width))
631  dist_field = 0.d0
632  mode_min = rank*width
633  DO i = 1, 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)
637  DO n = 1, nb_procs
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))
642  END DO
643  END DO
644  END DO
645 
646  CALL fft_par_real(communicator, dist_field, field_out_fft, opt_nb_plane=50)
647 
648  IF (SIZE(field_in,2)==2) THEN
649  dim = 1
650  ELSE IF (SIZE(field_in,2)==6) THEN
651  dim = 3
652  ELSE
653  CALL error_petsc('Bug in transfer_fourier_to_real: SIZE(field_in,2) not correct')
654  END IF
655  nb_angle = SIZE(field_out_fft,1)
656  pi = acos(-1.d0)
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))
662  count = 0
663  DO k = 1, nb_angle
664  DO i = 1, local_mesh%np
665  count = count + 1
666  field_out(:,count) = field_out_fft(k,:,i)
667  ENDDO
668  END DO
669  !===Create 3d mesh
670  ALLOCATE(mesh_3d%rr(dim,mesh_3d%np))
671  count = 0
672  DO k = 1, nb_angle
673  theta = (k-1)*dtheta
674  DO i = 1, local_mesh%np
675  count = count + 1
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)
679  ENDDO
680  END DO
681  ALLOCATE(mesh_3d%jj(6,mesh_3d%me)) !===P1 3d VTK wedge 13
682  mesh_3d%gauss%n_w = 6
683  count = 0
684  DO k = 1, nb_angle-1
685  DO m = 1, local_mesh%me
686  count = count + 1
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
689  END DO
690  END DO
691  DO m = 1, local_mesh%me
692  count = count + 1
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)
695  END DO
696  ELSE IF (mesh%gauss%n_w==6) THEN
697  ALLOCATE(virgin(local_mesh%np), renumber(local_mesh%np))
698  virgin = .true.
699  renumber = -1
700  count = 0
701  DO m = 1, local_mesh%me
702  DO n = 1, 3
703  i = local_mesh%jj(n,m)
704  IF (.NOT.virgin(i)) cycle
705  virgin(i) = .false.
706  count = count + 1 !===Count P1 nodes on local_mesh
707  renumber(i) = count
708  END DO
709  END DO
710  np1 = count
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)) !===P2 3d VTK wedge 26
715  mesh_3d%gauss%n_w = 15
716  count = 0
717  DO k = 1, nb_angle-1
718  DO m = 1, local_mesh%me
719  count = count + 1
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
730  END DO
731  END DO
732  k = nb_angle
733  DO m = 1, local_mesh%me
734  count = count + 1
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
745  END DO
746  !===Create 3d mesh
747  DO k = 1, nb_angle
748  theta = (k-1)*dtheta
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)
754  IF (virgin(i)) cycle !===This is a vertex
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)
759  ENDDO
760  END DO
761 
762  ALLOCATE(field_out(dim,mesh_3d%np))
763  DO k = 1, nb_angle
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)
767  END DO
768  END DO
769 
770  DO i = 1, local_mesh%np
771  IF (virgin(i)) cycle !===This is a vertex
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)
776  END DO
777  END DO
778  k = nb_angle - 1
779  DO i = 1, local_mesh%np
780  IF (virgin(i)) cycle !===This is a vertex
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)
784  END DO
785  k = nb_angle
786  DO i = 1, local_mesh%np
787  IF (virgin(i)) cycle !===This is a vertex
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)
791  END DO
792  DEALLOCATE(virgin, renumber)
793 
794  ELSE
795  CALL error_petsc('Bug in transfer_fourier_to_real: mesh%gauss%n_w is not correct')
796  END IF
797 
798  IF (PRESENT(opt_it)) THEN
799  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
800  field_out, name_of_field, what, opt_it=opt_it)
801  ELSE
802  CALL make_vtu_file_3d(mpi_comm_world, mesh_3d, header, &
803  field_out, name_of_field, what)
804  END IF
805 
806  !===Cleanup
807  DO n = 1, nb_procs
808  IF (ASSOCIATED(loc_to_glob(n)%DIL)) DEALLOCATE(loc_to_glob(n)%DIL)
809  END DO
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)
815  END SUBROUTINE transfer_fourier_to_real
816 
817  SUBROUTINE divide_mesh(communicator, mesh, local_mesh, loc_to_glob)
819  USE vtk_viz
820  IMPLICIT NONE
821  TYPE(mesh_type), INTENT(IN) :: mesh
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
831 
832  !===Deal with empty mesh
833  IF (mesh%me==0) THEN
834  local_mesh%me = 0
835  local_mesh%np = 0
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))
840  RETURN
841  END IF
842 
843  !===Continue with non empty mesh
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))
848 
849  !===Quotient and remainder
850  q = mesh%me/nb_procs
851  r = mesh%me - q*nb_procs
852 
853  DO loop_rank = 0, nb_procs-1
854  !===Compute me_b and me_b
855  IF (loop_rank+1.LE.r) THEN
856  me_b = loop_rank*(q+1) + 1
857  me_e = me_b + q
858  ELSE
859  me_b = loop_rank*q + r + 1
860  me_e = me_b + q - 1
861  END IF
862  me_loc = me_e - me_b + 1
863  !===Create local mesh
864  count = 0
865  virgin = .true.
866  DO m = me_b, me_e
867  DO n = 1, mesh%gauss%n_w
868  i = mesh%jj(n,m)
869  IF (.NOT.virgin(i)) cycle
870  virgin(i)=.false.
871  count = count + 1
872  END DO
873  END DO
874  np_loc = count
875  ALLOCATE(loc_to_glob(loop_rank+1)%DIL(np_loc))
876  count = 0
877  virgin = .true.
878  DO m = me_b, me_e
879  DO n = 1, mesh%gauss%n_w
880  i = mesh%jj(n,m)
881  IF (.NOT.virgin(i)) cycle
882  virgin(i)=.false.
883  count = count + 1
884  loc_to_glob(loop_rank+1)%DIL(count) = i
885  IF (loop_rank == rank) glob_to_loc(i) = count
886  END DO
887  END DO
888 
889  IF (loop_rank == rank) THEN
890  !===Create local mesh
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))
895  DO m = 1, me_loc
896  m_glob = me_b + m -1
897  local_mesh%jj(:,m) = glob_to_loc(mesh%jj(:,m_glob))
898  END DO
899  ALLOCATE(local_mesh%rr(2,local_mesh%np))
900  local_mesh%rr = mesh%rr(:,loc_to_glob(rank+1)%DIL)
901  END IF
902  END DO
903 
904  DEALLOCATE(glob_to_loc)
905  END SUBROUTINE divide_mesh
906 
907  SUBROUTINE compute_grad_phi(field_in, field_out)
909  USE my_util
910  USE fem_m_axi
911  USE st_matrix
912  USE solve_petsc
913  IMPLICIT NONE
914  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field_in
915  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: field_out
916 
917  INTEGER, POINTER, DIMENSION(:) :: ifrom
918  TYPE(solver_param) :: param_grad_phi
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
925 
926 
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.
933 
934  IF (if_first_grad_phi) THEN
935  if_first_grad_phi=.false.
936 
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)
942 
943  CALL create_local_petsc_matrix(vizu_grad_curl_comm(1),vizu_grad_phi_la , mat_grad_phi, clean=.true.)
944  CALL qs_diff_mass_scal_m (mesh_grad_phi, vizu_grad_phi_la, 0.d0, 1.d0, 0.d0, 0, mat_grad_phi)
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)
948 
949  END IF
950 
951  smb_grad_phi = 0
952  field_out = 0.d0
953  DO i = 1, SIZE(vizu_list_mode)
954  mode = vizu_list_mode(i)
955  smb_grad_phi=0.d0
956  DO m = 1, mesh_grad_phi%me
957  DO l = 1, mesh_grad_phi%gauss%l_G
958  ray = sum(mesh_grad_phi%rr(1,mesh_grad_phi%jj(:,m))*mesh_grad_phi%gauss%ww(:,l))
959  drdthetadz = mesh_grad_phi%gauss%rj(l,m)
960  grad_phi_loc = 0.d0
961  DO nj = 1,mesh_grad_phi%gauss%n_w
962  j_loc = mesh_grad_phi%jj(nj,m)
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
969  END DO
970  i_loc = mesh_grad_phi%jj(:,m)
971 
972  DO k = 1, 6
973  smb_grad_phi(i_loc,k) = smb_grad_phi(i_loc,k) + grad_phi_loc(k)* &
974  mesh_grad_phi%gauss%ww(:,l)*drdthetadz
975  END DO
976  END DO
977  END DO
978 
979  DO k = 1, 6
980  CALL veczeroentries(vb_phi, ierr)
981  CALL vecsetvalues(vb_phi, mesh_grad_phi%np, vizu_grad_phi_la%loc_to_glob(1,:)-1, smb_grad_phi(:,k), add_values, ierr)
982  CALL vecassemblybegin(vb_phi,ierr)
983  CALL vecassemblyend(vb_phi,ierr)
984 
985  CALL solver(ksp_grad_phi,vb_phi,vx_phi,reinit=.false.,verbose=.false.)
986 
987  CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
988  CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
989  IF (mesh_grad_phi%me/=0) THEN
990  CALL extract(vx_phi_ghost,1,1,vizu_grad_phi_la,field_out(:,k,i))
991  END IF
992  END DO
993 
994  IF (mode == 0) THEN
995  field_out(:,2,i) = 0.d0
996  field_out(:,4,i) = 0.d0
997  field_out(:,6,i) = 0.d0
998  END IF
999 
1000  END DO
1001 
1002  END SUBROUTINE compute_grad_phi
1003 
1004  SUBROUTINE compute_rot_h(field_in, field_out)
1006  USE my_util
1007  USE fem_m_axi
1008  USE st_matrix
1009  USE solve_petsc
1010  IMPLICIT NONE
1011  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field_in
1012  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: field_out
1013 
1014  INTEGER, POINTER, DIMENSION(:) :: ifrom
1015  TYPE(solver_param) :: param_rot_h
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
1022 
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.
1029 
1030  IF (if_first_rot_h) THEN
1031  if_first_rot_h=.false.
1032 
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)
1038 
1039  CALL create_local_petsc_matrix(vizu_grad_curl_comm(1),vizu_rot_h_la , mat_rot_h, clean=.true.)
1040  CALL qs_diff_mass_scal_m (mesh_rot_h, vizu_rot_h_la, 0.d0, 1.d0, 0.d0, 0, mat_rot_h)
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)
1044 
1045  END IF
1046 
1047  smb_rot_h = 0
1048  field_out = 0.d0
1049  DO i = 1, SIZE(vizu_list_mode)
1050  mode = vizu_list_mode(i)
1051  smb_rot_h=0.d0
1052  DO m = 1, mesh_rot_h%me
1053  DO l = 1, mesh_rot_h%gauss%l_G
1054  ray = sum(mesh_rot_h%rr(1,mesh_rot_h%jj(:,m))*mesh_rot_h%gauss%ww(:,l))
1055  drdthetadz = mesh_rot_h%gauss%rj(l,m)
1056  rot_h_loc = 0.d0
1057  DO nj = 1,mesh_rot_h%gauss%n_w
1058  j_loc = mesh_rot_h%jj(nj,m)
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)
1070 
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)
1074  END DO !=== sum on nj
1075  i_loc = mesh_rot_h%jj(:,m)
1076 
1077  DO k = 1, 6
1078  smb_rot_h(i_loc,k) = smb_rot_h(i_loc,k) + rot_h_loc(k)* &
1079  mesh_rot_h%gauss%ww(:,l)*drdthetadz
1080  END DO
1081  END DO
1082  END DO
1083 
1084  DO k = 1, 6
1085  CALL veczeroentries(vb_phi, ierr)
1086  CALL vecsetvalues(vb_phi, mesh_rot_h%np, vizu_rot_h_la%loc_to_glob(1,:)-1, smb_rot_h(:,k), add_values, ierr)
1087  CALL vecassemblybegin(vb_phi,ierr)
1088  CALL vecassemblyend(vb_phi,ierr)
1089 
1090  CALL solver(ksp_rot_h,vb_phi,vx_phi,reinit=.false.,verbose=.false.)
1091 
1092  CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr)
1093  CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
1094  IF (mesh_rot_h%me/=0) THEN
1095  CALL extract(vx_phi_ghost,1,1,vizu_rot_h_la,field_out(:,k,i))
1096  END IF
1097  END DO
1098 
1099  IF (mode == 0) THEN
1100  field_out(:,2,i) = 0.d0
1101  field_out(:,4,i) = 0.d0
1102  field_out(:,6,i) = 0.d0
1103  END IF
1104 
1105  END DO !===end loop on mode
1106 
1107  END SUBROUTINE compute_rot_h
1108 
1109  SUBROUTINE compute_rot_u(field_in, field_out, mesh_in)
1111  USE my_util
1112  USE fem_m_axi
1113  USE st_matrix
1114  USE solve_petsc
1115  IMPLICIT NONE
1116  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field_in
1117  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: field_out
1118  TYPE(mesh_type), INTENT(IN) :: mesh_in
1119 
1120  INTEGER, POINTER, DIMENSION(:) :: ifrom
1121  TYPE(solver_param) :: param_rot_u
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
1128 
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.
1135 
1136  IF (if_first_rot_u) THEN
1137  if_first_rot_u=.false.
1138 
1139  CALL create_my_ghost(mesh_in, vizu_rot_u_la, ifrom)
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)
1144 
1145  CALL create_local_petsc_matrix(vizu_grad_curl_comm(1),vizu_rot_u_la , mat_rot_u, clean=.true.)
1146  CALL qs_diff_mass_scal_m (mesh_in, vizu_rot_u_la, 0.d0, 1.d0, 0.d0, 0, mat_rot_u)
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)
1150  END IF
1151 
1152  smb_rot_u = 0
1153  field_out = 0.d0
1154  DO i = 1, SIZE(vizu_list_mode)
1155  mode = vizu_list_mode(i)
1156  smb_rot_u=0.d0
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)
1161  rot_u_loc = 0.d0
1162  DO nj = 1,mesh_in%gauss%n_w
1163  j_loc = mesh_in%jj(nj,m)
1164 
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
1169 
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
1174 
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)
1181  END DO !=== sum on nj
1182  i_loc = mesh_in%jj(:,m)
1183 
1184  DO k = 1, 6
1185  smb_rot_u(i_loc,k) = smb_rot_u(i_loc,k) + rot_u_loc(k)* &
1186  mesh_in%gauss%ww(:,l)*drdthetadz
1187  END DO
1188  END DO
1189  END DO
1190 
1191  DO k = 1, 6
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)
1196 
1197  CALL solver(ksp_rot_u,vb_phi,vx_phi,reinit=.false.,verbose=.false.)
1198 
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
1202  CALL extract(vx_phi_ghost,1,1,vizu_rot_u_la,field_out(:,k,i))
1203  END IF
1204  END DO
1205 
1206  IF (mode == 0) THEN
1207  field_out(:,2,i) = 0.d0
1208  field_out(:,4,i) = 0.d0
1209  field_out(:,6,i) = 0.d0
1210  END IF
1211 
1212  END DO !===end loop on mode
1213 
1214  END SUBROUTINE compute_rot_u
1215 
1216  SUBROUTINE divide_mesh_into_four_subcells(mesh_in,mesh_out)
1217  ! JLG + HF july 25, 2019
1218  ! Scheme of the sub_mesh partition
1219  ! n
1220  ! |\
1221  ! |n \
1222  ! | \
1223  ! | \n1+3
1224  ! n2+3- \
1225  ! | \
1226  ! |4*(m-1)+n \
1227  ! | \
1228  ! |n1 n+3 n2\ this is the element m (large mesh)
1229  ! |--------|--------|\ divided in 4*(m-1)+1, 4*(m-1)+2,
1230  ! |\n2 n+3 n1|n \ 4*(m-1)+3, 4*m in small mesh
1231  ! |n\ | \
1232  ! | \ 4m | \
1233  ! | \n1+3 | \n1+3
1234  ! | \ n2+3-n2+3 \
1235  ! n2+3- n1+3 \ | \
1236  ! | \ | \
1237  ! | \ | \
1238  ! | 4*(m-1)+n1 \n| 4*(m-1)+n2 \
1239  ! |n1 n2\|n1 n2\
1240  ! --------|--------------------|-----------
1241  ! n1 n+3 n+3 n2
1242  ! same edge touching (separed for clarity)
1243  ! n1op nop+3 nop+3 n2op
1244  ! ---------|--------------------|----------
1245  ! |n1op n2op /|n1op n2op /
1246  ! | 4(mop-1)+n1op / | 4(mop-1)+n2op /
1247  ! | /nop| /
1248  ! | / | /
1249  ! | / | /
1250  ! | / | /
1251  ! | / | /
1252  ! | / 4mop | /
1253  ! | / | /
1254  ! |/ | / this is the neigh element mop (large mesh)
1255  ! |-----------------|/ divided in 4*(mop-1)+1, 4*(mop-1)+2,
1256  ! | / 4*(mop-1)+3, 4*mop in small mesh
1257  ! | /
1258  ! | /
1259  ! | /
1260  ! | /
1261  ! | /
1262  ! | /
1263  ! | /
1264  ! |/
1265  ! nop
1266  IMPLICIT NONE
1267  TYPE(mesh_type), INTENT(IN) :: mesh_in
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, &
1274  j1, j2
1275  mesh_out%gauss%k_d = 2
1276  mesh_out%gauss%n_w = 6
1277  mesh_out%me = 4*mesh_in%me
1278 
1279  !===Nb of edges
1280  edge = 0
1281  DO m = 1, mesh_in%me
1282  DO n = 1, 3
1283  IF (mesh_in%neigh(n,m).LE.0) cycle !===ATTENTION, no dg mesh allowed here
1284  edge = edge + 1
1285  END DO
1286  END DO
1287  edge = edge/2
1288  DO m = 1, mesh_in%me
1289  DO n = 1, 3
1290  IF (mesh_in%neigh(n,m).GE.0) cycle !===ATTENTION, no dg mesh allowed here
1291  edge = edge + 1
1292  END DO
1293  END DO
1294  nb_edges = edge + mesh_in%mes
1295 
1296  IF (mesh_in%gauss%n_w==10) THEN
1297  nb_vertices = mesh_in%np - 2*nb_edges - mesh_in%me
1298  ELSE
1299  CALL error_petsc('BUG in divide_mesh_into_four_subcells: mesh to be divided is not P3')
1300  END IF
1301  mesh_out%np = nb_vertices + 3*nb_edges + 3*mesh_in%me !===Create P2-iso-P4 mesh
1302 
1303  !===Allocation
1304  ALLOCATE(mesh_out%jj(6,mesh_out%me), mesh_out%rr(2,mesh_out%np))
1305  v_virgin = .true.
1306  n_dof = 0
1307  DO m = 1, mesh_in%me
1308  DO n = 1, 3
1309  j = mesh_in%jj(n,m)
1310  IF (v_virgin(j)) THEN
1311  v_virgin(j) = .false.
1312  n_dof = n_dof + 1
1313  v_in_to_out(j) = n_dof
1314  END IF
1315  END DO
1316  END DO
1317 
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')
1321  END IF
1322 
1323  DO m = 1, mesh_in%me
1324  DO n = 1, 3 !===take care of vertices
1325  m_out = 4*(m-1)+n
1326  j = mesh_in%jj(n,m)
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)
1329  END DO
1330  END DO
1331 
1332  c_virgin = .true.
1333  n_dof = nb_vertices
1334  DO m = 1, mesh_in%me
1335  c_virgin(m)=.false.
1336  DO n = 1, 3 !===We take care of faces and interfaces
1337  m_op = mesh_in%neigh(n,m)
1338  IF (m_op.GT.0) THEN
1339  IF (.NOT.c_virgin(m_op)) cycle
1340  END IF
1341  IF (m_op.LE.0) THEN
1342  n_op = n !===Do not create the nodes on the opposite cell
1343  m_op = m
1344  ELSE
1345  DO n_op = 1, 3
1346  IF (mesh_in%neigh(n_op,m_op)==m) EXIT !===n_op is the local index of the opposite vertex
1347  END DO
1348  END IF
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
1354  nn = n1_op
1355  n1_op = n2_op
1356  n2_op = nn
1357  END IF
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')
1360  END IF
1361 
1362  n_dof = n_dof + 1 !===New vertex of P2 sub-mesh
1363  j1=mesh_in%jj(n1,m)
1364  j2=mesh_in%jj(n2,m)
1365  mesh_out%rr(:,n_dof) = (mesh_in%rr(:,j1)+mesh_in%rr(:,j2))/2
1366 
1367  m_out = 4*(m-1) +n1
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
1371 
1372  m_out = 4*(m-1) +n2
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
1376 
1377  m_out = 4*m
1378  m_out_op = 4*m_op
1379  mesh_out%jj(n,m_out) = n_dof
1380  mesh_out%jj(n_op,m_out_op) = n_dof
1381 
1382  n_dof = n_dof +1 !===New midpoint of P2 sub-mesh
1383  mesh_out%rr(:,n_dof) = mesh_in%rr(:,j1)*3/4 + mesh_in%rr(:,j2)/4
1384 
1385  m_out = 4*(m-1) +n1
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
1389 
1390  n_dof = n_dof +1 !===New midpoint of P2 sub-mesh
1391  mesh_out%rr(:,n_dof) = mesh_in%rr(:,j1)/4 + mesh_in%rr(:,j2)*3/4
1392 
1393  m_out = 4*(m-1) +n2
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
1397 
1398  END DO
1399 
1400  !===Now we take care of the dofs inside the cell m
1401  DO n = 1, 3
1402  n1 = modulo(n,3) + 1
1403  n2 = modulo(n+1,3) + 1
1404  j1=mesh_in%jj(n1,m)
1405  j2=mesh_in%jj(n2,m)
1406  j =mesh_in%jj(n ,m)
1407  n_dof = n_dof + 1
1408  mesh_out%rr(:,n_dof) = mesh_in%rr(:,j)/2 + mesh_in%rr(:,j1)/4 + mesh_in%rr(:,j2)/4
1409 
1410  m_out = 4*m
1411  mesh_out%jj(n+3,m_out) = n_dof
1412  m_out = 4*(m-1) + n
1413  mesh_out%jj(n+3,m_out) = n_dof
1414  END DO
1415  END DO
1416 
1417  END SUBROUTINE divide_mesh_into_four_subcells
1418 
1419  SUBROUTINE interpolate(field_in,mesh_in,field_out,mesh_out)
1420  IMPLICIT NONE
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
1444 
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))
1448  a = rr(:,2)-rr(:,1)
1449  b = rr(:,3)-rr(:,1)
1450  mesk = mesurek(a,b)
1451  DO nout = 1, 6
1452  rout = mesh_out%rr(:,mesh_out%jj(nout,m_out))
1453  DO nin = 1, 3
1454  n1in = modulo(nin,3) + 1
1455  n2in = modulo(nin+1,3) + 1
1456  a = rr(:,n1in)-rout
1457  b = rr(:,n2in)-rout
1458  mesin = mesurek(a,b)
1459  lambda(nin) = mesin/mesk
1460  END DO
1461  IF (abs(sum(lambda)-1.d0).GT.1.d-13) THEN
1462  CALL error_petsc('BUG in interpolate')
1463  END IF
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))
1475  END DO
1476  END DO
1477  CONTAINS
1478  FUNCTION mesurek(a,b) RESULT(mes)
1479  REAL(KIND=8), DIMENSION(2), INTENT(IN) :: a, b
1480  REAL(KIND=8) :: mes
1481  mes = abs(a(1)*b(2)-a(2)*b(1))
1482  END FUNCTION mesurek
1483 
1484  END SUBROUTINE interpolate
1485 
1486  SUBROUTINE make_vtu_file_2d(communicator, mesh_in, header, field_in, field_name, what, opt_it)
1488  USE my_util
1489  USE vtk_viz
1490  IMPLICIT NONE
1491  TYPE(mesh_type), INTENT(IN), TARGET :: mesh_in
1492  TYPE(mesh_type), POINTER :: mesh
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
1499  INTEGER :: j, it
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))
1508 
1509  IF (mesh_in%gauss%n_w==10 .AND. inputs%type_fe_velocity==3) THEN
1510  ALLOCATE(field_tmp(vv_mesh_p2_iso_p4%np,SIZE(field_in,2)))
1511  CALL interpolate(field_in,mesh_in,field_tmp,vv_mesh_p2_iso_p4)
1512  field => field_tmp
1513  mesh => vv_mesh_p2_iso_p4
1514  ELSE
1515  field => field_in
1516  mesh => mesh_in
1517  END IF
1518 
1519  IF (PRESENT(opt_it)) THEN
1520  it = opt_it
1521  WRITE(st_it,'(I3)') it
1522  DO j = 1, nb_procs
1523  WRITE(st_rank,'(I3)') j
1524  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))//&
1525  '_it_'//trim(adjustl(st_it))
1526  END DO
1527  ELSE
1528  DO j = 1, nb_procs
1529  WRITE(st_rank,'(I3)') j
1530  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
1531  END DO
1532  END IF
1533 
1534  CALL check_list(communicator, file_list, mesh%np)
1535  IF (rank==0) THEN
1536  IF (PRESENT(opt_it)) THEN
1537  it = opt_it
1538  ELSE
1539  it = 1
1540  END IF
1541  CALL create_pvd_file(file_list, trim(header), it, trim(what))
1542  END IF
1543 
1544  IF (SIZE(field,2) == 6) THEN
1545  CALL create_xml_vtu_vect_file(field, mesh, trim(adjustl(file_list(rank+1))), field_name)
1546  ELSE IF (SIZE(field,2) == 1) THEN
1547  CALL create_xml_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(rank+1))), field_name)
1548  ELSE IF (SIZE(field,2) .GT. 0) THEN
1549  CALL create_xml_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(rank+1))), field_name)
1550  ELSE
1551  CALL error_petsc('Bug in make_vtu_file_2D: field needs at least one component')
1552  END IF
1553 
1554  IF (ALLOCATED(field_tmp)) DEALLOCATE(field_tmp)
1555  END SUBROUTINE make_vtu_file_2d
1556 
1557 END MODULE fourier_to_real_for_vtu
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)
Definition: solver.f90:99
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)
Definition: st_csr.f90:34
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:15
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)
Definition: plot_vtk.f90:14
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)
Definition: solver.f90:147
subroutine error_petsc(string)
Definition: my_util.f90:16
type(my_data), public inputs
logical, save if_first_grad_phi
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
Definition: fem_M_axi.f90:130
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)
Definition: plot_vtk.f90:82
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)
Definition: solver.f90:12
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)
Definition: plot_vtk.f90:47
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)
Definition: plot_vtk.f90:546
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)
Definition: plot_vtk.f90:225
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