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