#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <string>
#include <sstream>
#include <cstdlib>

#include<stdlib.h>

#include <tuple>

#include <fstream> // for reading and writing files

#include <stdio.h>
#include <omp.h>

#include <chrono> // for timing of processes

// *****************************************************************************************************************************************
// Single-component multiphase LBM implementation for simulation of water retention behaviour in granular media
// Version 0: 06/2022
// *****************************************************************************************************************************************
// Author: Marius Milatz, Institute of Geotechnical Engineering and Construction Management, Hamburg University of Technology (TUHH)
// marius.milatz@tuhh.de
// *****************************************************************************************************************************************
// This multiphase Lattice Boltzmann model is based on the single-component Shan-Chen multiphase Lattice Boltzmann formulation
// as presented in the textbook by Huang et al.:

// [1] Huang, H., M. C. Sukop, and X.-Y. Lu (2015): Multiphase lattice Boltzmann methods: Theory and application. Chichester:
// Wiley Blackwell.

// This C++ implementation of the multiphase LBM is based on the Fortran code presented by Huang et al. (2015) in Chapter 2 and related
// Appendix of their textbook which is originally intended to model a "single component multiphase droplet contacting a wall".
// It has also been published online for own use and teaching purposes in the repository of the textbook [1].

// For the simulation of water retention behaviour in granular media, the Fortran code from [1] has been transferred to C++ with the
// possibility of parallelisation using OpenMP and it has been enhanced to consider simulation domains based on CT input data, to initialise
// density from a perturbated random field of densities leading to an initial phase separation step, to change the macroscopic degree of
// saturation after phase separation based on a density injection scheme as described in [2], to write relevant macroscopic data, such
// as macroscopic degree of saturation and capillary pressure, and to follow prescribed paths of degree of saturation.

// [2] Hosseini, R., K. Kumar, and J.-Y. Delenne (2021): Investigating the effect of porosity on the soil water retention curve using
// the multiphase Lattice Boltzmann Method. In: Proc. of Powders & Grains 2021 - 9th International Conference on Micromechanics
// on Granular Media. DOI: 10.1051/epjconf/202124909007.
// *****************************************************************************************************************************************

using real = double;

// Global variables:
int selected_number_of_threads = 28;		// Number of threads for multi-threaded calculations

real density_injection_rate = 0.005;		// Increment of density injection
int delta_t_density_injection = 1000;		// Time increment for density injection

int step_number = 0;						// Identifier for every hydraulic step

// In case of parallelisation, compile code with: g++ -o LBM LBM.cpp -fopenmp
const int Q = 19;
const bool output_info = true;				//output information on data structures etc.
bool debug = false;							// Output stuff for debugging (active if debug = true)
bool run_simulation = true;					// Switch for starting the main simulation

const int lx = 200, ly = 200, lz = 200;		// Computational domain size

const real rho_boundary_liquid = 0.2726; 	// Threshold density for liquid
const real rho_boundary_gas = 0.0484;		// Threshold density for gas


const int k_sel = 10;           // important: use k to select TTOW in EOS etc.


real rho_0 = 0.1;				// Initial density used in initial perturbation step
real RR;                        // bubble radius
real rho_w;                     // wall density
real tau;                       // tau
real visc; 
int t_max;                      // max. time increment
int Nwri;                       // interval for writing data

// Definition of D3Q19 velocity model: xc[ex], yc[ey], zc[ez]
const real xc[Q] = {0., 1., -1., 0., 0., 0., 0., 1., 1., -1., -1., 1., -1., 1., -1., 0., 0., 0., 0.};
const real yc[Q] = {0., 0., 0. , 1., -1., 0., 0., 1., -1., 1., -1., 0., 0., 0., 0., 1., 1., -1., -1.};
const real zc[Q] = {0., 0., 0., 0., 0., 1., -1., 0., 0., 0., 0., 1., 1., -1., -1., 1., -1., 1., -1.};

const signed int ex[Q] = {0, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, 1, -1, 1, -1, 0, 0, 0, 0};
const signed int ey[Q] = {0, 0, 0, 1, -1, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 1 ,1 ,-1, -1};
const signed int ez[Q] = {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, 1, -1, 1, -1};    

// The array opp gives the opposite direction for e_1, e_2, e_3, ..., e_18.
// It implements the simple bounce-back rule we use in the collision step for solid nodes (obst=1)
const int opp[Q] = {0, 2, 1, 4, 3, 6, 5, 10, 9, 8, 7, 14, 13, 12, 11, 18, 17, 16, 15}; //mm: added a 0 as entry 0 of array!!!

// Carnahan-Starling EOS
// RHW and RLW are the coexisting densities in the corresponding specified T/T0.
//
const real TT0W[12] ={0.975, 0.95, 0.925, 0.9, 0.875, 0.85, 0.825, 0.8, 0.775, 0.75, 0.7, 0.65};

//                  k  =  0     1      2     3       4      5      6      7       8     9      10     11
const real RHW[12]  = {0.16, 0.210, 0.23, 0.2470, 0.265, 0.279, 0.290, 0.3140, 0.30, 0.330, 0.358, 0.380};
const real RLW[12]  = {0.08, 0.067, 0.05, 0.0405, 0.038, 0.032, 0.025, 0.0245, 0.02, 0.015, 0.009, 0.006};


const real TT0   = TT0W[k_sel];
const real rho_h = RHW[k_sel];
const real rho_l = RLW[k_sel];


