SFEMaNS  version 5.3
Reference documentation for SFEMaNS
template_arpack_mhd.f90
Go to the documentation of this file.
1 MODULE arpack_mhd
2 
3  PUBLIC :: solver_arpack_mhd, post_proc_arpack
4  PRIVATE
5 
6 CONTAINS
7 
8  SUBROUTINE solver_arpack_mhd(communicator,H_mesh,phi_mesh,&
9  dt,list_mode,mu_h_field)
10  USE chaine_caractere
11  USE def_type_mesh
12  USE initialization
13  USE my_util
14  USE input_data
15 #include "petsc/finclude/petsc.h"
16  USE petsc
17  IMPLICIT NONE
18  TYPE(mesh_type), INTENT(IN) :: H_mesh, phi_mesh
19  REAL(KIND=8), INTENT(IN) :: dt
20  INTEGER, POINTER, DIMENSION(:),INTENT(IN) :: list_mode
21  REAL(KIND=8),DIMENSION(:), INTENT(IN) :: mu_H_field
22  !Arpack---------------------------------------------------------------------
23  CHARACTER(len=3) :: arpack_type
24  INTEGER :: code, rank
25  !---------------------------------------------------------------------------
26  mpi_comm, DIMENSION(2) :: communicator
27 
28  !------------------------ARPACK OR NOT? THAT IS THE QUESTION---------------
29  arpack_type = 'max' !Navier-Stokes not programmed yet
30  !---------------------------------------------------------------------------
31 
32  !------------------------ARPACK RESOLUTION----------------------------------
33  IF (inputs%if_arpack) THEN
34  IF (arpack_type=='nst') THEN
35  CALL arpack_navier_stokes
36  ELSE IF (arpack_type=='max') THEN
37  CALL arpack_maxwell_int_by_parts(communicator, h_mesh,phi_mesh, &
38  dt,list_mode,mu_h_field)
39  END IF
40  !----------------SAUVEGARDE-------------------------------------------------
41  CALL mpi_comm_rank(mpi_comm_world,rank,code)
42  CALL mpi_barrier(mpi_comm_world,code)
43  !save_run does not work here
44  !CALL save_run(0,1)
45  !---------------------------------------------------------------------------
46  END IF
47  !---------------------------------------------------------------------------
48 
49  RETURN
50  END SUBROUTINE solver_arpack_mhd
51 
52  !---------------------------------------------------------------------------
53  SUBROUTINE arpack_maxwell_int_by_parts(communicator, H_mesh, phi_mesh, dt, list_mode, &
54  mu_h_field)
55  USE initialization
56  USE def_type_mesh
57  USE chaine_caractere
58  USE vtk_viz
59  USE tn_axi
60  USE input_data
61  USE fourier_to_real_for_vtu
62 #include "petsc/finclude/petsc.h"
63  USE petsc
64  IMPLICIT NONE
65  TYPE(mesh_type) :: H_mesh, phi_mesh
66  REAL(KIND=8), INTENT(IN) :: dt
67  INTEGER, POINTER, DIMENSION(:),INTENT(IN) :: list_mode
68  REAL(KIND=8),DIMENSION(:), INTENT(IN):: mu_H_field
69  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: eigen
70  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: eigen_vect
71  REAL(KIND=8) :: rho, theta, rmax, norm, err
72  LOGICAL :: redo
73  INTEGER :: i, n, ndim, nconv, nmax, k, code, rank
74  CHARACTER(LEN=3) :: st_mode, st_eigen
75  REAL(KIND=8), ALLOCATABLE,DIMENSION(:,:,:):: vect, rot_Hn
76  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: ev_Hn, best_ev
77 !! 3D
78  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:):: ev_Hn_3d
79  CHARACTER(LEN=200) :: header
80  CHARACTER(LEN=3) :: what
81  INTEGER :: nb_procs_F, rank_F
82 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
83 !!$ REAL(KIND=8) :: time=0.
84 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
85 !!$!=============TEST FL, Apr. 11th, 2013
86 !!$ CHARACTER(LEN=200) :: name_best_ev
87 !!$ CHARACTER(LEN=3) :: tit_m
88 !!$!=============TEST FL, Apr. 11th, 2013
89  mpi_comm, DIMENSION(2) :: communicator
90  CALL mpi_comm_rank(mpi_comm_world,rank,code)
91  CALL mpi_comm_rank(communicator(2),rank_f, code)
92  CALL mpi_comm_size(communicator(2),nb_procs_f, code)
93 
94  redo = .false. !redo=.true. not yet programmed
95  ndim = 12*h_mesh%np
96  ALLOCATE(eigen(inputs%nb_eigenvalues,2),eigen_vect(ndim,inputs%nb_eigenvalues), &
97  ev_hn(h_mesh%np,6), rot_hn(h_mesh%np,6,1), ev_hn_3d(h_mesh%np,6,1),best_ev(h_mesh%np,6))
98  ALLOCATE(vect(6,h_mesh%np,1))
99 
100  ! ATTENTION mode_max_c must be equal to 1, meaning each processors is dealving with one Fourier mode.
101  DO i = 1, 1 !mode_max_c ! ATTENTION
102  ! ATTENTION
103  CALL arpack_not_sym(communicator, inputs%nb_eigenvalues, nconv, &
104  inputs%arpack_iter_max,inputs%arpack_tolerance, &
105  prodmat_maxwell_int_by_parts, eigen, eigen_vect, &
106  inputs%arpack_which, i, list_mode, redo)
107  !Postprocessing
108  WRITE(st_mode,'(I3)') list_mode(i)
109 
110  rmax = -1.d10
111  DO n = 1, min(nconv,inputs%nb_eigenvalues)
112  rho = sqrt(eigen(n,1)**2+eigen(n,2)**2)
113  theta = atan2(eigen(n,2),eigen(n,1))
114  IF (log(rho).GE.rmax) THEN
115  nmax = n
116  rmax = log(rho)
117  END IF
118  !===Visualization of fields. Create 2D vtu files.
119  IF (inputs%if_arpack_vtu_2d) THEN
120 !=============TEST FL, Feb. 11th, 2013
121  IF (n == 1) THEN
122  what='new'
123  ELSE
124  what='old'
125  END IF
126 !=============TEST FL, Feb. 11th, 2013
127  !===Magnetic field
128  DO k = 1, 6
129  ev_hn(:,k) = eigen_vect((k-1)*h_mesh%np+1:k*h_mesh%np,n)
130  ev_hn_3d(:,k,1) = eigen_vect((k-1)*h_mesh%np+1:k*h_mesh%np,n)
131  END DO
132  WRITE(st_eigen,'(I3)') n
133  header = 'Hn_Eigen_'//trim(adjustl(st_eigen))//'_mode_'//trim(adjustl(st_mode))
134  CALL make_vtu_file_2d(communicator(1), h_mesh, header, ev_hn, 'Hn', &
135  'new', opt_it=1)
136 !=============TEST FL, Feb. 11th, 2013
137  header='Hn_mode_'//trim(adjustl(st_mode))//'_theta_cos'
138  CALL make_vtu_file_arpack(communicator(1), h_mesh, header, ev_hn(:,3:3), 'Hn', &
139  what, n)
140 !! vtu_3D ne marche pas 19/03/2013, il faut ecrire par mode, 0 puis 1
141  !IF(n==nmax) THEN
142  ! CALL vtu_3d(ev_Hn_3d, 'H_mesh', 'MagField', 'mag', 'new', opt_it=1)
143  !ENDIF
144 !=============TEST FL, Feb. 11th, 2013
145 !=============TEST FL, Apr. 11th, 2013
146  IF (n==nmax) THEN
147  best_ev = ev_hn
148  ev_hn_3d(:,:,1) = ev_hn
149  END IF
150 !=============TEST FL, Apr. 11th, 2013
151 
152  !===Magnetic induction
153  DO k = 1, 6
154  ev_hn(:,k) = mu_h_field(:)*eigen_vect((k-1)*h_mesh%np+1:k*h_mesh%np,n)
155  END DO
156  header = 'Bn_Eigen_'//trim(adjustl(st_eigen))//'_mode_'//trim(adjustl(st_mode))
157  CALL make_vtu_file_2d(communicator(1), h_mesh, header, ev_hn, 'Bn', &
158  'new', opt_it=1)
159 !=============TEST FL, Feb. 11th, 2013
160  header='Bn_mode_'//trim(adjustl(st_mode))//'_theta_cos'
161  CALL make_vtu_file_arpack(communicator(1), h_mesh, header, ev_hn(:,3:3), 'Bn', &
162  what, n)
163 !=============TEST FL, Feb. 11th, 2013
164 !! vtu_3D ne marche pas 19/03/2013
165  !IF(n==nmax) THEN
166  ! CALL vtu_3d(ev_Hn_3d, 'H_mesh', 'IndField', 'ind', 'new', opt_it=1)
167  !ENDIF
168  END IF
169 
170 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
171  !===We could save things here
172  !CALL write_restart_maxwell(comm_one_d, H_mesh, phi_mesh, &
173  ! time, list_mode, Hn, Hn1, &
174  ! phin, phin1, inputs%file_name, it, freq_restart)
175 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
176 
177  END DO
178 
179  IF (nconv>0) THEN
180  rho = sqrt(eigen(nmax,1)**2+eigen(nmax,2)**2)
181  theta = atan2(eigen(nmax,2),eigen(nmax,1))
182  DO k = 1, 6
183  vect(k,:,1) = mu_h_field(:)*eigen_vect((k-1)*h_mesh%np+1:k*h_mesh%np,nmax)
184  END DO
185  norm = sqrt(norme_l2_champ_par(communicator(1), h_mesh, list_mode(i:i), vect)**2 &
186  + norme_h1_champ_par(communicator(1), h_mesh, list_mode(i:i), vect)**2)
187  err = norme_div_par(communicator(1), h_mesh, list_mode(i:i), vect)
188  WRITE(10+rank,*) log(rho)/dt, theta/dt, list_mode(i)
189  WRITE(10+rank,*) err/norm
190 
191  WRITE(*,'(A)') '==================================================='
192  WRITE(*,'(A,i2,A,i2,A,e12.5,A,e12.5,A)') '| mode',list_mode(i),', LR',nmax,&
193  ' = (',log(rho)/dt,' + I *',theta/dt,') |'
194  WRITE(*,'(A,e10.3,A)') '| |div(Bn)|_L2/|Bn|_H1 =', err/norm,' |'
195  WRITE(*,'(A)') '==================================================='
196  WRITE(10+rank,'(A)') '==================================================='
197  WRITE(10+rank,'(A,i2,A,i2,A,e12.5,A,e12.5,A)') '| mode',list_mode(i),', LR',nmax,&
198  ' = (',log(rho)/dt,' + I *',theta/dt,') |'
199  WRITE(10+rank,'(A,e10.3,A)') '| |div(Bn)|_L2/|Bn|_H1 =', err/norm,' |'
200  WRITE(10+rank,'(A)') '==================================================='
201  END IF
202 
203  DO n = 1, min(nconv,inputs%nb_eigenvalues)
204  rho = sqrt(eigen(n,1)**2+eigen(n,2)**2)
205  theta = atan2(eigen(n,2),eigen(n,1))
206  DO k = 1, 6
207  vect(k,:,1) = mu_h_field(:)*eigen_vect((k-1)*h_mesh%np+1:k*h_mesh%np,n)
208  END DO
209  norm = sqrt(norme_l2_champ_par(communicator(1), h_mesh, list_mode(i:i), vect)**2 &
210  + norme_h1_champ_par(communicator(1), h_mesh, list_mode(i:i), vect)**2)
211  err = norme_div_par(communicator(1), h_mesh, list_mode(i:i), vect)
212 
213  WRITE(*,'(A)') '==================================================='
214  WRITE(*,'(A,i2,A,i2,A,e12.5,A,e12.5,A)') '| mode',list_mode(i),', Eg',n,&
215  ' = (',log(rho)/dt,' + I *',theta/dt,') |'
216  WRITE(*,'(A,e10.3,A)') '| |div(Bn)|_L2/|Bn|_H1 =', err/norm,' |'
217  WRITE(*,'(A)') '==================================================='
218 
219  WRITE(10+rank,'(A)') '==================================================='
220  WRITE(10+rank,'(A,i2,A,i2,A,e12.5,A,e12.5,A)') '| mode',list_mode(i),', Eg',n,&
221  ' = (',log(rho)/dt,' + I *',theta/dt,') |'
222  WRITE(10+rank,'(A,e10.3,A)') '| |div(Bn)|_L2/|Bn|_H1 =', err/norm,' |'
223  WRITE(10+rank,'(A)') '==================================================='
224 
225  END do!===End postprocessing
226 
227 !!$!=============TEST FL, Apr. 11th, 2013
228 !!$ IF (inputs%if_arpack_vtu_2d) THEN
229 !!$ DO n = 1, nb_procs_F
230 !!$ IF (n==rank_F+1) THEN
231 !!$ ev_Hn_3d(:,:,1) = best_ev
232 !!$ ELSE
233 !!$ ev_Hn_3d = 0.d0
234 !!$ END IF
235 !!$ WRITE(tit_m,'(i3)') n-1
236 !!$ name_best_ev = 'Eigenvect_mode_'//trim(adjustl(tit_m))
237 !!$ CALL vtu_3d(ev_Hn_3d, 'H_mesh', name_best_ev, 'mag', 'new')
238 !!$ END DO
239 !!$ END IF
240 !!$!=============TEST FL, Apr. 11th, 2013
241 
242 
243  END DO !===End loop on Fourier modes
244 
245  CALL post_proc_arpack(communicator, h_mesh, phi_mesh, list_mode, mu_h_field, ev_hn_3d)
246 
247  CLOSE(10+rank)
248  DEALLOCATE(eigen,eigen_vect,vect)
249  END SUBROUTINE arpack_maxwell_int_by_parts
250  !-----------------------------------------------------------------
251 
252  !-----------------------------------------------------------------
253  SUBROUTINE arpack_navier_stokes
254  RETURN
255  END SUBROUTINE arpack_navier_stokes
256  !-----------------------------------------------------------------
257 
258  !-----------------------------------------------------------------
259  SUBROUTINE arpack_not_sym(communicator, nb_vp, nconv, iter_max, tol_arpack, prodmat, &
260  eigen, eigen_vect, which, imode, list_mode, redo)
261  USE my_util
262 #include "petsc/finclude/petsc.h"
263  USE petsc
264  IMPLICIT NONE
265  INTERFACE
266  SUBROUTINE prodmat(x,b,m,i)
267  INTEGER :: m, i
268  REAL(KIND=8), DIMENSION(m) :: x
269  REAL(KIND=8), DIMENSION(m) :: b
270  END SUBROUTINE prodmat
271  END INTERFACE
272  INTEGER, INTENT(IN) :: nb_vp, iter_max, imode
273  REAL(KIND=8), INTENT(INOUT) :: tol_arpack
274  REAL(KIND=8), DIMENSION(nb_vp,2), INTENT(OUT) :: eigen
275  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: eigen_vect
276  INTEGER, POINTER, DIMENSION(:), INTENT(IN) :: list_mode
277  LOGICAL, INTENT(IN) :: redo
278 
279  !ARPACK --------------------------------------------------------------------
280  INTEGER :: IDO, INFO, it=0
281  CHARACTER(len=1) :: BMAT
282  CHARACTER(len=2) :: WHICH
283  INTEGER, DIMENSION(11) :: IPARAM, IPNTR
284  INTEGER :: NDIM, NEV, NCV, LWORKL, NCONV
285  REAL(KIND=8) :: SIGMAR, sigmai
286 !!$ REAL(KIND=8), POINTER, DIMENSION(:) :: RESID, WORKL
287 !!$ REAL(KIND=8), POINTER, DIMENSION(:) :: WORKD
288 !!$ REAL(KIND=8), POINTER, DIMENSION(:) :: WORKEV
289 !!$ REAL(KIND=8), POINTER, DIMENSION(:,:) :: V, d
290 !+++++++++++++++++++TEST FL
291  REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: RESID, WORKL
292  REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: WORKD
293  REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: WORKEV
294  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: V, d
295  LOGICAL, ALLOCATABLE, DIMENSION(:) :: SEL
296 !+++++++++++++++++++TEST FL
297  INTEGER :: NP, NEV2, i, j, unit_save, &
298  ndim_r, nb_vp_r, code, rank, nb_procs, ido_loop
299  REAL(KIND=8) :: rnorm
300  INTEGER, ALLOCATABLE, DIMENSION(:) :: tableau_ido
301  LOGICAL :: bidon
302  !-------------END OF DECLARATIONS-------------------------------------------
303  mpi_comm, DIMENSION(2) :: communicator
304 
305  CALL mpi_comm_rank(mpi_comm_world,rank,code)
306 
307  CALL mpi_comm_size(communicator(2),nb_procs,code)
308  ALLOCATE(tableau_ido(nb_procs))
309 
310  nev = nb_vp
311  ndim = SIZE(eigen_vect,1)
312  ncv = 2*nev+1
313  !NCV = NEV + 2 ! Does not work with P2P2
314  lworkl = 3*ncv**2+6*ncv
315  ALLOCATE(v(ndim,ncv), workl(lworkl), resid(ndim))
316  ALLOCATE(sel(ncv), d(ncv,3), workev(3*ncv))
317  ALLOCATE(workd(3*ndim))
318 
319  IF (redo) THEN
320  unit_save=200+list_mode(imode)
321  rewind(unit_save)
322  READ(unit_save) ido, tol_arpack, iparam, ipntr, info, &
323  np, rnorm, nconv, nev2, ndim_r, bmat, nb_vp_r
324  IF (ndim_r/=ndim .OR. nb_vp_r/=nb_vp) THEN
325  CALL error_petsc('BUG in restart of arpack_mhd')
326  END IF
327  READ(unit_save) (resid(i), i = 1, ndim)
328  READ(unit_save) (workd(i), i = 1, 3*ndim)
329  READ(unit_save) (workl(i), i = 1, lworkl)
330  DO j=1,ncv
331  READ(unit_save) (v(i,j), i = 1, ndim)
332  ENDDO
333 
334  ELSE
335  ido=0
336  bmat = 'I'
337  info=0
338  iparam(1)=1
339  iparam(2)=1
340  iparam(4)=1
341  iparam(5)=1
342  iparam(6)=1
343  iparam(7)=1
344  iparam(8)=1
345  iparam(9)=1
346  iparam(10)=1
347  iparam(11)=1
348 
349  !----INITIALIZE SO THAT BC ARE ENFORCED CORRECTLY---------------------------
350  workd(1:ndim) = 1.d0
351  CALL prodmat(workd(1:ndim), resid, ndim, imode)
352  workd(ndim/2+1:ndim)=workd(1:ndim/2)
353  info = 1
354  !---------------------------------------------------------------------------
355  END IF
356 
357  iparam(3)=iter_max
358  bidon = .false.
359  ido_loop = -99
360  !---------------------------------------------------------------------------
361  DO WHILE(ido_loop.NE.99)
362  it = it+1
363  IF(rank==0 .AND. mod(it,100)==0) WRITE(*,'(a,I5)') ' Arpack iteration', it
364  IF (.NOT.bidon) THEN
365  CALL pdnaupd(communicator(1), ido,bmat,ndim,which,nev, tol_arpack, &
366  resid, ncv, v, ndim, iparam, ipntr,&
367  workd, workl, lworkl, info)
368  END IF
369 
370  !===We make sure that all the processors in communicator(2) call prodmat to avoid hanging
371  CALL mpi_allgather(ido, 1, mpi_integer, tableau_ido, 1, mpi_integer, &
372  communicator(2), code)
373 
374  IF (minval(abs(tableau_ido+2))==0) THEN
375  ido = -2
376  END IF
377 
378  IF (maxval(abs(tableau_ido-99))/=0) THEN
379  ido_loop = -99
380  ELSE
381  ido_loop = 99
382  END IF
383 
384  IF (ido==99) THEN
385  IF (.NOT. bidon) THEN
386  WRITE(*,'(A,i3,A,i5)') 'ARPACK OK for mode', list_mode(imode), &
387  ' Nb of Arpack iterations', iparam(3)
388  END IF
389  bidon=.true.
390 
391  CALL prodmat(workd(ipntr(1):), workd(ipntr(2):), ndim, imode)
392  ELSE IF(ido==-1 .OR. ido==1) THEN
393  CALL prodmat(workd(ipntr(1):), workd(ipntr(2):), ndim, imode)
394  ENDIF
395 
396 
397 !!$ IF (ido == -2) THEN
398 !!$ !
399 !!$ ! %---------------------------------------------------%
400 !!$ ! | After maxitr iterations without convergence, |
401 !!$ ! | output the computed quantities to the file state. |
402 !!$ ! %---------------------------------------------------%
403 !!$ !
404 !!$ PRINT *,'***** DEBUT SAUVEGARDE # ', list_mode(imode), imode
405 !!$ unit_save=200+list_mode(imode)
406 !!$ REWIND(unit_save)
407 !!$ WRITE(unit_save) ido, tol_arpack, iparam, ipntr, info, &
408 !!$ np, rnorm, nconv, nev2, ndim, bmat, nb_vp
409 !!$ WRITE(unit_save) (resid(i), i = 1, ndim)
410 !!$ WRITE(unit_save) (workd(i), i = 1, 3*ndim)
411 !!$ WRITE(unit_save) (workl(i), i = 1, lworkl)
412 !!$ DO j=1,ncv
413 !!$ WRITE(unit_save) (v(i,j), i = 1, ndim)
414 !!$ ENDDO
415 !!$ nconv = -1
416 !!$ PRINT *,'***** FIN SAUVEGARDE # ', list_mode(imode), imode
417 !!$ DEALLOCATE(tableau_ido)
418 !!$ RETURN
419 !!$ ENDIF
420 
421  END DO
422 
423  nconv = iparam(5)
424  IF (nconv/=nev) THEN
425  WRITE(*,*) ' Only ', nconv, 'eigenvalues have been computed instead of', nev
426  WRITE(*,*) ' it_max not large enough'
427  END IF
428 
429  IF(nconv.NE.0) THEN
430  CALL pdneupd (communicator(1), .true., 'All', sel, d, d(1,2), v, ndim, &
431  sigmar, sigmai, workev, bmat, ndim, which, nev, tol_arpack, &
432  resid, ncv, v, ndim, iparam, ipntr, workd, workl, lworkl, info)
433 
434  nev2=min(nb_vp, nconv)
435  eigen_vect(:,1:nev2) = v(:,1:nev2)
436  eigen(1:nev2,1) = d(1:nev2,1)
437  eigen(1:nev2,2) = d(1:nev2,2)
438  END IF
439 
440 !!$ IF (ALLOCATED(tableau_ido)) DEALLOCATE(tableau_ido)
441 !!$ IF (ASSOCIATED(RESID)) NULLIFY(RESID)
442 !!$ IF (ASSOCIATED(SEL)) NULLIFY(SEL)
443 !!$ IF (ASSOCIATED(d)) NULLIFY(d)
444 !!$ IF (ASSOCIATED(workev)) NULLIFY(workev)
445 !!$ IF (ASSOCIATED(WORKD)) NULLIFY(WORKD)
446 !!$ IF (ASSOCIATED(WORKL)) NULLIFY(WORKL)
447 !!$ IF (ASSOCIATED(V)) NULLIFY(V)
448 !+++++++++++++++++++TEST FL
449 
450 
451  END SUBROUTINE arpack_not_sym
452  !-----------------------------------------------------------------
453 
454  SUBROUTINE post_proc_arpack(communicator, H_mesh, phi_mesh, list_mode, mu_H_field, eigen_vect)
455  USE def_type_mesh
456  USE chaine_caractere
457  USE vtk_viz
458  USE tn_axi
459  USE input_data
460  USE fourier_to_real_for_vtu
461  USE restart
462 #include "petsc/finclude/petsc.h"
463  USE petsc
464  IMPLICIT NONE
465  TYPE(mesh_type), INTENT(IN) :: H_mesh, phi_mesh
466  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
467  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
468  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: eigen_vect
469 
470  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: field_H, field_phi
471  INTEGER :: nb_procs_F, rank_F, code, n, k
472  CHARACTER(len=3) :: tit_m
473  CHARACTER(len=200) :: name_best_ev
474  mpi_comm, DIMENSION(2) :: communicator
475 
476  CALL mpi_comm_rank(communicator(2),rank_f, code)
477  CALL mpi_comm_size(communicator(2),nb_procs_f, code)
478  ALLOCATE(field_h(h_mesh%np, 6, 1), field_phi(phi_mesh%np, 2, 1))
479  field_phi = 0.d0
480  IF (inputs%if_arpack_vtu_2d) THEN
481  DO n = 1, nb_procs_f
482  IF (n==rank_f+1) THEN
483  field_h(:,:,1) = eigen_vect(:,:,1)
484  ELSE
485  field_h(:,:,1) = 0.d0
486  END IF
487  WRITE(tit_m,'(i3)') n-1
488  name_best_ev = 'Eigenvect_H_mode_'//trim(adjustl(tit_m))
489  CALL vtu_3d(field_h, 'H_mesh', name_best_ev, 'mag', 'new')
490  DO k = 1, 6
491  field_h(:,k,1) = mu_h_field*field_h(:,k,1)
492  END DO
493  name_best_ev = 'Eigenvect_B_mode_'//trim(adjustl(tit_m))
494  CALL vtu_3d(field_h, 'H_mesh', name_best_ev, 'ind', 'new')
495  END DO
496  END IF
497 
498  DO k = 1, 6
499  field_h(:,k,1) = mu_h_field*eigen_vect(:,k,1)
500  END DO
501  CALL write_restart_maxwell(communicator, h_mesh, phi_mesh, 1.d0, list_mode, eigen_vect, eigen_vect, &
502  field_h, field_h, field_phi, field_phi, &
503  inputs%file_name, 1, 1)
504 
505  END SUBROUTINE post_proc_arpack
506 
507 
508 END MODULE arpack_mhd
arpack_mhd::arpack_not_sym
subroutine arpack_not_sym(communicator, nb_vp, nconv, iter_max, tol_arpack, prodmat, eigen, eigen_vect, which, imode, list_mode, redo)
Definition: template_arpack_mhd.f90:261
initialization::prodmat_maxwell_int_by_parts
subroutine, public prodmat_maxwell_int_by_parts(vect_in, vect_out, ndim, i)
Definition: initialization.f90:1502
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
fourier_to_real_for_vtu::make_vtu_file_2d
subroutine, public make_vtu_file_2d(communicator, mesh_in, header, field_in, field_name, what, opt_it)
Definition: fourier_real.f90:1487
my_util
Definition: my_util.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
vtk_viz::make_vtu_file_arpack
subroutine, public make_vtu_file_arpack(communicator, mesh, header, field, field_name, what, num_vp)
Definition: plot_vtk.f90:394
arpack_mhd::arpack_maxwell_int_by_parts
subroutine arpack_maxwell_int_by_parts(communicator, H_mesh, phi_mesh, dt, list_mode, mu_H_field)
Definition: template_arpack_mhd.f90:55
vtk_viz
Definition: plot_vtk.f90:1
fourier_to_real_for_vtu::vtu_3d
subroutine, public vtu_3d(field, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl, opt_2D, opt_mesh_in)
Definition: fourier_real.f90:272
tn_axi::norme_h1_champ_par
real(kind=8) function norme_h1_champ_par(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:206
tn_axi::norme_div_par
real(kind=8) function norme_div_par(communicator, H_mesh, list_mode, v)
Definition: tn_axi.f90:223
arpack_mhd::arpack_navier_stokes
subroutine arpack_navier_stokes
Definition: template_arpack_mhd.f90:254
tn_axi::norme_l2_champ_par
real(kind=8) function norme_l2_champ_par(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:185
tn_axi
Definition: tn_axi.f90:5
fourier_to_real_for_vtu
Definition: fourier_real.f90:1
restart::write_restart_maxwell
subroutine write_restart_maxwell(communicator, H_mesh, phi_mesh, time, list_mode, Hn, Hn1, Bn, Bn1, phin, phin1, filename, it, freq_restart, opt_mono, opt_dt)
Definition: restart.f90:677
arpack_mhd::post_proc_arpack
subroutine, public post_proc_arpack(communicator, H_mesh, phi_mesh, list_mode, mu_H_field, eigen_vect)
Definition: template_arpack_mhd.f90:455
input_data
Definition: read_my_data.f90:233
chaine_caractere
Definition: chaine_caractere.f90:4
def_type_mesh
Definition: def_type_mesh.f90:4
arpack_mhd
Definition: template_arpack_mhd.f90:1
arpack_mhd::solver_arpack_mhd
subroutine, public solver_arpack_mhd(communicator, H_mesh, phi_mesh, dt, list_mode, mu_H_field)
Definition: template_arpack_mhd.f90:10