1 subroutine gendat (y, ac, x, iBC, BC, 2 & iper, ilwork, 3 & shp, shgl, shpb, shglb, 4 & ifath, velbar, nsons ) 5c 6c---------------------------------------------------------------------- 7c 8c This routine inputs the geometry and the boundary conditions. 9c 10c 11c Zdenek Johan, Winter 1991. (Fortran 90) 12c---------------------------------------------------------------------- 13c 14 15 use readarrays ! used to acess nBC 16 use dtnmod 17 use pointer_data 18 use wallData ! give access to wnorm, findWallNOrm 19 include "common.h" 20 include "mpif.h" 21 22c 23c arrays in the following line are now dimensioned in readnblk 24c dimension nBC(nshg) 25c 26 dimension y(nshg,ndof), ac(nshg,ndof), 27 & x(numnp,nsd), iBC(nshg), 28 & BC(nshg,ndofBC), 29 & nodflx(numflx), ilwork(nlwork), 30 & iper(nshg) 31c 32c.... shape function declarations 33c 34 dimension shp(MAXTOP,maxsh,MAXQPT), 35 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 36 & shpb(MAXTOP,maxsh,MAXQPT), 37 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 38c 39c stuff for dynamic model s.w.avg and wall model 40c 41 dimension ifath(numnp), velbar(nfath,nflow), nsons(nfath) 42 43! Hack to get suction right on part boundaries 44! dimension BC3(numnp,5) 45 46c... Duct 47 integer iTurbWall(nshg) 48 integer nTopNodes 49 integer TopSurface(nshg) 50 51c... 52 53c.... start the timer 54c 55 56CAD call timer ('PrProces') 57 58c 59c put a barrier here so that all of the files are done reading 60c This SHOULD shield any mpi_profile information from io bottlenecks 61c 62 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 63 64c... updated, truly useful things for Duct......... 65 if(isetInitial_Duct.gt.0)then 66 call setInitial_Duct(x) ! in INIprofile.f 67 endif 68 69c 70c.... ----------------------------> Nodes <-------------------------- 71c 72c.... compute length scales 73c 74 call xyzbound(x) 75c 76c.... echo the coordinates 77c 78 if ((necho .lt. 2).and.(myrank.eq.master)) then 79 do n = 1, numnp 80 if (mod(n,50) .eq. 1) write (iecho,1000) ititle,(i,i=1,nsd) 81 write (iecho,1100) n, (x(n,i),i=1,nsd) 82 enddo 83 endif 84c 85c.... prepare periodic boundary conditions 86c 87 do i = 1,nshg 88 if (iper(i) .ne. 0) then 89 nshg0 = nshg0 - 1 90 else 91 iper(i) = i 92 endif 93 enddo 94c 95c.... ----------------------> Interior Elements <-------------------- 96c 97 ibound = 0 98c 99c.... generate the interior nodal mapping 100c 101 call genshp ( shp, shgl, nshape, nelblk) 102c 103c.... ---------------------> Boundary Conditions <------------------- 104c 105c.... read and generate the boundary condition codes (iBC array) 106c 107 call geniBC (iBC) 108c 109c.... read and generate the essential boundary conditions (BC array) 110c 111 call genBC (iBC, BC, point2x, 112 & point2ilwork, point2iper) 113 deallocate(nBC) 114 115c========================================================================================= 116c Yi Chen 117c Duct geometry8 118c... updated, truly useful things......... 119 120 121c===== specify outlet pressure 122 if(isetOutletID .gt. 0)then 123 call SetUniOutPres(BC) ! in BCprofile2.f 124 endif 125 126c==== specify inlet boundary conditions 127c if(isetInlet_Duct.gt.0)then 128c call setInlet_Duct(x,BC,iTurbWall) ! in BCprofile2.f 129c endif 130 131c==== specify blowing conditions 132 if(isetBlowing_Duct.gt.0)then 133 if (ifixBlowingVel_Duct.eq.0)then 134 call setBlowing_Duct(BC,iTurbWall) 135 else 136 call setBlowing_Duct3(x,BC,iTurbWall) ! in setBlowing_Duct3.f, fixed jet inlet velocity 137 endif 138 endif 139 140c====== specify wall conditions 141 call findTurbWall(iTurbWall) 142 143c==== apply suction patch on sides 144c suction is applied at the end so it will overwrite the velocity at any nodes shared by the no-slip walls 145 call findWallNorm(x,iBC,ilwork,iper) 146 if(isetSuctionID_Duct.gt.0)then 147 call setSuction_Duct3(x, BC, y, ilwork) 148 endif 149 150c==== 151 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 152c.... Yi Chen 153c========================================================================================== 154c 155c.... ----------------------> Boundary Elements <-------------------- 156c 157 ibound = 1 158 call gtnods 159c 160c We now take care of Direchlet to Neumann BC's. It had to move here 161c so that the IBC array was of size nshg and ready to be marked. 162c 163 164 if(nsclr.gt.0) then 165 call initDtN ! Dirichlet to Neumann module: 166 ! initialize this once only 167 do iblk = 1, nelblb ! number of blocks 168 iel = lcblkb(1,iblk) 169 npro = lcblkb(1,iblk+1) - iel 170c 171c for the DtN BC we need to mark all of the nodes that are involved. 172c 173 do i=1,npro 174c 175c if this element has the BCB AND it has not been found yet then mark it 176c 177 if(miBCB(iblk)%p(i,2).lt.0) then 178 idtn = 1 !set the flag for dtn bc's 179 do j=1,nshapeb 180 do isclr=1,nsclr 181 ignd=mienb(iblk)%p(i,j) 182 ifeature(ignd) = abs(miBCB(iblk)%p(i,2)) 183 iBC(ignd)=ior(iBC(ignd),2**13) 184 ! must mark this as a Neumann BC now 185 miBCB(iblk)%p(i,1)= 186 & ior(miBCB(iblk)%p(i,1),2**(4+isclr)) 187 end do 188 end do 189 endif 190 end do 191 end do 192 endif 193c 194c.... generate the boundary element shape functions 195c 196 call genshpb ( shpb, shglb, nshapeb, nelblb) 197c.... Evaluate the shape funcs. and their gradients at the desired quadrature 198c.... for filtering. Save these evaluations using a module 199c 200c KEJ moved them to this point because cdelsq now passed with module 201c and it is read in with velb.<stepnum>.<proc#> now 202c 203 if (iLES .gt. 0) then 204 205 call setfilt ! For setting quad. rule to use for integrating 206 call filtprep ! the hat filter. 207 if(iLES/10 .eq. 2) then 208 call setave ! For averaging cdelsq computed at quad pts 209 call aveprep(shp,x) 210 endif 211 endif 212c 213c User sets request pzero in solver.inp now 214c 215c call genpzero(iBC,iper) 216c 217 if((myrank.eq.master).and.(irscale.ge.0)) then 218 call setSPEBC(numnp,nsd) 219 call eqn_plane(point2x, iBC) 220 endif 221c 222c.... ---------------------> Initial Conditions <-------------------- 223c 224c.... generate the initial conditions and initialize time varying BC 225c 226 call genini (iBC, BC, y, 227 & ac, iper, 228 & ilwork, ifath, velbar, 229 & nsons, x, 230 & shp, shgl, shpb, shglb) 231c 232c.... close the geometry, boundary condition and material files 233c 234!MR CHANGE 235! close (igeom) 236!MR CHANGE END 237 close (ibndc) 238 if (mexist) close (imat) 239c 240c.... return 241c 242CAD call timer ('Back ') 243 return 244c 245c.... end of file error handling 246c 247999 call error ('gendat ','end file',igeom) 248c 2491000 format(a80,//, 250 & ' N o d a l C o o r d i n a t e s ',//, 251 & ' Node ',12x,3('x',i1,:,17x)) 2521100 format(1p,2x,i5,13x,3(1e12.5,7x)) 2532000 format(a80,//, 254 & ' B o u n d a r y F l u x N o d e s '//, 255 & ' index Node ') 2562100 format(1x,i5,5x,i10) 257c 258 end 259 260 261 subroutine xyzbound(x) 262 263 include "common.h" 264 include "mpif.h" 265 include "auxmpi.h" 266 267 dimension x(numnp,3) 268 269 real*8 Forout(3), Forin(3) 270 271 xlngth=maxval(x(:,1)) 272 ylngth=maxval(x(:,2)) 273 zlngth=maxval(x(:,3)) 274 if(numpe. gt. 1) then 275 Forin=(/xlngth,ylngth,zlngth/) 276 call MPI_ALLREDUCE (Forin, Forout, 3, 277 & MPI_DOUBLE_PRECISION,MPI_MAX, MPI_COMM_WORLD,ierr) 278 xmax = Forout(1) 279 ymax = Forout(2) 280 zmax = Forout(3) 281 else 282 xmax = xlngth 283 ymax = ylngth 284 zmax = zlngth 285 endif 286 xlngth=minval(x(:,1)) 287 ylngth=minval(x(:,2)) 288 zlngth=minval(x(:,3)) 289 if(numpe .gt. 1) then 290 Forin=(/xlngth,ylngth,zlngth/) 291 call MPI_ALLREDUCE (Forin, Forout, 3, 292 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 293 else 294 Forout(1) = xlngth 295 Forout(2) = ylngth 296 Forout(3) = zlngth 297 endif 298 299 xlngth = xmax-Forout(1) 300 ylngth = ymax-Forout(2) 301 zlngth = zmax-Forout(3) 302 303 if(myrank.eq.master) then 304 print 108, xlngth,ylngth,zlngth 305 endif 306 108 format(' Domain size (x,y,z):',2x,3f15.10) 307 return 308 end 309