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