// Speeds and weighting factors
const real cc = 1.;
const real c_squ = cc *cc/3.;
const real t_0 = 1./3.;
const real t_1 = 1./18.;
const real t_2 = 1./36.;

const real t_k[Q] = { t_0, t_1, t_1, t_1, t_1, t_1, t_1, t_2, t_2, t_2, t_2, t_2, t_2, t_2, t_2, t_2, t_2, t_2, t_2 };
  
int obst[lx][ly][lz] = {0};           // boolean array containing solid = 1 and fluid nodes = 0

real u_x[lx][ly][lz] = {0.};
real u_y[lx][ly][lz] = {0.};
real u_z[lx][ly][lz] = {0.};

real rho[lx][ly][lz] = {0.};

real p[lx][ly][lz]   = {0.};
real upx[lx][ly][lz] = {0.};
real upy[lx][ly][lz] = {0.};
real upz[lx][ly][lz] = {0.};

real Fx[lx][ly][lz]  = {0.};
real Fy[lx][ly][lz]  = {0.};
real Fz[lx][ly][lz]  = {0.};

real Sx[lx][ly][lz]  = {0.};
real Sy[lx][ly][lz]  = {0.};
real Sz[lx][ly][lz]  = {0.};

real psx[lx][ly][lz] = {0.};

real u_n[Q]          = {0.};
real fequi[Q]        = {0.};
real fequ[Q]         = {0.};
real ff[Q][lx][ly][lz] = {0.};
real f_hlp[Q][lx][ly][lz] = {0.};

int pore_count = 0;		// Pore cell counter
int solid_count = 0;	// Solid cell counter

int liquid_count = 0;	// Liquid cell counter
int gas_count = 0;		// Gas cell counter
int vapour_count = 0;	// Vapour cell counter

real mean_gas_pressure = 0.;		// Mean gas pressure
real mean_liquid_pressure = 0.;		// Mean liquid pressure
real mean_vapour_pressure = 0.;		// Mean vapour pressure

real saturation = 0.;				// Degree of saturation
real capillary_pressure = 0.;		// Capillary pressure


/// @brief Read simulation parameters...
std::tuple<real, real, real, int, int> read_parameters(){
    real value[5] = {0.};
    std::cout << "Reading simulation parameters..." << std::endl;
    std::fstream myfile("params.txt", std::ios_base::in);


    std::ifstream data("params.txt");
    if (!data.is_open())
    {
        std::exit(EXIT_FAILURE);
    }
    std::string str;
    std::getline(data, str); // skip the first line
    while (std::getline(data, str))
    {
        std::istringstream iss(str);
        std::string token;
        int i = 0;

        while (std::getline(iss, token, ','))
        {   
            value[i] = std::stof(token);    // stof converts string to float
            i++;
            // process each token
            std::cout << token << " ";
        }
        std::cout << std::endl;
        
    }

    for (int i = 0; i<5; i++){
        std::cout << "value(" << i << ") = " << value[i] << std::endl;
    }
        real RR = value[0];                        // bubble radius
        real rho_w = value[1];                     // wall density
        real tau = value[2];                       // tau
        int t_max = static_cast<int>(value[3]);     // max. time increment
        int Nwri  = static_cast<int>(value[4]);     // interval for writing data
     

    return std::make_tuple(RR, rho_w, tau, t_max, Nwri);
}

/// @brief Function to read geometry data from file...
void  read_geometry_data(const char* name){ 
std::ifstream file(name);
  if (file.fail()) {
    std::cerr << "Cannot open sample file " << name << std::endl;
  }

  std::vector<int> ids;
  std::cout << "Reading sample file: " << name << std::endl;


 if (file.is_open() && file.good()) {
    std::string line;


    ids.reserve(1000);
    // Read the sample
    unsigned k = 0;
    int x, y, z;
    
    while (std::getline(file, line)) {
      std::istringstream istream(line);
      int vid;

      istream >> vid;
      ids.emplace_back(vid);

      ++k;

    }    
 }

  std::cout << __FILE__ << __LINE__ << std::endl;
  file.close(); 
  std::cout << __FILE__ << __LINE__ << std::endl;  
  int ID;
  unsigned k = 0;
  for (int z = 0; z < lz; z++) {
    for (int y = 0; y < ly; y++) {
      for (int x = 0; x < lx; x++) {
        ID = ids.at(k);         // access vector at entry k
        obst[x][y][z] = ID; //!ID;  // ID of 1 is solid and 0 is fluid
        ++k;
      }
    }
  }
  std::cout << __FILE__ << __LINE__ << std::endl;

}

/// @brief Function to initialize the density...
void init_density(){

    real my_sign;
    real my_val;
    real my_randnum;

    real u_squ;

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                my_val = ((real)rand() / (real)RAND_MAX);
                    if (my_val <= 0.5){
                         my_sign = -1.0;
                    }
                    else {
                        my_sign = 1.0;
                    }
                my_randnum = rho_0 + 1e-3*((real)rand() / (real)RAND_MAX)*my_sign;

                rho[x][y][z] = my_randnum;


            }
        }
    }

     for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                u_squ = u_x[x][y][z]*u_x[x][y][z] + u_y[x][y][z]*u_y[x][y][z] + u_z[x][y][z]*u_z[x][y][z];

                for (int k = 0; k < Q; k++){

                    u_n[k] = xc[k]*u_x[x][y][z] + yc[k]*u_y[x][y][z] + zc[k]*u_z[x][y][z];

                    fequi[k] = t_k[k] * rho[x][y][z] * (cc*u_n[k] / c_squ + (u_n[k]*cc) * (u_n[k]*cc) / (2.0 * c_squ * c_squ) - u_squ / (2.0 * c_squ)) + t_k[k] * rho[x][y][z];
                    ff[k][x][y][z] = fequi[k];                    

                }

            }
        }
    }   

}

