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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
JamesEMcclure 2020-10-13 14:19:35 -04:00
commit c479402256
13 changed files with 721 additions and 629 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

259
analysis/GreyPhase.cpp Normal file
View File

@ -0,0 +1,259 @@
#include "analysis/GreyPhase.h"
// Constructor
GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr <Domain> dm):
Dm(dm)
{
Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0;
// Global arrays
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
Porosity.resize(Nx,Ny,Nz); Porosity.fill(0);
//PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0);
Rho_n.resize(Nx,Ny,Nz); Rho_n.fill(0);
Rho_w.resize(Nx,Ny,Nz); Rho_w.fill(0);
Pressure.resize(Nx,Ny,Nz); Pressure.fill(0);
//Phi.resize(Nx,Ny,Nz); Phi.fill(0);
//DelPhi.resize(Nx,Ny,Nz); DelPhi.fill(0);
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);
//.........................................
if (Dm->rank()==0){
bool WriteHeader=false;
TIMELOG = fopen("timelog.csv","r");
if (TIMELOG != NULL)
fclose(TIMELOG);
else
WriteHeader=true;
TIMELOG = fopen("timelog.csv","a+");
if (WriteHeader)
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"sw krw krn vw vn pw pn\n");
}
}
}
// Destructor
GreyPhaseAnalysis::~GreyPhaseAnalysis()
{
}
void GreyPhaseAnalysis::Write(int timestep)
{
}
void GreyPhaseAnalysis::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double B, double GreyPorosity)
{
Fx = force_x;
Fy = force_y;
Fz = force_z;
rho_n = rhoA;
rho_w = rhoB;
nu_n = (tauA-0.5)/3.f;
nu_w = (tauB-0.5)/3.f;
gamma_wn = 6.0*alpha;
beta = B;
grey_porosity = GreyPorosity;
}
void GreyPhaseAnalysis::Basic(){
int i,j,k,n,imin,jmin,kmin,kmax;
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
imin=jmin=1;
if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z;
if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z;
Water_local.reset();
Oil_local.reset();
double count_w = 0.0;
double count_n = 0.0;
for (k=kmin; k<kmax; k++){
for (j=jmin; j<Ny-1; j++){
for (i=imin; i<Nx-1; i++){
n = k*Nx*Ny + j*Nx + i;
// Compute volume averages
if ( Dm->id[n] > 0 ){
// compute density
double nA = Rho_n(n);
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);
if ( phi > 0.99 ){
Oil_local.p += Pressure(n);
//Oil_local.p += pressure*(rho_n*nA)/(rho_n*nA+rho_w*nB);
count_n += 1.0;
}
else if ( phi < -0.99 ){
Water_local.p += Pressure(n);
//Water_local.p += pressure*(rho_w*nB)/(rho_n*nA+rho_w*nB);
count_w += 1.0;
}
}
}
}
}
Oil.M=sumReduce( Dm->Comm, Oil_local.M);
Oil.Px=sumReduce( Dm->Comm, Oil_local.Px);
Oil.Py=sumReduce( Dm->Comm, Oil_local.Py);
Oil.Pz=sumReduce( Dm->Comm, Oil_local.Pz);
Water.M=sumReduce( Dm->Comm, Water_local.M);
Water.Px=sumReduce( Dm->Comm, Water_local.Px);
Water.Py=sumReduce( Dm->Comm, Water_local.Py);
Water.Pz=sumReduce( Dm->Comm, Water_local.Pz);
//Oil.p /= Oil.M;
//Water.p /= Water.M;
count_w=sumReduce( Dm->Comm, count_w);
count_n=sumReduce( Dm->Comm, count_n);
if (count_w > 0.0)
Water.p=sumReduce( Dm->Comm, Water_local.p) / count_w;
else
Water.p = 0.0;
if (count_n > 0.0)
Oil.p=sumReduce( Dm->Comm, Oil_local.p) / count_n;
else
Oil.p = 0.0;
// check for NaN
bool err=false;
if (Water.M != Water.M) err=true;
if (Water.p != Water.p) err=true;
if (Water.Px != Water.Px) err=true;
if (Water.Py != Water.Py) err=true;
if (Water.Pz != Water.Pz) err=true;
if (Oil.M != Oil.M) err=true;
if (Oil.p != Oil.p) err=true;
if (Oil.Px != Oil.Px) err=true;
if (Oil.Py != Oil.Py) err=true;
if (Oil.Pz != Oil.Pz) err=true;
if (Dm->rank() == 0){
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = 0.0;
double dir_y = 0.0;
double dir_z = 0.0;
if (force_mag > 0.0){
dir_x = Fx/force_mag;
dir_y = Fy/force_mag;
dir_z = Fz/force_mag;
}
else {
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
}
if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){
// compute the pressure drop
double pressure_drop = (Pressure(Nx*Ny + Nx + 1) - 1.0) / 3.0;
double length = ((Nz-2)*Dm->nprocz());
force_mag -= pressure_drop/length;
}
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
saturation=Water.M/(Water.M + Oil.M); // assume constant density
water_flow_rate=grey_porosity*saturation*(Water.Px*dir_x + Water.Py*dir_y + Water.Pz*dir_z)/Water.M;
oil_flow_rate =grey_porosity*(1.0-saturation)*(Oil.Px*dir_x + Oil.Py*dir_y + Oil.Pz*dir_z)/Oil.M;
double h = Dm->voxel_length;
//TODO check if need greyporosity or domain porosity ? - compare to analytical solution
double krn = h*h*nu_n*oil_flow_rate / force_mag ;
double krw = h*h*nu_w*water_flow_rate / force_mag;
//printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*water_flow_rate,h*oil_flow_rate, Water.p, Oil.p);
fflush(TIMELOG);
}
if (err==true){
// exception if simulation produceds NaN
printf("GreyPhaseAnalysis.cpp: NaN encountered, may need to check simulation parameters \n");
}
ASSERT(err==false);
}
/*
inline void InterfaceTransportMeasures( double beta, double rA, double rB, double nA, double nB,
double nx, double ny, double nz, double ux, double uy, double uz, interface &I){
double A1,A2,A3,A4,A5,A6;
double B1,B2,B3,B4,B5,B6;
double nAB,delta;
// Instantiate mass transport distributions
// Stationary value - distribution 0
nAB = 1.0/(nA+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;
//...............................................
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
A3 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
B3 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
A4 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
B4 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
//...............................................
// q = 4
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
A5 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
B5 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
A6 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
B6 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
double unx = (A1-A2);
double uny = (A3-A4);
double unz = (A5-A6);
double uwx = (B1-B2);
double uwy = (B3-B4);
double uwz = (B5-B6);
I.Mn += rA*nA;
I.Mw += rB*nB;
I.Pnx += rA*nA*unx;
I.Pny += rA*nA*uny;
I.Pnz += rA*nA*unz;
I.Pwx += rB*nB*uwx;
I.Pwy += rB*nB*uwy;
I.Pwz += rB*nB*uwz;
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
}
*/

