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 12 13 module readarrays 14 15 real*8, allocatable :: point2x(:,:) 16 real*8, allocatable :: qold(:,:) 17 real*8, allocatable :: uold(:,:) 18 real*8, allocatable :: acold(:,:) 19 integer, allocatable :: iBCtmp(:) 20 real*8, allocatable :: BCinp(:,:) 21 22 integer, allocatable :: point2ilwork(:) 23 integer, allocatable :: nBC(:) 24 integer, allocatable :: point2iper(:) 25 integer, allocatable :: point2ifath(:) 26 integer, allocatable :: point2nsons(:) 27 28 end module 29 30 31 32 subroutine readnblk 33c 34 use readarrays 35 include "common.h" 36c 37 real*8, allocatable :: xread(:,:), qread(:,:), acread(:,:) 38 real*8, allocatable :: uread(:,:) 39 real*8, allocatable :: BCinpread(:,:) 40 integer, allocatable :: iperread(:), iBCtmpread(:) 41 integer, allocatable :: ilworkread(:), nBCread(:) 42 character*10 cname2 43 character*8 mach2 44!MR CHANGE 45! character*20 fmt1 46 character*30 fmt1 47!MR CHANGE END 48 character*255 fname1,fnamer,fnamelr 49 character*255 warning 50 integer igeomBAK, ibndc, irstin, ierr 51 integer intfromfile(50) ! integers read from headers 52 53cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 54ccccccccccccccccccccccccccccccccccccccccc New PhastaIO Definition Part ccccccccccccccccccccccccccccccccccccccccc 55 56 integer :: descriptor, descriptorG, GPID, color, nfiles, nfields 57 integer :: numparts, nppf 58 integer :: ierr_io, numprocs, itmp, itmp2 59 character*255 fname2, temp2 60 character*64 temp1 61 62cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 63cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 64 65 66c 67c.... determine the step number to start with 68c 69 open(unit=72,file='numstart.dat',status='old') 70 read(72,*) irstart 71 close(72) 72 lstep=irstart ! in case restart files have no fields 73 74c 75 fname1='geombc.dat' 76 fname1= trim(fname1) // cname2(myrank+1) 77c fnamelr='restart.latest' 78 79 itmp=1 80 if (irstart .gt. 0) itmp = int(log10(float(irstart)))+1 81 write (fmt1,"('(''restart.'',i',i1,',1x)')") itmp 82 write (fnamer,fmt1) irstart 83 fnamer = trim(fnamer) // cname2(myrank+1) 84c fnamelr = trim(fnamelr) // cname2(myrank+1) 85 86c 87c.... open input files 88c 89 90 91c call openfile( fname1, 'read?', igeomBAK ); 92 93 94c 95c.... try opening restart.latest.proc before trying restart.stepno.proc 96c 97c call openfile( fnamelr, 'read?', irstin ); 98c if ( irstin .eq. 0 ) 99 100!MR CHANGE 101! call openfile( fnamer, 'read?', irstin ); 102!MR CHANGE END 103 104! either one will work 105c 106c.... input the geometry parameters 107c 108 109cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 110!MR CHANGE 111 112! nfiles = 2 113! nfields = 31 114! numparts = 8 115! nppp = numparts/numpe 116! nppf = numparts/nfiles 117 118 nfiles = nsynciofiles 119! nfields = nsynciofieldsreadgeombc 120 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 121 122! nppp = numparts/numpe 123! nppf = numparts/nfiles 124!MR CHANGE END 125 126c fnamer="/home/nliu/develop/PIG4/512-procs_case/geombc-dat" 127 128 color = int(myrank/(numparts/nfiles)) !Should call the color routine in SyncIO here 129 itmp2 = int(log10(float(color+1)))+1 130 write (temp2,"('(''geombc-dat.'',i',i1,')')") itmp2 131 write (fnamer,temp2) (color+1) 132 fnamer = trim(fnamer) // char(0) 133 134c fnamer="/home/nliu/develop/test-case/512-procs_case/geombc-dat" 135 136 itwo=2 137 ione=1 138 ieleven=11 139 itmp = int(log10(float(myrank+1)))+1 140 141 call queryphmpiio(fnamer, nfields, nppf); 142 if (myrank == 0) then 143 write(*,*) 'Number of fields in geombc-dat: ',nfields 144 write(*,*) 'Number of parts per file geombc-dat: ',nppf 145 endif 146 call initphmpiio( nfields, nppf, nfiles, igeom, 147 & 'read' // char(0)) 148 call openfile( fnamer, 'read' // char(0), igeom ) 149 150 write (temp1,"('(''number of nodes@'',i',i1,',A1)')") itmp 151 write (fname2,temp1) (myrank+1),'?' 152 call readheader(igeom,fname2 // char(0),numnp,ione, 153 & 'integer' // char(0), iotype) 154 155 write (temp1,"('(''number of modes@'',i',i1,',A1)')") itmp 156 write (fname2,temp1) (myrank+1),'?' 157 call readheader(igeom,fname2 // char(0),nshg,ione, 158 & 'integer' // char(0), iotype) 159 160 write (temp1,"('(''number of interior elements@'',i',i1,',A1)')") 161 & itmp 162 write (fname2,temp1) (myrank+1),'?' 163 call readheader(igeom,fname2 // char(0),numel,ione, 164 & 'integer' // char(0), iotype) 165 166 write (temp1,"('(''number of boundary elements@'',i',i1,',A1)')") 167 & itmp 168 write (fname2,temp1) (myrank+1),'?' 169 call readheader(igeom,fname2 // char(0),numelb,ione, 170 & 'integer' // char(0),iotype) 171 172 write (temp1, 173 & "('(''maximum number of element nodes@'',i',i1,',A1)')") itmp 174 write (fname2,temp1) (myrank+1),'?' 175 call readheader(igeom,fname2 // char(0),nen,ione, 176 &'integer' // char(0),iotype) 177 178 write (temp1,"('(''number of interior tpblocks@'',i',i1,',A1)')") 179 & itmp 180 write (fname2,temp1) (myrank+1),'?' 181 call readheader(igeom,fname2 // char(0),nelblk,ione, 182 & 'integer' // char(0) ,iotype) 183 184 write (temp1,"('(''number of boundary tpblocks@'',i',i1,',A1)')") 185 & itmp 186 write (fname2,temp1) (myrank+1),'?' 187 call readheader(igeom,fname2 // char(0),nelblb,ione, 188 & 'integer' // char(0), iotype) 189 190 write (temp1, 191 & "('(''number of nodes with Dirichlet BCs@'',i',i1,',A1)')") itmp 192 write (fname2,temp1) (myrank+1),'?' 193 call readheader(igeom,fname2 // char(0),numpbc,ione, 194 & 'integer' // char(0),iotype) 195 196 write (temp1,"('(''number of shape functions@'',i',i1,',A1)')") 197 & itmp 198 write (fname2,temp1) (myrank+1),'?' 199 call readheader(igeom,fname2 // char(0),ntopsh,ione, 200 & 'integer' // char(0),iotype) 201 202c call closefile( igeom, "read" ) 203c call finalizephmpiio( igeom ) 204 205! if(myrank==0) then 206! print *, "Reading Finished, ********* : ", trim(fname2) 207! endif 208 209 210ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 211 212c ieleven=11 213c ione=1 214c fname1='number of nodes?' 215c call readheader(igeomBAK,fname1,numnp,ione,'integer', iotype) 216c fname1='number of modes?' 217c call readheader(igeomBAK,fname1,nshg,ione,'integer', iotype) 218cc fname1='number of global modes?' 219cc call readheader(igeomBAK,fname1,nshgt,ione,'integer', iotype) 220c fname1='number of interior elements?' 221c call readheader(igeomBAK,fname1,numel,ione,'integer', iotype) 222c fname1='number of boundary elements?' 223c call readheader(igeomBAK,fname1,numelb,ione,'integer', iotype) 224c fname1='maximum number of element nodes?' 225c call readheader(igeomBAK,fname1,nen,ione,'integer', iotype) 226c fname1='number of interior tpblocks?' 227c call readheader(igeomBAK,fname1,nelblk,ione,'integer', iotype) 228c fname1='number of boundary tpblocks?' 229c call readheader(igeomBAK,fname1,nelblb,ione,'integer', iotype) 230c fname1='number of nodes with Dirichlet BCs?' 231c call readheader(igeomBAK,fname1,numpbc,ione,'integer', iotype) 232c fname1='number of shape functions?' 233c call readheader(igeomBAK,fname1,ntopsh,ione,'integer', iotype) 234 235c 236c.... calculate the maximum number of boundary element nodes 237c 238 nenb = 0 239 do i = 1, melCat 240 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 241 enddo 242c 243 if (myrank == master) then 244 if (nenb .eq. 0) call error ('input ','nen ',nen) 245 endif 246c 247c.... setup some useful constants 248c 249 I3nsd = nsd / 3 ! nsd=3 integer flag 250 E3nsd = float(I3nsd) ! nsd=3 real flag 251c 252 if(matflg(1,1).lt.0) then 253 nflow = nsd + 1 254 else 255 nflow = nsd + 2 256 endif 257 ndof = nsd + 2 258 nsclr=impl(1)/100 259 ndof=ndof+nsclr ! number of sclr transport equations to solve 260 261 ndofBC = ndof + I3nsd ! dimension of BC array 262 ndiBCB = 2 ! dimension of iBCB array 263 ndBCB = ndof + 1 ! dimension of BCB array 264c 265 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 266c 267c.... ----------------------> Communication tasks <-------------------- 268c 269 270cc still read in new 271 272 if(numpe > 1) then 273 274cc MR CHANGE 275 write (temp1,"('(''size of ilwork array@'',i',i1,',A1)')") itmp 276 write (fname2,temp1) (myrank+1),'?' 277 call readheader(igeom,fname2 // char(0),nlwork,ione, 278 & 'integer' // char(0) ,iotype) 279 280 write (temp1,"('(''ilwork@'',i',i1,',A1)')") itmp 281 write (fname2,temp1) (myrank+1),'?' 282 call readheader(igeom,fname2 //char(0) ,nlwork,ione, 283 & 'integer' // char(0) , iotype) 284 285 allocate( point2ilwork(nlwork) ) 286 allocate( ilworkread(nlwork) ) 287 call readdatablock(igeom,fname2 // char(0),ilworkread, 288 & nlwork,'integer' // char(0) , iotype) 289 290c call closefile( igeom, "read" ) 291c call finalizephmpiio( igeom ) 292 293c fname1='size of ilwork array?' 294c call readheader(igeomBAK,fname1,nlwork,ione,'integer', iotype) 295 296c ione=1 297c fname1='ilwork?' 298c call readheader(igeomBAK,fname1,nlwork,ione,'integer', iotype) 299 300c allocate( point2ilwork(nlwork) ) 301c allocate( ilworkread(nlwork) ) 302c call readdatablock(igeomBAK,fname1,ilworkread, 303c & nlwork,'integer', iotype) 304cc MR CHANGE 305 306 307 point2ilwork = ilworkread 308 call ctypes (point2ilwork) 309 else 310 nlwork=1 311 allocate( point2ilwork(1)) 312 nshg0 = nshg 313 endif 314 315cccccccccccccccccccccccccccccccccccccccccc 316 317 itwo=2 318 write (temp1,"('(''co-ordinates@'',i',i1,',A1)')") itmp 319 write (fname2,temp1) (myrank+1),'?' 320 321c print *, "fname2 is", fname2 322 323cc MR CHANGE 324c call initphmpiio(nfields,nppf,nfiles,igeom,'read') 325c call openfile( fnamer, 'read', igeom ) 326CC MR CHANGE 327 328 call readheader(igeom,fname2 // char(0),intfromfile,itwo, 329 & 'double' // char(0), iotype) 330 numnp=intfromfile(1) 331c print *, "read out @@@@@@ is ", numnp 332 allocate( point2x(numnp,nsd) ) 333 allocate( xread(numnp,nsd) ) 334 ixsiz=numnp*nsd 335 call readdatablock(igeom,fname2 // char(0),xread,ixsiz, 336 & 'double' // char(0), iotype) 337 point2x = xread 338 339! call closefile( igeom, "read" ) 340! call finalizephmpiio( igeom ) 341 342! print *, "Finished finalize" 343 344c deallocate (point2x) 345c deallocate (xread) 346 347cccccccccccccccccccccccccccccccccccccccccc 348 349c 350c.... read the node coordinates 351c 352 353c itwo=2 354c fname1='co-ordinates?' 355c call readheader(igeomBAK,fname1,intfromfile,itwo, 'double', iotype) 356c numnp=intfromfile(1) 357cc nsd=intfromfile(2) 358c allocate( point2x(numnp,nsd) ) 359c allocate( xread(numnp,nsd) ) 360c ixsiz=numnp*nsd 361c call readdatablock(igeomBAK,fname1,xread,ixsiz, 'double',iotype) 362c point2x = xread 363 364c 365c.... read in and block out the connectivity 366c 367 368! !MR CHANGE 369! This is not necessary but this avoids to have the geombc files opend two times. 370! A better way consists in pasisng the file handler to genblk or make it global or use igeomBAK instead of igeom 371! call closefile( igeom, "read" ) 372! call finalizephmpiio( igeom ) 373! !MR CHANGE END 374 375 call genblk (IBKSIZ) 376 377! !MR CHANGE 378! call initphmpiio( nfields, nppf, nfiles, igeom, 'read') 379! call openfile( fnamer, 'read', igeom ) 380! !MR CHANGE END 381 382c 383c.... read the boundary condition mapping array 384c 385 386cc MR CHANGE 387! call initphmpiio(nfields,nppf,nfiles,igeom, 'read') 388! call openfile( fnamer, 'read', igeom ) 389cc MR CHANGE 390 391 ione=1 392 write (temp1,"('(''bc mapping array@'',i',i1,',A1)')") itmp 393 write (fname2,temp1) (myrank+1),'?' 394 call readheader(igeom,fname2 // char(0),nshg,ione, 395 & 'integer' // char(0),iotype) 396 397c fname1='bc mapping array?' 398c call readheader(igeomBAK,fname1,nshg, 399c & ione,'integer', iotype) 400 401 allocate( nBC(nshg) ) 402 403 allocate( nBCread(nshg) ) 404 405c call readdatablock(igeomBAK,fname1,nBCread,nshg,'integer',iotype) 406 call readdatablock(igeom,fname2 // char(0),nBCread,nshg, 407 & 'integer' // char(0),iotype) 408 409 nBC=nBCread 410 411c 412c.... read the temporary iBC array 413c 414 ione=1 415 write (temp1,"('(''bc codes array@'',i',i1,',A1)')") itmp 416 write (fname2,temp1) (myrank+1),'?' 417 call readheader(igeom,fname2 // char(0) ,numpbc,ione, 418 & 'integer' // char(0),iotype) 419 420c ione = 1 421c fname1='bc codes array?' 422c call readheader(igeomBAK,fname1,numpbc, 423c & ione, 'integer', iotype) 424 425!MR CHANGE 426! if ( numpbc > 0 ) then 427! allocate( iBCtmp(numpbc) ) 428! allocate( iBCtmpread(numpbc) ) 429! c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 430! c & 'integer',iotype) 431! call readdatablock(igeom,fname2,iBCtmpread,numpbc, 432! & 'integer',iotype) 433! iBCtmp=iBCtmpread 434! else ! sometimes a partition has no BC's 435! allocate( iBCtmp(1) ) 436! iBCtmp=0 437! endif 438 439 if ( numpbc > 0 ) then 440 allocate( iBCtmp(numpbc) ) 441 allocate( iBCtmpread(numpbc) ) 442 else 443 allocate( iBCtmp(1) ) 444 allocate( iBCtmpread(1) ) 445 endif 446c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 447c & 'integer',iotype) 448 call readdatablock(igeom,fname2 // char(0),iBCtmpread,numpbc, 449 & 'integer' // char(0),iotype) 450 451 if ( numpbc > 0 ) then 452 iBCtmp=iBCtmpread 453 else ! sometimes a partition has no BC's 454 deallocate( iBCtmpread) 455 iBCtmp=0 456 endif 457!MR CHANGE END 458 459c 460c.... read boundary condition data 461c 462 463 ione=1 464 write (temp1,"('(''boundary condition array@'',i',i1,',A1)')") 465 & itmp 466 write (fname2,temp1) (myrank+1),'?' 467 468c ione=1 469c fname1='boundary condition array?' 470c call readheader(igeomBAK,fname1,intfromfile, 471c & ione, 'double', iotype) 472 call readheader(igeom,fname2 // char(0),intfromfile, 473 & ione, 'double' // char(0), iotype) 474c here intfromfile(1) contains (ndof+7)*numpbc 475!MR CHANGE 476! if ( numpbc > 0 ) then 477! if(intfromfile(1).ne.(ndof+7)*numpbc) then 478! warning='WARNING more data in BCinp than needed: keeping 1st' 479! write(*,*) warning, ndof+7 480! endif 481! allocate( BCinp(numpbc,ndof+7) ) 482! nsecondrank=intfromfile(1)/numpbc 483! allocate( BCinpread(numpbc,nsecondrank) ) 484! iBCinpsiz=intfromfile(1) 485! c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 486! c & 'double',iotype) 487! call readdatablock(igeom,fname2,BCinpread,iBCinpsiz, 488! & 'double',iotype) 489! BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 490! else ! sometimes a partition has no BC's 491! allocate( BCinp(1,ndof+7) ) 492! BCinp=0 493! endif 494 495 if ( numpbc > 0 ) then 496! if(intfromfile(1).ne.(ndof+7)*numpbc) then 497! warning='WARNING more data in BCinp than needed: keeping 1st' 498! write(*,*) warning, ndof+7 499! endif 500 allocate( BCinp(numpbc,ndof+7) ) 501 nsecondrank=intfromfile(1)/numpbc 502 allocate( BCinpread(numpbc,nsecondrank) ) 503 iBCinpsiz=intfromfile(1) 504 else 505 allocate( BCinp(1,ndof+7) ) 506 allocate( BCinpread(0,0) ) !dummy 507 iBCinpsiz=intfromfile(1) 508 endif 509c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 510c & 'double',iotype) 511 512 call readdatablock(igeom,fname2 // char(0),BCinpread,iBCinpsiz, 513 & 'double' // char(0) ,iotype) 514 515 if ( numpbc > 0 ) then 516 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 517 else ! sometimes a partition has no BC's 518 deallocate(BCinpread) 519 BCinp=0 520 endif 521!MR CHANGE END 522 523c 524c.... read periodic boundary conditions 525c 526 527 ione=1 528 write (temp1,"('(''periodic masters array@'',i',i1,',A1)')") itmp 529 write (fname2,temp1) (myrank+1),'?' 530 531c fname1='periodic masters array?' 532c call readheader(igeomBAK,fname1,nshg, 533c & ione, 'integer', iotype) 534 call readheader(igeom,fname2 // char(0) ,nshg, 535 & ione, 'integer' // char(0), iotype) 536 allocate( point2iper(nshg) ) 537 allocate( iperread(nshg) ) 538c call readdatablock(igeomBAK,fname1,iperread,nshg, 539c & 'integer',iotype) 540 call readdatablock(igeom,fname2 // char(0),iperread,nshg, 541 & 'integer' // char(0),iotype) 542 point2iper=iperread 543 544 545! !MR CHANGE 546! call closefile( igeom, "read" ) 547! call finalizephmpiio( igeom ) 548! !MR CHANGE END 549 550c 551c.... generate the boundary element blocks 552c 553 call genbkb (ibksiz) 554 555 556! !MR CHANGE 557! write(*,*) 'HELLO 12 from ', myrank 558! !MR CHANGE END 559 560c 561c Read in the nsons and ifath arrays if needed 562c 563c There is a fundamental shift in the meaning of ifath based on whether 564c there exist homogenous directions in the flow. 565c 566c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 567c points in the TOTAL mesh. That is to say that each partition keeps a 568c link to ALL inhomogenous points. This link is furthermore not to the 569c sms numbering but to the original structured grid numbering. These 570c inhomogenous points are thought of as fathers, with their sons being all 571c the points in the homogenous directions that have this father's 572c inhomogeneity. The array ifath takes as an arguement the sms numbering 573c and returns as a result the father. 574c 575c In this case nsons is the number of sons that each father has and ifath 576c is an array which tells the 577c 578c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 579c if we followed the above strategy since every partition would index its 580c points to the ENTIRE mesh. Furthermore, there would never be a need 581c to average to a node off processor since there is no spatial averaging. 582c Therefore, to properly account for this case we must recognize it and 583c inerrupt certain actions (i.e. assembly of the average across partitions). 584c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 585c father has any sons). Reiterating to be clear, in this case ifath does 586c not point to a global numbering but instead just points to itself. 587c 588 nfath=1 ! some architectures choke on a zero or undeclared 589 ! dimension variable. This sets it to a safe, small value. 590 if(((iLES .lt. 20) .and. (iLES.gt.0)) 591 & .or. (itwmod.gt.0) ) then ! don't forget same 592 ! conditional in proces.f 593 594c read (igeomBAK) nfath ! nfath already read in input.f, 595 ! needed for alloc 596 ione=1 597c call creadlist(igeomBAK,ione,nfath) 598c fname1='keyword sonfath?' 599 if(nohomog.gt.0) then 600 601 602! fname1='number of father-nodes?' 603! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 604 605 write (temp1,"('(''number of father-nodes@'',i',i1,',A1)')") 606 & itmp 607 write (fname2,temp1) (myrank+1),'?' 608 call readheader(igeom,fname2 // char(0),nfath,ione, 609 & 'integer' // char(0), iotype) 610 611c 612c fname1='keyword nsons?' 613! fname1='number of son-nodes for each father?' 614! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 615 616 write (temp1, 617 & "('(''number of son-nodes for each father@'',i',i1,',A1)')") itmp 618 write (fname2,temp1) (myrank+1),'?' 619 call readheader(igeom,fname2 // char(0),nfath,ione, 620 & 'integer' // char(0), iotype) 621 622 allocate (point2nsons(nfath)) 623 624! call readdatablock(igeomBAK,fname1,point2nsons,nfath, 625! & 'integer',iotype) 626 call readdatablock(igeom,fname2 // char(0),point2nsons, 627 & nfath,'integer' // char(0), iotype) 628 629c 630! fname1='keyword ifath?' 631! call readheader(igeomBAK,fname1,nshg,ione,'integer', iotype) 632 633 write (temp1,"('(''keyword ifath@'',i',i1,',A1)')") itmp 634 write (fname2,temp1) (myrank+1),'?' 635 call readheader(igeom,fname2 // char(0),nshg,ione, 636 & 'integer' // char(0), iotype) 637 638 allocate (point2ifath(nshg)) 639 640! call readdatablock(igeomBAK,fname1,point2ifath,nshg, 641! & 'integer',iotype) 642 call readdatablock(igeom,fname2 // char(0),point2ifath, 643 & nshg,'integer' // char(0) , iotype) 644 645c 646 nsonmax=maxval(point2nsons) 647c 648 else ! this is the case where there is no homogeneity 649 ! therefore ever node is a father (too itself). sonfath 650 ! (a routine in NSpre) will set this up but this gives 651 ! you an option to avoid that. 652 nfath=nshg 653 allocate (point2nsons(nfath)) 654 point2nsons=1 655 allocate (point2ifath(nshg)) 656 do i=1,nshg 657 point2ifath(i)=i 658 enddo 659 nsonmax=1 660c 661 endif 662 else 663 allocate (point2nsons(1)) 664 allocate (point2ifath(1)) 665 endif 666 667! !MR CHANGE 668 call closefile( igeom, "read" // char(0) ) 669 call finalizephmpiio( igeom ) 670! !MR CHANGE END 671 672! !MR CHANGE 673! write(*,*) 'HELLO 13 from ', myrank 674! !MR CHANGE END 675 676c 677c renumber the master partition for SPEBC 678c 679c if((myrank.eq.master).and.(irscale.ge.0)) then 680c call setSPEBC(numnp, nfath, nsonmax) 681c call renum(point2x,point2ifath,point2nsons) 682c endif 683c 684c.... Read restart files 685 686c$$$ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 687c 688c nfields = 11 689c numparts = 512 690c nppp = numparts/numpe 691c startpart = myrank * nppp +1 692c int endpart = startpart + nppp - 1 693c nppf = numparts/nfiles 694cc fnamer = "/users/nliu/PIG4/4-procs_case/restart-file" 695cc fnamer="./4-procs_case/restart-file" 696 697 itmp=1 698 if (irstart .gt. 0) itmp = int(log10(float(irstart+1)))+1 699 700 write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp 701 702 write (fnamer,fmt1) irstart 703 fnamer = trim(fnamer) // cname2(color+1) 704 705 call queryphmpiio(fnamer // char(0), nfields, nppf); 706 if (myrank == 0) then 707 write(*,*) 'Number of fields in restart-dat: ',nfields 708 write(*,*) 'Number of parts per file restart-dat: ',nppf 709 endif 710 call initphmpiio(nfields,nppf,nfiles,descriptor, 711 & 'read' // char(0)) 712 call openfile( fnamer // char(0) , 713 & 'read' // char(0), descriptor ) 714 715 ithree=3 716c call creadlist(irstin,ithree,nshg2,ndof2,lstep) 717 718 itmp = int(log10(float(myrank+1)))+1 719 write (temp1,"('(''solution@'',i',i1,',A1)')") itmp 720 write (fname1,temp1) (myrank+1),'?' 721 fname1 = trim(fname1) 722 723c print *, "Solution is : ", fname1 724 725 intfromfile=0 726 call readheader(descriptor,fname1 // char(0) ,intfromfile, 727 & ithree,'integer' // char(0), 728 & iotype) 729c 730c.... read the values of primitive variables into q 731c 732 733c print *, "intfromfile(1) is ", intfromfile(1) 734c print *, "intfromfile(2) is ", intfromfile(2) 735c print *, "intfromfile(3) is ", intfromfile(3) 736 737 allocate( qold(nshg,ndof) ) 738 if(intfromfile(1).ne.0) then 739 nshg2=intfromfile(1) 740 ndof2=intfromfile(2) 741 lstep=intfromfile(3) 742 if(ndof2.ne.ndof) then 743 744 endif 745c 746 if (nshg2 .ne. nshg) 747 & call error ('restar ', 'nshg ', nshg) 748 allocate( qread(nshg,ndof2) ) 749 iqsiz=nshg*ndof2 750 call readdatablock(descriptor,fname1 // char(0),qread,iqsiz, 751 & 'double' // char(0),iotype) 752 qold(:,1:ndof)=qread(:,1:ndof) 753 deallocate(qread) 754 else 755 if (myrank.eq.master) then 756 if (matflg(1,1).eq.0) then ! compressible 757 warning='Solution is set to zero (with p and T to one)' 758 else 759 warning='Solution is set to zero' 760 endif 761 write(*,*) warning 762 endif 763 qold=zero 764 if (matflg(1,1).eq.0) then ! compressible 765 qold(:,1)=one ! avoid zero pressure 766 qold(:,nflow)=one ! avoid zero temperature 767 endif 768 endif 769 770 771! !MR CHANGE 772! write(*,*) 'HELLO 16-8 from ', myrank 773! !MR CHANGE END 774 775c itmp=1 776c if (myrank .gt. 0) itmp = int(log10(float(myrank)))+1 777 write (temp1,"('(''time derivative of solution@'',i',i1,',A1)')") 778 & itmp 779 write (fname1,temp1) (myrank+1),'?' 780 fname1 = trim(fname1) 781 intfromfile=0 782 call readheader(descriptor,fname1 // char(0) ,intfromfile, 783 & ithree,'integer' // char(0), 784 & iotype) 785 allocate( acold(nshg,ndof) ) 786 if(intfromfile(1).ne.0) then 787 nshg2=intfromfile(1) 788 ndof2=intfromfile(2) 789 lstep=intfromfile(3) 790 791c print *, "intfromfile(1) is ", intfromfile(1) 792c print *, "intfromfile(2) is ", intfromfile(2) 793c print *, "intfromfile(3) is ", intfromfile(3) 794 795 if (nshg2 .ne. nshg) 796 & call error ('restar ', 'nshg ', nshg) 797c 798 allocate( acread(nshg,ndof2) ) 799 acread=zero 800 iacsiz=nshg*ndof2 801 call readdatablock(descriptor,fname1 // char(0),acread, 802 & iacsiz, 'double' // char(0),iotype) 803 acold(:,1:ndof)=acread(:,1:ndof) 804 deallocate(acread) 805 else 806 if (myrank.eq.master) then 807 warning='Time derivative of solution is set to zero (SAFE)' 808 write(*,*) warning 809 endif 810 acold=zero 811 endif 812 813cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 814cc 815c 816cc 817cc.... read the header and check it against the run data 818cc 819 820 821c ithree=3 822ccc call creadlist(irstin,ithree,nshg2,ndof2,lstep) 823c fname1='solution?' 824c intfromfile=0 825c call readheader(irstin,fname1,intfromfile, 826c & ithree,'integer', iotype) 827cc 828cc.... read the values of primitive variables into q 829cc 830c allocate( qold(nshg,ndof) ) 831c if(intfromfile(1).ne.0) then 832c nshg2=intfromfile(1) 833c ndof2=intfromfile(2) 834c lstep=intfromfile(3) 835c if(ndof2.ne.ndof) then 836c warning='WARNING more data in restart than needed: keeping 1st ' 837c write(*,*) warning , ndof 838c endif 839cc 840c if (nshg2 .ne. nshg) 841c & call error ('restar ', 'nshg ', nshg) 842c allocate( qread(nshg,ndof2) ) 843 844c iqsiz=nshg*ndof2 845c call readdatablock(irstin,fname1,qread,iqsiz, 846c & 'double',iotype) 847c qold(:,1:ndof)=qread(:,1:ndof) 848c deallocate(qread) 849c else 850c if (myrank.eq.master) then 851c if (matflg(1,1).eq.0) then ! compressible 852c warning='Solution is set to zero (with p and T to one)' 853c else 854c warning='Solution is set to zero' 855c endif 856c write(*,*) warning 857c endif 858c qold=zero 859c if (matflg(1,1).eq.0) then ! compressible 860c qold(:,1)=one ! avoid zero pressure 861c qold(:,nflow)=one ! avoid zero temperature 862c endif 863c endif 864cc 865c fname1='time derivative of solution?' 866c intfromfile=0 867c call readheader(irstin,fname1,intfromfile, 868c & ithree,'integer', iotype) 869c allocate( acold(nshg,ndof) ) 870c if(intfromfile(1).ne.0) then 871c nshg2=intfromfile(1) 872c ndof2=intfromfile(2) 873c lstep=intfromfile(3) 874c 875c if (nshg2 .ne. nshg) 876c & call error ('restar ', 'nshg ', nshg) 877cc 878c allocate( acread(nshg,ndof2) ) 879c acread=zero 880c 881c iacsiz=nshg*ndof2 882c call readdatablock(irstin,fname1,acread,iacsiz, 883c & 'double',iotype) 884c acold(:,1:ndof)=acread(:,1:ndof) 885c deallocate(acread) 886c else 887c if (myrank.eq.master) then 888c warning='Time derivative of solution is set to zero (SAFE)' 889c write(*,*) warning 890c endif 891c acold=zero 892c endif 893 894c call creadlist(irstin,ithree,nshg2,ndisp,lstep) 895 896 if (ideformwall.eq.1) then 897! fname1='displacement?' 898! call readheader(irstin,fname1,intfromfile, 899! & ithree,'integer', iotype) 900 901 write (temp1,"('(''displacement@'',i',i1,',A1)')") 902 & itmp 903 write (fname1,temp1) (myrank+1),'?' 904 fname1 = trim(fname1) 905 intfromfile=0 906 call readheader(descriptor,fname1 // char(0),intfromfile, 907 & ithree,'integer' // char(0),iotype) 908 909 nshg2=intfromfile(1) 910 ndisp=intfromfile(2) 911 lstep=intfromfile(3) 912 if(ndisp.ne.nsd) then 913 warning='WARNING ndisp not equal nsd' 914 write(*,*) warning , ndisp 915 endif 916c 917 if (nshg2 .ne. nshg) 918 & call error ('restar ', 'nshg ', nshg) 919c 920c.... read the values of primitive variables into uold 921c 922 923 allocate( uold(nshg,nsd) ) 924 allocate( uread(nshg,nsd) ) 925 926 iusiz=nshg*nsd 927 928! call readdatablock(irstin,fname1,uread,iusiz, 929! & 'double',iotype) 930 call readdatablock(descriptor,fname1 // char(0) ,uread,iusiz, 931 & 'double' // char(0),iotype) 932 933 uold(:,1:nsd)=uread(:,1:nsd) 934 else 935 allocate( uold(nshg,nsd) ) 936 uold(:,1:nsd) = zero 937 endif 938 939c 940c.... close c-binary files 941c 942!MR CHANGE 943! call closefile( irstin, "read" ) 944 945 call closefile( descriptor, "read" // char(0) ) 946 call finalizephmpiio( descriptor ) 947 948!MR CHANGE 949! call closefile( igeomBAK, "read" ) 950c 951 deallocate(xread) 952 if ( numpbc > 0 ) then 953 deallocate(bcinpread) 954 deallocate(ibctmpread) 955 endif 956 deallocate(iperread) 957 if(numpe.gt.1) 958 & deallocate(ilworkread) 959 deallocate(nbcread) 960 961 return 962c 963 994 call error ('input ','opening ', igeomBAK) 964 995 call error ('input ','opening ', igeomBAK) 965 997 call error ('input ','end file', igeomBAK) 966 998 call error ('input ','end file', igeomBAK) 967c 968 end 969 970c 971c No longer called but kept around in case.... 972c 973 subroutine genpzero(iBC) 974 975 use pointer_data 976c 977 include "common.h" 978 integer iBC(nshg) 979c 980c.... check to see if any of the nodes have a dirichlet pressure 981c 982 pzero=1 983 if (any(btest(iBC,2))) pzero=0 984c 985 do iblk = 1, nelblb 986 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 987 do i=1, npro 988 iBCB1=miBCB(iblk)%p(i,1) 989c 990c.... check to see if any of the nodes have a Neumann pressure 991c but not periodic (note that 992c 993 if(btest(iBCB1,1)) pzero=0 994 enddo 995c 996c.... share results with other processors 997c 998 pzl=pzero 999 if (numpe .gt. 1) 1000 & call MPI_ALLREDUCE (pzl, pzero, 1, 1001 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 1002 1003 enddo 1004c 1005c.... return 1006c 1007 return 1008c 1009 end 1010 1011