1 module MachControl 2 logical exMC !Mach Control enable / Does the MC file 3 !exist 4 5 integer :: MC_iPP !Probe Point index 6 integer :: MC_surfID !Surface ID to identify where to apply 7 ! the boundary condition 8 9 real*8 :: MC_dt !simulaiton time step 10 real*8 :: MC_MtTarget !Target throat Mach number 11 real*8 :: MC_dMtCO !M_throat Cut Off (frequency) (for 12 ! filtering the state) 13 real*8 :: MC_ulim(2) !u_inlet velocity limit 14 !u_inlet(1) lower limit 15 !u_inlet(2) upper limit 16 !States 17 real*8 :: MC_Mt !actual Mach number 18 real*8 :: MC_dMtdt !time derivative of MC_dMt 19 real*8 :: MC_dMt !filtered state deviation from target 20 real*8 :: MC_IdMt !integreal of MC_dMt 21 22 !Gains 23 real*8 :: MC_Kd !derivative gain 24 real*8 :: MC_Kp !proportional gain 25 real*8 :: MC_KI !integral gain 26 27 !Variables to apply the boundary condition 28 integer, allocatable, dimension(:) :: MC_imapInlet(:) 29 integer :: MC_icountInlet !number of inlet boundary points 30 real*8, allocatable, dimension(:) :: MC_vscale(:) !scaling factor for BL 31 real*8 :: MC_BLdelta !boundary layer thickness 32 real*8 :: MC_uInlet !x-velocity at the inlet 33 34 end module 35 36 subroutine MC_init(dt, tstep, BC) 37 use MachControl 38 use turbSA 39 include "common.h" 40 include "mpif.h" 41 42 dimension :: BC(nshg, ndofBC) 43 logical :: MC_existDir 44 real*8 :: dt 45 real*8 :: vBC, vBCg 46 integer :: tstep 47 vBCg = 0 48 49 if(myrank == master) then 50 inquire(file="./MachControl/.", exist=MC_existDir) 51 if(.not. MC_existDir) then 52 call system("mkdir ./MachControl") !Doesn't seem to work on BGQ. 53 endif 54 endif 55 56 call MC_readState(tstep) 57 58 if(exMC) then 59 !find the node associated with isetBlowerID_Duct and write the 60 !scale vector to linear ramp velocity near the walls. 61 allocate(MC_imapInlet(nshg)) 62 allocate(MC_vscale(nshg)) 63 call sfID2np(MC_surfID, MC_icountInlet, MC_imapInlet) 64 65 if(MC_BLdelta > 0) then 66 do i = 1, MC_icountInlet 67 imapped = MC_imapInlet(i) 68 69 MC_vscale(i) = min(one, d2wall(imapped)/MC_BLdelta) 70 !MC_vscale(i) = d2wall(imapped)/MC_BLdelta 71 !if(MC_vscale(i) > 1) MC_vscale(i) = 1 72 73 !Also look for the inlet velocity. This will get 74 !overwritten in the next time step, but it's nice to know 75 !what it started with. 76 vBCg = max(vBCg, BC(imapped, 3)) 77 enddo 78 endif 79 80 !Look for the maximum global inlet velocity among all parts. 81 if(numpe > 1) then 82 call MPI_ALLREDUCE(vBCg, MC_uInlet, 1, 83 & MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr) 84 endif 85 86 MC_dt = dt 87 endif 88 end subroutine 89 90 91 subroutine MC_updateState() 92 !Updates the derivative, proportional and integral values used 93 !when calculating MC_uInlet. 94 !Note: This function uses the most recent values in the TimeData 95 !buffer vartsbuff. TD_bufferData should be called prior to this 96 !subroutine. 97 98 use MachControl 99 use timedataC 100 include "common.h" 101 include "mpif.h" 102 103 real*8 :: u2, T 104 real*8 :: Mt, dMt, MC_alpha 105 logical :: bangbang 106 107 if(myrank == master) then !MC_* are not set on other processors 108 T = vartsbuff(MC_iPP, 5, ivartsBuff) 109 u2 = vartsbuff(MC_iPP, 2, ivartsBuff)**2 + 110 & vartsbuff(MC_iPP, 3, ivartsBuff)**2 + 111 & vartsbuff(MC_iPP, 4, ivartsBuff)**2 112 113 MC_Mt = sqrt(u2/(1.4*Rgas*T)) 114 dMt = MC_Mt - MC_MtTarget 115 116 !Note: MC_uInlet should always equal the previous time step; 117 !on the first time step of a run, MC_uInlet is computed based 118 !off of the maximum inlet velocity in the solution. 119 bangbang = MC_uInlet <= MC_uLim(1) .or. 120 & MC_uInlet >= MC_uLim(2) 121 122 !Low pass filter dMt 123 MC_alpha = MC_dt/(MC_dMtCO + MC_dt) 124 dMt = MC_alpha*dMt + (1 - MC_alpha)*MC_dMt 125 126 !Calculate derivatives and integrals 127 if(.not. bangbang) then 128 !Only integrate when not hitting the limits. This introduces 129 !a nonlinearity in the control scheme which minimizes the 130 !overshoot when starting from zero initial conditions 131 MC_IdMt = (dMt + MC_dMt)*MC_dt/2 + MC_IdMt !trapezoidal rule 132 endif 133 MC_dMtdt = (dMt - MC_dMt)/MC_dt !first order BE FD derivative 134 MC_dMt = dMt !update the actual state 135 endif 136 end subroutine 137 138 139 subroutine MC_applyBC(BC) 140 use MachControl 141 include "common.h" 142 include "mpif.h" 143 144 dimension :: BC(nshg, ndofBC) 145 real*8 :: vInlet 146 integer :: imapped 147 real*8 :: tmp 148 149 if(myrank .eq. master) then 150 MC_uInlet = MC_Kd*MC_dMtdt + 151 & MC_Kp*MC_dMt + 152 & MC_KI*MC_IdMt 153 154 !clip to [MC_uLim(1),MC_uLim(2)] 155 if(MC_uLim(1) < MC_ulim(2)) then 156 MC_uInlet = min(max(MC_uInlet, MC_uLim(1)), MC_uLim(2)) 157 endif 158 endif 159 if(numpe .gt. 1) 160 & call MPI_Bcast(MC_uInlet, 1, MPI_DOUBLE_PRECISION, 161 & master, MPI_COMM_WORLD, ierr) 162 163 if(MC_BLdelta .gt. 0) then !Apply a velocity profile with a BL 164 do i = 1, MC_icountInlet 165 imapped = MC_imapInlet(i) 166 BC(imapped,3) = MC_uInlet*MC_vscale(i) 167 enddo 168 else !Apply a uniform inlet velocity. 169 do i = 1, MC_icountInlet 170 BC(MC_imapInlet(i), 3) = MC_uInlet 171 enddo 172 endif 173 174 end subroutine 175 176 177 subroutine MC_printState() 178 !Prints the Mach Control states into the output 179 180 use MachControl 181 include "common.h" 182 183 if(myrank == master) then 184 write(*, "('M_throat: ',F8.6,' dM_throat: ',F8.6)") 185 & MC_Mt, MC_Mt - MC_MtTarget 186 write(*,"('dM/dt = ',E14.6,' M = ',E14.6,' I(M) =',E14.6)") 187 & MC_dMtdt, MC_dMt, MC_IdMt 188 write(*,"('Kd = ',E14.6,' Kp = ',E14.6,' KI =',E14.6)") 189 & MC_Kd, MC_Kp, MC_KI 190 write(*, "('u_inlet: ',F9.6)") MC_uInlet 191 endif 192 end subroutine 193 194 subroutine MC_readState(tstep) 195 !Reads the intput / state file for MachControl. 196 !Written by: Nicholas Mati 2014-04-19 197 !Revision History: 198 ! - 2014-04-19: created 199 ! 200 !The file is expected to be of the form: 201 ! MC_iPP MC_surfID MC_MtTarget 202 ! MC_dMtdt MC_dMt MC_IdMt 203 ! MC_Kd MC_Kp MC_KI 204 ! MC_dMtCO MC_BLdelta 205 206 use MachControl 207 include "common.h" 208 include "mpif.h" 209 210 integer :: tstep 211 character(len=60) :: fname 212 213 if(myrank == master) then 214 write(fname, "('./MachControl/MachControl.',I0,'.dat')") tstep 215 inquire(file=fname, exist=exMC) 216 217 if(exMC) then 218 open(unit=1002, file=fname, status='old') 219 220 read(1002, *) MC_iPP, MC_surfID 221 read(1002, *) MC_Mt, MC_MtTarget 222 read(1002, *) MC_dMtdt, MC_dMt, MC_IdMt 223 read(1002, *) MC_Kd, MC_Kp, MC_KI 224 read(1002, *) MC_dMtCO, MC_BLdelta 225 read(1002, *) MC_uLim(1), MC_uLim(2) 226 close(1002) 227 228 !If d2wall is not available, then set the inlet BL thickness 229 !to zero to disable wall scalling. 230 if(iRANS >= 0) MC_BLdelta = 0 231 endif 232 endif 233 234 !Only a few values are needed across other processors. 235 if(numpe > 1) then 236 call MPI_BCAST(MC_BLdelta, 1, MPI_DOUBLE_PRECISION, 237 & master, MPI_COMM_WORLD, ierr) 238 call MPI_BCAST(MC_surfID, 1, MPI_INTEGER, 239 & master, MPI_COMM_WORLD, ierr) 240 call MPI_BCAST(exMC, 1, MPI_INTEGER, 241 & master, MPI_COMM_WORLD, ierr) 242 endif 243 244 end subroutine 245 246 subroutine MC_writeState(tstep) 247 !Writes a file with the name MachControl.[ts].dat that can later 248 !be read in and restarted from. 249 250 use MachControl 251 include "common.h" 252 253 integer tstep 254 character(len=60) :: fname 255 256 if(myrank == master) then 257 write(fname, "('./MachControl/MachControl.',I0,'.dat')") tstep 258 open(unit=1002, file=fname, status='unknown') 259 260 write(1002, "(2(i9))") MC_iPP, MC_surfID 261 write(1002, "(2(F18.14))") MC_Mt, MC_MtTarget 262 write(1002, "(3(F18.14))") MC_dMtdt, MC_dMt, MC_IdMt 263 write(1002, "(3(F18.10))") MC_Kd, MC_Kp, MC_KI 264 write(1002, "(2(F18.14))") MC_dMtCO, MC_BLdelta 265 write(1002, "(2(F18.14))") MC_uLim(1), MC_uLim(2) 266 close(1002) 267 endif 268 end subroutine 269 270 subroutine MC_finalize() 271 272 use MachControl 273 274 deallocate(MC_imapInlet) 275 deallocate(MC_vscale) 276 end subroutine 277 278 279