71
analysis/GreyPhase.h Normal file
View File

@ -0,0 +1,71 @@
/*
* Sub-phase averaging tools
*/
#ifndef GreyPhase_INC
#define GreyPhase_INC
#include <vector>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
class GreyPhase{
public:
double p;
double M,Px,Py,Pz;
void reset(){
p=M=Px=Py=Pz=0.0;
}
private:
};
class GreyPhaseAnalysis{
public:
std::shared_ptr <Domain> Dm;
double Volume;
// input variables
double rho_n, rho_w;
double nu_n, nu_w;
double gamma_wn, beta;
double Fx, Fy, Fz;
double grey_porosity;
// outputs
double saturation,water_flow_rate, oil_flow_rate;
//simulation outputs (averaged values)
GreyPhase Water, Oil;
GreyPhase Water_local, Oil_local;
//...........................................................................
int Nx,Ny,Nz;
//IntArray PhaseID; // Phase ID array
DoubleArray SDs; // contains porosity map
DoubleArray Porosity; // contains porosity map
DoubleArray Rho_n; // density field
DoubleArray Rho_w; // density field
//DoubleArray Phi; // phase indicator field
//DoubleArray DelPhi; // Magnitude of Gradient of the phase indicator field
DoubleArray Pressure; // pressure field
DoubleArray Vel_x; // velocity field
DoubleArray Vel_y;
DoubleArray Vel_z;
GreyPhaseAnalysis(std::shared_ptr <Domain> Dm);
~GreyPhaseAnalysis();
void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double beta, double GreyPorosity);
void Basic();
void Write(int time);
private:
FILE *TIMELOG;
};
#endif

152
common/Domain.cpp
View File

