1 subroutine proces 2c 3c---------------------------------------------------------------------- 4c 5c This subroutine generates the problem data and calls the solution 6c driver. 7c 8c 9c Zdenek Johan, Winter 1991. (Fortran 90) 10c---------------------------------------------------------------------- 11c 12 use readarrays ! used to access x, iper, ilwork 13 use turbsa ! used to access d2wall 14 use dtnmod 15 use periodicity 16 use pvsQbi 17 include "common.h" 18 include "mpif.h" 19c 20c arrays in the following 2 lines are now dimensioned in readnblk 21c dimension x(numnp,nsd) 22c dimension iper(nshg), ilwork(nlwork) 23c 24 dimension y(nshg,ndof), 25 & iBC(nshg), 26 & BC(nshg,ndofBC), 27 & ac(nshg,ndof) 28c 29c.... shape function declarations 30c 31 dimension shp(MAXTOP,maxsh,MAXQPT), 32 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 33 & shpb(MAXTOP,maxsh,MAXQPT), 34 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 35c 36c stuff for dynamic model s.w.avg and wall model 37c 38c dimension ifath(numnp), velbar(nfath,nflow), 39c & nsons(nfath) 40 41 dimension velbar(nfath,nflow) 42c 43c stuff to interpolate profiles at inlet 44c 45 real*8 bcinterp(100,ndof+1),interp_mask(ndof) 46 logical exlog 47 48 !Duct 49 real*8 c0, c1, c2, c3, x1, x2 50 integer nn 51 52c if ((irscale .ge. 0) .and. (myrank.eq.master)) then 53c call setSPEBC(numnp, point2nfath, nsonmax) 54c endif 55c 56c.... generate the geometry and boundary conditions data 57c 58 call gendat (y, ac, point2x, 59 & iBC, BC, 60 & point2iper, point2ilwork, shp, 61 & shgl, shpb, shglb, 62 & point2ifath, velbar, point2nsons ) 63 call setper(nshg) 64 call perprep(iBC,point2iper,nshg) 65 if (iLES/10 .eq. 1) then 66 call keeplhsG ! Allocating space for the mass (Gram) matrix to be 67 ! used for projection filtering and reconstruction 68 ! of the strain rate tensor. 69 70 call setrls ! Allocating space for the resolved Leonard stresses 71c See bardmc.f 72 endif 73c 74c.... time averaged statistics 75c 76 if (ioform .eq. 2) then 77 call initStats(point2x, iBC, point2iper, point2ilwork) 78 endif 79c 80c.... RANS turbulence model 81c 82! if (iRANS .lt. 0) then 83 call initTurb( point2x ) 84! endif 85c 86c.... p vs. Q boundary 87c 88 call initNABI( point2x, shpb ) 89c 90c.... check for execution mode 91c 92 if (iexec .eq. 0) then 93 lstep = 0 94 call restar ('out ', y ,ac) 95 return 96 endif 97c 98c.... initialize AutoSponge 99c 100 if(matflg(5,1).ge.4) then ! cool case (sponge) 101 call initSponge( y,point2x) 102 endif 103c 104c 105c.... adjust BC's to interpolate from file 106c 107 108 inquire(file="inlet.dat",exist=exlog) 109 if(exlog) then 110 open (unit=654,file="inlet.dat",status="old") 111 read(654,*) ninterp,ibcmatch,(interp_mask(j),j=1,ndof) 112! with no surfID control, this will also get applied to isothermal walls but! usually ok and can be extended if not. 113 114 do i=1,ninterp 115 read(654,*) (bcinterp(i,j),j=1,ndof+1) ! distance to wall+ 116 ! ndof but note order of BC's 117 ! p,t,u,v,w,scalars. Also note that must be in 118 ! increasing distance from the wall order. 119 120 enddo 121 do i=1,nshg ! only correct for linears at this time 122 iupper=0 123 do j=2,ninterp 124 if(bcinterp(j,1).gt.d2wall(i)) then !bound found 125 xi=(d2wall(i)-bcinterp(j-1,1))/ 126 & (bcinterp(j,1)-bcinterp(j-1,1)) 127 iupper=j 128 exit 129 endif 130 enddo 131 if(iupper.eq.0) then 132 write(*,*) 'failure in finterp, ynint=',d2wall(i) 133 stop 134 endif 135 if(lstep.eq.0) then ! apply this inlet bl as an IC if on step 0 136 y(i,1:3)=(xi*bcinterp(iupper,4:6) 137 & +(one-xi)*bcinterp(iupper-1,4:6)) 138 y(i,4)=(xi*bcinterp(iupper,2) 139 & +(one-xi)*bcinterp(iupper-1,2)) 140 y(i,5)=(xi*bcinterp(iupper,3) 141 & +(one-xi)*bcinterp(iupper-1,3)) 142 endif 143 if(point2x(i,1).lt.1.0e-5 .and. mod(iBC(i),1024).eq.ibcmatch) then 144 do j=1,ndof 145 if(interp_mask(j).ne.zero) then 146 BC(i,j)=(xi*bcinterp(iupper,j+1) 147 & +(one-xi)*bcinterp(iupper-1,j+1)) 148 endif 149 enddo 150 endif 151 enddo 152 endif 153c$$$$$$$$$$$$$$$$$$$$ 154 155!====================================================================== 156!Modifications for Duct. Need to encapsulate in a function call. 157 !specify an initial eddy viscosity ramp 158 if(isetEVramp .gt. 0) then 159 if(myrank .eq. 0) then 160 write(*,*) "Setting eddy viscosity ramp with:" 161 write(*,*) " - ramp X min = ", EVrampXmin 162 write(*,*) " - ramp X max = ", EVrampXmax 163 write(*,*) " - EV min = ", EVrampMin 164 write(*,*) " - EV max = ", EVrampMax 165 endif 166 167 x1 = EVrampXmin !stuff in a shorter variable name to 168 x2 = EVrampXmax !make the formulas cleaner 169 !Newton Divided differences to generate a polynomial of 170 !the form P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3)) 171 !satisfying P(x1) = EVrampMin, P(x2) = EVrampMax, 172 ! P'(x1) = 0, and P'(x2) = 0 173 174 c0 = EVrampMin 175 c1 = 0 !zero derivative 176 c2 = (EVrampMax - EVrampMin)/(x2 - x1) 177 c3 = 0 !zero derivative 178 c3 = (c3 - c2)/(x2 - x1) 179 c2 = (c2 - c1)/(x2 - x1) 180 c3 = (c3 - c2)/(x2 - x1) 181 182 do nn = 1,nshg 183 if(y(nn,6) .eq. 0) cycle !don't change wall boundary conditions, should be iTurbWall == 1 184 185 if(point2x(nn,1) .gt. EVrampXmax) then !downstream of the ramp 186 y(nn,6) = EVrampMax 187 elseif(point2x(nn,1) .gt. EVrampXmin) then !and x(:,1) <= EVrampXmax 188 189 !P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3)) 190 ! = c0 + x*(c1 + x*(c2 - (x2 - x1)*c3 + x*c3 191 y(nn,6) = c0 + (point2x(nn,1) - x1)*( 192 & c1 + (point2x(nn,1) - x1)*( 193 & (c2 - (x2 - x1)*c3) + (point2x(nn,1) - x1)*c3)) 194 else 195 y(nn,6) = EVrampMin 196 endif 197 enddo 198 endif 199!End modifications for Duct 200!====================================================================== 201c 202c 203c.... call the semi-discrete predictor multi-corrector iterative driver 204c 205 call itrdrv (y, ac, 206 & uold, point2x, 207 & iBC, BC, 208 & point2iper, point2ilwork, shp, 209 & shgl, shpb, shglb, 210 & point2ifath, velbar, point2nsons ) 211c 212c.... return 213c 214c 215c.... stop CPU-timer 216c 217CAD call timer ('End ') 218c 219c.... close echo file 220c 221 222!MR CHANGE 223 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 224 if(myrank.eq.0) then 225 write(*,*) 'process - before closing iecho' 226 endif 227!MR CHANGE END 228 229 close (iecho) 230 231!MR CHANGE 232 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 233 if(myrank.eq.0) then 234 write(*,*) 'process - after closing iecho' 235 endif 236!MR CHANGE END 237 238 239c 240c.... end of the program 241c 242CAD write(6,*) 'Life: ', second(0) - ttim(100) 243 deallocate(point2iper) 244 if(numpe.gt.1) then 245 call Dctypes(point2ilwork(1)) 246 endif 247 deallocate(point2ilwork) 248 deallocate(point2x) 249 deallocate(point2nsons) 250 deallocate(point2ifath) 251 deallocate(uold) 252 deallocate(wnrm) 253 deallocate(otwn) 254 call finalizeDtN 255 call clearper 256 call finalizeNABI 257 258 return 259 end 260 261 262