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 43 dataInt = c_char_'integer'//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(handle, 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(handle, 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(handle, 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(igeom,'connectivity boundary1' // char(0), 167 & ientp,iientpsiz,'integer' // char(0),iotype) 168 169!MR CHANGE 170! write (temp1,"('connectivityBoundaryDatablock_',i1,'_',i1)") 171! & iblk,myrank 172! temp1=trim(temp1) 173 174! open(unit=14,file=temp1,status='unknown') 175! do i=1,neltp 176! do j=1,nshl 177! write(14,*) ientp(i,j) 178! enddo 179! enddo 180! write(14,*) iientpsiz 181! close(14) 182!MR CHANGE END 183 184 185c 186c.... Read the boundary flux codes 187c 188 189! print *,"connectivity [] is ", trim(fname2),ientp(0,0) 190 191 192ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 193 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 194 195ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 196 197 call phio_readheader(handle, 198 & c_char_'nbc codes1' // char(0), 199 & c_loc(intfromfile), ieight, dataInt, iotype) 200 iiBCBtpsiz=neltp*ndiBCB 201 call phio_readdatablock(igeom,'nbc codes1' // char(0) ,iBCBtp, 202 & iiBCBtpsiz,'integer' // char(0),iotype) 203 204!MR CHANGE 205! print *, "ndiBCB is ",ndiBCB 206! write (temp1,"('nbcCodesDatablock_',i1,'_',i1)") 207! & iblk,myrank 208! temp1=trim(temp1) 209! 210! open(unit=13,file=temp1,status='unknown') 211! do i=1,neltp 212! do j=1,ndiBCB 213! write(13,*) iBCBtp(i,j) 214! enddo 215! enddo 216! write(13,*) iiBCBtpsiz 217! close(13) 218!! iBCBtp(:,:) = 0 ! JUST FOR TEST 219!MR CHANGE END 220 221c 222c.... read the boundary condition data 223c 224 225ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 226 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 227ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 228 229 call phio_readheader(handle, 230 & c_char_'nbc values1' // char(0), 231 & c_loc(intfromfile), ieight, dataInt, iotype) 232 BCBtp = zero 233 iBCBtpsiz=neltp*ndBCB 234 call phio_readdatablock(igeom,'nbc values1' // char(0), 235 & BCBtp,iBCBtpsiz,'double' // char(0) ,iotype) 236 237!MR CHANGE 238! write (temp1,"('nbcValuesDatablock_',i1,'_',i1)") 239! & iblk,myrank 240! temp1=trim(temp1) 241! open(unit=12,file=temp1,status='unknown') 242! do i=1,neltp 243! do j=1,ndBCB 244! write(12,*) BCBtp(i,j) 245! enddo 246! enddo 247! write(12,*) iBCBtpsiz 248! close(12) 249!MR CHANGE END 250 251c 252c This is a temporary fix until NSpre properly zeros this array where it 253c is not set. DEC has indigestion with these arrays though the 254c result is never used (never effects solution). 255c 256 257 258!MR CHANGE 259 if(writeLock==0) then 260!MR CHANGE 261! print *,"In ASSIGN ASSIGN",myrank 262 where(.not.btest(iBCBtp(:,1),0)) BCBtp(:,1)=zero 263 where(.not.btest(iBCBtp(:,1),1)) BCBtp(:,2)=zero 264 where(.not.btest(iBCBtp(:,1),3)) BCBtp(:,6)=zero 265 if(ndBCB.gt.6) then 266 do i=6,ndof 267 where(.not.btest(iBCBtp(:,1),i-1)) BCBtp(:,i+1)=zero 268 enddo 269 endif 270 where(.not.btest(iBCBtp(:,1),2)) 271 BCBtp(:,3)=zero 272 BCBtp(:,4)=zero 273 BCBtp(:,5)=zero 274 endwhere 275 276 do n=1,neltp,ibksz 277 nelblb=nelblb+1 278 npro= min(IBKSZ, neltp - n + 1) 279c 280 lcblkb(1,nelblb) = iel 281c lcblkb(2,nelblb) = iopen ! available for later use 282 lcblkb(3,nelblb) = lcsyst 283 lcblkb(4,nelblb) = ipordl 284 lcblkb(5,nelblb) = nenl 285 lcblkb(6,nelblb) = nenbl 286 lcblkb(7,nelblb) = mattyp 287 lcblkb(8,nelblb) = ndofl 288 lcblkb(9,nelblb) = nshl 289 lcblkb(10,nelblb) = nshlb ! # of shape functions per elt 290c 291c.... save the element block 292c 293 n1=n 294 n2=n+npro-1 295 materb=1 ! all one material for now 296c 297c.... allocate memory for stack arrays 298c 299 300 allocate (mienb(nelblb)%p(npro,nshl)) 301c 302 allocate (miBCB(nelblb)%p(npro,ndiBCB)) 303c 304 allocate (mBCB(nelblb)%p(npro,nshlb,ndBCB)) 305c 306 allocate (mmatb(nelblb)%p(npro)) 307c 308c.... save the boundary element block 309c 310 call gensvb (ientp(n1:n2,1:nshl), 311 & iBCBtp(n1:n2,:), BCBtp(n1:n2,:), 312 & materb, mienb(nelblb)%p, 313 & miBCB(nelblb)%p, mBCB(nelblb)%p, 314 & mmatb(nelblb)%p) 315c 316 iel=iel+npro 317 enddo 318 319!MR CHANGE 320 endif 321!MR CHANGE 322 deallocate(ientp) 323 deallocate(iBCBtp) 324 deallocate(BCBtp) 325 326 enddo 327 lcblkb(1,nelblb+1) = iel 328 329! call closefile( igeom, "read" ) 330! call finalizephmpiio( igeom ) 331 332c 333c.... return 334c 335 return 336c 337c.... end of file error handling 338c 339 911 call error ('genbcb ','end file',igeomBAK) 340c 3411000 format(a80,//, 342 & ' 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',//, 343 & ' Elem BC codes',/, 344 & ' Number C P V H ',5x,27('Node',i1,:,2x)) 3451100 format(2x,i5,2x,4i2,3x,27i7) 346c$$$2000 format(a80,//, 347c$$$ & ' B o u n d a r y E l e m e n t B C D a t a ',//, 348c$$$ & ' Node ',3x,'mass',/, 349c$$$ & ' Number ',3x,'flux',6x,'Pressure',6x,'Heat',6x, 350c$$$ & 3('Viscous',i1,:,4x)) 3512100 format(2x,i5,1p,1x,6e12.4) 352c 353 end 354 355