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