hips.cc

Go to the documentation of this file.

 
 
#include "hips.h"
 
#ifdef REACTIONS_ENABLED
#include "batchReactor_cvode.h"
#include "batchReactor_cantera.h"
#endif
 
#include "randomGenerator.h"
 
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <cmath>
#include <vector>
 
using namespace std;
 
hips::hips(int nLevels_, 
           double domainLength_, 
           double tau0_, 
           double C_param_, 
           int forceTurb_,
           int nVar_,
           vector<double> &ScHips_,
           bool performReaction_,
           shared_ptr<void> vcantSol,
           int seed,
           int realization_): 
    nLevels(nLevels_), 
    domainLength(domainLength_), 
    tau0(tau0_),
    C_param(C_param_), 
    forceTurb(forceTurb_),       
    ScHips(ScHips_),   
    nVar(nVar_),                       
    LrandSet(true),              
    rand(seed),
    performReaction(performReaction_),
    realization(realization_){
 
    #ifdef REACTIONS_ENABLED
    if(performReaction) {
        shared_ptr<Cantera::Solution> cantSol = static_pointer_cast<Cantera::Solution>(vcantSol);
        gas = cantSol->thermo(); 
        nsp = gas->nSpecies();
 
        bRxr = make_shared<batchReactor_cvode>(cantSol);                                // By default, use batchReactor_cvode
 
        // Uncomment the following line to switch to batchReactor_cantera
        // bRxr = make_unique<batchReactor_cantera>(cantSol);
    }
    #endif
 
    // Resize vectors to the number of variables
    varData.resize(nVar);
    varName.resize(nVar); 
 
    set_tree(nLevels, domainLength, tau0, ScHips);
}
 
 
hips::hips(double C_param_, 
           int forceTurb_,
           int nVar_,
           bool performReaction_,
           shared_ptr<void> vcantSol,
           int seed,
           int realization_): 
    C_param(C_param_), 
    forceTurb(forceTurb_),       
    nVar(nVar_),                       
    LrandSet(true),              
    rand(seed),
    performReaction(performReaction_),
    realization(realization_){
 
    #ifdef REACTIONS_ENABLED
    // Initialize Cantera thermo phase and species count.
    if(performReaction) {
        shared_ptr<Cantera::Solution> cantSol = static_pointer_cast<Cantera::Solution>(vcantSol);
        gas = cantSol->thermo(); 
        nsp = gas->nSpecies();
 
        // Set up the default batch reactor (cvode).
       // bRxr = make_shared<batchReactor_cvode>(cantSol);
 
        // Uncomment the following line to switch to batchReactor_cantera.
         bRxr = make_shared<batchReactor_cantera>(cantSol);
    }
    #endif
 
    // Resize vectors to accommodate the number of variables.
    varData.resize(nVar);
    varName.resize(nVar);        
}
 
 
void hips::set_tree(int nLevels_, double domainLength_, double tau0_, vector<double> &ScHips_){    
 
    nLevels= nLevels_; 
    domainLength = domainLength_; 
    tau0 = tau0_; 
    ScHips = ScHips_; 
 
    if (nLevels == -1)  
        nLevels = nL; 
 
    iEta = nLevels - 3;                                // Kolmogorov level; if nLevels = 7, then 0, 1, 2, 3, (4), 5, 6; iEta=4 is the lowest swap level: swap grandchildren of iEta=4 at level 6.
           
    int maxSc = 1.0;
 
    for (int i=0; i<ScHips.size(); i++)
        maxSc = ScHips[i]>maxSc ? ScHips[i] : maxSc;
    
    if (maxSc > 1.0)
        nLevels += ceil(log(maxSc)/log(4));            // Changing number of levels!
    
    Nm1 = nLevels - 1;
    Nm2 = nLevels - 2;
    Nm3 = nLevels - 3;
    
    // -------------------------- 
    
    nparcels = static_cast<int>(pow(2, Nm1));
    parcelTimes.resize(nparcels,0);
    i_batchelor.resize(nVar,0);
    
    vector<double> levelLengths(nLevels);              // Including all levels, but last 2 don't count:
    vector<double> levelTaus(nLevels);                 // Smallest scale is 2 levels up from bottom
    levelRates   = vector<double>(nLevels);
 
    for (int i=0; i<nLevels; i++) {
        levelLengths[i] = domainLength * pow(Afac,i);
       //levelLengths[i] = domainLength * pow(Anew,i);
 
        levelTaus[i] = tau0 * pow(levelLengths[i]/domainLength, 2.0/3.0) / C_param;
        levelRates[i] = 1.0/levelTaus[i] * pow(2.0,i);
    }
 
    LScHips = ScHips.size() > 0 ? true : false;
    if (LScHips) {                                     // Ccorrect levels for high Sc (levels > Kolmogorov)
        for (int i=iEta+1; i<nLevels; i++) {
            levelTaus[i] = tau0 *
            pow(levelLengths[iEta]/domainLength, 2.0/3.0) /
            C_param;
            levelRates[i] = 1.0/levelTaus[i] * pow(2.0,i);
        }
    }
 
    //-------------------------------------------------
 
    eddyRate_total = 0.0;
    for (int i=0; i<=Nm3; i++)
        eddyRate_total += levelRates[i];
    
    eddyRate_inertial = 0.0;
    for (int i=0; i<=iEta; i++)
        eddyRate_inertial += levelRates[i];
    
    //-------------------
    
    i_plus.resize(nVar);
    
    for (int k=0; k<nVar; k++) {
        if (ScHips[k] < 1.0)
            i_batchelor[k] = iEta + 1.5*log(ScHips[k])/log(4);
        else if (ScHips[k] > 1.0)
            i_batchelor[k] = iEta + log(ScHips[k])/log(4);
        else
            i_batchelor[k] = iEta;
 
        i_plus[k] = ceil(i_batchelor[k]);
    }
    
    //------------------- Set the parcel addresses (index array)
 
    varRho.resize(nparcels);
    Temp.resize(nparcels);
 
    pLoc.resize(nparcels);
    for (int i=0; i<nparcels; i++)
        pLoc[i] = i;
} 
 
 
void hips::set_tree(double Re_, double domainLength_, double tau0_, std::vector<double> &ScHips_, std::string approach_) {
    Re = Re_;
    domainLength = domainLength_;
    tau0 = tau0_;
    ScHips = ScHips_;
    approach = approach_;
 
    double baseLevelEstimate = (3.0 / 4) * log(1 / Re) / log(Afac);                               // Calculate the base tree level estimate (non-integer)
    int baseLevel;
 
    if (approach == "rounding") {
        baseLevel = round(baseLevelEstimate);                                                     // Round the base level to the nearest integer
    } 
    else if (approach == "probability") {
        baseLevel = ceil(baseLevelEstimate);                                                      // Ceil the base level to the nearest integer
        int previousLevel = baseLevel - 1;
        Prob = baseLevelEstimate - previousLevel;                                                 // Calculate the probability
    } 
    else if (approach == "micromixing") {
        baseLevel = ceil(baseLevelEstimate);                                                      // Ceil the base level to the nearest integer
        lStar = std::pow(Re, -3.0 / 4);                                                           // Calculate lStar based on Re
    } 
    else if (approach == "dynamic_A") {
        baseLevel = round(baseLevelEstimate);                                                     // Round the base level to the nearest integer
        Anew = exp(-log(Re) / ((4.0 / 3.0) * baseLevel));                                         // Calculate the new value of parameter A
    } 
    else {
        throw std::invalid_argument("Invalid approach specified");                                // Handle invalid approach case if needed
    }
 
    nL = baseLevel + 3;                                                                           // Set the number of levels for the binary tree structure
 
    //----------------------------------------------------------------------
    nLevels = nL;
    iEta = nLevels - 3;  // Kolmogorov level 
 
    int maxSc = 1;
    for (const auto &sc : ScHips)
        maxSc = std::max(maxSc, static_cast<int>(sc));
    if (maxSc > 1.0)
        nLevels += ceil(log(maxSc) / log(4));                           
 
    Nm1 = nLevels - 1;
    Nm2 = nLevels - 2;
    Nm3 = nLevels - 3;
 
    nparcels = static_cast<int>(pow(2, Nm1));
    parcelTimes.resize(nparcels, 0);
    i_batchelor.resize(nVar, 0);
 
    std::vector<double> levelLengths(nLevels);                                // Including all levels, but last 2 don't count
    std::vector<double> levelTaus(nLevels);                                   // Smallest scale is 2 levels up from bottom
    levelRates.resize(nLevels);
 
    for (int i = 0; i < nLevels; ++i) {
        levelLengths[i] = domainLength * pow((approach == "dynamic_A" ? Anew : Afac), i);
 
        levelTaus[i] = tau0 * pow(levelLengths[i] / domainLength, 2.0 / 3.0) / C_param;
        levelRates[i] = 1.0 / levelTaus[i] * pow(2.0, i);
    }
 
    if (approach == "micromixing") {                                          // Adjust rates for micromixing model
        levelTaus[Nm3] = tau0 * pow(lStar / domainLength, 2.0 / 3.0) / C_param;
        levelRates[Nm3] = 1.0 / levelTaus[Nm3] * pow(2.0, Nm3);
    }
 
    if (approach == "probability") {                                          // Adjust final mixing rate based on probability
        levelRates[Nm3] = levelRates[nL - 3] * Prob;
    }
 
    LScHips = !ScHips.empty();                                               // Correct levels for high Sc (levels > Kolmogorov)
    if (LScHips) {
        for (int i = iEta + 1; i < nLevels; ++i) {
            levelTaus[i] = tau0 * pow(levelLengths[iEta] / domainLength, 2.0 / 3.0) / C_param;
            levelRates[i] = 1.0 / levelTaus[i] * pow(2.0, i);
        }
    }
 
    //-----------------------------------------------------
 
    eddyRate_total = 0.0;
    for (int i = 0; i <= Nm3; ++i)
        eddyRate_total += levelRates[i];
 
    eddyRate_inertial = 0.0;
    for (int i = 0; i <= iEta; ++i)
        eddyRate_inertial += levelRates[i];
 
    //-----------------------------------------------------
 
    i_plus.resize(nVar);
    for (int k = 0; k < nVar; ++k) {
        if (ScHips[k] < 1.0)
            i_batchelor[k] = iEta + 1.5 * log(ScHips[k]) / log(4);
        else if (ScHips[k] > 1.0)
            i_batchelor[k] = iEta + log(ScHips[k]) / log(4);
        else
            i_batchelor[k] = iEta;
        i_plus[k] = ceil(i_batchelor[k]);
    }
 
    //-----------------------------------------------------
 
    varRho.resize(nparcels);
    Temp.resize(nparcels);
    pLoc.resize(nparcels);
    for (int i = 0; i < nparcels; ++i)
        pLoc[i] = i;
}
 
 
void hips::set_varData(std::vector<double> &v, std::vector<double> &w, const std::string &varN) {  
    
    varData[currentIndex] = std::make_shared<std::vector<double>>(projection(v, w));
    varName[currentIndex] = varN;
 
    currentIndex++; 
}
 
 
void hips::set_varData(std::vector<double> &v, std::vector<double> &w, const std::string &varN, const std::vector<double> &rho) {
    
    std::pair<std::vector<double>, std::vector<double>> results = projection(v, w, rho);
    std::vector<double> vh = results.first;
    std::vector<double> rho_h = results.second;
 
    varData[currentIndex] = std::make_shared<std::vector<double>>(projection(v, w));            
    varRho = std::vector<double>(rho_h);
    varName[currentIndex] = varN;
 
    currentIndex++; 
}
 
 
std::vector<double> hips::projection(std::vector<double> &vcfd, std::vector<double> &weight) {
    
    xc = setGridCfd(weight);                               // Populate the physical domain for flow particles
    xh = setGridHips(nparcels);                            // Populate the physical domain for hips parcels
 
    int nc = xc.size() - 1;                 
    int nh = xh.size() - 1; 
 
    std::vector<double> vh(nh, 0.0);
    int jprev = 0;
 
    for(int i = 0; i < nh; i++) {
        for(int j = jprev + 1; j <= nc; ++j) {
            if(xc[j] <= xh[i + 1]) {
                double d1 = xc[j] - xc[j - 1];
                double d2 = xc[j] - xh[i];
                // calculation of shortest distance
                // handling if distance is zero but substarction gives comutational error
                double d  = std::min(d1, d2) < 1E-15 ? 0 : std::min(d1, d2);
 
                vh[i] += vcfd[j - 1] * d;
            } 
            else {
                double d1 = xh[i + 1] - xc[j - 1];
                double d2 = xh[i + 1] - xh[i];
                // calculation of shortest distance
                // handling if distance is zero but substarction gives comutational error
                double d  = std::min(d1, d2) < 1E-15 ? 0 : std::min(d1, d2);
 
                vh[i] += vcfd[j - 1] * d;
                jprev = j - 1;
                break;
            }
        }
        vh[i] /= (xh[i + 1] - xh[i]);   
    }
    return vh;
}
 
 
std::pair<std::vector<double>, std::vector<double>> hips::projection(std::vector<double> &vcfd, 
                                                                     std::vector<double> &weight, 
                                                                     const std::vector<double> &density) {
    
    xc = setGridCfd(weight);
    xh = setGridHips(nparcels);
 
    // Initialize variables
    int nc = xc.size() - 1;
    int nh = xh.size() - 1;
 
    std::vector<double> vh(nh, 0.0);
    std::vector<double> rho_h(nh, 0.0);
    int jprev = 0;
 
    for (int i = 0; i < nh; i++) {
        double total_dx = 0.0;
 
        for (int j = jprev + 1; j <= nc; ++j) {
            if (xc[j] <= xh[i + 1]) {
                double d = std::min(xc[j] - xc[j - 1], xc[j] - xh[i]);
                total_dx += d;
                rho_h[i] += density[j - 1] * d;
                vh[i] += density[j - 1] * vcfd[j - 1] * d;
            } 
            else {
                double d = std::min(xh[i + 1] - xc[j - 1], xh[i + 1] - xh[i]);
                total_dx += d;
                rho_h[i] += density[j - 1] * d;
                vh[i] += density[j - 1] * vcfd[j - 1] * d;
                jprev = j - 1;
                break;
            }
        }
 
        // ------------------------ Normalize results
 
        rho_h[i] /= (xh[i + 1] - xh[i]);
        vh[i] /= rho_h[i] * (xh[i + 1] - xh[i]);
    }
    return {vh, rho_h};
}
 
 
std::vector<double> hips::setGridCfd(std::vector<double> &w) {
    
    std::vector<double> pos;                               // Initializing a vector to hold the grid positions
    double posL = 0.0;                                     // Initializing the starting position
 
    int i = 0;
 
    while (i <= w.size()) {                               // Generate the grid positions based on the weights
        pos.push_back(posL);                              // Add the current position to the grid
        posL += w[i];                                     // Move to the next position by adding the corresponding weight
        i++;                                              
    }
    return pos;                                           // Return the generated grid positions
}
 
 
std::vector<double> hips::setGridHips(int N){
 
    std::vector<double> xh(N + 1);                               // Initialize a vector to hold the grid points
    double step = 1.0 / N;                                       // Calculate the step size
 
    for(int i = 0; i <= N; i++)                                 // Populate the grid with evenly spaced points
        xh[i] = i * step;
        
    return xh;                                                  // Return the generated grid
}
 
 
void hips::calculateSolution(const double tRun, bool shouldWriteData) {
    
    unsigned long long nEddies = 0;                               // Number of eddy events
    int fileCounter = 0;                                          // Number of data files written
    int iLevel;                                                   // Tree level of EE with top at iLevel=0
    int iTree;                                                    // One of two subtrees involved in swap at iLevel                                           
    time = 0.0;                                                   // Initialize simulation time
    int lastEddyOutput = 0;                                       // Track last eddy-based output event
 
    // Apply default values if user hasn't set them
    if (!useEddyBasedWriting && !useTimeBasedWriting) {
        outputIntervalEddy = DEFAULT_EDDY_INTERVAL;
        useEddyBasedWriting = true;  // Default to eddy-based writing
    }
 
    sample_hips_eddy(dtEE, iLevel);                               // Get first EE at time 0+dtEE
    nEddies++;
    eddyCounter = 0;                                              // Reset eddy counter at start
    lastOutputTime = 0.0;                                         // Reset last output time
 
    while (time + dtEE <= tRun) {
        time += dtEE;
        selectAndSwapTwoSubtrees(iLevel, iTree);
        advanceHips(iLevel, iTree);                              // Reaction and micromixing (if needed) to t=time
 
        sample_hips_eddy(dtEE, iLevel);
        nEddies++;
        eddyCounter++;  
 
        //  Only check the selected mode (set in example code or by default)
        bool writeByEddy = (useEddyBasedWriting && eddyCounter >= lastEddyOutput + outputIntervalEddy);
        bool writeByTime = (useTimeBasedWriting && time - lastOutputTime >= outputIntervalTime);
 
        if (shouldWriteData) {
            if (writeByEddy) {  //  Only write if using eddy-based writing
                writeData(realization, ++fileCounter, time);
                lastEddyOutput = eddyCounter;  //  Update last output event
            }
            else if (writeByTime) {  // Only write if using time-based writing
                 writeData(realization, ++fileCounter, time);
                lastOutputTime = time;
            }
        }
    }
 
    // Ensure the final time step completes
    time = tRun;
    iLevel = 0; 
    iTree = 0;
 
    if (performReaction)
        reactParcels_LevelTree(iLevel, iTree);                   // React all parcels up to end time
    saveAllParameters();
}
 
 
void hips::sample_hips_eddy(double &dtEE, int &iLevel) {
 
    static double c1 = 1.0 - pow(2.0, 5.0/3.0*(iEta+1));
    static double c2 = pow(2.0, Nm2) - pow(2.0, iEta+1);
    static double c3 = pow(2.0, iEta+1);
 
    //--------------- time to next eddy
 
    double r = rand.getRand();
    dtEE = -log(r)/eddyRate_total;
 
    //----------------- get eddy level
 
    r = rand.getRand();
 
    if ( r <= eddyRate_inertial/eddyRate_total) {     // Inertial region
        r = rand.getRand();
        iLevel = ceil(3.0/5.0*log2(1.0-r*c1) - 1.0);
        if (iLevel < 0)    iLevel = 0;
        if (iLevel > iEta) iLevel = iEta;
    }
 
    else {                                            // "Batchelor" region
        r = rand.getRand();
        iLevel = ceil(log2(r*c2 + c3) - 1.0);
        if (iLevel < iEta+1) iLevel = iEta+1;
        if (iLevel > Nm3) iLevel = Nm3;
    }
    return;
}
 
 
void hips::selectAndSwapTwoSubtrees(const int iLevel, int &iTree) {
 
    iTree = rand.getRandInt((1 << iLevel)-1);
    int zero_q = rand.getRandInt(1);                                    // 0q where q is 0 or 1
    int one_r  = 2 + rand.getRandInt(1);                                // 1r where r is 0 or 1
 
    int Qstart = (zero_q << (Nm3-iLevel)) + (iTree << (Nm1-iLevel));     // starting index of Q parcels
    int Rstart = (one_r  << (Nm3-iLevel)) + (iTree << (Nm1-iLevel));     // starting index of R parcels
    int nPswap = 1 << (Nm3-iLevel);                                      // number of parcels that will be swapped
 
    int Qend = Qstart + nPswap;                                          // inclusive indices are Qstart to Qend-1
    int Rend = Rstart + nPswap;                                          // inclusive indices are Rstart to Rend-1
    vector<int> aa(pLoc.begin()+Qstart, pLoc.begin()+Qend);
    copy(pLoc.begin()+Rstart, pLoc.begin()+Rend, pLoc.begin()+Qstart);   // python: pLoc[Qstart:Qend]=pLoc[Rstart:Rend]
    copy(aa.begin(), aa.end(), pLoc.begin()+Rstart);                     // python: pLoc[Rstart:Rend]=aa
}
 
 
void hips::advanceHips(const int iLevel, const int iTree) {
 
    if (forceTurb == 2 && iLevel == 0) {
        forceProfile();                                                  // Forcing for statistically stationary
    }
 
    bool rxnDone = false;                                               // React all variables once
    for (int k = 0; k < nVar; k++) {                                    // Upon finding first variable needing micromixing
        // Combined condition check with approach condition
        if ((iLevel >= i_plus[k]) || 
            (iLevel == i_plus[k] - 1 && rand.getRand() <= i_plus[k] - i_batchelor[k])) {
                if (!rxnDone && performReaction) {
                    reactParcels_LevelTree(iLevel, iTree);
                    rxnDone = true;
                }
                mixAcrossLevelTree(k, iLevel, iTree);
        }
    }
}
 
 
int hips::getVariableIndex(const std::string &varName) const {
 
    auto it = std::find(this->varName.begin(), this->varName.end(), varName);
    if (it == this->varName.end()) {
        throw std::runtime_error("Error: Variable name '" + varName + "' not found.");
    }
    return std::distance(this->varName.begin(), it);
}
 
 
void hips::reactParcels_LevelTree(const int iLevel, const int iTree) {
 
    #ifdef REACTIONS_ENABLED
 
    int enthalpyIdx = getVariableIndex("enthalpy");  // Dynamically find enthalpy index
    int nP = 1 << (Nm1 - iLevel);
    int istart = iTree * nP;
    int iend = istart + nP;
    int ime;
    double dt;
    double h;
    std::vector<double> y(nsp);
 
    for (int i = istart; i < iend; i++) {
        ime = pLoc[i];
        dt = time - parcelTimes[ime];
 
        // Access enthalpy using its index
        h = (*varData[enthalpyIdx])[ime];
 
          // Access species using their indices
        for (int k = 0; k < nsp; k++) {
            int speciesIdx = getVariableIndex(gas->speciesName(k));
            y[k] = (*varData[speciesIdx])[ime];
        }
        if (performReaction) {
            bRxr->react(h, y, dt);
            varRho[ime] = bRxr->getDensity();
            Temp[ime] = bRxr->temperature;
        }
        // Update enthalpy and species
        (*varData[enthalpyIdx])[ime] = h;
        for (int k = 0; k < nsp; k++) {
            int speciesIdx = getVariableIndex(gas->speciesName(k));
            (*varData[speciesIdx])[ime] = y[k];
        }
 
        parcelTimes[ime] = time;
    }
    #endif
}
 
 
void hips::mixAcrossLevelTree(int kVar, const int iLevel, const int iTree) {
    
    int istart;
    int iend;
 
    int nPmix = 1 << (nLevels - iLevel - 2);                   // Number of parcels mixed together
 
    int ime;
 
    //---------- Mix left branch of iTree
 
    istart = iTree << (Nm1-iLevel);  
    iend = istart + nPmix;
 
    double s = 0;                                               // Initialize sum to 0
    for (int i=istart; i<iend; i++) {
        ime = pLoc[i];
        s += (*varData[kVar])[ime];
    }
    for (int i=istart; i<iend; i++) {
        ime = pLoc[i];
        (*varData[kVar])[ime] = s / nPmix; 
    }
 
    //--------- Mix right branch of iTree
 
    istart = iend;
    iend = istart + nPmix;
 
    s = 0;                   // initialize sum to 0
    for (int i=istart; i<iend; i++) {
        ime = pLoc[i];
        s += (*varData[kVar])[ime];
    }
    for (int i=istart; i<iend; i++) {
        ime = pLoc[i];
        (*varData[kVar])[ime] = s / nPmix; 
    }
}
 
 
void hips::forceProfile() {
    // Loop through each variable in the HiPS profile
    for (int k = 0; k < varData.size(); k++) {
        double s;                                                    // Temporary variable for summation
 
        //---------- Force the left half of parcels to average 0 ----------
 
        for (int i = 0; i < nparcels >> 1; i++)
            s += (*varData[k])[pLoc[i]];                             // Calculate the sum of values in the left half of parcels
        
        s /= (nparcels >> 1); // Calculate the average of values in the left half of parcels
        
        for (int i = 0; i < nparcels >> 1; i++)
            (*varData[k])[pLoc[i]] += (-s - 0.0);                    // Adjust values in the left half of parcels to achieve an average of 0
 
        //---------- Force the right half of parcels to average 1 ----------
        s = 0.0;
 
        for (int i = nparcels >> 1; i < nparcels; i++)
            s += (*varData[k])[pLoc[i]];                             // Calculate the sum of values in the right half of parcels
        
        s /= (nparcels >> 1);                                        // Calculate the average of values in the right half of parcels
        
        for (int i = nparcels >> 1; i < nparcels; i++)
            (*varData[k])[pLoc[i]] += (-s + 1.0);                    // Adjust values in the right half of parcels to achieve an average of 1
    }
}
 
 
void hips::writeData(int real, const int ifile, const double outputTime) {
 
    stringstream ss1, ss2;
    string s1, s2;
 
    // Prepare directory path
    ss1 << "../data/rlz_" << setfill('0') << setw(5) << real;
    ss1 >> s1;
 
    // Create directories if they don't exist
    #ifdef _WIN32
        system(("mkdir " + s1).c_str());
    #else
        system(("mkdir -p " + s1).c_str());
    #endif
 
    ss2 << "Data_" << setfill('0') << setw(5) << ifile << ".dat";
    ss2 >> s2;
 
    string fname = s1 + "/" + s2;
    ofstream ofile(fname.c_str());
 
    // Check if file opened successfully
    if (!ofile) {
        cerr << "Error: Unable to open file " << fname << " for writing!" << endl;
        return;
    }
 
    // Write metadata (header information)
    ofile << "# time = " << outputTime << "\n";
    ofile << "# Grid Points = " << nparcels << "\n";
        
    // Write column names (include temperature if reactions are enabled)
    if(performReaction)
        ofile << setw(19) << "# Temp";  // Include temperature column if reactions are enabled
 
    for (const auto& varN : varName) {
        ofile << setw(19) << "# " << varN;
    }
    ofile << endl;  // End of the header line
 
    // Set scientific notation and precision
    ofile << scientific;
    ofile << setprecision(10);
 
    // Write data
    for (int i = 0; i < nparcels; i++) {
        // Write temperature first if reactions are enabled
        if(performReaction)
            ofile << setw(19) << Temp[pLoc[i]];
 
        // Write variables
        for (int k = 0; k < nVar; k++) {
            ofile << setw(19) << (*varData[k])[pLoc[i]];
        }
        ofile << endl;
    }
 
    ofile.close();
   // cout << "Data successfully written to: " << fname << endl;
}
   
 
std::vector<double> hips::projection_back(std::vector<double> &vh) {
 
    int nh = xh.size() - 1;
    int nc = xc.size() - 1;
 
    std::vector<double> vc(nc, 0.0);
    int jprev = 0;
 
    for (double val : vc)
        std::cout << val << " ";
    std::cout << std::endl;
 
    for (int i = 0; i < nc; ++i) {
        for (int j = jprev + 1; j <= nh; ++j) {
            if (xh[j] <= xc[i + 1]) {
                double d1 = xh[j] - xh[j - 1];
                double d2 = xh[j] - xc[i];
                // calculation of shortest distance
                // handling if distance is zero but substarction gives comutational error
                double d  = std::min(d1, d2) < 1E-15 ? 0 : std::min(d1, d2);
 
                vc[i] += vh[j - 1] * d;
            } else {
                double d1 = xc[i + 1] - xh[j - 1];
                double d2 = xc[i + 1] - xc[i];
                // calculation of shortest distance
                // handling if distance is zero but substarction gives comutational error
                double d  = std::min(d1, d2) < 1E-15 ? 0 : std::min(d1, d2);
 
                vc[i] += vh[j - 1] * d;
 
                jprev = j - 1;
                break;
            }
        }
        vc[i] /= (xc[i + 1] - xc[i]);
    }
    return vc;
}
 
 
std::pair<std::vector<double>, std::vector<double>> hips::projection_back_with_density(std::vector<double> &vh, 
                                                                                       std::vector<double> &rho_h) {
    int nh = xh.size() - 1;
    int nc = xc.size() - 1;
 
    std::vector<double> vc(nc, 0.0);
    std::vector<double> rho_c(nc, 0.0);
    int jprev = 0;
 
    for (int i = 0; i < nc; ++i) {
        double total_dx = 0.0;
 
        for (int j = jprev + 1; j <= nh; ++j) {
            if (xh[j] <= xc[i + 1]) {
                double d = std::min(xh[j] - xh[j - 1], xh[j] - xc[i]);
                total_dx += d;
                rho_c[i] += rho_h[j - 1] * d;
                vc[i] += vh[j - 1] * rho_h[j - 1] * d;
            } else {
                double d = std::min(xc[i + 1] - xh[j - 1], xc[i + 1] - xc[i]);
                total_dx += d;
                rho_c[i] += rho_h[j - 1] * d;
                vc[i] += vh[j - 1] * rho_h[j - 1] * d;
 
                jprev = j - 1;
                break;
            }
        }
 
        // Normalize the results
        vc[i] /= rho_c[i];
    }
    return {vc, rho_c};
}
 
 
std::vector<std::vector<double>> hips::get_varData(){
 
    std::vector<std::vector<double>> varDataProjections;               // Vector to store modified data projections
 
    for (int i = 0; i < varData.size(); i++)                           // Loop through each element of varData and project the data back
        varDataProjections.push_back(projection_back((*varData[i])));  // Project the data back and store it in vh
 
    return varDataProjections;                                         // Return the vector containing modified data projections
}
 
