Commit 28d33f28 authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

New Riemann solver added

parent 308700d4

Makefile

+2 −0
Original line number Diff line number Diff line
@@ -69,6 +69,8 @@ else ifeq ($(EXTRACTED_RIEMANN_SOLVER_SCHEME), HLLC+) 
        RIEMANN_SOLVER_SCHEME = HllcPlusApproximateRiemannSolver

else ifeq ($(EXTRACTED_RIEMANN_SOLVER_SCHEME), ECKEP)

        RIEMANN_SOLVER_SCHEME = EckepFluxApproximateRiemannSolver

else ifeq ($(EXTRACTED_RIEMANN_SOLVER_SCHEME), ECKEP-MOVERS-RH)

        RIEMANN_SOLVER_SCHEME = EckepMoversRhFluxApproximateRiemannSolver

endif

EXTRACTED_RUNGE_KUTTA_METHOD := $(shell grep "runge_kutta_time_scheme" $(CONFIG_FILE) | sed -n "s/.*'\([^']*\)'.*/\1/p")

EXPLICIT_RUNGE_KUTTA_METHOD = dummy

configuration_file.yaml

+1 −1
Original line number Diff line number Diff line
@@ -65,7 +65,7 @@ computational_parameters: 
   external_mesh: 'FALSE'                                               # Activate external mesh generation

   external_mesh_file: 'external_mesh_file.txt'                         # Name of external mesh file   

   CFL: 0.9                  					        # CFL coefficient

   # Riemann solver scheme: DIVERGENCE, MURMAN-ROE, KGP, SHIMA, HLL, HLLC, HLLC+, ECKEP

   # Riemann solver scheme: DIVERGENCE, MURMAN-ROE, KGP, SHIMA, HLL, HLLC, HLLC+, ECKEP, ECKEP-MOVERS-RH

   riemann_solver_scheme: 'HLLC'	        	  		# Riemann solver scheme

   # Runge-Kutta time scheme: RK1, SSP-RK2, SSP-RK3

   runge_kutta_time_scheme: 'SSP-RK3'		          		# Runge-Kutta time scheme

src/FlowSolverRHEA.cpp

+101 −0
Original line number Diff line number Diff line
@@ -74,6 +74,8 @@ FlowSolverRHEA::FlowSolverRHEA(const string name_configuration_file) : configura 
        riemann_solver = new HllcPlusApproximateRiemannSolver();

    } else if( riemann_solver_scheme == "ECKEP" ) {

        riemann_solver = new EckepFluxApproximateRiemannSolver();

    } else if( riemann_solver_scheme == "ECKEP-MOVERS-RH" ) {

        riemann_solver = new EckepMoversRHFluxApproximateRiemannSolver();	

    } else {

        cout << "Riemann solver not available!" << endl;

        MPI_Abort( MPI_COMM_WORLD, 1 );
@@ -4475,6 +4477,105 @@ double EckepFluxApproximateRiemannSolver::calculateIntercellFlux(const double &r 
};





////////// EckepMoversRHFluxApproximateRiemannSolver CLASS //////////



EckepMoversRHFluxApproximateRiemannSolver::EckepMoversRHFluxApproximateRiemannSolver() : BaseRiemannSolver() {};



EckepMoversRHFluxApproximateRiemannSolver::~EckepMoversRHFluxApproximateRiemannSolver() {};



double EckepMoversRHFluxApproximateRiemannSolver::calculateIntercellFlux(const double &rho_L, const double &rho_R, const double &u_L, const double &u_R, const double &v_L, const double &v_R, const double &w_L, const double &w_R, const double &E_L, const double &E_R, const double &s_L, const double &s_R, const double &P_L, const double &P_R, const double &T_L, const double &T_R, const double &a_L, const double &a_R, const int &var_type) {



    /// Method of Optimal Viscosity for Enhanced Resolution of Shocks (MOVERS) scheme:

    /// S. Jaisankar, S.V. Raghurama Rao.

    /// A central Rankine-Hugoniot solver for hyperbolic conservation laws.

    /// Journal of Computational Physics, 228, 770-798, 2009.



    /// Calculate eigenvalues

    double u = ( 1.0/2.0 )*( u_L + u_R );

    double a = ( 1.0/2.0 )*( a_L + a_R );

    double lambda_1 = abs( u - a );

    double lambda_2 = abs( u );

    double lambda_3 = abs( u + a );

    double lambda_min = min( lambda_1, min( lambda_2, lambda_3 ) );

    double lambda_max = max( lambda_1, max( lambda_2, lambda_3 ) );



    /// Conserved variable increment

    double deltaU_1 = rho_R - rho_L;          

    double deltaU_2 = rho_R*u_R - rho_L*u_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_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

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

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

    double theta = 0.1;									// !! theta needs to be larger than 0.0 !!

    alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta );  			// Prevent sonic glitch



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

    if( var_type == 0 ) {

        F *= 1.0 + 1.0;

        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;

    } else if ( var_type == 2 ) {

        F *= v_L + v_R;

    } else if ( var_type == 3 ) {

        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 );

        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 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 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 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  = ( bar_F_1*deltaV_1 + bar_F_2*deltaV_2 + bar_F_3*deltaV_3 - deltaPsi )/( deltaV_3*deltaV_3 + epsilon );

        double alpha_3  = ( 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;

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

    }



    return( F );



};





