13 TYPE(mesh_type),
POINTER :: pp_mesh, vv_mesh
14 REAL(KIND=8),
POINTER,
DIMENSION(:,:,:) :: un, pn
16 TYPE(mesh_type),
POINTER :: h_mesh, phi_mesh
18 REAL(KIND=8),
POINTER,
DIMENSION(:,:,:) :: hn, bn, phin, vel
19 REAL(KIND=8),
POINTER,
DIMENSION(:) :: sigma_field, mu_h_field
22 REAL(KIND=8),
POINTER,
DIMENSION(:,:,:) :: temperature
23 REAL(KIND=8),
POINTER,
DIMENSION(:) :: vol_heat_capacity_field
24 REAL(KIND=8),
POINTER,
DIMENSION(:) :: temperature_diffusivity_field
26 REAL(KIND=8),
POINTER,
DIMENSION(:,:,:,:) :: level_set
28 REAL(KIND=8),
POINTER,
DIMENSION(:,:,:) :: density
31 INTEGER,
POINTER,
DIMENSION(:) :: list_mode
36 REAL(KIND=8) :: tps, tploc, tploc_max=0.d0
38 #include "petsc/finclude/petsc.h"
39 petscerrorcode :: ierr
41 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d, comm_one_d_ns, comm_one_d_temp
44 CALL petscinitialize(petsc_null_character,ierr)
45 CALL mpi_comm_rank(petsc_comm_world,rank,ierr)
51 CALL
initial(vv_mesh, pp_mesh, h_mesh, phi_mesh, temp_mesh, &
52 interface_h_phi, interface_h_mu, list_mode, &
53 un, pn, hn, bn, phin, vel, &
54 vol_heat_capacity_field, temperature_diffusivity_field, mu_h_field, sigma_field, &
55 time, m_max_c, comm_one_d, comm_one_d_ns, comm_one_d_temp, temperature, level_set, density)
60 IF (inputs%if_just_processing)
THEN
69 IF (inputs%if_arpack)
THEN
72 inputs%dt,list_mode,mu_h_field)
74 IF (inputs%test_de_convergence)
THEN
75 CALL
post_proc_test(vv_mesh, pp_mesh, temp_mesh, h_mesh, phi_mesh, list_mode, &
76 un, pn, hn, bn, phin, temperature, level_set, mu_h_field, &
77 time, m_max_c, comm_one_d, comm_one_d_ns)
95 DO it = 1, inputs%nb_iteration
97 time = time + inputs%dt
102 IF (.NOT.inputs%test_de_convergence)
THEN
107 IF (mod(it, inputs%freq_restart) == 0)
THEN
108 CALL
save_run(it,inputs%freq_restart)
113 IF (it>1) tploc_max = tploc_max + tploc
122 IF (inputs%test_de_convergence)
THEN
123 CALL
post_proc_test(vv_mesh, pp_mesh, temp_mesh, h_mesh, phi_mesh, list_mode, &
124 un, pn, hn, bn, phin, temperature, level_set, mu_h_field, &
125 time, m_max_c, comm_one_d, comm_one_d_ns)
141 INTEGER,
INTENT(IN) :: it
142 REAL(KIND=8) :: err, norm
143 INTEGER :: i, it_plot
144 CHARACTER(LEN=3) :: what
145 INTEGER :: rank_s, rank_f
146 INTEGER :: rank_ns_s, rank_ns_f
147 REAL(KIND=8),
DIMENSION(vv_mesh%np, 2, SIZE(list_mode)) :: level_1_p2
148 REAL(KIND=8),
DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: level_1_p1
155 IF (vv_mesh%me /=0)
THEN
156 CALL mpi_comm_rank(comm_one_d_ns(1), rank_ns_s, ierr)
157 CALL mpi_comm_rank(comm_one_d_ns(2), rank_ns_f, ierr)
162 CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
163 CALL mpi_comm_rank(comm_one_d(2), rank_f, ierr)
166 IF (inputs%check_numerical_stability)
THEN
167 IF (inputs%type_pb==
'nst' .OR. inputs%type_pb==
'mhd' .OR. inputs%type_pb==
'fhd')
THEN
168 norm =
norm_sf(comm_one_d_ns,
'L2', vv_mesh, list_mode, un)
170 norm =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, hn)
173 CALL
error_petsc(
'From my_post_processing: numerical unstability')
178 IF (mod(it,inputs%freq_en) == 0)
THEN
183 IF (inputs%type_pb==
'nst' .OR. inputs%type_pb==
'mhd' .OR. inputs%type_pb==
'fhd')
THEN
184 IF (inputs%if_compute_momentum_pseudo_force)
THEN
191 err =
norm_sf(comm_one_d,
'div', vv_mesh, list_mode, un)
192 norm =
norm_sf(comm_one_d,
'H1', vv_mesh, list_mode, un)
195 WRITE(31,*) time, err/norm
197 DO i=1,
SIZE(list_mode)
198 norm =
norm_s(comm_one_d,
'L2', vv_mesh, list_mode(i:i), un(:,:,i:i))
199 IF (rank_ns_s == 0)
THEN
201 WRITE(100+list_mode(i),*) time, norm
205 err =
norm_sf(comm_one_d,
'L2', vv_mesh, list_mode, un)
206 norm =
norm_sf(comm_one_d,
'sH1', vv_mesh, list_mode, un)
208 WRITE(98,*) time, err
209 WRITE(*,*)
'norm L2 of velocity', time, err
210 WRITE(*,*)
'semi norm H1 of velocity', time, norm
213 err =
norm_sf(comm_one_d,
'L2', pp_mesh, list_mode, pn)
215 WRITE(*,*)
'norm L2 of pressure', time, err
218 IF (inputs%if_level_set)
THEN
221 IF (inputs%if_level_set_P2)
THEN
228 IF (inputs%if_temperature)
THEN
229 norm =
norm_sf(comm_one_d,
'L2', temp_mesh, list_mode, temperature)
231 WRITE(99,*)
'norm L2 of temperature', time, norm
232 WRITE(*,*)
'norm L2 of temperature', time, norm
237 IF (inputs%type_pb/=
'nst')
THEN
238 err =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, hn)
241 WRITE(41,*) time, err
243 err =
norm_sf(comm_one_d,
'div', h_mesh, list_mode, bn)
244 norm =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, bn)
247 WRITE(51,*) time, err, err/norm
248 WRITE(52,*) time, err, norm
249 WRITE(*,*)
'norm L2 of magnetic field', time, norm
251 DO i=1,
SIZE(list_mode)
252 norm =
norm_s(comm_one_d,
'L2', h_mesh, list_mode(i:i), hn(:,:,i:i))
253 IF (rank_s == 0)
THEN
255 WRITE(200+list_mode(i),*) time, norm
261 IF (mod(it,inputs%freq_plot) == 0)
THEN
263 IF (it==inputs%freq_plot)
THEN
268 it_plot = it/inputs%freq_plot
270 IF (inputs%if_level_set)
THEN
271 IF (inputs%if_level_set_P2)
THEN
272 level_1_p2=level_set(1,:,:,:)
273 CALL
vtu_3d(level_1_p2,
'vv_mesh',
'Level_1',
'level_1', what, opt_it=it_plot)
275 level_1_p1=level_set(1,:,:,:)
276 CALL
vtu_3d(level_1_p1,
'pp_mesh',
'Level_1',
'level_1', what, opt_it=it_plot)
278 CALL
vtu_3d(density,
'vv_mesh',
'Density',
'density', what, opt_it=it_plot)
280 IF (inputs%type_pb/=
'mxw' .AND. inputs%type_pb/=
'mxx')
THEN
281 CALL
vtu_3d(un,
'vv_mesh',
'Velocity',
'vel', what, opt_it=it_plot)
282 CALL
vtu_3d(pn,
'pp_mesh',
'Pressure',
'pre', what, opt_it=it_plot)
284 IF (inputs%if_temperature)
THEN
285 CALL
vtu_3d(temperature,
'temp_mesh',
'Temperature',
'temp', what, opt_it=it_plot)
287 IF (inputs%type_pb/=
'nst')
THEN
288 CALL
vtu_3d(hn,
'H_mesh',
'MagField',
'mag', what, opt_it=it_plot)
289 IF (inputs%nb_dom_phi>0)
THEN
290 CALL
vtu_3d(phin,
'phi_mesh',
'ScalPot',
'phi', what, opt_it=it_plot)
370 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
371 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: un
372 REAL(KIND=8),
INTENT(IN) :: time
373 REAL(KIND=8),
DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: vel_gauss, vel_gauss_penal
374 REAL(KIND=8),
DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
375 INTEGER,
DIMENSION(vv_mesh%gauss%n_w) :: j_loc
376 REAL(KIND=8) :: vel_torque, vel_torque_tot
377 INTEGER :: m, l , i, mode, index, type, nb_procs, m_max_pad, bloc_size
378 #include "petsc/finclude/petsc.h"
379 petscerrorcode :: ierr
380 mpi_comm,
DIMENSION(2) :: communicator
383 DO m = 1, vv_mesh%dom_me
384 j_loc = vv_mesh%jj(:,m)
385 DO l = 1, vv_mesh%gauss%l_G
387 rr_gauss(1,index) = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
388 rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
392 DO i = 1,
SIZE(list_mode)
395 DO m = 1, vv_mesh%dom_me
396 j_loc = vv_mesh%jj(:,m)
397 DO l = 1, vv_mesh%gauss%l_G
400 vel_gauss(index,type,i) = sum(un(j_loc,type,i)*vv_mesh%gauss%ww(:,l))*(3/(2*inputs%dt))
404 IF(inputs%if_impose_vel_in_solids)
THEN
411 CALL mpi_comm_size(communicator(2), nb_procs, ierr)
412 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
413 bloc_size =
SIZE(vel_gauss,1)/nb_procs+1
415 vel_gauss, vel_gauss_penal, nb_procs, bloc_size, m_max_pad, rr_gauss, time)
418 DO i = 1,
SIZE(list_mode)
424 DO m = 1, vv_mesh%dom_me
425 j_loc = vv_mesh%jj(:,m)
426 DO l = 1, vv_mesh%gauss%l_G
429 vel_torque = vel_torque + (vel_gauss_penal(index,5,i) - vel_gauss(index,5,i)) &
430 *rr_gauss(1,index)*vv_mesh%gauss%rj(l,m)
433 CALL mpi_allreduce(vel_torque, vel_torque_tot,1,mpi_double_precision, mpi_sum, communicator(1), ierr)
434 WRITE(*,*)
' FORCES_AND_MOMENTS ', time, 2*acos(-1.d0)*vel_torque_tot/(0.5d0*acos(-1.d0))
435 WRITE(12,*) time, 2*acos(-1.d0)*vel_torque_tot/(0.5d0*acos(-1.d0))
447 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
448 REAL(KIND=8),
DIMENSION(:,:,:,:),
INTENT(IN) :: level_set
449 REAL(KIND=8),
INTENT(IN) :: time
450 LOGICAL,
SAVE :: once_compute=.true.
451 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:),
SAVE :: volum_init
452 REAL(KIND=8) :: volum_init_loc, volum_init_f
453 REAL(KIND=8) :: inte_fft_loc, inte_fft_tot_f
454 REAL(KIND=8),
DIMENSION(inputs%nb_fluid-1) :: inte_fft_tot
455 REAL(KIND=8),
DIMENSION(mesh%np, 2, SIZE(list_mode)) :: level_posi_fft
457 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
458 INTEGER :: m, l , i, nb_inter
459 INTEGER :: my_petscworld_rank, code
460 #include "petsc/finclude/petsc.h"
461 mpi_comm,
DIMENSION(2) :: communicator
463 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
465 103
FORMAT(1500(e22.9,2x))
468 IF (once_compute)
THEN
469 once_compute = .false.
471 ALLOCATE(volum_init(
SIZE(level_set,1)))
473 DO nb_inter=1,
SIZE(level_set,1)
474 volum_init_loc = 0.d0
475 DO i = 1,
SIZE(list_mode)
476 IF (list_mode(i)==0)
THEN
479 DO l = 1, mesh%gauss%l_G
481 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
482 volum_init_loc = volum_init_loc + sum(
level_set_exact(nb_inter,1,mesh%rr(:,j_loc),list_mode(i),0.d0)* &
483 mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
488 volum_init_loc = volum_init_loc*2*acos(-1.d0)
489 CALL mpi_allreduce(volum_init_loc, volum_init_f, 1, mpi_double_precision, mpi_sum, &
490 communicator(2), code)
491 CALL mpi_allreduce(volum_init_f, volum_init(nb_inter), 1, mpi_double_precision, mpi_sum, &
492 communicator(1), code)
494 IF (my_petscworld_rank==0)
THEN
495 WRITE(*,*)
'mass initial = ', time, volum_init
500 DO nb_inter=1,
SIZE(level_set,1)
501 level_posi_fft = level_set(nb_inter,:,:,:)
503 DO i = 1,
SIZE(list_mode)
504 IF (list_mode(i)==0)
THEN
507 DO l = 1, mesh%gauss%l_G
509 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
510 inte_fft_loc = inte_fft_loc + sum(level_posi_fft(j_loc,1,i)*mesh%gauss%ww(:,l))* &
511 ray*mesh%gauss%rj(l,m)
516 inte_fft_loc = inte_fft_loc*2*acos(-1.d0)
517 CALL mpi_allreduce(inte_fft_loc, inte_fft_tot_f, 1, mpi_double_precision, mpi_sum, &
518 communicator(2), code)
519 CALL mpi_allreduce(inte_fft_tot_f, inte_fft_tot(nb_inter), 1, mpi_double_precision, mpi_sum, &
520 communicator(1), code)
522 IF (my_petscworld_rank==0)
THEN
523 WRITE(*,*)
'relative mass error of level set at t = ', &
524 time, abs(1-inte_fft_tot/(volum_init+1.d-14))
525 WRITE(97,103) time, abs(1-inte_fft_tot/(volum_init+1.d-14))
subroutine, public save_run(it, freq_restart)
subroutine, public vtu_3d(field_in, name_of_mesh, header, name_of_field, what, opt_it, opt_grad_curl)
subroutine, public solver_arpack_mhd(communicator, H_mesh, phi_mesh, dt, list_mode, mu_H_field)
subroutine, public write_verbose(rank, opt_tps, opt_tploc_max)
subroutine, public post_proc_test(vv_mesh, pp_mesh, temp_mesh, H_mesh, phi_mesh, list_mode, un, pn, Hn, Bn, phin, tempn, level_setn, mu_H_field, time, m_max_c, comm_one_d, comm_one_d_ns)
real(kind=8) function norm_s(communicator, norm_type, mesh, list_mode, v)
subroutine compute_level_set_conservation(time, mesh, communicator, list_mode, level_set)
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)
real(kind=8) function, dimension(nb_angles, ne-nb+1), public penal_in_real_space(mesh, rr_gauss, angles, nb_angles, nb, ne, time)
real(kind=8) function, dimension(size(rr, 2), 6), public imposed_velocity_by_penalty(rr, t)
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
subroutine my_post_processing(it)
subroutine, public fft_par_var_eta_prod_gauss_dcl(communicator, eta, H_mesh, c_in, c_out, nb_procs, bloc_size, m_max_pad, rr_gauss, time, temps)
real(kind=8) function user_time()
subroutine, public read_user_data(data_file)
subroutine error_petsc(string)
subroutine, public run_sfemans(time_in)
real(kind=8) function, dimension(size(rr, 2)), public level_set_exact(interface_nb, TYPE, rr, m, t)
subroutine forces_and_moments(time, vv_mesh, communicator, list_mode, un)