/// @brief Function for writing restart files...
void write_data_for_restart();				// Currently not implemented.

/// @brief Function for writing the results file for every time step...
void get_and_write_unsat_results(int sim_time){

    solid_count = 0;
    pore_count = 0;
    liquid_count = 0;
    gas_count = 0;
    vapour_count = 0;

    mean_gas_pressure = 0.;
    mean_liquid_pressure = 0.;
    saturation = 0.;
    capillary_pressure = 0.;

    real sum_liquid_pressure = 0.;
    real sum_gas_pressure = 0.;

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){
                
                if(obst[x][y][z] == 0){ // obst = 0: pore cells
                    pore_count += 1;
                    if(rho[x][y][z] > rho_boundary_liquid){     // liquid cells
                        liquid_count += 1;
                        sum_liquid_pressure += p[x][y][z];
                    }
                    else if(rho[x][y][z] < rho_boundary_gas){        // gas cells
                        gas_count += 1;
                        sum_gas_pressure += p[x][y][z];
                    }    
                    else if((rho[x][y][z] >= rho_boundary_gas) && (rho[x][y][z] <= rho_boundary_liquid)){    // vapour cells
                        vapour_count += 1;
                    }                                      
                }
                if(obst[x][y][z] == 1){ // obst = 1: solid cells
                    solid_count += 1;
                }

            }
        }
    }

    saturation = float(liquid_count) / float(pore_count);
    mean_liquid_pressure = sum_liquid_pressure / liquid_count;
    mean_gas_pressure = sum_gas_pressure / gas_count;

    capillary_pressure = mean_gas_pressure - mean_liquid_pressure;

    std::cout << "=========================================" << std::endl;
    std::cout << "Number of gas cells = " << gas_count << std::endl;
    std::cout << "Number of liquid cells = " << liquid_count << std::endl;
    std::cout << "Number of vapour cells = " << vapour_count << std::endl;
    std::cout << "Number of pore cells = " << pore_count << std::endl;
    std::cout << "Number of solid cells = " << solid_count << std::endl;    
    std::cout << "Current saturation = " << saturation << std::endl;
    std::cout << "Current capillary pressure = " << capillary_pressure << std::endl;
    std::cout << "=========================================" << std::endl;


    std::string my_time = std::to_string(sim_time);

    std::string my_filename = "unsat_results.txt";

    std::cout << "----------------------------------------------------------------------" << std::endl;
    std::cout << "Info: Writing unsat results to file " << my_filename <<  " at time = " << my_time << std::endl;
    std::cout << "----------------------------------------------------------------------" << std::endl;     

    std::ofstream myfile;


    myfile.open (my_filename, std::ios_base::app);

    if (sim_time == 1){
    myfile << "time,\t step_number,\t Sr,\t p_mean_gas,\t p_mean_liquid,\t pc,\t gas cells,\t vapour cells,\t liquid cells,\t pore cells,\t solid cells\n";

    myfile << sim_time << ",\t" << step_number << ",\t" << saturation << ",\t" << mean_gas_pressure
           << ",\t" << mean_liquid_pressure << ",\t" << capillary_pressure
           << ",\t" << gas_count << ",\t" << vapour_count << ",\t" << liquid_count
           << ",\t" << pore_count << ",\t" << solid_count << std::endl;
    }
    else{
    myfile << sim_time << ",\t" << step_number << ",\t" << saturation << ",\t" << mean_gas_pressure
           << ",\t" << mean_liquid_pressure << ",\t" << capillary_pressure
           << ",\t" << gas_count << ",\t" << vapour_count << ",\t" << liquid_count
           << ",\t" << pore_count << ",\t" << solid_count << std::endl;
    }    

    myfile.close();


}

/// @brief Function for writing spatial results to file at spaced intervals...
void write_results(int sim_time){

    std::string my_time = std::to_string(sim_time);
    
    std::string my_filename = my_time.append(".plt");

    my_time = std::to_string(sim_time);


    std::cout << "--------------------------------------------------" << std::endl;
    std::cout << "Info: Writing results to file " << my_filename <<  " at time = " << my_time << std::endl;
    std::cout << "--------------------------------------------------" << std::endl;    


    std::ofstream myfile;
    myfile.open (my_filename);

    myfile << "x, \t y, \t z, \t rho, \t upx, \t upy, \t upz, \t p \t obst\n";
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){
                myfile << x << "\t" << y << "\t" << z << "\t" << rho[x][y][z] << "\t"
                       << upx[x][y][z] << "\t" << upy[x][y][z] << "\t" << upz[x][y][z] << "\t"
                       << p[x][y][z] << "\t" << obst[x][y][z] << std::endl;
            }
        }
    }      

    myfile.close();

}

