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