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