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 include "common.h" 19 include "mpif.h" 20 21c 22c arrays in the following line are now dimensioned in readnblk 23c dimension nBC(nshg) 24c 25 dimension y(nshg,ndof), ac(nshg,ndof), 26 & x(numnp,nsd), iBC(nshg), 27 & BC(nshg,ndofBC), 28 & nodflx(numflx), ilwork(nlwork), 29 & iper(nshg) 30c 31c.... shape function declarations 32c 33 dimension shp(MAXTOP,maxsh,MAXQPT), 34 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 35 & shpb(MAXTOP,maxsh,MAXQPT), 36 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 37c 38c stuff for dynamic model s.w.avg and wall model 39c 40 dimension ifath(numnp), velbar(nfath,nflow), nsons(nfath) 41c 42c.... start the timer 43c 44 45CAD call timer ('PrProces') 46 47c 48c put a barrier here so that all of the files are done reading 49c This SHOULD shield any mpi_profile information from io bottlenecks 50c 51 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 52 53c 54c.... ----------------------------> Nodes <-------------------------- 55c 56c.... compute length scales 57c 58 call xyzbound(x) 59c 60c.... echo the coordinates 61c 62 if ((necho .lt. 2).and.(myrank.eq.master)) then 63 do n = 1, numnp 64 if (mod(n,50) .eq. 1) write (iecho,1000) ititle,(i,i=1,nsd) 65 write (iecho,1100) n, (x(n,i),i=1,nsd) 66 enddo 67 endif 68c 69c.... prepare periodic boundary conditions 70c 71 do i = 1,nshg 72 if (iper(i) .ne. 0) then 73 nshg0 = nshg0 - 1 74 else 75 iper(i) = i 76 endif 77 enddo 78c 79c.... ----------------------> Interior Elements <-------------------- 80c 81 ibound = 0 82c 83c.... generate the interior nodal mapping 84c 85 call genshp ( shp, shgl, nshape, nelblk) 86c 87c.... ---------------------> Boundary Conditions <------------------- 88c 89c.... read and generate the boundary condition codes (iBC array) 90c 91 call geniBC (iBC) 92c 93c.... read and generate the essential boundary conditions (BC array) 94c 95 call genBC (iBC, BC, point2x, 96 & point2ilwork, point2iper) 97 deallocate(nBC) 98c 99c.... ----------------------> Boundary Elements <-------------------- 100c 101 ibound = 1 102 call gtnods 103c 104c We now take care of Direchlet to Neumann BC's. It had to move here 105c so that the IBC array was of size nshg and ready to be marked. 106c 107 108 if(nsclr.gt.0) then 109 call initDtN ! Dirichlet to Neumann module: 110 ! initialize this once only 111 do iblk = 1, nelblb ! number of blocks 112 iel = lcblkb(1,iblk) 113 npro = lcblkb(1,iblk+1) - iel 114c 115c for the DtN BC we need to mark all of the nodes that are involved. 116c 117 do i=1,npro 118c 119c if this element has the BCB AND it has not been found yet then mark it 120c 121 if(miBCB(iblk)%p(i,2).lt.0) then 122 idtn = 1 !set the flag for dtn bc's 123 do j=1,nshapeb 124 do isclr=1,nsclr 125 ignd=mienb(iblk)%p(i,j) 126 ifeature(ignd) = abs(miBCB(iblk)%p(i,2)) 127 iBC(ignd)=ior(iBC(ignd),2**13) 128 ! must mark this as a Neumann BC now 129 miBCB(iblk)%p(i,1)= 130 & ior(miBCB(iblk)%p(i,1),2**(4+isclr)) 131 end do 132 end do 133 endif 134 end do 135 end do 136 endif 137c 138c.... generate the boundary element shape functions 139c 140 call genshpb ( shpb, shglb, nshapeb, nelblb) 141c.... Evaluate the shape funcs. and their gradients at the desired quadrature 142c.... for filtering. Save these evaluations using a module 143c 144c KEJ moved them to this point because cdelsq now passed with module 145c and it is read in with velb.<stepnum>.<proc#> now 146c 147 if (iLES .gt. 0) then 148 149 call setfilt ! For setting quad. rule to use for integrating 150 call filtprep ! the hat filter. 151 if(iLES/10 .eq. 2) then 152 call setave ! For averaging cdelsq computed at quad pts 153 call aveprep(shp,x) 154 endif 155 endif 156c 157c User sets request pzero in solver.inp now 158c 159c call genpzero(iBC,iper) 160c 161 if((myrank.eq.master).and.(irscale.ge.0)) then 162 call setSPEBC(numnp,nsd) 163 call eqn_plane(point2x, iBC) 164 endif 165c 166c.... ---------------------> Initial Conditions <-------------------- 167c 168c.... generate the initial conditions and initialize time varying BC 169c 170 call genini (iBC, BC, y, 171 & ac, iper, 172 & ilwork, ifath, velbar, 173 & nsons, x, 174 & shp, shgl, shpb, shglb) 175c 176c.... close the geometry, boundary condition and material files 177c 178!MR CHANGE 179! close (igeom) 180!MR CHANGE END 181 close (ibndc) 182 if (mexist) close (imat) 183c 184c.... return 185c 186CAD call timer ('Back ') 187 return 188c 189c.... end of file error handling 190c 191999 call error ('gendat ','end file',igeom) 192c 1931000 format(a80,//, 194 & ' N o d a l C o o r d i n a t e s ',//, 195 & ' Node ',12x,3('x',i1,:,17x)) 1961100 format(1p,2x,i5,13x,3(1e12.5,7x)) 1972000 format(a80,//, 198 & ' B o u n d a r y F l u x N o d e s '//, 199 & ' index Node ') 2002100 format(1x,i5,5x,i10) 201c 202 end 203 204 205 subroutine xyzbound(x) 206 207 include "common.h" 208 include "mpif.h" 209 include "auxmpi.h" 210 211 dimension x(numnp,3) 212 213 real*8 Forout(3), Forin(3) 214 215 xlngth=maxval(x(:,1)) 216 ylngth=maxval(x(:,2)) 217 zlngth=maxval(x(:,3)) 218 if(numpe. gt. 1) then 219 Forin=(/xlngth,ylngth,zlngth/) 220 call MPI_ALLREDUCE (Forin, Forout, 3, 221 & MPI_DOUBLE_PRECISION,MPI_MAX, MPI_COMM_WORLD,ierr) 222 xmax = Forout(1) 223 ymax = Forout(2) 224 zmax = Forout(3) 225 else 226 xmax = xlngth 227 ymax = ylngth 228 zmax = zlngth 229 endif 230 xlngth=minval(x(:,1)) 231 ylngth=minval(x(:,2)) 232 zlngth=minval(x(:,3)) 233 if(numpe .gt. 1) then 234 Forin=(/xlngth,ylngth,zlngth/) 235 call MPI_ALLREDUCE (Forin, Forout, 3, 236 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 237 else 238 Forout(1) = xlngth 239 Forout(2) = ylngth 240 Forout(3) = zlngth 241 endif 242 243 xlngth = xmax-Forout(1) 244 ylngth = ymax-Forout(2) 245 zlngth = zmax-Forout(3) 246 247 if(myrank.eq.master) then 248 print 108, xlngth,ylngth,zlngth 249 endif 250 108 format(' Domain size (x,y,z):',2x,3f15.10) 251 return 252 end 253