Commit 66e03346 authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

CBC improved

parent f68f2978

src/FlowSolverRHEA.cpp

+166 −58
Original line number Diff line number Diff line
@@ -1074,22 +1074,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_WEST_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_WEST_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->x[i+1] - mesh->x[i];
@@ -1237,22 +1255,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_EAST_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_EAST_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->x[i-1] - mesh->x[i];
@@ -1400,22 +1436,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_SOUTH_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

		            double drho_dy_g = ( rho_in - rho_g )/Delta_g;

		            double dv_dy_g   = (   v_in -   v_g )/Delta_g;

		            double dP_dy_g   = (   P_in -   P_g )/Delta_g;

                            double L_2_lambda_2_g = sos_in*sos_in*drho_dy_g - dP_dy_g;

                            double L_5_lambda_5_g = dP_dy_g + rho_in*sos_in*dv_dy_g;

		            double dQ_1_dy_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                            rho_g = rho_in - Delta_g*dQ_1_dy_in;

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_SOUTH_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->y[j+1] - mesh->y[j];
@@ -1563,22 +1617,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_NORTH_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

		            double drho_dy_g = ( rho_in - rho_g )/Delta_g;

		            double dv_dy_g   = (   v_in -   v_g )/Delta_g;

		            double dP_dy_g   = (   P_in -   P_g )/Delta_g;

                            double L_2_lambda_2_g = sos_in*sos_in*drho_dy_g - dP_dy_g;

                            double L_5_lambda_5_g = dP_dy_g + rho_in*sos_in*dv_dy_g;

		            double dQ_1_dy_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                            rho_g = rho_in - Delta_g*dQ_1_dy_in;

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_NORTH_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->y[j-1] - mesh->y[j];
@@ -1726,22 +1798,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_BACK_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

		            double drho_dz_g = ( rho_in - rho_g )/Delta_g;

		            double dw_dz_g   = (   w_in -   w_g )/Delta_g;

		            double dP_dz_g   = (   P_in -   P_g )/Delta_g;

                            double L_2_lambda_2_g = sos_in*sos_in*drho_dz_g - dP_dz_g;

                            double L_5_lambda_5_g = dP_dz_g + rho_in*sos_in*dw_dz_g;

		            double dQ_1_dz_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                            rho_g = rho_in - Delta_g*dQ_1_dz_in;

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_BACK_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->z[k+1] - mesh->z[k];
@@ -1889,22 +1979,40 @@ void FlowSolverRHEA::updateBoundaries() { 
                    rho_g = rho_field[I1D(i,j,k)];

		    T_g   = ( bocos_T[_FRONT_] - wg_in*T_in )/wg_g;

                    P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                    double rel_error = 1.0;

		    /// Aitken’s delta-squared process:

                    double x_0_Aitken = rho_g;

		    #pragma acc loop seq

                    for( int ite = 0; ite < max_iter; ite++ ) {

                        if( rel_error >= rel_tol ) {

		            double drho_dz_g = ( rho_in - rho_g )/Delta_g;

		            double dw_dz_g   = (   w_in -   w_g )/Delta_g;

		            double dP_dz_g   = (   P_in -   P_g )/Delta_g;

                            double L_2_lambda_2_g = sos_in*sos_in*drho_dz_g - dP_dz_g;

                            double L_5_lambda_5_g = dP_dz_g + rho_in*sos_in*dw_dz_g;

		            double dQ_1_dz_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

		            double rho_g_old = rho_g;

                            rho_g = rho_in - Delta_g*dQ_1_dz_in;

                        /// Aitken's x_1

                        double drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        double du_dx_g   = (   u_in -   u_g )/Delta_g;

                        double dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        double L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        double L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        double dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_1_Aitken = rho_g;

                        /// Aitken's x_2

                        drho_dx_g = ( rho_in - rho_g )/Delta_g;

                        du_dx_g   = (   u_in -   u_g )/Delta_g;

                        dP_dx_g   = (   P_in -   P_g )/Delta_g;

                        L_2_lambda_2_g = sos_in*sos_in*drho_dx_g - dP_dx_g;

                        L_5_lambda_5_g = dP_dx_g + rho_in*sos_in*du_dx_g;

                        dQ_1_dx_in = ( 1.0/( sos_in*sos_in ) )*( L_2_lambda_2_g + 0.5*( L_5_lambda_5_g + L_1_lambda_1_in_in ) );

                        rho_g = rho_in - Delta_g*dQ_1_dx_in;

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

                        double x_2_Aitken = rho_g;

                        /// Aitken's iteration

                        double denominator = x_2_Aitken - 2.0*x_1_Aitken + x_0_Aitken;

                        rho_g = x_2_Aitken - ( pow( x_2_Aitken - x_1_Aitken, 2.0 )/( denominator + epsilon ) );

                        P_g   = thermodynamics->calculatePressureFromTemperatureDensity( T_g, rho_g );

		            rel_error = abs( ( rho_g - rho_g_old )/rho_g_old );

		            //cout << ite << "  " << rho_g_old << "  " << rho_g << "  " << rel_error << endl;

                        ///cout << ite << "  " << rho_g << "  " << x_0_Aitken << "  " << x_1_Aitken << "  " << x_2_Aitken << endl;

                        /// Aitken's convergence

                        if( abs( (rho_g - x_2_Aitken )/rho_g ) < rel_tol ) {

                            break;		/// If the result is within tolerance, leave the loop!

			}

			x_0_Aitken = rho_g;	/// Otherwise, update x_0 to iterate again ...

		    }

		} else if( bocos_type[_FRONT_] == _SUBSONIC_OUTFLOW_ ) {

                    double Delta_g  = mesh->z[k-1] - mesh->z[k];