OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fv_up_switch.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!|| fv_up_switch ../engine/source/airbag/fv_up_switch.F
25!||--- called by ------------------------------------------------------
26!|| fvbag0 ../engine/source/airbag/fvbag0.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../engine/source/output/message/message.F
29!|| arret ../engine/source/system/arret.F
30!|| fvinjt6 ../engine/source/airbag/fvinjt6.F
31!|| fvinjt8 ../engine/source/airbag/fvinjt8.F
32!|| fvinjt8_1 ../engine/source/airbag/fvinjt8_1.F
33!|| fvtemp ../engine/source/airbag/fvtemp.F
34!|| fvvent0 ../engine/source/airbag/fvvent0.F
35!|| get_u_func ../engine/source/user_interface/ufunc.F
36!|| spmd_fvb_gath ../engine/source/mpi/airbags/spmd_fvb_gath.F
37!|| spmd_fvb_igath ../engine/source/mpi/airbags/spmd_fvb_igath.F
38!||--- uses -----------------------------------------------------
39!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
40!|| finter_mixed_mod ../engine/source/tools/finter_mixed.F90
41!|| fvbag_mod ../engine/share/modules/fvbag_mod.F
42!|| fvmbag_meshcontrol_mod ../common_source/modules/airbag/fvmbag_meshcontrol_mod.F
43!|| message_mod ../engine/share/message_module/message_mod.F
44!|| output_mod ../common_source/modules/output/output_mod.F90
45!|| python_funct_mod ../common_source/modules/python_mod.F90
46!|| sensor_mod ../common_source/modules/sensor_mod.F90
47!||====================================================================
48 SUBROUTINE fv_up_switch(OUTPUT,NN, NEL, ELEM, NJET, NPOLY,LENH,NBA,
49 2 IBAGJET, RBAGJET, NVENT , IBAGHOL, RBAGHOL ,
50 3 P , RHO , TK , U , SSPK ,
51 4 X , V , A , NSENSOR, SENSOR_TAB,
52 5 FSAV , NPC , TF , IVOLU , RVOLU ,
53 6 MPOLH , QPOLH , EPOLH ,
54 7 PPOLH , RPOLH , GPOLH ,
55 8 NPOLH , IFVNOD , RFVNOD , IFVTRI ,
56 9 IFVPOLY, IFVTADR, IFVPOLH ,
57 A IFVPADR, INFO , NNS , NNTR , IFV ,
58 B NPOLHA , DLH , CPAPOLH ,
59 C CPBPOLH, CPCPOLH, RMWPOLH ,
60 D ITAGEL , ELSINI , ICONTACT , IDPOLH ,
61 E ELFMASS, ELFVEL , IBUFA , ELEMA , TAGELA ,
62 F PA , RHOA , TKA , UA , BRNA ,
63 G NNA , NTGA , IBPOLH , DTPOLH , NNT ,
64 H NELT , XXXA , VVVA , NCONA , POROSITY,
65 I ITYP , IGEO , SSPKA ,
66 J GEO , PM , IPM , TPOLH , ELFEHPY ,
67 K CPDPOLH, CPEPOLH, CPFPOLH ,
68 L ELTG , IPARG , MATTG ,
69 M IGROUPTG,IGROUPC, ELBUF_TAB, CFL_COEF,
70 N PDISP_OLD, PDISP, WFEXT, PYTHON )
71C-----------------------------------------------
72C M O D U L E S
73C-----------------------------------------------
74 USE output_mod
75 USE fvbag_mod
76 USE message_mod
77 USE elbufdef_mod
78 USE fvmbag_meshcontrol_mod
79! USE OMP_LIB
80 USE sensor_mod
81 USE output_mod
82 USE python_funct_mod
83 use finter_mixed_mod
84C-----------------------------------------------
85C I M P L I C I T T Y P E S
86C-----------------------------------------------
87#include "implicit_f.inc"
88C-----------------------------------------------
89C C O M M O N B L O C K S
90C-----------------------------------------------
91#include "com01_c.inc"
92#include "com04_c.inc"
93#include "com06_c.inc"
94#include "com08_c.inc"
95#include "scr02_c.inc"
96#include "scr07_c.inc"
97#include "param_c.inc"
98#include "units_c.inc"
99#include "task_c.inc"
100#include "warn_c.inc"
101#include "tabsiz_c.inc"
102#include "mvsiz_p.inc"
103C-----------------------------------------------
104C D U M M Y A R G U M E N T S
105C-----------------------------------------------
106 TYPE(output_), INTENT(INOUT) :: OUTPUT
107 INTEGER, INTENT(IN) :: NPOLY, LENH, NPOLH, NBA,NSENSOR
108 INTEGER, INTENT(IN) :: NN, NEL, NELT, ELEM(3, NELT), NJET, IBAGJET(NIBJET, NJET),
109 . NVENT, NPC(SNPC), IFVNOD(3,NNS), IFVTRI(6,NNTR),
110 . IFVPOLY(NNTR),IFVTADR(NPOLY+1), IFVPOLH(LENH), IFVPADR(NPOLH + 1),
111 . NNS, NNTR, IFV, NPOLHA, ITAGEL(NELT), ICONTACT(SICONTACT),
112 . IDPOLH(NPOLH), IBUFA(NNA), ELEMA(3,NTGA), TAGELA(NTGA), BRNA(8,NBA),
113 . NNA, NTGA, IBPOLH(NPOLH), NNT, NCONA(16,NNA),
114 . ITYP, IGEO(NPROPGI,*), IPM(NPROPMI,NUMMAT), IPARG(NPARG,NGROUP)
115 INTEGER, INTENT(INOUT) :: IVOLU(NIMV), INFO, IBAGHOL(NIBHOL, NVENT)
116 INTEGER, INTENT(IN) :: ELTG(NELT), MATTG(NELT), IGROUPTG(NUMELTG), IGROUPC(NUMELC)
117
118 MY_REAL, INTENT(IN) ::
119 . X(3,NUMNOD), V(3,NUMNOD), A(3,NUMNOD),
120 . TF(STF), XXXA(3,NNA), VVVA(3,NNA),
121 . GEO(NPROPG,NUMGEO), PM(NPROPM,NUMMAT), CFL_COEF, PDISP_OLD
122 MY_REAL, INTENT(INOUT) :: POROSITY(NELT - NEL)
123 MY_REAL, INTENT(INOUT) :: RBAGJET(NRBJET,NJET)
124 MY_REAL, INTENT(INOUT) ::
125 . P(NNT), RHO(NNT), TK(NNT), U(3,NNT), SSPK(NNT), MPOLH(NPOLH), QPOLH(3,NPOLH),
126 . EPOLH(NPOLH), PPOLH(NPOLH), RPOLH(NPOLH), GPOLH(NPOLH), RFVNOD(2,NNS), DLH,
127 . CPAPOLH(NPOLH), CPBPOLH(NPOLH), CPCPOLH(NPOLH), RMWPOLH(NPOLH), ELSINI(NELT),
128 . elfmass(nelt), elfvel(nelt), pa(nna), rhoa(nna), tka(nna),sspka(nna), ua(3,nna),
129 . dtpolh(npolh), fsav(nthvki), tpolh(npolh), cpdpolh(npolh), cpepolh(npolh),
130 . cpfpolh(npolh), rbaghol(nrbhol,nvent), pdisp, elfehpy(nelt), rvolu(nrvolu)
131 TYPE(ELBUF_STRUCT_), DIMENSION(NGROUP), INTENT(IN) :: ELBUF_TAB
132 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
133 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
134 TYPE(PYTHON_),INTENT(IN) :: PYTHON
135C-----------------------------------------------
136C L O C A L V A R I A B L E S
137C-----------------------------------------------
138 INTEGER I, II, IINJ, IEL, N1, N2, N3, JEL,
139 . NN1, NN2, NN3, J, JJ, K, KK, ISENS, IVEL,
140 . itagp(npolh), nv, ilvout, ipri,
141 . ivent, idef, iporp,
142 . itvent, phdt, i1, i2, icont(nnt),ivdp,ittf, idpdef,
143 . idtpdef,itypl,
144 . prepare_anim,iventyp
145 INTEGER TITREVENT(20)
146
147 my_real
148 . nodarea(nnt), x1, y1, z1, x2, y2, z2, x3, y3, z3,
149 . x12, y12, z12, x13, y13, z13, nrx, nry, nrz, area2,
150 . elarea(nelt), norm(3,nelt), ksi, eta,
151 . pnorm(3,nntr), pvolu(npolh), area, nx, ny,
152 . nz, gamai, cpai, cpbi, cpci, rmwi, ti,
153 . cpdi, cpei, cpfi,
154 . rhoi, datainj(6,njet), scalt,
155 . fvel, tstart, tsg, injvel,
156 . cpa, cpb, cpc, efac,
157 . pu(3,npolh), gama, fel(nelt), dm(npolh),
158 . dq(3,npolh), de(npolh), dmi, rho1, re1, ux1,
159 . uy1, uz1, rho2, re2, ux2, uy2, uz2, vfx, vfy, vfz, vrel(nelt), ss(nelt), ss_,
160 . alpha, rhom, rem, ruxm, ruym, ruzm, p1, p2,
161 . uel(3,nelt), area1, area3, ff, temp, pext, volg, volph,
162 . areap, areael, qa, qb, dtx, al, dd, ad, ssp,
163 . qx, vv, dmini(npolh), dmcpa(npolh), dmcpb(npolh),
164 . dmcpc(npolh), dmrmw(npolh), msini,
165 . rmw, cp, cv, rhoel(nelt), tkel(nelt), sspel(nelt), gama1, gama2,
166 . qvisc(npolh), amtot, pmean, pmean2, elsout(nelt),tmp(nel),
167 . area_vent(nvent), pm_vent(nvent), scalp, scals,
168 . aout, deri, pdef, pvoltmp,
169 . dtpdefi, dtpdefc, tvent, poro(nelt), vola(nna),
170 . x23, y23, z23, x31, y31, z31, l12, l23, l31, h1,
171 . fac1, fac2, fac3, area_old, pp, dls(npolh), vmax, uu,
172 . pcrit, vx1, vy1, vz1, vx2, vy2, vz2,
173 . vx3, vy3, vz3, vvx, vvy, vvz, vel(nelt), fac,
174 . aoutot, svent(nvent), sventot, xxx(3,nnt),
175 . vvv(3,nnt), fvdp, rbid, ttf,
176 . dmcpd(npolh), dmcpe(npolh), dmcpf(npolh),
177 . cpd, cpe, cpf, temp0, pel(nelt), pstag, volno,tstope,
178 . enint, massflow, wfext0, eldmass(nel), eldehpy(nel), dmout, dhout,
179 . xx1, xx2, xx3, yy1, yy2, yy3, zz1, zz2, zz3, nnx, nny, nnz,
180 . pres_av, cp_av, cv_av, rgas_av, rho_av
181
182
183 LOGICAL :: EXIT_NEG_VOL,INJECTION_STARTED
184 LOGICAL :: UP_SWITCH
185C
186 MY_REAL
187 . GET_U_FUNC
188 EXTERNAL GET_U_FUNC
189 CHARACTER*20 VENTTITLE
190 INTEGER :: NODE_ID(5), INODE
191 MY_REAL :: TAB_PVOL(NELT), DOT_PROD
192 MY_REAL :: MOMENTUM_FLUX_X(NELT), MOMENTUM_FLUX_Y(NELT), MOMENTUM_FLUX_Z(NELT),
193 . MASS_FLUX(NELT), ENERGY_FLUX(NELT),
194 . CPA_FLUX(NELT), CPB_FLUX(NELT), CPC_FLUX(NELT), CPD_FLUX(NELT),
195 . CPE_FLUX(NELT), CPF_FLUX(NELT), RGAS_FLUX(NELT)
196 INTEGER(8) :: VEC_PTR_PLUS, VEC_PTR_MINUS
197 INTEGER :: COUNT, IDT_FVMBAG, PARAM_L_TYPE
198 MY_REAL :: CPAM, CPBM, CPCM, RGASM, CPDM, CPEM, CPFM
199 MY_REAL :: DTOLD, LAMBDA
200
201C-----------------------------------------------
202C C O M M E N T S
203C-----------------------------------------------
204C NEL : number of triangles (external surface only)
205C NELT : number of triangles (both internal & external surfaces)
206C NNS : number of nodes of TETRA mesh
207C NNA : number of 'additional' nodes
208C NNT : total number of nodes
209C NNTR : number of triangles (all polyhedra)
210C NPOLH : number of polyhedron
211C NPOLY : number of polygons
212
213C-----------------------------------------------
214C S O U R C E L I N E S
215C-----------------------------------------------
216 volph = zero
217 temp = zero
218
219 IF(nbgauge > 0 ) THEN
220 ALLOCATE(fvspmd(ifv)%GGG(3,nnt))
221 ALLOCATE(fvspmd(ifv)%GGA(3,nna))
222 ENDIF
223 ALLOCATE(fvspmd(ifv)%AAA(3,nnt))
224 fvspmd(ifv)%AAA(1:3,1:nnt) = zero
225
226! =================== !
227! PARITH/ON TREATMENT !
228! =================== !
229 IF (nspmd == 1) THEN
230 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
231 i1=fvspmd(ifv)%IBUF_L(1,i)
232 i2=fvspmd(ifv)%IBUF_L(2,i)
233 xxx(1,i1)=x(1,i2)
234 xxx(2,i1)=x(2,i2)
235 xxx(3,i1)=x(3,i2)
236 ENDDO
237
238 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
239 i1=fvspmd(ifv)%IBUF_L(1,i)
240 i2=fvspmd(ifv)%IBUF_L(2,i)
241 vvv(1,i1)=v(1,i2)
242 vvv(2,i1)=v(2,i2)
243 vvv(3,i1)=v(3,i2)
244 ENDDO
245
246 IF (intbag/=0.AND.nvent>0.AND.ivolu(39)/=0) THEN
247 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
248 i1=fvspmd(ifv)%IBUF_L(1,i)
249 i2=fvspmd(ifv)%IBUF_L(2,i)
250 icont(i1)=icontact(i2)
251 ENDDO
252 ENDIF
253
254 ELSE
255 CALL spmd_fvb_gath(ifv, x, xxx, rbid, rbid, 4)
256
257 CALL spmd_fvb_gath(ifv, v, vvv, rbid, rbid, 4)
258
259 IF (intbag/=0.AND.nvent>0.AND.ivolu(39)/=0)
260 . CALL spmd_fvb_igath(ifv, icontact, icont)
261
262
263 IF (ispmd/=fvspmd(ifv)%PMAIN-1) RETURN
264 END IF
265
266! ================================= !
267! REBUILD CONNECTIVITY AFTER SWITCH
268! ================================= !
269
270! INITIALIZATION PHASE, DONE ONLY ONCE
271! ------------------------------------
272
273! NB, FOR AN ELEMENT IEL, THE 'MINUS SIDE' IS BEHIND THE NORMAL, WHILE 'PLUS SIDE' IS IN THE DIRECTION OF
274! THE NORMAL
275
276 IF (ivolu(74) == -2) THEN
277 ALLOCATE(fvdata(ifv)%TRI_TO_ELEM(2, nelt))
278 fvdata(ifv)%TRI_TO_ELEM(:, :) = 0
279 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(npolh))
280 CALL intvector_create(vec_ptr_plus)
281 CALL intvector_create(vec_ptr_minus)
282 DO i = 1, npolh
283! CLEAR INT VECTORS
284 CALL intvector_clear(vec_ptr_plus)
285 CALL intvector_clear(vec_ptr_minus)
286 count = 0
287! LOOP ON POLYHEDRAL FACES
288 DO j = ifvpadr(i), ifvpadr(i+1)-1
289 jj = ifvpolh(j)
290! LOOP ON TRIANGLES OF THE POLYHEDRAL FACE
291 DO k = ifvtadr(jj), ifvtadr(jj+1)-1
292 kk = ifvpoly(k)
293 iel = ifvtri(4, kk)
294 n1 = ifvtri(1, kk)
295 IF (ifvnod(1, n1) /= 2) THEN
296 CALL abort()
297 ENDIF
298 n1 = ifvnod(3, n1)
299 x1 = xxxa(1, n1)
300 y1 = xxxa(2, n1)
301 z1 = xxxa(3, n1)
302 n2 = ifvtri(2, kk)
303 IF (ifvnod(1, n2) /= 2) THEN
304 CALL abort()
305 ENDIF
306 n2 = ifvnod(3, n2)
307 x2 = xxxa(1, n2)
308 y2 = xxxa(2, n2)
309 z2 = xxxa(3, n2)
310 n3 = ifvtri(3, kk)
311 IF (ifvnod(1, n3) /= 2) THEN
312 CALL abort()
313 ENDIF
314 n3 = ifvnod(3, n3)
315 x3 = xxxa(1, n3)
316 y3 = xxxa(2, n3)
317 z3 = xxxa(3, n3)
318 nx = (y2 - y1) * (z3 - z1) - (z2 - z1) * (y3 - y1)
319 ny = (z2 - z1) * (x3 - x1) - (x2 - x1) * (z3 - z1)
320 nz = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)
321 IF (iel /= 0) THEN
322 count = count + 1
323 nn1 = elem(1, abs(iel))
324 nn2 = elem(2, abs(iel))
325 nn3 = elem(3, abs(iel))
326 xx1 = xxx(1, nn1)
327 yy1 = xxx(2, nn1)
328 zz1 = xxx(3, nn1)
329 xx2 = xxx(1, nn2)
330 yy2 = xxx(2, nn2)
331 zz2 = xxx(3, nn2)
332 xx3 = xxx(1, nn3)
333 yy3 = xxx(2, nn3)
334 zz3 = xxx(3, nn3)
335 nnx = (yy2 - yy1) * (zz3 - zz1) - (zz2 - zz1) * (yy3 - yy1)
336 nny = (zz2 - zz1) * (xx3 - xx1) - (xx2 - xx1) * (zz3 - zz1)
337 nnz = (xx2 - xx1) * (yy3 - yy1) - (yy2 - yy1) * (xx3 - xx1)
338 dot_prod = nx * nnx + ny * nny + nz * nnz
339! EXTERNAL (IEL > 0) OR INTERNAL (IEL < 0) SURFACE ELEMENT
340 IF (iel < 0) THEN
341! INTERNAL, THERE IS A POTENTIAL NEIGHBOR
342 iel = -iel
343 i1 = ifvtri(5, kk)
344 i2 = ifvtri(6, kk)
345 IF (i1 == i2) THEN
346 CALL abort()
347 ENDIF
348 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i1
349 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i2
350 IF (i1 == i) THEN
351 IF (dot_prod >= zero) THEN
352 CALL intvector_push_back(vec_ptr_minus, iel)
353 ELSE
354 CALL intvector_push_back(vec_ptr_plus, iel)
355 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i2
356 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i1
357 ENDIF
358 ELSE
359 IF (dot_prod >= zero) THEN
360 CALL intvector_push_back(vec_ptr_plus, iel)
361 ELSE
362 CALL intvector_push_back(vec_ptr_minus, iel)
363 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i2
364 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i1
365 ENDIF
366 ENDIF
367 ELSE
368! EXTERNAL SURFACE : NO NEIGHBOR
369 IF (ifvtri(5, kk) /= ifvtri(6, kk)) THEN
370 CALL abort()
371 ENDIF
372 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i
373 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i
374 IF (dot_prod >= zero) THEN
375 CALL intvector_push_back(vec_ptr_minus, iel)
376 ELSE
377 CALL intvector_push_back(vec_ptr_plus, iel)
378 ENDIF
379 ENDIF
380 ELSE
381 CALL abort()
382 ENDIF
383 ENDDO
384 ENDDO
385 CALL intvector_get_size(vec_ptr_minus, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE)
386 CALL intvector_get_size(vec_ptr_plus, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE)
387 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE))
388 CALL intvector_copy_to(vec_ptr_minus, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(1))
389 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE))
390 CALL intvector_copy_to(vec_ptr_plus, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(1))
391 ENDDO
392 ivolu(74) = -20
393 CALL intvector_delete(vec_ptr_plus)
394 CALL intvector_delete(vec_ptr_minus)
395 ENDIF
396
397
398!$OMP PARALLEL PRIVATE(I,J)
399!$OMP+ PRIVATE(I1,I2,KSI,ETA,N1,N2,N3,IEL,X1,Y1,Z1,FAC,NN1,NN2,NN3)
400!$OMP+ PRIVATE(X2,Y2,Z2,X3,Y3,Z3,VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
401!$OMP+ PRIVATE(X12,Y12,Z12,X13,Y13,Z13,AREA2,NRX,NRY,NRZ,JJ,AREA,KK,K)
402!$OMP+ PRIVATE(NX,NY,NZ,PVOLTMP)
403
404! ================================================= !
405! SURFACE AND NORMAL COMPUTATION FOR FABRIC ELEMENT !
406! ================================================= !
407
408!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
409 DO i=1,nnt
410 nodarea(i)=zero
411 ENDDO
412!$OMP END DO
413!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
414 DO iel=1,nelt
415 n1=elem(1,iel)
416 n2=elem(2,iel)
417 n3=elem(3,iel)
418 x1=xxx(1,n1)
419 x2=xxx(1,n2)
420 x3=xxx(1,n3)
421 y1=xxx(2,n1)
422 y2=xxx(2,n2)
423 y3=xxx(2,n3)
424 z1=xxx(3,n1)
425 z2=xxx(3,n2)
426 z3=xxx(3,n3)
427 x12=x2-x1
428 y12=y2-y1
429 z12=z2-z1
430 x13=x3-x1
431 y13=y3-y1
432 z13=z3-z1
433 nrx=y12*z13-z12*y13
434 nry=z12*x13-x12*z13
435 nrz=x12*y13-y12*x13
436 area2=sqrt(nrx**2+nry**2+nrz**2)
437 elarea(iel)=area2
438 norm(1,iel)=nrx/area2
439 norm(2,iel)=nry/area2
440 norm(3,iel)=nrz/area2
441 IF(iel<=nel) THEN
442 tmp(iel) = one_over_6 * (x1*nrx+y1*nry+z1*nrz)
443 ENDIF
444 pel(iel)=zero
445 ENDDO
446!$OMP END DO
447
448!$OMP SINGLE
449 volg=zero
450 DO iel=1,nelt
451 n1=elem(1,iel)
452 n2=elem(2,iel)
453 n3=elem(3,iel)
454 area2 = elarea(iel)
455 elarea(iel) = half * area2
456 nodarea(n1)=nodarea(n1)+one_over_6*area2
457 nodarea(n2)=nodarea(n2)+one_over_6*area2
458 nodarea(n3)=nodarea(n3)+one_over_6*area2
459 IF(iel<=nel) THEN
460 volg=volg+tmp(iel)
461 ENDIF
462 ENDDO
463C
464 IF (dt1==zero.OR.ivolu(39)==0) THEN
465 DO iel=1,nelt
466 elfmass(iel)=zero
467 elfehpy(iel)=zero
468 ENDDO
469 fsav(1:nthvki)=zero
470 ENDIF
471!$OMP END SINGLE
472
473! ============================ !
474! VOLUME OF THE FINITE VOLUMES !
475! ============================ !
476
477!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
478 DO iel = 1, nelt
479 nx = norm(1, iel)
480 ny = norm(2, iel)
481 nz = norm(3, iel)
482 n1 = elem(1, iel)
483 x1 = xxx(1, n1)
484 x2 = xxx(2, n1)
485 x3 = xxx(3, n1)
486 area = elarea(iel)
487 pvoltmp = third * area * (nx * x1 + ny * x2 + nz * x3)
488 tab_pvol(iel) = pvoltmp
489 ENDDO
490!$OMP END DO
491
492!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
493 DO i = 1, npolh
494 pvolu(i) = zero
495 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
496 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
497! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
498 pvolu(i) = pvolu(i) + tab_pvol(iel)
499 ENDDO
500 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
501 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
502! finite volume i is on the right side of element iel, hence a - sign
503 pvolu(i) = pvolu(i) - tab_pvol(iel)
504 ENDDO
505 ENDDO
506!$OMP END DO
507!$OMP END PARALLEL
508
509
510C---------------------------------------------------
511C VERIFICATIONS
512C---------------------------------------------------
513 areael=zero
514 DO iel=1,nel
515 areael=areael+elarea(iel)
516 ENDDO
517 areap=areael
518 rvolu(18)=areael
519 ilvout=ivolu(44)
520 volph=zero
521 exit_neg_vol = .false.
522 DO i=1,npolh
523 volph=volph+pvolu(i)
524 IF(pvolu(i) <= 0) exit_neg_vol = .true.
525 ENDDO
526
527 ipri=mod(ncycle,iabs(ncpri))
528 IF(exit_neg_vol) THEN
529 DO i=1,npolh
530 IF (pvolu(i)<=zero) THEN
531 info=1
532 ierr=ierr+1
533 CALL ancmsg(msgid=185,anmode=aninfo,
534 . i1=idpolh(i),r1=pvolu(i),i3=i,i4=npolh)
535 CALL arret(2)
536 ENDIF
537 ENDDO
538 ENDIF
539
540
541C---------------------------------------------------
542C COMPUTATION OF QUANTITIES IN POLYEDRA
543C---------------------------------------------------
544 IF (dt1==zero) THEN
545 gamai=rvolu(1)
546 gpolh(1:npolh)=gamai
547 cpai =rvolu(7)
548 cpapolh(1:npolh)=cpai
549 cpbi =rvolu(8)
550 cpbpolh(1:npolh)=cpbi
551 cpci =rvolu(9)
552 cpcpolh(1:npolh)=cpci
553 rmwi =rvolu(10)
554 rmwpolh(1:npolh)=rmwi
555 ti =rvolu(13)
556 tpolh(1:npolh)=ti
557 cpdi =rvolu(56)
558 cpdpolh(1:npolh)=cpdi
559 cpei =rvolu(57)
560 cpepolh(1:npolh)=cpei
561 cpfi =rvolu(58)
562 cpfpolh(1:npolh)=cpfi
563 rhoi =rvolu(62)
564 rpolh(1:npolh)=rhoi
565 ENDIF
566 IF (dt1==zero.OR.ivolu(39)==0) THEN
567 rhoi =rvolu(62)
568 efac =rvolu(66)
569 mpolh(1:npolh)=rhoi*pvolu(1:npolh)
570 epolh(1:npolh)=mpolh(1:npolh)*efac
571 pu(1,1:npolh) = rvolu(67)
572 pu(2,1:npolh) = rvolu(68)
573 pu(3,1:npolh) = rvolu(69)
574 qpolh(1,1:npolh) = mpolh(1:npolh)*pu(1,1:npolh)
575 qpolh(2,1:npolh) = mpolh(1:npolh)*pu(2,1:npolh)
576 qpolh(3,1:npolh) = mpolh(1:npolh)*pu(3,1:npolh)
577 ENDIF
578
579 pext =rvolu(3)
580 IF(ivolu(39) /= 0) THEN
581
582C---------------------------------------------------
583C INJECTORS
584C---------------------------------------------------
585 datainj(1:6,1:njet)=zero
586 DO iel=1,nelt
587 IF (itagel(iel)>0) THEN
588 iinj=itagel(iel)
589 datainj(1,iinj)=datainj(1,iinj)+elarea(iel)
590 ENDIF
591 ENDDO
592C
593 scalt =rvolu(26)
594 IF(ityp==6) THEN
595 CALL fvinjt6(njet, ibagjet, rbagjet, npc, tf,nsensor,
596 . sensor_tab,scalt, datainj, python )
597 ELSEIF(ityp==8) THEN
598 CALL fvinjt8(njet, ibagjet, rbagjet, npc, tf,nsensor,
599 . sensor_tab,scalt, igeo, geo, pm, ivolu, datainj,python)
600 ENDIF
601C
602 DO iinj=1,njet
603 isens=ibagjet(4,iinj)
604 fvel=rbagjet(15,iinj)
605 ivel=ibagjet(11,iinj)
606 IF(isens==0)THEN
607 tstart=zero
608 ELSE
609 tstart=sensor_tab(isens)%TSTART
610 ENDIF
611C ITTF : 1 SHIFT EVENT TIME (CAN BE REMOVED IN A LATER VERSION)
612C ITTF : 2 SHIFT EVENT TIME + EVENT FUNCTION
613C ITTF : 1 OR 2 +10 1 FLAG HAS BEEN ACTIVATED AND RVOLU(60) IS FILLED
614 ittf=ivolu(17)
615 IF (tt>=tstart.AND.(ittf==1.OR.ittf==2.OR.ittf==3))THEN
616 ittf=ittf+10
617 rvolu(60)=tstart
618 ivolu(17)=ittf
619 END IF
620 IF (tt>=tstart.AND.dt1>zero)THEN
621 tsg=(tt-tstart)*scalt
622 IF (ivel>0) THEN
623 injvel=fvel*finter_mixed(python,nfunct,ivel,tsg,npc,tf)
624 ELSE
625 injvel = fvel
626 ENDIF
627 ELSE
628 injvel=zero
629 ENDIF
630 datainj(3,iinj)=injvel
631 ENDDO
632C---------------------------------------------------
633C VENT HOLES
634C---------------------------------------------------
635 IF (intbag==0) THEN
636 poro(1:nelt)=zero
637 ELSE
638!$OMP PARALLEL PRIVATE(IEL,N1,N2,N3)
639!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
640 DO iel=1,nelt
641 n1=elem(1,iel)
642 n2=elem(2,iel)
643 n3=elem(3,iel)
644 poro(iel)=zero
645 IF (icont(n1)/=0) poro(iel)=poro(iel)+third
646 IF (icont(n2)/=0) poro(iel)=poro(iel)+third
647 IF (icont(n3)/=0) poro(iel)=poro(iel)+third
648 ENDDO
649!$OMP END DO
650!$OMP END PARALLEL
651 ENDIF
652C
653 pext =rvolu(3)
654 scalt=rvolu(26)
655 scalp=rvolu(27)
656 scals=rvolu(28)
657 ttf =rvolu(60)
658C-----------------------------------------------------------
659C EFFECTIVE LEAKAGE SURFACE OF TRIANGLES
660C CORRECTION OF SCALE FACTORS ON THE VENT HOLES SURFACE
661C-----------------------------------------------------------
662 CALL fvvent0(
663 1 elsout ,aoutot ,nvent ,nelt ,ittf ,
664 2 elarea ,elsini ,elem ,itagel ,svent ,
665 3 ibaghol ,rvolu ,rbaghol ,poro ,p ,
666 4 eltg ,iparg ,mattg ,nel ,porosity,
667 5 ipm ,pm ,elbuf_tab,igroupc ,igrouptg)
668C---------------------------------------------------
669C FINITE VOLUME BALANCES
670C---------------------------------------------------
671!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,AREA,N1,N2,N3)
672!$OMP+ PRIVATE(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X12,Y12,Z12,X23,Y23,Z23,Z31,Y31,X31)
673!$OMP+ PRIVATE(L12,L23,L31,H1,NX,NY,NZ,IINJ,DMI,IVENT,IDEF,ITVENT,AOUT)
674!$OMP+ PRIVATE(P1,RHO1,RE1,GAMA1,UX1,UY1,UZ1,UU,RHO2,VV,ETA,KSI,PSTAG)
675!$OMP+ PRIVATE(PCRIT,P2,VMAX,IVDP,FVDP,GAMAI,RHOI,DMOUT,DHOUT,CPAI,CPBI,CPCI)
676!$OMP+ PRIVATE(RMWI,CPDI,CPEI,CPFI,VX1,VY1,VZ1,VX2,VY2,VZ2,VX3,VY3,VZ3,VVX,VVY,VVZ,I1,I2,GAMA2)
677!$OMP+ PRIVATE(VFX,VFY,VFZ,SS_,ALPHA,RHOM,RUXM,RUYM,RUZM,REM,MASSFLOW,CPAM,CPBM,CPCM,RGASM,CPDM,CPEM,CPFM)
678!$OMP SECTIONS
679!$OMP SECTION
680 fel(1:nelt) = zero
681 elfvel(1:nelt)=zero
682!$omp section
683 eldmass(1:nel)=zero
684 eldehpy(1:nel)=zero
685!$OMP SECTION
686 itagp(1:npolh)=0
687 dm(1:npolh)=zero
688!$OMP SECTION
689 dq(1,1:npolh)=zero
690 dq(2,1:npolh)=zero
691!$OMP SECTION
692 dq(3,1:npolh)=zero
693 de(1:npolh)=zero
694!$OMP SECTION
695 dmini(1:npolh)=zero
696 dmcpa(1:npolh)=zero
697!$OMP SECTION
698 dmcpb(1:npolh)=zero
699 dmcpc(1:npolh)=zero
700!$OMP SECTION
701 dmrmw(1:npolh)=zero
702 dmcpd(1:npolh)=zero
703!$OMP SECTION
704 dmcpe(1:npolh)=zero
705 dmcpf(1:npolh)=zero
706!$OMP SECTION
707C----------------------
708C TRANSPORTS
709C----------------------
710 pu(1,1:npolh)=qpolh(1,1:npolh)/mpolh(1:npolh)
711 pu(2,1:npolh)=qpolh(2,1:npolh)/mpolh(1:npolh)
712 pu(3,1:npolh)=qpolh(3,1:npolh)/mpolh(1:npolh)
713!$OMP END SECTIONS
714
715
716! ================== !
717! FLUXES COMPUTATION !
718! ================== !
719!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
720 DO iel = 1, nelt
721 mass_flux(iel) = zero
722 momentum_flux_x(iel) = zero
723 momentum_flux_y(iel) = zero
724 momentum_flux_z(iel) = zero
725 energy_flux(iel) = zero
726 cpa_flux(iel) = zero
727 cpb_flux(iel) = zero
728 cpc_flux(iel) = zero
729 rgas_flux(iel) = zero
730 IF (ityp == 8) THEN
731 cpd_flux(iel) = zero
732 cpe_flux(iel) = zero
733 cpf_flux(iel) = zero
734 ENDIF
735 elfmass(iel) = zero
736 elfehpy(iel) = zero
737 elfvel(iel) = zero
738 vel(iel) = zero
739! SURFACE
740 area = elarea(iel)
741 ss(iel) = zero
742 vrel(iel) = zero
743! NORMAL VECTOR
744 nx = norm(1, iel)
745 ny = norm(2, iel)
746 nz = norm(3, iel)
747 IF (itagel(iel) > 0) THEN
748! INJECTOR
749 iinj = itagel(iel)
750 dmi = datainj(2, iinj) * area / datainj(1, iinj)
751 mass_flux(iel) = dmi
752 momentum_flux_x(iel) = -dmi * nx * datainj(3, iinj)
753 momentum_flux_y(iel) = -dmi * ny * datainj(3, iinj)
754 momentum_flux_z(iel) = -dmi * nz * datainj(3, iinj)
755 energy_flux(iel) = dmi * datainj(4, iinj)
756 cpa_flux(iel) = dmi * rbagjet(2, iinj)
757 cpb_flux(iel) = dmi * rbagjet(3, iinj)
758 cpc_flux(iel) = dmi * rbagjet(4, iinj)
759 rgas_flux(iel) = dmi * rbagjet(1, iinj)
760 IF (ityp == 8) THEN
761 cpd_flux(iel) = dmi * rbagjet(16, iinj)
762 cpe_flux(iel) = dmi * rbagjet(17, iinj)
763 cpf_flux(iel) = dmi * rbagjet(18, iinj)
764 ENDIF
765 elfmass(iel) = elfmass(iel) + dmi * dt1
766 elfehpy(iel) = elfehpy(iel) + dmi * datainj(4, iinj) * dt1
767 elfvel(iel) = -datainj(3, iinj)
768 ELSE IF (itagel(iel) < 0) THEN
769! VENT HOLE
770 ivent = -itagel(iel)
771 idef = ibaghol(1, ivent)
772 itvent = ibaghol(10, ivent)
773 aout = elsout(iel)
774 IF (fvdata(ifv)%TRI_TO_ELEM(1, iel) /= fvdata(ifv)%TRI_TO_ELEM(2, iel)) THEN
775 CALL abort()
776 ENDIF
777 i = fvdata(ifv)%TRI_TO_ELEM(1, iel)
778 p1 = ppolh(i)
779 rho1 = rpolh(i)
780 re1 = epolh(i) / pvolu(i)
781 gama1 = gpolh(i)
782 ux1 = pu(1,i)
783 uy1 = pu(2,i)
784 uz1 = pu(3,i)
785 uu = zero
786! LOCAL VELOCITY
787 IF (itvent == 1) THEN
788 uu = nx * ux1 + ny * uy1 + nz * uz1
789 rho2 = rho1
790 IF (p1 < pext) uu = zero
791! BERNOULLI
792 ELSEIF ((itvent == 2 .AND. p1 > zero) .OR. (itvent == 4 .AND. p1 >= pext)) THEN
793 vv = nx * ux1 + ny * uy1 + nz * uz1
794 vv = max(vv, zero)
795 eta = (gama1 - one) / gama1
796 ksi = one / eta
797 pstag = p1 * (one + half * eta * rho1 * vv * vv / p1)**ksi
798 pcrit = pstag * (two / (gama1+one))**ksi
799 p2 = max(pext, pcrit)
800 rho2 = rho1 * (p2 / p1)**(one / gama1)
801 uu = two * p1 * (one - (p2 / p1)**eta) / (rho1 * eta)
802 uu = max(uu, zero)
803 uu = sqrt(uu)
804 vmax = half * ((p1 - pext) * pvolu(i) / (gama1 - one))
805 . / max(em20, rho2 * aout * dt1 * gama1 * re1 / rho1)
806 vmax = max(vmax, zero)
807 uu = min(uu, vmax)
808! CHEMKIN
809 ELSE IF (itvent == 3) THEN
810 p2 = p1 - pext
811 ivdp = ibaghol(9, ivent)
812 fvdp = rbaghol(13, ivent)
813 uu = fvdp * get_u_func(ivdp, p2 * scalp, deri)
814 IF (p2 < zero) THEN
815 rho2 = rvolu(62)
816 ELSE
817 rho2 = rho1
818 ENDIF
819! BERNOULLI WITH FLOW-IN
820 ELSEIF (itvent == 4 .AND. p1 < pext) THEN
821 gamai = rvolu(1)
822 rhoi = rvolu(62)
823 eta = (gamai - one) / gamai
824 pcrit = pext * (two / (gamai + one))**(one / eta)
825 p2 = max(p1, pcrit)
826 rho2 = rhoi * (p2 / pext)**(one / gamai)
827 uu = two * pext * (one - (p2 / pext)**eta)/(rhoi * eta)
828 uu = max(uu, zero)
829 uu = sqrt(uu)
830 vmax = half*((pext - p1) * pvolu(i) / (gama1 - one))
831 . / max(em20, rho2 * aout * dt1 * rvolu(63))
832 vmax = max(vmax, zero)
833 uu = min(uu, vmax)
834 uu = -uu
835! GRAEFE
836 ELSE IF (itvent == 5) THEN
837 p2 = p1 - pext
838 rho2 = rho1
839 uu = two * p2 / rho2
840 uu = max(uu, zero)
841 uu = sqrt(uu)
842 ENDIF
843
844 IF (uu > zero .AND. idef == 1) THEN
845 dmout = uu * rho2 * aout
846 mass_flux(iel) = -dmout
847 momentum_flux_x(iel) = -dmout * ux1
848 momentum_flux_y(iel) = -dmout * uy1
849 momentum_flux_z(iel) = -dmout * uz1
850 dhout = dmout * gama1 * re1 / rho1
851 energy_flux(iel) = -dhout
852 cpai = cpapolh(i)
853 cpbi = cpbpolh(i)
854 cpci = cpcpolh(i)
855 rmwi = rmwpolh(i)
856 cpa_flux(iel) = -dmout * cpai
857 cpb_flux(iel) = -dmout * cpbi
858 cpc_flux(iel) = -dmout * cpci
859 rgas_flux(iel) = -dmout * rmwi
860 IF(ityp == 8) THEN
861 cpdi = cpdpolh(i)
862 cpei = cpepolh(i)
863 cpfi = cpfpolh(i)
864 cpd_flux(iel) = -dmout * cpdi
865 cpe_flux(iel) = -dmout * cpei
866 cpf_flux(iel) = -dmout * cpfi
867 ENDIF
868 !/TODO
869 !DMINI(I)=DMINI(I)-DMOUT
870 elfmass(iel)=elfmass(iel)-dmout*dt1
871 elfehpy(iel)=elfehpy(iel)-dhout*dt1
872 elfvel(iel)=elfvel(iel)+uu*area/elarea(iel)
873 eldmass(iel)=-dmout
874 eldehpy(iel)=-dhout
875 ENDIF
876C
877 IF (uu < zero .AND. idef == 1) THEN
878 dmi = -uu * rho2 * aout
879 mass_flux(iel) = dmi
880 momentum_flux_x(iel) = dmi * ux1
881 momentum_flux_y(iel) = dmi * uy1
882 momentum_flux_z(iel) = dmi * uz1
883 cpai = rvolu(7)
884 cpbi = rvolu(8)
885 cpci = rvolu(9)
886 rmwi = rvolu(10)
887 cpa_flux(iel) = dmi * cpai
888 cpb_flux(iel) = dmi * cpbi
889 cpc_flux(iel) = dmi * cpci
890 rgas_flux(iel) = dmi * rmwi
891 IF (ityp == 8) THEN
892 cpdi = rvolu(56)
893 cpei = rvolu(57)
894 cpfi = rvolu(58)
895 cpd_flux(iel) = dmi * cpdi
896 cpe_flux(iel) = dmi * cpei
897 cpf_flux(iel) = dmi * cpfi
898 ENDIF
899 energy_flux(iel) = dmi * rvolu(63)
900 elfmass(iel) = elfmass(iel) + dmi * dt1
901 elfehpy(iel) = elfehpy(iel) + dmi * dt1 * rvolu(63)
902 elfvel(iel) = elfvel(iel) + uu * area / elarea(iel)
903 eldmass(iel) = dmi
904 eldehpy(iel) = dmi * rvolu(63)
905 ENDIF
906 ELSE IF (iel > nel) THEN
907! TRIANGLE ON INTERNAL SURFACE TAKE POROSITY INTO ACCOUNT
908 IF (porosity(iel - nel) /= zero) THEN
909 n1=elem(1,iel)
910 n2=elem(2,iel)
911 n3=elem(3,iel)
912 IF (tagela(iel) > 0) THEN
913 vx1=vvv(1,n1)
914 vy1=vvv(2,n1)
915 vz1=vvv(3,n1)
916 vx2=vvv(1,n2)
917 vy2=vvv(2,n2)
918 vz2=vvv(3,n2)
919 vx3=vvv(1,n3)
920 vy3=vvv(2,n3)
921 vz3=vvv(3,n3)
922 ELSE
923 vx1=vvva(1,n1)
924 vy1=vvva(2,n1)
925 vz1=vvva(3,n1)
926 vx2=vvva(1,n2)
927 vy2=vvva(2,n2)
928 vz2=vvva(3,n2)
929 vx3=vvva(1,n3)
930 vy3=vvva(2,n3)
931 vz3=vvva(3,n3)
932 ENDIF
933 !grid velocity
934 vvx=third*(vx1+vx2+vx3)
935 vvy=third*(vy1+vy2+vy3)
936 vvz=third*(vz1+vz2+vz3)
937 !normal vector
938 nx=norm(1,iel)
939 ny=norm(2,iel)
940 nz=norm(3,iel)
941 !normal velocity
942 vel(iel)=nx*vvx+ny*vvy+nz*vvz
943 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
944 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
945 IF (i1 == i2) THEN
946 CALL abort()
947 ENDIF
948 area = elarea(iel) * porosity(iel - nel)
949 rho1 = rpolh(i1)
950 re1 = epolh(i1) / pvolu(i1)
951 ux1 = pu(1, i1)
952 uy1 = pu(2, i1)
953 uz1 = pu(3, i1)
954 gama1 = gpolh(i1)
955 rho2 = rpolh(i2)
956 re2 = epolh(i2) / pvolu(i2)
957 ux2 = pu(1, i2)
958 uy2 = pu(2, i2)
959 uz2 = pu(3, i2)
960 gama2 = gpolh(i2)
961 vfx = half * (ux1 + ux2)
962 vfy = half * (uy1 + uy2)
963 vfz = half * (uz1 + uz2)
964 vrel(iel) = (vfx - vvx)*(vfx - vvx) + (vfy - vvy)*(vfy - vvy) + (vfz - vvz)*(vfz - vvz)
965 ss_= nx * (vfx - vvx) + ny * (vfy - vvy) + nz * (vfz - vvz)
966 ss(iel) = elarea(iel)*abs(ss_)
967 IF (ss_ <= zero) THEN
968 alpha = zero
969 ELSE
970 alpha = one
971 ENDIF
972 rhom = alpha * rho1 + (one - alpha) * rho2
973 ruxm = alpha * rho1 * ux1 + (one - alpha) * rho2 * ux2
974 ruym = alpha * rho1 * uy1 + (one - alpha) * rho2 * uy2
975 ruzm = alpha * rho1 * uz1 + (one - alpha) * rho2 * uz2
976 rem = alpha * gama1 * re1 + (one - alpha) * gama2 * re2
977
978 massflow = rhom * ss_ * area
979 mass_flux(iel) = -massflow
980 momentum_flux_x(iel) = -ruxm * ss_ * area
981 momentum_flux_y(iel) = -ruym * ss_ * area
982 momentum_flux_z(iel) = -ruzm * ss_ * area
983 energy_flux(iel) = -rem * ss_ * area
984
985 elfmass(iel) = elfmass(iel) + massflow * dt1
986 elfvel(iel) = elfvel(iel) + ss_ * area
987
988 cpam = alpha * cpapolh(i1) + (one - alpha) * cpapolh(i2)
989 cpbm = alpha * cpbpolh(i1) + (one - alpha) * cpbpolh(i2)
990 cpcm = alpha * cpcpolh(i1) + (one - alpha) * cpcpolh(i2)
991 rgasm = alpha * rmwpolh(i1) + (one - alpha) * rmwpolh(i2)
992 cpa_flux(iel) = -massflow * cpam
993 cpb_flux(iel) = -massflow * cpbm
994 cpc_flux(iel) = -massflow * cpcm
995 rgas_flux(iel) = -massflow * rgasm
996 IF (ityp == 8) THEN
997 cpdm = alpha * cpdpolh(i1) + (one - alpha) * cpdpolh(i2)
998 cpem = alpha * cpepolh(i1) + (one - alpha) * cpepolh(i2)
999 cpfm = alpha * cpfpolh(i1) + (one - alpha) * cpfpolh(i2)
1000 cpd_flux(iel) = -massflow * cpdm
1001 cpe_flux(iel) = -massflow * cpem
1002 cpf_flux(iel) = -massflow * cpfm
1003 ENDIF
1004 ENDIF
1005 ENDIF
1006 ENDDO
1007!$OMP END DO
1008
1009!$OMP END PARALLEL
1010
1011!$OMP PARALLEL PRIVATE(I,MSINI,CPA,CPB,CPC,CPD,CPE,CPF,EFAC,TEMP0)
1012!$OMP+ PRIVATE(TEMP,CP,CV,RMW)
1013!$OMP+ PRIVATE(AL,DD,AD,GAMA,SSP,QX,VV,SS_,AREA)
1014! --- PRIVATE(AREAMAX)
1015!$OMP+ PRIVATE(ITYPL,J,IEL,COUNT)
1016
1017!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1018 DO i = 1, npolh
1019 msini = mpolh(i)
1020 cpapolh(i) = msini * cpapolh(i)
1021 cpbpolh(i) = msini * cpbpolh(i)
1022 cpcpolh(i) = msini * cpcpolh(i)
1023 rmwpolh(i) = msini * rmwpolh(i)
1024 IF (ityp == 8) THEN
1025 cpdpolh(i) = msini * cpdpolh(i)
1026 cpepolh(i) = msini * cpepolh(i)
1027 cpfpolh(i) = msini * cpfpolh(i)
1028 ENDIF
1029 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1030 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1031! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1032 mpolh(i) = mpolh(i) + mass_flux(iel) * dt1
1033 qpolh(1, i) = qpolh(1, i) + momentum_flux_x(iel) * dt1
1034 qpolh(2, i) = qpolh(2, i) + momentum_flux_y(iel) * dt1
1035 qpolh(3, i) = qpolh(3, i) + momentum_flux_z(iel) * dt1
1036 epolh(i) = epolh(i) + energy_flux(iel) * dt1
1037 IF (mpolh(i) <= zero .OR. epolh(i) <= zero) THEN
1038 tpolh(i) = zero
1039 ELSE
1040 cpapolh(i) = cpapolh(i) + cpa_flux(iel) * dt1
1041 cpbpolh(i) = cpbpolh(i) + cpb_flux(iel) * dt1
1042 cpcpolh(i) = cpcpolh(i) + cpc_flux(iel) * dt1
1043 rmwpolh(i) = rmwpolh(i) + rgas_flux(iel) * dt1
1044 IF (ityp == 8) THEN
1045 cpdpolh(i) = cpdpolh(i) + cpd_flux(iel) * dt1
1046 cpepolh(i) = cpepolh(i) + cpe_flux(iel) * dt1
1047 cpfpolh(i) = cpfpolh(i) + cpf_flux(iel) * dt1
1048 ENDIF
1049 ENDIF
1050 ENDDO
1051 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1052 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1053! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1054 mpolh(i) = mpolh(i) - mass_flux(iel) * dt1
1055 qpolh(1, i) = qpolh(1, i) - momentum_flux_x(iel) * dt1
1056 qpolh(2, i) = qpolh(2, i) - momentum_flux_y(iel) * dt1
1057 qpolh(3, i) = qpolh(3, i) - momentum_flux_z(iel) * dt1
1058 epolh(i) = epolh(i) - energy_flux(iel) * dt1
1059 IF (mpolh(i) <= zero .OR. epolh(i) <= zero) THEN
1060 tpolh(i) = zero
1061 ELSE
1062 cpapolh(i) = cpapolh(i) - cpa_flux(iel) * dt1
1063 cpbpolh(i) = cpbpolh(i) - cpb_flux(iel) * dt1
1064 cpcpolh(i) = cpcpolh(i) - cpc_flux(iel) * dt1
1065 rmwpolh(i) = rmwpolh(i) - rgas_flux(iel) * dt1
1066 IF (ityp == 8) THEN
1067 cpdpolh(i) = cpdpolh(i) - cpd_flux(iel) * dt1
1068 cpepolh(i) = cpepolh(i) - cpe_flux(iel) * dt1
1069 cpfpolh(i) = cpfpolh(i) - cpf_flux(iel) * dt1
1070 ENDIF
1071 ENDIF
1072 ENDDO
1073! TEMPERATURE COMPUTATION
1074 efac = epolh(i) / mpolh(i)
1075 cpapolh(i) = cpapolh(i) / mpolh(i)
1076 cpbpolh(i) = cpbpolh(i) / mpolh(i)
1077 cpcpolh(i) = cpcpolh(i) / mpolh(i)
1078 rmwpolh(i) = rmwpolh(i) / mpolh(i)
1079 IF (ityp == 8) THEN
1080 cpdpolh(i) = cpdpolh(i) / mpolh(i)
1081 cpepolh(i) = cpepolh(i) / mpolh(i)
1082 cpfpolh(i) = cpfpolh(i) / mpolh(i)
1083 ENDIF
1084 cpa = cpapolh(i)
1085 cpb = cpbpolh(i)
1086 cpc = cpcpolh(i)
1087 rmw = rmwpolh(i)
1088 IF (ityp == 8) THEN
1089 cpd = cpdpolh(i)
1090 cpe = cpepolh(i)
1091 cpf = cpfpolh(i)
1092 ENDIF
1093 temp0 = rvolu(25)
1094 CALL fvtemp(ityp , efac , cpa , cpb , cpc ,
1095 . cpd , cpe , cpf , rmw , temp0,
1096 . temp )
1097 cp = cpa + cpb * temp + cpc * temp * temp
1098 IF (ityp == 8) THEN
1099 cp = cp + cpd * temp**3 + cpf * temp**4
1100 IF (cpe /= zero .AND. temp > zero) THEN
1101 cp = cp + cpe / (temp * temp)
1102 ENDIF
1103 ENDIF
1104 cv = cp - rmw
1105 rmwpolh(i) = rmw
1106 gpolh(i) = cp / cv
1107 tpolh(i) = temp
1108! PRESSURE AND MASS DENSITY
1109 rpolh(i) = mpolh(i) / pvolu(i)
1110 ppolh(i) = rpolh(i) * rmw * temp
1111 ENDDO
1112!$OMP END DO
1113
1114C---------------------------------------------------
1115C QA, QB AND STABILITY
1116C---------------------------------------------------
1117!$OMP SINGLE
1118
1119 ! /DT/FVMBAG/[ID_DT_OPTION]
1120 ! ID_DT_OPTION = 0 obsolete
1121 ! ID_DT_OPTION = 1 (legacy method), apply to all FVMBAG1
1122 ! ID_DT_OPTION = 2 (2022.3), apply to given FVMBAG1 or FVMBAG2.
1123 ! in this case characteristic length is calculated depending on L_TYPE flag
1124 ! and smoothing is applied : dt <- min(dt, dt_old*(1+lambda))
1125 !--------!-------------------------------------------------------------------!
1126 ! VALUE ! DESCRIPTION !
1127 !------------!--------!-------------------------------------------------------------------!
1128 ! ! 0,1 ! Constant (Starter minimum characteristic length) !
1129 ! L_TYPE ! 2 ! VOLUME^1/3 at current time !
1130 ! ! 3 ! VOLUME / SMAX at current time !
1131 !------------!--------!-------------------------------------------------------------------!
1132
1133
1134 !------------------------------!
1135 idt_fvmbag = fvdata(ifv)%ID_DT_OPTION !IDTMIN(52) ! /DT/FVMBAG/[IDT_FVMBAG]
1136
1137 qa=rvolu(46)
1138 qb=rvolu(47)
1139 dtx=ep30
1140 phdt=0
1141
1142!$OMP END SINGLE
1143
1144 !------------------------------!
1145 ! --- Characteristic Length ---!
1146 !------------------------------!
1147 idt_fvmbag = fvdata(ifv)%ID_DT_OPTION ! /dt/fvmbag/[idt_fvmbag]
1148
1149 SELECT CASE(idt_fvmbag)
1150
1151 ! /DT/FVMBAG/1
1152 CASE(0,1)
1153 ! ------------
1154 param_l_type = 0
1155!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1156 DO i=1,npolh
1157 IF (ibpolh(i) == 0) THEN
1158 dls(i)=dlh !POLYHEDRON CONSTANT CHARACTERISTIC LENGTH COMPUTED IN THE STARTER
1159 ELSE
1160 dls(i)=pvolu(i)**third !ADDITIONAL BRICKS
1161 ENDIF
1162 ENDDO
1163!$OMP END DO
1164
1165 ! /DT/FVMBAG/2
1166 CASE(2)
1167 ! ------------
1168 param_l_type = fvdata(ifv)%L_TYPE
1169
1170 SELECT CASE(param_l_type)
1171 CASE(0,1) !dlh = starter constant
1172 ! /DT/FVMBAG/2 + L_TYPE=0,1
1173!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1174 DO i=1,npolh
1175 dls(i)=dlh
1176 ENDDO
1177!$OMP END DO
1178 CASE(2) ! dl=VOL**1/3
1179 ! /DT/FVMBAG/2 + L_TYPE=2
1180!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1181 DO i=1,npolh
1182 dls(i)=pvolu(i)**third
1183 ENDDO
1184!$OMP END DO
1185 CASE(3) ! dl=VOL**1/3
1186 ! /DT/FVMBAG/2 + L_TYPE=3
1187!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1188 DO i=1,npolh
1189 dls(i)=pvolu(i)
1190 ENDDO
1191!$OMP END DO
1192 END select!PARAM_L_TYPE
1193
1194 END select!IDT_FVMBAG
1195
1196
1197 !------------------------------!
1198 ! --- TIME STEP ---!
1199 !------------------------------!
1200 SELECT CASE (idt_fvmbag)
1201
1202 CASE(0,1)
1203!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1204 DO i=1,npolh
1205 IF(mpolh(i)<=zero.OR.epolh(i)<=zero.OR.gpolh(i)<=one) THEN
1206 dtpolh(i)=ep30
1207 qvisc(i) =zero
1208 ELSE
1209C
1210 al=pvolu(i)**third
1211 dd=dm(i)/mpolh(i)
1212 ad=max(zero,dd)
1213 gama=gpolh(i)
1214 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1215 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1216C
1217 qx=qb*ssp+al*qa*qa*ad
1218 ssp=qx+sqrt(qx*qx+ssp*ssp)
1219 vv=sqrt(pu(1,i)**2+pu(2,i)**2+pu(3,i)**2)
1220 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1221 ENDIF
1222 ENDDO
1223!$OMP END DO
1224
1225 CASE(2) !time step computed from a characteristic length
1226
1227 param_l_type = fvdata(ifv)%L_TYPE !characteristic length formulation
1228
1229 IF(param_l_type /= 3)THEN !L_TYPE=0/1 & 2
1230!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1231 DO i=1,npolh
1232 IF(mpolh(i) <= zero.OR.epolh(i) <= zero.OR.gpolh(i) <= one) THEN
1233 dtpolh(i)=ep30
1234 qvisc(i) =zero
1235 ELSE
1236 ! 1.calculating Grid Velocity for each polyedron
1237 count = 0
1238 vv = zero
1239 !loop over polygons
1240 !AREAMAX=ZERO
1241 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1242 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1243 area = elarea(iel)
1244! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1245 !AREAMAX=MAX(AREAMAX,AREA)
1246 vv = vv+sqrt(vrel(iel))
1247 count=count+1
1248 ENDDO
1249 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1250 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1251! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1252 vv = vv+sqrt(vrel(iel))
1253 count=count+1
1254 ENDDO
1255 al=pvolu(i)**third
1256 dd=dm(i)/mpolh(i)
1257 ad=max(zero,dd)
1258 gama=gpolh(i)
1259 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1260 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1261 qx=qb*ssp+al*qa*qa*ad
1262 ssp=qx+sqrt(qx*qx+ssp*ssp)
1263 ! 3.relative velocty (PU:material velocity : PW:grid velocity)
1264 IF(count > 0)THEN
1265 vv = vv / float(count)
1266 ENDIF
1267 ! 4.time step
1268 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1269 ENDIF
1270 ENDDO
1271!$OMP END DO
1272
1273 ELSE !L_TYPE=3 ! check all faces and compute normal relative velocities
1274!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1275 DO i = 1, npolh
1276 ss_=zero
1277 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1278 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1279! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1280 ss_=max(ss_,ss(iel)) ! <|u-w|,n>
1281 ENDDO
1282 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1283 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1284! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1285 ss_=max(ss_,ss(iel)) ! <|u-w|,n>
1286 ENDDO
1287
1288 al=pvolu(i)**third
1289 dd=dm(i)/mpolh(i)
1290 ad=max(zero,dd)
1291 gama=gpolh(i)
1292 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1293 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1294 qx=qb*ssp+al*qa*qa*ad
1295 ssp=qx+sqrt(qx*qx+ssp*ssp)
1296 ! 3.relative velocty (PU:material velocity : PW:grid velocity)
1297 ! SS_ <- max( area*(u-w).n )
1298 ! u-w = 0 unless in case of porosity
1299 ! 4.time step
1300 dtpolh(i)=cfl_coef*dls(i)/(ssp+ss_)
1301 ENDDO
1302!$OMP END DO
1303 ENDIF
1304
1305 END SELECT
1306
1307
1308!$OMP END PARALLEL
1309
1310 DO i=1,npolh
1311 IF (dtpolh(i)<dtx) THEN
1312 dtx=dtpolh(i)
1313 phdt=i
1314 ENDIF
1315 ENDDO
1316
1317 !optional time step growth limitor (/DT/FVMBAG/2, 3rd parameter)
1318 IF(ncycle > 0 .AND. idt_fvmbag == 2)THEN
1319 lambda = fvdata(ifv)%LAMBDA
1320 IF(lambda > zero)THEN
1321 dtold = fvdata(ifv)%DTOLD
1322 IF(dtold > zero)THEN
1323 dtx = min( dtx, dtold*(one+lambda))
1324 ENDIF
1325 ENDIF
1326 ENDIF
1327 fvdata(ifv)%DTOLD = dtx
1328
1329 IF (ilvout >= 1.AND.ipri==0) THEN
1330 WRITE(iout,'(A,I10,A,E12.4,A,I10)')' ** MONITORED VOLUME ID: ',ivolu(1),' TIME STEP:',dtx,' FINITE VOLUME:',idpolh(phdt)
1331 ENDIF
1332
1333C
1334 IF (dtx<dt2) THEN
1335 dt2=dtx
1336 nelts=ivolu(1)
1337 itypts=52
1338 ENDIF
1339C---------------------------------------------------
1340C PRESSURE FORCES AND WORK COMPUTATION
1341C---------------------------------------------------
1342!$OMP PARALLEL PRIVATE(N1,N2,N3,IEL,NX,NY,NZ)
1343!$OMP+ PRIVATE(VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
1344!$OMP+ PRIVATE(VVX,VVY,VVZ,I1,I2,AREA,IVENT,IDEF,IVENTYP,AOUT,PP,P1,P2,PMEAN,AREA_OLD,AREA1,SS_)
1345!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1346 DO iel=1,nelt
1347 n1=elem(1,iel)
1348 n2=elem(2,iel)
1349 n3=elem(3,iel)
1350 vx1=vvv(1,n1)
1351 vy1=vvv(2,n1)
1352 vz1=vvv(3,n1)
1353 vx2=vvv(1,n2)
1354 vy2=vvv(2,n2)
1355 vz2=vvv(3,n2)
1356 vx3=vvv(1,n3)
1357 vy3=vvv(2,n3)
1358 vz3=vvv(3,n3)
1359 vvx=third*(vx1+vx2+vx3)
1360 vvy=third*(vy1+vy2+vy3)
1361 vvz=third*(vz1+vz2+vz3)
1362 nx=norm(1,iel)
1363 ny=norm(2,iel)
1364 nz=norm(3,iel)
1365 vel(iel)=nx*vvx+ny*vvy+nz*vvz
1366 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1367 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1368 area = elarea(iel)
1369 momentum_flux_x(iel) = zero
1370 momentum_flux_y(iel) = zero
1371 momentum_flux_z(iel) = zero
1372 energy_flux(iel) = zero
1373 fel(iel) = zero
1374 IF (i1 == i2) THEN
1375! EXTERNAL TRIANGLE
1376 IF (itagel(iel) < 0) THEN
1377! VENT HOLE
1378 ivent = -itagel(iel)
1379 idef = ibaghol(1, ivent)
1380 iventyp = ibaghol(13, ivent)
1381 IF (idef == 1 .AND. iventyp == 0) THEN
1382 aout = elsout(iel)
1383 area = area - aout
1384 pp = pext
1385 fel(iel) = pp * aout
1386 momentum_flux_x(iel) = -pp * aout * nx
1387 momentum_flux_y(iel) = -pp * aout * ny
1388 momentum_flux_z(iel) = -pp * aout * nz
1389 energy_flux(iel) = -pp * vel(iel) * aout
1390 ENDIF
1391 ENDIF
1392 pp = ppolh(i1) + qvisc(i1)
1393 fel(iel) = fel(iel) + pp * area
1394 momentum_flux_x(iel) = momentum_flux_x(iel) - pp * area * nx
1395 momentum_flux_y(iel) = momentum_flux_y(iel) - pp * area * ny
1396 momentum_flux_z(iel) = momentum_flux_z(iel) - pp * area * nz
1397 energy_flux(iel) = energy_flux(iel) - pp * vel(iel) * area
1398 energy_flux(iel) = energy_flux(iel) - rvolu(19) * area * (tpolh(i1) - rvolu(25))
1399 ELSE
1400! INTERNAL TRIANGLE, SEPARATES 2 FINITE VOLUMES
1401 p1 = ppolh(i1) + qvisc(i1)
1402 p2 = ppolh(i2) + qvisc(i2)
1403 pmean = half * (p1 + p2)
1404 area_old = area
1405 area = area * porosity(iel - nel)
1406 area1 = area_old - area
1407 ss_ = nx * vvx + ny * vvy + nz * vvz
1408 IF (itagel(iel) > 0) THEN
1409! INJECTOR
1410 fel(iel) = (p1 - p2) * area1
1411 momentum_flux_x(iel) = -pmean * nx * area
1412 momentum_flux_y(iel) = -pmean * ny * area
1413 momentum_flux_z(iel) = -pmean * nz * area
1414 energy_flux(iel) = -pmean * ss_ * area
1415 ELSE
1416 fel(iel) = (p1 - p2) * area1
1417 momentum_flux_x(iel) = -pmean * nx * area
1418 momentum_flux_y(iel) = -pmean * ny * area
1419 momentum_flux_z(iel) = -pmean * nz * area
1420 energy_flux(iel) = -pmean * ss_ * area
1421 ENDIF
1422 ENDIF
1423 ENDDO
1424!$OMP END DO
1425!$OMP END PARALLEL
1426
1427C---------------------------------------------------
1428C NEW QUANTITIES IN POLYEDRA
1429C---------------------------------------------------
1430!$OMP PARALLEL PRIVATE(I, J, IEL)
1431!$OMP+ PRIVATE(NX, NY, NZ, N1, N2, N3, VX1, VX2, VX3, VY1, VY2, VY3, VZ1, VZ2, VZ3)
1432!$OMP+ PRIVATE(VVX, VVY, VVZ, I1, I2, PP, AREA_OLD, AREA, AREA1, SS_)
1433!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1434 DO i = 1, npolh
1435 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1436 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1437 nx=norm(1,iel)
1438 ny=norm(2,iel)
1439 nz=norm(3,iel)
1440 n1=elem(1,iel)
1441 n2=elem(2,iel)
1442 n3=elem(3,iel)
1443 vx1=vvv(1,n1)
1444 vy1=vvv(2,n1)
1445 vz1=vvv(3,n1)
1446 vx2=vvv(1,n2)
1447 vy2=vvv(2,n2)
1448 vz2=vvv(3,n2)
1449 vx3=vvv(1,n3)
1450 vy3=vvv(2,n3)
1451 vz3=vvv(3,n3)
1452 vvx=third*(vx1+vx2+vx3)
1453 vvy=third*(vy1+vy2+vy3)
1454 vvz=third*(vz1+vz2+vz3)
1455 area = elarea(iel)
1456 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1457 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1458 IF (i1 /= i) THEN
1459 CALL abort()
1460 ENDIF
1461 pp = ppolh(i) + qvisc(i)
1462 IF (iel > nel) THEN
1463 area_old = area
1464 area = area * porosity(iel - nel)
1465 area1 = area_old - area
1466 ELSE
1467 area1 = zero
1468 ENDIF
1469
1470 ss_ = nx * vvx + ny * vvy + nz * vvz
1471
1472! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1473 qpolh(1, i) = qpolh(1, i) + (momentum_flux_x(iel) - pp * nx * area1) * dt1
1474 qpolh(2, i) = qpolh(2, i) + (momentum_flux_y(iel) - pp * ny * area1) * dt1
1475 qpolh(3, i) = qpolh(3, i) + (momentum_flux_z(iel) - pp * nz * area1) * dt1
1476 IF (itagel(iel) > 0) THEN
1477 epolh(i) = epolh(i) + (energy_flux(iel) + pp * vel(iel) * area1) * dt1
1478 ELSE
1479 epolh(i) = epolh(i) + (energy_flux(iel) - pp * ss_ * area1) * dt1
1480 ENDIF
1481 ENDDO
1482 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1483 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1484 nx=norm(1,iel)
1485 ny=norm(2,iel)
1486 nz=norm(3,iel)
1487 n1=elem(1,iel)
1488 n2=elem(2,iel)
1489 n3=elem(3,iel)
1490 vx1=vvv(1,n1)
1491 vy1=vvv(2,n1)
1492 vz1=vvv(3,n1)
1493 vx2=vvv(1,n2)
1494 vy2=vvv(2,n2)
1495 vz2=vvv(3,n2)
1496 vx3=vvv(1,n3)
1497 vy3=vvv(2,n3)
1498 vz3=vvv(3,n3)
1499 vvx=third*(vx1+vx2+vx3)
1500 vvy=third*(vy1+vy2+vy3)
1501 vvz=third*(vz1+vz2+vz3)
1502 area = elarea(iel)
1503 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1504 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1505 IF (i2 /= i) THEN
1506 CALL abort()
1507 ENDIF
1508 pp = ppolh(i) + qvisc(i)
1509 IF (iel > nel) THEN
1510 area_old = area
1511 area = area * porosity(iel - nel)
1512 area1 = area_old - area
1513 ELSE
1514 area1 = zero
1515 ENDIF
1516
1517 ss_ = nx * vvx + ny * vvy + nz * vvz
1518
1519
1520! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1521 qpolh(1, i) = qpolh(1, i) - (momentum_flux_x(iel) - pp * nx * area1) * dt1
1522 qpolh(2, i) = qpolh(2, i) - (momentum_flux_y(iel) - pp * ny * area1) * dt1
1523 qpolh(3, i) = qpolh(3, i) - (momentum_flux_z(iel) - pp * nz * area1) * dt1
1524 IF (itagel(iel) > 0) THEN
1525 epolh(i) = epolh(i) - (energy_flux(iel) - pp * vel(iel) * area1) * dt1
1526 ELSE
1527 epolh(i) = epolh(i) - (energy_flux(iel) - pp * ss_ * area1) * dt1
1528 ENDIF
1529 ENDDO
1530 ENDDO
1531!$OMP END DO
1532!$OMP END PARALLEL
1533
1534 ENDIF
1535
1536C---------------------------------------------------
1537C PRINT OUTS
1538C---------------------------------------------------
1539 250 IF (ilvout >= 1 .AND. ipri == 0) THEN
1540 WRITE(iout,*)
1541 WRITE(iout,'(A25,I10,A31)')' ** MONITORED VOLUME ID: ',ivolu(1),' - FINITE VOLUME MESH UPDATE **'
1542 WRITE(iout,'(A,I8)')' NUMBER OF FINITE VOLUMES : ',npolh
1543 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM VOLUME FINITE VOLUMES : ',volph,' (VOLUME AIRBAG: ',volg,')'
1544 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM AREA SURFACE POLYGONS : ',areap,' (AREA AIRBAG : ',areael,')'
1545 WRITE(iout,*)
1546 ENDIF
1547C---------------------------------------------------
1548C TIME HISTORY
1549C---------------------------------------------------
1550 amtot=zero
1551 pmean=zero
1552 enint=zero
1553 fsav(1:17)=zero
1554 DO i=1,npolh
1555 amtot=amtot+mpolh(i)
1556 enint=enint+epolh(i)
1557 pmean=pmean+ppolh(i)*pvolu(i)
1558 gama=gpolh(i)
1559 rmw=rmwpolh(i)
1560 temp=tpolh(i)
1561 cpa=cpapolh(i)
1562 cpb=cpbpolh(i)
1563 cpc=cpcpolh(i)
1564 cp=cpa+cpb*temp+cpc*temp*temp
1565 IF(ityp==8) THEN
1566 cpd=cpdpolh(i)
1567 cpe=cpepolh(i)
1568 cpf=cpfpolh(i)
1569 cp=cp+cpd*temp**3+cpf*temp**4
1570 IF(temp > zero) cp=cp+cpe/(temp*temp)
1571 ENDIF
1572 cv=cp-rmw
1573 fsav(5) =fsav(5) +mpolh(i)*temp
1574 fsav(10)=fsav(10)+mpolh(i)*cp
1575 fsav(11)=fsav(11)+mpolh(i)*cv
1576 fsav(12)=fsav(12)+mpolh(i)*gama
1577 ENDDO
1578C PRESSURE DISPERSION
1579 pdisp = zero
1580 pmean2 = pmean / volph
1581 DO i = 1, npolh
1582 pdisp = pdisp + pvolu(i) * (ppolh(i) - pmean2)**2
1583 ENDDO
1584 pdisp = sqrt(pdisp / volph) / pmean2
1585 fsav(19) = pdisp
1586 injection_started = .false.
1587 DO iinj = 1, njet
1588 injection_started = injection_started .OR. (datainj(2, iinj) > zero)
1589 ENDDO
1590 IF (.NOT. injection_started) THEN
1591 pdisp = zero
1592 ENDIF
1593
1594C
1595 fsav(1)=amtot
1596 fsav(2)=volph
1597 pres_av = pmean / volph
1598 fsav(3) = pres_av
1599 fsav(4)=areap
1600 fsav(7)=zero
1601 dmout =zero
1602 dhout =zero
1603 cp_av = fsav(10) / amtot
1604 fsav(10) = cp_av
1605 cv_av = fsav(11) / amtot
1606 fsav(11) = cv_av
1607 fsav(12) = cp_av / cv_av
1608 fsav(14) = npolh
1609 rgas_av = cp_av - cv_av
1610 rho_av = amtot / volph
1611 fsav(5) = fsav(5) / amtot ! sum(mass[i]*T[i])/total_mass
1612
1613 IF(ivolu(39) == 0) THEN
1614 fsav(6) =zero
1615 fsav(13)=zero
1616 ELSE
1617 fsav(6) =aoutot
1618 DO iel=1,nel
1619 IF (itagel(iel)<0) THEN
1620 ivent=-itagel(iel)
1621 rbaghol(18,ivent)=rbaghol(18,ivent)+elfvel(iel)*elsout(iel)
1622 rbaghol(19,ivent)=rbaghol(19,ivent)-elfmass(iel)
1623 rbaghol(20,ivent)=rbaghol(20,ivent)-elfehpy(iel)
1624 rbaghol(21,ivent)=rbaghol(21,ivent)+eldmass(iel)
1625 rbaghol(22,ivent)=rbaghol(22,ivent)+eldehpy(iel)
1626 ENDIF
1627 ENDDO
1628 DO ivent=1,nvent
1629 sventot=rbaghol(16,ivent)+rbaghol(17,ivent)
1630 IF(sventot <= zero) cycle
1631 rbaghol(18,ivent)=rbaghol(18,ivent)/sventot
1632 ENDDO
1633 DO ivent=1,nvent
1634 fsav(7)=fsav(7)+(rbaghol(16,ivent)+rbaghol(17,ivent))*rbaghol(18,ivent)
1635 dmout =dmout+rbaghol(21,ivent)
1636 dhout =dhout+rbaghol(22,ivent)
1637 ENDDO
1638 fsav(7) =fsav(7) /max(aoutot,em20)
1639 fsav(13)=dtx ! not DTPOLH(PHDT) since DTX may be updated with smoothing factor
1640
1641 IF(ityp == 8) THEN
1642 CALL fvinjt8_1(njet, ibagjet , rbagjet , igeo, geo, pm,
1643 . ivolu, rvolu, dmout, dhout)
1644 ttf =rvolu(60)
1645 IF (ivolu(74) == 0) THEN
1646 up_switch = tt-ttf >= rvolu(70) .OR. npolh <= ivolu(37)
1647 ENDIF
1648 IF (ivolu(74) == 1) THEN
1649C AUTOMATIC SWITCH TO UNIFORM PRESSURE WHEN DISPERSION OF PRESSURE IS LOW
1650 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1651 up_switch = up_switch .OR. tt-ttf >= rvolu(70)
1652 ENDIF
1653 IF (ivolu(74) == 2) THEN
1654C AUTOMATIC SWITCH TO UNIFORM PRESSURE WHEN DISPERSION OF PRESSURE IS LOW
1655 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1656 up_switch = up_switch .AND. tt-ttf >= rvolu(70)
1657 ENDIF
1658
1659 IF (up_switch) THEN
1660 rvolu(12)=fsav(3) ! PRESSION
1661 rvolu(13)=fsav(5) ! TEMPERATURE
1662 rvolu(16)=fsav(2) ! VOLUME
1663 rvolu(20)=fsav(1) ! MASSE
1664 ENDIF
1665 ENDIF
1666
1667 ENDIF
1668
1669 DO iel=1,nelt
1670 IF (itagel(iel) > 0) THEN
1671 fsav(15)=fsav(15)+elfmass(iel)
1672 fsav(16)=fsav(16)+elfehpy(iel)
1673 ENDIF
1674 ENDDO
1675 fsav(17)=enint
1676
1677
1678C---------------------------------------------------
1679C NODAL FLUID VELOCITIES
1680C---------------------------------------------------
1681 IF( (tt>=output%TANIM .AND. tt<=output%TANIM_STOP) .OR.tt>=toutp.OR. nbgauge > 0 .OR.
1682 . (manim>=4.AND.manim<=15) ) THEN
1683C COMPUTATIONS FOR ANIM, ALSO NEEDED IF THERE ARE GAUGES.
1684 prepare_anim = 1
1685 uel(1:3,1:nelt)=zero
1686 rhoel(1:nelt)=zero
1687 tkel(1:nelt) =zero
1688 sspel(1:nelt)=zero
1689 ELSE
1690 prepare_anim = 0
1691 ENDIF
1692
1693C PEL(1:NELT)=ZERO
1694
1695!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,II,NV)
1696!$OMP+ PRIVATE(NX,NY,NZ,NRX,NRY,NRZ)
1697!$OMP+ PRIVATE(FAC,AREA,I1,I2)
1698!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1699 DO iel = 1, nelt
1700 area = elarea(iel)
1701 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1702 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1703 pel(iel) = half * (ppolh(i1) + ppolh(i2))
1704 IF (prepare_anim == 1) THEN
1705 uel(1, iel) = half * (pu(1, i1) + pu(1, i2))
1706 uel(2, iel) = half * (pu(2, i1) + pu(2, i2))
1707 uel(3, iel) = half * (pu(3, i1) + pu(3, i2))
1708 rhoel(iel) = half * (rpolh(i1) + rpolh(i2))
1709 tkel(iel) = half * (tpolh(i1) + tpolh(i2))
1710 sspel(iel) = half * (sqrt(gpolh(i1)*ppolh(i1)/rpolh(i1)) + sqrt(gpolh(i2)*ppolh(i2)/rpolh(i2)))
1711 ENDIF
1712 ENDDO
1713!$OMP END DO
1714
1715C------------------
1716C VENT HOLES
1717C------------------
1718 IF(prepare_anim == 1) THEN
1719!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1720 DO iel=1,nel
1721 IF(itagel(iel)<0) THEN
1722 uel(1,iel)=elfvel(iel)*norm(1,iel)
1723 uel(2,iel)=elfvel(iel)*norm(2,iel)
1724 uel(3,iel)=elfvel(iel)*norm(3,iel)
1725 ENDIF
1726 ENDDO
1727!$OMP END DO NOWAIT
1728!$OMP SINGLE
1729 u(1:3,1:nnt)=zero
1730 rho(1:nnt) =zero
1731 tk(1:nnt) =zero
1732 sspk(1:nnt) =zero
1733 vola(1:nna)=zero
1734 pa(1:nna)=zero
1735 rhoa(1:nna)=zero
1736 tka(1:nna)=zero
1737 sspka(1:nna)=zero
1738 ua(1:3,1:nna)=zero
1739!$OMP END SINGLE
1740 ENDIF
1741
1742!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1743 DO i=1,nnt
1744 p(i)=zero
1745 ENDDO
1746!$OMP END DO
1747!$OMP END PARALLEL
1748
1749 DO iel=1,nelt
1750 n1=elem(1,iel)
1751 n2=elem(2,iel)
1752 n3=elem(3,iel)
1753 area1=nodarea(n1)
1754 area2=nodarea(n2)
1755 area3=nodarea(n3)
1756 fac1=third*elarea(iel)/area1
1757 fac2=third*elarea(iel)/area2
1758 fac3=third*elarea(iel)/area3
1759 p(n1) =p(n1) +fac1*pel(iel)
1760 p(n2) =p(n2) +fac2*pel(iel)
1761 p(n3) =p(n3) +fac3*pel(iel)
1762 ENDDO
1763
1764 IF(prepare_anim == 1) THEN
1765 DO iel=1,nelt
1766 n1=elem(1,iel)
1767 n2=elem(2,iel)
1768 n3=elem(3,iel)
1769 area1=nodarea(n1)
1770 area2=nodarea(n2)
1771 area3=nodarea(n3)
1772 fac1=third*elarea(iel)/area1
1773 fac2=third*elarea(iel)/area2
1774 fac3=third*elarea(iel)/area3
1775 u(1,n1)=u(1,n1)+fac1*uel(1,iel)
1776 u(2,n1)=u(2,n1)+fac1*uel(2,iel)
1777 u(3,n1)=u(3,n1)+fac1*uel(3,iel)
1778 u(1,n2)=u(1,n2)+fac2*uel(1,iel)
1779 u(2,n2)=u(2,n2)+fac2*uel(2,iel)
1780 u(3,n2)=u(3,n2)+fac2*uel(3,iel)
1781 u(1,n3)=u(1,n3)+fac3*uel(1,iel)
1782 u(2,n3)=u(2,n3)+fac3*uel(2,iel)
1783 u(3,n3)=u(3,n3)+fac3*uel(3,iel)
1784 rho(n1)=rho(n1)+fac1*rhoel(iel)
1785 rho(n2)=rho(n2)+fac2*rhoel(iel)
1786 rho(n3)=rho(n3)+fac3*rhoel(iel)
1787 tk(n1)=tk(n1)+fac1*tkel(iel)
1788 tk(n2)=tk(n2)+fac2*tkel(iel)
1789 tk(n3)=tk(n3)+fac3*tkel(iel)
1790 sspk(n1)=sspk(n1)+fac1*sspel(iel)
1791 sspk(n2)=sspk(n2)+fac2*sspel(iel)
1792 sspk(n3)=sspk(n3)+fac3*sspel(iel)
1793 ENDDO
1794
1795 DO i=1,npolh
1796 IF (ibpolh(i)/=0) THEN
1797 temp=tpolh(i)
1798 ii=iabs(ibpolh(i))
1799 IF(brna(1,ii)==brna(4,ii).AND.brna(5,ii)==brna(8,ii))THEN
1800C PENTAHEDRON
1801 DO k=1,6
1802 j=k
1803 IF(k>=4) j=j+1
1804 jj=brna(j,ii)
1805 vola(jj)=vola(jj)+one_over_6*pvolu(i)
1806 pa(jj)=pa(jj)+one_over_6*ppolh(i)*pvolu(i)
1807 rhoa(jj)=rhoa(jj)+one_over_6*rpolh(i)*pvolu(i)
1808 tka(jj)=tka(jj)+one_over_6*temp*pvolu(i)
1809 sspka(jj)=sspka(jj)+one_over_6*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1810 ua(1,jj)=ua(1,jj)+one_over_6*pu(1,i)*pvolu(i)
1811 ua(2,jj)=ua(2,jj)+one_over_6*pu(2,i)*pvolu(i)
1812 ua(3,jj)=ua(3,jj)+one_over_6*pu(3,i)*pvolu(i)
1813 ENDDO
1814 ELSEIF(brna(5,ii)==brna(8,ii).AND.
1815 . brna(6,ii)==brna(8,ii).AND.
1816 . brna(7,ii)==brna(8,ii))THEN
1817C PYRAMID
1818 DO j=1,5
1819 jj=brna(j,ii)
1820 vola(jj)=vola(jj)+one_fifth*pvolu(i)
1821 pa(jj)=pa(jj)+one_fifth*ppolh(i)*pvolu(i)
1822 rhoa(jj)=rhoa(jj)+one_fifth*rpolh(i)*pvolu(i)
1823 tka(jj)=tka(jj)+one_fifth*temp*pvolu(i)
1824 sspka(jj)=sspka(jj)+one_fifth*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1825 ua(1,jj)=ua(1,jj)+one_fifth*pu(1,i)*pvolu(i)
1826 ua(2,jj)=ua(2,jj)+one_fifth*pu(2,i)*pvolu(i)
1827 ua(3,jj)=ua(3,jj)+one_fifth*pu(3,i)*pvolu(i)
1828 ENDDO
1829 ELSE
1830C HEXAHEDRON OR TETRAEDRON
1831 DO j=1,8
1832 jj=brna(j,ii)
1833 vola(jj)=vola(jj)+one_over_8*pvolu(i)
1834 pa(jj)=pa(jj)+one_over_8*ppolh(i)*pvolu(i)
1835 rhoa(jj)=rhoa(jj)+one_over_8*rpolh(i)*pvolu(i)
1836 tka(jj)=tka(jj)+one_over_8*temp*pvolu(i)
1837 sspka(jj)=sspka(jj)+one_over_8*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1838 ua(1,jj)=ua(1,jj)+one_over_8*pu(1,i)*pvolu(i)
1839 ua(2,jj)=ua(2,jj)+one_over_8*pu(2,i)*pvolu(i)
1840 ua(3,jj)=ua(3,jj)+one_over_8*pu(3,i)*pvolu(i)
1841 ENDDO
1842 ENDIF
1843 ENDIF
1844 ENDDO
1845
1846 DO i=1,nna
1847 volno=vola(i)
1848 IF(volno > zero) THEN
1849 pa(i) =pa(i)/volno
1850 rhoa(i)=rhoa(i)/volno
1851 tka(i) =tka(i)/volno
1852 sspka(i) =sspka(i)/volno
1853 ua(1,i)=ua(1,i)/volno
1854 ua(2,i)=ua(2,i)/volno
1855 ua(3,i)=ua(3,i)/volno
1856 ENDIF
1857 ENDDO
1858 ENDIF ! PREPARE_ANIM
1859
1860
1861C---------------------------------------------------
1862C OPENING OF VENTHOLES
1863C---------------------------------------------------
1864 DO ivent=1,nvent
1865 area_vent(ivent)=zero
1866 pm_vent(ivent)=zero
1867 ENDDO
1868
1869C
1870 DO iel=1,nel
1871 IF (itagel(iel)<0) THEN
1872 ivent=-itagel(iel)
1873 area_vent(ivent)=area_vent(ivent)+elarea(iel)
1874 pm_vent(ivent)=pm_vent(ivent)+fel(iel)
1875 ENDIF
1876 ENDDO
1877
1878
1879
1880 DO ivent=1,nvent
1881 IF (area_vent(ivent)>zero) THEN
1882 pm_vent(ivent)=pm_vent(ivent)/area_vent(ivent)
1883 ELSE
1884 pm_vent(ivent)=zero
1885 ENDIF
1886C
1887 idef =ibaghol(1,ivent)
1888 idtpdef= ibaghol(11,ivent)
1889 idpdef = ibaghol(12,ivent)
1890C
1891 IF(idef==2) cycle
1892C
1893 pdef =rbaghol(1,ivent)
1894 dtpdefi=rbaghol(4,ivent)
1895 dtpdefc=rbaghol(5,ivent)
1896 tvent =rbaghol(3,ivent)
1897 tstope =rbaghol(14,ivent)
1898 IF (idtpdef==0) THEN
1899 IF(pm_vent(ivent)>pdef+pext) THEN
1900 dtpdefc = dtpdefc+dt1
1901 END IF
1902 ELSE IF (idtpdef==1) THEN
1903 IF (pm_vent(ivent)>pdef+pext) THEN
1904 idpdef = 1
1905 END IF
1906 IF (idpdef==1) THEN
1907 dtpdefc = dtpdefc+dt1
1908 END IF
1909 END IF
1910 ibaghol(12,ivent) = idpdef
1911 rbaghol(5,ivent) = dtpdefc
1912C
1913 ittf = ivolu(17)
1914 ttf =rvolu(60)
1915 DO k=1,20
1916 titrevent(k)=ibaghol(14+k,ivent)
1917 venttitle(k:k) = achar(titrevent(k))
1918 ENDDO
1919 IF(ittf==11.OR.ittf==12.OR.ittf==13) THEN
1920 IF (idef==0.AND.pm_vent(ivent)>pdef+pext
1921 . .AND.dtpdefc>dtpdefi
1922 . .AND.volph>em3*areap**three_half
1923 . .AND.tt<tstope+ttf) THEN
1924 idef=1
1925 WRITE(iout,'(A)')
1926 . ' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
1927 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1928 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1929 WRITE(istdo,'(2A)')
1930 . ' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
1931 ENDIF
1932 IF(idef==0 .AND. tt>tvent+ttf
1933 . .AND. tt<tstope+ttf) THEN
1934 idef=1
1935 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
1936 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1937 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1938 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
1939 ENDIF
1940 IF(idef==1 .AND. tt>=tstope+ttf) THEN
1941 idef=2
1942 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
1943 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1944 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1945 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
1946 ENDIF
1947C
1948 ELSE IF(ittf==0) THEN
1949 IF (idef==0.AND.pm_vent(ivent)>pdef+pext
1950 . .AND.dtpdefc>dtpdefi
1951 . .AND.volph>em3*areap**three_half
1952 . .AND.tt<tstope) THEN
1953 idef=1
1954 WRITE(iout,'(A)')
1955 . ' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
1956 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1957 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1958 WRITE(istdo,'(2A)')
1959 . ' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
1960 ENDIF
1961 IF(idef==0 .AND. tt>tvent
1962 . .AND. tt<tstope) THEN
1963 idef=1
1964 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
1965 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1966 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1967 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
1968 ENDIF
1969 IF(idef==1 .AND. tt>=tstope) THEN
1970 idef=2
1971 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
1972 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1973 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1974 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
1975 ENDIF
1976 END IF
1977C
1978 ibaghol(1,ivent)=idef
1979 ENDDO
1980C---------------------------------------------------
1981C FORCES ON AIRBAGS AND OUTPUTS
1982C---------------------------------------------------
1983C FVSPMD(IFV)%AAA(1:3,1:NNT) = ZERO
1984C
1985 IF(ivolu(39)==1) THEN
1986 wfext0=wfext
1987 DO iel=1,nelt
1988 n1=elem(1,iel)
1989 n2=elem(2,iel)
1990 n3=elem(3,iel)
1991 nrx=norm(1,iel)
1992 nry=norm(2,iel)
1993 nrz=norm(3,iel)
1994 IF(iel<=nel) THEN
1995 ff=fel(iel)-pext*elarea(iel)
1996 ELSE
1997 ff=fel(iel)
1998 ENDIF
1999 fvspmd(ifv)%AAA(1,n1)=fvspmd(ifv)%AAA(1,n1)+third*ff*nrx
2000 fvspmd(ifv)%AAA(2,n1)=fvspmd(ifv)%AAA(2,n1)+third*ff*nry
2001 fvspmd(ifv)%AAA(3,n1)=fvspmd(ifv)%AAA(3,n1)+third*ff*nrz
2002 fvspmd(ifv)%AAA(1,n2)=fvspmd(ifv)%AAA(1,n2)+third*ff*nrx
2003 fvspmd(ifv)%AAA(2,n2)=fvspmd(ifv)%AAA(2,n2)+third*ff*nry
2004 fvspmd(ifv)%AAA(3,n2)=fvspmd(ifv)%AAA(3,n2)+third*ff*nrz
2005 fvspmd(ifv)%AAA(1,n3)=fvspmd(ifv)%AAA(1,n3)+third*ff*nrx
2006 fvspmd(ifv)%AAA(2,n3)=fvspmd(ifv)%AAA(2,n3)+third*ff*nry
2007 fvspmd(ifv)%AAA(3,n3)=fvspmd(ifv)%AAA(3,n3)+third*ff*nrz
2008C
2009 wfext=wfext+ff*vel(iel)*dt1
2010 ENDDO
2011 fsav(18)=fsav(18)+wfext-wfext0
2012 ENDIF
2013C----------------------
2014C OUTPUT TH GAUGES
2015C----------------------
2016 IF (nbgauge > 0) THEN
2017 DO i=1,nnt
2018 fvspmd(ifv)%GGG(1,i)=p(i)
2019 fvspmd(ifv)%GGG(2,i)=rho(i)
2020 fvspmd(ifv)%GGG(3,i)=tk(i)
2021 ENDDO
2022
2023 DO i=1,nna
2024 fvspmd(ifv)%GGA(1,i)=pa(i)
2025 fvspmd(ifv)%GGA(2,i)=rhoa(i)
2026 fvspmd(ifv)%GGA(3,i)=tka(i)
2027 ENDDO
2028 ENDIF
2029
2030 RETURN
2031 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
fv_up_switch
subroutine fv_up_switch(output, nn, nel, elem, njet, npoly, lenh, nba, ibagjet, rbagjet, nvent, ibaghol, rbaghol, p, rho, tk, u, sspk, x, v, a, nsensor, sensor_tab, fsav, npc, tf, ivolu, rvolu, mpolh, qpolh, epolh, ppolh, rpolh, gpolh, npolh, ifvnod, rfvnod, ifvtri, ifvpoly, ifvtadr, ifvpolh, ifvpadr, info, nns, nntr, ifv, npolha, dlh, cpapolh, cpbpolh, cpcpolh, rmwpolh, itagel, elsini, icontact, idpolh, elfmass, elfvel, ibufa, elema, tagela, pa, rhoa, tka, ua, brna, nna, ntga, ibpolh, dtpolh, nnt, nelt, xxxa, vvva, ncona, porosity, ityp, igeo, sspka, geo, pm, ipm, tpolh, elfehpy, cpdpolh, cpepolh, cpfpolh, eltg, iparg, mattg, igrouptg, igroupc, elbuf_tab, cfl_coef, pdisp_old, pdisp, wfext, python)
Definition fv_up_switch.F:71
fvinjt6
subroutine fvinjt6(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, datainj, python)
Definition fvinjt6.F:35
fvinjt8
subroutine fvinjt8(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, igeo, geo, pm, ivolu, datainj, python)
Definition fvinjt8.F:36
fvinjt8_1
subroutine fvinjt8_1(njet, ibagjet, rbagjet, igeo, geo, pm, ivolu, rvolu, dmout, dhout)
Definition fvinjt8_1.F:32
fvtemp
subroutine fvtemp(ityp, efac_in, cpa_in, cpb_in, cpc_in, cpd_in, cpe_in, cpf_in, rmw_in, temp0_in, temp_in)
Definition fvtemp.F:35
fvvent0
subroutine fvvent0(elsout, aoutot, nvent, nelt, ittf, elarea, elsini, elem, itagel, svent, ibaghol, rvolu, rbaghol, poro, p, eltg, iparg, mattg, nel, porosity, ipm, pm, elbuf_tab, igroupc, igrouptg)
Definition fvvent0.F:45
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
fvbag_mod
Definition fvbag_mod.F:44
fvbag_mod::fvspmd
type(fvbag_spmd), dimension(:), allocatable fvspmd
Definition fvbag_mod.F:129
fvbag_mod::fvdata
type(fvbag_data), dimension(:), allocatable fvdata
Definition fvbag_mod.F:128
fvmbag_meshcontrol_mod
Definition fvmbag_meshcontrol_mod.F:46
message_mod
Definition message_mod.F:1257
poro
subroutine poro(geo, nodpor, ms, x, v, w, af, am, skew, weight, nporgeo)
Definition poro.F:40
section
subroutine section(nnod, n1, n2, n3, nstrf, x, v, vr, fsav, fopta, secfcum, ms, in, ifram, xsec)
Definition section.F:34
spmd_fvb_gath
subroutine spmd_fvb_gath(ifv, x, xxx, xxxa, xxxsa, ido)
Definition spmd_fvb_gath.F:40
spmd_fvb_igath
subroutine spmd_fvb_igath(ifv, itabs, itabg)
Definition spmd_fvb_igath.F:34
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895
arret
subroutine arret(nn)
Definition arret.F:86