Commit 4f6db54c authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

ECKEP & HES schemes improved

parent 23f4fa6a

src/FlowSolverRHEA.cpp

+67 −41
Original line number Diff line number Diff line
@@ -4589,29 +4589,43 @@ double EckepFluxApproximateRiemannSolver::calculateIntercellFlux( const double & 
        F *= w_L + w_R;

    } else if ( var_type == 4 ) {

        double bar_F_1  = ( 1.0/4.0 )*( rho_L + rho_R )*( u_L + u_R );

        double bar_F_2 = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) );

        double bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_R );

        F = bar_F_3;

        double ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( rho_R*T_R );

        double V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L;

        double V_3_R  = ( -1.0 )*rho_R*ds_drhoE_R;

        double deltaV_3 = V_3_R - V_3_L;

	if( abs( deltaV_3 ) > 1.0e-7 ) {

            double bar_F_2u = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) );

            double bar_F_2v = ( 1.0/2.0 )*bar_F_1*( v_L + v_R );

            double bar_F_2w = ( 1.0/2.0 )*bar_F_1*( w_L + w_R );

            double ds_drho_L  = ( u_L*u_L + v_L*v_L + w_L*w_L - E_L - P_L/rho_L )/( rho_L*T_L );

            double ds_drho_R  = ( u_R*u_R + v_R*v_R + w_R*w_R - E_R - P_R/rho_R )/( rho_R*T_R );

            double ds_drhou_L = ( -1.0 )*u_L/( rho_L*T_L );

            double ds_drhou_R = ( -1.0 )*u_R/( rho_R*T_R );

        double ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( rho_R*T_R );

            double ds_drhov_L = ( -1.0 )*v_L/( rho_L*T_L );

            double ds_drhov_R = ( -1.0 )*v_R/( rho_R*T_R );

            double ds_drhow_L = ( -1.0 )*w_L/( rho_L*T_L );

            double ds_drhow_R = ( -1.0 )*w_R/( rho_R*T_R );

            double V_1_L  = ( -1.0 )*( s_L + rho_L*ds_drho_L );

            double V_1_R  = ( -1.0 )*( s_R + rho_R*ds_drho_R );

        double V_2_L = ( -1.0 )*rho_L*ds_drhou_L;

        double V_2_R = ( -1.0 )*rho_R*ds_drhou_R;

        double V_3_L = ( -1.0 )*rho_L*ds_drhoE_L;

        double V_3_R = ( -1.0 )*rho_R*ds_drhoE_R;

            double V_2u_L = ( -1.0 )*rho_L*ds_drhou_L;

            double V_2u_R = ( -1.0 )*rho_R*ds_drhou_R;

            double V_2v_L = ( -1.0 )*rho_L*ds_drhov_L;

            double V_2v_R = ( -1.0 )*rho_R*ds_drhov_R;

            double V_2w_L = ( -1.0 )*rho_L*ds_drhow_L;

            double V_2w_R = ( -1.0 )*rho_R*ds_drhow_R;

            double deltaV_1 = V_1_R - V_1_L;

        double deltaV_2 = V_2_R - V_2_L;

        double deltaV_3 = V_3_R - V_3_L;

            double deltaV_2u = V_2u_R - V_2u_L;

            double deltaV_2v = V_2v_R - V_2v_L;

            double deltaV_2w = V_2w_R - V_2w_L;

            double psi_L = u_L*P_L/T_L;

            double psi_R = u_R*P_R/T_R;

            double deltaPsi = psi_R - psi_L;

        //double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon );

        double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + 1.0e-10 );	// ... modified for OpenACC

        F = bar_F_3 - ( 1.0/2.0 )*alpha_3*deltaV_3;

            double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2u*deltaV_2u + bar_F_2v*deltaV_2v + bar_F_2w*deltaV_2w + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 );

            F -= ( 1.0/2.0 )*alpha_3*deltaV_3;

	}

    }



    return( F );
