OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into electrokinetic

Browse Source
This commit is contained in:
Rex Zhe Li 2021-10-13 16:46:38 +11:00
commit 4e47b0ad90
16 changed files with 534 additions and 2578 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
analysis/GreyPhase.cpp
View File

@ -19,6 +19,7 @@ GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr <Domain> dm):
Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field 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_y.resize(Nx,Ny,Nz); Vel_y.fill(0);
Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0); Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0);
MobilityRatio.resize(Nx,Ny,Nz); MobilityRatio.fill(0);
//......................................... //.........................................
if (Dm->rank()==0){ if (Dm->rank()==0){
@ -89,14 +90,17 @@ void GreyPhaseAnalysis::Basic(){
double nB = Rho_w(n); double nB = Rho_w(n);
double phi = (nA-nB)/(nA+nB); double phi = (nA-nB)/(nA+nB);
double porosity = Porosity(n); double porosity = Porosity(n);
Water_local.M += rho_w*nB*porosity; double mobility_ratio = MobilityRatio(n);
Water_local.Px += porosity*rho_w*nB*Vel_x(n);
Water_local.Py += porosity*rho_w*nB*Vel_y(n); Water_local.M += nB*porosity;
Water_local.Pz += porosity*rho_w*nB*Vel_z(n); Water_local.Px += porosity*(nA+nB)*Vel_x(n)*0.5*(1.0-mobility_ratio);
Oil_local.M += rho_n*nA*porosity; Water_local.Py += porosity*(nA+nB)*Vel_y(n)*0.5*(1.0-mobility_ratio);
Oil_local.Px += porosity*rho_n*nA*Vel_x(n); Water_local.Pz += porosity*(nA+nB)*Vel_z(n)*0.5*(1.0-mobility_ratio);
Oil_local.Py += porosity*rho_n*nA*Vel_y(n);
Oil_local.Pz += porosity*rho_n*nA*Vel_z(n); 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 ){ if ( phi > 0.99 ){
Oil_local.p += Pressure(n); Oil_local.p += Pressure(n);

1
analysis/GreyPhase.h
View File

@ -71,6 +71,7 @@ public:
DoubleArray Vel_x; // velocity field DoubleArray Vel_x; // velocity field
DoubleArray Vel_y; DoubleArray Vel_y;
DoubleArray Vel_z; DoubleArray Vel_z;
DoubleArray MobilityRatio;
GreyPhaseAnalysis(std::shared_ptr <Domain> Dm); GreyPhaseAnalysis(std::shared_ptr <Domain> Dm);
~GreyPhaseAnalysis(); ~GreyPhaseAnalysis();

200
common/Domain.h
View File

@ -16,87 +16,98 @@
#include "common/MPI.h" #include "common/MPI.h"
#include "common/Communication.h" #include "common/Communication.h"
#include "common/Database.h" #include "common/Database.h"
/**
class Domain; * @file Domain.h
template<class TYPE> class PatchData; * \brief Parallel Domain data structures and helper functions
*/
//! Class to hold information about a box /**
* \class Box
*
* @details
* information about a box
*/
class Box { class Box {
public: public:
int ifirst[3]; int ifirst[3];
int ilast[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<class TYPE>
std::shared_ptr<PatchData<TYPE>> createPatchData( DataLocation location ) const;
private:
Box d_box;
int d_owner;
Domain *d_domain;
};
//! Class to hold domain info
class Domain{ class Domain{
public: public:
//! Default constructor /**
* \brief Constructor
* @param db input database
* @param Communicator MPI communicator
*/
Domain( std::shared_ptr<Database> db, const Utilities::MPI& Communicator); Domain( std::shared_ptr<Database> 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, Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
double lx, double ly, double lz, int BC); double lx, double ly, double lz, int BC);
//! Empty constructor /**
* \brief Empty constructor
*/
Domain() = delete; Domain() = delete;
//! Copy constructor /**
* \brief Copy constructor
*/
Domain( const Domain& ) = delete; Domain( const Domain& ) = delete;
//! Assignment operator /**
* \brief Assignment operator
*/
Domain& operator=( const Domain& ) = delete; Domain& operator=( const Domain& ) = delete;
//! Destructor /**
* \brief Destructor
*/
~Domain(); ~Domain();
//! Get the database /**
* \brief Get the database
*/
inline std::shared_ptr<const Database> getDatabase() const { return d_db; } inline std::shared_ptr<const Database> getDatabase() const { return d_db; }
//! Get the domain box /**
* \brief Get the domain box
*/
inline const Box& getBox() const { return d_box; } inline const Box& getBox() const { return d_box; }
//! Get local patch /**
* \brief Get local patch
*/
inline const Patch& getLocalPatch() const { return *d_localPatch; } inline const Patch& getLocalPatch() const { return *d_localPatch; }
//! Get all patches /**
* \brief Get all patches
*/
inline const std::vector<Patch>& getAllPatch() const { return d_patches; } inline const std::vector<Patch>& getAllPatch() const { return d_patches; }
private: private:
/**
* \brief initialize from database
*/
void initialize( std::shared_ptr<Database> db ); void initialize( std::shared_ptr<Database> db );
std::shared_ptr<Database> d_db; std::shared_ptr<Database> d_db;
@ -124,6 +135,9 @@ public: // Public variables (need to create accessors instead)
//********************************** //**********************************
// MPI ranks for all 18 neighbors // MPI ranks for all 18 neighbors
//********************************** //**********************************
/**
* \brief Compute the porosity based on the current domain id file
*/
inline double Porosity() const { return porosity; } inline double Porosity() const { return porosity; }
inline int iproc() const { return rank_info.ix; } inline int iproc() const { return rank_info.ix; }
inline int jproc() const { return rank_info.jy; } inline int jproc() const { return rank_info.jy; }
@ -165,22 +179,78 @@ public: // Public variables (need to create accessors instead)
// Solid indicator function // Solid indicator function
std::vector<signed char> id; std::vector<signed char> id;
/**
* \brief Read domain IDs from file
*/
void ReadIDs(); void ReadIDs();
/**
* \brief Compute the porosity
*/
void ComputePorosity(); void ComputePorosity();
/**
* \brief Read image and perform domain decomposition
* @param filename - name of file to read IDs
*/
void Decomp( const std::string& filename ); 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); 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(); void CommInit();
/**
* \brief Count number of pore nodes (labels > 1)
*/
int PoreCount(); 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); 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 ); 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 ); void AggregateLabels( const std::string& filename, DoubleArray &UserData );
private: 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); 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 UnpackID(int *list, int count, signed char *recvbuf, signed char *ID);
void CommHaloIDs();
//...................................................................................... //......................................................................................
MPI_Request req1[18], req2[18]; MPI_Request req1[18], req2[18];
@ -198,6 +268,44 @@ private:
}; };
template<class TYPE> 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<class TYPE>
std::shared_ptr<PatchData<TYPE>> createPatchData( DataLocation location ) const;
private:
Box d_box;
int d_owner;
Domain *d_domain;
};
// Class to hold data on a patch // Class to hold data on a patch
template<class TYPE> template<class TYPE>

