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 if(mod(iBC(i),1024).eq.ibcmatch) then 123 iupper=0 124 do j=2,ninterp 125 if(bcinterp(j,1).gt.d2wall(i)) then !bound found 126 xi=(d2wall(i)-bcinterp(j-1,1))/ 127 & (bcinterp(j,1)-bcinterp(j-1,1)) 128 iupper=j 129 exit 130 endif 131 enddo 132 if(iupper.eq.0) then 133 write(*,*) 'failure in finterp, ynint=',d2wall(i) 134 stop 135 endif 136 if(lstep.eq.0) then ! apply this inlet bl as an IC if on step 0 137 y(i,1:3)=(xi*bcinterp(iupper,4:6) 138 & +(one-xi)*bcinterp(iupper-1,4:6)) 139 y(i,4)=(xi*bcinterp(iupper,2) 140 & +(one-xi)*bcinterp(iupper-1,2)) 141 y(i,5)=(xi*bcinterp(iupper,3) 142 & +(one-xi)*bcinterp(iupper-1,3)) 143 endif 144 if(point2x(i,1).lt.1.0e-5 .and. mod(iBC(i),1024).eq.ibcmatch) then 145 do j=1,ndof 146 if(interp_mask(j).ne.zero) then 147 BC(i,j)=(xi*bcinterp(iupper,j+1) 148 & +(one-xi)*bcinterp(iupper-1,j+1)) 149 endif 150 enddo 151 endif 152 enddo 153 endif 154c$$$$$$$$$$$$$$$$$$$$ 155 156!====================================================================== 157!Modifications for Duct. Need to encapsulate in a function call. 158 !specify an initial eddy viscosity ramp 159 if(isetEVramp .gt. 0) then 160 if(myrank .eq. 0) then 161 write(*,*) "Setting eddy viscosity ramp with:" 162 write(*,*) " - ramp X min = ", EVrampXmin 163 write(*,*) " - ramp X max = ", EVrampXmax 164 write(*,*) " - EV min = ", EVrampMin 165 write(*,*) " - EV max = ", EVrampMax 166 endif 167 168 x1 = EVrampXmin !stuff in a shorter variable name to 169 x2 = EVrampXmax !make the formulas cleaner 170 !Newton Divided differences to generate a polynomial of 171 !the form P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3)) 172 !satisfying P(x1) = EVrampMin, P(x2) = EVrampMax, 173 ! P'(x1) = 0, and P'(x2) = 0 174 175 c0 = EVrampMin 176 c1 = 0 !zero derivative 177 c2 = (EVrampMax - EVrampMin)/(x2 - x1) 178 c3 = 0 !zero derivative 179 c3 = (c3 - c2)/(x2 - x1) 180 c2 = (c2 - c1)/(x2 - x1) 181 c3 = (c3 - c2)/(x2 - x1) 182 183 do nn = 1,nshg 184 if(y(nn,6) .eq. 0) cycle !don't change wall boundary conditions, should be iTurbWall == 1 185 186 if(point2x(nn,1) .gt. EVrampXmax) then !downstream of the ramp 187 y(nn,6) = EVrampMax 188 elseif(point2x(nn,1) .gt. EVrampXmin) then !and x(:,1) <= EVrampXmax 189 190 !P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3)) 191 ! = c0 + x*(c1 + x*(c2 - (x2 - x1)*c3 + x*c3 192 y(nn,6) = c0 + (point2x(nn,1) - x1)*( 193 & c1 + (point2x(nn,1) - x1)*( 194 & (c2 - (x2 - x1)*c3) + (point2x(nn,1) - x1)*c3)) 195 else 196 y(nn,6) = EVrampMin 197 endif 198 enddo 199 endif 200!End modifications for Duct 201!====================================================================== 202c 203c 204c.... call the semi-discrete predictor multi-corrector iterative driver 205c 206 call itrdrv (y, ac, 207 & uold, point2x, 208 & iBC, BC, 209 & point2iper, point2ilwork, shp, 210 & shgl, shpb, shglb, 211 & point2ifath, velbar, point2nsons ) 212c 213c.... return 214c 215c 216c.... stop CPU-timer 217c 218CAD call timer ('End ') 219c 220c.... close echo file 221c 222 223!MR CHANGE 224 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 225 if(myrank.eq.0) then 226 write(*,*) 'process - before closing iecho' 227 endif 228!MR CHANGE END 229 230 close (iecho) 231 232!MR CHANGE 233 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 234 if(myrank.eq.0) then 235 write(*,*) 'process - after closing iecho' 236 endif 237!MR CHANGE END 238 239 240c 241c.... end of the program 242c 243CAD write(6,*) 'Life: ', second(0) - ttim(100) 244 deallocate(point2iper) 245 if(numpe.gt.1) then 246 call Dctypes(point2ilwork(1)) 247 endif 248 deallocate(point2ilwork) 249 deallocate(point2x) 250 deallocate(point2nsons) 251 deallocate(point2ifath) 252 deallocate(uold) 253 deallocate(wnrm) 254 deallocate(otwn) 255 call finalizeDtN 256 call clearper 257 call finalizeNABI 258 259 return 260 end 261 262 263