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