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, 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 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( input_mode .eq. -1 ) then 88 call streamio_setup_read(fhandle, geomRestartStream) 89 else if( input_mode .eq. 0 ) then 90 call posixio_setup(fhandle, c_char_'r') 91 else if( input_mode .ge. 1 ) then 92 call syncio_setup_read(nsynciofiles, fhandle) 93 end if 94 call phio_constructName(fhandle, 95 & c_char_'geombc' // char(0), fname1) 96 call phio_openfile(fname1, 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(fhandle); 369c.... Read restart files 370 if( input_mode .eq. -1 ) then 371 call streamio_setup_read(fhandle, geomRestartStream) 372 else if( input_mode .eq. 0 ) then 373 call posixio_setup(fhandle, c_char_'r') 374 else if( input_mode .ge. 1 ) then 375 call syncio_setup_read(nsynciofiles, fhandle) 376 end if 377 call phio_constructName(fhandle, 378 & c_char_'restart' // char(0), fnamer) 379 call phio_appendInt(fnamer, irstart) 380 call phio_openfile(fnamer, fhandle); 381 382 ithree=3 383 384 itmp = int(log10(float(myrank+1)))+1 385 386 intfromfile=0 387 call phio_readheader(fhandle, 388 & c_char_'solution' // char(0), 389 & c_loc(intfromfile), ithree, dataInt, iotype) 390c 391c.... read the values of primitive variables into q 392c 393 allocate( qold(nshg,ndof) ) 394 if(intfromfile(1).ne.0) then 395 nshg2=intfromfile(1) 396 ndof2=intfromfile(2) 397 lstep=intfromfile(3) 398 if(ndof2.ne.ndof) then 399 400 endif 401 if (nshg2 .ne. nshg) 402 & call error ('restar ', 'nshg ', nshg) 403 allocate( qread(nshg,ndof2) ) 404 iqsiz=nshg*ndof2 405 call phio_readdatablock(fhandle, 406 & c_char_'solution' // char(0), 407 & c_loc(qread),iqsiz, dataDbl,iotype) 408 qold(:,1:ndof)=qread(:,1:ndof) 409 deallocate(qread) 410 else 411 if (myrank.eq.master) then 412 if (matflg(1,1).eq.0) then ! compressible 413 warning='Solution is set to zero (with p and T to one)' 414 else 415 warning='Solution is set to zero' 416 endif 417 write(*,*) warning 418 endif 419 qold=zero 420 if (matflg(1,1).eq.0) then ! compressible 421 qold(:,1)=one ! avoid zero pressure 422 qold(:,nflow)=one ! avoid zero temperature 423 endif 424 endif 425 426 intfromfile=0 427 call phio_readheader(fhandle, 428 & c_char_'time derivative of solution' // char(0), 429 & c_loc(intfromfile), ithree, dataInt, iotype) 430 allocate( acold(nshg,ndof) ) 431 if(intfromfile(1).ne.0) then 432 nshg2=intfromfile(1) 433 ndof2=intfromfile(2) 434 lstep=intfromfile(3) 435 436 if (nshg2 .ne. nshg) 437 & call error ('restar ', 'nshg ', nshg) 438 allocate( acread(nshg,ndof2) ) 439 acread=zero 440 iacsiz=nshg*ndof2 441 call phio_readdatablock(fhandle, 442 & c_char_'time derivative of solution' // char(0), 443 & c_loc(acread), iacsiz, dataDbl,iotype) 444 acold(:,1:ndof)=acread(:,1:ndof) 445 deallocate(acread) 446 else 447 if (myrank.eq.master) then 448 warning='Time derivative of solution is set to zero (SAFE)' 449 write(*,*) warning 450 endif 451 acold=zero 452 endif 453cc 454cc.... read the header and check it against the run data 455cc 456 if (ideformwall.eq.1) then 457 458 intfromfile=0 459 call phio_readheader(fhandle, 460 & c_char_'displacement' // char(0), 461 & c_loc(intfromfile), ithree, dataInt, iotype) 462 463 nshg2=intfromfile(1) 464 ndisp=intfromfile(2) 465 lstep=intfromfile(3) 466 if(ndisp.ne.nsd) then 467 warning='WARNING ndisp not equal nsd' 468 write(*,*) warning , ndisp 469 endif 470 if (nshg2 .ne. nshg) 471 & call error ('restar ', 'nshg ', nshg) 472c 473c.... read the values of primitive variables into uold 474c 475 allocate( uold(nshg,nsd) ) 476 allocate( uread(nshg,nsd) ) 477 478 iusiz=nshg*nsd 479 480 call phio_readdatablock(fhandle, 481 & c_char_'displacement' // char(0), 482 & c_loc(uread), iusiz, dataDbl, iotype) 483 484 uold(:,1:nsd)=uread(:,1:nsd) 485 else 486 allocate( uold(nshg,nsd) ) 487 uold(:,1:nsd) = zero 488 endif 489c 490c.... close c-binary files 491c 492 call phio_closefile(fhandle) 493 494 deallocate(xread) 495 if ( numpbc > 0 ) then 496 deallocate(bcinpread) 497 deallocate(ibctmpread) 498 endif 499 deallocate(iperread) 500 if(numpe.gt.1) 501 & deallocate(ilworkread) 502 deallocate(nbcread) 503 504 return 505 994 call error ('input ','opening ', igeomBAK) 506 995 call error ('input ','opening ', igeomBAK) 507 997 call error ('input ','end file', igeomBAK) 508 998 call error ('input ','end file', igeomBAK) 509 end 510c 511c No longer called but kept around in case.... 512c 513 subroutine genpzero(iBC) 514 515 use pointer_data 516 include "common.h" 517 integer iBC(nshg) 518c 519c.... check to see if any of the nodes have a dirichlet pressure 520c 521 pzero=1 522 if (any(btest(iBC,2))) pzero=0 523 do iblk = 1, nelblb 524 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 525 do i=1, npro 526 iBCB1=miBCB(iblk)%p(i,1) 527c 528c.... check to see if any of the nodes have a Neumann pressure 529c but not periodic (note that 530c 531 if(btest(iBCB1,1)) pzero=0 532 enddo 533c 534c.... share results with other processors 535c 536 pzl=pzero 537 if (numpe .gt. 1) 538 & call MPI_ALLREDUCE (pzl, pzero, 1, 539 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 540 enddo 541 return 542 end 543