diff --git a/analysis/GreyPhase.cpp b/analysis/GreyPhase.cpp index 74f4437a..011d9b6a 100644 --- a/analysis/GreyPhase.cpp +++ b/analysis/GreyPhase.cpp @@ -19,6 +19,7 @@ GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr dm): Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0); Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0); + MobilityRatio.resize(Nx,Ny,Nz); MobilityRatio.fill(0); //......................................... if (Dm->rank()==0){ @@ -89,14 +90,17 @@ void GreyPhaseAnalysis::Basic(){ double nB = Rho_w(n); double phi = (nA-nB)/(nA+nB); double porosity = Porosity(n); - Water_local.M += rho_w*nB*porosity; - Water_local.Px += porosity*rho_w*nB*Vel_x(n); - Water_local.Py += porosity*rho_w*nB*Vel_y(n); - Water_local.Pz += porosity*rho_w*nB*Vel_z(n); - Oil_local.M += rho_n*nA*porosity; - Oil_local.Px += porosity*rho_n*nA*Vel_x(n); - Oil_local.Py += porosity*rho_n*nA*Vel_y(n); - Oil_local.Pz += porosity*rho_n*nA*Vel_z(n); + double mobility_ratio = MobilityRatio(n); + + Water_local.M += nB*porosity; + Water_local.Px += porosity*(nA+nB)*Vel_x(n)*0.5*(1.0-mobility_ratio); + Water_local.Py += porosity*(nA+nB)*Vel_y(n)*0.5*(1.0-mobility_ratio); + Water_local.Pz += porosity*(nA+nB)*Vel_z(n)*0.5*(1.0-mobility_ratio); + + Oil_local.M += nA*porosity; + Oil_local.Px += porosity*(nA+nB)*Vel_x(n)*0.5*(1.0+mobility_ratio); + Oil_local.Py += porosity*(nA+nB)*Vel_y(n)*0.5*(1.0+mobility_ratio); + Oil_local.Pz += porosity*(nA+nB)*Vel_z(n)*0.5*(1.0+mobility_ratio); if ( phi > 0.99 ){ Oil_local.p += Pressure(n); diff --git a/analysis/GreyPhase.h b/analysis/GreyPhase.h index b2632a70..eb8724fb 100644 --- a/analysis/GreyPhase.h +++ b/analysis/GreyPhase.h @@ -71,6 +71,7 @@ public: DoubleArray Vel_x; // velocity field DoubleArray Vel_y; DoubleArray Vel_z; + DoubleArray MobilityRatio; GreyPhaseAnalysis(std::shared_ptr Dm); ~GreyPhaseAnalysis(); diff --git a/common/Domain.h b/common/Domain.h index 13eacd45..7313587e 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -16,87 +16,98 @@ #include "common/MPI.h" #include "common/Communication.h" #include "common/Database.h" - -class Domain; -template class PatchData; +/** + * @file Domain.h + * \brief Parallel Domain data structures and helper functions + */ -//! Class to hold information about a box +/** + * \class Box + * + * @details + * information about a box + */ class Box { public: int ifirst[3]; int ilast[3]; }; +class Patch; -enum class DataLocation { CPU, DEVICE }; +/** + * \class Domain + * + * @details + * the Domain class includes basic information to distribution 3D image data to multiple processes using MPI. + * A regular domain decomposision is performed, with each MPI process getting a [Nx,Ny,Nz] sub-domain. + * 8-bit image labels are retained internally. + * The domain class resides on the CPU and provides utilities to support CPU-based analysis. + * GPU-based data structures should be constructed separately but may utilize information that the Domain class provides. +*/ - -//! Class to hold information about a patch -class Patch { -public: - - //! Empty constructor - Patch() = delete; - - //! Copy constructor - Patch( const Patch& ) = delete; - - //! Assignment operator - Patch& operator=( const Patch& ) = delete; - - //! Return the box for the patch - inline const Box& getBox() const { return d_box; } - - //! Create patch data - template - std::shared_ptr> createPatchData( DataLocation location ) const; - -private: - Box d_box; - int d_owner; - Domain *d_domain; - -}; - - -//! Class to hold domain info class Domain{ public: - //! Default constructor + /** + * \brief Constructor + * @param db input database + * @param Communicator MPI communicator + */ Domain( std::shared_ptr db, const Utilities::MPI& Communicator); - //! Obsolete constructor + /** + * \brief Obsolete constructor + */ Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, double lx, double ly, double lz, int BC); - //! Empty constructor + /** + * \brief Empty constructor + */ Domain() = delete; - //! Copy constructor + /** + * \brief Copy constructor + */ Domain( const Domain& ) = delete; - //! Assignment operator + /** + * \brief Assignment operator + */ Domain& operator=( const Domain& ) = delete; - //! Destructor + /** + * \brief Destructor + */ ~Domain(); - //! Get the database + /** + * \brief Get the database + */ inline std::shared_ptr getDatabase() const { return d_db; } - //! Get the domain box + /** + * \brief Get the domain box + */ inline const Box& getBox() const { return d_box; } - //! Get local patch + /** + * \brief Get local patch + */ inline const Patch& getLocalPatch() const { return *d_localPatch; } - //! Get all patches + /** + * \brief Get all patches + */ inline const std::vector& getAllPatch() const { return d_patches; } private: + /** + * \brief initialize from database + */ void initialize( std::shared_ptr db ); std::shared_ptr d_db; @@ -124,6 +135,9 @@ public: // Public variables (need to create accessors instead) //********************************** // MPI ranks for all 18 neighbors //********************************** + /** + * \brief Compute the porosity based on the current domain id file + */ inline double Porosity() const { return porosity; } inline int iproc() const { return rank_info.ix; } inline int jproc() const { return rank_info.jy; } @@ -165,22 +179,78 @@ public: // Public variables (need to create accessors instead) // Solid indicator function std::vector id; + /** + * \brief Read domain IDs from file + */ void ReadIDs(); + + /** + * \brief Compute the porosity + */ void ComputePorosity(); + + /** + * \brief Read image and perform domain decomposition + * @param filename - name of file to read IDs + */ void Decomp( const std::string& filename ); + + /** + * \brief Perform a halo exchange using MPI + * @param Mesh - array data that holds scalar values + */ void CommunicateMeshHalo(DoubleArray &Mesh); + + /** + * \brief Initialize communication data structures within Domain object. + * This routine needs to be called before the communication functionality will work + */ void CommInit(); + + /** + * \brief Count number of pore nodes (labels > 1) + */ int PoreCount(); + /** + * \brief Read array data from a file and distribute to local arrays for each MPI process + * @param Filename - name of the file to read the data + * @param Datatype - data type to use + * @param UserData - Array to store the values that are read + */ void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData); + + /** + * \brief Aggregate labels from all MPI processes and write to a file + * @param filename - name of the file to write + */ void AggregateLabels( const std::string& filename ); + /** + * \brief Aggregate user provided array from all MPI processes and write to a single file + * @param filename - name of the file to write + * @param UserData - array data to aggregate and write + */ void AggregateLabels( const std::string& filename, DoubleArray &UserData ); private: + /** + * \brief Pack halo data for 8-bit integer + * @param list - list of values in the halo + * @param count - count of values in the halo + * @param sendbuf - memory buffer to use to pack values for MPI + * @param ID - 8-bit values on mesh [Nx, Ny, Nz] + */ void PackID(int *list, int count, signed char *sendbuf, signed char *ID); + + /** + * \brief Unpack halo data for 8-bit integer + * @param list - list of values in the halo + * @param count - count of values in the halo + * @param recvbuf - memory buffer containing values recieved by MPI + * @param ID - 8-bit values on mesh [Nx, Ny, Nz] + */ void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID); - void CommHaloIDs(); //...................................................................................... MPI_Request req1[18], req2[18]; @@ -198,6 +268,44 @@ private: }; +template class PatchData; + + +enum class DataLocation { CPU, DEVICE }; + + +/** + * \class Patch + * + * @details + * store patch data + */ +class Patch { +public: + + //! Empty constructor + Patch() = delete; + + //! Copy constructor + Patch( const Patch& ) = delete; + + //! Assignment operator + Patch& operator=( const Patch& ) = delete; + + //! Return the box for the patch + inline const Box& getBox() const { return d_box; } + + //! Create patch data + template + std::shared_ptr> createPatchData( DataLocation location ) const; + +private: + Box d_box; + int d_owner; + Domain *d_domain; + +}; + // Class to hold data on a patch template diff --git a/common/ScaLBL.h b/common/ScaLBL.h index ba253be6..816c1551 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -234,12 +234,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, double *Perm,double *Vel, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np); diff --git a/cpu/GreyscaleColor.cpp b/cpu/GreyscaleColor.cpp index 617cd293..332f8e14 100644 --- a/cpu/GreyscaleColor.cpp +++ b/cpu/GreyscaleColor.cpp @@ -1271,7 +1271,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); dist[16*Np+n] = fq; - // q = 17 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) @@ -1341,7 +1340,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl //CP: capillary penalty // also turn off recoloring for grey nodes extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, + double *Velocity, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ @@ -1378,7 +1378,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map /* Corey model parameters */ double Kn_grey,Kw_grey; double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB; - double GreyDiff; // grey diffusion + double GreyDiff=0.0e-4; // grey diffusion const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1399,7 +1399,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map nB = Den[Np + n]; porosity = Poros[n]; - GreyDiff = Perm[n]; + //GreyDiff = Perm[n]; perm = 1.0; W = GreySolidW[n]; Sn_grey = GreySn[n]; @@ -1423,21 +1423,33 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map Krw_grey = 0.0; if (nA/(nA+nB)=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero - Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero + Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero // recompute the effective permeability - perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5)); + perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5)); //mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); } else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){ perm = Kn_grey; + Krn_grey = Kn_grey; Swn = 1.0; - } - mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); + } + /* compute the mobility ratio */ + if (porosity != 1.0){ + mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5)); + } + else if (phi > 0.0){ + mobility_ratio = 1.0; + } + else { + mobility_ratio = -1.0; + } + MobilityRatio[n] = mobility_ratio; // Get the 1D index based on regular data layout ijk = Map[n]; @@ -2092,7 +2104,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio); } @@ -2122,7 +2134,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio); } @@ -2153,7 +2165,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio); } @@ -2181,7 +2193,8 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map //CP: capillary penalty // also turn off recoloring for grey nodes extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, + double *Velocity, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ @@ -2205,7 +2218,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do /* Corey model parameters */ double Kn_grey,Kw_grey; double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB; - double GreyDiff; // grey diffusion + double GreyDiff=0.0e-4; // grey diffusion //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) double porosity; @@ -2235,7 +2248,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do nB = Den[Np + n]; porosity = Poros[n]; - GreyDiff = Perm[n]; + //GreyDiff = Perm[n]; perm = 1.0; W = GreySolidW[n]; Sn_grey = GreySn[n]; @@ -2254,27 +2267,40 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity - + + mobility_ratio = 1.0; Krn_grey = 0.0; Krw_grey = 0.0; if (nA/(nA+nB)=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero - Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero + Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero // recompute the effective permeability - perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5)); + perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5)); //mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); } else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){ perm = Kn_grey; + Krn_grey = Kn_grey; Swn = 1.0; - } - mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); - + } + /* compute the mobility ratio */ + if (porosity != 1.0){ + mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5)); + } + else if (phi > 0.0){ + mobility_ratio = 1.0; + } + else { + mobility_ratio = -1.0; + } + MobilityRatio[n] = mobility_ratio; + // Get the 1D index based on regular data layout ijk = Map[n]; // COMPUTE THE COLOR GRADIENT @@ -2861,7 +2887,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio); } @@ -2887,7 +2913,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio); } @@ -2914,7 +2940,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ //delta = 0.0; - delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio); } @@ -2964,32 +2990,4 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl } } -//extern "C" void ScaLBL_D3Q19_GreyscaleColor_Init(double *dist, double *Porosity, int Np){ -// int n; -// double porosity; -// for (n=0; n=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero - Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero + Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero // recompute the effective permeability - perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5)); - mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); + perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5)); + //mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); } else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){ perm = Kn_grey; + Krn_grey = Kn_grey; Swn = 1.0; - } + } + /* compute the mobility ratio */ + if (porosity != 1.0){ + mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5)); + } + else if (phi > 0.0){ + mobility_ratio = 1.0; + } + else { + mobility_ratio = -1.0; + } + MobilityRatio[n] = mobility_ratio; + // Get the 1D index based on regular data layout ijk = Map[n]; // COMPUTE THE COLOR GRADIENT @@ -1606,27 +1625,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - //............Compute the Greyscale Potential Gradient..................... - // Fcpx = 0.0; - // Fcpy = 0.0; - // Fcpz = 0.0; - // if (porosity!=1.0){ - // //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14)); - // //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18)); - // //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18)); - // Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - // Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - // Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - // Fcpx *= alpha*W/sqrt(perm); - // Fcpy *= alpha*W/sqrt(perm); - // Fcpz *= alpha*W/sqrt(perm); - // //double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz); - // //double Fcp_mag = Fcp_mag_temp; - // //if (Fcp_mag_temp==0.0) Fcp_mag=1.0; - // //nx = Fcpx/Fcp_mag; - // //ny = Fcpy/Fcp_mag; - // //nz = Fcpz/Fcp_mag; - // } Fcpx = nx; Fcpy = ny; Fcpz = nz; @@ -2023,42 +2021,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int m18 = m18 + rlx_setB*( - m18); //----------------------------------------------------------------------// - //----------------With higher-order force ------------------------------// - //if (C == 0.0) nx = ny = nz = 0.0; - //m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1) - // + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity; - //m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2) - // + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity; - //jx = jx + Fx; - //m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); - //jy = jy + Fy; - //m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); - //jz = jz + Fz; - //m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); - //m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9) - // + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity; - ////m10 = m10 + rlx_setA*( - m10); - //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10) - // + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity; - //m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11) - // + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity; - ////m12 = m12 + rlx_setA*( - m12); - //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12) - // + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity; - //m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13); - // + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; - //m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14); - // + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; - //m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15); - // + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity; - //m16 = m16 + rlx_setB*( - m16); - //m17 = m17 + rlx_setB*( - m17); - //m18 = m18 + rlx_setB*( - m18); - //----------------------------------------------------------------------// - //.................inverse transformation...................................................... // q=0 fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; @@ -2194,8 +2156,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; - jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; + jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio); } if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; @@ -2223,7 +2186,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio); } @@ -2253,7 +2217,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio); } @@ -2282,7 +2247,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, + double *Perm, double *Velocity, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -2306,7 +2272,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis double Fcpx,Fcpy,Fcpz;//capillary penalty force double W;//greyscale wetting strength double Sn_grey,Sw_grey; - + double GreyDiff=0.0e-4; + /* Corey model parameters */ double Kn_grey,Kw_grey; double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB; @@ -2329,11 +2296,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis //........Get 1-D index for this thread.................... n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start; if (n=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero - Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero + Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero // recompute the effective permeability - perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5)); - mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); + perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5)); + //mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5)); } else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){ perm = Kn_grey; + Krn_grey = Kn_grey; Swn = 1.0; - } + } + /* compute the mobility ratio */ + if (porosity != 1.0){ + mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5)); + } + else if (phi > 0.0){ + mobility_ratio = 1.0; + } + else { + mobility_ratio = -1.0; + } + MobilityRatio[n] = mobility_ratio; + // Get the 1D index based on regular data layout ijk = Map[n]; // COMPUTE THE COLOR GRADIENT @@ -2433,27 +2418,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - //............Compute the Greyscale Potential Gradient..................... - // Fcpx = 0.0; - // Fcpy = 0.0; - // Fcpz = 0.0; - // if (porosity!=1.0){ - // //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14)); - // //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18)); - // //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18)); - // Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - // Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - // Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - // Fcpx *= alpha*W/sqrt(perm); - // Fcpy *= alpha*W/sqrt(perm); - // Fcpz *= alpha*W/sqrt(perm); - // double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz); - // double Fcp_mag = Fcp_mag_temp; - // if (Fcp_mag_temp==0.0) Fcp_mag=1.0; - // nx = Fcpx/Fcp_mag; - // ny = Fcpy/Fcp_mag; - // nz = Fcpz/Fcp_mag; - // } Fcpx = nx; Fcpy = ny; Fcpz = nz; @@ -2799,42 +2763,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis m18 = m18 + rlx_setB*( - m18); //----------------------------------------------------------------------// - //----------------With higher-order force ------------------------------// - //if (C == 0.0) nx = ny = nz = 0.0; - //m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1) - // + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity; - //m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2) - // + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity; - //jx = jx + Fx; - //m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); - //jy = jy + Fy; - //m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); - //jz = jz + Fz; - //m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8) - // + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); - //m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9) - // + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity; - ////m10 = m10 + rlx_setA*( - m10); - //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10) - // + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity; - //m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11) - // + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity; - ////m12 = m12 + rlx_setA*( - m12); - //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12) - // + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity; - //m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13); - // + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; - //m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14); - // + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; - //m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15); - // + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity; - //m16 = m16 + rlx_setB*( - m16); - //m17 = m17 + rlx_setB*( - m17); - //m18 = m18 + rlx_setB*( - m18); - //----------------------------------------------------------------------// - //.................inverse transformation...................................................... // q=0 fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; @@ -2954,7 +2882,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx; jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio); } @@ -2979,7 +2908,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny; jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio); } @@ -3005,7 +2935,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis if (!(nA*nB*nAB>0)) delta=0; //----------------newly added for better control of recoloring---------------// if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ - delta = 0.0; + //delta = 0.0; + delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz; jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio); } @@ -3023,6 +2954,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis Aq[6*Np+n] = a2; Bq[6*Np+n] = b2; //............................................... + } } } @@ -3059,1508 +2991,6 @@ __global__ void dvc_ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, d } } -//NOTE: so far it seems that we don't need this greyscale potental update; -// if we compute a grey-potential first, and take its gradient to work out the capillary penalty force, it is highly unstable; -// this is because the grey-potential is simply a scaling of the normal phase field, but such scaling create some artificial gradient at the open-grey interface -//__global__ void dvc_ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W, -// int start, int finish, int Np){ -// int idx,n; -// double phi,psi; -// double cap_penalty; -// double porosity,perm; -// -// int S = Np/NBLOCKS/NTHREADS + 1; -// for (int s=0; s1.0) t1 =((t1>0.0)-(t1<0.0))*(1.0-fabs(t1))+t1; -// //........................................................................ -// nn = ijk+1; // neighbor index (get convention) -// m2 = Phi[nn]; // get neighbor for phi - 2 -// t2 = m2+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t2)>1.0) t2 =((t2>0.0)-(t2<0.0))*(1.0-fabs(t2))+t2; -// //........................................................................ -// nn = ijk-strideY; // neighbor index (get convention) -// m3 = Phi[nn]; // get neighbor for phi - 3 -// t3 = m3+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t3)>1.0) t3 =((t3>0.0)-(t3<0.0))*(1.0-fabs(t3))+t3; -// //........................................................................ -// nn = ijk+strideY; // neighbor index (get convention) -// m4 = Phi[nn]; // get neighbor for phi - 4 -// t4 = m4+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t4)>1.0) t4 =((t4>0.0)-(t4<0.0))*(1.0-fabs(t4))+t4; -// //........................................................................ -// nn = ijk-strideZ; // neighbor index (get convention) -// m5 = Phi[nn]; // get neighbor for phi - 5 -// t5 = m5+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t5)>1.0) t5 =((t5>0.0)-(t5<0.0))*(1.0-fabs(t5))+t5; -// //........................................................................ -// nn = ijk+strideZ; // neighbor index (get convention) -// m6 = Phi[nn]; // get neighbor for phi - 6 -// t6 = m6+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t6)>1.0) t6 =((t6>0.0)-(t6<0.0))*(1.0-fabs(t6))+t6; -// //........................................................................ -// nn = ijk-strideY-1; // neighbor index (get convention) -// m7 = Phi[nn]; // get neighbor for phi - 7 -// t7 = m7+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t7)>1.0) t7 =((t7>0.0)-(t7<0.0))*(1.0-fabs(t7))+t7; -// //........................................................................ -// nn = ijk+strideY+1; // neighbor index (get convention) -// m8 = Phi[nn]; // get neighbor for phi - 8 -// t8 = m8+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t8)>1.0) t8 =((t8>0.0)-(t8<0.0))*(1.0-fabs(t8))+t8; -// //........................................................................ -// nn = ijk+strideY-1; // neighbor index (get convention) -// m9 = Phi[nn]; // get neighbor for phi - 9 -// t9 = m9+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t9)>1.0) t9 =((t9>0.0)-(t9<0.0))*(1.0-fabs(t9))+t9; -// //........................................................................ -// nn = ijk-strideY+1; // neighbor index (get convention) -// m10 = Phi[nn]; // get neighbor for phi - 10 -// t10 = m10+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t10)>1.0) t10 =((t10>0.0)-(t10<0.0))*(1.0-fabs(t10))+t10; -// //........................................................................ -// nn = ijk-strideZ-1; // neighbor index (get convention) -// m11 = Phi[nn]; // get neighbor for phi - 11 -// t11 = m11+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t11)>1.0) t11 =((t11>0.0)-(t11<0.0))*(1.0-fabs(t11))+t11; -// //........................................................................ -// nn = ijk+strideZ+1; // neighbor index (get convention) -// m12 = Phi[nn]; // get neighbor for phi - 12 -// t12 = m12+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t12)>1.0) t12 =((t12>0.0)-(t12<0.0))*(1.0-fabs(t12))+t12; -// //........................................................................ -// nn = ijk+strideZ-1; // neighbor index (get convention) -// m13 = Phi[nn]; // get neighbor for phi - 13 -// t13 = m13+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t13)>1.0) t13 =((t13>0.0)-(t13<0.0))*(1.0-fabs(t13))+t13; -// //........................................................................ -// nn = ijk-strideZ+1; // neighbor index (get convention) -// m14 = Phi[nn]; // get neighbor for phi - 14 -// t14 = m14+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t14)>1.0) t14 =((t14>0.0)-(t14<0.0))*(1.0-fabs(t14))+t14; -// //........................................................................ -// nn = ijk-strideZ-strideY; // neighbor index (get convention) -// m15 = Phi[nn]; // get neighbor for phi - 15 -// t15 = m15+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t15)>1.0) t15 =((t15>0.0)-(t15<0.0))*(1.0-fabs(t15))+t15; -// //........................................................................ -// nn = ijk+strideZ+strideY; // neighbor index (get convention) -// m16 = Phi[nn]; // get neighbor for phi - 16 -// t16 = m16+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t16)>1.0) t16 =((t16>0.0)-(t16<0.0))*(1.0-fabs(t16))+t16; -// //........................................................................ -// nn = ijk+strideZ-strideY; // neighbor index (get convention) -// m17 = Phi[nn]; // get neighbor for phi - 17 -// t17 = m17+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t17)>1.0) t17 =((t17>0.0)-(t17<0.0))*(1.0-fabs(t17))+t17; -// //........................................................................ -// nn = ijk-strideZ+strideY; // neighbor index (get convention) -// m18 = Phi[nn]; // get neighbor for phi - 18 -// t18 = m18+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t18)>1.0) t18 =((t18>0.0)-(t18<0.0))*(1.0-fabs(t18))+t18; -// //............Compute the Color Gradient................................... -// nx_phase = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); -// ny_phase = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); -// nz_phase = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); -// C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase); -// //correct the normal color gradient by considering the effect of grey solid -// nx = -(t1-t2+0.5*(t7-t8+t9-t10+t11-t12+t13-t14)); -// ny = -(t3-t4+0.5*(t7-t8-t9+t10+t15-t16+t17-t18)); -// nz = -(t5-t6+0.5*(t11-t12-t13+t14+t15-t16-t17+t18)); -// -// if (C_phase==0.0){//i.e. if in a bulk phase, there is no need for grey-solid correction -// nx = nx_phase; -// ny = ny_phase; -// nz = nz_phase; -// } -// -// //...........Normalize the Color Gradient................................. -// C = sqrt(nx*nx+ny*ny+nz*nz); -// double ColorMag = C; -// if (C==0.0) ColorMag=1.0; -// nx = nx/ColorMag; -// ny = ny/ColorMag; -// nz = nz/ColorMag; -// -// // q=0 -// fq = dist[n]; -// rho = fq; -// m1 = -30.0*fq; -// m2 = 12.0*fq; -// -// // q=1 -// //nread = neighborList[n]; // neighbor 2 -// //fq = dist[nread]; // reading the f1 data into register fq -// nr1 = neighborList[n]; -// fq = dist[nr1]; // reading the f1 data into register fq -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jx = fq; -// m4 = -4.0*fq; -// m9 = 2.0*fq; -// m10 = -4.0*fq; -// -// // f2 = dist[10*Np+n]; -// //nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) -// //fq = dist[nread]; // reading the f2 data into register fq -// nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) -// fq = dist[nr2]; // reading the f2 data into register fq -// rho += fq; -// m1 -= 11.0*(fq); -// m2 -= 4.0*(fq); -// jx -= fq; -// m4 += 4.0*(fq); -// m9 += 2.0*(fq); -// m10 -= 4.0*(fq); -// -// // q=3 -// //nread = neighborList[n+2*Np]; // neighbor 4 -// //fq = dist[nread]; -// nr3 = neighborList[n+2*Np]; // neighbor 4 -// fq = dist[nr3]; -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy = fq; -// m6 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 = fq; -// m12 = -2.0*fq; -// -// // q = 4 -// //nread = neighborList[n+3*Np]; // neighbor 3 -// //fq = dist[nread]; -// nr4 = neighborList[n+3*Np]; // neighbor 3 -// fq = dist[nr4]; -// rho+= fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy -= fq; -// m6 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 += fq; -// m12 -= 2.0*fq; -// -// // q=5 -// //nread = neighborList[n+4*Np]; -// //fq = dist[nread]; -// nr5 = neighborList[n+4*Np]; -// fq = dist[nr5]; -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz = fq; -// m8 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// -// // q = 6 -// //nread = neighborList[n+5*Np]; -// //fq = dist[nread]; -// nr6 = neighborList[n+5*Np]; -// fq = dist[nr6]; -// rho+= fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz -= fq; -// m8 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q=7 -// //nread = neighborList[n+6*Np]; -// //fq = dist[nread]; -// nr7 = neighborList[n+6*Np]; -// fq = dist[nr7]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 = fq; -// m16 = fq; -// m17 = -fq; -// -// // q = 8 -// //nread = neighborList[n+7*Np]; -// //fq = dist[nread]; -// nr8 = neighborList[n+7*Np]; -// fq = dist[nr8]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 += fq; -// m16 -= fq; -// m17 += fq; -// -// // q=9 -// //nread = neighborList[n+8*Np]; -// //fq = dist[nread]; -// nr9 = neighborList[n+8*Np]; -// fq = dist[nr9]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 += fq; -// m17 += fq; -// -// // q = 10 -// //nread = neighborList[n+9*Np]; -// //fq = dist[nread]; -// nr10 = neighborList[n+9*Np]; -// fq = dist[nr10]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 -= fq; -// m17 -= fq; -// -// // q=11 -// //nread = neighborList[n+10*Np]; -// //fq = dist[nread]; -// nr11 = neighborList[n+10*Np]; -// fq = dist[nr11]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 = fq; -// m16 -= fq; -// m18 = fq; -// -// // q=12 -// //nread = neighborList[n+11*Np]; -// //fq = dist[nread]; -// nr12 = neighborList[n+11*Np]; -// fq = dist[nr12]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 += fq; -// m16 += fq; -// m18 -= fq; -// -// // q=13 -// //nread = neighborList[n+12*Np]; -// //fq = dist[nread]; -// nr13 = neighborList[n+12*Np]; -// fq = dist[nr13]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 -= fq; -// m18 -= fq; -// -// // q=14 -// //nread = neighborList[n+13*Np]; -// //fq = dist[nread]; -// nr14 = neighborList[n+13*Np]; -// fq = dist[nr14]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 += fq; -// m18 += fq; -// -// // q=15 -// nread = neighborList[n+14*Np]; -// fq = dist[nread]; -// //fq = dist[17*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 = fq; -// m17 += fq; -// m18 -= fq; -// -// // q=16 -// nread = neighborList[n+15*Np]; -// fq = dist[nread]; -// //fq = dist[8*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 += fq; -// m17 -= fq; -// m18 += fq; -// -// // q=17 -// //fq = dist[18*Np+n]; -// nread = neighborList[n+16*Np]; -// fq = dist[nread]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 += fq; -// m18 += fq; -// -// // q=18 -// nread = neighborList[n+17*Np]; -// fq = dist[nread]; -// //fq = dist[9*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 -= fq; -// m18 -= fq; -// -// // Compute greyscale related parameters -// c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm); -// if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes -// //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity); -// c1 = porosity*0.5*GeoFun/sqrt(perm); -// if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes -// -// vx = jx/rho0+0.5*(porosity*Gx); -// vy = jy/rho0+0.5*(porosity*Gy); -// vz = jz/rho0+0.5*(porosity*Gz); -// v_mag=sqrt(vx*vx+vy*vy+vz*vz); -// ux = vx/(c0+sqrt(c0*c0+c1*v_mag)); -// uy = vy/(c0+sqrt(c0*c0+c1*v_mag)); -// uz = vz/(c0+sqrt(c0*c0+c1*v_mag)); -// u_mag=sqrt(ux*ux+uy*uy+uz*uz); -// -// //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium -// Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx); -// Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy); -// Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz); -// if (porosity==1.0){ -// Fx=rho0*(Gx); -// Fy=rho0*(Gy); -// Fz=rho0*(Gz); -// } -// -// // write the velocity -// Velocity[n] = ux; -// Velocity[Np+n] = uy; -// Velocity[2*Np+n] = uz; -// -// //........................................................................ -// //..............carry out relaxation process.............................. -// //..........Toelke, Fruediger et. al. 2006................................ -// if (C == 0.0) nx = ny = nz = 0.0; -// m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1); -// m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2); -// jx = jx + Fx; -// m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); -// jy = jy + Fy; -// m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); -// jz = jz + Fz; -// m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); -// m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); -// m10 = m10 + rlx_setA*( - m10); -// //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10); -// m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); -// m12 = m12 + rlx_setA*( - m12); -// //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12); -// m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13); -// m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14); -// m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15); -// m16 = m16 + rlx_setB*( - m16); -// m17 = m17 + rlx_setB*( - m17); -// m18 = m18 + rlx_setB*( - m18); -// -// //.................inverse transformation...................................................... -// // q=0 -// fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; -// dist[n] = fq; -// -// // q = 1 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); -// //nread = neighborList[n+Np]; -// dist[nr2] = fq; -// -// // q=2 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); -// //nread = neighborList[n]; -// dist[nr1] = fq; -// -// // q = 3 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// //nread = neighborList[n+3*Np]; -// dist[nr4] = fq; -// -// // q = 4 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// //nread = neighborList[n+2*Np]; -// dist[nr3] = fq; -// -// // q = 5 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// //nread = neighborList[n+5*Np]; -// dist[nr6] = fq; -// -// // q = 6 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// //nread = neighborList[n+4*Np]; -// dist[nr5] = fq; -// -// // q = 7 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); -// //nread = neighborList[n+7*Np]; -// dist[nr8] = fq; -// -// // q = 8 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 -// +mrt_V12*m12+0.25*m13+0.125*(m17-m16); -// //nread = neighborList[n+6*Np]; -// dist[nr7] = fq; -// -// // q = 9 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); -// //nread = neighborList[n+9*Np]; -// dist[nr10] = fq; -// -// // q = 10 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); -// //nread = neighborList[n+8*Np]; -// dist[nr9] = fq; -// -// // q = 11 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12+0.25*m15+0.125*(m18-m16); -// //nread = neighborList[n+11*Np]; -// dist[nr12] = fq; -// -// // q = 12 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ -// mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); -// //nread = neighborList[n+10*Np]; -// dist[nr11]= fq; -// -// // q = 13 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12-0.25*m15-0.125*(m16+m18); -// //nread = neighborList[n+13*Np]; -// dist[nr14] = fq; -// -// // q= 14 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12-0.25*m15+0.125*(m16+m18); -// //nread = neighborList[n+12*Np]; -// dist[nr13] = fq; -// -// -// // q = 15 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) -// -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); -// nread = neighborList[n+15*Np]; -// dist[nread] = fq; -// -// // q = 16 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) -// -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); -// nread = neighborList[n+14*Np]; -// dist[nread] = fq; -// -// -// // q = 17 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) -// -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); -// nread = neighborList[n+17*Np]; -// dist[nread] = fq; -// -// // q = 18 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) -// -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); -// nread = neighborList[n+16*Np]; -// dist[nread] = fq; -// //........................................................................ -// -// // Instantiate mass transport distributions -// // Stationary value - distribution 0 -// nAB = 1.0/(nA+nB); -// Aq[n] = 0.3333333333333333*nA; -// Bq[n] = 0.3333333333333333*nB; -// -// //............................................... -// // q = 0,2,4 -// // Cq = {1,0,0}, {0,1,0}, {0,0,1} -// delta = beta*nA*nB*nAB*0.1111111111111111*nx; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; -// -// // q = 1 -// //nread = neighborList[n+Np]; -// Aq[nr2] = a1; -// Bq[nr2] = b1; -// // q=2 -// //nread = neighborList[n]; -// Aq[nr1] = a2; -// Bq[nr1] = b2; -// -// //............................................... -// // Cq = {0,1,0} -// delta = beta*nA*nB*nAB*0.1111111111111111*ny; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; -// -// // q = 3 -// //nread = neighborList[n+3*Np]; -// Aq[nr4] = a1; -// Bq[nr4] = b1; -// // q = 4 -// //nread = neighborList[n+2*Np]; -// Aq[nr3] = a2; -// Bq[nr3] = b2; -// -// //............................................... -// // q = 4 -// // Cq = {0,0,1} -// delta = beta*nA*nB*nAB*0.1111111111111111*nz; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; -// -// // q = 5 -// //nread = neighborList[n+5*Np]; -// Aq[nr6] = a1; -// Bq[nr6] = b1; -// // q = 6 -// //nread = neighborList[n+4*Np]; -// Aq[nr5] = a2; -// Bq[nr5] = b2; -// //............................................... -// } -// } -//} -// -////Model-2&3 -//__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, -// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, -// double Gx, double Gy, double Gz, int strideY, int strideZ, int start, int finish, int Np){ -// int ijk,nn,n; -// double fq; -// // conserved momemnts -// double rho,jx,jy,jz; -// double vx,vy,vz,v_mag; -// double ux,uy,uz,u_mag; -// // non-conserved moments -// double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; -// double m3,m5,m7; -// double t1,t2,t4,t6,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18; -// double t3,t5,t7; -// double nA,nB; // number density -// double a1,b1,a2,b2,nAB,delta; -// double C,nx,ny,nz; //color gradient magnitude and direction -// double phi,tau,rho0,rlx_setA,rlx_setB; -// -// double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) -// double porosity; -// double perm;//voxel permeability -// double c0, c1; //Guo's model parameters -// double tau_eff; -// double mu_eff;//kinematic viscosity -// double nx_phase,ny_phase,nz_phase,C_phase; -// double Fx,Fy,Fz; -// -// const double mrt_V1=0.05263157894736842; -// const double mrt_V2=0.012531328320802; -// const double mrt_V3=0.04761904761904762; -// const double mrt_V4=0.004594820384294068; -// const double mrt_V5=0.01587301587301587; -// const double mrt_V6=0.0555555555555555555555555; -// const double mrt_V7=0.02777777777777778; -// const double mrt_V8=0.08333333333333333; -// const double mrt_V9=0.003341687552213868; -// const double mrt_V10=0.003968253968253968; -// const double mrt_V11=0.01388888888888889; -// const double mrt_V12=0.04166666666666666; -// -// int S = Np/NBLOCKS/NTHREADS + 1; -// for (int s=0; s1.0) t1 =((t1>0.0)-(t1<0.0))*(1.0-fabs(t1))+t1; -// //........................................................................ -// nn = ijk+1; // neighbor index (get convention) -// m2 = Phi[nn]; // get neighbor for phi - 2 -// t2 = m2+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t2)>1.0) t2 =((t2>0.0)-(t2<0.0))*(1.0-fabs(t2))+t2; -// //........................................................................ -// nn = ijk-strideY; // neighbor index (get convention) -// m3 = Phi[nn]; // get neighbor for phi - 3 -// t3 = m3+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t3)>1.0) t3 =((t3>0.0)-(t3<0.0))*(1.0-fabs(t3))+t3; -// //........................................................................ -// nn = ijk+strideY; // neighbor index (get convention) -// m4 = Phi[nn]; // get neighbor for phi - 4 -// t4 = m4+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t4)>1.0) t4 =((t4>0.0)-(t4<0.0))*(1.0-fabs(t4))+t4; -// //........................................................................ -// nn = ijk-strideZ; // neighbor index (get convention) -// m5 = Phi[nn]; // get neighbor for phi - 5 -// t5 = m5+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t5)>1.0) t5 =((t5>0.0)-(t5<0.0))*(1.0-fabs(t5))+t5; -// //........................................................................ -// nn = ijk+strideZ; // neighbor index (get convention) -// m6 = Phi[nn]; // get neighbor for phi - 6 -// t6 = m6+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t6)>1.0) t6 =((t6>0.0)-(t6<0.0))*(1.0-fabs(t6))+t6; -// //........................................................................ -// nn = ijk-strideY-1; // neighbor index (get convention) -// m7 = Phi[nn]; // get neighbor for phi - 7 -// t7 = m7+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t7)>1.0) t7 =((t7>0.0)-(t7<0.0))*(1.0-fabs(t7))+t7; -// //........................................................................ -// nn = ijk+strideY+1; // neighbor index (get convention) -// m8 = Phi[nn]; // get neighbor for phi - 8 -// t8 = m8+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t8)>1.0) t8 =((t8>0.0)-(t8<0.0))*(1.0-fabs(t8))+t8; -// //........................................................................ -// nn = ijk+strideY-1; // neighbor index (get convention) -// m9 = Phi[nn]; // get neighbor for phi - 9 -// t9 = m9+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t9)>1.0) t9 =((t9>0.0)-(t9<0.0))*(1.0-fabs(t9))+t9; -// //........................................................................ -// nn = ijk-strideY+1; // neighbor index (get convention) -// m10 = Phi[nn]; // get neighbor for phi - 10 -// t10 = m10+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t10)>1.0) t10 =((t10>0.0)-(t10<0.0))*(1.0-fabs(t10))+t10; -// //........................................................................ -// nn = ijk-strideZ-1; // neighbor index (get convention) -// m11 = Phi[nn]; // get neighbor for phi - 11 -// t11 = m11+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t11)>1.0) t11 =((t11>0.0)-(t11<0.0))*(1.0-fabs(t11))+t11; -// //........................................................................ -// nn = ijk+strideZ+1; // neighbor index (get convention) -// m12 = Phi[nn]; // get neighbor for phi - 12 -// t12 = m12+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t12)>1.0) t12 =((t12>0.0)-(t12<0.0))*(1.0-fabs(t12))+t12; -// //........................................................................ -// nn = ijk+strideZ-1; // neighbor index (get convention) -// m13 = Phi[nn]; // get neighbor for phi - 13 -// t13 = m13+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t13)>1.0) t13 =((t13>0.0)-(t13<0.0))*(1.0-fabs(t13))+t13; -// //........................................................................ -// nn = ijk-strideZ+1; // neighbor index (get convention) -// m14 = Phi[nn]; // get neighbor for phi - 14 -// t14 = m14+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t14)>1.0) t14 =((t14>0.0)-(t14<0.0))*(1.0-fabs(t14))+t14; -// //........................................................................ -// nn = ijk-strideZ-strideY; // neighbor index (get convention) -// m15 = Phi[nn]; // get neighbor for phi - 15 -// t15 = m15+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t15)>1.0) t15 =((t15>0.0)-(t15<0.0))*(1.0-fabs(t15))+t15; -// //........................................................................ -// nn = ijk+strideZ+strideY; // neighbor index (get convention) -// m16 = Phi[nn]; // get neighbor for phi - 16 -// t16 = m16+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t16)>1.0) t16 =((t16>0.0)-(t16<0.0))*(1.0-fabs(t16))+t16; -// //........................................................................ -// nn = ijk+strideZ-strideY; // neighbor index (get convention) -// m17 = Phi[nn]; // get neighbor for phi - 17 -// t17 = m17+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t17)>1.0) t17 =((t17>0.0)-(t17<0.0))*(1.0-fabs(t17))+t17; -// //........................................................................ -// nn = ijk-strideZ+strideY; // neighbor index (get convention) -// m18 = Phi[nn]; // get neighbor for phi - 18 -// t18 = m18+(1.0-porosity)*GreySolidGrad[nn]; -// if (fabs(t18)>1.0) t18 =((t18>0.0)-(t18<0.0))*(1.0-fabs(t18))+t18; -// //............Compute the Color Gradient................................... -// nx_phase = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); -// ny_phase = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); -// nz_phase = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); -// C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase); -// //correct the normal color gradient by considering the effect of grey solid -// nx = -(t1-t2+0.5*(t7-t8+t9-t10+t11-t12+t13-t14)); -// ny = -(t3-t4+0.5*(t7-t8-t9+t10+t15-t16+t17-t18)); -// nz = -(t5-t6+0.5*(t11-t12-t13+t14+t15-t16-t17+t18)); -// -// if (C_phase==0.0){ -// nx = nx_phase; -// ny = ny_phase; -// nz = nz_phase; -// } -// -// //...........Normalize the Color Gradient................................. -// C = sqrt(nx*nx+ny*ny+nz*nz); -// double ColorMag = C; -// if (C==0.0) ColorMag=1.0; -// nx = nx/ColorMag; -// ny = ny/ColorMag; -// nz = nz/ColorMag; -// -// // q=0 -// fq = dist[n]; -// rho = fq; -// m1 = -30.0*fq; -// m2 = 12.0*fq; -// -// // q=1 -// fq = dist[2*Np+n]; -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jx = fq; -// m4 = -4.0*fq; -// m9 = 2.0*fq; -// m10 = -4.0*fq; -// -// // f2 = dist[10*Np+n]; -// fq = dist[1*Np+n]; -// rho += fq; -// m1 -= 11.0*(fq); -// m2 -= 4.0*(fq); -// jx -= fq; -// m4 += 4.0*(fq); -// m9 += 2.0*(fq); -// m10 -= 4.0*(fq); -// -// // q=3 -// fq = dist[4*Np+n]; -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy = fq; -// m6 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 = fq; -// m12 = -2.0*fq; -// -// // q = 4 -// fq = dist[3*Np+n]; -// rho+= fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy -= fq; -// m6 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 += fq; -// m12 -= 2.0*fq; -// -// // q=5 -// fq = dist[6*Np+n]; -// rho += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz = fq; -// m8 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q = 6 -// fq = dist[5*Np+n]; -// rho+= fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz -= fq; -// m8 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q=7 -// fq = dist[8*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 = fq; -// m16 = fq; -// m17 = -fq; -// -// // q = 8 -// fq = dist[7*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 += fq; -// m16 -= fq; -// m17 += fq; -// -// // q=9 -// fq = dist[10*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 += fq; -// m17 += fq; -// -// // q = 10 -// fq = dist[9*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 -= fq; -// m17 -= fq; -// -// // q=11 -// fq = dist[12*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 = fq; -// m16 -= fq; -// m18 = fq; -// -// // q=12 -// fq = dist[11*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 += fq; -// m16 += fq; -// m18 -= fq; -// -// // q=13 -// fq = dist[14*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 -= fq; -// m18 -= fq; -// -// // q=14 -// fq = dist[13*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 += fq; -// m18 += fq; -// -// // q=15 -// fq = dist[16*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 = fq; -// m17 += fq; -// m18 -= fq; -// -// // q=16 -// fq = dist[15*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 += fq; -// m17 -= fq; -// m18 += fq; -// -// // q=17 -// fq = dist[18*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 += fq; -// m18 += fq; -// -// // q=18 -// fq = dist[17*Np+n]; -// rho += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 -= fq; -// m18 -= fq; -// -// // Compute greyscale related parameters -// c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm); -// if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes -// //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity); -// c1 = porosity*0.5*GeoFun/sqrt(perm); -// if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes -// -// vx = jx/rho0+0.5*(porosity*Gx); -// vy = jy/rho0+0.5*(porosity*Gy); -// vz = jz/rho0+0.5*(porosity*Gz); -// v_mag=sqrt(vx*vx+vy*vy+vz*vz); -// ux = vx/(c0+sqrt(c0*c0+c1*v_mag)); -// uy = vy/(c0+sqrt(c0*c0+c1*v_mag)); -// uz = vz/(c0+sqrt(c0*c0+c1*v_mag)); -// u_mag=sqrt(ux*ux+uy*uy+uz*uz); -// -// //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium -// Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx); -// Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy); -// Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz); -// if (porosity==1.0){ -// Fx=rho0*(Gx); -// Fy=rho0*(Gy); -// Fz=rho0*(Gz); -// } -// -// // write the velocity -// Velocity[n] = ux; -// Velocity[Np+n] = uy; -// Velocity[2*Np+n] = uz; -// -// //........................................................................ -// //..............carry out relaxation process.............................. -// //..........Toelke, Fruediger et. al. 2006................................ -// if (C == 0.0) nx = ny = nz = 0.0; -// m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1); -// m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2); -// jx = jx + Fx; -// m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); -// jy = jy + Fy; -// m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); -// jz = jz + Fz; -// m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); -// m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); -// m10 = m10 + rlx_setA*( - m10); -// //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10); -// m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); -// m12 = m12 + rlx_setA*( - m12); -// //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12); -// m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13); -// m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14); -// m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15); -// m16 = m16 + rlx_setB*( - m16); -// m17 = m17 + rlx_setB*( - m17); -// m18 = m18 + rlx_setB*( - m18); -// -// //.................inverse transformation...................................................... -// // q=0 -// fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; -// dist[n] = fq; -// -// // q = 1 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); -// dist[1*Np+n] = fq; -// -// // q=2 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); -// dist[2*Np+n] = fq; -// -// // q = 3 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// dist[3*Np+n] = fq; -// -// // q = 4 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// dist[4*Np+n] = fq; -// -// // q = 5 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// dist[5*Np+n] = fq; -// -// // q = 6 -// fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// dist[6*Np+n] = fq; -// -// // q = 7 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); -// dist[7*Np+n] = fq; -// -// -// // q = 8 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 -// +mrt_V12*m12+0.25*m13+0.125*(m17-m16); -// dist[8*Np+n] = fq; -// -// // q = 9 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); -// dist[9*Np+n] = fq; -// -// // q = 10 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ -// mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); -// dist[10*Np+n] = fq; -// -// -// // q = 11 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12+0.25*m15+0.125*(m18-m16); -// dist[11*Np+n] = fq; -// -// // q = 12 -// fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ -// mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); -// dist[12*Np+n] = fq; -// -// // q = 13 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12-0.25*m15-0.125*(m16+m18); -// dist[13*Np+n] = fq; -// -// // q= 14 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) -// +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -// -mrt_V12*m12-0.25*m15+0.125*(m16+m18); -// -// dist[14*Np+n] = fq; -// -// // q = 15 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) -// -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); -// dist[15*Np+n] = fq; -// -// // q = 16 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) -// -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); -// dist[16*Np+n] = fq; -// -// -// // q = 17 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) -// -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); -// dist[17*Np+n] = fq; -// -// // q = 18 -// fq = mrt_V1*rho+mrt_V9*m1 -// +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) -// -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); -// dist[18*Np+n] = fq; -// //........................................................................ -// -// // Instantiate mass transport distributions -// // Stationary value - distribution 0 -// nAB = 1.0/(nA+nB); -// Aq[n] = 0.3333333333333333*nA; -// Bq[n] = 0.3333333333333333*nB; -// -// //............................................... -// // q = 0,2,4 -// // Cq = {1,0,0}, {0,1,0}, {0,0,1} -// delta = beta*nA*nB*nAB*0.1111111111111111*nx; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; -// -// Aq[1*Np+n] = a1; -// Bq[1*Np+n] = b1; -// Aq[2*Np+n] = a2; -// Bq[2*Np+n] = b2; -// -// //............................................... -// // q = 2 -// // Cq = {0,1,0} -// delta = beta*nA*nB*nAB*0.1111111111111111*ny; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; -// -// Aq[3*Np+n] = a1; -// Bq[3*Np+n] = b1; -// Aq[4*Np+n] = a2; -// Bq[4*Np+n] = b2; -// //............................................... -// // q = 4 -// // Cq = {0,0,1} -// delta = beta*nA*nB*nAB*0.1111111111111111*nz; -// if (!(nA*nB*nAB>0)) delta=0; -// a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; -// b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; -// a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; -// b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; -// -// Aq[5*Np+n] = a1; -// Bq[5*Np+n] = b1; -// Aq[6*Np+n] = a2; -// Bq[6*Np+n] = b2; -// //............................................... -// -// } -// } -//} - -//__global__ void dvc_ScaLBL_D3Q19_GreyscaleColor_Init(double *dist, double *Porosity, int Np) -//{ -// int n; -// int S = Np/NBLOCKS/NTHREADS + 1; -// double porosity; -// for (int s=0; s>>(dist,Porosity,Np); -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q19_GreyscaleColor_Init: %s \n",cudaGetErrorString(err)); -// } -//} - //Model-1 & 4 extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, double *Pressure, @@ -4610,11 +3040,13 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, + double *Perm,double *Vel,double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, GreyKn, GreyKw, Poros, Perm, Vel, Pressure, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, + GreyKn, GreyKw, Poros, Perm, Vel, MobilityRatio, Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); @@ -4626,11 +3058,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, + double *Perm,double *Vel, double *MobilityRatio, double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, GreyKn, GreyKw, Poros, Perm,Vel,Pressure, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, + GreySw, GreyKn, GreyKw, Poros, Perm, Vel, MobilityRatio, Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); @@ -4638,52 +3072,3 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *M printf("CUDA error in ScaLBL_D3Q19_AAodd_GreyscaleColor_CP: %s \n",cudaGetErrorString(err)); } } - -//extern "C" void ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W, -// int start, int finish, int Np){ -// -// dvc_ScaLBL_Update_GreyscalePotential<<>>(Map, Phi, Psi, Poro, Perm, alpha, W, start, finish, Np); -// -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_Update_GreyscalePotential: %s \n",cudaGetErrorString(err)); -// } -//} - -////Model-2&3 -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -// -// //cudaProfilerStart(); -// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor, cudaFuncCachePreferL1); -// -// dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel, -// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColor: %s \n",cudaGetErrorString(err)); -// } -// //cudaProfilerStop(); -// -//} -// -////Model-2&3 -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -// -// //cudaProfilerStart(); -// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor, cudaFuncCachePreferL1); -// -// dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm,Vel, -// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); -// -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q19_AAodd_GreyscaleColor: %s \n",cudaGetErrorString(err)); -// } -// //cudaProfilerStop(); -//} diff --git a/docs/source/developerGuide/doxygen/Domain.rst b/docs/source/developerGuide/doxygen/Domain.rst new file mode 100644 index 00000000..0a1b0aac --- /dev/null +++ b/docs/source/developerGuide/doxygen/Domain.rst @@ -0,0 +1,7 @@ +============================================ +Domain +============================================ +.. doxygenfile:: Domain.h + :project: LBPM Doxygen + + diff --git a/docs/source/developerGuide/index.rst b/docs/source/developerGuide/index.rst index 0e5cbe11..2f373f21 100644 --- a/docs/source/developerGuide/index.rst +++ b/docs/source/developerGuide/index.rst @@ -17,6 +17,8 @@ into the framework. doxygen/ScaLBL.rst + doxygen/Domain.rst + doxygen/Models.rst doxygen/Analysis.rst diff --git a/docs/source/userGuide/models/greyscale/greyscale.rst b/docs/source/userGuide/models/greyscale/greyscale.rst index 47c0621f..52e5486e 100644 --- a/docs/source/userGuide/models/greyscale/greyscale.rst +++ b/docs/source/userGuide/models/greyscale/greyscale.rst @@ -1,7 +1,90 @@ ============================================= -Greyscale model model +Greyscale model ============================================= The LBPM greyscale lattice Boltzmann model is constructed to approximate the solution of the Darcy-Brinkman equations in grey regions, coupled to a Navier-Stokes -solution in open regions. +solution in open regions. To use the greyscale model, the input image should be segmented +into "grey" and open regions. Each particular "grey" label should be assigned both +a porosity and permeability value. + +A typical command to launch the LBPM color simulator is as follows + +``` +mpirun -np $NUMPROCS lbpm_greyscale_simulator input.db +``` + +where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is +the name of the input database that provides the simulation parameters. +Note that the specific syntax to launch MPI tasks may vary depending on your system. +For additional details please refer to your local system documentation. + +*************************** +Model parameters +*************************** + +The essential model parameters for the single phase greyscale model are + +- ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5` + +The kinematic viscosity is given by + +*************************** +Model formulation +*************************** + +A D3Q19 LBE is constructed to describe the momentum transport + +.. math:: + :nowrap: + + $$ + f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = + \sum^{Q-1}_{k=0} M^{-1}_{qk} S_{kk} (m_k^{eq}-m_k) + \sum^{Q-1}_{k=0} M^{-1}_{qk} (1-\frac{S_{kk}}{2}) \hat{F}_q\;, + $$ + + +The force is imposed based on the construction developed by Guo et al + +.. math:: + :nowrap: + + $$ + F_i = \rho_0 \omega_i \left[\frac{\bm{e}_i \cdot \bm{a}}{c_s^2} + + \frac{\bm{u} \bm{a}:(\bm{e}_i \bm{e}_i -c_s^2 \mathcal{I})}{\epsilon c_s^4} \right] , + $$ + + +where :math:`c_s^2 = 1/3` is the speed of sound for the D3Q19 lattice. +The acceleration includes contributions due to the external driving force :math:`\bm{g}` +as well as a drag force due to the permeability :math:`K` and flow rate :math:`\bm{u}` with the +porosity :math:`\epsilon` and viscosity :math:`\nu` determining the net forces acting within +a grey voxel + +.. math:: + :nowrap: + + $$ + \bm{a} = - \frac{\epsilon \nu}{K} \bm{u} + \bm{F}_{cp}/\rho_0 + \epsilon \bm{g}, + $$ + +The flow velocity is defined as + +.. math:: + :nowrap: + + $$ + \rho_0 \bm{u} = \sum_i \bm{e}_i f_i + \frac{\delta t}{2} \rho_0 \bm{a}. + $$ + +Combining the previous expressions, + +.. math:: + :nowrap: + + $$ + \bm{u} = \frac{\frac{1}{\rho_0}\sum_i \bm{e}_i f_i + \frac{\delta t}{2}\epsilon \bm{g} + + \frac{\delta t}{2} \frac{\bm{F}_{cp}}{\rho_0}}{1+ \frac{\delta t}{2} \frac{\epsilon \nu}{K}} + $$ + + diff --git a/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst b/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst new file mode 100644 index 00000000..8189ac29 --- /dev/null +++ b/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst @@ -0,0 +1,65 @@ +============================================= +Greyscale color model +============================================= + +The LBPM greyscale lattice Boltzmann model is constructed to approximate the +solution of the Darcy-Brinkman equations in grey regions coupled to a color model implementation +solution in open regions. A simple constitutive form is used to model the relative +permeability in the grey regions, + + + +A D3Q19 LBE is constructed to describe the momentum transport + +.. math:: + :nowrap: + + $$ + f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = + \sum^{Q-1}_{k=0} M^{-1}_{qk} S_{kk} (m_k^{eq}-m_k) + \sum^{Q-1}_{k=0} M^{-1}_{qk} (1-\frac{S_{kk}}{2}) \hat{F}_q\;, + $$ + + +The force is imposed based on the construction developed by Guo et al + +.. math:: + :nowrap: + + $$ + F_i = \rho_0 \omega_i \left[\frac{\bm{e}_i \cdot \bm{a}}{c_s^2} + + \frac{\bm{u} \bm{a}:(\bm{e}_i \bm{e}_i -c_s^2 \mathcal{I})}{\epsilon c_s^4} \right] , + $$ + + +The acceleration includes contributions due to the external driving force :math:`\bm{g}` +as well as a drag force due to the permeability :math:`K` and flow rate :math:`\bm{u}` with the +porosity :math:`\epsilon` and viscosity :math:`\nu` determining the net forces acting within +a grey voxel + +.. math:: + :nowrap: + + $$ + \bm{a} = - \frac{\epsilon \nu}{K} \bm{u} + \bm{F}_{cp}/\rho_0 + \epsilon \bm{g}, + $$ + +The flow velocity is defined as + +.. math:: + :nowrap: + + $$ + \rho_0 \bm{u} = \sum_i \bm{e}_i f_i + \frac{\delta t}{2} \rho_0 \bm{a}. + $$ + +Combining the previous expressions, + +.. math:: + :nowrap: + + $$ + \bm{u} = \frac{\frac{1}{\rho_0}\sum_i \bm{e}_i f_i + \frac{\delta t}{2}\epsilon \bm{g} + + \frac{\delta t}{2} \frac{\bm{F}_{cp}}{\rho_0}}{1+ \frac{\delta t}{2} \frac{\epsilon \nu}{K}} + $$ + + diff --git a/docs/source/userGuide/models/index.rst b/docs/source/userGuide/models/index.rst index 6a158aa6..dbae0914 100644 --- a/docs/source/userGuide/models/index.rst +++ b/docs/source/userGuide/models/index.rst @@ -18,6 +18,8 @@ Currently supported lattice Boltzmann models greyscale/* + greyscaleColor/* + freeEnergy/* diff --git a/docs/source/userGuide/models/mrt/mrt.rst b/docs/source/userGuide/models/mrt/mrt.rst index 02463480..ac5a4e5b 100644 --- a/docs/source/userGuide/models/mrt/mrt.rst +++ b/docs/source/userGuide/models/mrt/mrt.rst @@ -22,7 +22,7 @@ For additional details please refer to your local system documentation. Model parameters *************************** -The essential model parameters for the color model are +The essential model parameters for the single-phase MRT model are - ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5` diff --git a/doxygen/DoxygenMainpage.h b/doxygen/DoxygenMainpage.h index b88d8ba4..603f6c44 100644 --- a/doxygen/DoxygenMainpage.h +++ b/doxygen/DoxygenMainpage.h @@ -3,6 +3,7 @@ * C/C++ routines * * - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)" + * - \ref Domain.h "Domain structure" * - \ref models "Lattice Boltzmann models" * - \ref ScaLBL_ColorModel "Color model" * - \ref analysis "Analysis routines" diff --git a/models/ColorModel.h b/models/ColorModel.h index 6be45fe1..0c5b2bdd 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -111,7 +111,6 @@ public: std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; std::shared_ptr ScaLBL_Comm_Regular; - //std::shared_ptr Averages; std::shared_ptr Averages; // input database diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index c667eb30..94f64dab 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -356,8 +356,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt AFFINITY=AffinityList[idx]; Sn = SnList[idx]; Sw = SwList[idx]; - Kn = SnList[idx]; - Kw = SwList[idx]; + Kn = KnList[idx]; + Kw = KwList[idx]; idx = NLABELS; } } @@ -523,71 +523,6 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels() delete [] Permeability; } -//void ScaLBL_GreyscaleColorModel::AssignGreyscalePotential() -//{ -// double *psi;//greyscale potential -// psi = new double[N]; -// -// size_t NLABELS=0; -// signed char VALUE=0; -// double AFFINITY=0.f; -// -// auto LabelList = greyscaleColor_db->getVector( "ComponentLabels" ); -// auto AffinityList = greyscaleColor_db->getVector( "ComponentAffinity" ); -// NLABELS=LabelList.size(); -// -// //first, copy over normal phase field -// for (int k=0;kgetVector( "GreySolidLabels" ); -// auto PermeabilityList = greyscaleColor_db->getVector( "PermeabilityList" ); -// NLABELS=GreyLabelList.size(); -// -// for (int k=0;kvoxel_length/Dm->voxel_length); -// idx = NLABELS; -// } -// } -// //update greyscale potential -// psi[n] = psi[n]*Cap_Penalty; -// } -// } -// } -// -// ScaLBL_CopyToDevice(Psi, psi, N*sizeof(double)); -// ScaLBL_Comm->Barrier(); -// delete [] psi; -//} - void ScaLBL_GreyscaleColorModel::Create(){ /* * This function creates the variables needed to run a LBM @@ -635,7 +570,7 @@ void ScaLBL_GreyscaleColorModel::Create(){ //ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);//greyscale potential ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); - //ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &MobilityRatio, sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz); //ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np); @@ -686,8 +621,7 @@ void ScaLBL_GreyscaleColorModel::Create(){ AssignComponentLabels();//do open/black/grey nodes initialization AssignGreySolidLabels(); AssignGreyPoroPermLabels(); - //AssignGreyscalePotential(); - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity); + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity); ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);//porosity doesn't change over time } @@ -787,9 +721,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - int IMAGE_INDEX = 0; - int IMAGE_COUNT = 0; - std::vector ImageList; bool SET_CAPILLARY_NUMBER = false; bool RESCALE_FORCE = false; bool MORPH_ADAPT = false; @@ -845,16 +776,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ /* defaults for simulation protocols */ auto protocol = greyscaleColor_db->getWithDefault( "protocol", "none" ); - if (protocol == "image sequence"){ - // Get the list of images - USE_DIRECT = true; - ImageList = greyscaleColor_db->getVector( "image_sequence"); - IMAGE_INDEX = greyscaleColor_db->getWithDefault( "image_index", 0 ); - IMAGE_COUNT = ImageList.size(); - morph_interval = 10000; - USE_MORPH = true; - } - else if (protocol == "seed water"){ + if (protocol == "seed water"){ morph_delta = -0.05; seed_water = 0.01; USE_SEED = true; @@ -908,15 +830,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ if (rank==0){ printf("********************************************************\n"); - if (protocol == "image sequence"){ - printf(" using protocol = image sequence \n"); - printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); - printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); - printf(" tolerance = %f \n",tolerance); - std::string first_image = ImageList[IMAGE_INDEX]; - printf(" first image in sequence: %s ***\n", first_image.c_str()); - } - else if (protocol == "seed water"){ + if (protocol == "seed water"){ printf(" using protocol = seed water \n"); printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); @@ -963,18 +877,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ } // Halo exchange for phase field ScaLBL_Comm_Regular->SendHalo(Phi); - //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - //Model-1&4 - //ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ////Model-2&3 - //ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); @@ -992,18 +897,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - //Model-1&4 - //ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ////Model-2&3 - //ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); // *************EVEN TIMESTEP************* @@ -1025,18 +921,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - //Model-1&4 - //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ////Model-2&3 - //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); @@ -1054,18 +941,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - //Model-1&4 - //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ////Model-2&3 - //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, - // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, - // alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //************************************************************************ PROFILE_STOP("Update"); @@ -1121,11 +999,13 @@ void ScaLBL_GreyscaleColorModel::Run(){ if (timestep%analysis_interval == 0){ ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure); + ScaLBL_Comm->RegularLayout(Map,MobilityRatio,Averages->MobilityRatio); ScaLBL_Comm->RegularLayout(Map,&Den[0],Averages->Rho_n); ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w); ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x); ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y); ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z); + Averages->Basic(); } @@ -1212,43 +1092,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ double pA = Averages->Oil.p; double pB = Averages->Water.p; double pAB = (pA-pB)/(h*6.0*alpha); - - // -------- The following quantities may not make sense for greyscale simulation -----------// -// double pAc = Averages->gnc.p; -// double pBc = Averages->gwc.p; -// double pAB_connected = (pAc-pBc)/(h*6.0*alpha); -// // connected contribution -// double Vol_nc = Averages->gnc.V/Dm->Volume; -// double Vol_wc = Averages->gwc.V/Dm->Volume; -// double Vol_nd = Averages->gnd.V/Dm->Volume; -// double Vol_wd = Averages->gwd.V/Dm->Volume; -// double Mass_n = Averages->gnc.M + Averages->gnd.M; -// double Mass_w = Averages->gwc.M + Averages->gwd.M; -// double vAc_x = Averages->gnc.Px/Mass_n; -// double vAc_y = Averages->gnc.Py/Mass_n; -// double vAc_z = Averages->gnc.Pz/Mass_n; -// double vBc_x = Averages->gwc.Px/Mass_w; -// double vBc_y = Averages->gwc.Py/Mass_w; -// double vBc_z = Averages->gwc.Pz/Mass_w; -// // disconnected contribution -// double vAd_x = Averages->gnd.Px/Mass_n; -// double vAd_y = Averages->gnd.Py/Mass_n; -// double vAd_z = Averages->gnd.Pz/Mass_n; -// double vBd_x = Averages->gwd.Px/Mass_w; -// double vBd_y = Averages->gwd.Py/Mass_w; -// double vBd_z = Averages->gwd.Pz/Mass_w; -// -// double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z); -// double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z); -// double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z); -// double flow_rate_B_disconnected = (Vol_wd)*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z); -// -// double kAeff_connected = h*h*muA*flow_rate_A_connected/(force_mag); -// double kBeff_connected = h*h*muB*flow_rate_B_connected/(force_mag); -// -// double kAeff_disconnected = h*h*muA*flow_rate_A_disconnected/(force_mag); -// double kBeff_disconnected = h*h*muB*flow_rate_B_disconnected/(force_mag); -// //---------------------------------------------------------------------------------------// double kAeff = h*h*muA*(flow_rate_A)/(force_mag); double kBeff = h*h*muB*(flow_rate_B)/(force_mag); @@ -1300,22 +1143,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ if (MORPH_ADAPT ){ CURRENT_MORPH_TIMESTEPS += analysis_interval; - if (USE_DIRECT){ - // Use image sequence - IMAGE_INDEX++; - MORPH_ADAPT = false; - if (IMAGE_INDEX < IMAGE_COUNT){ - std::string next_image = ImageList[IMAGE_INDEX]; - if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX); - greyscaleColor_db->putScalar("image_index",IMAGE_INDEX); - ImageInit(next_image); - } - else{ - if (rank==0) printf("Finished simulating image sequence \n"); - timestep = timestepMax; - } - } - else if (USE_SEED){ + if (USE_SEED){ delta_volume = volA*Dm->Volume - initial_volume; CURRENT_MORPH_TIMESTEPS += analysis_interval; double massChange = SeedPhaseField(seed_water); @@ -1366,50 +1194,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ // ************************************************************************ } -void ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){ - if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str()); - Mask->Decomp(Filename); - for (int i=0; iid[i]; // save what was read - for (int i=0; iid[i] = Mask->id[i]; // save what was read - - AssignComponentLabels(); - AssignGreySolidLabels(); - AssignGreyPoroPermLabels(); - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity); - ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity); - - //NOTE in greyscale simulations, water may have multiple labels (e.g. 2, 21, 22, etc) - //so the saturaiton calculation is not that straightforward -// double Count = 0.0; -// double PoreCount = 0.0; -// for (int k=1; kComm.sumReduce( Count); -// PoreCount=Dm->Comm.sumReduce( PoreCount); -// if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount); - - ScaLBL_D3Q19_Init(fq, Np); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->Barrier(); - - //ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double)); - - //double saturation = Count/PoreCount; - //return saturation; - -} double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){ srand(time(NULL)); double mass_loss =0.f; @@ -1713,546 +1497,3 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){ */ } -//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-1 -//{ -// // ONLY initialize grey nodes -// // Key input parameters: -// // 1. GreySolidLabels -// // labels for grey nodes -// // 2. GreySolidAffinity -// // affinity ranges [-1,1] -// // oil-wet > 0 -// // water-wet < 0 -// // neutral = 0 -// double *SolidPotential_host = new double [Nx*Ny*Nz]; -// double *GreySolidGrad_host = new double [3*Np]; -// -// size_t NLABELS=0; -// signed char VALUE=0; -// double AFFINITY=0.f; -// -// auto LabelList = greyscaleColor_db->getVector( "GreySolidLabels" ); -// auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); -// -// NLABELS=LabelList.size(); -// if (NLABELS != AffinityList.size()){ -// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); -// } -// -// for (int k=0;kid[n] = 0; // set mask to zero since this is an immobile component -// } -// } -// SolidPotential_host[n] = AFFINITY; -// } -// } -// } -// -// // Calculate grey-solid color-gradient -// double *Dst; -// Dst = new double [3*3*3]; -// for (int kk=0; kk<3; kk++){ -// for (int jj=0; jj<3; jj++){ -// for (int ii=0; ii<3; ii++){ -// int index = kk*9+jj*3+ii; -// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1)); -// } -// } -// } -// double w_face = 1.f; -// double w_edge = 0.5; -// double w_corner = 0.f; -// //local -// Dst[13] = 0.f; -// //faces -// Dst[4] = w_face; -// Dst[10] = w_face; -// Dst[12] = w_face; -// Dst[14] = w_face; -// Dst[16] = w_face; -// Dst[22] = w_face; -// // corners -// Dst[0] = w_corner; -// Dst[2] = w_corner; -// Dst[6] = w_corner; -// Dst[8] = w_corner; -// Dst[18] = w_corner; -// Dst[20] = w_corner; -// Dst[24] = w_corner; -// Dst[26] = w_corner; -// // edges -// Dst[1] = w_edge; -// Dst[3] = w_edge; -// Dst[5] = w_edge; -// Dst[7] = w_edge; -// Dst[9] = w_edge; -// Dst[11] = w_edge; -// Dst[15] = w_edge; -// Dst[17] = w_edge; -// Dst[19] = w_edge; -// Dst[21] = w_edge; -// Dst[23] = w_edge; -// Dst[25] = w_edge; -// -// for (int k=1; kBarrier(); -// delete [] SolidPotential_host; -// delete [] GreySolidGrad_host; -// delete [] Dst; -//} -////----------------------------------------------------------------------------------------------------------// - - -//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-2 & Model-3 -//{ -// // ONLY initialize grey nodes -// // Key input parameters: -// // 1. GreySolidLabels -// // labels for grey nodes -// // 2. GreySolidAffinity -// // affinity ranges [-1,1] -// // oil-wet > 0 -// // water-wet < 0 -// // neutral = 0 -// -// double *GreySolidPhi_host = new double [Nx*Ny*Nz]; -// //initialize grey solid phase field -// for (int k=0;kgetVector( "GreySolidLabels" ); -// auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); -// -// size_t NLABELS=0; -// NLABELS=LabelList.size(); -// if (NLABELS != AffinityList.size()){ -// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); -// } -// -// double *Dst; -// Dst = new double [3*3*3]; -// for (int kk=0; kk<3; kk++){ -// for (int jj=0; jj<3; jj++){ -// for (int ii=0; ii<3; ii++){ -// int index = kk*9+jj*3+ii; -// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1)); -// } -// } -// } -// double w_face = 1.f; -// double w_edge = 1.f; -// double w_corner = 0.f; -// //local -// Dst[13] = 0.f; -// //faces -// Dst[4] = w_face; -// Dst[10] = w_face; -// Dst[12] = w_face; -// Dst[14] = w_face; -// Dst[16] = w_face; -// Dst[22] = w_face; -// // corners -// Dst[0] = w_corner; -// Dst[2] = w_corner; -// Dst[6] = w_corner; -// Dst[8] = w_corner; -// Dst[18] = w_corner; -// Dst[20] = w_corner; -// Dst[24] = w_corner; -// Dst[26] = w_corner; -// // edges -// Dst[1] = w_edge; -// Dst[3] = w_edge; -// Dst[5] = w_edge; -// Dst[7] = w_edge; -// Dst[9] = w_edge; -// Dst[11] = w_edge; -// Dst[15] = w_edge; -// Dst[17] = w_edge; -// Dst[19] = w_edge; -// Dst[21] = w_edge; -// Dst[23] = w_edge; -// Dst[25] = w_edge; -// -// for (int k=1; kid[n]; -// double AFFINITY=0.f; -// // Assign the affinity from the paired list -// for (unsigned int idx=0; idx < NLABELS; idx++){ -// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]); -// if (VALUE == LabelList[idx]){ -// AFFINITY=AffinityList[idx]; -// idx = NLABELS; -// //Mask->id[n] = 0; // set mask to zero since this is an immobile component -// } -// } -// -// if (VALUE>2){//i.e. a grey node -// double neighbor_counter = 0; -// for (int kk=0; kk<3; kk++){ -// for (int jj=0; jj<3; jj++){ -// for (int ii=0; ii<3; ii++){ -// -// int index = kk*9+jj*3+ii; -// double weight= Dst[index]; -// -// int idi=i+ii-1; -// int idj=j+jj-1; -// int idk=k+kk-1; -// -// if (idi < 0) idi=0; -// if (idj < 0) idj=0; -// if (idk < 0) idk=0; -// if (!(idi < Nx)) idi=Nx-1; -// if (!(idj < Ny)) idj=Ny-1; -// if (!(idk < Nz)) idk=Nz-1; -// -// int nn = idk*Nx*Ny + idj*Nx + idi; -// //if (Mask->id[nn] != VALUE){//Model-2:i.e. open nodes, impermeable solid nodes or any other type of greynodes -// if (Mask->id[nn] <=0){//Model-3:i.e. only impermeable solid nodes or any other type of greynodes -// neighbor_counter +=weight; -// } -// } -// } -// } -// if (neighbor_counter>0){ -// GreySolidPhi_host[n] = AFFINITY; -// } -// } -// } -// } -// } -// -// if (rank==0){ -// printf("Number of grey-solid labels: %lu \n",NLABELS); -// for (unsigned int idx=0; idxBarrier(); -// -// //debug -// //FILE *OUTFILE; -// //sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank); -// //OUTFILE = fopen(LocalRankFilename,"wb"); -// //fwrite(GreySolidPhi_host,8,N,OUTFILE); -// //fclose(OUTFILE); -// -// delete [] GreySolidPhi_host; -// delete [] Dst; -//} - -//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4 -//{ -// // ONLY initialize grey nodes -// // Key input parameters: -// // 1. GreySolidLabels -// // labels for grey nodes -// // 2. GreySolidAffinity -// // affinity ranges [-1,1] -// // oil-wet > 0 -// // water-wet < 0 -// // neutral = 0 -// double *SolidPotential_host = new double [Nx*Ny*Nz]; -// double *GreySolidGrad_host = new double [3*Np]; -// -// size_t NLABELS=0; -// signed char VALUE=0; -// double AFFINITY=0.f; -// -// auto LabelList = greyscaleColor_db->getVector( "GreySolidLabels" ); -// auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); -// -// NLABELS=LabelList.size(); -// if (NLABELS != AffinityList.size()){ -// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); -// } -// -// for (int k=0;kid[n] = 0; // set mask to zero since this is an immobile component -// } -// } -// SolidPotential_host[n] = AFFINITY; -// } -// } -// } -// -// // Calculate grey-solid color-gradient -// double *Dst; -// Dst = new double [3*3*3]; -// for (int kk=0; kk<3; kk++){ -// for (int jj=0; jj<3; jj++){ -// for (int ii=0; ii<3; ii++){ -// int index = kk*9+jj*3+ii; -// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1)); -// } -// } -// } -// double w_face = 1.f; -// double w_edge = 0.5; -// double w_corner = 0.f; -// //local -// Dst[13] = 0.f; -// //faces -// Dst[4] = w_face; -// Dst[10] = w_face; -// Dst[12] = w_face; -// Dst[14] = w_face; -// Dst[16] = w_face; -// Dst[22] = w_face; -// // corners -// Dst[0] = w_corner; -// Dst[2] = w_corner; -// Dst[6] = w_corner; -// Dst[8] = w_corner; -// Dst[18] = w_corner; -// Dst[20] = w_corner; -// Dst[24] = w_corner; -// Dst[26] = w_corner; -// // edges -// Dst[1] = w_edge; -// Dst[3] = w_edge; -// Dst[5] = w_edge; -// Dst[7] = w_edge; -// Dst[9] = w_edge; -// Dst[11] = w_edge; -// Dst[15] = w_edge; -// Dst[17] = w_edge; -// Dst[19] = w_edge; -// Dst[21] = w_edge; -// Dst[23] = w_edge; -// Dst[25] = w_edge; -// -// for (int k=1; kSDs(i,j,k)<2.0){ -// GreySolidGrad_host[idx+0*Np] = phi_x; -// GreySolidGrad_host[idx+1*Np] = phi_y; -// GreySolidGrad_host[idx+2*Np] = phi_z; -// } -// else{ -// GreySolidGrad_host[idx+0*Np] = 0.0; -// GreySolidGrad_host[idx+1*Np] = 0.0; -// GreySolidGrad_host[idx+2*Np] = 0.0; -// } -// } -// } -// } -// } -// -// -// if (rank==0){ -// printf("Number of Grey-solid labels: %lu \n",NLABELS); -// for (unsigned int idx=0; idxBarrier(); -// delete [] SolidPotential_host; -// delete [] GreySolidGrad_host; -// delete [] Dst; -//} - - -//--------- This is another old version of calculating greyscale-solid color-gradient modification-------// -// **not working effectively, to be deprecated -//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels() -//{ -// // ONLY initialize grey nodes -// // Key input parameters: -// // 1. GreySolidLabels -// // labels for grey nodes -// // 2. GreySolidAffinity -// // affinity ranges [-1,1] -// // oil-wet > 0 -// // water-wet < 0 -// // neutral = 0 -// -// //double *SolidPotential_host = new double [Nx*Ny*Nz]; -// double *GreySolidPhi_host = new double [Nx*Ny*Nz]; -// signed char VALUE=0; -// double AFFINITY=0.f; -// -// auto LabelList = greyscaleColor_db->getVector( "GreySolidLabels" ); -// auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); -// -// size_t NLABELS=0; -// NLABELS=LabelList.size(); -// if (NLABELS != AffinityList.size()){ -// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); -// } -// -// for (int k=0;kid[n] = 0; // set mask to zero since this is an immobile component -// } -// } -// GreySolidPhi_host[n] = AFFINITY; -// } -// } -// } -// -// if (rank==0){ -// printf("Number of grey-solid labels: %lu \n",NLABELS); -// for (unsigned int idx=0; idxBarrier(); -// -// //debug -// FILE *OUTFILE; -// sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank); -// OUTFILE = fopen(LocalRankFilename,"wb"); -// fwrite(GreySolidPhi_host,8,N,OUTFILE); -// fclose(OUTFILE); -// -// //delete [] SolidPotential_host; -// delete [] GreySolidPhi_host; -//} -//----------------------------------------------------------------------------------------------------------// diff --git a/models/GreyscaleColorModel.h b/models/GreyscaleColorModel.h index 008cb5f3..525eac07 100644 --- a/models/GreyscaleColorModel.h +++ b/models/GreyscaleColorModel.h @@ -15,19 +15,69 @@ Implementation of two-fluid greyscale color lattice boltzmann model #include "ProfilerApp.h" #include "threadpool/thread_pool.h" +/** + * \class ScaLBL_GreyscaleColorModel + * + * @details + * The ScaLBL_GreyscaleColorModel class extends the standard color model incorporate transport + * through sub-resolution "greyscale" regions. + * Momentum transport equations are described by a D3Q19 scheme + * Mass transport equations are described by D3Q7 scheme + */ + + class ScaLBL_GreyscaleColorModel{ public: + /** + * \brief Constructor + * @param RANK processor rank + * @param NP number of processors + * @param COMM MPI communicator + */ ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& COMM); ~ScaLBL_GreyscaleColorModel(); // functions in they should be run + /** + * \brief Read simulation parameters + * @param filename input database file that includes "Color" section + */ void ReadParams(string filename); + + /** + * \brief Read simulation parameters + * @param db0 input database that includes "Color" section + */ void ReadParams(std::shared_ptr db0); + + /** + * \brief Create domain data structures + */ void SetDomain(); + + /** + * \brief Read image data + */ void ReadInput(); + + /** + * \brief Create color model data structures + */ void Create(); + + /** + * \brief Initialize the simulation + */ void Initialize(); + + /** + * \brief Run the simulation + */ void Run(); + + /** + * \brief Debugging function to dump simulation state to disk + */ void WriteDebug(); bool Restart,pBC; @@ -72,7 +122,7 @@ public: double *GreySw; double *GreyKn; double *GreyKw; - //double *ColorGrad; + double *MobilityRatio; double *Velocity; double *Pressure; double *Porosity_dvc; @@ -91,14 +141,24 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); + + /** + * \brief Assign wetting affinity values + */ void AssignComponentLabels(); + + /** + * \brief Assign wetting affinity values in greyscale regions + */ void AssignGreySolidLabels(); + /** + * \brief Assign porosity and permeability in greyscale regions + */ void AssignGreyPoroPermLabels(); - //void AssignGreyscalePotential(); - void ImageInit(std::string filename); - double MorphInit(const double beta, const double morph_delta); + /** + * \brief Seed phase field + */ double SeedPhaseField(const double seed_water_in_oil); - double MorphOpenConnected(double target_volume_change); void WriteVisFiles(); };