1c----------------------------------------------------------------------- 2c 3c Initialize the predictor multicorrector (set up parameters) 4c 5c Modified by Alberto Figueroa. Winter 2004 6c----------------------------------------------------------------------- 7 subroutine itrSetup ( y, acold ) 8 9 include "common.h" 10 11 real*8 y(nshg,ndof), acold(nshg,ndof) 12 13c 14c Define the Hulbert parameters 15c second order if between (and including) 0 and 1 otherwise backward Euler 16c 17 if( rhoinf(itseq).lt.0.or.rhoinf(itseq).gt.1) then ! backward Euler 18 almi = one 19 alfi = one 20 gami = one 21 ipred = 1 22 else !second order family 23 almi = (three-rhoinf(itseq))/(one+rhoinf(itseq))/two 24 alfi = one/(one+rhoinf(itseq)) 25 gami = pt5+almi-alfi 26 if(ideformwall.eq.1) then 27 betai=1.0 28 else 29 betai=0.0 30 endif 31 endif 32c 33c.... set the jacobian type 34c 35 Jactyp=0 36 if (impl(itseq) .eq. 3) then 37 Jactyp = 1 38 impl(itseq) = 2 39 endif 40c 41c.... same_Dy predictor special case 42c 43 if (ipred.eq.4 .and. itseq .eq. 1 ) then 44 y=y-(one-alfi)*Delt(1)*acold 45 if ( rhoinf(itseq) .eq. 0.0 ) then 46 ipred = 3 47 endif 48 endif 49c 50c.... set the global time increment and the CFL data 51c 52 Dtgl = one / Delt(itseq) ! caution: inverse of time step 53 CFLfld = CFLfl(itseq) 54 CFLsld = CFLsl(itseq) 55 56 return 57 end 58 59 60c----------------------------------------------------------------------- 61c 62c Predict solution at time n+1 63c 64c----------------------------------------------------------------------- 65 subroutine itrPredict (yold, y, acold, ac, uold, u, iBC) 66 67 include "common.h" 68 69 real*8 y(nshg,ndof), ac(nshg,ndof), 70 & u(nshg,nsd), yold(nshg,ndof), 71 & acold(nshg,ndof), uold(nshg,nsd) 72 dimension iBC(nshg) 73 74 75 if ( ipred.eq.1) then ! yn+1_pred=yn 76 fct = (gami-one)/gami 77 y = yold 78 ac = acold * fct 79 if(ideformwall.eq.1) then 80 do isd = 1, nsd 81 where(btest(iBC,13)) 82 u(:,isd) = uold(:,isd) + Delt(itseq)*yold(:,isd) + 83 & pt5*((gami-two*betai)/gami)* 84 & Delt(itseq)*Delt(itseq)*acold(:,isd) 85 endwhere 86 enddo 87 endif 88 endif 89c 90 if ( ipred.eq.2) then ! an+1_pred=0 91 y = yold + (one - gami)/Dtgl * acold 92 ac = 0.0 93 endif 94c 95 if(ipred.eq.3 ) then ! an+1_pred=an 96 y = yold+alfi*Delt(itseq)*acold 97 ac = acold 98 endif 99c 100 if ( ipred.eq.4 ) then ! protect from DC=4 rho=0, same dV 101 fct1 = alfi/(one-alfi) 102 fct2 = one-almi/gami 103 fct3 = almi/gami/alfi*Dtgl 104 y = yold+fct1*(yold-y) 105 ac = acold*fct2+(y-yold)*fct3 106 endif 107c 108 109 return 110 end 111 112c----------------------------------------------------------------------- 113c 114c Correct solution at time n+1 115c 116c----------------------------------------------------------------------- 117 subroutine itrCorrect ( y, ac, u, solinc ,iBC) 118 119 include "common.h" 120 121 real*8 y(nshg,ndof), ac(nshg,ndof), 122 & u(nshg,nsd), solinc(nshg,4) 123 dimension iBC(nshg) 124 125 fct1 = gami*Delt(itseq) 126 fct2 = gami*alfi*Delt(itseq) 127 128 y(:,1:3) = y(:,1:3) + fct1 * solinc(:,1:3) 129 y(:,4 ) = y(:,4 ) + fct2 * solinc(:,4 ) 130 ac(:,1:3) = ac(:,1:3) + solinc(:,1:3) 131 if(ideformwall.eq.1) then 132 do isd = 1,nsd 133 where (btest(iBC,13)) 134 u(:,isd) = u(:,isd) + 135 & Delt(itseq)*Delt(itseq)*betai*solinc(:,isd) 136 endwhere 137 enddo 138 endif 139c 140 return 141 end 142 143c----------------------------------------------------------------------- 144c 145c Correct solution at time n+1 146c 147c----------------------------------------------------------------------- 148 subroutine itrCorrectSclr ( y, ac, solinc ) 149 150 include "common.h" 151 152 real*8 y(nshg,ndof), ac(nshg,ndof), 153 & solinc(nshg) 154 155 fct1 = gami*Delt(itseq) 156 is=5+isclr 157 y(:,is) = y(:,is) + fct1 * solinc(:) 158 ac(:,is) = ac(:,is) + solinc(:) 159c 160 return 161 end 162 163c----------------------------------------------------------------------- 164c 165c Correct solution at time n+1, protecting against negative values 166c 167c----------------------------------------------------------------------- 168 subroutine itrCorrectSclrPos ( y, ac, solinc ) 169 170 include "common.h" 171 172 real*8 y(nshg,ndof), ac(nshg,ndof), 173 & solinc(nshg), posinc(nshg) 174 175 fct1 = gami*Delt(itseq) 176 if(fct1.eq.0) then 177 call itrCorrectSclr( y, ac, solinc) 178 else 179 turbUpdFct = 1.0 180 updFct = turbUpdFct 181 updFct2 = 1.0 182c if(topHdTscFlag.and.PRM_TCS_PARTIAL_UPDATE)then 183c updFct2 = max(MTH_SQRT_EPS,min(topHdUpdFct,1.0)) 184c endif 185 fctCr = fct1*updFct*updFct2 186 c0 = 0.01 187 c1 = -(1.0-c0)/fctCr 188 is=5+isclr 189c if(any(updFct*updFct2*solinc(:).lt.c1*y(:,is))) then 190c write(*,*) 'INSTEAD OF GETTING NEGATIVE FIELD' 191c write(*,*) 'BROUGHT FIELD DOWN TO 1 PERCENT' 192c write(*,*) 'FOR SCALAR NUMBER ',isclr 193c write(*,*) '(SEE itrCorrectSclr in itrPC.f)' 194c endif 195 posinc = max(updFct*updFct2*solinc,c1*y(:,is)) 196 y(:,is) = y(:,is) + fct1 * posinc(:) 197 ac(:,is) = ac(:,is) + posinc(:) 198 endif 199c 200 return 201 end 202 203 204c----------------------------------------------------------------------- 205c 206c Compute solution and acceleration at n+alpha 207c 208c----------------------------------------------------------------------- 209 subroutine itrYAlpha ( uold, yold, acold, 210 & u, y, ac, 211 & uAlpha, yAlpha, acAlpha ) 212 213c use readarrays !reads in uold and acold 214 include "common.h" 215 216 real*8 yold(nshg,ndof), acold(nshg,ndof), 217 & y(nshg,ndof), ac(nshg,ndof), 218 & yAlpha(nshg,ndof), acAlpha(nshg,ndof), 219 & u(nshg,nsd), uold(nshg,nsd), 220 & uAlpha(nshg,nsd) 221 222 acAlpha(:,4) = zero !pressure acceleration is never used but.... 223 224 yAlpha (:,1:3) = yold(:,1:3) 225 & + alfi * (y(:,1:3) - yold(:,1:3)) 226 227 acAlpha(:,1:3) = acold(:,1:3) 228 & + almi * (ac(:,1:3) - acold(:,1:3)) 229 230 yAlpha (:,4 ) = y(:,4) 231 232 if(ideformwall.eq.1) uAlpha (:,1:3) = uold(:,1:3) 233 & + alfi * (u(:,1:3) - uold(:,1:3)) 234 235 if(ndof.ge.5) then 236c 237c Now take care of temperature, turbulence, what have you 238c 239 240 241 yAlpha (:,5:ndof ) = yold(:,5:ndof) 242 & + alfi * (y(:,5:ndof) - yold(:,5:ndof)) 243 acAlpha(:,5:ndof ) = acold(:,5:ndof) 244 & + almi * (ac(:,5:ndof) - acold(:,5:ndof)) 245 246 endif 247 return 248 end 249 250c----------------------------------------------------------------------- 251c 252c Update solution at end of time step 253c 254c----------------------------------------------------------------------- 255 subroutine itrUpdate( yold, acold, uold, 256 & y, ac, u ) 257 258c use readarrays !reads in uold and acold 259 260 include "common.h" 261 262 real*8 yold(nshg,ndof), acold(nshg,ndof), 263 & y(nshg,ndof), ac(nshg,ndof), 264 & u(nshg,nsd), uold(nshg,nsd) 265 266 267!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 268! Allow for redistance of the level set function, AD 5/8/00 269!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 270 271! if we are re-distancing the levelset function the we need 272! to update the levelset variable with the corrected value 273! stored in the "next" scalar position. 274 275 if(iLSet .eq.2 .and. isclr .eq. 2) then 276 277 y(:,6) = y(:,7) ! will need to also fix the yold ac etc.????? 278 ac(:,7)=zero 279 endif 280 281 282 283 yold = y 284 acold = ac 285 if(ideformwall.eq.1) uold = u 286 287 return 288 end 289 290