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

merge master into FOM_dev

Browse Source
This commit is contained in:
Rex Zhe Li
2021-05-05 21:46:27 -04:00
parent cd49bae2b4 ec108ef1d1
commit 8e1a4a4195
23 changed files with 8870 additions and 856 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
analysis/SubPhase.cpp
View File

@@ -93,7 +93,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"sw krw krn vw vn pw pn wet\n");
fprintf(TIMELOG,"sw krw krn vw vn force pw pn wet\n");
}
}
}
@@ -348,7 +348,7 @@ void SubPhase::Basic(){
double krn = h*h*nu_n*not_water_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,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, gwb.p, gnb.p, total_wetting_interaction_global);
fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global);
fflush(TIMELOG);
}
if (err==true){

39
analysis/morphology.cpp
View File

@@ -320,6 +320,8 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
DoubleArray phase(nx,ny,nz);
IntArray phase_label(nx,ny,nz);
Array<char> ID(nx,ny,nz);
fillHalo<char> fillChar(Dm->Comm,Dm->rank_info,{nx-2,ny-2,nz-2},{1,1,1},0,1);
int n;
double final_void_fraction;
@@ -337,10 +339,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
count += 1.0;
id[n] = 2;
}
ID(i,j,k) = id[n];
}
}
}
fillChar.fill(ID);
Dm->Comm.barrier();
// total Global is the number of nodes in the pore-space
@@ -351,7 +354,8 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
if (rank==0) printf("Volume fraction for morphological opening: %f \n",volume_fraction);
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
// Communication buffers
/* // Communication buffers
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
@@ -400,7 +404,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
*/
int ii,jj,kk;
int Nx = nx;
int Ny = ny;
@@ -455,9 +459,10 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (ii=imin; ii<imax; ii++){
int nn = kk*nx*ny+jj*nx+ii;
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (id[nn] == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
if (ID(ii,jj,kk) == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
LocalNumber+=1.0;
id[nn]=1;
//id[nn]=1;
ID(ii,jj,kk)=1;
}
}
}
@@ -468,8 +473,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
}
}
}
fillChar.fill(ID);
// Pack and send the updated ID values
PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
/* PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
PackID(Dm->sendList("X"), Dm->sendCount("X") ,sendID_X, id);
PackID(Dm->sendList("y"), Dm->sendCount("y") ,sendID_y, id);
PackID(Dm->sendList("Y"), Dm->sendCount("Y") ,sendID_Y, id);
@@ -527,12 +533,12 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
UnpackID(Dm->recvList("YZ"), Dm->recvCount("YZ") ,recvID_YZ, id);
//......................................................................................
// double GlobalNumber = Dm->Comm.sumReduce( LocalNumber );
*/
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;
if (id[n] == 1){
if (ID(i,j,k) == 1){
phase(i,j,k) = 1.0;
}
else
@@ -550,9 +556,10 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
if (id[n] == 1 && phase_label(i,j,k) > 1){
id[n] = 2;
if (ID(i,j,k) == 1 && phase_label(i,j,k) > 1){
ID(i,j,k) = 2;
}
id[n] = ID(i,j,k);
}
}
}
@@ -702,12 +709,14 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
if (rank == 0) printf(" delta=%f, growth=%f, max. displacement = %f \n",morph_delta, GrowthEstimate, MAX_DISPLACEMENT);
// Now adjust morph_delta
double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious);
GrowthPrevious = GrowthEstimate;
morph_delta_previous = morph_delta;
morph_delta += step_size;
if (fabs(GrowthEstimate - GrowthPrevious) > 0.0) {
double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious);
GrowthPrevious = GrowthEstimate;
morph_delta_previous = morph_delta;
morph_delta += step_size;
}
if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0;
//MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25);
if (morph_delta > 0.0 ){

133
analysis/runAnalysis.cpp
View File

@@ -706,6 +706,139 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
}
// Initialize the comms
for ( int i = 0; i < 1024; i++ )
d_comm_used[i] = false;
// Initialize the threads
int N_threads = db->getWithDefault<int>( "N_threads", 4 );
auto method = db->getWithDefault<std::string>( "load_balance", "default" );
createThreads( method, N_threads );
}
runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
/* std::shared_ptr<Database> input_db, const RankInfoStruct &rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> Dm, int Np,
bool Regular, IntArray Map )
: d_Np( Np ),
d_regular( Regular ),
d_rank_info( rank_info ),
d_Map( Map ),
d_comm( Dm->Comm.dup() ),
d_ScaLBL_Comm( ScaLBL_Comm )*/
{
d_comm = ColorModel.Dm->Comm.dup();
d_Np = ColorModel.Np;
bool Regular = false;
auto input_db = ColorModel.db;
auto db = input_db->getDatabase( "Analysis" );
auto vis_db = input_db->getDatabase( "Visualization" );
// Ids of work items to use for dependencies
ThreadPool::thread_id_t d_wait_blobID;
ThreadPool::thread_id_t d_wait_analysis;
ThreadPool::thread_id_t d_wait_vis;
ThreadPool::thread_id_t d_wait_restart;
ThreadPool::thread_id_t d_wait_subphase;
char rankString[20];
sprintf( rankString, "%05d", ColorModel.Dm->rank() );
d_n[0] = ColorModel.Dm->Nx - 2;
d_n[1] = ColorModel.Dm->Ny - 2;
d_n[2] = ColorModel.Dm->Nz - 2;
d_N[0] = ColorModel.Dm->Nx;
d_N[1] = ColorModel.Dm->Ny;
d_N[2] = ColorModel.Dm->Nz;
d_restart_interval = db->getScalar<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_subphase_analysis_interval = INT_MAX;
d_visualization_interval = INT_MAX;
d_blobid_interval = INT_MAX;
if ( db->keyExists( "blobid_interval" ) ) {
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
}
if ( db->keyExists( "visualization_interval" ) ) {
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
}
if ( db->keyExists( "subphase_analysis_interval" ) ) {
d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
}
auto restart_file = db->getScalar<std::string>( "restart_file" );
d_restartFile = restart_file + "." + rankString;
d_rank = d_comm.getRank();
writeIDMap( ID_map_struct(), 0, id_map_filename );
// Initialize IO for silo
IO::initialize( "", "silo", "false" );
// Create the MeshDataStruct
d_meshData.resize( 1 );
d_meshData[0].meshName = "domain";
d_meshData[0].mesh = std::make_shared<IO::DomainMesh>(
d_rank_info, d_n[0], d_n[1], d_n[2], ColorModel.Dm->Lx, ColorModel.Dm->Ly, ColorModel.Dm->Lz );
auto PhaseVar = std::make_shared<IO::Variable>();
auto PressVar = std::make_shared<IO::Variable>();
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
auto BlobIDVar = std::make_shared<IO::Variable>();
if ( vis_db->getWithDefault<bool>( "save_phase_field", true ) ) {
PhaseVar->name = "phase";
PhaseVar->type = IO::VariableType::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( PhaseVar );
}
if ( vis_db->getWithDefault<bool>( "save_pressure", false ) ) {
PressVar->name = "Pressure";
PressVar->type = IO::VariableType::VolumeVariable;
PressVar->dim = 1;
PressVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( PressVar );
}
if ( vis_db->getWithDefault<bool>( "save_velocity", false ) ) {
VxVar->name = "Velocity_x";
VxVar->type = IO::VariableType::VolumeVariable;
VxVar->dim = 1;
VxVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( VxVar );
VyVar->name = "Velocity_y";
VyVar->type = IO::VariableType::VolumeVariable;
VyVar->dim = 1;
VyVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( VyVar );
VzVar->name = "Velocity_z";
VzVar->type = IO::VariableType::VolumeVariable;
VzVar->dim = 1;
VzVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( VzVar );
}
if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( SignDistVar );
}
if ( vis_db->getWithDefault<bool>( "save_connected_components", false ) ) {
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VariableType::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( BlobIDVar );
}
// Initialize the comms
for ( int i = 0; i < 1024; i++ )
d_comm_used[i] = false;

