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 module readarrays 12 13 real*8, allocatable :: point2x(:,:) 14 real*8, allocatable :: qold(:,:) 15 real*8, allocatable :: uold(:,:) 16 real*8, allocatable :: acold(:,:) 17 integer, allocatable :: iBCtmp(:) 18 real*8, allocatable :: BCinp(:,:) 19 20 integer, allocatable :: point2ilwork(:) 21! integer, allocatable :: fncorp(:) 22 integer, allocatable :: twodncorp(:,:) 23 integer, allocatable :: nBC(:) 24 integer, allocatable :: point2iper(:) 25 integer, target, allocatable :: point2ifath(:) 26 integer, target, allocatable :: point2nsons(:) 27 28 end module 29 30 subroutine readnblk 31 use iso_c_binding 32 use readarrays 33 use fncorpmod 34 use phio 35 use phstr 36 use syncio 37 use posixio 38 use streamio 39 include "common.h" 40 41 real*8, target, allocatable :: xread(:,:), qread(:,:), acread(:,:) 42 real*8, target, allocatable :: uread(:,:) 43 real*8, target, allocatable :: BCinpread(:,:) 44 real*8 :: iotime 45 integer, target, allocatable :: iperread(:), iBCtmpread(:) 46 integer, target, allocatable :: ilworkread(:), nBCread(:) 47 integer, target, allocatable :: fncorpread(:) 48 integer fncorpsize 49 character*10 cname2, cname2nd 50 character*8 mach2 51 character*30 fmt1 52 character*255 fname1,fnamer,fnamelr 53 character*255 warning 54 integer igeomBAK, ibndc, irstin, ierr 55 integer, target :: intfromfile(50) ! integers read from headers 56 integer :: descriptor, descriptorG, GPID, color, nfields 57 integer :: numparts, nppf 58 integer :: ierr_io, numprocs, itmp, itmp2 59 integer :: ignored 60 integer :: fileFmt 61 character*255 fname2, temp2 62 character*64 temp1 63 type(c_ptr) :: handle 64 character(len=1024) :: dataInt, dataDbl 65 dataInt = c_char_'integer'//c_null_char 66 dataDbl = c_char_'double'//c_null_char 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 75 fname1='geombc.dat' 76 fname1= trim(fname1) // cname2(myrank+1) 77 78 itmp=1 79 if (irstart .gt. 0) itmp = int(log10(float(irstart)))+1 80 write (fmt1,"('(''restart.'',i',i1,',1x)')") itmp 81 write (fnamer,fmt1) irstart 82 fnamer = trim(fnamer) // cname2(myrank+1) 83c 84c.... open input files 85c.... input the geometry parameters 86c 87 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 88 89 itwo=2 90 ione=1 91 ieleven=11 92 itmp = int(log10(float(myrank+1)))+1 93 94 iotime = TMRC() 95 if( input_mode .eq. -1 ) then 96 call streamio_setup_read(fhandle, geomRestartStream) 97 else if( input_mode .eq. 0 ) then 98 call posixio_setup(fhandle, c_char_'r') 99 else if( input_mode .ge. 1 ) then 100 call syncio_setup_read(nsynciofiles, fhandle) 101 end if 102 call phio_constructName(fhandle, 103 & c_char_'geombc' // char(0), fname1) 104 call phio_openfile(fname1, fhandle); 105 106 call phio_readheader(fhandle,c_char_'number of nodes' // char(0), 107 & c_loc(numnp),ione, dataInt, iotype) 108 109 call phio_readheader(fhandle,c_char_'number of modes' // char(0), 110 & c_loc(nshg),ione, dataInt, iotype) 111 112 call phio_readheader(fhandle, 113 & c_char_'number of interior elements' // char(0), 114 & c_loc(numel),ione, dataInt, iotype) 115 116 call phio_readheader(fhandle, 117 & c_char_'number of boundary elements' // char(0), 118 & c_loc(numelb),ione, dataInt, iotype) 119 120 call phio_readheader(fhandle, 121 & c_char_'maximum number of element nodes' // char(0), 122 & c_loc(nen),ione, dataInt, iotype) 123 124 call phio_readheader(fhandle, 125 & c_char_'number of interior tpblocks' // char(0), 126 & c_loc(nelblk),ione, dataInt, iotype) 127 128 call phio_readheader(fhandle, 129 & c_char_'number of boundary tpblocks' // char(0), 130 & c_loc(nelblb),ione, dataInt, iotype) 131 132 call phio_readheader(fhandle, 133 & c_char_'number of nodes with Dirichlet BCs' // char(0), 134 & c_loc(numpbc),ione, dataInt, iotype) 135 136 call phio_readheader(fhandle, 137 & c_char_'number of shape functions' // char(0), 138 & c_loc(ntopsh),ione, dataInt, iotype) 139c 140c.... calculate the maximum number of boundary element nodes 141c 142 nenb = 0 143 do i = 1, melCat 144 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 145 enddo 146 if (myrank == master) then 147 if (nenb .eq. 0) call error ('input ','nen ',nen) 148 endif 149c 150c.... setup some useful constants 151c 152 I3nsd = nsd / 3 ! nsd=3 integer flag 153 E3nsd = float(I3nsd) ! nsd=3 real flag 154 if(matflg(1,1).lt.0) then 155 nflow = nsd + 1 156 else 157 nflow = nsd + 2 158 endif 159 ndof = nsd + 2 160 nsclr=impl(1)/100 161 ndof=ndof+nsclr ! number of sclr transport equations to solve 162 163 ndofBC = ndof + I3nsd ! dimension of BC array 164 ndiBCB = 2 ! dimension of iBCB array 165 ndBCB = ndof + 1 ! dimension of BCB array 166 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 167c 168c.... ----------------------> Communication tasks <-------------------- 169c 170 if(numpe > 1) then 171 call phio_readheader(fhandle, 172 & c_char_'size of ilwork array' // char(0), 173 & c_loc(nlwork),ione, dataInt, iotype) 174 175 call phio_readheader(fhandle, 176 & c_char_'ilwork' //char(0), 177 & c_loc(nlwork),ione, dataInt, iotype) 178 179 allocate( point2ilwork(nlwork) ) 180 allocate( ilworkread(nlwork) ) 181 call phio_readdatablock(fhandle, c_char_'ilwork' // char(0), 182 & c_loc(ilworkread), nlwork, dataInt , iotype) 183 184 point2ilwork = ilworkread 185 call ctypes (point2ilwork) 186 187 if((usingPETSc.eq.1).or.(svLSFlag.eq.1)) then 188 fncorpsize = nshg 189 allocate(fncorp(fncorpsize)) 190 call gen_ncorp(fncorp, ilworkread, nlwork, fncorpsize) 191! 192! the following code finds the global range of the owned nodes 193! 194 maxowned=0 195 minowned=maxval(fncorp) 196 do i = 1,nshg 197 if(fncorp(i).gt.0) then ! don't consider remote copies 198 maxowned=max(maxowned,fncorp(i)) 199 minowned=min(minowned,fncorp(i)) 200 endif 201 enddo 202! 203! end of global range code 204! 205 call commuInt(fncorp, point2ilwork, 1, 'out') 206 ncorpsize = fncorpsize 207 endif 208 if(svLSFlag.eq.1) then 209 allocate(ltg(ncorpsize)) 210 ltg(:)=fncorp(:) 211 endif 212 else 213 allocate(ltg(nshg)) 214 do i =1,nshg 215 ltg(i)=i 216 enddo 217 nlwork=1 218 allocate( point2ilwork(1)) 219 nshg0 = nshg 220 endif 221 222 223 itwo=2 224 225 call phio_readheader(fhandle, 226 & c_char_'co-ordinates' // char(0), 227 & c_loc(intfromfile),itwo, dataDbl, iotype) 228 numnp=intfromfile(1) 229 allocate( point2x(numnp,nsd) ) 230 allocate( xread(numnp,nsd) ) 231 ixsiz=numnp*nsd 232 call phio_readdatablock(fhandle, 233 & c_char_'co-ordinates' // char(0), 234 & c_loc(xread),ixsiz, dataDbl, iotype) 235 point2x = xread 236 237c..............................for Duct 238 if(istretchOutlet.eq.1)then 239 240c...geometry6 241 if(iDuctgeometryType .eq. 6) then 242 xmaxn = 1.276 243 xmaxo = 0.848 244 xmin = 0.42 245c...geometry8 246 elseif(iDuctgeometryType .eq. 8)then 247 xmaxn=1.6*4.5*0.0254+0.85*1.5 248 xmaxo=1.6*4.5*0.0254+0.85*1.0 249 xmin =1.6*4.5*0.0254+0.85*0.5 250 endif 251c... 252 alpha=(xmaxn-xmaxo)/(xmaxo-xmin)**2 253 where (point2x(:,1) .ge. xmin) 254c..... N=# of current elements from .42 to exit(~40) 255c..... (x_mx-x_mn)/N=.025 256c..... alpha=3 3*.025=.075 257 point2x(:,1)=point2x(:,1)+ 258 & alpha*(point2x(:,1)-xmin)**2 259c..... ftn to stretch x at exit 260 endwhere 261 endif 262 263c 264c.... read in and block out the connectivity 265c 266 call genblk (IBKSIZ) 267c 268c.... read the boundary condition mapping array 269c 270 ione=1 271 call phio_readheader(fhandle, 272 & c_char_'bc mapping array' // char(0), 273 & c_loc(nshg),ione, dataInt, iotype) 274 275 allocate( nBC(nshg) ) 276 277 allocate( nBCread(nshg) ) 278 279 call phio_readdatablock(fhandle, 280 & c_char_'bc mapping array' // char(0), 281 & c_loc(nBCread), nshg, dataInt, iotype) 282 283 nBC=nBCread 284c 285c.... read the temporary iBC array 286c 287 ione=1 288 call phio_readheader(fhandle, 289 & c_char_'bc codes array' // char(0), 290 & c_loc(numpbc),ione, dataInt, iotype) 291 292 if ( numpbc > 0 ) then 293 allocate( iBCtmp(numpbc) ) 294 allocate( iBCtmpread(numpbc) ) 295 else 296 allocate( iBCtmp(1) ) 297 allocate( iBCtmpread(1) ) 298 endif 299 call phio_readdatablock(fhandle, 300 & c_char_'bc codes array' // char(0), 301 & c_loc(iBCtmpread), numpbc, dataInt, iotype) 302 303 if ( numpbc > 0 ) then 304 iBCtmp=iBCtmpread 305 else ! sometimes a partition has no BC's 306 deallocate( iBCtmpread) 307 iBCtmp=0 308 endif 309c 310c.... read boundary condition data 311c 312 ione=1 313 call phio_readheader(fhandle, 314 & c_char_'boundary condition array' // char(0), 315 & c_loc(intfromfile),ione, dataDbl, iotype) 316 317 if ( numpbc > 0 ) then 318 allocate( BCinp(numpbc,ndof+7) ) 319 nsecondrank=intfromfile(1)/numpbc 320 allocate( BCinpread(numpbc,nsecondrank) ) 321 iBCinpsiz=intfromfile(1) 322 else 323 allocate( BCinp(1,ndof+7) ) 324 allocate( BCinpread(0,0) ) !dummy 325 iBCinpsiz=intfromfile(1) 326 endif 327 328 call phio_readdatablock(fhandle, 329 & c_char_'boundary condition array' // char(0), 330 & c_loc(BCinpread), iBCinpsiz, dataDbl, iotype) 331 332 if ( numpbc > 0 ) then 333 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 334 else ! sometimes a partition has no BC's 335 deallocate(BCinpread) 336 BCinp=0 337 endif 338c 339c.... read periodic boundary conditions 340c 341 ione=1 342 call phio_readheader(fhandle, 343 & c_char_'periodic masters array' // char(0), 344 & c_loc(nshg), ione, dataInt, iotype) 345 allocate( point2iper(nshg) ) 346 allocate( iperread(nshg) ) 347 call phio_readdatablock(fhandle, 348 & c_char_'periodic masters array' // char(0), 349 & c_loc(iperread), nshg, dataInt, iotype) 350 point2iper=iperread 351c 352c.... generate the boundary element blocks 353c 354 call genbkb (ibksiz) 355c 356c Read in the nsons and ifath arrays if needed 357c 358c There is a fundamental shift in the meaning of ifath based on whether 359c there exist homogenous directions in the flow. 360c 361c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 362c points in the TOTAL mesh. That is to say that each partition keeps a 363c link to ALL inhomogenous points. This link is furthermore not to the 364c sms numbering but to the original structured grid numbering. These 365c inhomogenous points are thought of as fathers, with their sons being all 366c the points in the homogenous directions that have this father's 367c inhomogeneity. The array ifath takes as an arguement the sms numbering 368c and returns as a result the father. 369c 370c In this case nsons is the number of sons that each father has and ifath 371c is an array which tells the 372c 373c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 374c if we followed the above strategy since every partition would index its 375c points to the ENTIRE mesh. Furthermore, there would never be a need 376c to average to a node off processor since there is no spatial averaging. 377c Therefore, to properly account for this case we must recognize it and 378c inerrupt certain actions (i.e. assembly of the average across partitions). 379c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 380c father has any sons). Reiterating to be clear, in this case ifath does 381c not point to a global numbering but instead just points to itself. 382c 383 nfath=1 ! some architectures choke on a zero or undeclared 384 ! dimension variable. This sets it to a safe, small value. 385 if(((iLES .lt. 20) .and. (iLES.gt.0)) 386 & .or. (itwmod.gt.0) ) then ! don't forget same 387 ! conditional in proces.f 388 ! needed for alloc 389 ione=1 390 if(nohomog.gt.0) then 391 call phio_readheader(fhandle, 392 & c_char_'number of father-nodes' // char(0), 393 & c_loc(nfath), ione, dataInt, iotype) 394 395 call phio_readheader(fhandle, 396 & c_char_'number of son-nodes for each father' // char(0), 397 & c_loc(nfath), ione, dataInt, iotype) 398 399 allocate (point2nsons(nfath)) 400 401 call phio_readdatablock(fhandle, 402 & c_char_'number of son-nodes for each father' // char(0), 403 & c_loc(point2nsons),nfath, dataInt, iotype) 404 405 call phio_readheader(fhandle, 406 & c_char_'keyword ifath' // char(0), 407 & c_loc(nshg), ione, dataInt, iotype); 408 409 allocate (point2ifath(nshg)) 410 411 call phio_readdatablock(fhandle, 412 & c_char_'keyword ifath' // char(0), 413 & c_loc(point2ifath), nshg, dataInt, iotype) 414 415 nsonmax=maxval(point2nsons) 416 else ! this is the case where there is no homogeneity 417 ! therefore ever node is a father (too itself). sonfath 418 ! (a routine in NSpre) will set this up but this gives 419 ! you an option to avoid that. 420 nfath=nshg 421 allocate (point2nsons(nfath)) 422 point2nsons=1 423 allocate (point2ifath(nshg)) 424 do i=1,nshg 425 point2ifath(i)=i 426 enddo 427 nsonmax=1 428 endif 429 else 430 allocate (point2nsons(1)) 431 allocate (point2ifath(1)) 432 endif 433 434 call phio_closefile(fhandle); 435 iotime = TMRC() - iotime 436 if (myrank.eq.master) then 437 write(*,*) 'time to read geombc (seconds)', iotime 438 endif 439 440c.... Read restart files 441 iotime = TMRC() 442 if( input_mode .eq. -1 ) then 443 call streamio_setup_read(fhandle, geomRestartStream) 444 else if( input_mode .eq. 0 ) then 445 call posixio_setup(fhandle, c_char_'r') 446 else if( input_mode .ge. 1 ) then 447 call syncio_setup_read(nsynciofiles, fhandle) 448 end if 449 call phio_constructName(fhandle, 450 & c_char_'restart' // char(0), fnamer) 451 call phstr_appendInt(fnamer, irstart) 452 call phstr_appendStr(fnamer, c_char_'.'//c_null_char) 453 call phio_openfile(fnamer, fhandle); 454 455 ithree=3 456 457 itmp = int(log10(float(myrank+1)))+1 458 459 intfromfile=0 460 call phio_readheader(fhandle, 461 & c_char_'solution' // char(0), 462 & c_loc(intfromfile), ithree, dataInt, iotype) 463c 464c.... read the values of primitive variables into q 465c 466 allocate( qold(nshg,ndof) ) 467 if(intfromfile(1).ne.0) then 468 nshg2=intfromfile(1) 469 ndof2=intfromfile(2) 470 lstep=intfromfile(3) 471 if(ndof2.ne.ndof) then 472 473 endif 474 if (nshg2 .ne. nshg) 475 & call error ('restar ', 'nshg ', nshg) 476 allocate( qread(nshg,ndof2) ) 477 iqsiz=nshg*ndof2 478 call phio_readdatablock(fhandle, 479 & c_char_'solution' // char(0), 480 & c_loc(qread),iqsiz, dataDbl,iotype) 481 qold(:,1:ndof)=qread(:,1:ndof) 482 deallocate(qread) 483 else 484 if (myrank.eq.master) then 485 if (matflg(1,1).eq.0) then ! compressible 486 warning='Solution is set to zero (with p and T to one)' 487 else 488 warning='Solution is set to zero' 489 endif 490 write(*,*) warning 491 endif 492 qold=zero 493 if (matflg(1,1).eq.0) then ! compressible 494 qold(:,1)=one ! avoid zero pressure 495 qold(:,nflow)=one ! avoid zero temperature 496 endif 497 endif 498 499 intfromfile=0 500 call phio_readheader(fhandle, 501 & c_char_'time derivative of solution' // char(0), 502 & c_loc(intfromfile), ithree, dataInt, iotype) 503 allocate( acold(nshg,ndof) ) 504 if(intfromfile(1).ne.0) then 505 nshg2=intfromfile(1) 506 ndof2=intfromfile(2) 507 lstep=intfromfile(3) 508 509 if (nshg2 .ne. nshg) 510 & call error ('restar ', 'nshg ', nshg) 511 allocate( acread(nshg,ndof2) ) 512 acread=zero 513 iacsiz=nshg*ndof2 514 call phio_readdatablock(fhandle, 515 & c_char_'time derivative of solution' // char(0), 516 & c_loc(acread), iacsiz, dataDbl,iotype) 517 acold(:,1:ndof)=acread(:,1:ndof) 518 deallocate(acread) 519 else 520 if (myrank.eq.master) then 521 warning='Time derivative of solution is set to zero (SAFE)' 522 write(*,*) warning 523 endif 524 acold=zero 525 endif 526cc 527cc.... read the header and check it against the run data 528cc 529 if (ideformwall.eq.1) then 530 531 intfromfile=0 532 call phio_readheader(fhandle, 533 & c_char_'displacement' // char(0), 534 & c_loc(intfromfile), ithree, dataInt, iotype) 535 536 nshg2=intfromfile(1) 537 ndisp=intfromfile(2) 538 lstep=intfromfile(3) 539 if(ndisp.ne.nsd) then 540 warning='WARNING ndisp not equal nsd' 541 write(*,*) warning , ndisp 542 endif 543 if (nshg2 .ne. nshg) 544 & call error ('restar ', 'nshg ', nshg) 545c 546c.... read the values of primitive variables into uold 547c 548 allocate( uold(nshg,nsd) ) 549 allocate( uread(nshg,nsd) ) 550 551 iusiz=nshg*nsd 552 553 call phio_readdatablock(fhandle, 554 & c_char_'displacement' // char(0), 555 & c_loc(uread), iusiz, dataDbl, iotype) 556 557 uold(:,1:nsd)=uread(:,1:nsd) 558 else 559 allocate( uold(nshg,nsd) ) 560 uold(:,1:nsd) = zero 561 endif 562c 563c.... close c-binary files 564c 565 call phio_closefile(fhandle) 566 iotime = TMRC() - iotime 567 if (myrank.eq.master) then 568 write(*,*) 'time to read restart (seconds)', iotime 569 endif 570 571 deallocate(xread) 572 if ( numpbc > 0 ) then 573 deallocate(bcinpread) 574 deallocate(ibctmpread) 575 endif 576 deallocate(iperread) 577 if(numpe.gt.1) 578 & deallocate(ilworkread) 579 deallocate(nbcread) 580 581 return 582 994 call error ('input ','opening ', igeomBAK) 583 995 call error ('input ','opening ', igeomBAK) 584 997 call error ('input ','end file', igeomBAK) 585 998 call error ('input ','end file', igeomBAK) 586 end 587c 588c No longer called but kept around in case.... 589c 590 subroutine genpzero(iBC) 591 592 use pointer_data 593 include "common.h" 594 integer iBC(nshg) 595c 596c.... check to see if any of the nodes have a dirichlet pressure 597c 598 pzero=1 599 if (any(btest(iBC,2))) pzero=0 600 do iblk = 1, nelblb 601 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 602 do i=1, npro 603 iBCB1=miBCB(iblk)%p(i,1) 604c 605c.... check to see if any of the nodes have a Neumann pressure 606c but not periodic (note that 607c 608 if(btest(iBCB1,1)) pzero=0 609 enddo 610c 611c.... share results with other processors 612c 613 pzl=pzero 614 if (numpe .gt. 1) 615 & call MPI_ALLREDUCE (pzl, pzero, 1, 616 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 617 enddo 618 return 619 end 620