Commit b645d76a authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

Implementation of HES flux scheme IV

parent 4f6db54c

src/FlowSolverRHEA.cpp

+20 −22
Original line number Diff line number Diff line
@@ -4576,7 +4576,9 @@ EckepFluxApproximateRiemannSolver::~EckepFluxApproximateRiemannSolver() {}; 
double EckepFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rho_LLL, const double &u_LLL, const double &v_LLL, const double &w_LLL, const double &E_LLL, const double &s_LLL, const double &P_LLL, const double &T_LLL, const double &a_LLL, const double &rho_LL, const double &u_LL, const double &v_LL, const double &w_LL, const double &E_LL, const double &s_LL, const double &P_LL, const double &T_LL, const double &a_LL, const double &rho_L, const double &u_L, const double &v_L, const double &w_L, const double &E_L, const double &s_L, const double &P_L, const double &T_L, const double &a_L, const double &rho_R, const double &u_R, const double &v_R, const double &w_R, const double &E_R, const double &s_R, const double &P_R, const double &T_R, const double &a_R, const double &rho_RR, const double &u_RR, const double &v_RR, const double &w_RR, const double &E_RR, const double &s_RR, const double &P_RR, const double &T_RR, const double &a_RR, const double &rho_RRR, const double &u_RRR, const double &v_RRR, const double &w_RRR, const double &E_RRR, const double &s_RRR, const double &P_RRR, const double &T_RRR, const double &a_RRR, const double &delta, const int &var_type ) {



    /// Entropy conservative and kinetic energy preserving (ECKEP) scheme:

    /// K. Bahuguna, R. Kolluru, S.V. R. Rao

    /// Structure-preserving schemes conserving entropy and kinetic energy.

    /// arXiv:2505.13374, 2025.



    double F = ( 1.0/8.0 )*( rho_L + rho_R )*( u_L + u_R );

    if( var_type == 0 ) {
@@ -4589,17 +4591,10 @@ 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_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 bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_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 );
@@ -4608,6 +4603,8 @@ double EckepFluxApproximateRiemannSolver::calculateIntercellFlux( const double & 
        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 ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( 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_2u_L = ( -1.0 )*rho_L*ds_drhou_L;
@@ -4616,16 +4613,19 @@ double EckepFluxApproximateRiemannSolver::calculateIntercellFlux( const double & 
        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 V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L;

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

        double deltaV_1 = V_1_R - V_1_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 deltaV_3 = V_3_R - V_3_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_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;

	}

        //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 + epsilon );

        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 + 1.0e-10 );	// ... modified for OpenACC

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

    }



    return( F );
@@ -4801,17 +4801,10 @@ double HesFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rh 
        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_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 bar_F_3  = ( 1.0/2.0 )*bar_F_1*( E_L + P_L/rho_L + E_R + P_R/rho_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 );
@@ -4820,6 +4813,8 @@ double HesFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rh 
        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 ds_drhoE_L = 1.0/( rho_L*T_L );

        double ds_drhoE_R = 1.0/( 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_2u_L = ( -1.0 )*rho_L*ds_drhou_L;
@@ -4828,16 +4823,19 @@ double HesFluxApproximateRiemannSolver::calculateIntercellFlux( const double &rh 
        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 V_3_L  = ( -1.0 )*rho_L*ds_drhoE_L;

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

        double deltaV_1 = V_1_R - V_1_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 deltaV_3 = V_3_R - V_3_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_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;

	}

        //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 + epsilon );

        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 + 1.0e-10 );	// ... modified for OpenACC

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

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

    }

stuff/rhea_flow_solver.py

+21 −20
Original line number Diff line number Diff line
@@ -985,6 +985,11 @@ def KGP_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_ 
@njit

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, P_L, P_R, T_L, T_R, a_L, a_R, var_type ):

   

    # Entropy conservative and kinetic energy preserving (ECKEP) scheme:

    # K. Bahuguna, R. Kolluru, S.V. R. Rao

    # Structure-preserving schemes conserving entropy and kinetic energy.

    # arXiv:2505.13374, 2025.



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

    if( var_type == 0 ):

        F *= 1.0 + 1.0
@@ -997,17 +1002,10 @@ 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_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 )

        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 )
@@ -1016,6 +1014,8 @@ 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, 
        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 )

        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_2u_L = ( -1.0 )*rho_L*ds_drhou_L
@@ -1024,15 +1024,18 @@ 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, 
        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

        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_2u = V_2u_R - V_2u_L

        deltaV_2v = V_2v_R - V_2v_L

        deltaV_2w = V_2w_R - V_2w_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  = 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

        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 + epsilon )

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



    return( F )
@@ -1126,17 +1129,10 @@ 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, 
        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_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 )

        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 )
@@ -1145,6 +1141,8 @@ 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, 
        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 )

        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_2u_L = ( -1.0 )*rho_L*ds_drhou_L
@@ -1153,15 +1151,18 @@ 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, 
        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

        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_2u = V_2u_R - V_2u_L

        deltaV_2v = V_2v_R - V_2v_L

        deltaV_2w = V_2w_R - V_2w_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  = 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

        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 + epsilon )

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

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

    

    return( F )