@@ -4643,33 +4657,29 @@ double HesFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rh 


    /// Conserved variable increment

    double deltaU_1  = rho_R - rho_L;          

    double deltaU_2 = rho_R*u_R - rho_L*u_L;

    double deltaU_2u = rho_R*u_R - rho_L*u_L;

    double deltaU_2v = rho_R*v_R - rho_L*v_L;

    double deltaU_2w = rho_R*w_R - rho_L*w_L;

    double deltaU_3  = rho_R*E_R - rho_L*E_L;



    /// Flux increment 

    double deltaF_1  = rho_R*u_R - rho_L*u_L;

    double deltaF_2 = rho_R*u_R*u_R + P_R - rho_L*u_L*u_L - P_L;

    double deltaF_2u = rho_R*u_R*u_R + P_R - rho_L*u_L*u_L - P_L;

    double deltaF_3  = rho_R*u_R*E_R + P_R*u_R - rho_L*u_L*E_L - P_L*u_L;

   

    /// Wave speed (limited)

    //double S_1 = deltaF_1 / ( deltaU_1 + epsilon );

    double S_1 = deltaF_1 / ( deltaU_1 + 1.0e-10 );					// ... modified for OpenACC

    //if( S_1 > lambda_max ) S_1 = ( S_1/( abs( S_1 ) + epsilon ) )*lambda_max;

    if( S_1 > lambda_max ) S_1 = ( S_1/( abs( S_1 ) + 1.0e-10 ) )*lambda_max;		// ... modified for OpenACC

    //if( S_1 < lambda_min ) S_1 = ( S_1/( abs( S_1 ) + epsilon ) )*lambda_min;

    if( S_1 < lambda_min ) S_1 = ( S_1/( abs( S_1 ) + 1.0e-10 ) )*lambda_min;		// ... modified for OpenACC

    //double S_2 = deltaF_2/( deltaU_2 + epsilon );

    double S_2 = deltaF_2/( deltaU_2 + 1.0e-10 );					// ... modified for OpenACC

    //if( S_2 > lambda_max ) S_2 = ( S_2/( abs( S_2 ) + epsilon ) )*lambda_max;

    if( S_2 > lambda_max ) S_2 = ( S_2/( abs( S_2 ) + 1.0e-10 ) )*lambda_max;		// ... modified for OpenACC

    //if( S_2 < lambda_min ) S_2 = ( S_2/( abs( S_2 ) + epsilon ) )*lambda_min;

    if( S_2 < lambda_min ) S_2 = ( S_2/( abs( S_2 ) + 1.0e-10 ) )*lambda_min;		// ... modified for OpenACC

    if( abs(S_1) > lambda_max ) S_1 = copysign( lambda_max, S_1 );

    if( abs(S_1) < lambda_min ) S_1 = copysign( lambda_min, S_1 );

    //double S_2 = deltaF_2u / ( deltaU_2u + epsilon );

    double S_2 = deltaF_2u / ( deltaU_2u + 1.0e-10 );					// ... modified for OpenACC

    if( abs(S_2) > lambda_max ) S_2 = copysign( lambda_max, S_2 );

    if( abs(S_2) < lambda_min ) S_2 = copysign( lambda_min, S_2 );

    //double S_3 = deltaF_3 / ( deltaU_3 + epsilon );

    double S_3 = deltaF_3 / ( deltaU_3 + 1.0e-10 );					// ... modified for OpenACC

    //if( S_3 > lambda_max ) S_3 = ( S_3/( abs( S_3 ) + epsilon ) )*lambda_max;

    if( S_3 > lambda_max ) S_3 = ( S_3/( abs( S_3 ) + 1.0e-10 ) )*lambda_max;		// ... modified for OpenACC

    //if( S_3 < lambda_min ) S_3 = ( S_3/( abs( S_3 ) + epsilon ) )*lambda_min;

    if( S_3 < lambda_min ) S_3 = ( S_3/( abs( S_3 ) + 1.0e-10 ) )*lambda_min;		// ... modified for OpenACC

    if( abs(S_3) > lambda_max ) S_3 = copysign( lambda_max, S_3 );

    if( abs(S_3) < lambda_min ) S_3 = copysign( lambda_min, S_3 );

    double alpha_S = min( abs(S_1), min( abs(S_2), abs(S_3) ) );



    /// Entropic variables: derivatives