/// @brief LBM-streaming step...
void stream(){

    int z_n, y_n, x_e, z_s, y_s, x_w;

#ifdef _OPENMP
#pragma omp parallel for private(z_n,y_n,x_e,z_s,y_s,x_w) schedule(runtime) num_threads(selected_number_of_threads)
#endif
    for (int z = 0; z < lz; z++){
       for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){
                z_n = (z % lz) + 1;
                y_n = (y % ly) + 1;
                x_e = (x % lx) + 1;
                z_s = lz - ((lz + 1 - z) % lz);
                y_s = ly - ((ly + 1 - y) % ly);
                x_w = lx - ((lx + 1 - x) % lx);

                if (x == 1){
                    x_w = 0;
                }
                if (y == 1){
                    y_s = 0;
                }
                if (z == 1){
                    z_s = 0;
                }


                if (z == lz - 1){
                    z_n = 0;
                }
                if (y == ly - 1){
                    y_n = 0;
                }
                if (x == lx - 1){
                    x_e = 0;
                }


                // Propagation:
                f_hlp[1][x_e][y][z  ] = ff[1][x][y][z];
                f_hlp[2][x_w][y][z  ] = ff[2][x][y][z];
                f_hlp[3][x  ][y_n][z  ] = ff[3][x][y][z];
                f_hlp[4][x  ][y_s][z  ] = ff[4][x][y][z];
                f_hlp[5][x  ][y  ][z_n] = ff[5][x][y][z];
                f_hlp[6][x  ][y  ][z_s] = ff[6][x][y][z];
                f_hlp[7][x_e][y_n][z  ] = ff[7][x][y][z];
                f_hlp[8][x_e][y_s][z  ] = ff[8][x][y][z];
                f_hlp[9][x_w][y_n][z  ] = ff[9][x][y][z];
                f_hlp[10][x_w][y_s][z ] = ff[10][x][y][z];
                f_hlp[11][x_e][y  ][z_n] = ff[11][x][y][z];
                f_hlp[12][x_w][y  ][z_n] = ff[12][x][y][z];
                f_hlp[13][x_e][y  ][z_s] = ff[13][x][y][z];
                f_hlp[14][x_w][y  ][z_s] = ff[14][x][y][z];
                f_hlp[15][x  ][y_n][z_n] = ff[15][x][y][z];
                f_hlp[16][x  ][y_n][z_s] = ff[16][x][y][z];
                f_hlp[17][x  ][y_s][z_n] = ff[17][x][y][z];
                f_hlp[18][x  ][y_s][z_s] = ff[18][x][y][z];
            }
        }
    }


#ifdef _OPENMP
#pragma omp parallel for schedule(runtime) num_threads(selected_number_of_threads)
#endif
    // Update the distribution function:
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){
                for (int k = 1; k < Q; k++){
                    ff[k][x][y][z] = f_hlp[k][x][y][z];
                }
            }
        }
    }
}

/// @brief Function for calculating the macro variables...
void getuv(){

#ifdef _OPENMP
#pragma omp parallel for schedule(runtime) num_threads(selected_number_of_threads)
#endif    
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){
                rho[x][y][z] = 0.;
                if (obst[x][y][z] == 0){
                    for (int k = 0; k < Q; k++){
                        rho[x][y][z] += ff[k][x][y][z];
                    }
                    if (rho[x][y][z] != 0.){
                        u_x[x][y][z] = (ff[1][x][y][z] + ff[7][x][y][z] + ff[8][x][y][z] + ff[11][x][y][z] + ff[13][x][y][z] - (ff[2][x][y][z] + ff[9][x][y][z] + ff[10][x][y][z] + ff[12][x][y][z] + ff[14][x][y][z]))/rho[x][y][z];
                        u_y[x][y][z] = (ff[3][x][y][z] + ff[7][x][y][z] + ff[9][x][y][z] + ff[15][x][y][z] + ff[16][x][y][z] - (ff[4][x][y][z] + ff[8][x][y][z] + ff[10][x][y][z] + ff[17][x][y][z] + ff[18][x][y][z]))/rho[x][y][z];
                        u_z[x][y][z] = (ff[5][x][y][z] + ff[11][x][y][z] + ff[12][x][y][z] + ff[15][x][y][z] + ff[17][x][y][z] - (ff[6][x][y][z] + ff[13][x][y][z] + ff[14][x][y][z] + ff[16][x][y][z] + ff[18][x][y][z]))/rho[x][y][z];

                    }

                }
            }
        }
    }
}

/// @brief Function for calculating the actual velocity...
void calcu_upr(){

#ifdef _OPENMP
#pragma omp parallel for schedule(runtime) num_threads(selected_number_of_threads)
#endif 
for (int z = 0; z < lz; z++){
    for (int y = 0; y < ly; y++){
        for (int x = 0; x < lx; x++){
            if (obst[x][y][z] == 0){

                upx[x][y][z] = u_x[x][y][z] + (Fx[x][y][z] + Sx[x][y][z])/2.0/rho[x][y][z];
                upy[x][y][z] = u_y[x][y][z] + (Fy[x][y][z] + Sy[x][y][z])/2.0/rho[x][y][z];
                upz[x][y][z] = u_z[x][y][z] + (Fz[x][y][z] + Sz[x][y][z])/2.0/rho[x][y][z];
            } else {
                upx[x][y][z] = u_x[x][y][z];
                upy[x][y][z] = u_y[x][y][z];
                upz[x][y][z] = u_z[x][y][z];               
            }

            
        }
    }
}

}

