1c----------------------------------------------------------------------- 2c 3c Initialize the predictor multicorrector (set up parameters) 4c 5c----------------------------------------------------------------------- 6 subroutine itrSetup ( y, acold ) 7 8 include "common.h" 9 10 real*8 y(nshg,ndof), acold(nshg,ndof) 11 12c 13c Define the Hulbert parameters 14c second order if between (and including) 0 and 1 otherwise backward Euler 15c 16 if( rhoinf(itseq).lt.0.or.rhoinf(itseq).gt.1) then ! backward Euler 17 almi = one 18 alfi = one 19 gami = one 20 ipred = 1 21 else !second order family 22 almi = (three-rhoinf(itseq))/(one+rhoinf(itseq))/two 23 alfi = one/(one+rhoinf(itseq)) 24 gami = pt5+almi-alfi 25 endif 26c 27c.... set the jacobian type 28c 29 Jactyp=0 30c 31 if(ipred.eq.4 .and. itseq .eq. 1 ) 32 & y=y-(one-alfi)*Delt(1)*CFLfl(1)*acold 33 34c 35c protect from ipred=4 and rhoi=0 36c 37 if(ipred.eq.4 .and. rhoinf(itseq) .eq. 0.0 ) ipred=3 38 39c 40c.... set the global time increment and the CFL data 41c 42 Dtgl = one / Delt(itseq) ! caution: inverse of time step 43 CFLfld = CFLfl(itseq) 44 CFLsld = CFLsl(itseq) 45 46 Dtgl = Dtgl / 1.0 ! check CFLfld 47 return 48 end 49 50 51c----------------------------------------------------------------------- 52c 53c Predict solution at time n+1 54c 55c----------------------------------------------------------------------- 56 subroutine itrPredict ( yold, acold, 57 & y, ac ) 58 59 include "common.h" 60 61 real*8 yold(nshg,ndof), acold(nshg,ndof), 62 & y(nshg,ndof), ac(nshg,ndof) 63 64c 65c.... set the solution at t_{n+alfi) and the acceleration at t_{n+almi) 66c 67 if(ipred.eq.1) then 68c 69c yn+1_pred=yn 70c 71 y = yold 72 ac = acold*(one-almi/gami) 73 endif 74c 75 if(ipred.eq.2) then 76c 77c an+1_pred=0 78c 79 y = yold+alfi/Dtgl*acold*(one-gami) 80 call itrBC (y, ac, iBC, BC, iper,ilwork) 81c Elaine-SPEBC 82 if((irscale.ge.0).and.(myrank.eq.master)) then 83 call genscale(y, x, iBC) 84 endif 85 ac = acold*(one-almi) 86 endif 87c 88 if(ipred.eq.3 ) then 89c 90c an+1_pred=an 91c 92 y = yold+alfi/Dtgl*acold 93 call itrBC (y, ac, iBC, BC, iper,ilwork) 94c Elaine-SPEBC 95 if((irscale.ge.0).and.(myrank.eq.master)) then 96 call genscale(y, x, iBC) 97 endif 98 ac = acold 99 endif 100c 101 if(ipred.eq.4 ) then ! protect from DC=4 rho=0 102c 103c same dV 104c 105 fct1=alfi/(one-alfi) 106 fct2=one-almi/gami 107 fct3=almi/gami/alfi*Dtgl 108 y = yold+fct1*(yold-y) 109 call itrBC (y, ac, iBC, BC, iper,ilwork) 110c Elaine-SPEBC 111 if((irscale.ge.0).and.(myrank.eq.master)) then 112 call genscale(y, x, iBC) 113 endif 114 ac = acold*fct2+(y-yold)*fct3 115 endif 116 117c 118 if (LCtime .eq. 0) time = time + Delt(itseq) 119c 120 return 121 end 122 123c----------------------------------------------------------------------- 124c 125c Correct solution at time n+1 126c 127c----------------------------------------------------------------------- 128 subroutine itrCorrect ( y, ac,yold, acold,Dy) 129 130 include "common.h" 131 132 real*8 y(nshg,ndof), ac(nshg,ndof), Dy(nshg,nflow) 133 real*8 yold(nshg,ndof), acold(nshg,ndof) 134c 135c 136c What we actually solved for here was delta y_{n+alpha_f). 137c This is the variable that we work with in the iteration loop 138c but we have another variable that must be updated 139c 140 y(:,1:3) = y(:,1:3) - Dy(:,2:4) 141 y(:,4) = y(:,4) - Dy(:,1) 142 y(:,5) = y(:,5) - Dy(:,5) 143c 144 fct1= (one-almi/gami) 145 fct2= almi*Dtgl/gami/alfi ! don't be confused by inverse dt=Dtgl 146 ac= acold*fct1+(y-yold)*fct2 147 return 148 end 149 150c----------------------------------------------------------------------- 151c 152c Correct solution at time n+1 153c 154c----------------------------------------------------------------------- 155 subroutine itrCorrectSclr (y, ac, yold, acold, Dyt) 156 157 include "common.h" 158 159 real*8 y(nshg,ndof), ac(nshg,ndof), Dyt(nshg) 160 real*8 yold(nshg,ndof), acold(nshg,ndof) 161c 162 is=5+isclr 163 y(:,is) = y(:,is) - Dyt(:) 164c 165 fct1= (one-almi/gami) 166 fct2= almi*Dtgl/gami/alfi ! don't be confused by inverse dt=Dtgl 167 ac(:,is)= acold(:,is)*fct1+(y(:,is)-yold(:,is))*fct2 168c 169 return 170 end 171c----------------------------------------------------------------------- 172c 173c Update solution at end of time step 174c 175c----------------------------------------------------------------------- 176 subroutine itrUpdate( yold, acold, 177 & y, ac ) 178 179 include "common.h" 180 181 real*8 yold(nshg,ndof), acold(nshg,ndof), 182 & y(nshg,ndof), ac(nshg,ndof) 183 184 185 if(iLset.eq.2 .and. isclr.eq.2) then 186c ... assigning the redistanced scalar to the first scalar 187c 188 y(:,6) = y(:,7) 189c 190c ... to calculate the right acceleration based on the updated value of the 191c ... scalar, Ref: Eq.25 in generalised alpha method paper by Prof.Jansen 192c 193 acold(:,6)= acold(:,6)*(1-1/gami) 194 & + (y(:,6)-yold(:,6))*Dtgl/gami 195 fct2=one/almi 196 fct3=one/alfi 197 acold(:,1:nflow) = acold(:,1:nflow) 198 & + (ac(:,1:nflow)-acold(:,1:nflow))*fct2 199 yold(:,1:nflow) = yold(:,1:nflow) 200 & + (y(:,1:nflow)-yold(:,1:nflow))*fct3 201 yold(:,6) = y(:,7) 202 acold(:,7)= zero 203 else 204 fct2=one/almi 205 fct3=one/alfi 206 acold = acold + (ac-acold)*fct2 207 yold = yold + (y-yold)*fct3 208 endif 209c 210 return 211 end 212 213 214 215 216 217 218 219