@ -181,6 +181,7 @@ Domain::~Domain()
delete [] recvData_yZ; delete [] recvData_Yz; delete [] recvData_YZ;
// Free id
delete [] id;
// Free the communicator
if ( Comm != MPI_COMM_WORLD && Comm != MPI_COMM_NULL ) {
MPI_Comm_free(&Comm);
@ -691,6 +692,7 @@ void Domain::Decomp( const std::string& Filename )
//.........................................................
}
void Domain::AggregateLabels( const std::string& filename ){
int nx = Nx;
@ -1247,7 +1249,7 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData);
}
// Ideally stuff below here should be moved somewhere else -- doesn't really belong here
// TODO Ideally stuff below here should be moved somewhere else -- doesn't really belong here
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np)
{
double value;
@ -1302,3 +1304,151 @@ void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
File.close();
}
void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData)
{
//........................................................................................
// Reading the user-defined input file
// NOTE: so far it only supports BC=0 (periodic) and BC=5 (mixed reflection)
// because if checkerboard or inlet/outlet buffer layers are added, the
// value of the void space is undefined.
// NOTE: if BC=5 is used, where the inlet and outlet layers of the domain are modified,
// user needs to modify the input file accordingly before LBPM simulator read
// the input file.
//........................................................................................
int rank_offset = 0;
int RANK = rank();
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k,n;
//TODO These offset we may still need them
int64_t xStart,yStart,zStart;
xStart=yStart=zStart=0;
// Read domain parameters
// TODO currently the size of the data is still read from Domain{};
// but user may have a user-specified size
auto size = database->getVector<int>( "n" );
auto SIZE = database->getVector<int>( "N" );
auto nproc = database->getVector<int>( "nproc" );
//TODO currently the funcationality "offset" is disabled as the user-defined input data may have a different size from that of the input domain
if (database->keyExists( "offset" )){
auto offset = database->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
global_Nx = SIZE[0];
global_Ny = SIZE[1];
global_Nz = SIZE[2];
nprocs=nprocx*nprocy*nprocz;
double *SegData = NULL;
if (RANK==0){
printf("User-defined input file: %s (data type: %s)\n",Filename.c_str(),Datatype.c_str());
printf("NOTE: currently only BC=0 or 5 supports user-defined input file!\n");
// Rank=0 reads the entire segmented data and distributes to worker processes
printf("Dimensions of the user-defined input file: %ld x %ld x %ld \n",global_Nx,global_Ny,global_Nz);
int64_t SIZE = global_Nx*global_Ny*global_Nz;
if (Datatype == "double"){
printf("Reading input data as double precision floating number\n");
SegData = new double[SIZE];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading user-defined file!\n");
size_t ReadSeg;
ReadSeg=fread(SegData,8,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading file: %s\n",Filename.c_str());
fclose(SEGDAT);
}
else{
ERROR("Error: User-defined input file only supports double-precision floating number!\n");
}
printf("Read file successfully from %s \n",Filename.c_str());
}
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);
// number of sites to use for periodic boundary condition transition zone
//int64_t z_transition_size = (nprocz*nz - (global_Nz - zStart))/2;
//if (z_transition_size < 0) z_transition_size=0;
int64_t z_transition_size = 0;
//char LocalRankFilename[1000];//just for debug
double *loc_id;
loc_id = new double [(nx+2)*(ny+2)*(nz+2)];
// Set up the sub-domains
if (RANK==0){
printf("Decomposing user-defined input file\n");
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
printf("Size of transition region: %ld \n", z_transition_size);
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*nprocx*nprocy + jp*nprocx + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int64_t x = xStart + ip*nx + i-1;
int64_t y = yStart + jp*ny + j-1;
// int64_t z = zStart + kp*nz + k-1;
int64_t z = zStart + kp*nz + k-1 - z_transition_size;
if (x<xStart) x=xStart;
if (!(x<global_Nx)) x=global_Nx-1;
if (y<yStart) y=yStart;
if (!(y<global_Ny)) y=global_Ny-1;
if (z<zStart) z=zStart;
if (!(z<global_Nz)) z=global_Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int64_t nglobal = z*global_Nx*global_Ny+y*global_Nx+x;
loc_id[nlocal] = SegData[nglobal];
}
}
}
if (rnk==0){
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
UserData[nlocal] = loc_id[nlocal];
}
}
}
}
else{
//printf("Sending data to process %i \n", rnk);
MPI_Send(loc_id,N,MPI_DOUBLE,rnk,15,Comm);
}
// Write the data for this rank data
// NOTE just for debug
//sprintf(LocalRankFilename,"%s.%05i",Filename.c_str(),rnk+rank_offset);
//FILE *ID = fopen(LocalRankFilename,"wb");
//fwrite(loc_id,8,(nx+2)*(ny+2)*(nz+2),ID);
//fclose(ID);
}
}
}
}
else{
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(UserData,N,MPI_DOUBLE,0,15,Comm,MPI_STATUS_IGNORE);
}
//Comm.barrier();
MPI_Barrier(Comm);
}

4
common/Domain.h
View File

@ -194,6 +194,7 @@ public: // Public variables (need to create accessors instead)
void ReadIDs();
void Decomp( const std::string& filename );
void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();
@ -260,6 +261,9 @@ private:
};
//void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
//void ReadFromFile(const std::string& Filename, DoubleArray &Mesh);
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np);
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, size_t Np);

4
common/ScaLBL.h
View File

@ -193,12 +193,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
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 *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,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, int strideY, int strideZ, int start, int finish, int Np);
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 *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,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, int strideY, int strideZ, int start, int finish, int Np);

4
cpu/Greyscale.cpp
View File

