Commit 269fe9df authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

New convection scheme modified II

parent f08a3197

src/FlowSolverRHEA.cpp

+10 −11
Original line number Diff line number Diff line
@@ -4539,24 +4539,23 @@ double EckepMoversRHFluxApproximateRiemannSolver::calculateIntercellFlux(const d 


    /// Pressure sensor: shocks have sharp pressure jump

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

    //double rel_P_jump = abs( P_R - P_L )/( P + epsilon );

    double rel_P_jump = abs( P_R - P_L )/( P + 1.0e-10 );				// ... modified for OpenACC

    //double sensor_P = abs( P_R - P_L )/( abs( P_L - P ) + abs( P_R - P ) + epsilon );

    double sensor_P = abs( P_R - P_L )/( abs( P_L - P ) + abs( P_R - P ) + 1.0e-10 );				// ... modified for OpenACC



    /// Apply sensor: if delta_u > 0 (rarefaction/expansion) or rel_P_jump < 5% --> alpha_s disabled

    if( delta_u > 0.0 || rel_P_jump < 0.05 ) alpha_S = 0.0;

    /// Apply sensor: if delta_u > 0 (rarefaction/expansion)

    if( delta_u > 0.0 ) alpha_S = 0.0;

    alpha_S *= sensor_P;



    /// -------------------------------///

    /// END: SHOCK SENSOR MODIFICATION ///

    /// -------------------------------///



    /// Prevent sonic/entropy glitch 

    ////if( alpha_S > epsilon ) {

    //if( alpha_S > 1.0e-10 ) {								// ... modified for OpenACC

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

    //    alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta ); 

    //}

    double delta = 0.1*lambda_max;

    if( abs( alpha_S ) < delta ) alpha_S = 0.5*( alpha_S*alpha_S/delta + delta );

    //if( alpha_S > epsilon ) {

    if( alpha_S > 1.0e-10 ) {						// ... modified for OpenACC

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

        alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta );

    }



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

    if( var_type == 0 ) {

stuff/rhea_flow_solver.py

+7 −9
Original line number Diff line number Diff line
@@ -1073,23 +1073,21 @@ 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, 


    # Pressure sensor: shocks have sharp pressure jump

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

    rel_P_jump = abs( P_R - P_L )/( P_avg + epsilon )

    sensor_P = abs( P_R - P_L )/( abs( P_L - P ) + abs( P_R - P ) + epsilon )



    # Apply sensor: if delta_u > 0 (rarefaction/expansion) or rel_P_jump < 5% --> alpha_s disabled

    if( delta_u > 0.0 or rel_P_jump < 0.05 ):

    # Apply sensor: if delta_u > 0 (rarefaction/expansion)

    if( delta_u > 0.0 ):

        alpha_S = 0.0

    alpha_S *= sensor_P;



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

    ### END: SHOCK SENSOR MODIFICATION ###

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



    # Prevent sonic/entropy glitch 

    #if( alpha_S > epsilon ):

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

    #    alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta )

    delta = 0.1*lambda_max

    if( abs( alpha_S ) < delta ):

       alpha_S = 0.5*( alpha_S*alpha_S/delta + delta )

    if( alpha_S > epsilon ):

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

        alpha_S = ( alpha_S*alpha_S + theta*theta )/( 2.0*theta )



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

    if( var_type == 0 ):

tests/3d_turbulent_boundary_layer/myRHEA.cpp

+4 −4
Original line number Diff line number Diff line
@@ -29,12 +29,12 @@ const double tau = delta_0_ast/U_infty; /// Characteristic time scale 
const double L_x        = 750*delta_0_ast;			/// Domain size in x-direction 

//const double L_y        = 40*delta_0_ast;			/// Domain size in y-direction

//const double L_z        = 34*delta_0_ast;			/// Domain size in z-direction

const double alpha_u    = 1.0e-1;                               /// Magnitude of velocity perturbations

const double alpha_u    = 1.0;                                  /// Magnitude of velocity perturbations

const double alpha_P    = 1.0e-3;                               /// Magnitude of pressure perturbations

const double x_0_domain = 0.0;					/// Domain origin in x-direction

const double y_0_domain = 0.0;					/// Domain origin in y-direction

const double L_x_forcing = 10*delta_0_ast;			/// Location of forcing downstream of inlet

const int time_iter_forcing = 100.0;				/// Forcing iteration interval

const int time_iter_forcing = 10.0;				/// Forcing iteration interval

const double delta_x    = L_x/1024.0;				/// Mesh size in x-direction

double random_number;						/// Random number
@@ -93,9 +93,9 @@ void myRHEA::calculateSourceTerms() { 
                f_rhou_field[I1D(i,j,k)] = 0.0;

                f_rhov_field[I1D(i,j,k)] = 0.0;

		x_position = x_field[I1D(i,j,k)];

		if( ( x_position > ( x_0_domain + L_x_forcing - delta_x ) ) && ( x_position < ( x_0_domain + L_x_forcing + delta_x ) ) ) {

		if( ( x_position > ( x_0_domain + L_x_forcing - 2.0*delta_x ) ) && ( x_position < ( x_0_domain + L_x_forcing + 2.0*delta_x ) ) ) {

		    y_position = y_field[I1D(i,j,k)];

	            if( y_position < 0.3*delta_0_ast ) {

	            if( y_position < 1.0*delta_0_ast ) {

                        f_rhov_field[I1D(i,j,k)] = rho_ref*( alpha_u*random_number*U_infty )/tau;

		    }

		}