4
common/ScaLBL.h
View File

@ -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); 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, 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 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); 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, 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 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); double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);

92
cpu/GreyscaleColor.cpp
View File

@ -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); -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
dist[16*Np+n] = fq; dist[16*Np+n] = fq;
// q = 17 // q = 17
fq = mrt_V1*rho+mrt_V9*m1 fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) +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 //CP: capillary penalty
// also turn off recoloring for grey nodes // 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, 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 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){ 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 */ /* Corey model parameters */
double Kn_grey,Kw_grey; double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB; 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_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802; 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]; nB = Den[Np + n];
porosity = Poros[n]; porosity = Poros[n];
GreyDiff = Perm[n]; //GreyDiff = Perm[n];
perm = 1.0; perm = 1.0;
W = GreySolidW[n]; W = GreySolidW[n];
Sn_grey = GreySn[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; Krw_grey = 0.0;
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){ if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey; perm = Kw_grey;
Krw_grey = Kw_grey;
Swn = 0.0; Swn = 0.0;
} }
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); 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 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 // 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)); //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){ else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey; perm = Kn_grey;
Krn_grey = Kn_grey;
Swn = 1.0; 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 // Get the 1D index based on regular data layout
ijk = Map[n]; 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---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 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---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 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---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 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 //CP: capillary penalty
// also turn off recoloring for grey nodes // 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, 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 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){ 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 */ /* Corey model parameters */
double Kn_grey,Kw_grey; double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB; 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 GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity; double porosity;
@ -2235,7 +2248,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
nB = Den[Np + n]; nB = Den[Np + n];
porosity = Poros[n]; porosity = Poros[n];
GreyDiff = Perm[n]; //GreyDiff = Perm[n];
perm = 1.0; perm = 1.0;
W = GreySolidW[n]; W = GreySolidW[n];
Sn_grey = GreySn[n]; Sn_grey = GreySn[n];
@ -2255,25 +2268,38 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
mobility_ratio = 1.0;
Krn_grey = 0.0; Krn_grey = 0.0;
Krw_grey = 0.0; Krw_grey = 0.0;
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){ if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey; perm = Kw_grey;
Krw_grey = Kw_grey;
Swn = 0.0; Swn = 0.0;
} }
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey); 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 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 // 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)); //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){ else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey; perm = Kn_grey;
Krn_grey = Kn_grey;
Swn = 1.0; 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 // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
@ -2861,7 +2887,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//----------------newly added for better control of recoloring---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 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---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 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---------------// //----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){ 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; delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio); jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 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<Np; n++){
// porosity = Porosity[n];
// if (porosity==0.0) porosity=1.f;
// dist[n] = 0.3333333333333333/porosity;
// dist[Np+n] = 0.055555555555555555/porosity; //double(100*n)+1.f;
// dist[2*Np+n] = 0.055555555555555555/porosity; //double(100*n)+2.f;
// dist[3*Np+n] = 0.055555555555555555/porosity; //double(100*n)+3.f;
// dist[4*Np+n] = 0.055555555555555555/porosity; //double(100*n)+4.f;
// dist[5*Np+n] = 0.055555555555555555/porosity; //double(100*n)+5.f;
// dist[6*Np+n] = 0.055555555555555555/porosity; //double(100*n)+6.f;
// dist[7*Np+n] = 0.0277777777777778/porosity; //double(100*n)+7.f;
// dist[8*Np+n] = 0.0277777777777778/porosity; //double(100*n)+8.f;
// dist[9*Np+n] = 0.0277777777777778/porosity; //double(100*n)+9.f;
// dist[10*Np+n] = 0.0277777777777778/porosity; //double(100*n)+10.f;
// dist[11*Np+n] = 0.0277777777777778/porosity; //double(100*n)+11.f;
// dist[12*Np+n] = 0.0277777777777778/porosity; //double(100*n)+12.f;
// dist[13*Np+n] = 0.0277777777777778/porosity; //double(100*n)+13.f;
// dist[14*Np+n] = 0.0277777777777778/porosity; //double(100*n)+14.f;
// dist[15*Np+n] = 0.0277777777777778/porosity; //double(100*n)+15.f;
// dist[16*Np+n] = 0.0277777777777778/porosity; //double(100*n)+16.f;
// dist[17*Np+n] = 0.0277777777777778/porosity; //double(100*n)+17.f;
// dist[18*Np+n] = 0.0277777777777778/porosity; //double(100*n)+18.f;
// }
//}

