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, 73 & 'total number of boundary tpblocks' // char(0), itpblktot,ione, 74 & 'integer' // char(0),iotype) 75 76! write (*,*) 'Rank: ',myrank,' boundary itpblktot intermediate:', 77! & itpblktot 78 79 if (itpblktot == -1) then 80 ! The field 'total number of different boundary tpblocks' was not found in the geombc file. 81 ! Scan all the geombc file for the 'connectivity interior' fields to get this information. 82 iblk=0 83 neltp=0 84 do while(neltp .ne. -1) 85 86 ! intfromfile is reinitialized to -1 every time. 87 ! If connectivity boundary@xxx is not found, then 88 ! readheader will return intfromfile unchanged 89 90 intfromfile(:)=-1 91 iblk = iblk+1 92 call phio_readheader(igeom,'connectivity boundary1' // char(0), 93 & intfromfile,ieight,'integer' // char(0),iotype) 94 neltp = intfromfile(1) ! -1 if fname2 was not found, >=0 otherwise 95 end do 96 itpblktot = iblk-1 97 end if 98 99 if (myrank == 0) then 100 write(*,*) 'Number of boundary topologies: ',itpblktot 101 endif 102! write (*,*) 'Rank: ',myrank,' boundary itpblktot final:', 103! & itpblktot 104 105!MR CHANGE END 106 nelblb=0 107 mattyp=0 108 ndofl = ndof 109!MR CHANGE 110! call initphmpiio( nfields, nppf, nfiles, igeom ) 111! call openfile( fnamer, 'read', igeom ) 112 113 do iblk = 1, itpblktot 114 writeLock=0; 115!MR CHANGE END 116 117 118! print *, "Loop ",iblk, myrank, itpblk, trim(fnamer) 119 120ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 121! write(*,*) 'rank, fname2',myrank, trim(adjustl(fname2)) 122 123ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 124 125 ! Synchronization for performance monitoring, as some parts do not include some topologies 126 call MPI_Barrier(MPI_COMM_WORLD,ierr) 127 call phio_readheader(igeom,'connectivity boundary1' // char(0), 128 & intfromfile,ieight,'integer' // char(0),iotype) 129 neltp =intfromfile(1) 130 nenl =intfromfile(2) 131 ipordl=intfromfile(3) 132 nshl =intfromfile(4) 133 nshlb =intfromfile(5) 134 nenbl =intfromfile(6) 135 lcsyst=intfromfile(7) 136 numnbc=intfromfile(8) 137 138!MR CHANGE 139! write (temp1,"('connectivityBoundaryHeader_',i1,'_',i1)") 140! & iblk,myrank 141! temp1=trim(temp1) 142! open(unit=14,file=temp1,status='unknown') 143! write(14,*) intfromfile(:) 144! close(14) 145!MR CHANGE END 146 147 allocate (ientp(neltp,nshl)) 148 allocate (iBCBtp(neltp,ndiBCB)) 149 allocate (BCBtp(neltp,ndBCB)) 150 iientpsiz=neltp*nshl 151 152 if (neltp==0) then 153 writeLock=1; 154 endif 155 156! print *, "neltp is ", neltp 157 158 call phio_readdatablock(igeom,'connectivity boundary1' // char(0), 159 & ientp,iientpsiz,'integer' // char(0),iotype) 160 161!MR CHANGE 162! write (temp1,"('connectivityBoundaryDatablock_',i1,'_',i1)") 163! & iblk,myrank 164! temp1=trim(temp1) 165 166! open(unit=14,file=temp1,status='unknown') 167! do i=1,neltp 168! do j=1,nshl 169! write(14,*) ientp(i,j) 170! enddo 171! enddo 172! write(14,*) iientpsiz 173! close(14) 174!MR CHANGE END 175 176 177c 178c.... Read the boundary flux codes 179c 180 181! print *,"connectivity [] is ", trim(fname2),ientp(0,0) 182 183 184ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 185 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 186 187ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 188 189 call phio_readheader(igeom,'nbc codes1' // char(0) ,intfromfile, 190 & ieight,'integer' // char(0),iotype) 191 iiBCBtpsiz=neltp*ndiBCB 192 call phio_readdatablock(igeom,'nbc codes1' // char(0) ,iBCBtp, 193 & iiBCBtpsiz,'integer' // char(0),iotype) 194 195!MR CHANGE 196! print *, "ndiBCB is ",ndiBCB 197! write (temp1,"('nbcCodesDatablock_',i1,'_',i1)") 198! & iblk,myrank 199! temp1=trim(temp1) 200! 201! open(unit=13,file=temp1,status='unknown') 202! do i=1,neltp 203! do j=1,ndiBCB 204! write(13,*) iBCBtp(i,j) 205! enddo 206! enddo 207! write(13,*) iiBCBtpsiz 208! close(13) 209!! iBCBtp(:,:) = 0 ! JUST FOR TEST 210!MR CHANGE END 211 212c 213c.... read the boundary condition data 214c 215 216ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 217 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 218ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 219 220 call phio_readheader(igeom,'nbc values1' // char(0) ,intfromfile, 221 & ieight,'integer' // char(0) ,iotype) 222 BCBtp = zero 223 iBCBtpsiz=neltp*ndBCB 224 call phio_readdatablock(igeom,'nbc values1' // char(0), 225 & BCBtp,iBCBtpsiz,'double' // char(0) ,iotype) 226 227!MR CHANGE 228! write (temp1,"('nbcValuesDatablock_',i1,'_',i1)") 229! & iblk,myrank 230! temp1=trim(temp1) 231! open(unit=12,file=temp1,status='unknown') 232! do i=1,neltp 233! do j=1,ndBCB 234! write(12,*) BCBtp(i,j) 235! enddo 236! enddo 237! write(12,*) iBCBtpsiz 238! close(12) 239!MR CHANGE END 240 241c 242c This is a temporary fix until NSpre properly zeros this array where it 243c is not set. DEC has indigestion with these arrays though the 244c result is never used (never effects solution). 245c 246 247 248!MR CHANGE 249 if(writeLock==0) then 250!MR CHANGE 251! print *,"In ASSIGN ASSIGN",myrank 252 where(.not.btest(iBCBtp(:,1),0)) BCBtp(:,1)=zero 253 where(.not.btest(iBCBtp(:,1),1)) BCBtp(:,2)=zero 254 where(.not.btest(iBCBtp(:,1),3)) BCBtp(:,6)=zero 255 if(ndBCB.gt.6) then 256 do i=6,ndof 257 where(.not.btest(iBCBtp(:,1),i-1)) BCBtp(:,i+1)=zero 258 enddo 259 endif 260 where(.not.btest(iBCBtp(:,1),2)) 261 BCBtp(:,3)=zero 262 BCBtp(:,4)=zero 263 BCBtp(:,5)=zero 264 endwhere 265 266 do n=1,neltp,ibksz 267 nelblb=nelblb+1 268 npro= min(IBKSZ, neltp - n + 1) 269c 270 lcblkb(1,nelblb) = iel 271c lcblkb(2,nelblb) = iopen ! available for later use 272 lcblkb(3,nelblb) = lcsyst 273 lcblkb(4,nelblb) = ipordl 274 lcblkb(5,nelblb) = nenl 275 lcblkb(6,nelblb) = nenbl 276 lcblkb(7,nelblb) = mattyp 277 lcblkb(8,nelblb) = ndofl 278 lcblkb(9,nelblb) = nshl 279 lcblkb(10,nelblb) = nshlb ! # of shape functions per elt 280c 281c.... save the element block 282c 283 n1=n 284 n2=n+npro-1 285 materb=1 ! all one material for now 286c 287c.... allocate memory for stack arrays 288c 289 290 allocate (mienb(nelblb)%p(npro,nshl)) 291c 292 allocate (miBCB(nelblb)%p(npro,ndiBCB)) 293c 294 allocate (mBCB(nelblb)%p(npro,nshlb,ndBCB)) 295c 296 allocate (mmatb(nelblb)%p(npro)) 297c 298c.... save the boundary element block 299c 300 call gensvb (ientp(n1:n2,1:nshl), 301 & iBCBtp(n1:n2,:), BCBtp(n1:n2,:), 302 & materb, mienb(nelblb)%p, 303 & miBCB(nelblb)%p, mBCB(nelblb)%p, 304 & mmatb(nelblb)%p) 305c 306 iel=iel+npro 307 enddo 308 309!MR CHANGE 310 endif 311!MR CHANGE 312 deallocate(ientp) 313 deallocate(iBCBtp) 314 deallocate(BCBtp) 315 316 enddo 317 lcblkb(1,nelblb+1) = iel 318 319! call closefile( igeom, "read" ) 320! call finalizephmpiio( igeom ) 321 322c 323c.... return 324c 325 return 326c 327c.... end of file error handling 328c 329 911 call error ('genbcb ','end file',igeomBAK) 330c 3311000 format(a80,//, 332 & ' 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',//, 333 & ' Elem BC codes',/, 334 & ' Number C P V H ',5x,27('Node',i1,:,2x)) 3351100 format(2x,i5,2x,4i2,3x,27i7) 336c$$$2000 format(a80,//, 337c$$$ & ' B o u n d a r y E l e m e n t B C D a t a ',//, 338c$$$ & ' Node ',3x,'mass',/, 339c$$$ & ' Number ',3x,'flux',6x,'Pressure',6x,'Heat',6x, 340c$$$ & 3('Viscous',i1,:,4x)) 3412100 format(2x,i5,1p,1x,6e12.4) 342c 343 end 344 345