3
analysis/runAnalysis.h
View File

@@ -7,6 +7,7 @@
#include "common/Communication.h"
#include "common/ScaLBL.h"
#include "threadpool/thread_pool.h"
#include "models/ColorModel.h"
#include <limits.h>
@@ -31,6 +32,8 @@ public:
runAnalysis( std::shared_ptr<Database> db, const RankInfoStruct &rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> dm, int Np,
bool Regular, IntArray Map );
runAnalysis( ScaLBL_ColorModel &ColorModel);
//! Destructor
~runAnalysis();

17
common/Domain.cpp
View File

@@ -138,7 +138,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
// Fill remaining variables
N = Nx*Ny*Nz;
Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
Volume = nx*ny*nz*nproc[0]*nproc[1]*nproc[2]*1.0;
if (myrank==0) printf("voxel length = %f micron \n", voxel_length);
@@ -620,12 +620,16 @@ void Domain::Decomp( const std::string& Filename )
Comm.recv(id.data(),N,0,15);
}
Comm.barrier();
ComputePorosity();
delete [] SegData;
}
void Domain::ComputePorosity(){
// Compute the porosity
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocx()*nprocy()*nprocz());
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
for (int j=1;j<Ny-1;j++){
@@ -640,8 +644,8 @@ void Domain::Decomp( const std::string& Filename )
sum = Comm.sumReduce(sum_local);
porosity = sum*iVol_global;
if (rank()==0) printf("Media porosity = %f \n",porosity);
//.........................................................
delete [] SegData;
//.........................................................
}
void Domain::AggregateLabels( const std::string& filename ){
@@ -1450,4 +1454,3 @@ void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData
Comm.barrier();
}

