OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_nrf_ebcs.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_nrf_ebcs ../engine/source/multifluid/multi_nrf_ebcs.F
25!||--- called by ------------------------------------------------------
26!|| multi_ebcs ../engine/source/multifluid/multi_ebcs.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| multi_submatlaw ../engine/source/multifluid/multi_submatlaw.F
30!||--- uses -----------------------------------------------------
31!|| ebcs_mod ../common_source/modules/boundary_conditions/ebcs_mod.f90
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
34!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
35!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.f
36!|| output_mod ../common_source/modules/output/output_mod.f90
37!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
38!||====================================================================
39 SUBROUTINE multi_nrf_ebcs(ITASK, EBCS_ID, MULTI_FVM, NELEM, ELEM_LIST, FACE_LIST,
40 . FVM_INLET_DATA, IXS, IXQ, IXTG, XGRID, WGRID, IPM, PM, FUNC_VALUE,
41 . EBCS,NPF,TF,FSAVSURF,TIMESTEP,MATPARAM,OUTPUT,PRED)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE multi_fvm_mod
46 USE ebcs_mod
47 USE th_surf_mod , ONLY : th_surf_num_channel
48 USE matparam_def_mod , ONLY : matparam_struct_
49 USE multi_submatlaw_mod , ONLY : multi_submatlaw
50 USE output_mod , ONLY : output_
51 use element_mod , only : nixs,nixq,nixtg
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59! NIXS
60#include "param_c.inc"
61! ISPMD
62#include "task_c.inc"
63! ALE
64#include "com01_c.inc"
65! NUMELS, NUMELQ, NUMELTG
66#include "com04_c.inc"
67! SNPC, STF
68 COMMON /tablesizf/ stf,snpc
69 INTEGER STF,SNPC
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 my_real,INTENT(INOUT) :: fsavsurf(th_surf_num_channel,nsurf)
74 INTEGER, INTENT(IN) :: EBCS_ID
75 TYPE(multi_fvm_struct), INTENT(INOUT) :: MULTI_FVM
76 INTEGER, INTENT(IN) :: ITASK, NELEM, ELEM_LIST(NELEM), FACE_LIST(NELEM)
77 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
78 my_real, INTENT(IN) :: xgrid(3, *), wgrid(3, *)
79 INTEGER, INTENT(IN) :: IPM(NPROPMI, *)
80 my_real, INTENT(IN) :: pm(npropm, *), func_value(*)
81 TYPE(fvm_inlet_data_struct), INTENT(IN) :: FVM_INLET_DATA
82 TYPE(t_ebcs_nrf), INTENT(INOUT) :: EBCS
83 INTEGER, INTENT(IN) :: NPF(SNPC)
84 my_real, INTENT(IN) :: tf(stf), timestep
85 TYPE(matparam_struct_),DIMENSION(NUMMAT), INTENT(IN) :: MATPARAM !material data structure
86 TYPE(output_),INTENT(INOUT) :: OUTPUT
87 LOGICAL,INTENT(IN) :: PRED
88C-----------------------------------------------
89C L o c a l V a r i a b l e s
90C-----------------------------------------------
91 INTEGER :: IELEM, ELEMID
92 INTEGER :: KFACE
93 my_real :: surf, nx, ny, nz
94 my_real :: sstar, pstar, sl, sr, wfac(3), vii(5), normalw
95 my_real :: fii(5), viistar(5), fiistar(5), pp(5)
96 INTEGER :: IMAT, NBMAT
97 INTEGER :: MATLAW(MULTI_FVM%NBMAT), LOCAL_MATID(MULTI_FVM%NBMAT)
98 my_real :: phase_rhoii(multi_fvm%NBMAT), phase_presii(multi_fvm%NBMAT),
99 . phase_eintii(multi_fvm%NBMAT),
100 . phase_alphaii(multi_fvm%NBMAT), phase_rhojj(multi_fvm%NBMAT),
101 . phase_presjj(multi_fvm%NBMAT), phase_eintjj(multi_fvm%NBMAT),
102 . phase_sspjj(multi_fvm%NBMAT), phase_alphajj(multi_fvm%NBMAT)
103 my_real :: dummy(6),dummy2(1), rhoii, pii, eintii, vxii, vyii, vzii, sspii, normal_velii, rhojj, sspjj,
104 . p_jj, normal_veljj, vxjj, vyjj, vzjj, velii2, alphaii, sub_rhoii, sub_rhoeintii, sub_viistar(3),
105 . sub_fiistar(3), alphastar, sub_rhostar, sub_pii, veljj2, sub_estar, eintjj
106 INTEGER :: IELEM_START, IELEM_END, ID_SURF
107 my_real :: TCAR_P, TCAR_VF,ALPHA,BETA,DP0,Vnew,Vold,Pvois,POld,MACH,RHOC2,ROC,PSURF
108 INTEGER :: VARTMP_EOS(1,128) !nel=1
109 INTEGER :: NVARTMP_EOS !upper bound
110 my_real :: aburn(1)
111 my_real :: m_in, m_out, e_in, e_out
112C-----------------------------------------------
113C B o d y
114C-----------------------------------------------
115 nvartmp_eos = 128
116 vartmp_eos(1,:) = 1
117
118 m_in = zero
119 m_out = zero
120 e_in = zero
121 e_out = zero
122
123 tcar_p = ebcs%TCAR_P
124 tcar_vf = ebcs%TCAR_VF
125 id_surf = ebcs%SURF_ID
126 aburn(:) = zero
127
128 alpha = timestep/tcar_p
129 beta = timestep/max(timestep,tcar_vf)
130 IF(tcar_vf>=ep20)beta=zero
131 IF(timestep == zero)THEN
132 alpha = one !zero
133 beta = one
134 ENDIF
135
136 nbmat = multi_fvm%NBMAT
137 ielem_start = 1 + itask * nelem / nthread
138 ielem_end = (1 + itask) * nelem / nthread
139 dummy(1:6) = zero
140 dummy2(:) = one
141 DO ielem = ielem_start, ielem_end
142 elemid = elem_list(ielem)
143 IF (elemid <= numels) THEN
144 DO imat = 1, nbmat
145 local_matid(imat) = ipm(20 + imat, ixs(1, elemid))
146 matlaw(imat) = ipm(2, local_matid(imat))
147 ENDDO
148 ELSE IF (elemid <= numels + numelq) THEN
149 DO imat = 1, nbmat
150 local_matid(imat) = ipm(20 + imat, ixq(1, elemid))
151 matlaw(imat) = ipm(2, local_matid(imat))
152 ENDDO
153 ELSE
154 DO imat = 1, nbmat
155 local_matid(imat) = ipm(20 + imat, ixtg(1, elemid))
156 matlaw(imat) = ipm(2, local_matid(imat))
157 ENDDO
158 ENDIF
159
160C Face of the element
161 kface = face_list(ielem)
162
163 nx = multi_fvm%FACE_DATA%NORMAL(1, kface, elemid)
164 ny = multi_fvm%FACE_DATA%NORMAL(2, kface, elemid)
165 nz = multi_fvm%FACE_DATA%NORMAL(3, kface, elemid)
166 surf = multi_fvm%FACE_DATA%SURF(kface, elemid)
167 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface, elemid)
168C Normal grid velocity
169 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
170
171C Current element
172 rhoii = multi_fvm%RHO(elemid)
173 pii = multi_fvm%PRES(elemid)
174 eintii = multi_fvm%EINT(elemid)
175 vxii = multi_fvm%VEL(1, elemid)
176 vyii = multi_fvm%VEL(2, elemid)
177 vzii = multi_fvm%VEL(3, elemid)
178 velii2 = vxii**2 + vyii**2 + vzii**2
179 IF (nbmat > 1) THEN
180 DO imat = 1, nbmat
181 phase_alphaii(imat) = multi_fvm%PHASE_ALPHA(imat, elemid)
182 phase_eintii(imat) = multi_fvm%PHASE_EINT(imat, elemid)
183 phase_rhoii(imat) = multi_fvm%PHASE_RHO(imat, elemid)
184 phase_presii(imat) = multi_fvm%PHASE_PRES(imat, elemid)
185 ENDDO
186 ELSE
187 phase_alphaii(1) = one
188 phase_eintii(1) = eintii
189 phase_rhoii(1) = rhoii
190 phase_presii(1) = pii
191 ENDIF
192 sspii = multi_fvm%SOUND_SPEED(elemid)
193 normal_velii = vxii * nx + vyii * ny + vzii * nz
194
195C Boundary "GHOST" element
196 rhojj = zero
197 DO imat = 1, nbmat
198 phase_rhojj(imat) = fvm_inlet_data%VAL_RHO(imat)
199 phase_rhojj(imat) = phase_rhoii(imat)
200 phase_alphajj(imat) = fvm_inlet_data%VAL_ALPHA(imat)
201 phase_alphajj(imat) = phase_alphaii(imat)
202 rhojj = rhojj + phase_rhojj(imat) * phase_alphajj(imat)
203 ENDDO
204
205C VE formulation
206 sspjj = zero
207 p_jj = zero
208 eintjj = zero
209 DO imat = 1, nbmat
210 phase_eintjj(imat) = fvm_inlet_data%VAL_PRES(imat)
211 phase_eintjj(imat) = phase_eintii(imat)
212 IF (phase_alphajj(imat) > zero) THEN
213 CALL multi_submatlaw(
214 1 0, matlaw(imat), local_matid(imat), 1,
215 2 phase_eintjj(imat), phase_presjj(imat), phase_rhojj(imat), phase_sspjj(imat),
216 3 dummy2, dummy, pm, ipm,
217 4 npropm, npropmi, dummy, dummy2,
218 5 dummy, multi_fvm%BFRAC(imat,elemid), multi_fvm%TBURN(elemid), dummy,
219 6 dummy(1), dummy, snpc, stf,
220 7 npf, tf, dummy, 1,
221 7 matparam(local_matid(imat)), nvartmp_eos, vartmp_eos, nummat, aburn)
222 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
223 p_jj = p_jj + phase_presjj(imat) * phase_alphajj(imat)
224 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
225 ENDIF
226 ENDDO
227
228 vnew = normal_velii
229 pvois = p_jj
230 mach = abs(normal_velii / sspii)
231 pold = ebcs%pold(ielem)
232 dp0 = ebcs%DP0(ielem)
233 IF(ncycle == 1)pold=pvois+dp0
234 vold = ebcs%vold(ielem)
235 IF(ncycle == 1) vold=normal_velii
236 ebcs%vold(ielem) = vnew
237
238 IF(mach >= one .AND. vnew > zero)THEN
239 !outgoing supersonic velocity: State = Neighboring State
240 pp = pvois
241 ELSE
242 rhoc2 = rhojj*sspii*sspii
243 roc = sqrt(rhojj*rhoc2)
244 p_jj = one/(one+alpha)*(pold+roc*(vnew-vold))+alpha*(pvois+dp0)/(alpha + one)
245 ENDIF
246
247 ebcs%pold(ielem) = p_jj
248
249 IF (sspjj / rhojj > zero) THEN
250 sspjj = sqrt(sspjj / rhojj)
251 ELSE
252 sspjj = multi_fvm%SOUND_SPEED(elemid)
253 ENDIF
254
255 !velocity continuity
256 vxjj = fvm_inlet_data%VAL_VEL(1)
257 vxjj = vxii
258 vyjj = fvm_inlet_data%VAL_VEL(2)
259 vyjj = vyii
260 vzjj = fvm_inlet_data%VAL_VEL(3)
261 vzjj = vzii
262 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
263
264 veljj2 = vxjj**2 + vyjj**2 + vzjj**2
265C HLL wave speed estimates
266 sl = min(normal_velii - sspii, normal_veljj - sspjj)
267 sr = max(normal_velii + sspii, normal_veljj + sspjj)
268
269C Intermediate wave speed
270 sstar = p_jj - pii + rhoii * normal_velii * (sl - normal_velii) - rhojj * normal_veljj * (sr - normal_veljj)
271 sstar = sstar / (rhoii * (sl - normal_velii) - rhojj * (sr - normal_veljj))
272
273 pstar = pii + rhoii * (sstar - normal_velii) * (sl - normal_velii)
274 pp(1) = zero
275 pp(2) = pstar * nx
276 pp(3) = pstar * ny
277 pp(4) = pstar * nz
278 pp(5) = sstar * pstar
279
280 IF (sl > normalw) THEN
281 vii(1) = rhoii
282 vii(2) = rhoii * vxii
283 vii(3) = rhoii * vyii
284 vii(4) = rhoii * vzii
285 vii(5) = eintii + half * rhoii * velii2
286! Normal physical flux current element
287 fii(1) = vii(1) * normal_velii
288 fii(2) = vii(2) * normal_velii + pii * nx
289 fii(3) = vii(3) * normal_velii + pii * ny
290 fii(4) = vii(4) * normal_velii + pii * nz
291 fii(5) = (vii(5) + pii) * normal_velii
292! /th/surf : pressure
293 psurf = pii
294C Take the fluxes of cell II
295C ===
296C Global fluxes
297 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
298 multi_fvm%FLUXES(6, kface, elemid) = normal_velii * surf
299C ===
300C Submaterial fluxes
301 IF (nbmat > 1) THEN
302 DO imat = 1, nbmat
303 alphaii = phase_alphaii(imat)
304 sub_rhoii = phase_rhoii(imat) ! ALPHA_RHO
305 sub_rhoeintii = phase_eintii(imat)
306 sub_viistar(1) = alphaii
307 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
308 sub_viistar(3) = alphaii * sub_rhoeintii
309 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
310 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
311 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
312 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
313 ENDDO
314 ENDIF
315C
316 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
317 vii(1) = rhoii
318 vii(2) = rhoii * vxii
319 vii(3) = rhoii * vyii
320 vii(4) = rhoii * vzii
321 vii(5) = eintii + half * rhoii * velii2
322! Normal physical flux current element
323 fii(1) = vii(1) * normal_velii
324 fii(2) = vii(2) * normal_velii + pii * nx
325 fii(3) = vii(3) * normal_velii + pii * ny
326 fii(4) = vii(4) * normal_velii + pii * nz
327 fii(5) = (vii(5) + pii) * normal_velii
328! /th/surf : pressure
329 psurf = pii
330C Take intermediate state flux (HLLC scheme)
331C ===
332C Global fluxes
333 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
334 viistar(1:5) = viistar(1:5) / (sstar - sl)
335 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
336 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
337 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
338C ===
339C Submaterial fluxes
340 IF (nbmat > 1) THEN
341 DO imat = 1, nbmat
342 matlaw(imat) = ipm(2, local_matid(imat))
343 alphastar = phase_alphaii(imat)
344 sub_rhostar = phase_rhoii(imat) * (normal_velii - sl) / (sstar - sl)
345 IF (alphastar > zero) THEN
346 sub_rhoii = phase_rhoii(imat)
347 sub_rhoeintii = phase_eintii(imat)
348 sub_pii = phase_presii(imat)
349 sub_estar = phase_eintii(imat) / phase_rhoii(imat) -
350 . phase_presii(imat) * (one / sub_rhostar - one / phase_rhoii(imat))
351
352 IF (sub_estar < zero) THEN
353 sub_estar = zero
354 ENDIF
355 ELSE
356 sub_estar = zero
357 ENDIF
358 sub_viistar(1) = alphastar
359 sub_viistar(2) = alphastar * sub_rhostar
360 sub_viistar(3) = alphastar * sub_rhostar * sub_estar
361
362 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
363 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
364 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
365 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
366 ENDDO
367 ENDIF
368
369 ELSE IF (sstar < normalw .AND. normalw <= sr) THEN
370 vii(1) = rhojj
371 vii(2) = rhojj * vxjj
372 vii(3) = rhojj * vyjj
373 vii(4) = rhojj * vzjj
374 vii(5) = eintjj + half * rhojj * veljj2
375! Normal physical flux current element
376 fii(1) = vii(1) * normal_veljj
377 fii(2) = vii(2) * normal_veljj + p_jj * nx
378 fii(3) = vii(3) * normal_veljj + p_jj * ny
379 fii(4) = vii(4) * normal_veljj + p_jj * nz
380 fii(5) = (vii(5) + p_jj) * normal_veljj
381! /th/surf : pressure
382 psurf = p_jj
383C Take intermediate state flux (HLLC scheme)
384C ===
385C Global fluxes
386 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
387 viistar(1:5) = viistar(1:5) / (sstar - sr)
388 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
389 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
390 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
391C ===
392C Submaterial fluxes
393 IF (nbmat > 1) THEN
394 DO imat = 1, nbmat
395 matlaw(imat) = ipm(2, local_matid(imat))
396 alphastar = phase_alphajj(imat)
397 sub_rhostar = phase_rhojj(imat) * (normal_veljj - sr) / (sstar - sr)
398 IF (alphastar > zero) THEN
399 sub_rhoii = phase_rhojj(imat)
400 sub_rhoeintii = phase_eintjj(imat)
401 sub_pii = phase_presjj(imat)
402 sub_estar = phase_eintjj(imat) / phase_rhojj(imat) -
403 . phase_presjj(imat) * (one / sub_rhostar - one / phase_rhojj(imat))
404 IF (sub_estar < zero) THEN
405 sub_estar = zero
406 ENDIF
407 ELSE
408 sub_estar = zero
409 ENDIF
410
411 sub_viistar(1) = phase_alphajj(imat)
412 sub_viistar(2) = phase_alphajj(imat) * sub_rhostar
413 sub_viistar(3) = phase_alphajj(imat) * sub_rhostar * sub_estar
414 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
415 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
416 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
417 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
418 ENDDO
419 ENDIF
420 ELSE IF (sr < normalw) THEN
421 vii(1) = rhojj
422 vii(2) = rhojj * vxjj
423 vii(3) = rhojj * vyjj
424 vii(4) = rhojj * vzjj
425 vii(5) = eintjj + half * rhojj * veljj2
426! Normal physical flux current element
427 fii(1) = vii(1) * normal_veljj
428 fii(2) = vii(2) * normal_veljj + p_jj * nx
429 fii(3) = vii(3) * normal_veljj + p_jj * ny
430 fii(4) = vii(4) * normal_veljj + p_jj * nz
431 fii(5) = (vii(5) + p_jj) * normal_veljj
432! /th/surf : pressure
433 psurf = p_jj
434C Take the fluxes of cell II
435C ===
436C Global fluxes
437 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
438 multi_fvm%FLUXES(6, kface, elemid) = normal_veljj * surf
439C ===
440C Submaterial fluxes
441 IF (nbmat > 1) THEN
442 DO imat = 1, nbmat
443 alphaii = phase_alphajj(imat)
444 sub_rhoii = phase_rhojj(imat) ! ALPHA_RHO
445 sub_rhoeintii = phase_eintjj(imat)
446 sub_viistar(1) = alphaii
447 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
448 sub_viistar(3) = alphaii * sub_rhoeintii
449 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
450 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
451 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
452 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
453 ENDDO
454 ENDIF
455
456 ELSE
457 CALL arret(2)
458C
459 ENDIF
460
461 IF(.NOT. pred)THEN
462 !muscl 'off', or muscl 'on' after 2nd step of prediction/correction
463
464 ! /TH/SURF (massflow,velocity,pressure,acceleration)
465 fsavsurf(2,id_surf) = fsavsurf(2,id_surf) + multi_fvm%FLUXES(1, kface, elemid) ! += rho.S.un
466 fsavsurf(3,id_surf) = fsavsurf(3,id_surf) + multi_fvm%FLUXES(1, kface, elemid) / vii(1) ! += S.un
467 fsavsurf(4,id_surf) = fsavsurf(4,id_surf) + surf*psurf ! += S.P
468 fsavsurf(6,id_surf) = fsavsurf(6,id_surf) + multi_fvm%FLUXES(1, kface, elemid) ! m <- m+dm (cumulative value)
469
470 ! FLUX > 0 if flow is exiting (OUTWARD NORMAL)
471 m_in = m_in - timestep * min(zero, multi_fvm%FLUXES(1, kface, elemid) )
472 m_out = m_out - timestep * max(zero, multi_fvm%FLUXES(1, kface, elemid) )
473 e_in = e_in - timestep * min(zero, multi_fvm%FLUXES(5, kface, elemid) )
474 e_out = e_out - timestep * max(zero, multi_fvm%FLUXES(5, kface, elemid) )
475 ENDIF
476 ENDDO !next IELEM
477
478 IF(.NOT. pred)THEN
479!$OMP CRITICAL
480 output%DATA%INOUT%DM_IN = output%DATA%INOUT%DM_IN + m_in
481 output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + m_out
482 output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + e_in
483 output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + e_out
484!$OMP END CRITICAL
485 ENDIF
486C-----------------------------------------------
487C E n d o f s u b r o u t i n e
488C-----------------------------------------------
489 END SUBROUTINE multi_nrf_ebcs
490
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_nrf_ebcs
subroutine multi_nrf_ebcs(itask, ebcs_id, multi_fvm, nelem, elem_list, face_list, fvm_inlet_data, ixs, ixq, ixtg, xgrid, wgrid, ipm, pm, func_value, ebcs, npf, tf, fsavsurf, timestep, matparam, output, pred)
Definition multi_nrf_ebcs.F:42
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
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:61
th_surf_mod::th_surf_num_channel
integer, parameter th_surf_num_channel
number of /TH/SURF channels : AREA, VELOCITY, MASSFLOW, P A, MASS
Definition th_surf_mod.F:100
arret
subroutine arret(nn)
Definition arret.F:86