/// @brief Function for calculating the interaction forces...
void calcu_Fxy(){

    const real R = 1.0;
    const real b = 4.0;
    const real a = 1.0;
    const real Tc = 0.3773*a/(b*R);
    const real TT = TT0*Tc;

    real sum_x;
    real sum_y;
    real sum_z;

    real G1;

    signed int xp;
    signed int yp;
    signed int zp;

    real psx_w = 0.;


#ifdef _OPENMP
#pragma omp parallel for schedule(runtime) num_threads(selected_number_of_threads)
#endif
    for (int k = 0; k < lz; k++){
        for (int j = 0; j < ly; j++){
            for (int i = 0; i < lx; i++){
                
                if ((obst[i][j][k] == 0) && (rho[i][j][k] != 0.)){

                    if ((R*TT*(1.0+(4.0*rho[i][j][k] - 2.0*rho[i][j][k]*rho[i][j][k])/pow((1.0-rho[i][j][k]),3.)) - a*rho[i][j][k] -1.0/3.0) > 0.){
                        G1 = 1.0/3.0;
                    }
                   
                    else{
                        G1 = -1.0/3.0;    
                    }  

                    psx[i][j][k] = sqrt(6.0*rho[i][j][k]*(R*TT*(1.0 + (4.0*rho[i][j][k] - 2.0*rho[i][j][k]*rho[i][j][k])/pow((1.0-rho[i][j][k]),3.)) - a*rho[i][j][k] - 1.0/3.0)/G1);

                    // Yuan Carnahan-Starling EOS:
                    p[i][j][k] = rho[i][j][k]/3.0 + G1/6.0 * psx[i][j][k]*psx[i][j][k];

                }
            }
        }
    }

    psx_w = sqrt(6.0*rho_w*(R*TT*(1.0 + (4.0*rho_w-2.0*rho_w*rho_w)/pow((1.0-rho_w),3.))-a*rho_w - 1.0/3.0)/G1);                                      


#ifdef _OPENMP
#pragma omp parallel for private(xp,yp,zp,sum_x,sum_y,sum_z) schedule(runtime) num_threads(selected_number_of_threads)
#endif
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                // Interaction between neighbouring with periodic boundaries:
                Fx[x][y][z] = 0.;
                Fy[x][y][z] = 0.;
                Fz[x][y][z] = 0.;

                if (obst[x][y][z] == 0){

                    sum_x = 0.;
                    sum_y = 0.;
                    sum_z = 0.;

                    for (int k = 1; k < Q; k++){

                        xp = x + ex[k];
                        yp = y + ey[k];
                        zp = z + ez[k];



                        if (xp < 0){
                            xp = lx-1;
                        }
                        if (xp > (lx-1)){
                            xp = 0;
                        }
                        if (yp < 0){
                            yp = ly-1;
                        }
                        if (yp > (ly-1)){
                            yp = 0;
                        }  
                        if (zp < 0){
                            zp = lz-1;
                        }
                        if (zp > (lz-1)){
                            zp = 0;
                        } 

                        

                        if (obst[xp][yp][zp] == 1){

                            // Interact with solid nodes (obst=1):
                            sum_x += t_k[k]*xc[k];
                            sum_y += t_k[k]*yc[k];
                            sum_z += t_k[k]*zc[k];
                        }
                        else{
                            
                            // Interact with fluid nodes (obst=0):
                            Fx[x][y][z] += t_k[k]*xc[k]*psx[xp][yp][zp];
                            Fy[x][y][z] += t_k[k]*yc[k]*psx[xp][yp][zp];
                            Fz[x][y][z] += t_k[k]*zc[k]*psx[xp][yp][zp];

                        }                                              

                    }

                    // Final wall-fluid interaction:
                    Sx[x][y][z] = -G1*sum_x*psx[x][y][z]*psx_w;
                    Sy[x][y][z] = -G1*sum_y*psx[x][y][z]*psx_w;
                    Sz[x][y][z] = -G1*sum_z*psx[x][y][z]*psx_w;    

                    // Final fluid-fluid interaction:
                    Fx[x][y][z] = -G1*psx[x][y][z]*Fx[x][y][z];
                    Fy[x][y][z] = -G1*psx[x][y][z]*Fy[x][y][z];
                    Fz[x][y][z] = -G1*psx[x][y][z]*Fz[x][y][z];                    

                }

            }     
        }
    }

}

/// @brief collision step...
void collision(){

real temp[Q] = {0.};

real ux = 0.;
real uy = 0.;
real uz = 0.;
real u_squ = 0.;

#ifdef _OPENMP
#pragma omp parallel for private(temp) schedule(runtime) num_threads(selected_number_of_threads)
#endif
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if (obst[x][y][z] == 1){

                    for (int k = 1; k < Q; k++){

                        temp[k] = ff[k][x][y][z];

                    }

                    for (int k = 1; k < Q; k++){

                        ff[opp[k]][x][y][z] = temp[k];

                    }                    

                }

            }
        }
    }


