Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

I'm having some issues with my code and I'm trying to rule out any syntax or programming mistakes. One of the things I'm trying to rule out is the different variable types used in my equations. I'm mostly operating with doubles, but I did try to ensure that all the integers in the equations were correct. Thank you in advance.

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <cmath>
#include <algorithm>

#define pi 3.14159265

using namespace std;

//Define Fluid Properties
double rho_L = 998; //Liquid Density
double rho_LG = 828.9; //Liquid-Gas Density Ratio
double mu_L = 0.000798; //Liquid Viscosity
double mu_LG = 40.24; //Liquid-Gas Viscosity Ratio
double sigma = 0.0712; //Surface Tension
double nu_G = (mu_L/mu_LG)/(rho_L/rho_LG);

//Define Injector Properties
double Uinj = 56.7; //Injection Velocity
double Dinj = 0.0998; //Injector Diameter
double theta = 45.0*pi/180.0; //Spray Cone Angle
double L = 100.0*Dinj; //Atomization Length
double Ls = L/Dinj; //Normalized Atomization Length

//Define Solver Parameters
double K = 26470; //Viscous Dissipation Coefficient
double Eps = pow(10,-5); //Residual Error
double dx = 0.01; //Step Size
double Ui = 15.0; //Initial Guess
//const int Z = static_cast<const int>((Ls-1)/dx + 1) + 1;
const int N = 9901;//Z;

double deriv (double U, double X, double delta, double m)
{
    double dudx;
    dudx = -(1.0/delta)*(1.0/U)*(U - sqrt(1.0 - U)/sqrt(m*X*X))*(U - sqrt(1.0 - U)/sqrt(m*X*X));
    return (dudx);
}

int main()
{
    //Declare Variables
    int max_step;
    double ERR;
    int step;
    double DEN;
    double SMD;
    double m;
    double Ug;
    double Re;
    double Cd;
    double delta;
    double K1;
    double K2;
    double K3;
    double K4;

    //Allocate Memory From Heap
    double *U = new double [N];
    double *X = new double [N];

    //Initialize Vectors and Variables
    DEN = 0.5*rho_L - (4.0/3.0)*K*(mu_L)/(Uinj*Dinj*Dinj)*L;
    m = 4.0/rho_LG*tan(theta)*tan(theta);

    for (int i = 0; i < N; i++)
    {
        X[i] = 1.0 + dx*i;
    }
    U[0] = 1.0;

    max_step = 10000;
    ERR = 1;
    step = 0;
    while(abs(ERR) > Eps && step < max_step)
    {

        //Calculate Ug
        Ug = sqrt(1.0 - (Ui/Uinj))/sqrt(m*Ls*Ls)*Uinj;

        //Calculate SMD
        SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

        //Calculate Re # and Drag Coefficient
        Re = abs(Ui-Ug)*SMD/nu_G;

        if(Re <= 0.01)
        {
            Cd = (0.1875) + (24.0/Re);
        }
        else if(Re > 0.01 && Re <= 260.0)
        {
            Cd = (24.0/Re)*(1.0 + 0.1315*pow(Re,(0.32 - 0.05*log10(Re))));
        }
        else
        {
            Cd = (24.0/Re)*(1.0 + 0.1935*pow(Re,0.6305));
        }

        //Determine New U
        delta = (4.0/3.0)*(1.0/Cd)*(rho_LG)*(SMD/Dinj);

        //RK4
        for (int i = 0; i < N-1; i++)
        {
            K1 = deriv(U[i],X[i],delta,m);
            K2 = deriv(U[i]+0.5*dx*K1,X[i]+0.5*dx,delta,m);
            K3 = deriv(U[i]+0.5*dx*K2,X[i]+0.5*dx,delta,m);
            K4 = deriv(U[i]+dx*K3,X[i+1],delta,m);
            U[i+1] = U[i] + dx/6.0*(K1 + 2.0*K2 + 2.0*K3 + K4);
        }

        ERR = abs(U[N-1]*Uinj - Ui)/Ui;

        Ui = U[N-1]*Uinj;

        step = step + 1;
    }

    SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

    /*cout << "U = " << Ui << endl;
    cout << "SMD = " << SMD << endl;
    cout << "DEN = " << DEN << endl;
    cout << "Ug = " << Ug << endl;
    cout << "m = " << m << endl;
    cout << "delta = " << delta << endl;
    cout << "Re = " << Re << endl;
    cout << "Cd = " << Cd << endl;
    cout << "U* = " << U[N-1] << endl;
    cout << "Error = " << ERR << endl;
    cout << "step = " << step << endl;*/

    //Output Data Into Text File
    ofstream outputdata("result.txt");
    for (int i = 0; i < N; i++)
    {
        outputdata << X[i] << " " << U[i] << '\n';
    }
    outputdata.close();

    delete [] U;
    delete [] X;

    return 0;
}
share|improve this question
    
This is likely off-topic as per the FAQ - 4. To the best of my knowledge, does the code work? –  Yuushi Feb 13 '13 at 1:43
    
The code is functioning, yes. The results are slightly off of what I expected, so I was hoping to localize the inaccuracies. –  user1562781 Feb 14 '13 at 3:41
add comment

1 Answer 1

There are no implicit type conversions here: the only variables that are not double are the step and max_step and they are not used with the other ones. Otherwise, you're using doubles everywhere so there is no issue.

Slight inaccuracies are possible when doing a lot of iterations, the most likely reason is that floating point numbers are simply not fully accurate. If you're seeing a difference when comparing this to another program, it's possible that the order of operations is the issue, too. On a computer, a * (b * c) is slightly different from (a * b) * c.

share|improve this answer
    
You probably meant a * (b * c) is different from (a * b) * c . –  Sylvain Defresne Feb 16 '13 at 20:54
    
Indeed. Fixed it. –  Quentin Pradet Feb 16 '13 at 21:06
    
What about all the different integers used in the equation (like the (4/3) factor)? I found that I get different results if I use (4/3) vs (4.0/3.0). –  user1562781 Feb 18 '13 at 1:45
    
@user1562781 Yes, indeed, 4/3 == 1, while 4.0/3.0 is very close to 1.25, which is what you want. Your question only contains doubles, either because the variable is double (double K = 25670), or because it is written as a double (4.0). –  Quentin Pradet Feb 18 '13 at 6:43
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.