OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_muscl_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_muscl_fluxes_computation_mod ../engine/source/multifluid/multi_muscl_fluxes_computation.F
25!||--- called by ------------------------------------------------------
26!|| multi_timeevolution ../engine/source/multifluid/multi_timeevolution.F
27!||====================================================================
28 MODULE multi_muscl_fluxes_computation_mod
29 CONTAINS
30!||====================================================================
31!|| multi_muscl_fluxes_computation ../engine/source/multifluid/multi_muscl_fluxes_computation.F
32!||--- called by ------------------------------------------------------
33!|| multi_timeevolution ../engine/source/multifluid/multi_timeevolution.F
34!||--- calls -----------------------------------------------------
35!|| arret ../engine/source/system/arret.F
36!|| multi_muscl_compute_pressure ../engine/source/multifluid/multi_muscl_compute_pressure.F90
37!|| my_flush ../engine/source/system/machine.F
38!||--- uses -----------------------------------------------------
39!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
40!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
41!|| element_mod ../common_source/modules/elements/element_mod.F90
42!|| initbuf_mod ../engine/share/resol/initbuf.F
43!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
44!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
45!|| multi_muscl_compute_pressure_mod ../engine/source/multifluid/multi_muscl_compute_pressure.f90
46!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.F
47!||====================================================================
48 SUBROUTINE multi_muscl_fluxes_computation(
49 . NG, ELBUF_TAB, IPARG, ITASK,
50 . IXS, IXQ, IXTG,
51 . PM, IPM, MULTI_FVM, ALE_CONNECTIVITY, WGRID, XGRID, ITAB, NBMAT, CURRENT_TIME, BUFMAT,
52 . ID_GLOBAL_VOIS,FACE_VOIS,NPF,TF,ISPMD,MATPARAM)
53C-----------------------------------------------
54!$COMMENT
55! MULTI_MUSCL_FLUXES_COMPUTATION description :
56! computation of fluxes with 2nd order algorithm
57! MULTI_MUSCL_FLUXES_COMPUTATION organization :
58! The parith/on is ensured by the same
59! order computation :
60! * check the user ID
61! * if user ID( element1) < user ID( element2)
62! --> element1 drives the computation
63! * if user ID( element1) > user ID( element2)
64! --> element2 drives the computation
65!
66!$ENDCOMMENT
67C-----------------------------------------------
68C M o d u l e s
69C-----------------------------------------------
70 USE initbuf_mod
71 USE elbufdef_mod
72 USE multi_fvm_mod
73 USE ale_connectivity_mod
74 USE matparam_def_mod, ONLY : matparam_struct_
75 USE multi_submatlaw_mod , ONLY : multi_submatlaw
76 USE multi_muscl_compute_pressure_mod , ONLY : multi_muscl_compute_pressure
77 use element_mod , only : nixs,nixq,nixtg
78C-----------------------------------------------
79C I m p l i c i t T y p e s
80C-----------------------------------------------
81#include "implicit_f.inc"
82C-----------------------------------------------
83C C o m m o n B l o c k s
84C-----------------------------------------------
85! NUMELS :
86#include "com01_c.inc"
87#include "com04_c.inc"
88#include "param_c.inc"
89#include "comlock.inc"
90! DTFAC1:
91#include "mvsiz_p.inc"
92#include "units_c.inc"
93 COMMON /tablesizf/ stf,snpc
94 INTEGER STF,SNPC
95C-----------------------------------------------
96C D u m m y A r g u m e n t s
97C-----------------------------------------------
98 TYPE(matparam_struct_),DIMENSION(NUMMAT),INTENT(IN) :: MATPARAM !< material buffer
99 INTEGER, INTENT(IN) :: NG, ISPMD
100 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
101 INTEGER, INTENT(IN) :: IPARG(NPARG, NGROUP)
102 INTEGER, INTENT(IN) :: ITASK ! SMP TASK
103 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
104 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT)
105 my_real, INTENT(IN) :: pm(npropm, nummat)
106 TYPE(multi_fvm_struct), INTENT(INOUT) :: MULTI_FVM
107 ! for parith/on : ID_GLOBAL_VOIS --> user id ; FACE_VOIS --> face of the remote element
108 INTEGER, INTENT(IN) :: ID_GLOBAL_VOIS(*),FACE_VOIS(*)
109 my_real, INTENT(IN) :: wgrid(3, *), xgrid(3, *), current_time
110 INTEGER, INTENT(IN) :: ITAB(*), NBMAT
111 my_real, INTENT(INOUT) :: BUFMAT(*)
112 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
113 INTEGER, INTENT(IN) :: NPF(SNPC)
114 my_real, INTENT(IN) :: tf(stf)
115C-----------------------------------------------
116C L o c a l V a r i a b l e s
117C-----------------------------------------------
118 INTEGER :: II, JJ, KFACE, I, J, KFACE2, KK, ID_OUT, IARRAY(2)
119 INTEGER :: ITY, NFT, NEL
120 INTEGER :: IMAT
121 INTEGER :: KFACE3,I_OLD,J_OLD,KFACE_OLD,KFACE2_OLD
122 INTEGER :: ID_GLOBAL_VOIS_1,ID_GLOBAL_VOIS_2
123 INTEGER, DIMENSION(MVSIZ) :: GLOBAL_ID_CURRENT_ELM
124 INTEGER :: IAD
125 INTEGER :: NB_FACE, ELEMID, NUMEL_SPMD
126 INTEGER :: LOCAL_MATID(NBMAT), MATLAW(NBMAT),IDX(6), NVAREOS(NBMAT),NVARTMP_EOS(NBMAT)
127 my_real :: fii(5), normal_velii, normal_veljj, vii(5), velii2, veljj2, surf
128 my_real :: sl, sr, sstar, pstar, fiistar(5), sub_fiistar(3), viistar(5), pp(5),sub_viistar(3)
129 my_real :: sspii(1), sspjj(1), rhoii(1), rhojj(1), pii(1), pjj(1), nx, ny, nz, normalw
130 my_real :: vxii, vyii, vzii
131 my_real :: vxjj, vyjj, vzjj, vx, vy, vz, normvel, ssp, rho, p
132 my_real :: machii, machjj, pstar2, theta, alphaii
133 my_real :: sub_rhoii, sub_rhoeintii,alphastar, sub_rhostar, sub_estar, sub_pii, wfac(3)
134 my_real :: xf(3), xk(3), xl(3)
135 my_real :: eintii(1), eintjj(1)
136 my_real :: phase_rhoii(nbmat), phase_presii(nbmat), phase_eintii(nbmat)
137 my_real :: phase_sspii(nbmat), phase_alphaii(nbmat)
138 my_real :: phase_rhojj(nbmat), phase_presjj(nbmat), phase_eintjj(nbmat)
139 my_real :: phase_sspjj(nbmat), phase_alphajj(nbmat), cumul_alpha, pshift
140 my_real :: direction
141 TYPE(lbuf_ptr) :: LBUFS(NBMAT)
142 TYPE(ebuf_ptr) :: EBUFS(NBMAT)
143 TYPE(g_bufel_), POINTER :: GBUF
144 TYPE(buf_eos_), POINTER :: EBUF
145C-----------------------------------------------
146C S o u r c e L i n e s
147C-----------------------------------------------
148 nb_face = -1
149
150 gbuf => elbuf_tab(ng)%GBUF
151 nel = iparg(2, ng)
152 nft = iparg(3, ng)
153 ity = iparg(5, ng)
154
155 pshift = multi_fvm%PRES_SHIFT
156
157 DO i=1,6
158 idx(i) = nel*(i-1)
159 ENDDO
160
161 IF (nbmat > 1) THEN
162!DIR$ NOVECTOR
163 DO imat = 1, nbmat
164 lbufs(imat)%LBUF => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
165 ebufs(imat)%EBUF => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
166 nvareos(imat) = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
167 nvartmp_eos(imat) = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
168 ENDDO
169 ELSE
170 DO imat = 1, nbmat
171 !GBUF will be used
172 NULLIFY( lbufs(imat)%LBUF, ebufs(imat)%EBUF )
173 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1,1,1)
174 nvareos(1) = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
175 nvartmp_eos(1) = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS
176 ENDDO
177 ENDIF
178C Normal and surface computation
179 SELECT CASE (multi_fvm%SYM)
180 CASE (0)
181 nb_face = 6
182C 3D
183 numel_spmd = numels
184 DO imat = 1, nbmat
185 local_matid(imat) = matparam(ixs(1, 1 + nft))%MULTIMAT%MID(imat)
186 ENDDO
187 DO ii = 1, nel
188 i = ii + nft
189 global_id_current_elm(ii) = ixs(nixs,i)
190 ENDDO
191 CASE (1, 2)
192 IF (ity == 2) THEN
193C QUADS
194 nb_face = 4
195 numel_spmd = numelq
196 DO imat = 1, nbmat
197 local_matid(imat) = matparam(ixq(1, 1 + nft))%MULTIMAT%MID(imat)
198 ENDDO
199 DO ii = 1, nel
200 i = ii + nft
201 global_id_current_elm(ii) = ixq(nixq,i)
202 ENDDO
203 ELSEIF (ity == 7) THEN
204C TRIANGLES
205 nb_face = 3
206 numel_spmd = numeltg
207 DO imat = 1, nbmat
208 local_matid(imat) = matparam(ixtg(1, 1 + nft))%MULTIMAT%MID(imat)
209 ENDDO
210 DO ii = 1, nel
211 i = ii + nft
212 global_id_current_elm(ii) = ixtg(nixtg,i)
213 ENDDO
214 ENDIF
215 CASE DEFAULT
216 CALL arret(2)
217 END SELECT
218
219C Computation flow
220 DO ii = 1, nel
221 i = ii + nft
222 iad = ale_connectivity%ee_connect%iad_connect(i)
223 nb_face = ale_connectivity%ee_connect%iad_connect(i+1) -
224 . ale_connectivity%ee_connect%iad_connect(i)
225 i_old = i
226C Face KFACE
227 DO kface3 = 1, nb_face
228 jj = ale_connectivity%ee_connect%connected(iad + kface3 - 1)
229 j_old = jj
230C Face normal
231 nx = multi_fvm%FACE_DATA%NORMAL(1, kface3, i_old)
232 ny = multi_fvm%FACE_DATA%NORMAL(2, kface3, i_old)
233 nz = multi_fvm%FACE_DATA%NORMAL(3, kface3, i_old)
234 surf = multi_fvm%FACE_DATA%SURF(kface3, i_old)
235 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface3, i_old)
236
237 id_global_vois_1 = global_id_current_elm(ii)
238 id_global_vois_2 = id_global_vois( nb_face * (i_old - 1) + kface3 )
239 kface_old = kface3
240 IF( j_old>0 ) THEN
241 IF (j_old <= numel_spmd) THEN
242 kface2_old = multi_fvm%FVM_CONNECTIVITY%KVOIS(nb_face * (i_old - 1) + kface3)
243 ELSE
244 kface2_old = face_vois( nb_face * (i_old - 1) + kface3 )
245 ENDIF
246 ENDIF
247
248! -----------------------------------------------
249 IF (j_old > 0 .AND. i_old < j_old) THEN
250 ! ----------------------
251
252 IF(id_global_vois_1<id_global_vois_2) THEN
253 direction = one
254 kface = kface_old
255 kface2 = kface2_old
256 i = i_old
257 jj = j_old
258 ELSE
259 direction = -one
260 kface2 = kface_old
261 kface = kface2_old
262 i = j_old
263 jj = i_old
264 ENDIF
265C Face normal
266 nx = nx *direction
267 ny = ny *direction
268 nz = nz *direction
269
270 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
271 ! ----------------------
272
273 j = jj
274! Conserved variable current element
275C ============
276C MUSCL
277C ============
278 xk(1:3) = multi_fvm%ELEM_DATA%CENTROID(1:3, i)
279 xl(1:3) = multi_fvm%ELEM_DATA%CENTROID(1:3, j)
280 xf(1:3) = multi_fvm%FACE_DATA%CENTROID(1:3, kface, i) ! <-- need to commnicate this value
281
282C Reconstructed velocity current element
283 vxii = multi_fvm%VEL(1, i)
284 vyii = multi_fvm%VEL(2, i)
285 vzii = multi_fvm%VEL(3, i)
286
287 IF (multi_fvm%MUSCL == 1) THEN
288 vxii = vxii + multi_fvm%GRAD_U(1, i) * (xf(1) - xk(1))
289 . + multi_fvm%GRAD_U(2, i) * (xf(2) - xk(2))
290 . + multi_fvm%GRAD_U(3, i) * (xf(3) - xk(3))
291 vyii = vyii + multi_fvm%GRAD_V(1, i) * (xf(1) - xk(1))
292 . + multi_fvm%GRAD_V(2, i) * (xf(2) - xk(2))
293 . + multi_fvm%GRAD_V(3, i) * (xf(3) - xk(3))
294 vzii = vzii + multi_fvm%GRAD_W(1, i) * (xf(1) - xk(1))
295 . + multi_fvm%GRAD_W(2, i) * (xf(2) - xk(2))
296 . + multi_fvm%GRAD_W(3, i) * (xf(3) - xk(3))
297 ENDIF
298
299C Squared velocity
300 velii2 = vxii * vxii + vyii * vyii + vzii * vzii
301C Normal velocity current element
302 normal_velii = vxii * nx + vyii * ny + vzii * nz
303C Reconstructed internal energy current element
304 IF (multi_fvm%NBMAT == 1) THEN
305C Reconstructed density current element
306 rhoii(1) = multi_fvm%RHO(i)
307 eintii(1) = multi_fvm%EINT(i)
308
309 IF (multi_fvm%MUSCL == 1) THEN
310 rhoii(1) = rhoii(1) + multi_fvm%GRAD_RHO(1, i) * (xf(1) - xk(1))
311 . + multi_fvm%GRAD_RHO(2, i) * (xf(2) - xk(2))
312 . + multi_fvm%GRAD_RHO(3, i) * (xf(3) - xk(3))
313 eintii(1) = eintii(1) + multi_fvm%GRAD_PRES(1, i) * (xf(1) - xk(1))
314 . + multi_fvm%GRAD_PRES(2, i) * (xf(2) - xk(2))
315 . + multi_fvm%GRAD_PRES(3, i) * (xf(3) - xk(3))
316 ENDIF
317 matlaw(1) = ipm(2, local_matid(1))
318 CALL multi_muscl_compute_pressure(matlaw(1), local_matid(1), pm, ipm, npropm, npropmi,
319 . eintii(1), rhoii(1), pii(1), sspii(1),
320 . multi_fvm%BFRAC(1, i), gbuf%TB(ii:), gbuf%DELTAX(ii:), current_time,
321 . bufmat, gbuf%OFF(ii:), gbuf%SIG(ii+idx(1):ii+idx(6)), snpc, stf, npf, tf, ebuf%VAR, nvareos(1),
322 . matparam(local_matid(1)) ,nvartmp_eos(1), ebuf%VARTMP,nummat,gbuf%ABURN(ii:))
323 sspii = sqrt(sspii)
324
325 ELSE
326 rhoii(1) = zero
327 pii(1) = zero
328 !EINTII = MULTI_FVM%EINT(I)
329 sspii(1) = zero
330 eintii(1) = zero
331 cumul_alpha = zero
332 DO imat = 1, nbmat
333 phase_alphaii(imat) = multi_fvm%PHASE_ALPHA(imat, i)
334 . + multi_fvm%PHASE_GRAD_ALPHA(1, imat, i) * (xf(1) - xk(1))
335 . + multi_fvm%PHASE_GRAD_ALPHA(2, imat, i) * (xf(2) - xk(2))
336 . + multi_fvm%PHASE_GRAD_ALPHA(3, imat, i) * (xf(3) - xk(3))
337 phase_alphaii(imat) = max(phase_alphaii(imat), zero)
338 cumul_alpha = cumul_alpha + phase_alphaii(imat)
339 ENDDO
340
341 DO imat = 1, nbmat
342 phase_alphaii(imat) = phase_alphaii(imat) / cumul_alpha
343 phase_rhoii(imat) = multi_fvm%PHASE_RHO(imat, i)
344 phase_eintii(imat) = multi_fvm%PHASE_EINT(imat, i)
345 IF (multi_fvm%MUSCL == 1) THEN
346 phase_rhoii(imat) = phase_rhoii(imat)
347 . + multi_fvm%PHASE_GRAD_RHO(1, imat, i) * (xf(1) - xk(1))
348 . + multi_fvm%PHASE_GRAD_RHO(2, imat, i) * (xf(2) - xk(2))
349 . + multi_fvm%PHASE_GRAD_RHO(3, imat, i) * (xf(3) - xk(3))
350 phase_eintii(imat) = phase_eintii(imat)
351 . + multi_fvm%PHASE_GRAD_PRES(1, imat, i) * (xf(1) - xk(1))
352 . + multi_fvm%PHASE_GRAD_PRES(2, imat, i) * (xf(2) - xk(2))
353 . + multi_fvm%PHASE_GRAD_PRES(3, imat, i) * (xf(3) - xk(3))
354 ENDIF
355C Global reconstructed density
356 rhoii(1) = rhoii(1) + phase_alphaii(imat) * phase_rhoii(imat)
357 matlaw(imat) = ipm(2, local_matid(imat))
358 IF (phase_alphaii(imat) > zero) THEN
359 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
360 CALL multi_muscl_compute_pressure(matlaw(imat), local_matid(imat), pm, ipm, npropm, npropmi,
361 . phase_eintii(imat), phase_rhoii(imat), phase_presii(imat), phase_sspii(imat),
362 . multi_fvm%BFRAC(imat, i), gbuf%TB(ii:), gbuf%DELTAX(ii:), current_time,
363 . bufmat, lbufs(imat)%LBUF%OFF(ii:),lbufs(imat)%LBUF%SIG(ii+idx(1):ii+idx(6)), snpc,stf,npf,tf,
364 . ebufs(imat)%EBUF%VAR,nvareos(imat), matparam(local_matid(imat)),
365 . nvartmp_eos(imat),ebufs(imat)%EBUF%VARTMP,nummat,gbuf%ABURN(ii:))
366C Global reconstructed pressure
367 pii(1) = pii(1) + phase_alphaii(imat) * phase_presii(imat)
368 eintii(1) = eintii(1) + phase_alphaii(imat) * phase_eintii(imat)
369 sspii(1) = sspii(1) + phase_alphaii(imat) * phase_rhoii(imat) * max(em20, phase_sspii(imat))
370 ENDIF
371 ENDDO
372 IF (sspii(1) / rhoii(1) > zero) THEN
373 sspii(1) = sqrt(sspii(1) / rhoii(1))
374 ELSE
375 sspii(1) = multi_fvm%SOUND_SPEED(i)
376 ENDIF
377 ENDIF
378
379C Reconstructed velocity adjacent element
380 vxjj = multi_fvm%VEL(1, j)
381 vyjj = multi_fvm%VEL(2, j)
382 vzjj = multi_fvm%VEL(3, j)
383 IF (multi_fvm%MUSCL == 1) THEN
384 vxjj = vxjj
385 . + multi_fvm%GRAD_U(1, j) * (xf(1) - xl(1))
386 . + multi_fvm%GRAD_U(2, j) * (xf(2) - xl(2))
387 . + multi_fvm%GRAD_U(3, j) * (xf(3) - xl(3))
388 vyjj = vyjj
389 . + multi_fvm%GRAD_V(1, j) * (xf(1) - xl(1))
390 . + multi_fvm%GRAD_V(2, j) * (xf(2) - xl(2))
391 . + multi_fvm%GRAD_V(3, j) * (xf(3) - xl(3))
392 vzjj = vzjj
393 . + multi_fvm%GRAD_W(1, j) * (xf(1) - xl(1))
394 . + multi_fvm%GRAD_W(2, j) * (xf(2) - xl(2))
395 . + multi_fvm%GRAD_W(3, j) * (xf(3) - xl(3))
396 ENDIF
397C Normal velocity adjacent element
398 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
399 veljj2 = vxjj * vxjj + vyjj * vyjj + vzjj * vzjj
400
401C Reconstructed internal energy adjacent element
402 IF (multi_fvm%NBMAT == 1) THEN
403C Reconstructed density adjacent element
404 rhojj(1) = multi_fvm%RHO(j)
405 eintjj(1) = multi_fvm%EINT(j)
406 IF (multi_fvm%MUSCL == 1) THEN
407 rhojj(1) = rhojj(1)
408 . + multi_fvm%GRAD_RHO(1, j) * (xf(1) - xl(1))
409 . + multi_fvm%GRAD_RHO(2, j) * (xf(2) - xl(2))
410 . + multi_fvm%GRAD_RHO(3, j) * (xf(3) - xl(3))
411 eintjj(1) = eintjj(1)
412 . + multi_fvm%GRAD_PRES(1, j) * (xf(1) - xl(1))
413 . + multi_fvm%GRAD_PRES(2, j) * (xf(2) - xl(2))
414 . + multi_fvm%GRAD_PRES(3, j) * (xf(3) - xl(3))
415 ENDIF
416 matlaw(1) = ipm(2, local_matid(1))
417 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
418 CALL multi_muscl_compute_pressure(matlaw(1), local_matid(1), pm, ipm, npropm, npropmi,
419 . eintjj(1), rhojj(1), pjj(1), sspjj(1),
420 . multi_fvm%BFRAC(1, j), gbuf%TB(ii:), gbuf%DELTAX(ii:), current_time,
421 . bufmat, gbuf%OFF(ii:),gbuf%SIG(ii+idx(1):ii+idx(6)), snpc,stf,npf,tf,ebuf%VAR,nvareos(1),
422 . matparam(local_matid(1)),
423 . nvartmp_eos(1),ebuf%VARTMP,nummat,gbuf%ABURN(ii:))
424 sspjj(1) = sqrt(sspjj(1))
425 ELSE
426 rhojj(1) = zero
427 pjj(1) = zero
428 eintjj(1) = zero
429 sspjj(1) = zero
430 cumul_alpha = zero
431 DO imat = 1, nbmat
432 phase_alphajj(imat) = multi_fvm%PHASE_ALPHA(imat, j)
433 . + multi_fvm%PHASE_GRAD_ALPHA(1, imat, j) * (xf(1) - xl(1))
434 . + multi_fvm%PHASE_GRAD_ALPHA(2, imat, j) * (xf(2) - xl(2))
435 . + multi_fvm%PHASE_GRAD_ALPHA(3, imat, j) * (xf(3) - xl(3))
436 phase_alphajj(imat) = max(phase_alphajj(imat), zero)
437 cumul_alpha = cumul_alpha + phase_alphajj(imat)
438 ENDDO
439
440 DO imat = 1, nbmat
441 phase_alphajj(imat) = phase_alphajj(imat) / cumul_alpha
442 phase_rhojj(imat) = multi_fvm%PHASE_RHO(imat, j)
443 phase_eintjj(imat) = multi_fvm%PHASE_EINT(imat, j)
444 IF (multi_fvm%MUSCL == 1) THEN
445 phase_rhojj(imat) = phase_rhojj(imat)
446 . + multi_fvm%PHASE_GRAD_RHO(1, imat, j) * (xf(1) - xl(1))
447 . + multi_fvm%PHASE_GRAD_RHO(2, imat, j) * (xf(2) - xl(2))
448 . + multi_fvm%PHASE_GRAD_RHO(3, imat, j) * (xf(3) - xl(3))
449 phase_eintjj(imat) = phase_eintjj(imat)
450 . + multi_fvm%PHASE_GRAD_PRES(1, imat, j) * (xf(1) - xl(1))
451 . + multi_fvm%PHASE_GRAD_PRES(2, imat, j) * (xf(2) - xl(2))
452 . + multi_fvm%PHASE_GRAD_PRES(3, imat, j) * (xf(3) - xl(3))
453 ENDIF
454C Global reconstructed density
455 rhojj(1) = rhojj(1) + phase_alphajj(imat) * phase_rhojj(imat)
456 matlaw(imat) = ipm(2, local_matid(imat))
457 IF (phase_alphajj(imat) > zero) THEN
458
459 CALL multi_muscl_compute_pressure(matlaw(imat), local_matid(imat), pm, ipm, npropm, npropmi,
460 . phase_eintjj(imat), phase_rhojj(imat), phase_presjj(imat), phase_sspjj(imat),
461 . multi_fvm%BFRAC(imat, j), gbuf%TB(ii:), gbuf%DELTAX(ii:), current_time,
462 . bufmat, lbufs(imat)%LBUF%OFF(ii:),lbufs(imat)%LBUF%SIG(ii+idx(1):ii+idx(6)),snpc,stf,npf,tf,
463 . ebufs(imat)%EBUF%VAR,nvareos(imat), matparam(local_matid(imat)),
464 . nvartmp_eos(imat),ebufs(imat)%EBUF%VARTMP,nummat,gbuf%ABURN(ii:))
465C Global reconstructed pressure
466 pjj(1) = pjj(1) + phase_alphajj(imat) * phase_presjj(imat)
467 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
468 sspjj(1) = sspjj(1) + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
469 ENDIF
470 ENDDO
471 IF (sspjj(1) / rhojj(1) > zero) THEN
472 sspjj(1) = sqrt(sspjj(1) / rhojj(1))
473 ELSE
474 sspjj(1) = multi_fvm%SOUND_SPEED(j)
475 ENDIF
476 ENDIF
477C Local Mach numbers
478 machii = abs(normal_velii) / sspii(1)
479 machjj = abs(normal_veljj) / sspjj(1)
480
481C HLL wave speed estimates
482 sl = min(normal_velii - sspii(1), normal_veljj - sspjj(1))
483 sr = max(normal_velii + sspii(1), normal_veljj + sspjj(1))
484
485C Intermediate wave speed
486 sstar = pjj(1) - pii(1) + rhoii(1) * normal_velii * (sl - normal_velii) -
487 . rhojj(1) * normal_veljj * (sr - normal_veljj)
488 sstar = sstar / (rhoii(1) * (sl - normal_velii) -
489 . rhojj(1) * (sr - normal_veljj))
490C Specific for Low Mach number corrections
491C Intermediate pressure
492 pstar2 = pii(1) + rhoii(1) * (sstar - normal_velii) * (sl - normal_velii)
493 IF (min(machii, machjj) < em01) THEN
494 theta = min(machii, machjj)
495 ELSE
496 theta = one
497 ENDIF
498 pstar = (one - theta) * half * (pii(1) + pjj(1)) + theta * pstar2
499 IF (multi_fvm%LOWMACH_OPT) THEN
500 pp(1) = zero
501 pp(2) = pstar * nx
502 pp(3) = pstar * ny
503 pp(4) = pstar * nz
504 pp(5) = sstar * (pstar + pshift)
505 ELSE
506 pp(1) = zero
507 pp(2) = pstar2 * nx
508 pp(3) = pstar2 * ny
509 pp(4) = pstar2 * nz
510 pp(5) = sstar * (pstar2 + pshift)
511 ENDIF
512
513 IF (sl > normalw) THEN
514 vii(1) = rhoii(1)
515 vii(2) = rhoii(1) * vxii
516 vii(3) = rhoii(1) * vyii
517 vii(4) = rhoii(1) * vzii
518 vii(5) = eintii(1) + half * rhoii(1) * velii2
519! Normal physical flux current element
520 fii(1) = vii(1) * normal_velii
521 fii(2) = vii(2) * normal_velii + pii(1) * nx
522 fii(3) = vii(3) * normal_velii + pii(1) * ny
523 fii(4) = vii(4) * normal_velii + pii(1) * nz
524 fii(5) = (vii(5) + pii(1) + pshift) * normal_velii
525C Take the fluxes of cell II
526C ===
527C Global fluxes
528 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction * (fii(1:5) - normalw * vii(1:5)) * surf
529 multi_fvm%FLUXES(6, kface_old, i_old) = direction * normal_velii * surf
530C ===
531C Submaterial fluxes
532 IF (nbmat > 1) THEN
533 DO imat = 1, nbmat
534 alphaii = phase_alphaii(imat)
535 sub_rhoii = phase_rhoii(imat) ! ALPHA_RHO
536 sub_rhoeintii = phase_eintii(imat)
537 sub_viistar(1) = alphaii
538 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
539 sub_viistar(3) = alphaii * sub_rhoeintii
540 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
541 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
542 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
543 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
544 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
545 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
546 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
547 ENDDO
548 ENDIF
549C
550 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
551 vii(1) = rhoii(1)
552 vii(2) = rhoii(1) * vxii
553 vii(3) = rhoii(1) * vyii
554 vii(4) = rhoii(1) * vzii
555 vii(5) = eintii(1) + half * rhoii(1) * velii2
556! Normal physical flux current element
557 fii(1) = vii(1) * normal_velii
558 fii(2) = vii(2) * normal_velii + pii(1) * nx
559 fii(3) = vii(3) * normal_velii + pii(1) * ny
560 fii(4) = vii(4) * normal_velii + pii(1) * nz
561 fii(5) = (vii(5) + pii(1) + pshift) * normal_velii
562C Take intermediate state flux (HLLC scheme)
563C ===
564C Global fluxes
565 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
566 viistar(1:5) = viistar(1:5) / (sstar - sl)
567 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
568 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction *
569 . (fiistar(1:5) - normalw * viistar(1:5)) * surf
570 multi_fvm%FLUXES(6, kface_old, i_old) = direction * sstar * surf
571C ===
572C Submaterial fluxes
573 IF (nbmat > 1) THEN
574 DO imat = 1, nbmat
575 matlaw(imat) = ipm(2, local_matid(imat))
576 alphastar = phase_alphaii(imat)
577 sub_rhostar = phase_rhoii(imat) * (normal_velii - sl) / (sstar - sl)
578 IF (alphastar > zero) THEN
579 sub_rhoii = phase_rhoii(imat)
580 sub_rhoeintii = phase_eintii(imat)
581 sub_pii = phase_presii(imat)
582 sub_estar = phase_eintii(imat) / phase_rhoii(imat) -
583 . phase_presii(imat) * (one / sub_rhostar - one / phase_rhoii(imat))
584
585 IF (sub_estar < zero) THEN
586 sub_estar = zero
587 ENDIF
588 ELSE
589 sub_estar = zero
590 ENDIF
591 sub_viistar(1) = alphastar
592 sub_viistar(2) = alphastar * sub_rhostar
593 sub_viistar(3) = alphastar * sub_rhostar * sub_estar
594
595 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
596 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
597 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
598 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
599 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
600 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
601 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
602 ENDDO
603 ENDIF
604
605 ELSE IF (sstar < normalw .AND. normalw <= sr) THEN
606 vii(1) = rhojj(1)
607 vii(2) = rhojj(1) * vxjj
608 vii(3) = rhojj(1) * vyjj
609 vii(4) = rhojj(1) * vzjj
610 vii(5) = eintjj(1) + half * rhojj(1) * veljj2
611! normal physical flux current element
612 fii(1) = vii(1) * normal_veljj
613 fii(2) = vii(2) * normal_veljj + pjj(1) * nx
614 fii(3) = vii(3) * normal_veljj + pjj(1) * ny
615 fii(4) = vii(4) * normal_veljj + pjj(1) * nz
616 fii(5) = (vii(5) + pjj(1) + pshift) * normal_veljj
617C Take intermediate state flux (HLLC scheme)
618C ===
619C Global fluxes
620 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
621 viistar(1:5) = viistar(1:5) / (sstar - sr)
622 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
623 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction *
624 . (fiistar(1:5) - normalw * viistar(1:5)) * surf
625 multi_fvm%FLUXES(6, kface_old, i_old) = direction * sstar * surf
626C ===
627C Submaterial fluxes
628 IF (nbmat > 1) THEN
629 DO imat = 1, nbmat
630 matlaw(imat) = ipm(2, local_matid(imat))
631 alphastar = phase_alphajj(imat)
632 sub_rhostar = phase_rhojj(imat) * (normal_veljj - sr) / (sstar - sr)
633 IF (alphastar > zero) THEN
634 sub_rhoii = phase_rhojj(imat)
635 sub_rhoeintii = phase_eintjj(imat)
636 sub_pii = phase_presjj(imat)
637 sub_estar = phase_eintjj(imat) / phase_rhojj(imat) -
638 . phase_presjj(imat) * (one / sub_rhostar - one / phase_rhojj(imat))
639 IF (sub_estar < zero) THEN
640 sub_estar = zero
641 ENDIF
642 ELSE
643 sub_estar = zero
644 ENDIF
645
646 sub_viistar(1) = phase_alphajj(imat)
647 sub_viistar(2) = phase_alphajj(imat) * sub_rhostar
648 sub_viistar(3) = phase_alphajj(imat) * sub_rhostar * sub_estar
649 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
650 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
651 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
652 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
653 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
654 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
655 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
656 ENDDO
657 ENDIF
658 ELSE IF (sr < normalw) THEN
659 vii(1) = rhojj(1)
660 vii(2) = rhojj(1) * vxjj
661 vii(3) = rhojj(1) * vyjj
662 vii(4) = rhojj(1) * vzjj
663 vii(5) = eintjj(1) + half * rhojj(1) * veljj2
664! Normal physical flux current element
665 fii(1) = vii(1) * normal_veljj
666 fii(2) = vii(2) * normal_veljj + pjj(1) * nx
667 fii(3) = vii(3) * normal_veljj + pjj(1) * ny
668 fii(4) = vii(4) * normal_veljj + pjj(1) * nz
669 fii(5) = (vii(5) + pjj(1) + pshift) * normal_veljj
670C Take the fluxes of cell II
671C ===
672C Global fluxes
673 multi_fvm%FLUXES(1:5, kface_old, i_old) = direction * (fii(1:5) - normalw * vii(1:5)) * surf
674 multi_fvm%FLUXES(6, kface_old, i_old) = direction * normal_veljj * surf
675C ===
676C Submaterial fluxes
677 IF (nbmat > 1) THEN
678 DO imat = 1, nbmat
679 alphaii = phase_alphajj(imat)
680 sub_rhoii = phase_rhojj(imat) ! ALPHA_RHO
681 sub_rhoeintii = phase_eintjj(imat)
682 sub_viistar(1) = alphaii
683 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
684 sub_viistar(3) = alphaii * sub_rhoeintii
685 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
686 multi_fvm%SUBVOL_FLUXES(imat, kface_old, i_old) = direction *
687 . (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
688 multi_fvm%SUBMASS_FLUXES(imat, kface_old, i_old) = direction *
689 . (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
690 multi_fvm%SUBENER_FLUXES(imat, kface_old, i_old) = direction *
691 . (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
692 ENDDO
693 ENDIF
694 ELSE
695
696 IF(ispmd == 0)THEN
697#include "lockon.inc"
698 iarray(1:2)=(/iout,istdo/)
699 DO kk=1,2
700 id_out=iarray(kk)
701 WRITE(id_out,*)" **error : MUSCL FLUXES CALCULATION"
702 WRITE(id_out,*)" potential instability with elem id = ", global_id_current_elm(ii)
703 WRITE(id_out,*)" face id = ", kface3
704 WRITE(id_out,*)" sound speed = ", sspii(1)
705 WRITE(id_out,*)" adjacent sound speed = ", sspjj(1)
706 WRITE(id_out,*)" normal velocity = ", normal_velii
707 WRITE(id_out,*)" adjacent normal velocity = ", normal_veljj
708 WRITE(id_out,*)""
709 WRITE(id_out,*)" recommendations :"
710 WRITE(id_out,*)" -1- check EoS parameters and unit system"
711 WRITE(id_out,*)" -2- check boundary conditions (type, parameters, and proximity from important gradients);"
712 WRITE(id_out,*)
713 . " -3- update beta parameters in /ale/muscl option(find compromise between precision and stability)"
714 ENDDO
715 CALL MY_FLUSH(IOUT)
716 CALL MY_FLUSH(ISTDO)
717 CALL ARRET(2)
718#include "lockoff.inc"
719 ENDIF!ISPMD==0
720
721 ENDIF
722
723 IF (J_OLD <= NUMEL_SPMD) THEN
724! KFACE2_OLD = MULTI_FVM%FVM_CONNECTIVITY%KVOIS(NB_FACE * (I_OLD - 1) + KFACE3)
725 MULTI_FVM%FLUXES(1:6, KFACE2_OLD, J_OLD) = -MULTI_FVM%FLUXES(1:6, KFACE_OLD, I_OLD)
726 IF (NBMAT > 1) THEN
727 DO IMAT = 1, NBMAT
728 MULTI_FVM%SUBVOL_FLUXES(IMAT, KFACE2_OLD, J_OLD) = -MULTI_FVM%SUBVOL_FLUXES(IMAT, KFACE_OLD, I_OLD)
729 MULTI_FVM%SUBMASS_FLUXES(IMAT, KFACE2_OLD, J_OLD) = -MULTI_FVM%SUBMASS_FLUXES(IMAT, KFACE_OLD, I_OLD)
730 MULTI_FVM%SUBENER_FLUXES(IMAT, KFACE2_OLD, J_OLD) = -MULTI_FVM%SUBENER_FLUXES(IMAT, KFACE_OLD, I_OLD)
731 ENDDO
732 ENDIF
733 ENDIF
734! -----------------------------------------------
735 ELSE IF (J_OLD <= 0) THEN
736 ELEMID = I_OLD
737 VX = MULTI_FVM%VEL(1, ELEMID)
738 VY = MULTI_FVM%VEL(2, ELEMID)
739 VZ = MULTI_FVM%VEL(3, ELEMID)
740 NORMVEL = VX * NX + VY * NY + VZ * NZ
741 NORMALW = WFAC(1) * NX + WFAC(2) * NY + WFAC(3) * NZ
742 SSP = MULTI_FVM%SOUND_SPEED(ELEMID)
743 RHO = MULTI_FVM%RHO(ELEMID)
744C HLL wave speed estimates
745 SL = MIN(NORMVEL - SSP, TWO * NORMALW - NORMVEL - SSP)
746 SR = MAX(NORMVEL + SSP, TWO * NORMALW - NORMVEL + SSP)
747
748C Intermediate wave speed
749 SSTAR = NORMALW
750
751 P = MULTI_FVM%PRES(ELEMID)
752 PSTAR = P + RHO * (SSTAR - NORMVEL) * (SL - NORMVEL)
753
754 PP(1) = ZERO
755 PP(2) = PSTAR * NX
756 PP(3) = PSTAR * NY
757 PP(4) = PSTAR * NZ
758 PP(5) = SSTAR * (PSTAR + PSHIFT)
759 MULTI_FVM%FLUXES(1:5, KFACE3, ELEMID) = SURF * PP(1:5)
760 MULTI_FVM%FLUXES(6, KFACE3, ELEMID) = SURF * NORMALW
761C ===
762C Submaterial fluxes
763 IF (NBMAT > 1) THEN
764 DO IMAT = 1, NBMAT
765 MULTI_FVM%SUBVOL_FLUXES(IMAT, KFACE3, I_OLD) = ZERO
766 MULTI_FVM%SUBMASS_FLUXES(IMAT, KFACE3, I_OLD) = ZERO
767 MULTI_FVM%SUBENER_FLUXES(IMAT, KFACE3, I_OLD) = ZERO
768 ENDDO
769 ENDIF
770 ENDIF
771! -----------------------------------------------
772 ENDDO !KFACE
773 ENDDO ! II = 1, NEL
774
775
776C----------------------
777C Boundary fluxes
778C----------------------
779 END SUBROUTINE MULTI_MUSCL_FLUXES_COMPUTATION
780 END MODULE MULTI_MUSCL_FLUXES_COMPUTATION_MOD
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
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
initbuf_mod
Definition initbuf.F:160
multi_muscl_fluxes_computation_mod
Definition multi_muscl_fluxes_computation.F:28
multi_muscl_fluxes_computation_mod::multi_muscl_fluxes_computation
subroutine multi_muscl_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, face_vois, npf, tf, ispmd, matparam)
Definition multi_muscl_fluxes_computation.F:53
multi_submatlaw_mod
Definition multi_submatlaw.F:33
multi_submatlaw_mod::multi_submatlaw
subroutine multi_submatlaw(iflag, matlaw, local_matid, nel, eint, pres, rho, ssp, vol, grun, pm, ipm, npropm, npropmi, bufmat, off, theta, burnfrac, burntime, deltax, current_time, sigold, snpc, stf, npf, tf, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos, nummat, aburn)
Definition multi_submatlaw.F:61
arret
subroutine arret(nn)
Definition arret.F:86
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256