@ -1974,7 +1974,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
}
//Calculate pressure for MRT model
pressure=rho/3.f;
pressure=rho/3.f/porosity;
//-------------------- MRT collison where body force has NO higher-order terms -------------//
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);
@ -2472,7 +2472,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int f
}
//Calculate pressure for Incompressible-MRT model
pressure=rho/3.f;
pressure=rho/3.f/porosity;
//-------------------- IMRT collison where body force has NO higher-order terms -------------//
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);

6
cpu/GreyscaleColor.cpp
View File

@ -1,7 +1,7 @@
#include <math.h>
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,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, int strideY, int strideZ, int start, int finish, int Np){
@ -494,6 +494,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Map, d
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
Pressure[n] = rho/3.f/porosity;
//........................................................................
//..............carry out relaxation process..............................
@ -709,7 +710,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Map, d
}
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 *Velocity,
double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,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, int strideY, int strideZ, int start, int finish, int Np){
@ -1148,6 +1149,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
Pressure[n] = rho/3.f/porosity;
//........................................................................
//..............carry out relaxation process..............................

4
gpu/Greyscale.cu
View File

@ -2006,7 +2006,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *
}
//Calculate pressure for MRT model
pressure=rho/3.f;
pressure=rho/3.f/porosity;
//-------------------- MRT collison where body force has NO higher-order terms -------------//
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);
@ -2512,7 +2512,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, i
}
//Calculate pressure for Incompressible-MRT model
pressure=rho/3.f;
pressure=rho/3.f/porosity;
//-------------------- IMRT collison where body force has NO higher-order terms -------------//
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);

15
gpu/GreyscaleColor.cu
View File

