SFEMaNS  version 5.3
Reference documentation for SFEMaNS
read_my_data.f90
Go to the documentation of this file.
1 MODULE my_data_module
2  USE def_type_mesh
3  USE solve_petsc
4  TYPE my_data
5  LOGICAL :: iformatted
6  CHARACTER(LEN=200) :: directory, file_name
7  LOGICAL :: if_read_partition
8  CHARACTER(LEN=3) :: type_pb
9  INTEGER, DIMENSION(2) :: ndim
10  LOGICAL :: select_mode
11  INTEGER :: m_max
12  INTEGER, DIMENSION(:), POINTER :: list_mode_lect
13  !===Data time iterations========================================================
14  INTEGER :: nb_iteration
15  REAL(KIND=8) :: dt
16  !===Data for LES================================================================
17  LOGICAL :: les
18  REAL(KIND=8) :: les_coeff1, les_coeff2, les_coeff3, les_coeff4
19  REAL(KIND=8) :: les_coeff1_mom
20  LOGICAL :: if_les_in_momentum
21  !===Data for precession=========================================================
22  REAL(KIND=8) :: taux_precession, angle_s_pi
23  LOGICAL :: precession
24  !===Data for NS penalization====================================================
25  LOGICAL :: if_ns_penalty
26  LOGICAL :: if_impose_vel_in_solids
27  LOGICAL :: if_compute_momentum_pseudo_force
28  REAL(KIND=8) :: div_stab_in_ns
29  !===Data for Navier-Stokes======================================================
30  LOGICAL :: if_navier_stokes_with_u
31  LOGICAL :: if_navier_stokes_art_comp
32  REAL(KIND=8) :: penal_coeff_art_comp
33  LOGICAL :: if_tensor_sym
34  LOGICAL :: if_moment_bdf2
35  LOGICAL :: if_level_bdf2
36  LOGICAL :: irestart_u
37  LOGICAL :: if_variable_visco
38  LOGICAL :: if_navier_stokes_with_taylor !===JLG july 20, 2019, p3 mesh
39  INTEGER :: taylor_order !===JLG Dec 2020
40  REAL(KIND=8) :: taylor_lambda!===JLG Dec 2020
41  INTEGER :: type_fe_velocity !===JLG july 20, 2019, p3 mesh
42  INTEGER :: nb_dom_ns
43  INTEGER, DIMENSION(:), POINTER :: list_dom_ns
44  REAL(KIND=8) :: re
45  TYPE(solver_param) :: my_par_vv, my_par_pp, my_par_mass
46  INTEGER :: pp_nb_dirichlet_sides
47  INTEGER, DIMENSION(:), POINTER :: pp_list_dirichlet_sides
48  INTEGER, DIMENSION(3) :: vv_nb_dirichlet_sides
49  TYPE(dyn_int_line), DIMENSION(3) :: vv_list_dirichlet_sides
50  INTEGER :: vv_nb_dirichlet
51  INTEGER :: vv_nb_dirichlet_normal_velocity
52  INTEGER, DIMENSION(:), POINTER :: vv_list_dirichlet_normal_velocity_sides
53  !===Data for Maxwell============================================================
54  LOGICAL :: if_maxwell_with_h
55  LOGICAL :: irestart_h
56  INTEGER :: nb_dom_h, nb_dom_phi, nb_inter, nb_inter_mu
57  INTEGER :: type_fe_h, type_fe_phi, nb_dirichlet_sides_h
58  INTEGER, DIMENSION(:), POINTER :: list_dom_h, list_dom_phi, list_inter_h_phi
59  INTEGER, DIMENSION(:), POINTER :: list_inter_mu, list_dirichlet_sides_h
60  REAL(KIND=8) :: rem, mu_phi
61  LOGICAL :: analytical_permeability
62  LOGICAL :: if_use_fem_integration_for_mu_bar
63  LOGICAL :: if_permeability_variable_in_theta
64  REAL(KIND=8), DIMENSION(:), POINTER :: mu_h, sigma
65  REAL(KIND=8), DIMENSION(3) :: stab
66  TYPE(solver_param) :: my_par_h_p_phi
67  INTEGER :: phi_nb_dirichlet_sides
68  INTEGER, DIMENSION(:), POINTER :: phi_list_dirichlet_sides
69  LOGICAL :: if_quasi_static_approx
70  LOGICAL :: if_steady_current_fhd ! MODIFICATION: for fhd, possibility of computing H only once in the case of steady current
71  !===Data for temperature========================================================
72  LOGICAL :: if_temperature
73  LOGICAL :: irestart_t
74  LOGICAL :: if_helmholtz_force
75  REAL(KIND=8) :: gravity_coefficient
76  REAL(KIND=8) :: mag_force_coefficient ! MODIFICATION: coefficient for magnetic force in ferrofluids
77  REAL(KIND=8), DIMENSION(:), POINTER :: temperature_diffusivity
78  REAL(KIND=8), DIMENSION(:), POINTER :: density, heat_capacity, vol_heat_capacity ! MODIFICATION: density and heat capacity found in the litterature more often than vol heat capacity
79  TYPE(solver_param) :: my_par_temperature
80  INTEGER :: temperature_nb_dirichlet_sides
81  INTEGER, DIMENSION(:), POINTER :: temperature_list_dirichlet_sides
82  INTEGER :: temperature_nb_robin_sides ! MODIFICATION: Robin
83  INTEGER, DIMENSION(:), POINTER :: temperature_list_robin_sides
84  REAL(KIND=8), DIMENSION(:), POINTER :: convection_coeff ! MODIFICATION: Robin
85  REAL(KIND=8), DIMENSION(:), POINTER :: exterior_temperature ! MODIFICATION: Robin
86  INTEGER :: nb_dom_temp
87  INTEGER, DIMENSION(:), POINTER :: list_dom_temp
88  INTEGER :: nb_inter_v_t
89  INTEGER, DIMENSION(:), POINTER :: list_inter_v_t
90  !===Data for level set==========================================================
91  LOGICAL :: if_level_set
92  LOGICAL :: if_level_set_fixed
93  TYPE(solver_param) :: my_par_level_set
94  INTEGER :: level_set_nb_dirichlet_sides
95  INTEGER, DIMENSION(:), POINTER :: level_set_list_dirichlet_sides
96  REAL(KIND=8) :: level_set_cmax
97  REAL(KIND=8) :: level_set_comp_factor
98  CHARACTER(LEN=3) :: level_set_reconstruction_type
99  REAL(KIND=8) :: level_set_tanh_coeff_reconstruction
100  INTEGER :: nb_fluid
101  REAL(KIND=8), DIMENSION(:), POINTER :: density_fluid
102  REAL(KIND=8), DIMENSION(:), POINTER :: dyna_visc_fluid
103  LOGICAL :: variation_sigma_fluid
104  REAL(KIND=8), DIMENSION(:), POINTER :: sigma_fluid
105  LOGICAL :: if_surface_tension
106  REAL(KIND=8), DIMENSION(:), POINTER :: coeff_surface
107  LOGICAL :: if_mass_correction
108  LOGICAL :: if_kill_overshoot
109  REAL(KIND=8) :: multiplier_for_h_min_distance
110  REAL(KIND=8) :: sigma_min, mu_min !LC 2017/01/27
111  !===Computed data for level set
112  REAL(KIND=8) :: h_min_distance
113  LOGICAL :: if_level_set_p2
114  LOGICAL :: if_compression_mthd_jlg
115  !===Data for periodicity========================================================
116  TYPE(periodic_data) :: my_periodic
117  !===Data for convergence tests==================================================
118  LOGICAL :: test_de_convergence
119  CHARACTER(len=200) :: data_directory_debug
120  INTEGER :: numero_du_test_debug
121  REAL(KIND=8), DIMENSION(4) :: norm_ref
122  !===Data for Arpack=============================================================
123  LOGICAL :: if_arpack
124  CHARACTER(len=2) :: arpack_which
125  INTEGER :: nb_eigenvalues, arpack_iter_max
126  REAL(KIND=8) :: arpack_tolerance
127  LOGICAL :: if_arpack_vtu_2d
128  !===Data for stress bc==========================================================
129  REAL(KIND=8) :: stab_bdy_ns
130  !===Data for postprocessing=====================================================
131  INTEGER :: number_of_planes_in_real_space
132  LOGICAL :: check_numerical_stability
133  LOGICAL :: is_mesh_symmetric
134  LOGICAL :: if_zero_out_modes
135  INTEGER :: nb_select_mode_ns, nb_select_mode_mxw
136  INTEGER, DIMENSION(:), POINTER :: list_select_mode_ns, list_select_mode_mxw
137  INTEGER :: freq_restart, freq_en, freq_plot
138  LOGICAL :: verbose_timing, verbose_divergence, verbose_cfl
139  LOGICAL :: if_just_processing
140  LOGICAL :: if_xml
141  CONTAINS
142  PROCEDURE, PUBLIC :: init
143  END TYPE my_data
144 CONTAINS
145  SUBROUTINE init(a)
146  CLASS(my_data), INTENT(INOUT) :: a
147  !===Logicals
148  a%iformatted=.false.
149  a%if_read_partition=.false.
150  a%select_mode=.false.
151  a%LES=.false.
152  a%if_LES_in_momentum=.false.
153  a%precession=.false.
154  a%if_ns_penalty=.false.
155  a%if_impose_vel_in_solids=.false.
156  a%if_compute_momentum_pseudo_force=.false.
157  a%if_navier_stokes_with_u=.false.
158  a%if_navier_stokes_art_comp=.false.
159  a%if_tensor_sym=.false.
160  a%if_moment_bdf2=.false.
161  a%if_level_bdf2=.false.
162  a%irestart_u=.false.
163  a%if_maxwell_with_H=.false.
164  a%irestart_h=.false.
165  a%irestart_T=.false.
166  a%if_helmholtz_force=.false.
167  a%analytical_permeability=.false.
168  a%if_use_fem_integration_for_mu_bar=.false.
169  a%if_permeability_variable_in_theta=.false.
170  a%if_quasi_static_approx=.false.
171  a%if_steady_current_fhd=.false. ! MODIFICATION
172  a%if_temperature=.false.
173  a%if_level_set=.false.
174  a%if_level_set_fixed=.false.
175  a%variation_sigma_fluid=.false.
176  a%if_surface_tension=.false.
177  a%if_mass_correction=.false.
178  a%if_kill_overshoot=.false.
179  a%if_level_set_P2=.false.
180  a%if_compression_mthd_JLG=.false.
181  a%if_variable_visco=.false.
182  a%if_xml = .true.
183  a%if_navier_stokes_with_taylor =.false.
184  !===Done in sfemansinitialize. Do not touch. a%test_de_convergence
185  a%if_arpack=.false.
186  a%if_arpack_vtu_2d=.false.
187  a%check_numerical_stability=.false.
188  a%is_mesh_symmetric=.false.
189  a%if_zero_out_modes=.false.
190  a%verbose_timing =.false.
191  a%verbose_divergence=.false.
192  a%verbose_CFL=.false.
193  a%if_just_processing=.false.
194  a%my_par_vv%verbose=.false.
195  a%my_par_pp%verbose=.false.
196  a%my_par_mass%verbose=.false.
197  a%my_par_H_p_phi%verbose=.false.
198  a%my_par_temperature%verbose=.false.
199  a%my_par_level_set%verbose=.false.
200  !===Reals
201  a%LES_coeff1=0.d0
202  a%LES_coeff2=0.d0
203  a%LES_coeff3=0.d0
204  a%LES_coeff4=0.d0
205  a%LES_coeff1_mom=0.d0
206  a%taux_precession=0.d0
207  a%angle_s_pi=0.d0
208  a%div_stab_in_ns=0.d0
209  a%stab=0.d0
210  a%gravity_coefficient=0.d0
211  a%mag_force_coefficient=0.d0
212  a%level_set_cmax=0.d0
213  a%level_set_comp_factor=0.d0
214  a%level_set_tanh_coeff_reconstruction=0.d0
215  a%norm_ref=0.d0
216  a%arpack_tolerance=0.d0
217  a%stab_bdy_ns=0.d0
218  a%h_min_distance=0.d0
219  a%multiplier_for_h_min_distance=0.d0
220  a%taylor_lambda = 1.d0
221  a%penal_coeff_art_comp=1.d0
222  !===Integers
223  a%vv_nb_dirichlet=0
224  a%vv_nb_dirichlet_normal_velocity=0
225  a%freq_restart = 10000000
226  a%freq_en = 10000000
227  a%freq_plot = 10000000
228  a%type_fe_velocity = 2
229  a%taylor_order = -1
230  END SUBROUTINE init
231 END MODULE my_data_module
232 
233 MODULE input_data
234  USE my_data_module
235  IMPLICIT NONE
236  PUBLIC :: read_my_data
237  TYPE(my_data), PUBLIC :: inputs
238  PRIVATE
239 
240 CONTAINS
241 
242  SUBROUTINE read_my_data(data_fichier)
243  USE chaine_caractere
244  USE my_util
245  IMPLICIT NONE
246  LOGICAL :: test, test_sigma_min
247  CHARACTER(len=200), INTENT(IN) :: data_fichier
248  INTEGER :: k
249  CHARACTER(len=8), DIMENSION(3) :: vel_component=(/'uradial','utheta ','uzaxis '/)
250 
251  !===Initialize data to zero and false by default===============================
252  CALL inputs%init
253 
254  !===Open data file=============================================================
255  OPEN(unit = 21, file = data_fichier, form = 'formatted', status = 'unknown')
256 
257  !===Location of mesh============================================================
258  CALL read_until(21, '===Is mesh file formatted (true/false)?')
259  READ(21,*) inputs%iformatted
260  CALL read_until(21, '===Directory and name of mesh file')
261  READ(21,*) inputs%directory, inputs%file_name
262 
263  !===Processor distribution======================================================
264  CALL read_until(21, '===Number of processors in meridian section')
265  READ (21,*) inputs%ndim(1)
266  CALL read_until(21, '===Number of processors in Fourier space')
267  READ (21,*) inputs%ndim(2)
268 
269  !===Fourier modes===============================================================
270  CALL read_until(21, '===Number of Fourier modes')
271  READ(21,*) inputs%m_max
272  IF (mod(inputs%m_max,inputs%ndim(2))/= 0) THEN
273  CALL error_petsc('BUG in read_my_data, MOD(nb_select_mode,nb_procs_F)/= 0')
274  END IF
275  CALL find_string(21,'===Select Fourier modes? (true/false)', test)
276  IF (test) THEN
277  READ (21,*) inputs%select_mode
278  ELSE
279  inputs%select_mode = .false.
280  END IF
281  IF (inputs%select_mode) THEN
282  ALLOCATE(inputs%list_mode_lect(inputs%m_max))
283  CALL read_until(21, '===List of Fourier modes (if select_mode=.TRUE.)')
284  READ(21,*) inputs%list_mode_lect
285  END IF
286 
287  !===Type of problem to be solved================================================
288  CALL read_until(21, '===Problem type: (nst, mxw, mhd, fhd)')
289  READ(21,*) inputs%type_pb
290  IF (inputs%type_pb/='nst' .AND. inputs%type_pb/='mhd' .AND. inputs%type_pb/='mxw' .AND. inputs%type_pb/='fhd') THEN
291  CALL error_petsc('BUG in read_my_data, type_pb of probleme not yet defined')
292  END IF
293 
294  !===Restarts====================================================================
295  CALL read_until(21, '===Restart on velocity (true/false)')
296  READ (21 ,*) inputs%irestart_u
297  IF (inputs%type_pb=='mhd' .OR. inputs%type_pb=='mxw' .OR. inputs%type_pb=='fhd') THEN
298  CALL read_until(21, '===Restart on magnetic field (true/false)')
299  READ (21 ,*) inputs%irestart_h
300  ELSE
301  inputs%irestart_h = .false.
302  END IF
303  CALL find_string(21, '===Restart on temperature (true/false)', test)
304  IF (test) THEN
305  READ (21 ,*) inputs%irestart_T
306  ELSE
307  inputs%irestart_T = .false.
308  END IF
309 
310  !===Mesh partitioning===========================================================
311  CALL find_string(21, '===Do we read metis partition? (true/false)', test)
312  IF (test) THEN
313  READ (21,*) inputs%if_read_partition
314  ELSE
315  inputs%if_read_partition = .false.
316  END IF
317  IF (.NOT.inputs%if_read_partition .AND. (inputs%irestart_u .OR. inputs%irestart_h .OR. inputs%irestart_T)) THEN
318  call error_petsc('Possibility of bug: set "===Do we read metis partition? (true/false) "' // &
319  'parameter to .true. when restarting a computation since number of procs, ' // &
320  'or machine, or type_pb may have changed. Make sure your mesh_part file is the correct one.')
321  END IF
322 
323  !===Time integration============================================================
324  CALL read_until(21, '===Time step and number of time iterations')
325  READ(21,*) inputs%dt, inputs%nb_iteration
326 
327  !===Check numerical stability===================================================
328  CALL find_string(21, '===Check numerical stability (true/false)', test)
329  IF (test) THEN
330  READ (21,*) inputs%check_numerical_stability
331  ELSE
332  inputs%check_numerical_stability = .false.
333  END IF
334 
335  !===Periodicity=================================================================
336  CALL find_string(21, '===How many pieces of periodic boundary?', test)
337  IF (test) THEN
338  READ (21,*) inputs%my_periodic%nb_periodic_pairs
339  ELSE
340  inputs%my_periodic%nb_periodic_pairs = 0
341  END IF
342  IF (inputs%my_periodic%nb_periodic_pairs.GE.1) THEN
343  ALLOCATE(inputs%my_periodic%list_periodic(2,inputs%my_periodic%nb_periodic_pairs))
344  ALLOCATE(inputs%my_periodic%vect_e(2,inputs%my_periodic%nb_periodic_pairs))
345  CALL read_until(21, '===Indices of periodic boundaries and corresponding vectors')
346  DO k = 1, inputs%my_periodic%nb_periodic_pairs
347  READ(21,*) inputs%my_periodic%list_periodic(:,k), inputs%my_periodic%vect_e(:,k)
348  END DO
349  END IF
350 
351  !===Mesh symmetry===============================================================
352  CALL find_string(21, '===Is the mesh symmetric (true/false)?', test)
353  IF (test) THEN
354  READ (21,*) inputs%is_mesh_symmetric
355  ELSE
356  inputs%is_mesh_symmetric = .false.
357  END IF
358 
359  !===Navier Stokes data==========================================================
360  IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd' .OR. inputs%irestart_u) THEN
361 
362  !==========Navier Stokes with artifical compression or projection-correction schemes===============!
363  CALL find_string(21, '===Solve Navier-Stokes with art comp scheme (true) or (false)?', test)
364  IF (test) THEN
365  READ(21,*) inputs%if_navier_stokes_art_comp
366  ELSE
367  inputs%if_navier_stokes_art_comp = .false.
368  END IF
369  IF (inputs%if_navier_stokes_art_comp) THEN
370  CALL find_string(21, '===Penalty coefficient for artifical compression', test)
371  IF (test) THEN
372  READ(21,*) inputs%penal_coeff_art_comp
373  ELSE
374  inputs%penal_coeff_art_comp=1.d0
375  END IF
376  END IF
377 
378  !==========Navier Stokes with u or m===============!
379  CALL find_string(21, '===Solve Navier-Stokes with u (true) or m (false)?', test)
380  IF (test) THEN
381  READ(21,*) inputs%if_navier_stokes_with_u
382  ELSE
383  inputs%if_navier_stokes_with_u = .true.
384  END IF
385  !==========Finite element type for velocity========!
386  CALL find_string(21, '===Type of finite element for velocity field', test)
387  IF (test) THEN
388  READ(21,*) inputs%type_fe_velocity
389  ELSE
390  inputs%type_fe_velocity=2
391  END IF
392 !!$ CALL find_string(21, '===Is tensor symmetric (true/false)?', test)
393 !!$ IF (test) THEN
394 !!$ READ(21,*) inputs%if_tensor_sym
395 !!$ ELSE
396 !!$ IF(inputs%if_navier_stokes_with_u) THEN
397 !!$ inputs%if_tensor_sym = .FALSE.
398 !!$ ELSE
399 !!$ inputs%if_tensor_sym = .TRUE.
400 !!$ END IF
401 !!$ END IF
402  !===Tensor symmetric always used.
403  inputs%if_tensor_sym = .true.
404 
405  IF (.NOT.inputs%if_navier_stokes_with_u) THEN
406  CALL find_string(21, '===Do we solve momentum with bdf2 (true/false)?', test)
407  IF (test) THEN
408  READ(21,*) inputs%if_moment_bdf2
409  ELSE
410  inputs%if_moment_bdf2=.false.
411  END IF
412  IF (inputs%if_moment_bdf2) THEN
413  inputs%if_level_bdf2=.true.
414  ELSE
415  CALL find_string(21, '===Do we solve level set with bdf2 (true/false)?', test)
416  IF (test) THEN
417  READ(21,*) inputs%if_level_bdf2
418  ELSE
419  inputs%if_level_bdf2=.false.
420  END IF
421  END IF
422  ELSE
423  !===JLG July 20, 2019, p3 mesh
424  CALL find_string(21, '===Do we use Taylor method to solve momentum (true/false)?', test)
425  IF (test) THEN
426  READ(21,*) inputs%if_navier_stokes_with_taylor
427  CALL read_until(21, '===Time accuracy of Taylor method (3 or 4)?')
428  READ(21,*) inputs%taylor_order
429  CALL read_until(21, '===Lambda parameter for Taylor method?')
430  READ(21,*) inputs%taylor_lambda
431  ELSE
432  inputs%if_navier_stokes_with_taylor = .false.
433  inputs%taylor_order = -1
434  END IF
435  !===JLG July 20, 2019, p3 mesh
436  END IF
437 
438  CALL read_until(21, '===Number of subdomains in Navier-Stokes mesh')
439  READ(21,*) inputs%nb_dom_ns
440  ALLOCATE(inputs%list_dom_ns(inputs%nb_dom_ns))
441  CALL read_until(21, '===List of subdomains for Navier-Stokes mesh')
442  READ(21,*) inputs%list_dom_ns
443 
444  !==========Dirichlet BCs for velocity==============!
445  DO k = 1, 3
446  CALL find_string(21, '===How many boundary pieces for Dirichlet BCs on '//trim(vel_component(k))//'?', test)
447  IF (test) THEN
448  CALL error_petsc( '===How many boundary pieces for Dirichlet BCs on '//trim(vel_component(k))//'? is disabled')
449  !READ(21,*) inputs%vv_nb_dirichlet_sides(k)
450  ELSE
451  inputs%vv_nb_dirichlet_sides(k) = 0
452  END IF
453  IF (inputs%vv_nb_dirichlet_sides(k)>0) THEN
454  ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(inputs%vv_nb_dirichlet_sides(k)))
455  CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on '//trim(vel_component(k)))
456  READ(21,*) inputs%vv_list_dirichlet_sides(k)%DIL
457  ELSE
458  ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(0))
459  END IF
460  END DO
461  CALL find_string(21, '===How many boundary pieces for full Dirichlet BCs on velocity?', test)
462  IF (test) THEN
463  READ(21,*) inputs%vv_nb_dirichlet
464  ELSE
465  inputs%vv_nb_dirichlet=0
466  END IF
467  IF (inputs%vv_nb_dirichlet>0) THEN
468  DO k = 1, 3
469  ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(inputs%vv_nb_dirichlet))
470  CALL read_until(21, '===List of boundary pieces for full Dirichlet BCs on velocity')
471  READ(21,*) inputs%vv_list_dirichlet_sides(k)%DIL
472  END DO
473  ELSE
474  DO k = 1, 3
475  ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(0))
476  END DO
477  END IF
478 
479  !===Homogeneous normal velocity====================!
480  CALL find_string(21, '===How many boundary pieces for homogeneous normal velocity?', test)
481  IF (test) THEN
482  READ(21,*) inputs%vv_nb_dirichlet_normal_velocity
483  ELSE
484  inputs%vv_nb_dirichlet_normal_velocity=0
485  END IF
486  IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
487  ALLOCATE(inputs%vv_list_dirichlet_normal_velocity_sides(inputs%vv_nb_dirichlet_normal_velocity))
488  CALL read_until(21, '===List of boundary pieces for homogeneous normal velocity')
489  READ(21,*) inputs%vv_list_dirichlet_normal_velocity_sides
490  CALL find_string(21, '===stab_bdy_ns', test)
491  IF (test) THEN
492  READ (21,*) inputs%stab_bdy_ns
493  ELSE
494  inputs%stab_bdy_ns=1.d0
495  END IF
496  DO k = 1, inputs%vv_nb_dirichlet_normal_velocity
497  IF (minval(abs(inputs%vv_list_dirichlet_normal_velocity_sides(k) &
498  - inputs%vv_list_dirichlet_sides(1)%DIL))==0) THEN
499  CALL error_petsc('Boundary conditions are mixed up')
500  END IF
501  END DO
502  END IF
503 
504  !===Stress boundary conditions=====================! !Disabled
505  CALL find_string(21, '===Stress boundary conditions? (true/false)', test)
506  IF (test) THEN
507  CALL error_petsc('===Stress boundary conditions Does not work any more')
508  END IF
509 
510  !==========Dirichlet BCs for pressure==============! !Disabled
511  !CALL find_string(21, '===How many boundary pieces for Dirichlet BCs on pressure?', test)
512  test=.false.
513  IF (test) THEN
514  READ(21,*) inputs%pp_nb_dirichlet_sides
515  ELSE
516  inputs%pp_nb_dirichlet_sides = 0
517  END IF
518  IF (inputs%pp_nb_dirichlet_sides>0) THEN
519  ALLOCATE(inputs%pp_list_dirichlet_sides(inputs%pp_nb_dirichlet_sides))
520  CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on pressure')
521  READ(21,*) inputs%pp_list_dirichlet_sides
522  ELSE
523  ALLOCATE(inputs%pp_list_dirichlet_sides(0))
524  END IF
525 
526  !==========Reynolds number=========================!
527  CALL find_string(21, '===Reynolds number',test)
528  IF (test) THEN
529  READ(21,*) inputs%Re
530  ELSE
531  CALL read_until(21, '===Kinematic viscosity')
532  READ(21,*) inputs%Re
533  inputs%Re = 1.d0 / inputs%Re
534  END IF
535 
536  !==========Variable viscosity======================!
537  CALL find_string(21, '===Variable viscosity (true/false)?',test)
538  IF (test) THEN
539  READ(21,*) inputs%if_variable_visco
540  ELSE
541  inputs%if_variable_visco = .false.
542  END IF
543 
544  !==========DIV penalty=============================!
545  IF (inputs%if_navier_stokes_with_u) THEN
546  CALL find_string(21, '===Coefficient for penalty of divergence in NS?', test)
547  IF (test) THEN
548  READ(21,*) inputs%div_stab_in_ns
549  ELSE
550  inputs%div_stab_in_ns = 0.d0
551  END IF
552  END IF
553 
554  !==========Precession==============================!
555  CALL find_string(21, '===Is there a precession term (true/false)?', test)
556  IF (test) THEN
557  READ(21,*) inputs%precession
558  ELSE
559  inputs%precession = .false.
560  END IF
561  IF (inputs%precession) THEN
562  CALL read_until(21, '===Precession rate')
563  READ(21,*) inputs%taux_precession
564  CALL read_until(21, '===Precession angle over pi')
565  READ(21,*) inputs%angle_s_pi
566  ELSE
567  inputs%taux_precession = 0.d0
568  inputs%angle_s_pi = 0.d0
569  END IF
570 
571  !==========NS penalty==============================!
572  CALL find_string(21, '===Use penalty in NS domain (true/false)?', test)
573  IF (test) THEN
574  READ(21,*) inputs%if_ns_penalty
575  ELSE
576  inputs%if_ns_penalty = .false.
577  END IF
578  IF (inputs%if_ns_penalty) THEN
579  CALL find_string(21, '===Use nonzero velocity in solids (true/false)?', test)
580  IF (test) THEN
581  READ(21,*) inputs%if_impose_vel_in_solids
582  ELSE
583  inputs%if_impose_vel_in_solids = .false.
584  END IF
585  CALL find_string(21, '===Compute z momentum (true/false)?', test)
586  IF (test) THEN
587  READ(21,*) inputs%if_compute_momentum_pseudo_force
588  ELSE
589  inputs%if_compute_momentum_pseudo_force = .false.
590  END IF
591  ELSE
592  inputs%if_impose_vel_in_solids = .false.
593  inputs%if_compute_momentum_pseudo_force = .false.
594  END IF
595 
596  !==========Solver parameters for velocity==========!
597  CALL read_until(21, '===Maximum number of iterations for velocity solver')
598  READ(21,*) inputs%my_par_vv%it_max
599  CALL read_until(21, '===Relative tolerance for velocity solver')
600  READ(21,*) inputs%my_par_vv%rel_tol
601  CALL read_until(21, '===Absolute tolerance for velocity solver')
602  READ(21,*) inputs%my_par_vv%abs_tol
603  CALL find_string(21, '===Velocity solver verbose? (true/false)', test)
604  IF (test) THEN
605  READ(21,*) inputs%my_par_vv%verbose
606  END IF
607  CALL read_until(21, '===Solver type for velocity (FGMRES, CG, ...)')
608  READ(21,*) inputs%my_par_vv%solver
609  CALL read_until(21, '===Preconditionner type for velocity solver (HYPRE, JACOBI, MUMPS...)')
610  READ(21,*) inputs%my_par_vv%precond
611 
612  !==========Solver parameters for pressure==========!
613  CALL read_until(21, '===Maximum number of iterations for pressure solver')
614  READ(21,*) inputs%my_par_pp%it_max
615  CALL read_until(21, '===Relative tolerance for pressure solver')
616  READ(21,*) inputs%my_par_pp%rel_tol
617  CALL read_until(21, '===Absolute tolerance for pressure solver')
618  READ(21,*) inputs%my_par_pp%abs_tol
619  CALL find_string(21, '===Pressure solver verbose? (true/false)',test)
620  IF (test) THEN
621  READ(21,*) inputs%my_par_pp%verbose
622  END IF
623  CALL read_until(21, '===Solver type for pressure (FGMRES, CG, ...)')
624  READ(21,*) inputs%my_par_pp%solver
625  CALL read_until(21, '===Preconditionner type for pressure solver (HYPRE, JACOBI, MUMPS...)')
626  READ(21,*) inputs%my_par_pp%precond
627 
628  !==========Solver parameters for mass matrix=======!
629  CALL read_until(21, '===Maximum number of iterations for mass matrix solver')
630  READ(21,*) inputs%my_par_mass%it_max
631  CALL read_until(21, '===Relative tolerance for mass matrix solver')
632  READ(21,*) inputs%my_par_mass%rel_tol
633  CALL read_until(21, '===Absolute tolerance for mass matrix solver')
634  READ(21,*) inputs%my_par_mass%abs_tol
635  CALL find_string(21, '===Mass matrix solver verbose? (true/false)',test)
636  IF (test) THEN
637  READ(21,*) inputs%my_par_mass%verbose
638  END IF
639  CALL read_until(21, '===Solver type for mass matrix (FGMRES, CG, ...)')
640  READ(21,*) inputs%my_par_mass%solver
641  CALL read_until(21, '===Preconditionner type for mass matrix solver (HYPRE, JACOBI, MUMPS...)')
642  READ(21,*) inputs%my_par_mass%precond
643 
644  !==========LES coefficients========================!
645  CALL find_string(21, '===Use LES? (true/false)', test)
646  IF (test) THEN
647  READ(21,*) inputs%LES
648  ELSE
649  inputs%LES = .false.
650  END IF
651 
652  IF (inputs%LES) THEN
653  CALL find_string(21, '===Coefficients for LES', test)
654  IF (test) THEN
655  CALL error_petsc('===Coefficients for LES is disabled')
656  END IF
657  !===LES_coeff4 is fixed and coeff3 is disabled
658  inputs%LES_coeff4 = 0.125d0
659  inputs%LES_coeff3 = 0.d0
660  !===LES_coeff2 in [0.15,1.0] for viscosity entropy
661  !===LES_coeff2 is equal to 1.d10 for test with first order viscosity
662  CALL read_until(21, '===Coefficient multiplying residual')
663  READ(21,*) inputs%LES_coeff2
664  !===LES_coeff1 is equal to LES_coeff4*MAX(velocity)=0.125d0*MAX(velocity)
665  CALL find_string(21, '===Coefficient for explicit LES', test)
666  IF (test) THEN
667  READ(21,*) inputs%LES_coeff1
668  ELSE
669  inputs%LES_coeff1 = 0.d0
670  END IF
671  ELSE
672  inputs%LES_coeff1 = 0.d0
673  inputs%LES_coeff2 = 0.d0
674  inputs%LES_coeff3 = 0.d0
675  inputs%LES_coeff4 = 0.d0
676  END IF
677  END IF
678  IF (.NOT.inputs%if_navier_stokes_with_u) THEN
679  CALL find_string(21, '===Use LES in momentum? (true/false)', test)
680  IF (test) THEN
681  READ(21,*) inputs%if_LES_in_momentum
682  ELSE
683  inputs%if_LES_in_momentum=.true.
684  END IF
685  IF (inputs%if_LES_in_momentum) THEN
686  inputs%LES_coeff1_mom=inputs%LES_coeff1
687  ELSE
688  inputs%LES_coeff1_mom=0.d0
689  END IF
690  END IF
691 
692  !===Maxwell data================================================================
693  IF (inputs%type_pb=='mxw' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
694  !==========Maxwell with H or B=====================!
695  CALL find_string(21, '===Solve Maxwell with H (true) or B (false)?', test)
696  IF (test) THEN
697  READ(21,*) inputs%if_maxwell_with_H
698  ELSE
699  inputs%if_maxwell_with_H = .false.
700  END IF
701 
702  !==========Quasi-static approximation==============!
703  CALL find_string(21, '===Quasi-static approximation (true) or (false)?', test)
704  IF (test) THEN
705  READ(21,*) inputs%if_quasi_static_approx
706  ELSE
707  inputs%if_quasi_static_approx = .false.
708  END IF
709 
710  !==========Steady current in fhd===================!
711  IF (inputs%type_pb=='fhd') THEN
712  CALL find_string(21, '===Steady current for fhd (true or false)?', test)
713  IF (test) THEN
714  READ(21,*) inputs%if_steady_current_fhd
715  ELSE
716  inputs%if_steady_current_fhd = .false.
717  END IF
718  END IF
719 
720  !==========Subdomains for H========================!
721  CALL read_until(21, '===Number of subdomains in magnetic field (H) mesh')
722  READ(21,*) inputs%nb_dom_H ! number of sub_domains for H
723  ALLOCATE(inputs%list_dom_H(inputs%nb_dom_H))
724  CALL read_until(21, '===List of subdomains for magnetic field (H) mesh')
725  READ(21,*) inputs%list_dom_H
726 
727  !==========Interfaces H/H==========================!
728  CALL read_until(21, '===Number of interfaces in H mesh')
729  READ(21,*) inputs%nb_inter_mu
730  ALLOCATE(inputs%list_inter_mu(inputs%nb_inter_mu))
731  IF (inputs%nb_inter_mu>0) THEN
732  CALL read_until(21, '===List of interfaces in H mesh')
733  READ(21,*) inputs%list_inter_mu
734  END IF
735 
736  !==========Dirichlet BCs for H=====================!
737  CALL read_until(21, '===Number of Dirichlet sides for Hxn')
738  READ(21,*) inputs%nb_dirichlet_sides_H
739  ALLOCATE(inputs%list_dirichlet_sides_H(inputs%nb_dirichlet_sides_H))
740  IF (inputs%nb_dirichlet_sides_H>0) THEN
741  CALL read_until(21, '===List of Dirichlet sides for Hxn')
742  READ(21,*) inputs%list_dirichlet_sides_H
743  END IF
744 
745  !==========Permeability for H======================!
746  CALL find_string(21, '===Is permeability defined analytically (true/false)?', test)
747  IF (test) THEN
748  READ(21,*) inputs%analytical_permeability
749  ELSE
750  inputs%analytical_permeability = .false.
751  END IF
752  IF (.NOT.inputs%analytical_permeability) THEN
753  CALL read_until(21, '===Permeability in the conductive part (1:nb_dom_H)')
754  ALLOCATE(inputs%mu_H(inputs%nb_dom_H))
755  READ(21,*) inputs%mu_H
756  END IF
757  inputs%if_permeability_variable_in_theta = .false.
758  IF (inputs%analytical_permeability) THEN
759  CALL find_string(21, '===Is permeability variable in theta (true/false)?', test)
760  IF (test) THEN
761  READ(21,*) inputs%if_permeability_variable_in_theta
762  END IF
763  END IF
764 
765  !==========FEM or Gaussian integration for mu_bar==!
766  inputs%if_use_fem_integration_for_mu_bar = .true.
767  IF (inputs%analytical_permeability) THEN
768  CALL find_string(21, '===Use FEM Interpolation for magnetic permeability (true/false)?', test)
769  IF (test) THEN
770  READ(21,*) inputs%if_use_fem_integration_for_mu_bar
771  END IF
772  END IF
773 
774  !==========Conductivity for H======================!
775  CALL read_until(21, '===Conductivity in the conductive part (1:nb_dom_H)')
776  ALLOCATE(inputs%sigma(inputs%nb_dom_H))
777  READ(21,*) inputs%sigma
778 
779  !LC 2017/01/27
780  !==========Minimum of conductivity=================!
781  CALL find_string(21,'===Minimum of conductivity in the whole domain', test)
782  IF (test) THEN
783  READ(21,*) inputs%sigma_min
784  ELSE
785  test_sigma_min = .true.
786  inputs%sigma_min=minval(inputs%sigma) !JLG (Feb 23/2017) !LC 1.d0
787  END IF
788 
789  !==========Minimum of permeability=================!
790  CALL find_string(21,'===Minimum of permeability in the whole domain', test)
791  IF (test) THEN
792  READ(21,*) inputs%mu_min
793  ELSE
794  IF (.NOT.inputs%analytical_permeability) THEN ! JLG (FEB 23, 2017), begin
795  inputs%mu_min = minval(inputs%mu_H)
796  ELSE
797  inputs%mu_min=1.d0 ! LC
798  END IF ! JLG (FEB 23, 2017), end
799  END IF
800 
801  !==========Finite element type=====================!
802  CALL read_until(21, '===Type of finite element for magnetic field')
803  READ(21,*) inputs%type_fe_H
804 
805  !==========Magnetic Reynolds number================!
806  CALL read_until(21, '===Magnetic Reynolds number')
807  READ(21,*) inputs%Rem
808 
809  !==========Stabilization parameters================!
810  CALL read_until(21, '===Stabilization coefficient (divergence)')
811  READ(21,*) inputs%stab(1)
812  IF (inputs%nb_inter_mu>0 .OR. inputs%nb_dirichlet_sides_H>0) THEN
813  CALL read_until(21, '===Stabilization coefficient for Dirichlet H and/or interface H/H')
814  READ(21,*) inputs%stab(3)
815  ELSE
816  inputs%stab(3) = 0.d0
817  END IF
818 
819  !==========Subdomains for phi======================!
820  CALL find_string(21, '===Number of subdomains in magnetic potential (phi) mesh', test)
821  IF (test) THEN
822  READ (21,*) inputs%nb_dom_phi
823  ELSE
824  inputs%nb_dom_phi = 0
825  END IF
826  ALLOCATE(inputs%list_dom_phi(inputs%nb_dom_phi))
827  IF (inputs%nb_dom_phi>0) THEN
828  !==========List of subdomains for phi======================!
829  CALL read_until(21, '===List of subdomains for magnetic potential (phi) mesh')
830  READ(21,*) inputs%list_dom_phi
831 
832  !==========Dirichlet BCs on phi====================!
833  CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on phi?')
834  READ(21,*) inputs%phi_nb_dirichlet_sides
835  IF (inputs%phi_nb_dirichlet_sides>0) THEN
836  ALLOCATE(inputs%phi_list_dirichlet_sides(inputs%phi_nb_dirichlet_sides))
837  CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on phi')
838  READ(21,*) inputs%phi_list_dirichlet_sides
839  ELSE
840  ALLOCATE(inputs%phi_list_dirichlet_sides(0))
841  END IF
842 
843  !==========H/phi interface=========================!
844  CALL read_until(21, '===Number of interfaces between H and phi')
845  READ(21,*) inputs%nb_inter ! number of interfaces between H and phi
846  ALLOCATE(inputs%list_inter_H_phi(inputs%nb_inter))
847  CALL read_until(21, '===List of interfaces between H and phi')
848  IF (inputs%nb_inter>0) THEN
849  READ(21,*) inputs%list_inter_H_phi
850  END IF
851 
852  !==========Permeability in vacuum==================!
853  CALL read_until(21, '===Permeability in vacuum')
854  READ(21,*) inputs%mu_phi
855 
856  !==========Finite element type=====================!
857  CALL read_until(21, '===Type of finite element for scalar potential')
858  READ(21,*) inputs%type_fe_phi
859 
860  !==========Stabilization parameters================!
861  CALL read_until(21, '===Stabilization coefficient (interface H/phi)')
862  READ(21,*) inputs%stab(2)
863  END IF
864 
865  !==========Solver parameters=======================!
866  CALL read_until(21, '===Maximum number of iterations for Maxwell solver')
867  READ(21,*) inputs%my_par_H_p_phi%it_max
868  CALL read_until(21, '===Relative tolerance for Maxwell solver')
869  READ(21,*) inputs%my_par_H_p_phi%rel_tol
870  CALL read_until(21, '===Absolute tolerance for Maxwell solver')
871  READ(21,*) inputs%my_par_H_p_phi%abs_tol
872  CALL find_string(21, '===Maxwell solver verbose? (true/false)',test)
873  IF (test) THEN
874  READ(21,*) inputs%my_par_H_p_phi%verbose
875  END IF
876  CALL read_until(21, '===Solver type for Maxwell (FGMRES, CG, ...)')
877  READ(21,*) inputs%my_par_H_p_phi%solver
878  CALL read_until(21, '===Preconditionner type for Maxwell solver (HYPRE, JACOBI, MUMPS...)')
879  READ(21,*) inputs%my_par_H_p_phi%precond
880  END IF
881 
882  !===Data temperature==============================================================
883  CALL find_string(21, '===Is there a temperature field?', test)
884  IF (test) THEN
885  READ (21,*) inputs%if_temperature
886  ELSE
887  inputs%if_temperature = .false.
888  END IF
889  IF (inputs%if_temperature) THEN
890  !==========Gravity coefficient for temperature=====!
891  CALL find_string(21, '===Non-dimensional gravity coefficient', test)
892  IF (test) THEN
893  READ (21,*) inputs%gravity_coefficient
894  ELSE
895  inputs%gravity_coefficient = 0.d0
896  END IF
897  !==========Magnetic force coefficient==============!
898  CALL find_string(21, '===Helmholtz force for ferrohydrodynamics? (true/false)', test)
899  IF (test) THEN
900  READ (21,*) inputs%if_helmholtz_force
901  END IF
902  CALL find_string(21, '===Non-dimensional magnetic force coefficient for ferrohydrodynamics', test) ! MODIFICATION: coefficient for magnetic force
903  IF (test) THEN
904  READ (21,*) inputs%mag_force_coefficient
905  END IF
906  !==========Subdomains for temp=====================!
907  CALL read_until(21, '===Number of subdomains in temperature mesh')
908  READ(21,*) inputs%nb_dom_temp ! number of sub_domains for temp
909  ALLOCATE(inputs%list_dom_temp(inputs%nb_dom_temp))
910  CALL read_until(21, '===List of subdomains for temperature mesh')
911  READ(21,*) inputs%list_dom_temp
912  ALLOCATE(inputs%vol_heat_capacity(inputs%nb_dom_temp)) ! Two choices for the user, indicate 1) the heat capacity and the conductivity (necessary if non uniform heat capacity) or 2) the diffusivity only
913  ALLOCATE(inputs%temperature_diffusivity(inputs%nb_dom_temp)) ! In both cases, temperature_diffusivity is used, it contains the conductivities in 1) and the diffusivities in 2)
914  CALL find_string(21, '===Density (1:nb_dom_temp)', test)
915  IF (test) THEN ! Case 1)
916  ALLOCATE(inputs%density(inputs%nb_dom_temp))
917  ALLOCATE(inputs%heat_capacity(inputs%nb_dom_temp))
918  READ(21,*) inputs%density
919  CALL read_until(21, '===Heat capacity (1:nb_dom_temp)')
920  READ(21,*) inputs%heat_capacity
921  inputs%vol_heat_capacity = inputs%density * inputs%heat_capacity
922  DEALLOCATE (inputs%density)
923  DEALLOCATE (inputs%heat_capacity)
924  CALL read_until(21, '===Thermal conductivity (1:nb_dom_temp)')
925  READ(21,*) inputs%temperature_diffusivity
926  ELSE ! Case 1) bis
927  CALL find_string(21, '===Volumetric heat capacity (1:nb_dom_temp)', test)
928  IF (test) THEN
929  READ(21,*) inputs%vol_heat_capacity
930  CALL read_until(21, '===Thermal conductivity (1:nb_dom_temp)')
931  READ(21,*) inputs%temperature_diffusivity
932  ELSE ! Case 2)
933  inputs%vol_heat_capacity = 1.d0 ! Heat capacity is equalized to one so that it does not impact the calculus
934  CALL read_until(21, '===Diffusivity coefficient for temperature (1:nb_dom_temp)')
935  READ(21,*) inputs%temperature_diffusivity
936  END IF
937  END IF
938  !==========Dirichlet BCs on temperature============!
939  CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on temperature?')
940  READ(21,*) inputs%temperature_nb_dirichlet_sides
941  IF (inputs%temperature_nb_dirichlet_sides>0) THEN
942  ALLOCATE(inputs%temperature_list_dirichlet_sides(inputs%temperature_nb_dirichlet_sides))
943  CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on temperature')
944  READ(21,*) inputs%temperature_list_dirichlet_sides
945  ELSE
946  ALLOCATE(inputs%temperature_list_dirichlet_sides(0))
947  END IF
948  !==========Robin BCs on temperature================!
949  CALL find_string(21, '===How many boundary pieces for Robin BCs on temperature?', test)
950  IF (test) THEN
951  READ(21,*) inputs%temperature_nb_robin_sides
952  ELSE
953  inputs%temperature_nb_robin_sides = 0
954  END IF
955  IF (test .AND. (inputs%temperature_nb_robin_sides>0)) THEN
956  ALLOCATE(inputs%temperature_list_robin_sides(inputs%temperature_nb_robin_sides))
957  CALL read_until(21, '===List of boundary pieces for Robin BCs on temperature')
958  READ(21,*) inputs%temperature_list_robin_sides
959  ALLOCATE(inputs%convection_coeff(inputs%temperature_nb_robin_sides))
960  CALL read_until(21, '===Convection heat transfert coefficient (1:temperature_nb_robin_sides)')
961  READ(21,*) inputs%convection_coeff
962  ALLOCATE(inputs%exterior_temperature(inputs%temperature_nb_robin_sides))
963  CALL read_until(21, '===Exterior temperature (1:temperature_nb_robin_sides)')
964  READ(21,*) inputs%exterior_temperature
965  ELSE
966  ALLOCATE(inputs%temperature_list_robin_sides(0))
967  END IF
968  !==========Interfaces between vel and temp=========!
969  CALL find_string(21, '===Number of interfaces between velocity and temperature only domains (for nst applications)', test)
970  IF (test) THEN
971  READ(21,*) inputs%nb_inter_v_T
972  ALLOCATE(inputs%list_inter_v_T(inputs%nb_inter_v_T))
973  CALL read_until(21, '===List of interfaces between velocity and temperature only domains (for nst applications)')
974  READ(21,*) inputs%list_inter_v_T
975  ELSE
976  ALLOCATE(inputs%list_inter_v_T(0))
977  END IF
978  !==========Solver parameters=======================!
979  CALL read_until(21, '===Maximum number of iterations for temperature solver')
980  READ(21,*) inputs%my_par_temperature%it_max
981  CALL read_until(21, '===Relative tolerance for temperature solver')
982  READ(21,*) inputs%my_par_temperature%rel_tol
983  CALL read_until(21, '===Absolute tolerance for temperature solver')
984  READ(21,*) inputs%my_par_temperature%abs_tol
985  CALL find_string(21, '===Temperature solver verbose? (true/false)',test)
986  IF (test) THEN
987  READ(21,*) inputs%my_par_temperature%verbose
988  END IF
989  CALL read_until(21, '===Solver type for temperature (FGMRES, CG, ...)')
990  READ(21,*) inputs%my_par_temperature%solver
991  CALL read_until(21, '===Preconditionner type for temperature solver (HYPRE, JACOBI, MUMPS...)')
992  READ(21,*) inputs%my_par_temperature%precond
993  END IF
994 
995  !===Data Level set================================================================
996  CALL find_string(21, '===Is there a level set?', test)
997  IF (test) THEN
998  READ (21,*) inputs%if_level_set
999  ELSE
1000  inputs%if_level_set = .false.
1001  END IF
1002  IF (inputs%if_level_set) THEN
1003  !==========Fix level set to test code==============!
1004  CALL find_string(21, '===Do we fix level set? (true/false)', test)
1005  IF (test) THEN
1006  READ (21, *) inputs%if_level_set_fixed
1007  ELSE
1008  inputs%if_level_set_fixed = .false.
1009  END IF
1010  !==========Number of fluids========================!
1011  CALL read_until(21, '===How many fluids?')
1012  READ(21,*) inputs%nb_fluid
1013  !==========Level_set multiplier for h_min==========!
1014  CALL read_until(21, '===multiplier for h_min for level set')
1015  READ(21,*) inputs%multiplier_for_h_min_distance
1016  !==========Level_set c_max=========================!
1017  !===Disabled functionality
1018  inputs%level_set_cmax=0.d0
1019  !==========Level_set compression factor============!
1020  CALL read_until(21, '===Compression factor for level set')
1021  READ(21,*) inputs%level_set_comp_factor
1022  !==========Densities of fluids=====================!
1023  ALLOCATE(inputs%density_fluid( inputs%nb_fluid))
1024  CALL read_until(21, '===Density of fluid 0, fluid 1, ...')
1025  READ(21,*) inputs%density_fluid
1026  !==========Dynamic viscosities of fluids===========!
1027  ALLOCATE(inputs%dyna_visc_fluid( inputs%nb_fluid))
1028  CALL read_until(21, '===Dynamic viscosity of fluid 0, fluid 1, ...')
1029  READ(21,*) inputs%dyna_visc_fluid
1030  !==========Conductivities of fluids================!
1031  IF (inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
1032  CALL find_string(21, '===Conductivity of fluid 0, fluid 1, ...', test)
1033  IF (test) THEN
1034  inputs%variation_sigma_fluid =.true.
1035  ALLOCATE(inputs%sigma_fluid(inputs%nb_fluid))
1036  READ(21,*) inputs%sigma_fluid
1037  !===sigma_min has not been been read=========!
1038  IF(test_sigma_min) THEN
1039  inputs%sigma_min=minval(inputs%sigma_fluid) !===JLG May 5,
1040  END IF
1041  ELSE
1042  inputs%variation_sigma_fluid =.false.
1043  END IF
1044  END IF
1045  !==========Surface tension=========================!
1046  CALL find_string(21, '===Is there a surface tension?', test)
1047  IF (test) THEN
1048  READ (21,*) inputs%if_surface_tension
1049  ELSE
1050  inputs%if_surface_tension = .false.
1051  END IF
1052  IF (inputs%if_surface_tension) THEN
1053  !==========Coefficient for surface tension======!
1054  CALL read_until(21, '===Coefficients of surface tension for level set 0, level set 1, ...')
1055  ALLOCATE(inputs%coeff_surface(inputs%nb_fluid-1))
1056  READ(21,*) inputs%coeff_surface
1057  END IF
1058  !==========Mass correction=========================!
1059  CALL find_string(21, '===Do we apply mass correction? (true/false)', test)
1060  IF (test) THEN
1061  READ (21,*) inputs%if_mass_correction
1062  ELSE
1063  inputs%if_mass_correction=.true.
1064  END IF
1065  !==========Dirichlet BCs on level_set==============!
1066  CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on level set?')
1067  READ(21,*) inputs%level_set_nb_dirichlet_sides
1068  IF (inputs%level_set_nb_dirichlet_sides>0) THEN
1069  ALLOCATE(inputs%level_set_list_dirichlet_sides(inputs%level_set_nb_dirichlet_sides))
1070  CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on level set')
1071  READ(21,*) inputs%level_set_list_dirichlet_sides
1072  ELSE
1073  ALLOCATE(inputs%level_set_list_dirichlet_sides(0))
1074  END IF
1075  !==========Solver parameters=======================!
1076  CALL read_until(21, '===Maximum number of iterations for level set solver')
1077  READ(21,*) inputs%my_par_level_set%it_max
1078  CALL read_until(21, '===Relative tolerance for level set solver')
1079  READ(21,*) inputs%my_par_level_set%rel_tol
1080  CALL read_until(21, '===Absolute tolerance for level set solver')
1081  READ(21,*) inputs%my_par_level_set%abs_tol
1082  CALL find_string(21, '===Level set solver verbose? (true/false)',test)
1083  IF (test) THEN
1084  READ(21,*) inputs%my_par_level_set%verbose
1085  END IF
1086  CALL read_until(21, '===Solver type for level set (FGMRES, CG, ...)')
1087  READ(21,*) inputs%my_par_level_set%solver
1088  CALL read_until(21, '===Preconditionner type for level set solver (HYPRE, JACOBI, MUMPS...)')
1089  READ(21,*) inputs%my_par_level_set%precond
1090 
1091  !==========Reconstruction parameters===============!
1092  CALL find_string(21, '===How are the variables reconstructed from the level set function? (lin, reg)', test)
1093  IF (test) THEN
1094  READ (21,*) inputs%level_set_reconstruction_type
1095  ELSE
1096  inputs%level_set_reconstruction_type = 'reg'
1097  END IF
1098  IF (inputs%level_set_reconstruction_type=='reg') THEN
1099  CALL find_string(21, '===Value of the regularization coefficient in (0,0.5]', test)
1100  IF (test) THEN
1101  READ (21,*) inputs%level_set_tanh_coeff_reconstruction
1102  IF (inputs%level_set_tanh_coeff_reconstruction.LE.1d-2 .OR. &
1103  .5d0 .LE. inputs%level_set_tanh_coeff_reconstruction) THEN
1104  inputs%level_set_tanh_coeff_reconstruction = 0.45d0
1105  END IF
1106  ELSE
1107  inputs%level_set_tanh_coeff_reconstruction = 0.45d0
1108  END IF
1109  ELSE IF (inputs%level_set_reconstruction_type/='lin') THEN
1110  CALL error_petsc('BUG in read_my_data: variable reconstruction type not correct')
1111  END IF
1112  !==========Level set overshoot=====================!
1113  IF(inputs%level_set_reconstruction_type/='reg') THEN
1114  CALL find_string(21, '===Do we kill level set overshoot? (true/false)', test)
1115  IF (test) THEN
1116  READ (21,*) inputs%if_kill_overshoot
1117  ELSE
1118  inputs%if_kill_overshoot=.false.
1119  END IF
1120  ELSE
1121  inputs%if_kill_overshoot=.false.
1122  END IF
1123 
1124  !==========Level set type finite element===========!
1125  CALL find_string(21, '===Do we use P2 finite element for level_set? (true/false)', test)
1126  IF (test) THEN
1127  READ(21,*) inputs%if_level_set_P2
1128  ELSE
1129  inputs%if_level_set_P2=.true.
1130  END IF
1131 
1132  !==========Level set compression method============!
1133  CALL find_string(21, '===Do we use JLG compression method for level_set? (true/false)', test)
1134  IF (test) THEN
1135  READ(21,*) inputs%if_compression_mthd_JLG
1136  ELSE
1137  inputs%if_compression_mthd_JLG=.true.
1138  END IF
1139 
1140  !==========Check coherence=========================!
1141  IF (inputs%type_pb=='mxw') THEN
1142  CALL error_petsc('BUG in read_my_data: Level_set without Navier-Stokes')
1143  END IF
1144  END IF
1145 
1146  !===Data for arpack (eigenvalue problems)==========!
1147  !==========Frequency parameters====================!
1148  IF (inputs%type_pb=='mxw') THEN
1149  CALL find_string(21, '===Do we use Arpack?', test)
1150  IF (test) THEN
1151  READ (21,*) inputs%if_arpack
1152  ELSE
1153  inputs%if_arpack = .false.
1154  END IF
1155  IF (inputs%if_arpack)THEN
1156  CALL read_until(21, '===Number of eigenvalues to compute')
1157  READ (21,*) inputs%nb_eigenvalues
1158  CALL read_until(21, '===Maximum number of Arpack iterations')
1159  READ (21,*) inputs%arpack_iter_max
1160  CALL read_until(21, '===Tolerance for Arpack')
1161  READ (21,*) inputs%arpack_tolerance
1162  CALL read_until(21, &
1163  '===Which eigenvalues (''LM'', ''SM'', ''SR'', ''LR'' ''LI'', ''SI'')')
1164  READ (21,*) inputs%arpack_which
1165  CALL find_string(21, '===Create 2D vtu files for Arpack? (true/false)', test)
1166  IF (test) THEN
1167  READ (21,*) inputs%if_arpack_vtu_2d
1168  ELSE
1169  inputs%if_arpack_vtu_2d = .false.
1170  END IF
1171  END IF
1172  ELSE
1173  !==========Arpack currently works only for Maxwell=!
1174  inputs%if_arpack = .false.
1175  END IF
1176 
1177  !===Format for paraview============================!
1178  CALL find_string(21, '===Vtu files in xml format? (true=xml/false=ascii)', test)
1179  IF (test) THEN
1180  READ (21,*) inputs%if_xml
1181  ELSE
1182  inputs%if_xml=.true.
1183  END IF
1184 
1185  !===Data for post processing=======================!
1186  CALL find_string(21, '===Number of planes in real space for Visualization', test)
1187  IF (test) THEN
1188  READ (21,*) inputs%number_of_planes_in_real_space
1189  ELSE
1190  inputs%number_of_planes_in_real_space = 10
1191  END IF
1192 
1193  !==========Frequency parameters====================!
1194  CALL find_string(21, '===Frequency to write restart file', test)
1195  IF (test) THEN
1196  READ (21,*) inputs%freq_restart
1197  ELSE
1198  inputs%freq_restart = 100000000
1199  END IF
1200  CALL find_string(21, '===Frequency to write energies', test)
1201  IF (test) THEN
1202  READ (21,*) inputs%freq_en
1203  ELSE
1204  inputs%freq_en = 100000000
1205  END IF
1206  CALL find_string(21, '===Frequency to create plots', test)
1207  IF (test) THEN
1208  READ (21,*) inputs%freq_plot
1209  ELSE
1210  inputs%freq_plot = 100000000
1211  END IF
1212  CALL find_string(21, '===Just postprocessing without computing? (true/false)', test)
1213  IF (test) THEN
1214  READ (21,*) inputs%if_just_processing
1215  ELSE
1216  inputs%if_just_processing = .false.
1217  END IF
1218 
1219 
1220  !==========Modes to be zeroed out==================!
1221  IF (inputs%type_pb=='mhd') THEN
1222  CALL find_string(21, '===Should some modes be zeroed out?', test)
1223  IF (test) THEN
1224  READ (21,*) inputs%if_zero_out_modes
1225  ELSE
1226  inputs%if_zero_out_modes = .false.
1227  END IF
1228  IF (inputs%if_zero_out_modes) THEN
1229  CALL read_until(21, '===How many Navier-Stokes modes to zero out?')
1230  READ(21,*) inputs%nb_select_mode_ns
1231  ALLOCATE(inputs%list_select_mode_ns(inputs%nb_select_mode_ns))
1232  CALL read_until(21, '===List of Navier-Stokes modes to zero out?')
1233  READ(21,*) inputs%list_select_mode_ns
1234  CALL read_until(21, '===How Maxwell modes to zero out?')
1235  READ(21,*) inputs%nb_select_mode_mxw
1236  ALLOCATE(inputs%list_select_mode_mxw(inputs%nb_select_mode_mxw))
1237  CALL read_until(21, '===List of Maxwell modes to zero out?')
1238  READ(21,*) inputs%list_select_mode_mxw
1239  END IF
1240  END IF
1241 
1242  !==========Verbose=================================!
1243  CALL find_string(21, '===Verbose timing? (true/false)', test)
1244  IF (test) THEN
1245  READ (21,*) inputs%verbose_timing
1246  ELSE
1247  inputs%verbose_timing = .false.
1248  END IF
1249  CALL find_string(21, '===Verbose divergence? (true/false)', test)
1250  IF (test) THEN
1251  READ (21,*) inputs%verbose_divergence
1252  ELSE
1253  inputs%verbose_divergence = .false.
1254  END IF
1255  CALL find_string(21, '===Verbose CFL? (true/false)', test)
1256  IF (test) THEN
1257  READ (21,*) inputs%verbose_CFL
1258  ELSE
1259  inputs%verbose_CFL = .false.
1260  END IF
1261 
1262  !===Norms for reference tests===================================================
1263  IF (inputs%test_de_convergence) THEN
1264  CALL read_until(21, '===Reference results')
1265  DO k = 1, 4
1266  READ(21,*) inputs%norm_ref(k)
1267  END DO
1268  ELSE
1269  inputs%norm_ref = 1.d0
1270  END IF
1271 
1272  !===mxw becomes mxx if restart_u and mxw========================================
1273  IF (inputs%type_pb == 'mxw' .AND. inputs%irestart_u) THEN
1274  inputs%type_pb = 'mxx'
1275  END IF
1276 
1277  CLOSE(21)
1278 
1279  !===Check coherence of data=====================================================
1280  CALL check_coherence_of_data
1281  RETURN
1282 
1283  END SUBROUTINE read_my_data
1284 
1285  SUBROUTINE check_coherence_of_data
1286  USE my_util
1287  IMPLICIT NONE
1288  INTEGER :: k, mode_min, mode_max, Delta_mode
1289  LOGICAL :: test
1290 
1291  !===Dirichlet BCs for Navier-Stokes=============================================
1292  IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
1293  IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd' .OR. inputs%irestart_u) THEN ! MODIFICATION: fhd added
1294  !==========Velocity================================!
1295  DO k = 1, 3
1296  IF (inputs%vv_nb_dirichlet_sides(k)<1) cycle
1297  test = check_coherence_with_periodic_bcs(inputs%vv_list_dirichlet_sides(k)%DIL)
1298  IF (test) THEN
1299  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1300  ' and periodic BCs on velocity')
1301  END IF
1302  END DO
1303  !==========Pressure================================!
1304  IF (inputs%pp_nb_dirichlet_sides>0) THEN
1305  test = check_coherence_with_periodic_bcs(inputs%pp_list_dirichlet_sides)
1306  IF (test) THEN
1307  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1308  ' and periodic BCs on pressure')
1309  END IF
1310  END IF
1311  END IF
1312  END IF
1313 
1314  !===No penalty with momentum formulation========================================
1315  IF (.NOT.inputs%if_navier_stokes_with_u .AND. inputs%if_ns_penalty) THEN
1316  CALL error_petsc('BUG in read_my_data: No penalty with momentum formulation')
1317  END IF
1318 
1319  !===Dirichlet BCs for temperature for Navier-Stokes=============================
1320  IF (inputs%if_temperature.AND. inputs%my_periodic%nb_periodic_pairs>0) THEN
1321  IF (inputs%temperature_nb_dirichlet_sides>0) THEN
1322  test = check_coherence_with_periodic_bcs(inputs%temperature_list_dirichlet_sides)
1323  IF (test) THEN
1324  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1325  ' and periodic BCs on temperature')
1326  END IF
1327  END IF
1328  END IF
1329 
1330  !===Dirichlet BCs for Level_set for Navier-Stokes=============================
1331  IF (inputs%if_level_set.AND. inputs%my_periodic%nb_periodic_pairs>0) THEN
1332  IF (inputs%level_set_nb_dirichlet_sides>0) THEN
1333  test = check_coherence_with_periodic_bcs(inputs%level_set_list_dirichlet_sides)
1334  IF (test) THEN
1335  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1336  ' and periodic BCs on level_set')
1337  END IF
1338  END IF
1339  END IF
1340 
1341  !===Robin and Dirichlet BCs===================================================== ! MODIFICATION/ The user should not specify sides that are Robin and Dirichlet
1342  IF (inputs%temperature_nb_robin_sides>0) THEN
1343  DO k = 1, inputs%temperature_nb_robin_sides
1344  IF (minval(abs(inputs%temperature_list_dirichlet_sides - inputs%temperature_list_robin_sides(k))) == 0) THEN
1345  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1346  ' and Robin BCs for temperature')
1347  END IF
1348  END DO
1349  END IF
1350 
1351  !===Dirichlet BCs magnetic field for Maxwell====================================
1352  IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
1353  IF (inputs%type_pb=='mxw' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
1354  !==========Magnetic field==========================!
1355  IF (inputs%nb_dirichlet_sides_H>0) THEN
1356  test = check_coherence_with_periodic_bcs(inputs%list_dirichlet_sides_H)
1357  IF (test) THEN
1358  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1359  ' and periodic BCs on magnetic field')
1360  END IF
1361  END IF
1362 
1363  !==========Scalar potential========================!
1364  IF (inputs%nb_dom_phi>0 .AND. inputs%phi_nb_dirichlet_sides>0) THEN
1365  test = check_coherence_with_periodic_bcs(inputs%phi_list_dirichlet_sides)
1366  IF (test) THEN
1367  CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
1368  ' and periodic BCs on scalar potential')
1369  END IF
1370  END IF
1371  END IF
1372  END IF
1373 
1374  !===Check temperature with Maxwell==============================================
1375  IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
1376  IF (inputs%if_temperature .AND. inputs%type_pb=='mxw') THEN
1377  CALL error_petsc('Bug in read_my_data: incompatible temperature with maxwell')
1378  END IF
1379  END IF
1380 
1381  !===Check temperature with Ferrohydrodynamics===================================
1382  IF ((.NOT. inputs%if_temperature) .AND. inputs%type_pb=='fhd') THEN
1383  CALL error_petsc('Bug in read_my_data: ferrohydrodynamics but no temperature')
1384  END IF
1385 
1386  !===Check Arpack================================================================
1387  IF (inputs%if_arpack) THEN
1388  IF (inputs%ndim(2) /= inputs%m_max) THEN
1389  CALL error_petsc('Bug in read_my_data: #Fourier modes'// &
1390  ' not equal to #processors in Fourier direction')
1391  END IF
1392  END IF
1393 
1394  !===Check Fourier modes=========================================================
1395  IF (inputs%select_mode .AND. .NOT.inputs%if_arpack) THEN
1396  IF (SIZE(inputs%list_mode_lect)/=1) THEN
1397  mode_max = maxval(inputs%list_mode_lect)
1398  mode_min = minval(inputs%list_mode_lect)
1399  IF (mod(mode_max-mode_min,SIZE(inputs%list_mode_lect)-1)/=0) THEN
1400  CALL error_petsc('Bug in read_my_data: Fourier modes not equally spaced ')
1401  END IF
1402  delta_mode = (mode_max-mode_min)/(SIZE(inputs%list_mode_lect)-1)
1403  DO k = 0, SIZE(inputs%list_mode_lect)-1
1404  IF (minval(abs(inputs%list_mode_lect-(delta_mode*k+mode_min)))/=0) THEN
1405  CALL error_petsc('Bug in read_my_data: Fourier modes not equally spaced ')
1406  END IF
1407  END DO
1408  DO k = 1, SIZE(inputs%list_mode_lect)-1
1409  IF (inputs%list_mode_lect(k+1).LE.inputs%list_mode_lect(k)) THEN
1410  CALL error_petsc('Bug in read_my_data: Fourier modes not in increasing order ')
1411  END IF
1412  END DO
1413  END IF
1414  END IF
1415 
1416  !===Irestart for postprocessing=================================================
1417  IF (inputs%if_just_processing) THEN
1418  inputs%irestart_u = .false.
1419  inputs%irestart_h = .false.
1420  IF (inputs%type_pb/='mxw') inputs%irestart_u = .true.
1421  IF (inputs%type_pb/='nst') inputs%irestart_h = .true.
1422  END IF
1423 
1424  !===Allocate dummy list_dom_* for new partitioning==============================
1425  IF (.NOT. ASSOCIATED(inputs%list_dom_ns)) ALLOCATE(inputs%list_dom_ns(0))
1426  IF (.NOT. ASSOCIATED(inputs%list_dom_H)) ALLOCATE(inputs%list_dom_H(0))
1427  IF (.NOT. ASSOCIATED(inputs%list_dom_phi)) ALLOCATE(inputs%list_dom_phi(0))
1428  IF (.NOT. ASSOCIATED(inputs%list_dom_temp)) ALLOCATE(inputs%list_dom_temp(0))
1429 
1430  END SUBROUTINE check_coherence_of_data
1431 
1432  FUNCTION check_coherence_with_periodic_bcs(list) RESULT(test)
1433  IMPLICIT NONE
1434  INTEGER, DIMENSION(:) :: list
1435  LOGICAL :: test
1436  INTEGER :: k, n
1437 
1438  test = .false.
1439  DO k = 1, inputs%my_periodic%nb_periodic_pairs
1440  DO n = 1, SIZE(list)
1441  IF (minval(abs(list(n)-inputs%my_periodic%list_periodic(:,k)))==0) THEN
1442  test = .true.
1443  RETURN
1444  END IF
1445  END DO
1446  END DO
1447  END FUNCTION check_coherence_with_periodic_bcs
1448 
1449 END MODULE input_data
1450 
input_data::read_my_data
subroutine, public read_my_data(data_fichier)
Definition: read_my_data.f90:243
chaine_caractere::find_string
subroutine find_string(unit, string, okay)
Definition: chaine_caractere.f90:119
def_type_mesh::periodic_data
Definition: def_type_mesh.f90:12
my_util
Definition: my_util.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
chaine_caractere::read_until
subroutine read_until(unit, string, error)
Definition: chaine_caractere.f90:96
my_data_module
Definition: read_my_data.f90:1
dyn_line::dyn_int_line
Definition: dyn_line.f90:5
my_data_module::my_data
Definition: read_my_data.f90:4
input_data::check_coherence_of_data
subroutine check_coherence_of_data
Definition: read_my_data.f90:1286
my_data_module::init
subroutine init(a)
Definition: read_my_data.f90:146
solve_petsc
Definition: solver.f90:1
input_data::check_coherence_with_periodic_bcs
logical function check_coherence_with_periodic_bcs(list)
Definition: read_my_data.f90:1433
solve_petsc::solver_param
Definition: solver.f90:3
input_data
Definition: read_my_data.f90:233
chaine_caractere
Definition: chaine_caractere.f90:4
def_type_mesh
Definition: def_type_mesh.f90:4
form
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