39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
67 USE elbufdef_mod
69
70
71
72#include "implicit_f.inc"
73
74
75
76#include "mvsiz_p.inc"
77
78
79
80#include "param_c.inc"
81#include "inter22.inc"
82#include "task_c.inc"
83#include "com01_c.inc"
84#include "com08_c.inc"
85
86
87
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
92
93
94 INTEGER MAT, MLW, NC(8), I, IE,J, K, Kv,L,ITASK, IDm, IDV
95 INTEGER IDvm, IBvm, NGvm, IEvm, IEV, Jm, IMAT, IALEFVM_FLG
96 INTEGER IB,IBv,NIN, NBCUT,ICELL,NCELL,NGm
99 INTEGER :: NBF,NBL, MCELL,iNOD,ICELLv,numnod, numnodV
100 TYPE(G_BUFEL_) ,POINTER :: GBUF
101 my_real :: face , z(3), zadj(3), zzadj_, cf(3), zcf(3),zzadj(3)
103 INTEGER :: NUM, NADJ, IADJ, JV, NG, IE_M, IEm, IBm
104 my_real :: facev, ddvol, valel(6), valvois(6,6,5), sr1, sr2
105 LOGICAL debug_outp
106
107
108
109 IF(int22 == 0)RETURN
110
111
112
113
114
115
116
117 ibvm = 0
118 vf = zero
119 nin = 1
120 nbf = 1+itask*
nb/nthread
121 nbl = (itask+1)*
nb/nthread
123
124
125
126
127
128 debug_outp = .false.
130 debug_outp = .false.
132 do ib=nbf,nbl
135 enddo
137 debug_outp = .true.
138 endif
139 endif
140 if(debug_outp .AND. itask==0)then
141 print *, " |---------eflux3_int22_fvm.F---------|"
142 print *, " | THREAD INFORMATION |"
143 print *, " |------------------------------------|"
144 print *, " NCYCLE =", ncycle
145 endif
146
147
148
149
150
151
152
153
154 DO ib=nbf,nbl
159 icell = 0
160
161
162
163
164 flux(1:6,ie) = zero
165 flu1(ie) = zero
166
167
168
169
170
171
172
176 DO j=1,nv46
177 lnorm(j) = sqrt(
norm(1,j)**2 +
norm(2,j)**2 +
norm(3,j)**2 )
178 norm(1:3,j) =
norm(1:3,j) / lnorm(j)
179 ENDDO
180
181
182
183
184 DO WHILE (icell<=ncell)
185 icell = icell +1
186 IF (icell>ncell .AND. ncell/=0)icell=9
188 jm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
189 iem =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
190 ibm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
193 gbuf =>elbuf_tab(ngm)%GBUF
200 sr1 = sqrt(valel(4))
201 DO j=1,nv46
202 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
206 cellflux(j,icell,ib,1:5) = zero
207 DO iadj=1,nadj
208 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
209 IF(icellv==0)THEN
210 valvois(j,1:3,iadj) = -valel(1:3)
211 valvois(j,4,iadj) = valel(4)
212 ELSE
213
214 IF(ibv/=0)THEN
215 ievm =
brick_list(nin,ibv)%POLY(icellv)%WhereIsMain
216 ibvm =
brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
225 ELSE
226
233 ENDIF
234 ENDIF
235 sr2 = sqrt(valvois(j,4,iadj))
236
237
238
239 imat = ixs(1,ie)
240 ialefvm_flg = ipm(251,imat)
242
243
245 CASE(1)
246
247 vf(1) = half * (valel(1)/valel(4)+valvois(j,1,iadj)/valvois(j,4,iadj))
248 vf(2) = half * (valel(2)/valel(4)+valvois(j,2,iadj)/valvois(j
249 vf(3) = half * (valel(3)/valel(4)+valvois(j,3,iadj)/valvois(j,4,iadj))
250 CASE(2)
251
252 vf(1) = (valel(1)+valvois(j,1,iadj))/(valel(4)+valvois(j,4,iadj))
253 vf(2) = (valel(2)+valvois(j,2,iadj))/(valel(4)+valvois(j,4,iadj))
254 vf(3) = (valel(3)+valvois(j,3,iadj))/(valel(4)+valvois(j,4,iadj))
255 CASE(3)
256
257 vf(1) = (valel(1)/sr1+valvois(j,1,iadj)/sr2)/(sr1+sr2)
258 vf(2) = (valel(2)/sr1+valvois(j,2,iadj)/sr2)/(sr1+sr2)
259 vf(3) = (valel(3)/sr1+valvois(j,3,iadj)/sr2)/(sr1+sr2)
260 CASE(4)
261 IF(dt1==zero)THEN
262 vf(1) = (valel(1) +valvois(j,1,iadj) )
263 . /(valel(4) +valvois(j,4,iadj) )
264 vf(2) = (valel(2) +valvois(j,2,iadj) )
265 . /(valel(4) +valvois(j,4,iadj
266 vf(3) = (valel(3) +valvois(j,3,iadj) )
267 . /(valel(4) +valvois(j,4,iadj
268 ELSE
269
270
271 vf(1) = (valel(1)*valel(5)+valvois(j,1,iadj)*valvois(j,5,iadj))
272 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
273 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois(j,5,iadj))
274 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
275 vf(3) = (valel(3)*valel(5)+valvois(j,3,iadj)*valvois(j,5,iadj))
276 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
277
278 ENDIF
279
280 CASE(5)
281
282 IF(dt1==zero)THEN
283 vf(1) = (valel(1) +valvois(j,1,iadj)
284 . /(valel(4) +valvois
285 vf(2) = (valel(2) +valvois(j,2,iadj)
286 . /(valel(4) +valvois(j,4,iadj
287 vf(3) = (valel(3) +valvois(j,3,iadj) )
288 . /(valel(4) +valvois(j,4,iadj) )
289 ELSE
290 term2 = ( valel(6)-valvois(j,6,iadj) )/ (valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj))
291 vf(1) = (valel(1)*valel(5)+valvois(j,1,iadj)*valvois(j,5,iadj))
292 . /(valel(4)*valel(5)+valvois(j,4
293 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois(j,5,iadj))
294 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *
norm(2,j)
295 vf(3) = (valel(3)*valel(5)+valvois(j,3,iadj)*valvois(j,5,iadj))
296 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *
norm(3,j)
297 ENDIF
298 CASE(6)
300 IF(ibv /=0)THEN
301 zadj(1:3) =
brick_list(nin,ibvm)%ScellCenter(1:3)
302 ELSE
303 nc(1) = ixs(2,iev)
304 nc(2) = ixs(3,iev)
305 nc(3) = ixs(4,iev)
306 nc(4) = ixs(5,iev)
307 nc(5) = ixs(6,iev)
308 nc(6) = ixs(7,iev)
309 nc(7) = ixs(8,iev)
310 nc(8) = ixs(9,iev)
311 zadj(1) = one_over_8*sum(x(1,nc(1:8
312 zadj(2) = one_over_8*sum(x(2,nc(1:8)))
313 zadj(3) = one_over_8*sum(x(3,nc(1:8)))
314 ENDIF
315
316 cf(1:3) =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
317 IF(ibv /= 0)THEN
318 face =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
319 facev =
brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
321 ENDIF
322 zzadj(1) = zadj(1)-z(1)
323 zzadj(2) = zadj(2)-z(2)
324 zzadj(3) = zadj(3)-z(3)
325 zcf(1) = cf(1) - z(1)
326 zcf(2) = cf(2) - z(2)
327 zcf(3) = cf(3) - z(3)
328 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj
329 zzadj_ = zzadj(1)**2 + zzadj(2)**2 + zzadj(3)**2
330 lambda = ps /
max(em20,zzadj_)
331 IF(lambda<zero .OR. lambda >one)then
332 print *, "labmda=", lambda
333
334 endif
335 lambda =
min(
max(zero,lambda) , one)
336 lambda = sin(half*3.14159265358979d00*lambda)
337 lambda = lambda * lambda
338
339 sr1 = valel(4)
340 sr2 = valvois(j,4,iadj)
341 srf = sr1 + lambda*(sr2-sr1)
342
343 vf(1) = valel(1) + lambda*(valvois(j,1,iadj)-valel(1))
344 vf(2) = valel(2) + lambda*(valvois(j,2,iadj)-valel(2))
345 vf(3) = valel(3) + lambda*(valvois(j,3,iadj)-valel(3))
346
347 vf(1) = vf(1) / srf
348 vf(2) = vf(2) / srf
349 vf(3) = vf(3) / srf
350
351
352
353
354 END SELECT
355
356
357 cellflux(j,icell,ib,iadj) = (vf(1)*
norm(1,j) + vf(2)*
norm(2,j) + vf(3)*
norm(3,j))
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(1) = vf(1)
375 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(2) = vf(2)
376 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(3) = vf(3)
377
378
379
380
381
382
383 numnod =
brick_list(nin,ib)%POLY(icell)%NumNOD
384 k = 0
385 DO l=1,numnod
386 inod =
brick_list(nin,ib)%POLY(icell)%ListNodID(l)
387 IF(
int22_buf%IsNodeOnFace(inod,j))k = k +1
388 ENDDO
389 face =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
390 IF(ibv==0)THEN
391
392 facev = face
393 numnodv = 8
394 kv = 1
395 ELSE
396 facev =
brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
397 numnod =
brick_list(nin,ibv)%POLY(icellv)%NumNOD
398 kv = 0
399 DO l=1,numnod
400 inod =
brick_list(nin,ibv)%POLY(icellv)%ListNodID(l)
401 IF(
int22_buf%IsNodeOnFace(inod,jv))kv= kv +1
402 ENDDO
403 ENDIF
404 face =
min(face,facev)
405 If(k==0 .OR. kv==0) face = zero
406 IF(ibv/=0 .AND. ibm==ibvm) face=zero
407 cellflux(j,icell,ib,iadj) = face * cellflux(j,icell,ib,iadj)
408
409 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_FLUX(iadj) = cellflux(j,icell,ib,iadj)
410 enddo
411 enddo
412 enddo
413 enddo
414
415
416
417
418 if(debug_outp)then
420 if(itask==0)then
422 print *, " |------e22flux3_int22_fvm.F------|"
426 icell = 0
429 nbcut = ncell
430 DO WHILE (icell<=ncell)
431 icell = icell +1
432 IF (icell>ncell .AND. ncell/=0)icell=9
433
434
435
436
437
438
439 do i=1,6
440 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(i)%NAdjCell
441
442
443
444
445
446!
write (*,fmt=
'(A,I10,A,5E26.14)')
" phi(1:5,",i,
")=",
brick_list(nin,ib)%POLY(icell)%FACE(2)%Adjacent_FLUX(1:nadj)
447
448!
write (*,fmt=
'(A,I10,A,5E26.14)')
" phi(1:5,",i,
")=",
brick_list(nin,ib)%POLY(icell)%FACE(4)%Adjacent_FLUX(1:nadj)
449
450
451
452
453
454
455!
write (*,fmt=
'(A,6E26.14)')
" N-x =",
brick_list(nin,ib)%N(1:6,1)
456
457
458 enddo
459 enddo
460 enddo
461 print *, " ------------------------"
462 endif
463 endif
464 endif
465
466
467
468
469
471
472
473
474
475
476
477 DO ib=nbf,nbl
482 icell = 0
483 DO WHILE (icell<=ncell)
484 icell = icell +1
485 IF (icell>ncell .AND. ncell/=0)icell=9
486 DO j=1,6
487 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (1) = zero
488 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (2) = zero
489 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (3) = zero
490 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (4) = zero
491 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (5) = zero
492 enddo
493
494 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 = zero
495
496
497
498 ie_m =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
499 mat = ixs(1,ie_m)
500 upwl(1:6) = pm(16,mat)
501 reduc = pm(92,mat)
502 ddvol = zero
503 DO j=1,6
504 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
506 IF(iclose==1) nadj=0
507 DO iadj = 1,nadj
511 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
512 IF(idv==0)THEN
513 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*reduc
514 ELSEIF(idv>0)THEN
516 isilent = iparg(64,ng)
517 IF(isilent==1)THEN
518 upwl(j)=one
519 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*pm(92,ixs(1,idv))
520 ENDIF
521 ENDIF
522 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX(iadj) =
523 . cellflux(j,icell,ib,iadj)-upwl(j)*abs(cellflux(j,icell,ib,iadj))
524 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 =
525 .
brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 + cellflux(j,icell,ib,iadj)+upwl(j)*abs(cellflux(j,icell,ib,iadj))
526
527 ddvol = ddvol + cellflux(j,icell,ib,iadj)
528
529 enddo
530 enddo
531
533
534
535 enddo
536 enddo
537
538
540
541
542
543
544
545
546 IF(int22>0)THEN
547 nin = 1
548 DO ib=nbf,nbl
552 ddvol = zero
553 DO k=1,num
555 icellv =
brick_list(nin,ib)%SecndList%ICELLv(k)
556 ddvol = ddvol +
brick_list(nin,ibv)%POLY(icellv)%DDVOL
557 ENDDO
558 ddvol = ddvol +
brick_list(nin,ib)%POLY(mcell)%DDVOL
559
561 enddo
562 ENDIF
563
564
566
567
568 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
type(alefvm_buffer_), target alefvm_buffer
type(alefvm_param_), target alefvm_param
type(brick_entity), dimension(:,:), allocatable, target brick_list
type(int22_buf_t) int22_buf