1
common/Domain.h
View File

@@ -166,6 +166,7 @@ public: // Public variables (need to create accessors instead)
std::vector<signed char> id;
void ReadIDs();
void ComputePorosity();
void Decomp( const std::string& filename );
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();

13
common/ScaLBL.h
View File

@@ -87,6 +87,19 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
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_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, 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, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *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, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
//extern "C" void ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W,
// int start, int finish, int Np);
// ION TRANSPORT MODEL
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np);

1456
cpu/GreyscaleColor.cpp
View File

File diff suppressed because it is too large Load Diff

2
cpu/MixedGradient.cpp
View File

@@ -7,7 +7,7 @@ extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradie
{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
int i,j,k,n,N;
int i,j,k,n;
int np,np2,nm; // neighbors
double v,vp,vp2,vm; // values at neighbors
double grad;

1396
cuda/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

1542
cuda/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

16
example/Bubble/input.db
View File

@@ -12,6 +12,22 @@ Color {
ComponentAffinity = -1.0, 1.0, -1.0
}
FreeLee {
tauA = 1.0;
tauB = 1.0;
tauM = 1.0;//relaxation parameter for the phase field
rhoA = 1.0;
rhoB = 1.0;
gamma = 1.0e-4;//surface tension parameter in Lee model
W = 3.0; //theoretical interfacial thickness in Lee model; unit:[voxel]
F = 0, 0, 0
Restart = false
timestepMax = 1000
flux = 0.0
ComponentLabels = 0
ComponentAffinity = -1.0
}
Domain {
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 80, 80, 80 // Size of local domain (Nx,Ny,Nz)

2813
hip/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

1616
hip/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

128
models/ColorModel.cpp
View File

@@ -531,6 +531,121 @@ void ScaLBL_ColorModel::Initialize(){
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
}
double ScaLBL_ColorModel::Run(int returntime){
int nprocs=nprocx*nprocy*nprocz;
//************ MAIN ITERATION LOOP ***************************************/
comm.barrier();
PROFILE_START("Loop");
//std::shared_ptr<Database> analysis_db;
bool Regular = false;
auto current_db = db->cloneDatabase();
auto t1 = std::chrono::system_clock::now();
int START_TIMESTEP = timestep;
int EXIT_TIMESTEP = min(timestepMax,returntime);
while (timestep < EXIT_TIMESTEP ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
// *************ODD TIMESTEP*************
timestep++;
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
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);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set BCs
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
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);
ScaLBL_Comm->Barrier();
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
// Halo exchange for phase field
if (BoundaryCondition > 0 && BoundaryCondition < 5){
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_Regular->SendHalo(Phi);
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);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
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);
ScaLBL_Comm->Barrier();
//************************************************************************
}
PROFILE_STOP("Update");
PROFILE_STOP("Loop");
PROFILE_SAVE("lbpm_color_simulator",1);
//************************************************************************
// Compute the walltime per timestep
auto t2 = std::chrono::system_clock::now();
double cputime = std::chrono::duration<double>( t2 - t1 ).count() / (timestep - START_TIMESTEP);
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
return(MLUPS);
MLUPS *= nprocs;
}
void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
@@ -580,7 +695,6 @@ void ScaLBL_ColorModel::Run(){
if (color_db->keyExists( "krA_morph_factor" )){
KRA_MORPH_FACTOR = color_db->getScalar<double>( "krA_morph_factor" );
}
/* defaults for simulation protocols */
auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){
@@ -625,7 +739,7 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "seed_water" )){
seed_water = analysis_db->getScalar<double>( "seed_water" );
if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water);
USE_SEED = true;
ASSERT(protocol == "seed water");
}
if (analysis_db->keyExists( "morph_delta" )){
morph_delta = analysis_db->getScalar<double>( "morph_delta" );
@@ -656,7 +770,6 @@ void ScaLBL_ColorModel::Run(){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
@@ -1320,7 +1433,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
double vF = 0.f;
double vS = 0.f;
double delta_volume;
double WallFactor = 0.0;
double WallFactor = 1.0;
bool USE_CONNECTED_NWP = false;
DoubleArray phase(Nx,Ny,Nz);
@@ -1343,6 +1456,11 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
}
}
double volume_initial = Dm->Comm.sumReduce( count);
double PoreVolume = Dm->Volume*Dm->Porosity();
/*ensure target isn't an absurdly small fraction of pore volume */
if (volume_initial < target_delta_volume*PoreVolume){
volume_initial = target_delta_volume*PoreVolume;
}
/*
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb");
@@ -1492,7 +1610,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
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)));
if (rank == 0) printf(" new saturation = %f \n", volume_final/(Mask->Porosity()*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");

6
models/ColorModel.h
View File

@@ -16,6 +16,10 @@ Implementation of color lattice boltzmann model
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
#ifndef ScaLBL_ColorModel_INC
#define ScaLBL_ColorModel_INC
class ScaLBL_ColorModel{
public:
ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM);
@@ -29,6 +33,7 @@ public:
void Create();
void Initialize();
void Run();
double Run(int returntime);
void WriteDebug();
void getPhaseField(DoubleArray &f);
@@ -99,4 +104,5 @@ private:
int timestep;
int timestep_previous;
};
#endif

453
models/GreyscaleColorModel.cpp
View File

@@ -17,7 +17,7 @@ void DeleteArray( const TYPE *p )
ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& 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),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),GreyPorosity(0),RecoloringOff(0),
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;
@@ -43,6 +43,8 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
Restart=false;
din=dout=1.0;
flux=0.0;
RecoloringOff = false;
//W=1.0;
// Color Model parameters
if (greyscaleColor_db->keyExists( "timestepMax" )){
@@ -85,6 +87,9 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
if (greyscaleColor_db->keyExists( "flux" )){
flux = greyscaleColor_db->getScalar<double>( "flux" );
}
if (greyscaleColor_db->keyExists( "RecoloringOff" )){
RecoloringOff = greyscaleColor_db->getScalar<bool>( "RecoloringOff" );
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
@@ -212,6 +217,7 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
}
void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
{
// Initialize impermeability solid nodes and grey nodes
@@ -256,6 +262,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
label_count[idx] += 1.0;
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
@@ -290,19 +297,18 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
delete [] phase;
}
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalty wetting strength W
{
// 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];
// ranges [-1,1]
// water-wet > 0
// oil-wet < 0
// neutral = 0 (i.e. no penalty)
double *GreySolidW_host = new double [Np];
size_t NLABELS=0;
signed char VALUE=0;
@@ -324,117 +330,19 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
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);
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;
}
}
GreySolidW_host[idx] = AFFINITY;
}
}
}
}
if (rank==0){
printf("Number of Grey-solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
@@ -442,14 +350,13 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
AFFINITY=AffinityList[idx];
printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
}
printf("NOTE: grey-solid affinity>0: water-wet || grey-solid affinity<0: oil-wet \n");
}
ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] SolidPotential_host;
delete [] GreySolidGrad_host;
delete [] Dst;
delete [] GreySolidW_host;
}
////----------------------------------------------------------------------------------------------------------//
@@ -575,6 +482,71 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability;
}
//void ScaLBL_GreyscaleColorModel::AssignGreyscalePotential()
//{
// double *psi;//greyscale potential
// psi = new double[N];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<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(){
/*
* This function creates the variables needed to run a LBM
@@ -619,11 +591,13 @@ void ScaLBL_GreyscaleColorModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);
//ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);//greyscale potential
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &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 **) &Porosity_dvc, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
//...........................................................................
@@ -667,6 +641,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
AssignComponentLabels();//do open/black/grey nodes initialization
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
//AssignGreyscalePotential();
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);//porosity doesn't change over time
}
@@ -931,9 +906,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,0,ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@@ -943,10 +920,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,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);
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,
@@ -968,10 +949,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,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);
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,
@@ -983,9 +968,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,0,ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
@@ -995,10 +982,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,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);
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,
@@ -1020,10 +1011,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,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);
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,
@@ -1575,6 +1570,13 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
//ScaLBL_CopyToHost(PhaseField.data(), Psi, sizeof(double)*N);
//FILE *PSIFILE;
//sprintf(LocalRankFilename,"Psi.%05i.raw",rank);
//PSIFILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,PSIFILE);
//fclose(PSIFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank);
@@ -1631,26 +1633,26 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,PERM_FILE);
fclose(PERM_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[0],PhaseField);
FILE *GreySG_X_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_X.%05i.raw",rank);
GreySG_X_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_X_FILE);
fclose(GreySG_X_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[0],PhaseField);
//FILE *GreySG_X_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_X.%05i.raw",rank);
//GreySG_X_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_X_FILE);
//fclose(GreySG_X_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[Np],PhaseField);
FILE *GreySG_Y_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_Y.%05i.raw",rank);
GreySG_Y_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_Y_FILE);
fclose(GreySG_Y_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[Np],PhaseField);
//FILE *GreySG_Y_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_Y.%05i.raw",rank);
//GreySG_Y_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_Y_FILE);
//fclose(GreySG_Y_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[2*Np],PhaseField);
FILE *GreySG_Z_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_Z.%05i.raw",rank);
GreySG_Z_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_Z_FILE);
fclose(GreySG_Z_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[2*Np],PhaseField);
//FILE *GreySG_Z_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_Z.%05i.raw",rank);
//GreySG_Z_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_Z_FILE);
//fclose(GreySG_Z_FILE);
/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
FILE *CGX_FILE;
@@ -1984,6 +1986,169 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
// 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()

9
models/GreyscaleColorModel.h
View File

@@ -39,6 +39,8 @@ public:
double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB;
double GreyPorosity;
bool RecoloringOff;//recoloring can be turn off for grey nodes if this is true
//double W;//wetting strength paramter for capillary pressure penalty for grey nodes
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@@ -48,7 +50,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<GreyPhaseAnalysis> Averages;
std::shared_ptr<GreyPhaseAnalysis> Averages;
// input database
std::shared_ptr<Database> db;
@@ -64,12 +66,14 @@ public:
double *fq, *Aq, *Bq;
double *Den, *Phi;
//double *GreySolidPhi; //Model 2 & 3
double *GreySolidGrad;//Model 1 & 4
//double *GreySolidGrad;//Model 1 & 4
double *GreySolidW;
//double *ColorGrad;
double *Velocity;
double *Pressure;
double *Porosity_dvc;
double *Permeability_dvc;
//double *Psi;
private:
Utilities::MPI comm;
@@ -86,6 +90,7 @@ private:
void AssignComponentLabels();
void AssignGreySolidLabels();
void AssignGreyPoroPermLabels();
//void AssignGreyscalePotential();
void ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);

34
sample_scripts/configure_lonestar Executable file
View File

@@ -0,0 +1,34 @@
# load the module for cmake
#module load cmake
#source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
module load cmake gcc
module load cuda
export HDF5_DIR=$HOME/local/hdf5/1.8.12/
export SILO_DIR=$HOME/local/silo/4.10.2/
export NETCDF_DIR=$HOME/local/netcdf/4.6.1
# configure
rm -rf CMake*
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_CXX_STANDARD=14 \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_70 -Xptxas=-v -Xptxas -dlcm=cg -lineinfo" \
-D CMAKE_CUDA_HOST_COMPILER="/opt/apps/gcc/7.3.0/bin/gcc" \
-D USE_HDF5=1 \
-D HDF5_DIRECTORY="$HDF5_DIR" \
-D HDF5_LIB="$HDF5_DIR/lib/libhdf5.a" \
-D USE_SILO=1 \
-D SILO_LIB="$SILO_DIR/lib/libsiloh5.a" \
-D SILO_DIRECTORY="$SILO_DIR" \
-D USE_NETCDF=0 \
-D NETCDF_DIRECTORY="$NETCDF_DIR" \
-D USE_DOXYGEN:BOOL=false \
-D USE_TIMER=0 \
~/src/LBPM
make VERBOSE=1 -j1 && make install

4
sample_scripts/configure_summit
View File

@@ -4,7 +4,7 @@
#source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
module load cmake gcc
module load cuda
#/ccs/proj/csc380/mcclurej
export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1
@@ -28,7 +28,7 @@ cmake \
-D USE_SILO=1 \
-D SILO_LIB="$SILO_DIR/lib/libsiloh5.a" \
-D SILO_DIRECTORY="$SILO_DIR" \
-D USE_NETCDF=1 \
-D USE_NETCDF=0 \
-D NETCDF_DIRECTORY="$NETCDF_DIR" \
-D USE_DOXYGEN:BOOL=false \
-D USE_TIMER=0 \

2
tests/CMakeLists.txt
View File

@@ -62,7 +62,7 @@ ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestWideHalo )
ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db)
ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
#ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
#ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db)

35
tests/lbpm_color_simulator.cpp
View File

@@ -27,19 +27,24 @@ int main( int argc, char **argv )
// Initialize
Utilities::startup( argc, argv );
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
std::string SimulationMode = "production";
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
if (argc > 2) {
SimulationMode = "development";
}
if ( rank == 0 ) {
printf( "********************************************************\n" );
printf( "Running Color LBM \n" );
printf( "********************************************************\n" );
if (SimulationMode == "development")
printf("**** DEVELOPMENT MODE ENABLED *************\n");
}
// Initialize compute device
int device = ScaLBL_SetDevice( rank );
@@ -62,8 +67,28 @@ int main( int argc, char **argv )
ColorModel.Create(); // creating the model will create data structure to match the pore
// structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
ColorModel.Run();
// ColorModel.WriteDebug();
if (SimulationMode == "development"){
double MLUPS=0.0;
int timestep = 0;
int analysis_interval = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "" )){
analysis_interval = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
}
FlowAdaptor Adapt(ColorModel);
runAnalysis analysis(ColorModel);
while (ColorModel.timestep < ColorModel.timestepMax){
timestep += analysis_interval;
MLUPS = ColorModel.Run(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
Adapt.MoveInterface(ColorModel);
}
ColorModel.WriteDebug();
} //Analysis.WriteVis(LeeModel,LeeModel.db, timestep);
else
ColorModel.Run();
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );

4
tests/lbpm_freelee_simulator.cpp
View File

@@ -62,8 +62,8 @@ int main( int argc, char **argv )
double MLUPS=0.0;
int timestep = 0;
int visualization_time = LeeModel.timestepMax;
if (LeeModel.vis_db->keyExists( "visualizataion_interval" )){
visualization_time = LeeModel.vis_db->getScalar<int>( "visualizataion_interval" );
if (LeeModel.vis_db->keyExists( "visualization_interval" )){
visualization_time = LeeModel.vis_db->getScalar<int>( "visualization_interval" );
timestep += visualization_time;
}
while (LeeModel.timestep < LeeModel.timestepMax){