SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
initialization.f90
Go to the documentation of this file.
1 !
2 !Authors Jean-Luc Guermond, Caroline Nore, Copyrights 2005
3 !Revised June 2008, Jean-Luc Guermond
4 !Revised for PETSC, Jean-Luc Guermond, Francky Luddens, January 2011
5 !Revised July 9th 2013, JLG, Loic Cappanera, Remi Menard, Daniel Castanon
6 !Revised July 25th 2016, JLG, Loic Cappanera, Raphael Zanella
8  USE def_type_mesh
10  USE input_data
11  IMPLICIT NONE
12  PUBLIC:: initial, save_run, run_sfemans
14  PRIVATE
15 
16  !Logicals for equations-----------------------------------------------------
17  LOGICAL :: if_momentum, if_mass, if_induction, if_energy
18 
19  !Fields for Navier-Stokes---------------------------------------------------
20  TYPE(mesh_type), TARGET :: pp_mesh, vv_mesh
21  TYPE(petsc_csr_la) :: vv_1_LA, vv_2_LA, pp_1_LA
22  TYPE(petsc_csr_la) :: vv_3_LA ! for stress bc
23  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: un, un_m1
24  ! (noeuds,type,mode) composante du champ de vitesse a deux instants sur vv_mesh
25  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: pn, pn_m1
26  ! (noeuds,type,mode) composante du champ de pression a deux instants sur pp_mesh
27  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: incpn, incpn_m1
28  !---------------------------------------------------------------------------
29 
30  !Fields for level sets in Navier-Stokes-------------------------------------
31  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:,:) :: level_set, level_set_m1
32  !Fields for density---------------------------------------------------------
33  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: density, density_m1, density_m2
34  !Entropy viscosity for level-set--------------------------------------------
35  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: visc_entro_level
36  !Maximum of velocity--------------------------------------------------------
37  REAL(KIND=8) :: max_vel
38 
39  !Fields for temperature-----------------------------------------------------
40  ! (noeuds,type,mode) composante du champ de phase a deux instants sur vv_mesh
41  TYPE(mesh_type), TARGET :: temp_mesh
42  TYPE(petsc_csr_la) :: temp_1_LA
43  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: tempn, tempn_m1
44  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:) :: vol_heat_capacity_field
45  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:) :: temperature_diffusivity_field
46  !---------------------------------------------------------------------------
47 
48  !Fields for Maxwell---------------------------------------------------------
49  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: Hn, Hn1, Hext, phin, phin1
50  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: Bn, Bn1, Bext
51  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:) :: sigma_field, mu_H_field
52  TYPE(mesh_type), TARGET :: H_mesh, phi_mesh, pmag_mesh
53  TYPE(petsc_csr_la) :: LA_H, LA_pmag, LA_phi, LA_mhd
54  TYPE(interface_type), TARGET :: interface_H_mu, interface_H_phi
55  !---------------------------------------------------------------------------
56 
57  !Periodic structures--------------------------------------------------------
58  TYPE(periodic_type) :: H_phi_per
59  TYPE(periodic_type) :: vvrt_per
60  TYPE(periodic_type) :: vvrtz_per
61  TYPE(periodic_type) :: vvz_per
62  TYPE(periodic_type) :: pp_per
63  TYPE(periodic_type) :: temp_per
64  TYPE(periodic_type) :: level_set_per
65  !---------------------------------------------------------------------------
66 
67  !Coupling variables---------------------------------------------------------
68  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: v_to_Max
69  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: H_to_NS
70  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: B_to_NS
71  !October 7, 2008, JLG
72  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: T_to_NS
73  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: v_to_energy
74  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: T_to_Max
75  REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) :: H_to_energy
76  !---------------------------------------------------------------------------
77 
78  !Connectivity structures----------------------------------------------------
79  !October 7, 2008, JLG
80  INTEGER, ALLOCATABLE, DIMENSION(:) :: jj_v_to_H
81  INTEGER, ALLOCATABLE, DIMENSION(:) :: jj_v_to_temp
82  INTEGER, ALLOCATABLE, DIMENSION(:) :: jj_temp_to_H
83  !---------------------------------------------------------------------------
84 
85  !Modes list-----------------------------------------------------------------
86  INTEGER, TARGET, ALLOCATABLE, DIMENSION(:) :: list_mode
87  INTEGER :: m_max_c
88  !---------------------------------------------------------------------------
89 
90  !Names already used---------------------------------------------------------
91  REAL(KIND=8) :: time
92  REAL(KIND=8) :: R_fourier
93  INTEGER :: index_fourier
94  !---------------------------------------------------------------------------
95 
96  !Communicators for Petsc, in space and Fourier------------------------------
97 #include "petsc/finclude/petsc.h"
98  mpi_comm :: comm_cart
99  mpi_comm, DIMENSION(:), POINTER :: comm_one_d, comm_one_d_ns, comm_one_d_temp, coord_cart
100  !---------------------------------------------------------------------------
101 
102  !-------------END OF DECLARATIONS------------------------------------------
103 CONTAINS
104  !---------------------------------------------------------------------------
105 
106  !---------------------------------------------------------------------------
107  SUBROUTINE initial(vv_mesh_out, pp_mesh_out, H_mesh_out, phi_mesh_out, temp_mesh_out, &
108  interface_h_phi_out, interface_h_mu_out, list_mode_out, &
109  un_out, pn_out, hn_out, bn_out, phin_out, v_to_max_out, &
110  vol_heat_capacity_field_out, temperature_diffusivity_field_out, &
111  mu_h_field_out, sigma_field_out, &
112  time_out, m_max_c_out, comm_one_d_out, comm_one_d_ns_out, comm_one_d_temp_out, &
113  tempn_out, level_set_out, density_out)
115  IMPLICIT NONE
116  TYPE(mesh_type), POINTER :: pp_mesh_out, vv_mesh_out
117  TYPE(mesh_type), POINTER :: h_mesh_out, phi_mesh_out
118  TYPE(mesh_type), POINTER :: temp_mesh_out
119  TYPE(interface_type), POINTER :: interface_h_mu_out, interface_h_phi_out
120  INTEGER, POINTER, DIMENSION(:) :: list_mode_out
121  REAL(KIND=8), POINTER, DIMENSION(:,:,:) :: un_out, pn_out, hn_out, bn_out, phin_out, v_to_max_out, tempn_out, density_out
122  REAL(KIND=8), POINTER, DIMENSION(:,:,:,:):: level_set_out
123  REAL(KIND=8), POINTER, DIMENSION(:) :: sigma_field_out, mu_h_field_out
124  REAL(KIND=8), POINTER, DIMENSION(:) :: vol_heat_capacity_field_out, temperature_diffusivity_field_out
125  REAL(KIND=8) :: time_out
126  INTEGER :: m_max_c_out
127 #include "petsc/finclude/petsc.h"
128  mpi_comm, DIMENSION(:), POINTER :: comm_one_d_out, comm_one_d_ns_out, comm_one_d_temp_out
129 
130  CALL init
131 
132  !===Initialize meshes for vtu post processing
133  CALL sfemans_initialize_postprocessing(comm_one_d, vv_mesh, pp_mesh, h_mesh, phi_mesh, temp_mesh, &
134  list_mode, inputs%number_of_planes_in_real_space)
135 
136  vv_mesh_out => vv_mesh
137  pp_mesh_out => pp_mesh
138  h_mesh_out => h_mesh
139  phi_mesh_out => phi_mesh
140  temp_mesh_out => temp_mesh
141  interface_h_mu_out => interface_h_mu
142  interface_h_phi_out => interface_h_phi
143  list_mode_out => list_mode
144  un_out => un
145  pn_out => pn
146  tempn_out => tempn
147  level_set_out => level_set
148  density_out => density
149  hn_out => hn
150  bn_out => bn
151  phin_out => phin
152  v_to_max_out => v_to_max
153  vol_heat_capacity_field_out => vol_heat_capacity_field
154  temperature_diffusivity_field_out => temperature_diffusivity_field
155  mu_h_field_out => mu_h_field
156  sigma_field_out => sigma_field
157  time_out = time
158  m_max_c_out = m_max_c
159  comm_one_d_out => comm_one_d
160  comm_one_d_ns_out => comm_one_d_ns
161  comm_one_d_temp_out => comm_one_d_temp
162  END SUBROUTINE initial
163  !---------------------------------------------------------------------------
164 
165  !---------------------------------------------------------------------------
166  SUBROUTINE run_sfemans(time_in)
167  USE subroutine_mass
170  USE update_maxwell
171  USE input_data
172  IMPLICIT NONE
173  REAL(KIND=8), INTENT(in) :: time_in
174  REAL(KIND=8), DIMENSION(vv_mesh%np,2,SIZE(list_mode)) :: visco_dyn_m1
175  REAL(KIND=8), DIMENSION(vv_mesh%np,2,SIZE(list_mode)) :: sigma_ns
176 
177  CALL zero_out_modes
178 
179  time = time_in
180 
181  IF (if_mass) THEN
182  CALL three_level_mass(comm_one_d_ns, time, pp_1_la, vv_1_la, list_mode, pp_mesh, vv_mesh, &
183  2*un-un_m1, max_vel, level_set_per, density_m2, density_m1, density, level_set_m1, level_set,&
184  visc_entro_level)
185  CALL reconstruct_variable(comm_one_d_ns, list_mode, pp_mesh, vv_mesh, level_set_m1, &
186  inputs%dyna_visc_fluid, visco_dyn_m1)
187  IF (inputs%variation_sigma_fluid) THEN
188  CALL reconstruct_variable(comm_one_d_ns, list_mode, pp_mesh, vv_mesh, level_set, &
189  inputs%sigma_fluid, sigma_ns)
190  ELSE
191  sigma_ns = 0.d0
192  END IF
193  END IF
194 
195  IF (if_energy) THEN
196  IF (if_momentum) THEN
197  CALL projection_velocity(temp_mesh, 2*un-un_m1, jj_v_to_temp, v_to_energy)
198  END IF
199  IF (if_induction) THEN
200  CALL projection_mag_field(temp_mesh, 2*hn-hn1, jj_temp_to_h, h_to_energy)
201  END IF
202  CALL three_level_temperature(comm_one_d_temp, time, temp_1_la, inputs%dt, list_mode, &
203  temp_mesh, tempn_m1, tempn, v_to_energy, &
204  vol_heat_capacity_field, temperature_diffusivity_field, &
205  inputs%my_par_temperature, inputs%temperature_list_dirichlet_sides, temp_per)
206  END IF
207 
208  IF (if_momentum) THEN
209  IF (if_energy) THEN
210  CALL projection_temperature(vv_mesh, tempn, jj_v_to_temp, .true., t_to_ns)
211  END IF
212  IF (if_induction) THEN
213  CALL projection_mag_field(vv_mesh, 2*hn-hn1, jj_v_to_h, h_to_ns)
214  CALL projection_mag_field(vv_mesh, 2*bn-bn1, jj_v_to_h, b_to_ns)
215  END IF
216  CALL navier_stokes_decouple(comm_one_d_ns,time, vv_3_la, pp_1_la, &
217  list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, &
218  pn_m1, pn, un_m1, un, vvz_per, pp_per, h_to_ns, b_to_ns, &
219  density_m2, density_m1, density, visco_dyn_m1, t_to_ns, level_set_m1, level_set, &
220  visc_entro_level)
221  END IF
222 
223  IF (if_induction) THEN
224  IF (if_energy) THEN
225  CALL projection_temperature(h_mesh, tempn, jj_temp_to_h, .false., t_to_max)
226  END IF
227  IF (if_momentum) THEN
228  CALL projection_velocity(h_mesh, un, jj_v_to_h, v_to_max)
229  END IF
230  CALL maxwell_decouple(comm_one_d, h_mesh, pmag_mesh, phi_mesh, &
231  interface_h_phi, interface_h_mu, hn, bn, phin, hn1, bn1, phin1, v_to_max, &
232  inputs%stab, sigma_field, r_fourier, index_fourier, mu_h_field, inputs%mu_phi, &
233  time, inputs%dt, inputs%Rem, list_mode, h_phi_per, la_h, la_pmag, la_phi, la_mhd, sigma_ns, jj_v_to_h)
234  END IF
235 
236  END SUBROUTINE run_sfemans
237  !---------------------------------------------------------------------------
238 
239  !---------------------------------------------------------------------------
240  SUBROUTINE zero_out_modes
241  USE input_data
242  IMPLICIT NONE
243  LOGICAL, SAVE :: once_zero_out_mode=.true.
244  INTEGER, POINTER, DIMENSION(:), SAVE :: select_mode_ns, select_mode_mxw
245 
246  IF (once_zero_out_mode) THEN
247  once_zero_out_mode=.false.
248  IF (inputs%if_zero_out_modes) THEN
249  CALL prepare_zero_out_modes(list_mode, inputs%list_select_mode_ns, select_mode_ns)
250  CALL prepare_zero_out_modes(list_mode, inputs%list_select_mode_mxw, select_mode_mxw)
251  END IF
252  END IF
253  IF (inputs%if_zero_out_modes) THEN
254  IF (h_mesh%me /=0 .AND. SIZE(select_mode_mxw)>0) THEN
255  hn(:,:,select_mode_mxw) = 0.d0
256  hn1(:,:,select_mode_mxw) = 0.d0
257  END IF
258  IF (phi_mesh%me /= 0 .AND. SIZE(select_mode_mxw)>0) THEN
259  phin(:,:,select_mode_mxw) = 0.d0
260  phin1(:,:,select_mode_mxw) = 0.d0
261  END IF
262  IF (vv_mesh%me /=0 .AND. SIZE(select_mode_ns)>0) THEN
263  un(:,:,select_mode_ns) = 0.d0
264  pn(:,:,select_mode_ns) = 0.d0
265  incpn(:,:,select_mode_ns) = 0.d0
266  un_m1(:,:,select_mode_ns) = 0.d0
267  pn_m1(:,:,select_mode_ns) = 0.d0
268  incpn_m1(:,:,select_mode_ns) = 0.d0
269  END IF
270  END IF
271  END SUBROUTINE zero_out_modes
272  !---------------------------------------------------------------------------
273 
274  !---------------------------------------------------------------------------
275  SUBROUTINE prepare_zero_out_modes(list_mode, list_mode_to_zero_out, select_mode)
276  USE input_data
277  IMPLICIT NONE
278  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode_to_zero_out
279  INTEGER, DIMENSION(:) :: list_mode
280  INTEGER, POINTER, DIMENSION(:) :: select_mode
281  INTEGER :: i, inc
282  INTEGER, DIMENSION(1) :: kloc
283  inc = 0
284  DO i = 1, SIZE(list_mode_to_zero_out)
285  IF (minval(abs(list_mode-list_mode_to_zero_out(i)))==0) THEN
286  inc = inc + 1
287  END IF
288  END DO
289  ALLOCATE(select_mode(inc))
290  inc = 0
291  DO i = 1, SIZE(list_mode_to_zero_out)
292  IF (minval(abs(list_mode-list_mode_to_zero_out(i)))==0) THEN
293  inc = inc + 1
294  kloc = minloc(abs(list_mode-list_mode_to_zero_out(i)))
295  select_mode(inc) = kloc(1)
296  END IF
297  END DO
298  END SUBROUTINE prepare_zero_out_modes
299  !---------------------------------------------------------------------------
300 
301  !---------------------------------------------------------------------------
302  SUBROUTINE projection_velocity(mesh, vn, connectivity_structure, coupling_variable)
303  USE boundary
304  IMPLICIT NONE
305  TYPE(mesh_type), INTENT(IN) :: mesh
306  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vn
307  INTEGER, DIMENSION(:), INTENT(IN) :: connectivity_structure
308  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: coupling_variable
309  INTEGER :: i, j, k
310 
311  IF (mesh%np>0) THEN !we check that ns is a subset of temp before constructing the extended vv field
312  DO i = 1, m_max_c
313  DO k= 1, 6 !The user has to code extension_vel
314  coupling_variable(:,k,i) = extension_velocity(k, mesh, list_mode(i), time, 1)
315  END DO
316  END DO
317  END IF
318  IF (mesh%me /=0) THEN !construction of the extended field
319  DO j = 1, SIZE(connectivity_structure)
320  IF (connectivity_structure(j) == -1) cycle
321  coupling_variable(j,:,:) = vn(connectivity_structure(j),:,:)
322  END DO
323  END IF
324  END SUBROUTINE projection_velocity
325  !---------------------------------------------------------------------------
326 
327  !---------------------------------------------------------------------------
328  SUBROUTINE projection_temperature(mesh, vn, connectivity_structure, if_restriction, coupling_variable) !projection of temp on another mesh subroutine added
329  IMPLICIT NONE
330  TYPE(mesh_type), INTENT(IN) :: mesh
331  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vn
332  INTEGER, DIMENSION(:), INTENT(IN) :: connectivity_structure
333  LOGICAL, INTENT(IN) :: if_restriction
334  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: coupling_variable
335  INTEGER :: j
336 
337  coupling_variable = 0.d0
338 
339  IF (mesh%me /=0 ) THEN !construction of the restricted temperature field
340  DO j = 1, SIZE(connectivity_structure)
341  IF (connectivity_structure(j) == -1) cycle
342  IF (if_restriction) THEN
343  coupling_variable(connectivity_structure(j),:,:) = vn(j,:,:) ! restriction of T on vv_mesh
344  ELSE
345  coupling_variable(j,:,:) = vn(connectivity_structure(j),:,:) ! extension of T on H_mesh
346  END IF
347  END DO
348  END IF
349  END SUBROUTINE projection_temperature
350  !---------------------------------------------------------------------------
351 
352  !---------------------------------------------------------------------------
353  SUBROUTINE projection_mag_field(mesh, vn, connectivity_structure, coupling_variable) !projection of mag field on another mesh subroutine added
354  USE my_util
355  IMPLICIT NONE
356  TYPE(mesh_type), INTENT(IN) :: mesh
357  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: vn
358  INTEGER, DIMENSION(:), INTENT(IN) :: connectivity_structure
359  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: coupling_variable
360  INTEGER :: j, ms, ms1, ms2, m1, m2, ns, j1, j2, m
361 
362  DO j = 1, SIZE(connectivity_structure)
363  IF (connectivity_structure(j) == -1) cycle
364  coupling_variable(connectivity_structure(j),:,:) = vn(j,:,:)
365  END DO
366  DO ms = 1, interface_h_mu%mes
367  ms2 = interface_h_mu%mesh2(ms)
368  ms1 = interface_h_mu%mesh1(ms)
369  m2 = h_mesh%neighs(ms2)
370  m1 = h_mesh%neighs(ms1)
371  IF (m1> mesh%me .OR. m2> mesh%me) cycle
372  DO ns = 1, h_mesh%gauss%n_ws
373  j1 = interface_h_mu%jjs1(ns,ms)
374  j2 = interface_h_mu%jjs2(ns,ms)
375  IF (connectivity_structure(j1)/=connectivity_structure(j2)) THEN
376  WRITE(*,*) connectivity_structure(j1), connectivity_structure(j2)
377  CALL error_petsc(' BUG in mhd: connectivity_structure(j1)/=connectivity_structure(j2) ')
378  END IF
379  coupling_variable(connectivity_structure(j1),:,:) = (vn(j1,:,:)+vn(j2,:,:))/2
380  END DO
381  END DO
382  IF (h_mesh%gauss%n_w/=mesh%gauss%n_w) THEN
383  DO m = 1, mesh%me
384  coupling_variable(mesh%jj(4,m), :,:) = (coupling_variable(mesh%jj(2,m),:,:) + coupling_variable(mesh%jj(3,m),:,:))/2
385  coupling_variable(mesh%jj(5,m), :,:) = (coupling_variable(mesh%jj(3,m),:,:) + coupling_variable(mesh%jj(1,m),:,:))/2
386  coupling_variable(mesh%jj(6,m), :,:) = (coupling_variable(mesh%jj(1,m),:,:) + coupling_variable(mesh%jj(2,m),:,:))/2
387  END DO
388  END IF
389 
390  END SUBROUTINE projection_mag_field
391  !---------------------------------------------------------------------------
392 
393  !----------------SAVE RUN---------------------------------------------------
394  SUBROUTINE save_run(it, freq_restart)
395  USE restart
396  IMPLICIT NONE
397  INTEGER, INTENT(IN) :: it, freq_restart
398 
399  IF (if_momentum) THEN
400  IF (pp_mesh%me /= 0) THEN
401  IF (.NOT. if_mass) THEN
402  CALL write_restart_ns(comm_one_d_ns, vv_mesh, pp_mesh, time, &
403  list_mode, un, un_m1, pn, pn_m1, &
404  incpn, incpn_m1, inputs%file_name, it, freq_restart)
405  ELSE
406  CALL write_restart_ns(comm_one_d_ns, vv_mesh, pp_mesh, time, &
407  list_mode, un, un_m1, pn, pn_m1, &
408  incpn, incpn_m1, inputs%file_name, it, freq_restart, &
409  opt_level_set=level_set, opt_level_set_m1=level_set_m1,opt_max_vel=max_vel)
410  END IF
411  END IF
412  END IF
413 
414  IF (if_induction) THEN
415  CALL write_restart_maxwell(comm_one_d, h_mesh, phi_mesh, &
416  time, list_mode, hn, hn1, bn, bn1, &
417  phin, phin1, inputs%file_name, it, freq_restart)
418  END IF
419 
420  IF (if_energy) THEN
421  CALL write_restart_temp(comm_one_d_temp, temp_mesh, time, &
422  list_mode, tempn, tempn_m1, inputs%file_name, it, freq_restart)
423  END IF
424 
425  END SUBROUTINE save_run
426  !---------------------------------------------------------------------------
427 
428  SUBROUTINE init
429  !==================
430  USE my_util
431  USE chaine_caractere
432  USE periodic
433  USE prep_maill
435  USE restart
436  USE boundary
437  USE sub_plot
438  USE def_type_mesh
439  USE create_comm
440  USE metis_sfemans
441  USE generation_gauss
442  USE st_matrix
443  USE st_csr_mhd
444  USE matrix_type
445  USE symmetric_field
446  USE tn_axi
448  USE subroutine_mass
449  USE sft_parallele
450  IMPLICIT NONE
451  TYPE(mesh_type) :: vv_mesh_glob, pp_mesh_glob
452  TYPE(mesh_type) :: h_mesh_glob, phi_mesh_glob, pmag_mesh_glob, temp_mesh_glob
453  TYPE(mesh_type) :: p1_mesh_glob, p2_mesh_glob, p1_c0_mesh_glob, p2_c0_mesh_glob_temp
454  TYPE(interface_type) :: interface_h_phi_glob, interface_h_mu_glob
455  INTEGER, DIMENSION(:), ALLOCATABLE :: list_dom_h, list_dom_temp
456  INTEGER, DIMENSION(:), ALLOCATABLE :: list_dom, list_inter, part, list_dummy, list_inter_temp
457  INTEGER, DIMENSION(:), ALLOCATABLE :: h_in_to_new, temp_in_to_new
458  CHARACTER(len=200) :: data_file
459  CHARACTER(len=200) :: data_directory
460  CHARACTER(len=200) :: tit_part, mesh_part_name
461  CHARACTER(len=200) :: data_fichier
462  INTEGER :: nsize
463  INTEGER :: k, kp, m, n, i, j
464  INTEGER :: code, rank, rank_s, nb_procs, petsc_rank, bloc_size, m_max_pad
465  REAL(KIND=8) :: time_u, time_h, time_t, error, max_vel_s
466  LOGICAL :: ns_periodic, mxw_periodic, temp_periodic
467  CHARACTER(len=2) :: tit
468  INTEGER :: ms, ns, j1, j2
469 #include "petsc/finclude/petsc.h"
470  petscmpiint :: nb_procs_f, nb_procs_s
471 
472  !===Get numbers of processors===================================================
473  CALL mpi_comm_size(petsc_comm_world,nb_procs,code)
474  CALL mpi_comm_rank(petsc_comm_world,petsc_rank,code)
475 
476  !===Decide whether debugging or not=============================================
477  CALL sfemansinitialize
478  IF (inputs%test_de_convergence) THEN
479  IF (inputs%numero_du_test_debug<1 .OR. inputs%numero_du_test_debug>34) THEN
480  CALL error_petsc('BUG in INIT: debug_test_number is not in the correct range')
481  END IF
482  WRITE(tit,'(i2)') inputs%numero_du_test_debug
483  data_file = 'data_'//trim(adjustl(tit))
484  data_directory = inputs%data_directory_debug
485  data_fichier = trim(adjustl(data_directory))//'/debug_'//trim(adjustl(data_file))
486  ELSE
487  data_directory = '.'
488  data_file='data'
489  data_fichier = trim(adjustl(data_directory))//'/'//trim(adjustl(data_file))
490  END IF
491 
492  !===Assign correct user functions in boundary module
493  call assign_boundary
494 
495  !===Read data file==============================================================
496  CALL read_my_data(data_fichier)
497 
498  !===Debugging===================================================================
499  IF (inputs%test_de_convergence) THEN
500  inputs%directory = data_directory
501  END IF
502 
503  !===Initialization for empty vacuum=============================================
504  IF (inputs%nb_dom_phi==0) THEN
505  inputs%phi_nb_dirichlet_sides = 0
506  inputs%nb_inter = 0
507  inputs%mu_phi = 1.d0
508  inputs%type_fe_phi = -1
509  inputs%stab(2) = 0.d0
510  IF (ASSOCIATED(inputs%phi_list_dirichlet_sides)) DEALLOCATE(inputs%phi_list_dirichlet_sides)
511  ALLOCATE(inputs%phi_list_dirichlet_sides(0))
512  IF (ASSOCIATED(inputs%list_inter_H_phi)) DEALLOCATE(inputs%list_inter_H_phi)
513  ALLOCATE(inputs%list_inter_H_phi(0))
514  END IF
515 
516  !===Control inputs==============================================================
517  nb_procs_s = inputs%ndim(1)
518  nb_procs_f = inputs%ndim(2)
519  IF (inputs%ndim(1)*inputs%ndim(2)/=nb_procs) THEN
520  CALL error_petsc('BUG in INIT, nb_proc_space*nb_proc_fourier/=nb_procs')
521  END IF
522 
523  !===Create communicators========================================================
524  CALL create_cart_comm(inputs%ndim,comm_cart,comm_one_d,coord_cart)
525  CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
526  CALL mpi_comm_rank(comm_one_d(2),rank,code)
527  IF (nb_procs_f/=nb_procs) THEN
528  CALL error_petsc('BUG in INIT, nb_procs_F/=nb_procs')
529  END IF
530 
531  !===Sort out fourier modes======================================================
532  m_max_c = inputs%m_max/nb_procs_f
533  ALLOCATE(list_mode(m_max_c))
534  IF (m_max_c==0) THEN
535  CALL error_petsc('BUG in INIT, m_max_c==0')
536  END IF
537  IF (inputs%select_mode) THEN
538  DO i = 1, m_max_c
539  list_mode(i) = inputs%list_mode_lect(i + rank*m_max_c)
540  END DO
541  ELSE
542  DO i = 1, m_max_c
543  list_mode(i) = i + rank*m_max_c - 1
544  END DO
545  END IF
546 
547  !===Check periodicity===========================================================
548  IF (inputs%my_periodic%nb_periodic_pairs< 1) THEN
549  ns_periodic=.false.; mxw_periodic=.false.; temp_periodic = .false.
550  vvrtz_per%n_bord = 0; vvrt_per%n_bord = 0; vvz_per%n_bord = 0; pp_per%n_bord = 0
551  h_phi_per%n_bord = 0 ; temp_per%n_bord = 0
552  level_set_per%n_bord = 0
553  ELSE
554  ns_periodic=.true.; mxw_periodic=.true.; temp_periodic=.true.
555  END IF
556 
557  !===Creation of logicals for equations==========================================
558  if_mass = inputs%if_level_set
559  if_momentum = inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd'
560  if_induction = inputs%type_pb=='mxw' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='mxx' .OR. inputs%type_pb=='fhd'
561  if_energy = inputs%if_temperature
562 
563  !===Check mesh that vv_mesh is a subset of H_mesh===============================
564  IF (if_induction) THEN
565  ALLOCATE(list_dom_h(inputs%nb_dom_H), h_in_to_new(inputs%nb_dom_H)) ! JLG/AR Nov 17 2008
566  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
567  IF (SIZE(list_dom_h) < SIZE(inputs%list_dom_ns)) THEN
568  CALL error_petsc(' BUG: NS must be a subset of Maxwell ')
569  END IF
570  DO k = 1, inputs%nb_dom_ns
571  ! JLG/AR Nov 17 2008
572  IF (minval(abs(inputs%list_dom_H - inputs%list_dom_ns(k))) /= 0) THEN
573  CALL error_petsc(' BUG: NS must be a subset of Maxwell ')
574  END IF
575  DO kp = 1, inputs%nb_dom_H
576  IF (inputs%list_dom_H(kp) == inputs%list_dom_ns(k)) EXIT
577  END DO
578  h_in_to_new(k) = kp
579  ! JLG/AR Nov 17 2008
580  list_dom_h(k) = inputs%list_dom_ns(k)
581  END DO
582  m = inputs%nb_dom_ns
583  DO k = 1, inputs%nb_dom_H
584  IF (minval(abs(inputs%list_dom_H(k) - inputs%list_dom_ns)) == 0) cycle
585  m = m + 1
586  ! JLG/AR Nov 17 2008
587  h_in_to_new(m) = k
588  ! JLG/AR Nov 17 2008
589  list_dom_h(m) = inputs%list_dom_H(k)
590  END DO
591  IF (m/=inputs%nb_dom_H) THEN
592  CALL error_petsc(' BUG: m/=inputs%nb_dom_H ')
593  END IF
594  ELSE
595  ! JLG/AR Nov 17 2008
596  DO k = 1, inputs%nb_dom_H
597  h_in_to_new(k) = k
598  END DO
599  ! JLG/AR Nov 17 2008
600  list_dom_h = inputs%list_dom_H
601  END IF
602  END IF
603 
604  !===Check mesh that vv_mesh is a subset of temp_mesh============================
605  IF (if_energy) THEN
606  ALLOCATE(list_dom_temp(inputs%nb_dom_temp), temp_in_to_new(inputs%nb_dom_temp))
607  IF (SIZE(list_dom_temp) < SIZE(inputs%list_dom_ns)) THEN
608  CALL error_petsc(' BUG: NS must be a subset of temp ')
609  END IF
610  DO k = 1, inputs%nb_dom_ns
611  IF (minval(abs(inputs%list_dom_temp - inputs%list_dom_ns(k))) /= 0) THEN
612  CALL error_petsc(' BUG: NS must be a subset of temp ')
613  END IF
614  DO kp = 1, inputs%nb_dom_temp
615  IF (inputs%list_dom_temp(kp) == inputs%list_dom_ns(k)) EXIT
616  END DO
617  temp_in_to_new(k) = kp
618  ! JLG/AR Nov 17 2008
619  list_dom_temp(k) = inputs%list_dom_ns(k)
620  END DO
621  m = inputs%nb_dom_ns
622  DO k = 1, inputs%nb_dom_temp
623  IF (minval(abs(inputs%list_dom_temp(k) - inputs%list_dom_ns)) == 0) cycle
624  m = m + 1
625  ! JLG/AR Nov 17 2008
626  temp_in_to_new(m) = k
627  ! JLG/AR Nov 17 2008
628  list_dom_temp(m) = inputs%list_dom_temp(k)
629  END DO
630  IF (m/=inputs%nb_dom_temp) THEN
631  CALL error_petsc(' BUG: m/=inputs%nb_dom_temp ')
632  END IF
633  END IF
634 
635  !===Check mesh that temp_mesh is a subset of H_mesh=============================
636  IF ((if_energy).AND.(if_induction)) THEN
637  IF (SIZE(list_dom_h) < SIZE(list_dom_temp)) THEN
638  CALL error_petsc(' BUG: temp must be a subset of H ')
639  END IF
640  DO k = 1, inputs%nb_dom_temp
641  IF (minval(abs(list_dom_h - list_dom_temp(k))) /= 0) THEN
642  CALL error_petsc(' BUG: temp must be a subset of H ')
643  END IF
644  END DO
645  END IF
646 
647  !===Create interfaces in meshes=================================================
648  IF (if_momentum .AND. (.NOT. if_induction)) THEN
649  IF (if_energy) THEN
650  nsize = SIZE(list_dom_temp)
651  ALLOCATE(list_dom(nsize))
652  list_dom = list_dom_temp
653  ALLOCATE(list_inter(SIZE(inputs%list_inter_v_T)))
654  list_inter = inputs%list_inter_v_T
655  ELSE
656  nsize = SIZE(inputs%list_dom_ns)
657  ALLOCATE(list_dom(nsize))
658  list_dom = inputs%list_dom_ns
659  ALLOCATE(list_inter(0))
660  END IF
661  ELSE
662  nsize = SIZE(list_dom_h)+SIZE(inputs%list_dom_phi)
663  ALLOCATE(list_dom(nsize))
664  list_dom(1:SIZE(list_dom_h)) = list_dom_h
665  IF (SIZE(inputs%list_dom_phi)>0) THEN
666  list_dom(SIZE(list_dom_h)+1:) = inputs%list_dom_phi
667  END IF
668  nsize = SIZE(inputs%list_inter_mu)+SIZE(inputs%list_inter_H_phi)
669  ALLOCATE(list_inter(nsize))
670  IF (SIZE(inputs%list_inter_mu)>0) THEN
671  list_inter(1:SIZE(inputs%list_inter_mu)) = inputs%list_inter_mu
672  END IF
673  IF (SIZE(inputs%list_inter_H_phi)>0) THEN
674  list_inter(SIZE(inputs%list_inter_mu)+1:) = inputs%list_inter_H_phi
675  END IF
676  END IF
677  !WRITE(*,*) 'before load_dg_mesh'
678  IF (if_energy) THEN
679  ALLOCATE(list_inter_temp(0))
680  END IF
681 
682  !===Create meshes===============================================================
683  CALL load_dg_mesh_free_format(inputs%directory, inputs%file_name, list_dom, &
684  list_inter, 1, p1_mesh_glob, inputs%iformatted)
685  CALL load_dg_mesh_free_format(inputs%directory, inputs%file_name, list_dom, &
686  list_inter, 2, p2_mesh_glob, inputs%iformatted)
687  IF (if_energy) THEN
688  CALL load_dg_mesh_free_format(inputs%directory, inputs%file_name, list_dom, &
689  list_inter_temp, 2, p2_c0_mesh_glob_temp, inputs%iformatted)
690  END IF
691  IF (if_induction) THEN
692  ALLOCATE(list_dummy(0))
693  CALL load_dg_mesh_free_format(inputs%directory, inputs%file_name, list_dom, &
694  inputs%list_inter_H_phi, 1, p1_c0_mesh_glob, inputs%iformatted)
695  END IF
696 
697  !===Start Metis mesh generation=================================================
698  ALLOCATE(part(p1_mesh_glob%me))
699  WRITE(tit_part,'(i4)') inputs%ndim(1)
700  mesh_part_name='mesh_part_S'//trim(adjustl(tit_part))//'.'//trim(adjustl(inputs%file_name))
701  IF (inputs%if_read_partition) THEN
702  OPEN(unit=51, file=mesh_part_name, status='unknown', form='formatted')
703  READ(51,*) part
704  CLOSE(51)
705  WRITE(*,*) 'read partition'
706  ELSE
707  WRITE(*,*) 'create partition'
708  CALL part_mesh_m_t_h_phi(nb_procs_s, inputs%list_dom_ns, inputs%list_dom_temp, inputs%list_dom_H, &
709  inputs%list_dom_phi, p1_mesh_glob,list_inter,part, inputs%my_periodic)
710  IF (petsc_rank==0) THEN
711  OPEN(unit=51, file=mesh_part_name, status='replace', form='formatted')
712  WRITE(51,*) part
713  CLOSE(51)
714  END IF
715  END IF
716 
717  !===Extract local meshes from global meshes=====================================
718  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
719  CALL extract_mesh(comm_one_d(1),nb_procs_s,p1_mesh_glob,part,inputs%list_dom_ns,pp_mesh_glob,pp_mesh)
720  CALL extract_mesh(comm_one_d(1),nb_procs_s,p2_mesh_glob,part,inputs%list_dom_ns,vv_mesh_glob,vv_mesh)
721  ALLOCATE(comm_one_d_ns(2))
722  !Correction FL-ID, 4/4/13
723 !!$ comm_one_d_ns(2) = comm_one_d(2)
724  CALL mpi_comm_dup(comm_one_d(2), comm_one_d_ns(2), code)
725  !Correction FL-ID, 4/4/13
726  CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
727  IF (pp_mesh%me/=0) THEN
728  CALL mpi_comm_split(comm_one_d(1),1,rank_s,comm_one_d_ns(1),code)
729  ELSE
730  CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_ns(1),code)
731  END IF
732 
733  !===Test whether pp_mesh and vv_mesh meshes coincide=========================
734  DO m = 1, vv_mesh%me
735  DO n = 1, 3
736  IF (minval(abs(vv_mesh%rr(1,vv_mesh%jj(n,m)) &
737  -pp_mesh%rr(1,pp_mesh%jj(:,m))))>1.d-16 &
738  .OR. minval(abs(vv_mesh%rr(2,vv_mesh%jj(n,m))&
739  -pp_mesh%rr(2,pp_mesh%jj(:,m))))>1.d-16) THEN
740  CALL error_petsc('BUG in INIT, vv and pp global meshes are different')
741  END IF
742  END DO
743  END DO
744 
745  !===Create periodic structures===============================================
746  IF (ns_periodic) THEN
747  CALL prep_periodic(inputs%my_periodic, pp_mesh, pp_per)
748  CALL prep_periodic(inputs%my_periodic, vv_mesh, vvz_per)
749  CALL prep_periodic_bloc(inputs%my_periodic, vv_mesh, vvrt_per, 2)
750  CALL prep_periodic_bloc(inputs%my_periodic, vv_mesh, vvrtz_per, 3)
751  IF (inputs%if_level_set_P2) THEN
752  CALL prep_periodic(inputs%my_periodic, vv_mesh, level_set_per)
753  ELSE
754  CALL prep_periodic(inputs%my_periodic, pp_mesh, level_set_per)
755  END IF
756  END IF
757 
758  !===Create global csr structure==============================================
759  IF (pp_mesh%me/=0) THEN
760  CALL st_aij_csr_glob_block(comm_one_d_ns(1),1,vv_mesh_glob,vv_mesh,vv_1_la, opt_per=vvz_per)
761  CALL st_aij_csr_glob_block(comm_one_d_ns(1),2,vv_mesh_glob,vv_mesh,vv_2_la, opt_per=vvrt_per)
762  CALL st_aij_csr_glob_block(comm_one_d_ns(1),3,vv_mesh_glob,vv_mesh,vv_3_la, opt_per=vvrtz_per)
763  CALL st_aij_csr_glob_block(comm_one_d_ns(1),1,pp_mesh_glob,pp_mesh,pp_1_la, opt_per=pp_per)
764  END IF
765 
766  !===Create symmetric points==================================================
767  IF (inputs%is_mesh_symmetric) THEN
768  ALLOCATE(vv_mz_la(vv_mesh%np))
769  CALL symmetric_points(vv_mesh, vv_mesh_glob, vv_mz_la)
770  END IF
771 
772  !===Deallocate global meshes=================================================
773  CALL free_mesh(vv_mesh_glob)
774  CALL free_mesh(pp_mesh_glob)
775 
776  !===Start Gauss points generation============================================
777  CALL gen_gauss(vv_mesh,edge_stab=.false.)
778  CALL gen_gauss(pp_mesh,edge_stab=.false.)
779 
780  END IF
781 
782  !===Extract local meshes from global meshes for Maxwell=========================
783  IF (if_induction) THEN
784  CALL extract_mesh(comm_one_d(1),nb_procs_s,p1_c0_mesh_glob,part,list_dom_h,pmag_mesh_glob,pmag_mesh)
785  IF (inputs%type_fe_H==1) THEN
786  CALL extract_mesh(comm_one_d(1),nb_procs_s,p1_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
787  ELSE
788  CALL extract_mesh(comm_one_d(1),nb_procs_s,p2_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
789  END IF
790  IF (inputs%type_fe_phi==1) THEN
791  CALL extract_mesh(comm_one_d(1),nb_procs_s,p1_mesh_glob,part,inputs%list_dom_phi,phi_mesh_glob,phi_mesh)
792  ELSE
793  CALL extract_mesh(comm_one_d(1),nb_procs_s,p2_mesh_glob,part,inputs%list_dom_phi,phi_mesh_glob,phi_mesh)
794  END IF
795  END IF
796 
797  !===Extract local meshes from global meshes for temperature=====================
798  IF (if_energy) THEN
799  CALL extract_mesh(comm_one_d(1),nb_procs_s,p2_c0_mesh_glob_temp,part,list_dom_temp,temp_mesh_glob,temp_mesh)
800  ALLOCATE(comm_one_d_temp(2))
801  CALL mpi_comm_dup(comm_one_d(2), comm_one_d_temp(2), code)
802  CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
803  IF (temp_mesh%me/=0) THEN
804  CALL mpi_comm_split(comm_one_d(1),1,rank_s,comm_one_d_temp(1),code)
805  ELSE
806  CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_temp(1),code)
807  END IF
808  END IF
809 
810  !===Deallocate global meshes====================================================
811  CALL free_mesh(p1_mesh_glob)
812  CALL free_mesh(p2_mesh_glob)
813  IF (if_induction) THEN
814  DEALLOCATE(list_dummy)
815  CALL free_mesh(p1_c0_mesh_glob)
816  END IF
817  IF (if_energy) THEN
818  DEALLOCATE(list_inter_temp)
819  CALL free_mesh(p2_c0_mesh_glob_temp)
820  END IF
821 
822  !===Specific to induction equation==============================================
823  IF (if_induction) THEN
824 
825  !===Verify that pmag_mesh and H_mesh coincide================================
826  IF (pmag_mesh%me/=0) THEN
827  error = 0.d0
828  DO k = 1, 2
829  DO n = 1, SIZE(pmag_mesh%jj,1)
830  error = error + maxval(abs(pmag_mesh%rr(k,pmag_mesh%jj(n,:))-h_mesh%rr(k,h_mesh%jj(n,1:pmag_mesh%me))))
831  END DO
832  END DO
833  IF (error/maxval(abs(h_mesh%rr(1,1) -h_mesh%rr(1,:))) .GE. 5.d-14) THEN
834  CALL error_petsc(.GE.'BUG in INIT, (error/MAXVAL(ABS(H_mesh%rr(1,1) -H_mesh%rr(1,:))) 5.d-14')
835  END IF
836  END IF
837 
838  !===Verify if the two meshes coincide on the NS domain=======================
839  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
840  IF (vv_mesh%me/=0) THEN
841  error = 0.d0
842  DO k = 1, 2
843  DO n = 1, SIZE(h_mesh%jj,1)
844  error = error + maxval(abs(vv_mesh%rr(k,vv_mesh%jj(n,:))-h_mesh%rr(k,h_mesh%jj(n,1:vv_mesh%me))))
845  END DO
846  END DO
847  IF (error/maxval(abs(h_mesh%rr(1,1) -h_mesh%rr(1,:))) .GE. 5.d-14) THEN
848  CALL error_petsc(.GE.'BUG in INIT, (error/MAXVAL(ABS(H_mesh%rr(1,1) -H_mesh%rr(1,:))) 5.d-14')
849  END IF
850 
851  error = error + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(4,1:vv_mesh%me)) &
852  -(h_mesh%rr(1,h_mesh%jj(2,1:vv_mesh%me))+h_mesh%rr(1,h_mesh%jj(3,1:vv_mesh%me)))/2))&
853  + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(5,:)) &
854  -(h_mesh%rr(1,h_mesh%jj(3,1:vv_mesh%me))+h_mesh%rr(1,h_mesh%jj(1,1:vv_mesh%me)))/2))&
855  + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(6,:)) &
856  -(h_mesh%rr(1,h_mesh%jj(1,1:vv_mesh%me))+h_mesh%rr(1,h_mesh%jj(2,1:vv_mesh%me)))/2))&
857  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(4,:)) &
858  -(h_mesh%rr(2,h_mesh%jj(2,1:vv_mesh%me))+h_mesh%rr(2,h_mesh%jj(3,1:vv_mesh%me)))/2))&
859  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(5,:)) &
860  -(h_mesh%rr(2,h_mesh%jj(3,1:vv_mesh%me))+h_mesh%rr(2,h_mesh%jj(1,1:vv_mesh%me)))/2))&
861  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(6,:)) &
862  -(h_mesh%rr(2,h_mesh%jj(1,1:vv_mesh%me))+h_mesh%rr(2,h_mesh%jj(2,1:vv_mesh%me)))/2))
863  IF (error/maxval(abs(h_mesh%rr(1,1) -h_mesh%rr(1,:))) .GE. 5.d-14) THEN
864  WRITE(*,*) ' WARNING: vv_mesh and H_mesh do not coincide on the NS domain.'
865  WRITE(*,*) ' WARNING: Either you use curved elements P2 elements or BUG, ', &
866  error/maxval(abs(h_mesh%rr(1,1) -h_mesh%rr(1,:)))
867  END IF
868 
869  error=0.d0
870  DO k = 1, vv_mesh%me
871  DO n = 1, 2
872  error = error+ maxval(abs(vv_mesh%rr(n,vv_mesh%jj(1:3,k))-h_mesh%rr(n,h_mesh%jj(1:3,k))))
873  END DO
874  END DO
875  IF (error/maxval(abs(h_mesh%rr(1,1) -h_mesh%rr(1,:))) .GE. 5.d-14) THEN
876  CALL error_petsc('vv_mesh and H_mesh do NOT have the same P1 nodes')
877  END IF
878 
879  END IF
880  END IF
881 
882  !===Create interface structures==============================================
883  CALL load_interface(h_mesh_glob, h_mesh_glob, inputs%list_inter_mu, interface_h_mu_glob, .false.)
884  CALL load_interface(h_mesh_glob, phi_mesh_glob, inputs%list_inter_H_phi, interface_h_phi_glob, .true.)
885 
886  IF (h_mesh%me /=0) THEN
887  CALL load_interface(h_mesh, h_mesh, inputs%list_inter_mu, interface_h_mu, .false.)
888  ELSE
889  interface_h_mu%mes = 0
890  END IF
891  IF (h_mesh%me * phi_mesh%me /=0) THEN
892  CALL load_interface(h_mesh, phi_mesh, inputs%list_inter_H_phi, interface_h_phi, .true.)
893  ELSE
894  interface_h_phi%mes = 0
895  END IF
896 
897  !===Create periodic structures for Maxwell===================================
898  IF (mxw_periodic) THEN
899  CALL prep_periodic_h_p_phi_bc(inputs%my_periodic, h_mesh, pmag_mesh, phi_mesh, h_phi_per)
900  WRITE(*,*) 'periodic MHD done'
901  END IF
902 
903  !===Create global csr structure==============================================
904  CALL st_scr_maxwell_mu_h_p_phi(comm_one_d(1), h_mesh_glob, h_mesh, pmag_mesh_glob, pmag_mesh, &
905  phi_mesh_glob, phi_mesh, interface_h_phi_glob, interface_h_mu_glob, &
906  la_h, la_pmag, la_phi, la_mhd, opt_per=h_phi_per)
907 
908  !===Prepare csr structure for post processing grad phi=======================+++
909  CALL st_aij_csr_glob_block(comm_one_d(1),1,phi_mesh_glob,phi_mesh, vizu_grad_phi_la)
910 
911  IF (inputs%is_mesh_symmetric) THEN
912  ALLOCATE(h_mz_la(h_mesh%np))
913  CALL symmetric_points(h_mesh, h_mesh_glob, h_mz_la)
914  END IF
915 
916  !===Deallocate global meshes=================================================
917  CALL free_mesh(h_mesh_glob)
918  CALL free_mesh(pmag_mesh_glob)
919  CALL free_mesh(phi_mesh_glob)
920  CALL free_interface(interface_h_mu_glob)
921  CALL free_interface(interface_h_phi_glob)
922 
923  !===Start Gauss points generation============================================
924  CALL gen_gauss(h_mesh,edge_stab=.false.)
925  CALL gen_gauss(pmag_mesh,edge_stab=.false.)
926  CALL gen_gauss(phi_mesh,edge_stab=.false.)
927 
928  !===Create sigma_field=======================================================
929  ALLOCATE(sigma_field(h_mesh%me)) !sigma field is P0 and defined on H_mesh
930  DO m = 1, h_mesh%me
931  DO k=1, inputs%nb_dom_H
932  IF (h_mesh%i_d(m) == list_dom_h(k)) THEN
933  sigma_field(m) = inputs%sigma(h_in_to_new(k))
934  ENDIF
935  ENDDO
936  END DO
937 
938  !===Create mu_H_field========================================================
939  ALLOCATE(mu_h_field(h_mesh%np)) !mu_H_field is defined at nodes of H_mesh
940  IF (inputs%analytical_permeability) THEN !JLG + DCQ, June 26 2013
941  mu_h_field = mu_bar_in_fourier_space(h_mesh,1,h_mesh%np)
942  ELSE
943  DO m = 1, h_mesh%me
944  DO k=1, inputs%nb_dom_H
945  IF (h_mesh%i_d(m) == list_dom_h(k)) THEN
946  mu_h_field(h_mesh%jj(:,m)) = inputs%mu_H(h_in_to_new(k))
947  ENDIF
948  ENDDO
949  END DO
950  END IF
951  !===Create mu_H_field========================================================
952  !===Artificial boundary condition on phi on sphere of radius R
953  !===d(phi)/dR + (1/R)*phi = 0. Assumes that phi=C/r at infinity
954  !===Feature is currently disabled.
955  r_fourier=-1.d0 !Negative radius disables the boundary condition
956  index_fourier=0 !Index of spherical boundary where Fourier BC enforced
957 
958  END IF
959 
960  !===Specific to temperature=====================================================
961  IF (if_energy) THEN
962  !===Verify if the two meshes coincide on the NS domain=======================
963  IF (vv_mesh%me/=0) THEN
964  error = 0.d0
965  DO k = 1, 2
966  DO n = 1, SIZE(temp_mesh%jj,1)
967  error = error + maxval(abs(vv_mesh%rr(k,vv_mesh%jj(n,:))-temp_mesh%rr(k,temp_mesh%jj(n,1:vv_mesh%me))))
968  END DO
969  END DO
970  IF (error/maxval(abs(temp_mesh%rr(1,1) -temp_mesh%rr(1,:))) .GE. 5.d-14) THEN
971  CALL error_petsc(.GE.'BUG in INIT, (error/MAXVAL(ABS(temp_mesh%rr(1,1) -temp_mesh%rr(1,:))) 5.d-14')
972  END IF
973 
974  error = error + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(4,1:vv_mesh%me)) &
975  -(temp_mesh%rr(1,temp_mesh%jj(2,1:vv_mesh%me))+temp_mesh%rr(1,temp_mesh%jj(3,1:vv_mesh%me)))/2))&
976  + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(5,:)) &
977  -(temp_mesh%rr(1,temp_mesh%jj(3,1:vv_mesh%me))+temp_mesh%rr(1,temp_mesh%jj(1,1:vv_mesh%me)))/2))&
978  + maxval(abs(vv_mesh%rr(1,vv_mesh%jj(6,:)) &
979  -(temp_mesh%rr(1,temp_mesh%jj(1,1:vv_mesh%me))+temp_mesh%rr(1,temp_mesh%jj(2,1:vv_mesh%me)))/2))&
980  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(4,:)) &
981  -(temp_mesh%rr(2,temp_mesh%jj(2,1:vv_mesh%me))+temp_mesh%rr(2,temp_mesh%jj(3,1:vv_mesh%me)))/2))&
982  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(5,:)) &
983  -(temp_mesh%rr(2,temp_mesh%jj(3,1:vv_mesh%me))+temp_mesh%rr(2,temp_mesh%jj(1,1:vv_mesh%me)))/2))&
984  + maxval(abs(vv_mesh%rr(2,vv_mesh%jj(6,:)) &
985  -(temp_mesh%rr(2,temp_mesh%jj(1,1:vv_mesh%me))+temp_mesh%rr(2,temp_mesh%jj(2,1:vv_mesh%me)))/2))
986  IF (error/maxval(abs(temp_mesh%rr(1,1) -temp_mesh%rr(1,:))) .GE. 5.d-14) THEN
987  WRITE(*,*) ' WARNING: vv_mesh and temp_mesh do not coincide on the NS domain.'
988  WRITE(*,*) ' WARNING: Either you use curved elements P2 elements or BUG, ', &
989  error/maxval(abs(temp_mesh%rr(1,1) -temp_mesh%rr(1,:)))
990  END IF
991 
992  error=0.d0
993  DO k = 1, vv_mesh%me
994  DO n = 1, 2
995  error = error+ maxval(abs(vv_mesh%rr(n,vv_mesh%jj(1:3,k))-temp_mesh%rr(n,temp_mesh%jj(1:3,k))))
996  END DO
997  END DO
998  IF (error/maxval(abs(temp_mesh%rr(1,1) -temp_mesh%rr(1,:))) .GE. 5.d-14) THEN
999  CALL error_petsc('vv_mesh and temp_mesh do NOT have the same P1 nodes')
1000  END IF
1001 
1002  END IF
1003 
1004  !===Create periodic structures for temperature===============================
1005  IF (temp_periodic) THEN
1006  CALL prep_periodic(inputs%my_periodic, temp_mesh, temp_per)
1007  END IF
1008 
1009  !===Create csr structure for temperature=====================================
1010  CALL st_aij_csr_glob_block(comm_one_d_temp(1),1,temp_mesh_glob,temp_mesh,temp_1_la, opt_per=temp_per)
1011 
1012  !===Deallocate global meshes=================================================
1013  CALL free_mesh(temp_mesh_glob)
1014 
1015  !===Start Gauss points generation============================================
1016  CALL gen_gauss(temp_mesh,edge_stab=.false.)
1017 
1018  !===Create temperature_diffusivity_field=====================================
1019  ALLOCATE(temperature_diffusivity_field(temp_mesh%me))
1020  DO m = 1, temp_mesh%me
1021  DO k=1, inputs%nb_dom_temp
1022  IF (temp_mesh%i_d(m) == list_dom_temp(k)) THEN
1023  temperature_diffusivity_field(m) = inputs%temperature_diffusivity(temp_in_to_new(k))
1024  END IF
1025  END DO
1026  END DO
1027 
1028  !===Create vol_heat_capacity_field===========================================
1029  ALLOCATE(vol_heat_capacity_field(temp_mesh%me))
1030  DO m = 1, temp_mesh%me
1031  DO k=1, inputs%nb_dom_temp
1032  IF (temp_mesh%i_d(m) == list_dom_temp(k)) THEN
1033  vol_heat_capacity_field(m) = inputs%vol_heat_capacity(temp_in_to_new(k))
1034  END IF
1035  END DO
1036  END DO
1037  END IF
1038 
1039  !===Check coherence of vv_mesh and H_mesh=======================================
1040  IF ((if_induction .AND. if_momentum) .OR. inputs%type_pb=='mxx') THEN
1041  IF (vv_mesh%me /=0) THEN
1042  DO m = 1, vv_mesh%me
1043  IF (maxval(abs(h_mesh%rr(:,h_mesh%jj(1:3,m))-vv_mesh%rr(:,vv_mesh%jj(1:3,m))))/=0.d0) THEN
1044  CALL error_petsc( ' BUG in init: H_mesh and vv_mesh do not coincide ')
1045  END IF
1046  END DO
1047  END IF
1048  END IF
1049  IF (ALLOCATED(list_dom_h)) DEALLOCATE(list_dom_h)
1050  IF (ASSOCIATED(inputs%list_dom_ns)) DEALLOCATE(inputs%list_dom_ns)
1051 
1052  !===Check coherence of vv_mesh and temp_mesh====================================
1053  IF (if_energy) THEN
1054  IF (vv_mesh%me /=0) THEN
1055  DO m = 1, vv_mesh%me
1056  IF (maxval(abs(temp_mesh%rr(:,temp_mesh%jj(1:3,m))-vv_mesh%rr(:,vv_mesh%jj(1:3,m))))/=0.d0) THEN
1057  CALL error_petsc( ' BUG in init: temp_mesh and vv_mesh do not coincide ')
1058  END IF
1059  END DO
1060  END IF
1061  END IF
1062  IF (ALLOCATED(list_dom_temp)) DEALLOCATE(list_dom_temp)
1063 
1064  !===Compute local mesh size for stabilization================================
1065  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
1066  CALL compute_local_mesh_size(pp_mesh)
1067  CALL compute_local_mesh_size(vv_mesh)
1069  END IF
1070  IF (if_induction) THEN
1071  CALL compute_local_mesh_size(pmag_mesh)
1072  CALL compute_local_mesh_size(phi_mesh)
1073  CALL compute_local_mesh_size(h_mesh)
1074  END IF
1075  IF (if_energy) THEN
1076  CALL compute_local_mesh_size(temp_mesh)
1077  END IF
1078 
1079  !===Allocate arrays for Navier-Stokes===========================================
1080  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
1081  ALLOCATE(un_m1(vv_mesh%np, 6, m_max_c))
1082  ALLOCATE(un(vv_mesh%np, 6, m_max_c))
1083  ALLOCATE(pn_m1(pp_mesh%np, 2, m_max_c))
1084  ALLOCATE(pn(pp_mesh%np, 2, m_max_c))
1085  ALLOCATE(incpn_m1(pp_mesh%np, 2, m_max_c))
1086  ALLOCATE(incpn(pp_mesh%np, 2, m_max_c))
1087  ALLOCATE(density_m2(vv_mesh%np, 2, m_max_c))
1088  ALLOCATE(density_m1(vv_mesh%np, 2, m_max_c))
1089  ALLOCATE(density(vv_mesh%np, 2, m_max_c))
1090  !===Allocate arrays for Level sets===========================================
1091  IF (if_mass) THEN
1092  IF (inputs%if_level_set_P2) THEN
1093  ALLOCATE(level_set_m1(inputs%nb_fluid-1, vv_mesh%np, 2, m_max_c))
1094  ALLOCATE(level_set(inputs%nb_fluid-1, vv_mesh%np, 2, m_max_c))
1095  CALL mpi_comm_size(comm_one_d_ns(2), nb_procs, code)
1096  bloc_size = vv_mesh%gauss%l_G*vv_mesh%dom_me/nb_procs+1
1097  bloc_size = vv_mesh%gauss%l_G*(bloc_size/vv_mesh%gauss%l_G)+vv_mesh%gauss%l_G
1098  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1099  ALLOCATE(visc_entro_level(2*m_max_pad-1,bloc_size))
1100  ELSE
1101  ALLOCATE(level_set_m1(inputs%nb_fluid-1, pp_mesh%np, 2, m_max_c))
1102  ALLOCATE(level_set(inputs%nb_fluid-1, pp_mesh%np, 2, m_max_c))
1103  CALL mpi_comm_size(comm_one_d_ns(2), nb_procs, code)
1104  bloc_size = pp_mesh%gauss%l_G*pp_mesh%dom_me/nb_procs+1
1105  bloc_size = pp_mesh%gauss%l_G*(bloc_size/pp_mesh%gauss%l_G)+pp_mesh%gauss%l_G
1106  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1107  ALLOCATE(visc_entro_level(2*m_max_pad-1,bloc_size))
1108  END IF
1109  END IF
1110  END IF
1111 
1112  !===Allocate arrays for Maxwell=================================================
1113  IF (if_induction) THEN
1114  ALLOCATE(hn1(h_mesh%np, 6, m_max_c))
1115  ALLOCATE(hn(h_mesh%np, 6, m_max_c))
1116  ALLOCATE(bn1(h_mesh%np, 6, m_max_c))
1117  ALLOCATE(bn(h_mesh%np, 6, m_max_c))
1118  ALLOCATE(phin1(phi_mesh%np,2, m_max_c))
1119  ALLOCATE(phin(phi_mesh%np,2, m_max_c))
1120  END IF
1121 
1122  !===Allocate arrays for temperature=============================================
1123  IF (if_energy) THEN
1124  ALLOCATE(tempn_m1(temp_mesh%np, 2, m_max_c))
1125  ALLOCATE(tempn(temp_mesh%np, 2, m_max_c))
1126  END IF
1127 
1128  !===Create data structure jj_v_to_H=============================================
1129  IF (if_induction) THEN
1130  ALLOCATE(v_to_max(h_mesh%np, 6, m_max_c))
1131  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
1132  ALLOCATE(jj_v_to_h(h_mesh%np))
1133  jj_v_to_h = -1
1134  DO m = 1, vv_mesh%me
1135  jj_v_to_h(h_mesh%jj(:,m)) = vv_mesh%jj(1:h_mesh%gauss%n_w,m)
1136  END DO
1137 
1138  DO ms = 1, interface_h_mu%mes
1139  DO ns = 1, h_mesh%gauss%n_ws
1140  j1 = interface_h_mu%jjs1(ns,ms)
1141  j2 = interface_h_mu%jjs2(ns,ms)
1142  IF (jj_v_to_h(j1) == -1) THEN
1143  jj_v_to_h(j1) = jj_v_to_h(j2)
1144  ELSE
1145  jj_v_to_h(j2) = jj_v_to_h(j1)
1146  END IF
1147  END DO
1148  END DO
1149  END IF
1150  END IF
1151  IF (if_momentum) THEN
1152  ALLOCATE(h_to_ns(vv_mesh%np, 6, m_max_c))
1153  ALLOCATE(b_to_ns(vv_mesh%np, 6, m_max_c))
1154  ALLOCATE(hext(h_mesh%np, 6, m_max_c))
1155  ALLOCATE(bext(h_mesh%np, 6, m_max_c))
1156  ELSE
1157  ALLOCATE(h_to_ns(1, 1, 1))
1158  ALLOCATE(b_to_ns(1, 1, 1))
1159  END IF
1160 
1161  !===Create data structure jj_v_to_temp==========================================
1162  IF (if_energy) THEN
1163  ALLOCATE(v_to_energy(temp_mesh%np, 6, m_max_c))
1164  ALLOCATE(jj_v_to_temp(temp_mesh%np))
1165  jj_v_to_temp = -1
1166  IF (vv_mesh%me/=0) THEN
1167  DO m = 1, vv_mesh%me
1168  jj_v_to_temp(temp_mesh%jj(:,m)) = vv_mesh%jj(1:temp_mesh%gauss%n_w,m)
1169  END DO
1170  ALLOCATE(t_to_ns(vv_mesh%np, 2, m_max_c))
1171  ELSE
1172  ALLOCATE(t_to_ns(1, 1, 1))
1173  END IF
1174  END IF
1175 
1176  !===Create data structure jj_temp_to_H==========================================
1177  IF ((if_energy).AND.(if_induction)) THEN
1178  ALLOCATE(t_to_max(h_mesh%np, 2, m_max_c))
1179  ALLOCATE(h_to_energy(temp_mesh%np, 6, m_max_c))
1180  ALLOCATE(jj_temp_to_h(h_mesh%np))
1181  jj_temp_to_h = -1
1182  DO m = 1, temp_mesh%me
1183  jj_temp_to_h(h_mesh%jj(:,m)) = temp_mesh%jj(1:h_mesh%gauss%n_w,m)
1184  END DO
1185  DO ms = 1, interface_h_mu%mes
1186  DO ns = 1, h_mesh%gauss%n_ws
1187  j1 = interface_h_mu%jjs1(ns,ms)
1188  j2 = interface_h_mu%jjs2(ns,ms)
1189  IF (jj_temp_to_h(j1) == -1) THEN
1190  jj_temp_to_h(j1) = jj_temp_to_h(j2)
1191  ELSE
1192  jj_temp_to_h(j2) = jj_temp_to_h(j1)
1193  END IF
1194  END DO
1195  END DO
1196  END IF
1197 
1198  !===Initialize Navier-Stokes====================================================
1199  time_u = 0.d0
1200  IF (if_momentum .OR. inputs%type_pb=='mxx') THEN
1201  IF (vv_mesh%me/=0) THEN
1202  IF (inputs%irestart_u) THEN
1203  IF (.NOT. if_mass) THEN
1204  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un, un_m1, pn, pn_m1, &
1205  incpn, incpn_m1, inputs%file_name)
1206  ELSE
1207  CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un, un_m1, pn, pn_m1, &
1208  incpn, incpn_m1, inputs%file_name, opt_level_set=level_set, &
1209  opt_level_set_m1=level_set_m1, opt_max_vel=max_vel)
1210  END IF
1211  ELSE
1212  CALL init_velocity_pressure(vv_mesh, pp_mesh, time_u, &
1213  inputs%dt, list_mode, un_m1, un, pn_m1, pn, incpn_m1, incpn)
1214  IF (if_mass) THEN
1215  IF (inputs%if_level_set_P2) THEN
1216  CALL init_level_set(vv_mesh, time_u, &
1217  inputs%dt, list_mode, level_set_m1, level_set)
1218  ELSE
1219  CALL init_level_set(pp_mesh, time_u, &
1220  inputs%dt, list_mode, level_set_m1, level_set)
1221  END IF
1222 
1223  bloc_size = SIZE(un,1)/inputs%ndim(2) + 1
1224  m_max_pad = 3*SIZE(list_mode)*inputs%ndim(2)/2
1225  CALL fft_max_vel_dcl(comm_one_d_ns(2), 2*un-un_m1, max_vel_s, inputs%ndim(2), bloc_size, m_max_pad)
1226  CALL mpi_allreduce(max_vel_s, max_vel, 1, mpi_double_precision, &
1227  mpi_max, comm_one_d_ns(1), code)
1228  max_vel = max(1.1d0*max_vel, 0.1d0)
1229  END IF
1230  END IF
1231 
1232  IF (if_mass) THEN
1233  CALL reconstruct_variable(comm_one_d_ns, list_mode, pp_mesh, vv_mesh, level_set_m1, &
1234  inputs%density_fluid, density_m1)
1235  CALL reconstruct_variable(comm_one_d_ns, list_mode, pp_mesh, vv_mesh, level_set, &
1236  inputs%density_fluid, density)
1237  visc_entro_level = 0.d0
1238  END IF
1239 
1240  !===Force sine coefficients of zero mode to zero==========================
1241  DO i = 1, m_max_c
1242  IF (list_mode(i) == 0) THEN
1243  DO k= 1, 3
1244  un(:,2*k,i) = 0.d0
1245  un_m1(:,2*k,i) = 0.d0
1246  END DO
1247  pn(:,2,i) = 0.d0
1248  incpn(:,2,i) = 0.d0
1249  density(:,2,i) = 0.d0
1250  pn_m1(:,2,i) = 0.d0
1251  incpn_m1(:,2,i) = 0.d0
1252  density_m1(:,2,i) = 0.d0
1253  IF (if_mass) THEN
1254  level_set(:,:,2,i) = 0.d0
1255  level_set_m1(:,:,2,i) = 0.d0
1256  END IF
1257  END IF
1258  END DO
1259  !===End force sine coefficients of zero mode to zero======================
1260  density_m2 = density_m1
1261 
1262  END IF
1263  END IF
1264 
1265  !===Initialize velocity (time-independent) if mxw===============================
1266  IF ( (inputs%type_pb=='mxw') .AND. (h_mesh%me/=0) ) THEN
1267  DO i = 1, m_max_c !Initialization of vel
1268  v_to_max(:,:,i) = vexact(list_mode(i), h_mesh)
1269  END DO
1270  ENDIF
1271 
1272  !===Initialize velocity using un (time-independent) if mxx======================
1273  IF (inputs%type_pb=='mxx') THEN
1274  !===Use extension_velocity===================================================
1275  IF (h_mesh%np>0) THEN !We extend v_to_Max
1276  DO i = 1, m_max_c
1277  DO k= 1, 6 !The user has to code extension_vel
1278  v_to_max(:,k,i) = extension_velocity(k, h_mesh, list_mode(i), time_u, 1)
1279  END DO
1280  END DO
1281  END IF
1282 
1283  !===Use extension_velocity===================================================
1284  IF (vv_mesh%me /=0) THEN
1285  DO j = 1, SIZE(jj_v_to_h)
1286  IF (jj_v_to_h(j) == -1) cycle
1287  v_to_max(j,:,:) = un(jj_v_to_h(j),:,:)
1288  END DO
1289  END IF
1290 
1291  !===Cleanup==================================================================
1292  DEALLOCATE(un_m1)
1293  DEALLOCATE(un)
1294  DEALLOCATE(pn_m1)
1295  DEALLOCATE(pn)
1296  DEALLOCATE(incpn_m1)
1297  DEALLOCATE(incpn)
1298  DEALLOCATE(density_m2)
1299  DEALLOCATE(density_m1)
1300  DEALLOCATE(density)
1301  IF (if_mass) THEN
1302  IF (ALLOCATED(level_set)) DEALLOCATE(level_set,level_set_m1)
1303  END IF
1304  ENDIF
1305 
1306  !===Initialize Maxwell==========================================================
1307  time_h = 0.d0
1308  IF (if_induction) THEN
1309  IF (inputs%irestart_h) THEN
1310  CALL read_restart_maxwell(comm_one_d, h_mesh, phi_mesh, time_h, list_mode, hn, hn1, bn, bn1, &
1311  phin, phin1, inputs%file_name)
1312  ELSE
1313  CALL init_maxwell(h_mesh,phi_mesh,time_h,inputs%dt,mu_h_field,inputs%mu_phi,list_mode,&
1314  hn1,hn,phin1,phin)
1315  !===Initialize Bn and Bn1
1316  IF (h_mesh%me/=0) THEN
1317  IF (inputs%if_permeability_variable_in_theta) THEN
1318  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
1319  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
1320  bloc_size = SIZE(bn,1)/nb_procs+1
1321  CALL fft_par_var_eta_prod_t_dcl(comm_one_d(2), mu_in_real_space, &
1322  h_mesh, hn, bn, nb_procs, bloc_size, m_max_pad, time)
1323  CALL fft_par_var_eta_prod_t_dcl(comm_one_d(2), mu_in_real_space, &
1324  h_mesh, hn1, bn1, nb_procs, bloc_size, m_max_pad, time)
1325  ELSE
1326  DO i = 1, m_max_c
1327  DO k = 1, 6
1328  bn(:,k,i) = mu_h_field*hn(:,k,i)
1329  bn1(:,k,i) = mu_h_field*hn1(:,k,i)
1330  END DO
1331  END DO
1332  END IF
1333  END IF
1334  END IF
1335  !===Force sine coefficients of zero mode to zero=============================
1336  DO i = 1, m_max_c
1337  IF (list_mode(i) == 0) THEN
1338  IF (h_mesh%me/=0) THEN
1339  DO k = 1, 3
1340  hn(:,2*k,i) = 0.d0
1341  hn1(:,2*k,i) = 0.d0
1342  END DO
1343  END IF
1344  IF (phi_mesh%me/=0) THEN
1345  phin(:,2,i) = 0.d0
1346  phin1(:,2,i) = 0.d0
1347  END IF
1348  END IF
1349  END DO
1350 
1351  END IF
1352 
1353  !===Initialize temperature======================================================
1354  time_t = 0.d0
1355  IF (if_energy) THEN
1356  IF (temp_mesh%me/=0) THEN
1357  IF (inputs%irestart_T) THEN
1358  CALL read_restart_temp(comm_one_d_temp, time_t, list_mode, tempn, tempn_m1, &
1359  inputs%file_name)
1360  ELSE
1361  CALL init_temperature(temp_mesh, time_t, inputs%dt, list_mode, tempn_m1, tempn)
1362  END IF
1363  !===Force sine coefficients of zero mode to zero==========================
1364  DO i = 1, m_max_c
1365  IF (list_mode(i) == 0) THEN
1366  tempn(:,2,i) = 0.d0
1367  tempn_m1(:,2,i) = 0.d0
1368  END IF
1369  END DO
1370 
1371  END IF
1372  END IF
1373 
1374  !===Initialize time=============================================================
1375  IF (inputs%irestart_h .OR. inputs%irestart_u .OR. inputs%irestart_T) THEN
1376  time = max(time_u,time_h,time_t)
1377  ELSE
1378  time = 0.d0
1379  END IF
1380 
1381  !===Guardrail===================================================================
1382  IF (if_mass.AND.inputs%variation_sigma_fluid) THEN
1383  IF (inputs%analytical_permeability.OR.inputs%nb_inter_mu>0.OR.inputs%nb_dom_phi>0) THEN
1384  CALL error_petsc(' BUG in INIT : sigma reconstruct via level set not implemented with vacuum or variable mu')
1385  END IF
1386  END IF
1387 
1388  END SUBROUTINE init
1389 
1390  SUBROUTINE prodmat_maxwell_int_by_parts(vect_in, vect_out, ndim, i)
1391  USE update_maxwell
1392  IMPLICIT NONE
1393  INTEGER :: ndim
1394  REAL(KIND=8), DIMENSION(ndim) :: vect_in
1395  REAL(KIND=8), DIMENSION(ndim) :: vect_out
1396  INTEGER :: i
1397  INTEGER :: type, i_deb, i_fin
1398  REAL(KIND=8), DIMENSION(vv_mesh%np,2,SIZE(list_mode)) :: sigma_ns
1399 
1400  time = 0.d0
1401  DO TYPE = 1, 6
1402  i_deb = (type-1)*h_mesh%np+1
1403  i_fin = i_deb + h_mesh%np -1
1404  IF (modulo(type,2)==0 .AND. list_mode(i)==0) THEN
1405  hn(:,type,i) = 0.d0
1406  ELSE
1407  hn(:,type,i) = vect_in(i_deb:i_fin)
1408  END IF
1409  END DO
1410  DO TYPE = 1, 2
1411  phin(:,type,i) =0.d0
1412  END DO
1413 
1414  DO TYPE = 1, 6
1415  i_deb = 6*h_mesh%np + (type-1)*h_mesh%np+1
1416  i_fin = i_deb + h_mesh%np -1
1417  IF (modulo(type,2)==0 .AND. list_mode(i)==0) THEN
1418  hn1(:,type,i) = 0.d0
1419  ELSE
1420  hn1(:,type,i) = vect_in(i_deb:i_fin)
1421  END IF
1422  END DO
1423  DO TYPE = 1, 2
1424  phin1(:,type,i) =0.d0
1425  END DO
1426 
1427  CALL maxwell_decouple(comm_one_d, h_mesh, pmag_mesh, phi_mesh, interface_h_phi, interface_h_mu, &
1428  hn, bn, phin, hn1, bn1, phin1, v_to_max, inputs%stab, sigma_field, r_fourier, index_fourier, mu_h_field, inputs%mu_phi, &
1429  time, inputs%dt, inputs%Rem, list_mode, h_phi_per, la_h, la_pmag, la_phi, la_mhd, sigma_ns, jj_v_to_h)
1430 
1431  DO TYPE = 1, 6
1432  i_deb = (type-1)*h_mesh%np+1
1433  i_fin = i_deb + h_mesh%np -1
1434  vect_out(i_deb:i_fin) = hn(:,type,i)
1435  END DO
1436 
1437  DO TYPE = 1, 6
1438  i_deb = 6*h_mesh%np + (type-1)*h_mesh%np+1
1439  i_fin = i_deb + h_mesh%np -1
1440  vect_out(i_deb:i_fin) = hn1(:,type,i)
1441  END DO
1442 
1443  END SUBROUTINE prodmat_maxwell_int_by_parts
1444 
1445 
1447  IMPLICIT NONE
1448  INTEGER :: narg, i
1449  CHARACTER(len=200),DIMENSION(:,:), ALLOCATABLE :: inline
1450  LOGICAL :: ok
1451  CHARACTER(len=3) :: tit
1452 
1453  narg = 0
1454  ok = .true.
1455 
1456  DO WHILE (ok)
1457  CALL getarg(narg+1,tit)
1458  IF (tit == ' ') THEN
1459  ok = .false.
1460  ELSE
1461  narg = narg+1
1462  END IF
1463  END DO
1464 
1465  narg = narg/2
1466  ALLOCATE(inline(2,narg))
1467 
1468  DO i = 1, narg
1469  CALL getarg(2*(i-1)+1,inline(1,i))
1470  CALL getarg(2*i ,inline(2,i))
1471  END DO
1472 
1473  inputs%test_de_convergence = .false.
1474  inputs%numero_du_test_debug = 0
1475  inputs%data_directory_debug = '.'
1476  DO i = 1, narg
1477  IF (trim(adjustl(inline(1,i))) == 'debug') THEN
1478  inputs%test_de_convergence = .true.
1479  READ(inline(2,i),*) inputs%numero_du_test_debug
1480  ELSE IF (trim(adjustl(inline(1,i))) == 'debug_dir') THEN
1481  inputs%data_directory_debug=inline(2,i)
1482  END IF
1483  END DO
1484 
1485  END SUBROUTINE sfemansinitialize
1486 
1487  SUBROUTINE compute_local_mesh_size(mesh)
1488  USE def_type_mesh
1489  TYPE(mesh_type) :: mesh
1490  INTEGER :: m, type_fe, index, l
1491  IF (mesh%gauss%n_w==3) THEN
1492  type_fe = 1
1493  ELSE
1494  type_fe = 2
1495  END IF
1496  ALLOCATE(mesh%hloc(mesh%dom_me))
1497  ALLOCATE(mesh%hloc_gauss(mesh%gauss%l_G*mesh%dom_me))
1498 
1499  index = 0
1500  DO m = 1, mesh%dom_me
1501  mesh%hloc(m) = sqrt(sum(mesh%gauss%rj(:,m)))/type_fe
1502  DO l = 1, mesh%gauss%l_G
1503  index = index + 1
1504  mesh%hloc_gauss(index) = mesh%hloc(m)
1505  END DO
1506  END DO
1507  END SUBROUTINE compute_local_mesh_size
1508 
1510  REAL(KIND=8) :: h_min, h_min_f
1511  INTEGER :: code
1512  !===Compute h_min
1513  IF (inputs%if_level_set_P2) THEN
1514  h_min_f=minval(vv_mesh%hloc_gauss)
1515  ELSE
1516  h_min_f=minval(pp_mesh%hloc_gauss)
1517  END IF
1518  CALL mpi_allreduce(h_min_f, h_min, 1, mpi_double_precision, &
1519  mpi_min, comm_one_d_ns(1), code)
1520  inputs%h_min_distance = inputs%multiplier_for_h_min_distance*h_min
1521  !===End compute h_min
1522  END SUBROUTINE compute_local_mesh_size_level_set
1523 
1524 END MODULE initialization
subroutine, public save_run(it, freq_restart)
Definition: tn_axi.f90:5
subroutine, public reconstruct_variable(comm_one_d, list_mode, mesh_P1, mesh_P2, level_set, values, variable)
Definition: sub_mass.f90:70
subroutine prepare_zero_out_modes(list_mode, list_mode_to_zero_out, select_mode)
subroutine zero_out_modes
real(kind=8) function, dimension(h_mesh%np), public extension_velocity(TYPE, H_mesh, mode, t, n_start)
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
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)
Definition: restart.f90:8
subroutine, public free_interface(interf)
subroutine, public free_mesh(mesh)
subroutine, public init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, list_mode, Hn1, Hn, phin1, phin)
subroutine projection_velocity(mesh, vn, connectivity_structure, coupling_variable)
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 compute_local_mesh_size_level_set
subroutine write_restart_temp(communicator, temp_mesh, time, list_mode, tempn, tempn_m1, filename, it, freq_restart, opt_mono)
Definition: restart.f90:534
subroutine, public load_interface(mesh_master, mesh_slave, list_inter, mesh_INTERFACE, disjoint)
subroutine, public load_dg_mesh_free_format(dir, fil, list_dom, list_inter, type_fe, mesh, mesh_formatted)
Definition: prep_mesh.f90:16
subroutine, public fft_max_vel_dcl(communicator, V1_in, V_out, nb_procs, bloc_size, m_max_pad)
subroutine, public navier_stokes_decouple(comm_one_d_ns, time, vv_3_LA, pp_1_LA, list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, vvz_per, pp_per, Hn_p2, Bn_p2, density_m2, density_m1, density, visco_dyn, tempn, level_set_m1, level_set, visc_entro_level)
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:611
subroutine, public st_scr_maxwell_mu_h_p_phi(communicator, H_mesh_glob, H_mesh, pmag_mesh_glob, pmag_mesh, phi_mesh_glob, phi_mesh, interface_glob, interface_H_mu_glob, LA_H, LA_pmag, LA_phi, LA_mhd, opt_per)
Definition: st_csr_mhd.f90:9
subroutine sfemansinitialize
subroutine, public fft_par_var_eta_prod_t_dcl(communicator, eta, H_mesh, c_in, c_out, nb_procs, bloc_size, m_max_pad, time, temps)
subroutine, public initial(vv_mesh_out, pp_mesh_out, H_mesh_out, phi_mesh_out, temp_mesh_out, interface_H_phi_out, interface_H_mu_out, list_mode_out, un_out, pn_out, Hn_out, Bn_out, phin_out, v_to_Max_out, vol_heat_capacity_field_out, temperature_diffusivity_field_out, mu_H_field_out, sigma_field_out, time_out, m_max_c_out, comm_one_d_out, comm_one_d_ns_out, comm_one_d_temp_out, tempn_out, level_set_out, density_out)
subroutine projection_temperature(mesh, vn, connectivity_structure, if_restriction, coupling_variable)
real(kind=8) function, dimension(h_mesh%np, 6), public vexact(m, H_mesh)
subroutine, public read_my_data(data_fichier)
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)
Definition: restart.f90:99
subroutine assign_boundary
Definition: boundary.f90:40
subroutine, public st_aij_csr_glob_block(communicator, kmax, mesh_glob, mesh, LA, opt_per)
Definition: st_csr.f90:163
subroutine, public prep_periodic_bloc(my_periodic, mesh, periodic, nb_bloc)
real(kind=8) function, dimension(nb_angles, ne-nb+1), public mu_in_real_space(H_mesh, angles, nb_angles, nb, ne, time)
subroutine, public prodmat_maxwell_int_by_parts(vect_in, vect_out, ndim, i)
subroutine, public create_cart_comm(ndim, comm_cart, comm_one_d, coord_cart)
subroutine, public three_level_temperature(comm_one_d, time, temp_1_LA, dt, list_mode, temp_mesh, tempn_m1, tempn, chmp_vit, vol_heat_capacity, temp_diffusivity, my_par_cc, temp_list_dirichlet_sides, temp_per)
subroutine, public sfemans_initialize_postprocessing(comm_one_d, vv_mesh, pp_mesh, H_mesh, phi_mesh, temp_mesh, list_mode_in, opt_nb_plane)
subroutine, public symmetric_points(mesh_loc, mesh_glob, ltg_LA)
Definition: symmetry.f90:18
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
subroutine compute_local_mesh_size(mesh)
subroutine, public gen_gauss(mesh, edge_stab)
Definition: gen_gauss.f90:6
real(kind=8) function, dimension(ne-nb+1), public mu_bar_in_fourier_space(H_mesh, nb, ne, pts, pts_ids)
subroutine, public prep_periodic_h_p_phi_bc(my_periodic, H_mesh, pmag_mesh, phi_mesh, H_p_phi_per)
subroutine error_petsc(string)
Definition: my_util.f90:15
subroutine projection_mag_field(mesh, vn, connectivity_structure, coupling_variable)
subroutine, public run_sfemans(time_in)
subroutine read_restart_maxwell(communicator, H_mesh, phi_mesh, time, list_mode, Hn, Hn1, Bn, Bn1, phin, phin1, filename, val_init, interpol, opt_mono)
Definition: restart.f90:374
subroutine, public three_level_mass(comm_one_d, time, level_set_LA_P1, level_set_LA_P2, list_mode, mesh_P1, mesh_P2, chmp_vit_P2, max_vel, level_set_per, density_m2, density_m1, density, level_set_m1, level_set, visc_entro_level)
Definition: sub_mass.f90:9
subroutine, public init_level_set(vv_mesh, time, dt, list_mode, level_set_m1, level_set)
subroutine, public init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
subroutine, public prep_periodic(my_periodic, mesh, periodic)
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:193
subroutine, public maxwell_decouple(comm_one_d, H_mesh, pmag_mesh, phi_mesh, interface_H_phi, interface_H_mu, Hn, Bn, phin, Hn1, Bn1, phin1, vel, stab_in, sigma_in, R_fourier, index_fourier, mu_H_field, mu_phi, time, dt, Rem, list_mode, H_phi_per, LA_H, LA_pmag, LA_phi, LA_mhd, sigma_ns, jj_v_to_H)
Definition: sub_maxwell.f90:7