1 subroutine INIprofile(BC,y,x) 2 3 include "common.h" 4 real*8 BC(nshg,ndofBC) 5 real*8 y(nshg,ndof) 6 real*8 x(nshg,nsd) 7 8 if(iI2Binlet .gt. 0)then ! useless for geometry with contraction 9 call TakeBCfromIC_ScaleInlVelX(BC,y,x) 10 endif 11 12 if(isetInitial .gt.0)then ! obscure 13 call setInitial(x,y) 14 endif 15 16 return 17 end 18c.............................................................................. 19 subroutine TakeBCfromIC_ScaleInlVelX(BC,y,x) 20 21 use BCsfIDmap 22 include "common.h" 23 include "mpif.h" 24 include "auxmpi.h" 25 26 integer inlf_tmp(nshg) 27 real*8 BC(nshg,ndofBC) 28 real*8 y(nshg,ndof) 29 real*8 x(nshg,nsd) 30 real*8 u_max,u_max_all 31 32 call sfID2np(iI2Binlet,ninlet,inlf_tmp) 33 if(ninlet.gt.0) then ! only for process including inlet nodes 34 allocate(inlf(ninlet)) 35 inlf(1:ninlet)=inlf_tmp(1:ninlet) ! generate map mapping inlet nodes to nshg nodes 36 do i=1,ninlet 37 nn=inlf(i) 38 if(y(nn,1).gt.u_max)u_max=y(nn,1) 39 enddo 40 endif 41 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 42 call MPI_ALLREDUCE( u_max, u_max_all, 1, 43 & MPI_REAL8, MPI_MAX, MPI_COMM_WORLD,ierr) 44 if(myrank.eq.0)then 45 write(*,*)'' 46 write(*,556)'u_max_all=',u_max_all ! find the largest x velocity at inlet 47 endif 48c apply IC to BC, must be run before applying BC to IC later in genini 49 if(ninlet.gt.0) then 50 BC(inlf(:),2) = y(inlf(:),5) ! set temp 51 BC(inlf(:),3) = y(inlf(:),1)/u_max_all*inletVelX ! set and scale x velocity 52 BC(inlf(:),4) = 0 ! set y velocity 53 BC(inlf(:),5) = 0 ! set z velocity 54 if(nsclr.eq.1)then 55 BC(inlf(:),7) = y(inlf(:),6) ! set scalar_1 56 endif 57 endif 58555 format(I3,F10.3) 59556 format(A,F10.3) 60 61 return 62 end 63 64c.............................................................................. 65c... this subroutine gives the formular of initial conditions, user may want to 66c... change this file and recompile the code, just like other commercial code which 67c... supports User Defined Function, you need to recompile and relink the code 68 subroutine setInitial(x,y) 69 70 use turbSA ! this gives us d2wall(1:nshg) nodal distance to the closest wall 71 use BCsfIDmap 72 include "common.h" 73 include "mpif.h" 74 include "auxmpi.h" 75 76 real*8 x(nshg,nsd) 77 real*8 y(nshg,ndof) 78 if(isetinitial.eq.1) then 79 80c... user gives the formular of initial conditions from here 81 y(:,1)= xvel_ini 82 y(:,4)= pres_ini 83 y(:,2)=yvel_ini 84 y(:,3)=zvel_ini 85 y(:,5)=temp_ini 86 if(nsclr.eq.1)then 87 y(:,6)=evis_ini 88 endif 89!moving previously coded spatially varying pressure to 10 so that 1 is simple 90 elseif (isetInitial.eq.10) then 91 if(ninlet.gt.0)then 92 xinlet=x(inlf(1),1) 93 else 94 xinlet=1e18 95 endif 96 97 if(noutlet.gt.0)then 98 xoutlet=x(outf(1),1) 99 else 100 xoutlet=-1e18 101 endif 102 103 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 104 call MPI_ALLREDUCE(xinlet, xinlet_all, 1, 105 & MPI_REAL8, MPI_MIN, MPI_COMM_WORLD,ierr) 106 call MPI_ALLREDUCE(xoutlet, xoutlet_all, 1, 107 & MPI_REAL8, MPI_MAX, MPI_COMM_WORLD,ierr) 108 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 109 do nn=1,nshg 110 xdel=(x(nn,1)-xinlet_all)/(xoutlet_all-xinlet_all) 111 xdel=1-tanh(7*xdel) 112 y(nn,1)= xvel_ini + (inletVelX-xvel_ini)*xdel 113 y(nn,4)= outPres1 + (pres_ini-outPres1)*xdel 114 enddo 115 y(:,2)=yvel_ini 116 y(:,3)=zvel_ini 117 y(:,5)=temp_ini 118 if(nsclr.eq.1)then 119 y(:,6)=evis_ini 120 endif 121 elseif (isetInitial.eq.2) then 122 !---------------------------------- 123 ! set initial condition for blower 124 !---------------------------------- 125 126! blthickness=1.0e-4 127! yfloor=-0.0333375 - 1e-7 128! xthroat=-0.3048 + 1e-5 129! xSlit = 0.025 !x-coordinate of the blower slit 130! 131! p0 = 121800 !total pressure in blower 132! T0 = 305 !total temperature 133! 134! do nn=1,nshg 135! xcoor = x(nn,1) 136! 137! if((xcoor .gt. xthroat) .and. (x(nn,2).lt.yfloor)) then 138! y(nn,3:4)=0 !velocity 139! y(nn,6)=1e-3 !eddy viscosity 140! 141! if(xcoor .le. xSlit) then !test whether the point is infront of or behind the slit. If infront, use M = 0.8 142! !4th order polynomial best fit for an exit Mach number of 0.8 143! MEuler=0.0806961 + xcoor*(2.96302 + xcoor*(72.5101 + xcoor*(13914.6 + xcoor*1.09778e6))) 144! else 145! MEuler = 0.8 146! endif 147! 148! pEuler = p0/((1 + 0.2*M*M)**3.5) 149! TEuler = T0/(1 + 0.2*M*M) 150! uEuler = MEuler*sqrt(1.4*287.06*TEuler) 151! 152! bldamper=min(one, d2wall(nn)/blthickness) 153! y(nn,2)=uEuler*bldamper 154! y(nn,1)=pEuler 155! y(nn,5)=TEuler 156! endif 157! enddo 158 endif 159 160 return 161 end 162 163c========================================================================= 164c Set the initial condition based on coordiantes 165c========================================================================= 166 167 subroutine setInitial_Duct(x) 168 169 use readarrays 170 include "common.h" 171 include "mpif.h" 172 include "auxmpi.h" 173 174 real*8 x(nshg,nsd) 175 real*8 xcoor,ycoor,zcoor 176 integer nn 177 real*8 xVel, pres, Temp 178!p u v w T s(calar_1) in array qold 179!x y z in array x 180 181 182 183 a1 = 329.7 184 b1 = -0.07732 185 c1 = 9.82 186 a2 = 6.431 187 b2 = 1.621 188 c2 = 0.8998 189 a3 = -4.637 190 b3 = 1.728 191 c3 = 0.2168 192 193 194 if (iDuctgeometryType .eq. 6) then 195 xContraStart = -74*0.0254; 196 elseif (iDuctgeometryType .eq. 8) then 197 xContraStart = -84*0.0254; 198 endif 199 xContraEnd = xContraStart + 72*0.0254 200 201 do nn=1,nshg 202 203 xcoor=x(nn,1) 204 ycoor=x(nn,2) 205 zcoor=x(nn,3) 206c.... [xVel,pres,Temp] are functions of [xcoor,ycoor,zcoor] 207 208c if(ycoor<(-0.587375+1.0e-4))then 209c ry=(ycoor+0.587375)/1.0e-4 210c elseif(ycoor>(0.587375+1.0e-4))then 211c ry=(0.587375-ycoor)/1.0e-4 212c else 213c ry=1.0 214c endif 215 216c if(zcoor<(-0.587375+1.0e-4))then 217c rz=(zcoor+0.587375)/1.0e-4 218c elseif(zcoor>(0.587375+1.0e-4))then 219c rz=(0.587375-zcoor)/1.0e-4 220c else 221c rz=1.0 222c endif 223 224c ry=max(0.0,ry) 225c rz=max(0.0,rz) 226 227c if(xcoor.le.0.0.and.xcoor.gt.-12*0.0254)then 228c xVel=120 229c elseif(xcoor.lt.-12*0.0254)then 230c xVel=(1.013158+120)/2+(120-1.013158)/2* 231c & tanh((xcoor-(-84*0.0254-12*0.0254)/2)*4) 232c else 233c xVel=max(0.0,120*(1-(xcoor/(1.6*4.5*0.0254+0.85))**2)) 234c endif 235 236c Temp=317+(-tanh(4*(xcoor+0.35+0.254))* 237c & (13.578+5)*0.5+(13.578-5)*0.5)*ry*rz 238c pres=-tanh(4*(xcoor+0.35+0.254))* 239c & (116600-97000)*0.5+(116600+97000)*0.5 240c... above, Onkar's method 241 242c if(xcoor.lt.0.0 .and. 243c & xcoor.gt.12*0.0254 .and. 244c & ycoor .lt. -0.04446)then 245c xVel=0 246c pres=97800 247c Temp=317 248c else 249c if(xcoor .ge. -0.0254)then 250c xVel=158 251c pres=97800 252c Temp=318 253c elseif(xcoor .le. -1.0163)then 254c Temp=330 255c pres=111896 256c xVel=1.0132 257c else 258c pres=97800+0.5*(111896-97800)* 259c & (tanh(-3.4*((xcoor+0.52085)/0.49545))+1) 260c Temp=318+0.5*(330-318)* 261c & (tanh(-3.4*((xcoor+0.52085)/0.49545))+1) 262c xVel=158+0.5*(1.0132-158)* 263c & (tanh(-3.4*((xcoor+0.52085)/0.49545))+1) 264c endif 265c endif 266c.................... 267 268 269c pres= 97000 270c xvel= 0.1 271 272c if (xcoor .gt. xContraEnd)then 273c Temp = 320 274c elseif (xcoor . le. xContraEnd)then 275c xSim = xcoor - xContraStart 276c Temp = a1*exp(-((xSim-b1)/c1)**2) + 277c & a2*exp(-((xSim-b2)/c2)**2) + 278c & a3*exp(-((xSim-b3)/c3)**2) 279c endif 280c varSA = 0.0 281 282c qold(nn,1)=pres 283c qold(nn,2)=xvel 284c qold(nn,3)=0 285c qold(nn,4)=0 286c qold(nn,5)=Temp 287c qold(nn,6)=varSA 288 289c... above Yi Chen's Method 290 291 if(xcoor .lt. 0.0 .and. 292 & xcoor .gt. (-12*0.0254) .and. 293 & ycoor .lt. -0.044449)then ! this is the jet device 294 qold(nn,1)=97000 295 qold(nn,2)=0 296 qold(nn,3)=0 297 qold(nn,4)=0 298 qold(nn,5)=320 299 qold(nn,6)=1.825e-5 300 endif 301 302c... above, used for restart from M2M 303 304 enddo ! end of loop over all nshg points 305c........................... 306 307 return 308 end 309 310 311