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*293/((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*293/((0.0168910950375201 31 & -0.01054110785157094)*(2*0.041275)) ! bottom FC 32 ctrl_factor(3) = 2.25*Rgas*293/(((-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=zero 51 rmaxvbot=zero 52 rmaxvtop=zero 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 callled form init, this is going to be used to setc the yold value on the first step. Consequently, we will need to set istep toc -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 toc zero after this call 82c 83 istep=-1 !offset the istep+1 increment so that it gets t0 right 84 call BCprofileScale(vbc_prof,BC,yold) 85 istep=0 !set back to correct value 86 return 87 end 88 89 subroutine BCprofileScale(vbc_prof,BC,yold) 90 use rampBC 91 use timedata 92 include "common.h" 93 include "mpif.h" 94 include "auxmpi.h" 95! assuming three profiles to control (inlet, bottom FC and top FC) 96 real*8 vbc_prof(nshg,3), loc_ctrl_pr(3), ctrl_pr(3) 97 real*8 BC(nshg,ndofBC),yold(nshg,ndof),mdotNow(3) 98c 99c in some way find mdot for this way 100c 101 call rampedMdot(mdotNow) 102c call flowControl(mdotNow) 103 104c 105c compute the pressure from the current solution at the peak of profile 106c 107 108 loc_ctrl_pr(1) = yold(kmaxinf,4) ! inlet 109 loc_ctrl_pr(2) = yold(kmaxbot,4) ! bottom FC 110 loc_ctrl_pr(3) = yold(kmaxtop,4) ! top FC 111c 112c communicate it so that we are sure all agree 113c (in case the profile is split across procs) 114c 115 if(numpe > 1) then 116 call MPI_ALLREDUCE ( loc_ctrl_pr, ctrl_pr, 3, 117 & MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr ) 118 else 119 ctrl_pr = loc_ctrl_pr 120 endif 121 122 123 ctrl_a(1) = ctrl_factor(1)*mdotNow(1)/ctrl_pr(1) ! inlet 124 ctrl_a(2) = ctrl_factor(2)*mdotNow(2)/ctrl_pr(2) ! bottom FC 125 ctrl_a(3) = ctrl_factor(3)*mdotNow(3)/ctrl_pr(3) ! top FC 126 if(max(ninflow,max(nfctop,nfcbot)).gt. 0) then 127 do kk=1,numnp 128 if(vbc_prof(kk,1).ne.zero) then 129 if(vbc_prof(kk,2).ne.zero) then ! both of these true means top FC 130 BC(kk,3)=ctrl_a(3)*vbc_prof(kk,1) 131 BC(kk,4)=ctrl_a(3)*vbc_prof(kk,2) 132 else ! only x true means inflow 133 BC(kk,3)=ctrl_a(1)*vbc_prof(kk,1) 134 endif 135 else if(vbc_prof(kk,2).ne.zero) then ! y true means low FC 136 BC(kk,4)=ctrl_a(2)*vbc_prof(kk,2) 137 endif 138 enddo 139 if(ninflow.gt.0) write(myrank+2000,2000) 140 & BC(kmaxinf,3),BC(kmaxinf,4),ctrl_a(1),ctrl_pr(1) 141 if(nfcbot .gt.0) write(myrank+2000,2000) 142 & BC(kmaxbot,3),BC(kmaxbot,4),ctrl_a(2),ctrl_pr(2) 143 if(nfctop .gt.0) write(myrank+2000,2000) 144 & BC(kmaxtop,3),BC(kmaxtop,4),ctrl_a(3),ctrl_pr(3) 145 endif 146 2000 format(4(e14.7,2x)) 147 return 148 end 149 150 151 subroutine rampedMdot(mdotNow) 152 include "common.h" !gives us rampmndot and nstep 153 154! assuming three profiles to control (inlet, bottom FC and top FC) 155 real*8 mdotNow(3),xi,cstep,rnstep 156 157 rnstep=nstep(1) 158 cstep=istep+1 159 xi=cstep/rnstep 160 161 mdotNow(:)=rampmdot(1,:)+xi*(rampmdot(2,:)-rampmdot(1,:)) 162 return 163 end 164