////////// BaseExplicitRungeKuttaMethod CLASS //////////



BaseExplicitRungeKuttaMethod::BaseExplicitRungeKuttaMethod() {};

src/FlowSolverRHEA.hpp

+28 −0
Original line number Diff line number Diff line
@@ -36,6 +36,7 @@ class HllApproximateRiemannSolver; /// HLL approximate Riemann solver 
class HllcApproximateRiemannSolver;			/// HLLC approximate Riemann solver

class HllcPlusApproximateRiemannSolver;			/// HLLC+ approximate Riemann solver

class EckepFluxApproximateRiemannSolver;		/// ECKEP scheme approximate Riemann solver

class EckepMoversRHFluxApproximateRiemannSolver;	/// ECKEP MOVERS Rankine-Hugoniot scheme approximate Riemann solver



class BaseExplicitRungeKuttaMethod;			/// Base explicit Runge-Kutta method

class RungeKutta1Method;				/// Runge-Kutta 1 (RK1) method
@@ -732,6 +733,33 @@ class EckepFluxApproximateRiemannSolver : public BaseRiemannSolver { 


};



////////// EckepMoversRHFluxApproximateRiemannSolver CLASS //////////

class EckepMoversRHFluxApproximateRiemannSolver : public BaseRiemannSolver {

   

    public:



        ////////// CONSTRUCTORS & DESTRUCTOR //////////

        EckepMoversRHFluxApproximateRiemannSolver();						/// Default constructor

        virtual ~EckepMoversRHFluxApproximateRiemannSolver();					/// Destructor



	////////// GET FUNCTIONS //////////



	////////// SET FUNCTIONS //////////



	////////// METHODS //////////

       

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

	#pragma acc routine

        double calculateIntercellFlux(const double &rho_L, const double &rho_R, const double &u_L, const double &u_R, const double &v_L, const double &v_R, const double &w_L, const double &w_R, const double &E_L, const double &E_R, const double &s_L, const double &s_R, const double &P_L, const double &P_R, const double &T_L, const double &T_R, const double &a_L, const double &a_R, const int &var_type);



    protected:



        ////////// PARAMETERS //////////



    private:



};



////////// BaseExplicitRungeKuttaMethod CLASS //////////

class BaseExplicitRungeKuttaMethod {

stuff/rhea_flow_solver.py

+88 −46
Original line number Diff line number Diff line
@@ -977,62 +977,103 @@ 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_ 




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

#@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, 2025.

#    

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

#    if( var_type == 0 ):

#        F *= 1.0

#    elif ( var_type == 1 ):

#        F *= ( 1.0/2.0 )*( u_L + u_R )

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

#    elif ( var_type == 2 ):

#        F *= ( 1.0/2.0 )*( v_L + v_R )

#    elif ( var_type == 3 ):

#        F *= ( 1.0/2.0 )*( w_L + w_R )

#    elif ( var_type == 4 ):

#        bar_F_1 = ( 1.0/2.0 )*( rho_L*u_L + rho_R*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 )*( rho_L*u_L*E_L + P_L*u_L + rho_R*u_R*E_R + P_R*u_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 )

#        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

#    

#    return( F )

@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 ):

    

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

    if( var_type == 0 ):

        F *= 1.0 + 1.0

    elif ( var_type == 1 ):

        F *= u_L + u_R

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

    elif ( var_type == 2 ):

        F *= v_L + v_R

    elif ( var_type == 3 ):

        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 )

        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 )

        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

    

    return( F )



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



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

@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 ):

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

  

    # Method of Optimal Viscosity for Enhanced Resolution of Shocks (MOVERS) scheme:

    # S. Jaisankar, S.V. Raghurama Rao.

    # A central Rankine-Hugoniot solver for hyperbolic conservation laws.

    # Journal of Computational Physics, 228, 770-798, 2009.    



    # Calculate eigenvalues

    u = ( 1.0/2.0 )*( u_L + u_R )

    a = ( 1.0/2.0 )*( a_L + a_R )

    lambda_1 = abs( u - a )

    lambda_2 = abs( u )

    lambda_3 = abs( u + a )

    lambda_min = min( lambda_1, lambda_2, lambda_3 )

    lambda_max = max( lambda_1, lambda_2, lambda_3 )



    # Conserved variable increment

    deltaU_1 = rho_R - rho_L            

    deltaU_2 = rho_R*u_R - rho_L*u_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_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

    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 )

    theta = 0.1                                                 # !! theta needs to be larger than 0.0 !!

    alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta )   # Prevent sonic glitch



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

    if( var_type == 0 ):

        F *= 1.0 + 1.0

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

    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

    elif ( var_type == 2 ):

        F *= v_L + v_R

    elif ( var_type == 3 ):
@@ -1061,6 +1102,7 @@ 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, 
        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 )