@ -6,7 +6,7 @@
//Model-1 & 4
__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, 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, int strideY, int strideZ, int start, int finish, int Np){
@ -512,6 +512,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Ma
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
Pressure[n] = rho/3.f/porosity;
//........................................................................
//..............carry out relaxation process..............................
@ -767,7 +768,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Ma
//Model-1 & 4
__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 *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, 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, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
@ -1217,6 +1218,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
Pressure[n] = rho/3.f/porosity;
//........................................................................
//..............carry out relaxation process..............................
@ -2918,14 +2920,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
//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 *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, 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, int strideY, int strideZ, int start, int finish, int Np){
//cudaProfilerStart();
//cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor, cudaFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel,
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel, Pressure,
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){
@ -2937,14 +2939,15 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
//Model-1 & 4
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 *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,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, int strideY, int strideZ, int start, int finish, int Np){
//cudaProfilerStart();
//cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor, cudaFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm,Vel,
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm,Vel,Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
cudaError_t err = cudaGetLastError();

756
models/GreyscaleColorModel.cpp
View File

@ -12,50 +12,13 @@ Two-fluid greyscale color lattice boltzmann model
ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauA_eff(0),tauB_eff(0),rhoA(0),rhoB(0),alpha(0),beta(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),GreyPorosity(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),greyMode(0),Lx(0),Ly(0),Lz(0),comm(COMM)
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
REVERSE_FLOW_DIRECTION = false;
}
ScaLBL_GreyscaleColorModel::~ScaLBL_GreyscaleColorModel(){
}
/*void ScaLBL_GreyscaleColorModel::WriteCheckpoint(const char *FILENAME, const double *cPhi, const double *cfq, int Np)
{
int q,n;
double value;
ofstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
value = cPhi[n];
File.write((char*) &value, sizeof(value));
// Write the even distributions
for (q=0; q<19; q++){
value = cfq[q*Np+n];
File.write((char*) &value, sizeof(value));
}
}
File.close();
}
void ScaLBL_GreyscaleColorModel::ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np)
{
int q=0, n=0;
double value=0;
ifstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
File.read((char*) &value, sizeof(value));
cPhi[n] = value;
// Read the distributions
for (q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value;
}
}
File.close();
}
*/
void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
@ -74,7 +37,6 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
Restart=false;
din=dout=1.0;
flux=0.0;
greyMode=false;
// Color Model parameters
if (greyscaleColor_db->keyExists( "timestepMax" )){
@ -117,12 +79,6 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
if (greyscaleColor_db->keyExists( "flux" )){
flux = greyscaleColor_db->getScalar<double>( "flux" );
}
if (greyscaleColor_db->keyExists("GreySolidLabels")&&
greyscaleColor_db->keyExists("GreySolidAffinity")&&
greyscaleColor_db->keyExists("PorosityList")&&
greyscaleColor_db->keyExists("PermeabilityList")){
greyMode = true;
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
@ -172,8 +128,7 @@ void ScaLBL_GreyscaleColorModel::SetDomain(){
N = Nx*Ny*Nz;
id = new signed char [N];
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
Averages = std::shared_ptr<SubPhase> ( new SubPhase(Dm) ); // TwoPhase analysis object
Averages = std::shared_ptr<GreyPhaseAnalysis> ( new GreyPhaseAnalysis(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
Dm->CommInit();
MPI_Barrier(comm);
@ -249,7 +204,6 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
if (rank == 0) cout << "Domain set." << endl;
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
}
void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
@ -662,13 +616,11 @@ void ScaLBL_GreyscaleColorModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
if (greyMode==true){
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
}
//ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
@ -708,10 +660,10 @@ void ScaLBL_GreyscaleColorModel::Create(){
// initialize phi based on PhaseLabel (include solid component labels)
AssignComponentLabels();//do open/black/grey nodes initialization
if (greyMode==true){
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
}
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
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
}
void ScaLBL_GreyscaleColorModel::Initialize(){
@ -799,7 +751,7 @@ void ScaLBL_GreyscaleColorModel::Initialize(){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
//ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
}
void ScaLBL_GreyscaleColorModel::Run(){
@ -869,15 +821,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
USE_SEED = true;
USE_MORPH = true;
}
else if (protocol == "open connected oil"){
morph_delta = -0.05;
USE_MORPH = true;
USE_MORPHOPEN_OIL = true;
}
else if (protocol == "shell aggregation"){
morph_delta = -0.05;
USE_MORPH = true;
}
if (greyscaleColor_db->keyExists( "capillary_number" )){
capillary_number = greyscaleColor_db->getScalar<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
@ -906,11 +850,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
morph_interval = analysis_db->getScalar<int>( "morph_interval" );
USE_MORPH = true;
}
if (analysis_db->keyExists( "use_morphopen_oil" )){
USE_MORPHOPEN_OIL = analysis_db->getScalar<bool>( "use_morphopen_oil" );
if (rank == 0 && USE_MORPHOPEN_OIL) printf("Volume change by morphological opening \n");
USE_MORPH = true;
}
if (analysis_db->keyExists( "tolerance" )){
tolerance = analysis_db->getScalar<double>( "tolerance" );
}
@ -946,20 +885,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
printf(" morph_delta = %f \n",morph_delta);
printf(" seed_water = %f \n",seed_water);
}
else if (protocol == "open connected oil"){
printf(" using protocol = open connected oil \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
printf(" tolerance = %f \n",tolerance);
printf(" morph_delta = %f \n",morph_delta);
}
else if (protocol == "shell aggregation"){
printf(" using protocol = shell aggregation \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
printf(" tolerance = %f \n",tolerance);
printf(" morph_delta = %f \n",morph_delta);
}
printf("No. of timesteps: %i \n", timestepMax);
fflush(stdout);
}
@ -976,7 +901,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
//std::shared_ptr<Database> analysis_db;
bool Regular = false;
auto current_db = db->cloneDatabase();
runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
//runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
//analysis.createThreads( analysis_method, 4 );
while (timestep < timestepMax ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
@ -999,20 +924,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,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);
////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);
}
else{
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 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_DeviceBarrier();
@ -1029,20 +948,15 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,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);
////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);
}
else{
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 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_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
@ -1063,20 +977,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_Regular->SendHalo(Phi);
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,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);
////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);
}
else{
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 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_DeviceBarrier();
@ -1093,47 +1001,39 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,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);
////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);
}
else{
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 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_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
//************************************************************************
PROFILE_STOP("Update");
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){
printf("%i %f \n",timestep,din);
printf("%i %.5g \n",timestep,din);
}
// Run the analysis
analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
if (timestep%analysis_interval == 0){
ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure);
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();
}
// allow initial ramp-up to get closer to steady state
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
analysis.finish();
//analysis.finish();
CURRENT_STEADY_TIMESTEPS += analysis_interval;
double volB = Averages->gwb.V;
double volA = Averages->gnb.V;
volA /= Dm->Volume;
volB /= Dm->Volume;;
//initial_volume = volA*Dm->Volume;
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
double vB_x = Averages->gwb.Px/Averages->gwb.M;
double vB_y = Averages->gwb.Py/Averages->gwb.M;
double vB_z = Averages->gwb.Pz/Averages->gwb.M;
double muA = rhoA*(tauA-0.5)/3.f;
double muB = rhoB*(tauB-0.5)/3.f;
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
@ -1147,10 +1047,12 @@ void ScaLBL_GreyscaleColorModel::Run(){
dir_z = 1.0;
force_mag = 1.0;
}
double current_saturation = volB/(volA+volB);
double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
double current_saturation = Averages->saturation;
double volA = current_saturation*GreyPorosity;
double volB = (1.0-current_saturation)*GreyPorosity;
double flow_rate_A = Averages->oil_flow_rate;
double flow_rate_B = Averages->water_flow_rate;
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(6.0*alpha);
if ( morph_timesteps > morph_interval ){
@ -1173,8 +1075,8 @@ void ScaLBL_GreyscaleColorModel::Run(){
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
if (rank == 0) printf(" -- adjust force by factor %.5g \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
greyscaleColor_db->putVector<double>("F",{Fx,Fy,Fz});
}
if ( isSteady ){
@ -1195,61 +1097,65 @@ void ScaLBL_GreyscaleColorModel::Run(){
slope_krA_volume = (log_krA - log_krA_prev)/(Dm->Volume*(volA - volA_prev));
delta_volume_target=min(delta_volume_target,Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume));
if (rank==0){
printf(" Enabling endpoint adaptation: krA = %f, krB = %f \n",krA_TMP,krB_TMP);
printf(" log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
printf(" Enabling endpoint adaptation: krA = %.5g, krB = %.5g \n",krA_TMP,krB_TMP);
printf(" log(kr)=%.5g, volume=%.5g, TARGET log(kr)=%.5g, volume change=%.5g \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
}
}
log_krA_prev = log_krA;
volA_prev = volA;
//******************************** **/
/** compute averages & write data **/
Averages->Full();
/*Averages->Full();
Averages->Write(timestep);
analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.finish();
*/
if (rank==0){
printf("** WRITE STEADY POINT *** ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
printf("Ca = %.5g, (previous = %.5g) \n",Ca,Ca_previous);
double h = Dm->voxel_length;
// pressures
double pA = Averages->gnb.p;
double pB = Averages->gwb.p;
double pAc = Averages->gnc.p;
double pBc = Averages->gwc.p;
double pA = Averages->Oil.p;
double pB = Averages->Water.p;
double pAB = (pA-pB)/(h*6.0*alpha);
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);
// -------- 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);
@ -1265,12 +1171,13 @@ void ScaLBL_GreyscaleColorModel::Run(){
WriteHeader=true;
kr_log_file = fopen("relperm.csv","a");
if (WriteHeader)
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water cap.pressure.norm pressure.drop Ca M\n");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",
CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility);
fclose(kr_log_file);
printf(" Measured capillary number %f \n ",Ca);
printf(" Measured capillary number %.5g \n ",Ca);
}
if (SET_CAPILLARY_NUMBER ){
Fx *= capillary_number / Ca;
@ -1281,8 +1188,8 @@ void ScaLBL_GreyscaleColorModel::Run(){
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
if (rank == 0) printf(" -- adjust force by factor %.5g \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
greyscaleColor_db->putVector<double>("F",{Fx,Fy,Fz});
}
@ -1291,7 +1198,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
else{
if (rank==0){
printf("** Continue to simulate steady *** \n ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
printf("Ca = %.5g, (previous = %.5g) \n",Ca,Ca_previous);
}
}
morph_timesteps=0;
@ -1319,17 +1226,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
double massChange = SeedPhaseField(seed_water);
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target);
}
else if (USE_MORPHOPEN_OIL){
delta_volume = volA*Dm->Volume - initial_volume;
if (rank==0) printf("***Morphological opening of connected oil, target volume change %f ***\n", delta_volume_target);
MorphOpenConnected(delta_volume_target);
}
else {
if (rank==0) printf("***Shell aggregation, target volume change %f ***\n", delta_volume_target);
//double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A
delta_volume += MorphInit(beta,delta_volume_target-delta_volume);
if (rank==0) printf("***Seed water in oil %.5g, volume change %.5g / %.5g ***\n", massChange, delta_volume, delta_volume_target);
}
if ( (delta_volume - delta_volume_target)/delta_volume_target > 0.0 ){
@ -1354,7 +1251,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
}
analysis.finish();
//analysis.finish();
PROFILE_STOP("Loop");
PROFILE_SAVE("lbpm_color_simulator",1);
//************************************************************************
@ -1377,190 +1274,50 @@ void ScaLBL_GreyscaleColorModel::Run(){
// ************************************************************************
}
double ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
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);
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=sumReduce( Dm->Comm, Count);
PoreCount=sumReduce( Dm->Comm, PoreCount);
if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
//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=sumReduce( Dm->Comm, Count);
// PoreCount=sumReduce( Dm->Comm, 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);
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
//ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
double saturation = Count/PoreCount;
return saturation;
//double saturation = Count/PoreCount;
//return saturation;
}
double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_change){
int nx = Nx;
int ny = Ny;
int nz = Nz;
int n;
int N = nx*ny*nz;
double volume_change=0.0;
if (target_volume_change < 0.0){
Array<char> id_solid(nx,ny,nz);
Array<int> phase_label(nx,ny,nz);
DoubleArray distance(Nx,Ny,Nz);
DoubleArray phase(nx,ny,nz);
signed char *id_connected;
id_connected = new signed char [nx*ny*nz];
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
// Extract only the connected part of NWP
BlobIDstruct new_index;
double vF=0.0; double vS=0.0;
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
MPI_Barrier(Dm->Comm);
long long count_connected=0;
long long count_porespace=0;
long long count_water=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++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( phase_label(i,j,k) == 0){
count_connected++;
}
if (id[n] > 0){
count_porespace++;
}
if (id[n] == 2){
count_water++;
}
}
}
}
count_connected=sumReduce( Dm->Comm, count_connected);
count_porespace=sumReduce( Dm->Comm, count_porespace);
count_water=sumReduce( Dm->Comm, count_water);
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( phase_label(i,j,k) == 0){
id_solid(i,j,k) = 1;
id_connected[n] = 2;
id[n] = 2;
/* delete the connected component */
phase(i,j,k) = -1.0;
}
else{
id_solid(i,j,k) = 0;
id_connected[n] = 0;
}
}
}
}
CalcDist(distance,id_solid,*Dm);
signed char water=2;
signed char notwater=1;
double SW=-(target_volume_change)/count_connected;
MorphOpen(distance, id_connected, Dm, SW, water, notwater);
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( id_connected[n] == 1){
id_solid(i,j,k) = 0;
}
else{
id_solid(i,j,k) = 1;
}
}
}
}
CalcDist(distance,id_solid,*Dm);
// re-initialize
double beta = 0.95;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
double d = distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){
if (d < 3.f){
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
int count_morphopen=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++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( id_connected[n] == 1){
count_morphopen++;
}
}
}
}
count_morphopen=sumReduce( Dm->Comm, count_morphopen);
volume_change = double(count_morphopen - count_connected);
if (rank==0) printf(" opening of connected oil %f \n",volume_change/count_connected);
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
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);
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
}
return(volume_change);
}
double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
@ -1626,7 +1383,7 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
count= sumReduce( Dm->Comm, count);
mass_loss= sumReduce( Dm->Comm, mass_loss);
if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
if (rank == 0) printf("Remove mass %.5g from %.5g voxels \n",mass_loss,count);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
@ -1636,219 +1393,6 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
return(mass_loss);
}
double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double target_delta_volume){
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
double vF = 0.f;
double vS = 0.f;
double delta_volume;
double WallFactor = 0.0;
bool USE_CONNECTED_NWP = false;
DoubleArray phase(Nx,Ny,Nz);
IntArray phase_label(Nx,Ny,Nz);;
DoubleArray phase_distance(Nx,Ny,Nz);
Array<char> phase_id(Nx,Ny,Nz);
fillHalo<double> fillDouble(Dm->Comm,Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
// Basic algorithm to
// 1. Copy phase field to CPU
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
double count = 0.f;
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 (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
double volume_initial = sumReduce( Dm->Comm, count);
/*
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb");
fwrite(phase.data(),8,N,INPUT);
fclose(INPUT);
*/
// 2. Identify connected components of phase field -> phase_label
double volume_connected = 0.0;
double second_biggest = 0.0;
if (USE_CONNECTED_NWP){
BlobIDstruct new_index;
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
MPI_Barrier(Dm->Comm);
// only operate on component "0"
count = 0.0;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int label = phase_label(i,j,k);
if (label == 0 ){
phase_id(i,j,k) = 0;
count += 1.0;
}
else
phase_id(i,j,k) = 1;
if (label == 1 ){
second_biggest += 1.0;
}
}
}
}
volume_connected = sumReduce( Dm->Comm, count);
second_biggest = sumReduce( Dm->Comm, second_biggest);
}
else {
// use the whole NWP
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (Averages->SDs(i,j,k) > 0.f){
if (phase(i,j,k) > 0.f ){
phase_id(i,j,k) = 0;
}
else {
phase_id(i,j,k) = 1;
}
}
else {
phase_id(i,j,k) = 1;
}
}
}
}
}
/*int reach_x, reach_y, reach_z;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
}
}
}*/
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*Dm);
double temp,value;
double factor=0.5/beta;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 3.f ){
value = phase(i,j,k);
if (value > 1.f) value=1.f;
if (value < -1.f) value=-1.f;
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
temp = -factor*log((1.0+value)/(1.0-value));
/// use this approximation close to the object
if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
// erase the original object
phase(i,j,k) = -1.0;
}
}
}
}
if (USE_CONNECTED_NWP){
if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
// if connected volume is less than 2% just delete the whole thing
if (rank==0) printf("Connected region has shrunk! \n");
REVERSE_FLOW_DIRECTION = true;
}
/* else{*/
if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
}
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
double target_delta_volume_incremental = target_delta_volume;
if (fabs(target_delta_volume) > 0.01*volume_initial)
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
else phase_id(i,j,k) = 1;
//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
}
}
}
CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
// 5. Update phase indicator field based on new distnace
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double d = phase_distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){
if (d < 3.f){
//phase(i,j,k) = -1.0;
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
fillDouble.fill(phase);
//}
count = 0.f;
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 (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f){
count+=1.f;
}
}
}
}
double volume_final= sumReduce( Dm->Comm, count);
delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
if (rank == 0) printf(" new saturation = %f \n", volume_final/(0.238323*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs)));
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
/*
sprintf(LocalRankFilename,"dist_final.%05i.raw",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(phase_distance.data(),8,N,DIST);
fclose(DIST);
sprintf(LocalRankFilename,"phi_final.%05i.raw",rank);
FILE *PHI = fopen(LocalRankFilename,"wb");
fwrite(phase.data(),8,N,PHI);
fclose(PHI);
*/
// 7. Re-initialize phase field and density
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);
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
return delta_volume;
}
void ScaLBL_GreyscaleColorModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);

