subroutine outgoing_branch_cold_circulator(me,mdot,VarIn,deriv_wr_both_br,ssOut,FluxOut)
type(circulator_t), intent(inout) :: me
real(dp), intent(in) :: mdot,VarIn
type(derivatives_involving_both_branches_t), intent(in) :: deriv_wr_both_br
real(dp), intent(out) :: ssOut
type(flux_and_derivatives_for_cold_circulator_t), intent(out) :: FluxOut
real(dp) :: Ror,Pr,Ur,Cr,Zr,pmin,pmax,f_pmin,f_pmax,dpp
real(dp) :: fct,Utest,Rotest,Ptest,Ttest
real(dp) :: Pstar,Ustar,Rostar,R_Correction_Star,estar,Tstar,Cstar
real(dp) :: dcdro,dcdp,DerhdPstar,BB,dTdp_Ro,dTdRo_p,dRstardRo,dRstardT,Rstar
real(dp), dimension(Nb_VarC) :: dc,Cons_star,DerBBdUr,DerhDur,DerRcorrstardUr,DerRostardUr,DerUstardUr,DerhdUoth
real(dp), dimension(Nb_VarC) :: DerRcorrstardUr_Glob,DerRostardUr_Glob,DerUstardUr_Glob,DerPstardUr_Glob
real(dp), dimension(Nb_VarC) :: dRodUoth,dRdUoth
real(dp), dimension(Nb_VarC) :: DerRcorrstardUoth_Glob,DerRostardUoth_Glob,DerUstardUoth_Glob,DerPstardUoth_Glob
real(dp), dimension(Nb_VarC,Nb_VarC) :: Jacob_star,DerConstardUr,DerConstardUoth
real(dp) :: kk, KKK, dedT, dPdT, HH, dPde
real(dp), dimension(Nb_VarP) :: PPP
real(dp), dimension(Nb_VarP) :: PrimStar
integer :: i,nbpoints
Ror=me%br(2)%p%Prim(Pri_ro)
Pr=me%br(2)%p%Prim(Pri_P)
Ur=me%br(2)%p%Prim(Pri_u)
Cr=me%br(2)%p%Prim(Pri_c)
Zr=Ror*Cr
me%dynCC%mdot=mdot
me%dynCC%VarIn=VarIn
me%dynCC%area=me%br(2)%p%Area
me%dynCC%Pr=Pr
me%dynCC%Ur=Ur
me%dynCC%Zr=Zr
! Preparatory work for the solution of F to be found, knowing the fact that F is diminishing
fct=-1.0_dp
nbpoints=1000
pmin=Pr/1.5_dp
pmax=Pr*1.5_dp
dpp=(pmax-pmin)/(nbpoints-1)
do i=1,nbpoints
Ptest=Pmin+dpp*(i-1)
Utest=Ur+(Ptest-Pr)/Zr
Rotest=mdot/(Utest*me%br(2)%p%Area)
if(Rotest>1.0e-3_dp .and. Rotest<290.0_dp) then
Ttest=T_roP(Rotest,Ptest)
if(Ttest>2.5_dp .and. Ttest<500.0_dp) then
fct=me%dynCC%f(Ptest)
if(fct>0.0_dp) exit
endif
endif
enddo
if(fct<=0.0_dp .or. i==nbpoints) then
fct=-1.0_dp
nbpoints=1000000
pmin=Pr/1.5_dp
pmax=Pr*1.5_dp
dpp=(pmax-pmin)/(nbpoints-1)
do i=1,nbpoints
Ptest=Pmin+dpp*(i-1)
Utest=Ur+(Ptest-Pr)/Zr
Rotest=mdot/(Utest*me%br(2)%p%Area)
if(Rotest>1.0e-3_dp .and. Rotest<290.0_dp) then
Ttest=T_roP(Rotest,Ptest)
if(Ttest>2.5_dp .and. Ttest<500.0_dp) then
fct=me%dynCC%f(Ptest)
if(fct>0.0_dp) exit
endif
endif
enddo
if(fct<=0.0_dp .or. i==nbpoints) then
call set_error('initial bracket search failed in outgoing_branch_cold_circulator')
return
endif
endif
! Dichotomy
pmin=Ptest
f_pmin=me%dynCC%f(pmin)
pmax=Ptest*1.1_dp
f_pmax=me%dynCC%f(pmax)
if(f_pmin*f_pmax>0.0_dp) then
call set_error('initial range does not contain solution in outgoing_branch_cold_circulator')
return
else
!Brent's method
Pstar=zero(me%dynCC,pmin,pmax,1.0e-15_dp,1.0e-2_dp)
endif
Ustar=Ur+(Pstar-Pr)/Zr
Rostar=mdot/(Ustar*me%br(2)%p%Area)
call state_roP(Rostar, Pstar, R_Correction_Star, estar, Tstar, Cstar)
if(Tstar<2.0_dp .or. Tstar>50.0_dp) then
call set_error('unphysical temperature in outgoing_branch_cold_circulator')
return
endif
FluxOut%flx(Con_Mas)=Rostar*Ustar
FluxOut%flx(Con_Qdm)=Rostar*Ustar*Ustar+Pstar
FluxOut%flx(Con_Ene)=(Rostar*(estar+0.5_dp*Ustar*Ustar)+Pstar)*Ustar
FluxOut%flx(Con_R)=0.0_dp
if(R_Correction) FluxOut%flx(Con_R)=R_Correction_Star*Ustar
FluxOut%vit=Ustar
!------------------------------------------------------------------!
! Derivatives !
!------------------------------------------------------------------!
call dc2_roT(Ror, me%br(2)%p%Prim(Pri_T), dcdro, dcdp)
dcdro = dcdro/(2.0_dp*Cr)
dcdp = dcdp /(2.0_dp*Cr)
if(R_Correction) then
dc(Con_Mas)=dcdro-0.5_dp*dcdp*(me%br(2)%p%Cons(Con_Qdm)**2)/(me%br(2)%p%Cons(Con_Mas)**2)
dc(Con_Qdm)=dcdp*me%br(2)%p%Cons(Con_Qdm)/me%br(2)%p%Cons(Con_Mas)
dc(Con_Ene)=-dcdp
dc(Con_R)=dcdp
Cons_star(Con_Mas)=Rostar
Cons_star(Con_Qdm)=Rostar*Ustar
Cons_star(Con_Ene)=Rostar*(estar+0.5_dp*Ustar*Ustar)
Cons_star(Con_R)=R_Correction_Star
call jacobian_star(Cons_star,Jacob_star)
DerhdPstar=me%dynCC%df(Pstar)
BB=pstar-(me%br(2)%p%Cons(Con_R)-me%br(2)%p%Cons(Con_Ene)+0.5_dp*(me%br(2)%p%Cons(Con_Qdm)**2)/me%br(2)%p%Cons(Con_Mas))
DerBBdUr(1)=0.5_dp*(me%br(2)%p%Cons(Con_Qdm)**2)/(me%br(2)%p%Cons(Con_Mas)**2)
DerBBdUr(2)=-me%br(2)%p%Cons(Con_Qdm)/me%br(2)%p%Cons(Con_Mas)
DerBBdUr(3)=1.0_dp
DerBBdUr(4)=-1.0_dp
DerUstardUr(1)=-1.0_dp/(me%br(2)%p%Cons(Con_Mas)**2)*(me%br(2)%p%Cons(Con_Qdm)+BB/Cr)+&
+((DerBBdUr(1)*Cr-BB*dc(1))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr(2)=(1.0_dp+(DerBBdUr(2)*Cr-BB*dc(2))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr(3)=((DerBBdUr(3)*Cr-BB*dc(3))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr(4)=((DerBBdUr(4)*Cr-BB*dc(4))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerRostardUr(:)=(deriv_wr_both_br%DerMdotDerUout(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUr(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
DerRcorrstardUr(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUr(:)
DerhDur(:)=-((DerRcorrstardUr(:)*rostar-R_Correction_Star*DerRostardUr(:))/(rostar**2)+Ustar*DerUstardUr(:)-&
deriv_wr_both_br%DerVarIndUout(:))/1.0e3_dp
DerPstardUr_Glob(:)=-DerhDur(:)/DerhdPstar
BB=pstar-(me%br(2)%p%Cons(Con_R)-me%br(2)%p%Cons(Con_Ene)+0.5_dp*(me%br(2)%p%Cons(Con_Qdm)**2)/me%br(2)%p%Cons(Con_Mas))
DerBBdUr(1)=DerPstardUr_Glob(1)+0.5_dp*(me%br(2)%p%Cons(Con_Qdm)**2)/(me%br(2)%p%Cons(Con_Mas)**2)
DerBBdUr(2)=DerPstardUr_Glob(2)-me%br(2)%p%Cons(Con_Qdm)/me%br(2)%p%Cons(Con_Mas)
DerBBdUr(3)=DerPstardUr_Glob(3)+1.0_dp
DerBBdUr(4)=DerPstardUr_Glob(4)-1.0_dp
DerUstardUr_Glob(1)=-1.0_dp/(me%br(2)%p%Cons(Con_Mas)**2)*(me%br(2)%p%Cons(Con_Qdm)+BB/Cr)+&
+((DerBBdUr(1)*Cr-BB*dc(1))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr_Glob(2)=(1.0_dp+(DerBBdUr(2)*Cr-BB*dc(2))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr_Glob(3)=((DerBBdUr(3)*Cr-BB*dc(3))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerUstardUr_Glob(4)=((DerBBdUr(4)*Cr-BB*dc(4))/(Cr**2))/me%br(2)%p%Cons(Con_Mas)
DerRostardUr_Glob(:)=(deriv_wr_both_br%DerMdotDerUout(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUr_Glob(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
DerRcorrstardUr_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUr_Glob(:)+dRstardT*dTdp_Ro*DerPstardUr_Glob(:)
DerConstardUr(:,1)=DerRostardUr_Glob(:)
DerConstardUr(:,2)=Rostar*DerUstardUr_Glob(:)+DerRostardUr_Glob(:)*Ustar
DerConstardUr(:,3)=DerRcorrstardUr_Glob(:)+0.5_dp*(DerRostardUr_Glob(:)*Ustar*Ustar+&
2.0_dp*Rostar*Ustar*DerUstardUr_Glob(:))-DerPstardUr_Glob(:)
DerConstardUr(:,4)=DerRcorrstardUr_Glob(:)
FluxOut%DerFlxdBr_b(:,:)=matmul(DerConstardUr(:,:),Jacob_star(:,:))
FluxOut%DerVitdBr_b(:)=DerUstardUr_Glob(:)
! Contribution of the cold circulator inlet
dRodUoth(:)=deriv_wr_both_br%DerMdotDerUin(:)/(Ustar*me%br(2)%p%Area)
dRdUoth(:)=(dRstardRo+dRstardT*dTdRo_p)*dRodUoth(:)
DerhdUoth(:)=-((dRdUoth(:)*Rostar-R_Correction_Star*dRodUoth(:))/(Rostar**2)-deriv_wr_both_br%DerVarIndUin(:))/1.0e3_dp
DerPstardUoth_Glob(:)=-DerhDuOth(:)/DerhdPstar
DerUstardUoth_Glob(:)=DerPstardUoth_Glob(:)/Zr
DerRostardUoth_Glob(:)=(deriv_wr_both_br%DerMdotDerUin(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUoth_Glob(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
DerRcorrstardUoth_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUoth_Glob(:)+dRstardT*dTdp_Ro*DerPstardUoth_Glob(:)
DerConstardUoth(:,1)=DerRostardUoth_Glob(:)
DerConstardUoth(:,2)=Rostar*DerUstardUoth_Glob(:)+DerRostardUoth_Glob(:)*Ustar
DerConstardUoth(:,3)=DerRcorrstardUoth_Glob(:)+0.5_dp*(DerRostardUoth_Glob(:)*Ustar*Ustar+&
2.0_dp*Rostar*Ustar*DerUstardUoth_Glob(:))-DerPstardUoth_Glob(:)
DerConstardUoth(:,4)=DerRcorrstardUoth_Glob(:)
FluxOut%DerFlxdBr_a(:,:)=matmul(DerConstardUoth(:,:),Jacob_star(:,:))
FluxOut%DerVitdBr_a(:)=DerUstardUoth_Glob(:)
else
PrimStar(Pri_ro)=Rostar
PrimStar(Pri_u)=Ustar
PrimStar(Pri_p)=pstar
PrimStar(Pri_e)=estar
PrimStar(Pri_T)=Tstar
PrimStar(Pri_c)=Cstar
call jacobian_star_Req_rem(PrimStar,Jacob_star)
DerhdPstar=me%dynCC%df(Pstar)
PPP(:)=me%br(2)%p%Prim(:)
HH=PPP(Pri_e)+0.5_dp*PPP(Pri_u)**2+PPP(Pri_p)/PPP(Pri_ro)
call jacobian_roT(PPP(Pri_ro), PPP(Pri_T), dedT, dPdT)
dPde=dPdT/dedT
kk=dPde/PPP(Pri_ro)
KKK=PPP(Pri_c)**2+kk*(PPP(Pri_u)**2-HH)
BB=Pstar-Pr
DerBBdUr(1)=-KKK
DerBBdUr(2)=kk*PPP(Pri_u)
DerBBdUr(3)=-kk
DerBBdUr(4)=0.0_dp
dc(Con_Mas)=dcdro+dcdp*KKK
dc(Con_Qdm)=-dcdp*kk*PPP(Pri_u)
dc(Con_Ene)=dcdp*kk
dc(Con_R)=0.0_dp
DerUstardUr(1)=-PPP(Pri_u)/PPP(Pri_ro)+(DerBBdUr(1)*PPP(Pri_ro)*Cr-BB*(Cr+PPP(Pri_ro)*dc(1)))/((PPP(Pri_ro)*Cr)**2)
DerUstardUr(2)=(1.0_dp/PPP(Pri_ro))*(1.0_dp+(DerBBdUr(2)*Cr-BB*dc(2))/(Cr**2))
DerUstardUr(3)=(1.0_dp/PPP(Pri_ro))*(DerBBdUr(3)*Cr-BB*dc(3))/(Cr**2)
DerUstardUr(4)=0.0_dp
DerRostardUr(:)=(deriv_wr_both_br%DerMdotDerUout(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUr(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
DerRcorrstardUr(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUr(:)
DerhDur(:)=-((DerRcorrstardUr(:)*rostar-R_Correction_Star*DerRostardUr(:))/(rostar**2)+Ustar*DerUstardUr(:)-&
deriv_wr_both_br%DerVarIndUout(:))/1.0e3_dp
DerPstardUr_Glob(:)=-DerhDur(:)/DerhdPstar
BB=Pstar-Pr
DerBBdUr(1)=DerPstardUr_Glob(1)-KKK
DerBBdUr(2)=DerPstardUr_Glob(2)+kk*PPP(Pri_u)
DerBBdUr(3)=DerPstardUr_Glob(3)-kk
DerBBdUr(4)=0.0_dp
DerUstardUr_Glob(1)=-PPP(Pri_u)/PPP(Pri_ro)+(DerBBdUr(1)*PPP(Pri_ro)*Cr-BB*(Cr+PPP(Pri_ro)*dc(1)))/((PPP(Pri_ro)*Cr)**2)
DerUstardUr_Glob(2)=(1.0_dp/PPP(Pri_ro))*(1.0_dp+(DerBBdUr(2)*Cr-BB*dc(2))/(Cr**2))
DerUstardUr_Glob(3)=(1.0_dp/PPP(Pri_ro))*(DerBBdUr(3)*Cr-BB*dc(3))/(Cr**2)
DerUstardUr_Glob(4)=0.0_dp
DerRostardUr_Glob(:)=(deriv_wr_both_br%DerMdotDerUout(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUr_Glob(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
DerRcorrstardUr_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUr_Glob(:)+dRstardT*dTdp_Ro*DerPstardUr_Glob(:)
DerConstardUr(:,1)=DerRostardUr_Glob(:)
DerConstardUr(:,2)=Rostar*DerUstardUr_Glob(:)+DerRostardUr_Glob(:)*Ustar
DerConstardUr(:,3)=DerRcorrstardUr_Glob(:)+0.5_dp*(DerRostardUr_Glob(:)*Ustar*Ustar+&
2.0_dp*Rostar*Ustar*DerUstardUr_Glob(:))-DerPstardUr_Glob(:)
DerConstardUr(:,4)=0.0_dp
FluxOut%DerFlxdBr_b(:,:)=matmul(DerConstardUr(:,:),Jacob_star(:,:))
FluxOut%DerFlxdBr_b(:,4)=0.0_dp
FluxOut%DerVitdBr_b(:)=DerUstardUr_Glob(:)
! Contribution of the cold circulator inlet
dRodUoth(:)=deriv_wr_both_br%DerMdotDerUin(:)/(Ustar*me%br(2)%p%Area)
dRdUoth(:)=(dRstardRo+dRstardT*dTdRo_p)*dRodUoth(:)
DerhdUoth(:)=-((dRdUoth(:)*Rostar-R_Correction_Star*dRodUoth(:))/(Rostar**2)-deriv_wr_both_br%DerVarIndUin(:))/1.0e3_dp
DerPstardUoth_Glob(:)=-DerhDuOth(:)/DerhdPstar
DerUstardUoth_Glob(:)=DerPstardUoth_Glob(:)/Zr
DerRostardUoth_Glob(:)=(deriv_wr_both_br%DerMdotDerUin(:)*Ustar*me%br(2)%p%Area-mdot*DerUstardUoth_Glob(:)*me%br(2)%p%Area)/&
((Ustar*me%br(2)%p%Area)**2)
DerRcorrstardUoth_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUoth_Glob(:)+dRstardT*dTdp_Ro*DerPstardUoth_Glob(:)
DerConstardUoth(:,1)=DerRostardUoth_Glob(:)
DerConstardUoth(:,2)=Rostar*DerUstardUoth_Glob(:)+DerRostardUoth_Glob(:)*Ustar
DerConstardUoth(:,3)=DerRcorrstardUoth_Glob(:)+0.5_dp*(DerRostardUoth_Glob(:)*Ustar*Ustar+&
2.0_dp*Rostar*Ustar*DerUstardUoth_Glob(:))-DerPstardUoth_Glob(:)
DerConstardUoth(:,4)=0.0_dp
FluxOut%DerFlxdBr_a(:,:)=matmul(DerConstardUoth(:,:),Jacob_star(:,:))
FluxOut%DerFlxdBr_a(:,4)=0.0_dp
FluxOut%DerVitdBr_a(:)=DerUstardUoth_Glob(:)
endif
ssOut=abs(Ur)+Cr
end subroutine outgoing_branch_cold_circulator