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