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 allocate(ltg(nshg)) 214 do i =1,nshg 215 ltg(i)=i 216 enddo 217 nlwork=1 218 allocate( point2ilwork(1)) 219 nshg0 = nshg 220 endif 221 222 223 itwo=2 224 225 call phio_readheader(fhandle, 226 & c_char_'co-ordinates' // char(0), 227 & c_loc(intfromfile),itwo, dataDbl, iotype) 228 numnp=intfromfile(1) 229 allocate( point2x(numnp,nsd) ) 230 allocate( xread(numnp,nsd) ) 231 ixsiz=numnp*nsd 232 call phio_readdatablock(fhandle, 233 & c_char_'co-ordinates' // char(0), 234 & c_loc(xread),ixsiz, dataDbl, iotype) 235 point2x = xread 236 237c..............................for Duct 238 if(istretchOutlet.eq.1)then 239 240c...geometry6 241 if(iDuctgeometryType .eq. 6) then 242 xmaxn = 1.276 243 xmaxo = 0.848 244 xmin = 0.42 245c...geometry8 246 elseif(iDuctgeometryType .eq. 8)then 247 xmaxn=1.6*4.5*0.0254+0.85*1.5 248 xmaxo=1.6*4.5*0.0254+0.85*1.0 249 xmin =1.6*4.5*0.0254+0.85*0.5 250 endif 251c... 252 alpha=(xmaxn-xmaxo)/(xmaxo-xmin)**2 253 where (point2x(:,1) .ge. xmin) 254c..... N=# of current elements from .42 to exit(~40) 255c..... (x_mx-x_mn)/N=.025 256c..... alpha=3 3*.025=.075 257 point2x(:,1)=point2x(:,1)+ 258 & alpha*(point2x(:,1)-xmin)**2 259c..... ftn to stretch x at exit 260 endwhere 261 endif 262 263c 264c.... read in and block out the connectivity 265c 266 if( input_mode .eq. -1 ) then 267 call genblk (IBKSIZ) 268 else if( input_mode .eq. 0 ) then 269 call genblkPosix (IBKSIZ) 270 else if( input_mode .ge. 1 ) then 271 call genblkSyncIO (IBKSIZ) 272 end if 273c 274c.... read the boundary condition mapping array 275c 276 ione=1 277 call phio_readheader(fhandle, 278 & c_char_'bc mapping array' // char(0), 279 & c_loc(nshg),ione, dataInt, iotype) 280 281 allocate( nBC(nshg) ) 282 283 allocate( nBCread(nshg) ) 284 285 call phio_readdatablock(fhandle, 286 & c_char_'bc mapping array' // char(0), 287 & c_loc(nBCread), nshg, dataInt, iotype) 288 289 nBC=nBCread 290c 291c.... read the temporary iBC array 292c 293 ione=1 294 call phio_readheader(fhandle, 295 & c_char_'bc codes array' // char(0), 296 & c_loc(numpbc),ione, dataInt, iotype) 297 298 if ( numpbc > 0 ) then 299 allocate( iBCtmp(numpbc) ) 300 allocate( iBCtmpread(numpbc) ) 301 else 302 allocate( iBCtmp(1) ) 303 allocate( iBCtmpread(1) ) 304 endif 305 call phio_readdatablock(fhandle, 306 & c_char_'bc codes array' // char(0), 307 & c_loc(iBCtmpread), numpbc, dataInt, iotype) 308 309 if ( numpbc > 0 ) then 310 iBCtmp=iBCtmpread 311 else ! sometimes a partition has no BC's 312 deallocate( iBCtmpread) 313 iBCtmp=0 314 endif 315c 316c.... read boundary condition data 317c 318 ione=1 319 call phio_readheader(fhandle, 320 & c_char_'boundary condition array' // char(0), 321 & c_loc(intfromfile),ione, dataDbl, iotype) 322 323 if ( numpbc > 0 ) then 324 allocate( BCinp(numpbc,ndof+7) ) 325 nsecondrank=intfromfile(1)/numpbc 326 allocate( BCinpread(numpbc,nsecondrank) ) 327 iBCinpsiz=intfromfile(1) 328 else 329 allocate( BCinp(1,ndof+7) ) 330 allocate( BCinpread(0,0) ) !dummy 331 iBCinpsiz=intfromfile(1) 332 endif 333 334 call phio_readdatablock(fhandle, 335 & c_char_'boundary condition array' // char(0), 336 & c_loc(BCinpread), iBCinpsiz, dataDbl, iotype) 337 338 if ( numpbc > 0 ) then 339 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 340 else ! sometimes a partition has no BC's 341 deallocate(BCinpread) 342 BCinp=0 343 endif 344c 345c.... read periodic boundary conditions 346c 347 ione=1 348 call phio_readheader(fhandle, 349 & c_char_'periodic masters array' // char(0), 350 & c_loc(nshg), ione, dataInt, iotype) 351 allocate( point2iper(nshg) ) 352 allocate( iperread(nshg) ) 353 call phio_readdatablock(fhandle, 354 & c_char_'periodic masters array' // char(0), 355 & c_loc(iperread), nshg, dataInt, iotype) 356 point2iper=iperread 357c 358c.... generate the boundary element blocks 359c 360 if( input_mode .eq. -1 ) then 361 call genbkb (IBKSIZ) 362 else if( input_mode .eq. 0 ) then 363 call genbkbPosix (IBKSIZ) 364 else if( input_mode .ge. 1 ) then 365 call genbkbSyncIO (IBKSIZ) 366 end if 367c 368c Read in the nsons and ifath arrays if needed 369c 370c There is a fundamental shift in the meaning of ifath based on whether 371c there exist homogenous directions in the flow. 372c 373c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 374c points in the TOTAL mesh. That is to say that each partition keeps a 375c link to ALL inhomogenous points. This link is furthermore not to the 376c sms numbering but to the original structured grid numbering. These 377c inhomogenous points are thought of as fathers, with their sons being all 378c the points in the homogenous directions that have this father's 379c inhomogeneity. The array ifath takes as an arguement the sms numbering 380c and returns as a result the father. 381c 382c In this case nsons is the number of sons that each father has and ifath 383c is an array which tells the 384c 385c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 386c if we followed the above strategy since every partition would index its 387c points to the ENTIRE mesh. Furthermore, there would never be a need 388c to average to a node off processor since there is no spatial averaging. 389c Therefore, to properly account for this case we must recognize it and 390c inerrupt certain actions (i.e. assembly of the average across partitions). 391c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 392c father has any sons). Reiterating to be clear, in this case ifath does 393c not point to a global numbering but instead just points to itself. 394c 395 nfath=1 ! some architectures choke on a zero or undeclared 396 ! dimension variable. This sets it to a safe, small value. 397 if(((iLES .lt. 20) .and. (iLES.gt.0)) 398 & .or. (itwmod.gt.0) ) then ! don't forget same 399 ! conditional in proces.f 400 ! needed for alloc 401 ione=1 402 if(nohomog.gt.0) then 403 call phio_readheader(fhandle, 404 & c_char_'number of father-nodes' // char(0), 405 & c_loc(nfath), ione, dataInt, iotype) 406 407 call phio_readheader(fhandle, 408 & c_char_'number of son-nodes for each father' // char(0), 409 & c_loc(nfath), ione, dataInt, iotype) 410 411 allocate (point2nsons(nfath)) 412 413 call phio_readdatablock(fhandle, 414 & c_char_'number of son-nodes for each father' // char(0), 415 & c_loc(point2nsons),nfath, dataInt, iotype) 416 417 call phio_readheader(fhandle, 418 & c_char_'keyword ifath' // char(0), 419 & c_loc(nshg), ione, dataInt, iotype); 420 421 allocate (point2ifath(nshg)) 422 423 call phio_readdatablock(fhandle, 424 & c_char_'keyword ifath' // char(0), 425 & c_loc(point2ifath), nshg, dataInt, iotype) 426 427 nsonmax=maxval(point2nsons) 428 else ! this is the case where there is no homogeneity 429 ! therefore ever node is a father (too itself). sonfath 430 ! (a routine in NSpre) will set this up but this gives 431 ! you an option to avoid that. 432 nfath=nshg 433 allocate (point2nsons(nfath)) 434 point2nsons=1 435 allocate (point2ifath(nshg)) 436 do i=1,nshg 437 point2ifath(i)=i 438 enddo 439 nsonmax=1 440 endif 441 else 442 allocate (point2nsons(1)) 443 allocate (point2ifath(1)) 444 endif 445 446 call phio_closefile(fhandle); 447 iotime = TMRC() - iotime 448 if (myrank.eq.master) then 449 write(*,*) 'time to read geombc (seconds)', iotime 450 endif 451 452c.... Read restart files 453 iotime = TMRC() 454 if( input_mode .eq. -1 ) then 455 call streamio_setup_read(fhandle, geomRestartStream) 456 else if( input_mode .eq. 0 ) then 457 call posixio_setup(fhandle, c_char_'r') 458 else if( input_mode .ge. 1 ) then 459 call syncio_setup_read(nsynciofiles, fhandle) 460 end if 461 call phio_constructName(fhandle, 462 & c_char_'restart' // char(0), fnamer) 463 call phstr_appendInt(fnamer, irstart) 464 call phstr_appendStr(fnamer, c_char_'.'//c_null_char) 465 call phio_openfile(fnamer, fhandle); 466 467 ithree=3 468 469 itmp = int(log10(float(myrank+1)))+1 470 471 intfromfile=0 472 call phio_readheader(fhandle, 473 & c_char_'solution' // char(0), 474 & c_loc(intfromfile), ithree, dataInt, iotype) 475c 476c.... read the values of primitive variables into q 477c 478 allocate( qold(nshg,ndof) ) 479 if(intfromfile(1).ne.0) then 480 nshg2=intfromfile(1) 481 ndof2=intfromfile(2) 482 lstep=intfromfile(3) 483 if(ndof2.ne.ndof) then 484 485 endif 486 if (nshg2 .ne. nshg) 487 & call error ('restar ', 'nshg ', nshg) 488 allocate( qread(nshg,ndof2) ) 489 iqsiz=nshg*ndof2 490 call phio_readdatablock(fhandle, 491 & c_char_'solution' // char(0), 492 & c_loc(qread),iqsiz, dataDbl,iotype) 493 qold(:,1:ndof)=qread(:,1:ndof) 494 deallocate(qread) 495 else 496 if (myrank.eq.master) then 497 if (matflg(1,1).eq.0) then ! compressible 498 warning='Solution is set to zero (with p and T to one)' 499 else 500 warning='Solution is set to zero' 501 endif 502 write(*,*) warning 503 endif 504 qold=zero 505 if (matflg(1,1).eq.0) then ! compressible 506 qold(:,1)=one ! avoid zero pressure 507 qold(:,nflow)=one ! avoid zero temperature 508 endif 509 endif 510 511 intfromfile=0 512 call phio_readheader(fhandle, 513 & c_char_'time derivative of solution' // char(0), 514 & c_loc(intfromfile), ithree, dataInt, iotype) 515 allocate( acold(nshg,ndof) ) 516 if(intfromfile(1).ne.0) then 517 nshg2=intfromfile(1) 518 ndof2=intfromfile(2) 519 lstep=intfromfile(3) 520 521 if (nshg2 .ne. nshg) 522 & call error ('restar ', 'nshg ', nshg) 523 allocate( acread(nshg,ndof2) ) 524 acread=zero 525 iacsiz=nshg*ndof2 526 call phio_readdatablock(fhandle, 527 & c_char_'time derivative of solution' // char(0), 528 & c_loc(acread), iacsiz, dataDbl,iotype) 529 acold(:,1:ndof)=acread(:,1:ndof) 530 deallocate(acread) 531 else 532 if (myrank.eq.master) then 533 warning='Time derivative of solution is set to zero (SAFE)' 534 write(*,*) warning 535 endif 536 acold=zero 537 endif 538cc 539cc.... read the header and check it against the run data 540cc 541 if (ideformwall.eq.1) then 542 543 intfromfile=0 544 call phio_readheader(fhandle, 545 & c_char_'displacement' // char(0), 546 & c_loc(intfromfile), ithree, dataInt, iotype) 547 548 nshg2=intfromfile(1) 549 ndisp=intfromfile(2) 550 lstep=intfromfile(3) 551 if(ndisp.ne.nsd) then 552 warning='WARNING ndisp not equal nsd' 553 write(*,*) warning , ndisp 554 endif 555 if (nshg2 .ne. nshg) 556 & call error ('restar ', 'nshg ', nshg) 557c 558c.... read the values of primitive variables into uold 559c 560 allocate( uold(nshg,nsd) ) 561 allocate( uread(nshg,nsd) ) 562 563 iusiz=nshg*nsd 564 565 call phio_readdatablock(fhandle, 566 & c_char_'displacement' // char(0), 567 & c_loc(uread), iusiz, dataDbl, iotype) 568 569 uold(:,1:nsd)=uread(:,1:nsd) 570 else 571 allocate( uold(nshg,nsd) ) 572 uold(:,1:nsd) = zero 573 endif 574c 575c.... close c-binary files 576c 577 call phio_closefile(fhandle) 578 iotime = TMRC() - iotime 579 if (myrank.eq.master) then 580 write(*,*) 'time to read restart (seconds)', iotime 581 endif 582 583 deallocate(xread) 584 if ( numpbc > 0 ) then 585 deallocate(bcinpread) 586 deallocate(ibctmpread) 587 endif 588 deallocate(iperread) 589 if(numpe.gt.1) 590 & deallocate(ilworkread) 591 deallocate(nbcread) 592 593 return 594 994 call error ('input ','opening ', igeomBAK) 595 995 call error ('input ','opening ', igeomBAK) 596 997 call error ('input ','end file', igeomBAK) 597 998 call error ('input ','end file', igeomBAK) 598 end 599c 600c No longer called but kept around in case.... 601c 602 subroutine genpzero(iBC) 603 604 use pointer_data 605 include "common.h" 606 integer iBC(nshg) 607c 608c.... check to see if any of the nodes have a dirichlet pressure 609c 610 pzero=1 611 if (any(btest(iBC,2))) pzero=0 612 do iblk = 1, nelblb 613 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 614 do i=1, npro 615 iBCB1=miBCB(iblk)%p(i,1) 616c 617c.... check to see if any of the nodes have a Neumann pressure 618c but not periodic (note that 619c 620 if(btest(iBCB1,1)) pzero=0 621 enddo 622c 623c.... share results with other processors 624c 625 pzl=pzero 626 if (numpe .gt. 1) 627 & call MPI_ALLREDUCE (pzl, pzero, 1, 628 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 629 enddo 630 return 631 end 632