solve_boundary_MT Subroutine

public subroutine solve_boundary_MT(me, ss)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Arguments

Type IntentOptional Attributes Name
type(boundary_t), intent(inout) :: me
real(kind=dp), intent(out) :: ss

Calls

proc~~solve_boundary_mt~~CallsGraph proc~solve_boundary_mt solve_boundary_MT proc~dc2_rot dc2_roT proc~solve_boundary_mt->proc~dc2_rot proc~function_of_pstar_mt function_of_pstar_MT proc~solve_boundary_mt->proc~function_of_pstar_mt proc~jacobian_rot jacobian_roT proc~solve_boundary_mt->proc~jacobian_rot proc~jacobian_star_req_rem~2 jacobian_star_Req_rem proc~solve_boundary_mt->proc~jacobian_star_req_rem~2 proc~jacobian_star~2 jacobian_star proc~solve_boundary_mt->proc~jacobian_star~2 proc~set_error set_error proc~solve_boundary_mt->proc~set_error proc~signal_v0d signal_t%signal_v0d proc~solve_boundary_mt->proc~signal_v0d proc~state_rop state_roP proc~solve_boundary_mt->proc~state_rop proc~t_rop T_roP proc~solve_boundary_mt->proc~t_rop proc~eos_terms eos_terms proc~dc2_rot->proc~eos_terms proc~function_of_pstar_mt->proc~t_rop proc~jacobian_rot->proc~eos_terms proc~jacobian_star_req_rem~2->proc~jacobian_rot proc~signal_v1d signal_t%signal_v1d proc~signal_v0d->proc~signal_v1d proc~state_rop->proc~jacobian_rot proc~state_rop->proc~t_rop proc~brent brent proc~t_rop->proc~brent proc~fill_f_terms fill_f_terms proc~t_rop->proc~fill_f_terms proc~fill_g_dreg_terms fill_g_Dreg_terms proc~t_rop->proc~fill_g_dreg_terms proc~brent->proc~set_error f f proc~brent->f proc~zero zero proc~brent->proc~zero proc~eos_terms->proc~fill_f_terms proc~eos_e_terms eos_e_terms proc~eos_terms->proc~eos_e_terms proc~fill_e_tpart fill_e_Tpart proc~eos_e_terms->proc~fill_e_tpart proc~fill_ff_terms fill_ff_terms proc~eos_e_terms->proc~fill_ff_terms proc~zero->f proc~fill_e_tpart->proc~fill_g_dreg_terms

Called by

proc~~solve_boundary_mt~~CalledByGraph proc~solve_boundary_mt solve_boundary_MT proc~boundary_resolution_from_and_to_ports boundary_resolution_from_and_to_ports proc~boundary_resolution_from_and_to_ports->proc~solve_boundary_mt proc~main_loop main_loop proc~main_loop->proc~boundary_resolution_from_and_to_ports program~reims_p reims_p program~reims_p->proc~main_loop

Source Code

