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
65
68 USE elbufdef_mod
70 use element_mod , only : nixs
71
72
73
74#include "implicit_f.inc"
75
76
77
78#include "mvsiz_p.inc"
79
80
81
82#include "param_c.inc"
83#include "inter22.inc"
84#include "task_c.inc"
85#include "com01_c.inc"
86#include "com08_c.inc"
87
88
89
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
94
95
96 INTEGER MAT, MLW, NC(8), I, IE,J, K, Kv,L,ITASK, IDm, IDV
97 INTEGER , IBvm, NGvm, IEvm, IEV, Jm, IMAT, IALEFVM_FLG
98 INTEGER IB,IBv,NIN, NBCUT,ICELL,NCELL,NGm
101 INTEGER :: NBF,NBL, MCELL,iNOD,ICELLv,numnod, numnodV
102 TYPE(G_BUFEL_) ,POINTER :: GBUF
103 my_real :: face , z(3), zadj(3), zzadj_, cf(3), zcf(3),zzadj(3)
105 INTEGER :: NUM, NADJ, IADJ, JV, NG, IE_M, IEm, IBm
106 my_real :: facev, ddvol, valel(6), valvois(6,6,5), sr1, sr2, srf
107 LOGICAL debug_outp
108
109
110
111 IF(int22 == 0)RETURN
112
113
114
115
116
117
118
119 ibvm = 0
120 vf = zero
121 nin = 1
122 nbf = 1+itask*
nb/nthread
123 nbl = (itask+1)*
nb/nthread
125
126
127
128
129
130 debug_outp = .false.
132 debug_outp = .false.
134 do ib=nbf,nbl
137 enddo
139 debug_outp = .true.
140 endif
141 endif
142 if(debug_outp .AND. itask==0)then
143 print *, " |---------eflux3_int22_fvm.F---------|"
144 print *, " | THREAD INFORMATION |"
145 print *, " |------------------------------------|"
146 print *, " NCYCLE =", ncycle
147 endif
148
149
150
151
152
153
154 ! a.1 :
main & secnd cells fluxes
155
156 DO ib=nbf,nbl
161 icell = 0
162
163
164
165
166 flux(1:6,ie) = zero
167 flu1(ie) = zero
168
169
170
171
172
173
174
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
185
186 DO WHILE (icell<=ncell)
187 icell = icell +1
188 IF (icell>ncell .AND. ncell/=0)icell=9
190 jm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
191 iem =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
192 ibm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
195 gbuf =>elbuf_tab(ngm)%GBUF
198
202 sr1 = sqrt(valel(4))
203 DO j=1,nv46
204 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
208 cellflux(j,icell,ib,1:5) = zero
209 DO iadj=1,nadj
210 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
211 IF(icellv==0)THEN
212 valvois(j,1:3,iadj) = -valel(1:3)
213 valvois(j,4,iadj) = valel(4)
214 ELSE
215
216 IF(ibv/=0)THEN
217 ievm =
brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
218 ibvm =
brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
227 ELSE
228
235 ENDIF
236 ENDIF
237 sr2 = sqrt(valvois(j,4,iadj))
238
239
240
241 imat = ixs(1,ie)
242 ialefvm_flg = ipm(251,imat)
244
245
247 CASE(1)
248
249 vf(1) = half * (valel(1)/valel(4)+valvois(j,1,iadj
250 vf(2) = half * (valel
251 vf(3) = half * (valel(3)/valel(4)+valvois(j,3,iadj
252 CASE(2)
253
254 vf(1) = (valel(1)+valvois(j,1,iadj
255 vf(2) = (valel(2)+valvois(j,2,iadj
256 vf(3) = (valel(3)+valvois
257 CASE(3)
258
259 vf(1) = (valel(1)/sr1
260 vf(2) = (valel(2)/sr1+valvois(j,2,iadj)/sr2)/(sr1+sr2)
261 vf(3) = (valel(3)/sr1+valvois(j
262 CASE(4)
263 IF(dt1==zero)THEN
264 vf(1) = (valel(1) +valvois(j
265 . /(valel(4)
266 vf(2) = (valel(2) +valvois
267 . /(valel(4) +valvois(j,4,iadj) )
268 vf(3)
269 . /(valel(4) +valvois(j,4,iadj)
270 ELSE
271
272
273 vf(1) = (valel(1)*valel(5)+valvois
274 . /(valel(4)
275 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois
276 . /(valel(4)*valel
277 vf(3) = (valel(3)*valel(5)
278 . /(valel(4)*valel
279
280 ENDIF
281
282 CASE(5)
283
284 IF(dt1==zeroTHEN
285 vf(1) = (valel(1) +valvois
286 . /(valel(4) +valvois(j,4,iadj) )
287 vf(2) = (valel(2) +valvois(j,2,iadj
288 .
289
290 . /(valel(4) +valvois
291 ELSE
292 term2 = ( valel(6)-valvois(j,6,iadj
293 vf(1) = (valel(1)*valel(5)
294 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *
norm(1,j)
295 vf(2) = (valel(2)*valel(5)+valvois(j,2,iadj)*valvois(j,5,iadj))
296 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *
norm(2,j)
297 vf(3) = (valel(3)*valel(5)+valvois(j,3,iadj)*valvois(j,5,iadj))
298 . /(valel(4)*valel(5)+valvois(j,4,iadj)*valvois(j,5,iadj)) + term2 *
norm(3,j)
299 ENDIF
300 CASE(6)
302 IF(ibv /=0)THEN
303 zadj(1:3) =
brick_list(nin,ibvm)%ScellCenter(1:3)
304 ELSE
305 nc(1) = ixs(2,iev)
306 nc(2) = ixs(3,iev)
307 nc(3) = ixs(4,iev)
308 nc(4) = ixs(5,iev)
309 nc(5) = ixs(6,iev)
310 nc(6) = ixs(7,iev)
311 nc(7) = ixs(8,iev)
312 nc(8) = ixs(9,iev)
313 zadj(1) = one_over_8*sum(x(1,nc(1:8)))
314 zadj(2) = one_over_8*sum(x(2,nc(1:8)))
315 zadj(3) = one_over_8*sum(x(3,nc(1:8)))
316 ENDIF
317
318 cf(1:3) =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
319 IF(ibv /= 0)THEN
320 face =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
321 facev =
brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
322 IF(facev<face) cf(1:3) =
brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Center(1:3)
323 ENDIF
324 zzadj(1) = zadj(1)-z(1)
325 zzadj(2) = zadj(2)-z(2)
326 zzadj(3) = zadj(3)-z(3)
327 zcf(1) = cf(1) - z(1)
328 zcf(2) = cf(2) - z(2)
329 zcf(3) = cf(3) - z(3)
330 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj(3)
331 zzadj_ = zzadj(1)**2 + zzadj
332 lambda = ps /
max(em20,zzadj_)
333 IF(lambda<zero .OR. lambda >one)then
334 print *, "labmda=", lambda
335
336 endif
337 lambda =
min(
max(zero,lambda) , one)
338 lambda = sin(half*3.14159265358979d00*lambda)
339 lambda = lambda * lambda
340
341 sr1 = valel(4)
342 sr2 = valvois(j
343 srf = sr1 + lambda*(sr2-sr1)
344
345 vf(1) = valel(1) + lambda*(valvois(j,1,iadj)-valel(1))
346 vf(2) = valel(2) + lambda*(valvois(j,2,iadj)-valel(2))
347 vf(3) = valel(3) + lambda*(valvois(j,3,iadj)-valel(3))
348
349 vf(1) = vf(1) / srf
350 vf(2) = vf(2) / srf
351 vf(3) = vf(3) / srf
352
353
354
355
356 END SELECT
357
358
359 cellflux(j,icell,ib,iadj) = (vf(1)*
norm(1,j) + vf(2)*
norm(2,j) + vf(3)*
norm(3,j))
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(1) = vf(1)
377 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(2) = vf(2)
378 brick_list(nin,ib)%POLY(icell)%FACE(j)%Vel(3) = vf(3)
379
380
381
382
383
384
385 numnod =
brick_list(nin,ib)%POLY(icell)%NumNOD
386 k = 0
387 DO l=1,numnod
388 inod =
brick_list(nin,ib)%POLY(icell)%ListNodID(l)
389 IF(
int22_buf%IsNodeOnFace(inod,j))k = k +1
390 ENDDO
392 IF(ibv==0)THEN
393
394 facev = face
395 numnodv = 8
396 kv = 1
397 ELSE
398 facev =
brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf
399 numnod =
brick_list(nin,ibv)%POLY(icellv)%NumNOD
400 kv = 0
401 DO l=1,numnod
402 inod =
brick_list(nin,ibv)%POLY(icellv)%ListNodID(l)
403 IF(
int22_buf%IsNodeOnFace(inod,jv))kv= kv +1
404 ENDDO
405 ENDIF
406 face =
min(face,facev)
407 If(k==0 .OR. kv==0) face = zero
408 IF(ibv/=0 .AND. ibm==ibvm) face=zero
409 cellflux(j,icell,ib,iadj) = face * cellflux(j,icell,ib,iadj)
410
411 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_FLUX(iadj) = cellflux(j,icell,ib,iadj)
412 enddo
413 enddo
414 enddo
415 enddo
416
417
418
419
420 if(debug_outp)then
422 if(itask==0)then
424 print *, " |------e22flux3_int22_fvm.F------|"
428 icell = 0
431 nbcut = ncell
432 DO WHILE (icell<=ncell)
433 icell = icell +1
434 IF (icell>ncell .AND. ncell/=0)icell=9
435
436
437
438
439
440
441 do i=1,6
442 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(i)%NAdjCell
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460 enddo
461 enddo
462 enddo
463 print *, " ------------------------"
464 endif
465 endif
466 endif
467
468
469
470
471
473
474
475
476
477
478
479 DO ib=nbf,nbl
484 icell = 0
485 DO WHILE (icell<=ncell)
486 icell = icell +1
487 IF (icell>ncell .AND. ncell/=0)icell=9
488 DO j=1,6
489 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (1) = zero
490 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (2) = zero
491 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (3) = zero
492 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (4) = zero
493 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX (5) = zero
494 enddo
495
496 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 = zero
497
498
499
500 ie_m =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(3)
501 mat = ixs(1,ie_m)
502 upwl(1:6) = pm(16,mat)
503 reduc = pm(92,mat)
504 ddvol = zero
505 DO j=1,6
506 nadj =
brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
508 IF(iclose==1) nadj=0
509 DO iadj = 1,nadj
513 icellv =
brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
514 IF(idv==0)THEN
515 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*reduc
516 ELSEIF(idv>0)THEN
518 isilent = iparg(64,ng)
519 IF(isilent==1)THEN
520 upwl(j)=one
521 cellflux(j,icell,ib,iadj)=cellflux(j,icell,ib,iadj)*pm(92,ixs(1,idv))
522 ENDIF
523 ENDIF
524 brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_upwFLUX(iadj) =
525 . cellflux(j,icell,ib,iadj)-upwl(j)*abs(cellflux(j,icell,ib,iadj))
526 brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 =
527 .
brick_list(nin,ib)%POLY(icell)%Adjacent_FLU1 + cellflux(j,icell,ib,iadj)+upwl(j)*abs(cellflux(j,icell,ib,iadj))
528
529 ddvol = ddvol + cellflux(j,icell,ib,iadj)
530
531 enddo
532 enddo
533
535
536
537 enddo
538 enddo
539
540
542
543
544
545
546
547
548 IF(int22>0)THEN
549 nin = 1
550 DO ib=nbf,nbl
554 ddvol = zero
555 DO k=1,num
557 icellv =
brick_list(nin,ib)%SecndList%ICELLv(k)
558 ddvol = ddvol +
brick_list(nin,ibv)%POLY(icellv)%DDVOL
559 ENDDO
560 ddvol = ddvol +
brick_list(nin,ib)%POLY(mcell)%DDVOL
561
563 enddo
564 ENDIF
565
566
568
569
570 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
for(i8=*sizetab-1;i8 >=0;i8--)
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
int main(int argc, char *argv[])