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