SFEMaNS  version 5.3
Reference documentation for SFEMaNS
computing_rhs_gauss.f90
Go to the documentation of this file.
1 MODULE rhs_gauss_computing
2  PUBLIC:: rhs_ns_gauss_3x3, rhs_residual_ns_gauss_3x3, rhs_residual_ns_gauss_3x3_mom,&
3  rhs_ns_gauss_3x3_art_comp_mom
4  PRIVATE
5 CONTAINS
6 
7  SUBROUTINE rhs_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, pn_inc, rotv_v, &
8  rhs_gauss, opt_tempn)
9  !=================================
10  !RHS for Navier-Stokes
11  USE def_type_mesh
12  USE my_util
13  USE input_data
14  USE fem_rhs_axi
15  USE sft_parallele
16  USE boundary
17  IMPLICIT NONE
18  TYPE(mesh_type) :: vv_mesh, pp_mesh
19  REAL(KIND=8), INTENT(IN) :: time
20  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotv_v
21  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: V1m
22  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn, pn_inc
23  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
24  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: opt_tempn
25  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: rhs_gauss
26  REAL(KIND=8), DIMENSION(6) :: fs, ft, fp_inc
27  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
28  REAL(KIND=8), DIMENSION(vv_mesh%np,6) :: ff, imposed_vel
29  REAL(KIND=8), DIMENSION(vv_mesh%np,2) :: P, P_inc
30  REAL(KIND=8), DIMENSION(vv_mesh%gauss%k_d,vv_mesh%gauss%n_w) :: dw_loc
31  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6) :: imposed_vel_gauss
32  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6,SIZE(list_mode)) :: fp, rhs_gauss_penal
33  REAL(KIND=8), DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
34  REAL(KIND=8) :: ray
35  INTEGER :: m, l , i, k, index, TYPE
36  INTEGER :: nb_procs, m_max_pad, bloc_size
37 #include "petsc/finclude/petsc.h"
38  petscerrorcode :: ierr
39  mpi_comm :: communicator
40 
41  DO i = 1, SIZE(list_mode)
42  DO k= 1, 6
43  !IF (PRESENT(opt_tempn)) THEN
44  !TEST LC-CN 15/12/2016
45  IF (inputs%if_temperature) THEN
46  !TEST LC-CN 15/12/2016
47  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
48  opt_tempn=opt_tempn)
49  ELSE
50  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns')
51  END IF
52  END DO
53  DO k = 1, 2
54  !===JLG+HF Dec 10 2019
55  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn(:,k,i), p(:,k))
56  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn_inc(:,k,i), p_inc(:,k))
57  !CALL inject(pp_mesh%jj, vv_mesh%jj, pn(:,k,i), P(:,k))
58  !CALL inject(pp_mesh%jj, vv_mesh%jj, pn_inc(:,k,i), P_inc(:,k))
59  ENDDO
60 
61  index = 0
62  DO m = 1, vv_mesh%dom_me
63  j_loc = vv_mesh%jj(:,m)
64 
65  DO l = 1, vv_mesh%gauss%l_G
66  index = index +1
67  dw_loc = vv_mesh%gauss%dw(:,:,l,m)
68 
69  !===Compute radius of Gauss point
70  ray = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
71  rr_gauss(1,index) = ray
72  rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
73 
74  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
75  fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l))
76  ft(1) = sum(v1m(j_loc,1,i) * vv_mesh%gauss%ww(:,l))
77  fp(index,1,i) = -sum(p(j_loc,1)*dw_loc(1,:))
78  fp_inc(1) = -sum(p_inc(j_loc,1)*dw_loc(1,:))
79  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
80  fs(2) = sum(ff(j_loc,2) * vv_mesh%gauss%ww(:,l))
81  ft(2) = sum(v1m(j_loc,2,i) * vv_mesh%gauss%ww(:,l))
82  fp(index,2,i) = -sum(p(j_loc,2)*dw_loc(1,:))
83  fp_inc(2) = -sum(p_inc(j_loc,2)*dw_loc(1,:))
84  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
85  fs(3) = sum(ff(j_loc,3) * vv_mesh%gauss%ww(:,l))
86  ft(3) = sum(v1m(j_loc,3,i) * vv_mesh%gauss%ww(:,l))
87  fp(index,3,i) = -sum(p(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
88  fp_inc(3) = -sum(p_inc(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
89  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
90  fs(4) = sum(ff(j_loc,4) * vv_mesh%gauss%ww(:,l))
91  ft(4) = sum(v1m(j_loc,4,i) * vv_mesh%gauss%ww(:,l))
92  fp(index,4,i) = sum(p(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
93  fp_inc(4) = sum(p_inc(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
94  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
95  fs(5) = sum(ff(j_loc,5) * vv_mesh%gauss%ww(:,l))
96  ft(5) = sum(v1m(j_loc,5,i) * vv_mesh%gauss%ww(:,l))
97  fp(index,5,i) = -sum(p(j_loc,1)*dw_loc(2,:))
98  fp_inc(5) = -sum(p_inc(j_loc,1)*dw_loc(2,:))
99  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
100  fs(6) = sum(ff(j_loc,6) * vv_mesh%gauss%ww(:,l))
101  ft(6) = sum(v1m(j_loc,6,i) * vv_mesh%gauss%ww(:,l))
102  fp(index,6,i) = -sum(p(j_loc,2)*dw_loc(2,:))
103  fp_inc(6) = -sum(p_inc(j_loc,2)*dw_loc(2,:))
104 
105  rhs_gauss(index,:,i) = (ft+fp_inc+fs-rotv_v(index,:,i))
106 
107  ENDDO
108  ENDDO
109  IF (inputs%if_ns_penalty) THEN
110  IF(inputs%if_impose_vel_in_solids) THEN
111  IF (list_mode(i)==0) THEN
112  !Velocity imposed by penalty in solids (using BDF2).
113  imposed_vel(:,:) = 3.d0*imposed_velocity_by_penalty(vv_mesh%rr,time)/(2*inputs%dt)
114  index = 0
115  DO m = 1, vv_mesh%dom_me
116  j_loc = vv_mesh%jj(:,m)
117  DO l = 1, vv_mesh%gauss%l_G
118  index = index +1
119  DO TYPE = 1, 6
120  imposed_vel_gauss(index,type) = sum(imposed_vel(j_loc,type) &
121  * vv_mesh%gauss%ww(:,l))
122  END DO
123  END DO
124  END DO
125  rhs_gauss(:,:,i) = rhs_gauss(:,:,i) - imposed_vel_gauss(:,:)
126  ENDIF
127  END IF
128  END IF
129  END DO
130 
131  IF (inputs%if_ns_penalty) THEN
132  CALL mpi_comm_size(communicator, nb_procs, ierr)
133  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
134  bloc_size = SIZE(rhs_gauss,1)/nb_procs+1
135  CALL fft_par_var_eta_prod_gauss_dcl(communicator, penal_in_real_space, vv_mesh, &
136  rhs_gauss, rhs_gauss_penal, nb_procs, bloc_size, m_max_pad, rr_gauss, time)
137 
138  rhs_gauss = rhs_gauss_penal
139 
140  IF(inputs%if_impose_vel_in_solids) THEN
141  DO i = 1, SIZE(list_mode)
142  IF (list_mode(i)==0) THEN
143  rhs_gauss(:,:,i) = rhs_gauss_penal(:,:,i) + imposed_vel_gauss(:,:)
144  END IF
145  END DO
146  END IF
147  END IF
148 
149  rhs_gauss = rhs_gauss + fp
150 
151  END SUBROUTINE rhs_ns_gauss_3x3
152 
153  SUBROUTINE rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, du_dt,&
154  pn, rotv_v, rhs_gauss, opt_tempn)
155  !=================================
156  !RHS for Residual of Navier-Stokes
157  USE def_type_mesh
158  USE my_util
159  USE input_data
160  USE fem_rhs_axi
161  USE sft_parallele
162  USE boundary
163  IMPLICIT NONE
164  TYPE(mesh_type) :: vv_mesh, pp_mesh
165  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
166  REAL(KIND=8), INTENT(IN) :: time
167  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: du_dt
168  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn
169  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotv_v
170  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: opt_tempn
171  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: rhs_gauss
172  REAL(KIND=8), DIMENSION(6) :: fs
173  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
174  REAL(KIND=8), DIMENSION(vv_mesh%np,6) :: ff
175  REAL(KIND=8), DIMENSION(vv_mesh%np,2) :: P
176  REAL(KIND=8), DIMENSION(vv_mesh%gauss%k_d,vv_mesh%gauss%n_w) :: dw_loc
177  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6,SIZE(list_mode)) :: rhs_gauss_penal, fp, ft
178  REAL(KIND=8), DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
179  REAL(KIND=8) :: ray
180  INTEGER :: m, l , i, k, index
181  INTEGER :: nb_procs, m_max_pad, bloc_size
182 #include "petsc/finclude/petsc.h"
183  petscerrorcode :: ierr
184  mpi_comm :: communicator
185 
186  DO i = 1, SIZE(list_mode)
187  DO k = 1, 6
188  IF (PRESENT(opt_tempn)) THEN
189  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
190  opt_tempn=opt_tempn)
191  ELSE
192  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns')
193  END IF
194  END DO
195  DO k = 1, 2
196  !===JLG+HF Dec 10 2019
197  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn(:,k,i), p(:,k))
198  !CALL inject(pp_mesh%jj, vv_mesh%jj, pn(:,k,i), P(:,k))
199  ENDDO
200 
201  index = 0
202  DO m = 1, vv_mesh%dom_me
203  j_loc = vv_mesh%jj(:,m)
204 
205  DO l = 1, vv_mesh%gauss%l_G
206  index = index +1
207  dw_loc = vv_mesh%gauss%dw(:,:,l,m)
208 
209  !===Compute radius of Gauss point
210  ray = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
211  rr_gauss(1,index) = ray
212  rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
213 
214  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
215  fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l))
216  ft(index,1,i) = sum(du_dt(j_loc,1,i) * vv_mesh%gauss%ww(:,l))
217  fp(index,1,i) = sum(p(j_loc,1)*dw_loc(1,:))
218  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
219  fs(2) = sum(ff(j_loc,2) * vv_mesh%gauss%ww(:,l))
220  ft(index,2,i) = sum(du_dt(j_loc,2,i) * vv_mesh%gauss%ww(:,l))
221  fp(index,2,i) = sum(p(j_loc,2)*dw_loc(1,:))
222  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
223  fs(3) = sum(ff(j_loc,3) * vv_mesh%gauss%ww(:,l))
224  ft(index,3,i) = sum(du_dt(j_loc,3,i) * vv_mesh%gauss%ww(:,l))
225  fp(index,3,i) = sum(p(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
226  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
227  fs(4) = sum(ff(j_loc,4) * vv_mesh%gauss%ww(:,l))
228  ft(index,4,i) = sum(du_dt(j_loc,4,i) * vv_mesh%gauss%ww(:,l))
229  fp(index,4,i) = -sum(p(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
230  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
231  fs(5) = sum(ff(j_loc,5) * vv_mesh%gauss%ww(:,l))
232  ft(index,5,i) = sum(du_dt(j_loc,5,i) * vv_mesh%gauss%ww(:,l))
233  fp(index,5,i) = sum(p(j_loc,1)*dw_loc(2,:))
234  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
235  fs(6) = sum(ff(j_loc,6) * vv_mesh%gauss%ww(:,l))
236  ft(index,6,i) = sum(du_dt(j_loc,6,i) * vv_mesh%gauss%ww(:,l))
237  fp(index,6,i) = sum(p(j_loc,2)*dw_loc(2,:))
238 
239  rhs_gauss(index,:,i) = -fs+rotv_v(index,:,i)
240  ENDDO
241  ENDDO
242  END DO
243 
244  IF (inputs%if_ns_penalty) THEN
245  CALL mpi_comm_size(communicator, nb_procs, ierr)
246  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
247  bloc_size = SIZE(rhs_gauss,1)/nb_procs+1
248  CALL fft_par_var_eta_prod_gauss_dcl(communicator, penal_in_real_space, vv_mesh, &
249  rhs_gauss, rhs_gauss_penal, nb_procs, bloc_size, m_max_pad, rr_gauss, time)
250  rhs_gauss = rhs_gauss_penal
251  END IF
252 
253  !===Add term time derivative and grad pressure
254  rhs_gauss = rhs_gauss + ft + fp
255 
256  END SUBROUTINE rhs_residual_ns_gauss_3x3
257 
258  SUBROUTINE rhs_residual_ns_gauss_3x3_mom(vv_mesh, pp_mesh, list_mode, time, du_dt, pn, &
259  density, rotb_b, rhs_gauss, opt_tempn)
260  !=================================
261  !RHS for Navier-Stokes
262  USE def_type_mesh
263  USE my_util
264  USE input_data
265  USE fem_rhs_axi
266  USE sft_parallele
267  USE boundary
268  IMPLICIT NONE
269  TYPE(mesh_type) :: vv_mesh, pp_mesh
270  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
271  REAL(KIND=8), INTENT(IN) :: time
272  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: du_dt
273  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn
274  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: density
275  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotb_b
276  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: opt_tempn
277  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: rhs_gauss
278  REAL(KIND=8), DIMENSION(6) :: fs, ft, fp
279  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
280  REAL(KIND=8), DIMENSION(vv_mesh%np,6) :: ff
281  REAL(KIND=8), DIMENSION(vv_mesh%np,2) :: P
282  REAL(KIND=8), DIMENSION(vv_mesh%gauss%k_d,vv_mesh%gauss%n_w) :: dw_loc
283  REAL(KIND=8), DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
284  REAL(KIND=8) :: ray
285  INTEGER :: m, l , i, k, index
286 
287  DO i = 1, SIZE(list_mode)
288  DO k = 1, 6
289  IF (inputs%if_temperature) THEN
290  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
291  opt_density=density, opt_tempn=opt_tempn)
292  ELSE
293  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
294  opt_density=density)
295  END IF
296  END DO
297  DO k = 1, 2
298  !===JLG+HF Dec 10 2019
299  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn(:,k,i), p(:,k))
300  !CALL inject(pp_mesh%jj, vv_mesh%jj, pn(:,k,i), P(:,k))
301  ENDDO
302 
303  index = 0
304  DO m = 1, vv_mesh%dom_me
305  j_loc = vv_mesh%jj(:,m)
306 
307  DO l = 1, vv_mesh%gauss%l_G
308  index = index +1
309  dw_loc = vv_mesh%gauss%dw(:,:,l,m)
310 
311  !===Compute radius of Gauss point
312  ray = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
313  rr_gauss(1,index) = ray
314  rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
315 
316  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
317  fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l))
318  ft(1) = sum(du_dt(j_loc,1,i) * vv_mesh%gauss%ww(:,l))
319  fp(1) = sum(p(j_loc,1)*dw_loc(1,:))
320  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
321  fs(2) = sum(ff(j_loc,2) * vv_mesh%gauss%ww(:,l))
322  ft(2) = sum(du_dt(j_loc,2,i) * vv_mesh%gauss%ww(:,l))
323  fp(2) = sum(p(j_loc,2)*dw_loc(1,:))
324  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
325  fs(3) = sum(ff(j_loc,3) * vv_mesh%gauss%ww(:,l))
326  ft(3) = sum(du_dt(j_loc,3,i) * vv_mesh%gauss%ww(:,l))
327  fp(3) = sum(p(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
328  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
329  fs(4) = sum(ff(j_loc,4) * vv_mesh%gauss%ww(:,l))
330  ft(4) = sum(du_dt(j_loc,4,i) * vv_mesh%gauss%ww(:,l))
331  fp(4) = -sum(p(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
332  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
333  fs(5) = sum(ff(j_loc,5) * vv_mesh%gauss%ww(:,l))
334  ft(5) = sum(du_dt(j_loc,5,i) * vv_mesh%gauss%ww(:,l))
335  fp(5) = sum(p(j_loc,1)*dw_loc(2,:))
336  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
337  fs(6) = sum(ff(j_loc,6) * vv_mesh%gauss%ww(:,l))
338  ft(6) = sum(du_dt(j_loc,6,i) * vv_mesh%gauss%ww(:,l))
339  fp(6) = sum(p(j_loc,2)*dw_loc(2,:))
340 
341  rhs_gauss(index,:,i) = ft+fp-fs-rotb_b(index,:,i)
342  ENDDO
343  ENDDO
344  END DO
345 
346  END SUBROUTINE rhs_residual_ns_gauss_3x3_mom
347 
348  SUBROUTINE rhs_ns_gauss_3x3_art_comp_mom(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, rotv_v, &
349  rhs_gauss, tempn, density)
350  !=================================
351  !RHS for Navier-Stokes
352  USE def_type_mesh
353  USE my_util
354  USE input_data
355  USE fem_rhs_axi
356  USE sft_parallele
357  USE boundary
358  IMPLICIT NONE
359  TYPE(mesh_type) :: vv_mesh, pp_mesh
360  REAL(KIND=8), INTENT(IN) :: time
361  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: rotv_v
362  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: V1m
363  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: pn
364  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
365  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tempn
366  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: density
367  REAL(KIND=8), DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)), INTENT(OUT) :: rhs_gauss
368  REAL(KIND=8), DIMENSION(6) :: fs, ft
369  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
370  REAL(KIND=8), DIMENSION(vv_mesh%np,6) :: ff
371  REAL(KIND=8), DIMENSION(vv_mesh%np,2) :: P
372  REAL(KIND=8), DIMENSION(vv_mesh%gauss%k_d,vv_mesh%gauss%n_w) :: dw_loc
373  REAL(KIND=8), DIMENSION(vv_mesh%dom_me*vv_mesh%gauss%l_G,6,SIZE(list_mode)) :: fp
374  REAL(KIND=8), DIMENSION(2,vv_mesh%gauss%l_G*vv_mesh%dom_me) :: rr_gauss
375  REAL(KIND=8) :: ray
376  INTEGER :: m, l , i, k, index
377 #include "petsc/finclude/petsc.h"
378  petscerrorcode :: ierr
379  petscmpiint :: rank
380  mpi_comm :: communicator
381  CALL mpi_comm_rank(communicator,rank,ierr)
382  DO i = 1, SIZE(list_mode)
383  DO k = 1, 6
384  IF (inputs%if_temperature) THEN
385  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
386  opt_density=density, opt_tempn=tempn)
387  ELSE
388  ff(:,k) = source_in_ns_momentum(k, vv_mesh%rr, list_mode(i), i, time, inputs%Re, 'ns', &
389  opt_density=density)
390  END IF
391  END DO
392  DO k = 1, 2
393  CALL inject_generic(inputs%type_fe_velocity, pp_mesh%jj, vv_mesh%jj, pn(:,k,i), p(:,k))
394  !CALL inject(pp_mesh%jj, vv_mesh%jj, pn(:,k,i), P(:,k))
395  ENDDO
396 
397  index = 0
398  DO m = 1, vv_mesh%dom_me
399  j_loc = vv_mesh%jj(:,m)
400 
401  DO l = 1, vv_mesh%gauss%l_G
402  index = index +1
403  dw_loc = vv_mesh%gauss%dw(:,:,l,m)
404 
405  !===Compute radius of Gauss point
406  ray = sum(vv_mesh%rr(1,j_loc)*vv_mesh%gauss%ww(:,l))
407  rr_gauss(1,index) = ray
408  rr_gauss(2,index) = sum(vv_mesh%rr(2,j_loc)*vv_mesh%gauss%ww(:,l))
409 
410  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
411  fs(1) = sum(ff(j_loc,1) * vv_mesh%gauss%ww(:,l))
412  ft(1) = sum(v1m(j_loc,1,i) * vv_mesh%gauss%ww(:,l))
413  fp(index,1,i) = -sum(p(j_loc,1)*dw_loc(1,:))
414  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
415  fs(2) = sum(ff(j_loc,2) * vv_mesh%gauss%ww(:,l))
416  ft(2) = sum(v1m(j_loc,2,i) * vv_mesh%gauss%ww(:,l))
417  fp(index,2,i) = -sum(p(j_loc,2)*dw_loc(1,:))
418  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
419  fs(3) = sum(ff(j_loc,3) * vv_mesh%gauss%ww(:,l))
420  ft(3) = sum(v1m(j_loc,3,i) * vv_mesh%gauss%ww(:,l))
421  fp(index,3,i) = -sum(p(j_loc,2)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
422  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
423  fs(4) = sum(ff(j_loc,4) * vv_mesh%gauss%ww(:,l))
424  ft(4) = sum(v1m(j_loc,4,i) * vv_mesh%gauss%ww(:,l))
425  fp(index,4,i) = sum(p(j_loc,1)*vv_mesh%gauss%ww(:,l))/ray*list_mode(i)
426  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
427  fs(5) = sum(ff(j_loc,5) * vv_mesh%gauss%ww(:,l))
428  ft(5) = sum(v1m(j_loc,5,i) * vv_mesh%gauss%ww(:,l))
429  fp(index,5,i) = -sum(p(j_loc,1)*dw_loc(2,:))
430  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
431  fs(6) = sum(ff(j_loc,6) * vv_mesh%gauss%ww(:,l))
432  ft(6) = sum(v1m(j_loc,6,i) * vv_mesh%gauss%ww(:,l))
433  fp(index,6,i) = -sum(p(j_loc,2)*dw_loc(2,:))
434 
435  rhs_gauss(index,:,i) = (ft+fs-rotv_v(index,:,i))
436 
437  ENDDO
438  ENDDO
439  END DO
440 
441  rhs_gauss = rhs_gauss + fp
442 
443  END SUBROUTINE rhs_ns_gauss_3x3_art_comp_mom
444 
445  !===Local subroutine
446 !!$ SUBROUTINE inject(jj_c, jj_f, pp_c, pp_f)
447 !!$ IMPLICIT NONE
448 !!$ INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_c, jj_f
449 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_c
450 !!$ REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_f
451 !!$ REAL(KIND=8) :: half = 0.5
452 !!$ INTEGER:: m
453 !!$ IF (SIZE(jj_c,1)==3) THEN
454 !!$ DO m = 1, SIZE(jj_f,2)
455 !!$ pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
456 !!$ pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
457 !!$ pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
458 !!$ pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
459 !!$ END DO
460 !!$
461 !!$ ELSE
462 !!$ DO m = 1, SIZE(jj_f,2)
463 !!$ pp_f(jj_f(1:4,m)) = pp_c(jj_c(:,m))
464 !!$ END DO
465 !!$ pp_f(jj_f(5,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(4,:)))*half
466 !!$ pp_f(jj_f(6,:)) = (pp_c(jj_c(4,:)) + pp_c(jj_c(2,:)))*half
467 !!$ pp_f(jj_f(7,:)) = (pp_c(jj_c(2,:)) + pp_c(jj_c(3,:)))*half
468 !!$ pp_f(jj_f(8,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(4,:)))*half
469 !!$ pp_f(jj_f(9,:)) = (pp_c(jj_c(3,:)) + pp_c(jj_c(1,:)))*half
470 !!$ pp_f(jj_f(10,:)) = (pp_c(jj_c(1,:)) + pp_c(jj_c(2,:)))*half
471 !!$
472 !!$ END IF
473 !!$
474 !!$ END SUBROUTINE inject
475 END MODULE rhs_gauss_computing
rhs_gauss_computing::rhs_ns_gauss_3x3_art_comp_mom
subroutine, public rhs_ns_gauss_3x3_art_comp_mom(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, rotv_v, rhs_gauss, tempn, density)
Definition: computing_rhs_gauss.f90:350
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
sft_parallele::fft_par_var_eta_prod_gauss_dcl
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)
Definition: fft_parallel.f90:1851
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
rhs_gauss_computing::rhs_ns_gauss_3x3
subroutine, public rhs_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, pn_inc, rotv_v, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:9
sft_parallele
Definition: fft_parallel.f90:4
rhs_gauss_computing
Definition: computing_rhs_gauss.f90:1
rhs_gauss_computing::rhs_residual_ns_gauss_3x3_mom
subroutine, public rhs_residual_ns_gauss_3x3_mom(vv_mesh, pp_mesh, list_mode, time, du_dt, pn, density, rotb_b, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:260
rhs_gauss_computing::rhs_residual_ns_gauss_3x3
subroutine, public rhs_residual_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, du_dt, pn, rotv_v, rhs_gauss, opt_tempn)
Definition: computing_rhs_gauss.f90:155
input_data
Definition: read_my_data.f90:233
fem_rhs_axi
Definition: fem_rhs_axi.f90:1
def_type_mesh
Definition: def_type_mesh.f90:4
fem_rhs_axi::inject_generic
subroutine inject_generic(type_fe_velocity, jj_c, jj_f, pp_c, pp_f)
Definition: fem_rhs_axi.f90:753