subroutine solve_boundary_MT(me,SS)     
  type(boundary_t), intent(inout) :: me
  real(dp), intent(out) :: ss
 
  real(dp) :: Ttank, ro_LR, p_LR, u_LR, c_LR, z_LR, R_Correction_Star,R_Phys_Star
  real(dp) :: Pstar, Tstar, Ustar, Rostar, estar, Cstar
  real(dp) :: DerhdPstar,dcdro,dcdp,BB,dTdp_Ro,dTdRo_p,dRstardRo,dRstardT,Rstar,T_temp,T_temp_plus
  real(dp), dimension(Nb_VarC) :: DerBBdU_LR,DerUstardU_LR_Glob,Cs_cell,Cons_star
  real(dp), dimension(Nb_VarC) :: DerPstardU_LR_Glob,DerRostardU_LR_Glob,DerRcorrstardU_LR_Glob,DerRphysstardU_LR_Glob
  real(dp), dimension(Nb_VarC) :: DerUstardU_LR,DerhdU_LR,dc,DerRostardU_LR
  real(dp), dimension(Nb_VarC,Nb_VarC) :: Jacob_star,DerConstardU_LR
  real(dp) :: kk, KKK, dedT, dPdT, HH, dPde
  real(dp), dimension(Nb_VarP) :: PrimStar
  real(dp), dimension(Nb_VarP) :: Pr_cell
  integer :: NbreMaxIte, ite, i
  real(dp) :: pmin,pmid,pmax,f_pmin,f_pmid,f_pmax,AreaPP,mdot,sgnBC,Nbcel,dpp,rostarLoc,FctLoc

  Pr_cell(:)=me%port%Prim(:)
  Cs_cell(:)=me%port%Cons(:)
  AreaPP=me%port%Area
  sgnBC=-me%port%Sgn4j  ! sgnBC=1 if 1-port, -1 if 2-port

  Ttank=me%imposed_T%v0d()
  mdot=me%imposed_PorM%v0d()

  ro_LR=Pr_cell(Pri_ro)
  p_LR=Pr_cell(Pri_P)
  u_LR=Pr_cell(Pri_u)
  c_LR=Pr_cell(Pri_c)
  z_LR=ro_LR*c_LR

  ! Analysis of increasing function
  pmin=p_LR/2.0_dp
  pmax=p_LR*2.0_dp
  Nbcel=100000.0_dp
  dpp=(pmax-pmin)/Nbcel
  do i=1,Nbcel
    rostarLoc=mdot/((u_LR+sgnBC*(pmin+dpp*(i-1)-p_LR)/z_LR)*AreaPP)
    if(rostarLoc<=0.0_dp .or. rostarLoc>200.0_dp) cycle
    FctLoc=function_of_pstar_MT(sgnBC,pmin+dpp*(i-1),Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)/1.0e3_dp
    pmin=pmin+dpp*(i-1)
    exit
  enddo

  ! Dichotomy
  f_pmin=function_of_pstar_MT(sgnBC,pmin,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)/1.0e3_dp

  f_pmax=function_of_pstar_MT(sgnBC,pmax,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)/1.0e3_dp  

  if(f_pmin*f_pmax>0.0_dp) then
    call set_error('initial range does not contain solution in solve_boundary_MT')
    return
  else
    NbreMaxIte=200
    do ite=1,NbreMaxIte
      pmid=0.5_dp*(pmin+pmax)
      f_pmid=function_of_pstar_MT(sgnBC,pmid,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)/1.0e3_dp
      if(f_pmin*f_pmid<0.0_dp) then
        pmax=pmid
      else
        pmin=pmid
        f_pmin=f_pmid
      endif
      if(abs(pmax-pmin)<1.0e-2_dp) exit
      if(ite==NbreMaxIte) then
        call set_error('convergence failure in solve_boundary_MT')
        return
      endif
    enddo  
  endif
  Pstar=0.5_dp*(pmin+pmax)

  Ustar=u_LR+sgnBC*(Pstar-p_LR)/z_LR
  Rostar=mdot/(Ustar*AreaPP)    ! Ustar and mdot always share the same sign
  call state_roP(Rostar, Pstar, R_Phys_Star, estar, Tstar, Cstar)
  R_Correction_Star=1.0_dp
  if(R_Correction) R_Correction_Star=R_Phys_Star

  me%port%flx(Con_Mas)=Rostar*Ustar
  me%port%flx(Con_Qdm)=Rostar*Ustar*Ustar+Pstar
  me%port%flx(Con_Ene)=(Rostar*(estar+0.5_dp*Ustar*Ustar)+Pstar)*Ustar
  me%port%flx(Con_R)=0.0_dp
  if(R_Correction) me%port%flx(Con_R)=R_Correction_Star*Ustar

  me%port%vit=Ustar

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                            Derivatives                           !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  call dc2_roT(ro_LR, Pr_cell(Pri_T), dcdro, dcdp)
  dcdro = dcdro/(2.0_dp*c_LR)
  dcdp  = dcdp /(2.0_dp*c_LR)

  dc(Con_Mas)=dcdro-0.5_dp*dcdp*(Cs_cell(Con_Qdm)**2)/(Cs_cell(Con_Mas)**2)
  dc(Con_Qdm)=dcdp*Cs_cell(Con_Qdm)/Cs_cell(Con_Mas)
  dc(Con_Ene)=-dcdp
  dc(Con_R)=dcdp   

  if(R_Correction) then

    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=(function_of_pstar_MT(sgnBC,Pstar+1.0e-3_dp,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)-&
                function_of_pstar_MT(sgnBC,Pstar,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR))    

    BB=pstar-(Cs_cell(Con_R)-Cs_cell(Con_Ene)+0.5_dp*(Cs_cell(Con_Qdm)**2)/Cs_cell(Con_Mas))
    DerBBdU_LR(1)=0.5_dp*(Cs_cell(Con_Qdm)**2)/(Cs_cell(Con_Mas)**2)
    DerBBdU_LR(2)=-Cs_cell(Con_Qdm)/Cs_cell(Con_Mas)
    DerBBdU_LR(3)=1.0_dp
    DerBBdU_LR(4)=-1.0_dp            
    
    DerUstardU_LR(1)=-1.0_dp/(Cs_cell(Con_Mas)**2)*(Cs_cell(Con_Qdm)+sgnBC*BB/c_LR)+&
                   +(sgnBC*(DerBBdU_LR(1)*c_LR-BB*dc(1))/(c_LR**2))/Cs_cell(Con_Mas)
                   DerUstardU_LR(2)=(1.0_dp+sgnBC*(DerBBdU_LR(2)*c_LR-BB*dc(2))/(c_LR**2))/Cs_cell(Con_Mas)
    DerUstardU_LR(3)=(sgnBC*(DerBBdU_LR(3)*c_LR-BB*dc(3))/(c_LR**2))/Cs_cell(Con_Mas)
    DerUstardU_LR(4)=(sgnBC*(DerBBdU_LR(4)*c_LR-BB*dc(4))/(c_LR**2))/Cs_cell(Con_Mas)

    DerRostardU_LR(:)=-(mdot/AreaPP)*DerUstardU_LR(:)/(Ustar**2)

    T_temp=T_roP(Rostar,Pstar)
    T_temp_plus=T_roP(Rostar+1.0e-3_dp,Pstar)
    DerhdU_LR(:)=DerRostardU_LR(:)*(T_temp_plus-T_temp)
    DerPstardU_LR_Glob(:)=-DerhdU_LR(:)/DerhdPstar    

    BB=pstar-(Cs_cell(Con_R)-Cs_cell(Con_Ene)+0.5_dp*(Cs_cell(Con_Qdm)**2)/Cs_cell(Con_Mas))
    DerBBdU_LR(1)=DerPstardU_LR_Glob(1)+0.5_dp*(Cs_cell(Con_Qdm)**2)/(Cs_cell(Con_Mas)**2)
    DerBBdU_LR(2)=DerPstardU_LR_Glob(2)-Cs_cell(Con_Qdm)/Cs_cell(Con_Mas)
    DerBBdU_LR(3)=DerPstardU_LR_Glob(3)+1.0_dp
    DerBBdU_LR(4)=DerPstardU_LR_Glob(4)-1.0_dp    

    DerUstardU_LR_Glob(1)=-1.0_dp/(Cs_cell(Con_Mas)**2)*(Cs_cell(Con_Qdm)+sgnBC*BB/c_LR)+&
                        +(sgnBC*(DerBBdU_LR(1)*c_LR-BB*dc(1))/(c_LR**2))/Cs_cell(Con_Mas)
    DerUstardU_LR_Glob(2)=(1.0_dp+sgnBC*(DerBBdU_LR(2)*c_LR-BB*dc(2))/(c_LR**2))/Cs_cell(Con_Mas)
    DerUstardU_LR_Glob(3)=(sgnBC*(DerBBdU_LR(3)*c_LR-BB*dc(3))/(c_LR**2))/Cs_cell(Con_Mas)
    DerUstardU_LR_Glob(4)=(sgnBC*(DerBBdU_LR(4)*c_LR-BB*dc(4))/(c_LR**2))/Cs_cell(Con_Mas)    

    DerRostardU_LR_Glob(:)=-(mdot/AreaPP)*DerUstardU_LR_Glob(:)/(Ustar**2)

    call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
    DerRcorrstardU_LR_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardU_LR_Glob(:)+dRstardT*dTdp_Ro*DerPstardU_LR_Glob(:) 

    DerConstardU_LR(:,1)=DerRostardU_LR_Glob(:)
    DerConstardU_LR(:,2)=Rostar*DerUstardU_LR_Glob(:)+DerRostardU_LR_Glob(:)*Ustar
    DerConstardU_LR(:,3)=DerRcorrstardU_LR_Glob(:)+0.5_dp*(DerRostardU_LR_Glob(:)*Ustar*Ustar+&
                         2.0_dp*Rostar*Ustar*DerUstardU_LR_Glob(:))-DerPstardU_LR_Glob(:)
    DerConstardU_LR(:,4)=DerRcorrstardU_LR_Glob(:)

    me%port%derFlx_derCon(:,:)=matmul(DerConstardU_LR(:,:),Jacob_star(:,:))
    
    me%port%derVit_derCon(:)=DerUstardU_LR_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
    PrimStar(Pri_R)=R_Correction_Star
    call jacobian_star_Req_rem(PrimStar,Jacob_star)

    DerhdPstar=(function_of_pstar_MT(sgnBC,Pstar+1.0e-3_dp,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR)-&
                function_of_pstar_MT(sgnBC,Pstar,Ttank,mdot,AreaPP,p_LR,u_LR,z_LR))

    HH=Pr_cell(Pri_e)+0.5_dp*Pr_cell(Pri_u)**2+Pr_cell(Pri_p)/Pr_cell(Pri_ro)
    call jacobian_roT(Pr_cell(Pri_ro), Pr_cell(Pri_T), dedT, dPdT)
    dPde=dPdT/dedT
    kk=dPde/Pr_cell(Pri_ro)
    KKK=Pr_cell(Pri_c)**2+kk*(Pr_cell(Pri_u)**2-HH)

    BB=Pstar-p_LR
    DerBBdU_LR(1)=-KKK
    DerBBdU_LR(2)=kk*Pr_cell(Pri_u)
    DerBBdU_LR(3)=-kk
    DerBBdU_LR(4)=0.0_dp

    dc(Con_Mas)=dcdro+dcdp*KKK
    dc(Con_Qdm)=-dcdp*kk*Pr_cell(Pri_u)
    dc(Con_Ene)=dcdp*kk
    dc(Con_R)=0.0_dp

    DerUstardU_LR(1)=-Pr_cell(Pri_u)/Pr_cell(Pri_ro)+sgnBC*(DerBBdU_LR(1)*Pr_cell(Pri_ro)*c_LR-BB*(c_LR+Pr_cell(Pri_ro)*dc(1)))/&
                    ((Pr_cell(Pri_ro)*c_LR)**2)
    DerUstardU_LR(2)=(1.0_dp/Pr_cell(Pri_ro))*(1.0_dp+sgnBC*(DerBBdU_LR(2)*c_LR-BB*dc(2))/(c_LR**2))
    DerUstardU_LR(3)=sgnBC*(1.0_dp/Pr_cell(Pri_ro))*(DerBBdU_LR(3)*c_LR-BB*dc(3))/(c_LR**2)
    DerUstardU_LR(4)=0.0_dp

    DerRostardU_LR(:)=-(mdot/AreaPP)*DerUstardU_LR(:)/(Ustar**2)
     
    T_temp=T_roP(Rostar,Pstar)
    T_temp_plus=T_roP(Rostar+1.0e-3_dp,Pstar)
    DerhdU_LR(:)=DerRostardU_LR(:)*(T_temp_plus-T_temp)
    DerPstardU_LR_Glob(:)=-DerhdU_LR(:)/DerhdPstar

    BB=Pstar-p_LR
    DerBBdU_LR(1)=DerPstardU_LR_Glob(1)-KKK
    DerBBdU_LR(2)=DerPstardU_LR_Glob(2)+kk*Pr_cell(Pri_u)
    DerBBdU_LR(3)=DerPstardU_LR_Glob(3)-kk
    DerBBdU_LR(4)=0.0_dp

    DerUstardU_LR_Glob(1)=-Pr_cell(Pri_u)/Pr_cell(Pri_ro)+sgnBC*(DerBBdU_LR(1)*Pr_cell(Pri_ro)*c_LR-BB*(c_LR+Pr_cell(Pri_ro)*&
                          dc(1)))/((Pr_cell(Pri_ro)*c_LR)**2)
    DerUstardU_LR_Glob(2)=(1.0_dp/Pr_cell(Pri_ro))*(1.0_dp+sgnBC*(DerBBdU_LR(2)*c_LR-BB*dc(2))/(c_LR**2))
    DerUstardU_LR_Glob(3)=sgnBC*(1.0_dp/Pr_cell(Pri_ro))*(DerBBdU_LR(3)*c_LR-BB*dc(3))/(c_LR**2)
    DerUstardU_LR_Glob(4)=0.0_dp

    DerRostardU_LR_Glob(:)=-(mdot/AreaPP)*DerUstardU_LR_Glob(:)/(Ustar**2)

    call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
    DerRphysstardU_LR_Glob(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardU_LR_Glob(:)+dRstardT*dTdp_Ro*DerPstardU_LR_Glob(:)

    DerConstardU_LR(:,1)=DerRostardU_LR_Glob(:)
    DerConstardU_LR(:,2)=Rostar*DerUstardU_LR_Glob(:)+DerRostardU_LR_Glob(:)*Ustar
    DerConstardU_LR(:,3)=DerRphysstardU_LR_Glob(:)+0.5_dp*(DerRostardU_LR_Glob(:)*Ustar*Ustar+&
                       2.0_dp*Rostar*Ustar*DerUstardU_LR_Glob(:))-DerPstardU_LR_Glob(:)
    DerConstardU_LR(:,4)=0.0_dp     

    me%port%derFlx_derCon(:,:)=matmul(DerConstardU_LR(:,:),Jacob_star(:,:))
    me%port%derFlx_derCon(:,4)=0.0_dp

    me%port%derVit_derCon(:)=DerUstardU_LR_Glob(:)

  endif

  ss=abs(u_LR)+c_LR
end subroutine solve_boundary_MT