@@ -4782,36 +4792,52 @@ double HesFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rh 
        F -= ( 1.0/2.0 )*alpha_S*deltaU_1;

    } else if ( var_type == 1 ) {

        F *= u_L + u_R; F += ( 1.0/2.0 )*( P_L + P_R );

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2;

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2u;

    } else if ( var_type == 2 ) {

        F *= v_L + v_R;

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2v;

    } else if ( var_type == 3 ) {

        F *= w_L + w_R;

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2w;

    } else if ( var_type == 4 ) {

        double bar_F_1  = ( 1.0/4.0 )*( rho_L + rho_R )*( u_L + u_R );

        double bar_F_2 = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) );

        double bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_R );

        F = bar_F_3;

        double ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( rho_R*T_R );

        double V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L;

        double V_3_R  = ( -1.0 )*rho_R*ds_drhoE_R;

        double deltaV_3 = V_3_R - V_3_L;

	if( abs( deltaV_3 ) > 1.0e-7 ) {

            double bar_F_2u = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) );

            double bar_F_2v = ( 1.0/2.0 )*bar_F_1*( v_L + v_R );

            double bar_F_2w = ( 1.0/2.0 )*bar_F_1*( w_L + w_R );

            double ds_drho_L  = ( u_L*u_L + v_L*v_L + w_L*w_L - E_L - P_L/rho_L )/( rho_L*T_L );

            double ds_drho_R  = ( u_R*u_R + v_R*v_R + w_R*w_R - E_R - P_R/rho_R )/( rho_R*T_R );

            double ds_drhou_L = ( -1.0 )*u_L/( rho_L*T_L );

            double ds_drhou_R = ( -1.0 )*u_R/( rho_R*T_R );

        double ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( rho_R*T_R );

            double ds_drhov_L = ( -1.0 )*v_L/( rho_L*T_L );

            double ds_drhov_R = ( -1.0 )*v_R/( rho_R*T_R );

            double ds_drhow_L = ( -1.0 )*w_L/( rho_L*T_L );

            double ds_drhow_R = ( -1.0 )*w_R/( rho_R*T_R );

            double V_1_L  = ( -1.0 )*( s_L + rho_L*ds_drho_L );

            double V_1_R  = ( -1.0 )*( s_R + rho_R*ds_drho_R );

        double V_2_L = ( -1.0 )*rho_L*ds_drhou_L;

        double V_2_R = ( -1.0 )*rho_R*ds_drhou_R;

        double V_3_L = ( -1.0 )*rho_L*ds_drhoE_L;

        double V_3_R = ( -1.0 )*rho_R*ds_drhoE_R;

            double V_2u_L = ( -1.0 )*rho_L*ds_drhou_L;

            double V_2u_R = ( -1.0 )*rho_R*ds_drhou_R;

            double V_2v_L = ( -1.0 )*rho_L*ds_drhov_L;

            double V_2v_R = ( -1.0 )*rho_R*ds_drhov_R;

            double V_2w_L = ( -1.0 )*rho_L*ds_drhow_L;

            double V_2w_R = ( -1.0 )*rho_R*ds_drhow_R;

            double deltaV_1 = V_1_R - V_1_L;

        double deltaV_2 = V_2_R - V_2_L;

        double deltaV_3 = V_3_R - V_3_L;

            double deltaV_2u = V_2u_R - V_2u_L;

            double deltaV_2v = V_2v_R - V_2v_L;

            double deltaV_2w = V_2w_R - V_2w_L;

            double psi_L = u_L*P_L/T_L;

            double psi_R = u_R*P_R/T_R;

            double deltaPsi = psi_R - psi_L;

        //double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon );

        double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + 1.0e-10 );	// ... modified for OpenACC

        F = bar_F_3 - ( 1.0/2.0 )*alpha_3*deltaV_3;

            double alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2u*deltaV_2u + bar_F_2v*deltaV_2v + bar_F_2w*deltaV_2w + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 );

            F -= ( 1.0/2.0 )*alpha_3*deltaV_3;

	}

        F -= ( 1.0/2.0 )*alpha_S*deltaU_3;

    }

