SFEMaNS  version 5.3
Reference documentation for SFEMaNS
mesh_interpolation.f90
Go to the documentation of this file.
2 
3  PUBLIC :: mesh_interpol
4  PRIVATE
5 CONTAINS
6  SUBROUTINE mesh_interpol
7 
9  USE def_type_mesh
10  USE my_util
11  USE create_comm
13  USE metis_sfemans
14  USE prep_maill
15  USE sub_plot
16  USE restart
17 #include "petsc/finclude/petsc.h"
18  USE petsc
19  IMPLICIT NONE
20 
21  TYPE(mesh_type) :: p1_mesh_glob, p2_mesh_glob, p2_c0_mesh_glob_temp
22 
23  TYPE(mesh_type), TARGET :: H_mesh, phi_mesh
24  TYPE(mesh_type), TARGET :: H_mesh_glob, phi_mesh_glob
25  TYPE(mesh_type), POINTER :: H_mesh_in, H_mesh_out, phi_mesh_in, phi_mesh_out
26 
27  TYPE(mesh_type), TARGET :: vv_mesh, pp_mesh
28  TYPE(mesh_type), TARGET :: vv_mesh_glob, pp_mesh_glob
29  TYPE(mesh_type), POINTER :: vv_mesh_in, vv_mesh_out, pp_mesh_in, pp_mesh_out
30 
31  TYPE(mesh_type), TARGET :: temp_mesh
32  TYPE(mesh_type), TARGET :: temp_mesh_glob
33  TYPE(mesh_type), POINTER :: temp_mesh_in, temp_mesh_out
34 
35  CHARACTER(len=200) :: old_filename, old_directory
36  CHARACTER(len=200) :: new_filename, new_directory
37  CHARACTER(len=200) :: file_name, directory
38  CHARACTER(len=200) :: file_name_m, directory_m
39  LOGICAL :: if_momentum, if_mass, if_induction, if_energy, &
40  old_is_form, new_is_form, iformatted, is_form_m, mono_in, mono_out, is_in, &
41  dom_H_larger_dom_ns, dom_temp_larger_dom_ns, &
42  if_temperature_in, if_temperature_out, if_temperature, if_level_set, inter_mesh, &
43  rw_ns, rw_temp, rw_mxw, check_plt, test, if_select
44 
45  INTEGER :: nb_S, nb_F, nb_procs, rank, type_fe_H, type_fe_phi, &
46  nb_dom_ns, nb_dom_temp, nb_dom_H, nb_dom_phi, nsize, code, m, &
47  m_max_c, rank_S, nb_inter, nb_inter_mu, nb_inter_v_T, &
48  k, kp, i, nb_mode, rang_ns_S, rang_temp_S, rang_S, l, lblank, nb_fluid
49  INTEGER :: nb_fic, index_start
50  REAL(KIND=8) :: time_h, time_u, time_T, max_vel
51  TYPE(periodic_data) :: my_periodic
52 
53  INTEGER, DIMENSION(:), ALLOCATABLE :: list_inter, list_inter_temp, list_inter_v_T, &
54  list_inter_mu, list_inter_H_phi, list_dom_H, &
55  list_dom_ns, list_dom_temp, list_dom_temp_in, list_dom_phi, list_dom_H_in, list_dom, part, &
56  list_mode, controle_H, controle_phi, &
57  controle_vv, controle_pp, controle_temp, temp_in_to_new, H_in_to_new, &
58  l_t_g_vv, l_t_g_pp, l_t_g_H, l_t_g_phi, l_t_g_temp
59 
60  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: Hn, Hn1, Bn, Bn1, phin, phin1
61  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: Hn_glob, Hn1_glob, Bn_glob, Bn1_glob, phin_glob, phin1_glob
62  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: un, un_m1, pn, pn_m1
63  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: incpn, incpn_m1, tempn, tempn_m1
64  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: level_setn, level_setn_m1
65  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: un_glob, un_m1_glob, pn_glob, pn_m1_glob
66  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: incpn_glob, incpn_m1_glob, tempn_glob, tempn_m1_glob
67  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: level_setn_glob, level_setn_m1_glob
68 
69  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: Hn_in, Hn1_in, Bn_in, Bn1_in, phin_in, phin1_in
70  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: Hn_out, Hn1_out, Bn_out, Bn1_out, phin_out, phin1_out
71  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: un_in, un_m1_in, pn_in, pn_m1_in
72  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: incpn_in, incpn_m1_in, tempn_in, tempn_m1_in
73  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: level_setn_in, level_setn_m1_in
74  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: un_out, un_m1_out, pn_out, pn_m1_out
75  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: incpn_out, incpn_m1_out, tempn_out, tempn_m1_out
76  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: level_setn_out, level_setn_m1_out
77 
78  CHARACTER(len=3) :: type_pb, tit_m, tit_s, tit
79  CHARACTER(len=4) :: tit_part
80  CHARACTER(len=200) :: mesh_part_name
81  INTEGER :: petsc_rank
82  LOGICAL :: if_read_partition
83  LOGICAL :: if_level_set_P2
84  mpi_comm :: comm_cart
85  mpi_comm, DIMENSION(:), POINTER :: comm_one_d, comm_one_d_ns, &
86  comm_one_d_temp, coord_cart
87  petscerrorcode :: ierr
88 
89  CALL petscinitialize(petsc_null_character,ierr)
90  CALL mpi_comm_rank(petsc_comm_world, petsc_rank, code)
91 
92  OPEN(unit=22, file='data_interpol', status='unknown', access='sequential')
93 
94  !===Compute is_in (is_in=.TRUE. if Petsc -> NonPetsc; is_in=.FALSE. if NonPetsc -> Petsc)
95  CALL which_pb(is_in)
96 
97  IF (is_in) THEN
98  CALL read_until(22, '===Number of processors in meridian section (Input)')
99  READ(22,*) nb_s
100  CALL read_until(22, '===Problem type (Input): (nst, mxw, mhd, fhd)')
101  READ(22,*) type_pb
102  CALL find_string(22, '===Is there an input temperature field?', test)
103  IF (test) THEN
104  READ (22, *) if_temperature_in
105  ELSE
106  if_temperature_in = .false.
107  END IF
108  if_temperature = if_temperature_in
109  inter_mesh = .false.
110  if_read_partition = .true.
111  ELSE
112  CALL read_until(22, '===Number of processors in meridian section (Output)')
113  READ(22,*) nb_s
114  CALL read_until(22, '===Problem type (Output): (nst, mxw, mhd, fhd)')
115  READ(22,*) type_pb
116  CALL find_string(22, '===Is there an output temperature field?', test)
117  IF (test) THEN
118  READ (22, *) if_temperature_out
119  ELSE
120  if_temperature_out = .false.
121  END IF
122  if_temperature = if_temperature_out
123  CALL read_until(22, '===Should data be interpolated on new mesh? (True/False)')
124  READ(22,*) inter_mesh
125  if_read_partition = .false.
126  END IF
127 
128  !===Set rw_ns, rw_mxw and rw_temp (function of Problem type (Input))
129  CALL read_until(22, '===Problem type (Input): (nst, mxw, mhd, fhd)')
130  READ(22,*) tit
131  IF (tit=='nst') THEN
132  rw_ns = .true.
133  rw_mxw = .false.
134  ELSE IF ((tit=='mhd').OR.(tit=='fhd')) THEN
135  rw_ns = .true.
136  rw_mxw = .true.
137  ELSE IF (tit=='mxw') THEN
138  rw_ns = .false.
139  rw_mxw = .true.
140  ELSE
141  CALL error_petsc('BUG in interpol, type_pb (Input) not correct')
142  END IF
143  CALL find_string(22, '===Is there an input temperature field?', test)
144  IF (test) THEN
145  READ (22, *) rw_temp
146  ELSE
147  rw_temp = .false.
148  END IF
149 
150  !===Default parameters
151  nb_f = 1
152 
153  !===Check construction with plotmtv
154  CALL find_string(22, '===Check construction with plotmtv? (True/False)', test)
155  IF (test) THEN
156  READ (22, *) check_plt
157  ELSE
158  check_plt = .false.
159  END IF
160 
161  !===Number of time steps to process
162  CALL read_until(22, '===How many files (times steps) should be converted?')
163  READ(22, *) nb_fic
164 
165  !===Starting index
166  CALL find_string(22, '===What is the starting index? (suite*_I*index.mesh*), index is an integer in [1,999]', test)
167  IF (test) THEN
168  READ (22, *) index_start
169  ELSE
170  index_start = 1
171  END IF
172 
173  !===Number of Fourier modes
174  CALL read_until(22, '===Number of Fourier modes')
175  READ(22, *) nb_mode
176  ALLOCATE(list_mode(nb_mode))
177  CALL find_string(22, '===Select Fourier modes? (true/false)',test)
178  IF (test) THEN
179  READ(22,*) if_select
180  IF (if_select) THEN
181  CALL read_until(22, '===List of Fourier modes (if select_mode=.TRUE.)')
182  READ(22,*) list_mode
183  ELSE
184  list_mode = (/ (m, m=0, nb_mode-1) /)
185  END IF
186  ELSE
187  list_mode = (/ (m, m=0, nb_mode-1) /)
188  END IF
189 
190  !===Input and Output mesh files
191  CALL read_until(22, '===Directory and name of input mesh file')
192  READ(22,*) old_directory, old_filename
193  CALL read_until(22, '===Is input mesh file formatted (true/false)?')
194  READ(22,*) old_is_form
195  IF (is_in) THEN !Interpolation is done at second step only
196  new_directory = old_directory
197  new_filename = old_filename
198  new_is_form = old_is_form
199  ELSE
200  CALL read_until(22, '===Directory and name of output mesh file')
201  READ(22,*) new_directory, new_filename
202  CALL read_until(22, '===Is output mesh file formatted (true/false)?')
203  READ(22,*) new_is_form
204  END IF
205 
206  !===Data for NS
207  IF ( (type_pb == 'nst') .OR. (type_pb == 'mhd') .OR. (type_pb == 'fhd')) THEN
208  CALL read_until(22, '===Number of subdomains in Navier-Stokes mesh')
209  READ(22,*) nb_dom_ns
210  ALLOCATE(list_dom_ns(nb_dom_ns))
211  CALL read_until(22, '===List of subdomains for Navier-Stokes mesh')
212  READ (22, *) list_dom_ns
213  CALL find_string(22, '===Is there a level set?', test)
214  IF (test) THEN
215  READ (22, *) if_level_set
216  IF (if_level_set) THEN
217  CALL read_until(22, '===How many fluids?')
218  READ(22,*) nb_fluid
219 
220  CALL find_string(22, '===Do we use P2 finite element for level_set? (true/false)', test)
221  IF (test) THEN
222  READ (22, *) if_level_set_p2
223  ELSE
224  if_level_set_p2=.true.
225  END IF
226  END IF
227  ELSE
228  if_level_set = .false.
229  END IF
230  IF (type_pb == 'nst') ALLOCATE(list_dom_h(0), list_dom_phi(0))
231  END IF
232 
233  !===Data for temperature
234  IF (if_temperature) THEN
235  CALL read_until(22, '===Number of subdomains in temperature mesh')
236  READ(22,*) nb_dom_temp
237  ALLOCATE(list_dom_temp_in(nb_dom_temp), list_dom_temp(nb_dom_temp), temp_in_to_new(nb_dom_temp))
238  CALL read_until(22, '===List of subdomains for temperature mesh')
239  READ (22, *) list_dom_temp_in
240  CALL find_string(22, '===Number of interfaces between velocity and temperature only domains (for nst applications)', test)
241  IF (test) THEN
242  READ(22,*) nb_inter_v_t
243  ALLOCATE(list_inter_v_t(nb_inter_v_t))
244  CALL read_until(22, '===List of interfaces between velocity and temperature only domains (for nst applications)')
245  READ(22,*) list_inter_v_t
246  ELSE
247  ALLOCATE(list_inter_v_t(0))
248  END IF
249  ELSE
250  ALLOCATE(list_dom_temp(0))
251  END IF
252 
253  !===Data for Maxwell
254  IF ( (type_pb == 'mhd') .OR. (type_pb == 'mxw') .OR. (type_pb == 'fhd')) THEN
255  CALL read_until(22, '===Type of finite element for magnetic field')
256  READ (22, *) type_fe_h
257  CALL read_until(22, '===Number of subdomains in magnetic field (H) mesh')
258  READ(22,*) nb_dom_h ! number of sub_domains for H
259  ALLOCATE(list_dom_h_in(nb_dom_h), list_dom_h(nb_dom_h), h_in_to_new(nb_dom_h)) ! JLG/AR Nov 17 2008
260  CALL read_until(22, '===List of subdomains for magnetic field (H) mesh')
261  READ(22,*) list_dom_h_in
262 
263  !==========Interfaces H/H==========================!
264  CALL read_until(22, '===Number of interfaces in H mesh')
265  READ(22,*) nb_inter_mu
266  ALLOCATE(list_inter_mu(nb_inter_mu))
267  IF (nb_inter_mu>0) THEN
268  CALL read_until(22, '===List of interfaces in H mesh')
269  READ(22,*) list_inter_mu
270  END IF
271 
272  !==========Subdomains for phi======================!
273  CALL find_string(22, '===Number of subdomains in magnetic potential (phi) mesh', test)
274  IF (test) THEN
275  READ (22, *) nb_dom_phi
276  ELSE
277  nb_dom_phi = 0
278  END IF
279  ALLOCATE(list_dom_phi(nb_dom_phi))
280  IF (nb_dom_phi>0) THEN
281  !==========List of subdomains for phi======================!
282  CALL read_until(22, '===List of subdomains for magnetic potential (phi) mesh')
283  READ(22,*) list_dom_phi
284 
285  !==========H/phi interface=========================!
286  CALL read_until(22, '===Number of interfaces between H and phi')
287  READ(22,*) nb_inter ! number of interfaces between H and phi
288  ALLOCATE(list_inter_h_phi(nb_inter))
289  IF (nb_inter>0) THEN
290  CALL read_until(22, '===List of interfaces between H and phi')
291  READ(22,*) list_inter_h_phi
292  END IF
293 
294  !==========Finite element type=====================!
295  CALL read_until(22, '===Type of finite element for scalar potential')
296  READ(22,*) type_fe_phi
297  ELSE
298  ALLOCATE(list_inter_h_phi(0))
299  type_fe_phi = -1
300  END IF
301  IF (type_pb == 'mxw') THEN
302  ALLOCATE(list_dom_ns(0))
303  END IF
304  END IF
305 
306  CALL find_string(22, '===How many pieces of periodic boundary?', test)
307  IF (test) THEN
308  READ (22, *) my_periodic%nb_periodic_pairs
309  ELSE
310  my_periodic%nb_periodic_pairs = 0
311  END IF
312 
313  IF (my_periodic%nb_periodic_pairs.GE.1) THEN
314  ALLOCATE(my_periodic%list_periodic(2,my_periodic%nb_periodic_pairs))
315  ALLOCATE(my_periodic%vect_e(2,my_periodic%nb_periodic_pairs))
316  CALL read_until(22, '===Indices of periodic boundaries and corresponding vectors')
317  DO k = 1, my_periodic%nb_periodic_pairs
318  READ(22,*) my_periodic%list_periodic(:,k), my_periodic%vect_e(:,k)
319  END DO
320  END IF
321  CLOSE(22)
322 
323  !===Creation of logicals for equations
324  if_mass = if_level_set
325  if_momentum = type_pb=='nst' .OR. type_pb=='mhd' .OR. type_pb=='fhd'
326  if_induction = type_pb=='mxw' .OR. type_pb=='mhd' .OR. type_pb=='fhd' ! No mxx: for mxx, put mhd
327  if_energy = if_temperature
328 
329  !===Check that vv_mesh is a subset of H_mesh
330  IF (if_induction) THEN
331  IF (if_momentum) THEN
332  IF (SIZE(list_dom_h) < SIZE(list_dom_ns)) THEN
333  WRITE(*,*) ' BUG: NS must be a subset of Maxwell '
334  stop
335  END IF
336  DO k = 1, nb_dom_ns
337  ! JLG/AR Nov 17 2008
338  IF (minval(abs(list_dom_h_in - list_dom_ns(k))) /= 0) THEN
339  WRITE(*,*) ' BUG : NS must be a subset of Maxwell '
340  stop
341  END IF
342  DO kp = 1, nb_dom_h
343  IF (list_dom_h_in(kp) == list_dom_ns(k)) EXIT
344  END DO
345  h_in_to_new(k) = kp
346  ! JLG/AR Nov 17 2008
347  list_dom_h(k) = list_dom_ns(k)
348  END DO
349  m = nb_dom_ns
350  DO k = 1, nb_dom_h
351  IF (minval(abs(list_dom_h_in(k) - list_dom_ns)) == 0) cycle
352  m = m + 1
353  ! JLG/AR Nov 17 2008
354  h_in_to_new(m) = k
355  ! JLG/AR Nov 17 2008
356  list_dom_h(m) = list_dom_h_in(k)
357  END DO
358  IF (m/=nb_dom_h) THEN
359  WRITE(*,*) ' BUG : m/=nb_dom_H '
360  stop
361  END IF
362  IF (nb_dom_h > nb_dom_ns) THEN
363  dom_h_larger_dom_ns = .true.
364  ELSE
365  dom_h_larger_dom_ns = .false.
366  END IF
367  ELSE
368  ! JLG/AR Nov 17 2008
369  DO k = 1, nb_dom_h
370  h_in_to_new(k) = k
371  END DO
372  ! JLG/AR Nov 17 2008
373  list_dom_h = list_dom_h_in
374  END IF
375  END IF
376 
377  !===Check that vv_mesh is a subset of temp_mesh
378  IF (if_energy) THEN
379  IF (SIZE(list_dom_temp) < SIZE(list_dom_ns)) THEN
380  WRITE(*,*) ' BUG: NS must be a subset of temp '
381  stop
382  END IF
383  DO k = 1, nb_dom_ns
384  IF (minval(abs(list_dom_temp_in - list_dom_ns(k))) /= 0) THEN
385  WRITE(*,*) ' BUG : NS must be a subset of temp '
386  stop
387  END IF
388  DO kp = 1, nb_dom_temp
389  IF (list_dom_temp_in(kp) == list_dom_ns(k)) EXIT
390  END DO
391  temp_in_to_new(k) = kp
392  list_dom_temp(k) = list_dom_ns(k)
393  END DO
394  m = nb_dom_ns
395  DO k = 1, nb_dom_temp
396  IF (minval(abs(list_dom_temp_in(k) - list_dom_ns)) == 0) cycle
397  m = m + 1
398  temp_in_to_new(m) = k
399  list_dom_temp(m) = list_dom_temp_in(k)
400  END DO
401  IF (m/=nb_dom_temp) THEN
402  WRITE(*,*) ' BUG : m/=nb_dom_temp '
403  stop
404  END IF
405  IF (nb_dom_temp > nb_dom_ns) THEN
406  dom_temp_larger_dom_ns = .true.
407  ELSE
408  dom_temp_larger_dom_ns = .false.
409  END IF
410  END IF
411 
412  !===Check that temp_mesh is a subset of H_mesh
413  IF (if_induction .AND. if_energy) THEN
414  IF (SIZE(list_dom_h) < SIZE(list_dom_temp)) THEN
415  WRITE(*,*) ' BUG: temp must be a subset of Maxwell '
416  stop
417  END IF
418  DO k = 1, nb_dom_temp
419  IF (minval(abs(list_dom_h - list_dom_temp(k))) /= 0) THEN
420  WRITE(*,*) ' BUG: temp must be a subset of Maxwell '
421  stop
422  END IF
423  END DO
424  END IF
425 
426  !===Create interfaces in meshes
427  IF (if_momentum .AND. (.NOT. if_induction)) THEN
428  IF (if_energy) THEN
429  nsize = SIZE(list_dom_temp)
430  ALLOCATE(list_dom(nsize))
431  list_dom = list_dom_temp
432  ALLOCATE(list_inter(SIZE(list_inter_v_t)))
433  list_inter = list_inter_v_t
434  ELSE
435  nsize = SIZE(list_dom_ns)
436  ALLOCATE(list_dom(nsize))
437  list_dom = list_dom_ns
438  ALLOCATE(list_inter(0))
439  END IF
440  ELSE
441  nsize = SIZE(list_dom_h)+SIZE(list_dom_phi)
442  ALLOCATE(list_dom(nsize))
443  list_dom(1:SIZE(list_dom_h)) = list_dom_h
444  IF (SIZE(list_dom_phi)>0) THEN
445  list_dom(SIZE(list_dom_h)+1:) = list_dom_phi
446  END IF
447  nsize = SIZE(list_inter_mu)+SIZE(list_inter_h_phi)
448  ALLOCATE(list_inter(nsize))
449  IF (SIZE(list_inter_mu)>0) THEN
450  list_inter(1:SIZE(list_inter_mu)) = list_inter_mu
451  END IF
452  IF (SIZE(list_inter_h_phi)>0) THEN
453  list_inter(SIZE(list_inter_mu)+1:) = list_inter_h_phi
454  END IF
455  END IF
456  IF (if_energy) THEN
457  ALLOCATE(list_inter_temp(0))
458  END IF
459 
460  !===Directory, file name and format
461  IF (is_in) THEN
462  directory = old_directory
463  file_name = old_filename
464  iformatted = old_is_form
465  ELSE
466  directory = new_directory
467  file_name = new_filename
468  iformatted = new_is_form
469  END IF
470 
471  IF (is_in) THEN
472  directory_m = new_directory
473  file_name_m = new_filename
474  is_form_m = new_is_form
475  ELSE
476  directory_m = old_directory
477  file_name_m = old_filename
478  is_form_m = old_is_form
479  END IF
480 
481  !===Check coherence
482  CALL mpi_comm_size(petsc_comm_world, nb_procs, ierr)
483  IF ( nb_s*nb_f /= nb_procs ) THEN
484  CALL error_petsc('BUG: nb_S * nb_F /= nb_procs')
485  END IF
486 
487  CALL create_cart_comm((/nb_s, nb_f/),comm_cart,comm_one_d,coord_cart)
488  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
489  CALL mpi_comm_rank(comm_one_d(2),rank,code)
490  IF (nb_f/=nb_procs) THEN
491  CALL error_petsc('BUG in INIT, nb_procs_F/=nb_procs')
492  END IF
493 
494  IF ( (rw_ns) .AND. (type_pb=='mxw')) THEN
495  CALL error_petsc('WARNING: type_pb and/or rw_ns might be wrong')
496  END IF
497  IF ( (rw_mxw) .AND. ((type_pb=='nst')) ) THEN
498  CALL error_petsc('WARNING: type_pb and/or rw_mxw might be wrong')
499  END IF
500 
501  !===Prepare meshes and pointers
502  CALL load_dg_mesh_free_format(directory, file_name, list_dom, list_inter, 1, p1_mesh_glob, iformatted)
503  CALL load_dg_mesh_free_format(directory, file_name, list_dom, list_inter, 2, p2_mesh_glob, iformatted)
504  IF (if_energy) THEN
505  CALL load_dg_mesh_free_format(directory, file_name, list_dom, &
506  list_inter_temp, 2, p2_c0_mesh_glob_temp, iformatted)
507  END IF
508 
509  !===Start Metis mesh generation=================================================
510  ALLOCATE(part(p1_mesh_glob%me))
511  WRITE(tit_part,'(i4)') nb_s
512  mesh_part_name='mesh_part_S'//trim(adjustl(tit_part))//'.'//trim(adjustl(file_name))
513  IF (if_read_partition) THEN
514  OPEN(unit=51, file=mesh_part_name, status='unknown', form='formatted')
515  READ(51,*) part
516  CLOSE(51)
517  WRITE(*,*) 'read partition'
518  ELSE
519  WRITE(*,*) 'create partition'
520  CALL part_mesh_m_t_h_phi(nb_s, list_dom_ns, list_dom_temp, list_dom_h, &
521  list_dom_phi, p1_mesh_glob, list_inter, part, my_periodic)
522  IF (petsc_rank==0) THEN
523  OPEN(unit=51, file=mesh_part_name, status='replace', form='formatted')
524  WRITE(51,*) part
525  CLOSE(51)
526  END IF
527  END IF
528 
529  !===Extract local meshes from global meshes
530  IF (if_momentum) THEN
531  CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_ns,pp_mesh_glob,pp_mesh)
532  CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_ns,vv_mesh_glob,vv_mesh)
533 
534  ALLOCATE(comm_one_d_ns(2))
535  comm_one_d_ns(2) = comm_one_d(2)
536  CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
537  IF (pp_mesh%me/=0) THEN
538  CALL mpi_comm_split (comm_one_d(1),1,rank_s,comm_one_d_ns(1),code)
539  ELSE
540  CALL mpi_comm_split (comm_one_d(1),mpi_undefined,rank_s,comm_one_d_ns(1),code)
541  END IF
542  END IF
543 
544  IF (if_energy) THEN
545  CALL extract_mesh(comm_one_d(1),nb_s,p2_c0_mesh_glob_temp,part,list_dom_temp,temp_mesh_glob,temp_mesh)
546  ALLOCATE(comm_one_d_temp(2))
547  CALL mpi_comm_dup(comm_one_d(2), comm_one_d_temp(2), code)
548  CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
549  IF (temp_mesh%me/=0) THEN
550  CALL mpi_comm_split (comm_one_d(1),1,rank_s,comm_one_d_temp(1),code)
551  ELSE
552  CALL mpi_comm_split (comm_one_d(1),mpi_undefined,rank_s,comm_one_d_temp(1),code)
553  END IF
554  END IF
555 
556  IF (if_induction) THEN
557  IF (type_fe_h==1) THEN
558  CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
559  ELSE
560  CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
561  END IF
562  IF (type_fe_phi==1) THEN
563  CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_phi,phi_mesh_glob,phi_mesh)
564  ELSE
565  CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_phi,phi_mesh_glob,phi_mesh)
566  END IF
567  END IF
568 
569  !===Cleanup
570  CALL free_mesh(p1_mesh_glob)
571  CALL free_mesh(p2_mesh_glob)
572  IF (if_energy) THEN
573  DEALLOCATE(list_inter_temp)
574  CALL free_mesh(p2_c0_mesh_glob_temp)
575  END IF
576 
577  m_max_c = nb_mode/nb_f
578 
579  !===Load meshes for monoproc
580  IF (if_momentum) THEN
581  CALL free_mesh(vv_mesh_glob)
582  CALL free_mesh(pp_mesh_glob)
583  CALL load_mesh_free_format(directory_m, file_name_m, list_dom_ns, 2, vv_mesh_glob, is_form_m)
584  CALL load_mesh_free_format(directory_m, file_name_m, list_dom_ns, 1, pp_mesh_glob, is_form_m)
585  IF (check_plt) THEN
586  CALL plot_const_p1_label(vv_mesh_glob%jj, vv_mesh_glob%rr, 1.d0*vv_mesh_glob%i_d, 'vv.plt')
587  END IF
588  END IF
589  IF (if_energy) THEN
590  CALL free_mesh(temp_mesh_glob)
591  CALL load_mesh_free_format(directory_m, file_name_m, list_dom_temp, 2, temp_mesh_glob, is_form_m)
592  IF (check_plt) THEN
593  CALL plot_const_p1_label(temp_mesh_glob%jj, temp_mesh_glob%rr, 1.d0*temp_mesh_glob%i_d, 'temp.plt')
594  END IF
595  END IF
596  IF (if_induction) THEN
597  CALL free_mesh(h_mesh_glob)
598  CALL free_mesh(phi_mesh_glob)
599  CALL load_dg_mesh_free_format(directory_m, file_name_m, list_dom_h, list_inter_mu, type_fe_h, h_mesh_glob, is_form_m)
600  CALL load_mesh_free_format(directory_m, file_name_m, list_dom_phi, type_fe_phi, phi_mesh_glob, is_form_m)
601  IF (check_plt) THEN
602  CALL plot_const_p1_label(h_mesh_glob%jj, h_mesh_glob%rr, 1.d0*h_mesh_glob%i_d, 'HH.plt')
603  CALL plot_const_p1_label(phi_mesh_glob%jj, phi_mesh_glob%rr, 1.d0*phi_mesh_glob%i_d, 'phi.plt')
604  END IF
605  END IF
606 
607  !===Array allocation
608  IF (if_momentum) THEN
609  ALLOCATE(un_glob(vv_mesh_glob%np, 6, m_max_c))
610  ALLOCATE(un_m1_glob(vv_mesh_glob%np, 6, m_max_c))
611  ALLOCATE(pn_glob(pp_mesh_glob%np, 2, m_max_c))
612  ALLOCATE(pn_m1_glob(pp_mesh_glob%np, 2, m_max_c))
613  ALLOCATE(incpn_m1_glob(pp_mesh_glob%np, 2, m_max_c))
614  ALLOCATE(incpn_glob(pp_mesh_glob%np, 2, m_max_c))
615  ALLOCATE(un(vv_mesh%np, 6, m_max_c))
616  ALLOCATE(un_m1(vv_mesh%np, 6, m_max_c))
617  ALLOCATE(pn(pp_mesh%np, 2, m_max_c))
618  ALLOCATE(pn_m1(pp_mesh%np, 2, m_max_c))
619  ALLOCATE(incpn_m1(pp_mesh%np, 2, m_max_c))
620  ALLOCATE(incpn(pp_mesh%np, 2, m_max_c))
621  un_glob = 0.d0
622  un_m1_glob = 0.d0
623  pn_glob = 0.d0
624  pn_m1_glob = 0.d0
625  incpn_m1_glob = 0.d0
626  incpn_glob = 0.d0
627  un = 0.d0
628  un_m1 = 0.d0
629  pn = 0.d0
630  pn_m1 = 0.d0
631  incpn_m1= 0.d0
632  incpn = 0.d0
633  IF (if_mass) THEN
634  IF (if_level_set_p2) THEN
635  ALLOCATE(level_setn_m1_glob(nb_fluid-1,vv_mesh_glob%np, 2, m_max_c))
636  ALLOCATE(level_setn_glob(nb_fluid-1,vv_mesh_glob%np, 2, m_max_c))
637  ALLOCATE(level_setn_m1(nb_fluid-1,vv_mesh%np, 2, m_max_c))
638  ALLOCATE(level_setn(nb_fluid-1,vv_mesh%np, 2, m_max_c))
639  ELSE
640  ALLOCATE(level_setn_m1_glob(nb_fluid-1,pp_mesh_glob%np, 2, m_max_c))
641  ALLOCATE(level_setn_glob(nb_fluid-1,pp_mesh_glob%np, 2, m_max_c))
642  ALLOCATE(level_setn_m1(nb_fluid-1,pp_mesh%np, 2, m_max_c))
643  ALLOCATE(level_setn(nb_fluid-1,pp_mesh%np, 2, m_max_c))
644  END IF
645  level_setn_m1_glob = 0.d0
646  level_setn_glob = 0.d0
647  level_setn_m1 = 0.d0
648  level_setn = 0.d0
649  END IF
650  END IF
651 
652  IF (if_energy) THEN
653  ALLOCATE(tempn_m1_glob(temp_mesh_glob%np, 2, m_max_c))
654  ALLOCATE(tempn_glob(temp_mesh_glob%np, 2, m_max_c))
655  ALLOCATE(tempn_m1(temp_mesh%np, 2, m_max_c))
656  ALLOCATE(tempn(temp_mesh%np, 2, m_max_c))
657  tempn_m1_glob = 0.d0
658  tempn_glob = 0.d0
659  tempn_m1 = 0.d0
660  tempn = 0.d0
661  END IF
662 
663  IF (if_induction) THEN
664  ALLOCATE(hn1(h_mesh%np, 6, m_max_c))
665  ALLOCATE(hn(h_mesh%np, 6, m_max_c))
666  ALLOCATE(bn1(h_mesh%np, 6, m_max_c))
667  ALLOCATE(bn(h_mesh%np, 6, m_max_c))
668  ALLOCATE(phin1(phi_mesh%np,2, m_max_c))
669  ALLOCATE(phin(phi_mesh%np,2, m_max_c))
670  ALLOCATE(hn1_glob(h_mesh_glob%np, 6, m_max_c))
671  ALLOCATE(hn_glob(h_mesh_glob%np, 6, m_max_c))
672  ALLOCATE(bn1_glob(h_mesh_glob%np, 6, m_max_c))
673  ALLOCATE(bn_glob(h_mesh_glob%np, 6, m_max_c))
674  ALLOCATE(phin1_glob(phi_mesh_glob%np,2, m_max_c))
675  ALLOCATE(phin_glob(phi_mesh_glob%np,2, m_max_c))
676  hn1 = 0.d0
677  hn = 0.d0
678  bn1 = 0.d0
679  bn = 0.d0
680  phin1 = 0.d0
681  phin = 0.d0
682  hn1_glob = 0.d0
683  hn_glob = 0.d0
684  bn1_glob = 0.d0
685  bn_glob = 0.d0
686  phin1_glob = 0.d0
687  phin_glob = 0.d0
688  END IF
689 
690  !===Pointers
691  IF (is_in) THEN
692  mono_in = .false.
693  mono_out = .true.
694  IF (if_induction) THEN
695  hn_in => hn
696  hn1_in => hn1
697  bn_in => bn
698  bn1_in => bn1
699  phin_in => phin
700  phin1_in => phin1
701  h_mesh_in => h_mesh
702  phi_mesh_in => phi_mesh
703  hn_out => hn_glob
704  hn1_out => hn1_glob
705  bn_out => bn_glob
706  bn1_out => bn1_glob
707  phin_out => phin_glob
708  phin1_out => phin1_glob
709  h_mesh_out => h_mesh_glob
710  phi_mesh_out => phi_mesh_glob
711  END IF
712  IF (if_momentum) THEN
713  un_in => un
714  un_m1_in => un_m1
715  pn_in => pn
716  pn_m1_in => pn_m1
717  incpn_in => incpn
718  incpn_m1_in=> incpn_m1
719  vv_mesh_in => vv_mesh
720  pp_mesh_in => pp_mesh
721  un_out => un_glob
722  un_m1_out => un_m1_glob
723  pn_out => pn_glob
724  pn_m1_out => pn_m1_glob
725  incpn_out => incpn_glob
726  incpn_m1_out=> incpn_m1_glob
727  vv_mesh_out => vv_mesh_glob
728  pp_mesh_out => pp_mesh_glob
729  IF (if_mass) THEN
730  level_setn_in => level_setn
731  level_setn_m1_in => level_setn_m1
732  level_setn_out => level_setn_glob
733  level_setn_m1_out=> level_setn_m1_glob
734  END IF
735  END IF
736  IF (if_energy) THEN
737  tempn_in => tempn
738  tempn_m1_in => tempn_m1
739  tempn_out => tempn_glob
740  tempn_m1_out => tempn_m1_glob
741  temp_mesh_in => temp_mesh
742  temp_mesh_out => temp_mesh_glob
743  END IF
744 
745  ELSE
746  mono_in = .true.
747  mono_out = .false.
748  IF (if_induction) THEN
749  hn_in => hn_glob
750  hn1_in => hn1_glob
751  bn_in => bn_glob
752  bn1_in => bn1_glob
753  phin_in => phin_glob
754  phin1_in => phin1_glob
755  h_mesh_in => h_mesh_glob
756  phi_mesh_in => phi_mesh_glob
757  hn_out => hn
758  hn1_out => hn1
759  bn_out => bn
760  bn1_out => bn1
761  phin_out => phin
762  phin1_out => phin1
763  h_mesh_out => h_mesh
764  phi_mesh_out => phi_mesh
765  END IF
766  IF (if_momentum) THEN
767  un_in => un_glob
768  un_m1_in => un_m1_glob
769  pn_in => pn_glob
770  pn_m1_in => pn_m1_glob
771  incpn_in => incpn_glob
772  incpn_m1_in=> incpn_m1_glob
773  vv_mesh_in => vv_mesh_glob
774  pp_mesh_in => pp_mesh_glob
775  un_out => un
776  un_m1_out => un_m1
777  pn_out => pn
778  pn_m1_out => pn_m1
779  incpn_out => incpn
780  incpn_m1_out=> incpn_m1
781  vv_mesh_out => vv_mesh
782  pp_mesh_out => pp_mesh
783  IF (if_mass) THEN
784  level_setn_in => level_setn_glob
785  level_setn_m1_in => level_setn_m1_glob
786  level_setn_out => level_setn
787  level_setn_m1_out=> level_setn_m1
788  END IF
789  END IF
790  IF (if_energy) THEN
791  tempn_in => tempn_glob
792  tempn_m1_in => tempn_m1_glob
793  tempn_out => tempn
794  tempn_m1_out=> tempn_m1
795  temp_mesh_in => temp_mesh_glob
796  temp_mesh_out => temp_mesh
797  END IF
798  END IF
799 
800  !===Interpolation for Maxwell
801  IF (rw_mxw) THEN
802  IF (rank==0) WRITE(*,*) 'Start interpolation Maxwell'
803  IF (inter_mesh) THEN
804  ALLOCATE(controle_h(h_mesh_out%np), controle_phi(phi_mesh_out%np))
805  DO m = index_start, index_start+nb_fic-1
806  hn1 = 0.d0
807  hn = 0.d0
808  bn1 = 0.d0
809  bn = 0.d0
810  phin1 = 0.d0
811  phin = 0.d0
812  hn1_glob = 0.d0
813  hn_glob = 0.d0
814  bn1_glob = 0.d0
815  bn_glob = 0.d0
816  phin1_glob = 0.d0
817  phin_glob = 0.d0
818  WRITE(tit, '(i3)') m
819  lblank = eval_blank(3,tit)
820  DO l = 1, lblank - 1
821  tit(l:l) = '0'
822  END DO
823 
824  IF (petsc_rank==0) CALL system('mv suite_maxwell_I'//tit//'.'//old_filename//'suite_maxwell.'//old_filename)
825  CALL mpi_barrier( mpi_comm_world, code)
826  CALL read_restart_maxwell(comm_one_d, h_mesh_in, phi_mesh_in, time_h, list_mode, hn_in, hn1_in, &
827  bn_in, bn1_in, phin_in, phin1_in, old_filename, interpol=.false., opt_mono=mono_in)
828 
829  !===Controle_H and controle_phi are initialized to zero in interp_mesh
830  CALL interp_mesh(h_mesh_in, h_mesh_out, hn_in, hn_out, controle_h, type_fe_h)
831  CALL interp_mesh(h_mesh_in, h_mesh_out, hn1_in, hn1_out, controle_h, type_fe_h)
832  CALL interp_mesh(h_mesh_in, h_mesh_out, bn_in, bn_out, controle_h, type_fe_h)
833  CALL interp_mesh(h_mesh_in, h_mesh_out, bn1_in, bn1_out, controle_h, type_fe_h)
834  CALL interp_mesh(phi_mesh_in, phi_mesh_out, phin_in, phin_out, controle_phi, type_fe_phi)
835  CALL interp_mesh(phi_mesh_in, phi_mesh_out, phin1_in, phin1_out, controle_phi, type_fe_phi)
836 
837  IF ( min(minval(controle_h),minval(controle_phi)) == 0) THEN
838  CALL error_petsc('certains points non trouve H/phi 2')
839  END IF
840 
841  CALL write_restart_maxwell(comm_one_d, h_mesh_out, phi_mesh_out, time_h, list_mode, hn_out, hn1_out, &
842  bn_out, bn1_out, phin_out, phin1_out, new_filename, m, 1, opt_mono=mono_out)
843  CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
844  END DO
845  DEALLOCATE(controle_h, controle_phi)
846  ELSE
847  ALLOCATE(l_t_g_h(h_mesh%np), l_t_g_phi(phi_mesh%np))
848  l_t_g_h = 0
849  l_t_g_phi = 0
850  CALL loc_to_glob(h_mesh, h_mesh_glob, l_t_g_h)
851  CALL loc_to_glob(phi_mesh, phi_mesh_glob, l_t_g_phi)
852  DO m = index_start, index_start+nb_fic-1
853  hn1 = 0.d0
854  hn = 0.d0
855  bn1 = 0.d0
856  bn = 0.d0
857  phin1 = 0.d0
858  phin = 0.d0
859  hn1_glob = 0.d0
860  hn_glob = 0.d0
861  bn1_glob = 0.d0
862  bn_glob = 0.d0
863  phin1_glob = 0.d0
864  phin_glob = 0.d0
865  WRITE(tit, '(i3)') m
866  lblank = eval_blank(3,tit)
867  DO l = 1, lblank - 1
868  tit(l:l) = '0'
869  END DO
870 
871  IF (is_in) THEN
872  CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
873  WRITE(tit_s,'(i3)') rank_s
874  lblank = eval_blank(3,tit_s)
875  DO l = 1, lblank - 1
876  tit_s(l:l) = '0'
877  END DO
878  CALL system('mv suite_maxwell_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_maxwell_S'//tit_s//'.'//old_filename)
879  ELSE
880  IF (petsc_rank==0) CALL system('mv suite_maxwell_I'//tit//'.'//old_filename//'suite_maxwell.'//old_filename)
881  CALL mpi_barrier( mpi_comm_world, code)
882  END IF
883  CALL read_restart_maxwell(comm_one_d, h_mesh_in, phi_mesh_in, time_h, list_mode, hn_in, hn1_in, &
884  bn_in, bn1_in, phin_in, phin1_in, old_filename, interpol=.false., opt_mono=mono_in)
885 
886  CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, hn_in, hn_out, l_t_g_h, is_in, comm_one_d(1))
887  CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, hn1_in, hn1_out, l_t_g_h, is_in, comm_one_d(1))
888  CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, bn_in, bn_out, l_t_g_h, is_in, comm_one_d(1))
889  CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, bn1_in, bn1_out, l_t_g_h, is_in, comm_one_d(1))
890  CALL inter_mesh_loc_to_glob(phi_mesh_in, phi_mesh_out, phin_in, phin_out, l_t_g_phi, is_in, comm_one_d(1))
891  CALL inter_mesh_loc_to_glob(phi_mesh_in, phi_mesh_out, phin1_in, phin1_out, l_t_g_phi, is_in, comm_one_d(1))
892 
893  CALL write_restart_maxwell(comm_one_d, h_mesh_out, phi_mesh_out, time_h, list_mode, hn_out, hn1_out, &
894  bn_out, bn1_out, phin_out, phin1_out, new_filename, m, 1, opt_mono=mono_out)
895  CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
896  END DO
897 
898  END IF
899 
900  IF (check_plt) THEN
901  IF (h_mesh%me /= 0) THEN
902  CALL mpi_comm_rank(comm_one_d(1), rang_s, ierr)
903  WRITE(tit_s,'(i3)') rang_s
904  lblank = eval_blank(3,tit_s)
905  DO l = 1, lblank - 1
906  tit_s(l:l) = '0'
907  END DO
908  DO i = 1, SIZE(list_mode)
909  WRITE(tit_m,'(i3)') list_mode(i)
910  lblank = eval_blank(3,tit_m)
911  DO l = 1, lblank - 1
912  tit_m(l:l) = '0'
913  END DO
914  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,1,i), 'H_r_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
915  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,2,i), 'H_r_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
916  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,3,i), 'H_t_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
917  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,4,i), 'H_t_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
918  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,5,i), 'H_z_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
919  CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,6,i), 'H_z_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
920 
921  END DO
922 
923  IF (rang_s == 0) THEN
924  DO i = 1, SIZE(list_mode)
925  WRITE(tit_m,'(i3)') list_mode(i)
926  lblank = eval_blank(3,tit_m)
927  DO l = 1, lblank - 1
928  tit_m(l:l) = '0'
929  END DO
930  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,1,i), 'gH_r_cos_m='//tit_m//'_999.plt' )
931  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,2,i), 'gH_r_sin_m='//tit_m//'_999.plt' )
932  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,3,i), 'gH_t_cos_m='//tit_m//'_999.plt' )
933  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,4,i), 'gH_t_sin_m='//tit_m//'_999.plt' )
934  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,5,i), 'gH_z_cos_m='//tit_m//'_999.plt' )
935  CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,6,i), 'gH_z_sin_m='//tit_m//'_999.plt' )
936  END DO
937  END IF
938  END IF
939  END IF
940  IF (rank==0) WRITE(*,*) 'End interpolation Maxwell'
941 
942  END IF
943 
944  !===Interpolation for NS
945  IF (rw_ns) THEN
946  IF (rank==0) WRITE(*,*) 'Start interpolation NS'
947 
948  IF (inter_mesh) THEN
949 
950  ALLOCATE(controle_vv(vv_mesh_out%np), controle_pp(pp_mesh_out%np))
951  DO m = index_start, index_start+nb_fic-1
952  un_glob = 0.d0
953  un_m1_glob = 0.d0
954  pn_glob = 0.d0
955  pn_m1_glob = 0.d0
956  incpn_m1_glob = 0.d0
957  incpn_glob = 0.d0
958  un = 0.d0
959  un_m1 = 0.d0
960  pn = 0.d0
961  pn_m1 = 0.d0
962  incpn_m1= 0.d0
963  incpn = 0.d0
964  IF (if_level_set) THEN
965  level_setn_m1_glob = 0.d0
966  level_setn_glob = 0.d0
967  level_setn_m1 = 0.d0
968  level_setn = 0.d0
969  END IF
970 
971  WRITE(tit, '(i3)') m
972  lblank = eval_blank(3,tit)
973  DO l = 1, lblank - 1
974  tit(l:l) = '0'
975  END DO
976 
977  IF (petsc_rank==0) THEN
978  CALL system('mv suite_ns_I'//tit//'.'//old_filename//'suite_ns.'//old_filename)
979  END IF
980  CALL mpi_barrier( mpi_comm_world, code)
981 
982  IF (vv_mesh%me/=0) THEN
983  IF (if_level_set)THEN
984  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
985  pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename,opt_level_set=level_setn_in, &
986  opt_level_set_m1=level_setn_m1_in, opt_max_vel=max_vel, opt_mono = mono_in)
987 
988  ELSE
989  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
990  pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename, opt_mono = mono_in)
991  END IF
992  END IF
993  controle_vv = 0
994  controle_pp = 0
995  CALL interp_mesh(vv_mesh_in, vv_mesh_out, un_in, un_out, controle_vv, 2)
996  CALL interp_mesh(vv_mesh_in, vv_mesh_out, un_m1_in, un_m1_out, controle_vv, 2)
997  CALL interp_mesh(pp_mesh_in, pp_mesh_out, pn_in, pn_out, controle_pp, 1)
998  CALL interp_mesh(pp_mesh_in, pp_mesh_out, pn_m1_in, pn_m1_out, controle_pp, 1)
999  CALL interp_mesh(pp_mesh_in, pp_mesh_out, incpn_in, incpn_out, controle_pp, 1)
1000  CALL interp_mesh(pp_mesh_in, pp_mesh_out, incpn_m1_in, incpn_m1_out, controle_pp, 1)
1001  IF (if_level_set) THEN
1002  IF (if_level_set_p2) THEN
1003  DO k = 1, nb_fluid-1
1004  CALL interp_mesh(vv_mesh_in, vv_mesh_out, level_setn_in(k,:,:,:),&
1005  level_setn_out(k,:,:,:), controle_vv, 2)
1006  CALL interp_mesh(vv_mesh_in, vv_mesh_out, level_setn_m1_in(k,:,:,:),&
1007  level_setn_m1_out(k,:,:,:), controle_vv, 2)
1008  END DO
1009  ELSE
1010  DO k = 1, nb_fluid-1
1011  CALL interp_mesh(pp_mesh_in, pp_mesh_out, level_setn_in(k,:,:,:),&
1012  level_setn_out(k,:,:,:), controle_pp, 1)
1013  CALL interp_mesh(pp_mesh_in, pp_mesh_out, level_setn_m1_in(k,:,:,:),&
1014  level_setn_m1_out(k,:,:,:), controle_pp, 1)
1015  END DO
1016  END IF
1017  END IF
1018 
1019  IF ( min(minval(controle_vv),minval(controle_pp)) == 0) THEN
1020  CALL error_petsc('certains points non trouve')
1021  END IF
1022 
1023  IF (pp_mesh%me /=0) THEN
1024  IF (if_level_set) THEN
1025  CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
1026  pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_level_set=level_setn_out, &
1027  opt_level_set_m1=level_setn_m1_out, opt_max_vel=max_vel, opt_mono=mono_out)
1028  ELSE
1029  CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
1030  pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_mono=mono_out)
1031  END IF
1032  CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
1033  END IF
1034  END DO
1035  DEALLOCATE(controle_vv, controle_pp)
1036 
1037  ELSE
1038 
1039  ALLOCATE(l_t_g_vv(vv_mesh%np), l_t_g_pp(pp_mesh%np))
1040  l_t_g_vv = 0
1041  l_t_g_pp = 0
1042  CALL loc_to_glob(vv_mesh, vv_mesh_glob, l_t_g_vv)
1043  CALL loc_to_glob(pp_mesh, pp_mesh_glob, l_t_g_pp)
1044  IF (pp_mesh%me /=0) THEN
1045  DO m = index_start, index_start+nb_fic-1
1046  un_glob = 0.d0
1047  un_m1_glob = 0.d0
1048  pn_glob = 0.d0
1049  pn_m1_glob = 0.d0
1050  incpn_m1_glob = 0.d0
1051  incpn_glob = 0.d0
1052  un = 0.d0
1053  un_m1 = 0.d0
1054  pn = 0.d0
1055  pn_m1 = 0.d0
1056  incpn_m1= 0.d0
1057  incpn = 0.d0
1058  IF (if_level_set) THEN
1059  level_setn_m1_glob = 0.d0
1060  level_setn_glob = 0.d0
1061  level_setn_m1 = 0.d0
1062  level_setn = 0.d0
1063  END IF
1064 
1065  WRITE(tit, '(i3)') m
1066  lblank = eval_blank(3,tit)
1067  DO l = 1, lblank - 1
1068  tit(l:l) = '0'
1069  END DO
1070  IF (is_in) THEN
1071  CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
1072  WRITE(tit_s,'(i3)') rang_ns_s
1073  lblank = eval_blank(3,tit_s)
1074  DO l = 1, lblank - 1
1075  tit_s(l:l) = '0'
1076  END DO
1077 
1078  CALL system('mv suite_ns_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_ns_S'//tit_s//'.'//old_filename)
1079 
1080  ELSE
1081  IF (petsc_rank==0) CALL system('mv suite_ns_I'//tit//'.'//old_filename//'suite_ns.'//old_filename)
1082  CALL mpi_barrier( mpi_comm_world, code)
1083  END IF
1084 
1085  IF (if_level_set) THEN
1086  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
1087  pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename, opt_level_set=level_setn_in, &
1088  opt_level_set_m1=level_setn_m1_in, opt_max_vel=max_vel, opt_mono = mono_in)
1089  ELSE
1090  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
1091  pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename, opt_mono = mono_in)
1092  END IF
1093 
1094  CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, un_in, un_out, l_t_g_vv, is_in, &
1095  comm_one_d_ns(1))
1096  CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, un_m1_in, un_m1_out, l_t_g_vv, is_in, &
1097  comm_one_d_ns(1))
1098  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, pn_in, pn_out, l_t_g_pp, is_in, &
1099  comm_one_d_ns(1))
1100  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, pn_m1_in, pn_m1_out, l_t_g_pp, is_in, &
1101  comm_one_d_ns(1))
1102  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, incpn_in, incpn_out, l_t_g_pp, is_in, &
1103  comm_one_d_ns(1))
1104  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, incpn_m1_in, incpn_m1_out, l_t_g_pp, is_in, &
1105  comm_one_d_ns(1))
1106  IF (if_level_set) THEN
1107  IF (if_level_set_p2) THEN
1108  DO k = 1, nb_fluid-1
1109  CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, level_setn_in(k,:,:,:),&
1110  level_setn_out(k,:,:,:), l_t_g_vv, is_in, comm_one_d_ns(1))
1111  CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, level_setn_m1_in(k,:,:,:),&
1112  level_setn_m1_out(k,:,:,:), l_t_g_vv, is_in, comm_one_d_ns(1))
1113  END DO
1114  ELSE
1115  DO k = 1, nb_fluid-1
1116  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, level_setn_in(k,:,:,:),&
1117  level_setn_out(k,:,:,:), l_t_g_pp, is_in, comm_one_d_ns(1))
1118  CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, level_setn_m1_in(k,:,:,:),&
1119  level_setn_m1_out(k,:,:,:), l_t_g_pp, is_in, comm_one_d_ns(1))
1120  END DO
1121  END IF
1122  END IF
1123 
1124  IF (if_level_set) THEN
1125  CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
1126  pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_level_set=level_setn_out, &
1127  opt_level_set_m1=level_setn_m1_out, opt_max_vel=max_vel, opt_mono=mono_out)
1128  ELSE
1129  CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
1130  pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_mono=mono_out)
1131  END IF
1132  CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
1133  END DO
1134  END IF
1135  END IF
1136 
1137  IF (check_plt) THEN
1138  IF (pp_mesh%me /= 0) THEN
1139  CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
1140  WRITE(tit_s,'(i3)') rang_ns_s
1141  lblank = eval_blank(3,tit_s)
1142  DO l = 1, lblank - 1
1143  tit_s(l:l) = '0'
1144  END DO
1145  DO i = 1, SIZE(list_mode)
1146  WRITE(tit_m,'(i3)') list_mode(i)
1147  lblank = eval_blank(3,tit_m)
1148  DO l = 1, lblank - 1
1149  tit_m(l:l) = '0'
1150  END DO
1151  CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,1,i), 'u_r_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
1152  CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,3,i), 'u_t_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
1153  CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,4,i), 'u_t_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
1154  CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,5,i), 'u_z_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
1155 
1156  END DO
1157 
1158  IF (rang_ns_s == 0) THEN
1159  DO i = 1, SIZE(list_mode)
1160  WRITE(tit_m,'(i3)') list_mode(i)
1161  lblank = eval_blank(3,tit_m)
1162  DO l = 1, lblank - 1
1163  tit_m(l:l) = '0'
1164  END DO
1165  CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,1,i), 'gu_r_cos_m='//tit_m//'_999.plt' )
1166  CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,3,i), 'gu_t_cos_m='//tit_m//'_999.plt' )
1167  CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,4,i), 'gu_t_sin_m='//tit_m//'_999.plt' )
1168  CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,5,i), 'gu_z_cos_m='//tit_m//'_999.plt' )
1169  END DO
1170  END IF
1171  END IF
1172  END IF
1173  IF (rank==0) WRITE(*,*) 'End interpolation NS'
1174  END IF
1175 
1176  !===Interpolation for temperature
1177  IF (rw_temp) THEN
1178  IF (rank==0) WRITE(*,*) 'Start interpolation temperature'
1179 
1180  IF (inter_mesh) THEN
1181  ALLOCATE(controle_temp(temp_mesh_out%np))
1182  DO m = index_start, index_start+nb_fic-1
1183  tempn_m1_glob = 0.d0
1184  tempn_glob = 0.d0
1185  tempn_m1 = 0.d0
1186  tempn = 0.d0
1187  WRITE(tit, '(i3)') m
1188  lblank = eval_blank(3,tit)
1189  DO l = 1, lblank - 1
1190  tit(l:l) = '0'
1191  END DO
1192 
1193  IF (petsc_rank==0) THEN
1194  CALL system('mv suite_temp_I'//tit//'.'//old_filename//'suite_temp.'//old_filename)
1195  END IF
1196  CALL mpi_barrier( mpi_comm_world, code)
1197 
1198  IF (temp_mesh%me/=0) THEN
1199  CALL read_restart_temp(comm_one_d_temp, time_t, list_mode, tempn_in, tempn_m1_in, old_filename, &
1200  opt_mono = mono_in)
1201  END IF
1202  controle_temp = 0
1203  CALL interp_mesh(temp_mesh_in, temp_mesh_out, tempn_in, tempn_out, controle_temp, 2)
1204  CALL interp_mesh(temp_mesh_in, temp_mesh_out, tempn_m1_in, tempn_m1_out, controle_temp, 2)
1205 
1206  IF (temp_mesh%me /= 0) THEN
1207  CALL write_restart_temp(comm_one_d_temp, temp_mesh_out, time_t, &
1208  list_mode, tempn_out, tempn_m1_out, new_filename, m, 1, opt_mono = mono_out)
1209  CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
1210  END IF
1211  END DO
1212  DEALLOCATE(controle_temp)
1213 
1214  ELSE
1215 
1216  ALLOCATE(l_t_g_temp(temp_mesh%np))
1217  l_t_g_temp = 0
1218  CALL loc_to_glob(temp_mesh, temp_mesh_glob, l_t_g_temp)
1219  IF (temp_mesh%me /=0) THEN
1220  DO m = index_start, index_start+nb_fic-1
1221  tempn_m1_glob = 0.d0
1222  tempn_glob = 0.d0
1223  tempn_m1 = 0.d0
1224  tempn = 0.d0
1225 
1226  WRITE(tit, '(i3)') m
1227  lblank = eval_blank(3,tit)
1228  DO l = 1, lblank - 1
1229  tit(l:l) = '0'
1230  END DO
1231  IF (is_in) THEN
1232  CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
1233  WRITE(tit_s,'(i3)') rang_temp_s
1234  lblank = eval_blank(3,tit_s)
1235  DO l = 1, lblank - 1
1236  tit_s(l:l) = '0'
1237  END DO
1238 
1239  CALL system('mv suite_temp_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_temp_S'//tit_s//'.'//old_filename)
1240  ELSE
1241  IF (petsc_rank==0) CALL system('mv suite_temp_I'//tit//'.'//old_filename//'suite_temp.'//old_filename)
1242  CALL mpi_barrier( mpi_comm_world, code)
1243  END IF
1244 
1245  CALL read_restart_temp(comm_one_d_temp, time_t, list_mode, tempn_in, tempn_m1_in, old_filename, &
1246  opt_mono = mono_in)
1247 
1248  CALL inter_mesh_loc_to_glob(temp_mesh_in, temp_mesh_out, tempn_in, tempn_out, l_t_g_temp, is_in, &
1249  comm_one_d_temp(1))
1250  CALL inter_mesh_loc_to_glob(temp_mesh_in, temp_mesh_out, tempn_m1_in, tempn_m1_out, l_t_g_temp, is_in, &
1251  comm_one_d_temp(1))
1252 
1253  CALL write_restart_temp(comm_one_d_temp, temp_mesh_out, time_t, list_mode, tempn_out, tempn_m1_out, &
1254  new_filename, m, 1, opt_mono = mono_out)
1255  CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
1256  END DO
1257  END IF
1258  END IF
1259 
1260  IF (check_plt) THEN
1261  IF (temp_mesh%me /= 0) THEN
1262  CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
1263  WRITE(tit_s,'(i3)') rang_temp_s
1264  lblank = eval_blank(3,tit_s)
1265  DO l = 1, lblank - 1
1266  tit_s(l:l) = '0'
1267  END DO
1268  DO i = 1, SIZE(list_mode)
1269  WRITE(tit_m,'(i3)') list_mode(i)
1270  lblank = eval_blank(3,tit_m)
1271  DO l = 1, lblank - 1
1272  tit_m(l:l) = '0'
1273  END DO
1274  CALL plot_scalar_field(temp_mesh%jj, temp_mesh%rr, tempn(:,1,i), 'temp_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
1275  CALL plot_scalar_field(temp_mesh%jj, temp_mesh%rr, tempn(:,2,i), 'temp_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
1276  END DO
1277  IF (rang_temp_s == 0) THEN
1278  DO i = 1, SIZE(list_mode)
1279  WRITE(tit_m,'(i3)') list_mode(i)
1280  lblank = eval_blank(3,tit_m)
1281  DO l = 1, lblank - 1
1282  tit_m(l:l) = '0'
1283  END DO
1284  CALL plot_scalar_field(temp_mesh_glob%jj, temp_mesh_glob%rr, tempn(:,1,i), &
1285  'gtemp_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
1286  CALL plot_scalar_field(temp_mesh_glob%jj, temp_mesh_glob%rr, tempn(:,2,i), &
1287  'gtemp_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
1288  END DO
1289  END IF
1290  END IF
1291  END IF
1292  IF (rank==0) WRITE(*,*) 'End interpolation temperature'
1293  END IF
1294 
1295  IF (is_in) THEN
1296  IF (petsc_rank==0 .AND. rw_mxw) CALL system('rm -rf suite_maxwell_S*')
1297  IF (petsc_rank==0 .AND. rw_ns) CALL system('rm -rf suite_ns_S*')
1298  IF (petsc_rank==0 .AND. rw_temp) CALL system('rm -rf suite_temp_S*')
1299  ELSE
1300  IF (petsc_rank==0 .AND. rw_mxw) CALL system('rm -rf suite_maxwell.*')
1301  IF (petsc_rank==0 .AND. rw_ns) CALL system('rm -rf suite_ns.*')
1302  IF (petsc_rank==0 .AND. rw_temp) CALL system('rm -rf suite_temp.*')
1303  END IF
1304 
1305  IF (check_plt .AND. petsc_rank==0) THEN
1306  CALL system('mkdir -p PLT')
1307  CALL system('mv *.plt PLT/')
1308  END IF
1309  CALL mpi_barrier( mpi_comm_world, code)
1310 
1311  CALL petscfinalize(ierr)
1312 
1313  END SUBROUTINE mesh_interpol
1314 
1315  SUBROUTINE which_pb(is_in)
1317  IMPLICIT NONE
1318  LOGICAL, INTENT(OUT) :: is_in
1319  CHARACTER(len=200) :: inline
1320  CALL getarg(1,inline)
1321  IF (trim(adjustl(inline)) == 'petsc_nonpetsc') THEN
1322  is_in = .true.
1323  ELSE IF (trim(adjustl(inline)) == 'nonpetsc_petsc') THEN
1324  is_in = .false.
1325  ELSE
1326  CALL error_petsc('BUG in which_pb')
1327  END IF
1328  END SUBROUTINE which_pb
1329 
1330 END MODULE mesh_interpolation
subroutine write_restart_temp(communicator, temp_mesh, time, list_mode, tempn, tempn_m1, filename, it, freq_restart, opt_mono)
Definition: restart.f90:1004
subroutine find_string(unit, string, okay)
subroutine, public load_dg_mesh_free_format(dir, fil, list_dom, list_inter, type_fe, mesh, mesh_formatted)
Definition: prep_mesh.f90:276
subroutine read_restart_ns(communicator, time, list_mode, un, un_m1, pn, pn_m1, incpn, incpn_m1, filename, val_init, interpol, opt_level_set, opt_level_set_m1, opt_max_vel, opt_mono, opt_it, opt_dt)
Definition: restart.f90:110
subroutine plot_scalar_field(jj, rr, uu, file_name)
Definition: sub_plot.f90:704
subroutine error_petsc(string)
Definition: my_util.f90:16
subroutine read_restart_maxwell(communicator, H_mesh, phi_mesh, time, list_mode, Hn, Hn1, Bn, Bn1, phin, phin1, filename, val_init, interpol, opt_mono, opt_it, opt_dt)
Definition: restart.f90:781
integer function eval_blank(len_str, string)
subroutine read_until(unit, string, error)
subroutine plot_const_p1_label(jj, rr, uu, file_name)
Definition: sub_plot.f90:265
subroutine, public mesh_interpol
subroutine, public part_mesh_m_t_h_phi(nb_proc, list_u, list_T_in, list_h_in, list_phi, mesh, list_of_interfaces, part, my_periodic)
subroutine, public free_mesh(mesh)
subroutine which_pb(is_in)
subroutine interp_mesh(mesh_in, mesh_out, in_field, out_field, controle, type_fe)
subroutine inter_mesh_loc_to_glob(mesh_in, mesh_out, in_field, out_field, l_t_g, is_in, comm)
subroutine, public load_mesh_free_format(dir, fil, list_dom, type_fe, mesh, mesh_formatted, edge_stab)
Definition: prep_mesh.f90:914
subroutine, public create_cart_comm(ndim, comm_cart, comm_one_d, coord_cart)
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
subroutine write_restart_ns(communicator, vv_mesh, pp_mesh, time, list_mode, un, un_m1, pn, pn_m1, incpn, incpn_m1, filename, it, freq_restart, opt_level_set, opt_level_set_m1, opt_max_vel, opt_mono, opt_dt)
Definition: restart.f90:11
subroutine, public extract_mesh(communicator, nb_proc, mesh_glob, part, list_dom, mesh, mesh_loc)
subroutine read_restart_temp(communicator, time, list_mode, tempn, tempn_m1, filename, val_init, interpol, opt_mono)
Definition: restart.f90:1081
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 form
Definition: doc_intro.h:190