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 integer, target, allocatable :: iperread(:), iBCtmpread(:) 42 integer, target, allocatable :: ilworkread(:), nBCread(:) 43 character*10 cname2 44 character*8 mach2 45 character*30 fmt1 46 character*255 fname1,fnamer,fnamelr 47 character*255 warning 48 integer igeomBAK, ibndc, irstin, ierr 49 integer, target :: intfromfile(50) ! integers read from headers 50 integer :: descriptor, descriptorG, GPID, color, nfields 51 integer :: numparts, nppf 52 integer :: ierr_io, numprocs, itmp, itmp2 53 integer :: ignored 54 integer :: fileFmt 55 character*255 fname2, temp2 56 character*64 temp1 57 type(c_ptr) :: handle 58 character(len=1024) :: dataInt, dataDbl 59 dataInt = c_char_'integer'//c_null_char 60 dataDbl = c_char_'double'//c_null_char 61c 62c.... determine the step number to start with 63c 64 open(unit=72,file='numstart.dat',status='old') 65 read(72,*) irstart 66 close(72) 67 lstep=irstart ! in case restart files have no fields 68 69 fname1='geombc.dat' 70 fname1= trim(fname1) // cname2(myrank+1) 71 72 itmp=1 73 if (irstart .gt. 0) itmp = int(log10(float(irstart)))+1 74 write (fmt1,"('(''restart.'',i',i1,',1x)')") itmp 75 write (fnamer,fmt1) irstart 76 fnamer = trim(fnamer) // cname2(myrank+1) 77c 78c.... open input files 79c.... input the geometry parameters 80c 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( input_mode .eq. -1 ) then 89 call streamio_setup_read(fhandle, geomRestartStream) 90 else if( input_mode .eq. 0 ) then 91 call posixio_setup(fhandle, c_char_'r') 92 else if( input_mode .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( input_mode .eq. -1 ) then 372 call streamio_setup_read(fhandle, geomRestartStream) 373 else if( input_mode .eq. 0 ) then 374 call posixio_setup(fhandle, c_char_'r') 375 else if( input_mode .ge. 1 ) then 376 call syncio_setup_read(nsynciofiles, fhandle) 377 end if 378 call phio_constructName(fhandle, 379 & c_char_'restart' // char(0), fnamer) 380 call phstr_appendInt(fnamer, irstart) 381 call phstr_appendStr(fnamer, c_char_'.'//c_null_char) 382 call phio_openfile(fnamer, fhandle); 383 384 ithree=3 385 386 itmp = int(log10(float(myrank+1)))+1 387 388 intfromfile=0 389 call phio_readheader(fhandle, 390 & c_char_'solution' // char(0), 391 & c_loc(intfromfile), ithree, dataInt, iotype) 392c 393c.... read the values of primitive variables into q 394c 395 allocate( qold(nshg,ndof) ) 396 if(intfromfile(1).ne.0) then 397 nshg2=intfromfile(1) 398 ndof2=intfromfile(2) 399 lstep=intfromfile(3) 400 if(ndof2.ne.ndof) then 401 402 endif 403 if (nshg2 .ne. nshg) 404 & call error ('restar ', 'nshg ', nshg) 405 allocate( qread(nshg,ndof2) ) 406 iqsiz=nshg*ndof2 407 call phio_readdatablock(fhandle, 408 & c_char_'solution' // char(0), 409 & c_loc(qread),iqsiz, dataDbl,iotype) 410 qold(:,1:ndof)=qread(:,1:ndof) 411 deallocate(qread) 412 else 413 if (myrank.eq.master) then 414 if (matflg(1,1).eq.0) then ! compressible 415 warning='Solution is set to zero (with p and T to one)' 416 else 417 warning='Solution is set to zero' 418 endif 419 write(*,*) warning 420 endif 421 qold=zero 422 if (matflg(1,1).eq.0) then ! compressible 423 qold(:,1)=one ! avoid zero pressure 424 qold(:,nflow)=one ! avoid zero temperature 425 endif 426 endif 427 428 intfromfile=0 429 call phio_readheader(fhandle, 430 & c_char_'time derivative of solution' // char(0), 431 & c_loc(intfromfile), ithree, dataInt, iotype) 432 allocate( acold(nshg,ndof) ) 433 if(intfromfile(1).ne.0) then 434 nshg2=intfromfile(1) 435 ndof2=intfromfile(2) 436 lstep=intfromfile(3) 437 438 if (nshg2 .ne. nshg) 439 & call error ('restar ', 'nshg ', nshg) 440 allocate( acread(nshg,ndof2) ) 441 acread=zero 442 iacsiz=nshg*ndof2 443 call phio_readdatablock(fhandle, 444 & c_char_'time derivative of solution' // char(0), 445 & c_loc(acread), iacsiz, dataDbl,iotype) 446 acold(:,1:ndof)=acread(:,1:ndof) 447 deallocate(acread) 448 else 449 if (myrank.eq.master) then 450 warning='Time derivative of solution is set to zero (SAFE)' 451 write(*,*) warning 452 endif 453 acold=zero 454 endif 455cc 456cc.... read the header and check it against the run data 457cc 458 if (ideformwall.eq.1) then 459 460 intfromfile=0 461 call phio_readheader(fhandle, 462 & c_char_'displacement' // char(0), 463 & c_loc(intfromfile), ithree, dataInt, iotype) 464 465 nshg2=intfromfile(1) 466 ndisp=intfromfile(2) 467 lstep=intfromfile(3) 468 if(ndisp.ne.nsd) then 469 warning='WARNING ndisp not equal nsd' 470 write(*,*) warning , ndisp 471 endif 472 if (nshg2 .ne. nshg) 473 & call error ('restar ', 'nshg ', nshg) 474c 475c.... read the values of primitive variables into uold 476c 477 allocate( uold(nshg,nsd) ) 478 allocate( uread(nshg,nsd) ) 479 480 iusiz=nshg*nsd 481 482 call phio_readdatablock(fhandle, 483 & c_char_'displacement' // char(0), 484 & c_loc(uread), iusiz, dataDbl, iotype) 485 486 uold(:,1:nsd)=uread(:,1:nsd) 487 else 488 allocate( uold(nshg,nsd) ) 489 uold(:,1:nsd) = zero 490 endif 491c 492c.... close c-binary files 493c 494 call phio_closefile(fhandle) 495 496 deallocate(xread) 497 if ( numpbc > 0 ) then 498 deallocate(bcinpread) 499 deallocate(ibctmpread) 500 endif 501 deallocate(iperread) 502 if(numpe.gt.1) 503 & deallocate(ilworkread) 504 deallocate(nbcread) 505 506 return 507 994 call error ('input ','opening ', igeomBAK) 508 995 call error ('input ','opening ', igeomBAK) 509 997 call error ('input ','end file', igeomBAK) 510 998 call error ('input ','end file', igeomBAK) 511 end 512c 513c No longer called but kept around in case.... 514c 515 subroutine genpzero(iBC) 516 517 use pointer_data 518 include "common.h" 519 integer iBC(nshg) 520c 521c.... check to see if any of the nodes have a dirichlet pressure 522c 523 pzero=1 524 if (any(btest(iBC,2))) pzero=0 525 do iblk = 1, nelblb 526 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 527 do i=1, npro 528 iBCB1=miBCB(iblk)%p(i,1) 529c 530c.... check to see if any of the nodes have a Neumann pressure 531c but not periodic (note that 532c 533 if(btest(iBCB1,1)) pzero=0 534 enddo 535c 536c.... share results with other processors 537c 538 pzl=pzero 539 if (numpe .gt. 1) 540 & call MPI_ALLREDUCE (pzl, pzero, 1, 541 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 542 enddo 543 return 544 end 545