stuff/rhea_flow_solver.py

+67 −129
Original line number Diff line number Diff line
@@ -997,28 +997,42 @@ def ECKEP_flux( rho_L, rho_R, u_L, u_R, v_L, v_R, w_L, w_R, E_L, E_R, s_L, s_R, 
        F *= w_L + w_R

    elif ( var_type == 4 ):

        bar_F_1  = ( 1.0/4.0 )*( rho_L + rho_R )*( u_L + u_R )

        bar_F_2 = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) )

        bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_R )

        F = bar_F_3

        ds_drhoE_L = 1.0/( rho_L*T_L )

        ds_drhoE_R = 1.0/( rho_R*T_R )

        V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L

        V_3_R  = ( -1.0 )*rho_R*ds_drhoE_R

        deltaV_3 = V_3_R - V_3_L

	    if( abs( deltaV_3 ) > 1.0e-7 ):

            bar_F_2u = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) )

            bar_F_2v = ( 1.0/2.0 )*bar_F_1*( v_L + v_R )

            bar_F_2w = ( 1.0/2.0 )*bar_F_1*( w_L + w_R )

            ds_drho_L  = ( u_L*u_L + v_L*v_L + w_L*w_L - E_L - P_L/rho_L )/( rho_L*T_L )

            ds_drho_R  = ( u_R*u_R + v_R*v_R + w_R*w_R - E_R - P_R/rho_R )/( rho_R*T_R )

            ds_drhou_L = ( -1.0 )*u_L/( rho_L*T_L )

            ds_drhou_R = ( -1.0 )*u_R/( rho_R*T_R )

        ds_drhoE_L = 1.0/( rho_L*T_L )

        ds_drhoE_R = 1.0/( rho_R*T_R )

            ds_drhov_L = ( -1.0 )*v_L/( rho_L*T_L )

            ds_drhov_R = ( -1.0 )*v_R/( rho_R*T_R )

            ds_drhow_L = ( -1.0 )*w_L/( rho_L*T_L )

            ds_drhow_R = ( -1.0 )*w_R/( rho_R*T_R )

            V_1_L  = ( -1.0 )*( s_L + rho_L*ds_drho_L )

            V_1_R  = ( -1.0 )*( s_R + rho_R*ds_drho_R )

        V_2_L = ( -1.0 )*rho_L*ds_drhou_L

        V_2_R = ( -1.0 )*rho_R*ds_drhou_R

        V_3_L = ( -1.0 )*rho_L*ds_drhoE_L

        V_3_R = ( -1.0 )*rho_R*ds_drhoE_R

            V_2u_L = ( -1.0 )*rho_L*ds_drhou_L

            V_2u_R = ( -1.0 )*rho_R*ds_drhou_R

            V_2v_L = ( -1.0 )*rho_L*ds_drhov_L

            V_2v_R = ( -1.0 )*rho_R*ds_drhov_R

            V_2w_L = ( -1.0 )*rho_L*ds_drhow_L

            V_2w_R = ( -1.0 )*rho_R*ds_drhow_R

            deltaV_1 = V_1_R - V_1_L

        deltaV_2 = V_2_R - V_2_L

        deltaV_3 = V_3_R - V_3_L

            deltaV_2u = V_2u_R - V_2u_L

            deltaV_2v = V_2v_R - V_2v_L

            deltaV_2w = V_2w_R - V_2w_L

            psi_L = u_L*P_L/T_L

            psi_R = u_R*P_R/T_R

            deltaPsi = psi_R - psi_L

        alpha_3  = ( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon )

        F = bar_F_3 - ( 1.0/2.0 )*alpha_3*deltaV_3

            alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2u*deltaV_2u + bar_F_2v*deltaV_2v + bar_F_2w*deltaV_2w + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 )

            F -= ( 1.0/2.0 )*alpha_3*deltaV_3



    return( F )
