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