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