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>
#include <algorithm>
 
using namespace std;
 
hips::hips(int nLevels_, 
           double domainLength_, 
           double tau0_, 
           double C_param_, 
           bool 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_),                                     
    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);
 
        if(forceTurb)
            throw std::runtime_error("Error: forceTurb should be false if preformReaction is true");
    }
    #endif
 
    // Resize vectors to the number of variables
    varData.resize(nVar);
    varName.resize(nVar); 
 
    set_tree(nLevels, domainLength, tau0);
}
 
 
hips::hips(double C_param_, 
           bool forceTurb_,
           int nVar_,
           vector<double> &ScHips_,
           bool performReaction_,
           shared_ptr<void> vcantSol,
           int seed,
           int realization_): 
    C_param(C_param_), 
    forceTurb(forceTurb_),       
    nVar(nVar_),                       
    ScHips(ScHips_),                 
    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_){    
 
    nLevels= nLevels_; 
    domainLength = domainLength_; 
    tau0 = tau0_; 
 
    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);
    wPar.assign(nparcels, 1.0 / nparcels);
 
    pLoc.resize(nparcels);
    for (int i=0; i<nparcels; i++)
        pLoc[i] = i;
 
    currentIndex = 0.0;
} 
 
 
void hips::set_tree(double Re_, double domainLength_, double tau0_, std::string ReApproach_) {
    Re = Re_;
    domainLength = domainLength_;
    tau0 = tau0_;
    ReApproach = ReApproach_;
 
    double baseLevelEstimate = (3.0 / 4) * log(1 / Re) / log(Afac);                               // Calculate the base tree level estimate (non-integer)
    int baseLevel;
 
    if (ReApproach == "rounding") {
        baseLevel = round(baseLevelEstimate);                                                     // Round the base level to the nearest integer
    } 
    else if (ReApproach == "probability") {
        baseLevel = ceil(baseLevelEstimate);                                                      // Ceil the base level to the nearest integer
        int previousLevel = baseLevel - 1;
        Prob = baseLevelEstimate - previousLevel;                                                 // Calculate the probability
    } 
    else if (ReApproach == "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 (ReApproach == "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 ReApproach 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((ReApproach == "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 (ReApproach == "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 (ReApproach == "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);
    wPar.assign(nparcels, 1.0 / nparcels);
    pLoc.resize(nparcels);
    for (int i = 0; i < nparcels; ++i)
        pLoc[i] = i;
 
    currentIndex = 0.0;
}
 
 
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);
 
    varData[currentIndex] = std::make_shared<std::vector<double>>(results.first);
    varRho = results.second;
    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) {
    // Build CFD and HiPS grids
    xc = setGridCfd(weight);
    xh = setGridHips(nparcels);
 
    int nc = xc.size() - 1; // CFD cell count
    int nh = xh.size() - 1; // HiPS parcel count
 
    std::vector<double> vh(nh, 0.0);     // parcel-averaged phi
    std::vector<double> rho_h(nh, 0.0);  // parcel-averaged density
 
    int jprev = 0; // remember where we left off in CFD cells
 
    for (int i = 0; i < nh; i++) 
    {
        double M = 0.0;     // total mass in this parcel
        double Mphi = 0.0;  // total (mass * phi) in this parcel
 
        double parcel_length = xh[i + 1] - xh[i];
 
        for (int j = jprev + 1; j <= nc; ++j) 
        {
            // Find geometric overlap between CFD cell and HiPS parcel
            double overlap_start = std::max(xh[i], xc[j - 1]);
            double overlap_end   = std::min(xh[i + 1], xc[j]);
 
            double overlap_len = overlap_end - overlap_start;
            if (overlap_len <= 0.0) continue;
 
            // Effective length scaled by wPar[i]
            double effective_len = overlap_len * (wPar[i] / parcel_length);
 
            // Accumulate mass and mass*phi
            double rho = density[j - 1];
            double phi = vcfd[j - 1];
            M    += rho * effective_len;
            Mphi += rho * phi * effective_len;
 
            // If CFD cell ends after parcel end  move to next parcel
            if (xc[j] >= xh[i + 1]) {
                jprev = j - 1;
                break;
            }
        }
 
        // Convert total mass  average density and phi
        if (wPar[i] > 0.0) rho_h[i] = M / wPar[i];
        if (M > 0.0)       vh[i]    = Mphi / M;
    }
 
    return {vh, rho_h};
}
 
 
std::vector<double> hips::setGridCfd(std::vector<double> &w) {
 
    double sumw = 0.0;
    for(int i=0; i<w.size(); i++)
        sumw += w[i];
   
    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]/sumw;                                // 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 && 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
    // ---- cache indices once
    const int enthalpyIdx = getVariableIndex("enthalpy");
    std::vector<int> yIdx(nsp);
    for (int k = 0; k < nsp; ++k) {
        yIdx[k] = getVariableIndex(gas->speciesName(k)); // map species name -> var index
    }
 
    const int nP     = 1 << (Nm1 - iLevel);
    const int istart = iTree * nP;
    const int iend   = istart + nP;
 
    std::vector<double> y(nsp);
 
    for (int i = istart; i < iend; ++i) {
        const int ime = pLoc[i];
        const double dt = time - parcelTimes[ime];
 
        // Pull current state
        double h = (*varData[enthalpyIdx])[ime];
        for (int k = 0; k < nsp; ++k) {
            y[k] = (*varData[yIdx[k]])[ime];
        }
 
        // ---- store old density BEFORE chemistry
        const double rho_old = varRho[ime];
 
        if (performReaction) {
            // Advance chemistry; bRxr updates its state internally
            bRxr->react(h, y, dt);
 
            // Get new density from the reactor/EOS
            const double rho_new = bRxr->getDensity();
 
            // Keep per-parcel mass m = rho * wPar * V_tot constant
            if (rho_new > 0.0) {
                wPar[ime] *= (rho_old / rho_new);
                varRho[ime] = rho_new;
            }
 
        }
 
        // Write back enthalpy and species (post-reaction)
        (*varData[enthalpyIdx])[ime] = h;
        for (int k = 0; k < nsp; ++k) {
            (*varData[yIdx[k]])[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.0;  // sum (value or value*mass)
    double msum = 0.0;  // mass sum (only used if weighted)
 
    for (int i = istart; i < iend; i++) {
        ime = pLoc[i];
        if (performReaction) {
            double m = varRho[ime] * wPar[ime];
            s    += (*varData[kVar])[ime] * m;   
            msum += m;
        } else {
            s += (*varData[kVar])[ime];          
        }
    }
 
    double avg = performReaction
               ? ((msum > 0.0) ? (s / msum) : 0.0)   
               : (s / nPmix);
 
    for (int i = istart; i < iend; i++) {
        ime = pLoc[i];
        (*varData[kVar])[ime] = avg;              
    }
 
    //---------- Mix right branch of iTree ----------
    istart = iend;
    iend   = istart + nPmix;
 
    s    = 0.0;
    msum = 0.0;
 
    for (int i = istart; i < iend; i++) {
        ime = pLoc[i];
        if (performReaction) {
            double m = varRho[ime] * wPar[ime];
            s    += (*varData[kVar])[ime] * m;       
            msum += m;
        } else {
            s += (*varData[kVar])[ime];              
        }
    }
 
    avg = performReaction
        ? ((msum > 0.0) ? (s / msum) : 0.0)
        : (s / nPmix);
 
    for (int i = istart; i < iend; i++) {
        ime = pLoc[i];
        (*varData[kVar])[ime] = avg;                 
    }
}
 
 
void hips::forceProfile() {
    // Loop through each variable in the HiPS profile
    for (int k = 0; k < varData.size(); k++) {
        double s=0;                                                  // 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());
    cout << endl << "writing data for time " << outputTime << " to file: " << 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
        #ifdef REACTIONS_ENABLED
        if(performReaction) {
            vector<double> yy(nsp);
            for(int k=0; k<nsp; k++)
                yy[k] = (*varData[k+1])[pLoc[i]];
            gas->setMassFractions(yy.data());
            gas->setState_HP((*varData[0])[pLoc[i]], gas->pressure());
            ofile << setw(19) << gas->temperature();
        }
        #endif
 
        // 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;
}
   
 
//  cfd
//  i=    0         1         2         3
//   |    *    |    *    |    *    |    *    |
 
//  hips
//   | * | * | * | * | * | * | * | * | * | * |
//  j= 0   1   2   3   4   5   6   7   8   9
 
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 (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::vector<double> hips::projection_back_with_density(std::vector<double> &vh, 
                                                       std::vector<double> &rho_h,
                                                       std::vector<double> &rho_c) {
    const int nh = static_cast<int>(xh.size()) - 1;  // # HiPS parcels
    const int nc = static_cast<int>(xc.size()) - 1;  // # CFD cells
 
    std::vector<double> phi_c(nc, 0.0);   // CFD-side variable (output)
    std::vector<double> M_cfd(nc, 0.0);   // mass in each CFD cell
    std::vector<double> Mphi_cfd(nc, 0.0);// mass*phi in each CFD cell
 
    int jprev = 0;
    for (int i = 0; i < nc; ++i) {
        // cell geometry (normalized length, since xc is 0..1)
        const double cell_len = xc[i + 1] - xc[i];
 
        for (int j = jprev + 1; j <= nh; ++j) {
            // geometric overlap between parcel j-1 and CFD cell i
            const double overlap_start = std::max(xh[j - 1], xc[i]);
            const double overlap_end   = std::min(xh[j],     xc[i + 1]);
            const double overlap_len   = overlap_end - overlap_start;
            if (overlap_len <= 0.0) continue;
 
            // scale overlap by current parcel volume fraction wPar
            const double parcel_len    = xh[j] - xh[j - 1]; // = 1.0/nh
            const double effective_len = overlap_len * (wPar[pLoc[j - 1]] / parcel_len);
 
            // accumulate mass and mass*phi into CFD cell i
            const double rhoP = rho_h[j - 1];
            const double phiP = vh[j - 1];
            M_cfd[i]    += rhoP * effective_len;
            Mphi_cfd[i] += rhoP * phiP * effective_len;
 
            // advance parcel index if we reached end of this CFD cell
            if (xc[i + 1] <= xh[j]) { jprev = j - 1; break; }
        }
 
        // recover CFD density and variable
        rho_c[i] = (cell_len > 0.0)        ? (M_cfd[i] / cell_len) : 0.0;
        phi_c[i] = (M_cfd[i] > 1e-300)     ? (Mphi_cfd[i] / M_cfd[i]) : 0.0;
    }
    return phi_c;
}
 
 
std::vector<std::vector<double>> hips::get_varData() {
    std::vector<std::vector<double>> varDataProjections;
 
    for (int i = 0; i < varData.size(); i++) {
 
        // Reorder using pLoc
        std::vector<double> vh(nparcels);
        for (int j = 0; j < nparcels; j++) {
            vh[j] = (*varData[i])[pLoc[j]];
        }
 
        std::vector<double> vc = projection_back(vh);
        varDataProjections.push_back(vc);
    }
 
    return varDataProjections;
}
 
 
std::pair<std::vector<std::vector<double>>, std::vector<double>> hips::get_varData_with_density() {
    std::vector<std::vector<double>> varDataProjections;
    
    // Reorder varRho based on pLoc
    std::vector<double> rho_h(nparcels);
    for (int i = 0; i < nparcels; i++) {
        rho_h[i] = varRho[pLoc[i]];
    }
 
    std::vector<double> rho_c = projection_back(rho_h);
 
    for (int i = 0; i < varData.size(); i++) {
 
        // Reorder variable data using pLoc
        std::vector<double> vh(nparcels);
        for (int j = 0; j < nparcels; j++) {
            vh[j] = (*varData[i])[pLoc[j]];
        }
 
        auto vc = projection_back_with_density(vh, rho_h, rho_c);
        varDataProjections.push_back(vc);
    }
 
    return {varDataProjections, rho_c};
}
 
 
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();
    //cout << endl << "All parameters saved in: " << filepath << std::endl;
}