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