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