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