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