solve_boundary_PT

public subroutine solve_boundary_PT(me, ss)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Arguments

Type	Intent	Optional	Attributes		Name
type(boundary_t),	intent(inout)			::	me
real(kind=dp),	intent(out)			::	ss
Called by

Help
Source Code

subroutine solve_boundary_PT(me,ss)
  type(boundary_t), intent(inout) :: me
  real(dp), intent(out) :: ss
 
  real(dp) :: kk, KKK, dedT, dPdT, HH, dPde, sgnBC
  real(dp), dimension(Nb_VarP) :: PrimStar
  real(dp) :: Ptank, Ttank, Rotank, Ctank, Ro_LR, p_LR, u_LR, Einttank, c_LR, z_LR, R_Correction_Star
  real(dp) :: Pstar, Tstar, Ustar, Rostar, r_Corr_LR, estar, Cstar, Funct, DerFunct
  real(dp) :: DerhdPstar,dcdro,dcdp,BB,dTdp_Ro,dTdRo_p,dRstardRo,dRstardT,Rstar,dRdP,T_temp
  real(dp) :: Rout,Pout,eout,cout,dR_dRo,R_Phys_Star
  real(dp), dimension(Nb_VarP) :: Pr_cell
  real(dp), dimension(Nb_VarC) :: Cs_cell,Cons_star,DerBBdU_LR,DerUstardU_LR_Glob
  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,DerRcorrstardU_LR,DerRphystardU_LR
  real(dp), dimension(Nb_VarC,Nb_VarC) :: Jacob_star,DerConstardU_LR
  integer :: NbreMaxIte, ite
  logical :: Inflow, use_hugoniot
  type(hug_star_t) :: st
  real(dp), dimension(Nb_VarC) :: dvstar_dU, dPsi_dU
  real(dp) :: v_LR, DeltaP_hug, Deltav_hug

  st%ok = .false.
  Ptank = me%imposed_PorM%v0d()
  Ttank = me%imposed_T%v0d()
  Rotank = ro_pT(Ptank,Ttank)
  call state_roT(Rotank, Ttank, Rout, Pout, eout, cout)
  Ctank = cout
  Einttank = eout

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

  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
  r_Corr_LR=Pr_cell(Pri_R)
  use_hugoniot = me%hugoniot_bc .and. &
                 abs(Ptank - p_LR) >= 1.0e-7_dp * Ro_LR * c_LR**2

  ! Inflow/outflow criterion: Hugoniot-consistent when use_hugoniot, acoustic otherwise
  if (use_hugoniot) then
    call hugoniot_star_state(sgnBC, Ptank, Ro_LR, p_LR, u_LR, c_LR, st)
    if (.not. st%ok) then; call set_error('hugoniot criterion: solver failed'); return; end if
    Inflow = sgnBC * st%ustar > 0.0_dp
  else
    Inflow = Ptank > p_LR - sgnBC*z_LR*u_LR
  end if

  if (Inflow) then
    Pstar=0.5_dp*(Ptank+p_LR)
    NbreMaxIte=200

    ! Newton method for Pstar calculation
    do ite=1,NbreMaxIte

      if (use_hugoniot) then
        call hugoniot_star_state(sgnBC, Pstar, Ro_LR, p_LR, u_LR, c_LR, st)
        if (.not. st%ok) then; call set_error('hugoniot inflow: solver failed'); return; end if
        Rostar=Rotank+(Pstar-Ptank)/(Ctank**2)
        call state_roP(Rostar, Pstar, R_Phys_Star, estar, Tstar, Cstar)
        if (ieee_is_nan(Tstar)) then; call set_error('convergence failure in trop (hugoniot inflow)'); return; end if
        Funct=estar+Pstar/Rostar+0.5_dp*(st%ustar**2)-Einttank-Ptank/Rotank
        call jacobian_roT(Rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
        dRdP=(dRstardRo+dRstardT*dTdRo_p)/(Ctank**2)+dRstardT*dTdp_Ro
        DerFunct=(dRdP*Rostar-Rstar/(Ctank**2))/(Rostar**2)+st%ustar*st%dustar_dPstar
      else
        Funct=function_of_pstar_PT(sgnBC,Pstar,Rotank,Ctank,Ptank,Einttank,p_LR,u_LR,z_LR)
        DerFunct=derivative_function_of_pstar_PT(sgnBC,Pstar,Rotank,Ctank,Ptank,p_LR,u_LR,z_LR)
      end if

      if(abs(Funct)>1.0e-8_dp) then
        Pstar=Pstar-Funct/DerFunct
        Pstar=max(Pstar, 1.0e-6_dp*p_LR)
      else
        if (use_hugoniot) then
          Ustar=st%ustar
        else
          Ustar=u_LR+sgnBC*(Pstar-p_LR)/z_LR
        end if
        Rostar=Rotank+(Pstar-Ptank)/(Ctank**2)
        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
        exit
      endif

      if(ite==NbreMaxIte) then
        call set_error('convergence failure in solve_boundary_condition_pipe_left')
        return
      endif

    enddo
    ! End of Newton method

  else

  ! Outflow
    Pstar=Ptank
    if(R_Correction) then
      if (use_hugoniot) then
        Rostar = 1.0_dp / st%vstar
        Ustar  = st%ustar
      else
        Ustar=u_LR+sgnBC*(Pstar-p_LR)/z_LR
        Rostar=(Pstar-p_LR)/(c_LR**2)+Ro_LR
      end if
      R_Correction_Star=r_Corr_LR*Rostar/Ro_LR+(c_LR**2)*(Rostar-Ro_LR)
      call state_roP_withR(Rostar, Pstar, R_Correction_Star, estar, Tstar, Cstar)
    else
      if (use_hugoniot) then
        Rostar = 1.0_dp / st%vstar
        Ustar  = st%ustar
      else
        Ustar=u_LR+sgnBC*(Pstar-p_LR)/z_LR
        Rostar=(Pstar-p_LR)/(c_LR**2)+Ro_LR
      end if
      call state_roP(Rostar, Pstar, R_Phys_Star, estar, Tstar, Cstar)
      if (ieee_is_nan(Tstar)) then; call set_error('convergence failure in T_roP (boundary)'); return; end if
      R_Correction_Star=r_Corr_LR
    endif
  endif
   
  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

  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(Inflow) then

    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)

      if (use_hugoniot) then
        DerhdPstar = DerFunct
      else
        DerhdPstar=derivative_function_of_pstar_PT(sgnBC,Pstar,Rotank,Ctank,Ptank,p_LR,u_LR,z_LR)
      end if

      if (use_hugoniot) then
        v_LR       = 1.0_dp / Ro_LR
        DeltaP_hug = Pstar - p_LR
        Deltav_hug = v_LR - st%vstar

        dvstar_dU(1) = st%dvstar_dv_LR*(-v_LR**2) + st%dvstar_dp_LR*(-u_LR**2/2.0_dp)
        dvstar_dU(2) = st%dvstar_dp_LR*u_LR
        dvstar_dU(3) = st%dvstar_dp_LR*(-1.0_dp)
        dvstar_dU(4) = st%dvstar_dp_LR*1.0_dp

        dPsi_dU(1) = ((-v_LR**2 - dvstar_dU(1))*DeltaP_hug + Deltav_hug*(-u_LR**2/2.0_dp)) / DeltaP_hug**2
        dPsi_dU(2) = (            -dvstar_dU(2) *DeltaP_hug + Deltav_hug*u_LR)              / DeltaP_hug**2
        dPsi_dU(3) = (            -dvstar_dU(3) *DeltaP_hug + Deltav_hug*(-1.0_dp))         / DeltaP_hug**2
        dPsi_dU(4) = (            -dvstar_dU(4) *DeltaP_hug + Deltav_hug*1.0_dp)            / DeltaP_hug**2

        DerUstardU_LR(1) = -u_LR*v_LR + sgnBC*st%K*(u_LR**2/2.0_dp) + sgnBC*DeltaP_hug*dPsi_dU(1)/(2.0_dp*st%K)
        DerUstardU_LR(2) =  v_LR      - sgnBC*st%K*u_LR              + sgnBC*DeltaP_hug*dPsi_dU(2)/(2.0_dp*st%K)
        DerUstardU_LR(3) =              sgnBC*st%K                    + sgnBC*DeltaP_hug*dPsi_dU(3)/(2.0_dp*st%K)
        DerUstardU_LR(4) =            - sgnBC*st%K                    + sgnBC*DeltaP_hug*dPsi_dU(4)/(2.0_dp*st%K)

        DerhdU_LR(:) = Ustar * DerUstardU_LR(:)
        DerPstardU_LR_Glob(:) = -DerhdU_LR(:)/DerhdPstar
        DerRostardU_LR_Glob(:) = DerPstardU_LR_Glob(:)/(Ctank**2)
        DerUstardU_LR_Glob(:) = DerUstardU_LR(:) + st%dustar_dPstar * DerPstardU_LR_Glob(:)

      else

        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR)
        DerhdU_LR(:)=Ustar*DerUstardU_LR(:)

        DerPstardU_LR_Glob(:)=-DerhdU_LR(:)/DerhdPstar

        DerRostardU_LR_Glob(:)=DerPstardU_LR_Glob(:)/(Ctank**2)

        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR_Glob)

      end if

      call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
      dRdP=(dRstardRo+dRstardT*dTdRo_p)/(Ctank**2)+dRstardT*dTdp_Ro
      DerRcorrstardU_LR_Glob(:)=dRdP*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)

      if (use_hugoniot) then
        DerhdPstar=DerFunct
      else
        DerhdPstar=derivative_function_of_pstar_PT(sgnBC,Pstar,Rotank,Ctank,Ptank,p_LR,u_LR,z_LR)
      end if

      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)

      if (use_hugoniot) then
        v_LR       = 1.0_dp / Ro_LR
        DeltaP_hug = Pstar - p_LR
        Deltav_hug = v_LR - st%vstar

        dvstar_dU(1) = st%dvstar_dv_LR*(-v_LR**2) + st%dvstar_dp_LR*KKK
        dvstar_dU(2) = st%dvstar_dp_LR*(-kk*u_LR)
        dvstar_dU(3) = st%dvstar_dp_LR*kk
        dvstar_dU(4) = 0.0_dp

        dPsi_dU(1) = ((-v_LR**2 - dvstar_dU(1))*DeltaP_hug + Deltav_hug*KKK)       / DeltaP_hug**2
        dPsi_dU(2) = (           -dvstar_dU(2) *DeltaP_hug + Deltav_hug*(-kk*u_LR)) / DeltaP_hug**2
        dPsi_dU(3) = (           -dvstar_dU(3) *DeltaP_hug + Deltav_hug*kk)         / DeltaP_hug**2
        dPsi_dU(4) = 0.0_dp

        DerUstardU_LR(1) = -u_LR*v_LR - sgnBC*st%K*KKK     + sgnBC*DeltaP_hug*dPsi_dU(1)/(2.0_dp*st%K)
        DerUstardU_LR(2) =  v_LR      + sgnBC*st%K*kk*u_LR + sgnBC*DeltaP_hug*dPsi_dU(2)/(2.0_dp*st%K)
        DerUstardU_LR(3) =            - sgnBC*st%K*kk       + sgnBC*DeltaP_hug*dPsi_dU(3)/(2.0_dp*st%K)
        DerUstardU_LR(4) = 0.0_dp

        DerhdU_LR(:) = Ustar*DerUstardU_LR(:)
        DerPstardU_LR_Glob(:) = -DerhdU_LR(:)/DerhdPstar
        DerRostardU_LR_Glob(:) = DerPstardU_LR_Glob(:)/(Ctank**2)

        DerUstardU_LR_Glob(:) = DerUstardU_LR(:) + st%dustar_dPstar*DerPstardU_LR_Glob(:)

      else
        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR)
        DerhdU_LR(:)=Ustar*DerUstardU_LR(:)

        DerPstardU_LR_Glob(:)=-DerhdU_LR(:)/DerhdPstar

        DerRostardU_LR_Glob(:)=DerPstardU_LR_Glob(:)/(Ctank**2)

        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR_Glob)
      end if

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

      DerRcorrstardU_LR_Glob(:)=0.0_dp

      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! DerRcorrstardU_LR_Glob(:)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      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

  else

    if(R_Correction) then

      if (use_hugoniot) then
        v_LR       = 1.0_dp / Ro_LR
        DeltaP_hug = Pstar - p_LR
        Deltav_hug = v_LR - st%vstar

        dvstar_dU(1) = st%dvstar_dv_LR*(-v_LR**2) + st%dvstar_dp_LR*(-u_LR**2/2.0_dp)
        dvstar_dU(2) = st%dvstar_dp_LR*u_LR
        dvstar_dU(3) = st%dvstar_dp_LR*(-1.0_dp)
        dvstar_dU(4) = st%dvstar_dp_LR*1.0_dp

        dPsi_dU(1) = ((-v_LR**2 - dvstar_dU(1))*DeltaP_hug + Deltav_hug*(-u_LR**2/2.0_dp)) / DeltaP_hug**2
        dPsi_dU(2) = (            -dvstar_dU(2) *DeltaP_hug + Deltav_hug*u_LR)              / DeltaP_hug**2
        dPsi_dU(3) = (            -dvstar_dU(3) *DeltaP_hug + Deltav_hug*(-1.0_dp))         / DeltaP_hug**2
        dPsi_dU(4) = (            -dvstar_dU(4) *DeltaP_hug + Deltav_hug*1.0_dp)            / DeltaP_hug**2

        DerUstardU_LR(1) = -u_LR*v_LR + sgnBC*st%K*(u_LR**2/2.0_dp) + sgnBC*DeltaP_hug*dPsi_dU(1)/(2.0_dp*st%K)
        DerUstardU_LR(2) =  v_LR      - sgnBC*st%K*u_LR              + sgnBC*DeltaP_hug*dPsi_dU(2)/(2.0_dp*st%K)
        DerUstardU_LR(3) =              sgnBC*st%K                    + sgnBC*DeltaP_hug*dPsi_dU(3)/(2.0_dp*st%K)
        DerUstardU_LR(4) =            - sgnBC*st%K                    + sgnBC*DeltaP_hug*dPsi_dU(4)/(2.0_dp*st%K)

        DerRostardU_LR(:)    = -dvstar_dU(:) * Rostar**2
        DerRcorrstardU_LR(:) = (estar - st%ev_v*st%vstar) * DerRostardU_LR(:)

      else

        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR)

        DerRostardU_LR(1)=(DerBBdU_LR(1)*(c_LR**2)-2.0_dp*BB*c_LR*dc(1))/(c_LR**4)+1.0_dp
        DerRostardU_LR(2:4)=(DerBBdU_LR(2:4)*(c_LR**2)-2.0_dp*BB*c_LR*dc(2:4))/(c_LR**4)

        DerRcorrstardU_LR(1)=-Cs_cell(Con_R)*Rostar/(Cs_cell(Con_Mas)**2)+Cs_cell(Con_R)*DerRostardU_LR(1)/Cs_cell(Con_Mas)+&
                           2.0_dp*c_LR*dc(1)*(Rostar-Cs_cell(Con_Mas))+(c_LR**2)*(DerRostardU_LR(1)-1.0_dp)
        DerRcorrstardU_LR(2:3)=Cs_cell(Con_R)*DerRostardU_LR(2:3)/Cs_cell(Con_Mas)+2.0_dp*c_LR*dc(2:3)*(Rostar-Cs_cell(Con_Mas))+&
                            (c_LR**2)*DerRostardU_LR(2:3)
        DerRcorrstardU_LR(4)=Rostar/Cs_cell(Con_Mas)+Cs_cell(Con_R)*DerRostardU_LR(4)/Cs_cell(Con_Mas)+2.0_dp*c_LR*dc(4)*&
                           (Rostar-Cs_cell(Con_Mas))+(c_LR**2)*DerRostardU_LR(4)

      end if

      DerConstardU_LR(:,1)=DerRostardU_LR(:)
      DerConstardU_LR(:,2)=Rostar*DerUstardU_LR(:)+DerRostardU_LR(:)*Ustar
      DerConstardU_LR(:,3)=DerRcorrstardU_LR(:)+0.5_dp*(DerRostardU_LR(:)*Ustar*Ustar+2.0_dp*Rostar*Ustar*DerUstardU_LR(:))
      DerConstardU_LR(:,4)=DerRcorrstardU_LR(:)
    
      me%port%derFlx_derCon(:,1)=DerConstardU_LR(:,2)
      me%port%derFlx_derCon(:,2)=DerConstardU_LR(:,2)*Ustar+Rostar*Ustar*DerUstardU_LR(:)
      me%port%derFlx_derCon(:,3)=DerConstardU_LR(:,3)*Ustar+(Rostar*(estar+0.5_dp*Ustar*Ustar)+pstar)*DerUstardU_LR(:)
      me%port%derFlx_derCon(:,4)=DerRcorrstardU_LR(:)*Ustar+R_Correction_Star*DerUstardU_LR(:)

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

    else

      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)

      if (use_hugoniot) then

        v_LR       = 1.0_dp / Ro_LR
        DeltaP_hug = Pstar - p_LR
        Deltav_hug = v_LR - st%vstar

        dvstar_dU(1) = st%dvstar_dv_LR*(-v_LR**2) + st%dvstar_dp_LR*KKK
        dvstar_dU(2) = st%dvstar_dp_LR*(-kk*u_LR)
        dvstar_dU(3) = st%dvstar_dp_LR*kk
        dvstar_dU(4) = 0.0_dp

        dPsi_dU(1) = ((-v_LR**2 - dvstar_dU(1))*DeltaP_hug + Deltav_hug*KKK)       / DeltaP_hug**2
        dPsi_dU(2) = (           -dvstar_dU(2) *DeltaP_hug + Deltav_hug*(-kk*u_LR)) / DeltaP_hug**2
        dPsi_dU(3) = (           -dvstar_dU(3) *DeltaP_hug + Deltav_hug*kk)         / DeltaP_hug**2
        dPsi_dU(4) = 0.0_dp

        DerUstardU_LR(1) = -u_LR*v_LR - sgnBC*st%K*KKK     + sgnBC*DeltaP_hug*dPsi_dU(1)/(2.0_dp*st%K)
        DerUstardU_LR(2) =  v_LR      + sgnBC*st%K*kk*u_LR + sgnBC*DeltaP_hug*dPsi_dU(2)/(2.0_dp*st%K)
        DerUstardU_LR(3) =            - sgnBC*st%K*kk       + sgnBC*DeltaP_hug*dPsi_dU(3)/(2.0_dp*st%K)
        DerUstardU_LR(4) = 0.0_dp

        DerRostardU_LR(:) = -dvstar_dU(:) * Rostar**2

      else

        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

        call derUstar_acoustic(BB, DerBBdU_LR, dc, sgnBC, c_LR, Ro_LR, u_LR, DerUstardU_LR)

        DerRostardU_LR(1)=(DerBBdU_LR(1)*(c_LR**2)-2.0_dp*BB*c_LR*dc(1))/(c_LR**4)+1.0_dp
        DerRostardU_LR(2)=(DerBBdU_LR(2)*(c_LR**2)-2.0_dp*BB*c_LR*dc(2))/(c_LR**4)
        DerRostardU_LR(3)=(DerBBdU_LR(3)*(c_LR**2)-2.0_dp*BB*c_LR*dc(3))/(c_LR**4)
        DerRostardU_LR(4)=0.0_dp

      end if

      call jacobian_roT(rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
      dR_dRo=dRstardRo+dRstardT*dTdRo_p
      DerRphystardU_LR(:)=dR_dRo*DerRostardU_LR(:)

      DerRcorrstardU_LR(:)=0.0_dp
      DerRcorrstardU_LR(4)=1.0_dp

      DerConstardU_LR(:,1)=DerRostardU_LR(:)
      DerConstardU_LR(:,2)=Rostar*DerUstardU_LR(:)+DerRostardU_LR(:)*Ustar
      DerConstardU_LR(:,3)=DerRphystardU_LR(:)+0.5_dp*(DerRostardU_LR(:)*Ustar*Ustar+2.0_dp*Rostar*Ustar*DerUstardU_LR(:))
      DerConstardU_LR(:,4)=DerRcorrstardU_LR(:)

      me%port%derFlx_derCon(:,1)=DerConstardU_LR(:,2)
      me%port%derFlx_derCon(:,2)=DerConstardU_LR(:,2)*Ustar+Rostar*Ustar*DerUstardU_LR(:)
      me%port%derFlx_derCon(:,3)=DerConstardU_LR(:,3)*Ustar+(Rostar*(estar+0.5_dp*Ustar*Ustar)+pstar)*DerUstardU_LR(:)
      me%port%derFlx_derCon(:,4)=0.0_dp

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

    endif

  endif

  if (use_hugoniot) then
    ss = abs(u_LR) + max(c_LR, st%W / Ro_LR)
  else
    ss = abs(u_LR) + c_LR
  end if
end subroutine solve_boundary_PT
solve_boundary_PT Subroutine

Contents

Source Code

public subroutine solve_boundary_PT(me, ss)

Arguments

Calls

Called by

Source Code
