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