OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_muscl_fluxes_computation_mod Module Reference

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ multi_muscl_fluxes_computation()

subroutine multi_muscl_fluxes_computation_mod::multi_muscl_fluxes_computation ( integer, intent(in) ng,
type(elbuf_struct_), dimension(ngroup), target elbuf_tab,
integer, dimension(nparg, ngroup), intent(in) iparg,
integer, intent(in) itask,
integer, dimension(nixs, *), intent(in) ixs,
integer, dimension(nixq, *), intent(in) ixq,
integer, dimension(nixtg, *), intent(in) ixtg,
dimension(npropm, nummat), intent(in) pm,
integer, dimension(npropmi, nummat), intent(in) ipm,
type(multi_fvm_struct), intent(inout) multi_fvm,
type(t_ale_connectivity), intent(in) ale_connectivity,
dimension(3, *), intent(in) wgrid,
dimension(3, *), intent(in) xgrid,
integer, dimension(*), intent(in) itab,
integer, intent(in) nbmat,
intent(in) current_time,
dimension(*), intent(inout) bufmat,
integer, dimension(*), intent(in) id_global_vois,
integer, dimension(*), intent(in) face_vois,
integer, dimension(snpc), intent(in) npf,
dimension(stf), intent(in) tf,
integer, intent(in) ispmd,
type(matparam_struct_), dimension(nummat), intent(in) matparam )
Parameters
[in]matparammaterial buffer

Definition at line 47 of file multi_muscl_fluxes_computation.F.

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