1c----------------------------------------------------------------------- 2c 3c This module conveys ramped BC data. 4c 5c----------------------------------------------------------------------- 6 module rampBC 7 integer kmaxinf, kmaxbot,kmaxtop,nfctop,nfcbot,ninflow 8 real*8 ctrl_a(3) 9 real*8 ctrl_factor(3) 10 end module 11 12 subroutine initBCprofileScale(vbc_prof,BC,yold,x) 13 use rampBC 14 include "common.h" 15 include "mpif.h" 16 include "auxmpi.h" 17 18! assuming three profiles to control (inlet, bottom FC and top FC) 19 real*8 vbc_prof(nshg,3), loc_ctrl_pr(3), ctrl_pr(3) 20 real*8 BC(nshg,ndofBC),yold(nshg,ndof),x(nshg,3) 21 nstp=nstep(1) 22! set constant factors for each profile ((V_peak/V_bar)*R*T/A) 23c 24c ALL COMPUTABLE WITH READ DATA BUT FOR NOW USER MUST COMPUTE FOR 25c EACH GEOMETRY + PROFILE 26C 27 ctrl_factor(1) = 1.00*Rgas*330.578/((2*0.587375)**2) ! inlet 28c ctrl_factor(1) = 1.00*Rgas*293/((0.0889*0.1143)) ! inlet 29c no chamber ctrl_factor(2) = 2.25*Rgas*293/((0.02441024589128274-0.008949756886432915)*(2*0.0508)) ! bottom FC 30 ctrl_factor(2) = 2.25*Rgas*317/((0.0168910950375201 31 & -0.01054110785157094)*(2*0.041275)) ! bottom FC 32 ctrl_factor(3) = 2.25*Rgas*317/(((-0.126612) 33 & -(-0.156330))*(2*0.0508)) ! top FC 34c ctrl_factor(2) = 2.25*Rgas*293/(lbot *wbot) 35c ctrl_factor(1) = CfactInl ! inlet 36c ctrl_factor(2) = CfactBot ! bottom FC 37c ctrl_factor(3) = CfactTop ! top FC 38c2D ctrl_factor(1) = 1.0*Rgas*273/0.046990500 ! inlet 39c2D ctrl_factor(2) = 1.5*Rgas*273/0.000618413 ! bottom FC 40c2D ctrl_factor(3) = 1.5*Rgas*273/0.001188690 ! top FC 41 42!set the first first step 43 ninflow=0 44 nfctop=0 45 nfcbot=0 46!assuming max at center line 47 kmaxinf=0 48 kmaxbot=0 49 kmaxtop=0 50 rmaxvinflow=-one 51 rmaxvbot=-one 52 rmaxvtop=-one 53c 54c find the peak value (and node number) of each inlet 55c SWITCH THIS TO USE SURFID's 56c 57 do kk=1,numnp 58 if(vbc_prof(kk,1).ne.zero) then 59 if(vbc_prof(kk,2).ne.zero) then ! both of these true means top FC 60 if(abs(vbc_prof(kk,1)).gt.rmaxvtop) then 61 kmaxtop=kk 62 rmaxvtop=abs(vbc_prof(kk,1)) 63 endif 64 nfctop=nfctop+1 65 else ! only x true means inflow 66 if(vbc_prof(kk,1).gt.rmaxvinflow) then 67 kmaxinf=kk 68 rmaxvinflow=vbc_prof(kk,1) 69 endif 70 ninflow=ninflow+1 71 endif 72 else if(vbc_prof(kk,2).ne.zero) then ! y true means low FC 73 if(vbc_prof(kk,2).gt.rmaxvbot) then 74 kmaxbot=kk 75 rmaxvbot=vbc_prof(kk,2) 76 endif 77 nfcbot=nfcbot+1 78 endif 79 enddo 80c 81c call BCprofileScale. Since called form init, this is going to be used to set 82c the yold value on the first step. Consequently, we will need to set istep to 83c -1 temporarily so that when it looks to the target value at the "end" of the c step it will actually be looking to time zero. Note we will reset it back to 84c zero after this call 85c 86 istep=-1 !offset the istep+1 increment so that it gets t0 right 87 call BCprofileScale(vbc_prof,BC,yold) 88 istep=0 !set back to correct value 89 return 90 end 91 92 subroutine BCprofileScale(vbc_prof,BC,yold) 93 use rampBC 94 use timedataC 95 include "common.h" 96 include "mpif.h" 97 include "auxmpi.h" 98! assuming three profiles to control (inlet, bottom FC and top FC) 99 real*8 vbc_prof(nshg,3), loc_ctrl_pr(3), ctrl_pr(3) 100 real*8 BC(nshg,ndofBC),yold(nshg,ndof),mdotNow(3) 101c 102c in some way find mdot for this way 103c 104 call rampedMdot(mdotNow) 105c call flowControl(mdotNow) 106 107c 108c compute the pressure from the current solution at the peak of profile 109c 110 111 loc_ctrl_pr(1) = yold(kmaxinf,4) ! inlet 112 loc_ctrl_pr(2) = yold(kmaxbot,4) ! bottom FC 113 loc_ctrl_pr(3) = yold(kmaxtop,4) ! top FC 114c 115c communicate it so that we are sure all agree 116c (in case the profile is split across procs) 117c 118 if(numpe > 1) then 119 call MPI_ALLREDUCE ( loc_ctrl_pr, ctrl_pr, 3, 120 & MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr ) 121 else 122 ctrl_pr = loc_ctrl_pr 123 endif 124 125 126 ctrl_a(1) = ctrl_factor(1)*mdotNow(1)/ctrl_pr(1) ! inlet 127 ctrl_a(2) = ctrl_factor(2)*mdotNow(2)/ctrl_pr(2) ! bottom FC 128 ctrl_a(3) = ctrl_factor(3)*mdotNow(3)/ctrl_pr(3) ! top FC 129 if(max(ninflow,max(nfctop,nfcbot)).gt. 0) then 130 do kk=1,numnp 131 if(vbc_prof(kk,1).ne.zero) then 132 if(vbc_prof(kk,2).ne.zero) then ! both of these true means top FC 133 if(ctrl_a(3).ge.zero) then 134 BC(kk,3)=ctrl_a(3)*vbc_prof(kk,1) 135 BC(kk,4)=ctrl_a(3)*vbc_prof(kk,2) 136 endif 137 else ! only x true means inflow 138 if(ctrl_a(1).ge.zero) then 139 BC(kk,3)=ctrl_a(1)*vbc_prof(kk,1) 140 endif 141 if(ctrl_a(1).eq.zero) then 142 BC(kk,7) = 0.0d0 ! scalar_1 set to zero for zero flow 143 endif 144 endif 145 else if(vbc_prof(kk,2).ne.zero) then ! y true means low FC 146 if(ctrl_a(2).ge.zero) then 147 BC(kk,4)=ctrl_a(2)*vbc_prof(kk,2) 148 endif 149 if(ctrl_a(2).eq.zero) then 150 BC(kk,7) = 0.0d0 ! scalar_1 set to zero for zero flow 151 endif 152 endif 153 enddo 154c if(ninflow.gt.0) write(myrank+2000,2000), lstep,kmaxinf,one, 155c & BC(kmaxinf,3),BC(kmaxinf,4),ctrl_a(1),ctrl_pr(1) 156c if(nfcbot .gt.0) write(myrank+2000,2000), lstep,kmaxbot,two, 157c & BC(kmaxbot,3),BC(kmaxbot,4),ctrl_a(2),ctrl_pr(2) 158c if(nfctop .gt.0) write(myrank+2000,2000), lstep,kmaxtop,three, 159c & BC(kmaxtop,3),BC(kmaxtop,4),ctrl_a(3),ctrl_pr(3) 160 endif 161 2000 format(i6,i6,5(e14.7,2x)) 162 return 163 end 164 165 166 subroutine rampedMdot(mdotNow) 167 include "common.h" !gives us rampmndot and nstep 168 169! assuming three profiles to control (inlet, bottom FC and top FC) 170 real*8 mdotNow(3),xi,cstep,rnstep 171 172 rnstep=nstep(1) 173 cstep=istep+1 174 xi=cstep/rnstep 175 176 mdotNow(:)=rampmdot(1,:)+xi*(rampmdot(2,:)-rampmdot(1,:)) 177 return 178 end 179