Junction_resolution_from_and_to_ports Subroutine

public subroutine Junction_resolution_from_and_to_ports(me)

Called at the beginning of main loop

Arguments

Type IntentOptional Attributes Name
type(junction_t), intent(inout) :: me

Calls

proc~~junction_resolution_from_and_to_ports~~CallsGraph proc~junction_resolution_from_and_to_ports Junction_resolution_from_and_to_ports proc~set_error set_error proc~junction_resolution_from_and_to_ports->proc~set_error proc~solve_junction solve_junction proc~junction_resolution_from_and_to_ports->proc~solve_junction proc~solve_junction->proc~set_error proc~droeint_drop droeint_droP proc~solve_junction->proc~droeint_drop proc~fsolve fsolve proc~solve_junction->proc~fsolve proc~inv inv proc~solve_junction->proc~inv proc~jacobian_rot jacobian_roT proc~solve_junction->proc~jacobian_rot proc~jacobian_star_loc jacobian_star_loc proc~solve_junction->proc~jacobian_star_loc proc~jacobian_star_req_rem_loc jacobian_star_Req_rem_loc proc~solve_junction->proc~jacobian_star_req_rem_loc proc~r_rot r_roT proc~solve_junction->proc~r_rot proc~state_rop state_roP proc~solve_junction->proc~state_rop proc~t_roe T_roE proc~solve_junction->proc~t_roe proc~eos_terms eos_terms proc~droeint_drop->proc~eos_terms proc~t_rop T_roP proc~droeint_drop->proc~t_rop proc~hybrd hybrd proc~fsolve->proc~hybrd proc~inv->proc~set_error dgetrf dgetrf proc~inv->dgetrf dgetri dgetri proc~inv->dgetri proc~jacobian_rot->proc~eos_terms proc~jacobian_star_req_rem_loc->proc~jacobian_rot proc~fill_f_terms fill_f_terms proc~r_rot->proc~fill_f_terms proc~fill_ff_terms fill_ff_terms proc~r_rot->proc~fill_ff_terms proc~fill_g_dreg_terms fill_g_Dreg_terms proc~r_rot->proc~fill_g_dreg_terms proc~state_rop->proc~jacobian_rot proc~state_rop->proc~t_rop proc~brent brent proc~t_roe->proc~brent proc~fill_e_tpart fill_e_Tpart proc~t_roe->proc~fill_e_tpart proc~t_roe->proc~fill_ff_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->proc~fill_g_dreg_terms fcn fcn proc~hybrd->fcn proc~dogleg dogleg proc~hybrd->proc~dogleg proc~enorm eNorm proc~hybrd->proc~enorm proc~fdjac1 fdJac1 proc~hybrd->proc~fdjac1 proc~qform qform proc~hybrd->proc~qform proc~qrfac qrFac proc~hybrd->proc~qrfac proc~r1mpyq r1mpyq proc~hybrd->proc~r1mpyq proc~r1updt r1updt proc~hybrd->proc~r1updt proc~t_rop->proc~brent proc~t_rop->proc~fill_f_terms proc~t_rop->proc~fill_g_dreg_terms proc~dogleg->proc~enorm proc~eos_e_terms->proc~fill_e_tpart proc~eos_e_terms->proc~fill_ff_terms proc~fdjac1->fcn proc~qrfac->proc~enorm proc~zero->f

Called by

proc~~junction_resolution_from_and_to_ports~~CalledByGraph proc~junction_resolution_from_and_to_ports Junction_resolution_from_and_to_ports proc~main_loop main_loop proc~main_loop->proc~junction_resolution_from_and_to_ports program~reims_p reims_p program~reims_p->proc~main_loop

Source Code