std::vector<std::pair<std::vector<double>, std::vector<double>>> hips::get_varData_with_density() {
 
    std::vector<std::pair<std::vector<double>, std::vector<double>>> varDataProjections;  // Vector to store data projections with densities
 
    for (int i = 0; i < varData.size(); i++) {
        // Extract the value and density data
        std::vector<double> vh = (*varData[i]);  // Assuming varData[i][0] holds the values
        std::vector<double> rho_h = varRho;      // Assuming varData[i][1] holds the densities
 
        // Project the data back with density and store the result
        varDataProjections.push_back(projection_back_with_density(vh, rho_h));
    }
 
    return varDataProjections;  // Return the vector containing modified data projections with densities
}
 
 
void hips::setOutputIntervalEddy(int interval) {
 
    outputIntervalEddy = interval;
    useEddyBasedWriting = true;  
    useTimeBasedWriting = false; 
}
 
 
void hips::setOutputIntervalTime(double interval) {
 
    outputIntervalTime = interval;
    useTimeBasedWriting = true;  
    useEddyBasedWriting = false; 
}
 
 
void hips::saveAllParameters() {
 
    std::string filepath = "../post/parameters.dat";  
    std::ofstream file(filepath);
 
    if (!file) {
        std::cerr << "Error: Could not open " << filepath << " for writing!\n";
        return;
    }
 
    // Write user-defined input parameters
    file << "nLevels " << nLevels << "\n";            
    file << "domainLength " << domainLength << "\n";  
    file << "tau0 " << tau0 << "\n";                  
    file << "C_param " << C_param << "\n";            
    file << "forceTurb " << forceTurb << "\n";        
    file << "nVar " << nVar << "\n";                  
    file << "performReaction " << performReaction << "\n";  
    file << "realization " << realization << "\n";    
 
    // Write variable names
    if (!varName.empty()) {
        file << "varName ";
        for (const std::string &name : varName) {
            file << name << " ";  
        }
        file << "\n";
    } else {
        file << "varName (undefined)\n";
    }
 
    // Write i_batchelor vector if it exists and has the correct size
    if (!i_batchelor.empty() && i_batchelor.size() == nVar) {
        file << "i_batchelor ";
        for (const auto &val : i_batchelor) {
            file << val << " ";  
        }
        file << "\n";
    } else {
        file << "i_batchelor (undefined or size mismatch)\n";
    }
 
    file.close();
    std::cout << "All parameters saved in: " << filepath << std::endl;
}