#ifdef _OPENMP
#pragma omp parallel for private(ux,uy,uz,u_squ,u_n,fequ) schedule(runtime) num_threads(selected_number_of_threads)
#endif
    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if (obst[x][y][z] == 0){

                    ux = u_x[x][y][z] + tau * (Fx[x][y][z] + Sx[x][y][z])/rho[x][y][z];
                    uy = u_y[x][y][z] + tau * (Fy[x][y][z] + Sy[x][y][z])/rho[x][y][z];
                    uz = u_z[x][y][z] + tau * (Fz[x][y][z] + Sz[x][y][z])/rho[x][y][z];

                    u_squ = ux*ux + uy*uy + uz*uz;

                    for (int k = 0; k < Q; k++){

                        // Equilibrium distribution function:
                        u_n[k]  = xc[k]*ux + yc[k]*uy + zc[k]*uz;
                        fequ[k] = t_k[k]*rho[x][y][z]*(cc*u_n[k]/c_squ + (u_n[k]*cc)*(u_n[k]*cc)/(2.0 * c_squ*c_squ) - u_squ/(2.0*c_squ)) + t_k[k]*rho[x][y][z];
                        // Collision step:
                        ff[k][x][y][z] = fequ[k] + (1.0-1.0/tau)*(ff[k][x][y][z] - fequ[k]);                        

                    }

                }

            }
        }
    }

}

/// @brief Function for decreasing density...
void decrease_density(){

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if (obst[x][y][z] == 0){

                    for (int k = 0; k < Q; k++){

                        // Decrease density:
                        ff[k][x][y][z] = ff[k][x][y][z] - density_injection_rate/19.0;

                        if (ff[k][x][y][z] <= 0.0){

                            ff[k][x][y][z] = 0.009/19.0;

                        }

                    }

                }

            }
        }
    }

}

/// @brief Function for decreasing density in liquid nodes only...
void decrease_density_in_liquid_nodes_only(){

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if ((obst[x][y][z] == 0) && (rho[x][y][z] > rho_boundary_liquid)){

                    for (int k = 0; k < Q; k++){

                        // Decrease density:
                        ff[k][x][y][z] = ff[k][x][y][z] - density_injection_rate/19.0;

                    }

                }

            }
        }
    }

}

/// @brief Function for increasing density...
void increase_density(){

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if (obst[x][y][z] == 0){

                    for (int k = 0; k < Q; k++){

                        // Increase density:
                        ff[k][x][y][z] = ff[k][x][y][z] + density_injection_rate/19.0;

                    }

                }

            }
        }
    }


}

/// @brief Function for increasing density in liquid nodes only...
void increase_density_in_liquid_nodes_only(){

    for (int z = 0; z < lz; z++){
        for (int y = 0; y < ly; y++){
            for (int x = 0; x < lx; x++){

                if ((obst[x][y][z] == 0) && (rho[x][y][z] > rho_boundary_liquid)){

                    for (int k = 0; k < Q; k++){

                        // Increase density:
                        ff[k][x][y][z] = ff[k][x][y][z] + density_injection_rate/19.0;

                    }

                }

            }
        }
    }

}