subroutine Junction_resolution_from_and_to_ports(me)
    !! Called at the beginning of main loop
    type(junction_t), intent(inout) :: me

    real(dp), parameter :: epsM = 1.0e-12_dp
    real(dp) :: m_dot(me%NbTotBr), angle(me%NbTotBr)
    integer :: i, i_ref, i_conv(me%NbTotBr), i_range(me%NbTotBr)
    logical :: in(me%NbTotBr), out(me%NbTotBr)

    ! calculation of the mass flow rate per branch to determine inlets, outlets
    in = .false.
    do i = 1, size(me%br); associate(prt => me%br(i)%p) 
        m_dot(i) = prt%Sgn4j * prt%Prim(Pri_u) * prt%Prim(Pri_ro) * prt%Area
        if (m_dot(i) > epsM) in(i) = .true. ! Identification of the incoming branches
    end associate; enddo
    if (any(ieee_is_nan(m_dot))) then; call set_error('NaN in junction m_dot (junction)'); return; end if
    i_ref = maxLoc(m_dot, 1) ! reference branch should be with highest mass flow rate
    if(m_dot(i_ref) <= epsM) in = .true. ! no inlet branches so make them all inlet

    ! calculate conversion index array
    out = .not.in
    in(i_ref) = .false. ! exclude reference branch from incoming
    i_range = [(i, i = 1, size(me%br))]
    i_conv = [pack(i_range, out), i_ref, pack(i_range, in)]

    angle = abs(me%br(i_ref)%angle - me%br%angle) / 180.0_dp * Pi_value
    
    me%dynamic%NbTotBr = size(me%br)
    me%dynamic%NbIn = count(in)+1
    me%dynamic%NbOut = count(out)
     me%dynamic%kappa = me%kappa
    me%dynamic%idxKappaRef = findLoc(i_conv, 1, 1)   
    do i = 1, size(me%br); associate(dyn => me%dynamic%br(i), prt => me%br(i_conv(i))%p)
        dyn%rho0 = prt%Prim(Pri_ro)
        dyn%vit0 = -prt%Sgn4j * prt%Prim(Pri_u)
        dyn%p0 = prt%Prim(Pri_p)
        dyn%ETot0 = prt%Prim(Pri_e) + 0.5_dp * prt%Prim(Pri_u) * prt%Prim(Pri_u)
        dyn%rCorr0 = prt%Prim(Pri_R)
        dyn%SpeedS0 = abs(prt%Prim(Pri_u)) + prt%Prim(Pri_c)
        dyn%CSound0 = prt%Prim(Pri_c)
        dyn%T0 = prt%Prim(Pri_T)
        dyn%AreaJGB = prt%Area
        dyn%ThetaJGB = angle(i_conv(i))
        dyn%SgnJGB = -prt%Sgn4j
        dyn%U1 = prt%Cons(Con_Mas)
        dyn%U2 = -prt%Sgn4j * prt%Cons(Con_Qdm)
        dyn%U3 = prt%Cons(Con_Ene)
        dyn%U4 = prt%Cons(Con_R)
        dyn%q0 = 0.0_dp
        dyn%qBar = 0.0_dp
        dyn%dx0 = prt%dxLoc
        dyn%Frc      = -prt%Sgn4j * prt%Su_fric
        dyn%DerFric1 = -prt%Sgn4j * prt%DerSu_fric(Con_Mas)
        dyn%DerFric2 =              prt%DerSu_fric(Con_Qdm)
        dyn%DerFric3 = -prt%Sgn4j * prt%DerSu_fric(Con_Ene)
        dyn%DerFric4 = -prt%Sgn4j * prt%DerSu_fric(Con_R)
        if(R_Correction) then
            dyn%FrcR      = prt%Sr_fric
            dyn%DerFricR1 = prt%DerSr_fric(Con_Mas)
            dyn%DerFricR2 = -prt%Sgn4j * prt%DerSr_fric(Con_Qdm)
            dyn%DerFricR3 = prt%DerSr_fric(Con_Ene)
            dyn%DerFricR4 = prt%DerSr_fric(Con_R)
        endif
    end associate; enddo

    call solve_junction(me%dynamic,me%FlJ)
    if (sim_error > 0) return

    do i = 1, size(me%br)
        me%DerFlx_dBr(:,:,i_conv,i_conv(i)) = me%FlJ(i)%DerFlx
        me%DerVel_dBr(:,  i_conv,i_conv(i)) = me%FlJ(i)%DerVit
    enddo

    do i = 1, size(me%br)
        me%br(i_conv(i))%p%flx = me%FlJ(i)%Flx
        me%br(i_conv(i))%p%vit = me%FlJ(i)%Vit
        me%br(i)%p%derFlx_derCon = me%DerFlx_dBr(:,:,i,i)
        me%br(i)%p%derVit_derCon = me%DerVel_dBr(:,i,i)
    enddo

    me%wave_time = min(minVal(me%dynamic%br%dx0 / me%dynamic%br%SpeedS0), 1e30_dp)
end subroutine Junction_resolution_from_and_to_ports