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