14 temp_mesh, tempn_m1, tempn, chmp_vit, chmp_mag, chmp_pdt_h, vol_heat_capacity, &
15 temp_diffusivity, my_par_cc, temp_list_dirichlet_sides, &
16 temp_list_robin_sides, convection_coeff, exterior_temperature, temp_per)
32 #include "petsc/finclude/petsc.h" 35 REAL(KIND=8) :: time, dt
36 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
41 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: tempn_m1, tempn
42 INTEGER,
DIMENSION(:),
INTENT(IN) :: temp_list_dirichlet_sides
43 INTEGER,
DIMENSION(:),
INTENT(IN) :: temp_list_robin_sides
44 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vol_heat_capacity, temp_diffusivity
45 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: convection_coeff, exterior_temperature
46 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: chmp_vit, chmp_mag, chmp_pdt_H
47 LOGICAL,
SAVE :: once = .true.
48 INTEGER,
SAVE :: m_max_c
49 INTEGER,
DIMENSION(:),
POINTER,
SAVE :: temp_js_D
50 INTEGER,
SAVE :: my_petscworld_rank
51 REAL(KIND=8),
SAVE :: mass0, hmoy
52 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: temp_global_D
53 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: temp_mode_global_js_D
57 INTEGER,
POINTER,
DIMENSION(:) :: temp_1_ifrom
58 INTEGER :: i, m, n, l, index
61 REAL(KIND=8),
DIMENSION(temp_mesh%np) :: ff
62 REAL(KIND=8),
DIMENSION(temp_mesh%np, 2) :: tempn_p1
63 REAL(KIND=8),
DIMENSION(temp_mesh%gauss%l_G*temp_mesh%me,2, SIZE(list_mode)) :: ff_conv, pyromag_term
64 REAL(KIND=8) ::tps, tps_tot, tps_cumul
65 REAL(KIND=8) :: one, zero, three
66 DATA zero, one, three/0.d0,1.d0,3.d0/
67 REAL(KIND=8),
DIMENSION(2,temp_mesh%gauss%l_G*temp_mesh%me) :: rr_gauss
68 INTEGER,
DIMENSION(temp_mesh%gauss%n_w) :: j_loc
70 petscerrorcode :: ierr
71 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
72 mat,
DIMENSION(:),
POINTER,
SAVE :: temp_mat
73 vec,
SAVE :: cb_1, cb_2, cx_1, cx_1_ghost
74 ksp,
DIMENSION(:),
POINTER,
SAVE :: temp_ksp
81 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
86 CALL veccreateghost(comm_one_d(1), n, &
87 petsc_determine,
SIZE(temp_1_ifrom), temp_1_ifrom, cx_1, ierr)
88 CALL vecghostgetlocalform(cx_1, cx_1_ghost, ierr)
89 CALL vecduplicate(cx_1, cb_1, ierr)
90 CALL vecduplicate(cx_1, cb_2, ierr)
94 m_max_c =
SIZE(list_mode)
100 ALLOCATE(temp_global_d(m_max_c))
102 ALLOCATE(temp_global_d(i)%DRL(
SIZE(temp_mode_global_js_d(i)%DIL)))
108 DO m = 1, temp_mesh%dom_me
109 hmoy = hmoy + sqrt(sum(temp_mesh%gauss%rj(:,m)))/2
111 hmoy = hmoy/temp_mesh%dom_me
116 CALL mass_tot(comm_one_d(1),temp_mesh, tempn(:,1,i), mass0)
122 ALLOCATE(temp_mat(m_max_c),temp_ksp(m_max_c))
130 1.5d0/dt, temp_list_robin_sides, convection_coeff, zero, mode, temp_mat(i))
132 IF (temp_per%n_bord/=0)
THEN 138 CALL init_solver(my_par_cc,temp_ksp(i),temp_mat(i),comm_one_d(1),&
139 solver=my_par_cc%solver,precond=my_par_cc%precond)
148 CALL smb_ugradc_gauss_fft_par(comm_one_d(2),temp_mesh,list_mode,vol_heat_capacity,chmp_vit,2*tempn-tempn_m1,ff_conv)
149 IF (
inputs%type_pb==
'fhd')
THEN 151 CALL smb_pyromag_gauss_fft_par(comm_one_d(2),temp_mesh,list_mode,2*tempn-tempn_m1,chmp_vit,chmp_mag,chmp_pdt_h,pyromag_term)
152 ff_conv = ff_conv + pyromag_term
154 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
158 DO m = 1, temp_mesh%me
159 j_loc = temp_mesh%jj(:,m)
160 DO l = 1, temp_mesh%gauss%l_G
162 rr_gauss(1,index) = sum(temp_mesh%rr(1,j_loc)*temp_mesh%gauss%ww(:,l))
163 rr_gauss(2,index) = sum(temp_mesh%rr(2,j_loc)*temp_mesh%gauss%ww(:,l))
171 ff = (2d0/dt)*tempn(:,1,i) - (1d0/(2*dt))*tempn_m1(:,1,i)
172 CALL qs_00_gauss (temp_mesh, temp_1_la, vol_heat_capacity, ff, &
173 -ff_conv(:,1,i) + source_in_temperature(1, rr_gauss, mode, time), cb_1)
175 ff = (2d0/dt)*tempn(:,2,i) - (1d0/(2*dt))*tempn_m1(:,2,i)
176 CALL qs_00_gauss (temp_mesh, temp_1_la, vol_heat_capacity, ff, &
177 -ff_conv(:,2,i) + source_in_temperature(2, rr_gauss, mode, time), cb_2)
181 CALL qs_00_gauss_surface(temp_mesh, temp_1_la, temp_list_robin_sides, convection_coeff, exterior_temperature, cb_1)
185 IF (temp_per%n_bord/=0)
THEN 186 CALL periodic_rhs_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, cb_1, temp_1_la)
187 CALL periodic_rhs_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, cb_2, temp_1_la)
192 temp_global_d(i)%DRL(n+1:) = 0.d0
193 temp_global_d(i)%DRL(1:n) = temperature_exact(1,temp_mesh%rr(:,temp_js_d), mode, time)
194 CALL dirichlet_rhs(temp_mode_global_js_d(i)%DIL-1,temp_global_d(i)%DRL,cb_1)
195 temp_global_d(i)%DRL(1:n) = temperature_exact(2,temp_mesh%rr(:,temp_js_d), mode, time)
196 CALL dirichlet_rhs(temp_mode_global_js_d(i)%DIL-1,temp_global_d(i)%DRL,cb_2)
198 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
205 CALL solver(temp_ksp(i),cb_1,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
206 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
207 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
208 CALL extract(cx_1_ghost,1,1,temp_1_la,tempn_p1(:,1))
211 CALL solver(temp_ksp(i),cb_2,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
212 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
213 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
214 CALL extract(cx_1_ghost,1,1,temp_1_la,tempn_p1(:,2))
215 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
229 tempn_m1(:,:,i) = tempn(:,:,i)
230 tempn(:,:,i) = tempn_p1
245 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
263 #include "petsc/finclude/petsc.h" 267 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
268 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: heat_capa_in
269 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: V_in, c_in
270 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: c_out
271 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: Gradc, W
272 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: Div, Cgauss, cint
273 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
274 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
275 INTEGER :: m, l , i, mode, index, k
276 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: Vs
277 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: cs
278 INTEGER,
DIMENSION(:,:),
POINTER :: jj
279 INTEGER,
POINTER :: me
280 REAL(KIND=8) :: ray, tps
281 REAL(KIND=8),
DIMENSION(3) :: temps
282 INTEGER :: code, m_max_pad, bloc_size, nb_procs
283 mpi_comm :: communicator
290 DO i = 1,
SIZE(list_mode)
293 DO m = 1, mesh%dom_me
296 vs(:,k) = v_in(j_loc,k,i)
299 cs(:,k) = c_in(j_loc,k,i)
306 ray = sum(mesh%rr(1,j_loc)*
ww(:,l))
309 w(index,1,i) = sum(vs(:,1)*
ww(:,l))
310 w(index,3,i) = sum(vs(:,3)*
ww(:,l))
311 w(index,5,i) = sum(vs(:,5)*
ww(:,l))
313 w(index,2,i) = sum(vs(:,2)*
ww(:,l))
314 w(index,4,i) = sum(vs(:,4)*
ww(:,l))
315 w(index,6,i) = sum(vs(:,6)*
ww(:,l))
317 div(index,1,i) = sum(vs(:,1)*dw_loc(1,:)) + sum(vs(:,1)*
ww(:,l))/ray &
318 + (mode/ray)*sum(vs(:,4)*
ww(:,l)) + sum(vs(:,5)*dw_loc(2,:))
319 div(index,2,i) = sum(vs(:,2)*dw_loc(1,:)) + sum(vs(:,2)*
ww(:,l))/ray &
320 - (mode/ray)*sum(vs(:,3)*
ww(:,l)) + sum(vs(:,6)*dw_loc(2,:))
324 gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
325 gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
326 gradc(index,3,i) = mode/ray*sum(cs(:,2)*
ww(:,l))
327 gradc(index,4,i) = -mode/ray*sum(cs(:,1)*
ww(:,l))
328 gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
329 gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))
331 gradc(index,:,i) = heat_capa_in(m) * gradc(index,:,i)
333 cgauss(index,1,i) = sum(cs(:,1)*
ww(:,l))
334 cgauss(index,2,i) = sum(cs(:,2)*
ww(:,l))
336 cgauss(index,:,i) = heat_capa_in(m) * cgauss(index,:,i)
347 CALL mpi_comm_size(communicator, nb_procs, code)
348 bloc_size =
SIZE(gradc,1)/nb_procs+1
349 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
351 bloc_size =
SIZE(div,1)/nb_procs+1
352 CALL fft_par_prod_dcl(communicator, div, cgauss, cint, nb_procs, bloc_size, m_max_pad, temps)
362 SUBROUTINE mass_tot(communicator,mesh,tempn,RESLT)
366 #include "petsc/finclude/petsc.h" 370 REAL(KIND=8),
DIMENSION(:) ,
INTENT(IN) :: tempn
371 REAL(KIND=8) ,
INTENT(OUT) :: RESLT
372 REAL(KIND=8) :: r_loc, r_out
373 INTEGER :: m, l , i , ni, code
374 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
376 mpi_comm :: communicator
379 DO m = 1, mesh%dom_me
381 DO l = 1, mesh%gauss%l_G
383 DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
384 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
387 r_loc = r_loc + sum(tempn(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
392 CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
402 #include "petsc/finclude/petsc.h" 406 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
407 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: scal_in, vect_1_in, vect_2_in, vect_3_in
408 REAL(KIND=8),
DIMENSION(:,:,:) :: scal_out
409 REAL(KIND=8),
DIMENSION(mesh%np,2,SIZE(list_mode)) :: T_dchi_dT_coeff, vect_2_in_square, v2_dot_v3
410 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: T_dchi_dT_coeff_gauss, pdtH2, ugradH2, DH2_Dt
411 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: vect_1_in_gauss, grad_vect_2_in_square
412 INTEGER :: i, mode, index, m, k, l
413 INTEGER,
DIMENSION(:,:),
POINTER :: jj
414 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
417 INTEGER :: nb_procs, bloc_size, m_max_pad, code
418 mpi_comm :: communicator
423 CALL mpi_comm_size(communicator, nb_procs, code)
424 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
427 t_dchi_dt_coeff = scal_in
428 bloc_size =
SIZE(t_dchi_dt_coeff,1)/nb_procs+1
429 CALL fft_par_scal_funct(communicator, t_dchi_dt_coeff, t_dchi_dt_coeff_law, nb_procs, bloc_size, m_max_pad)
432 DO i = 1,
SIZE(list_mode)
434 DO m = 1, mesh%dom_me
435 IF (minval(abs(mesh%i_d(m) -
inputs%list_dom_ns)) == 0)
THEN 439 t_dchi_dt_coeff_gauss(index,k,i) = sum(t_dchi_dt_coeff(jj(:,m),k,i)*
ww(:,l))
443 t_dchi_dt_coeff_gauss(index+1:index+
l_g,:,i) = 0.d0
450 IF (
inputs%if_steady_current_fhd)
THEN 453 bloc_size =
SIZE(vect_2_in,1)/nb_procs+1
454 CALL fft_par_dot_prod_dcl(communicator, vect_2_in, vect_3_in, v2_dot_v3, nb_procs, bloc_size, m_max_pad)
455 DO i = 1,
SIZE(list_mode)
457 DO m = 1, mesh%dom_me
458 IF (minval(abs(mesh%i_d(m) -
inputs%list_dom_ns)) == 0)
THEN 462 pdth2(index,k,i) = 2*sum(v2_dot_v3(jj(:,m),k,i)*
ww(:,l))
466 pdth2(index+1:index+
l_g,:,i) = 0.d0
474 DO i = 1,
SIZE(list_mode)
476 DO m = 1, mesh%dom_me
477 IF (minval(abs(mesh%i_d(m) -
inputs%list_dom_ns)) == 0)
THEN 481 vect_1_in_gauss(index,k,i) = sum(vect_1_in(jj(:,m),k,i)*
ww(:,l))
485 vect_1_in_gauss(index+1:index+
l_g,:,i) = 0.d0
492 bloc_size =
SIZE(vect_2_in,1)/nb_procs+1
493 CALL fft_par_dot_prod_dcl(communicator, vect_2_in, vect_2_in, vect_2_in_square, nb_procs, bloc_size, m_max_pad)
494 DO i = 1,
SIZE(list_mode)
497 DO m = 1, mesh%dom_me
498 IF (minval(abs(mesh%i_d(m) -
inputs%list_dom_ns)) == 0)
THEN 502 rad = sum(mesh%rr(1,jj(:,m))*
ww(:,l))
503 grad_vect_2_in_square(index,1,i) = sum(vect_2_in_square(jj(:,m),1,i)*dw_loc(1,:))
504 grad_vect_2_in_square(index,2,i) = sum(vect_2_in_square(jj(:,m),2,i)*dw_loc(1,:))
505 grad_vect_2_in_square(index,3,i) = mode/rad * sum(vect_2_in_square(jj(:,m),2,i)*
ww(:,l))
506 grad_vect_2_in_square(index,4,i) = - mode/rad * sum(vect_2_in_square(jj(:,m),1,i)*
ww(:,l))
507 grad_vect_2_in_square(index,5,i) = sum(vect_2_in_square(jj(:,m),1,i)*dw_loc(2,:))
508 grad_vect_2_in_square(index,6,i) = sum(vect_2_in_square(jj(:,m),2,i)*dw_loc(2,:))
511 grad_vect_2_in_square(index+1:index+
l_g,:,i) = 0.d0
518 bloc_size =
SIZE(vect_1_in_gauss,1)/nb_procs+1
519 CALL fft_par_dot_prod_dcl(communicator, vect_1_in_gauss, grad_vect_2_in_square, ugradh2, nb_procs, bloc_size, m_max_pad)
522 dh2_dt = pdth2 + ugradh2
525 bloc_size =
SIZE(t_dchi_dt_coeff_gauss,1)/nb_procs+1
526 CALL fft_par_prod_dcl(communicator, t_dchi_dt_coeff_gauss, 0.5*dh2_dt, scal_out, nb_procs, bloc_size, m_max_pad)
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine, public three_level_temperature(comm_one_d, time, temp_1_LA, dt, list_mode, temp_mesh, tempn_m1, tempn, chmp_vit, chmp_mag, chmp_pdt_H, vol_heat_capacity, temp_diffusivity, my_par_cc, temp_list_dirichlet_sides, temp_list_robin_sides, convection_coeff, exterior_temperature, temp_per)
subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, temp_list_robin_sides, convection_coeff, stab, mode, matrix)
subroutine, public fft_par_scal_funct(communicator, c1_inout, funct, nb_procs, bloc_size, m_max_pad, temps)
subroutine qs_00_gauss_surface(mesh, vv_1_LA, temp_list_robin_sides, convection_coeff, exterior_temperature, cb)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
real(kind=8) function user_time()
subroutine smb_pyromag_gauss_fft_par(communicator, mesh, list_mode, scal_in, vect_1_in, vect_2_in, vect_3_in, scal_out)
subroutine, public fft_par_prod_dcl(communicator, c1_in, c2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine scalar_with_bc_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_js_D, pp_mode_global_js_D)
subroutine qs_00_gauss(mesh, LA, heat_capa, ff, ff_gauss, vect)
subroutine smb_ugradc_gauss_fft_par(communicator, mesh, list_mode, heat_capa_in, V_in, c_in, c_out)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine dirichlet_rhs(js_D, bs_D, b)
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine mass_tot(communicator, mesh, tempn, RESLT)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
subroutine, public fft_par_dot_prod_dcl(communicator, V1_in, V2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
real(kind=8), dimension(:,:,:,:), pointer dw
real(kind=8), dimension(:,:), pointer ww
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)