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 36 of file aflux3_int22_fvm.F.

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