37
38
39
40 USE elbufdef_mod
44 use element_mod , only : nixs
45
46
47
48
49
50
51
52
53
54
55
56
57
58#include "implicit_f.inc"
59
60
61
62#include "spmd_c.inc"
63#include "com01_c.inc"
64#include "com04_c.inc"
65#include "com08_c.inc"
66#include "param_c.inc"
67#include "task_c.inc"
68
69
70
71 INTEGER IXS(NIXS,*),NV46,ITRIMAT,ITASK
73 . flux(nv46,*),vfrac(*)
74 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
75
76
77
78 INTEGER K,KK,J1,J2,J
80 . vol0,av0,uav0,alphi,ualphi,aaa,ff(nv46,5),udt,phi0
81 INTEGER :: IE, MLW, IADJv, NADJv, IB, NBF, NBL, ICELL,ICELLM, MCELL, , IBM,NG,IDLOC,NADJ,IADJ
82 INTEGER :: NIN,NCELL,IBV,IFV,ICELLv, IEV
83 my_real :: volg, alph, alphv(6,5), tmpflux(nv46,5)
85 LOGICAL :: debug_outp
86
87
88
89 IF(trimat==0)RETURN
90
91
92
93
94 IF(dt1>zero)THEN
95 udt = one/dt1
96 ELSE
97 udt = zero
98 ENDIF
99
100 nin = 1
101 nbf = 1+itask*
nb/nthread
102 nbl = (itask+1)*
nb/nthread
104
105
106
107 debug_outp = .false.
109 debug_outp = .false.
111 do ib=nbf,nbl
115 if(mlw==51)then
116 debug_outp=.true.
117 endif
118 endif
119 enddo
121 debug_outp = .true.
122 kk = 1
123 do ib=nbf,nbl
125 if(mlw/=51)then
126 kk = 0
127 endif
128 enddo
129 if (kk==0)debug_outp=.false.
130 endif
132 endif
133 if(debug_outp)then
134 print *, " |----ale51_antidiff3_int22.F-----|"
135 print *, " | THREAD INFORMATION |"
136 print *, " |--------------------------------|"
137 print *, " NCYCLE =", ncycle
138 print *, " ITRIMAT =", itrimat
139 endif
140
141
142
143 DO ib=nbf,nbl
148 icell = 0
149 DO WHILE (icell<=ncell)
150 icell = icell +1
151 IF (icell>ncell .AND. ncell/=0)icell=9
152
153 j =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
154 icellm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(2)
155 IF(j==0)THEN
156 ie_m = ie
157 ibm = ib
158 icellm = mcell
159 ELSEIF(j<=nv46)THEN
162 ELSE
163 j1 = j/10
164 j2 = mod(j,10)
165 ibv =
brick_list(nin,ib )%Adjacent_Brick(j1,4)
166 ie_m =
brick_list(nin,ibv)%Adjacent_Brick(j2,1)
167 ibm =
brick_list(nin,ibv)%Adjacent_Brick(j2,4)
168 ENDIF
172 IF(mlw/=51)cycle
173 alph =
brick_list(nin,ibm)%POLY(icellm)%VFRACm(itrimat)
174 volg = elbuf_tab(ng)%GBUF%VOL(idloc)
175 vol0 = volg*udt
176 av0 = alph * vol0
177 uav0 = vol0 - av0
178 alphi = zero
179 ualphi = zero
180 phi0 = zero
181
182
183
184
185 DO k=1,nv46
186 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
187 DO iadj=1,nadj
188 tmpflux(k,iadj) =
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_UpwFLUX(iadj)
189 IF(tmpflux(k,iadj)>zero)THEN
193 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_Cell(iadj)
194 IF(icellv==0) THEN
195 alphv(k,iadj) = alph
196 ELSE
197 IF(ibv==0)THEN
198 IF(iev==0)print *, "inter22 : potential material leakage, Check domain boundaries..."
199 alphv(k,iadj) = vfrac(iev)
200 ELSE
201 alphv(k,iadj) =
brick_list(nin,ibv)%POLY(icellv)%VFRACm(itrimat)
202 ENDIF
203 ENDIF
204 ff(k,iadj)= alphv(k,iadj) * tmpflux(k,iadj)
205
206 alphi = alphi + ff(k,iadj)
207
208 phi0 = phi0 + tmpflux(k,iadj)
209 ENDIF
210 enddo
211 enddo
212
213 ualphi = phi0 - alphi
214
215
216
217 IF(alphi>av0.AND.av0>zero)THEN
218
219
220
221 aaa = av0 / alphi
222 DO k=1,nv46
223 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
224 DO iadj=1,nadj
225 IF(tmpflux(k,iadj)>zero)THEN
226 ff(k,iadj) = ff(k,iadj) * aaa
227 ENDIF
228 enddo
229 enddo
230 ELSEIF(ualphi>uav0.AND.uav0>zero)THEN
231
232
233
234 aaa = uav0/ualphi
235 DO k=1,nv46
236 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
237 DO iadj=1,nadj
238 IF(tmpflux(k,iadj)>zero)THEN
239 ff(k,iadj) = tmpflux(k,iadj) + (ff(k,iadj)-tmpflux(k,iadj))*aaa
240 ENDIF
241 enddo
242 enddo
243 ENDIF
244
245
246
247 DO k=1,nv46
248 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(k)%NAdjCell
249 DO iadj=1,nadj
253 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_Cell(iadj)
254 IF(tmpflux(k,iadj)>zero)THEN
255 ff(k,iadj) = half * ( ff(k,iadj)*(one-
ale%UPWIND%UPWSM)+alph*tmpflux(k,iadj)*(one+
ale%UPWIND%UPWSM) )
256
257
258
259
260 if(debug_outp)then
263
264 print *, " brique =", ixs(11,ie)
265 print *, " icell =", icell
266 print *, " FACE =", k
267 print *, " ALPH =", alph
268 print *, " ALPHv =", alphv(k,iadj)
269 write (*,fmt=
'(A,6E26.14)')
" WAS Flux(J) =",
brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_upwFLUX(iadj)
270 write (*,fmt='(A,6E26.14)')" IS Flux(J) =", ff(k,iadj)
271 print *, " ------------------------"
272 endif
273 endif
274
275
276
277 brick_list(nin,ib)%POLY(icell)%FACE(k)%Adjacent_UpwFLUX(iadj) = ff(k,iadj)
278
279
280 IF(icellv>0)THEN
281
282
283 IF(ibv>0)THEN
284
285 nadjv =
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%NAdjCell
286 DO iadjv=1,nadjv
287 IF(
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_Cell(iadjv)==icell)
EXIT
288
289 ENDDO
290 debug_tmp =
brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_UpwFLUX(iadjv)
291 brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%Adjacent_UpwFLUX(iadjv) = -ff(k,iadj)
292 ELSE
293
294 debug_tmp = flux(ifv,iev)
295 flux(ifv,iev) = -ff(k,iadj)
296 ENDIF
297
298
299 if(debug_outp)then
301 print *, " => Setting adjacent flux consequently :"
302 print *, " brique.V =", ixs(11,iev)
303 print *, " icell.V =", icellv
304 print *, " FACE.V =", ifv
305 write (*,fmt='(A,6E26.14)')
306 . " WAS Flux(J) =", debug_tmp
307 write (*,fmt='(A,6E26.14)')
308 . " IS Flux(J) =", -ff(k,iadj)
309 print *, " ---"
310 endif
311 endif
312
313 ELSE
314
315 ENDIF
316 ENDIF
317 enddo
318 enddo
319
320 enddo
321 enddo
322
323
324 RETURN
type(brick_entity), dimension(:,:), allocatable, target brick_list