@@ -1043,31 +1057,33 @@ def ECKEP_MOVERS_RH_flux( rho_L, rho_R, u_L, u_R, v_L, v_R, w_L, w_R, E_L, E_R, 


    # Conserved variable increment

    deltaU_1  = rho_R - rho_L            

    deltaU_2 = rho_R*u_R - rho_L*u_L

    deltaU_2u = rho_R*u_R - rho_L*u_L

    deltaU_2v = rho_R*v_R - rho_L*v_L

    deltaU_2w = rho_R*w_R - rho_L*w_L

    deltaU_3  = rho_R*E_R - rho_L*E_L



    # Flux increment 

    deltaF_1  = rho_R*u_R - rho_L*u_L

    deltaF_2 = rho_R*u_R*u_R + P_R - rho_L*u_L*u_L - P_L

    deltaF_2u = rho_R*u_R*u_R + P_R - rho_L*u_L*u_L - P_L

    deltaF_3  = rho_R*u_R*E_R + P_R*u_R - rho_L*u_L*E_L - P_L*u_L

   

    # Wave speed (limited)

    S_1 = deltaF_1 / (deltaU_1 + epsilon)

    if( S_1 > lambda_max ):

        S_1 = ( S_1/( abs( S_1 ) + epsilon ) )*lambda_max

    if( S_1 < lambda_min ):

        S_1 = ( S_1/( abs( S_1 ) + epsilon ) )*lambda_min

    S_2 = deltaF_2/( deltaU_2 + epsilon )

    if( S_2 > lambda_max ):

        S_2 = ( S_2/( abs( S_2 ) + epsilon ) )*lambda_max

    if( S_2 < lambda_min ):

        S_2 = ( S_2/( abs( S_2 ) + epsilon ) )*lambda_min

    if abs(S_1) > lambda_max:

        S_1 = math.copysign(lambda_max, S_1)

    if abs(S_1) < lambda_min:

        S_1 = math.copysign(lambda_min, S_1)

    S_2 = deltaF_2u / (deltaU_2 + epsilon)

    if abs(S_2) > lambda_max:

        S_2 = math.copysign(lambda_max, S_2)

    if abs(S_2) < lambda_min:

        S_2 = math.copysign(lambda_min, S_2)

    S_3 = deltaF_3 / (deltaU_3 + epsilon)

    if( S_3 > lambda_max ):

        S_3 = ( S_3/( abs( S_3 ) + epsilon ) )*lambda_max

    if( S_3 < lambda_min ):

        S_3 = ( S_3/( abs( S_3 ) + epsilon ) )*lambda_min

    alpha_S = min( S_1, S_2, S_3 )

    if abs(S_3) > lambda_max:

        S_3 = math.copysign(lambda_max, S_3)

    if abs(S_3) < lambda_min:

        S_3 = math.copysign(lambda_min, S_3)

    alpha_S = min(abs(S_1), abs(S_2), abs(S_3))



    ### ---------------------------------###

    ### START: SHOCK SENSOR MODIFICATION ###
@@ -1101,131 +1117,53 @@ def ECKEP_MOVERS_RH_flux( rho_L, rho_R, u_L, u_R, v_L, v_R, w_L, w_R, E_L, E_R, 
    elif ( var_type == 1 ):

        F *= u_L + u_R

        F += ( 1.0/2.0 )*( P_L + P_R )

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2u

    elif ( var_type == 2 ):

        F *= v_L + v_R

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2v

    elif ( var_type == 3 ):

        F *= w_L + w_R

        F -= ( 1.0/2.0 )*alpha_S*deltaU_2w

    elif ( var_type == 4 ):

        bar_F_1  = ( 1.0/4.0 )*( rho_L + rho_R )*( u_L + u_R )

        bar_F_2 = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) )

        bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_R )

        ds_drho_L  = ( u_L*u_L + v_L*v_L + w_L*w_L - E_L - P_L/rho_L )/( rho_L*T_L )

        ds_drho_R  = ( u_R*u_R + v_R*v_R + w_R*w_R - E_R - P_R/rho_R )/( rho_R*T_R )

        ds_drhou_L = ( -1.0 )*u_L/( rho_L*T_L )

        ds_drhou_R = ( -1.0 )*u_R/( rho_R*T_R )

        F = bar_F_3

        ds_drhoE_L = 1.0/( rho_L*T_L )

        ds_drhoE_R = 1.0/( rho_R*T_R )

        V_1_L = ( -1.0 )*( s_L + rho_L*ds_drho_L ) 

        V_1_R = ( -1.0 )*( s_R + rho_R*ds_drho_R ) 

        V_2_L = ( -1.0 )*rho_L*ds_drhou_L

        V_2_R = ( -1.0 )*rho_R*ds_drhou_R

        V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L

        V_3_R  = ( -1.0 )*rho_R*ds_drhoE_R

        deltaV_1 = V_1_R - V_1_L

        deltaV_2 = V_2_R - V_2_L

        deltaV_3 = V_3_R - V_3_L

        psi_L = u_L*P_L/T_L

        psi_R = u_R*P_R/T_R

        deltaPsi = psi_R - psi_L

        alpha_3  = ( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon )

        F = bar_F_3 - ( 1.0/2.0 )*alpha_3*deltaV_3

        F -= ( 1.0/2.0 )*alpha_S*deltaU_3

    

    return( F )