1755
cuda/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

7
docs/source/developerGuide/doxygen/Domain.rst Normal file
View File

@ -0,0 +1,7 @@
============================================
Domain
============================================
.. doxygenfile:: Domain.h
:project: LBPM Doxygen

2
docs/source/developerGuide/index.rst
View File

@ -17,6 +17,8 @@ into the framework.
doxygen/ScaLBL.rst doxygen/ScaLBL.rst
doxygen/Domain.rst
doxygen/Models.rst doxygen/Models.rst
doxygen/Analysis.rst doxygen/Analysis.rst

87
docs/source/userGuide/models/greyscale/greyscale.rst
View File

@ -1,7 +1,90 @@
============================================= =============================================
Greyscale model model Greyscale model
============================================= =============================================
The LBPM greyscale lattice Boltzmann model is constructed to approximate the 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 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}}
$$

65
docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst Normal file
View File

@ -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}}
$$

2
docs/source/userGuide/models/index.rst
View File

@ -18,6 +18,8 @@ Currently supported lattice Boltzmann models
greyscale/* greyscale/*
greyscaleColor/*
freeEnergy/* freeEnergy/*

2
docs/source/userGuide/models/mrt/mrt.rst
View File

@ -22,7 +22,7 @@ For additional details please refer to your local system documentation.
Model parameters 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` - ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5`

1
doxygen/DoxygenMainpage.h
View File

@ -3,6 +3,7 @@
* C/C++ routines * C/C++ routines
* *
* - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)" * - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)"
* - \ref Domain.h "Domain structure"
* - \ref models "Lattice Boltzmann models" * - \ref models "Lattice Boltzmann models"
* - \ref ScaLBL_ColorModel "Color model" * - \ref ScaLBL_ColorModel "Color model"
* - \ref analysis "Analysis routines" * - \ref analysis "Analysis routines"

1
models/ColorModel.h
View File

@ -111,7 +111,6 @@ public:
std::shared_ptr<Domain> Mask; // this domain is for lbm std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm; std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular; std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
//std::shared_ptr<TwoPhase> Averages;
std::shared_ptr<SubPhase> Averages; std::shared_ptr<SubPhase> Averages;
// input database // input database

783
models/GreyscaleColorModel.cpp
View File

@ -356,8 +356,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
AFFINITY=AffinityList[idx]; AFFINITY=AffinityList[idx];
Sn = SnList[idx]; Sn = SnList[idx];
Sw = SwList[idx]; Sw = SwList[idx];
Kn = SnList[idx]; Kn = KnList[idx];
Kw = SwList[idx]; Kw = KwList[idx];
idx = NLABELS; idx = NLABELS;
} }
} }
@ -523,71 +523,6 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability; 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<int>( "ComponentLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "ComponentAffinity" );
// NLABELS=LabelList.size();
//
// //first, copy over normal phase field
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// // 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;
// }
// }
// // fluid labels are reserved
// if (VALUE == 1) AFFINITY=1.0;
// else if (VALUE == 2) AFFINITY=-1.0;
// psi[n] = AFFINITY;
// }
// }
// }
//
// //second, scale the phase field for grey nodes
// double Cap_Penalty=1.f;
// auto GreyLabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto PermeabilityList = greyscaleColor_db->getVector<double>( "PermeabilityList" );
// NLABELS=GreyLabelList.size();
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// Cap_Penalty=1.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// if (VALUE == GreyLabelList[idx]){
// Cap_Penalty=alpha*W/sqrt(PermeabilityList[idx]/Dm->voxel_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(){ void ScaLBL_GreyscaleColorModel::Create(){
/* /*
* This function creates the variables needed to run a LBM * 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 **) &Psi, sizeof(double)*Nx*Ny*Nz);//greyscale potential
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*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 **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np);
@ -686,7 +621,6 @@ void ScaLBL_GreyscaleColorModel::Create(){
AssignComponentLabels();//do open/black/grey nodes initialization AssignComponentLabels();//do open/black/grey nodes initialization
AssignGreySolidLabels(); AssignGreySolidLabels();
AssignGreyPoroPermLabels(); 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 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; int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
bool SET_CAPILLARY_NUMBER = false; bool SET_CAPILLARY_NUMBER = false;
bool RESCALE_FORCE = false; bool RESCALE_FORCE = false;
bool MORPH_ADAPT = false; bool MORPH_ADAPT = false;
@ -845,16 +776,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
/* defaults for simulation protocols */ /* defaults for simulation protocols */
auto protocol = greyscaleColor_db->getWithDefault<std::string>( "protocol", "none" ); auto protocol = greyscaleColor_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){ if (protocol == "seed water"){
// Get the list of images
USE_DIRECT = true;
ImageList = greyscaleColor_db->getVector<std::string>( "image_sequence");
IMAGE_INDEX = greyscaleColor_db->getWithDefault<int>( "image_index", 0 );
IMAGE_COUNT = ImageList.size();
morph_interval = 10000;
USE_MORPH = true;
}
else if (protocol == "seed water"){
morph_delta = -0.05; morph_delta = -0.05;
seed_water = 0.01; seed_water = 0.01;
USE_SEED = true; USE_SEED = true;
@ -908,15 +830,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (rank==0){ if (rank==0){
printf("********************************************************\n"); printf("********************************************************\n");
if (protocol == "image sequence"){ if (protocol == "seed water"){
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"){
printf(" using protocol = seed water \n"); printf(" using protocol = seed water \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
@ -963,18 +877,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
} }
// Halo exchange for phase field // Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi); 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,MobilityRatio,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); 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_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
@ -992,18 +897,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); 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,MobilityRatio,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); 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(); ScaLBL_Comm->Barrier();
// *************EVEN TIMESTEP************* // *************EVEN TIMESTEP*************
@ -1025,18 +921,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
ScaLBL_Comm_Regular->SendHalo(Phi); 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,MobilityRatio,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); 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_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
@ -1054,18 +941,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); 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,MobilityRatio,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); 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(); ScaLBL_Comm->Barrier();
//************************************************************************ //************************************************************************
PROFILE_STOP("Update"); PROFILE_STOP("Update");
@ -1121,11 +999,13 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (timestep%analysis_interval == 0){ if (timestep%analysis_interval == 0){
ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure); 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[0],Averages->Rho_n);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w); ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x); ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y); ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z); ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z);
Averages->Basic(); Averages->Basic();
} }
@ -1213,43 +1093,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
double pB = Averages->Water.p; double pB = Averages->Water.p;
double pAB = (pA-pB)/(h*6.0*alpha); 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 kAeff = h*h*muA*(flow_rate_A)/(force_mag);
double kBeff = h*h*muB*(flow_rate_B)/(force_mag); double kBeff = h*h*muB*(flow_rate_B)/(force_mag);
@ -1300,22 +1143,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (MORPH_ADAPT ){ if (MORPH_ADAPT ){
CURRENT_MORPH_TIMESTEPS += analysis_interval; CURRENT_MORPH_TIMESTEPS += analysis_interval;
if (USE_DIRECT){ if (USE_SEED){
// 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<int>("image_index",IMAGE_INDEX);
ImageInit(next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
timestep = timestepMax;
}
}
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume; delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval; CURRENT_MORPH_TIMESTEPS += analysis_interval;
double massChange = SeedPhaseField(seed_water); 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; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[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; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// if (id[Nx*Ny*k+Nx*j+i] == 2){
// PoreCount++;
// Count++;
// }
// else if (id[Nx*Ny*k+Nx*j+i] == 1){
// PoreCount++;
// }
// }
// }
// }
// Count=Dm->Comm.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){ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL)); srand(time(NULL));
double mass_loss =0.f; 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<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // 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
// }
// }
// 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; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// int idx=Map(i,j,k);
// if (!(idx < 0)){
// double phi_x = 0.f;
// double phi_y = 0.f;
// double phi_z = 0.f;
// 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;
// double vec_x = double(ii-1);
// double vec_y = double(jj-1);
// double vec_z = double(kk-1);
// double GWNS=SolidPotential_host[nn];
// phi_x += GWNS*weight*vec_x;
// phi_y += GWNS*weight*vec_y;
// phi_z += GWNS*weight*vec_z;
// }
// }
// }
// GreySolidGrad_host[idx+0*Np] = phi_x;
// GreySolidGrad_host[idx+1*Np] = phi_y;
// GreySolidGrad_host[idx+2*Np] = phi_z;
// }
// }
// }
// }
//
// if (rank==0){
// printf("Number of Grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_Comm->Barrier();
// 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;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// GreySolidPhi_host[n]=0.f;
// }
// }
// }
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "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; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
//
// int n = k*Nx*Ny+j*Nx+i;
// signed char VALUE=Mask->id[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; idx<NLABELS; idx++){
// signed char VALUE=LabelList[idx];
// double AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_Comm->Barrier();
//
// //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<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // 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
// }
// }
// 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; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// int idx=Map(i,j,k);
// if (!(idx < 0)){
// double phi_x = 0.f;
// double phi_y = 0.f;
// double phi_z = 0.f;
// 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;
// double vec_x = double(ii-1);
// double vec_y = double(jj-1);
// double vec_z = double(kk-1);
// double GWNS=SolidPotential_host[nn];
// phi_x += GWNS*weight*vec_x;
// phi_y += GWNS*weight*vec_y;
// phi_z += GWNS*weight*vec_z;
// }
// }
// }
// if (Averages->SDs(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; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_Comm->Barrier();
// 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<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "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;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // 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
// }
// }
// GreySolidPhi_host[n] = AFFINITY;
// }
// }
// }
//
// if (rank==0){
// printf("Number of grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, solid-affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_Comm->Barrier();
//
// //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;
//}
//----------------------------------------------------------------------------------------------------------//

70
models/GreyscaleColorModel.h
View File

@ -15,19 +15,69 @@ Implementation of two-fluid greyscale color lattice boltzmann model
#include "ProfilerApp.h" #include "ProfilerApp.h"
#include "threadpool/thread_pool.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{ class ScaLBL_GreyscaleColorModel{
public: 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(int RANK, int NP, const Utilities::MPI& COMM);
~ScaLBL_GreyscaleColorModel(); ~ScaLBL_GreyscaleColorModel();
// functions in they should be run // functions in they should be run
/**
* \brief Read simulation parameters
* @param filename input database file that includes "Color" section
*/
void ReadParams(string filename); void ReadParams(string filename);
/**
* \brief Read simulation parameters
* @param db0 input database that includes "Color" section
*/
void ReadParams(std::shared_ptr<Database> db0); void ReadParams(std::shared_ptr<Database> db0);
/**
* \brief Create domain data structures
*/
void SetDomain(); void SetDomain();
/**
* \brief Read image data
*/
void ReadInput(); void ReadInput();
/**
* \brief Create color model data structures
*/
void Create(); void Create();
/**
* \brief Initialize the simulation
*/
void Initialize(); void Initialize();
/**
* \brief Run the simulation
*/
void Run(); void Run();
/**
* \brief Debugging function to dump simulation state to disk
*/
void WriteDebug(); void WriteDebug();
bool Restart,pBC; bool Restart,pBC;
@ -72,7 +122,7 @@ public:
double *GreySw; double *GreySw;
double *GreyKn; double *GreyKn;
double *GreyKw; double *GreyKw;
//double *ColorGrad; double *MobilityRatio;
double *Velocity; double *Velocity;
double *Pressure; double *Pressure;
double *Porosity_dvc; double *Porosity_dvc;
@ -91,14 +141,24 @@ private:
//int rank,nprocs; //int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0); void LoadParams(std::shared_ptr<Database> db0);
/**
* \brief Assign wetting affinity values
*/
void AssignComponentLabels(); void AssignComponentLabels();
/**
* \brief Assign wetting affinity values in greyscale regions
*/
void AssignGreySolidLabels(); void AssignGreySolidLabels();
/**
* \brief Assign porosity and permeability in greyscale regions
*/
void AssignGreyPoroPermLabels(); void AssignGreyPoroPermLabels();
//void AssignGreyscalePotential(); /**
void ImageInit(std::string filename); * \brief Seed phase field
double MorphInit(const double beta, const double morph_delta); */
double SeedPhaseField(const double seed_water_in_oil); double SeedPhaseField(const double seed_water_in_oil);
double MorphOpenConnected(double target_volume_change);
void WriteVisFiles(); void WriteVisFiles();
}; };