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