Are my implicit type conversions correct?

Question

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

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.