1 subroutine genbkb (ibksz) 2c 3c---------------------------------------------------------------------- 4c 5c This routine reads the boundary elements, reorders them and 6c generates traces for the gather/scatter operations. 7c 8c Zdenek Johan, Fall 1991. 9c---------------------------------------------------------------------- 10c 11 use dtnmod 12 use pointer_data 13c 14 include "common.h" 15!MR CHANGE 16 include "mpif.h" !Required to determine the max for itpblk 17!MR CHANGE END 18c 19 20 integer, allocatable :: ientp(:,:),iBCBtp(:,:) 21 real*8, allocatable :: BCBtp(:,:) 22 integer materb(ibksz) 23 integer intfromfile(50) ! integers read from headers 24 character*255 fname1 25 26cccccccccccccc New Phasta IO starts here cccccccccccccccccccccccccccccc 27 28 integer :: descriptor, descriptorG, GPID, color, nfiles, nfields 29 integer :: numparts, nppf, nppp, nprocs, writeLock 30 integer :: ierr_io, numprocs, itmp, itmp2 31! integer :: num_local_loop, num_global_loop 32!MR CHANGE 33 integer :: itpblktot,ierr 34!MR CHANGE END 35 36 character*255 fnamer, fname2, temp2 37 character*64 temp1, temp3 38 39 nfiles = nsynciofiles 40 nfields = nsynciofieldsreadgeombc 41 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 42 43 nppp = numparts/numpe 44 nppf = numparts/nfiles 45 46 color = int(myrank/(numparts/nfiles)) 47 itmp2 = int(log10(float(color+1)))+1 48 write (temp2,"('(''geombc-dat.'',i',i1,')')") itmp2 49 temp2=trim(temp2) 50 write (fnamer,temp2) (color+1) 51 fnamer=trim(fnamer) 52 53 ione=1 54 itwo=2 55 ieight=8 56 ieleven=11 57 itmp = int(log10(float(myrank+1)))+1 58 59! num_global_loop = 4 60! num_local_loop = nelblb 61 62cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 63 64 iel=1 65 itpblk=nelblb 66 67 68!MR CHANGE 69 ! Get the total number of different interior topologies in the whole domain. 70 ! Try to read from a field. If the field does not exist, scan the geombc file. 71 itpblktot=-1 72 call phio_readheader(igeom,'total number of boundary tpblocks' // char(0), 73 & itpblktot,ione,'integer' // char(0),iotype) 74 75! write (*,*) 'Rank: ',myrank,' boundary itpblktot intermediate:', 76! & itpblktot 77 78 if (itpblktot == -1) then 79 ! The field 'total number of different boundary tpblocks' was not found in the geombc file. 80 ! Scan all the geombc file for the 'connectivity interior' fields to get this information. 81 iblk=0 82 neltp=0 83 do while(neltp .ne. -1) 84 85 ! intfromfile is reinitialized to -1 every time. 86 ! If connectivity boundary@xxx is not found, then 87 ! readheader will return intfromfile unchanged 88 89 intfromfile(:)=-1 90 iblk = iblk+1 91 call phio_readheader(igeom,'connectivity boundary1' // char(0), 92 & intfromfile,ieight,'integer' // char(0),iotype) 93 neltp = intfromfile(1) ! -1 if fname2 was not found, >=0 otherwise 94 end do 95 itpblktot = iblk-1 96 end if 97 98 if (myrank == 0) then 99 write(*,*) 'Number of boundary topologies: ',itpblktot 100 endif 101! write (*,*) 'Rank: ',myrank,' boundary itpblktot final:', 102! & itpblktot 103 104!MR CHANGE END 105 nelblb=0 106 mattyp=0 107 ndofl = ndof 108!MR CHANGE 109! call initphmpiio( nfields, nppf, nfiles, igeom ) 110! call openfile( fnamer, 'read', igeom ) 111 112 do iblk = 1, itpblktot 113 writeLock=0; 114!MR CHANGE END 115 116 117! print *, "Loop ",iblk, myrank, itpblk, trim(fnamer) 118 119ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 120! write(*,*) 'rank, fname2',myrank, trim(adjustl(fname2)) 121 122ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 123 124 ! Synchronization for performance monitoring, as some parts do not include some topologies 125 call MPI_Barrier(MPI_COMM_WORLD,ierr) 126 call phio_readheader(igeom,'connectivity boundary1' // char(0), 127 & intfromfile,ieight,'integer' // char(0),iotype) 128 neltp =intfromfile(1) 129 nenl =intfromfile(2) 130 ipordl=intfromfile(3) 131 nshl =intfromfile(4) 132 nshlb =intfromfile(5) 133 nenbl =intfromfile(6) 134 lcsyst=intfromfile(7) 135 numnbc=intfromfile(8) 136 137!MR CHANGE 138! write (temp1,"('connectivityBoundaryHeader_',i1,'_',i1)") 139! & iblk,myrank 140! temp1=trim(temp1) 141! open(unit=14,file=temp1,status='unknown') 142! write(14,*) intfromfile(:) 143! close(14) 144!MR CHANGE END 145 146 allocate (ientp(neltp,nshl)) 147 allocate (iBCBtp(neltp,ndiBCB)) 148 allocate (BCBtp(neltp,ndBCB)) 149 iientpsiz=neltp*nshl 150 151 if (neltp==0) then 152 writeLock=1; 153 endif 154 155! print *, "neltp is ", neltp 156 157 call readdatablock(igeom,'connectivity boundary1' // char(0), 158 & ientp,iientpsiz,'integer' // char(0),iotype) 159 160!MR CHANGE 161! write (temp1,"('connectivityBoundaryDatablock_',i1,'_',i1)") 162! & iblk,myrank 163! temp1=trim(temp1) 164 165! open(unit=14,file=temp1,status='unknown') 166! do i=1,neltp 167! do j=1,nshl 168! write(14,*) ientp(i,j) 169! enddo 170! enddo 171! write(14,*) iientpsiz 172! close(14) 173!MR CHANGE END 174 175 176c 177c.... Read the boundary flux codes 178c 179 180! print *,"connectivity [] is ", trim(fname2),ientp(0,0) 181 182 183ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 184 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 185 186ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 187 188 call phio_readheader(igeom,'nbc codes1' // char(0) ,intfromfile, 189 & ieight,'integer' // char(0),iotype) 190 iiBCBtpsiz=neltp*ndiBCB 191 call readdatablock(igeom,'nbc codes1' // char(0) ,iBCBtp, 192 & iiBCBtpsiz,'integer' // char(0),iotype) 193 194!MR CHANGE 195! print *, "ndiBCB is ",ndiBCB 196! write (temp1,"('nbcCodesDatablock_',i1,'_',i1)") 197! & iblk,myrank 198! temp1=trim(temp1) 199! 200! open(unit=13,file=temp1,status='unknown') 201! do i=1,neltp 202! do j=1,ndiBCB 203! write(13,*) iBCBtp(i,j) 204! enddo 205! enddo 206! write(13,*) iiBCBtpsiz 207! close(13) 208!! iBCBtp(:,:) = 0 ! JUST FOR TEST 209!MR CHANGE END 210 211c 212c.... read the boundary condition data 213c 214 215ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 216 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 217ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 218 219 call phio_readheader(igeom,'nbc values1' // char(0) ,intfromfile, 220 & ieight,'integer' // char(0) ,iotype) 221 BCBtp = zero 222 iBCBtpsiz=neltp*ndBCB 223 call readdatablock(igeom,'nbc values1' // char(0), 224 & BCBtp,iBCBtpsiz,'double' // char(0) ,iotype) 225 226!MR CHANGE 227! write (temp1,"('nbcValuesDatablock_',i1,'_',i1)") 228! & iblk,myrank 229! temp1=trim(temp1) 230! open(unit=12,file=temp1,status='unknown') 231! do i=1,neltp 232! do j=1,ndBCB 233! write(12,*) BCBtp(i,j) 234! enddo 235! enddo 236! write(12,*) iBCBtpsiz 237! close(12) 238!MR CHANGE END 239 240c 241c This is a temporary fix until NSpre properly zeros this array where it 242c is not set. DEC has indigestion with these arrays though the 243c result is never used (never effects solution). 244c 245 246 247!MR CHANGE 248 if(writeLock==0) then 249!MR CHANGE 250! print *,"In ASSIGN ASSIGN",myrank 251 where(.not.btest(iBCBtp(:,1),0)) BCBtp(:,1)=zero 252 where(.not.btest(iBCBtp(:,1),1)) BCBtp(:,2)=zero 253 where(.not.btest(iBCBtp(:,1),3)) BCBtp(:,6)=zero 254 if(ndBCB.gt.6) then 255 do i=6,ndof 256 where(.not.btest(iBCBtp(:,1),i-1)) BCBtp(:,i+1)=zero 257 enddo 258 endif 259 where(.not.btest(iBCBtp(:,1),2)) 260 BCBtp(:,3)=zero 261 BCBtp(:,4)=zero 262 BCBtp(:,5)=zero 263 endwhere 264 265 do n=1,neltp,ibksz 266 nelblb=nelblb+1 267 npro= min(IBKSZ, neltp - n + 1) 268c 269 lcblkb(1,nelblb) = iel 270c lcblkb(2,nelblb) = iopen ! available for later use 271 lcblkb(3,nelblb) = lcsyst 272 lcblkb(4,nelblb) = ipordl 273 lcblkb(5,nelblb) = nenl 274 lcblkb(6,nelblb) = nenbl 275 lcblkb(7,nelblb) = mattyp 276 lcblkb(8,nelblb) = ndofl 277 lcblkb(9,nelblb) = nshl 278 lcblkb(10,nelblb) = nshlb ! # of shape functions per elt 279c 280c.... save the element block 281c 282 n1=n 283 n2=n+npro-1 284 materb=1 ! all one material for now 285c 286c.... allocate memory for stack arrays 287c 288 289 allocate (mienb(nelblb)%p(npro,nshl)) 290c 291 allocate (miBCB(nelblb)%p(npro,ndiBCB)) 292c 293 allocate (mBCB(nelblb)%p(npro,nshlb,ndBCB)) 294c 295 allocate (mmatb(nelblb)%p(npro)) 296c 297c.... save the boundary element block 298c 299 call gensvb (ientp(n1:n2,1:nshl), 300 & iBCBtp(n1:n2,:), BCBtp(n1:n2,:), 301 & materb, mienb(nelblb)%p, 302 & miBCB(nelblb)%p, mBCB(nelblb)%p, 303 & mmatb(nelblb)%p) 304c 305 iel=iel+npro 306 enddo 307 308!MR CHANGE 309 endif 310!MR CHANGE 311 deallocate(ientp) 312 deallocate(iBCBtp) 313 deallocate(BCBtp) 314 315 enddo 316 lcblkb(1,nelblb+1) = iel 317 318! call closefile( igeom, "read" ) 319! call finalizephmpiio( igeom ) 320 321c 322c.... return 323c 324 return 325c 326c.... end of file error handling 327c 328 911 call error ('genbcb ','end file',igeomBAK) 329c 3301000 format(a80,//, 331 & ' B o u n d a r y E l e m e n t C o n n e c t i v i t y',//, 332 & ' Elem BC codes',/, 333 & ' Number C P V H ',5x,27('Node',i1,:,2x)) 3341100 format(2x,i5,2x,4i2,3x,27i7) 335c$$$2000 format(a80,//, 336c$$$ & ' B o u n d a r y E l e m e n t B C D a t a ',//, 337c$$$ & ' Node ',3x,'mass',/, 338c$$$ & ' Number ',3x,'flux',6x,'Pressure',6x,'Heat',6x, 339c$$$ & 3('Viscous',i1,:,4x)) 3402100 format(2x,i5,1p,1x,6e12.4) 341c 342 end 343 344