OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
aflux3_int22_fvm.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "inter22.inc"
#include "task_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine aflux3_int22_fvm (pm, ixs, flux, flu1, iparg, elbuf_tab, itask, nv46, ipm, x, w)

Function/Subroutine Documentation

◆ aflux3_int22_fvm()

subroutine aflux3_int22_fvm ( pm,
integer, dimension(nixs,*) ixs,
flux,
flu1,
integer, dimension(nparg,*) iparg,
type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
integer itask,
integer nv46,
integer, dimension(npropmi,*) ipm,
x,
dimension(3,numnod), intent(in) w )

Definition at line 35 of file aflux3_int22_fvm.F.

38C-----------------------------------------------
39C D e s c r i p t i o n
40C-----------------------------------------------
41C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
42C which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
43C This cut cell method is not completed, abandoned, and is not an official option.
44C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
45C
46C This subroutine is computing fluxes for cut cells.
47C Result are stored in cut cell buffer (currently called BRICK_LIST)
48C Main Cell flux is already computed with the usual
49C subroutine (EFLUX3)
50C All Cell not in cut cell buffer only have uncut adjacent cells. Since
51C a cut cell has all its adjacent cells in the cut cell buffer.
52C Cells in cut cell buffer might have adjacent cut cells.
53 !IPM(251,MAT(1)) = I_ALE_SOLVER
54 ! 0: Default = 1 expect if /ALE/MAT or /EULER/MAT has IFROM flag defined.
55 ! 1 : FEM
56 ! 2 : FVM U average
57 ! 3 : FVM rho.U average
58 ! 4 : FVM rho.c.U average
59 ! 5 : Godunov Acoustic
60 ! 6 : experimental
61C-----------------------------------------------
62C M o d u l e s
63C-----------------------------------------------
64 USE i22bufbric_mod
65 USE i22tri_mod
66 USE elbufdef_mod
67 USE alefvm_mod
68C-----------------------------------------------
69C I m p l i c i t T y p e s
70C-----------------------------------------------
71#include "implicit_f.inc"
72C-----------------------------------------------
73C G l o b a l P a r a m e t e r s
74C-----------------------------------------------
75#include "mvsiz_p.inc"
76C-----------------------------------------------
77C C o m m o n B l o c k s
78C-----------------------------------------------
79#include "param_c.inc"
80#include "inter22.inc"
81#include "task_c.inc"
82#include "com01_c.inc"
83#include "com04_c.inc"
84#include "com08_c.inc"
85C-----------------------------------------------
86C D u m m y A r g u m e n t s
87C-----------------------------------------------
88 INTEGER :: IXS(NIXS,*), IPARG(NPARG,*),ISILENT, NV46,IPM(NPROPMI,*)
89 my_real :: pm(npropm,*),flux(6,*), flu1(*),x(3,*)
90 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP),TARGET :: ELBUF_TAB
91 my_real, INTENT(IN) :: w(3,numnod)
92C-----------------------------------------------
93C L o c a l V a r i a b l e s
94C-----------------------------------------------
95 INTEGER MAT, MLW, NC(8), I, IE,J, K, Kv,L,ITASK, IDm, IDV, II
96 INTEGER IDvm, IBvm, NGvm, IEvm, IEV, Jm, IMAT, IALEFVM_FLG
97 INTEGER IB,IBv, NIN, NBCUT, ICELL,NCELL,NGm
98 my_real :: cellflux(6,9,nb,5),reduc,upwl(6)
99 my_real :: vf(3), norm(3,6), lnorm(6),term2
100 INTEGER :: NBF,NBL, MCELL,iNOD,ICELLv, numnod_, numnod_V
101 TYPE(G_BUFEL_) ,POINTER :: GBUF
102 my_real :: face , z(3), zadj(3), zzadj_, cf(3), zcf(3),zzadj(3)
103 my_real :: ps, lambda
104 INTEGER :: NUM, NADJ, IADJ, JV, NG, IE_M, IEm, IBm
105 my_real :: facev, ddvol, valel(6), valvois(6,6,5), sr1, sr2, srf
106 my_real :: wface(3), wfacev(3), wfaceb(3,6)
107 LOGICAL debug_outp
108 INTEGER :: NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8
109C-----------------------------------------------
110C P r e - C o n d i t i o n
111C-----------------------------------------------
112 IF(int22 == 0)RETURN
113C-----------------------------------------------
114C S o u r c e L i n e s
115C-----------------------------------------------
116 ibvm = 0
117 !======================================================!
118 ! INITIALIZATIONS !
119 !======================================================!
120 nin = 1
121 nbf = 1+itask*nb/nthread
122 nbl = (itask+1)*nb/nthread
123 nbl = min(nbl,nb)
124
125 !======================================================!
126 ! DEBUG OUTPUT !
127 !======================================================!
128 !INTERFACE 22 ONLY - OUTPUT---------------!
129 debug_outp = .false.
130 if(ibug22_flux22/=0)then
131 debug_outp = .false.
132 if(ibug22_flux22>0)then
133 do ib=nbf,nbl
134 ie = brick_list(nin,ib)%id
135 if(ixs(11,ie)==ibug22_flux22)debug_outp=.true.
136 enddo
137 elseif(ibug22_flux22==-1)then
138 debug_outp = .true.
139 endif
140 endif
141 if(debug_outp)then
142 print *, " |---------eflux3_int22_fvm.F---------|"
143 print *, " | THREAD INFORMATION |"
144 print *, " |------------------------------------|"
145 print *, " NCYCLE =", ncycle
146 endif
147 !INTERFACE 22 ONLY - OUTPUT---------------!
148
149
150
151 !======================================================!
152 ! STEP A : CUT CELL FLUXES !
153 ! A.1 : main & SECND CELLS FLUXES !
154 !======================================================!
155 DO ib=nbf,nbl
156 ie = brick_list(nin,ib)%ID
157 nbcut = brick_list(nin,ib)%NBCUT
158 mlw = brick_list(nin,ib)%MLW
159 ncell = brick_list(nin,ib)%NBCUT
160 mcell = brick_list(nin,ib)%MainID
161 icell = 0
162
163 !======================================================!
164 ! RESET FLUXES !
165 !======================================================!
166 flux(1:6,ie) = zero
167 flu1(ie) = zero
168
169 !======================================================!
170 ! NORMAL VECTORS ON EACH DACE !
171 ! 2S[n] = [diag1] x [diag2] !
172 ! where !
173 ! [n] : unitary normal vector on face !
174 !======================================================!
175 norm(1,1:6) = brick_list(nin,ib)%N(1:6,1) !VEUL(14:19,IE)
176 norm(2,1:6) = brick_list(nin,ib)%N(1:6,2) !VEUL(20:25,IE)
177 norm(3,1:6) = brick_list(nin,ib)%N(1:6,3) !VEUL(26:31,IE)
178 DO j=1,nv46
179 lnorm(j) = sqrt( norm(1,j)**2 + norm(2,j)**2 + norm(3,j)**2 )
180 norm(1:3,j) = norm(1:3,j) / lnorm(j)
181 ENDDO
182
183 !======================================================!
184 ! FACE VELOCITIES !
185 !======================================================!
186
187 ii = ie
188 nc1 = ixs(2,ii)
189 nc2 = ixs(3,ii)
190 nc3 = ixs(4,ii)
191 nc4 = ixs(5,ii)
192 nc5 = ixs(6,ii)
193 nc6 = ixs(7,ii)
194 nc7 = ixs(8,ii)
195 nc8 = ixs(9,ii)
196
197 wfaceb(1,1) = fourth*(w(1,nc1)+w(1,nc2)+w(1,nc3)+w(1,nc4))
198 wfaceb(2,1) = fourth*(w(2,nc1)+w(2,nc2)+w(2,nc3)+w(2,nc4))
199 wfaceb(3,1) = fourth*(w(3,nc1)+w(3,nc2)+w(3,nc3)+w(3,nc4))
200
201 wfaceb(1,2) = fourth*(w(1,nc3)+w(1,nc4)+w(1,nc7)+w(1,nc8))
202 wfaceb(2,2) = fourth*(w(2,nc3)+w(2,nc4)+w(2,nc7)+w(2,nc8))
203 wfaceb(3,2) = fourth*(w(3,nc3)+w(3,nc4)+w(3,nc7)+w(3,nc8))
204
205
206 wfaceb(1,3) = fourth*(w(1,nc5)+w(1,nc6)+w(1,nc7)+w(1,nc8))
207 wfaceb(2,3) = fourth*(w(2,nc5)+w(2,nc6)+w(2,nc7)+w(2,nc8))
208 wfaceb(3,3) = fourth*(w(3,nc5)+w(3,nc6)+w(3,nc7)+w(3,nc8))
209
210 wfaceb(1,4) = fourth*(w(1,nc1)+w(1,nc2)+w(1,nc5)+w(1,nc6))
211 wfaceb(2,4) = fourth*(w(2,nc1)+w(2,nc2)+w(2,nc5)+w(2,nc6))
212 wfaceb(3,4) = fourth*(w(3,nc1)+w(3,nc2)+w(3,nc5)+w(3,nc6))
213
214 wfaceb(1,5) = fourth*(w(1,nc2)+w(1,nc3)+w(1,nc6)+w(1,nc7))
215 wfaceb(2,5) = fourth*(w(2,nc2)+w(2,nc3)+w(2,nc6)+w(2,nc7))
216 wfaceb(3,5) = fourth*(w(3,nc2)+w(3,nc3)+w(3,nc6)+w(3,nc7))
217
218 wfaceb(1,6) = fourth*(w(1,nc1)+w(1,nc4)+w(1,nc5)+w(1,nc8))
219 wfaceb(2,6) = fourth*(w(2,nc1)+w(2,nc4)+w(2,nc5)+w(2,nc8))
220 wfaceb(3,6) = fourth*(w(3,nc1)+w(3,nc4)+w(3,nc5)+w(3,nc8))
221
222 IF(nbcut==0)THEN
223 brick_list(nin,ib)%POLY(1)%FACE(1)%W(1:3) = wfaceb(1:3,1)
224 brick_list(nin,ib)%POLY(1)%FACE(2)%W(1:3) = wfaceb(1:3,2)
225 brick_list(nin,ib)%POLY(1)%FACE(3)%W(1:3) = wfaceb(1:3,3)
226 brick_list(nin,ib)%POLY(1)%FACE(4)%W(1:3) = wfaceb(1:3,4)
227 brick_list(nin,ib)%POLY(1)%FACE(5)%W(1:3) = wfaceb(1:3,5)
228 brick_list(nin,ib)%POLY(1)%FACE(6)%W(1:3) = wfaceb(1:3,6)
229 ENDIF
230
231 !======================================================!
232 !COMPUTES LOCAL(VALEL) AND ADJACENT (VALVOIS) MOMENTUMS!
233 !======================================================!
234 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
235 icell = icell +1
236 IF (icell>ncell .AND. ncell/=0)icell=9
237 nbcut = brick_list(nin,ib)%NBCUT
238 jm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
239 iem = brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
240 ibm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
241 idm = brick_list(nin,ibm)%ID
242 ngm = brick_list(nin,ibm)%NG
243 gbuf =>elbuf_tab(ngm)%GBUF
244 valel(1) = alefvm_buffer%FCELL(1,iem) !rho.ux
245 valel(2) = alefvm_buffer%FCELL(2,iem) !rho.uy
246 valel(3) = alefvm_buffer%FCELL(3,iem) !rho.uz
247 valel(4) = alefvm_buffer%FCELL(4,iem) !rho
248 valel(5) = alefvm_buffer%FCELL(5,iem) !ssp
249 valel(6) = alefvm_buffer%FCELL(6,iem) !P
250 sr1 = sqrt(valel(4))
251 DO j=1,nv46
252 nadj = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
253 iev = brick_list(nin,ib)%Adjacent_Brick(j,1)
254 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
255 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
256 cellflux(j,icell,ib,1:5) = zero
257 DO iadj=1,nadj
258 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
259 IF(icellv==0)THEN
260 valvois(j,1:3,iadj) = -valel(1:3)
261 valvois(j,4,iadj) = valel(4)
262 ELSE
263 !IF INSIDE CUT CELL
264 IF(ibv/=0)THEN
265 ievm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
266 ibvm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
267 ngvm = brick_list(nin,ibvm)%NG
268 idvm = brick_list(nin,ibvm)%IDLOC
269 valvois(j,1,iadj) = alefvm_buffer%FCELL(1,ievm)
270 valvois(j,2,iadj) = alefvm_buffer%FCELL(2,ievm)
271 valvois(j,3,iadj) = alefvm_buffer%FCELL(3,ievm)
272 valvois(j,4,iadj) = alefvm_buffer%FCELL(4,ievm)
273 valvois(j,5,iadj) = alefvm_buffer%FCELL(5,ievm)
274 valvois(j,6,iadj) = alefvm_buffer%FCELL(6,ievm)
275 ELSE
276 !NOT IN CUT CELL (FRONTIER TO USUAL DOMAIN)
277 valvois(j,1,iadj) = alefvm_buffer%FCELL(1,iev)
278 valvois(j,2,iadj) = alefvm_buffer%FCELL(2,iev)
279 valvois(j,3,iadj) = alefvm_buffer%FCELL(3,iev)
280 valvois(j,4,iadj) = alefvm_buffer%FCELL(4,iev)
281 valvois(j,5,iadj) = alefvm_buffer%FCELL(5,iev)
282 valvois(j,6,iadj) = alefvm_buffer%FCELL(6,iev)
283 ENDIF
284 ENDIF
285 sr2 = sqrt(valvois(j,4,iadj))
286 !======================================================!
287 ! FLUXES CALCULATION ON EACH FACE(J) AND EACH ADJ CELL !
288 !======================================================!
289 imat = ixs(1,ie)
290 ialefvm_flg = ipm(251,imat)
291 IF(alefvm_param%IFORM/=0)ialefvm_flg = alefvm_param%IFORM !for debug purpose if set in alefvm_init.F
292
293 vf = zero
294 SELECT CASE(alefvm_param%ISOLVER)
295 CASE(1)
296 !CENTERED FVM - U average
297 vf(1) = half * (valel(1)/valel(4)+valvois(j,1,iadj)/valvois(j,4,iadj))
298 vf(2) = half * (valel(2)/valel(4)+valvois(j,2,iadj)/valvois(j,4,iadj))
299 vf(3) = half * (valel(3)/valel(4)+valvois(j,3,iadj)/valvois(j,4,iadj))
300 vf(1) = vf(1)
301 vf(2) = vf(2)
302 vf(3) = vf(3)
303 CASE(2)
304 !CENTERED FVM - rho.U average
305 vf(1) = (valel(1)+valvois(j,1,iadj))/(valel(4)+valvois(j,4,iadj))
306 vf(2) = (valel(2)+valvois(j,2,iadj))/(valel(4)+valvois(j,4,iadj))
307 vf(3) = (valel(3)+valvois(j,3,iadj))/(valel(4)+valvois(j,4,iadj))
308 vf(1) = vf(1)
309 vf(2) = vf(2)
310 vf(3) = vf(3)
311 CASE(3)
312 !CENTERED FVM - Roe-average
313 vf(1) = (valel(1)/sr1+valvois(j,1,iadj)/sr2)/(sr1+sr2)
314 vf(2) = (valel(2)/sr1+valvois(j,2,iadj)/sr2)/(sr1+sr2)
315 vf(3) = (valel(3)/sr1+valvois(j,3,iadj)/sr2)/(sr1+sr2)
316 vf(1) = vf(1)
317 vf(2) = vf(2)
318 vf(3) = vf(3)
319 CASE(4)
320 IF(dt1==zero)THEN
321 vf(1) = (valel(1) +valvois(j,1,iadj) )
322 . /(valel(4) +valvois(j,4,iadj) )
323 vf(2) = (valel(2) +valvois(j,2,iadj) )
324 . /(valel(4) +valvois(j,4,iadj) )
325 vf(3) = (valel(3) +valvois(j,3,iadj) )
326 . /(valel(4) +valvois(j,4,iadj) )
327 ELSE
328
329 !CENTERED FVM - rho.c.U average
330 vf(1) = (valel(1)*valel(5)+valvois(j,1,iadj)*valvois(j,5,iadj))
331 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
332 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois(j,5,iadj))
333 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
334 vf(3) = (valel(3)*valel(5)+valvois(j,3,iadj)*valvois(j,5,iadj))
335 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
336
337 ENDIF
338 vf(1) = vf(1)
339 vf(2) = vf(2)
340 vf(3) = vf(3)
341
342 CASE(5)
343 !Godunoc/Riemann Acoustic
344 IF(dt1==zero)THEN
345 vf(1) = (valel(1) +valvois(j,1,iadj) )
346 . /(valel(4) +valvois(j,4,iadj) )
347 vf(2) = (valel(2) +valvois(j,2,iadj) )
348 . /(valel(4) +valvois(j,4,iadj) )
349 vf(3) = (valel(3) +valvois(j,3,iadj) )
350 . /(valel(4) +valvois(j,4,iadj) )
351 ELSE
352 term2 = ( valel(6)-valvois(j,6,iadj) )/ (valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
353 vf(1) = (valel(1)*valel(5)+valvois(j,1,iadj)*valvois(j,5,iadj))
354 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *norm(1,j)
355 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois(j,5,iadj))
356 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *norm(2,j)
357 vf(3) = (valel(3)*valel(5)+valvois(j,3,iadj)*valvois(j,5,iadj))
358 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *norm(3,j)
359 ENDIF
360 vf(1) = vf(1)
361 vf(2) = vf(2)
362 vf(3) = vf(3)
363
364 CASE(6)
365 z(1:3) = brick_list(nin,ibm)%SCellCenter(1:3)
366 IF(ibv /=0)THEN
367 zadj(1:3) = brick_list(nin,ibvm)%ScellCenter(1:3)
368 ELSE
369 nc(1) = ixs(2,iev)
370 nc(2) = ixs(3,iev)
371 nc(3) = ixs(4,iev)
372 nc(4) = ixs(5,iev)
373 nc(5) = ixs(6,iev)
374 nc(6) = ixs(7,iev)
375 nc(7) = ixs(8,iev)
376 nc(8) = ixs(9,iev)
377 zadj(1) = one_over_8*sum(x(1,nc(1:8)))
378 zadj(2) = one_over_8*sum(x(2,nc(1:8)))
379 zadj(3) = one_over_8*sum(x(3,nc(1:8)))
380 ENDIF
381
382 cf(1:3) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
383 IF(ibv /= 0)THEN
384 face = brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
385 facev = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
386 IF(facev<face) cf(1:3) = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Center(1:3)
387 ENDIF
388 zzadj(1) = zadj(1)-z(1)
389 zzadj(2) = zadj(2)-z(2)
390 zzadj(3) = zadj(3)-z(3)
391 zcf(1) = cf(1) - z(1)
392 zcf(2) = cf(2) - z(2)
393 zcf(3) = cf(3) - z(3)
394 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj(3)
395 zzadj_ = zzadj(1)**2 + zzadj(2)**2 + zzadj(3)**2
396 lambda = ps / max(em20,zzadj_)
397C IF(lambda<ZERO .OR. lambda >ONE)then
398C print *, "labmda=", lambda
399C pause
400C endif
401 lambda = min(max(zero,lambda) , one)
402 lambda = sin(half*3.14159265358979d00*lambda)
403 lambda = lambda * lambda
404 !face density
405 sr1 = valel(4)
406 sr2 = valvois(j,4,iadj)
407 srf = sr1 + lambda*(sr2-sr1)
408 !face momentum density
409 vf(1) = valel(1) + lambda*(valvois(j,1,iadj)-valel(1))
410 vf(2) = valel(2) + lambda*(valvois(j,2,iadj)-valel(2))
411 vf(3) = valel(3) + lambda*(valvois(j,3,iadj)-valel(3))
412 !face velocity
413 vf(1) = vf(1) / srf
414 vf(2) = vf(2) / srf
415 vf(3) = vf(3) / srf
416 !CASE(2) gives
417 !VF(1) = (VALEL(1)+VALVOIS(J,1,IADj))/(VALEL(4)+VALVOIS(J,4,IADj))
418 !VF(2) = (VALEL(2)+VALVOIS(J,2,IADj))/(VALEL(4)+VALVOIS(J,4,IADj))
419 !VF(3) = (VALEL(3)+VALVOIS(J,3,IADj))/(VALEL(4)+VALVOIS(J,4,IADj))
420 vf(1) = vf(1)
421 vf(2) = vf(2)
422 vf(3) = vf(3)
423 END SELECT
424
425 cellflux(j,icell,ib,iadj) = (vf(1)*norm(1,j) + vf(2)*norm(2,j) + vf(3)*norm(3,j))
426
427 wface(1) = brick_list(nin,ib)%POLY(icell)%FACE(j)%W(1)
428 wface(2) = brick_list(nin,ib)%POLY(icell)%FACE(j)%W(2)
429 wface(3) = brick_list(nin,ib)%POLY(icell)%FACE(j)%W(3)
430
431 IF(icellv /=0)THEN
432 wfacev(1) = brick_list(nin,ib)%POLY(icellv)%FACE(jv)%W(1)
433 wfacev(2) = brick_list(nin,ib)%POLY(icellv)%FACE(jv)%W(2)
434 wfacev(3) = brick_list(nin,ib)%POLY(icellv)%FACE(jv)%W(3)
435 ELSE
436 wfacev(1) = wface(1)
437 wfacev(2) = wface(2)
438 wfacev(3) = wface(3)
439 ENDIF
440
441 wface(1) = half*(wface(1)+wfacev(1))
442 wface(2) = half*(wface(2)+wfacev(2))
443 wface(3) = half*(wface(3)+wfacev(3))
444
445
446 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(1) = vf(1)-wface(1)
447 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(2) = vf(2)-wface(2)
448 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(3) = vf(3)-wface(3)
449
450
451 !======================================================!
452 ! REDUCTION FACTOR FOR POLYHEDRA FACES !
453 !======================================================!
454 !facette non communicante
455 numnod_ = brick_list(nin,ib)%POLY(icell)%NumNOD
456 k = 0
457 DO l=1,numnod_
458 inod = brick_list(nin,ib)%POLY(icell)%ListNodID(l)
459 IF(int22_buf%IsNodeOnFace(inod,j))k = k +1
460 ENDDO
461 face = brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
462 IF(ibv==0)THEN
463 !calcul de facette
464 facev = face
465 numnod_v = 8
466 kv = 1
467 ELSE
468 facev = brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
469 numnod_ = brick_list(nin,ibv)%POLY(icellv)%NumNOD
470 kv = 0
471 DO l=1,numnod_
472 inod = brick_list(nin,ibv)%POLY(icellv)%ListNodID(l)
473 IF(int22_buf%IsNodeOnFace(inod,jv))kv= kv +1
474 ENDDO
475 ENDIF
476 face = min(face,facev)
477 If(k==0 .OR. kv==0) face = zero !facette non communicante si pas de noeud de brique.
478 IF(ibv/=0 .AND. ibm==ibvm) face=zero !blocage flux entre le main et son escalve
479 cellflux(j,icell,ib,iadj) = face * cellflux(j,icell,ib,iadj)
480 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_FLUX(iadj) = cellflux(j,icell,ib,iadj)
481 enddo!next IADj
482 enddo!next J
483 enddo!next ICELL
484 enddo!next IB
485
486 !======================================================!
487 ! debug output !
488 !======================================================!
489 if(debug_outp)then
490 if(ibug22_flux22>0 .OR. ibug22_flux22==-1)then
491 print *, " |------e22flux3_int22_fvm.F------|"
492 do ib=nbf,nbl
493 ie = brick_list(nin,ib)%Id
494 if (ibug22_flux22>0 .and. ibug22_flux22/= ixs(11,ie))cycle
495 icell = 0
496 ncell = brick_list(nin,ib)%nbcut
497 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
498 icell = icell +1
499 IF (icell>ncell .AND. ncell/=0)icell=9
500
501 print *, " brique =", ixs(11,ie)
502 print *, " NCYCLE =", ncycle
503 print *, " icell =", icell
504 print *, " mcell =", mcell
505 print *, " nbcut =", nbcut
506 do i=1,6
507 nadj = brick_list(nin,ib)%POLY(icell)%FACE(i)%NAdjCell
508 write (*,fmt='(A,I10,A,5E26.14)') " phi(1:5,",i,")=", cellflux(i,icell,ib,1:nadj)
509 enddo
510 enddo! next ICELL
511 enddo!next IB
512 print *, " ------------------------"
513 endif
514 endif
515
516 !======================================================!
517
518
519 !==============!
520 CALL my_barrier
521 !==============!
522
523 !======================================================!
524 ! STEP B : NON CONFORM MESH !
525 ! USE CONSISTENT FLUX !
526 !======================================================!
527 DO ib=nbf,nbl
528 ie = brick_list(nin,ib)%ID
529 mlw = brick_list(nin,ib)%MLW
530 ncell = brick_list(nin,ib)%NBCUT
531 mcell = brick_list(nin,ib)%MainID
532 icell = 0
533 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
534 icell = icell +1
535 IF (icell>ncell .AND. ncell/=0)icell=9
536 DO j=1,6
537 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (1) = zero
538 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (2) = zero
539 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (3) = zero
540 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (4) = zero
541 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (5) = zero
542 enddo!next J
543
544 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 = zero
545 !======================================================!
546 ! MONOMATERIAL UPWIND TREATMENT !
547 !======================================================!
548 ie_m = brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
549 mat = ixs(1,ie_m)
550 upwl(1:6) = pm(16,mat)
551 reduc = pm(92,mat)
552 ddvol = zero
553 DO j=1,6
554 nadj = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
555 iclose = brick_list(nin,ib)%ClosedSurf(j) !pour l'hypothese de surface fermee : double coupe sur la facette et pas de partitionnement de l'autre cote => ON FERME.
556 IF(iclose==1) nadj=0
557 DO iadj = 1,nadj
558 idv = brick_list(nin,ib)%Adjacent_Brick(j,1)
559 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
560 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
561 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
562 IF(idv==0)THEN
563 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*reduc
564 ELSEIF(idv>0)THEN
565 ng = brick_list(nin,ib)%NG
566 isilent = iparg(64,ng)
567 IF(isilent==1)THEN
568 upwl(j)=one
569 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*pm(92,ixs(1,idv))
570 ENDIF
571 ENDIF
572 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX(iadj) =
573 . cellflux(j,icell,ib,iadj)-upwl(j)*abs(cellflux(j,icell,ib,iadj))
574 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 =
575 . brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 + cellflux(j,icell,ib,iadj)+upwl(j)*abs(cellflux(j,icell,ib,iadj))
576 !-----DDVOL
577 ddvol = ddvol + cellflux(j,icell,ib,iadj)
578 !-----DDVOL*DT IS SUM FOR INCOMING AND OUTCOMING VOLUMES
579 enddo!next IADJ
580 enddo!next J
581
582 brick_list(nin,ib)%POLY(icell)%DDVOL = ddvol
583
584 !======================================================!
585 enddo!next ICELL
586 enddo!next IB
587
588 !==============!
589 CALL my_barrier
590 !==============!
591
592 !---------------------------------------------------------!
593 ! SECND CELLS STACK !
594 !---------------------------------------------------------!
595 !STACK Secnd cells values from ones connected to current main cell
596 IF(int22>0)THEN
597 nin = 1
598 DO ib=nbf,nbl
599 mlw = brick_list(nin,ib)%MLW
600 num = brick_list(nin,ib)%SecndList%Num
601 mcell = brick_list(nin,ib)%mainID
602 ddvol = zero
603 DO k=1,num
604 ibv = brick_list(nin,ib)%SecndList%IBV(k)
605 icellv = brick_list(nin,ib)%SecndList%ICELLv(k)
606 ddvol = ddvol + brick_list(nin,ibv)%POLY(icellv)%DDVOL
607 ENDDO
608 ddvol = ddvol + brick_list(nin,ib)%POLY(mcell)%DDVOL
609 !updating ddvol cumulative stack
610 brick_list(nin,ib)%POLY(mcell)%DDVOL = ddvol
611 enddo!next IB
612 ENDIF
613
614 !==============!
615 CALL my_barrier
616 !==============!
617
618 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22bufbric_mod::int22_buf
type(int22_buf_t) int22_buf
Definition cut-cell-search_mod.F:130
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
i22tri_mod::ibug22_flux22
integer ibug22_flux22
Definition cut-cell-search_mod.F:296
my_barrier
subroutine my_barrier
Definition machine.F:31