1c readnblk.f (pronounce "Reed and Block Dot Eff") contains: 2c 3c module readarrays ("Red Arrays") -- contains the arrays that 4c are read in from binary files but not immediately blocked 5c through pointers. 6c 7c subroutine readnblk ("Reed and Block") -- allocates space for 8c and reads data to be contained in module readarrays. Reads 9c all remaining data and blocks them with pointers. 10c 11 module readarrays 12 13 real*8, allocatable :: point2x(:,:) 14 real*8, allocatable :: qold(:,:) 15 real*8, allocatable :: uold(:,:) 16 real*8, allocatable :: acold(:,:) 17 integer, allocatable :: iBCtmp(:) 18 real*8, allocatable :: BCinp(:,:) 19 20 integer, allocatable :: point2ilwork(:) 21! integer, allocatable :: fncorp(:) 22 integer, allocatable :: twodncorp(:,:) 23 integer, allocatable :: nBC(:) 24 integer, allocatable :: point2iper(:) 25 integer, target, allocatable :: point2ifath(:) 26 integer, target, allocatable :: point2nsons(:) 27 28 end module 29 30 subroutine readnblk 31 use iso_c_binding 32 use readarrays 33 use fncorpmod 34 use phio 35 use phstr 36 use syncio 37 use posixio 38 use streamio 39 include "common.h" 40 41 real*8, target, allocatable :: xread(:,:), qread(:,:), acread(:,:) 42 real*8, target, allocatable :: uread(:,:) 43 real*8, target, allocatable :: BCinpread(:,:) 44 real*8 :: iotime 45 integer, target, allocatable :: iperread(:), iBCtmpread(:) 46 integer, target, allocatable :: ilworkread(:), nBCread(:) 47 integer, target, allocatable :: fncorpread(:) 48 integer fncorpsize 49 character*10 cname2, cname2nd 50 character*8 mach2 51 character*30 fmt1 52 character*255 fname1,fnamer,fnamelr 53 character*255 warning 54 integer igeomBAK, ibndc, irstin, ierr 55 integer, target :: intfromfile(50) ! integers read from headers 56 integer :: descriptor, descriptorG, GPID, color, nfields 57 integer :: numparts, nppf 58 integer :: ierr_io, numprocs, itmp, itmp2 59 integer :: ignored 60 integer :: fileFmt 61 character*255 fname2, temp2 62 character*64 temp1 63 type(c_ptr) :: handle 64 character(len=1024) :: dataInt, dataDbl 65 dataInt = c_char_'integer'//c_null_char 66 dataDbl = c_char_'double'//c_null_char 67c 68c.... determine the step number to start with 69c 70 open(unit=72,file='numstart.dat',status='old') 71 read(72,*) irstart 72 close(72) 73 lstep=irstart ! in case restart files have no fields 74 75 fname1='geombc.dat' 76 fname1= trim(fname1) // cname2(myrank+1) 77 78 itmp=1 79 if (irstart .gt. 0) itmp = int(log10(float(irstart)))+1 80 write (fmt1,"('(''restart.'',i',i1,',1x)')") itmp 81 write (fnamer,fmt1) irstart 82 fnamer = trim(fnamer) // cname2(myrank+1) 83c 84c.... open input files 85c.... input the geometry parameters 86c 87 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 88 89 itwo=2 90 ione=1 91 ieleven=11 92 itmp = int(log10(float(myrank+1)))+1 93 94 iotime = TMRC() 95 if( input_mode .eq. -1 ) then 96 call streamio_setup_read(fhandle, geomRestartStream) 97 else if( input_mode .eq. 0 ) then 98 call posixio_setup(fhandle, c_char_'r') 99 else if( input_mode .ge. 1 ) then 100 call syncio_setup_read(nsynciofiles, fhandle) 101 end if 102 call phio_constructName(fhandle, 103 & c_char_'geombc' // char(0), fname1) 104 call phio_openfile(fname1, fhandle); 105 106 call phio_readheader(fhandle,c_char_'number of nodes' // char(0), 107 & c_loc(numnp),ione, dataInt, iotype) 108 109 call phio_readheader(fhandle,c_char_'number of modes' // char(0), 110 & c_loc(nshg),ione, dataInt, iotype) 111 112 call phio_readheader(fhandle, 113 & c_char_'number of interior elements' // char(0), 114 & c_loc(numel),ione, dataInt, iotype) 115 116 call phio_readheader(fhandle, 117 & c_char_'number of boundary elements' // char(0), 118 & c_loc(numelb),ione, dataInt, iotype) 119 120 call phio_readheader(fhandle, 121 & c_char_'maximum number of element nodes' // char(0), 122 & c_loc(nen),ione, dataInt, iotype) 123 124 call phio_readheader(fhandle, 125 & c_char_'number of interior tpblocks' // char(0), 126 & c_loc(nelblk),ione, dataInt, iotype) 127 128 call phio_readheader(fhandle, 129 & c_char_'number of boundary tpblocks' // char(0), 130 & c_loc(nelblb),ione, dataInt, iotype) 131 132 call phio_readheader(fhandle, 133 & c_char_'number of nodes with Dirichlet BCs' // char(0), 134 & c_loc(numpbc),ione, dataInt, iotype) 135 136 call phio_readheader(fhandle, 137 & c_char_'number of shape functions' // char(0), 138 & c_loc(ntopsh),ione, dataInt, iotype) 139c 140c.... calculate the maximum number of boundary element nodes 141c 142 nenb = 0 143 do i = 1, melCat 144 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 145 enddo 146 if (myrank == master) then 147 if (nenb .eq. 0) call error ('input ','nen ',nen) 148 endif 149c 150c.... setup some useful constants 151c 152 I3nsd = nsd / 3 ! nsd=3 integer flag 153 E3nsd = float(I3nsd) ! nsd=3 real flag 154 if(matflg(1,1).lt.0) then 155 nflow = nsd + 1 156 else 157 nflow = nsd + 2 158 endif 159 ndof = nsd + 2 160 nsclr=impl(1)/100 161 ndof=ndof+nsclr ! number of sclr transport equations to solve 162 163 ndofBC = ndof + I3nsd ! dimension of BC array 164 ndiBCB = 2 ! dimension of iBCB array 165 ndBCB = ndof + 1 ! dimension of BCB array 166 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 167c 168c.... ----------------------> Communication tasks <-------------------- 169c 170 if(numpe > 1) then 171 call phio_readheader(fhandle, 172 & c_char_'size of ilwork array' // char(0), 173 & c_loc(nlwork),ione, dataInt, iotype) 174 175 call phio_readheader(fhandle, 176 & c_char_'ilwork' //char(0), 177 & c_loc(nlwork),ione, dataInt, iotype) 178 179 allocate( point2ilwork(nlwork) ) 180 allocate( ilworkread(nlwork) ) 181 call phio_readdatablock(fhandle, c_char_'ilwork' // char(0), 182 & c_loc(ilworkread), nlwork, dataInt , iotype) 183 184 point2ilwork = ilworkread 185 call ctypes (point2ilwork) 186 187 if((usingPETSc.eq.1).or.(svLSFlag.eq.1)) then 188 fncorpsize = nshg 189 allocate(fncorp(fncorpsize)) 190 call gen_ncorp(fncorp, ilworkread, nlwork, fncorpsize) 191! 192! the following code finds the global range of the owned nodes 193! 194 maxowned=0 195 minowned=maxval(fncorp) 196 do i = 1,nshg 197 if(fncorp(i).gt.0) then ! don't consider remote copies 198 maxowned=max(maxowned,fncorp(i)) 199 minowned=min(minowned,fncorp(i)) 200 endif 201 enddo 202! 203! end of global range code 204! 205 call commuInt(fncorp, point2ilwork, 1, 'out') 206 ncorpsize = fncorpsize 207 endif 208 if(svLSFlag.eq.1) then 209 allocate(ltg(ncorpsize)) 210 ltg(:)=fncorp(:) 211 endif 212 else 213 if((usingPETSc.eq.1)) then 214 fncorpsize = nshg 215 allocate(fncorp(fncorpsize)) 216 maxowned=nshg 217 minowned=1 218 iownnodes=nshg 219 do i =1,nshg 220 fncorp(i)=i 221 enddo 222 endif 223 if((svLSFlag.eq.1)) then 224 allocate(ltg(nshg)) 225 do i =1,nshg 226 ltg(i)=i 227 enddo 228 endif 229 nlwork=1 230 allocate( point2ilwork(1)) 231 nshg0 = nshg 232 endif 233 234 235 itwo=2 236 237 call phio_readheader(fhandle, 238 & c_char_'co-ordinates' // char(0), 239 & c_loc(intfromfile),itwo, dataDbl, iotype) 240 numnp=intfromfile(1) 241 allocate( point2x(numnp,nsd) ) 242 allocate( xread(numnp,nsd) ) 243 ixsiz=numnp*nsd 244 call phio_readdatablock(fhandle, 245 & c_char_'co-ordinates' // char(0), 246 & c_loc(xread),ixsiz, dataDbl, iotype) 247 point2x = xread 248 249c..............................for Duct 250 if(istretchOutlet.eq.1)then 251 252c...geometry6 253 if(iDuctgeometryType .eq. 6) then 254 xmaxn = 1.276 255 xmaxo = 0.848 256 xmin = 0.42 257c...geometry8 258 elseif(iDuctgeometryType .eq. 8)then 259 xmaxn=1.6*4.5*0.0254+0.85*1.5 260 xmaxo=1.6*4.5*0.0254+0.85*1.0 261 xmin =1.6*4.5*0.0254+0.85*0.5 262 endif 263c... 264 alpha=(xmaxn-xmaxo)/(xmaxo-xmin)**2 265 where (point2x(:,1) .ge. xmin) 266c..... N=# of current elements from .42 to exit(~40) 267c..... (x_mx-x_mn)/N=.025 268c..... alpha=3 3*.025=.075 269 point2x(:,1)=point2x(:,1)+ 270 & alpha*(point2x(:,1)-xmin)**2 271c..... ftn to stretch x at exit 272 endwhere 273 endif 274 275c 276c.... read in and block out the connectivity 277c 278 if( input_mode .eq. -1 ) then 279 call genblk (IBKSIZ) 280 else if( input_mode .eq. 0 ) then 281 call genblkPosix (IBKSIZ) 282 else if( input_mode .ge. 1 ) then 283 call genblkSyncIO (IBKSIZ) 284 end if 285c 286c.... read the boundary condition mapping array 287c 288 ione=1 289 call phio_readheader(fhandle, 290 & c_char_'bc mapping array' // char(0), 291 & c_loc(nshg),ione, dataInt, iotype) 292 293 allocate( nBC(nshg) ) 294 295 allocate( nBCread(nshg) ) 296 297 call phio_readdatablock(fhandle, 298 & c_char_'bc mapping array' // char(0), 299 & c_loc(nBCread), nshg, dataInt, iotype) 300 301 nBC=nBCread 302c 303c.... read the temporary iBC array 304c 305 ione=1 306 call phio_readheader(fhandle, 307 & c_char_'bc codes array' // char(0), 308 & c_loc(numpbc),ione, dataInt, iotype) 309 310 if ( numpbc > 0 ) then 311 allocate( iBCtmp(numpbc) ) 312 allocate( iBCtmpread(numpbc) ) 313 else 314 allocate( iBCtmp(1) ) 315 allocate( iBCtmpread(1) ) 316 endif 317 call phio_readdatablock(fhandle, 318 & c_char_'bc codes array' // char(0), 319 & c_loc(iBCtmpread), numpbc, dataInt, iotype) 320 321 if ( numpbc > 0 ) then 322 iBCtmp=iBCtmpread 323 else ! sometimes a partition has no BC's 324 deallocate( iBCtmpread) 325 iBCtmp=0 326 endif 327c 328c.... read boundary condition data 329c 330 ione=1 331 call phio_readheader(fhandle, 332 & c_char_'boundary condition array' // char(0), 333 & c_loc(intfromfile),ione, dataDbl, iotype) 334 335 if ( numpbc > 0 ) then 336 allocate( BCinp(numpbc,ndof+7) ) 337 nsecondrank=intfromfile(1)/numpbc 338 allocate( BCinpread(numpbc,nsecondrank) ) 339 iBCinpsiz=intfromfile(1) 340 else 341 allocate( BCinp(1,ndof+7) ) 342 allocate( BCinpread(0,0) ) !dummy 343 iBCinpsiz=intfromfile(1) 344 endif 345 346 call phio_readdatablock(fhandle, 347 & c_char_'boundary condition array' // char(0), 348 & c_loc(BCinpread), iBCinpsiz, dataDbl, iotype) 349 350 if ( numpbc > 0 ) then 351 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 352 else ! sometimes a partition has no BC's 353 deallocate(BCinpread) 354 BCinp=0 355 endif 356c 357c.... read periodic boundary conditions 358c 359 ione=1 360 call phio_readheader(fhandle, 361 & c_char_'periodic masters array' // char(0), 362 & c_loc(nshg), ione, dataInt, iotype) 363 allocate( point2iper(nshg) ) 364 allocate( iperread(nshg) ) 365 call phio_readdatablock(fhandle, 366 & c_char_'periodic masters array' // char(0), 367 & c_loc(iperread), nshg, dataInt, iotype) 368 point2iper=iperread 369c 370c.... generate the boundary element blocks 371c 372 if( input_mode .eq. -1 ) then 373 call genbkb (IBKSIZ) 374 else if( input_mode .eq. 0 ) then 375 call genbkbPosix (IBKSIZ) 376 else if( input_mode .ge. 1 ) then 377 call genbkbSyncIO (IBKSIZ) 378 end if 379c 380c Read in the nsons and ifath arrays if needed 381c 382c There is a fundamental shift in the meaning of ifath based on whether 383c there exist homogenous directions in the flow. 384c 385c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 386c points in the TOTAL mesh. That is to say that each partition keeps a 387c link to ALL inhomogenous points. This link is furthermore not to the 388c sms numbering but to the original structured grid numbering. These 389c inhomogenous points are thought of as fathers, with their sons being all 390c the points in the homogenous directions that have this father's 391c inhomogeneity. The array ifath takes as an arguement the sms numbering 392c and returns as a result the father. 393c 394c In this case nsons is the number of sons that each father has and ifath 395c is an array which tells the 396c 397c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 398c if we followed the above strategy since every partition would index its 399c points to the ENTIRE mesh. Furthermore, there would never be a need 400c to average to a node off processor since there is no spatial averaging. 401c Therefore, to properly account for this case we must recognize it and 402c inerrupt certain actions (i.e. assembly of the average across partitions). 403c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 404c father has any sons). Reiterating to be clear, in this case ifath does 405c not point to a global numbering but instead just points to itself. 406c 407 nfath=1 ! some architectures choke on a zero or undeclared 408 ! dimension variable. This sets it to a safe, small value. 409 if(((iLES .lt. 20) .and. (iLES.gt.0)) 410 & .or. (itwmod.gt.0) ) then ! don't forget same 411 ! conditional in proces.f 412 ! needed for alloc 413 ione=1 414 if(nohomog.gt.0) then 415 call phio_readheader(fhandle, 416 & c_char_'number of father-nodes' // char(0), 417 & c_loc(nfath), ione, dataInt, iotype) 418 419 call phio_readheader(fhandle, 420 & c_char_'number of son-nodes for each father' // char(0), 421 & c_loc(nfath), ione, dataInt, iotype) 422 423 allocate (point2nsons(nfath)) 424 425 call phio_readdatablock(fhandle, 426 & c_char_'number of son-nodes for each father' // char(0), 427 & c_loc(point2nsons),nfath, dataInt, iotype) 428 429 call phio_readheader(fhandle, 430 & c_char_'keyword ifath' // char(0), 431 & c_loc(nshg), ione, dataInt, iotype); 432 433 allocate (point2ifath(nshg)) 434 435 call phio_readdatablock(fhandle, 436 & c_char_'keyword ifath' // char(0), 437 & c_loc(point2ifath), nshg, dataInt, iotype) 438 439 nsonmax=maxval(point2nsons) 440 else ! this is the case where there is no homogeneity 441 ! therefore ever node is a father (too itself). sonfath 442 ! (a routine in NSpre) will set this up but this gives 443 ! you an option to avoid that. 444 nfath=nshg 445 allocate (point2nsons(nfath)) 446 point2nsons=1 447 allocate (point2ifath(nshg)) 448 do i=1,nshg 449 point2ifath(i)=i 450 enddo 451 nsonmax=1 452 endif 453 else 454 allocate (point2nsons(1)) 455 allocate (point2ifath(1)) 456 endif 457 458 call phio_closefile(fhandle); 459 iotime = TMRC() - iotime 460 if (myrank.eq.master) then 461 write(*,*) 'time to read geombc (seconds)', iotime 462 endif 463 464c.... Read restart files 465 iotime = TMRC() 466 if( input_mode .eq. -1 ) then 467 call streamio_setup_read(fhandle, geomRestartStream) 468 else if( input_mode .eq. 0 ) then 469 call posixio_setup(fhandle, c_char_'r') 470 else if( input_mode .ge. 1 ) then 471 call syncio_setup_read(nsynciofiles, fhandle) 472 end if 473 call phio_constructName(fhandle, 474 & c_char_'restart' // char(0), fnamer) 475 call phstr_appendInt(fnamer, irstart) 476 call phstr_appendStr(fnamer, c_char_'.'//c_null_char) 477 call phio_openfile(fnamer, fhandle); 478 479 ithree=3 480 481 itmp = int(log10(float(myrank+1)))+1 482 483 intfromfile=0 484 call phio_readheader(fhandle, 485 & c_char_'solution' // char(0), 486 & c_loc(intfromfile), ithree, dataInt, iotype) 487c 488c.... read the values of primitive variables into q 489c 490 allocate( qold(nshg,ndof) ) 491 if(intfromfile(1).ne.0) then 492 nshg2=intfromfile(1) 493 ndof2=intfromfile(2) 494 lstep=intfromfile(3) 495 if(ndof2.ne.ndof) then 496 497 endif 498 if (nshg2 .ne. nshg) 499 & call error ('restar ', 'nshg ', nshg) 500 allocate( qread(nshg,ndof2) ) 501 iqsiz=nshg*ndof2 502 call phio_readdatablock(fhandle, 503 & c_char_'solution' // char(0), 504 & c_loc(qread),iqsiz, dataDbl,iotype) 505 qold(:,1:ndof)=qread(:,1:ndof) 506 deallocate(qread) 507 else 508 if (myrank.eq.master) then 509 if (matflg(1,1).eq.0) then ! compressible 510 warning='Solution is set to zero (with p and T to one)' 511 else 512 warning='Solution is set to zero' 513 endif 514 write(*,*) warning 515 endif 516 qold=zero 517 if (matflg(1,1).eq.0) then ! compressible 518 qold(:,1)=one ! avoid zero pressure 519 qold(:,nflow)=one ! avoid zero temperature 520 endif 521 endif 522 523 intfromfile=0 524 call phio_readheader(fhandle, 525 & c_char_'time derivative of solution' // char(0), 526 & c_loc(intfromfile), ithree, dataInt, iotype) 527 allocate( acold(nshg,ndof) ) 528 if(intfromfile(1).ne.0) then 529 nshg2=intfromfile(1) 530 ndof2=intfromfile(2) 531 lstep=intfromfile(3) 532 533 if (nshg2 .ne. nshg) 534 & call error ('restar ', 'nshg ', nshg) 535 allocate( acread(nshg,ndof2) ) 536 acread=zero 537 iacsiz=nshg*ndof2 538 call phio_readdatablock(fhandle, 539 & c_char_'time derivative of solution' // char(0), 540 & c_loc(acread), iacsiz, dataDbl,iotype) 541 acold(:,1:ndof)=acread(:,1:ndof) 542 deallocate(acread) 543 else 544 if (myrank.eq.master) then 545 warning='Time derivative of solution is set to zero (SAFE)' 546 write(*,*) warning 547 endif 548 acold=zero 549 endif 550cc 551cc.... read the header and check it against the run data 552cc 553 if (ideformwall.eq.1) then 554 555 intfromfile=0 556 call phio_readheader(fhandle, 557 & c_char_'displacement' // char(0), 558 & c_loc(intfromfile), ithree, dataInt, iotype) 559 560 nshg2=intfromfile(1) 561 ndisp=intfromfile(2) 562 lstep=intfromfile(3) 563 if(ndisp.ne.nsd) then 564 warning='WARNING ndisp not equal nsd' 565 write(*,*) warning , ndisp 566 endif 567 if (nshg2 .ne. nshg) 568 & call error ('restar ', 'nshg ', nshg) 569c 570c.... read the values of primitive variables into uold 571c 572 allocate( uold(nshg,nsd) ) 573 allocate( uread(nshg,nsd) ) 574 575 iusiz=nshg*nsd 576 577 call phio_readdatablock(fhandle, 578 & c_char_'displacement' // char(0), 579 & c_loc(uread), iusiz, dataDbl, iotype) 580 581 uold(:,1:nsd)=uread(:,1:nsd) 582 else 583 allocate( uold(nshg,nsd) ) 584 uold(:,1:nsd) = zero 585 endif 586c 587c.... close c-binary files 588c 589 call phio_closefile(fhandle) 590 iotime = TMRC() - iotime 591 if (myrank.eq.master) then 592 write(*,*) 'time to read restart (seconds)', iotime 593 endif 594 595 deallocate(xread) 596 if ( numpbc > 0 ) then 597 deallocate(bcinpread) 598 deallocate(ibctmpread) 599 endif 600 deallocate(iperread) 601 if(numpe.gt.1) 602 & deallocate(ilworkread) 603 deallocate(nbcread) 604 605 return 606 994 call error ('input ','opening ', igeomBAK) 607 995 call error ('input ','opening ', igeomBAK) 608 997 call error ('input ','end file', igeomBAK) 609 998 call error ('input ','end file', igeomBAK) 610 end 611c 612c No longer called but kept around in case.... 613c 614 subroutine genpzero(iBC) 615 616 use pointer_data 617 include "common.h" 618 integer iBC(nshg) 619c 620c.... check to see if any of the nodes have a dirichlet pressure 621c 622 pzero=1 623 if (any(btest(iBC,2))) pzero=0 624 do iblk = 1, nelblb 625 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 626 do i=1, npro 627 iBCB1=miBCB(iblk)%p(i,1) 628c 629c.... check to see if any of the nodes have a Neumann pressure 630c but not periodic (note that 631c 632 if(btest(iBCB1,1)) pzero=0 633 enddo 634c 635c.... share results with other processors 636c 637 pzl=pzero 638 if (numpe .gt. 1) 639 & call MPI_ALLREDUCE (pzl, pzero, 1, 640 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 641 enddo 642 return 643 end 644