17 LOGICAL :: if_momentum, if_mass, if_induction, if_energy
23 REAL(KIND=8),
TARGET,
ALLOCATABLE,
DIMENSION(:,:,:) :: un, un_m1
25 REAL(KIND=8),
TARGET,
ALLOCATABLE,
DIMENSION(:,:,:) :: pn, pn_m1
27 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:,:,:) :: incpn, incpn_m1
31 REAL(KIND=8),
TARGET,
ALLOCATABLE,
DIMENSION(:,:,:,:) :: level_set, level_set_m1
33 REAL(KIND=8),
TARGET,
ALLOCATABLE,
DIMENSION(:,:,:) :: density, density_m1, density_m2
35 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:,:) :: visc_entro_level
37 REAL(KIND=8) :: max_vel
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
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
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
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
80 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: jj_v_to_H
81 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: jj_v_to_temp
82 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: jj_temp_to_H
86 INTEGER,
TARGET,
ALLOCATABLE,
DIMENSION(:) :: list_mode
92 REAL(KIND=8) :: R_fourier
93 INTEGER :: index_fourier
97 #include "petsc/finclude/petsc.h"
99 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d, comm_one_d_ns, comm_one_d_temp, coord_cart
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)
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
134 list_mode, inputs%number_of_planes_in_real_space)
136 vv_mesh_out => vv_mesh
137 pp_mesh_out => pp_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
147 level_set_out => level_set
148 density_out => density
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
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
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
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,&
186 inputs%dyna_visc_fluid, visco_dyn_m1)
187 IF (inputs%variation_sigma_fluid)
THEN
189 inputs%sigma_fluid, sigma_ns)
196 IF (if_momentum)
THEN
199 IF (if_induction)
THEN
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)
208 IF (if_momentum)
THEN
212 IF (if_induction)
THEN
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, &
223 IF (if_induction)
THEN
227 IF (if_momentum)
THEN
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)
243 LOGICAL,
SAVE :: once_zero_out_mode=.true.
244 INTEGER,
POINTER,
DIMENSION(:),
SAVE :: select_mode_ns, select_mode_mxw
246 IF (once_zero_out_mode)
THEN
247 once_zero_out_mode=.false.
248 IF (inputs%if_zero_out_modes)
THEN
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
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
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
278 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode_to_zero_out
279 INTEGER,
DIMENSION(:) :: list_mode
280 INTEGER,
POINTER,
DIMENSION(:) :: select_mode
282 INTEGER,
DIMENSION(1) :: kloc
284 DO i = 1,
SIZE(list_mode_to_zero_out)
285 IF (minval(abs(list_mode-list_mode_to_zero_out(i)))==0)
THEN
289 ALLOCATE(select_mode(inc))
291 DO i = 1,
SIZE(list_mode_to_zero_out)
292 IF (minval(abs(list_mode-list_mode_to_zero_out(i)))==0)
THEN
294 kloc = minloc(abs(list_mode-list_mode_to_zero_out(i)))
295 select_mode(inc) = kloc(1)
306 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: vn
307 INTEGER,
DIMENSION(:),
INTENT(IN) :: connectivity_structure
308 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: coupling_variable
318 IF (mesh%me /=0)
THEN
319 DO j = 1,
SIZE(connectivity_structure)
320 IF (connectivity_structure(j) == -1) cycle
321 coupling_variable(j,:,:) = vn(connectivity_structure(j),:,:)
328 SUBROUTINE projection_temperature(mesh, vn, connectivity_structure, if_restriction, coupling_variable) !projection of temp on another mesh subroutine added
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
337 coupling_variable = 0.d0
339 IF (mesh%me /=0 )
THEN
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,:,:)
345 coupling_variable(j,:,:) = vn(connectivity_structure(j),:,:)
353 SUBROUTINE projection_mag_field(mesh, vn, connectivity_structure, coupling_variable) !projection of mag field on another mesh subroutine added
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
362 DO j = 1,
SIZE(connectivity_structure)
363 IF (connectivity_structure(j) == -1) cycle
364 coupling_variable(connectivity_structure(j),:,:) = vn(j,:,:)
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) ')
379 coupling_variable(connectivity_structure(j1),:,:) = (vn(j1,:,:)+vn(j2,:,:))/2
382 IF (h_mesh%gauss%n_w/=mesh%gauss%n_w)
THEN
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
397 INTEGER,
INTENT(IN) :: it, freq_restart
399 IF (if_momentum)
THEN
400 IF (pp_mesh%me /= 0)
THEN
401 IF (.NOT. if_mass)
THEN
403 list_mode, un, un_m1, pn, pn_m1, &
404 incpn, incpn_m1, inputs%file_name, it, freq_restart)
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)
414 IF (if_induction)
THEN
416 time, list_mode, hn, hn1, bn, bn1, &
417 phin, phin1, inputs%file_name, it, freq_restart)
422 list_mode, tempn, tempn_m1, inputs%file_name, it, freq_restart)
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
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
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
473 CALL mpi_comm_size(petsc_comm_world,nb_procs,code)
474 CALL mpi_comm_rank(petsc_comm_world,petsc_rank,code)
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')
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))
489 data_fichier = trim(adjustl(data_directory))//
'/'//trim(adjustl(data_file))
499 IF (inputs%test_de_convergence)
THEN
500 inputs%directory = data_directory
504 IF (inputs%nb_dom_phi==0)
THEN
505 inputs%phi_nb_dirichlet_sides = 0
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))
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')
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')
532 m_max_c = inputs%m_max/nb_procs_f
533 ALLOCATE(list_mode(m_max_c))
537 IF (inputs%select_mode)
THEN
539 list_mode(i) = inputs%list_mode_lect(i + rank*m_max_c)
543 list_mode(i) = i + rank*m_max_c - 1
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
554 ns_periodic=.true.; mxw_periodic=.true.; temp_periodic=.true.
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
564 IF (if_induction)
THEN
565 ALLOCATE(list_dom_h(inputs%nb_dom_H), h_in_to_new(inputs%nb_dom_H))
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 ')
570 DO k = 1, inputs%nb_dom_ns
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 ')
575 DO kp = 1, inputs%nb_dom_H
576 IF (inputs%list_dom_H(kp) == inputs%list_dom_ns(k))
EXIT
580 list_dom_h(k) = inputs%list_dom_ns(k)
583 DO k = 1, inputs%nb_dom_H
584 IF (minval(abs(inputs%list_dom_H(k) - inputs%list_dom_ns)) == 0) cycle
589 list_dom_h(m) = inputs%list_dom_H(k)
591 IF (m/=inputs%nb_dom_H)
THEN
596 DO k = 1, inputs%nb_dom_H
600 list_dom_h = inputs%list_dom_H
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 ')
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 ')
614 DO kp = 1, inputs%nb_dom_temp
615 IF (inputs%list_dom_temp(kp) == inputs%list_dom_ns(k))
EXIT
617 temp_in_to_new(k) = kp
619 list_dom_temp(k) = inputs%list_dom_ns(k)
622 DO k = 1, inputs%nb_dom_temp
623 IF (minval(abs(inputs%list_dom_temp(k) - inputs%list_dom_ns)) == 0) cycle
626 temp_in_to_new(m) = k
628 list_dom_temp(m) = inputs%list_dom_temp(k)
630 IF (m/=inputs%nb_dom_temp)
THEN
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 ')
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 ')
648 IF (if_momentum .AND. (.NOT. if_induction))
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
656 nsize =
SIZE(inputs%list_dom_ns)
657 ALLOCATE(list_dom(nsize))
658 list_dom = inputs%list_dom_ns
659 ALLOCATE(list_inter(0))
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
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
673 IF (
SIZE(inputs%list_inter_H_phi)>0)
THEN
674 list_inter(
SIZE(inputs%list_inter_mu)+1:) = inputs%list_inter_H_phi
679 ALLOCATE(list_inter_temp(0))
684 list_inter, 1, p1_mesh_glob, inputs%iformatted)
686 list_inter, 2, p2_mesh_glob, inputs%iformatted)
689 list_inter_temp, 2, p2_c0_mesh_glob_temp, inputs%iformatted)
691 IF (if_induction)
THEN
692 ALLOCATE(list_dummy(0))
694 inputs%list_inter_H_phi, 1, p1_c0_mesh_glob, inputs%iformatted)
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')
705 WRITE(*,*)
'read partition'
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')
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))
724 CALL mpi_comm_dup(comm_one_d(2), comm_one_d_ns(2), code)
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)
730 CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_ns(1),code)
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')
746 IF (ns_periodic)
THEN
751 IF (inputs%if_level_set_P2)
THEN
752 CALL
prep_periodic(inputs%my_periodic, vv_mesh, level_set_per)
754 CALL
prep_periodic(inputs%my_periodic, pp_mesh, level_set_per)
759 IF (pp_mesh%me/=0)
THEN
767 IF (inputs%is_mesh_symmetric)
THEN
768 ALLOCATE(vv_mz_la(vv_mesh%np))
777 CALL
gen_gauss(vv_mesh,edge_stab=.false.)
778 CALL
gen_gauss(pp_mesh,edge_stab=.false.)
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)
788 CALL
extract_mesh(comm_one_d(1),nb_procs_s,p2_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
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)
793 CALL
extract_mesh(comm_one_d(1),nb_procs_s,p2_mesh_glob,part,inputs%list_dom_phi,phi_mesh_glob,phi_mesh)
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)
806 CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_temp(1),code)
813 IF (if_induction)
THEN
814 DEALLOCATE(list_dummy)
818 DEALLOCATE(list_inter_temp)
823 IF (if_induction)
THEN
826 IF (pmag_mesh%me/=0)
THEN
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))))
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')
839 IF (if_momentum .OR. inputs%type_pb==
'mxx')
THEN
840 IF (vv_mesh%me/=0)
THEN
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))))
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')
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,:)))
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))))
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')
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.)
886 IF (h_mesh%me /=0)
THEN
887 CALL
load_interface(h_mesh, h_mesh, inputs%list_inter_mu, interface_h_mu, .false.)
889 interface_h_mu%mes = 0
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.)
894 interface_h_phi%mes = 0
898 IF (mxw_periodic)
THEN
900 WRITE(*,*)
'periodic MHD done'
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)
911 IF (inputs%is_mesh_symmetric)
THEN
912 ALLOCATE(h_mz_la(h_mesh%np))
925 CALL
gen_gauss(pmag_mesh,edge_stab=.false.)
926 CALL
gen_gauss(phi_mesh,edge_stab=.false.)
929 ALLOCATE(sigma_field(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))
939 ALLOCATE(mu_h_field(h_mesh%np))
940 IF (inputs%analytical_permeability)
THEN
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))
963 IF (vv_mesh%me/=0)
THEN
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))))
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')
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,:)))
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))))
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')
1005 IF (temp_periodic)
THEN
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)
1016 CALL
gen_gauss(temp_mesh,edge_stab=.false.)
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))
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))
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 ')
1049 IF (
ALLOCATED(list_dom_h))
DEALLOCATE(list_dom_h)
1050 IF (
ASSOCIATED(inputs%list_dom_ns))
DEALLOCATE(inputs%list_dom_ns)
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 ')
1062 IF (
ALLOCATED(list_dom_temp))
DEALLOCATE(list_dom_temp)
1065 IF (if_momentum .OR. inputs%type_pb==
'mxx')
THEN
1070 IF (if_induction)
THEN
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))
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))
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))
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))
1124 ALLOCATE(tempn_m1(temp_mesh%np, 2, m_max_c))
1125 ALLOCATE(tempn(temp_mesh%np, 2, m_max_c))
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))
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)
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)
1145 jj_v_to_h(j2) = jj_v_to_h(j1)
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))
1157 ALLOCATE(h_to_ns(1, 1, 1))
1158 ALLOCATE(b_to_ns(1, 1, 1))
1163 ALLOCATE(v_to_energy(temp_mesh%np, 6, m_max_c))
1164 ALLOCATE(jj_v_to_temp(temp_mesh%np))
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)
1170 ALLOCATE(t_to_ns(vv_mesh%np, 2, m_max_c))
1172 ALLOCATE(t_to_ns(1, 1, 1))
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))
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)
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)
1192 jj_temp_to_h(j2) = jj_temp_to_h(j1)
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)
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)
1213 inputs%dt, list_mode, un_m1, un, pn_m1, pn, incpn_m1, incpn)
1215 IF (inputs%if_level_set_P2)
THEN
1217 inputs%dt, list_mode, level_set_m1, level_set)
1220 inputs%dt, list_mode, level_set_m1, level_set)
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)
1234 inputs%density_fluid, density_m1)
1236 inputs%density_fluid, density)
1237 visc_entro_level = 0.d0
1242 IF (list_mode(i) == 0)
THEN
1245 un_m1(:,2*k,i) = 0.d0
1249 density(:,2,i) = 0.d0
1251 incpn_m1(:,2,i) = 0.d0
1252 density_m1(:,2,i) = 0.d0
1254 level_set(:,:,2,i) = 0.d0
1255 level_set_m1(:,:,2,i) = 0.d0
1260 density_m2 = density_m1
1266 IF ( (inputs%type_pb==
'mxw') .AND. (h_mesh%me/=0) )
THEN
1268 v_to_max(:,:,i) =
vexact(list_mode(i), h_mesh)
1273 IF (inputs%type_pb==
'mxx')
THEN
1275 IF (h_mesh%np>0)
THEN
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),:,:)
1296 DEALLOCATE(incpn_m1)
1298 DEALLOCATE(density_m2)
1299 DEALLOCATE(density_m1)
1302 IF (
ALLOCATED(level_set))
DEALLOCATE(level_set,level_set_m1)
1308 IF (if_induction)
THEN
1309 IF (inputs%irestart_h)
THEN
1311 phin, phin1, inputs%file_name)
1313 CALL
init_maxwell(h_mesh,phi_mesh,time_h,inputs%dt,mu_h_field,inputs%mu_phi,list_mode,&
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
1322 h_mesh, hn, bn, nb_procs, bloc_size, m_max_pad, time)
1324 h_mesh, hn1, bn1, nb_procs, bloc_size, m_max_pad, time)
1328 bn(:,k,i) = mu_h_field*hn(:,k,i)
1329 bn1(:,k,i) = mu_h_field*hn1(:,k,i)
1337 IF (list_mode(i) == 0)
THEN
1338 IF (h_mesh%me/=0)
THEN
1344 IF (phi_mesh%me/=0)
THEN
1356 IF (temp_mesh%me/=0)
THEN
1357 IF (inputs%irestart_T)
THEN
1361 CALL
init_temperature(temp_mesh, time_t, inputs%dt, list_mode, tempn_m1, tempn)
1365 IF (list_mode(i) == 0)
THEN
1367 tempn_m1(:,2,i) = 0.d0
1375 IF (inputs%irestart_h .OR. inputs%irestart_u .OR. inputs%irestart_T)
THEN
1376 time = max(time_u,time_h,time_t)
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')
1394 REAL(KIND=8),
DIMENSION(ndim) :: vect_in
1395 REAL(KIND=8),
DIMENSION(ndim) :: vect_out
1397 INTEGER :: type, i_deb, i_fin
1398 REAL(KIND=8),
DIMENSION(vv_mesh%np,2,SIZE(list_mode)) :: sigma_ns
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
1407 hn(:,type,i) = vect_in(i_deb:i_fin)
1411 phin(:,type,i) =0.d0
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
1420 hn1(:,type,i) = vect_in(i_deb:i_fin)
1424 phin1(:,type,i) =0.d0
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)
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)
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)
1449 CHARACTER(len=200),
DIMENSION(:,:),
ALLOCATABLE :: inline
1451 CHARACTER(len=3) :: tit
1457 CALL getarg(narg+1,tit)
1458 IF (tit ==
' ')
THEN
1466 ALLOCATE(inline(2,narg))
1469 CALL getarg(2*(i-1)+1,inline(1,i))
1470 CALL getarg(2*i ,inline(2,i))
1473 inputs%test_de_convergence = .false.
1474 inputs%numero_du_test_debug = 0
1475 inputs%data_directory_debug =
'.'
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)
1490 INTEGER :: m, type_fe, index, l
1491 IF (mesh%gauss%n_w==3)
THEN
1496 ALLOCATE(mesh%hloc(mesh%dom_me))
1497 ALLOCATE(mesh%hloc_gauss(mesh%gauss%l_G*mesh%dom_me))
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
1504 mesh%hloc_gauss(index) = mesh%hloc(m)
1510 REAL(KIND=8) :: h_min, h_min_f
1513 IF (inputs%if_level_set_P2)
THEN
1514 h_min_f=minval(vv_mesh%hloc_gauss)
1516 h_min_f=minval(pp_mesh%hloc_gauss)
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
subroutine, public save_run(it, freq_restart)
subroutine, public reconstruct_variable(comm_one_d, list_mode, mesh_P1, mesh_P2, level_set, values, variable)
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)
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)
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)
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)
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)
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)
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 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)
subroutine assign_boundary
subroutine, public st_aij_csr_glob_block(communicator, kmax, mesh_glob, mesh, LA, opt_per)
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)
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)
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)
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)
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)
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
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)