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