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 2021-04-04 21:02:42 -04:00
commit fde6af1dc0
14 changed files with 9021 additions and 272 deletions
Show Stats 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){

30
common/ScaLBL.h
View File

@ -103,6 +103,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 *Psi, 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, bool RecoloringOff, double W, 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 *Psi, 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, bool RecoloringOff, double W, 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);
@ -209,11 +222,12 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double
extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad,
double rhonA, double rhoB, double tauM, double W, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
// double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi,
// double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np);
@ -230,6 +244,14 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np);

1276
cpu/FreeLee.cpp
View File

File diff suppressed because it is too large Load Diff

1456
cpu/GreyscaleColor.cpp
View File

File diff suppressed because it is too large Load Diff

1590
cuda/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

1614
cuda/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

1480
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

74
models/FreeLeeModel.cpp
View File

@ -77,8 +77,10 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
Fx = Fy = Fz = 0.0;
gamma=1e-3;//surface tension
W=5.0;//interfacial thickness
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
//beta = 12.0*gamma/W;
//kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
beta = 0.75*gamma/W;
kappa = 0.375*gamma*W;//beta and kappa are related to surface tension \gamma
Restart=false;
din=dout=1.0;
flux=0.0;
@ -136,8 +138,10 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
outletA=0.f;
outletB=1.f;
//update secondary parameters
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
//beta = 12.0*gamma/W;
//kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
beta = 0.75*gamma/W;
kappa = 0.375*gamma*W;//beta and kappa are related to surface tension \gamma
//if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
BoundaryCondition = 0;
@ -795,27 +799,31 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
// Compute the Phase indicator field
// Read for hq happens in this routine (requires communication)
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
// Halo exchange for phase field
ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi);
//ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){
//TODO to be revised
// Need to add BC for hq!!!
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set BCs
@ -832,7 +840,9 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar);
}
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
//ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
@ -841,24 +851,29 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMA
ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
// Halo exchange for phase field
ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi);
//ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM 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);
}
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set boundary conditions
@ -874,7 +889,9 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar);
}
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
//ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
//************************************************************************
@ -1055,6 +1072,13 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
fwrite(PhaseField.data(),8,N,PFILE);
fclose(PFILE);
ScaLBL_Comm->RegularLayout(Map,mu_phi,PhaseField);
FILE *CHEMFILE;
sprintf(LocalRankFilename,"ChemPotential.%05i.raw",rank);
CHEMFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CHEMFILE);
fclose(CHEMFILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);

135
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),W(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,12 @@ 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" );
}
if (greyscaleColor_db->keyExists( "W" )){
W = greyscaleColor_db->getScalar<double>( "W" );
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
@ -212,6 +220,7 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
}
void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
{
// Initialize impermeability solid nodes and grey nodes
@ -256,6 +265,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
@ -575,6 +585,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
@ -595,6 +670,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
ScaLBL_Comm_Regular = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
ScaLBL_Comm_Regular_2 = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N);
@ -619,6 +695,7 @@ 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);
@ -667,6 +744,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 +1009,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,14 +1023,20 @@ 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,
ScaLBL_Comm_Regular_2->SendHalo(Psi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Psi, 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);
alpha, beta, Fx, Fy, Fz, RecoloringOff, W, 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_2->RecvHalo(Psi);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
@ -968,10 +1054,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, Psi, 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);
alpha, beta, Fx, Fy, Fz, RecoloringOff, W, 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 +1073,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,14 +1087,20 @@ 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,
ScaLBL_Comm_Regular_2->SendHalo(Psi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, Psi, 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);
alpha, beta, Fx, Fy, Fz, RecoloringOff, W, 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_2->RecvHalo(Psi);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
@ -1020,10 +1118,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, Psi, 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);
alpha, beta, Fx, Fy, Fz, RecoloringOff, W, 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 +1677,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);

7
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,8 @@ 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<ScaLBL_Communicator> ScaLBL_Comm_Regular_2;
std::shared_ptr<GreyPhaseAnalysis> Averages;
// input database
std::shared_ptr<Database> db;
@ -70,6 +73,7 @@ public:
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);

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 \

5
tests/lbpm_color_simulator.cpp
View File

@ -84,12 +84,11 @@ int main( int argc, char **argv )
Adapt.MoveInterface(ColorModel);
}
ColorModel.WriteDebug();
} //Analysis.WriteVis(LeeModel,LeeModel.db, timestep);
else
ColorModel.Run();
ColorModel.WriteDebug();
ColorModel.Run();
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );

2
tests/lbpm_freelee_simulator.cpp
View File

@ -72,7 +72,7 @@ int main( int argc, char **argv )
Analysis.WriteVis(LeeModel,LeeModel.db, timestep);
timestep += visualization_time;
}
//LeeModel.WriteDebug_TwoFluid();
LeeModel.WriteDebug_TwoFluid();
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;