### Calculate ECKEP_HLLC flux ... var_type corresponds to: 0 for rho, 1-3 for rhouvw, 4 for rhoE 

@njit

def ECKEP_HLLC_flux( rho_L, rho_R, u_L, u_R, v_L, v_R, w_L, w_R, E_L, E_R, s_L, s_R, P_L, P_R, T_L, T_R, a_L, a_R, var_type ):



    # Calculate waves speed

    S_L, S_R = waves_speed( rho_L, rho_R, u_L, u_R, P_L, P_R, a_L, a_R )

    S_star   = ( P_R - P_L + rho_L*u_L*( S_L - u_L ) - rho_R*u_R*( S_R - u_R ) )/( rho_L*( S_L - u_L ) - rho_R*( S_R - u_R ) + epsilon )

    

    # Calculate ECKEP flux + HLLC L-R fluxes & states

    F_ECKEP = ( 1.0/8.0 )*( rho_L + rho_R )*( u_L + u_R )

    F_L = rho_L*u_L

    F_R = rho_R*u_R

    U_L = rho_L

    U_R = rho_R

    U_star_L = rho_L*( ( S_L - u_L )/( S_L - S_star + epsilon ) )

    U_star_R = rho_R*( ( S_R - u_R )/( S_R - S_star + epsilon ) )

    if( var_type == 0 ):

        F_ECKEP *= 1.0 + 1.0

        F_L *= 1.0

        F_R *= 1.0        

        U_L *= 1.0

        U_R *= 1.0

        U_star_L *= 1.0

        U_star_R *= 1.0        

    elif ( var_type == 1 ):

        F_ECKEP *= u_L + u_R

        F_ECKEP += ( 1.0/2.0 )*( P_L + P_R )

        F_L *= u_L; F_L += P_L

        F_R *= u_R; F_R += P_R        

        U_L *= u_L

        U_R *= u_R

        U_star_L *= S_star

        U_star_R *= S_star        

    elif ( var_type == 2 ):

        F_ECKEP *= v_L + v_R

        F_L *= v_L

        F_R *= v_R        

        U_L *= v_L

        U_R *= v_R

        U_star_L *= v_L

        U_star_R *= v_R        

    elif ( var_type == 3 ):

        F_ECKEP *= w_L + w_R

        F_L *= w_L

        F_R *= w_R        

        U_L *= w_L

        U_R *= w_R

        U_star_L *= w_L

        U_star_R *= w_R        

    elif ( var_type == 4 ):

        bar_F_1 = ( 1.0/4.0 )*( rho_L + rho_R )*( u_L + u_R )

        bar_F_2 = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) )

        bar_F_3 = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_R )

	    if( abs( deltaV_3 ) > 1.0e-7 ):

            bar_F_2u = ( 1.0/2.0 )*( bar_F_1*( u_L + u_R ) + ( P_L + P_R ) )

            bar_F_2v = ( 1.0/2.0 )*bar_F_1*( v_L + v_R )

            bar_F_2w = ( 1.0/2.0 )*bar_F_1*( w_L + w_R )

            ds_drho_L  = ( u_L*u_L + v_L*v_L + w_L*w_L - E_L - P_L/rho_L )/( rho_L*T_L )

            ds_drho_R  = ( u_R*u_R + v_R*v_R + w_R*w_R - E_R - P_R/rho_R )/( rho_R*T_R )

            ds_drhou_L = ( -1.0 )*u_L/( rho_L*T_L )

            ds_drhou_R = ( -1.0 )*u_R/( rho_R*T_R )

        ds_drhoE_L = 1.0/( rho_L*T_L )

        ds_drhoE_R = 1.0/( rho_R*T_R )

            ds_drhov_L = ( -1.0 )*v_L/( rho_L*T_L )

            ds_drhov_R = ( -1.0 )*v_R/( rho_R*T_R )

            ds_drhow_L = ( -1.0 )*w_L/( rho_L*T_L )

            ds_drhow_R = ( -1.0 )*w_R/( rho_R*T_R )

            V_1_L  = ( -1.0 )*( s_L + rho_L*ds_drho_L )

            V_1_R  = ( -1.0 )*( s_R + rho_R*ds_drho_R )

        V_2_L = ( -1.0 )*rho_L*ds_drhou_L

        V_2_R = ( -1.0 )*rho_R*ds_drhou_R

        V_3_L = ( -1.0 )*rho_L*ds_drhoE_L

        V_3_R = ( -1.0 )*rho_R*ds_drhoE_R

            V_2u_L = ( -1.0 )*rho_L*ds_drhou_L

            V_2u_R = ( -1.0 )*rho_R*ds_drhou_R

            V_2v_L = ( -1.0 )*rho_L*ds_drhov_L

            V_2v_R = ( -1.0 )*rho_R*ds_drhov_R

            V_2w_L = ( -1.0 )*rho_L*ds_drhow_L

            V_2w_R = ( -1.0 )*rho_R*ds_drhow_R

            deltaV_1 = V_1_R - V_1_L

        deltaV_2 = V_2_R - V_2_L

        deltaV_3 = V_3_R - V_3_L

            deltaV_2u = V_2u_R - V_2u_L

            deltaV_2v = V_2v_R - V_2v_L

            deltaV_2w = V_2w_R - V_2w_L

            psi_L = u_L*P_L/T_L

            psi_R = u_R*P_R/T_R

            deltaPsi = psi_R - psi_L

        alpha_3  = ( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon )

        F_ECKEP = bar_F_3 - ( 1.0/2.0 )*alpha_3*deltaV_3

        F_L *= E_L; F_L += u_L*P_L

        F_R *= E_R; F_R += u_R*P_R

        U_L *= E_L

        U_R *= E_R

        U_star_L *= ( E_L + ( S_star - u_L )*( S_star + P_L/( rho_L*( S_L - u_L ) + epsilon ) ) )

        U_star_R *= ( E_R + ( S_star - u_R )*( S_star + P_R/( rho_R*( S_R - u_R ) + epsilon ) ) )        

       

    # Calculate HLLC flux

    F_HLLC = (1.0/2.0)*( ( 1.0 + np.sign( S_star ) )*( F_L + S_L*( U_star_L - U_L ) ) + ( 1.0 - np.sign( S_star ) )*( F_R + S_R*( U_star_R - U_R ) ) ) 



    # Pressure-based shock indicator: shocks have sharp pressure jump

    P = ( 1.0/2.0 )*( P_L + P_R )

    phi = min( 1.0, 50.0*abs( P_R - P_L )/( P + epsilon ) )     # Factor: 10, 50, 100



    # Calculate ECKEP-HLLC flux

    F = ( 1.0 - phi )*F_ECKEP + phi*F_HLLC 

            alpha_3  = 2.0*( bar_F_1*deltaV_1 + bar_F_2u*deltaV_2u + bar_F_2v*deltaV_2v + bar_F_2w*deltaV_2w + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 )

            F -= ( 1.0/2.0 )*alpha_3*deltaV_3

        F -= ( 1.0/2.0 )*alpha_S*deltaU_3

    

    # Return F value

    return( F )