9
models/GreyscaleColorModel.h
View File

@ -10,8 +10,7 @@ Implementation of two-fluid greyscale color lattice boltzmann model
#include <fstream>
#include "common/Communication.h"
#include "analysis/TwoPhase.h"
#include "analysis/runAnalysis.h"
#include "analysis/GreyPhase.h"
#include "common/MPI_Helpers.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
@ -40,7 +39,6 @@ public:
double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB;
double GreyPorosity;
bool greyMode;//run greyColor model if true
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@ -50,8 +48,7 @@ public:
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
//std::shared_ptr<TwoPhase> Averages;
std::shared_ptr<SubPhase> Averages;
std::shared_ptr<GreyPhaseAnalysis> Averages;
// input database
std::shared_ptr<Database> db;
@ -89,7 +86,7 @@ private:
void AssignComponentLabels();
void AssignGreySolidLabels();
void AssignGreyPoroPermLabels();
double ImageInit(std::string filename);
void ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);
double MorphOpenConnected(double target_volume_change);

64
models/GreyscaleModel.cpp
View File

@ -303,6 +303,56 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
}
}
void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity,double *Permeability,const vector<std::string> &File_poro,const vector<std::string> &File_perm)
{
double *Porosity_host, *Permeability_host;
Porosity_host = new double[N];
Permeability_host = new double[N];
double POROSITY=0.f;
double PERMEABILITY=0.f;
//Initialize a weighted porosity after considering grey voxels
double GreyPorosity_loc=0.0;
GreyPorosity=0.0;
//double label_count_loc = 0.0;
//double label_count_glb = 0.0;
Mask->ReadFromFile(File_poro[0],File_poro[1],Porosity_host);
Mask->ReadFromFile(File_perm[0],File_perm[1],Permeability_host);
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int idx = Map(i,j,k);
if (!(idx < 0)){
int n = k*Nx*Ny+j*Nx+i;
POROSITY = Porosity_host[n];
PERMEABILITY = Permeability_host[n];
if (POROSITY<=0.0){
ERROR("Error: Porosity for grey voxels must be 0.0 < Porosity <= 1.0 !\n");
}
else if (PERMEABILITY<=0.0){
ERROR("Error: Permeability for grey voxel must be > 0.0 ! \n");
}
else{
Porosity[idx] = POROSITY;
Permeability[idx] = PERMEABILITY;
GreyPorosity_loc += POROSITY;
//label_count_loc += 1.0;
}
}
}
}
}
GreyPorosity = sumReduce( Dm->Comm, GreyPorosity_loc);
GreyPorosity = GreyPorosity/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (rank==0){
printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
printf("The weighted porosity, considering both open and grey voxels, is %.3g\n",GreyPorosity);
}
delete [] Porosity_host;
delete [] Permeability_host;
}
void ScaLBL_GreyscaleModel::Create(){
/*
@ -356,7 +406,19 @@ void ScaLBL_GreyscaleModel::Create(){
double *Poros, *Perm;
Poros = new double[Np];
Perm = new double[Np];
AssignComponentLabels(Poros,Perm);
if (greyscale_db->keyExists("FileVoxelPorosityMap")){
//NOTE: FileVoxel**Map is a vector, including "file_name, datatype"
auto File_poro = greyscale_db->getVector<std::string>( "FileVoxelPorosityMap" );
auto File_perm = greyscale_db->getVector<std::string>( "FileVoxelPermeabilityMap" );
AssignComponentLabels(Poros,Perm,File_poro,File_perm);
}
else if (greyscale_db->keyExists("PorosityList")){
//initialize voxel porosity and perm from the input list
AssignComponentLabels(Poros,Perm);
}
else {
ERROR("Error: PorosityList or FilenameVoxelPorosityMap cannot be found! \n");
}
ScaLBL_CopyToDevice(Porosity, Poros, Np*sizeof(double));
ScaLBL_CopyToDevice(Permeability, Perm, Np*sizeof(double));
delete [] Poros;

2
models/GreyscaleModel.h
View File

@ -101,6 +101,6 @@ private:
char LocalRestartFile[40];
void AssignComponentLabels(double *Porosity, double *Permeablity);
void AssignComponentLabels(double *Porosity,double *Permeability,const vector<std::string> &File_poro,const vector<std::string> &File_perm);
};