SFEMaNS  version 5.3
Reference documentation for SFEMaNS
sub_temperature.f90
Go to the documentation of this file.
1 !
2 !Authors Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyrights 2005
3 !Revised for PETSC, Jean-Luc Guermond, Franky Luddens, January 2011
4 !
6  USE my_util
7  USE boundary
8 
10  PRIVATE
11 CONTAINS
12 
13  SUBROUTINE three_level_temperature(comm_one_d,time, temp_1_LA, dt, list_mode, &
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) ! MODIFICATION: robin
17  !==============================
18  USE def_type_mesh
19  USE fem_m_axi
20  USE fem_rhs_axi
21  USE fem_tn_axi
22  USE dir_nodes_petsc
23  USE periodic
24  USE st_matrix
25  USE solve_petsc
26  USE dyn_line
28  USE sub_plot
29  USE st_matrix
30  USE sft_parallele
31  USE input_data
32 #include "petsc/finclude/petsc.h"
33  USE petsc
34  IMPLICIT NONE
35  REAL(KIND=8) :: time, dt
36  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
37  TYPE(mesh_type), INTENT(IN) :: temp_mesh
38  type(petsc_csr_la) :: temp_1_LA
39  TYPE(periodic_type), INTENT(IN) :: temp_per
40  TYPE(solver_param), INTENT(IN) :: my_par_cc
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 ! MODIFICATION: robin
44  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vol_heat_capacity, temp_diffusivity
45  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: convection_coeff, exterior_temperature ! MODIFICATION: robin
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 ! Dirichlet nodes
50  INTEGER, SAVE :: my_petscworld_rank
51  REAL(KIND=8), SAVE :: mass0, hmoy
52  TYPE(dyn_real_line),DIMENSION(:), ALLOCATABLE, SAVE :: temp_global_D ! MODIFICATION: axis BC
53  TYPE(dyn_int_line), DIMENSION(:), POINTER, SAVE :: temp_mode_global_js_D ! MODIFICATION: axis BC
54  !----------FIN SAVE--------------------------------------------------------------------
55 
56  !----------Declaration sans save-------------------------------------------------------
57  INTEGER, POINTER, DIMENSION(:) :: temp_1_ifrom
58  INTEGER :: i, m, n, l, index
59  INTEGER :: code, mode
60  !allocations des variables locales
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 ! MODIFICATION: pyromagnetic term for fhd
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
69  !Communicators for Petsc, in space and Fourier------------------------------
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
75  !------------------------------END OF DECLARATION--------------------------------------
76 
77  IF (once) THEN
78 
79  once = .false.
80 
81  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
82 
83  !-----CREATE PETSC VECTORS AND GHOSTS-----------------------------------------
84  CALL create_my_ghost(temp_mesh,temp_1_la,temp_1_ifrom)
85  n = temp_mesh%dom_np
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)
91  !------------------------------------------------------------------------------
92 
93  !-------------DIMENSIONS-------------------------------------------------------
94  m_max_c = SIZE(list_mode)
95  !------------------------------------------------------------------------------
96 
97  !---------PREPARE pp_js_D ARRAY FOR TEMPERATURE--------------------------------------
98  CALL dirichlet_nodes_parallel(temp_mesh, temp_list_dirichlet_sides, temp_js_d)
99  CALL scalar_with_bc_glob_js_d(temp_mesh, list_mode, temp_1_la, temp_js_d, temp_mode_global_js_d) ! MODIFICATION: axis BC
100  ALLOCATE(temp_global_d(m_max_c))
101  DO i = 1, m_max_c
102  ALLOCATE(temp_global_d(i)%DRL(SIZE(temp_mode_global_js_d(i)%DIL)))
103  END DO
104  !------------------------------------------------------------------------------
105 
106  !--------------------------------------------------------------------------------
107  hmoy = 0
108  DO m = 1, temp_mesh%dom_me
109  hmoy = hmoy + sqrt(sum(temp_mesh%gauss%rj(:,m)))/2
110  END DO
111  hmoy = hmoy/temp_mesh%dom_me
112  mass0 = 0.d0
113  DO i = 1, m_max_c
114  mode = list_mode(i)
115  IF (mode == 0) THEN
116  CALL mass_tot(comm_one_d(1),temp_mesh, tempn(:,1,i), mass0)
117  ENDIF
118  END DO
119  !--------------------------------------------------------------------------------
120 
121  !-------------ASSEMBLE TEMPERATURE MATRICES--------------------------------------------
122  ALLOCATE(temp_mat(m_max_c),temp_ksp(m_max_c))
123 
124  DO i = 1, m_max_c
125  mode = list_mode(i)
126 
127  !---TEMPERATURE MATRIX
128  CALL create_local_petsc_matrix(comm_one_d(1), temp_1_la, temp_mat(i), clean=.false.)
129  CALL qs_diff_mass_scal_m_variant(temp_mesh, temp_1_la, vol_heat_capacity, temp_diffusivity, &
130  1.5d0/dt, temp_list_robin_sides, convection_coeff, zero, mode, temp_mat(i)) ! MODIFICATION: double precision 1.5d0/dt instead of simple precision 1.5/dt & Robin sides and coefficients added
131 
132  IF (temp_per%n_bord/=0) THEN
133  CALL periodic_matrix_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, temp_mat(i), temp_1_la)
134  END IF
135 
136  CALL dirichlet_m_parallel(temp_mat(i),temp_mode_global_js_d(i)%DIL) ! MODIFICATION: axis BC
137 
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)
140  ENDDO
141 
142  ENDIF
143  tps_tot = user_time()
144  tps_cumul = 0
145 
146  !===Compute rhs by FFT at Gauss points
147  tps = user_time()
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
150  !===Compute pyromagnetic term if fhd
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
153  END IF
154  tps = user_time() - tps; tps_cumul=tps_cumul+tps
155  !WRITE(*,*) ' Tps fft vitesse', tps
156  !------------CONSTRUCTION OF rr_gauss------------------
157  index = 0
158  DO m = 1, temp_mesh%me
159  j_loc = temp_mesh%jj(:,m)
160  DO l = 1, temp_mesh%gauss%l_G
161  index = index + 1
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))
164  END DO
165  END DO
166  !------------DEBUT BOUCLE SUR LES MODES----------------
167  DO i = 1, m_max_c
168  mode = list_mode(i)
169 
170  !===RHS temperature
171  ff = (2d0/dt)*tempn(:,1,i) - (1d0/(2*dt))*tempn_m1(:,1,i) ! MODIFICATION: double precision (2 -> 2d0, 1 -> 1d0)
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)
174 
175  ff = (2d0/dt)*tempn(:,2,i) - (1d0/(2*dt))*tempn_m1(:,2,i) ! MODIFICATION: double precision (2 -> 2d0, 1 -> 1d0)
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)
178 
179  !===RHS Robins BCs
180  IF (mode == 0) THEN ! exterior temperature = constant
181  CALL qs_00_gauss_surface(temp_mesh, temp_1_la, temp_list_robin_sides, convection_coeff, exterior_temperature, cb_1) ! MODIFICATION: implementation of the term int_(partial Omega) h*Text*v, with h the convection coefficient
182  END IF
183 
184  !===RHS periodicity
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)
188  END IF
189  !----------------------------------------------------------
190 
191  n = SIZE(temp_js_d) ! MODIFICATION: axis BC
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)
197 
198  tps = user_time() - tps; tps_cumul=tps_cumul+tps
199  !WRITE(*,*) ' Tps second membre vitesse', tps
200  !-------------------------------------------------------------------------------------
201 
202  !--------------------INVERSION DES OPERATEURS--------------
203  tps = user_time()
204  !Solve system temp_c
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))
209 
210  !Solve system temp_s
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
216  !WRITE(*,*) ' Tps solution des pb de vitesse', tps, 'for mode ', mode
217  !-------------------------------------------------------------------------------------
218 
219  !---------------UPDATES-----------------------
220  tps = user_time()
221 
222  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
223  IF (mode==0) THEN
224  tempn_p1(:,2) = 0.d0
225  END IF
226  !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
227 
228 
229  tempn_m1(:,:,i) = tempn(:,:,i)
230  tempn(:,:,i) = tempn_p1
231 
232  ! Reconstruct dtempndt on Gauss points
233  ! WARNING FL (1/2/13) If reactivated, declare l and index
234  !tempn_p1 = (3*tempn_p1 - 4*tempn(:,:,i) + tempn_m1(:,:,i))/(2*dt)
235  !index = 0
236  !DO m = 1, temp_mesh%me
237  ! DO l = 1, temp_mesh%gauss%l_G
238  ! index = index + 1
239  ! dtempndt(index,1,i) = SUM(tempn_p1(temp_mesh%jj(:,m),1)*temp_mesh%gauss%ww(:,l)) + ff_conv(index,1,i)
240  ! dtempndt(index,2,i) = SUM(tempn_p1(temp_mesh%jj(:,m),2)*temp_mesh%gauss%ww(:,l)) + ff_conv(index,2,i)
241  ! END DO
242  !END DO
243  ! WARNING FL (1/2/13) If reactivated, declare l and index
244 
245  tps = user_time() - tps; tps_cumul=tps_cumul+tps
246  !WRITE(*,*) ' Tps des updates', tps
247  !-------------------------------------------------------------------------------------
248  ENDDO
249 
250  tps_tot = user_time() - tps_tot
251  !WRITE(*,'(A,2(f13.3,2x))') ' Tps boucle en temps Navier_stokes', tps_tot, tps_cumul
252  !WRITE(*,*) ' TIME = ', time, '========================================'
253 
254  END SUBROUTINE three_level_temperature
255  !============================================
256 
257  SUBROUTINE smb_ugradc_gauss_fft_par(communicator,mesh,list_mode,heat_capa_in,V_in,c_in,c_out)
258  !=================================
259  USE gauss_points
260  USE sft_parallele
261  USE chaine_caractere
262  USE boundary
263 #include "petsc/finclude/petsc.h"
264  USE petsc
265  IMPLICIT NONE
266  TYPE(mesh_type), TARGET :: mesh
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
284 
285  CALL gauss(mesh)
286  jj => mesh%jj
287  me => mesh%me
288 
289  tps = user_time()
290  DO i = 1, SIZE(list_mode)
291  mode = list_mode(i)
292  index = 0
293  DO m = 1, mesh%dom_me
294  j_loc = jj(:,m)
295  DO k = 1, 6
296  vs(:,k) = v_in(j_loc,k,i)
297  END DO
298  DO k = 1, 2
299  cs(:,k) = c_in(j_loc,k,i)
300  END DO
301  DO l = 1, l_g
302  index = index + 1
303  dw_loc = dw(:,:,l,m)
304 
305  !===Compute radius of Gauss point
306  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
307 
308  !-----------------vitesse sur les points de Gauss---------------------------
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))
312 
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))
316 
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,:))
321 
322  !-----------------gradient de c sur les points de Gauss---------------------------
323  !coeff sur les cosinus et sinus
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,:))
330 
331  gradc(index,:,i) = heat_capa_in(m) * gradc(index,:,i)
332 
333  cgauss(index,1,i) = sum(cs(:,1)*ww(:,l))
334  cgauss(index,2,i) = sum(cs(:,2)*ww(:,l))
335 
336  cgauss(index,:,i) = heat_capa_in(m) * cgauss(index,:,i)
337 
338  ENDDO
339  ENDDO
340  END DO
341 
342  !tps = user_time() - tps
343  !WRITE(*,*) ' Tps dans la grande boucle', tps
344  !tps = user_time()
345  temps = 0
346 
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
350  CALL fft_par_dot_prod_dcl(communicator, gradc, w, c_out, nb_procs, bloc_size, m_max_pad, temps)
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)
353  c_out = c_out + cint
354  tps = user_time() - tps
355  !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
356  !write(*,*) ' Temps de Comm ', temps(1)
357  !write(*,*) ' Temps de Calc ', temps(2)
358  !write(*,*) ' Temps de Chan ', temps(3)
359 
360  END SUBROUTINE smb_ugradc_gauss_fft_par
361 
362  SUBROUTINE mass_tot(communicator,mesh,tempn,RESLT)
363  !===========================
364  !moyenne
365  USE def_type_mesh
366 #include "petsc/finclude/petsc.h"
367  USE petsc
368  IMPLICIT NONE
369  TYPE(mesh_type) :: mesh
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
375  REAL(KIND=8) :: ray
376  mpi_comm :: communicator
377  r_loc = 0.d0
378 
379  DO m = 1, mesh%dom_me
380  j_loc = mesh%jj(:,m)
381  DO l = 1, mesh%gauss%l_G
382  ray = 0
383  DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
384  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
385  END DO
386 
387  r_loc = r_loc + sum(tempn(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
388 
389  ENDDO
390  ENDDO
391 
392  CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
393  reslt = r_out
394 
395  END SUBROUTINE mass_tot
396 
397  SUBROUTINE smb_pyromag_gauss_fft_par(communicator,mesh,list_mode,scal_in,vect_1_in,vect_2_in,vect_3_in,scal_out)
399  USE sft_parallele
400  USE def_type_mesh
401  USE input_data
402 #include "petsc/finclude/petsc.h"
403  USE petsc
404  IMPLICIT NONE
405  TYPE(mesh_type), INTENT(IN) :: mesh
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
415  REAL(KIND=8) :: rad
416  !===FOR FFT
417  INTEGER :: nb_procs, bloc_size, m_max_pad, code
418  mpi_comm :: communicator
419 
420  ! We compute the pyromagnetic coefficient term: T_dchi_dT_coeff(T) * 1/2 * D(H**2)/Dt
421 
422  !===nb_procs and m_max_pad calculus for FFT
423  CALL mpi_comm_size(communicator, nb_procs, code)
424  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
425 
426  !===T_dchi_dT_coeff(T) on Gauss nodes computation
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)
430  CALL gauss(mesh)
431  jj => mesh%jj
432  DO i = 1, SIZE(list_mode)
433  index = 0 ! global index of Gauss node
434  DO m = 1, mesh%dom_me
435  IF (minval(abs(mesh%i_d(m) - inputs%list_dom_ns)) == 0) THEN ! dchi/dT non null inside the (ferro)fluid domain only
436  DO l=1, l_g
437  index = index + 1
438  DO k = 1, 2
439  t_dchi_dt_coeff_gauss(index,k,i) = sum(t_dchi_dt_coeff(jj(:,m),k,i)*ww(:,l))
440  END DO
441  END DO
442  ELSE
443  t_dchi_dt_coeff_gauss(index+1:index+l_g,:,i) = 0.d0
444  index = index + l_g
445  END IF
446  END DO
447  END DO
448 
449  !===pdt(H**2) on Gauss nodes computation
450  IF (inputs%if_steady_current_fhd) THEN
451  pdth2 = 0.d0
452  ELSE
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)
456  index = 0 ! global index of Gauss node
457  DO m = 1, mesh%dom_me
458  IF (minval(abs(mesh%i_d(m) - inputs%list_dom_ns)) == 0) THEN ! dchi/dT non null inside the (ferro)fluid domain only
459  DO l=1, l_g
460  index = index + 1
461  DO k = 1, 2
462  pdth2(index,k,i) = 2*sum(v2_dot_v3(jj(:,m),k,i)*ww(:,l))
463  END DO
464  END DO
465  ELSE
466  pdth2(index+1:index+l_g,:,i) = 0.d0
467  index = index + l_g
468  END IF
469  END DO
470  END DO
471  END IF
472 
473  !===u on Gauss nodes computation
474  DO i = 1, SIZE(list_mode)
475  index = 0 ! global index of Gauss node
476  DO m = 1, mesh%dom_me
477  IF (minval(abs(mesh%i_d(m) - inputs%list_dom_ns)) == 0) THEN ! dchi/dT non null inside the (ferro)fluid domain only
478  DO l=1, l_g
479  index = index + 1
480  DO k = 1, 6
481  vect_1_in_gauss(index,k,i) = sum(vect_1_in(jj(:,m),k,i)*ww(:,l))
482  END DO
483  END DO
484  ELSE
485  vect_1_in_gauss(index+1:index+l_g,:,i) = 0.d0
486  index = index + l_g
487  END IF
488  END DO
489  END DO
490 
491  !===grad(H**2) on Gauss nodes computation
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)
495  mode = list_mode(i)
496  index = 0
497  DO m = 1, mesh%dom_me
498  IF (minval(abs(mesh%i_d(m) - inputs%list_dom_ns)) == 0) THEN ! dchi/dT non null inside the (ferro)fluid domain only
499  DO l=1, l_g
500  index = index + 1
501  dw_loc = dw(:,:,l,m)
502  rad = sum(mesh%rr(1,jj(:,m))*ww(:,l)) ! radius of Gauss node
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,:))
509  END DO
510  ELSE
511  grad_vect_2_in_square(index+1:index+l_g,:,i) = 0.d0
512  index = index + l_g
513  END IF
514  END DO
515  END DO
516 
517  !===u.grad(H**2) on Gauss nodes computation
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)
520 
521  !===D(H**2)/Dt on Gauss nodes computation
522  dh2_dt = pdth2 + ugradh2
523 
524  !===T_dchi_dT_coeff(T) * 1/2 * D(H**2)/Dt
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)
527 
528  END SUBROUTINE smb_pyromag_gauss_fft_par
529 
530 END MODULE subroutine_temperature
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:99
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:34
integer, public l_g
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:15
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)
Definition: fem_M_axi.f90:202
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)
Definition: solver.f90:147
type(my_data), public inputs
real(kind=8) function user_time()
Definition: my_util.f90:8
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)
Definition: solver.f90:12
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 gauss(mesh)
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)