/// @brief Main function...
int main(int argc, char const* argv[]) {

    std::cout << "=================" << std::endl;
    std::cout << "3D Multiphase LBM" << std::endl;
    std::cout << "=================" << std::endl;

    if (output_info == true) {
        std::cout << "Dimensions of calculation:" << std::endl;
        std::cout << "lx = " << lx << std::endl;
        std::cout << "ly = " << ly << std::endl;
        std::cout << "lz = " << lz << std::endl;
    }



    if (debug == true) {
        std::cout << "D3Q19 velocity model:" << std::endl;
        std::cout << "xc = ";
        for (int i = 0; i < Q; i++){   
            std::cout <<  xc[i] << ", " ; //<< std::endl;
        }
        std::cout << std::endl;
        std::cout << "yc = ";
        for (int i = 0; i < Q; i++){   
            std::cout <<  yc[i] << ", " ; //<< std::endl;
        }
        std::cout << std::endl;
        std::cout << "zc = ";
        for (int i = 0; i < Q; i++){   
            std::cout <<  zc[i] << ", " ; //<< std::endl;
        }
        std::cout << std::endl;

        std::cout << "t_k = ";
        for (int i = 0; i < Q; i++){   
            std::cout <<  t_k[i] << ", " ; //<< std::endl;
        }
        std::cout << std::endl;
    }

    if (output_info == true){
        std::cout << "-----------------" << std::endl;
        std::cout << "TT0 = " << TT0 << std::endl;
        std::cout << "rho_l = " << rho_l << std::endl;
        std::cout << "rho_h = " << rho_h << std::endl;
        std::cout << "-----------------" << std::endl;    
    }
// Read the simulation parameters from file:

    auto data = read_parameters();  // data is a tuple containing several return values!!!

    RR = std::get<0>(data);
    rho_w = std::get<1>(data);
    tau = std::get<2>(data);
    t_max = std::get<3>(data);
    Nwri = std::get<4>(data);

    real visc = c_squ*(tau-0.5);

    if (output_info == true){
        std::cout << "-----------------" << std::endl;
        std::cout << "RR = "    << RR    << std::endl;
        std::cout << "rho_w = " << rho_w << std::endl;
        std::cout << "tau = "   << tau   << std::endl;
        std::cout << "t_max = " << t_max << std::endl;
        std::cout << "Nwri = "  << Nwri  << std::endl;
        std::cout << "visc = "  << visc  << std::endl;
        std::cout << "-----------------" << std::endl;       
    } 

// Read the geometry of the simulation from file:
    std::cout << "Reading obstacles..." << std::endl;
    read_geometry_data("geometry_data.txt"); 
    std::cout << "Done." << std::endl;


    if (debug == true){
        std::cout << "Check obst entries..." << std::endl;
        std::cout << "obst[0][0][0] = " << obst[0][0][0] << std::endl;
        std::cout << "obst[47][0][0] = " << obst[47][0][0] << std::endl; 
        std::cout << "obst[48][0][0] = " << obst[48][0][0] << std::endl;         
        std::cout << "obst[74][0][0] = " << obst[74][0][0] << std::endl; 
        std::cout << "obst[79][0][0] = " << obst[79][0][0] << std::endl;    
        std::cout << "obst[80][0][0] = " << obst[80][0][0] << std::endl;    
        std::cout << "obst[81][0][0] = " << obst[81][0][0] << std::endl;  
        std::cout << "obst[82][0][0] = " << obst[82][0][0] << std::endl;  
        std::cout << "obst[83][0][0] = " << obst[83][0][0] << std::endl;     
        std::cout << "obst[84][0][0] = " << obst[84][0][0] << std::endl;   
        std::cout << "obst[85][0][0] = " << obst[85][0][0] << std::endl;  
        std::cout << "obst[86][0][0] = " << obst[86][0][0] << std::endl; 
        std::cout << "obst[87][0][0] = " << obst[87][0][0] << std::endl;    
        std::cout << "obst[100][0][0] = " << obst[100][0][0] << std::endl;
        std::cout << "obst[0][1][0] = " << obst[0][1][0] << std::endl;                                 

        int x, y, z;
        for (int z = 0; z < 5; z++) {
            for (int y = 0; y < 5; y++) {
                for (int x = 0; x < 5; x++) {

                    std::cout << x << '\t' << y << '\t' << z << '\t' << obst[x][y][z] << std::endl;;

                }
            }
        }
    }


  // Initialize density:
    std::cout << "Initializing density..." << std::endl;
    init_density(); 
    std::cout << "Done." << std::endl; 


    if (debug == true){
        std::cout << "Check rho entries..." << std::endl;
        std::cout << "rho[0][0][0] = " << rho[0][0][0] << std::endl;
        std::cout << "rho[47][0][0] = " << rho[47][0][0] << std::endl; 
        std::cout << "rho[48][0][0] = " << rho[48][0][0] << std::endl;         
        std::cout << "rho[74][0][0] = " << rho[74][0][0] << std::endl; 
        std::cout << "rho[79][0][0] = " << rho[79][0][0] << std::endl;    
        std::cout << "rho[80][0][0] = " << rho[80][0][0] << std::endl;    
        std::cout << "rho[81][0][0] = " << rho[81][0][0] << std::endl;  
        std::cout << "rho[82][0][0] = " << rho[82][0][0] << std::endl;  
        std::cout << "rho[83][0][0] = " << rho[83][0][0] << std::endl;     
        std::cout << "rho[84][0][0] = " << rho[84][0][0] << std::endl;   
        std::cout << "rho[85][0][0] = " << rho[85][0][0] << std::endl;  
        std::cout << "rho[86][0][0] = " << rho[86][0][0] << std::endl; 
        std::cout << "rho[87][0][0] = " << rho[87][0][0] << std::endl;    
        std::cout << "rho[100][0][0] = " << rho[100][0][0] << std::endl;
        std::cout << "rho[0][1][0] = " << rho[0][1][0] << std::endl;                                 

        int x, y, z;
         for (int z = 0; z < 5; z++) {
            for (int y = 0; y < 5; y++) {
                for (int x = 0; x < 5; x++) {

                    std::cout << x << '\t' << y << '\t' << z << '\t' << rho[x][y][z] << std::endl;;

                }
            }
        }
    }  



    if (debug == true){
        std::cout << "Calling stream()..." << std::endl;
        auto t1 = std::chrono::high_resolution_clock::now();
        stream(); 
        auto t2 = std::chrono::high_resolution_clock::now();
        std::cout << "Done." << std::endl;

        std::cout << "stream() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t2-t1).count()
                  << " milliseconds\n";

        std::cout << "Calling getuv()..." << std::endl;
        auto t3 = std::chrono::high_resolution_clock::now();
        getuv();
        auto t4 = std::chrono::high_resolution_clock::now();
        std::cout<< "Done." << std::endl;

        std::cout << "getuv() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t4-t3).count()
                  << " milliseconds\n";

        std::cout << "Calling calcu_upr()..." << std::endl;
        auto t5 = std::chrono::high_resolution_clock::now();
        calcu_upr();
        auto t6 = std::chrono::high_resolution_clock::now();
        std::cout<< "Done." << std::endl;

        std::cout << "calcu_upr() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t6-t5).count()
                  << " milliseconds\n";

        std::cout << "Calling calcu_Fxy()..." << std::endl;
        auto t7 = std::chrono::high_resolution_clock::now();
        calcu_Fxy();
        auto t8 = std::chrono::high_resolution_clock::now();
        std::cout<< "Done." << std::endl;

        std::cout << "calcu_Fxy() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t8-t7).count()
                  << " milliseconds\n";              

        std::cout << "Calling collision()..." << std::endl;
        auto t9 = std::chrono::high_resolution_clock::now();
        collision();
        auto t10 = std::chrono::high_resolution_clock::now();
        std::cout<< "Done." << std::endl;

        std::cout << "collision() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t10-t9).count()
                  << " milliseconds\n";   

        int test_time = 250;          

        std::cout << "Calling get_and_write_unsat_results()..." << std::endl;
        auto t11 = std::chrono::high_resolution_clock::now();
        get_and_write_unsat_results(test_time);
        auto t12 = std::chrono::high_resolution_clock::now();
        std::cout<< "Done." << std::endl;

        std::cout << "get_and_write_unsat_results() took "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(t12-t11).count()
                  << " milliseconds\n"; 

    }              

    if (run_simulation ==true){

        std::cout << "===================================" << std::endl;
        std::cout << "Main calculation loop" << std::endl;
        std::cout << "===================================" << std::endl;

        #ifdef _OPENMP
        std::cout << "----------" << std::endl;
        std::cout << "Running simulation in parallel using " << (selected_number_of_threads) << " threads." << std::endl;
        std::cout << "----------" << std::endl;        
        #endif

        // Timing of simulation:
        auto t13 = std::chrono::high_resolution_clock::now();

        for (int time = 1; time < t_max; time++){

            if ((time % Nwri == 0) | (time == 1)){
                write_results(time);
                std::cout << "Results written to file at time = " << time << std::endl;
            }
           
            get_and_write_unsat_results(time);

            std::cout << "***" << std::endl;
            std::cout << "**********************" << std::endl;
            std::cout << "Progress: " << (float(time-1)/float(t_max)*100.0) << " %" << std::endl;
            std::cout << "**********************" << std::endl;            
            std::cout << "***" << std::endl;

            stream();

            getuv();

            calcu_upr();

            calcu_Fxy();

            collision();

            // Hydraulic steps after 5000 time steps:
            if (time==5000){
                step_number = 1;
            }

            // Check of degree of saturation:
            if ((saturation >= 0.98) && (step_number == 1)){
                step_number = 2;
            }

            if ((saturation <= 0.3) && (step_number == 2)){
                step_number = 3;
            }   

            if ((saturation >= 0.857) && (step_number == 3)){
                step_number = 4;
            }   

            if ((saturation <= 0.393) && (step_number == 4)){
                step_number = 5;
            }   
            
            if ((saturation >= 0.764) && (step_number == 5)){
                step_number = 6;
            }       

            if ((saturation <= 0.486) && (step_number == 6)){
                step_number = 7;
            }     

            if ((saturation >= 0.671) && (step_number == 7)){
                step_number = 8;
            }   

            if ((saturation <= 0.579) && (step_number == 8)){
                step_number = 9;
            }                                                                         


            // Hydraulic step 1: Imbibition (step_number=1)
            if ( (time % 1000 == 0) && (step_number == 1)){

                std::cout << "Increasing fluid node density in liquid nodes only by +0.005 after collision step at time = "
                    << time << std::endl; 

            increase_density_in_liquid_nodes_only();

            }   

            // Hydraulic step 2: Drainage (step_number=2)
            if ( (time % 1000 == 0) && (step_number == 2)){
        
                std::cout << "Decreasing fluid node density in liquid nodes only by -0.005 after collision step at time = "
                    << time << std::endl; 

            decrease_density_in_liquid_nodes_only();

            }

            // Hydraulic step 3: Imbibition (step_number=3)
            if ( (time % 1000 == 0) && (step_number == 3)){

                std::cout << "Increasing fluid node density in liquid nodes only by +0.005 after collision step at time = "
                    << time << std::endl; 

            increase_density_in_liquid_nodes_only();

            }

            // Hydraulic step 4: Drainage (step_number=4)
            if ( (time % 1000 == 0) && (step_number == 4)){
        
                std::cout << "Decreasing fluid node density in liquid nodes only by -0.005 after collision step at time = "
                    << time << std::endl; 

            decrease_density_in_liquid_nodes_only();

            }    

            // Hydraulic step 5: Imbibition (step_number=5)
            if ( (time % 1000 == 0) && (step_number == 5)){

                std::cout << "Increasing fluid node density in liquid nodes only by +0.005 after collision step at time = "
                    << time << std::endl; 

            increase_density_in_liquid_nodes_only();

            }

            // Hydraulic step 6: Drainage (step_number=6)
            if ( (time % 1000 == 0) && (step_number == 6)){
        
                std::cout << "Decreasing fluid node density in liquid nodes only by -0.005 after collision step at time = "
                    << time << std::endl; 

            decrease_density_in_liquid_nodes_only();

            }

            // Hydraulic step 7: Imbibition (step_number=7)
            if ( (time % 1000 == 0) && (step_number == 7)){

                std::cout << "Increasing fluid node density in liquid nodes only by +0.005 after collision step at time = "
                    << time << std::endl; 

            increase_density_in_liquid_nodes_only();

            } 

            // Hydraulic step 8: Drainage (step_number=8)
            if ( (time % 1000 == 0) && (step_number == 8)){
        
                std::cout << "Decreasing fluid node density in liquid nodes only by -0.005 after collision step at time = "
                    << time << std::endl; 

            decrease_density_in_liquid_nodes_only();

            }  

            // Step 9: Exit (step_number=9)
            if ( (time % 1000 == 0) && (step_number == 9)){

                std::cout << "*** The maximum step number is reached. ***"
                    << std::endl; 

            break;

            }                                                                       

        }

        auto t14 = std::chrono::high_resolution_clock::now();

        std::cout << "=============================" << std::endl;
        std::cout << "End of main calculation loop." << std::endl;
        std::cout << "=============================" << std::endl;


        std::cout << "Runtime of simulation is "
                  << std::chrono::duration_cast<std::chrono::seconds>(t14-t13).count()
                  << " seconds\n";         

    }

  return 0;  
}