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