OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_fluxes_computation.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| multi_fluxes_computation ../engine/source/multifluid/multi_fluxes_computation.f
25!||--- called by ------------------------------------------------------
26!|| multi_timeevolution ../engine/source/multifluid/multi_timeevolution.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!||--- uses -----------------------------------------------------
30!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
31!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| initbuf_mod ../engine/share/resol/initbuf.F
34!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
35!||====================================================================
36 SUBROUTINE multi_fluxes_computation(NG, ELBUF_TAB, IPARG, ITASK, IXS, IXQ, IXTG,
37 . PM, IPM, MULTI_FVM, ALE_CONNECTIVITY, WGRID, XGRID, ITAB, NBMAT, CURRENT_TIME, BUFMAT,
38 . ID_GLOBAL_VOIS,NPF,TF)
39!$COMMENT
40! MULTI_FLUXES_COMPUTATION description :
41! computation of fluxes with 1rst order algorithm
42! MULTI_FLUXES_COMPUTATION organization :
43! The parith/on is ensured by the same
44! order computation :
45! * check the user ID
46! * if user ID( element1) < user ID( element2)
47! --> element1 drives the computation
48! * if user ID( element1) > user ID( element2)
49! --> element2 drives the computation
50!
51!$ENDCOMMENT
52C-----------------------------------------------
53C M o d u l e s
54C-----------------------------------------------
55 USE initbuf_mod
56 USE elbufdef_mod
57 USE multi_fvm_mod
58 USE ale_connectivity_mod
59 use element_mod , only : nixs,nixq,nixtg
60C-----------------------------------------------
61C I m p l i c i t T y p e s
62C-----------------------------------------------
63#include "implicit_f.inc"
64C-----------------------------------------------
65C C o m m o n B l o c k s
66C-----------------------------------------------
67#include "com01_c.inc"
68! NUMELS
69#include "com04_c.inc"
70#include "param_c.inc"
71#include "mvsiz_p.inc"
72! DTFAC1
73 COMMON /tablesizf/ stf,snpc
74 INTEGER STF,SNPC
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER, INTENT(IN) :: NG
79 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
80 INTEGER, INTENT(IN) :: IPARG(NPARG, *)
81 INTEGER, INTENT(IN) :: ITASK ! SMP TASK
82 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
83 INTEGER, INTENT(IN) :: IPM(NPROPMI, *)
84 my_real, INTENT(IN) :: pm(npropm, *)
85 TYPE(multi_fvm_struct), INTENT(INOUT) :: MULTI_FVM
86 ! for parith/on : ID_GLOBAL_VOIS --> user id
87 INTEGER, INTENT(IN) :: ID_GLOBAL_VOIS(*)
88 my_real, INTENT(IN) :: WGRID(3, *), XGRID(3, *), CURRENT_TIME
89 INTEGER, INTENT(IN) :: ITAB(*), NBMAT
90 my_real, INTENT(INOUT) :: bufmat(*)
91 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
92 INTEGER, INTENT(IN) :: NPF(SNPC)
93 my_real, INTENT(IN) :: tf(stf)
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 TYPE(g_bufel_), POINTER :: GBUF, GBUFJJ
98 TYPE(L_BUFEL_), POINTER :: LBUF
99
100 INTEGER :: II, JJ, KFACE, I, J, KFACE2, NGJJ, NFTJJ, NELJJ, IFACE
101 INTEGER :: NEL, ISOLNOD, ITY, NFT
102 INTEGER :: IMAT
103 my_real :: x1(3), x2(3), x3(3), x4(3)
104 my_real :: lambdaii, lambdaf, normuii, normujj
105 my_real :: fii(5), sub_fii(3), fjj(5), normal_velii, normal_veljj, vii(5), sub_vii(3), vjj(5),
106 . velii2, veljj2, surf
107 my_real :: sl, sr, sstar, pstar, fiistar(5), sub_fiistar(3), fjjstar(5), viistar(5), vjjstar(5), pp(5),
108 . sub_viistar(3)
109 my_real :: sspii, sspjj, rhoii, rhojj, pii, pjj, nx, ny, nz, normalw
110 INTEGER :: NODE1, NODE2, NODE3, NODE4
111 my_real :: vxii, vyii, vzii
112 my_real :: vxjj, vyjj, vzjj, vx, vy, vz, normvel, ssp, rho, p
113 INTEGER :: NB_FACE, ELEMID, NUMEL_SPMD
114 my_real :: machii, machjj, pstar2, theta, alphaii, sub_rhoii, sub_rhoeintii,
115 . alphastar, sub_rhostar, sub_estar, sub_pii, wfac(3)
116 TYPE(lbuf_ptr) :: LBUFS(NBMAT)
117 INTEGER :: LOCAL_MATID(NBMAT), MATLAW(NBMAT)
118 my_real :: xf(3), xk(3), xl(3)
119 my_real :: massflux1, massflux2, pshift
120
121 INTEGER :: KFACE3,I_OLD,J_OLD,KFACE_OLD,KFACE2_OLD
122 INTEGER :: ID_GLOBAL_VOIS_1,ID_GLOBAL_VOIS_2,IJK
123 INTEGER, DIMENSION(MVSIZ) :: GLOBAL_ID_CURRENT_ELM
124 my_real :: direction
125 INTEGER :: IAD
126
127 nb_face = -1
128
129 gbuf => elbuf_tab(ng)%GBUF
130 nel = iparg(2, ng)
131 nft = iparg(3, ng)
132 ity = iparg(5, ng)
133 isolnod = iparg(28, ng)
134
135 pshift = multi_fvm%PRES_SHIFT
136
137 IF (nbmat > 1) THEN
138!DIR$ NOVECTOR
139 DO imat = 1, nbmat
140 lbufs(imat)%LBUF => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
141 ENDDO
142 ENDIF
143C Normal and surface computation
144 SELECT CASE (multi_fvm%SYM)
145 CASE (0)
146 nb_face = 6
147C 3D
148 numel_spmd = numels
149 DO imat = 1, nbmat
150 local_matid(imat) = ipm(20 + imat, ixs(1, 1 + nft))
151 ENDDO
152 DO ii = 1, nel
153 i = ii + nft
154 global_id_current_elm(ii) = ixs(nixs,i)
155 ENDDO
156 CASE (1, 2)
157 IF (ity == 2) THEN
158C QUADS
159 nb_face = 4
160 numel_spmd = numelq
161 DO imat = 1, nbmat
162 local_matid(imat) = ipm(20 + imat, ixq(1, 1 + nft))
163 ENDDO
164 DO ii = 1, nel
165 i = ii + nft
166 global_id_current_elm(ii) = ixq(nixq,i)
167 ENDDO
168 ELSEIF (ity == 7) THEN
169C TRIANGLES
170 nb_face = 3
171 numel_spmd = numeltg
172 DO imat = 1, nbmat
173 local_matid(imat) = ipm(20 + imat, ixtg(1, 1 + nft))
174 ENDDO
175 DO ii = 1, nel
176 i = ii + nft
177 global_id_current_elm(ii) = ixtg(nixtg,i)
178 ENDDO
179 ENDIF
180 CASE DEFAULT
181 CALL arret(2)
182 END SELECT
183
184C Computation flow
185 DO ii = 1, nel
186 i = ii + nft
187 iad = ale_connectivity%ee_connect%iad_connect(i)
188 nb_face =ale_connectivity%ee_connect%iad_connect(i+1) -
189 . ale_connectivity%ee_connect%iad_connect(i)
190 i_old = i
191C Face KFACE
192 DO kface3 = 1, nb_face
193 jj = ale_connectivity%ee_connect%connected(iad + kface3 - 1)
194 j_old = jj
195C Face normal
196 nx = multi_fvm%FACE_DATA%NORMAL(1, kface3, i_old)
197 ny = multi_fvm%FACE_DATA%NORMAL(2, kface3, i_old)
198 nz = multi_fvm%FACE_DATA%NORMAL(3, kface3, i_old)
199 surf = multi_fvm%FACE_DATA%SURF(kface3, i_old)
200 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface3, i_old)
201
202 id_global_vois_1 = global_id_current_elm(ii)
203 id_global_vois_2 = id_global_vois( nb_face * (i_old - 1) + kface3 )
204! -----------------------------------------------
205 IF (j_old > 0 .AND. i_old < j_old) THEN
206 ! ----------------------
207 kface_old = kface3
208 kface2_old = multi_fvm%FVM_CONNECTIVITY%KVOIS(nb_face * (i_old - 1) + kface3)
209
210 IF(id_global_vois_1<id_global_vois_2) THEN
211 direction = one
212 kface = kface3
213 kface2 = multi_fvm%FVM_CONNECTIVITY%KVOIS(nb_face * (i_old - 1) + kface)
214 i = i_old
215 jj = j_old
216 ELSE
217 kface2 = kface3
218 kface = multi_fvm%FVM_CONNECTIVITY%KVOIS(nb_face * (i_old - 1) + kface3)
219
220 ijk = i
221 i = j_old
222 jj = i_old
223 direction = -one
224 ENDIF
225C Face normal
226 nx = nx *direction
227 ny = ny *direction
228 nz = nz *direction
229
230 ! ----------------------
231 j = jj
232! Conserved variable current element
233 rhoii = multi_fvm%RHO(i)
234 vxii = multi_fvm%VEL(1, i)
235 vyii = multi_fvm%VEL(2, i)
236 vzii = multi_fvm%VEL(3, i)
237 velii2 = vxii * vxii + vyii * vyii + vzii * vzii
238
239! Pressure current element
240 pii = multi_fvm%PRES(i)
241! Sound speed current element
242 sspii = multi_fvm%SOUND_SPEED(i)
243
244
245C Densities
246 rhojj = multi_fvm%RHO(j)
247C Pressures
248 pjj = multi_fvm%PRES(j)
249C Normal velocity current element
250 normal_velii = vxii * nx + vyii * ny + vzii * nz
251
252! Normal velocity adjacent element
253 vxjj = multi_fvm%VEL(1, j)
254 vyjj = multi_fvm%VEL(2, j)
255 vzjj = multi_fvm%VEL(3, j)
256
257 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
258 veljj2 = vxjj * vxjj + vyjj * vyjj + vzjj * vzjj
259C Sound speeds adjacent element
260 sspjj = multi_fvm%SOUND_SPEED(j)
261C Normal grid velocity
262 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
263
264C Local Mach numbers
265 machii = abs(normal_velii) / sspii
266 machjj = abs(normal_veljj) / sspjj
267
268C HLL wave speed estimates
269 sl = min(normal_velii - sspii, normal_veljj - sspjj)
270 sr = max(normal_velii + sspii, normal_veljj + sspjj)
271
272C Intermediate wave speed
273 sstar = pjj - pii + rhoii * normal_velii * (sl - normal_velii) -
274 . rhojj * normal_veljj * (sr - normal_veljj)
275
276 sstar = sstar / (rhoii * (sl - normal_velii) -
277 . rhojj * (sr - normal_veljj))
278
279
280C Specific for Low Mach number corrections
281C Intermediate pressure
282 pstar2 = pii + rhoii * (sstar - normal_velii) * (sl - normal_velii)
283C PSTAR2 = HALF * (PSTAR + PJJ + RHOJJ * (SSTAR - NORMAL_VELJJ) * (SR - NORMAL_VELJJ))
284 IF (min(machii, machjj) < em01) THEN
285 theta = min(machii, machjj)
286 ELSE
287 theta = one
288 ENDIF
289 pstar = (one - theta) * half * (pii + pjj) + theta * pstar2
290 IF (multi_fvm%LOWMACH_OPT) THEN
291 pp(1) = zero
292 pp(2) = pstar * nx
293 pp(3) = pstar * ny
294 pp(4) = pstar * nz
295 pp(5) = sstar * (pstar + pshift)
296 ELSE
297 pp(1) = zero
298 pp(2) = pstar2 * nx
299 pp(3) = pstar2 * ny
300 pp(4) = pstar2 * nz
301 pp(5) = sstar * (pstar2 + pshift)
302 ENDIF
303
304 IF (sl > normalw) THEN
305C Take the fluxes of cell II
306C ===
307C Global fluxes
308 vii(1) = rhoii
309 vii(2) = rhoii * vxii
310 vii(3) = rhoii * vyii
311 vii(4) = rhoii * vzii
312 vii(5) = multi_fvm%EINT(i) + half * rhoii * velii2
313C Normal physical flux current element
314 fii(1) = vii(1) * normal_velii
315 fii(2) = vii(2) * normal_velii + pii * nx
316 fii(3) = vii(3) * normal_velii + pii * ny
317 fii(4) = vii(4) * normal_velii + pii * nz
318 fii(5) = (vii(5) + pii + pshift) * normal_velii
319 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction * (fii(1:5) - normalw * vii(1:5)) * surf
320 multi_fvm%FLUXES(6, kface_old, i_old) = direction * normal_velii * surf
321 massflux1 = multi_fvm%FLUXES(1, kface_old, i_old)
322C ===
323C Submaterial fluxes
324 IF (nbmat > 1) THEN
325 massflux2 = zero
326 DO imat = 1, nbmat
327 alphaii = multi_fvm%PHASE_ALPHA(imat, i)
328 sub_rhoii = multi_fvm%PHASE_RHO(imat, i)
329 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, i)
330 sub_viistar(1) = alphaii
331 sub_viistar(2) = alphaii * sub_rhoii
332 sub_viistar(3) = alphaii * sub_rhoeintii
333 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
334 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
335 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
336 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
337 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
338 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
339 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
340 massflux2 = massflux2 + multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old)
341 ENDDO
342 ENDIF
343C
344 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
345C Global fluxes
346 vii(1) = rhoii
347 vii(2) = rhoii * vxii
348 vii(3) = rhoii * vyii
349 vii(4) = rhoii * vzii
350 vii(5) = multi_fvm%EINT(i) + half * rhoii * velii2
351C Normal physical flux current element
352 fii(1) = vii(1) * normal_velii
353 fii(2) = vii(2) * normal_velii + pii * nx
354 fii(3) = vii(3) * normal_velii + pii * ny
355 fii(4) = vii(4) * normal_velii + pii * nz
356 fii(5) = (vii(5) + pii + pshift) * normal_velii
357C Take intermediate state flux (HLLC scheme)
358C ===
359C Global fluxes
360 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
361 viistar(1:5) = viistar(1:5) / (sstar - sl)
362 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
363 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction *
364 . (fiistar(1:5) - normalw * viistar(1:5)) * surf
365 multi_fvm%FLUXES(6, kface_old, i_old) = direction * sstar * surf
366 massflux1 = multi_fvm%FLUXES(1, kface_old, i_old)
367C ===
368C Submaterial fluxes
369 IF (nbmat > 1) THEN
370 massflux2 = zero
371 DO imat = 1, nbmat
372 matlaw(imat) = ipm(2, local_matid(imat))
373 alphaii = multi_fvm%PHASE_ALPHA(imat, i)
374 sub_rhoii = multi_fvm%PHASE_RHO(imat, i)
375 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, i)
376 alphastar = alphaii
377 sub_rhostar = sub_rhoii * (normal_velii - sl) / (sstar - sl)
378 sub_pii = multi_fvm%PHASE_PRES(imat, i)
379 IF (alphastar > zero) THEN
380c$$$ CALL ENERGY_JUMP(MATLAW(IMAT), LOCAL_MATID(IMAT), PM, IPM, NPROPM, NPROPMI,
381c$$$ . SUB_RHOEINTII / SUB_RHOII, SUB_RHOSTAR, SUB_PII, (SSTAR - NORMAL_VELII) / (SSTAR - SL),
382c$$$ . LBUFS(IMAT)%LBUF%BFRAC(II), GBUF%TB(II), GBUF%DELTAX(II), CURRENT_TIME,
383c$$$ . SUB_ESTAR, BUFMAT)
384 sub_estar = sub_rhoeintii / sub_rhoii -
385 . sub_pii * (one / sub_rhostar - one / sub_rhoii)
386 ELSE
387 sub_estar = zero
388 ENDIF
389 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, i)
390 sub_viistar(1) = alphastar
391 sub_viistar(2) = alphastar * sub_rhostar
392 sub_viistar(3) = alphaii * sub_rhostar * sub_estar
393 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
394 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
395 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
396 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
397 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
398 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
399 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
400 massflux2 = massflux2 + multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old)
401 ENDDO
402 ENDIF
403 ELSEIF (sstar < normalw .AND. normalw <= sr) THEN
404C Global fluxes
405 vii(1) = rhojj
406 vii(2) = rhojj * vxjj
407 vii(3) = rhojj * vyjj
408 vii(4) = rhojj * vzjj
409 vii(5) = multi_fvm%EINT(j) + half * rhojj * veljj2
410C Normal physical flux current element
411 fii(1) = vii(1) * normal_veljj
412 fii(2) = vii(2) * normal_veljj + pjj * nx
413 fii(3) = vii(3) * normal_veljj + pjj * ny
414 fii(4) = vii(4) * normal_veljj + pjj* nz
415 fii(5) = (vii(5) + pjj + pshift) * normal_veljj
416C Take intermediate state flux (HLLC scheme)
417C ===
418C Global fluxes
419 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
420 viistar(1:5) = viistar(1:5) / (sstar - sr)
421 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
422 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction *
423 . (fiistar(1:5) - normalw * viistar(1:5)) * surf
424 multi_fvm%FLUXES(6, kface_old, i_old) = direction * sstar * surf
425 massflux1 = multi_fvm%FLUXES(1, kface_old, i_old)
426C ===
427C Submaterial fluxes
428 IF (nbmat > 1) THEN
429 massflux2 = zero
430 DO imat = 1, nbmat
431 matlaw(imat) = ipm(2, local_matid(imat))
432 alphaii = multi_fvm%PHASE_ALPHA(imat, j)
433 sub_rhoii = multi_fvm%PHASE_RHO(imat, j)
434 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, j)
435 alphastar = alphaii
436 sub_rhostar = sub_rhoii * (normal_veljj - sr) / (sstar - sr)
437 sub_pii = multi_fvm%PHASE_PRES(imat, j)
438 IF (alphastar > zero) THEN
439c$$$ CALL ENERGY_JUMP(MATLAW(IMAT), LOCAL_MATID(IMAT), PM, IPM, NPROPM, NPROPMI,
440c$$$ . SUB_RHOEINTII / SUB_RHOII, SUB_RHOSTAR, SUB_PII, (SSTAR - NORMAL_VELJJ) / (SSTAR - SR),
441c$$$ . LBUFS(IMAT)%LBUF%BFRAC(II), GBUF%TB(II), GBUF%DELTAX(II), CURRENT_TIME,
442c$$$ . SUB_ESTAR, BUFMAT)
443 sub_estar = sub_rhoeintii / sub_rhoii -
444 . sub_pii * (one / sub_rhostar - one / sub_rhoii)
445
446 ELSE
447 sub_estar = zero
448 ENDIF
449 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, j)
450 sub_viistar(1) = alphastar
451 sub_viistar(2) = alphastar * sub_rhostar
452 sub_viistar(3) = alphaii * sub_rhostar * sub_estar
453 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
454 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
455 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
456 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
457 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
458 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
459 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
460 massflux2 = massflux2 + multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old)
461 ENDDO
462 ENDIF
463 ELSE IF (sr < normalw) THEN
464C Take the fluxes of cell JJ
465C ===
466C Global fluxes
467 vii(1) = rhojj
468 vii(2) = rhojj * vxjj
469 vii(3) = rhojj * vyjj
470 vii(4) = rhojj * vzjj
471 vii(5) = multi_fvm%EINT(j) + half * rhojj * veljj2
472C Normal physical flux current element
473 fii(1) = vii(1) * normal_veljj
474 fii(2) = vii(2) * normal_veljj + pjj * nx
475 fii(3) = vii(3) * normal_veljj + pjj * ny
476 fii(4) = vii(4) * normal_veljj + pjj * nz
477 fii(5) = (vii(5) + pjj + pshift) * normal_veljj
478 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction * (fii(1:5) - normalw * vii(1:5)) * surf
479 multi_fvm%FLUXES(6, kface_old, i_old) = direction * normal_veljj * surf
480 massflux1 = multi_fvm%FLUXES(1, kface_old, i_old)
481C ===
482C Submaterial fluxes
483 IF (nbmat > 1) THEN
484 massflux2 = zero
485 DO imat = 1, nbmat
486 alphaii = multi_fvm%PHASE_ALPHA(imat, j)
487 sub_rhoii = multi_fvm%PHASE_RHO(imat, j)
488 sub_rhoeintii = multi_fvm%PHASE_EINT(imat, j)
489 sub_viistar(1) = alphaii
490 sub_viistar(2) = alphaii * sub_rhoii
491 sub_viistar(3) = alphaii * sub_rhoeintii
492 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
493 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
494 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
495 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
496 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
497 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
498 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
499 massflux2 = massflux2 + multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old)
500 ENDDO
501 ENDIF
502 ELSE
503 CALL arret(2)
504C
505 ENDIF
506
507 IF (j_old <= numel_spmd) THEN
508 multi_fvm%FLUXES(1:6, kface2_old, j_old) = -multi_fvm%FLUXES(1:6, kface_old, i_old)
509 IF (nbmat > 1) THEN
510 DO imat = 1, nbmat
511 multi_fvm%SUBVOL_FLUXES(imat, kface2_old, j_old) = -multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old)
512 multi_fvm%SUBMASS_FLUXES(imat, kface2_old, j_old) = -multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old)
513 multi_fvm%SUBENER_FLUXES(imat, kface2_old, j_old) = -multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old)
514 ENDDO
515 ENDIF
516 ELSE
517 !PRINT*, "VOISIN DOMAINE"
518 ENDIF
519! -----------------------------------------------
520 ELSE IF(j_old <= 0) THEN
521 elemid = i_old
522 vx = multi_fvm%VEL(1, elemid)
523 vy = multi_fvm%VEL(2, elemid)
524 vz = multi_fvm%VEL(3, elemid)
525 normvel = vx * nx + vy * ny + vz * nz
526 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
527 ssp = multi_fvm%SOUND_SPEED(elemid)
528 rho = multi_fvm%RHO(elemid)
529C HLL wave speed estimates
530 sl = min(normvel - ssp, two * normalw - normvel - ssp)
531 sr = max(normvel + ssp, two * normalw - normvel + ssp)
532
533C Intermediate wave speed
534 sstar = normalw
535
536 p = multi_fvm%PRES(elemid)
537 pstar = p + rho * (sstar - normvel) * (sl - normvel)
538
539 pp(1) = zero
540 pp(2) = pstar * nx
541 pp(3) = pstar * ny
542 pp(4) = pstar * nz
543 pp(5) = sstar * (pstar + pshift)
544 multi_fvm%FLUXES(1:5, kface3, elemid) = surf * pp(1:5)
545 multi_fvm%FLUXES(6, kface3, elemid) = surf * normalw
546C ===
547C Submaterial fluxes
548 IF (nbmat > 1) THEN
549 DO imat = 1, nbmat
550 multi_fvm%SUBVOL_FLUXES(imat, kface3, i_old) = zero
551 multi_fvm%SUBMASS_FLUXES(imat, kface3, i_old) = zero
552 multi_fvm%SUBENER_FLUXES(imat, kface3, i_old) = zero
553 ENDDO
554 ENDIF
555 ENDIF
556! -----------------------------------------------
557 ENDDO !KFACE3
558 ENDDO ! II = 1, NEL
559
560
561C----------------------
562C Boundary fluxes
563C----------------------
564 END SUBROUTINE multi_fluxes_computation
565
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
multi_fluxes_computation
subroutine multi_fluxes_computation(ng, elbuf_tab, iparg, itask, ixs, ixq, ixtg, pm, ipm, multi_fvm, ale_connectivity, wgrid, xgrid, itab, nbmat, current_time, bufmat, id_global_vois, npf, tf)
Definition multi_fluxes_computation.F:39
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
initbuf_mod
Definition initbuf.F:160
arret
subroutine arret(nn)
Definition arret.F:86
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256