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