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